Commit 9f559b95 authored by Bryan Cazabonne's avatar Bryan Cazabonne

Added orbit determination tutorial using SLR measurements.

parent fe5873c3
/* Copyright 2002-2020 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.tutorials.estimation;
import java.io.File;
import java.io.IOException;
import java.net.URISyntaxException;
import java.util.NoSuchElementException;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.bodies.CelestialBody;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
import org.orekit.forces.ForceModel;
import org.orekit.forces.drag.DragForce;
import org.orekit.forces.drag.DragSensitive;
import org.orekit.forces.empirical.AccelerationModel;
import org.orekit.forces.empirical.ParametricAcceleration;
import org.orekit.forces.empirical.PolynomialAccelerationModel;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.OceanTides;
import org.orekit.forces.gravity.Relativity;
import org.orekit.forces.gravity.SolidTides;
import org.orekit.forces.gravity.ThirdBodyAttraction;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.forces.radiation.RadiationSensitive;
import org.orekit.forces.radiation.SolarRadiationPressure;
import org.orekit.models.earth.atmosphere.Atmosphere;
import org.orekit.models.earth.atmosphere.NRLMSISE00;
import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherData;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
import org.orekit.propagation.conversion.ODEIntegratorBuilder;
import org.orekit.time.TimeScalesFactory;
import org.orekit.tutorials.estimation.common.AbstractOrbitDetermination;
import org.orekit.tutorials.estimation.common.TutorialOrbitDetermination;
import org.orekit.tutorials.yaml.TutorialForceModel.TutorialGravity;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.ParameterDriver;
/** Orekit tutorial for orbit determination based on laser ranging data.
* @author Bryan Cazabonne
*/
public class LaserRangingOrbitDetermination extends AbstractOrbitDetermination<NumericalPropagatorBuilder> {
/** Gravity field. */
private NormalizedSphericalHarmonicsProvider gravityField;
/** Earth body. */
private OneAxisEllipsoid earth;
/** Program entry point.
* @param args program arguments (unused here)
*/
public static void main(final String[] args) {
try {
// input in tutorial resources directory
final String inputPath = NumericalOrbitDetermination.class.getClassLoader().
getResource("slr-od/slr-orbit-determination.yaml").toURI().getPath();
new LaserRangingOrbitDetermination().run(new File(inputPath));
} catch (URISyntaxException | IOException | OrekitException e) {
System.err.println(e.getLocalizedMessage());
System.exit(1);
}
}
/** {@inheritDoc} */
@Override
protected void createGravityField(final TutorialOrbitDetermination inputData)
throws NoSuchElementException {
final TutorialGravity gravity = inputData.getPropagator().getForceModels().getGravity();
final int degree = gravity.getDegree();
final int order = FastMath.min(degree, gravity.getOrder());
gravityField = GravityFieldFactory.getNormalizedProvider(degree, order);
}
/** {@inheritDoc} */
@Override
protected double getMu() {
return gravityField.getMu();
}
/** {@inheritDoc} */
@Override
protected NumericalPropagatorBuilder createPropagatorBuilder(final Orbit referenceOrbit,
final ODEIntegratorBuilder builder,
final double positionScale) {
return new NumericalPropagatorBuilder(referenceOrbit, builder, PositionAngle.MEAN,
positionScale);
}
/** {@inheritDoc} */
@Override
protected void setMass(final NumericalPropagatorBuilder propagatorBuilder,
final double mass) {
propagatorBuilder.setMass(mass);
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setGravity(final NumericalPropagatorBuilder propagatorBuilder,
final OneAxisEllipsoid body) {
earth = body;
final ForceModel gravityModel = new HolmesFeatherstoneAttractionModel(body.getBodyFrame(), gravityField);
propagatorBuilder.addForceModel(gravityModel);
return gravityModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setOceanTides(final NumericalPropagatorBuilder propagatorBuilder,
final IERSConventions conventions,
final OneAxisEllipsoid body,
final int degree, final int order) {
final ForceModel tidesModel = new OceanTides(body.getBodyFrame(),
gravityField.getAe(), gravityField.getMu(),
degree, order, conventions,
TimeScalesFactory.getUT1(conventions, true));
propagatorBuilder.addForceModel(tidesModel);
return tidesModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setSolidTides(final NumericalPropagatorBuilder propagatorBuilder,
final IERSConventions conventions,
final OneAxisEllipsoid body,
final CelestialBody[] solidTidesBodies) {
final ForceModel tidesModel = new SolidTides(body.getBodyFrame(),
gravityField.getAe(), gravityField.getMu(),
gravityField.getTideSystem(), conventions,
TimeScalesFactory.getUT1(conventions, true),
solidTidesBodies);
propagatorBuilder.addForceModel(tidesModel);
return tidesModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setThirdBody(final NumericalPropagatorBuilder propagatorBuilder,
final CelestialBody thirdBody) {
final ForceModel thirdBodyModel = new ThirdBodyAttraction(thirdBody);
propagatorBuilder.addForceModel(thirdBodyModel);
return thirdBodyModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setDrag(final NumericalPropagatorBuilder propagatorBuilder,
final Atmosphere atmosphere, final DragSensitive spacecraft) {
final Atmosphere atm = new NRLMSISE00(new CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt"), DataContext.getDefault().getCelestialBodies().getSun(), earth);
final ForceModel dragModel = new DragForce(atm, spacecraft);
propagatorBuilder.addForceModel(dragModel);
return dragModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setSolarRadiationPressure(final NumericalPropagatorBuilder propagatorBuilder, final CelestialBody sun,
final double equatorialRadius, final RadiationSensitive spacecraft) {
final ForceModel srpModel = new SolarRadiationPressure(sun, equatorialRadius, spacecraft);
propagatorBuilder.addForceModel(srpModel);
return srpModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setRelativity(final NumericalPropagatorBuilder propagatorBuilder) {
final ForceModel relativityModel = new Relativity(gravityField.getMu());
propagatorBuilder.addForceModel(relativityModel);
return relativityModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected ParameterDriver[] setPolynomialAcceleration(final NumericalPropagatorBuilder propagatorBuilder,
final String name, final Vector3D direction, final int degree) {
final AccelerationModel accModel = new PolynomialAccelerationModel(name, null, degree);
final ForceModel polynomialModel = new ParametricAcceleration(direction, true, accModel); propagatorBuilder.addForceModel(polynomialModel);
return polynomialModel.getParametersDrivers();
}
/** {@inheritDoc} */
@Override
protected void setAttitudeProvider(final NumericalPropagatorBuilder propagatorBuilder,
final AttitudeProvider attitudeProvider) {
propagatorBuilder.setAttitudeProvider(attitudeProvider);
}
/** {@inheritDoc} */
@Override
protected void compareWithReference(final Orbit estimatedOrbit) {
// Reference results
final Vector3D refPos = new Vector3D(7526994.072, -9646309.832, 1464110.239);
final Vector3D refVel = new Vector3D(3033.794, 1715.265, -4447.659);
// Print results
System.out.println("Comparison with reference orbit: ");
System.out.println("ΔP(m) = " + Vector3D.distance(refPos, estimatedOrbit.getPVCoordinates().getPosition()));
System.out.println("ΔV(m/s) = " + Vector3D.distance(refVel, estimatedOrbit.getPVCoordinates().getVelocity()));
}
}
......@@ -30,6 +30,9 @@ public class TutorialMeasurements {
/** List of measurement files. */
private List<String> measurementFiles;
/** Sinex files containing station positions and eccentricities. */
private TutorialSinex sinex;
/** Flag for station position estimation. */
private boolean withStationPositionEstimated;
......@@ -85,6 +88,23 @@ public class TutorialMeasurements {
this.measurementFiles = measurementFiles;
}
/**
* Get the sinex files containing station positions and eccentricities.
* @return sinex files
*/
public TutorialSinex getSinex() {
return sinex;
}
/**
* Set the sinex files containing station positions and eccentricities.
* @param sinex the sinex file object to set
*/
public void setSinex(final TutorialSinex sinex) {
this.sinex = sinex;
}
/**
* Get the flag for station position estimation.
* @return true if station position is estimated
......
/* Copyright 2002-2020 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.tutorials.yaml;
/**
* Container for SINEX data.
* @author Bryan Cazabonne
*/
public class TutorialSinex {
/** Path for sinex station positions file. */
private String stationPositions;
/** Path for sinex station eccentricities file. */
private String stationEccentricities;
/**
* Get the path for sinex station positions file.
* @return path for sinex station positions file
*/
public String getStationPositions() {
return stationPositions;
}
/**
* Set the path for sinex station positions file.
* @param stationPositions the path to set
*/
public void setStationPositions(final String stationPositions) {
this.stationPositions = stationPositions;
}
/**
* Get the path for sinex station eccentricities file.
* @return the path for sinex station eccentricities file.
*/
public String getStationEccentricities() {
return stationEccentricities;
}
/**
* Set the path for sinex station eccentricities file.
* @param stationEccentricities the path to set
*/
public void setStationEccentricities(final String stationEccentricities) {
this.stationEccentricities = stationEccentricities;
}
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
# input file for laser ranging orbit determination
# Lageos 2 orbit determination using laser ranging data
# Laser measurements are read from CRD files provided by ILRS
# Lines starting with '#' are silently ignored
# Estimated parameters (clock offsets, atmospheric parameters, etc.)
# are given by [initialValue, minValue, maxValue]
# Orbit definition (position in meters, velocity in meter per seconds, angles in degrees, date in ISO 8601)
orbit:
date: "2016-02-13T16:00:00.000"
frameName: "EME2000"
orbitType:
name: "CARTESIAN"
cartesian:
x: 7526990.0
y: -9646310.0
z: 1464110.0
vx: 3033.0
vy: 1715.0
vz: -4447.0
# Spacecraft definition (mass in kilogrammes)
spacecraft:
mass: 405.380
bias:
# !!! Here the satellite bias corresponds to the center of mass offset !!!
values: [0.251, -1000.0, 1000.0]
isEstimated: false
# Propagator definition
propagator:
# Numerical integrator (min step (s), max step (s) and position error (m))
integrator:
minStep: 0.001
maxStep: 300
positionError: 10.0
# Force models used by the propagator (only the ones used !!!)
forceModels:
# Central body gravity
gravity:
degree: 30
order: 30
# 3rd body attraction
thirdBody:
- name: "Sun"
withSolidTides: true
- name: "Moon"
withSolidTides: true
- name: "Venus"
withSolidTides: false
# Post-Newtonian correction force due to general relativity ("isUsed: true" to not have an empty value)
relativity:
isUsed: true
# Oean tides:
oceanTides:
degree: 6
order: 6
# Polynomial acceleration
polynomialAcceleration:
- name: "Z-bias"
directions: [0.0, 0.0, 1.0]
coefficients: [0.0]
isEstimated: true
# Body definition (According to ILRS recommendations, the computation models follow the IERS 2003 Conventions as closely as possible)
body:
iersConventionYear: 2003
frameName: "CIO/2003-based ITRF accurate EOP"
equatorialRadius: 6378137.0
inverseFlattening: 298.257223563
# Measurements definition
measurements:
measurementFiles: ["lageos2_20160214.npt"]
sinex:
stationPositions: "SLRF2014_POS+VEL_2030.0_200428.snx"
stationEccentricities: "ecc_une.snx"
range:
sigma: 20.0
weight: 1.0
outlierRejectionMultiplier: 6
outlierRejectionStartingIteration: 2
withStationPositionEstimated: false
# Troposphere: According to ILRS recommendations, Mendes-Pavlis model must be used
# Therefore, no need to configure the model here since it is done directly in the tutorial
# Shapiro correction on range measurements
withShapiro: true
# correction of ground stations displacements if remove_permanent_deformation is true,
# the station coordinates are considered *mean tide* and already include the permanent
# deformation, hence it should be removed from the displacement to avoid considering it
# twice if false, the station coordinates are considered *conventional tide free*
# so the permanent deformation must be included in the displacement
stationDisplacement:
withTidalCorrection: true
removePermanentDeformation: false
withOceanLoadingCorrection: false
# Stations
# Station names correspond to their identifier in Sinex and CRD files
# Station coordinates are read in a SINEX file
stations:
- name: "7090"
rangeBias:
values: [0.0, -1000.0, 1000.0]
isEstimated: true
observationTypes:
- name: ""
clockOffset:
values: [0.0, -0.001, 0.001]
isEstimated: false
- name: "7825"
rangeBias:
values: [0.0, -1000.0, 1000.0]
isEstimated: true
observationTypes:
- name: ""
clockOffset:
values: [0.0, -0.001, 0.001]
isEstimated: false
- name: "7119"
rangeBias:
values: [0.0, -1000.0, 1000.0]
isEstimated: true
observationTypes:
- name: ""
clockOffset:
values: [0.0, -0.001, 0.001]
isEstimated: false
- name: "7941"
rangeBias:
values: [0.0, -1000.0, 1000.0]
isEstimated: true
observationTypes:
- name: ""
clockOffset:
values: [0.0, -0.001, 0.001]
isEstimated: false
# Estimator definition
estimator:
# normalized parameters p are computed from physical parameters x as p = (x - x0) / sx
# where x0 is the reference value for physical parameter x and sx is the associated
# scaling factor for orbital parameters normalization (m)
orbitalParametersPositionScale: 100.0
# Levenberg-Marquardt or a Gauss-Newton
optimizationEngine:
# levenbergMarquardt or gaussNewton
levenbergMarquardt:
initialStep: 1.0e6
maxIterations: 20
maxEvaluations: 25
# convergence is reached when max|p(k+1) - p(k)| < ε for each normalized estimated
# parameters p and iterations k and k+1 so the ε threshold
# Normalized values are computed as (x - x0) / sx, so convergence is reached when
# the following condition holds for all estimated parameters: |x[i+1] - x[i]| <= ε * sx
# So the convergence threshold specified here can be considered as a multiplication
# factor applied to scale. Since for all parameters the scale is often small,
# then the threshold should not be too small. A value of 10⁻³ is often quite accurate.
convergenceThreshold: 1.0e-3
# base name of the output files (log and residuals), no files created if empty
outputBaseName: ""
\ No newline at end of file
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