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

Example 66 with RandomGenerator

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

the class JacobianPropagatorConverterTest method doTestDerivatives.

private void doTestDerivatives(double tolP, double tolV, String... names) throws OrekitException {
    // we use a fixed step integrator on purpose
    // as the test is based on external differentiation using finite differences,
    // an adaptive step size integrator would introduce *lots* of numerical noise
    NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit), new LutherIntegratorBuilder(10.0), PositionAngle.TRUE, dP);
    builder.setMass(200.0);
    builder.addForceModel(drag);
    builder.addForceModel(gravity);
    // retrieve a state slightly different from the initial state,
    // using normalized values different from 0.0 for the sake of generality
    RandomGenerator random = new Well19937a(0xe67f19c1a678d037l);
    List<ParameterDriver> all = new ArrayList<ParameterDriver>();
    for (final ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
        all.add(driver);
    }
    for (final ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
        all.add(driver);
    }
    double[] normalized = new double[names.length];
    List<ParameterDriver> selected = new ArrayList<ParameterDriver>(names.length);
    int index = 0;
    for (final ParameterDriver driver : all) {
        boolean found = false;
        for (final String name : names) {
            if (name.equals(driver.getName())) {
                found = true;
                normalized[index++] = driver.getNormalizedValue() + (2 * random.nextDouble() - 1);
                selected.add(driver);
            }
        }
        driver.setSelected(found);
    }
    // create a one hour sample that starts 10 minutes after initial state
    // the 10 minutes offset implies even the first point is influenced by model parameters
    final List<SpacecraftState> sample = new ArrayList<SpacecraftState>();
    Propagator propagator = builder.buildPropagator(normalized);
    propagator.setMasterMode(60.0, new OrekitFixedStepHandler() {

        @Override
        public void handleStep(SpacecraftState currentState, boolean isLast) {
            sample.add(currentState);
        }
    });
    propagator.propagate(orbit.getDate().shiftedBy(600.0), orbit.getDate().shiftedBy(4200.0));
    JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder, 1.0e-3, 5000);
    try {
        Method setSample = AbstractPropagatorConverter.class.getDeclaredMethod("setSample", List.class);
        setSample.setAccessible(true);
        setSample.invoke(fitter, sample);
    } catch (NoSuchMethodException | SecurityException | IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
        Assert.fail(e.getLocalizedMessage());
    }
    MultivariateVectorFunction f = fitter.getObjectiveFunction();
    Pair<RealVector, RealMatrix> p = fitter.getModel().value(new ArrayRealVector(normalized));
    // check derivatives
    // a h offset on normalized parameter represents a physical offset of h * scale
    RealMatrix m = p.getSecond();
    double h = 10.0;
    double[] shifted = normalized.clone();
    double maxErrorP = 0;
    double maxErrorV = 0;
    for (int j = 0; j < selected.size(); ++j) {
        shifted[j] = normalized[j] + 2.0 * h;
        double[] valueP2 = f.value(shifted);
        shifted[j] = normalized[j] + 1.0 * h;
        double[] valueP1 = f.value(shifted);
        shifted[j] = normalized[j] - 1.0 * h;
        double[] valueM1 = f.value(shifted);
        shifted[j] = normalized[j] - 2.0 * h;
        double[] valueM2 = f.value(shifted);
        shifted[j] = normalized[j];
        for (int i = 0; i < valueP2.length; ++i) {
            double d = (8 * (valueP1[i] - valueM1[i]) - (valueP2[i] - valueM2[i])) / (12 * h);
            if (i % 6 < 3) {
                // position
                maxErrorP = FastMath.max(maxErrorP, FastMath.abs(m.getEntry(i, j) - d));
            } else {
                // velocity
                maxErrorV = FastMath.max(maxErrorV, FastMath.abs(m.getEntry(i, j) - d));
            }
        }
    }
    Assert.assertEquals(0.0, maxErrorP, tolP);
    Assert.assertEquals(0.0, maxErrorV, tolV);
}

Example 67 with RandomGenerator

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

the class TabulatedLofOffsetTest method testConstantOffset.

