Examples with NormalizedSphericalHarmonics org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics used on opensource projects

Example 1 with NormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModelTest method gradientHessian.

private GradientHessian gradientHessian(final HolmesFeatherstoneAttractionModel hfModel, final AbsoluteDate date, final Vector3D position) throws OrekitException {
    try {
        java.lang.reflect.Field providerField = HolmesFeatherstoneAttractionModel.class.getDeclaredField("provider");
        providerField.setAccessible(true);
        NormalizedSphericalHarmonicsProvider provider = (NormalizedSphericalHarmonicsProvider) providerField.get(hfModel);
        java.lang.reflect.Method createDistancePowersArrayMethod = HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("createDistancePowersArray", Double.TYPE);
        createDistancePowersArrayMethod.setAccessible(true);
        java.lang.reflect.Method createCosSinArraysMethod = HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("createCosSinArrays", Double.TYPE, Double.TYPE);
        createCosSinArraysMethod.setAccessible(true);
        java.lang.reflect.Method computeTesseralMethod = HolmesFeatherstoneAttractionModel.class.getDeclaredMethod("computeTesseral", Integer.TYPE, Integer.TYPE, Integer.TYPE, Double.TYPE, Double.TYPE, Double.TYPE, double[].class, double[].class, double[].class, double[].class, double[].class, double[].class);
        computeTesseralMethod.setAccessible(true);
        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];
        double[] pnm1Plus1 = new double[degree + 1];
        double[] pnm1 = new double[degree + 1];
        final double[] pnm2 = 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 = (double[]) createDistancePowersArrayMethod.invoke(hfModel, provider.getAe() / r);
        // compute longitude cosines/sines
        final double[][] cosSinLambda = (double[][]) createCosSinArraysMethod.invoke(hfModel, position.getX() / rho, position.getY() / rho);
        // outer summation over order
        int index = 0;
        double value = 0;
        final double[] gradient = new double[3];
        final double[][] hessian = new double[3][3];
        for (int m = degree; m >= 0; --m) {
            // compute tesseral terms
            index = ((Integer) computeTesseralMethod.invoke(hfModel, m, degree, index, t, u, tOu, pnm0Plus2, pnm0Plus1, pnm1Plus1, pnm0, pnm1, pnm2)).intValue();
            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;
                double d2SumDegreeSdRdR = 0;
                double d2SumDegreeSdRdTheta = 0;
                double d2SumDegreeSdThetadTheta = 0;
                double d2SumDegreeCdRdR = 0;
                double d2SumDegreeCdRdTheta = 0;
                double d2SumDegreeCdThetadTheta = 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 nnP1Or2 = nOr * (n + 1) / 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;
                    final double s2 = pnm2[n] * qSnm;
                    final double c2 = pnm2[n] * qCnm;
                    sumDegreeS += s0;
                    sumDegreeC += c0;
                    dSumDegreeSdR -= nOr * s0;
                    dSumDegreeCdR -= nOr * c0;
                    dSumDegreeSdTheta += s1;
                    dSumDegreeCdTheta += c1;
                    d2SumDegreeSdRdR += nnP1Or2 * s0;
                    d2SumDegreeSdRdTheta -= nOr * s1;
                    d2SumDegreeSdThetadTheta += s2;
                    d2SumDegreeCdRdR += nnP1Or2 * c0;
                    d2SumDegreeCdRdTheta -= nOr * c1;
                    d2SumDegreeCdThetadTheta += c2;
                }
                // contribution to outer summation over order
                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;
                hessian[0][0] = hessian[0][0] * u + sML * d2SumDegreeSdRdR + cML * d2SumDegreeCdRdR;
                hessian[1][0] = hessian[1][0] * u + m * (cML * dSumDegreeSdR - sML * dSumDegreeCdR);
                hessian[2][0] = hessian[2][0] * u + sML * d2SumDegreeSdRdTheta + cML * d2SumDegreeCdRdTheta;
                hessian[1][1] = hessian[1][1] * u - m * m * (sML * sumDegreeS + cML * sumDegreeC);
                hessian[2][1] = hessian[2][1] * u + m * (cML * dSumDegreeSdTheta - sML * dSumDegreeCdTheta);
                hessian[2][2] = hessian[2][2] * u + sML * d2SumDegreeSdThetadTheta + cML * d2SumDegreeCdThetadTheta;
            }
            // rotate the recursion arrays
            final double[] tmp0 = pnm0Plus2;
            pnm0Plus2 = pnm0Plus1;
            pnm0Plus1 = pnm0;
            pnm0 = tmp0;
            final double[] tmp1 = pnm1Plus1;
            pnm1Plus1 = pnm1;
            pnm1 = tmp1;
        }
        // scale back
        value = FastMath.scalb(value, SCALING);
        for (int i = 0; i < 3; ++i) {
            gradient[i] = FastMath.scalb(gradient[i], SCALING);
            for (int j = 0; j <= i; ++j) {
                hessian[i][j] = FastMath.scalb(hessian[i][j], SCALING);
            }
        }
        // apply the global mu/r factor
        final double muOr = provider.getMu() / r;
        value *= muOr;
        gradient[0] = muOr * gradient[0] - value / r;
        gradient[1] *= muOr;
        gradient[2] *= muOr;
        hessian[0][0] = muOr * hessian[0][0] - 2 * gradient[0] / r;
        hessian[1][0] = muOr * hessian[1][0] - gradient[1] / r;
        hessian[2][0] = muOr * hessian[2][0] - gradient[2] / r;
        hessian[1][1] *= muOr;
        hessian[2][1] *= muOr;
        hessian[2][2] *= muOr;
        // convert gradient and Hessian from spherical to Cartesian
        final SphericalCoordinates sc = new SphericalCoordinates(position);
        return new GradientHessian(sc.toCartesianGradient(gradient), sc.toCartesianHessian(hessian, gradient));
    } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException | InvocationTargetException | NoSuchMethodException e) {
        return null;
    }
}

