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

Example 16 with Well19937a

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

the class TransformTest method testShiftDerivatives.

@Test
public void testShiftDerivatives() {
    RandomGenerator random = new Well19937a(0x5acda4f605aadce7l);
    for (int i = 0; i < 10; ++i) {
        Transform t = randomTransform(random);
        for (double dt = -10.0; dt < 10.0; dt += 0.125) {
            Transform t0 = t.shiftedBy(dt);
            double v = t0.getVelocity().getNorm();
            double a = t0.getAcceleration().getNorm();
            double omega = t0.getRotationRate().getNorm();
            double omegaDot = t0.getRotationAcceleration().getNorm();
            // numerical derivatives
            double h = 0.01 / omega;
            Transform tm4h = t.shiftedBy(dt - 4 * h);
            Transform tm3h = t.shiftedBy(dt - 3 * h);
            Transform tm2h = t.shiftedBy(dt - 2 * h);
            Transform tm1h = t.shiftedBy(dt - 1 * h);
            Transform tp1h = t.shiftedBy(dt + 1 * h);
            Transform tp2h = t.shiftedBy(dt + 2 * h);
            Transform tp3h = t.shiftedBy(dt + 3 * h);
            Transform tp4h = t.shiftedBy(dt + 4 * h);
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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());
            double 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
            double theXDot = t0.getVelocity().getX();
            double theYDot = t0.getVelocity().getY();
            double theZDot = t0.getVelocity().getZ();
            double theXDot2 = t0.getAcceleration().getX();
            double theYDot2 = t0.getAcceleration().getY();
            double theZDot2 = t0.getAcceleration().getZ();
            Rotation r0 = t0.getRotation();
            Vector3D w = t0.getRotationRate();
            Vector3D q = new Vector3D(r0.getQ1(), r0.getQ2(), r0.getQ3());
            Vector3D qw = Vector3D.crossProduct(q, w);
            double theQ0Dot = -0.5 * Vector3D.dotProduct(q, w);
            double theQ1Dot = 0.5 * (r0.getQ0() * w.getX() + qw.getX());
            double theQ2Dot = 0.5 * (r0.getQ0() * w.getY() + qw.getY());
            double theQ3Dot = 0.5 * (r0.getQ0() * w.getZ() + qw.getZ());
            double theOxDot2 = t0.getRotationAcceleration().getX();
            double theOyDot2 = t0.getRotationAcceleration().getY();
            double theOzDot2 = t0.getRotationAcceleration().getZ();
            // check consistency
            Assert.assertEquals(theXDot, numXDot, 1.0e-13 * v);
            Assert.assertEquals(theYDot, numYDot, 1.0e-13 * v);
            Assert.assertEquals(theZDot, numZDot, 1.0e-13 * v);
            Assert.assertEquals(theXDot2, numXDot2, 1.0e-13 * a);
            Assert.assertEquals(theYDot2, numYDot2, 1.0e-13 * a);
            Assert.assertEquals(theZDot2, numZDot2, 1.0e-13 * a);
            Assert.assertEquals(theQ0Dot, numQ0Dot, 1.0e-13 * omega);
            Assert.assertEquals(theQ1Dot, numQ1Dot, 1.0e-13 * omega);
            Assert.assertEquals(theQ2Dot, numQ2Dot, 1.0e-13 * omega);
            Assert.assertEquals(theQ3Dot, numQ3Dot, 1.0e-13 * omega);
            Assert.assertEquals(theOxDot2, numOxDot, 1.0e-12 * omegaDot);
            Assert.assertEquals(theOyDot2, numOyDot, 1.0e-12 * omegaDot);
            Assert.assertEquals(theOzDot2, numOzDot, 1.0e-12 * omegaDot);
        }
    }
}

Example 17 with Well19937a

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

the class TransformTest method testTranslation.

