use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.
the class HolmesFeatherstoneAttractionModel method gradient.
/**
* Compute the gradient of 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 gradient of the non-central part of the gravity field
* @exception OrekitException if position cannot be converted to central body frame
*/
public double[] gradient(final AbsoluteDate date, final Vector3D position, final double mu) throws OrekitException {
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];
final double[] pnm1 = 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 rho2 = x2 + y2;
final double rho = FastMath.sqrt(rho2);
// cos(theta), where theta is the polar angle
final double t = z / r;
// sin(theta), where theta is the polar angle
final double u = rho / r;
final double tOu = z / rho;
// compute distance powers
final double[] aOrN = createDistancePowersArray(provider.getAe() / r);
// compute longitude cosines/sines
final double[][] cosSinLambda = createCosSinArrays(position.getX() / rho, position.getY() / rho);
// outer summation over order
int index = 0;
double value = 0;
final double[] gradient = new double[3];
for (int m = degree; m >= 0; --m) {
// compute tesseral terms with derivatives
index = computeTesseral(m, degree, index, t, u, tOu, pnm0Plus2, pnm0Plus1, null, pnm0, pnm1, null);
if (m <= order) {
// compute contribution of current order to field (equation 5 of the paper)
// inner summation over degree, for fixed order
double sumDegreeS = 0;
double sumDegreeC = 0;
double dSumDegreeSdR = 0;
double dSumDegreeCdR = 0;
double dSumDegreeSdTheta = 0;
double dSumDegreeCdTheta = 0;
for (int n = FastMath.max(2, m); n <= degree; ++n) {
final double qSnm = aOrN[n] * harmonics.getNormalizedSnm(n, m);
final double qCnm = aOrN[n] * harmonics.getNormalizedCnm(n, m);
final double nOr = n / r;
final double s0 = pnm0[n] * qSnm;
final double c0 = pnm0[n] * qCnm;
final double s1 = pnm1[n] * qSnm;
final double c1 = pnm1[n] * qCnm;
sumDegreeS += s0;
sumDegreeC += c0;
dSumDegreeSdR -= nOr * s0;
dSumDegreeCdR -= nOr * c0;
dSumDegreeSdTheta += s1;
dSumDegreeCdTheta += c1;
}
// contribution to outer summation over order
// beware that we need to order gradient using the mathematical conventions
// compliant with the SphericalCoordinates class, so our lambda is its theta
// (and hence at index 1) and our theta is its phi (and hence at index 2)
final double sML = cosSinLambda[1][m];
final double cML = cosSinLambda[0][m];
value = value * u + sML * sumDegreeS + cML * sumDegreeC;
gradient[0] = gradient[0] * u + sML * dSumDegreeSdR + cML * dSumDegreeCdR;
gradient[1] = gradient[1] * u + m * (cML * sumDegreeS - sML * sumDegreeC);
gradient[2] = gradient[2] * u + sML * dSumDegreeSdTheta + cML * dSumDegreeCdTheta;
}
// rotate the recursion arrays
final double[] tmp = pnm0Plus2;
pnm0Plus2 = pnm0Plus1;
pnm0Plus1 = pnm0;
pnm0 = tmp;
}
// scale back
value = FastMath.scalb(value, SCALING);
gradient[0] = FastMath.scalb(gradient[0], SCALING);
gradient[1] = FastMath.scalb(gradient[1], SCALING);
gradient[2] = FastMath.scalb(gradient[2], SCALING);
// apply the global mu/r factor
final double muOr = mu / r;
value *= muOr;
gradient[0] = muOr * gradient[0] - value / r;
gradient[1] *= muOr;
gradient[2] *= muOr;
// convert gradient from spherical to Cartesian
return new SphericalCoordinates(position).toCartesianGradient(gradient);
}
use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.
the class SHMFormatReaderTest method testRegular03cNormalized.
@Test
public void testRegular03cNormalized() throws OrekitException {
Utils.setDataRoot("potential");
GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("eigen_cg03c_coef", false));
NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
Assert.assertEquals(TideSystem.TIDE_FREE, provider.getTideSystem());
AbsoluteDate refDate = new AbsoluteDate("1997-01-01T12:00:00", TimeScalesFactory.getTT());
Assert.assertEquals(refDate, provider.getReferenceDate());
AbsoluteDate date = new AbsoluteDate("2011-05-01T01:02:03", TimeScalesFactory.getTT());
Assert.assertEquals(date.durationFrom(refDate), provider.getOffset(date), Precision.SAFE_MIN);
NormalizedSphericalHarmonics harmonics = provider.onDate(date);
double offset = date.durationFrom(refDate);
double offsetYear = offset / Constants.JULIAN_YEAR;
Assert.assertEquals(0.957201462136e-06 + offsetYear * 0.490000000000e-11, harmonics.getNormalizedCnm(3, 0), 1.0e-15);
Assert.assertEquals(0.174786174485e-06, harmonics.getNormalizedCnm(5, 5), 1.0e-15);
Assert.assertEquals(0.0, harmonics.getNormalizedSnm(4, 0), 1.0e-15);
Assert.assertEquals(0.308834784975e-06, harmonics.getNormalizedSnm(4, 4), 1.0e-15);
Assert.assertEquals(0.3986004415E+15, provider.getMu(), 0);
Assert.assertEquals(0.6378136460E+07, provider.getAe(), 0);
}
use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.
the class OceanTidesFieldTest method testDeltaCnmSnm.
@Test
public void testDeltaCnmSnm() throws OrekitException {
// this is an arbitrarily truncated model, limited to 4x4 and with only a few waves
List<OceanTidesWave> waves = getWaves(4, 4, 55565, 56554, 85455, 135655, 273555);
UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate date = new AbsoluteDate(2003, 5, 6, 13, 43, 32.125, utc);
OceanTidesField tidesField = new OceanTidesField(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, Constants.EIGEN5C_EARTH_MU, waves, IERSConventions.IERS_2010.getNutationArguments(ut1), null);
NormalizedSphericalHarmonics harmonics = tidesField.onDate(date);
double[][] refDeltaCnm = new double[][] { { 0.0, 0.0, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0, 0.0, 0.0 }, { -4.812565797928061E-11, -4.1748378190052583E-11, 7.013273986245356E-11, 0.0, 0.0 }, { -2.5341227608443308E-11, 9.76515813742254E-11, -1.21931214469994E-10, 1.3179722429471184E-10, 0.0 }, { -2.7496974839179478E-11, 8.419627031293907E-11, 6.56546217101275E-11, -3.375298928713117E-11, -7.588006744166988E-11 } };
double[][] refDeltaSnm = new double[][] { { 0.0, 0.0, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0, 0.0, 0.0 }, { -1.168129177701461E-10, 5.646187590518608E-12, 1.742233297668071E-10, 0.0, 0.0 }, { -6.586546350227345E-11, -8.032186864783105E-11, -3.118910148495339E-11, 1.0566857199592183E-10, 0.0 }, { 7.665313525684617E-11, 7.37884528812169E-11, -1.3085142873419844E-10, -1.5813709543115768E-10, 1.770903634801541E-10 } };
for (int n = 0; n < refDeltaCnm.length; ++n) {
double threshold = 4.0e-17;
for (int m = 0; m <= n; ++m) {
Assert.assertEquals(refDeltaCnm[n][m], harmonics.getNormalizedCnm(n, m), threshold);
Assert.assertEquals(refDeltaSnm[n][m], harmonics.getNormalizedSnm(n, m), threshold);
}
}
}
use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.
the class CachedNormalizedSphericalHarmonicsProviderTest method testReverseEntryGeneration.
@Test
public void testReverseEntryGeneration() throws OrekitException {
// setup
// generate points on grid with date as the origin
cache.onDate(date);
// sample before the current cached values
AbsoluteDate sampleDate = date.shiftedBy(-step * 3);
NormalizedSphericalHarmonics expected = raw.onDate(sampleDate);
// action
NormalizedSphericalHarmonics actual = cache.onDate(sampleDate);
// verify
double tol = Precision.EPSILON;
for (int n = 0; n < raw.getMaxDegree(); n++) {
for (int m = 0; m < n; m++) {
Assert.assertEquals(expected.getNormalizedCnm(n, m), actual.getNormalizedCnm(n, m), tol);
Assert.assertEquals(expected.getNormalizedSnm(n, m), actual.getNormalizedSnm(n, m), tol);
}
}
}
use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.
the class GravityFieldFactoryTest method testNormalizer.
@Test
public void testNormalizer() throws OrekitException {
Utils.setDataRoot("potential/icgem-format");
final double shift = 1.23456e8;
NormalizedSphericalHarmonicsProvider ref = GravityFieldFactory.getNormalizedProvider(5, 5);
NormalizedSphericalHarmonics refHarmonics = ref.onDate(ref.getReferenceDate().shiftedBy(shift));
UnnormalizedSphericalHarmonicsProvider unnormalized = GravityFieldFactory.getUnnormalizedProvider(5, 5);
NormalizedSphericalHarmonicsProvider normalized = GravityFieldFactory.getNormalizedProvider(unnormalized);
NormalizedSphericalHarmonics normalizedHarmonics = normalized.onDate(normalized.getReferenceDate().shiftedBy(shift));
Assert.assertEquals(ref.getMaxDegree(), normalized.getMaxDegree());
Assert.assertEquals(ref.getMaxOrder(), normalized.getMaxOrder());
Assert.assertEquals(ref.getReferenceDate(), normalized.getReferenceDate());
Assert.assertEquals(ref.getAe(), normalized.getAe(), FastMath.ulp(ref.getAe()));
Assert.assertEquals(ref.getMu(), normalized.getMu(), FastMath.ulp(ref.getMu()));
Assert.assertEquals(ref.getOffset(AbsoluteDate.GPS_EPOCH), normalized.getOffset(AbsoluteDate.GPS_EPOCH), FastMath.ulp(ref.getOffset(AbsoluteDate.GPS_EPOCH)));
for (int i = 0; i <= 5; ++i) {
for (int j = 0; j <= i; ++j) {
double cRef = refHarmonics.getNormalizedCnm(i, j);
double cTest = normalizedHarmonics.getNormalizedCnm(i, j);
Assert.assertEquals(cRef, cTest, FastMath.ulp(cRef));
double sRef = refHarmonics.getNormalizedSnm(i, j);
double sTest = normalizedHarmonics.getNormalizedSnm(i, j);
Assert.assertEquals(sRef, sTest, FastMath.ulp(sRef));
}
}
}
Aggregations