Example 2 with NormalizedSphericalHarmonics

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
 * @param <T> type of field used
 * @return gradient of the non-central part of the gravity field
 * @exception OrekitException if position cannot be converted to central body frame
 */
public <T extends RealFieldElement<T>> T[] gradient(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position, final T mu) throws OrekitException {
    final int degree = provider.getMaxDegree();
    final int order = provider.getMaxOrder();
    final NormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());
    final T zero = date.getField().getZero();
    // allocate the columns for recursion
    T[] pnm0Plus2 = MathArrays.buildArray(date.getField(), degree + 1);
    T[] pnm0Plus1 = MathArrays.buildArray(date.getField(), degree + 1);
    T[] pnm0 = MathArrays.buildArray(date.getField(), degree + 1);
    final T[] pnm1 = MathArrays.buildArray(date.getField(), degree + 1);
    // compute polar coordinates
    final T x = position.getX();
    final T y = position.getY();
    final T z = position.getZ();
    final T x2 = x.multiply(x);
    final T y2 = y.multiply(y);
    final T z2 = z.multiply(z);
    final T r2 = x2.add(y2).add(z2);
    final T r = r2.sqrt();
    final T rho2 = x2.add(y2);
    final T rho = rho2.sqrt();
    // cos(theta), where theta is the polar angle
    final T t = z.divide(r);
    // sin(theta), where theta is the polar angle
    final T u = rho.divide(r);
    final T tOu = z.divide(rho);
    // compute distance powers
    final T[] aOrN = createDistancePowersArray(r.reciprocal().multiply(provider.getAe()));
    // compute longitude cosines/sines
    final T[][] cosSinLambda = createCosSinArrays(rho.reciprocal().multiply(position.getX()), rho.reciprocal().multiply(position.getY()));
    // outer summation over order
    int index = 0;
    T value = zero;
    final T[] gradient = MathArrays.buildArray(zero.getField(), 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
            T sumDegreeS = zero;
            T sumDegreeC = zero;
            T dSumDegreeSdR = zero;
            T dSumDegreeCdR = zero;
            T dSumDegreeSdTheta = zero;
            T dSumDegreeCdTheta = zero;
            for (int n = FastMath.max(2, m); n <= degree; ++n) {
                final T qSnm = aOrN[n].multiply(harmonics.getNormalizedSnm(n, m));
                final T qCnm = aOrN[n].multiply(harmonics.getNormalizedCnm(n, m));
                final T nOr = r.reciprocal().multiply(n);
                final T s0 = pnm0[n].multiply(qSnm);
                final T c0 = pnm0[n].multiply(qCnm);
                final T s1 = pnm1[n].multiply(qSnm);
                final T c1 = pnm1[n].multiply(qCnm);
                sumDegreeS = sumDegreeS.add(s0);
                sumDegreeC = sumDegreeC.add(c0);
                dSumDegreeSdR = dSumDegreeSdR.subtract(nOr.multiply(s0));
                dSumDegreeCdR = dSumDegreeCdR.subtract(nOr.multiply(c0));
                dSumDegreeSdTheta = dSumDegreeSdTheta.add(s1);
                dSumDegreeCdTheta = dSumDegreeCdTheta.add(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 T sML = cosSinLambda[1][m];
            final T cML = cosSinLambda[0][m];
            value = value.multiply(u).add(sML.multiply(sumDegreeS)).add(cML.multiply(sumDegreeC));
            gradient[0] = gradient[0].multiply(u).add(sML.multiply(dSumDegreeSdR)).add(cML.multiply(dSumDegreeCdR));
            gradient[1] = gradient[1].multiply(u).add(cML.multiply(sumDegreeS).subtract(sML.multiply(sumDegreeC)).multiply(m));
            gradient[2] = gradient[2].multiply(u).add(sML.multiply(dSumDegreeSdTheta)).add(cML.multiply(dSumDegreeCdTheta));
        }
        // rotate the recursion arrays
        final T[] tmp = pnm0Plus2;
        pnm0Plus2 = pnm0Plus1;
        pnm0Plus1 = pnm0;
        pnm0 = tmp;
    }
    // scale back
    value = value.scalb(SCALING);
    gradient[0] = gradient[0].scalb(SCALING);
    gradient[1] = gradient[1].scalb(SCALING);
    gradient[2] = gradient[2].scalb(SCALING);
    // apply the global mu/r factor
    final T muOr = r.reciprocal().multiply(mu);
    value = value.multiply(muOr);
    gradient[0] = muOr.multiply(gradient[0]).subtract(value.divide(r));
    gradient[1] = gradient[1].multiply(muOr);
    gradient[2] = gradient[2].multiply(muOr);
    // convert gradient from spherical to Cartesian
    // Cartesian coordinates
    // remaining spherical coordinates
    final T rPos = position.getNorm();
    // intermediate variables
    final T xPos = position.getX();
    final T yPos = position.getY();
    final T zPos = position.getZ();
    final T rho2Pos = x.multiply(x).add(y.multiply(y));
    final T rhoPos = rho2.sqrt();
    final T r2Pos = rho2.add(z.multiply(z));
    final T[][] jacobianPos = MathArrays.buildArray(zero.getField(), 3, 3);
    // row representing the gradient of r
    jacobianPos[0][0] = xPos.divide(rPos);
    jacobianPos[0][1] = yPos.divide(rPos);
    jacobianPos[0][2] = zPos.divide(rPos);
    // row representing the gradient of theta
    jacobianPos[1][0] = yPos.negate().divide(rho2Pos);
    jacobianPos[1][1] = xPos.divide(rho2Pos);
    // jacobian[1][2] is already set to 0 at allocation time
    // row representing the gradient of phi
    jacobianPos[2][0] = xPos.multiply(zPos).divide(rhoPos.multiply(r2Pos));
    jacobianPos[2][1] = yPos.multiply(zPos).divide(rhoPos.multiply(r2Pos));
    jacobianPos[2][2] = rhoPos.negate().divide(r2Pos);
    final T[] cartGradPos = MathArrays.buildArray(zero.getField(), 3);
    cartGradPos[0] = gradient[0].multiply(jacobianPos[0][0]).add(gradient[1].multiply(jacobianPos[1][0])).add(gradient[2].multiply(jacobianPos[2][0]));
    cartGradPos[1] = gradient[0].multiply(jacobianPos[0][1]).add(gradient[1].multiply(jacobianPos[1][1])).add(gradient[2].multiply(jacobianPos[2][1]));
    cartGradPos[2] = gradient[0].multiply(jacobianPos[0][2]).add(gradient[2].multiply(jacobianPos[2][2]));
    return cartGradPos;
}

