Commit e2c5883f authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Fixed harmonic accelerations multipliers.

parent 03f77b23
......@@ -60,9 +60,6 @@ import org.orekit.propagation.events.FieldEventDetector;
*/
public abstract class AbstractParametricAcceleration extends AbstractForceModel {
/** Suffix for parameter name. */
public static final String PARAMETER_SUFFIX = " parametric acceleration";
/** Direction of the acceleration in defining frame. */
private final Vector3D direction;
......
......@@ -64,15 +64,15 @@ public class HarmonicParametricAcceleration extends AbstractParametricAccelerati
/** Reference date for computing phase. */
private AbsoluteDate referenceDate;
/** Angular frequency ω = 2π/(kT). */
/** Angular frequency ω = 2kπ/T. */
private final double omega;
/** Simple constructor.
* <p>
* The signed amplitude of the acceleration is a sin[2π(t-t₀)/(kT) + φ], where
* a is parameter {@code 0} and represents the full amplitude, t is current
* The signed amplitude of the acceleration is γ sin[2kπ(t-t₀)/T + φ], where
* γ is parameter {@code 0} and represents the full amplitude, t is current
* date, t₀ is reference date, {@code T} is fundamental period, {@code k} is
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* The value t-t₀ is in seconds.
* </p>
* <p>
......@@ -112,10 +112,10 @@ public class HarmonicParametricAcceleration extends AbstractParametricAccelerati
/** Simple constructor.
* <p>
* The signed amplitude of the acceleration is a sin[2π(t-t₀)/(kT) + φ], where
* a is parameter {@code 0} and represents the full amplitude, t is current
* The signed amplitude of the acceleration is γ sin[2kπ(t-t₀)/T + φ], where
* γ is parameter {@code 0} and represents the full amplitude, t is current
* date, t₀ is reference date, {@code T} is fundamental period, {@code k} is
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* The value t-t₀ is in seconds.
* </p>
* <p>
......@@ -143,8 +143,8 @@ public class HarmonicParametricAcceleration extends AbstractParametricAccelerati
* fundamental period)
*/
public HarmonicParametricAcceleration(final Vector3D direction, final AttitudeProvider attitudeOverride,
final String prefix, final AbsoluteDate referenceDate,
final double fundamentalPeriod, final int harmonicMultiplier) {
final String prefix, final AbsoluteDate referenceDate,
final double fundamentalPeriod, final int harmonicMultiplier) {
this(direction, false, attitudeOverride, prefix, referenceDate,
fundamentalPeriod, harmonicMultiplier);
}
......@@ -168,17 +168,17 @@ public class HarmonicParametricAcceleration extends AbstractParametricAccelerati
* fundamental period)
*/
private HarmonicParametricAcceleration(final Vector3D direction, final boolean isInertial,
final AttitudeProvider attitudeOverride,
final String prefix, final AbsoluteDate referenceDate,
final double fundamentalPeriod, final int harmonicMultiplier) {
final AttitudeProvider attitudeOverride,
final String prefix, final AbsoluteDate referenceDate,
final double fundamentalPeriod, final int harmonicMultiplier) {
super(direction, isInertial, attitudeOverride);
this.referenceDate = referenceDate;
this.omega = MathUtils.TWO_PI / (harmonicMultiplier * fundamentalPeriod);
this.omega = harmonicMultiplier * MathUtils.TWO_PI / fundamentalPeriod;
try {
drivers = new ParameterDriver[] {
new ParameterDriver(prefix + " amplitude" + PARAMETER_SUFFIX,
0.0, AMPLITUDE_SCALE, 0.0, Double.POSITIVE_INFINITY),
new ParameterDriver(prefix + " phase" + PARAMETER_SUFFIX,
new ParameterDriver(prefix + " γ",
0.0, AMPLITUDE_SCALE, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY),
new ParameterDriver(prefix + " φ",
0.0, PHASE_SCALE, -MathUtils.TWO_PI, MathUtils.TWO_PI),
};
} catch (OrekitException oe) {
......@@ -197,10 +197,10 @@ public class HarmonicParametricAcceleration extends AbstractParametricAccelerati
}
/** {@inheritDoc}.
* The signed amplitude of the acceleration is a sin[2π(t-t₀)/(kT) + φ], where
* a is parameter {@code 0} and represents the full amplitude, t is current
* The signed amplitude of the acceleration is γ sin[2kπ(t-t₀)/T + φ], where
* γ is parameter {@code 0} and represents the full amplitude, t is current
* date, t₀ is reference date, {@code T} is fundamental period, {@code k} is
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* The value t-t₀ is in seconds.
*/
@Override
......@@ -210,10 +210,10 @@ public class HarmonicParametricAcceleration extends AbstractParametricAccelerati
}
/** {@inheritDoc}
* The signed amplitude of the acceleration is a sin[2π(t-t₀)/(kT) + φ], where
* a is parameter {@code 0} and represents the full amplitude, t is current
* The signed amplitude of the acceleration is γ sin[2kπ(t-t₀)/T + φ], where
* γ is parameter {@code 0} and represents the full amplitude, t is current
* date, t₀ is reference date, {@code T} is fundamental period, {@code k} is
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* harmonic multiplier, and φ is parameter {@code 1} and represents phase at t₀.
* The value t-t₀ is in seconds.
*/
@Override
......
......@@ -133,8 +133,7 @@ public class PolynomialParametricAcceleration extends AbstractParametricAccelera
this.drivers = new ParameterDriver[degree + 1];
try {
for (int i = 0; i < drivers.length; ++i) {
drivers[i] = new ParameterDriver(prefix + "[" + i + "]" + PARAMETER_SUFFIX,
0.0, ACCELERATION_SCALE,
drivers[i] = new ParameterDriver(prefix + "[" + i + "]", 0.0, ACCELERATION_SCALE,
Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
}
} catch (OrekitException oe) {
......
......@@ -17,7 +17,9 @@
package org.orekit.forces;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
......@@ -28,6 +30,7 @@ import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
import org.hipparchus.util.Decimal64Field;
import org.hipparchus.util.FastMath;
import org.junit.Assert;
......@@ -40,10 +43,14 @@ import org.orekit.attitudes.InertialProvider;
import org.orekit.attitudes.LofOffset;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.errors.OrekitException;
import org.orekit.estimation.leastsquares.BatchLSEstimator;
import org.orekit.estimation.measurements.ObservedMeasurement;
import org.orekit.estimation.measurements.PV;
import org.orekit.forces.maneuvers.ConstantThrustManeuver;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.LOFType;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.CircularOrbit;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngle;
......@@ -51,6 +58,8 @@ import org.orekit.propagation.FieldBoundedPropagator;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.PropagatorsParallelizer;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder;
import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
import org.orekit.propagation.numerical.FieldNumericalPropagator;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.propagation.sampling.MultiSatStepHandler;
......@@ -61,6 +70,7 @@ import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.ParameterDriver;
public class HarmonicParametricAccelerationTest extends AbstractForceModelTest {
......@@ -299,12 +309,12 @@ public class HarmonicParametricAccelerationTest extends AbstractForceModelTest {
@Test
public void testParameterDerivative2T() throws OrekitException {
doTestParameterDerivative(2, 7.0e-14, 5.0e-12);
doTestParameterDerivative(2, 3.0e-14, 7.0e-12);
}
@Test
public void testParameterDerivative3T() throws OrekitException {
doTestParameterDerivative(3, 3.0e-14, 2.0e-12);
doTestParameterDerivative(3, 2.0e-14, 2.0e-11);
}
private void doTestParameterDerivative(final int harmonicMultiplier,
......@@ -328,11 +338,124 @@ public class HarmonicParametricAccelerationTest extends AbstractForceModelTest {
hpa.init(state, state.getDate().shiftedBy(3600.0));
hpa.getParametersDrivers()[0].setValue(0.00001);
hpa.getParametersDrivers()[1].setValue(0.00002);
checkParameterDerivative(state, hpa, "kT amplitude parametric acceleration", 1.0e-3, amplitudeDerivativeTolerance);
checkParameterDerivative(state, hpa, "kT phase parametric acceleration", 1.0e-3, phaseDerivativeTolerance);
checkParameterDerivative(state, hpa, "kT γ", 1.0e-3, amplitudeDerivativeTolerance);
checkParameterDerivative(state, hpa, "kT φ", 1.0e-3, phaseDerivativeTolerance);
}
@Test
public void testCoefficientsDetermination() throws OrekitException {
final double mass = 2500;
final Orbit orbit = new CircularOrbit(7500000.0, 1.0e-4, 1.0e-3, 1.7, 0.3, 0.5, PositionAngle.TRUE,
FramesFactory.getEME2000(),
new AbsoluteDate(new DateComponents(2004, 2, 3), TimeComponents.H00,
TimeScalesFactory.getUTC()),
Constants.EIGEN5C_EARTH_MU);
final double period = orbit.getKeplerianPeriod();
AttitudeProvider maneuverLaw = new LofOffset(orbit.getFrame(), LOFType.VNC);
SpacecraftState initialState = new SpacecraftState(orbit,
maneuverLaw.getAttitude(orbit,
orbit.getDate(),
orbit.getFrame()),
mass);
double dP = 10.0;
double minStep = 0.001;
double maxStep = 100;
double[][] tolerance = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
// generate PV measurements corresponding to a tangential maneuver
AdaptiveStepsizeIntegrator integrator0 =
new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
integrator0.setInitialStepSize(60);
final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
propagator0.setInitialState(initialState);
propagator0.setAttitudeProvider(maneuverLaw);
ForceModel hpaRefX1 = new HarmonicParametricAcceleration(Vector3D.PLUS_I, true, "refX1", null, period, 1);
ForceModel hpaRefY1 = new HarmonicParametricAcceleration(Vector3D.PLUS_J, true, "refY1", null, period, 1);
ForceModel hpaRefZ2 = new HarmonicParametricAcceleration(Vector3D.PLUS_K, true, "refZ2", null, period, 2);
hpaRefX1.getParametersDrivers()[0].setValue(2.4e-2);
hpaRefX1.getParametersDrivers()[1].setValue(3.1);
hpaRefY1.getParametersDrivers()[0].setValue(4.0e-2);
hpaRefY1.getParametersDrivers()[1].setValue(0.3);
hpaRefZ2.getParametersDrivers()[0].setValue(1.0e-2);
hpaRefZ2.getParametersDrivers()[1].setValue(1.8);
propagator0.addForceModel(hpaRefX1);
propagator0.addForceModel(hpaRefY1);
propagator0.addForceModel(hpaRefZ2);
final List<ObservedMeasurement<?>> measurements = new ArrayList<>();
propagator0.setMasterMode(10.0,
(state, isLast) ->
measurements.add(new PV(state.getDate(),
state.getPVCoordinates().getPosition(), state.getPVCoordinates().getVelocity(),
1.0e-3, 1.0e-6, 1.0)));
propagator0.propagate(orbit.getDate().shiftedBy(900));
// set up an estimator to retrieve the maneuver as several harmonic accelerations in inertial frame
final NumericalPropagatorBuilder propagatorBuilder =
new NumericalPropagatorBuilder(orbit,
new DormandPrince853IntegratorBuilder(minStep, maxStep, dP),
PositionAngle.TRUE, dP);
propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_I, true, "X1", null, period, 1));
propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_J, true, "Y1", null, period, 1));
propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_K, true, "Z2", null, period, 2));
final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
estimator.setParametersConvergenceThreshold(1.0e-2);
estimator.setMaxIterations(20);
estimator.setMaxEvaluations(100);
for (final ObservedMeasurement<?> measurement : measurements) {
estimator.addMeasurement(measurement);
}
// we will estimate only the force model parameters, not the orbit
for (final ParameterDriver d : estimator.getOrbitalParametersDrivers(false).getDrivers()) {
d.setSelected(false);
}
setParameter(estimator, "X1 γ", 1.0e-2);
setParameter(estimator, "X1 φ", 4.0);
setParameter(estimator, "Y1 γ", 1.0e-2);
setParameter(estimator, "Y1 φ", 0.0);
setParameter(estimator, "Z2 γ", 1.0e-2);
setParameter(estimator, "Z2 φ", 1.0);
estimator.estimate();
Assert.assertTrue(estimator.getIterationsCount() < 15);
Assert.assertTrue(estimator.getEvaluationsCount() < 15);
Assert.assertEquals(0.0, estimator.getOptimum().getRMS(), 1.0e-5);
Assert.assertEquals(hpaRefX1.getParametersDrivers()[0].getValue(), getParameter(estimator, "X1 γ"), 1.e-12);
Assert.assertEquals(hpaRefX1.getParametersDrivers()[1].getValue(), getParameter(estimator, "X1 φ"), 1.e-12);
Assert.assertEquals(hpaRefY1.getParametersDrivers()[0].getValue(), getParameter(estimator, "Y1 γ"), 1.e-12);
Assert.assertEquals(hpaRefY1.getParametersDrivers()[1].getValue(), getParameter(estimator, "Y1 φ"), 1.e-12);
Assert.assertEquals(hpaRefZ2.getParametersDrivers()[0].getValue(), getParameter(estimator, "Z2 γ"), 1.e-12);
Assert.assertEquals(hpaRefZ2.getParametersDrivers()[1].getValue(), getParameter(estimator, "Z2 φ"), 1.e-12);
}
private void setParameter(BatchLSEstimator estimator, String name, double value)
throws OrekitException {
for (final ParameterDriver driver : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
if (driver.getName().equals(name)) {
driver.setSelected(true);
driver.setValue(value);
return;
}
}
Assert.fail("unknown parameter " + name);
}
private double getParameter(BatchLSEstimator estimator, String name)
throws OrekitException {
for (final ParameterDriver driver : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
if (driver.getName().equals(name)) {
return driver.getValue();
}
}
Assert.fail("unknown parameter " + name);
return Double.NaN;
}
@Before
public void setUp() {
try {
......
......@@ -302,10 +302,10 @@ public class PolynomialParametricAccelerationTest extends AbstractForceModelTest
ppa.getParametersDrivers()[2].setValue(0.00003);
ppa.getParametersDrivers()[3].setValue(0.00004);
checkParameterDerivative(state, ppa, "ppa[0] parametric acceleration", 1.0e-3, 7.0e-12);
checkParameterDerivative(state, ppa, "ppa[1] parametric acceleration", 1.0e-3, 9.0e-13);
checkParameterDerivative(state, ppa, "ppa[2] parametric acceleration", 1.0e-3, 2.0e-12);
checkParameterDerivative(state, ppa, "ppa[3] parametric acceleration", 1.0e-3, 3.0e-12);
checkParameterDerivative(state, ppa, "ppa[0]", 1.0e-3, 7.0e-12);
checkParameterDerivative(state, ppa, "ppa[1]", 1.0e-3, 9.0e-13);
checkParameterDerivative(state, ppa, "ppa[2]", 1.0e-3, 2.0e-12);
checkParameterDerivative(state, ppa, "ppa[3]", 1.0e-3, 3.0e-12);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment