Commit 714ff57f authored by Francesco Rocca's avatar Francesco Rocca Committed by Luc Maisonobe

Added a way to fit Two-Line Elements against a spacecraft states sample.

parent 1cf76731
......@@ -112,6 +112,9 @@
<contributor>
<name>Silvia R&#237;os Berganti&#241;os</name>
</contributor>
<contributor>
<name>Francesco Rocca</name>
</contributor>
<contributor>
<name>Mathieu Rom&#233;ro</name>
</contributor>
......
/* Copyright 2002-2012 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.tle;
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.util.FastMath;
import org.orekit.errors.OrekitException;
/** Orbit converter for Two-Lines Elements using differential algorithm.
* @author Rocca
* @since 6.0
*/
public class DifferentialOrbitConverter extends AbstractTLEFitter {
/** Maximum number of iterations for fitting. */
private final int maxIterations;
/** Simple constructor.
* @param maxIterations maximum number of iterations for fitting
* @param satelliteNumber satellite number
* @param classification classification (U for unclassified)
* @param launchYear launch year (all digits)
* @param launchNumber launch number
* @param launchPiece launch piece
* @param elementNumber element number
* @param revolutionNumberAtEpoch revolution number at epoch
*/
protected DifferentialOrbitConverter(final int maxIterations,
final int satelliteNumber, final char classification,
final int launchYear, final int launchNumber, final String launchPiece,
final int elementNumber, final int revolutionNumberAtEpoch) {
super(satelliteNumber, classification,
launchYear, launchNumber, launchPiece, elementNumber, revolutionNumberAtEpoch);
this.maxIterations = maxIterations;
}
/** {@inheritDoc} */
protected double[] fit(final double[] initial) throws OrekitException, MaxCountExceededException {
final double[] w = getWeight();
for (int i = 0; i < w.length; ++i) {
w[i] = FastMath.sqrt(w[i]);
}
final MultivariateMatrixFunction jacobian = getPVFunction().jacobian();
final double[] result = initial.clone();
double previousRMS = Double.NaN;
for (int iterations = 0; iterations < maxIterations; ++iterations) {
final RealMatrix A = new Array2DRowRealMatrix(jacobian.value(result));
for (int i = 0; i < A.getRowDimension(); i++) {
for (int j = 0; j < A.getColumnDimension(); j++) {
A.multiplyEntry(i, j, w[i]);
}
}
final double[] residuals = getResiduals(result);
final RealVector y = new ArrayRealVector(residuals.length);
for (int i = 0; i < y.getDimension(); i++) {
y.setEntry(i, residuals[i] * w[i]);
}
final RealVector dx = new QRDecomposition(A).getSolver().solve(y);
for (int i = 0; i < result.length; i++) {
result[i] = result[i] + dx.getEntry(i);
}
final double rms = getRMS(result);
if (iterations > 0 && FastMath.abs(rms - previousRMS) <= getPositionTolerance()) {
return result;
}
previousRMS = rms;
}
throw new MaxCountExceededException(maxIterations);
}
}
/* Copyright 2002-2012 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.tle;
import org.apache.commons.math3.optimization.ConvergenceChecker;
import org.apache.commons.math3.optimization.PointVectorValuePair;
import org.apache.commons.math3.optimization.SimpleVectorValueChecker;
import org.apache.commons.math3.optimization.general.LevenbergMarquardtOptimizer;
import org.orekit.errors.OrekitException;
/** Orbit converter for Two-Lines Elements using differential algorithm.
* @author Rocca
* @since 6.0
*/
public class LevenbergMarquardtOrbitConverter extends AbstractTLEFitter {
/** Maximum number of iterations for fitting. */
private final int maxIterations;
/** Simple constructor.
* @param maxIterations maximum number of iterations for fitting
* @param satelliteNumber satellite number
* @param classification classification (U for unclassified)
* @param launchYear launch year (all digits)
* @param launchNumber launch number
* @param launchPiece launch piece
* @param elementNumber element number
* @param revolutionNumberAtEpoch revolution number at epoch
*/
protected LevenbergMarquardtOrbitConverter(final int maxIterations,
final int satelliteNumber, final char classification,
final int launchYear, final int launchNumber, final String launchPiece,
final int elementNumber, final int revolutionNumberAtEpoch) {
super(satelliteNumber, classification,
launchYear, launchNumber, launchPiece, elementNumber, revolutionNumberAtEpoch);
this.maxIterations = maxIterations;
}
/** {@inheritDoc} */
protected double[] fit(final double[] initial) throws OrekitException {
final ConvergenceChecker<PointVectorValuePair> checker =
new SimpleVectorValueChecker(-1.0, getPositionTolerance());
final LevenbergMarquardtOptimizer optimizer =
new LevenbergMarquardtOptimizer(checker);
final PointVectorValuePair optimum =
optimizer.optimize(maxIterations, getPVFunction(),
getTarget(), getWeight(), initial);
return optimum.getPointRef();
}
}
......@@ -22,6 +22,10 @@
<body>
<release version="6.0" date="TBD"
description="TBD.">
<action dev="luc" type="add" due-to="Francesco Rocca">
Added two different fitters to fit Two-Line Elements against a sample of
spacecraft states.
</action>
<action dev="luc" type="fix">
Improved testing of error messages.
</action>
......
/* Copyright 2002-2012 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.tle;
import java.util.ArrayList;
import java.util.List;
import junit.framework.Assert;
import org.junit.Before;
import org.orekit.Utils;
import org.orekit.errors.OrekitException;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.SpacecraftState;
public abstract class AbstractTLEFitterTest {
private TLE geoTLE;
private TLE leoTLE;
protected TLE getGeoTLE() {
return geoTLE;
}
protected TLE getLeoTLE() {
return leoTLE;
}
protected abstract AbstractTLEFitter getFitter(TLE tle);
protected void checkFit(final TLE tle, final double duration, final double stepSize,
final double positionTolerance, final boolean positionOnly,
final boolean withBStar, final double expectedRMS)
throws OrekitException {
Propagator p = TLEPropagator.selectExtrapolator(tle);
List<SpacecraftState> sample = new ArrayList<SpacecraftState>();
for (double dt = 0; dt < duration; dt += stepSize) {
sample.add(p.propagate(tle.getDate().shiftedBy(dt)));
}
AbstractTLEFitter fitter = getFitter(tle);
fitter.toTLE(sample, positionTolerance, positionOnly, withBStar);
Assert.assertEquals(expectedRMS, fitter.getRMS(), 0.01 * expectedRMS);
TLE fitted = fitter.getTLE();
Assert.assertEquals(tle.getSatelliteNumber(), fitted.getSatelliteNumber());
Assert.assertEquals(tle.getClassification(), fitted.getClassification());
Assert.assertEquals(tle.getLaunchYear(), fitted.getLaunchYear());
Assert.assertEquals(tle.getLaunchNumber(), fitted.getLaunchNumber());
Assert.assertEquals(tle.getLaunchPiece(), fitted.getLaunchPiece());
Assert.assertEquals(tle.getElementNumber(), fitted.getElementNumber());
Assert.assertEquals(tle.getRevolutionNumberAtEpoch(), fitted.getRevolutionNumberAtEpoch());
final double eps = 1.0e-5;
Assert.assertEquals(tle.getMeanMotion(), fitted.getMeanMotion(), eps * tle.getMeanMotion());
Assert.assertEquals(tle.getE(), fitted.getE(), eps * tle.getE());
Assert.assertEquals(tle.getI(), fitted.getI(), eps * tle.getI());
Assert.assertEquals(tle.getPerigeeArgument(), fitted.getPerigeeArgument(), eps * tle.getPerigeeArgument());
Assert.assertEquals(tle.getRaan(), fitted.getRaan(), eps * tle.getRaan());
Assert.assertEquals(tle.getMeanAnomaly(), fitted.getMeanAnomaly(), eps * tle.getMeanAnomaly());
}
@Before
public void setUp() throws OrekitException {
Utils.setDataRoot("regular-data");
geoTLE =
new TLE("1 27508U 02040A 12021.25695307 -.00000113 00000-0 10000-3 0 7326",
"2 27508 0.0571 356.7800 0005033 344.4621 218.7816 1.00271798 34501");
leoTLE =
new TLE("1 31135U 07013A 11003.00000000 .00000816 00000+0 47577-4 0 11",
"2 31135 2.4656 183.9084 0021119 236.4164 60.4567 15.10546832 15");
}
}
/* Copyright 2002-2012 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.tle;
import org.junit.Test;
import org.orekit.errors.OrekitException;
public class DifferentialOrbitConverterTest extends AbstractTLEFitterTest {
@Test
public void testConversionGeoPositionVelocity() throws OrekitException {
checkFit(getGeoTLE(), 86400, 300, 1.0e-3, false, false, 1.234e-7);
}
@Test
public void testConversionGeoPositionOnly() throws OrekitException {
checkFit(getGeoTLE(), 86400, 300, 1.0e-3, true, false, 1.213e-7);
}
@Test
public void testConversionLeoPositionVelocityWithoutBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, false, false, 10.77);
}
@Test
public void testConversionLeoPositionOnlyWithoutBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, true, false, 15.23);
}
@Test
public void testConversionLeoPositionVelocityWithBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, false, true, 2.154e-5);
}
@Test
public void testConversionLeoPositionOnlyWithBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, true, true, 1.147e-6);
}
protected AbstractTLEFitter getFitter(final TLE tle) {
return new DifferentialOrbitConverter(1000,
tle.getSatelliteNumber(), tle.getClassification(),
tle.getLaunchYear(), tle.getLaunchNumber(), tle.getLaunchPiece(),
tle.getElementNumber(), tle.getRevolutionNumberAtEpoch());
}
}
/* Copyright 2002-2012 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.tle;
import org.junit.Test;
import org.orekit.errors.OrekitException;
public class LevenbergMarquardtOrbitConverterTest extends AbstractTLEFitterTest {
@Test
public void testConversionGeoPositionVelocity() throws OrekitException {
checkFit(getGeoTLE(), 86400, 300, 1.0e-3, false, false, 8.600e-8);
}
@Test
public void testConversionGeoPositionOnly() throws OrekitException {
checkFit(getGeoTLE(), 86400, 300, 1.0e-3, true, false, 1.059e-7);
}
@Test
public void testConversionLeoPositionVelocityWithoutBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, false, false, 10.77);
}
@Test
public void testConversionLeoPositionOnlyWithoutBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, true, false, 15.23);
}
@Test
public void testConversionLeoPositionVelocityWithBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, false, true, 9.920e-7);
}
@Test
public void testConversionLeoPositionOnlyWithBStar() throws OrekitException {
checkFit(getLeoTLE(), 86400, 300, 1.0e-3, true, true, 1.147e-6);
}
protected AbstractTLEFitter getFitter(final TLE tle) {
return new LevenbergMarquardtOrbitConverter(1000,
tle.getSatelliteNumber(), tle.getClassification(),
tle.getLaunchYear(), tle.getLaunchNumber(), tle.getLaunchPiece(),
tle.getElementNumber(), tle.getRevolutionNumberAtEpoch());
}
}
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