Examples with UnnormalizedSphericalHarmonics org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics used on opensource projects

Example 1 with UnnormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics in project Orekit by CS-SI.

the class DSSTTesseral method computeNSum.

/**
 * Compute the n-SUM for potential derivatives components.
 *  @param date current date
 *  @param j resonant index <i>j</i>
 *  @param m resonant order <i>m</i>
 *  @param s d'Alembert characteristic <i>s</i>
 *  @param maxN maximum possible value for <i>n</i> index
 *  @param roaPow powers of R/a up to degree <i>n</i>
 *  @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials
 *  @param gammaMNS &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function
 *  @return Components of U<sub>n</sub> derivatives for fixed j, m, s
 * @throws OrekitException if some error occurred
 */
private double[][] computeNSum(final AbsoluteDate date, final int j, final int m, final int s, final int maxN, final double[] roaPow, final GHmsjPolynomials ghMSJ, final GammaMnsFunction gammaMNS) throws OrekitException {
    // spherical harmonics
    final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
    // Potential derivatives components
    double dUdaCos = 0.;
    double dUdaSin = 0.;
    double dUdhCos = 0.;
    double dUdhSin = 0.;
    double dUdkCos = 0.;
    double dUdkSin = 0.;
    double dUdlCos = 0.;
    double dUdlSin = 0.;
    double dUdAlCos = 0.;
    double dUdAlSin = 0.;
    double dUdBeCos = 0.;
    double dUdBeSin = 0.;
    double dUdGaCos = 0.;
    double dUdGaSin = 0.;
    // I^m
    @SuppressWarnings("unused") final int Im = I > 0 ? 1 : (m % 2 == 0 ? 1 : -1);
    // jacobi v, w, indices from 2.7.1-(15)
    final int v = FastMath.abs(m - s);
    final int w = FastMath.abs(m + s);
    // Initialise lower degree nmin = (Max(2, m, |s|)) for summation over n
    final int nmin = FastMath.max(FastMath.max(2, m), FastMath.abs(s));
    // Get the corresponding Hansen object
    final int sIndex = maxDegree + (j < 0 ? -s : s);
    final int jIndex = FastMath.abs(j);
    final HansenTesseralLinear hans = this.hansenObjects[sIndex][jIndex];
    // n-SUM from nmin to N
    for (int n = nmin; n <= maxN; n++) {
        // If (n - s) is odd, the contribution is null because of Vmns
        if ((n - s) % 2 == 0) {
            // Vmns coefficient
            final double vMNS = CoefficientsFactory.getVmns(m, n, s);
            // Inclination function Gamma and derivative
            final double gaMNS = gammaMNS.getValue(m, n, s);
            final double dGaMNS = gammaMNS.getDerivative(m, n, s);
            // Hansen kernel value and derivative
            final double kJNS = hans.getValue(-n - 1, chi);
            final double dkJNS = hans.getDerivative(-n - 1, chi);
            // Gjms, Hjms polynomials and derivatives
            final double gMSJ = ghMSJ.getGmsj(m, s, j);
            final double hMSJ = ghMSJ.getHmsj(m, s, j);
            final double dGdh = ghMSJ.getdGmsdh(m, s, j);
            final double dGdk = ghMSJ.getdGmsdk(m, s, j);
            final double dGdA = ghMSJ.getdGmsdAlpha(m, s, j);
            final double dGdB = ghMSJ.getdGmsdBeta(m, s, j);
            final double dHdh = ghMSJ.getdHmsdh(m, s, j);
            final double dHdk = ghMSJ.getdHmsdk(m, s, j);
            final double dHdA = ghMSJ.getdHmsdAlpha(m, s, j);
            final double dHdB = ghMSJ.getdHmsdBeta(m, s, j);
            // Jacobi l-index from 2.7.1-(15)
            final int l = FastMath.min(n - m, n - FastMath.abs(s));
            // Jacobi polynomial and derivative
            final DerivativeStructure jacobi = JacobiPolynomials.getValue(l, v, w, factory.variable(0, gamma));
            // Geopotential coefficients
            final double cnm = harmonics.getUnnormalizedCnm(n, m);
            final double snm = harmonics.getUnnormalizedSnm(n, m);
            // Common factors from expansion of equations 3.3-4
            final double cf_0 = roaPow[n] * Im * vMNS;
            final double cf_1 = cf_0 * gaMNS * jacobi.getValue();
            final double cf_2 = cf_1 * kJNS;
            final double gcPhs = gMSJ * cnm + hMSJ * snm;
            final double gsMhc = gMSJ * snm - hMSJ * cnm;
            final double dKgcPhsx2 = 2. * dkJNS * gcPhs;
            final double dKgsMhcx2 = 2. * dkJNS * gsMhc;
            final double dUdaCoef = (n + 1) * cf_2;
            final double dUdlCoef = j * cf_2;
            final double dUdGaCoef = cf_0 * kJNS * (jacobi.getValue() * dGaMNS + gaMNS * jacobi.getPartialDerivative(1));
            // dU / da components
            dUdaCos += dUdaCoef * gcPhs;
            dUdaSin += dUdaCoef * gsMhc;
            // dU / dh components
            dUdhCos += cf_1 * (kJNS * (cnm * dGdh + snm * dHdh) + h * dKgcPhsx2);
            dUdhSin += cf_1 * (kJNS * (snm * dGdh - cnm * dHdh) + h * dKgsMhcx2);
            // dU / dk components
            dUdkCos += cf_1 * (kJNS * (cnm * dGdk + snm * dHdk) + k * dKgcPhsx2);
            dUdkSin += cf_1 * (kJNS * (snm * dGdk - cnm * dHdk) + k * dKgsMhcx2);
            // dU / dLambda components
            dUdlCos += dUdlCoef * gsMhc;
            dUdlSin += -dUdlCoef * gcPhs;
            // dU / alpha components
            dUdAlCos += cf_2 * (dGdA * cnm + dHdA * snm);
            dUdAlSin += cf_2 * (dGdA * snm - dHdA * cnm);
            // dU / dBeta components
            dUdBeCos += cf_2 * (dGdB * cnm + dHdB * snm);
            dUdBeSin += cf_2 * (dGdB * snm - dHdB * cnm);
            // dU / dGamma components
            dUdGaCos += dUdGaCoef * gcPhs;
            dUdGaSin += dUdGaCoef * gsMhc;
        }
    }
    return new double[][] { { dUdaCos, dUdaSin }, { dUdhCos, dUdhSin }, { dUdkCos, dUdkSin }, { dUdlCos, dUdlSin }, { dUdAlCos, dUdAlSin }, { dUdBeCos, dUdBeSin }, { dUdGaCos, dUdGaSin } };
}

