Examples with Well19937a org.hipparchus.random.Well19937a used on opensource projects

Example 61 with Well19937a

use of org.hipparchus.random.Well19937a in project Orekit by CS-SI.

the class FieldTransformTest method doTestTranslationDouble.

private <T extends RealFieldElement<T>> void doTestTranslationDouble(Field<T> field) {
    RandomGenerator rnd = new Well19937a(0x7e9d737ba4147787l);
    for (int i = 0; i < 10; ++i) {
        FieldVector3D<T> delta = randomVector(field, 1.0e3, rnd);
        FieldTransform<T> transform = new FieldTransform<>(FieldAbsoluteDate.getJ2000Epoch(field), delta);
        for (int j = 0; j < 10; ++j) {
            Vector3D a = createVector(field, rnd.nextDouble(), rnd.nextDouble(), rnd.nextDouble()).toVector3D();
            FieldVector3D<T> b = transform.transformVector(a);
            Assert.assertEquals(0, b.subtract(a).getNorm().getReal(), 1.0e-15);
            FieldVector3D<T> c = transform.transformPosition(a);
            Assert.assertEquals(0, c.subtract(a).subtract(delta).getNorm().getReal(), 1.0e-14);
        }
    }
}

Example 62 with Well19937a

use of org.hipparchus.random.Well19937a in project Orekit by CS-SI.

the class FieldTransformTest method doTestLine.

private <T extends RealFieldElement<T>> void doTestLine(Field<T> field) {
    RandomGenerator random = new Well19937a(0x4a5ff67426c5731fl);
    for (int i = 0; i < 100; ++i) {
        FieldTransform<T> transform = randomTransform(field, random);
        for (int j = 0; j < 20; ++j) {
            FieldVector3D<T> p0 = randomVector(field, 1.0e3, random);
            FieldVector3D<T> p1 = randomVector(field, 1.0e3, random);
            FieldLine<T> l = new FieldLine<>(p0, p1, 1.0e-10);
            FieldLine<T> transformed = transform.transformLine(l);
            for (int k = 0; k < 10; ++k) {
                FieldVector3D<T> p = l.pointAt(random.nextDouble() * 1.0e6);
                Assert.assertEquals(0.0, transformed.distance(transform.transformPosition(p)).getReal(), 1.0e-9);
            }
        }
    }
}

Example 63 with Well19937a

use of org.hipparchus.random.Well19937a in project Orekit by CS-SI.

the class FieldTransformTest method doTestJacobianPVA.