Example 3 with NormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModel method nonCentralPart.

/**
 * 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
 * @exception OrekitException if position cannot be converted to central body frame
 */
public double nonCentralPart(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];
    // 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);
    // 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;
    for (int m = degree; m >= 0; --m) {
        // compute tesseral terms without derivatives
        index = computeTesseral(m, degree, index, t, u, tOu, pnm0Plus2, pnm0Plus1, null, pnm0, null, 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;
            for (int n = FastMath.max(2, m); n <= degree; ++n) {
                sumDegreeS += pnm0[n] * aOrN[n] * harmonics.getNormalizedSnm(n, m);
                sumDegreeC += pnm0[n] * aOrN[n] * harmonics.getNormalizedCnm(n, m);
            }
            // contribution to outer summation over order
            value = value * u + cosSinLambda[1][m] * sumDegreeS + cosSinLambda[0][m] * sumDegreeC;
        }
        // rotate the recursion arrays
        final double[] tmp = pnm0Plus2;
        pnm0Plus2 = pnm0Plus1;
        pnm0Plus1 = pnm0;
        pnm0 = tmp;
    }
    // scale back
    value = FastMath.scalb(value, SCALING);
    // apply the global mu/r factor
    return mu * value / r;
}

Example 4 with NormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.

the class SolidTidesFieldTest method testInterpolationAccuracy.