Example 2 with UnnormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics in project Orekit by CS-SI.

the class DSSTZonal method computeMeanElementsTruncations.

/**
 * Compute indices truncations for mean elements computations.
 * @param aux auxiliary elements
 * @throws OrekitException if an error occurs
 */
private void computeMeanElementsTruncations(final AuxiliaryElements aux) throws OrekitException {
    // Compute the max eccentricity power for the mean element rate expansion
    if (maxDegree == 2) {
        maxEccPowMeanElements = 0;
    } else {
        // Initializes specific parameters.
        initializeStep(aux);
        final UnnormalizedSphericalHarmonics harmonics = provider.onDate(aux.getDate());
        // Utilities for truncation
        final double ax2or = 2. * a / provider.getAe();
        double xmuran = provider.getMu() / a;
        // Set a lower bound for eccentricity
        final double eo2 = FastMath.max(0.0025, ecc / 2.);
        final double x2o2 = XX / 2.;
        final double[] eccPwr = new double[maxDegree + 1];
        final double[] chiPwr = new double[maxDegree + 1];
        final double[] hafPwr = new double[maxDegree + 1];
        eccPwr[0] = 1.;
        chiPwr[0] = X;
        hafPwr[0] = 1.;
        for (int i = 1; i <= maxDegree; i++) {
            eccPwr[i] = eccPwr[i - 1] * eo2;
            chiPwr[i] = chiPwr[i - 1] * x2o2;
            hafPwr[i] = hafPwr[i - 1] * 0.5;
            xmuran /= ax2or;
        }
        // Set highest power of e and degree of current spherical harmonic.
        maxEccPowMeanElements = 0;
        int n = maxDegree;
        // Loop over n
        do {
            // Set order of current spherical harmonic.
            int m = 0;
            // Loop over m
            do {
                // Compute magnitude of current spherical harmonic coefficient.
                final double cnm = harmonics.getUnnormalizedCnm(n, m);
                final double snm = harmonics.getUnnormalizedSnm(n, m);
                final double csnm = FastMath.hypot(cnm, snm);
                if (csnm == 0.)
                    break;
                // Set magnitude of last spherical harmonic term.
                double lastTerm = 0.;
                // Set current power of e and related indices.
                int nsld2 = (n - maxEccPowMeanElements - 1) / 2;
                int l = n - 2 * nsld2;
                // Loop over l
                double term = 0.;
                do {
                    // Compute magnitude of current spherical harmonic term.
                    if (m < l) {
                        term = csnm * xmuran * (fact[n - l] / (fact[n - m])) * (fact[n + l] / (fact[nsld2] * fact[nsld2 + l])) * eccPwr[l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) * (UpperBounds.getRnml(gamma, n, l, m, 1, I) + UpperBounds.getRnml(gamma, n, l, m, -1, I));
                    } else {
                        term = csnm * xmuran * (fact[n + m] / (fact[nsld2] * fact[nsld2 + l])) * eccPwr[l] * hafPwr[m - l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) * (UpperBounds.getRnml(gamma, n, m, l, 1, I) + UpperBounds.getRnml(gamma, n, m, l, -1, I));
                    }
                    // Is the current spherical harmonic term bigger than the truncation tolerance ?
                    if (term >= TRUNCATION_TOLERANCE) {
                        maxEccPowMeanElements = l;
                    } else {
                        // Is the current term smaller than the last term ?
                        if (term < lastTerm) {
                            break;
                        }
                    }
                    // Proceed to next power of e.
                    lastTerm = term;
                    l += 2;
                    nsld2--;
                } while (l < n);
                // Is the current spherical harmonic term bigger than the truncation tolerance ?
                if (term >= TRUNCATION_TOLERANCE) {
                    maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
                    return;
                }
                // Proceed to next order.
                m++;
            } while (m <= FastMath.min(n, maxOrder));
            // Proceed to next degree.
            xmuran *= ax2or;
            n--;
        } while (n > maxEccPowMeanElements + 2);
        maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
    }
}

