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

Example 31 with Well19937a

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

the class FieldTransformTest method doTestLinear.

private <T extends RealFieldElement<T>> void doTestLinear(final Field<T> field) {
    RandomGenerator random = new Well19937a(0x14f6411217b148d8l);
    for (int n = 0; n < 100; ++n) {
        FieldTransform<T> t = randomTransform(field, random);
        // build an equivalent linear transform by extracting raw translation/rotation
        FieldMatrix<T> linearA = MatrixUtils.createFieldMatrix(field, 3, 4);
        linearA.setSubMatrix(t.getRotation().getMatrix(), 0, 0);
        FieldVector3D<T> 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
        FieldMatrix<T> linearB = MatrixUtils.createFieldMatrix(field, 3, 4);
        FieldVector3D<T> p0 = t.transformPosition(FieldVector3D.getZero(field));
        FieldVector3D<T> pI = t.transformPosition(FieldVector3D.getPlusI(field)).subtract(p0);
        FieldVector3D<T> pJ = t.transformPosition(FieldVector3D.getPlusJ(field)).subtract(p0);
        FieldVector3D<T> pK = t.transformPosition(FieldVector3D.getPlusK(field)).subtract(p0);
        linearB.setEntry(0, 0, pI.getX());
        linearB.setEntry(1, 0, pI.getY());
        linearB.setEntry(2, 0, pI.getZ());
        linearB.setEntry(0, 1, pJ.getX());
        linearB.setEntry(1, 1, pJ.getY());
        linearB.setEntry(2, 1, pJ.getZ());
        linearB.setEntry(0, 2, pK.getX());
        linearB.setEntry(1, 2, pK.getY());
        linearB.setEntry(2, 2, pK.getZ());
        linearB.setEntry(0, 3, p0.getX());
        linearB.setEntry(1, 3, p0.getY());
        linearB.setEntry(2, 3, p0.getZ());
        // both linear transforms should be equal
        FieldMatrix<T> sub = linearB.subtract(linearA);
        double refMax = 0;
        double diffMax = 0;
        for (int i = 0; i < linearA.getRowDimension(); ++i) {
            for (int j = 0; j < linearA.getColumnDimension(); ++j) {
                refMax = FastMath.max(linearA.getEntry(i, j).getReal(), refMax);
                diffMax = FastMath.max(sub.getEntry(i, j).getReal(), diffMax);
            }
        }
        Assert.assertEquals(0.0, diffMax, 2.0e-12 * refMax);
        for (int i = 0; i < 100; ++i) {
            FieldVector3D<T> p = randomVector(field, 1.0e3, random);
            FieldVector3D<T> q = t.transformPosition(p);
            T[] pField = MathArrays.buildArray(field, 4);
            pField[0] = p.getX();
            pField[1] = p.getY();
            pField[2] = p.getZ();
            pField[3] = field.getOne();
            T[] qA = linearA.operate(pField);
            Assert.assertEquals(q.getX().getReal(), qA[0].getReal(), 1.0e-13 * p.getNorm().getReal());
            Assert.assertEquals(q.getY().getReal(), qA[1].getReal(), 1.0e-13 * p.getNorm().getReal());
            Assert.assertEquals(q.getZ().getReal(), qA[2].getReal(), 1.0e-13 * p.getNorm().getReal());
            T[] qB = linearB.operate(pField);
            Assert.assertEquals(q.getX().getReal(), qB[0].getReal(), 1.0e-10 * p.getNorm().getReal());
            Assert.assertEquals(q.getY().getReal(), qB[1].getReal(), 1.0e-10 * p.getNorm().getReal());
            Assert.assertEquals(q.getZ().getReal(), qB[2].getReal(), 1.0e-10 * p.getNorm().getReal());
        }
    }
}

Example 32 with Well19937a

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

the class FieldTransformTest method doTestTransPV.

