import yaml
from typing import Union, Dict
from datetime import datetime, timezone
from math import pi
import libs.orekit_setup
from orekit.pyhelpers import datetime_to_absolutedate
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.models.earth.atmosphere import NRLMSISE00
from org.orekit.models.earth.atmosphere.data import MarshallSolarActivityFutureEstimation
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.propagation import SpacecraftState
Osculating Propagator
from org.orekit.propagation.numerical import NumericalPropagator
Mean Propagator
from org.orekit.propagation.semianalytical.dsst import DSSTPropagator
from org.orekit.propagation.semianalytical.dsst.forces import DSSTAtmosphericDrag, DSSTSolarRadiationPressure,
DSSTThirdBody, DSSTZonal, DSSTTesseral
from org.orekit.propagation import PropagationType
from org.orekit.orbits import KeplerianOrbit, PositionAngle, OrbitType
from org.hipparchus.geometry.euclidean.threed import RotationOrder
from org.orekit.frames import FramesFactory, LOFType
from org.orekit.bodies import CelestialBodyFactory
from org.orekit.utils import Constants, IERSConventions
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.attitudes import LofOffset
Constants
sun = CelestialBodyFactory.getSun()
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(IERSConventions.IERS_2010, True))
moon = CelestialBodyFactory.getMoon()
R_E = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
w_E = Constants.WGS84_EARTH_ANGULAR_VELOCITY
mu = Constants.WGS84_EARTH_MU
velocity_pointing = LofOffset(FramesFactory.getEME2000(), LOFType.TNW, RotationOrder.YXY, pi/2, 0.0, 0.0)
spacecraft_data = {‘info’: {‘platform’: ‘Ciaociao’, ‘name’: ‘Ciaociao’, ‘COSPAR_ID’: ‘2027-001G’, ‘NORAD_ID’: ‘88888’},
‘disturbances’: {‘SRP’: {‘crossSection’: 1.5, ‘reflectivityCoefficient’: 1.3}, ‘drag’: {‘crossSection’: 1.5, ‘dragCoefficient’: 2.2}, ‘geometry’: {‘spacecraft_size_x’: 1.0509, ‘spacecraft_size_y’: 1.19269, ‘spacecraft_size_z’: 1.1}},
‘platform’: {‘mass_dry_empty’: 277.486, ‘cog_dry_empty’: [-0.00363, -0.00201, 0.45725], ‘I_dry_empty’: [64.8, 0.311, -0.5, 0.311, 68.9, 0.219, -0.502, 0.219, 22.8]},
‘deployables’: {‘1’: {‘mass’: 50.4, ‘cog’: [-0.00135, -0.00115, 0.05935], ‘cgeom’: [0.25425, -0.10375, 0.87176], ‘I’: [0.0113, 0.004, 0.0113], ‘disturbances’: {‘SRP’: {‘crossSection’: 0.03, ‘reflectivityCoefficient’: 1.8}, ‘drag’: {‘crossSection’: 0.03, ‘dragCoefficient’: 2.2}}}},
‘propellant’: {‘reference_date’: datetime(2022, 11, 17, 0, 0, tzinfo=timezone.utc), ‘mass_ox’: 28.8, ‘mass_fu’: 6.9, ‘cog_ox’: [0, 0, 0.34353], ‘cog_fu’: [0, 0, 0.35067], ‘I_ox’: [4.14, 0, 0, 0, 4.14, 0, 0, 0, 0.931], ‘I_fu’: [2.14, 0, 0, 0, 0.861, 0, 0, 0, 1.31]}
}
propagation_data = {‘computeMeanState’: {‘convergence_threshold’: ‘1E-12’, ‘max_iteration’: ‘1E13’},
‘integrator’: {‘minStep’: 1e-07, ‘maxStep’: 2700.0, ‘absTol’: 1e-11, ‘relTol’: 1e-11},
‘gravity’: {‘harmonics_degree’: 64, ‘harmonics_order’: 64}, ‘drag’: {‘solar_activity_level’: ‘AVERAGE’}, ‘elevationDetector’: {‘maxcheck’: 60.0, ‘threshold’: 0.001},
‘attitude’: {‘mode’: ‘NONE’, ‘details’: {‘mode’: ‘NONE’, ‘details’: {‘sunPointingAxis’: ‘+X’, ‘phasingAxis’: ‘+Z’}}},
‘propagator_type’: {‘cubesat_deployment_analysis’: ‘OSCULATING’, ‘orbit_transfer_analysis’: ‘OSCULATING’}}
def _add_dsst_force_models(dsst_propagator: DSSTPropagator, spacecraft_data, propagation_data) → DSSTPropagator:
“”"
Adds semi-analytical force models to DSST Propagator. Drag, SRP and Gravity models are configured according to
spacecraft data and propagation configurations
:param dsst_propagator: DSST propagator where forces will be added
:param spacecraft_data: dictionary or path to yml containing ‘drag’ (‘dragCoefficient’ and ‘crossSection’) and ‘SRP’
(‘reflectivityCoefficient’ and ‘crossSection’) data
:param propagation_data: dictionary or path to yml containing ‘gravity’ (‘harmonics_degree’ and ‘harmonics_order’) data
:return: DSST Propagator with Drag, SRP, Zonal Harmonics, Moon and Sun Third-Body force models
“”"
propagation_conf = propagation_data
spacecraft_conf = spacecraft_data
cd = spacecraft_conf[‘disturbances’][‘drag’][‘dragCoefficient’]
cr = spacecraft_conf[‘disturbances’][‘SRP’][‘reflectivityCoefficient’]
drag_area = spacecraft_conf[‘disturbances’][‘drag’][‘crossSection’]
srp_area = spacecraft_conf[‘disturbances’][‘SRP’][‘crossSection’]
harmonics_degree = propagation_conf[‘gravity’][‘harmonics_degree’]
harmonics_order = propagation_conf[‘gravity’][‘harmonics_order’]
if propagation_conf[‘drag’][‘solar_activity_level’] == ‘AVERAGE’:
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
elif propagation_conf[‘drag’][‘solar_activity_level’] == ‘STRONG’:
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.STRONG
elif propagation_conf[‘drag’][‘solar_activity_level’] == ‘WEAK’:
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.WEAK
else:
print(“Unrecognized setting for solar_activity_level in propagation_conf. Setting AVERAGE.”)
strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
# DRAG
supported_names = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
input_data = MarshallSolarActivityFutureEstimation(supported_names, strength_level)
dsst_atmosphere = NRLMSISE00(input_data, sun, earth)
dsst_drag = DSSTAtmosphericDrag(dsst_atmosphere, cd, drag_area, mu)
dsst_propagator.addForceModel(dsst_drag)
# SRP
dsst_srp = DSSTSolarRadiationPressure(cr, srp_area, sun, R_E, mu)
dsst_propagator.addForceModel(dsst_srp)
# Third Body
dsst_moon = DSSTThirdBody(moon, mu)
dsst_sun = DSSTThirdBody(sun, mu)
dsst_propagator.addForceModel(dsst_moon)
dsst_propagator.addForceModel(dsst_sun)
# Zonal and Tesseral Harmonics
provider = GravityFieldFactory.getUnnormalizedProvider(harmonics_degree, harmonics_order)
dsst_zonal = DSSTZonal(provider)
dsst_propagator.addForceModel(dsst_zonal)
dsst_tesseral = DSSTTesseral(FramesFactory.getITRF(IERSConventions.IERS_2010, True), w_E, provider)
dsst_propagator.addForceModel(dsst_tesseral)
return dsst_propagator
minStep = 0.0000001 # [s] smallest integration step allowed
maxStep = 2700.0 # [s] longest integration step allowed
absTol = 9.999999999999999e-12 # absolute tolerance on integration error
relTol = 9.999999999999999e-12 # relative tolerance on difference between solutions
integrator = DormandPrince853Integrator(minStep, maxStep, absTol, relTol)
dsst_propagator_temp = DSSTPropagator(integrator, PropagationType.MEAN)
propagator = _add_dsst_force_models(dsst_propagator_temp, spacecraft_data, propagation_data)
orbit = KeplerianOrbit(float(6958000), float(0.001), float(0.090757), float(0.01), float(0.01), float(0.01), PositionAngle.TRUE,
FramesFactory.getEME2000(), datetime_to_absolutedate(datetime.strptime(“2027-01-01T00:05:00.000Z”,“%Y-%m-%dT%H:%M:%S.%fz”)), mu)
orbit = OrbitType.convertType(OrbitType.EQUINOCTIAL,orbit)
state = SpacecraftState(orbit)
mean_state = propagator.computeMeanState(state, velocity_pointing, propagator.getAllForceModels(),
float(propagation_data[‘computeMeanState’][‘convergence_threshold’]),
int(float(propagation_data[‘computeMeanState’][‘max_iteration’])))