I am writing a code that computes battery charge generated during Sun contact events. To do this, I add the event detectors, propagate my satellite, and then loop through each event from start to end every 60 seconds to compute necessary parameters to calculate the charge generated. I am attaching a summarized version of my code.
My issue happens when computing the true anomaly of my orbit. I need this angle in order to compute a rotation matrix to convert the normal direction of my solar panel to the ECI frame, according to the body frame I defined. Every 60 seconds, based on my circular orbit, the true anomaly should vary by about 4 degrees, however it makes very large jumps seemingly randomly, throwing off my calculation.
I am not entirely sure what the error is, but I also noticed my argument of perigee changes substantially every 60 seconds, even though I am using the KeplerianPropagator propagator, which I believe is a two-body problem. Is there any problems with my code? Note that the variable propTempSun is a new propagator created so that the events dont get repeated when I access PVCoordinates.
Thank you!
double a = 6978140;
double ecc = 0;
double inc = FastMath.toRadians(90);
double RAAN = FastMath.toRadians(45);
double omega = FastMath.toRadians(0);
double true_anomaly= FastMath.toRadians(0);
Orbit orbit = new KeplerianOrbit(a, ecc, inc, omega, RAAN, true_anomaly, PositionAngleType.TRUE, inertialFrame, start_date, mu);
KeplerianPropagator Prop = new KeplerianPropagator(orbit);
Prop.setAttitudeProvider(new NadirPointing(inertialFrame, earth));
// skipping event detector part //
for (AbsoluteDate extrapDate = start_date; extrapDate.compareTo(final_date) <= 0; extrapDate = extrapDate.shiftedBy(stepT)) {
Prop.getPVCoordinates(extrapDate, inertialFrame);
}
// a for loop then goes through every event, and a while loop checks every 60 seconds for each event
// the following code is inside the while loop
KeplerianOrbit newOrb = new KeplerianOrbit(propTempSun.getPVCoordinates(extrapDate, inertialFrame),inertialFrame,extrapDate,mu);
double newTrueAnomaly = newOrb.getTrueAnomaly();
extrapDate = extrapDate.shiftedBy(dt);
Computing Keplerian anomaly (be it mean, eccentric or true anomaly) on an orbit with exactly 0 eccentricity is error-prone. You are at the singularity point of Keplerian elements.
You should use CircularOrbit.
Also for computing extra parameters like this, I suggest you include the computation right within the propagation using additional parameters, by letting the propagator know about the equations that drive battery charge. In fact, the additional parameters were added years ago exactly for that purpose: computing battery charging.
I would like to make my algorithm ss generic as possible, so Id like for it to run with any input. What do you recommend I do in this case? Just an if statement to check if the eccentricity is 0?
Regarding battery, what exactly do you mean by adding the equations to the propagation?
I don’t get why you need the anomaly itself. Can’t you retrieve the direction of satellite axes, including solar panel from the attitude embedded in the SpacecraftState?
Also if you want your code to be as generic as possible, perhaps EquinoctialOrbit or even CartesianOrbit would avoid most singularities, but they won’t allow you to get the anomaly itself, if it is really what you need.
For context, I have as an input the direction of the solar panel in terms of a body frame. This body frame I have defined as: z parallel to position vector in ECI, x parallel to orbit angular momentum, and y perpendicular to both. In order to compute the charge, I need the direction of the solar panel in the ECI frame so I can then get the angle between it and the Sun vector from the spacecraft.
This means I need the rotation matrix to convert my body frame to the ECI. So I need the three orientation orbital elements, as well as the true anomaly, and an axis flip matrix. I haven’t defined the axes in Propagation, and I propagate attitude using Nadir.
No, you only need the position and velocity in ECI, which is something all orbits provide thanks to PVCoordinates pv = orbit.getPVCoordinates(). Then you can create a Rotation using Rotation r = new Rotation(Vector3D.PLUS_K, Vector3D.PLUS_J, pv.getPosition(), pv.getMomentum()), this gives a rotation that transforms vectors in body frame into vectors in ECI, so you apply this rotation to your solar panel direction in body frame and you get it in ECI. I did not check everything, but you have the rough algorithm.