From 040880ec6e4bf4738a28a175e771a7081004a7f4 Mon Sep 17 00:00:00 2001 From: Jorge Martinez Date: Thu, 29 Jan 2026 11:40:17 +0100 Subject: [PATCH 01/18] doc: add geostationary phasing example --- examples/phasing-geo-orbits.py | 336 +++++++++++++++++++++++++++++++++ 1 file changed, 336 insertions(+) create mode 100644 examples/phasing-geo-orbits.py diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py new file mode 100644 index 0000000000..138eb287ff --- /dev/null +++ b/examples/phasing-geo-orbits.py @@ -0,0 +1,336 @@ +# # GEO Orbit Phasing +# +# Once we get the intended orbit, a change of the true anomaly is sometime needed to target a well-defined orbit slot. This could be the case of constellations that have more than one satellite in the same orbital plane and a certain relative geometry has to be respected or a geostationary satellite that is moved from its original location to another one for operational purposes. +# +# In both cases we can solve for a reference (or target) satellite in Astrogator, whose position over time is taken as reference by a maneuvering (or chaser) satellite using a targeted MCS. +# +# ## STK Scenario +# +# This lesson focuses on the design of phasing maneuver to be executed by a geostationary satellite who needs to change its position along the geostationary belt. +# +#
Drawing
+ +# ## Import Libraries + +# + +import math + +from ansys.stk.core.stkdesktop import STKDesktop +from ansys.stk.core.stkobjects import * +from ansys.stk.core.stkobjects.astrogator import * +from ansys.stk.core.utilities.colors import Colors + + +# - + +# ### Create a New Scenario and Configure Satellites +# +# We now create a new scenario containing 2 satellites: the *Target* and the *Chaser*. Their initial states are not defined yet, but they will be on the geostationary orbit as specified in the next code cell. + +# + +stk = STKDesktop.start_application() +app = STKDesktop.start_application(visible=True) +root = app.root +root.new_scenario("Phasing_orbit") +scen = Scenario(root.current_scenario) + +# configure the target satellite +targetSat = Satellite(scen.children.new(STKObjectType.SATELLITE, "Target")) +targetSat.set_propagator_type(PropagatorType.ASTROGATOR) +tDriver = MCSDriver(targetSat.propagator) +tDriver.options.draw_trajectory_in_3d = False +tInitState = MCSInitialState(tDriver.main_sequence["Initial State"]) +tInitState.set_element_type(ElementSetType.KEPLERIAN) +tInitState.initial_state.epoch = scen.start_time +tPropagator = MCSPropagate(tDriver.main_sequence["Propagate"]) +tPropagator.propagator_name = "Earth point mass" +tPropagator.properties.color = Colors.Red +tKep = tInitState.element + +# configure the chaser satellite +chaserSat = Satellite(scen.children.new(STKObjectType.SATELLITE, "Chaser")) +chaserSat.set_propagator_type(PropagatorType.ASTROGATOR) +cDriver = MCSDriver(chaserSat.propagator) +cDriver.options.draw_trajectory_in_3d = False +cInitState = MCSInitialState(cDriver.main_sequence["Initial State"]) +cInitState.set_element_type(ElementSetType.KEPLERIAN) +cInitState.initial_state.epoch = scen.start_time +cPropagator = MCSPropagate(cDriver.main_sequence["Propagate"]) +cPropagator.propagator_name = "Earth point mass" +cPropagator.stopping_conditions.add("Periapsis") +cPropagator.stopping_conditions.remove("Duration") +cKep = cInitState.element + +root.execute_command("VO */Satellite/Chaser Pass3D OrbitLead None") +root.execute_command("VO */Satellite/Chaser Pass3D OrbitTrail All") +root.execute_command("VO */Satellite/Target Pass3D OrbitLead None") +root.execute_command("VO */Satellite/Target Pass3D OrbitTrail All") +# - + +# ### Define the Input Data and Run the Analysis +# +# The cell below defines and runs the MCS for both satellites. +# +# The **phase angle** is the angular displacement between the two satellites at initial time. In this context, it is defined as always positive and between 0 and 360 degrees. To make the rendezvous between the satellites, it has to be reduced to 0 in a specified number of orbits (called **phasing orbits**), so each phase orbit will reduce the angular gap by just a fraction of the overall value. +# +# #### $$ \theta_{gap}=\frac{\theta}{n_{orbits}}$$ +# +# , where $\theta$ is the phase angle, and $\theta_{gap}$ is the angular dispacement to recover after each phasing orbit. +# +# In the code below the **Kepler's equation** is used to calculate the period and semimajor axis of the phasing orbit, given the number of revolution around the phasing orbit itself. As first, the **eccentric anomaly E** is calculated from $\theta_{gap}$ : +# +# #### $$\tan\frac{E}{2} = \sqrt{\frac{1 - e}{1 + e}}\tan\frac{\theta_{gap}}{2}$$ +# +# ...and then the period of the phasing orbit is derived: +# +# #### $$ T_{phasing}=\frac{T_{chaser}}{2\pi}\left ( E-e \sin E \right )$$ +# +# Having those data available, a differential corrector is configured to let the two satellites be in the same position at the end of the transfer + +# + +########################################## INPUT DATA ########################################################## +true_anom = 70 # deg - chaser true anomaly (between 0 and 360) +target_ta = 20 # deg - target true anomaly (between 0 and 360) +nOrbits = 5 # number of phasing orbits +prop_time = 16 # days +################################################################################################################ + +try: + root.execute_command( + 'VectorTool * CentralBody/Earth Create Angle Phase "Dihedral Angle"' + ) + root.execute_command( + 'VectorTool * CentralBody/Earth Modify Angle Phase "Dihedral Angle" "Satellite/Chaser Position" "Satellite/Target Position" "Satellite/Target Orbit_Normal" 0-360 Positive' + ) + root.execute_command( + 'VO * SetVectorGeometry Add "CentralBody/Earth Phase Angle" Color #fcba03 Show On' + ) + root.execute_command( + 'VO * SetVectorGeometry Add "Satellite/Target Position Vector"' + ) + root.execute_command( + 'VO * SetVectorGeometry Modify "Satellite/Target Position Vector" Color red Thickness 5 Show On' + ) + root.execute_command( + 'VO * SetVectorGeometry Add "Satellite/Chaser Position Vector"' + ) + root.execute_command( + 'VO * SetVectorGeometry Modify "Satellite/Chaser Position Vector" Color green Thickness 5 Show On' + ) +except: + print("") + + +# fixed keplerian parameters (geostationary orbit) +sma = 42164 +ecc = 0.0001 +inc = 0.0 +raan = 0.0 + +# change the scenario duration +scen.stop_time = "+" + str(prop_time) + " day" +try: + cDriver.main_sequence.remove("Start Phasing") + cDriver.main_sequence.remove("Circularization") +except: + print("") +finally: + # set the chaser initial state + cKep.semimajor_axis = sma + cKep.eccentricity = ecc + cKep.inclination = inc + cKep.raan = raan + cKep.true_anomaly = true_anom + cKep.arg_of_periapsis = 0.0 + + # set the target initial state + tKep.semimajor_axis = sma + tKep.eccentricity = ecc + tKep.inclination = inc + tKep.raan = raan + tKep.true_anomaly = target_ta + tKep.arg_of_periapsis = 0.0 + + tPropagator.stopping_conditions["Duration"].properties.trip = prop_time * 86400 + + cDriver.run_mcs() + tDriver.run_mcs() + + # get the mean period and SMA of the initial orbit + meanKepDp = chaserSat.data_providers["Kozai-Izsak Mean"].group["ICRF"] + meanKep = meanKepDp.execute(scen.start_time, scen.start_time, 60) + meanKepDf = meanKep.data_sets.to_pandas_dataframe() + meanPeriod = float(meanKepDf.at[0, "mean nodal period"]) + meanSma = float(meanKepDf.at[0, "mean semi-major axis"]) + + # calculate the phase angle + phaseAngle = math.radians(target_ta - true_anom) + + if phaseAngle < 0: + phaseAngle = 2 * math.pi + phaseAngle + + print("#################### Initial geometry #####################") + print("Phase angle = " + str(math.degrees(phaseAngle)) + " deg") + print("") + + # recalculate the phase angle accordingly with the number of phasing orbits + phaseAngle = phaseAngle / nOrbits + + # calculate the eccentric anomaly + E = 2 * math.atan(math.sqrt((1 - ecc) / (1 + ecc)) * math.tan(phaseAngle / 2)) + + # calculate the time elapsed to cover phase angle in original orbit (Kepler's equation) + tPhaseAngle = (meanPeriod / (2 * math.pi)) * (E - ecc * math.sin(E)) + + # calculate the period of the phasing orbit + tPhasingOrbit = meanPeriod - tPhaseAngle + + # calculate the SMA of the phasing orbit + mu = 398600.44 + transferSma = math.pow((tPhasingOrbit * math.sqrt(mu)) / (2 * math.pi), 2 / 3) + + ## calculate eccentricity and radii + kepDp = chaserSat.data_providers["Astrogator Values"].group["Keplerian Elems"] + kep = kepDp.execute(scen.start_time, scen.start_time, 60) + kepDf = kep.data_sets.to_pandas_dataframe() + initRadius = float(kepDf.at[0, "radius_of_periapsis"]) + + if transferSma > sma: + periRadius = initRadius + apoRadius = 2 * transferSma - periRadius + else: + apoRadius = initRadius + periRadius = 2 * transferSma - apoRadius + + # estimate the needed Delta V for phasing + vPeriInitial = math.sqrt(mu * ((2 / initRadius) - (1 / meanSma))) + vPeriFinal = math.sqrt(mu * ((2 / initRadius) - (1 / transferSma))) + deltaV = vPeriFinal - vPeriInitial + + # add the first Target Sequence segment + ts1 = MCSTargetSequence( + cDriver.main_sequence.insert(SegmentType.TARGET_SEQUENCE, "Start Phasing", "-") + ) + ts1.action = TargetSequenceAction.RUN_ACTIVE_PROFILES + + ## add a Maneuver segment + dv1 = MCSManeuver(ts1.segments.insert(SegmentType.MANEUVER, "DV1", "-")) + dv1.set_maneuver_type(ManeuverType.IMPULSIVE) + dv1.enable_control_parameter( + ControlManeuver.IMPULSIVE_CARTESIAN_X + ) # control parameter + + #### maneuver attitude definition + dv1.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) + thrustVector = dv1.maneuver.attitude_control + thrustVector.thrust_axes_name = "Satellite/Chaser VNC(Earth)" + thrustVector.x = deltaV * 1000 + + ## add a Propagate segment + phasingOrbit = MCSPropagate( + ts1.segments.insert(SegmentType.PROPAGATE, "Phasing Orbit", "-") + ) + phasingOrbit.properties.color = Colors.Orange + phasingOrbit.propagator_name = "Earth point mass" + if transferSma > sma: + phasingOrbit.stopping_conditions.add("Periapsis") + else: + phasingOrbit.stopping_conditions.add("Apoapsis") + phasingOrbit.stopping_conditions.remove("Duration") + phasingOrbit.stopping_conditions.item(0).properties.repeat_count = nOrbits + phasingOrbit.results.add("Vector/Angle Between Vectors") + phasingOrbit.results[0].vector1_name = "Satellite/Chaser Position" + phasingOrbit.results[0].vector2_name = "Satellite/Target Position" + + # customize the Differential Corrector + dc1 = ProfileDifferentialCorrector(ts1.profiles["Differential Corrector"]) + dc1.mode = ProfileMode.ITERATE + dc1.max_iterations = 50 + + # set Control Parameters and Results + xControlParam1 = dc1.control_parameters.get_control_by_paths( + "DV1", "ImpulsiveMnvr.Cartesian.X" + ) + xControlParam1.enable = True + xControlParam1.max_step = 0.001 + + roaResult = dc1.results.get_result_by_paths( + "Phasing Orbit", "Angle_Between_Vectors" + ) + roaResult.enable = True + roaResult.desired_value = 0.0 + roaResult.tolerance = 0.1 + + # add the second Target Sequence segment + ts2 = MCSTargetSequence( + cDriver.main_sequence.insert( + SegmentType.TARGET_SEQUENCE, "Circularization", "-" + ) + ) + ts2.action = TargetSequenceAction.RUN_ACTIVE_PROFILES + + ## add a Maneuver segment + dv2 = MCSManeuver(ts2.segments.insert(SegmentType.MANEUVER, "DV2", "-")) + dv2.set_maneuver_type(ManeuverType.IMPULSIVE) + dv2.enable_control_parameter( + ControlManeuver.IMPULSIVE_CARTESIAN_X + ) # control parameter + + ## add a Propagate segment + finalOrbit = MCSPropagate( + ts2.segments.insert(SegmentType.PROPAGATE, "Final Orbit", "-") + ) + finalOrbit.properties.color = Colors.Yellow + finalOrbit.propagator_name = "Earth point mass" + finalOrbit.stopping_conditions["Duration"].properties.trip = 86400 + finalOrbit.results.add("Keplerian Elems/Eccentricity") + + # customize the Differential Corrector + dc2 = ProfileDifferentialCorrector(ts2.profiles["Differential Corrector"]) + dc2.mode = ProfileMode.ITERATE + dc2.max_iterations = 50 + + # set Control Parameters and Results + xControlParam2 = dc2.control_parameters.get_control_by_paths( + "DV2", "ImpulsiveMnvr.Cartesian.X" + ) + xControlParam2.enable = True + xControlParam2.max_step = 0.01 + + roaResult = dc2.results.get_result_by_paths("Final Orbit", "Eccentricity") + roaResult.enable = True + roaResult.desired_value = 0.0001 + roaResult.tolerance = 0.00001 + cDriver.run_mcs() + tDriver.run_mcs() + root.rewind() + + # get the maneuver data providers + manDp = chaserSat.data_providers["Maneuver Summary"] + man = manDp.execute(scen.start_time, scen.stop_time) + manDf = man.data_sets.to_pandas_dataframe() + deltaV = manDf.at[0, "delta v"] + + # change orbit system + try: + root.execute_command( + 'VO */Satellite/Chaser OrbitSystem Modify System "InertialByWindow" Show Off' + ) + root.execute_command( + 'VO */Satellite/Chaser OrbitSystem Add System "Satellite/Target Body"' + ) + root.execute_command( + 'VO */Satellite/Chaser OrbitSystem Modify System "TargetBody" Show On' + ) + except: + print("") + finally: + print("################### Transfer orbit data ###################") + print("N phasing orbits = " + str(nOrbits)) + print("Period = " + str(tPhasingOrbit) + " sec") + print("SMA = " + str(transferSma) + " km") + print("Perigee radius = " + str(periRadius) + " km") + print("Apogee radius = " + str(apoRadius) + " km") + print("Delta V = " + str(deltaV) + " m/sec") + print("Transfer time = " + str(tPhasingOrbit * nOrbits / 86400) + " days") +# - From d83083fc395916a6b192f0c24b928bd163ac5c81 Mon Sep 17 00:00:00 2001 From: Jorge Martinez Date: Thu, 29 Jan 2026 13:27:23 +0100 Subject: [PATCH 02/18] fix: refactor --- examples/phasing-geo-orbits.py | 594 +++++++++++++++++---------------- 1 file changed, 309 insertions(+), 285 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 138eb287ff..b79480fd0c 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -1,336 +1,360 @@ # # GEO Orbit Phasing # -# Once we get the intended orbit, a change of the true anomaly is sometime needed to target a well-defined orbit slot. This could be the case of constellations that have more than one satellite in the same orbital plane and a certain relative geometry has to be respected or a geostationary satellite that is moved from its original location to another one for operational purposes. +# This tutorial provides a practical example on how to solve a phasing maneuver problem using Python. # -# In both cases we can solve for a reference (or target) satellite in Astrogator, whose position over time is taken as reference by a maneuvering (or chaser) satellite using a targeted MCS. +# ## What is orbit phasing? # -# ## STK Scenario +# Orbit phasing is an orbital maneuver used to adjust the position of a spacecraft within its orbit without changing the orbit's shape or orientation. This is particularly important for satellites in geostationary orbit (GEO) that need to move along the geostationary belt to a different orbital slot, or for constellation satellites that must maintain specific relative positions within the same orbital plane. # -# This lesson focuses on the design of phasing maneuver to be executed by a geostationary satellite who needs to change its position along the geostationary belt. +# The phasing maneuver involves temporarily changing the orbital period by adjusting the semi-major axis. By placing the satellite into a slightly different orbit (the phasing orbit), it will drift relative to its target position. After a specified number of orbits, the satellite returns to the original orbit at the desired location. # -#
Drawing
+# The key parameters for a phasing maneuver are: +# - **Phase angle**: The angular displacement between the chaser and target satellites +# - **Number of phasing orbits**: How many orbits are used to close the phase angle +# - **Phasing orbit**: A temporary orbit with a different period to achieve the desired drift +# +# ## Problem statement +# +# Two geostationary satellites occupy the same orbit but at different angular positions. The target satellite is at a true anomaly of 20 degrees, while the chaser satellite is at 70 degrees. Design a phasing maneuver for the chaser satellite to rendezvous with the target satellite using 5 phasing orbits over a 16-day period. +# +# Both satellites have the following initial orbital parameters: +# - Semi-major axis: 42164 km (GEO altitude) +# - Eccentricity: 0.0001 +# - Inclination: 0.0 degrees (equatorial) +# - RAAN: 0.0 degrees +# +# Compute the required $\Delta v$ for the phasing maneuvers and determine the phasing orbit parameters. -# ## Import Libraries +# ## Launch a new STK instance +# +# Start by launching a new STK instance. In this example, ``STKEngine`` is used with graphics (``no_graphics`` mode set to ``False``). This means that the graphic user interface (GUI) of the product is not launched but 2D and 3D visualization is still available through the STK Engine controls: # + -import math - -from ansys.stk.core.stkdesktop import STKDesktop -from ansys.stk.core.stkobjects import * -from ansys.stk.core.stkobjects.astrogator import * -from ansys.stk.core.utilities.colors import Colors +from ansys.stk.core.stkengine import STKEngine +stk = STKEngine.start_application(no_graphics=False) +print(f"Using {stk.version}") # - -# ### Create a New Scenario and Configure Satellites +# ## Create a new scenario # -# We now create a new scenario containing 2 satellites: the *Target* and the *Chaser*. Their initial states are not defined yet, but they will be on the geostationary orbit as specified in the next code cell. +# Start by creating a new scenario in STK: # + -stk = STKDesktop.start_application() -app = STKDesktop.start_application(visible=True) -root = app.root +from ansys.stk.core.stkobjects import PropagatorType, STKObjectType + + +root = stk.new_object_root() root.new_scenario("Phasing_orbit") -scen = Scenario(root.current_scenario) +scen = root.current_scenario -# configure the target satellite -targetSat = Satellite(scen.children.new(STKObjectType.SATELLITE, "Target")) +# - + +# ## Configure the target satellite +# +# Create the target satellite which will remain at a fixed position in geostationary orbit. This satellite serves as the reference point for the phasing maneuver: + +# + +from ansys.stk.core.stkobjects.astrogator import ElementSetType +from ansys.stk.core.utilities.colors import Colors + + +targetSat = scen.children.new(STKObjectType.SATELLITE, "Target") targetSat.set_propagator_type(PropagatorType.ASTROGATOR) -tDriver = MCSDriver(targetSat.propagator) +tDriver = targetSat.propagator tDriver.options.draw_trajectory_in_3d = False -tInitState = MCSInitialState(tDriver.main_sequence["Initial State"]) +tInitState = tDriver.main_sequence["Initial State"] tInitState.set_element_type(ElementSetType.KEPLERIAN) tInitState.initial_state.epoch = scen.start_time -tPropagator = MCSPropagate(tDriver.main_sequence["Propagate"]) +tPropagator = tDriver.main_sequence["Propagate"] tPropagator.propagator_name = "Earth point mass" tPropagator.properties.color = Colors.Red tKep = tInitState.element +# - + +# ## Configure the chaser satellite +# +# Create the chaser satellite which will perform the phasing maneuver to rendezvous with the target: -# configure the chaser satellite -chaserSat = Satellite(scen.children.new(STKObjectType.SATELLITE, "Chaser")) +# + +chaserSat = scen.children.new(STKObjectType.SATELLITE, "Chaser") chaserSat.set_propagator_type(PropagatorType.ASTROGATOR) -cDriver = MCSDriver(chaserSat.propagator) +cDriver = chaserSat.propagator cDriver.options.draw_trajectory_in_3d = False -cInitState = MCSInitialState(cDriver.main_sequence["Initial State"]) +cInitState = cDriver.main_sequence["Initial State"] cInitState.set_element_type(ElementSetType.KEPLERIAN) cInitState.initial_state.epoch = scen.start_time -cPropagator = MCSPropagate(cDriver.main_sequence["Propagate"]) +cPropagator = cDriver.main_sequence["Propagate"] cPropagator.propagator_name = "Earth point mass" cPropagator.stopping_conditions.add("Periapsis") cPropagator.stopping_conditions.remove("Duration") cKep = cInitState.element - -root.execute_command("VO */Satellite/Chaser Pass3D OrbitLead None") -root.execute_command("VO */Satellite/Chaser Pass3D OrbitTrail All") -root.execute_command("VO */Satellite/Target Pass3D OrbitLead None") -root.execute_command("VO */Satellite/Target Pass3D OrbitTrail All") # - -# ### Define the Input Data and Run the Analysis -# -# The cell below defines and runs the MCS for both satellites. -# -# The **phase angle** is the angular displacement between the two satellites at initial time. In this context, it is defined as always positive and between 0 and 360 degrees. To make the rendezvous between the satellites, it has to be reduced to 0 in a specified number of orbits (called **phasing orbits**), so each phase orbit will reduce the angular gap by just a fraction of the overall value. -# -# #### $$ \theta_{gap}=\frac{\theta}{n_{orbits}}$$ -# -# , where $\theta$ is the phase angle, and $\theta_{gap}$ is the angular dispacement to recover after each phasing orbit. -# -# In the code below the **Kepler's equation** is used to calculate the period and semimajor axis of the phasing orbit, given the number of revolution around the phasing orbit itself. As first, the **eccentric anomaly E** is calculated from $\theta_{gap}$ : -# -# #### $$\tan\frac{E}{2} = \sqrt{\frac{1 - e}{1 + e}}\tan\frac{\theta_{gap}}{2}$$ -# -# ...and then the period of the phasing orbit is derived: +# ## Define the input parameters # -# #### $$ T_{phasing}=\frac{T_{chaser}}{2\pi}\left ( E-e \sin E \right )$$ -# -# Having those data available, a differential corrector is configured to let the two satellites be in the same position at the end of the transfer +# The cell below defines the input parameters for the phasing maneuver analysis: # + +import math + + ########################################## INPUT DATA ########################################################## true_anom = 70 # deg - chaser true anomaly (between 0 and 360) target_ta = 20 # deg - target true anomaly (between 0 and 360) nOrbits = 5 # number of phasing orbits -prop_time = 16 # days +prop_time = 16 # days - total propagation time ################################################################################################################ -try: - root.execute_command( - 'VectorTool * CentralBody/Earth Create Angle Phase "Dihedral Angle"' - ) - root.execute_command( - 'VectorTool * CentralBody/Earth Modify Angle Phase "Dihedral Angle" "Satellite/Chaser Position" "Satellite/Target Position" "Satellite/Target Orbit_Normal" 0-360 Positive' - ) - root.execute_command( - 'VO * SetVectorGeometry Add "CentralBody/Earth Phase Angle" Color #fcba03 Show On' - ) - root.execute_command( - 'VO * SetVectorGeometry Add "Satellite/Target Position Vector"' - ) - root.execute_command( - 'VO * SetVectorGeometry Modify "Satellite/Target Position Vector" Color red Thickness 5 Show On' - ) - root.execute_command( - 'VO * SetVectorGeometry Add "Satellite/Chaser Position Vector"' - ) - root.execute_command( - 'VO * SetVectorGeometry Modify "Satellite/Chaser Position Vector" Color green Thickness 5 Show On' - ) -except: - print("") - - -# fixed keplerian parameters (geostationary orbit) -sma = 42164 +# Fixed Keplerian parameters (geostationary orbit) +sma = 42164 # km ecc = 0.0001 inc = 0.0 raan = 0.0 +# - -# change the scenario duration +# ## Configure initial orbital states +# +# Set the initial orbital elements for both satellites and run an initial propagation: +# + +# Change the scenario duration scen.stop_time = "+" + str(prop_time) + " day" -try: - cDriver.main_sequence.remove("Start Phasing") - cDriver.main_sequence.remove("Circularization") -except: - print("") -finally: - # set the chaser initial state - cKep.semimajor_axis = sma - cKep.eccentricity = ecc - cKep.inclination = inc - cKep.raan = raan - cKep.true_anomaly = true_anom - cKep.arg_of_periapsis = 0.0 - - # set the target initial state - tKep.semimajor_axis = sma - tKep.eccentricity = ecc - tKep.inclination = inc - tKep.raan = raan - tKep.true_anomaly = target_ta - tKep.arg_of_periapsis = 0.0 - - tPropagator.stopping_conditions["Duration"].properties.trip = prop_time * 86400 - - cDriver.run_mcs() - tDriver.run_mcs() - - # get the mean period and SMA of the initial orbit - meanKepDp = chaserSat.data_providers["Kozai-Izsak Mean"].group["ICRF"] - meanKep = meanKepDp.execute(scen.start_time, scen.start_time, 60) - meanKepDf = meanKep.data_sets.to_pandas_dataframe() - meanPeriod = float(meanKepDf.at[0, "mean nodal period"]) - meanSma = float(meanKepDf.at[0, "mean semi-major axis"]) - - # calculate the phase angle - phaseAngle = math.radians(target_ta - true_anom) - - if phaseAngle < 0: - phaseAngle = 2 * math.pi + phaseAngle - - print("#################### Initial geometry #####################") - print("Phase angle = " + str(math.degrees(phaseAngle)) + " deg") - print("") - - # recalculate the phase angle accordingly with the number of phasing orbits - phaseAngle = phaseAngle / nOrbits - - # calculate the eccentric anomaly - E = 2 * math.atan(math.sqrt((1 - ecc) / (1 + ecc)) * math.tan(phaseAngle / 2)) - - # calculate the time elapsed to cover phase angle in original orbit (Kepler's equation) - tPhaseAngle = (meanPeriod / (2 * math.pi)) * (E - ecc * math.sin(E)) - - # calculate the period of the phasing orbit - tPhasingOrbit = meanPeriod - tPhaseAngle - - # calculate the SMA of the phasing orbit - mu = 398600.44 - transferSma = math.pow((tPhasingOrbit * math.sqrt(mu)) / (2 * math.pi), 2 / 3) - - ## calculate eccentricity and radii - kepDp = chaserSat.data_providers["Astrogator Values"].group["Keplerian Elems"] - kep = kepDp.execute(scen.start_time, scen.start_time, 60) - kepDf = kep.data_sets.to_pandas_dataframe() - initRadius = float(kepDf.at[0, "radius_of_periapsis"]) - - if transferSma > sma: - periRadius = initRadius - apoRadius = 2 * transferSma - periRadius - else: - apoRadius = initRadius - periRadius = 2 * transferSma - apoRadius - - # estimate the needed Delta V for phasing - vPeriInitial = math.sqrt(mu * ((2 / initRadius) - (1 / meanSma))) - vPeriFinal = math.sqrt(mu * ((2 / initRadius) - (1 / transferSma))) - deltaV = vPeriFinal - vPeriInitial - - # add the first Target Sequence segment - ts1 = MCSTargetSequence( - cDriver.main_sequence.insert(SegmentType.TARGET_SEQUENCE, "Start Phasing", "-") - ) - ts1.action = TargetSequenceAction.RUN_ACTIVE_PROFILES - - ## add a Maneuver segment - dv1 = MCSManeuver(ts1.segments.insert(SegmentType.MANEUVER, "DV1", "-")) - dv1.set_maneuver_type(ManeuverType.IMPULSIVE) - dv1.enable_control_parameter( - ControlManeuver.IMPULSIVE_CARTESIAN_X - ) # control parameter - - #### maneuver attitude definition - dv1.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) - thrustVector = dv1.maneuver.attitude_control - thrustVector.thrust_axes_name = "Satellite/Chaser VNC(Earth)" - thrustVector.x = deltaV * 1000 - - ## add a Propagate segment - phasingOrbit = MCSPropagate( - ts1.segments.insert(SegmentType.PROPAGATE, "Phasing Orbit", "-") - ) - phasingOrbit.properties.color = Colors.Orange - phasingOrbit.propagator_name = "Earth point mass" - if transferSma > sma: - phasingOrbit.stopping_conditions.add("Periapsis") - else: - phasingOrbit.stopping_conditions.add("Apoapsis") - phasingOrbit.stopping_conditions.remove("Duration") - phasingOrbit.stopping_conditions.item(0).properties.repeat_count = nOrbits - phasingOrbit.results.add("Vector/Angle Between Vectors") - phasingOrbit.results[0].vector1_name = "Satellite/Chaser Position" - phasingOrbit.results[0].vector2_name = "Satellite/Target Position" - - # customize the Differential Corrector - dc1 = ProfileDifferentialCorrector(ts1.profiles["Differential Corrector"]) - dc1.mode = ProfileMode.ITERATE - dc1.max_iterations = 50 - - # set Control Parameters and Results - xControlParam1 = dc1.control_parameters.get_control_by_paths( - "DV1", "ImpulsiveMnvr.Cartesian.X" - ) - xControlParam1.enable = True - xControlParam1.max_step = 0.001 - - roaResult = dc1.results.get_result_by_paths( - "Phasing Orbit", "Angle_Between_Vectors" - ) - roaResult.enable = True - roaResult.desired_value = 0.0 - roaResult.tolerance = 0.1 - - # add the second Target Sequence segment - ts2 = MCSTargetSequence( - cDriver.main_sequence.insert( - SegmentType.TARGET_SEQUENCE, "Circularization", "-" - ) - ) - ts2.action = TargetSequenceAction.RUN_ACTIVE_PROFILES - - ## add a Maneuver segment - dv2 = MCSManeuver(ts2.segments.insert(SegmentType.MANEUVER, "DV2", "-")) - dv2.set_maneuver_type(ManeuverType.IMPULSIVE) - dv2.enable_control_parameter( - ControlManeuver.IMPULSIVE_CARTESIAN_X - ) # control parameter - - ## add a Propagate segment - finalOrbit = MCSPropagate( - ts2.segments.insert(SegmentType.PROPAGATE, "Final Orbit", "-") - ) - finalOrbit.properties.color = Colors.Yellow - finalOrbit.propagator_name = "Earth point mass" - finalOrbit.stopping_conditions["Duration"].properties.trip = 86400 - finalOrbit.results.add("Keplerian Elems/Eccentricity") - - # customize the Differential Corrector - dc2 = ProfileDifferentialCorrector(ts2.profiles["Differential Corrector"]) - dc2.mode = ProfileMode.ITERATE - dc2.max_iterations = 50 - - # set Control Parameters and Results - xControlParam2 = dc2.control_parameters.get_control_by_paths( - "DV2", "ImpulsiveMnvr.Cartesian.X" - ) - xControlParam2.enable = True - xControlParam2.max_step = 0.01 - - roaResult = dc2.results.get_result_by_paths("Final Orbit", "Eccentricity") - roaResult.enable = True - roaResult.desired_value = 0.0001 - roaResult.tolerance = 0.00001 - cDriver.run_mcs() - tDriver.run_mcs() - root.rewind() - - # get the maneuver data providers - manDp = chaserSat.data_providers["Maneuver Summary"] - man = manDp.execute(scen.start_time, scen.stop_time) - manDf = man.data_sets.to_pandas_dataframe() - deltaV = manDf.at[0, "delta v"] - - # change orbit system - try: - root.execute_command( - 'VO */Satellite/Chaser OrbitSystem Modify System "InertialByWindow" Show Off' - ) - root.execute_command( - 'VO */Satellite/Chaser OrbitSystem Add System "Satellite/Target Body"' - ) - root.execute_command( - 'VO */Satellite/Chaser OrbitSystem Modify System "TargetBody" Show On' - ) - except: - print("") - finally: - print("################### Transfer orbit data ###################") - print("N phasing orbits = " + str(nOrbits)) - print("Period = " + str(tPhasingOrbit) + " sec") - print("SMA = " + str(transferSma) + " km") - print("Perigee radius = " + str(periRadius) + " km") - print("Apogee radius = " + str(apoRadius) + " km") - print("Delta V = " + str(deltaV) + " m/sec") - print("Transfer time = " + str(tPhasingOrbit * nOrbits / 86400) + " days") + +# Set the chaser initial state +cKep.semimajor_axis = sma +cKep.eccentricity = ecc +cKep.inclination = inc +cKep.raan = raan +cKep.true_anomaly = true_anom +cKep.arg_of_periapsis = 0.0 + +# Set the target initial state +tKep.semimajor_axis = sma +tKep.eccentricity = ecc +tKep.inclination = inc +tKep.raan = raan +tKep.true_anomaly = target_ta +tKep.arg_of_periapsis = 0.0 + +# Configure target propagation duration +tPropagator.stopping_conditions["Duration"].properties.trip = prop_time * 86400 + +# Run initial propagation for both satellites +cDriver.run_mcs() +tDriver.run_mcs() +# - + +# ## Calculate phasing orbit parameters +# +# The **phase angle** is the angular displacement between the two satellites at initial time. To achieve rendezvous, this angle must be reduced to 0 over a specified number of phasing orbits. +# +# The angular gap per orbit is: +# $$ \theta_{gap}=\frac{\theta}{n_{orbits}}$$ +# +# where $\theta$ is the phase angle, and $\theta_{gap}$ is the angular displacement to recover after each phasing orbit. +# +# Using **Kepler's equation**, we calculate the period and semi-major axis of the phasing orbit. First, the **eccentric anomaly E** is calculated from $\theta_{gap}$: +# +# $$\tan\frac{E}{2} = \sqrt{\frac{1 - e}{1 + e}}\tan\frac{\theta_{gap}}{2}$$ +# +# Then the period of the phasing orbit is derived: +# +# $$ T_{phasing}=\frac{T_{chaser}}{2\pi}\left ( E-e \sin E \right )$$ + +# + +# Get the mean period and SMA of the initial orbit +meanKepDp = chaserSat.data_providers["Kozai-Izsak Mean"].group["ICRF"] +meanKep = meanKepDp.execute(scen.start_time, scen.start_time, 60) +meanKepDf = meanKep.data_sets.to_pandas_dataframe() +meanPeriod = float(meanKepDf.at[0, "mean nodal period"]) +meanSma = float(meanKepDf.at[0, "mean semi-major axis"]) + +# Calculate the phase angle +phaseAngle = math.radians(target_ta - true_anom) + +if phaseAngle < 0: + phaseAngle = 2 * math.pi + phaseAngle + +print("#################### Initial geometry #####################") +print("Phase angle = " + str(math.degrees(phaseAngle)) + " deg") +print("") + +# Recalculate the phase angle accordingly with the number of phasing orbits +phaseAngle = phaseAngle / nOrbits + +# Calculate the eccentric anomaly +E = 2 * math.atan(math.sqrt((1 - ecc) / (1 + ecc)) * math.tan(phaseAngle / 2)) + +# Calculate the time elapsed to cover phase angle in original orbit (Kepler's equation) +tPhaseAngle = (meanPeriod / (2 * math.pi)) * (E - ecc * math.sin(E)) + +# Calculate the period of the phasing orbit +tPhasingOrbit = meanPeriod - tPhaseAngle + +# Calculate the SMA of the phasing orbit +mu = 398600.44 +transferSma = math.pow((tPhasingOrbit * math.sqrt(mu)) / (2 * math.pi), 2 / 3) + +# Calculate eccentricity and radii +kepDp = chaserSat.data_providers["Astrogator Values"].group["Keplerian Elems"] +kep = kepDp.execute(scen.start_time, scen.start_time, 60) +kepDf = kep.data_sets.to_pandas_dataframe() +initRadius = float(kepDf.at[0, "radius_of_periapsis"]) + +if transferSma > sma: + periRadius = initRadius + apoRadius = 2 * transferSma - periRadius +else: + apoRadius = initRadius + periRadius = 2 * transferSma - apoRadius + +# Estimate the needed Delta V for phasing +vPeriInitial = math.sqrt(mu * ((2 / initRadius) - (1 / meanSma))) +vPeriFinal = math.sqrt(mu * ((2 / initRadius) - (1 / transferSma))) +deltaV = vPeriFinal - vPeriInitial +# - + +# ## Set up the phasing maneuver sequence +# +# Add a target sequence to perform the first delta-V maneuver and propagate through the phasing orbits: + +# + +from ansys.stk.core.stkobjects.astrogator import ( + AttitudeControl, + ControlManeuver, + ManeuverType, + ProfileMode, + SegmentType, + TargetSequenceAction, +) + + +# Add the first Target Sequence segment +ts1 = cDriver.main_sequence.insert(SegmentType.TARGET_SEQUENCE, "Start Phasing", "-") +ts1.action = TargetSequenceAction.RUN_ACTIVE_PROFILES + +# Add a Maneuver segment +dv1 = ts1.segments.insert(SegmentType.MANEUVER, "DV1", "-") +dv1.set_maneuver_type(ManeuverType.IMPULSIVE) +dv1.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) + +# Maneuver attitude definition +dv1.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) +thrustVector = dv1.maneuver.attitude_control +thrustVector.thrust_axes_name = "Satellite/Chaser VNC(Earth)" +thrustVector.x = deltaV * 1000 + +# Add a Propagate segment +phasingOrbit = ts1.segments.insert(SegmentType.PROPAGATE, "Phasing Orbit", "-") +phasingOrbit.properties.color = Colors.Orange +phasingOrbit.propagator_name = "Earth point mass" +if transferSma > sma: + phasingOrbit.stopping_conditions.add("Periapsis") +else: + phasingOrbit.stopping_conditions.add("Apoapsis") +phasingOrbit.stopping_conditions.remove("Duration") +phasingOrbit.stopping_conditions.item(0).properties.repeat_count = nOrbits +phasingOrbit.results.add("Vector/Angle Between Vectors") +phasingOrbit.results[0].vector1_name = "Satellite/Chaser Position" +phasingOrbit.results[0].vector2_name = "Satellite/Target Position" +# - + +# ## Configure the differential corrector for phasing +# +# Set up the differential corrector to adjust the first delta-V to achieve zero phase angle at the end of the phasing orbits: + +# + +# Customize the Differential Corrector +dc1 = ts1.profiles["Differential Corrector"] +dc1.mode = ProfileMode.ITERATE +dc1.max_iterations = 50 + +# Set Control Parameters and Results +xControlParam1 = dc1.control_parameters.get_control_by_paths( + "DV1", "ImpulsiveMnvr.Cartesian.X" +) +xControlParam1.enable = True +xControlParam1.max_step = 0.001 + +roaResult = dc1.results.get_result_by_paths("Phasing Orbit", "Angle_Between_Vectors") +roaResult.enable = True +roaResult.desired_value = 0.0 +roaResult.tolerance = 0.1 +# - + +# ## Set up the circularization maneuver +# +# Add a second target sequence to perform the circularization delta-V to return to the original circular orbit: + +# + +# Add the second Target Sequence segment +ts2 = cDriver.main_sequence.insert(SegmentType.TARGET_SEQUENCE, "Circularization", "-") +ts2.action = TargetSequenceAction.RUN_ACTIVE_PROFILES + +# Add a Maneuver segment +dv2 = ts2.segments.insert(SegmentType.MANEUVER, "DV2", "-") +dv2.set_maneuver_type(ManeuverType.IMPULSIVE) +dv2.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) + +# Add a Propagate segment +finalOrbit = ts2.segments.insert(SegmentType.PROPAGATE, "Final Orbit", "-") +finalOrbit.properties.color = Colors.Yellow +finalOrbit.propagator_name = "Earth point mass" +finalOrbit.stopping_conditions["Duration"].properties.trip = 86400 +finalOrbit.results.add("Keplerian Elems/Eccentricity") +# - + +# ## Configure the differential corrector for circularization +# +# Set up the differential corrector to achieve a circular orbit with the desired eccentricity: + +# + +# Customize the Differential Corrector +dc2 = ts2.profiles["Differential Corrector"] +dc2.mode = ProfileMode.ITERATE +dc2.max_iterations = 50 + +# Set Control Parameters and Results +xControlParam2 = dc2.control_parameters.get_control_by_paths( + "DV2", "ImpulsiveMnvr.Cartesian.X" +) +xControlParam2.enable = True +xControlParam2.max_step = 0.01 + +roaResult = dc2.results.get_result_by_paths("Final Orbit", "Eccentricity") +roaResult.enable = True +roaResult.desired_value = 0.0001 +roaResult.tolerance = 0.00001 +# - + +# ## Run the main control sequence +# +# Execute the mission control sequence to solve for the phasing maneuver: + +# + +cDriver.run_mcs() +tDriver.run_mcs() +root.rewind() +# - + +# ## Retrieve the results +# +# Once the analysis has been performed, retrieve the delta-V values and phasing orbit parameters: + +# + +# Get the maneuver data providers +manDp = chaserSat.data_providers["Maneuver Summary"] +man = manDp.execute(scen.start_time, scen.stop_time) +manDf = man.data_sets.to_pandas_dataframe() +deltaV = manDf.at[0, "delta v"] + +print("################### Transfer orbit data ###################") +print("N phasing orbits = " + str(nOrbits)) +print("Period = " + str(tPhasingOrbit) + " sec") +print("SMA = " + str(transferSma) + " km") +print("Perigee radius = " + str(periRadius) + " km") +print("Apogee radius = " + str(apoRadius) + " km") +print("Delta V = " + str(deltaV) + " m/sec") +print("Transfer time = " + str(tPhasingOrbit * nOrbits / 86400) + " days") # - From 6ca6c3e726654373cbb7071d4af25507d7936f47 Mon Sep 17 00:00:00 2001 From: Jorge Martinez Date: Thu, 29 Jan 2026 14:02:38 +0100 Subject: [PATCH 03/18] fix: refactor --- examples/phasing-geo-orbits.py | 360 ++++++++++++++++++--------------- 1 file changed, 196 insertions(+), 164 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index b79480fd0c..cae828b1a9 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -55,95 +55,89 @@ # # Create the target satellite which will remain at a fixed position in geostationary orbit. This satellite serves as the reference point for the phasing maneuver: -# + from ansys.stk.core.stkobjects.astrogator import ElementSetType from ansys.stk.core.utilities.colors import Colors -targetSat = scen.children.new(STKObjectType.SATELLITE, "Target") -targetSat.set_propagator_type(PropagatorType.ASTROGATOR) -tDriver = targetSat.propagator -tDriver.options.draw_trajectory_in_3d = False -tInitState = tDriver.main_sequence["Initial State"] -tInitState.set_element_type(ElementSetType.KEPLERIAN) -tInitState.initial_state.epoch = scen.start_time -tPropagator = tDriver.main_sequence["Propagate"] -tPropagator.propagator_name = "Earth point mass" -tPropagator.properties.color = Colors.Red -tKep = tInitState.element -# - +target_satellite = scen.children.new(STKObjectType.SATELLITE, "Target") +target_satellite.set_propagator_type(PropagatorType.ASTROGATOR) +target_propagator = target_satellite.propagator +target_propagator.options.draw_trajectory_in_3d = False +target_initial_state = target_propagator.main_sequence["Initial State"] +target_initial_state.set_element_type(ElementSetType.KEPLERIAN) +target_initial_state.initial_state.epoch = scen.start_time +target_propagate_segment = target_propagator.main_sequence["Propagate"] +target_propagate_segment.propagator_name = "Earth point mass" +target_propagate_segment.properties.color = Colors.Red +target_keplerian_elements = target_initial_state.element # ## Configure the chaser satellite # # Create the chaser satellite which will perform the phasing maneuver to rendezvous with the target: -# + -chaserSat = scen.children.new(STKObjectType.SATELLITE, "Chaser") -chaserSat.set_propagator_type(PropagatorType.ASTROGATOR) -cDriver = chaserSat.propagator -cDriver.options.draw_trajectory_in_3d = False -cInitState = cDriver.main_sequence["Initial State"] -cInitState.set_element_type(ElementSetType.KEPLERIAN) -cInitState.initial_state.epoch = scen.start_time -cPropagator = cDriver.main_sequence["Propagate"] -cPropagator.propagator_name = "Earth point mass" -cPropagator.stopping_conditions.add("Periapsis") -cPropagator.stopping_conditions.remove("Duration") -cKep = cInitState.element -# - +chaser_satellite = scen.children.new(STKObjectType.SATELLITE, "Chaser") +chaser_satellite.set_propagator_type(PropagatorType.ASTROGATOR) +chaser_propagator = chaser_satellite.propagator +chaser_propagator.options.draw_trajectory_in_3d = False +chaser_initial_state = chaser_propagator.main_sequence["Initial State"] +chaser_initial_state.set_element_type(ElementSetType.KEPLERIAN) +chaser_initial_state.initial_state.epoch = scen.start_time +chaser_propagate_segment = chaser_propagator.main_sequence["Propagate"] +chaser_propagate_segment.propagator_name = "Earth point mass" +chaser_propagate_segment.stopping_conditions.add("Periapsis") +chaser_propagate_segment.stopping_conditions.remove("Duration") +chaser_keplerian_elements = chaser_initial_state.element # ## Define the input parameters # # The cell below defines the input parameters for the phasing maneuver analysis: -# + import math -########################################## INPUT DATA ########################################################## -true_anom = 70 # deg - chaser true anomaly (between 0 and 360) -target_ta = 20 # deg - target true anomaly (between 0 and 360) -nOrbits = 5 # number of phasing orbits -prop_time = 16 # days - total propagation time -################################################################################################################ +# Input parameters for phasing maneuver +chaser_true_anomaly = 70 # deg - chaser true anomaly (between 0 and 360) +target_true_anomaly = 20 # deg - target true anomaly (between 0 and 360) +number_of_phasing_orbits = 5 # number of phasing orbits +propagation_time = 16 # days - total propagation time # Fixed Keplerian parameters (geostationary orbit) -sma = 42164 # km -ecc = 0.0001 -inc = 0.0 -raan = 0.0 -# - +semi_major_axis = 42164 # km +eccentricity = 0.0001 +inclination = 0.0 +right_ascension_ascending_node = 0.0 # ## Configure initial orbital states # # Set the initial orbital elements for both satellites and run an initial propagation: -# + + # Change the scenario duration -scen.stop_time = "+" + str(prop_time) + " day" +scen.stop_time = "+" + str(propagation_time) + " day" # Set the chaser initial state -cKep.semimajor_axis = sma -cKep.eccentricity = ecc -cKep.inclination = inc -cKep.raan = raan -cKep.true_anomaly = true_anom -cKep.arg_of_periapsis = 0.0 +chaser_keplerian_elements.semimajor_axis = semi_major_axis +chaser_keplerian_elements.eccentricity = eccentricity +chaser_keplerian_elements.inclination = inclination +chaser_keplerian_elements.raan = right_ascension_ascending_node +chaser_keplerian_elements.true_anomaly = chaser_true_anomaly +chaser_keplerian_elements.arg_of_periapsis = 0.0 # Set the target initial state -tKep.semimajor_axis = sma -tKep.eccentricity = ecc -tKep.inclination = inc -tKep.raan = raan -tKep.true_anomaly = target_ta -tKep.arg_of_periapsis = 0.0 +target_keplerian_elements.semimajor_axis = semi_major_axis +target_keplerian_elements.eccentricity = eccentricity +target_keplerian_elements.inclination = inclination +target_keplerian_elements.raan = right_ascension_ascending_node +target_keplerian_elements.true_anomaly = target_true_anomaly +target_keplerian_elements.arg_of_periapsis = 0.0 # Configure target propagation duration -tPropagator.stopping_conditions["Duration"].properties.trip = prop_time * 86400 +target_propagate_segment.stopping_conditions["Duration"].properties.trip = ( + propagation_time * 86400 +) # Run initial propagation for both satellites -cDriver.run_mcs() -tDriver.run_mcs() -# - +chaser_propagator.run_mcs() +target_propagator.run_mcs() # ## Calculate phasing orbit parameters # @@ -162,64 +156,81 @@ # # $$ T_{phasing}=\frac{T_{chaser}}{2\pi}\left ( E-e \sin E \right )$$ -# + # Get the mean period and SMA of the initial orbit -meanKepDp = chaserSat.data_providers["Kozai-Izsak Mean"].group["ICRF"] -meanKep = meanKepDp.execute(scen.start_time, scen.start_time, 60) -meanKepDf = meanKep.data_sets.to_pandas_dataframe() -meanPeriod = float(meanKepDf.at[0, "mean nodal period"]) -meanSma = float(meanKepDf.at[0, "mean semi-major axis"]) +mean_keplerian_data_provider = chaser_satellite.data_providers[ + "Kozai-Izsak Mean" +].group["ICRF"] +mean_keplerian_result = mean_keplerian_data_provider.execute( + scen.start_time, scen.start_time, 60 +) +mean_keplerian_dataframe = mean_keplerian_result.data_sets.to_pandas_dataframe() +mean_orbital_period = float(mean_keplerian_dataframe.at[0, "mean nodal period"]) +mean_semi_major_axis = float(mean_keplerian_dataframe.at[0, "mean semi-major axis"]) # Calculate the phase angle -phaseAngle = math.radians(target_ta - true_anom) +phase_angle = math.radians(target_true_anomaly - chaser_true_anomaly) -if phaseAngle < 0: - phaseAngle = 2 * math.pi + phaseAngle +if phase_angle < 0: + phase_angle = 2 * math.pi + phase_angle print("#################### Initial geometry #####################") -print("Phase angle = " + str(math.degrees(phaseAngle)) + " deg") +print("Phase angle = " + str(math.degrees(phase_angle)) + " deg") print("") # Recalculate the phase angle accordingly with the number of phasing orbits -phaseAngle = phaseAngle / nOrbits +phase_angle = phase_angle / number_of_phasing_orbits # Calculate the eccentric anomaly -E = 2 * math.atan(math.sqrt((1 - ecc) / (1 + ecc)) * math.tan(phaseAngle / 2)) +eccentric_anomaly = 2 * math.atan( + math.sqrt((1 - eccentricity) / (1 + eccentricity)) * math.tan(phase_angle / 2) +) # Calculate the time elapsed to cover phase angle in original orbit (Kepler's equation) -tPhaseAngle = (meanPeriod / (2 * math.pi)) * (E - ecc * math.sin(E)) +time_to_cover_phase_angle = (mean_orbital_period / (2 * math.pi)) * ( + eccentric_anomaly - eccentricity * math.sin(eccentric_anomaly) +) # Calculate the period of the phasing orbit -tPhasingOrbit = meanPeriod - tPhaseAngle +phasing_orbit_period = mean_orbital_period - time_to_cover_phase_angle # Calculate the SMA of the phasing orbit -mu = 398600.44 -transferSma = math.pow((tPhasingOrbit * math.sqrt(mu)) / (2 * math.pi), 2 / 3) +earth_gravitational_parameter = 398600.44 +transfer_semi_major_axis = math.pow( + (phasing_orbit_period * math.sqrt(earth_gravitational_parameter)) / (2 * math.pi), + 2 / 3, +) # Calculate eccentricity and radii -kepDp = chaserSat.data_providers["Astrogator Values"].group["Keplerian Elems"] -kep = kepDp.execute(scen.start_time, scen.start_time, 60) -kepDf = kep.data_sets.to_pandas_dataframe() -initRadius = float(kepDf.at[0, "radius_of_periapsis"]) - -if transferSma > sma: - periRadius = initRadius - apoRadius = 2 * transferSma - periRadius +keplerian_elements_data_provider = chaser_satellite.data_providers[ + "Astrogator Values" +].group["Keplerian Elems"] +keplerian_result = keplerian_elements_data_provider.execute( + scen.start_time, scen.start_time, 60 +) +keplerian_dataframe = keplerian_result.data_sets.to_pandas_dataframe() +initial_radius = float(keplerian_dataframe.at[0, "radius_of_periapsis"]) + +if transfer_semi_major_axis > semi_major_axis: + periapsis_radius = initial_radius + apoapsis_radius = 2 * transfer_semi_major_axis - periapsis_radius else: - apoRadius = initRadius - periRadius = 2 * transferSma - apoRadius + apoapsis_radius = initial_radius + periapsis_radius = 2 * transfer_semi_major_axis - apoapsis_radius # Estimate the needed Delta V for phasing -vPeriInitial = math.sqrt(mu * ((2 / initRadius) - (1 / meanSma))) -vPeriFinal = math.sqrt(mu * ((2 / initRadius) - (1 / transferSma))) -deltaV = vPeriFinal - vPeriInitial -# - +velocity_periapsis_initial = math.sqrt( + earth_gravitational_parameter * ((2 / initial_radius) - (1 / mean_semi_major_axis)) +) +velocity_periapsis_final = math.sqrt( + earth_gravitational_parameter + * ((2 / initial_radius) - (1 / transfer_semi_major_axis)) +) +delta_v_estimate = velocity_periapsis_final - velocity_periapsis_initial # ## Set up the phasing maneuver sequence # # Add a target sequence to perform the first delta-V maneuver and propagate through the phasing orbits: -# + from ansys.stk.core.stkobjects.astrogator import ( AttitudeControl, ControlManeuver, @@ -231,130 +242,151 @@ # Add the first Target Sequence segment -ts1 = cDriver.main_sequence.insert(SegmentType.TARGET_SEQUENCE, "Start Phasing", "-") -ts1.action = TargetSequenceAction.RUN_ACTIVE_PROFILES +phasing_start_sequence = chaser_propagator.main_sequence.insert( + SegmentType.TARGET_SEQUENCE, "Start Phasing", "-" +) +phasing_start_sequence.action = TargetSequenceAction.RUN_ACTIVE_PROFILES # Add a Maneuver segment -dv1 = ts1.segments.insert(SegmentType.MANEUVER, "DV1", "-") -dv1.set_maneuver_type(ManeuverType.IMPULSIVE) -dv1.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) +first_delta_v_maneuver = phasing_start_sequence.segments.insert( + SegmentType.MANEUVER, "DV1", "-" +) +first_delta_v_maneuver.set_maneuver_type(ManeuverType.IMPULSIVE) +first_delta_v_maneuver.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) # Maneuver attitude definition -dv1.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) -thrustVector = dv1.maneuver.attitude_control -thrustVector.thrust_axes_name = "Satellite/Chaser VNC(Earth)" -thrustVector.x = deltaV * 1000 +first_delta_v_maneuver.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) +thrust_vector = first_delta_v_maneuver.maneuver.attitude_control +thrust_vector.thrust_axes_name = "Satellite/Chaser VNC(Earth)" +thrust_vector.x = delta_v_estimate * 1000 # Add a Propagate segment -phasingOrbit = ts1.segments.insert(SegmentType.PROPAGATE, "Phasing Orbit", "-") -phasingOrbit.properties.color = Colors.Orange -phasingOrbit.propagator_name = "Earth point mass" -if transferSma > sma: - phasingOrbit.stopping_conditions.add("Periapsis") +phasing_orbit_propagate = phasing_start_sequence.segments.insert( + SegmentType.PROPAGATE, "Phasing Orbit", "-" +) +phasing_orbit_propagate.properties.color = Colors.Orange +phasing_orbit_propagate.propagator_name = "Earth point mass" +if transfer_semi_major_axis > semi_major_axis: + phasing_orbit_propagate.stopping_conditions.add("Periapsis") else: - phasingOrbit.stopping_conditions.add("Apoapsis") -phasingOrbit.stopping_conditions.remove("Duration") -phasingOrbit.stopping_conditions.item(0).properties.repeat_count = nOrbits -phasingOrbit.results.add("Vector/Angle Between Vectors") -phasingOrbit.results[0].vector1_name = "Satellite/Chaser Position" -phasingOrbit.results[0].vector2_name = "Satellite/Target Position" -# - + phasing_orbit_propagate.stopping_conditions.add("Apoapsis") +phasing_orbit_propagate.stopping_conditions.remove("Duration") +phasing_orbit_propagate.stopping_conditions.item( + 0 +).properties.repeat_count = number_of_phasing_orbits +phasing_orbit_propagate.results.add("Vector/Angle Between Vectors") +phasing_orbit_propagate.results[0].vector1_name = "Satellite/Chaser Position" +phasing_orbit_propagate.results[0].vector2_name = "Satellite/Target Position" # ## Configure the differential corrector for phasing # # Set up the differential corrector to adjust the first delta-V to achieve zero phase angle at the end of the phasing orbits: -# + # Customize the Differential Corrector -dc1 = ts1.profiles["Differential Corrector"] -dc1.mode = ProfileMode.ITERATE -dc1.max_iterations = 50 +phasing_differential_corrector = phasing_start_sequence.profiles[ + "Differential Corrector" +] +phasing_differential_corrector.mode = ProfileMode.ITERATE +phasing_differential_corrector.max_iterations = 50 # Set Control Parameters and Results -xControlParam1 = dc1.control_parameters.get_control_by_paths( - "DV1", "ImpulsiveMnvr.Cartesian.X" +phasing_x_control_parameter = ( + phasing_differential_corrector.control_parameters.get_control_by_paths( + "DV1", "ImpulsiveMnvr.Cartesian.X" + ) ) -xControlParam1.enable = True -xControlParam1.max_step = 0.001 +phasing_x_control_parameter.enable = True +phasing_x_control_parameter.max_step = 0.001 -roaResult = dc1.results.get_result_by_paths("Phasing Orbit", "Angle_Between_Vectors") -roaResult.enable = True -roaResult.desired_value = 0.0 -roaResult.tolerance = 0.1 -# - +phasing_angle_result = phasing_differential_corrector.results.get_result_by_paths( + "Phasing Orbit", "Angle_Between_Vectors" +) +phasing_angle_result.enable = True +phasing_angle_result.desired_value = 0.0 +phasing_angle_result.tolerance = 0.1 # ## Set up the circularization maneuver # # Add a second target sequence to perform the circularization delta-V to return to the original circular orbit: -# + # Add the second Target Sequence segment -ts2 = cDriver.main_sequence.insert(SegmentType.TARGET_SEQUENCE, "Circularization", "-") -ts2.action = TargetSequenceAction.RUN_ACTIVE_PROFILES +circularization_sequence = chaser_propagator.main_sequence.insert( + SegmentType.TARGET_SEQUENCE, "Circularization", "-" +) +circularization_sequence.action = TargetSequenceAction.RUN_ACTIVE_PROFILES # Add a Maneuver segment -dv2 = ts2.segments.insert(SegmentType.MANEUVER, "DV2", "-") -dv2.set_maneuver_type(ManeuverType.IMPULSIVE) -dv2.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) +second_delta_v_maneuver = circularization_sequence.segments.insert( + SegmentType.MANEUVER, "DV2", "-" +) +second_delta_v_maneuver.set_maneuver_type(ManeuverType.IMPULSIVE) +second_delta_v_maneuver.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) # Add a Propagate segment -finalOrbit = ts2.segments.insert(SegmentType.PROPAGATE, "Final Orbit", "-") -finalOrbit.properties.color = Colors.Yellow -finalOrbit.propagator_name = "Earth point mass" -finalOrbit.stopping_conditions["Duration"].properties.trip = 86400 -finalOrbit.results.add("Keplerian Elems/Eccentricity") -# - +final_orbit_propagate = circularization_sequence.segments.insert( + SegmentType.PROPAGATE, "Final Orbit", "-" +) +final_orbit_propagate.properties.color = Colors.Yellow +final_orbit_propagate.propagator_name = "Earth point mass" +final_orbit_propagate.stopping_conditions["Duration"].properties.trip = 86400 +final_orbit_propagate.results.add("Keplerian Elems/Eccentricity") # ## Configure the differential corrector for circularization # # Set up the differential corrector to achieve a circular orbit with the desired eccentricity: -# + # Customize the Differential Corrector -dc2 = ts2.profiles["Differential Corrector"] -dc2.mode = ProfileMode.ITERATE -dc2.max_iterations = 50 +circularization_differential_corrector = circularization_sequence.profiles[ + "Differential Corrector" +] +circularization_differential_corrector.mode = ProfileMode.ITERATE +circularization_differential_corrector.max_iterations = 50 # Set Control Parameters and Results -xControlParam2 = dc2.control_parameters.get_control_by_paths( - "DV2", "ImpulsiveMnvr.Cartesian.X" +circularization_x_control_parameter = ( + circularization_differential_corrector.control_parameters.get_control_by_paths( + "DV2", "ImpulsiveMnvr.Cartesian.X" + ) ) -xControlParam2.enable = True -xControlParam2.max_step = 0.01 +circularization_x_control_parameter.enable = True +circularization_x_control_parameter.max_step = 0.01 -roaResult = dc2.results.get_result_by_paths("Final Orbit", "Eccentricity") -roaResult.enable = True -roaResult.desired_value = 0.0001 -roaResult.tolerance = 0.00001 -# - +circularization_eccentricity_result = ( + circularization_differential_corrector.results.get_result_by_paths( + "Final Orbit", "Eccentricity" + ) +) +circularization_eccentricity_result.enable = True +circularization_eccentricity_result.desired_value = 0.0001 +circularization_eccentricity_result.tolerance = 0.00001 # ## Run the main control sequence # # Execute the mission control sequence to solve for the phasing maneuver: -# + -cDriver.run_mcs() -tDriver.run_mcs() +chaser_propagator.run_mcs() +target_propagator.run_mcs() root.rewind() -# - # ## Retrieve the results # # Once the analysis has been performed, retrieve the delta-V values and phasing orbit parameters: -# + # Get the maneuver data providers -manDp = chaserSat.data_providers["Maneuver Summary"] -man = manDp.execute(scen.start_time, scen.stop_time) -manDf = man.data_sets.to_pandas_dataframe() -deltaV = manDf.at[0, "delta v"] +maneuver_data_provider = chaser_satellite.data_providers["Maneuver Summary"] +maneuver_result = maneuver_data_provider.execute(scen.start_time, scen.stop_time) +maneuver_dataframe = maneuver_result.data_sets.to_pandas_dataframe() +delta_v_actual = maneuver_dataframe.at[0, "delta v"] print("################### Transfer orbit data ###################") -print("N phasing orbits = " + str(nOrbits)) -print("Period = " + str(tPhasingOrbit) + " sec") -print("SMA = " + str(transferSma) + " km") -print("Perigee radius = " + str(periRadius) + " km") -print("Apogee radius = " + str(apoRadius) + " km") -print("Delta V = " + str(deltaV) + " m/sec") -print("Transfer time = " + str(tPhasingOrbit * nOrbits / 86400) + " days") -# - +print("N phasing orbits = " + str(number_of_phasing_orbits)) +print("Period = " + str(phasing_orbit_period) + " sec") +print("SMA = " + str(transfer_semi_major_axis) + " km") +print("Perigee radius = " + str(periapsis_radius) + " km") +print("Apogee radius = " + str(apoapsis_radius) + " km") +print("Delta V = " + str(delta_v_actual) + " m/sec") +print( + "Transfer time = " + + str(phasing_orbit_period * number_of_phasing_orbits / 86400) + + " days" +) From fcfa44b9fd4cb65991987cd5e12ea98afbc751e5 Mon Sep 17 00:00:00 2001 From: Jorge Martinez Date: Thu, 29 Jan 2026 14:10:23 +0100 Subject: [PATCH 04/18] fix: refactor --- examples/phasing-geo-orbits.py | 43 ++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 10 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index cae828b1a9..1787b6f84c 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -55,21 +55,44 @@ # # Create the target satellite which will remain at a fixed position in geostationary orbit. This satellite serves as the reference point for the phasing maneuver: -from ansys.stk.core.stkobjects.astrogator import ElementSetType -from ansys.stk.core.utilities.colors import Colors +# + +target_satellite = scen.children.new(STKObjectType.SATELLITE, "Target") +# - +# Then, declare the type of orbit propagator used for the satellite: -target_satellite = scen.children.new(STKObjectType.SATELLITE, "Target") +# + target_satellite.set_propagator_type(PropagatorType.ASTROGATOR) -target_propagator = target_satellite.propagator -target_propagator.options.draw_trajectory_in_3d = False -target_initial_state = target_propagator.main_sequence["Initial State"] +# - + +# Initialize the propagator and configure additional settings: + +target_satellite.propagator.options.draw_trajectory_in_3d = False + +# ## Set up the initial state of the target satellite +# +# Access the existing initial state segment in the main sequence and configure the element type: + +# + +from ansys.stk.core.stkobjects.astrogator import ElementSetType + + +target_initial_state = target_satellite.propagator.main_sequence["Initial State"] target_initial_state.set_element_type(ElementSetType.KEPLERIAN) target_initial_state.initial_state.epoch = scen.start_time -target_propagate_segment = target_propagator.main_sequence["Propagate"] +target_keplerian_elements = target_initial_state.element +# - + +# Configure the propagation segment: + +# + +from ansys.stk.core.utilities.colors import Colors + + +target_propagate_segment = target_satellite.propagator.main_sequence["Propagate"] target_propagate_segment.propagator_name = "Earth point mass" target_propagate_segment.properties.color = Colors.Red -target_keplerian_elements = target_initial_state.element +# - # ## Configure the chaser satellite # @@ -137,7 +160,7 @@ # Run initial propagation for both satellites chaser_propagator.run_mcs() -target_propagator.run_mcs() +target_satellite.propagator.run_mcs() # ## Calculate phasing orbit parameters # @@ -365,7 +388,7 @@ # Execute the mission control sequence to solve for the phasing maneuver: chaser_propagator.run_mcs() -target_propagator.run_mcs() +target_satellite.propagator.run_mcs() root.rewind() # ## Retrieve the results From 93c81c4e885fd5e01d909adcc84525f7dca29a2e Mon Sep 17 00:00:00 2001 From: Jorge Martinez Garrido Date: Fri, 20 Feb 2026 10:01:18 +0100 Subject: [PATCH 05/18] wip: example --- .github/agents/astrogator.agent.md | 286 ++++++++++++++++++++++ examples/phasing-geo-orbits.py | 375 +++++++++++++++++++++-------- 2 files changed, 555 insertions(+), 106 deletions(-) create mode 100644 .github/agents/astrogator.agent.md diff --git a/.github/agents/astrogator.agent.md b/.github/agents/astrogator.agent.md new file mode 100644 index 0000000000..d1e647a2e8 --- /dev/null +++ b/.github/agents/astrogator.agent.md @@ -0,0 +1,286 @@ +--- +name: astrogator +description: Focuses on astrogator usage +--- + +You are a software developer with strong background in digital mission engineering and astrodynamics. You know how to use Ansys Systems Tool Kit (STK), and have a strong foundations about the Astrogator propagator. + +## Your Scope + +You are responsible for: +1. Creating new examples for project documentation using Astrogator +2. Ensuring examples have a consistent layout and structure +3. Reviewing existing examples for correctness and adherence to guidelines +4. Providing guidance on how to use Astrogator segments and profiles + +You do NOT: +- Review source code outside of example files +- Review CI/CD configuration or workflows +- Suggest improvements beyond example writing and correctness + +## Astrogator Overview + +Astrogator is a powerful orbit propagator in STK that allows for mission design and analysis. It uses a sequence-based approach where trajectories are built using different segment types. + +### Key Components + +1. **Main Control Sequence (MCS)**: The primary sequence containing all trajectory segments +2. **Segments**: Building blocks of a trajectory (Initial State, Propagate, Maneuver, etc.) +3. **Profiles**: Solvers that can optimize parameters within target sequences (Differential Corrector, Lambert Profile, etc.) +4. **Stopping Conditions**: Criteria that determine when propagation segments end +5. **Control Parameters**: Variables that can be adjusted by solvers +6. **Results**: Target values that solvers attempt to achieve + +## Segment Types + +### INITIAL_STATE +Defines the starting state of a spacecraft. Use element sets like: +- **Keplerian**: periapsis_radius_size, eccentricity, inclination, raan, arg_of_periapsis, true_anomaly +- **Cartesian**: x, y, z, vx, vy, vz + +```python +from ansys.stk.core.stkobjects.astrogator import SegmentType, ElementSetType + +initial_state = satellite.propagator.main_sequence.insert( + SegmentType.INITIAL_STATE, "Initial State", "-" +) +initial_state.set_element_type(ElementSetType.KEPLERIAN) +initial_state.element.periapsis_radius_size = 6700.00 +initial_state.element.eccentricity = 0.00 +initial_state.element.inclination = 0.00 +``` + +### PROPAGATE +Propagates the spacecraft forward in time using a specified propagator and stopping conditions. + +```python +propagate = satellite.propagator.main_sequence.insert( + SegmentType.PROPAGATE, "Propagate Parking Orbit", "-" +) +propagate.propagator_name = "Earth Point Mass" +propagate.stopping_conditions["Duration"].properties.trip = 7200 # seconds +propagate.properties.color = Colors.Blue +``` + +Common stopping conditions: +- `Duration`: Time-based propagation +- `Apoapsis`: Stop at orbit apoapsis +- `Periapsis`: Stop at orbit periapsis +- `Altitude`: Stop at specific altitude + +### MANEUVER +Applies velocity changes to the spacecraft. Two main types: +- **IMPULSIVE**: Instantaneous velocity change +- **FINITE**: Continuous thrust over time + +```python +from ansys.stk.core.stkobjects.astrogator import ManeuverType, AttitudeControl + +maneuver = sequence.segments.insert(SegmentType.MANEUVER, "First Impulse", "-") +maneuver.set_maneuver_type(ManeuverType.IMPULSIVE) +maneuver.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) +``` + +### SEQUENCE +Groups multiple segments together for organization. + +```python +transfer_sequence = satellite.propagator.main_sequence.insert( + SegmentType.SEQUENCE, "Hohmann Transfer", "-" +) +``` + +### TARGET_SEQUENCE +A special sequence that contains a solver profile for optimization. Used when you need to achieve specific results by adjusting control parameters. + +```python +from ansys.stk.core.stkobjects.astrogator import TargetSequenceAction + +target_seq = satellite.propagator.main_sequence.insert( + SegmentType.TARGET_SEQUENCE, "Targeting Sequence", "-" +) +target_seq.action = TargetSequenceAction.RUN_ACTIVE_PROFILES +``` + +## Profiles + +### Differential Corrector +Iteratively adjusts control parameters to achieve desired results. + +```python +from ansys.stk.core.stkobjects.astrogator import ProfileMode, ControlManeuver + +# Enable control parameter on a maneuver +maneuver.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) +maneuver.results.add("Keplerian Elems/Radius of Apoapsis") + +# Configure the differential corrector +dc = target_seq.profiles["Differential Corrector"] +dc.mode = ProfileMode.ITERATE +dc.max_iterations = 50 + +# Configure control parameter +control = dc.control_parameters.get_control_by_paths( + "First Impulse", "ImpulsiveMnvr.Cartesian.X" +) +control.enable = True +control.max_step = 0.30 + +# Configure result +result = dc.results.get_result_by_paths( + "First Impulse", "Radius Of Apoapsis" +) +result.enable = True +result.desired_value = 42238.00 +result.tolerance = 0.10 +``` + +### Lambert Profile +Solves Lambert's problem: finding the orbit connecting two position vectors in a given time. + +```python +from ansys.stk.core.stkobjects.astrogator import ( + LambertTargetCoordinateType, + LambertDirectionOfMotionType, + LambertSolutionOptionType +) + +lambert = target_seq.profiles.add("Lambert Profile") +lambert.coord_system_name = "CentralBody/Sun Inertial" +lambert.set_target_coord_type(LambertTargetCoordinateType.CARTESIAN) + +# Set target position and velocity +lambert.target_position_x = final_position[0] * 1000 +lambert.target_position_y = final_position[1] * 1000 +lambert.target_position_z = final_position[2] * 1000 + +lambert.target_velocity_x = final_velocity[0] * 1000 +lambert.target_velocity_y = final_velocity[1] * 1000 +lambert.target_velocity_z = final_velocity[2] * 1000 + +# Configure time of flight and solution parameters +lambert.time_of_flight = tof +lambert.solution_option = LambertSolutionOptionType.FIXED_TIME +lambert.revolutions = 0 +lambert.direction_of_motion = LambertDirectionOfMotionType.SHORT + +# Enable writing to segments +lambert.enable_write_to_first_maneuver = True +lambert.first_maneuver_segment = first_impulse.name +lambert.enable_write_duration_to_propagate = True +lambert.propagate_segment = propagate.name +lambert.enable_write_to_second_maneuver = True +lambert.second_maneuver_segment = last_impulse.name +``` + +## Common Transfer Patterns + +### Hohmann Transfer +Two-impulse transfer between circular orbits: +1. Initial State segment with parking orbit parameters +2. Propagate parking orbit +3. Target Sequence containing: + - First impulse to raise apoapsis + - Propagate to apoapsis +4. Target Sequence containing: + - Last impulse to circularize orbit +5. Propagate final orbit + +### Lambert Transfer +Transfers between two position vectors in specified time: +1. Initial State segment with departure state +2. Target Sequence containing: + - First impulse (optimized by Lambert profile) + - Propagate for time of flight + - Last impulse (optimized by Lambert profile) +3. Lambert Profile configured with target position, velocity, and time of flight + +### Bi-elliptic Transfer +Three-impulse transfer for large orbit changes: +1. Initial State segment with parking orbit parameters +2. Propagate parking orbit +3. Target Sequence containing: + - First impulse to raise apoapsis to intermediate value + - Propagate to apoapsis +4. Target Sequence containing: + - Second impulse to change periapsis + - Propagate to periapsis +5. Target Sequence containing: + - Last impulse to circularize orbit +6. Propagate final orbit + +## Running and Applying Results + +After configuring all segments: + +```python +from ansys.stk.core.stkobjects.astrogator import ProfilesFinish + +# Configure target sequence behavior +target_seq.when_profiles_finish = ProfilesFinish.RUN_TO_RETURN_AND_CONTINUE +target_seq.continue_on_failure = False +target_seq.reset_inner_targeters = False + +# Run the main control sequence +satellite.propagator.run_mcs() + +# Apply computed results to the trajectory +satellite.propagator.apply_all_profile_changes() +``` + +## Retrieving Results + +Access computed values from solvers: + +```python +# For differential corrector control parameters +delta_v = control_param.final_value +print(f"ΔV = {delta_v:.5f} km/s") + +# For maneuver magnitude +delta_v = maneuver.maneuver.attitude_control.magnitude +print(f"Maneuver magnitude: {delta_v:.2f} km/s") +``` + +## Example Creation Guidelines + +When creating Astrogator examples: + +1. **Setup**: Initialize STK, create scenario, add satellite with Astrogator propagator +2. **Clear sequence**: Remove all existing segments from main sequence +3. **Initial state**: Define starting orbital parameters +4. **Parking orbit**: Optional propagation before maneuvers +5. **Transfer sequences**: Build maneuver sequences with appropriate segment types +6. **Configure solvers**: Set up profiles with control parameters and results +7. **Execute**: Run MCS and apply results +8. **Retrieve**: Extract and display key results (ΔV, orbital elements) +9. **Visualize**: Update camera and show 3D trajectory + +Always enable 3D trajectory drawing: +```python +satellite.propagator.options.draw_trajectory_in_3d = True +``` + +Use color coding for different phases: +```python +from ansys.stk.core.utilities.colors import Colors + +parking_orbit.properties.color = Colors.Blue +transfer_orbit.properties.color = Colors.Red +final_orbit.properties.color = Colors.Green +``` + +## Review Criteria + +When reviewing Astrogator examples, check: + +1. **Correct segment types**: Appropriate use of INITIAL_STATE, PROPAGATE, MANEUVER, SEQUENCE, TARGET_SEQUENCE +2. **Proper element sets**: Keplerian or Cartesian elements set correctly +3. **Stopping conditions**: Correctly configured (Duration, Apoapsis, Periapsis, etc.) +4. **Solver configuration**: Differential Corrector or Lambert Profile properly set up +5. **Control parameters**: Enabled on correct maneuvers with reasonable max_step values +6. **Results**: Specified with appropriate desired values and tolerances +7. **Execution order**: run_mcs() followed by apply_all_profile_changes() +8. **Color coding**: Different phases use distinct colors for clarity +9. **Units consistency**: All values in correct units (km, km/s, seconds, degrees) +10. **Documentation**: Each section clearly explained with Markdown headings diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 1787b6f84c..9e7c6bd391 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -13,18 +13,49 @@ # - **Number of phasing orbits**: How many orbits are used to close the phase angle # - **Phasing orbit**: A temporary orbit with a different period to achieve the desired drift # +# +# The **phase angle** is the angular displacement between the two satellites at initial time. In this context, it is defined as always positive and between 0 and 360 degrees. To make the rendezvous between the satellites, it has to be reduced to 0 in a specified number of orbits (called **phasing orbits**), so each phase orbit will reduce the angular gap by just a fraction of the overall value. +# +# $$ \theta_{gap}=\frac{\theta}{n_{orbits}}$$ +# +# , where $\theta$ is the phase angle, and $\theta_{gap}$ is the angular dispacement to recover after each phasing orbit. +# +# In the code below the **Kepler's equation** is used to calculate the period and semimajor axis of the phasing orbit, given the number of revolution around the phasing orbit itself. As first, the **eccentric anomaly E** is calculated from $\theta_{gap}$ : +# +# $$\tan\left(\frac{E}{2}\right)=\sqrt{\frac{1 - e}{1 + e}}\tan\left(\frac{\theta_{gap}}{2}\right)$$ +# +# ...and then the period of the phasing orbit is derived: +# +# $$T_{phasing}=\frac{T_{chaser}}{2\pi}\left ( E-e \sin E \right )$$ +# # ## Problem statement # # Two geostationary satellites occupy the same orbit but at different angular positions. The target satellite is at a true anomaly of 20 degrees, while the chaser satellite is at 70 degrees. Design a phasing maneuver for the chaser satellite to rendezvous with the target satellite using 5 phasing orbits over a 16-day period. # # Both satellites have the following initial orbital parameters: +# # - Semi-major axis: 42164 km (GEO altitude) # - Eccentricity: 0.0001 # - Inclination: 0.0 degrees (equatorial) # - RAAN: 0.0 degrees +# - Argument of periapsis: 0.0 degrees # # Compute the required $\Delta v$ for the phasing maneuvers and determine the phasing orbit parameters. +# + +semimajor_axis = 42164.0 +eccentricity = 0.0001 +inclination = 0.0 +raan = 0.0 +arg_of_periapsis = 0.0 + +true_anomaly_target = 20.0 +true_anomaly_chaser = 70.0 + +propagation_time = 16 * 86400 +phasing_orbits = 5 +# - + # ## Launch a new STK instance # # Start by launching a new STK instance. In this example, ``STKEngine`` is used with graphics (``no_graphics`` mode set to ``False``). This means that the graphic user interface (GUI) of the product is not launched but 2D and 3D visualization is still available through the STK Engine controls: @@ -47,161 +78,197 @@ root = stk.new_object_root() root.new_scenario("Phasing_orbit") -scen = root.current_scenario +scenario = root.current_scenario +# - + +# Configure the scenario time: + +start_time, stop_time = "today", "+16 day" +scenario.set_time_period(start_time, stop_time) + +# Once created, it is possible to show a 3D graphics window by running: + +# + +from ansys.stk.core.stkengine.experimental.jupyterwidgets import GlobeWidget + +plotter = GlobeWidget(root, 640, 480) +plotter.show() # - -# ## Configure the target satellite +# ## Set up the target satellite # # Create the target satellite which will remain at a fixed position in geostationary orbit. This satellite serves as the reference point for the phasing maneuver: -# + -target_satellite = scen.children.new(STKObjectType.SATELLITE, "Target") -# - +target_satellite = scenario.children.new(STKObjectType.SATELLITE, "Target") -# Then, declare the type of orbit propagator used for the satellite: +# Then, declare the type of orbit propagator used for the satellite. Ensure a clean main sequence: -# + target_satellite.set_propagator_type(PropagatorType.ASTROGATOR) -# - +target_satellite.propagator.main_sequence.remove_all() -# Initialize the propagator and configure additional settings: +# Configure graphics settings: -target_satellite.propagator.options.draw_trajectory_in_3d = False +target_satellite.propagator.options.draw_trajectory_in_3d = True # ## Set up the initial state of the target satellite # # Access the existing initial state segment in the main sequence and configure the element type: # + -from ansys.stk.core.stkobjects.astrogator import ElementSetType +from ansys.stk.core.stkobjects.astrogator import ElementSetType, SegmentType -target_initial_state = target_satellite.propagator.main_sequence["Initial State"] +target_initial_state = target_satellite.propagator.main_sequence.insert( + SegmentType.INITIAL_STATE, "Initial State", "-" +) target_initial_state.set_element_type(ElementSetType.KEPLERIAN) -target_initial_state.initial_state.epoch = scen.start_time +target_initial_state.initial_state.epoch = scenario.start_time +# - + +# Delcare the Keplerian elements for the initial state: + +# + target_keplerian_elements = target_initial_state.element + +target_keplerian_elements.semimajor_axis = semimajor_axis +target_keplerian_elements.eccentricity = eccentricity +target_keplerian_elements.inclination = inclination +target_keplerian_elements.raan = raan +target_keplerian_elements.true_anomaly = true_anomaly_target +target_keplerian_elements.arg_of_periapsis = arg_of_periapsis # - -# Configure the propagation segment: +# Configure the propagation segment. Propagate the orbit of the target for 16 days: # + from ansys.stk.core.utilities.colors import Colors - -target_propagate_segment = target_satellite.propagator.main_sequence["Propagate"] +target_propagate_segment = target_satellite.propagator.main_sequence.insert( + SegmentType.PROPAGATE, "Propagate", "-" +) target_propagate_segment.propagator_name = "Earth point mass" + +target_propagate_segment.stopping_conditions["Duration"].properties.trip = ( + propagation_time +) +# - + +# Red color is used to identify the target satellite: + target_propagate_segment.properties.color = Colors.Red + +# Propagate the target satellite and show it: + +# + +from ansys.stk.core.stkobjects.astrogator import RunCode + + +run_code = target_satellite.propagator.run_mcs2() +if run_code != RunCode.MARCHING: + raise ValueError("Could not propagate target satellite orbit.") # - -# ## Configure the chaser satellite -# -# Create the chaser satellite which will perform the phasing maneuver to rendezvous with the target: +# Finally, show the orbit of the satellite: -chaser_satellite = scen.children.new(STKObjectType.SATELLITE, "Chaser") -chaser_satellite.set_propagator_type(PropagatorType.ASTROGATOR) -chaser_propagator = chaser_satellite.propagator -chaser_propagator.options.draw_trajectory_in_3d = False -chaser_initial_state = chaser_propagator.main_sequence["Initial State"] -chaser_initial_state.set_element_type(ElementSetType.KEPLERIAN) -chaser_initial_state.initial_state.epoch = scen.start_time -chaser_propagate_segment = chaser_propagator.main_sequence["Propagate"] -chaser_propagate_segment.propagator_name = "Earth point mass" -chaser_propagate_segment.stopping_conditions.add("Periapsis") -chaser_propagate_segment.stopping_conditions.remove("Duration") -chaser_keplerian_elements = chaser_initial_state.element +plotter.show() -# ## Define the input parameters +# ## Setup the chaser satellite # -# The cell below defines the input parameters for the phasing maneuver analysis: +# Create the chaser satellite which will get closer to the target satellite: -import math +chaser_satellite = scenario.children.new(STKObjectType.SATELLITE, "Chaser") +# Then, declare the type of orbit propagator used for the satellite. Ensure a clean main sequence: -# Input parameters for phasing maneuver -chaser_true_anomaly = 70 # deg - chaser true anomaly (between 0 and 360) -target_true_anomaly = 20 # deg - target true anomaly (between 0 and 360) -number_of_phasing_orbits = 5 # number of phasing orbits -propagation_time = 16 # days - total propagation time +chaser_satellite.set_propagator_type(PropagatorType.ASTROGATOR) +chaser_satellite.propagator.main_sequence.remove_all() -# Fixed Keplerian parameters (geostationary orbit) -semi_major_axis = 42164 # km -eccentricity = 0.0001 -inclination = 0.0 -right_ascension_ascending_node = 0.0 +# Configure graphics settings: + +chaser_satellite.propagator.options.draw_trajectory_in_3d = True -# ## Configure initial orbital states +# ## Setup the initial state of the chaser satellite # -# Set the initial orbital elements for both satellites and run an initial propagation: +# Access the existing initial state segment in the main sequence and configure the element type: -# Change the scenario duration -scen.stop_time = "+" + str(propagation_time) + " day" +# + +from ansys.stk.core.stkobjects.astrogator import ElementSetType, SegmentType + + +chaser_initial_state = chaser_satellite.propagator.main_sequence.insert( + SegmentType.INITIAL_STATE, "Initial State", "-" +) +chaser_initial_state.set_element_type(ElementSetType.KEPLERIAN) +chaser_initial_state.initial_state.epoch = scenario.start_time +# - + +# Declare the keplerian elements: + +# + +chaser_keplerian_elements = chaser_initial_state.element -# Set the chaser initial state -chaser_keplerian_elements.semimajor_axis = semi_major_axis +chaser_keplerian_elements.semimajor_axis = semimajor_axis chaser_keplerian_elements.eccentricity = eccentricity chaser_keplerian_elements.inclination = inclination -chaser_keplerian_elements.raan = right_ascension_ascending_node -chaser_keplerian_elements.true_anomaly = chaser_true_anomaly -chaser_keplerian_elements.arg_of_periapsis = 0.0 +chaser_keplerian_elements.raan = raan +chaser_keplerian_elements.true_anomaly = true_anomaly_chaser +chaser_keplerian_elements.arg_of_periapsis = arg_of_periapsis +# - -# Set the target initial state -target_keplerian_elements.semimajor_axis = semi_major_axis -target_keplerian_elements.eccentricity = eccentricity -target_keplerian_elements.inclination = inclination -target_keplerian_elements.raan = right_ascension_ascending_node -target_keplerian_elements.true_anomaly = target_true_anomaly -target_keplerian_elements.arg_of_periapsis = 0.0 +# Configure the propagation segment. Propagate the orbit of the chaser for 16 days: -# Configure target propagation duration -target_propagate_segment.stopping_conditions["Duration"].properties.trip = ( - propagation_time * 86400 +# + +chaser_propagate_segment = chaser_satellite.propagator.main_sequence.insert( + SegmentType.PROPAGATE, "Propagate", "-" ) +chaser_propagate_segment.propagator_name = "Earth point mass" -# Run initial propagation for both satellites -chaser_propagator.run_mcs() -target_satellite.propagator.run_mcs() -# ## Calculate phasing orbit parameters -# -# The **phase angle** is the angular displacement between the two satellites at initial time. To achieve rendezvous, this angle must be reduced to 0 over a specified number of phasing orbits. -# -# The angular gap per orbit is: -# $$ \theta_{gap}=\frac{\theta}{n_{orbits}}$$ -# -# where $\theta$ is the phase angle, and $\theta_{gap}$ is the angular displacement to recover after each phasing orbit. -# -# Using **Kepler's equation**, we calculate the period and semi-major axis of the phasing orbit. First, the **eccentric anomaly E** is calculated from $\theta_{gap}$: -# -# $$\tan\frac{E}{2} = \sqrt{\frac{1 - e}{1 + e}}\tan\frac{\theta_{gap}}{2}$$ -# -# Then the period of the phasing orbit is derived: +chaser_propagate_segment.stopping_conditions.add("Periapsis") +chaser_propagate_segment.stopping_conditions.remove("Duration") +# - + +# Run the mission control sequence of the chaser: + +run_code = chaser_satellite.propagator.run_mcs2() +if run_code != RunCode.MARCHING: + raise ValueError("Could not propagate chaser satellite orbit.") + +# Finally, show the orbit of the satellite: + +plotter.show() + +# ## Calculate phasing parameters # -# $$ T_{phasing}=\frac{T_{chaser}}{2\pi}\left ( E-e \sin E \right )$$ +# Now calculate the parameters for the phasing maneuver using Kepler's equation and orbital mechanics principles: + +# + +import math # Get the mean period and SMA of the initial orbit mean_keplerian_data_provider = chaser_satellite.data_providers[ "Kozai-Izsak Mean" ].group["ICRF"] mean_keplerian_result = mean_keplerian_data_provider.execute( - scen.start_time, scen.start_time, 60 + scenario.start_time, scenario.start_time, 60 ) mean_keplerian_dataframe = mean_keplerian_result.data_sets.to_pandas_dataframe() mean_orbital_period = float(mean_keplerian_dataframe.at[0, "mean nodal period"]) mean_semi_major_axis = float(mean_keplerian_dataframe.at[0, "mean semi-major axis"]) # Calculate the phase angle -phase_angle = math.radians(target_true_anomaly - chaser_true_anomaly) +phase_angle = math.radians(true_anomaly_target - true_anomaly_chaser) if phase_angle < 0: phase_angle = 2 * math.pi + phase_angle print("#################### Initial geometry #####################") -print("Phase angle = " + str(math.degrees(phase_angle)) + " deg") +print(f"Phase angle = {math.degrees(phase_angle):.1f} deg") print("") # Recalculate the phase angle accordingly with the number of phasing orbits -phase_angle = phase_angle / number_of_phasing_orbits +phase_angle = phase_angle / phasing_orbits # Calculate the eccentric anomaly eccentric_anomaly = 2 * math.atan( @@ -228,12 +295,12 @@ "Astrogator Values" ].group["Keplerian Elems"] keplerian_result = keplerian_elements_data_provider.execute( - scen.start_time, scen.start_time, 60 + scenario.start_time, scenario.start_time, 60 ) keplerian_dataframe = keplerian_result.data_sets.to_pandas_dataframe() initial_radius = float(keplerian_dataframe.at[0, "radius_of_periapsis"]) -if transfer_semi_major_axis > semi_major_axis: +if transfer_semi_major_axis > semimajor_axis: periapsis_radius = initial_radius apoapsis_radius = 2 * transfer_semi_major_axis - periapsis_radius else: @@ -249,23 +316,23 @@ * ((2 / initial_radius) - (1 / transfer_semi_major_axis)) ) delta_v_estimate = velocity_periapsis_final - velocity_periapsis_initial +# - # ## Set up the phasing maneuver sequence # # Add a target sequence to perform the first delta-V maneuver and propagate through the phasing orbits: +# + from ansys.stk.core.stkobjects.astrogator import ( AttitudeControl, ControlManeuver, ManeuverType, ProfileMode, - SegmentType, TargetSequenceAction, ) - # Add the first Target Sequence segment -phasing_start_sequence = chaser_propagator.main_sequence.insert( +phasing_start_sequence = chaser_satellite.propagator.main_sequence.insert( SegmentType.TARGET_SEQUENCE, "Start Phasing", "-" ) phasing_start_sequence.action = TargetSequenceAction.RUN_ACTIVE_PROFILES @@ -289,17 +356,18 @@ ) phasing_orbit_propagate.properties.color = Colors.Orange phasing_orbit_propagate.propagator_name = "Earth point mass" -if transfer_semi_major_axis > semi_major_axis: +if transfer_semi_major_axis > semimajor_axis: phasing_orbit_propagate.stopping_conditions.add("Periapsis") else: phasing_orbit_propagate.stopping_conditions.add("Apoapsis") phasing_orbit_propagate.stopping_conditions.remove("Duration") phasing_orbit_propagate.stopping_conditions.item( 0 -).properties.repeat_count = number_of_phasing_orbits +).properties.repeat_count = phasing_orbits phasing_orbit_propagate.results.add("Vector/Angle Between Vectors") phasing_orbit_propagate.results[0].vector1_name = "Satellite/Chaser Position" phasing_orbit_propagate.results[0].vector2_name = "Satellite/Target Position" +# - # ## Configure the differential corrector for phasing # @@ -333,7 +401,7 @@ # Add a second target sequence to perform the circularization delta-V to return to the original circular orbit: # Add the second Target Sequence segment -circularization_sequence = chaser_propagator.main_sequence.insert( +circularization_sequence = chaser_satellite.propagator.main_sequence.insert( SegmentType.TARGET_SEQUENCE, "Circularization", "-" ) circularization_sequence.action = TargetSequenceAction.RUN_ACTIVE_PROFILES @@ -358,6 +426,7 @@ # # Set up the differential corrector to achieve a circular orbit with the desired eccentricity: +# + # Customize the Differential Corrector circularization_differential_corrector = circularization_sequence.profiles[ "Differential Corrector" @@ -382,12 +451,13 @@ circularization_eccentricity_result.enable = True circularization_eccentricity_result.desired_value = 0.0001 circularization_eccentricity_result.tolerance = 0.00001 +# - # ## Run the main control sequence # # Execute the mission control sequence to solve for the phasing maneuver: -chaser_propagator.run_mcs() +chaser_satellite.propagator.run_mcs() target_satellite.propagator.run_mcs() root.rewind() @@ -395,21 +465,114 @@ # # Once the analysis has been performed, retrieve the delta-V values and phasing orbit parameters: +# + # Get the maneuver data providers maneuver_data_provider = chaser_satellite.data_providers["Maneuver Summary"] -maneuver_result = maneuver_data_provider.execute(scen.start_time, scen.stop_time) +maneuver_result = maneuver_data_provider.execute(scenario.start_time, scenario.stop_time) maneuver_dataframe = maneuver_result.data_sets.to_pandas_dataframe() delta_v_actual = maneuver_dataframe.at[0, "delta v"] +print("") print("################### Transfer orbit data ###################") -print("N phasing orbits = " + str(number_of_phasing_orbits)) -print("Period = " + str(phasing_orbit_period) + " sec") -print("SMA = " + str(transfer_semi_major_axis) + " km") -print("Perigee radius = " + str(periapsis_radius) + " km") -print("Apogee radius = " + str(apoapsis_radius) + " km") -print("Delta V = " + str(delta_v_actual) + " m/sec") -print( - "Transfer time = " - + str(phasing_orbit_period * number_of_phasing_orbits / 86400) - + " days" -) +print(f"N phasing orbits = {phasing_orbits}") +print(f"Period = {phasing_orbit_period:.2f} sec") +print(f"SMA = {transfer_semi_major_axis:.2f} km") +print(f"Perigee radius = {periapsis_radius:.2f} km") +print(f"Apogee radius = {apoapsis_radius:.2f} km") +print(f"Delta V = {delta_v_actual} m/sec") +print(f"Transfer time = {phasing_orbit_period * phasing_orbits / 86400:.4f} days") +print("") +# - + +# Finally, show the complete orbit trajectory: + +plotter.show() + +# ## Plot the phase angle history +# +# Retrieve and plot the phase angle between the chaser and target satellites over time: + +# + +import matplotlib.pyplot as plt +import numpy as np + +# Get the phase angle data over time +angle_data_provider = chaser_satellite.data_providers["Vector"].group["Angle Between Vectors"] +angle_result = angle_data_provider.execute(scenario.start_time, scenario.stop_time, 3600) +angle_dataframe = angle_result.data_sets.to_pandas_dataframe() + +# Convert time to days +time_seconds = angle_dataframe["Time"].values +time_days = time_seconds / 86400.0 + +# Get the angle values +phase_angles = angle_dataframe["Angle"].values + +# Create the plot +plt.figure(figsize=(10, 6)) +plt.plot(time_days, phase_angles, 'b-', linewidth=2) +plt.xlabel('Time (days)', fontsize=12) +plt.ylabel('Phase Angle (degrees)', fontsize=12) +plt.title('Phase Angle History During Phasing Maneuver', fontsize=14, fontweight='bold') +plt.grid(True, alpha=0.3) +plt.xlim([0, propagation_time / 86400]) +plt.tight_layout() +plt.show() +# - + +# ## Plot the orbital altitude history +# +# Retrieve and plot the altitude of both satellites to visualize the phasing orbit: + +# + +# Get altitude data for chaser satellite +chaser_altitude_provider = chaser_satellite.data_providers["Cartesian Position"].group["Fixed"] +chaser_altitude_result = chaser_altitude_provider.execute(scenario.start_time, scenario.stop_time, 3600) +chaser_altitude_df = chaser_altitude_result.data_sets.to_pandas_dataframe() + +# Calculate magnitude of position vector and subtract Earth radius +earth_radius = 6378.137 # km +chaser_x = chaser_altitude_df["x"].values +chaser_y = chaser_altitude_df["y"].values +chaser_z = chaser_altitude_df["z"].values +chaser_altitude = np.sqrt(chaser_x**2 + chaser_y**2 + chaser_z**2) - earth_radius + +# Get altitude data for target satellite +target_altitude_provider = target_satellite.data_providers["Cartesian Position"].group["Fixed"] +target_altitude_result = target_altitude_provider.execute(scenario.start_time, scenario.stop_time, 3600) +target_altitude_df = target_altitude_result.data_sets.to_pandas_dataframe() + +target_x = target_altitude_df["x"].values +target_y = target_altitude_df["y"].values +target_z = target_altitude_df["z"].values +target_altitude = np.sqrt(target_x**2 + target_y**2 + target_z**2) - earth_radius + +# Convert time to days +chaser_time_days = chaser_altitude_df["Time"].values / 86400.0 +target_time_days = target_altitude_df["Time"].values / 86400.0 + +# Create the plot +plt.figure(figsize=(10, 6)) +plt.plot(chaser_time_days, chaser_altitude, 'b-', linewidth=2, label='Chaser') +plt.plot(target_time_days, target_altitude, 'r--', linewidth=2, label='Target') +plt.xlabel('Time (days)', fontsize=12) +plt.ylabel('Altitude (km)', fontsize=12) +plt.title('Altitude History During Phasing Maneuver', fontsize=14, fontweight='bold') +plt.legend(fontsize=11) +plt.grid(True, alpha=0.3) +plt.xlim([0, propagation_time / 86400]) +plt.tight_layout() +plt.show() +# - + +# ## Summary +# +# This example demonstrates how to: +# +# 1. Set up two geostationary satellites with different true anomalies +# 2. Calculate the required phasing orbit parameters using Kepler's equation +# 3. Design a phasing maneuver sequence with differential corrector profiles +# 4. Execute the maneuver and retrieve the results +# 5. Visualize the phase angle and altitude histories +# +# The phasing maneuver successfully brings the chaser satellite to the same orbital position as the target satellite after the specified number of phasing orbits. From feb762801d7ace55e9cf785f0faa2bb936617ae7 Mon Sep 17 00:00:00 2001 From: pyansys-ci-bot <92810346+pyansys-ci-bot@users.noreply.github.com> Date: Fri, 20 Feb 2026 09:02:46 +0000 Subject: [PATCH 06/18] chore: adding changelog file 950.documentation.md [dependabot-skip] --- doc/source/changelog/950.documentation.md | 1 + 1 file changed, 1 insertion(+) create mode 100644 doc/source/changelog/950.documentation.md diff --git a/doc/source/changelog/950.documentation.md b/doc/source/changelog/950.documentation.md new file mode 100644 index 0000000000..2b49ebb301 --- /dev/null +++ b/doc/source/changelog/950.documentation.md @@ -0,0 +1 @@ +Geosynchronous phasing example From f55361917c511b287eb767617cd956da26a7a6f4 Mon Sep 17 00:00:00 2001 From: Jorge Martinez Date: Fri, 27 Feb 2026 10:48:31 +0100 Subject: [PATCH 07/18] dbg: check push --- file | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 file diff --git a/file b/file new file mode 100644 index 0000000000..e69de29bb2 From 90cafc17f34f3047b5ddf3a51a4af008ad6d9875 Mon Sep 17 00:00:00 2001 From: Dupont Date: Wed, 8 Jul 2026 16:40:57 -0400 Subject: [PATCH 08/18] doc: working on example --- .github/agents/astrogator.agent.md | 286 ------------------------ examples/phasing-geo-orbits.py | 344 ++++++++++++++++++++--------- 2 files changed, 236 insertions(+), 394 deletions(-) delete mode 100644 .github/agents/astrogator.agent.md diff --git a/.github/agents/astrogator.agent.md b/.github/agents/astrogator.agent.md deleted file mode 100644 index d1e647a2e8..0000000000 --- a/.github/agents/astrogator.agent.md +++ /dev/null @@ -1,286 +0,0 @@ ---- -name: astrogator -description: Focuses on astrogator usage ---- - -You are a software developer with strong background in digital mission engineering and astrodynamics. You know how to use Ansys Systems Tool Kit (STK), and have a strong foundations about the Astrogator propagator. - -## Your Scope - -You are responsible for: -1. Creating new examples for project documentation using Astrogator -2. Ensuring examples have a consistent layout and structure -3. Reviewing existing examples for correctness and adherence to guidelines -4. Providing guidance on how to use Astrogator segments and profiles - -You do NOT: -- Review source code outside of example files -- Review CI/CD configuration or workflows -- Suggest improvements beyond example writing and correctness - -## Astrogator Overview - -Astrogator is a powerful orbit propagator in STK that allows for mission design and analysis. It uses a sequence-based approach where trajectories are built using different segment types. - -### Key Components - -1. **Main Control Sequence (MCS)**: The primary sequence containing all trajectory segments -2. **Segments**: Building blocks of a trajectory (Initial State, Propagate, Maneuver, etc.) -3. **Profiles**: Solvers that can optimize parameters within target sequences (Differential Corrector, Lambert Profile, etc.) -4. **Stopping Conditions**: Criteria that determine when propagation segments end -5. **Control Parameters**: Variables that can be adjusted by solvers -6. **Results**: Target values that solvers attempt to achieve - -## Segment Types - -### INITIAL_STATE -Defines the starting state of a spacecraft. Use element sets like: -- **Keplerian**: periapsis_radius_size, eccentricity, inclination, raan, arg_of_periapsis, true_anomaly -- **Cartesian**: x, y, z, vx, vy, vz - -```python -from ansys.stk.core.stkobjects.astrogator import SegmentType, ElementSetType - -initial_state = satellite.propagator.main_sequence.insert( - SegmentType.INITIAL_STATE, "Initial State", "-" -) -initial_state.set_element_type(ElementSetType.KEPLERIAN) -initial_state.element.periapsis_radius_size = 6700.00 -initial_state.element.eccentricity = 0.00 -initial_state.element.inclination = 0.00 -``` - -### PROPAGATE -Propagates the spacecraft forward in time using a specified propagator and stopping conditions. - -```python -propagate = satellite.propagator.main_sequence.insert( - SegmentType.PROPAGATE, "Propagate Parking Orbit", "-" -) -propagate.propagator_name = "Earth Point Mass" -propagate.stopping_conditions["Duration"].properties.trip = 7200 # seconds -propagate.properties.color = Colors.Blue -``` - -Common stopping conditions: -- `Duration`: Time-based propagation -- `Apoapsis`: Stop at orbit apoapsis -- `Periapsis`: Stop at orbit periapsis -- `Altitude`: Stop at specific altitude - -### MANEUVER -Applies velocity changes to the spacecraft. Two main types: -- **IMPULSIVE**: Instantaneous velocity change -- **FINITE**: Continuous thrust over time - -```python -from ansys.stk.core.stkobjects.astrogator import ManeuverType, AttitudeControl - -maneuver = sequence.segments.insert(SegmentType.MANEUVER, "First Impulse", "-") -maneuver.set_maneuver_type(ManeuverType.IMPULSIVE) -maneuver.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) -``` - -### SEQUENCE -Groups multiple segments together for organization. - -```python -transfer_sequence = satellite.propagator.main_sequence.insert( - SegmentType.SEQUENCE, "Hohmann Transfer", "-" -) -``` - -### TARGET_SEQUENCE -A special sequence that contains a solver profile for optimization. Used when you need to achieve specific results by adjusting control parameters. - -```python -from ansys.stk.core.stkobjects.astrogator import TargetSequenceAction - -target_seq = satellite.propagator.main_sequence.insert( - SegmentType.TARGET_SEQUENCE, "Targeting Sequence", "-" -) -target_seq.action = TargetSequenceAction.RUN_ACTIVE_PROFILES -``` - -## Profiles - -### Differential Corrector -Iteratively adjusts control parameters to achieve desired results. - -```python -from ansys.stk.core.stkobjects.astrogator import ProfileMode, ControlManeuver - -# Enable control parameter on a maneuver -maneuver.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) -maneuver.results.add("Keplerian Elems/Radius of Apoapsis") - -# Configure the differential corrector -dc = target_seq.profiles["Differential Corrector"] -dc.mode = ProfileMode.ITERATE -dc.max_iterations = 50 - -# Configure control parameter -control = dc.control_parameters.get_control_by_paths( - "First Impulse", "ImpulsiveMnvr.Cartesian.X" -) -control.enable = True -control.max_step = 0.30 - -# Configure result -result = dc.results.get_result_by_paths( - "First Impulse", "Radius Of Apoapsis" -) -result.enable = True -result.desired_value = 42238.00 -result.tolerance = 0.10 -``` - -### Lambert Profile -Solves Lambert's problem: finding the orbit connecting two position vectors in a given time. - -```python -from ansys.stk.core.stkobjects.astrogator import ( - LambertTargetCoordinateType, - LambertDirectionOfMotionType, - LambertSolutionOptionType -) - -lambert = target_seq.profiles.add("Lambert Profile") -lambert.coord_system_name = "CentralBody/Sun Inertial" -lambert.set_target_coord_type(LambertTargetCoordinateType.CARTESIAN) - -# Set target position and velocity -lambert.target_position_x = final_position[0] * 1000 -lambert.target_position_y = final_position[1] * 1000 -lambert.target_position_z = final_position[2] * 1000 - -lambert.target_velocity_x = final_velocity[0] * 1000 -lambert.target_velocity_y = final_velocity[1] * 1000 -lambert.target_velocity_z = final_velocity[2] * 1000 - -# Configure time of flight and solution parameters -lambert.time_of_flight = tof -lambert.solution_option = LambertSolutionOptionType.FIXED_TIME -lambert.revolutions = 0 -lambert.direction_of_motion = LambertDirectionOfMotionType.SHORT - -# Enable writing to segments -lambert.enable_write_to_first_maneuver = True -lambert.first_maneuver_segment = first_impulse.name -lambert.enable_write_duration_to_propagate = True -lambert.propagate_segment = propagate.name -lambert.enable_write_to_second_maneuver = True -lambert.second_maneuver_segment = last_impulse.name -``` - -## Common Transfer Patterns - -### Hohmann Transfer -Two-impulse transfer between circular orbits: -1. Initial State segment with parking orbit parameters -2. Propagate parking orbit -3. Target Sequence containing: - - First impulse to raise apoapsis - - Propagate to apoapsis -4. Target Sequence containing: - - Last impulse to circularize orbit -5. Propagate final orbit - -### Lambert Transfer -Transfers between two position vectors in specified time: -1. Initial State segment with departure state -2. Target Sequence containing: - - First impulse (optimized by Lambert profile) - - Propagate for time of flight - - Last impulse (optimized by Lambert profile) -3. Lambert Profile configured with target position, velocity, and time of flight - -### Bi-elliptic Transfer -Three-impulse transfer for large orbit changes: -1. Initial State segment with parking orbit parameters -2. Propagate parking orbit -3. Target Sequence containing: - - First impulse to raise apoapsis to intermediate value - - Propagate to apoapsis -4. Target Sequence containing: - - Second impulse to change periapsis - - Propagate to periapsis -5. Target Sequence containing: - - Last impulse to circularize orbit -6. Propagate final orbit - -## Running and Applying Results - -After configuring all segments: - -```python -from ansys.stk.core.stkobjects.astrogator import ProfilesFinish - -# Configure target sequence behavior -target_seq.when_profiles_finish = ProfilesFinish.RUN_TO_RETURN_AND_CONTINUE -target_seq.continue_on_failure = False -target_seq.reset_inner_targeters = False - -# Run the main control sequence -satellite.propagator.run_mcs() - -# Apply computed results to the trajectory -satellite.propagator.apply_all_profile_changes() -``` - -## Retrieving Results - -Access computed values from solvers: - -```python -# For differential corrector control parameters -delta_v = control_param.final_value -print(f"ΔV = {delta_v:.5f} km/s") - -# For maneuver magnitude -delta_v = maneuver.maneuver.attitude_control.magnitude -print(f"Maneuver magnitude: {delta_v:.2f} km/s") -``` - -## Example Creation Guidelines - -When creating Astrogator examples: - -1. **Setup**: Initialize STK, create scenario, add satellite with Astrogator propagator -2. **Clear sequence**: Remove all existing segments from main sequence -3. **Initial state**: Define starting orbital parameters -4. **Parking orbit**: Optional propagation before maneuvers -5. **Transfer sequences**: Build maneuver sequences with appropriate segment types -6. **Configure solvers**: Set up profiles with control parameters and results -7. **Execute**: Run MCS and apply results -8. **Retrieve**: Extract and display key results (ΔV, orbital elements) -9. **Visualize**: Update camera and show 3D trajectory - -Always enable 3D trajectory drawing: -```python -satellite.propagator.options.draw_trajectory_in_3d = True -``` - -Use color coding for different phases: -```python -from ansys.stk.core.utilities.colors import Colors - -parking_orbit.properties.color = Colors.Blue -transfer_orbit.properties.color = Colors.Red -final_orbit.properties.color = Colors.Green -``` - -## Review Criteria - -When reviewing Astrogator examples, check: - -1. **Correct segment types**: Appropriate use of INITIAL_STATE, PROPAGATE, MANEUVER, SEQUENCE, TARGET_SEQUENCE -2. **Proper element sets**: Keplerian or Cartesian elements set correctly -3. **Stopping conditions**: Correctly configured (Duration, Apoapsis, Periapsis, etc.) -4. **Solver configuration**: Differential Corrector or Lambert Profile properly set up -5. **Control parameters**: Enabled on correct maneuvers with reasonable max_step values -6. **Results**: Specified with appropriate desired values and tolerances -7. **Execution order**: run_mcs() followed by apply_all_profile_changes() -8. **Color coding**: Different phases use distinct colors for clarity -9. **Units consistency**: All values in correct units (km, km/s, seconds, degrees) -10. **Documentation**: Each section clearly explained with Markdown headings diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 9e7c6bd391..b0be215bcb 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -42,12 +42,14 @@ # # Compute the required $\Delta v$ for the phasing maneuvers and determine the phasing orbit parameters. +# First, define the initial orbital parameters for use during calculations: + # + semimajor_axis = 42164.0 eccentricity = 0.0001 inclination = 0.0 raan = 0.0 -arg_of_periapsis = 0.0 +arg_of_periapsis = 0.0 true_anomaly_target = 20.0 true_anomaly_chaser = 70.0 @@ -73,7 +75,7 @@ # Start by creating a new scenario in STK: # + -from ansys.stk.core.stkobjects import PropagatorType, STKObjectType +from ansys.stk.core.stkobjects import STKObjectType root = stk.new_object_root() @@ -81,12 +83,12 @@ scenario = root.current_scenario # - -# Configure the scenario time: +# Configure the scenario time to cover a 16 day period starting July 1, 2026: -start_time, stop_time = "today", "+16 day" +start_time, stop_time = "1 Jul 2026 16:00:00.000", "+16 day" scenario.set_time_period(start_time, stop_time) -# Once created, it is possible to show a 3D graphics window by running: +# It is now possible to show a 3D graphics window by running: # + from ansys.stk.core.stkengine.experimental.jupyterwidgets import GlobeWidget @@ -96,24 +98,87 @@ plotter.show() # - -# ## Set up the target satellite -# +# ## Create the satellites + # Create the target satellite which will remain at a fixed position in geostationary orbit. This satellite serves as the reference point for the phasing maneuver: target_satellite = scenario.children.new(STKObjectType.SATELLITE, "Target") -# Then, declare the type of orbit propagator used for the satellite. Ensure a clean main sequence: +# Create the chaser satellite, which will need to approach the target satellite: + +chaser_satellite = scenario.children.new(STKObjectType.SATELLITE, "Chaser") + +# ## Set up vectors in analysis workbench + +# Analysis workbench (AWB) is used in this example to calculate the phase angle between the two satellites. To begin using analysis workbench, first retrieve the AWB component root: + +awb_root = root.analysis_workbench_components_root + +# The phase angle component will be associated with the Earth central body, so retrieve the factory that allows creating AWB components on this body: + +earth_angle_factory = awb_root.get_provider("CentralBody/Earth").angles.factory + +# Create a dihedral angle component, which is a type of angle between two vectors and about a vector: + +# + +from ansys.stk.core.analysis_workbench import AngleType + + +phase_angle_awb = earth_angle_factory.create("Phase", "Phase", AngleType.DIHEDRAL_ANGLE) +# - + +# Next, retrieve all of the vectors associated with the Earth, the target satellite, and the chaser satellite: + +target_vectors = awb_root.get_provider(target_satellite.class_name+"/"+target_satellite.instance_name).vectors +chaser_vectors = awb_root.get_provider(chaser_satellite.class_name+"/"+chaser_satellite.instance_name).vectors + +# Retrieve the target satellite's position vector, which is automatically created on STK objects when they are created. Then, designate this vector as the vector that the dihedral angle sweeps to: + +target_position_vector = target_vectors.item("Position") +phase_angle_awb.to_vector.set_vector(target_position_vector) + +# Retrieve the chaser satellite's position vector and set it as the vector that the dihedral angle begins from: + +chaser_position_vector = chaser_vectors.item("Position") +phase_angle_awb.from_vector.set_vector(chaser_position_vector) + +# Then, retrieve the target satellite's orbit normal vector, which is perpendicular to its orbit with respect to the Earth's center. Set this vector as the vector that the angle sweeps about: + +about_vector = target_vectors.item("Orbit_Normal") +phase_angle_awb.pole_about.set_vector(about_vector) + +# Configure the angle to have a signed value: + +phase_angle_awb.signed_angle = True + +# Finally, create a scalar calculation using the dihedral angle as input in order to record the angle values over time: + +earth_calc_scalar_factory = awb_root.get_provider("CentralBody/Earth").calculation_scalars.factory +phase_angle_calc = earth_calc_scalar_factory.create_angle("PhaseAngle", "PhaseAngle") +phase_angle_calc.input_angle = phase_angle_awb + +# ## Set up the target satellite + +# Then, declare the type of orbit propagator used for the satellite. Set the propagator to an Astrogator type: + +# + +from ansys.stk.core.stkobjects import PropagatorType + target_satellite.set_propagator_type(PropagatorType.ASTROGATOR) +# - + +# Ensure a clean main sequence: + target_satellite.propagator.main_sequence.remove_all() -# Configure graphics settings: +# Then, configure graphics settings to allow 3D visualization of the satellite's trajectory: target_satellite.propagator.options.draw_trajectory_in_3d = True # ## Set up the initial state of the target satellite # -# Access the existing initial state segment in the main sequence and configure the element type: +# Access the existing initial state segment in the main sequence and configure the element type to a Keplerian element: # + from ansys.stk.core.stkobjects.astrogator import ElementSetType, SegmentType @@ -126,7 +191,7 @@ target_initial_state.initial_state.epoch = scenario.start_time # - -# Delcare the Keplerian elements for the initial state: +# Delcare the Keplerian elements for the initial state using the previously defined parameters: # + target_keplerian_elements = target_initial_state.element @@ -142,8 +207,6 @@ # Configure the propagation segment. Propagate the orbit of the target for 16 days: # + -from ansys.stk.core.utilities.colors import Colors - target_propagate_segment = target_satellite.propagator.main_sequence.insert( SegmentType.PROPAGATE, "Propagate", "-" ) @@ -154,11 +217,16 @@ ) # - -# Red color is used to identify the target satellite: +# Use the color red to identify the target satellite: + +# + +from ansys.stk.core.utilities.colors import Colors + target_propagate_segment.properties.color = Colors.Red +# - -# Propagate the target satellite and show it: +# Propagate the target satellite and verify that the propagation worked correctly: # + from ansys.stk.core.stkobjects.astrogator import RunCode @@ -169,41 +237,32 @@ raise ValueError("Could not propagate target satellite orbit.") # - -# Finally, show the orbit of the satellite: +# Finally, show the orbit of the satellite in the 3D graphics window: plotter.show() -# ## Setup the chaser satellite -# -# Create the chaser satellite which will get closer to the target satellite: - -chaser_satellite = scenario.children.new(STKObjectType.SATELLITE, "Chaser") +# ## Set up the chaser satellite -# Then, declare the type of orbit propagator used for the satellite. Ensure a clean main sequence: +# Then, declare the type of orbit propagator used for the satellite to an Astrogator type and ensure a clean main sequence: chaser_satellite.set_propagator_type(PropagatorType.ASTROGATOR) chaser_satellite.propagator.main_sequence.remove_all() -# Configure graphics settings: +# Configure graphics settings to allow 3D trajectory visualization: -chaser_satellite.propagator.options.draw_trajectory_in_3d = True +chaser_satellite.propagator.options.draw_trajectory_in_3d = False -# ## Setup the initial state of the chaser satellite +# ## Set up the initial state of the chaser satellite # -# Access the existing initial state segment in the main sequence and configure the element type: - -# + -from ansys.stk.core.stkobjects.astrogator import ElementSetType, SegmentType - +# Access the existing initial state segment in the main sequence and configure the element type to a Keplerian type: chaser_initial_state = chaser_satellite.propagator.main_sequence.insert( SegmentType.INITIAL_STATE, "Initial State", "-" ) chaser_initial_state.set_element_type(ElementSetType.KEPLERIAN) chaser_initial_state.initial_state.epoch = scenario.start_time -# - -# Declare the keplerian elements: +# Declare the keplerian elements using the previously defined parameters: # + chaser_keplerian_elements = chaser_initial_state.element @@ -216,81 +275,92 @@ chaser_keplerian_elements.arg_of_periapsis = arg_of_periapsis # - -# Configure the propagation segment. Propagate the orbit of the chaser for 16 days: +# Configure the propagation segment to propagate the orbit of the chaser for 16 days: -# + chaser_propagate_segment = chaser_satellite.propagator.main_sequence.insert( SegmentType.PROPAGATE, "Propagate", "-" ) chaser_propagate_segment.propagator_name = "Earth point mass" - - chaser_propagate_segment.stopping_conditions.add("Periapsis") chaser_propagate_segment.stopping_conditions.remove("Duration") -# - -# Run the mission control sequence of the chaser: +# Run the mission control sequence of the chaser and verify that the propagation was successful: run_code = chaser_satellite.propagator.run_mcs2() if run_code != RunCode.MARCHING: raise ValueError("Could not propagate chaser satellite orbit.") -# Finally, show the orbit of the satellite: +# Finally, show the orbit of the satellite using the 3D graphics window: plotter.show() # ## Calculate phasing parameters # -# Now calculate the parameters for the phasing maneuver using Kepler's equation and orbital mechanics principles: +# Now, calculate the parameters for the phasing maneuver using Kepler's equation and orbital mechanics principles. -# + -import math +# First, retrieve the mean period and semi-major axis (SMA) of the chaser satellite's initial orbit. To do so, use the Kozai-Izsak Mean data provider associated with the satellite. Use an ICRF (International Celestial Reference Frame) coordinate system: -# Get the mean period and SMA of the initial orbit mean_keplerian_data_provider = chaser_satellite.data_providers[ "Kozai-Izsak Mean" ].group["ICRF"] + +# Compute the data every 60 seconds over the scenario, then convert the result to a Pandas dataframe: + mean_keplerian_result = mean_keplerian_data_provider.execute( scenario.start_time, scenario.start_time, 60 ) mean_keplerian_dataframe = mean_keplerian_result.data_sets.to_pandas_dataframe() + +# Finally, retrieve the initial mean period and SMA: + mean_orbital_period = float(mean_keplerian_dataframe.at[0, "mean nodal period"]) mean_semi_major_axis = float(mean_keplerian_dataframe.at[0, "mean semi-major axis"]) -# Calculate the phase angle +# Next, calculate the phase angle, which is defined as the difference between the target and chaser orbit anomalies: + +# + +import math + + phase_angle = math.radians(true_anomaly_target - true_anomaly_chaser) if phase_angle < 0: phase_angle = 2 * math.pi + phase_angle -print("#################### Initial geometry #####################") -print(f"Phase angle = {math.degrees(phase_angle):.1f} deg") +print(f"The initial phase angle is {math.degrees(phase_angle):.1f} deg.") print("") +# - + +# Then, calculate the phase angle per phasing orbit: -# Recalculate the phase angle accordingly with the number of phasing orbits phase_angle = phase_angle / phasing_orbits -# Calculate the eccentric anomaly +# Calculate the eccentric anomaly from the phase angle: + eccentric_anomaly = 2 * math.atan( math.sqrt((1 - eccentricity) / (1 + eccentricity)) * math.tan(phase_angle / 2) ) -# Calculate the time elapsed to cover phase angle in original orbit (Kepler's equation) +# Use Kepler's equation to calculate the time needed to cover the phase angle in the original orbit: + time_to_cover_phase_angle = (mean_orbital_period / (2 * math.pi)) * ( eccentric_anomaly - eccentricity * math.sin(eccentric_anomaly) ) -# Calculate the period of the phasing orbit +# Calculate the period of the phasing orbit: + phasing_orbit_period = mean_orbital_period - time_to_cover_phase_angle -# Calculate the SMA of the phasing orbit +# Then, calculate the SMA of the phasing orbit using the Earth's gravitational parameter, which is approximately 398600.44 $km^3s^{-2} $: + earth_gravitational_parameter = 398600.44 transfer_semi_major_axis = math.pow( (phasing_orbit_period * math.sqrt(earth_gravitational_parameter)) / (2 * math.pi), 2 / 3, ) -# Calculate eccentricity and radii +# Next, calculate the eccentricity and the radii of the initial orbit. Start by using the chaser satellite's Astrogator data provider to retrieve the Keplerian elements of the orbit. Retrieve this data during the entire scenario at 60 second intervals and convert the data to a Pandas dataframe: + keplerian_elements_data_provider = chaser_satellite.data_providers[ "Astrogator Values" ].group["Keplerian Elems"] @@ -298,8 +368,10 @@ scenario.start_time, scenario.start_time, 60 ) keplerian_dataframe = keplerian_result.data_sets.to_pandas_dataframe() -initial_radius = float(keplerian_dataframe.at[0, "radius_of_periapsis"]) +# Retrieve the initial radius of periapsis from the dataframe then use it to determine the periapsis and apoapsis radii: + +initial_radius = float(keplerian_dataframe.at[0, "radius_of_periapsis"]) if transfer_semi_major_axis > semimajor_axis: periapsis_radius = initial_radius apoapsis_radius = 2 * transfer_semi_major_axis - periapsis_radius @@ -307,80 +379,113 @@ apoapsis_radius = initial_radius periapsis_radius = 2 * transfer_semi_major_axis - apoapsis_radius -# Estimate the needed Delta V for phasing +# Then, calculate the velocity of the satellite at the initial periapsis: + velocity_periapsis_initial = math.sqrt( earth_gravitational_parameter * ((2 / initial_radius) - (1 / mean_semi_major_axis)) ) + +# Calculate the estimated desired velocity of the satellite at periapsis after phasing: + velocity_periapsis_final = math.sqrt( earth_gravitational_parameter * ((2 / initial_radius) - (1 / transfer_semi_major_axis)) ) + +# Finally, estimate the velocity change (Delta V) needed for phasing: + delta_v_estimate = velocity_periapsis_final - velocity_periapsis_initial -# - # ## Set up the phasing maneuver sequence # -# Add a target sequence to perform the first delta-V maneuver and propagate through the phasing orbits: +# Add a target sequence to perform the first delta-V maneuver and propagate through the phasing orbits. Add the first target sequence segment and configure it to use the active profiles: # + -from ansys.stk.core.stkobjects.astrogator import ( - AttitudeControl, - ControlManeuver, - ManeuverType, - ProfileMode, - TargetSequenceAction, -) +from ansys.stk.core.stkobjects.astrogator import TargetSequenceAction + -# Add the first Target Sequence segment phasing_start_sequence = chaser_satellite.propagator.main_sequence.insert( SegmentType.TARGET_SEQUENCE, "Start Phasing", "-" ) phasing_start_sequence.action = TargetSequenceAction.RUN_ACTIVE_PROFILES +# - + +# Next, add an impulsive maneuver segment, through which Astrogator calculates the new state of the satellite by adding a Delta-V vector to the final state velocity of the previous segment. Configure the maneuver to use a Cartesian X control parameter: + +# + +from ansys.stk.core.stkobjects.astrogator import ManeuverType, ControlManeuver + -# Add a Maneuver segment first_delta_v_maneuver = phasing_start_sequence.segments.insert( SegmentType.MANEUVER, "DV1", "-" ) first_delta_v_maneuver.set_maneuver_type(ManeuverType.IMPULSIVE) first_delta_v_maneuver.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) +# - + +# Configure the maneuver attitude to use a thrust vector with an x component 1000 times the estimated delta V: + +# + +from ansys.stk.core.stkobjects.astrogator import AttitudeControl + -# Maneuver attitude definition first_delta_v_maneuver.maneuver.set_attitude_control_type(AttitudeControl.THRUST_VECTOR) thrust_vector = first_delta_v_maneuver.maneuver.attitude_control thrust_vector.thrust_axes_name = "Satellite/Chaser VNC(Earth)" thrust_vector.x = delta_v_estimate * 1000 +# - + +# Finally, add a propagate segment using an Earth Point Mass propagator: -# Add a Propagate segment phasing_orbit_propagate = phasing_start_sequence.segments.insert( SegmentType.PROPAGATE, "Phasing Orbit", "-" ) -phasing_orbit_propagate.properties.color = Colors.Orange phasing_orbit_propagate.propagator_name = "Earth point mass" + +# Configure the segment to display in orange: + +phasing_orbit_propagate.properties.color = Colors.Orange + +# Configure the stopping conditions of the segment to use either apoapsis or periapsis based on the difference in SMAs between the orbits: + if transfer_semi_major_axis > semimajor_axis: phasing_orbit_propagate.stopping_conditions.add("Periapsis") else: phasing_orbit_propagate.stopping_conditions.add("Apoapsis") + +# Remove the default stopping condition of the segment: + phasing_orbit_propagate.stopping_conditions.remove("Duration") + +# Then, designate that the stopping condition must be satisfied for each phasing orbit before the overall condition is met: + phasing_orbit_propagate.stopping_conditions.item( 0 ).properties.repeat_count = phasing_orbits + +# Finally, designate that the stopping condition should be calculated for the angles between the satellite and the chaser and the satellite and target: + phasing_orbit_propagate.results.add("Vector/Angle Between Vectors") phasing_orbit_propagate.results[0].vector1_name = "Satellite/Chaser Position" phasing_orbit_propagate.results[0].vector2_name = "Satellite/Target Position" -# - # ## Configure the differential corrector for phasing # # Set up the differential corrector to adjust the first delta-V to achieve zero phase angle at the end of the phasing orbits: -# Customize the Differential Corrector +# + +from ansys.stk.core.stkobjects.astrogator import ProfileMode + + phasing_differential_corrector = phasing_start_sequence.profiles[ "Differential Corrector" ] phasing_differential_corrector.mode = ProfileMode.ITERATE phasing_differential_corrector.max_iterations = 50 +# - + +# Then, set the control parameter: -# Set Control Parameters and Results phasing_x_control_parameter = ( phasing_differential_corrector.control_parameters.get_control_by_paths( "DV1", "ImpulsiveMnvr.Cartesian.X" @@ -389,6 +494,8 @@ phasing_x_control_parameter.enable = True phasing_x_control_parameter.max_step = 0.001 +# Finally, configure the result of the corrector: + phasing_angle_result = phasing_differential_corrector.results.get_result_by_paths( "Phasing Orbit", "Angle_Between_Vectors" ) @@ -398,43 +505,51 @@ # ## Set up the circularization maneuver # -# Add a second target sequence to perform the circularization delta-V to return to the original circular orbit: +# Add a second target sequence to perform the circularization delta-V to return to the original circular orbit. + +# First, add a second target sequence segment to the chaser satellite's main sequence and configure it to use the active profiles: -# Add the second Target Sequence segment circularization_sequence = chaser_satellite.propagator.main_sequence.insert( SegmentType.TARGET_SEQUENCE, "Circularization", "-" ) circularization_sequence.action = TargetSequenceAction.RUN_ACTIVE_PROFILES -# Add a Maneuver segment +# Next, add a maneuver segment and configure it to use an impulsive control maneuver: + second_delta_v_maneuver = circularization_sequence.segments.insert( SegmentType.MANEUVER, "DV2", "-" ) second_delta_v_maneuver.set_maneuver_type(ManeuverType.IMPULSIVE) second_delta_v_maneuver.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) -# Add a Propagate segment +# Finally, add a propagate segment to move the satellite to the final orbit. Configure it to use an Earth Point Mass propagator and stop after a duration of one day: + final_orbit_propagate = circularization_sequence.segments.insert( SegmentType.PROPAGATE, "Final Orbit", "-" ) -final_orbit_propagate.properties.color = Colors.Yellow final_orbit_propagate.propagator_name = "Earth point mass" final_orbit_propagate.stopping_conditions["Duration"].properties.trip = 86400 + +# Record the eccentricity of the orbit during the segment as a result: + final_orbit_propagate.results.add("Keplerian Elems/Eccentricity") +# Set the segment to display in yellow: + +final_orbit_propagate.properties.color = Colors.Yellow + # ## Configure the differential corrector for circularization # # Set up the differential corrector to achieve a circular orbit with the desired eccentricity: -# + -# Customize the Differential Corrector circularization_differential_corrector = circularization_sequence.profiles[ "Differential Corrector" ] circularization_differential_corrector.mode = ProfileMode.ITERATE circularization_differential_corrector.max_iterations = 50 -# Set Control Parameters and Results +# Configure the corrector's control parameter: + circularization_x_control_parameter = ( circularization_differential_corrector.control_parameters.get_control_by_paths( "DV2", "ImpulsiveMnvr.Cartesian.X" @@ -443,6 +558,8 @@ circularization_x_control_parameter.enable = True circularization_x_control_parameter.max_step = 0.01 +# Set the result of the corrector to record the eccentricity of the orbit during the correction: + circularization_eccentricity_result = ( circularization_differential_corrector.results.get_result_by_paths( "Final Orbit", "Eccentricity" @@ -451,7 +568,6 @@ circularization_eccentricity_result.enable = True circularization_eccentricity_result.desired_value = 0.0001 circularization_eccentricity_result.tolerance = 0.00001 -# - # ## Run the main control sequence # @@ -463,15 +579,21 @@ # ## Retrieve the results # -# Once the analysis has been performed, retrieve the delta-V values and phasing orbit parameters: +# Once the analysis has been performed, retrieve the delta-V values and phasing orbit parameters using the chaser satellite's Maneuver Summary data provider: -# + -# Get the maneuver data providers maneuver_data_provider = chaser_satellite.data_providers["Maneuver Summary"] + +# Retrieve the data for this provider over the entire scenario and convert the data to a Pandas dataframe: + maneuver_result = maneuver_data_provider.execute(scenario.start_time, scenario.stop_time) maneuver_dataframe = maneuver_result.data_sets.to_pandas_dataframe() + +# Retrieve the actual delta v: + delta_v_actual = maneuver_dataframe.at[0, "delta v"] +# Display the transfer orbit data: + print("") print("################### Transfer orbit data ###################") print(f"N phasing orbits = {phasing_orbits}") @@ -482,7 +604,6 @@ print(f"Delta V = {delta_v_actual} m/sec") print(f"Transfer time = {phasing_orbit_period * phasing_orbits / 86400:.4f} days") print("") -# - # Finally, show the complete orbit trajectory: @@ -490,32 +611,34 @@ # ## Plot the phase angle history # -# Retrieve and plot the phase angle between the chaser and target satellites over time: +# Retrieve and plot the phase angle between the chaser and target satellites over time. To do so, use the scalar calculation AWB component representing the value of the phase angle: # + import matplotlib.pyplot as plt import numpy as np +import pandas as pd + -# Get the phase angle data over time -angle_data_provider = chaser_satellite.data_providers["Vector"].group["Angle Between Vectors"] -angle_result = angle_data_provider.execute(scenario.start_time, scenario.stop_time, 3600) +phase_angle_data = scenario.data_providers.item("Scalar Calculations").group.item("Earth PhaseAngle") +angle_result = phase_angle_data.execute(scenario.start_time, scenario.stop_time, 3600) angle_dataframe = angle_result.data_sets.to_pandas_dataframe() +# Convert columns to correct types +angle_dataframe["time"] = pd.to_datetime(angle_dataframe["time"]) +angle_dataframe["scalar"] = angle_dataframe["scalar"].apply(pd.to_numeric) -# Convert time to days -time_seconds = angle_dataframe["Time"].values -time_days = time_seconds / 86400.0 +# Get the time values +time = angle_dataframe["time"] # Get the angle values -phase_angles = angle_dataframe["Angle"].values +phase_angles = angle_dataframe["scalar"] # Create the plot plt.figure(figsize=(10, 6)) -plt.plot(time_days, phase_angles, 'b-', linewidth=2) -plt.xlabel('Time (days)', fontsize=12) +plt.plot(time, phase_angles, '-', linewidth=2, color="dodgerblue") +plt.xlabel('Time', fontsize=12) plt.ylabel('Phase Angle (degrees)', fontsize=12) -plt.title('Phase Angle History During Phasing Maneuver', fontsize=14, fontweight='bold') +plt.title('Phase Angle History During Phasing Maneuver', fontsize=14) plt.grid(True, alpha=0.3) -plt.xlim([0, propagation_time / 86400]) plt.tight_layout() plt.show() # - @@ -525,6 +648,9 @@ # Retrieve and plot the altitude of both satellites to visualize the phasing orbit: # + +from datetime import datetime + + # Get altitude data for chaser satellite chaser_altitude_provider = chaser_satellite.data_providers["Cartesian Position"].group["Fixed"] chaser_altitude_result = chaser_altitude_provider.execute(scenario.start_time, scenario.stop_time, 3600) @@ -532,9 +658,9 @@ # Calculate magnitude of position vector and subtract Earth radius earth_radius = 6378.137 # km -chaser_x = chaser_altitude_df["x"].values -chaser_y = chaser_altitude_df["y"].values -chaser_z = chaser_altitude_df["z"].values +chaser_x = chaser_altitude_df["x"].values.astype(float) +chaser_y = chaser_altitude_df["y"].values.astype(float) +chaser_z = chaser_altitude_df["z"].values.astype(float) chaser_altitude = np.sqrt(chaser_x**2 + chaser_y**2 + chaser_z**2) - earth_radius # Get altitude data for target satellite @@ -542,25 +668,27 @@ target_altitude_result = target_altitude_provider.execute(scenario.start_time, scenario.stop_time, 3600) target_altitude_df = target_altitude_result.data_sets.to_pandas_dataframe() -target_x = target_altitude_df["x"].values -target_y = target_altitude_df["y"].values -target_z = target_altitude_df["z"].values +target_x = target_altitude_df["x"].values.astype(float) +target_y = target_altitude_df["y"].values.astype(float) +target_z = target_altitude_df["z"].values.astype(float) target_altitude = np.sqrt(target_x**2 + target_y**2 + target_z**2) - earth_radius # Convert time to days -chaser_time_days = chaser_altitude_df["Time"].values / 86400.0 -target_time_days = target_altitude_df["Time"].values / 86400.0 +chaser_altitude_df["time"] = pd.to_datetime(chaser_altitude_df["time"]) +target_altitude_df["time"] = pd.to_datetime(target_altitude_df["time"]) +start = np.datetime64(datetime.strptime(start_time, "%d %b %Y %H:%M:%S.%f")) +chaser_time_days = (chaser_altitude_df["time"].values - start) / 86400.0 +target_time_days = (target_altitude_df["time"].values - start) / 86400.0 # Create the plot plt.figure(figsize=(10, 6)) -plt.plot(chaser_time_days, chaser_altitude, 'b-', linewidth=2, label='Chaser') -plt.plot(target_time_days, target_altitude, 'r--', linewidth=2, label='Target') +plt.plot(chaser_time_days, chaser_altitude, '-', linewidth=2, label='Chaser', color="dodgerblue") +plt.plot(target_time_days, target_altitude, '--', linewidth=2, label='Target', color="coral") plt.xlabel('Time (days)', fontsize=12) plt.ylabel('Altitude (km)', fontsize=12) -plt.title('Altitude History During Phasing Maneuver', fontsize=14, fontweight='bold') +plt.title('Altitude History During Phasing Maneuver', fontsize=14) plt.legend(fontsize=11) plt.grid(True, alpha=0.3) -plt.xlim([0, propagation_time / 86400]) plt.tight_layout() plt.show() # - From df49bc2ae47447692fa228591324e99cfa5ead49 Mon Sep 17 00:00:00 2001 From: Dupont Date: Wed, 8 Jul 2026 16:56:09 -0400 Subject: [PATCH 09/18] doc: working on example code style --- examples/phasing-geo-orbits.py | 71 +++++++++++++++++++++++----------- 1 file changed, 49 insertions(+), 22 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index b0be215bcb..72db503fb9 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -49,7 +49,7 @@ eccentricity = 0.0001 inclination = 0.0 raan = 0.0 -arg_of_periapsis = 0.0 +arg_of_periapsis = 0.0 true_anomaly_target = 20.0 true_anomaly_chaser = 70.0 @@ -129,8 +129,12 @@ # Next, retrieve all of the vectors associated with the Earth, the target satellite, and the chaser satellite: -target_vectors = awb_root.get_provider(target_satellite.class_name+"/"+target_satellite.instance_name).vectors -chaser_vectors = awb_root.get_provider(chaser_satellite.class_name+"/"+chaser_satellite.instance_name).vectors +target_vectors = awb_root.get_provider( ++ target_satellite.class_name + "/" + target_satellite.instance_name ++).vectors ++chaser_vectors = awb_root.get_provider( ++ chaser_satellite.class_name + "/" + chaser_satellite.instance_name ++).vectors # Retrieve the target satellite's position vector, which is automatically created on STK objects when they are created. Then, designate this vector as the vector that the dihedral angle sweeps to: @@ -153,7 +157,9 @@ # Finally, create a scalar calculation using the dihedral angle as input in order to record the angle values over time: -earth_calc_scalar_factory = awb_root.get_provider("CentralBody/Earth").calculation_scalars.factory +earth_calc_scalar_factory = awb_root.get_provider( + "CentralBody/Earth" +).calculation_scalars.factory phase_angle_calc = earth_calc_scalar_factory.create_angle("PhaseAngle", "PhaseAngle") phase_angle_calc.input_angle = phase_angle_awb @@ -212,9 +218,9 @@ ) target_propagate_segment.propagator_name = "Earth point mass" -target_propagate_segment.stopping_conditions["Duration"].properties.trip = ( - propagation_time -) +target_propagate_segment.stopping_conditions[ + "Duration" +].properties.trip = propagation_time # - # Use the color red to identify the target satellite: @@ -413,7 +419,7 @@ # Next, add an impulsive maneuver segment, through which Astrogator calculates the new state of the satellite by adding a Delta-V vector to the final state velocity of the previous segment. Configure the maneuver to use a Cartesian X control parameter: # + -from ansys.stk.core.stkobjects.astrogator import ManeuverType, ControlManeuver +from ansys.stk.core.stkobjects.astrogator import ControlManeuver, ManeuverType first_delta_v_maneuver = phasing_start_sequence.segments.insert( @@ -585,7 +591,9 @@ # Retrieve the data for this provider over the entire scenario and convert the data to a Pandas dataframe: -maneuver_result = maneuver_data_provider.execute(scenario.start_time, scenario.stop_time) +maneuver_result = maneuver_data_provider.execute( + scenario.start_time, scenario.stop_time +) maneuver_dataframe = maneuver_result.data_sets.to_pandas_dataframe() # Retrieve the actual delta v: @@ -619,7 +627,9 @@ import pandas as pd -phase_angle_data = scenario.data_providers.item("Scalar Calculations").group.item("Earth PhaseAngle") +phase_angle_data = scenario.data_providers.item("Scalar Calculations").group.item( + "Earth PhaseAngle" +) angle_result = phase_angle_data.execute(scenario.start_time, scenario.stop_time, 3600) angle_dataframe = angle_result.data_sets.to_pandas_dataframe() # Convert columns to correct types @@ -635,9 +645,9 @@ # Create the plot plt.figure(figsize=(10, 6)) plt.plot(time, phase_angles, '-', linewidth=2, color="dodgerblue") -plt.xlabel('Time', fontsize=12) -plt.ylabel('Phase Angle (degrees)', fontsize=12) -plt.title('Phase Angle History During Phasing Maneuver', fontsize=14) +plt.xlabel("Time", fontsize=12) +plt.ylabel("Phase Angle (degrees)", fontsize=12) +plt.title("Phase Angle History During Phasing Maneuver", fontsize=14) plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() @@ -652,8 +662,12 @@ # Get altitude data for chaser satellite -chaser_altitude_provider = chaser_satellite.data_providers["Cartesian Position"].group["Fixed"] -chaser_altitude_result = chaser_altitude_provider.execute(scenario.start_time, scenario.stop_time, 3600) +chaser_altitude_provider = chaser_satellite.data_providers["Cartesian Position"].group[ + "Fixed" +] +chaser_altitude_result = chaser_altitude_provider.execute( + scenario.start_time, scenario.stop_time, 3600 +) chaser_altitude_df = chaser_altitude_result.data_sets.to_pandas_dataframe() # Calculate magnitude of position vector and subtract Earth radius @@ -664,8 +678,12 @@ chaser_altitude = np.sqrt(chaser_x**2 + chaser_y**2 + chaser_z**2) - earth_radius # Get altitude data for target satellite -target_altitude_provider = target_satellite.data_providers["Cartesian Position"].group["Fixed"] -target_altitude_result = target_altitude_provider.execute(scenario.start_time, scenario.stop_time, 3600) +target_altitude_provider = target_satellite.data_providers["Cartesian Position"].group[ + "Fixed" +] +target_altitude_result = target_altitude_provider.execute( + scenario.start_time, scenario.stop_time, 3600 +) target_altitude_df = target_altitude_result.data_sets.to_pandas_dataframe() target_x = target_altitude_df["x"].values.astype(float) @@ -682,11 +700,20 @@ # Create the plot plt.figure(figsize=(10, 6)) -plt.plot(chaser_time_days, chaser_altitude, '-', linewidth=2, label='Chaser', color="dodgerblue") -plt.plot(target_time_days, target_altitude, '--', linewidth=2, label='Target', color="coral") -plt.xlabel('Time (days)', fontsize=12) -plt.ylabel('Altitude (km)', fontsize=12) -plt.title('Altitude History During Phasing Maneuver', fontsize=14) +plt.plot( + chaser_time_days, + chaser_altitude, + "-", + linewidth=2, + label="Chaser", + color="dodgerblue", +) +plt.plot( + target_time_days, target_altitude, "--", linewidth=2, label="Target", color="coral" +) +plt.xlabel("Time (days)", fontsize=12) +plt.ylabel("Altitude (km)", fontsize=12) +plt.title("Altitude History During Phasing Maneuver", fontsize=14) plt.legend(fontsize=11) plt.grid(True, alpha=0.3) plt.tight_layout() From 25da7499126d5455e51536d8ddad5062fe584e9e Mon Sep 17 00:00:00 2001 From: Dupont Date: Wed, 8 Jul 2026 17:03:16 -0400 Subject: [PATCH 10/18] doc: working on example code style --- examples/phasing-geo-orbits.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 72db503fb9..dff40a8ae7 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -130,11 +130,11 @@ # Next, retrieve all of the vectors associated with the Earth, the target satellite, and the chaser satellite: target_vectors = awb_root.get_provider( -+ target_satellite.class_name + "/" + target_satellite.instance_name -+).vectors -+chaser_vectors = awb_root.get_provider( -+ chaser_satellite.class_name + "/" + chaser_satellite.instance_name -+).vectors + target_satellite.class_name + "/" + target_satellite.instance_name +).vectors +chaser_vectors = awb_root.get_provider( + chaser_satellite.class_name + "/" + chaser_satellite.instance_name +).vectors # Retrieve the target satellite's position vector, which is automatically created on STK objects when they are created. Then, designate this vector as the vector that the dihedral angle sweeps to: From 03ab3c59c31cff100e08f357377394a865c2e1f2 Mon Sep 17 00:00:00 2001 From: Dupont Date: Thu, 9 Jul 2026 10:28:54 -0400 Subject: [PATCH 11/18] doc: finish example --- examples/phasing-geo-orbits.py | 54 +++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 21 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index dff40a8ae7..039546cf76 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -40,7 +40,7 @@ # - RAAN: 0.0 degrees # - Argument of periapsis: 0.0 degrees # -# Compute the required $\Delta v$ for the phasing maneuvers and determine the phasing orbit parameters. +# Compute the required $\Delta V$ for the phasing maneuvers and determine the phasing orbit parameters. # First, define the initial orbital parameters for use during calculations: @@ -98,7 +98,7 @@ plotter.show() # - -# ## Create the satellites +# ## Add the satellites # Create the target satellite which will remain at a fixed position in geostationary orbit. This satellite serves as the reference point for the phasing maneuver: @@ -114,7 +114,7 @@ awb_root = root.analysis_workbench_components_root -# The phase angle component will be associated with the Earth central body, so retrieve the factory that allows creating AWB components on this body: +# The phase angle component is associated with the Earth central body, so retrieve the factory that allows creating AWB components on this body: earth_angle_factory = awb_root.get_provider("CentralBody/Earth").angles.factory @@ -127,7 +127,7 @@ phase_angle_awb = earth_angle_factory.create("Phase", "Phase", AngleType.DIHEDRAL_ANGLE) # - -# Next, retrieve all of the vectors associated with the Earth, the target satellite, and the chaser satellite: +# Next, retrieve all of the vectors associated with the target satellite and the chaser satellite: target_vectors = awb_root.get_provider( target_satellite.class_name + "/" + target_satellite.instance_name @@ -136,7 +136,7 @@ chaser_satellite.class_name + "/" + chaser_satellite.instance_name ).vectors -# Retrieve the target satellite's position vector, which is automatically created on STK objects when they are created. Then, designate this vector as the vector that the dihedral angle sweeps to: +# Retrieve the target satellite's position vector. Position vectors are automatically created on STK objects when they are added to a scenario. Then, designate this vector as the vector that the dihedral angle sweeps to: target_position_vector = target_vectors.item("Position") phase_angle_awb.to_vector.set_vector(target_position_vector) @@ -197,7 +197,7 @@ target_initial_state.initial_state.epoch = scenario.start_time # - -# Delcare the Keplerian elements for the initial state using the previously defined parameters: +# Declare the Keplerian elements for the initial state using the previously defined parameters: # + target_keplerian_elements = target_initial_state.element @@ -256,7 +256,7 @@ # Configure graphics settings to allow 3D trajectory visualization: -chaser_satellite.propagator.options.draw_trajectory_in_3d = False +chaser_satellite.propagator.options.draw_trajectory_in_3d = True # ## Set up the initial state of the chaser satellite # @@ -398,13 +398,13 @@ * ((2 / initial_radius) - (1 / transfer_semi_major_axis)) ) -# Finally, estimate the velocity change (Delta V) needed for phasing: +# Finally, estimate the velocity change (Delta-V) needed for phasing: delta_v_estimate = velocity_periapsis_final - velocity_periapsis_initial # ## Set up the phasing maneuver sequence # -# Add a target sequence to perform the first delta-V maneuver and propagate through the phasing orbits. Add the first target sequence segment and configure it to use the active profiles: +# Add a target sequence to perform the first Delta-V maneuver and propagate through the phasing orbits. Add the first target sequence segment and configure it to use the active profiles: # + from ansys.stk.core.stkobjects.astrogator import TargetSequenceAction @@ -429,7 +429,7 @@ first_delta_v_maneuver.enable_control_parameter(ControlManeuver.IMPULSIVE_CARTESIAN_X) # - -# Configure the maneuver attitude to use a thrust vector with an x component 1000 times the estimated delta V: +# Configure the maneuver attitude to use a thrust vector with an x component 1000 times the estimated Delta-V: # + from ansys.stk.core.stkobjects.astrogator import AttitudeControl @@ -477,7 +477,7 @@ # ## Configure the differential corrector for phasing # -# Set up the differential corrector to adjust the first delta-V to achieve zero phase angle at the end of the phasing orbits: +# Set up the differential corrector to adjust the first Delta-V to achieve zero phase angle at the end of the phasing orbits: # + from ansys.stk.core.stkobjects.astrogator import ProfileMode @@ -511,7 +511,7 @@ # ## Set up the circularization maneuver # -# Add a second target sequence to perform the circularization delta-V to return to the original circular orbit. +# Add a second target sequence to perform the circularization Delta-V to return to the original circular orbit. # First, add a second target sequence segment to the chaser satellite's main sequence and configure it to use the active profiles: @@ -585,7 +585,7 @@ # ## Retrieve the results # -# Once the analysis has been performed, retrieve the delta-V values and phasing orbit parameters using the chaser satellite's Maneuver Summary data provider: +# Once the analysis has been performed, retrieve the Delta-V values and phasing orbit parameters using the chaser satellite's Maneuver Summary data provider: maneuver_data_provider = chaser_satellite.data_providers["Maneuver Summary"] @@ -596,7 +596,7 @@ ) maneuver_dataframe = maneuver_result.data_sets.to_pandas_dataframe() -# Retrieve the actual delta v: +# Retrieve the actual Delta-V: delta_v_actual = maneuver_dataframe.at[0, "delta v"] @@ -609,7 +609,7 @@ print(f"SMA = {transfer_semi_major_axis:.2f} km") print(f"Perigee radius = {periapsis_radius:.2f} km") print(f"Apogee radius = {apoapsis_radius:.2f} km") -print(f"Delta V = {delta_v_actual} m/sec") +print(f"Delta-V = {delta_v_actual} m/sec") print(f"Transfer time = {phasing_orbit_period * phasing_orbits / 86400:.4f} days") print("") @@ -625,6 +625,7 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd +from datetime import datetime phase_angle_data = scenario.data_providers.item("Scalar Calculations").group.item( @@ -639,17 +640,26 @@ # Get the time values time = angle_dataframe["time"] +# Convert time to days +start = np.datetime64(datetime.strptime(start_time, "%d %b %Y %H:%M:%S.%f")) +time_days = (time.values - start) / 86400.0 + # Get the angle values phase_angles = angle_dataframe["scalar"] # Create the plot plt.figure(figsize=(10, 6)) -plt.plot(time, phase_angles, '-', linewidth=2, color="dodgerblue") -plt.xlabel("Time", fontsize=12) +plt.plot(time_days, phase_angles, "-", linewidth=2, color="dodgerblue") +plt.xlabel("Time (days from start of maneuver)", fontsize=12) plt.ylabel("Phase Angle (degrees)", fontsize=12) plt.title("Phase Angle History During Phasing Maneuver", fontsize=14) plt.grid(True, alpha=0.3) plt.tight_layout() + +# Configure style +plt.gca().set_facecolor("whitesmoke") +plt.gca().grid(visible=True, which="both", linestyle="--") + plt.show() # - @@ -658,9 +668,6 @@ # Retrieve and plot the altitude of both satellites to visualize the phasing orbit: # + -from datetime import datetime - - # Get altitude data for chaser satellite chaser_altitude_provider = chaser_satellite.data_providers["Cartesian Position"].group[ "Fixed" @@ -711,12 +718,17 @@ plt.plot( target_time_days, target_altitude, "--", linewidth=2, label="Target", color="coral" ) -plt.xlabel("Time (days)", fontsize=12) +plt.xlabel("Time (days from start of maneuver)", fontsize=12) plt.ylabel("Altitude (km)", fontsize=12) plt.title("Altitude History During Phasing Maneuver", fontsize=14) plt.legend(fontsize=11) plt.grid(True, alpha=0.3) plt.tight_layout() + +# Configure style +plt.gca().set_facecolor("whitesmoke") +plt.gca().grid(visible=True, which="both", linestyle="--") + plt.show() # - From 301b46db524c72f8f1f2bd77d3ecc7b072b6b47e Mon Sep 17 00:00:00 2001 From: Dupont Date: Thu, 9 Jul 2026 10:35:23 -0400 Subject: [PATCH 12/18] code style fix --- examples/phasing-geo-orbits.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 039546cf76..9a371e1a87 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -622,10 +622,11 @@ # Retrieve and plot the phase angle between the chaser and target satellites over time. To do so, use the scalar calculation AWB component representing the value of the phase angle: # + +from datetime import datetime + import matplotlib.pyplot as plt import numpy as np import pandas as pd -from datetime import datetime phase_angle_data = scenario.data_providers.item("Scalar Calculations").group.item( From 0080a5d5b0234a85e19f8be7015cf3ba9347d09e Mon Sep 17 00:00:00 2001 From: Dupont Date: Thu, 9 Jul 2026 15:23:55 -0400 Subject: [PATCH 13/18] fix doc build --- examples/phasing-geo-orbits.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 9a371e1a87..bc6c97c13f 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -74,14 +74,9 @@ # # Start by creating a new scenario in STK: -# + -from ansys.stk.core.stkobjects import STKObjectType - - root = stk.new_object_root() root.new_scenario("Phasing_orbit") scenario = root.current_scenario -# - # Configure the scenario time to cover a 16 day period starting July 1, 2026: @@ -91,7 +86,7 @@ # It is now possible to show a 3D graphics window by running: # + -from ansys.stk.core.stkengine.experimental.jupyterwidgets import GlobeWidget +from ansys.stk.core.experimental.jupyterwidgets import GlobeWidget plotter = GlobeWidget(root, 640, 480) @@ -102,7 +97,12 @@ # Create the target satellite which will remain at a fixed position in geostationary orbit. This satellite serves as the reference point for the phasing maneuver: +# + +from ansys.stk.core.stkobjects import STKObjectType + + target_satellite = scenario.children.new(STKObjectType.SATELLITE, "Target") +# - # Create the chaser satellite, which will need to approach the target satellite: From 0c45fa11843e0f72aa78ea904d2165901fe166a1 Mon Sep 17 00:00:00 2001 From: Dupont Date: Fri, 10 Jul 2026 10:21:51 -0400 Subject: [PATCH 14/18] fix mcs to mcs --- examples/phasing-geo-orbits.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index bc6c97c13f..589c087d13 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -238,7 +238,7 @@ from ansys.stk.core.stkobjects.astrogator import RunCode -run_code = target_satellite.propagator.run_mcs2() +run_code = target_satellite.propagator.run_mcs() if run_code != RunCode.MARCHING: raise ValueError("Could not propagate target satellite orbit.") # - @@ -292,7 +292,7 @@ # Run the mission control sequence of the chaser and verify that the propagation was successful: -run_code = chaser_satellite.propagator.run_mcs2() +run_code = chaser_satellite.propagator.run_mcs() if run_code != RunCode.MARCHING: raise ValueError("Could not propagate chaser satellite orbit.") From 5114af162f74f604764dbfd87362bc09eaf485e3 Mon Sep 17 00:00:00 2001 From: Dupont Date: Sat, 11 Jul 2026 23:34:56 -0400 Subject: [PATCH 15/18] fix toctree --- doc/source/examples.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/source/examples.rst b/doc/source/examples.rst index a3df2336b4..9b0fce15cd 100644 --- a/doc/source/examples.rst +++ b/doc/source/examples.rst @@ -51,6 +51,7 @@ Examples examples/hohmann-transfer examples/bielliptic-transfer examples/lambert-transfer + examples/phasing-geo-orbits .. nbgallery:: :caption: Mission analysis From dd7e9ffcc9c47665cc9e95b2652249d49f9ca532 Mon Sep 17 00:00:00 2001 From: Dupont Date: Wed, 15 Jul 2026 09:49:47 -0400 Subject: [PATCH 16/18] code review comments --- examples/phasing-geo-orbits.py | 19 +++++++++---------- file | 0 2 files changed, 9 insertions(+), 10 deletions(-) delete mode 100644 file diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 589c087d13..7677a0c7f6 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -18,9 +18,9 @@ # # $$ \theta_{gap}=\frac{\theta}{n_{orbits}}$$ # -# , where $\theta$ is the phase angle, and $\theta_{gap}$ is the angular dispacement to recover after each phasing orbit. +# where $\theta$ is the phase angle, and $\theta_{gap}$ is the angular displacement to recover after each phasing orbit. # -# In the code below the **Kepler's equation** is used to calculate the period and semimajor axis of the phasing orbit, given the number of revolution around the phasing orbit itself. As first, the **eccentric anomaly E** is calculated from $\theta_{gap}$ : +# In the code below, **Kepler's equation** is used to calculate the period and semimajor axis of the phasing orbit, given the number of revolutions around the phasing orbit itself. As first, the **eccentric anomaly E** is calculated from $\theta_{gap}$ : # # $$\tan\left(\frac{E}{2}\right)=\sqrt{\frac{1 - e}{1 + e}}\tan\left(\frac{\theta_{gap}}{2}\right)$$ # @@ -60,7 +60,7 @@ # ## Launch a new STK instance # -# Start by launching a new STK instance. In this example, ``STKEngine`` is used with graphics (``no_graphics`` mode set to ``False``). This means that the graphic user interface (GUI) of the product is not launched but 2D and 3D visualization is still available through the STK Engine controls: +# Start by launching a new STK instance. In this example, ``STKEngine`` is used with graphics (``no_graphics`` mode set to ``False``). This means that the STK graphic user interface (GUI) of the product is not launched but 2D and 3D visualization is still available through the STK Engine controls: # + from ansys.stk.core.stkengine import STKEngine @@ -75,7 +75,7 @@ # Start by creating a new scenario in STK: root = stk.new_object_root() -root.new_scenario("Phasing_orbit") +root.new_scenario("PhasingOrbit") scenario = root.current_scenario # Configure the scenario time to cover a 16 day period starting July 1, 2026: @@ -165,7 +165,7 @@ # ## Set up the target satellite -# Then, declare the type of orbit propagator used for the satellite. Set the propagator to an Astrogator type: +# Declare the type of orbit propagator used for the chaser satellite. Set the propagator to an Astrogator type: # + from ansys.stk.core.stkobjects import PropagatorType @@ -249,7 +249,7 @@ # ## Set up the chaser satellite -# Then, declare the type of orbit propagator used for the satellite to an Astrogator type and ensure a clean main sequence: +# Declare the type of orbit propagator used for the target satellite to an Astrogator type and ensure a clean main sequence: chaser_satellite.set_propagator_type(PropagatorType.ASTROGATOR) chaser_satellite.propagator.main_sequence.remove_all() @@ -334,17 +334,16 @@ phase_angle = 2 * math.pi + phase_angle print(f"The initial phase angle is {math.degrees(phase_angle):.1f} deg.") -print("") # - # Then, calculate the phase angle per phasing orbit: -phase_angle = phase_angle / phasing_orbits +phase_angle_per_orbit = phase_angle / phasing_orbits # Calculate the eccentric anomaly from the phase angle: eccentric_anomaly = 2 * math.atan( - math.sqrt((1 - eccentricity) / (1 + eccentricity)) * math.tan(phase_angle / 2) + math.sqrt((1 - eccentricity) / (1 + eccentricity)) * math.tan(phase_angle_per_orbit / 2) ) # Use Kepler's equation to calculate the time needed to cover the phase angle in the original orbit: @@ -469,7 +468,7 @@ 0 ).properties.repeat_count = phasing_orbits -# Finally, designate that the stopping condition should be calculated for the angles between the satellite and the chaser and the satellite and target: +# Finally, designate that the stopping condition should be calculated for the angle between the chaser's position vector and the target's position vector: phasing_orbit_propagate.results.add("Vector/Angle Between Vectors") phasing_orbit_propagate.results[0].vector1_name = "Satellite/Chaser Position" diff --git a/file b/file deleted file mode 100644 index e69de29bb2..0000000000 From b4ac14acf44ddbc13896ffb18708c58d6de69bd8 Mon Sep 17 00:00:00 2001 From: Dupont Date: Wed, 15 Jul 2026 10:00:55 -0400 Subject: [PATCH 17/18] code style --- examples/phasing-geo-orbits.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 7677a0c7f6..4722dcb4fd 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -249,7 +249,7 @@ # ## Set up the chaser satellite -# Declare the type of orbit propagator used for the target satellite to an Astrogator type and ensure a clean main sequence: +# Declare the type of orbit propagator used for the target satellite as the Astrogator type and ensure a clean main sequence: chaser_satellite.set_propagator_type(PropagatorType.ASTROGATOR) chaser_satellite.propagator.main_sequence.remove_all() @@ -343,7 +343,8 @@ # Calculate the eccentric anomaly from the phase angle: eccentric_anomaly = 2 * math.atan( - math.sqrt((1 - eccentricity) / (1 + eccentricity)) * math.tan(phase_angle_per_orbit / 2) + math.sqrt((1 - eccentricity) / (1 + eccentricity)) + * math.tan(phase_angle_per_orbit / 2) ) # Use Kepler's equation to calculate the time needed to cover the phase angle in the original orbit: From 75cda617d2b9c1396ca89a4f6c01571cb6ad9396 Mon Sep 17 00:00:00 2001 From: Dupont Date: Wed, 15 Jul 2026 10:06:42 -0400 Subject: [PATCH 18/18] trying to fix unit --- examples/phasing-geo-orbits.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/phasing-geo-orbits.py b/examples/phasing-geo-orbits.py index 4722dcb4fd..2a6bfd164c 100644 --- a/examples/phasing-geo-orbits.py +++ b/examples/phasing-geo-orbits.py @@ -357,7 +357,7 @@ phasing_orbit_period = mean_orbital_period - time_to_cover_phase_angle -# Then, calculate the SMA of the phasing orbit using the Earth's gravitational parameter, which is approximately 398600.44 $km^3s^{-2} $: +# Then, calculate the SMA of the phasing orbit using the Earth's gravitational parameter, which is approximately 398600.44 $km^3 s^{-2} $: earth_gravitational_parameter = 398600.44 transfer_semi_major_axis = math.pow(