@Test
public void testTranslation() {
    RandomGenerator rnd = new Well19937a(0x7e9d737ba4147787l);
    for (int i = 0; i < 10; ++i) {
        Vector3D delta = randomVector(1.0e3, rnd);
        Transform transform = new Transform(AbsoluteDate.J2000_EPOCH, delta);
        for (int j = 0; j < 10; ++j) {
            Vector3D a = new Vector3D(rnd.nextDouble(), rnd.nextDouble(), rnd.nextDouble());
            Vector3D b = transform.transformVector(a);
            Assert.assertEquals(0, b.subtract(a).getNorm(), 1.0e-15);
            Vector3D c = transform.transformPosition(a);
            Assert.assertEquals(0, c.subtract(a).subtract(delta).getNorm(), 1.0e-14);
        }
    }
}

Example 18 with Well19937a

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

the class TransformTest method testJacobianPVA.

@Test
public void testJacobianPVA() {
    // base directions for finite differences
    PVCoordinates[] directions = new PVCoordinates[] { new PVCoordinates(Vector3D.PLUS_I, Vector3D.ZERO, Vector3D.ZERO), new PVCoordinates(Vector3D.PLUS_J, Vector3D.ZERO, Vector3D.ZERO), new PVCoordinates(Vector3D.PLUS_K, Vector3D.ZERO, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.PLUS_I, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.PLUS_J, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.PLUS_K, Vector3D.ZERO), new PVCoordinates(Vector3D.ZERO, Vector3D.ZERO, Vector3D.PLUS_I), new PVCoordinates(Vector3D.ZERO, Vector3D.ZERO, Vector3D.PLUS_J), new PVCoordinates(Vector3D.ZERO, Vector3D.ZERO, Vector3D.PLUS_K) };
    double h = 0.01;
    RandomGenerator random = new Well19937a(0xd223e88b6232198fl);
    for (int i = 0; i < 20; ++i) {
        // generate a random transform
        Transform combined = randomTransform(random);
        // compute Jacobian
        double[][] jacobian = new double[9][9];
        for (int l = 0; l < jacobian.length; ++l) {
            for (int c = 0; c < jacobian[l].length; ++c) {
                jacobian[l][c] = l + 0.1 * c;
            }
        }
        combined.getJacobian(CartesianDerivativesFilter.USE_PVA, jacobian);
        for (int j = 0; j < 100; ++j) {
            PVCoordinates pv0 = new PVCoordinates(randomVector(1e3, random), randomVector(1.0, random), randomVector(1.0e-3, random));
            double epsilonP = 2.0e-12 * pv0.getPosition().getNorm();
            double epsilonV = 6.0e-11 * pv0.getVelocity().getNorm();
            double epsilonA = 2.0e-9 * pv0.getAcceleration().getNorm();
            for (int c = 0; c < directions.length; ++c) {
                // eight points finite differences estimation of a Jacobian column
                PVCoordinates pvm4h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, -4 * h, directions[c]));
                PVCoordinates pvm3h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, -3 * h, directions[c]));
                PVCoordinates pvm2h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, -2 * h, directions[c]));
                PVCoordinates pvm1h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, -1 * h, directions[c]));
                PVCoordinates pvp1h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, +1 * h, directions[c]));
                PVCoordinates pvp2h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, +2 * h, directions[c]));
                PVCoordinates pvp3h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, +3 * h, directions[c]));
                PVCoordinates pvp4h = combined.transformPVCoordinates(new PVCoordinates(1.0, pv0, +4 * h, directions[c]));
                PVCoordinates d4 = new PVCoordinates(pvm4h, pvp4h);
                PVCoordinates d3 = new PVCoordinates(pvm3h, pvp3h);
                PVCoordinates d2 = new PVCoordinates(pvm2h, pvp2h);
                PVCoordinates d1 = new PVCoordinates(pvm1h, pvp1h);
                double d = 1.0 / (840 * h);
                PVCoordinates estimatedColumn = new PVCoordinates(-3 * d, d4, 32 * d, d3, -168 * d, d2, 672 * d, d1);
                // check analytical Jacobian against finite difference reference
                Assert.assertEquals(estimatedColumn.getPosition().getX(), jacobian[0][c], epsilonP);
                Assert.assertEquals(estimatedColumn.getPosition().getY(), jacobian[1][c], epsilonP);
                Assert.assertEquals(estimatedColumn.getPosition().getZ(), jacobian[2][c], epsilonP);
                Assert.assertEquals(estimatedColumn.getVelocity().getX(), jacobian[3][c], epsilonV);
                Assert.assertEquals(estimatedColumn.getVelocity().getY(), jacobian[4][c], epsilonV);
                Assert.assertEquals(estimatedColumn.getVelocity().getZ(), jacobian[5][c], epsilonV);
                Assert.assertEquals(estimatedColumn.getAcceleration().getX(), jacobian[6][c], epsilonA);
                Assert.assertEquals(estimatedColumn.getAcceleration().getY(), jacobian[7][c], epsilonA);
                Assert.assertEquals(estimatedColumn.getAcceleration().getZ(), jacobian[8][c], epsilonA);
            }
        }
    }
}