@Test
public void testConstantOffset() throws OrekitException {
    RandomGenerator random = new Well19937a(0x1199d4bb8f53d2b6l);
    for (LOFType type : LOFType.values()) {
        for (int i = 0; i < 100; ++i) {
            double a1 = random.nextDouble() * MathUtils.TWO_PI;
            double a2 = random.nextDouble() * MathUtils.TWO_PI;
            double a3 = random.nextDouble() * MathUtils.TWO_PI;
            LofOffset law = new LofOffset(orbit.getFrame(), type, RotationOrder.XYZ, a1, a2, a3);
            Rotation offsetAtt = law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()).getRotation();
            LofOffset aligned = new LofOffset(orbit.getFrame(), type);
            Rotation alignedAtt = aligned.getAttitude(orbit, orbit.getDate(), orbit.getFrame()).getRotation();
            Rotation offsetProper = offsetAtt.compose(alignedAtt.revert(), RotationConvention.VECTOR_OPERATOR);
            TabulatedLofOffset tabulated = new TabulatedLofOffset(orbit.getFrame(), type, Arrays.asList(new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(-10), offsetProper, Vector3D.ZERO, Vector3D.ZERO), new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(0), offsetProper, Vector3D.ZERO, Vector3D.ZERO), new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(+10), offsetProper, Vector3D.ZERO, Vector3D.ZERO)), 2, AngularDerivativesFilter.USE_R);
            Rotation rebuilt = tabulated.getAttitude(orbit, orbit.getDate(), orbit.getFrame()).getRotation();
            Assert.assertEquals(0.0, Rotation.distance(offsetAtt, rebuilt), 1.2e-15);
        }
    }
}

Example 68 with RandomGenerator

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

the class EllipsoidTest method testRandomNormalSections.

@Test
public void testRandomNormalSections() throws IOException {
    RandomGenerator random = new Well19937a(0x573c54d152aeafe4l);
    for (int i = 0; i < 100; ++i) {
        double a = 10 * random.nextDouble();
        double b = 10 * random.nextDouble();
        double c = 10 * random.nextDouble();
        double size = FastMath.max(FastMath.max(a, b), c);
        final Ellipsoid ellipsoid = new Ellipsoid(FramesFactory.getEME2000(), a, b, c);
        for (int j = 0; j < 50; ++j) {
            double phi = FastMath.PI * (random.nextDouble() - 0.5);
            double lambda = 2 * FastMath.PI * random.nextDouble();
            double cPhi = FastMath.cos(phi);
            double sPhi = FastMath.sin(phi);
            double cLambda = FastMath.cos(lambda);
            double sLambda = FastMath.sin(lambda);
            Vector3D surfacePoint = new Vector3D(ellipsoid.getA() * cPhi * cLambda, ellipsoid.getB() * cPhi * sLambda, ellipsoid.getC() * sPhi);
            Vector3D t1 = new Vector3D(-ellipsoid.getA() * cPhi * sLambda, ellipsoid.getB() * cPhi * cLambda, 0).normalize();
            Vector3D t2 = new Vector3D(-ellipsoid.getA() * sPhi * cLambda, -ellipsoid.getB() * sPhi * sLambda, ellipsoid.getC() * cPhi).normalize();
            final double azimuth = 2 * FastMath.PI * random.nextDouble();
            double cAzimuth = FastMath.cos(azimuth);
            double sAzimuth = FastMath.sin(azimuth);
            Vector3D tAz = new Vector3D(cAzimuth, t1, sAzimuth, t2);
            final Ellipse ps = ellipsoid.getPlaneSection(surfacePoint, tAz);
            Assert.assertEquals(0.0, errorOnEllipsoid(ps, ellipsoid), 1.0e-12 * size);
            Assert.assertEquals(0.0, errorOnPlane(ps, surfacePoint, tAz), 1.0e-10 * size);
            double cos = Vector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getU()) / ps.getA();
            double sin = Vector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getV()) / ps.getB();
            final Vector3D rebuilt = ps.pointAt(FastMath.atan2(sin, cos));
            Assert.assertEquals(0, Vector3D.distance(surfacePoint, rebuilt), 1.0e-11 * size);
        }
    }
}

Example 69 with RandomGenerator

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

the class NutationCodecTest method testKeySymmetry.

@Test
public void testKeySymmetry() {
    RandomGenerator random = new Well1024a(0x8fef7f6f99ad5d56l);
    int[] multipliers = new int[15];
    for (int i = 0; i < 100000; ++i) {
        Arrays.fill(multipliers, 0);
        int nb = 1 + random.nextInt(7);
        for (int k = 0; k < nb; ++k) {
            int index = random.nextInt(15);
            while (multipliers[index] == 0) {
                multipliers[index] = random.nextInt(128) - 64;
            }
        }
        long key = NutationCodec.encode(multipliers);
        int[] rebuilt = NutationCodec.decode(key);
        for (int k = 0; k < multipliers.length; ++k) {
            Assert.assertEquals(multipliers[k], rebuilt[k]);
        }
    }
}

Example 70 with RandomGenerator

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

the class FieldTransformTest method doTestJacobianPV.

private <T extends RealFieldElement<T>> void doTestJacobianPV(Field<T> field) {
    // base directions for finite differences
    @SuppressWarnings("unchecked") FieldPVCoordinates<T>[] directions = (FieldPVCoordinates<T>[]) Array.newInstance(FieldPVCoordinates.class, 6);
    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));
    double h = 0.01;
    RandomGenerator random = new Well19937a(0xce2bfddfbb9796bel);
    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_PV, 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();
            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);
                // check the rest of the matrix remains untouched
                for (int l = 6; 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);
                }
            }
        }
    }
}