private <T extends RealFieldElement<T>> void doTestJacobianPVA(Field<T> field) {
    // base directions for finite differences
    @SuppressWarnings("unchecked") FieldPVCoordinates<T>[] directions = (FieldPVCoordinates<T>[]) Array.newInstance(FieldPVCoordinates.class, 9);
    directions[0] = new FieldPVCoordinates<>(FieldVector3D.getPlusI(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
    directions[1] = new FieldPVCoordinates<>(FieldVector3D.getPlusJ(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
    directions[2] = new FieldPVCoordinates<>(FieldVector3D.getPlusK(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
    directions[3] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getPlusI(field), FieldVector3D.getZero(field));
    directions[4] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getPlusJ(field), FieldVector3D.getZero(field));
    directions[5] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getPlusK(field), FieldVector3D.getZero(field));
    directions[6] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getZero(field), FieldVector3D.getPlusI(field));
    directions[7] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getZero(field), FieldVector3D.getPlusJ(field));
    directions[8] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getZero(field), FieldVector3D.getPlusK(field));
    double h = 0.01;
    RandomGenerator random = new Well19937a(0xd223e88b6232198fl);
    for (int i = 0; i < 20; ++i) {
        // generate a random transform
        FieldTransform<T> combined = randomTransform(field, random);
        // compute Jacobian
        T[][] jacobian = MathArrays.buildArray(field, 9, 9);
        for (int l = 0; l < jacobian.length; ++l) {
            for (int c = 0; c < jacobian[l].length; ++c) {
                jacobian[l][c] = field.getZero().add(1 + 0.1 * c);
            }
        }
        combined.getJacobian(CartesianDerivativesFilter.USE_PVA, jacobian);
        for (int j = 0; j < 100; ++j) {
            FieldPVCoordinates<T> pv0 = new FieldPVCoordinates<>(randomVector(field, 1e3, random), randomVector(field, 1.0, random), randomVector(field, 1.0e-3, random));
            double epsilonP = 2.0e-12 * pv0.getPosition().getNorm().getReal();
            double epsilonV = 6.0e-11 * pv0.getVelocity().getNorm().getReal();
            double epsilonA = 2.0e-9 * pv0.getAcceleration().getNorm().getReal();
            for (int c = 0; c < directions.length; ++c) {
                // eight points finite differences estimation of a Jacobian column
                FieldPVCoordinates<T> pvm4h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -4 * h, directions[c]));
                FieldPVCoordinates<T> pvm3h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -3 * h, directions[c]));
                FieldPVCoordinates<T> pvm2h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -2 * h, directions[c]));
                FieldPVCoordinates<T> pvm1h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -1 * h, directions[c]));
                FieldPVCoordinates<T> pvp1h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +1 * h, directions[c]));
                FieldPVCoordinates<T> pvp2h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +2 * h, directions[c]));
                FieldPVCoordinates<T> pvp3h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +3 * h, directions[c]));
                FieldPVCoordinates<T> pvp4h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +4 * h, directions[c]));
                FieldPVCoordinates<T> d4 = new FieldPVCoordinates<>(pvm4h, pvp4h);
                FieldPVCoordinates<T> d3 = new FieldPVCoordinates<>(pvm3h, pvp3h);
                FieldPVCoordinates<T> d2 = new FieldPVCoordinates<>(pvm2h, pvp2h);
                FieldPVCoordinates<T> d1 = new FieldPVCoordinates<>(pvm1h, pvp1h);
                double d = 1.0 / (840 * h);
                FieldPVCoordinates<T> estimatedColumn = new FieldPVCoordinates<>(-3 * d, d4, 32 * d, d3, -168 * d, d2, 672 * d, d1);
                // check analytical Jacobian against finite difference reference
                Assert.assertEquals(estimatedColumn.getPosition().getX().getReal(), jacobian[0][c].getReal(), epsilonP);
                Assert.assertEquals(estimatedColumn.getPosition().getY().getReal(), jacobian[1][c].getReal(), epsilonP);
                Assert.assertEquals(estimatedColumn.getPosition().getZ().getReal(), jacobian[2][c].getReal(), epsilonP);
                Assert.assertEquals(estimatedColumn.getVelocity().getX().getReal(), jacobian[3][c].getReal(), epsilonV);
                Assert.assertEquals(estimatedColumn.getVelocity().getY().getReal(), jacobian[4][c].getReal(), epsilonV);
                Assert.assertEquals(estimatedColumn.getVelocity().getZ().getReal(), jacobian[5][c].getReal(), epsilonV);
                Assert.assertEquals(estimatedColumn.getAcceleration().getX().getReal(), jacobian[6][c].getReal(), epsilonA);
                Assert.assertEquals(estimatedColumn.getAcceleration().getY().getReal(), jacobian[7][c].getReal(), epsilonA);
                Assert.assertEquals(estimatedColumn.getAcceleration().getZ().getReal(), jacobian[8][c].getReal(), epsilonA);
            }
        }
    }
}

Example 64 with Well19937a

use of org.hipparchus.random.Well19937a in project Orekit by CS-SI.

the class FieldTransformTest method doTestShiftDerivatives.

