use of org.hipparchus.random.RandomGenerator 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);
}
}
}
use of org.hipparchus.random.RandomGenerator 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);
}
}
}
use of org.hipparchus.random.RandomGenerator 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);
}
}
}
}
use of org.hipparchus.random.RandomGenerator 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);
}
}
use of org.hipparchus.random.RandomGenerator 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());
}
}
}
Aggregations