Does the current implementation of the spherical harmonics force through the Holmes-Featherstone method allow for accelerations to be calculated above the poles, i.e, with x = y = 0 and z = r? The cited paper in the code mentions that this is not possible with the formulation presented, and the code seems to not take this into account.
Hi rjpower4,
I believe you are correct. I looked at this issue once several years ago and implemented an alternative formulation that does not have a singularity at the poles. The alternative was significantly slower and was only different from Holmes-Featherstone very close to the poles. So close that i regarded the alternative to be of no practical use. I didn’t evaluate derivatives, perhaps they break down earlier.
What is your use case? There are several papers that describe alternative formulations to remove the singularity at the poles.
Regards,
Evan
Hey Evan,
Thanks for the reply. I am looking at targeting low-thrust transfers to low polar orbits. Therefore, I will (likely) at some point be passing over the poles – or at least very very close to the poles. Consequently, I want to make sure I have a high confidence that my non-spherical effects won’t break down during the transfer.
RJ
Hi @rjpower4,
As long as you’re above the Brillouin sphere (circumscribing sphere of your Earth ellipsoid) there should be no problem.
With a classical ellipsoid like WGS84 it means staying above 20km over the poles.
Maxime
I am familiar with that limitation of the convergence of the series, however, I don’t think that the algorithm as implemented can handle any points directly above the poles.
Here is the code for the first several lines of nonCentralPart in the class HolmesFeatherstoneAttractionModel:
/** Compute the non-central part of the gravity field.
* @param date current date
* @param position position at which gravity field is desired in body frame
* @param mu central attraction coefficient to use
* @return value of the non-central part of the gravity field
*/
public double nonCentralPart(final AbsoluteDate date, final Vector3D position, final double mu) {
final int degree = provider.getMaxDegree();
final int order = provider.getMaxOrder();
final NormalizedSphericalHarmonics harmonics = provider.onDate(date);
// allocate the columns for recursion
double[] pnm0Plus2 = new double[degree + 1];
double[] pnm0Plus1 = new double[degree + 1];
double[] pnm0 = new double[degree + 1];
// compute polar coordinates
final double x = position.getX();
final double y = position.getY();
final double z = position.getZ();
final double x2 = x * x;
final double y2 = y * y;
final double z2 = z * z;
final double r2 = x2 + y2 + z2;
final double r = FastMath.sqrt (r2);
final double rho = FastMath.sqrt(x2 + y2);
final double t = z / r; // cos(theta), where theta is the polar angle
final double u = rho / r; // sin(theta), where theta is the polar angle
final double tOu = z / rho;
// .... more code
As you can see, if the position of the spacecraft is such that there are no components in the x and y directions but any component in the z direction (even above the Brillouin sphere radius) rho will be zero and thus tOu will be the result of a divide-by-zero. This is a limitation in the methodology used by Holmes and Featherstone in their paper; and, while they mention circumventing measures are available, they do not provide these measures.
I wasn’t aware of this limitation thank you.
It seems I missed the point on this question, sorry about that.