Example 19 with Well19937a

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

the class TransformTest method testReverse.

@Test
public void testReverse() {
    RandomGenerator random = new Well19937a(0x9f82ba2b2c98dac5l);
    for (int i = 0; i < 20; ++i) {
        Transform combined = randomTransform(random);
        checkNoTransform(new Transform(AbsoluteDate.J2000_EPOCH, combined, combined.getInverse()), random);
    }
}

Example 20 with Well19937a

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

the class TransformTest method testLinear.

@Test
public void testLinear() {
    RandomGenerator random = new Well19937a(0x14f6411217b148d8l);
    for (int n = 0; n < 100; ++n) {
        Transform t = randomTransform(random);
        // build an equivalent linear transform by extracting raw translation/rotation
        RealMatrix linearA = MatrixUtils.createRealMatrix(3, 4);
        linearA.setSubMatrix(t.getRotation().getMatrix(), 0, 0);
        Vector3D rt = t.getRotation().applyTo(t.getTranslation());
        linearA.setEntry(0, 3, rt.getX());
        linearA.setEntry(1, 3, rt.getY());
        linearA.setEntry(2, 3, rt.getZ());
        // build an equivalent linear transform by observing transformed points
        RealMatrix linearB = MatrixUtils.createRealMatrix(3, 4);
        Vector3D p0 = t.transformPosition(Vector3D.ZERO);
        Vector3D pI = t.transformPosition(Vector3D.PLUS_I).subtract(p0);
        Vector3D pJ = t.transformPosition(Vector3D.PLUS_J).subtract(p0);
        Vector3D pK = t.transformPosition(Vector3D.PLUS_K).subtract(p0);
        linearB.setColumn(0, new double[] { pI.getX(), pI.getY(), pI.getZ() });
        linearB.setColumn(1, new double[] { pJ.getX(), pJ.getY(), pJ.getZ() });
        linearB.setColumn(2, new double[] { pK.getX(), pK.getY(), pK.getZ() });
        linearB.setColumn(3, new double[] { p0.getX(), p0.getY(), p0.getZ() });
        // both linear transforms should be equal
        Assert.assertEquals(0.0, linearB.subtract(linearA).getNorm(), 1.0e-15 * linearA.getNorm());
        for (int i = 0; i < 100; ++i) {
            Vector3D p = randomVector(1.0e3, random);
            Vector3D q = t.transformPosition(p);
            double[] qA = linearA.operate(new double[] { p.getX(), p.getY(), p.getZ(), 1.0 });
            Assert.assertEquals(q.getX(), qA[0], 1.0e-13 * p.getNorm());
            Assert.assertEquals(q.getY(), qA[1], 1.0e-13 * p.getNorm());
            Assert.assertEquals(q.getZ(), qA[2], 1.0e-13 * p.getNorm());
            double[] qB = linearB.operate(new double[] { p.getX(), p.getY(), p.getZ(), 1.0 });
            Assert.assertEquals(q.getX(), qB[0], 1.0e-10 * p.getNorm());
            Assert.assertEquals(q.getY(), qB[1], 1.0e-10 * p.getNorm());
            Assert.assertEquals(q.getZ(), qB[2], 1.0e-10 * p.getNorm());
        }
    }
}