Commit db0d964c authored by Bryan Cazabonne's avatar Bryan Cazabonne
Browse files

Added computation of state transition matrix using keplerian model.

parent d39253e8
Pipeline #1666 failed with stages
in 21 minutes and 30 seconds
/* Copyright 2002-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.analytical;
import java.util.Collections;
import java.util.List;
import org.hipparchus.analysis.differentiation.Gradient;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.utils.ParameterDriver;
/**
* Converter for Keplerian propagator.
*
* @author Bryan Cazabonne
* @since 11.1
*/
class KeplerianGradientConverter extends AbstractAnalyticalGradientConverter {
/** Fixed dimension of the state. */
public static final int FREE_STATE_PARAMETERS = 6;
/** Orbit propagator. */
private final KeplerianPropagator propagator;
/** Simple constructor.
* @param propagator orbit propagator used to access initial orbit
*/
KeplerianGradientConverter(final KeplerianPropagator propagator) {
super(propagator, propagator.getInitialState().getMu(), FREE_STATE_PARAMETERS);
// Initialize fields
this.propagator = propagator;
}
/** {@inheritDoc} */
@Override
public FieldKeplerianPropagator<Gradient> getPropagator(final FieldSpacecraftState<Gradient> state,
final Gradient[] parameters) {
// Zero
final Gradient zero = state.getA().getField().getZero();
// Model parameters
final AttitudeProvider provider = propagator.getAttitudeProvider();
// Central attraction coefficient
final Gradient mu = zero.add(propagator.getInitialState().getMu());
// Return the "Field" propagator
return new FieldKeplerianPropagator<Gradient>(state.getOrbit(), provider, mu);
}
/** {@inheritDoc} */
@Override
public List<ParameterDriver> getParametersDrivers() {
return Collections.emptyList();
}
}
/* Copyright 2002-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.analytical;
import org.hipparchus.linear.RealMatrix;
import org.orekit.utils.DoubleArrayDictionary;
/**
* Harvester between two-dimensional Jacobian matrices and
* one-dimensional {@link KeplerianPropagator}.
*
* @author Bryan Cazabonne
* @since 11.1
*/
class KeplerianHarvester extends AbstractAnalyticalMatricesHarvester {
/** Propagator bound to this harvester. */
private final KeplerianPropagator propagator;
/** Simple constructor.
* <p>
* The arguments for initial matrices <em>must</em> be compatible with the
* {@link org.orekit.orbits.OrbitType orbit type}
* and {@link org.orekit.orbits.PositionAngle position angle} that will be used by propagator
* </p>
* @param propagator propagator bound to this harvester
* @param stmName State Transition Matrix state name
* @param initialStm initial State Transition Matrix ∂Y/∂Y₀,
* if null (which is the most frequent case), assumed to be 6x6 identity
* @param initialJacobianColumns initial columns of the Jacobians matrix with respect to parameters,
* if null or if some selected parameters are missing from the dictionary, the corresponding
* initial column is assumed to be 0
*/
KeplerianHarvester(final KeplerianPropagator propagator, final String stmName,
final RealMatrix initialStm, final DoubleArrayDictionary initialJacobianColumns) {
super(propagator, stmName, initialStm, initialJacobianColumns);
this.propagator = propagator;
}
/** {@inheritDoc} */
@Override
public AbstractAnalyticalGradientConverter getGradientConverter() {
return new KeplerianGradientConverter(propagator);
}
}
......@@ -16,12 +16,14 @@
*/
package org.orekit.propagation.analytical;
import org.hipparchus.linear.RealMatrix;
import org.orekit.attitudes.Attitude;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.InertialProvider;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.AbstractMatricesHarvester;
import org.orekit.propagation.SpacecraftState;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.DoubleArrayDictionary;
......@@ -192,4 +194,11 @@ public class KeplerianPropagator extends AbstractAnalyticalPropagator {
return states.get(date).getMass();
}
/** {@inheritDoc} */
@Override
protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
final DoubleArrayDictionary initialJacobianColumns) {
return new KeplerianHarvester(this, stmName, initialStm, initialJacobianColumns);
}
}
......@@ -16,18 +16,28 @@
*/
package org.orekit.propagation.conversion;
import java.util.List;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.InertialProvider;
import org.orekit.estimation.leastsquares.AbstractBatchLSModel;
import org.orekit.estimation.leastsquares.AnalyticalBatchLSModel;
import org.orekit.estimation.leastsquares.ModelObserver;
import org.orekit.estimation.measurements.ObservedMeasurement;
import org.orekit.estimation.sequential.AbstractKalmanModel;
import org.orekit.estimation.sequential.CovarianceMatrixProvider;
import org.orekit.estimation.sequential.KalmanModel;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.utils.ParameterDriversList;
/** Builder for Keplerian propagator.
* @author Pascal Parraud
* @since 6.0
*/
public class KeplerianPropagatorBuilder extends AbstractPropagatorBuilder {
public class KeplerianPropagatorBuilder extends AbstractPropagatorBuilder implements OrbitDeterminationPropagatorBuilder {
/** Build a new instance.
* <p>
......@@ -81,4 +91,22 @@ public class KeplerianPropagatorBuilder extends AbstractPropagatorBuilder {
return new KeplerianPropagator(createInitialOrbit(), getAttitudeProvider());
}
/** {@inheritDoc} */
@Override
public AbstractBatchLSModel buildLSModel(final OrbitDeterminationPropagatorBuilder[] builders,
final List<ObservedMeasurement<?>> measurements,
final ParameterDriversList estimatedMeasurementsParameters,
final ModelObserver observer) {
return new AnalyticalBatchLSModel(builders, measurements, estimatedMeasurementsParameters, observer);
}
/** {@inheritDoc} */
@Override
public AbstractKalmanModel buildKalmanModel(final List<OrbitDeterminationPropagatorBuilder> propagatorBuilders,
final List<CovarianceMatrixProvider> covarianceMatricesProviders,
final ParameterDriversList estimatedMeasurementsParameters,
final CovarianceMatrixProvider measurementProcessNoiseMatrix) {
return new KalmanModel(propagatorBuilders, covarianceMatricesProviders, estimatedMeasurementsParameters, measurementProcessNoiseMatrix);
}
}
......@@ -126,7 +126,8 @@ want to use the estimated covariance and state vector of an orbit determination
the results of the daily orbit determination.
The two batch least squares estimator implementations in Orekit are compatible to work with a `NumericalPropagator`, a
`DSSTPropagator`, a `TLEPropagator`, an `EcksteinHechlerPropagator`, or a `BrouwerLyddanePropagator`.
`DSSTPropagator`, a `TLEPropagator`, an `EcksteinHechlerPropagator`, a `BrouwerLyddanePropagator`,
or a `KeplerianPropagator`.
### Kalman filter
......
......@@ -206,7 +206,7 @@
* measurements parameters estimation (biases, satellite clock offset, station clock offset,
station position, pole motion and rate, prime meridian correction and rate, total zenith
delay in tropospheric correction)
* can be used with numerical, DSST, SDP4/SGP4, Eckstein-Hechler, or Brouwer-Lyddane propagators
* can be used with numerical, DSST, SDP4/SGP4, Eckstein-Hechler, Brouwer-Lyddane, or Keplerian propagators
* multi-satellites orbit determination
* initial orbit determination methods (Gibbs, Gooding, Lambert and Laplace)
* ground stations displacements due to solid tides
......
package org.orekit.propagation.analytical;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.util.FastMath;
import org.junit.Assert;
import org.junit.Before;
import org.junit.Test;
import org.orekit.Utils;
import org.orekit.attitudes.Attitude;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.EquinoctialOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.MatricesHarvester;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.PVCoordinates;
public class KeplerianStateTransitionMatrixTest {
@Before
public void setUp() {
Utils.setDataRoot("regular-data");
}
@Test(expected=OrekitException.class)
public void testNullStmName() {
// Definition of initial conditions with position and velocity
//------------------------------------------------------------
Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
FramesFactory.getEME2000(), initDate, Constants.WGS84_EARTH_MU);
// Extrapolator definition
// -----------------------
KeplerianPropagator extrapolator = new KeplerianPropagator(initialOrbit);
extrapolator.setupMatricesComputation(null, null, null);
}
@Test
public void testStateJacobian() {
// Definition of initial conditions with position and velocity
//------------------------------------------------------------
Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
FramesFactory.getEME2000(), initDate, Constants.WGS84_EARTH_MU);
// compute state Jacobian using PartialDerivatives
KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit);
final SpacecraftState initialState = propagator.getInitialState();
final double[] stateVector = new double[6];
OrbitType.CARTESIAN.mapOrbitToArray(initialState.getOrbit(), PositionAngle.MEAN, stateVector, null);
final AbsoluteDate target = initialState.getDate().shiftedBy(initialState.getKeplerianPeriod());
MatricesHarvester harvester = propagator.setupMatricesComputation("stm", null, null);
final SpacecraftState finalState = propagator.propagate(target);
RealMatrix dYdY0 = harvester.getStateTransitionMatrix(finalState);
// compute reference state Jacobian using finite differences
double[][] dYdY0Ref = new double[6][6];
KeplerianPropagator propagator2;
double[] steps = NumericalPropagator.tolerances(10, initialState.getOrbit(), OrbitType.CARTESIAN)[0];
for (int i = 0; i < 6; ++i) {
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, -4 * steps[i], i).getOrbit());
SpacecraftState sM4h = propagator2.propagate(target);
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, -3 * steps[i], i).getOrbit());
SpacecraftState sM3h = propagator2.propagate(target);
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, -2 * steps[i], i).getOrbit());
SpacecraftState sM2h = propagator2.propagate(target);
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, -1 * steps[i], i).getOrbit());
SpacecraftState sM1h = propagator2.propagate(target);
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, +1 * steps[i], i).getOrbit());
SpacecraftState sP1h = propagator2.propagate(target);
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, +2 * steps[i], i).getOrbit());
SpacecraftState sP2h = propagator2.propagate(target);
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, +3 * steps[i], i).getOrbit());
SpacecraftState sP3h = propagator2.propagate(target);
propagator2 = new KeplerianPropagator(shiftState(initialState, OrbitType.CARTESIAN, +4 * steps[i], i).getOrbit());
SpacecraftState sP4h = propagator2.propagate(target);
fillJacobianColumn(dYdY0Ref, i, OrbitType.CARTESIAN, steps[i],
sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
}
// Verify
for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) {
if (stateVector[i] != 0) {
double error = FastMath.abs((dYdY0.getEntry(i, j) - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
Assert.assertEquals(0, error, 4.51e-14);
}
}
}
}
private void fillJacobianColumn(double[][] jacobian, int column,
OrbitType orbitType, double h,
SpacecraftState sM4h, SpacecraftState sM3h,
SpacecraftState sM2h, SpacecraftState sM1h,
SpacecraftState sP1h, SpacecraftState sP2h,
SpacecraftState sP3h, SpacecraftState sP4h) {
double[] aM4h = stateToArray(sM4h, orbitType)[0];
double[] aM3h = stateToArray(sM3h, orbitType)[0];
double[] aM2h = stateToArray(sM2h, orbitType)[0];
double[] aM1h = stateToArray(sM1h, orbitType)[0];
double[] aP1h = stateToArray(sP1h, orbitType)[0];
double[] aP2h = stateToArray(sP2h, orbitType)[0];
double[] aP3h = stateToArray(sP3h, orbitType)[0];
double[] aP4h = stateToArray(sP4h, orbitType)[0];
for (int i = 0; i < jacobian.length; ++i) {
jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
32 * (aP3h[i] - aM3h[i]) -
168 * (aP2h[i] - aM2h[i]) +
672 * (aP1h[i] - aM1h[i])) / (840 * h);
}
}
private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
double delta, int column) {
double[][] array = stateToArray(state, orbitType);
array[0][column] += delta;
return arrayToState(array, state.getFrame(), state.getDate(),
state.getMu(), state.getAttitude());
}
private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
double[][] array = new double[2][6];
orbitType.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, array[0], array[1]);
return array;
}
private SpacecraftState arrayToState(double[][] array,
Frame frame, AbsoluteDate date, double mu,
Attitude attitude) {
CartesianOrbit orbit = (CartesianOrbit) OrbitType.CARTESIAN.mapArrayToOrbit(array[0], array[1], PositionAngle.MEAN, date, mu, frame);
return new SpacecraftState(orbit, attitude);
}
}
Markdown is supported
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