Spinning satellite determination

Orekit Python Wrapper usage
1

Hello everyone,

Finishing my topic on station keeping maneuver last week, I’m now hadling a new project quite different this time.

I would like to study the impact of spinning satellite (without ADCS) and the impact on power generation through solar pannel. An uncontrolled spin, combined with specific orbits has in fact an effet on the power, generated by the solar pannels. (For a few weeks the solar pannels could be approximately align with the sun but it can be the opposite a few weeks after.)

The aims of the study is to :

  • Propagate this state and obtain a rotation rate through time (To see how the natural spin is evolving)
  • Propagate this state and obtain the Beta angle between a specific spacecraft axis and the sun. If the beta angle is greater than a threshold, it would mean that the solar pannels aren’t generating enough power, which could be dramatic for the spacecraft survival. The idea would be to register the date of such event.

Do you have any idea on how to handle it, and if it’s possible ?I’m gonna try on my own to but I’m struggling for the moment.

Thank you for your time !

Cheers,

Emilien

1 Like
2

You should look at the TorqueFree attitude provider which was designed to represent attitudes for non-controlled objects.

1 Like
3

Hello Luc,

Thank you for your advice, I tried to propagate the spin by using this class. Even if I obtained outputs by configuring the spin through TorqueFree, I’m not sure of the way to add this attitude to the propagator.
You can find below the code I used.

from imports import *
from Plotting import Handler_For_Date_Definition
from org.hipparchus.geometry.euclidean.threed import Rotation, Vector3D
from org.orekit.attitudes import TorqueFree, Attitude, Inertia, InertiaAxis


######################## Environment Definition #########################

def Env (area):
    GravityProvider = GravityFieldFactory.getNormalizedProvider(10, 10)
    earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                            Constants.WGS84_EARTH_FLATTENING,
                            FramesFactory.getITRF(IERSConventions.IERS_2010, True))

    inertialFrame = FramesFactory.getEME2000()
    sun = CelestialBodyFactory.getSun()
    itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
    earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf)
    dragCoef = 2.2
    drag = IsotropicDrag(area, dragCoef)
    strengthLevel = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
    supportedNames = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
    inputParameters_Drag = MarshallSolarActivityFutureEstimation(supportedNames, strengthLevel)
    atmosphere = NRLMSISE00(inputParameters_Drag, sun, earth)
    orbitType = OrbitType.CARTESIAN
    isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(area, 1.8)
    solarRadiationPressure = SolarRadiationPressure(sun, earth,
                                                isotropicRadiationSingleCoeff)
    return earth, orbitType, inertialFrame, atmosphere, drag, GravityProvider, solarRadiationPressure

######################## Propagator Definition #########################
Altitude_reentry = 200.0 * 1000.0
def propagator_definition(initialOrbit,orbitType, inertialFrame,satellite_mass, earth,atmosphere, drag, GravityProvider, solarRadiationPressure):
    minStep = 0.0000000001
    maxStep = 100.0
    initStep = 60.0
    positionTolerance = 1.0
    tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit,orbitType)
    integrator = DormandPrince853Integrator(minStep, maxStep,
                                            JArray_double.cast_(tolerances[0]), JArray_double.cast_(tolerances[1]))

    integrator.setInitialStepSize(initStep)
    initialState = SpacecraftState(initialOrbit, satellite_mass)
    propagator = NumericalPropagator(integrator)
    propagator.setOrbitType(orbitType)
    propagator.setInitialState(initialState)
 
    propagator.setAttitudeProvider(LofOffset(inertialFrame, LOFType.VNC))
    altitudeDetect =  AltitudeDetector(Altitude_reentry, earth).withHandler(StopOnEvent())
    propagator.addEventDetector(altitudeDetect)
    propagator.addForceModel(DragForce(atmosphere, drag))
    propagator.addForceModel(solarRadiationPressure)
    #propagator.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), GravityProvider))

    return propagator

############################ Date Definition ######################################
##########################################s##########################################

simDur = 3600 *24*60.0            # Simulation Duration in seconds (Duration of the propagation)
# Here 60 days


initial_date = AbsoluteDate(2024, 7, 11, 0, 0, 0.0, TimeScalesFactory.getUTC())


############################ Orbit Definition ######################################
####################################################################################

ra = 580.0 * 1000  # Apogee in meters
rp = 580.0 * 1000  # Perigee in meters
i = 62.0  # Inclination (degres)
omega = 20.0  # Perigee argument in radians
raan = 10.0    # Right ascension of ascending node in radians
lv = 0.0     # True anomaly in radians
a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0    
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a

mu = Constants.EGM96_EARTH_MU 

# Setting the initial orbit : 
initialOrbit = KeplerianOrbit(a, e, math.radians(i), math.radians(omega), math.radians(raan), math.radians(lv),
                              PositionAngleType.TRUE, FramesFactory.getEME2000(), initial_date, mu)


############################ Satellite Definition ######################################
########################################################################################

satellite_mass = 3.9        
area = 0.0635

##################################################################################################################

earth, orbitType, inertialFrame, atmosphere, drag, GravityProvider,solarRadiationPressure = Env(area)


############################ Integrator and propagator Definition ################################################


propagator = propagator_definition(initialOrbit,orbitType, inertialFrame, satellite_mass, earth, atmosphere, 
                                   drag, GravityProvider,solarRadiationPressure)
# Setting up the propagator

inertia_matrix = np.array([[0.0639, -0.00001, -0.0041],
                           [-0.001, 0.0573, 0.00004],
                           [-0.001, 0.00004, 0.0253]])