private <T extends RealFieldElement<T>> void doTestShiftDerivatives(Field<T> field) {
    RandomGenerator random = new Well19937a(0x5acda4f605aadce7l);
    for (int i = 0; i < 10; ++i) {
        FieldTransform<T> t = randomTransform(field, random);
        for (double dtD = -10.0; dtD < 10.0; dtD += 0.125) {
            T dt = field.getZero().add(dtD);
            FieldTransform<T> t0 = t.shiftedBy(dt);
            T v = t0.getVelocity().getNorm();
            T a = t0.getAcceleration().getNorm();
            T omega = t0.getRotationRate().getNorm();
            T omegaDot = t0.getRotationAcceleration().getNorm();
            // numerical derivatives
            T h = omega.reciprocal().multiply(0.01);
            FieldTransform<T> tm4h = t.shiftedBy(dt.subtract(h.multiply(4)));
            FieldTransform<T> tm3h = t.shiftedBy(dt.subtract(h.multiply(3)));
            FieldTransform<T> tm2h = t.shiftedBy(dt.subtract(h.multiply(2)));
            FieldTransform<T> tm1h = t.shiftedBy(dt.subtract(h.multiply(1)));
            FieldTransform<T> tp1h = t.shiftedBy(dt.add(h.multiply(1)));
            FieldTransform<T> tp2h = t.shiftedBy(dt.add(h.multiply(2)));
            FieldTransform<T> tp3h = t.shiftedBy(dt.add(h.multiply(3)));
            FieldTransform<T> tp4h = t.shiftedBy(dt.add(h.multiply(4)));
            T numXDot = derivative(h, tm4h.getTranslation().getX(), tm3h.getTranslation().getX(), tm2h.getTranslation().getX(), tm1h.getTranslation().getX(), tp1h.getTranslation().getX(), tp2h.getTranslation().getX(), tp3h.getTranslation().getX(), tp4h.getTranslation().getX());
            T numYDot = derivative(h, tm4h.getTranslation().getY(), tm3h.getTranslation().getY(), tm2h.getTranslation().getY(), tm1h.getTranslation().getY(), tp1h.getTranslation().getY(), tp2h.getTranslation().getY(), tp3h.getTranslation().getY(), tp4h.getTranslation().getY());
            T numZDot = derivative(h, tm4h.getTranslation().getZ(), tm3h.getTranslation().getZ(), tm2h.getTranslation().getZ(), tm1h.getTranslation().getZ(), tp1h.getTranslation().getZ(), tp2h.getTranslation().getZ(), tp3h.getTranslation().getZ(), tp4h.getTranslation().getZ());
            T numXDot2 = derivative(h, tm4h.getVelocity().getX(), tm3h.getVelocity().getX(), tm2h.getVelocity().getX(), tm1h.getVelocity().getX(), tp1h.getVelocity().getX(), tp2h.getVelocity().getX(), tp3h.getVelocity().getX(), tp4h.getVelocity().getX());
            T numYDot2 = derivative(h, tm4h.getVelocity().getY(), tm3h.getVelocity().getY(), tm2h.getVelocity().getY(), tm1h.getVelocity().getY(), tp1h.getVelocity().getY(), tp2h.getVelocity().getY(), tp3h.getVelocity().getY(), tp4h.getVelocity().getY());
            T numZDot2 = derivative(h, tm4h.getVelocity().getZ(), tm3h.getVelocity().getZ(), tm2h.getVelocity().getZ(), tm1h.getVelocity().getZ(), tp1h.getVelocity().getZ(), tp2h.getVelocity().getZ(), tp3h.getVelocity().getZ(), tp4h.getVelocity().getZ());
            T numQ0Dot = derivative(h, tm4h.getRotation().getQ0(), tm3h.getRotation().getQ0(), tm2h.getRotation().getQ0(), tm1h.getRotation().getQ0(), tp1h.getRotation().getQ0(), tp2h.getRotation().getQ0(), tp3h.getRotation().getQ0(), tp4h.getRotation().getQ0());
            T numQ1Dot = derivative(h, tm4h.getRotation().getQ1(), tm3h.getRotation().getQ1(), tm2h.getRotation().getQ1(), tm1h.getRotation().getQ1(), tp1h.getRotation().getQ1(), tp2h.getRotation().getQ1(), tp3h.getRotation().getQ1(), tp4h.getRotation().getQ1());
            T numQ2Dot = derivative(h, tm4h.getRotation().getQ2(), tm3h.getRotation().getQ2(), tm2h.getRotation().getQ2(), tm1h.getRotation().getQ2(), tp1h.getRotation().getQ2(), tp2h.getRotation().getQ2(), tp3h.getRotation().getQ2(), tp4h.getRotation().getQ2());
            T numQ3Dot = derivative(h, tm4h.getRotation().getQ3(), tm3h.getRotation().getQ3(), tm2h.getRotation().getQ3(), tm1h.getRotation().getQ3(), tp1h.getRotation().getQ3(), tp2h.getRotation().getQ3(), tp3h.getRotation().getQ3(), tp4h.getRotation().getQ3());
            T numOxDot = derivative(h, tm4h.getRotationRate().getX(), tm3h.getRotationRate().getX(), tm2h.getRotationRate().getX(), tm1h.getRotationRate().getX(), tp1h.getRotationRate().getX(), tp2h.getRotationRate().getX(), tp3h.getRotationRate().getX(), tp4h.getRotationRate().getX());
            T numOyDot = derivative(h, tm4h.getRotationRate().getY(), tm3h.getRotationRate().getY(), tm2h.getRotationRate().getY(), tm1h.getRotationRate().getY(), tp1h.getRotationRate().getY(), tp2h.getRotationRate().getY(), tp3h.getRotationRate().getY(), tp4h.getRotationRate().getY());
            T numOzDot = derivative(h, tm4h.getRotationRate().getZ(), tm3h.getRotationRate().getZ(), tm2h.getRotationRate().getZ(), tm1h.getRotationRate().getZ(), tp1h.getRotationRate().getZ(), tp2h.getRotationRate().getZ(), tp3h.getRotationRate().getZ(), tp4h.getRotationRate().getZ());
            // theoretical derivatives
            T theXDot = t0.getVelocity().getX();
            T theYDot = t0.getVelocity().getY();
            T theZDot = t0.getVelocity().getZ();
            T theXDot2 = t0.getAcceleration().getX();
            T theYDot2 = t0.getAcceleration().getY();
            T theZDot2 = t0.getAcceleration().getZ();
            FieldRotation<T> r0 = t0.getRotation();
            FieldVector3D<T> w = t0.getRotationRate();
            FieldVector3D<T> q = new FieldVector3D<>(r0.getQ1(), r0.getQ2(), r0.getQ3());
            FieldVector3D<T> qw = FieldVector3D.crossProduct(q, w);
            T theQ0Dot = FieldVector3D.dotProduct(q, w).multiply(-0.5);
            T theQ1Dot = r0.getQ0().multiply(w.getX()).add(qw.getX()).multiply(0.5);
            T theQ2Dot = r0.getQ0().multiply(w.getY()).add(qw.getY()).multiply(0.5);
            T theQ3Dot = r0.getQ0().multiply(w.getZ()).add(qw.getZ()).multiply(0.5);
            T theOxDot2 = t0.getRotationAcceleration().getX();
            T theOyDot2 = t0.getRotationAcceleration().getY();
            T theOzDot2 = t0.getRotationAcceleration().getZ();
            // check consistency
            Assert.assertEquals(theXDot.getReal(), numXDot.getReal(), 1.0e-13 * v.getReal());
            Assert.assertEquals(theYDot.getReal(), numYDot.getReal(), 1.0e-13 * v.getReal());
            Assert.assertEquals(theZDot.getReal(), numZDot.getReal(), 1.0e-13 * v.getReal());
            Assert.assertEquals(theXDot2.getReal(), numXDot2.getReal(), 1.0e-13 * a.getReal());
            Assert.assertEquals(theYDot2.getReal(), numYDot2.getReal(), 1.0e-13 * a.getReal());
            Assert.assertEquals(theZDot2.getReal(), numZDot2.getReal(), 1.0e-13 * a.getReal());
            Assert.assertEquals(theQ0Dot.getReal(), numQ0Dot.getReal(), 1.0e-13 * omega.getReal());
            Assert.assertEquals(theQ1Dot.getReal(), numQ1Dot.getReal(), 1.0e-13 * omega.getReal());
            Assert.assertEquals(theQ2Dot.getReal(), numQ2Dot.getReal(), 1.0e-13 * omega.getReal());
            Assert.assertEquals(theQ3Dot.getReal(), numQ3Dot.getReal(), 1.0e-13 * omega.getReal());
            Assert.assertEquals(theOxDot2.getReal(), numOxDot.getReal(), 1.0e-12 * omegaDot.getReal());
            Assert.assertEquals(theOyDot2.getReal(), numOyDot.getReal(), 1.0e-12 * omegaDot.getReal());
            Assert.assertEquals(theOzDot2.getReal(), numOzDot.getReal(), 1.0e-12 * omegaDot.getReal());
        }
    }
}