private <T extends RealFieldElement<T>> void doTestTransPV(Field<T> field) {
    RandomGenerator rnd = new Well19937a(0x73d5554d99427af0l);
    for (int i = 0; i < 10; ++i) {
        // random position, velocity and acceleration
        FieldVector3D<T> pos = randomVector(field, 1.0e3, rnd);
        FieldVector3D<T> vel = randomVector(field, 1.0, rnd);
        FieldVector3D<T> acc = randomVector(field, 1.0e-3, rnd);
        FieldPVCoordinates<T> pvOne = new FieldPVCoordinates<>(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);
        double dt = 1;
        // we should obtain
        FieldVector3D<T> good = tr.transformPosition(pos.add(new FieldVector3D<>(dt, vel))).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, vel, 1.0e-15);
    }
}

Example 33 with Well19937a

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

the class FieldTransformTest method doTestJacobianP.

private <T extends RealFieldElement<T>> void doTestJacobianP(Field<T> field) {
    // base directions for finite differences
    @SuppressWarnings("unchecked") FieldPVCoordinates<T>[] directions = (FieldPVCoordinates<T>[]) Array.newInstance(FieldPVCoordinates.class, 3);
    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));
    double h = 0.01;
    RandomGenerator random = new Well19937a(0x47fd0d6809f4b173l);
    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(l + 0.1 * c);
            }
        }
        combined.getJacobian(CartesianDerivativesFilter.USE_P, 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();
            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);
                // check the rest of the matrix remains untouched
                for (int l = 3; l < jacobian.length; ++l) {
                    Assert.assertEquals(l + 0.1 * c, jacobian[l][c].getReal(), 1.0e-15);
                }
            }
            // check the rest of the matrix remains untouched
            for (int c = directions.length; c < jacobian[0].length; ++c) {
                for (int l = 0; l < jacobian.length; ++l) {
                    Assert.assertEquals(l + 0.1 * c, jacobian[l][c].getReal(), 1.0e-15);
                }
            }
        }
    }
}

Example 34 with Well19937a

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

the class FieldTransformTest method doTestIdentityLine.

private <T extends RealFieldElement<T>> void doTestIdentityLine(Field<T> field) {
    RandomGenerator random = new Well19937a(0x98603025df70db7cl);
    FieldVector3D<T> p1 = randomVector(field, 100.0, random);
    FieldVector3D<T> p2 = randomVector(field, 100.0, random);
    FieldLine<T> line = new FieldLine<>(p1, p2, 1.0e-6);
    FieldLine<T> transformed = FieldTransform.getIdentity(field).transformLine(line);
    Assert.assertSame(line, transformed);
}

Example 35 with Well19937a

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

the class GroundStationTest method doTestCartesianDerivatives.