Example 3 with UnnormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics in project Orekit by CS-SI.

the class DSSTZonal method computeUDerivatives.

/**
 * Compute the derivatives of the gravitational potential U [Eq. 3.1-(6)].
 *  <p>
 *  The result is the array
 *  [dU/da, dU/dk, dU/dh, dU/dα, dU/dβ, dU/dγ]
 *  </p>
 *  @param date current date
 *  @return potential derivatives
 *  @throws OrekitException if an error occurs in hansen computation
 */
private double[] computeUDerivatives(final AbsoluteDate date) throws OrekitException {
    final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
    // Reset U
    U = 0.;
    // Gs and Hs coefficients
    final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPowMeanElements);
    // Qns coefficients
    final double[][] Qns = CoefficientsFactory.computeQns(gamma, maxDegree, maxEccPowMeanElements);
    final double[] roaPow = new double[maxDegree + 1];
    roaPow[0] = 1.;
    for (int i = 1; i <= maxDegree; i++) {
        roaPow[i] = roa * roaPow[i - 1];
    }
    // Potential derivatives
    double dUda = 0.;
    double dUdk = 0.;
    double dUdh = 0.;
    double dUdAl = 0.;
    double dUdBe = 0.;
    double dUdGa = 0.;
    for (int s = 0; s <= maxEccPowMeanElements; s++) {
        // Initialize the Hansen roots
        this.hansenObjects[s].computeInitValues(X);
        // Get the current Gs coefficient
        final double gs = GsHs[0][s];
        // Compute Gs partial derivatives from 3.1-(9)
        double dGsdh = 0.;
        double dGsdk = 0.;
        double dGsdAl = 0.;
        double dGsdBe = 0.;
        if (s > 0) {
            // First get the G(s-1) and the H(s-1) coefficients
            final double sxgsm1 = s * GsHs[0][s - 1];
            final double sxhsm1 = s * GsHs[1][s - 1];
            // Then compute derivatives
            dGsdh = beta * sxgsm1 - alpha * sxhsm1;
            dGsdk = alpha * sxgsm1 + beta * sxhsm1;
            dGsdAl = k * sxgsm1 - h * sxhsm1;
            dGsdBe = h * sxgsm1 + k * sxhsm1;
        }
        // Kronecker symbol (2 - delta(0,s))
        final double d0s = (s == 0) ? 1 : 2;
        for (int n = s + 2; n <= maxDegree; n++) {
            // (n - s) must be even
            if ((n - s) % 2 == 0) {
                // Extract data from previous computation :
                final double kns = this.hansenObjects[s].getValue(-n - 1, X);
                final double dkns = this.hansenObjects[s].getDerivative(-n - 1, X);
                final double vns = Vns.get(new NSKey(n, s));
                final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
                final double coef1 = coef0 * Qns[n][s];
                final double coef2 = coef1 * kns;
                final double coef3 = coef2 * gs;
                // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
                final double dqns = Qns[n][s + 1];
                // Compute U
                U += coef3;
                // Compute dU / da :
                dUda += coef3 * (n + 1);
                // Compute dU / dEx
                dUdk += coef1 * (kns * dGsdk + k * XXX * gs * dkns);
                // Compute dU / dEy
                dUdh += coef1 * (kns * dGsdh + h * XXX * gs * dkns);
                // Compute dU / dAlpha
                dUdAl += coef2 * dGsdAl;
                // Compute dU / dBeta
                dUdBe += coef2 * dGsdBe;
                // Compute dU / dGamma
                dUdGa += coef0 * kns * dqns * gs;
            }
        }
    }
    // Multiply by -(μ / a)
    U *= -muoa;
    return new double[] { dUda * muoa / a, dUdk * -muoa, dUdh * -muoa, dUdAl * -muoa, dUdBe * -muoa, dUdGa * -muoa };
}