Example 65 with Well19937a

use of org.hipparchus.random.Well19937a in project Orekit by CS-SI.

the class FieldTransformTest method doTestTransPVDouble.

private <T extends RealFieldElement<T>> void doTestTransPVDouble(Field<T> field) {
    RandomGenerator rnd = new Well19937a(0x73d5554d99427af0l);
    for (int i = 0; i < 10; ++i) {
        // random position, velocity and acceleration
        Vector3D pos = randomVector(field, 1.0e3, rnd).toVector3D();
        Vector3D vel = randomVector(field, 1.0, rnd).toVector3D();
        Vector3D acc = randomVector(field, 1.0e-3, rnd).toVector3D();
        PVCoordinates pvOne = new PVCoordinates(pos, vel, acc);
        // random transform
        FieldVector3D<T> transPos = randomVector(field, 1.0e3, rnd);
        FieldVector3D<T> transVel = randomVector(field, 1.0, rnd);
        FieldVector3D<T> transAcc = randomVector(field, 1.0e-3, rnd);
        FieldTransform<T> tr = new FieldTransform<>(FieldAbsoluteDate.getJ2000Epoch(field), transPos, transVel, transAcc);
        T dt = field.getZero().add(1);
        // we should obtain
        FieldVector3D<T> good = tr.transformPosition(new FieldVector3D<>(dt, vel).add(pos)).add(new FieldVector3D<>(dt, transVel));
        // we have
        FieldPVCoordinates<T> pvTwo = tr.transformPVCoordinates(pvOne);
        FieldVector3D<T> result = pvTwo.getPosition().add(new FieldVector3D<>(dt, pvTwo.getVelocity()));
        checkVector(good, result, 1.0e-15);
        // test inverse
        FieldVector3D<T> resultvel = tr.getInverse().transformPVCoordinates(pvTwo).getVelocity();
        checkVector(resultvel, new FieldVector3D<>(field, vel), 1.0e-15);
    }
}