/* Copyright 2002-2017 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.estimation;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import org.hipparchus.RealFieldElement;
import org.hipparchus.geometry.euclidean.threed.FieldRotation;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
import org.hipparchus.util.FastMath;
import org.junit.Assert;
import org.orekit.Utils;
import org.orekit.bodies.CelestialBodyFactory;
//import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataProvidersManager;
import org.orekit.errors.OrekitException;
import org.orekit.estimation.leastsquares.BatchLSEstimator;
import org.orekit.estimation.measurements.EstimatedMeasurement;
import org.orekit.estimation.measurements.GroundStation;
import org.orekit.estimation.measurements.ObservedMeasurement;
import org.orekit.estimation.measurements.MeasurementCreator;
import org.orekit.forces.drag.IsotropicDrag;
import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
import org.orekit.forces.gravity.potential.GRGSFormatReader;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
import org.orekit.forces.radiation.IsotropicRadiationClassicalConvention;
import org.orekit.frames.FieldTransform;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.Transform;
import org.orekit.frames.TransformProvider;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.conversion.PropagatorBuilder;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
/** Utility class for orbit determination tests. */
public class EstimationTestUtils {
public static Context eccentricContext() throws OrekitException {
Utils.setDataRoot("regular-data:potential:tides");
Context context = new Context();
context.conventions = IERSConventions.IERS_2010;
context.earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(context.conventions, true));
context.sun = CelestialBodyFactory.getSun();
context.moon = CelestialBodyFactory.getMoon();
context.radiationSensitive = new IsotropicRadiationClassicalConvention(2.0, 0.2, 0.8);
context.dragSensitive = new IsotropicDrag(2.0, 1.2);
context.utc = TimeScalesFactory.getUTC();
context.ut1 = TimeScalesFactory.getUT1(context.conventions, true);
GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
AstronomicalAmplitudeReader aaReader =
new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
0.01, FastMath.toRadians(1.0),
OceanLoadDeformationCoefficients.IERS_2010,
map));
context.gravity = GravityFieldFactory.getNormalizedProvider(20, 20);
context.initialOrbit = new KeplerianOrbit(15000000.0, 0.125, 1.25,
0.250, 1.375, 0.0625, PositionAngle.TRUE,
FramesFactory.getEME2000(),
new AbsoluteDate(2000, 2, 24, 11, 35,47.0, context.utc),
context.gravity.getMu());
context.stations = Arrays.asList(//context.createStation(-18.59146, -173.98363, 76.0, "Leimatu`a"),
context.createStation(-53.05388, -75.01551, 1750.0, "Isla Desolación"),
context.createStation( 62.29639, -7.01250, 880.0, "Slættaratindur")
//context.createStation( -4.01583, 103.12833, 3173.0, "Gunung Dempo")
);
// Turn-around range stations
// Map entry = master station
// Map value = slave station associated
context.TARstations = new HashMap<GroundStation, GroundStation>();
context.TARstations.put(context.createStation(-53.05388, -75.01551, 1750.0, "Isla Desolación"),
context.createStation(-54.815833, -68.317778, 6.0, "Ushuaïa"));
context.TARstations.put(context.createStation( 62.29639, -7.01250, 880.0, "Slættaratindur"),
context.createStation( 61.405833, -6.705278, 470.0, "Sumba"));
return context;
}
public static Context geoStationnaryContext() throws OrekitException {
Utils.setDataRoot("regular-data:potential:tides");
Context context = new Context();
context.conventions = IERSConventions.IERS_2010;
context.utc = TimeScalesFactory.getUTC();
context.ut1 = TimeScalesFactory.getUT1(context.conventions, true);
String Myframename = "MyEarthFrame";
final AbsoluteDate datedef = new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc);
final double omega = Constants.WGS84_EARTH_ANGULAR_VELOCITY;
final Vector3D rotationRate = new Vector3D(0.0, 0.0, omega);
TransformProvider MyEarthFrame = new TransformProvider() {
private static final long serialVersionUID = 1L;
public Transform getTransform(final AbsoluteDate date) {
final double rotationduration = date.durationFrom(datedef);
final Vector3D alpharot = new Vector3D(rotationduration, rotationRate);
final Rotation rotation = new Rotation(Vector3D.PLUS_K, -alpharot.getZ(),
RotationConvention.VECTOR_OPERATOR);
return new Transform(date, rotation, rotationRate);
}
public <T extends RealFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
final T rotationduration = date.durationFrom(datedef);
final FieldVector3D<T> alpharot = new FieldVector3D<>(rotationduration, rotationRate);
final FieldRotation<T> rotation = new FieldRotation<>(FieldVector3D.getPlusK(date.getField()),
alpharot.getZ().negate(),
RotationConvention.VECTOR_OPERATOR);
return new FieldTransform<>(date, rotation, new FieldVector3D<>(date.getField(), rotationRate));
}
};
Frame FrameTest = new Frame(FramesFactory.getEME2000(), MyEarthFrame, Myframename, true);
// Earth is spherical, rotating in one sidereal day
context.earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0.0, FrameTest);
context.sun = CelestialBodyFactory.getSun();
context.moon = CelestialBodyFactory.getMoon();
context.radiationSensitive = new IsotropicRadiationClassicalConvention(2.0, 0.2, 0.8);
context.dragSensitive = new IsotropicDrag(2.0, 1.2);
GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
AstronomicalAmplitudeReader aaReader =
new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
0.01, FastMath.toRadians(1.0),
OceanLoadDeformationCoefficients.IERS_2010,
map));
context.gravity = GravityFieldFactory.getNormalizedProvider(20, 20);
// semimajor axis for a geostationnary satellite
double da = FastMath.cbrt(context.gravity.getMu() / (omega * omega));
//context.stations = Arrays.asList(context.createStation( 0.0, 0.0, 0.0, "Lat0_Long0"),
// context.createStation( 62.29639, -7.01250, 880.0, "Slættaratindur")
// );
context.stations = Arrays.asList(context.createStation(0.0, 0.0, 0.0, "Lat0_Long0") );
// Station position & velocity in EME2000
final Vector3D geovelocity = new Vector3D (0., 0., 0.);
// Compute the frames transformation from station frame to EME2000
Transform topoToEME =
context.stations.get(0).getOffsetFrame().getTransformTo(FramesFactory.getEME2000(),new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc));
// Station position in EME2000 at reference date
Vector3D stationPositionEME = topoToEME.transformPosition(Vector3D.ZERO);
// Satellite position and velocity in Station Frame
final Vector3D sat_pos = new Vector3D(0., 0., da-stationPositionEME.getNorm());
final Vector3D acceleration = new Vector3D(-context.gravity.getMu(), sat_pos);
final PVCoordinates pv_sat_topo = new PVCoordinates(sat_pos, geovelocity, acceleration);
// satellite position in EME2000
final PVCoordinates pv_sat_iner = topoToEME.transformPVCoordinates(pv_sat_topo);
// Geo-stationary Satellite Orbit, tightly above the station (l0-L0)
context.initialOrbit = new KeplerianOrbit(pv_sat_iner,
FramesFactory.getEME2000(),
new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc),
context.gravity.getMu());
context.stations = Arrays.asList(context.createStation(10.0, 45.0, 0.0, "Lat10_Long45") );
// Turn-around range stations
// Map entry = master station
// Map value = slave station associated
context.TARstations = new HashMap<GroundStation, GroundStation>();
context.TARstations.put(context.createStation( 41.977, 13.600, 671.354, "Fucino"),
context.createStation( 43.604, 1.444, 263.0 , "Toulouse"));
context.TARstations.put(context.createStation( 49.867, 8.65 , 144.0 , "Darmstadt"),
context.createStation( -25.885, 27.707, 1566.633, "Pretoria"));
return context;
}
public static Propagator createPropagator(final Orbit initialOrbit,
final PropagatorBuilder propagatorBuilder)
throws OrekitException {
// override orbital parameters
double[] orbitArray = new double[6];
propagatorBuilder.getOrbitType().mapOrbitToArray(initialOrbit,
propagatorBuilder.getPositionAngle(),
orbitArray);
for (int i = 0; i < orbitArray.length; ++i) {
propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(i).setValue(orbitArray[i]);
}
return propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
}
public static List<ObservedMeasurement<?>> createMeasurements(final Propagator propagator,
final MeasurementCreator creator,
final double startPeriod, final double endPeriod,
final double step)
throws OrekitException {
propagator.setMasterMode(step, creator);
final double period = propagator.getInitialState().getKeplerianPeriod();
final AbsoluteDate start = propagator.getInitialState().getDate().shiftedBy(startPeriod * period);
final AbsoluteDate end = propagator.getInitialState().getDate().shiftedBy(endPeriod * period);
propagator.propagate(start, end);
return creator.getMeasurements();
}
public static void checkFit(final Context context, final BatchLSEstimator estimator,
final int iterations, final int evaluations,
final double expectedRMS, final double rmsEps,
final double expectedMax, final double maxEps,
final double expectedDeltaPos, final double posEps,
final double expectedDeltaVel, final double velEps)
throws OrekitException {
final Orbit estimatedOrbit = estimator.estimate().getInitialState().getOrbit();
final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition();
final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity();
Assert.assertEquals(iterations, estimator.getIterationsCount());
Assert.assertEquals(evaluations, estimator.getEvaluationsCount());
Optimum optimum = estimator.getOptimum();
Assert.assertEquals(iterations, optimum.getIterations());
Assert.assertEquals(evaluations, optimum.getEvaluations());
int k = 0;
double sum = 0;
double max = 0;
for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry :
estimator.getLastEstimations().entrySet()) {
final ObservedMeasurement<?> m = entry.getKey();
final EstimatedMeasurement<?> e = entry.getValue();
final double[] weight = m.getBaseWeight();
final double[] sigma = m.getTheoreticalStandardDeviation();
final double[] observed = m.getObservedValue();
final double[] theoretical = e.getEstimatedValue();
for (int i = 0; i < m.getDimension(); ++i) {
final double weightedResidual = weight[i] * (theoretical[i] - observed[i]) / sigma[i];
++k;
sum += weightedResidual * weightedResidual;
max = FastMath.max(max, FastMath.abs(weightedResidual));
}
}
Assert.assertEquals(expectedRMS,
FastMath.sqrt(sum / k),
rmsEps);
Assert.assertEquals(expectedMax,
max,
maxEps);
Assert.assertEquals(expectedDeltaPos,
Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), estimatedPosition),
posEps);
Assert.assertEquals(expectedDeltaVel,
Vector3D.distance(context.initialOrbit.getPVCoordinates().getVelocity(), estimatedVelocity),
velEps);
}
}