@Test
public void testInterpolationAccuracy() throws OrekitException {
    // The shortest periods are slightly below one half day for the tidal waves
    // considered here. This implies the sampling rate should be fast enough.
    // The tuning parameters we have finally settled correspond to a two hours
    // sample containing 12 points (i.e. one new point is computed every 10 minutes).
    // The observed relative interpolation error with these settings are essentially
    // due to Runge phenomenon at points sampling rate. Plotting the errors shows
    // singular peaks pointing out of merely numerical noise.
    final IERSConventions conventions = IERSConventions.IERS_2010;
    Frame itrf = FramesFactory.getITRF(conventions, true);
    TimeScale utc = TimeScalesFactory.getUTC();
    UT1Scale ut1 = TimeScalesFactory.getUT1(conventions, true);
    NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getConstantNormalizedProvider(5, 5);
    SolidTidesField raw = new SolidTidesField(conventions.getLoveNumbers(), conventions.getTideFrequencyDependenceFunction(ut1), conventions.getPermanentTide(), conventions.getSolidPoleTide(ut1.getEOPHistory()), itrf, gravityField.getAe(), gravityField.getMu(), gravityField.getTideSystem(), CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
    int step = 600;
    int nbPoints = 12;
    CachedNormalizedSphericalHarmonicsProvider interpolated = new CachedNormalizedSphericalHarmonicsProvider(raw, step, nbPoints, OrekitConfiguration.getCacheSlotsNumber(), 7 * Constants.JULIAN_DAY, 0.5 * Constants.JULIAN_DAY);
    // the following time range is located around the maximal observed error
    AbsoluteDate start = new AbsoluteDate(2003, 6, 12, utc);
    AbsoluteDate end = start.shiftedBy(3 * Constants.JULIAN_DAY);
    StreamingStatistics stat = new StreamingStatistics();
    for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(60)) {
        NormalizedSphericalHarmonics rawHarmonics = raw.onDate(date);
        NormalizedSphericalHarmonics interpolatedHarmonics = interpolated.onDate(date);
        for (int n = 2; n < 5; ++n) {
            for (int m = 0; m <= n; ++m) {
                if (n < 4 || m < 3) {
                    double cnmRaw = rawHarmonics.getNormalizedCnm(n, m);
                    double cnmInterp = interpolatedHarmonics.getNormalizedCnm(n, m);
                    double errorC = (cnmInterp - cnmRaw) / FastMath.abs(cnmRaw);
                    stat.addValue(errorC);
                    if (m > 0) {
                        double snmRaw = rawHarmonics.getNormalizedSnm(n, m);
                        double snmInterp = interpolatedHarmonics.getNormalizedSnm(n, m);
                        double errorS = (snmInterp - snmRaw) / FastMath.abs(snmRaw);
                        stat.addValue(errorS);
                    }
                }
            }
        }
    }
    Assert.assertEquals(0.0, stat.getMean(), 2.0e-12);
    Assert.assertTrue(stat.getStandardDeviation() < 2.0e-9);
    Assert.assertTrue(stat.getMin() > -9.0e-8);
    Assert.assertTrue(stat.getMax() < 2.2e-7);
}

Example 5 with NormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics in project Orekit by CS-SI.

the class SolidTidesFieldTest method testDeltaCnmSnm.

@Test
public void testDeltaCnmSnm() throws OrekitException {
    NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getConstantNormalizedProvider(8, 8);
    UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
    TimeScale utc = TimeScalesFactory.getUTC();
    AbsoluteDate date = new AbsoluteDate(2003, 5, 6, 13, 43, 32.125, utc);
    SolidTidesField tidesField = new SolidTidesField(IERSConventions.IERS_2010.getLoveNumbers(), IERSConventions.IERS_2010.getTideFrequencyDependenceFunction(ut1), IERSConventions.IERS_2010.getPermanentTide(), null, FramesFactory.getITRF(IERSConventions.IERS_2010, true), gravityField.getAe(), gravityField.getMu(), TideSystem.TIDE_FREE, CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
    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 }, { -2.6732289327355114E-9, 4.9078992447259636E-9, 3.5894110538262888E-9, 0.0, 0.0 }, // { -2.6598001259383122E-9,   4.907899244804072E-9,   3.5894110542365972E-9,         0.0 ,                 0.0  },
    { -1.290639603871307E-11, -9.287425756410472E-14, 8.356574033404024E-12, -2.2644465207860626E-12, 0.0 }, { 7.888138856951149E-12, -1.4422209452877158E-11, -6.815519349970944E-12, 0.0, 0.0 } };
    double[][] refDeltaSnm = new double[][] { { 0.0, 0.0, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0, 0.0, 0.0 }, { 0.0, 1.599927449004677E-9, 2.1815888169727694E-9, 0.0, 0.0 }, { 0.0, -4.6129961143785774E-14, 1.8097527720906976E-11, 1.633889224766215E-11, 0.0 }, { 0.0, -4.897228975221076E-12, -4.1034042689652575E-12, 0.0, 0.0 } };
    for (int n = 0; n < refDeltaCnm.length; ++n) {
        double threshold = (n == 2) ? 1.3e-17 : 1.0e-24;
        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);
        }
    }
}