Example 4 with UnnormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics in project Orekit by CS-SI.

the class EGMFormatReaderTest method testReadUnnormalized.

@Test
public void testReadUnnormalized() throws OrekitException {
    Utils.setDataRoot("potential");
    GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("egm96_to5.ascii", true));
    UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
    UnnormalizedSphericalHarmonics harmonics = provider.onDate(AbsoluteDate.FUTURE_INFINITY);
    Assert.assertEquals(TideSystem.TIDE_FREE, provider.getTideSystem());
    int maxUlps = 1;
    checkValue(harmonics.getUnnormalizedCnm(3, 0), 3, 0, 0.957254173792E-06, maxUlps);
    checkValue(harmonics.getUnnormalizedCnm(5, 5), 5, 5, 0.174971983203E-06, maxUlps);
    checkValue(harmonics.getUnnormalizedSnm(4, 0), 4, 0, 0.0, maxUlps);
    checkValue(harmonics.getUnnormalizedSnm(4, 4), 4, 4, 0.308853169333E-06, maxUlps);
    double a = (-0.295301647654E-06);
    double b = 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2;
    double c = 2 * 11 / b;
    double result = a * FastMath.sqrt(c);
    Assert.assertEquals(result, harmonics.getUnnormalizedCnm(5, 4), 1.0e-20);
    a = -0.188560802735E-06;
    b = 8 * 7 * 6 * 5 * 4 * 3 * 2;
    c = 2 * 9 / b;
    result = a * FastMath.sqrt(c);
    Assert.assertEquals(result, harmonics.getUnnormalizedCnm(4, 4), 1.0e-20);
    Assert.assertEquals(1.0826266835531513e-3, -harmonics.getUnnormalizedCnm(2, 0), 1.0e-20);
    Assert.assertNull(provider.getReferenceDate());
    Assert.assertEquals(0, provider.getOffset(AbsoluteDate.J2000_EPOCH), Precision.SAFE_MIN);
    Assert.assertEquals(0, provider.getOffset(AbsoluteDate.MODIFIED_JULIAN_EPOCH), Precision.SAFE_MIN);
}

Example 5 with UnnormalizedSphericalHarmonics

use of org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics in project Orekit by CS-SI.

the class GRGSFormatReaderTest method testAdditionalColumn.

@Test
public void testAdditionalColumn() throws OrekitException {
    GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim5-c1.txt", true));
    UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
    Assert.assertEquals(TideSystem.UNKNOWN, 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);
    UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
    int maxUlps = 2;
    checkValue(harmonics.getUnnormalizedCnm(3, 0), date, 3, 0, 1997, 1, 1, 0.95857491635129E-06, 0.28175700027753E-11, maxUlps);
    checkValue(harmonics.getUnnormalizedCnm(5, 5), date, 5, 5, 1997, 1, 1, 0.17481512311600E-06, 0.0, maxUlps);
    checkValue(harmonics.getUnnormalizedSnm(4, 0), date, 4, 0, 1997, 1, 1, 0, 0, maxUlps);
    checkValue(harmonics.getUnnormalizedSnm(4, 4), date, 4, 4, 1997, 1, 1, 0.30882755318300E-06, 0, maxUlps);
    Assert.assertEquals(0.3986004415E+15, provider.getMu(), 0);
    Assert.assertEquals(0.6378136460E+07, provider.getAe(), 0);
}