ConstantThrustManeuver implementation

[daiana]

Hello!!

maneuverStartDate = maneuverStartDates[0]
maneuverStopDate = maneuverStopDates[-1]
generator = propagator_num.getEphemerisGenerator()
propagator_num.propagate(maneuverStartDate, maneuverStopDate)
ephemeris = generator.getGeneratedEphemeris()

lvlh = LocalOrbitalFrame(FramesFactory.getGCRF(), LOFType.LVLH, ephemeris, 'LVLH')
#LVLH X axis aligned with position, Z axis aligned with orbital momentum
#GCRF The X-axis is aligned with the mean equinox of Earth at 12:00 Terrestrial Time on the 1st of January, 2000, and the Z-axis is aligned with the Earth´s rotation axis.

propagator_num = NumericalPropagator(integrator)
for i in range(len(maneuverStartDates)): 
    maneuverStartDate = maneuverStartDates[i]
    maneuverDuration = maneuverDurations[i]
    maneuverStopDate = maneuverStopDates[i]
    delta_Ax = delta_Axs[i]
    delta_Ay = delta_Ays[i]
    delta_Az = delta_Azs[i]
    delta_Vx = delta_Ax*maneuverDuration
    delta_Vy = delta_Ay*maneuverDuration
    delta_Vz = delta_Az*maneuverDuration
    delta_V_abs = (delta_Vx**2+delta_Vy**2+delta_Vz**2)**0.5
    delta_A = (delta_Ax**2+delta_Ay**2+delta_Az**2)**0.5
    
    delta_V = Vector3D(float(delta_Vx), float(delta_Vy), float(delta_Vz))
    #the acceleration direction in satellite frame.
    delta_V_lvlh = inertialFrame.getStaticTransformTo(lvlh, maneuverStartDate).transformVector(delta_V)
    
    specificImpulse = 4400.0   #sec (depends on type of propellant used for maneuver)   #ion propulsion   #efficiency 
    delta_m = -1700*(1-np.exp(delta_V_abs/specificImpulse/9.80665))   #g0 = 9.80665 - sea level standard acceleration of gravity
    thrustForce = delta_A*1700.0
    maneuverDirection = Vector3D(float(delta_V_lvlh.getX()/delta_V_abs), float(delta_V_lvlh.getY()/delta_V_abs), float(delta_V_lvlh.getZ()/delta_V_abs))
    print('maneuver data', thrustForce, delta_A, delta_V, delta_m, maneuverDirection)
    propagator_num.addForceModel(ConstantThrustManeuver(maneuverStartDate, float(maneuverDuration), 
                                                        float(thrustForce), float(specificImpulse), maneuverDirection))

step_time = 1800
duration = 86400 * 6.0 + step_time

propagator_num.setOrbitType(OrbitType.CARTESIAN)
propagator_num.setInitialState(initialState)
propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(ITRF, gravityProvider))
propagator_num.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getSun()))
propagator_num.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getMoon()))
propagator_num.addForceModel(SolarRadiationPressure(sun, earth, isotropicRadiationSingleCoeff))

t = [initialDate.shiftedBy(float(dt)) \
        for dt in np.arange(0, duration, step_time)]
pos1 = [propagator_num.propagate(tt).getPVCoordinates().getPosition() for tt in t]
vel1 = [propagator_num.propagate(tt).getPVCoordinates().getVelocity() for tt in t]

I have input data for maneuvers in GCRF frame, so i need to translate it to spacecraft frame (Do i need to use specific one??). To do so, it’s necessary to use getEphemerisGenerator() and propagate the orbit. I wanted to know do I have to create new propagator and add all of the disturbances again, i.e. to make it twice (the first propagator without including maneuvers and second one with maneuvers)??

[MaximeJ]

Hi @daiana,

@daiana

I have input data for maneuvers in GCRF frame, so i need to translate it to spacecraft frame (Do i need to use specific one??).

It’s up to you which Local Orbital Frame (LOF) you want to use
I tend to prefer LOFType.TNW since the X axis is aligned with velocity and well adapted for tangential thrust.
But you will find usage of LOFType.QSW a lot in the literature.

@daiana

I wanted to know do I have to create new propagator and add all of the disturbances again, i.e. to make it twice (the first propagator without including maneuvers and second one with maneuvers)??

You could use a NumericalPropagatorBuilder that was made exactly for this (generate a new propagator using the method buildPropagator).
Now if you just want to propagate the first one then the second one you can just:

• Propagate a first time without maneuvers
• Add the maneuvers’ force models
• Re-initialize the state using resetInitialState
• Propagate a second time with maneuvers

Hope this helps,
Maxime

[daiana]

Thank you, Maxime!!

[daiana]

hello again!! I’ve realized something. I can’t use resetInitialState() because the calculated state is not the state on initial date. I mean we have to take all of the main perturbations (earths oblateness, third body attraction, srp) into account for whole time span. To translate the deltaV vector to LVLH frame we have to do it on the time duration that covers moments of maneuvers, otherwise the error will occur. Thus we can’t use another propagator using buildPropagator() method to add maneuvers as perturbations. With using reselinitialState() the divergences at the beggining will raise, because propagating the second time we use not right state. It only means that we can use propagatorBuilder to propagate for the first time without maneuvers and then create a new one and add all of the forces once again including ConstantThrustManeuver this time and propagate the second time to avoid this divergences. I’m adding the graph to make it more understandable. If i’m not mistaken the green graph is the case without maneuvers at all (propagating just once), orange and blue graphs are those ones that i’ve got using resetInitialState() without and with maneuvers correspondingly. Green vertical lines are marking the starts of the maneuvers (but that doesn’t matter).

Translating maneuverDirection vector does not make any difference in results