Examples with Dfp org.hipparchus.dfp.Dfp used on opensource projects

Example 1 with Dfp

use of org.hipparchus.dfp.Dfp in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModelTest method testRelativeNumericPrecision.

@Test
public void testRelativeNumericPrecision() throws OrekitException {
    // this test is similar in spirit to section 4.2 of Holmes and Featherstone paper,
    // but reduced to lower degree since our reference implementation is MUCH slower
    // than the one used in the paper (Clenshaw method)
    int max = 50;
    NormalizedSphericalHarmonicsProvider provider = new GleasonProvider(max, max);
    HolmesFeatherstoneAttractionModel model = new HolmesFeatherstoneAttractionModel(itrf, provider);
    Assert.assertTrue(model.dependsOnPositionOnly());
    // Note that despite it uses adjustable high accuracy, the reference model
    // uses unstable formulas and hence loses lots of digits near poles.
    // This implies that if we still want about 16 digits near the poles, we
    // need to ask for at least 30 digits in computation. Setting for example
    // the following to 28 digits makes the test fail as the relative errors
    // raise to about 10^-12 near North pole and near 10^-11 near South pole.
    // The reason for this is that the reference becomes less accurate than
    // the model we are testing!
    int digits = 30;
    ReferenceFieldModel refModel = new ReferenceFieldModel(provider, digits);
    double r = 1.25;
    for (double theta = 0.01; theta < 3.14; theta += 0.1) {
        Vector3D position = new Vector3D(r * FastMath.sin(theta), 0.0, r * FastMath.cos(theta));
        Dfp refValue = refModel.nonCentralPart(null, position);
        double value = model.nonCentralPart(null, position, model.getMu());
        double relativeError = error(refValue, value).divide(refValue).toDouble();
        Assert.assertEquals(0, relativeError, 7.0e-15);
    }
}

Example 2 with Dfp

use of org.hipparchus.dfp.Dfp in project Orekit by CS-SI.

the class ReferenceFieldModel method nonCentralPart.

public Dfp nonCentralPart(final AbsoluteDate date, final Vector3D position) throws OrekitException {
    int degree = provider.getMaxDegree();
    int order = provider.getMaxOrder();
    // use coefficients without caring if they are the correct type
    final RawSphericalHarmonics harmonics = raw(provider).onDate(date);
    Dfp x = dfpField.newDfp(position.getX());
    Dfp y = dfpField.newDfp(position.getY());
    Dfp z = dfpField.newDfp(position.getZ());
    Dfp rho2 = x.multiply(x).add(y.multiply(y));
    Dfp rho = rho2.sqrt();
    Dfp r2 = rho2.add(z.multiply(z));
    Dfp r = r2.sqrt();
    Dfp aOr = dfpField.newDfp(provider.getAe()).divide(r);
    Dfp lambda = position.getX() > 0 ? dfpField.getTwo().multiply(DfpMath.atan(y.divide(rho.add(x)))) : dfpField.getPi().subtract(dfpField.getTwo().multiply(DfpMath.atan(y.divide(rho.subtract(x)))));
    Dfp cosTheta = z.divide(r);
    Dfp value = dfpField.getZero();
    Dfp aOrN = aOr;
    for (int n = 2; n <= degree; ++n) {
        Dfp sum = dfpField.getZero();
        for (int m = 0; m <= FastMath.min(n, order); ++m) {
            double cnm = harmonics.getRawCnm(n, m);
            double snm = harmonics.getRawSnm(n, m);
            Dfp mLambda = lambda.multiply(m);
            Dfp c = DfpMath.cos(mLambda).multiply(dfpField.newDfp(cnm));
            Dfp s = DfpMath.sin(mLambda).multiply(dfpField.newDfp(snm));
            Dfp pnm = alf[n][m].value(cosTheta);
            sum = sum.add(pnm.multiply(c.add(s)));
        }
        aOrN = aOrN.multiply(aOr);
        value = value.add(aOrN.multiply(sum));
    }
    return value.multiply(dfpField.newDfp(provider.getMu())).divide(r);
}

Example 3 with Dfp

use of org.hipparchus.dfp.Dfp in project Orekit by CS-SI.

the class AssociatedLegendreFunction method value.

public Dfp value(Dfp t) {
    Dfp y1 = polynomial[polynomial.length - 1];
    for (int j = polynomial.length - 2; j >= 0; j--) {
        y1 = y1.multiply(t).add(polynomial[j]);
    }
    Dfp oneMinusT2 = t.getField().getOne().subtract(t.multiply(t));
    Dfp y2 = t.getField().getOne();
    for (int j = 0; j < m; ++j) {
        y2 = y2.multiply(oneMinusT2);
    }
    y2 = y2.sqrt();
    return y1.multiply(y2).multiply(normalization);
}

Example 4 with Dfp

use of org.hipparchus.dfp.Dfp in project Orekit by CS-SI.

the class AssociatedLegendreFunction method getLegendrePolynomial.

private Dfp[] getLegendrePolynomial(int n, DfpField dfpField) {
    // get (or create) the list of polynomials for the specified field
    List<Dfp[]> list = LEGENDRE_POLYNOMIALS.get(dfpField.getRadixDigits());
    if (list == null) {
        list = new ArrayList<Dfp[]>();
        list.add(new Dfp[] { // P0(X) = 1
        dfpField.getOne() });
        list.add(new Dfp[] { // P1(X) = 0 + 1 * X
        dfpField.getZero(), // P1(X) = 0 + 1 * X
        dfpField.getOne() });
    }
    while (list.size() <= n) {
        // build polynomial Pk+1 using recursion formula
        // (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
        int k = list.size() - 1;
        Dfp kDfp = dfpField.newDfp(k);
        Dfp kP1Dfp = kDfp.add(dfpField.getOne());
        Dfp ckDfp = kP1Dfp.add(kDfp).divide(kP1Dfp);
        Dfp[] pk = list.get(k);
        Dfp ckM1Dfp = kDfp.divide(kP1Dfp).negate();
        Dfp[] pkM1 = list.get(k - 1);
        Dfp[] pkP1 = new Dfp[k + 2];
        pkP1[0] = ckM1Dfp.multiply(pkM1[0]);
        for (int i = 0; i < k; ++i) {
            if ((k - i) % 2 == 1) {
                pkP1[i + 1] = dfpField.getZero();
            } else {
                pkP1[i + 1] = ckDfp.multiply(pk[i]).add(ckM1Dfp.multiply(pkM1[i + 1]));
            }
        }
        pkP1[k + 1] = ckDfp.multiply(pk[k]);
        list.add(pkP1);
    }
    // retrieve degree n polynomial
    return list.get(n);
}