private void doTestCartesianDerivatives(double latitude, double longitude, double altitude, double stepFactor, double relativeTolerancePositionValue, double relativeTolerancePositionDerivative, double relativeToleranceVelocityValue, double relativeToleranceVelocityDerivative, String... parameterPattern) throws OrekitException {
    Utils.setDataRoot("regular-data");
    final Frame eme2000 = FramesFactory.getEME2000();
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
    final AbsoluteDate date0 = date.shiftedBy(50000);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    final GroundStation station = new GroundStation(new TopocentricFrame(earth, new GeodeticPoint(latitude, longitude, altitude), "dummy"));
    final DSFactory factory = new DSFactory(parameterPattern.length, 1);
    final FieldAbsoluteDate<DerivativeStructure> dateDS = new FieldAbsoluteDate<>(factory.getDerivativeField(), date);
    ParameterDriver[] selectedDrivers = new ParameterDriver[parameterPattern.length];
    UnivariateDifferentiableVectorFunction[] dFCartesian = new UnivariateDifferentiableVectorFunction[parameterPattern.length];
    final ParameterDriver[] allDrivers = selectAllDrivers(station);
    for (ParameterDriver driver : allDrivers) {
        driver.setReferenceDate(date0);
    }
    Map<String, Integer> indices = new HashMap<>();
    for (int k = 0; k < dFCartesian.length; ++k) {
        for (int i = 0; i < allDrivers.length; ++i) {
            if (allDrivers[i].getName().matches(parameterPattern[k])) {
                selectedDrivers[k] = allDrivers[i];
                dFCartesian[k] = differentiatedStationPV(station, eme2000, date, selectedDrivers[k], stepFactor);
                indices.put(selectedDrivers[k].getName(), k);
            }
        }
    }
    ;
    DSFactory factory11 = new DSFactory(1, 1);
    RandomGenerator generator = new Well19937a(0x084d58a19c498a54l);
    double maxPositionValueRelativeError = 0;
    double maxPositionDerivativeRelativeError = 0;
    double maxVelocityValueRelativeError = 0;
    double maxVelocityDerivativeRelativeError = 0;
    for (int i = 0; i < 1000; ++i) {
        // randomly change one parameter
        ParameterDriver changed = allDrivers[generator.nextInt(allDrivers.length)];
        changed.setNormalizedValue(2 * generator.nextDouble() - 1);
        // transform to check
        FieldTransform<DerivativeStructure> t = station.getOffsetToInertial(eme2000, dateDS, factory, indices);
        FieldPVCoordinates<DerivativeStructure> pv = t.transformPVCoordinates(FieldPVCoordinates.getZero(factory.getDerivativeField()));
        for (int k = 0; k < dFCartesian.length; ++k) {
            // reference values and derivatives computed using finite differences
            DerivativeStructure[] refCartesian = dFCartesian[k].value(factory11.variable(0, selectedDrivers[k].getValue()));
            // position
            final Vector3D refP = new Vector3D(refCartesian[0].getValue(), refCartesian[1].getValue(), refCartesian[2].getValue());
            final Vector3D resP = new Vector3D(pv.getPosition().getX().getValue(), pv.getPosition().getY().getValue(), pv.getPosition().getZ().getValue());
            maxPositionValueRelativeError = FastMath.max(maxPositionValueRelativeError, Vector3D.distance(refP, resP) / refP.getNorm());
            final Vector3D refPD = new Vector3D(refCartesian[0].getPartialDerivative(1), refCartesian[1].getPartialDerivative(1), refCartesian[2].getPartialDerivative(1));
            final Vector3D resPD = new Vector3D(pv.getPosition().getX().getAllDerivatives()[k + 1], pv.getPosition().getY().getAllDerivatives()[k + 1], pv.getPosition().getZ().getAllDerivatives()[k + 1]);
            maxPositionDerivativeRelativeError = FastMath.max(maxPositionDerivativeRelativeError, Vector3D.distance(refPD, resPD) / refPD.getNorm());
            // velocity
            final Vector3D refV = new Vector3D(refCartesian[3].getValue(), refCartesian[4].getValue(), refCartesian[5].getValue());
            final Vector3D resV = new Vector3D(pv.getVelocity().getX().getValue(), pv.getVelocity().getY().getValue(), pv.getVelocity().getZ().getValue());
            maxVelocityValueRelativeError = FastMath.max(maxVelocityValueRelativeError, Vector3D.distance(refV, resV) / refV.getNorm());
            final Vector3D refVD = new Vector3D(refCartesian[3].getPartialDerivative(1), refCartesian[4].getPartialDerivative(1), refCartesian[5].getPartialDerivative(1));
            final Vector3D resVD = new Vector3D(pv.getVelocity().getX().getAllDerivatives()[k + 1], pv.getVelocity().getY().getAllDerivatives()[k + 1], pv.getVelocity().getZ().getAllDerivatives()[k + 1]);
            maxVelocityDerivativeRelativeError = FastMath.max(maxVelocityDerivativeRelativeError, Vector3D.distance(refVD, resVD) / refVD.getNorm());
        }
    }
    Assert.assertEquals(0.0, maxPositionValueRelativeError, relativeTolerancePositionValue);
    Assert.assertEquals(0.0, maxPositionDerivativeRelativeError, relativeTolerancePositionDerivative);
    Assert.assertEquals(0.0, maxVelocityValueRelativeError, relativeToleranceVelocityValue);
    Assert.assertEquals(0.0, maxVelocityDerivativeRelativeError, relativeToleranceVelocityDerivative);
}