i_xx, i_yy, i_zz = 0.0639, 0.0573, 0.0253


a1 = Vector3D(1.0, 0.0, 0.0)  # Direction of the x axis
a2 = Vector3D(0.0, 1.0, 0.0)  # Direction of the y axis
a3 = Vector3D(0.0, 0.0, 1.0)  # Direction of the z axis


iA1 = InertiaAxis(i_xx, a1)
iA2 = InertiaAxis(i_yy, a2)
iA3 = InertiaAxis(i_zz, a3)

inertia = Inertia(iA1, iA2, iA3)

initial_spin = Vector3D(math.radians(130.0), math.radians(-10.0), math.radians(1.0))
# An inital spin of 130°/s around x   /  -10°/s around y and 1.0°/s around z is considered

initial_rotation = Rotation.IDENTITY    
initialAttitude = Attitude(initial_date, FramesFactory.getEME2000(), initial_rotation, initial_spin, Vector3D.ZERO)


torque_free = TorqueFree(initialAttitude, inertia)




time_step = 30.0
Handler = Handler_For_Date_Definition(initial_date, simDur, time_step, satellite_mass)     # The handler is used is another file and used for user convenience purpose.
Propagator.cast_(propagator).setStepHandler(time_step,Handler)


############################ Propagate #################################################################
########################################################################################################


final_state = propagator.propagate(initial_date, initial_date.shiftedBy(simDur))

future_date = initial_date.shiftedBy(simDur)

                                            # Get the attitude at the future date
pvProv = None                               # Position-velocity provider (you may need to define this)
frame = FramesFactory.getEME2000()          # Reference frame for attitude computation
future_attitude = torque_free.getAttitude(pvProv, future_date, frame)
future_spin_vector = future_attitude.getSpin()
                                            # Convert spin to degrees per second for better readability
future_spin_degrees = Vector3D(
    math.degrees(future_spin_vector.getX()),
    math.degrees(future_spin_vector.getY()),
    math.degrees(future_spin_vector.getZ())
)

print(f"At the end of the propagation, the spin is : {future_spin_degrees}")

With this code, after a propagation of 60 days, and for an initial spin of (130°/s, -10°/s, 1°/s) around x,y,z, I obtained the following result :
At the end of the propagation, the spin is : {130,1302503837; -7,3771398375; -4,3185620109}

Passing from 1°/s for z axis to -4.3°/s seems a bit strange to me.

Do you have any ideas and see any problems ?

Thank you !

Emilien

4

TorqueFree is an attitude provider, so you could directly register it using propagator.setAttitudeProvider.

1 Like
5

Hi,

I updated the code by implementing propagator.setAttitudeprovider, as you suggested.

Here is the updated code :


propagator = propagator_definition(initialOrbit,orbitType, inertialFrame, satellite_mass, earth, atmosphere, 
                                   drag, GravityProvider,solarRadiationPressure)
# Setting up the propagator

inertia_matrix = np.array([[0.0639, -0.00001, -0.0041],
                           [-0.001, 0.0573, 0.00004],
                           [-0.001, 0.00004, 0.0253]])

i_xx, i_yy, i_zz = 0.0639, 0.0573, 0.0253


a1 = Vector3D(1.0, 0.0, 0.0)  # Direction of the x axis
a2 = Vector3D(0.0, 1.0, 0.0)  # Direction of the y axis
a3 = Vector3D(0.0, 0.0, 1.0)  # Direction of the z axis


iA1 = InertiaAxis(i_xx, a1)
iA2 = InertiaAxis(i_yy, a2)
iA3 = InertiaAxis(i_zz, a3)

inertia = Inertia(iA1, iA2, iA3)


initAttitude = Attitude(initial_date, FramesFactory.getEME2000(), Rotation(RotationOrder.XYZ, RotationConvention.FRAME_TRANSFORM, 0.0, 0.0, 0.0),
                        Vector3D(math.radians(130.0), math.radians(-10.0), 30 * math.radians(1.0)), Vector3D(0.0, 0.0, 0.0))

propagator.setAttitudeProvider(TorqueFree(initAttitude, inertia))


time_step = 30.0
Handler = Handler_For_Date_Definition(initial_date, simDur, time_step, satellite_mass)     # The handler is used is another file and used for user convenience purpose.
Propagator.cast_(propagator).setStepHandler(time_step,Handler)


############################ Propagate #################################################################
########################################################################################################


final_state = propagator.propagate(initial_date, initial_date.shiftedBy(simDur))

Final_attitude = final_state.getAttitude()
print("first rot angle = ", math.degrees(Final_attitude.getRotation().getAngles(RotationOrder.XYZ, RotationConvention.FRAME_TRANSFORM)[0]))
print("second rot angle = ", math.degrees(Final_attitude.getRotation().getAngles(RotationOrder.XYZ, RotationConvention.FRAME_TRANSFORM)[1]))
print("third rot angle = ", math.degrees(Final_attitude.getRotation().getAngles(RotationOrder.XYZ, RotationConvention.FRAME_TRANSFORM)[2]))

Which leads to the following outputs :
first rot angle = 131.15493458054863 °/s
second rot angle = 5.22361515876958 °/s
third rot angle = 15.484435657499523 °/s

which is quite different. However, I have data that shows the decreasing rate of the spin around the x-axis, which is obviously not the case in here. I think the frame definition could be an issue in my code but I’m not sure.

Thanks for your help !

Cheers,

Emilien