/* Copyright 2002-2019 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 fr.cs.examples;
import java.io.File;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.Locale;
import java.util.SortedSet;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.linear.MatrixUtils;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.ode.ODEIntegrator;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.random.CorrelatedRandomVectorGenerator;
import org.hipparchus.random.GaussianRandomGenerator;
import org.hipparchus.random.RandomGenerator;
import org.hipparchus.random.Well19937c;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.estimation.measurements.GroundStation;
import org.orekit.estimation.measurements.ObservableSatellite;
import org.orekit.estimation.measurements.ObservedMeasurement;
import org.orekit.estimation.measurements.RangeRate;
import org.orekit.estimation.measurements.generation.EventBasedScheduler;
import org.orekit.estimation.measurements.generation.Generator;
import org.orekit.estimation.measurements.generation.RangeRateBuilder;
import org.orekit.estimation.measurements.generation.Scheduler;
import org.orekit.estimation.measurements.generation.SignSemantic;
import org.orekit.estimation.measurements.modifiers.Bias;
import org.orekit.estimation.measurements.modifiers.RangeRateTroposphericDelayModifier;
import org.orekit.forces.drag.DragForce;
import org.orekit.forces.drag.IsotropicDrag;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.ThirdBodyAttraction;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.TopocentricFrame;
import org.orekit.models.AtmosphericRefractionModel;
import org.orekit.models.earth.EarthITU453AtmosphereRefraction;
import org.orekit.models.earth.atmosphere.Atmosphere;
import org.orekit.models.earth.atmosphere.DTM2000;
import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
import org.orekit.models.earth.troposphere.DiscreteTroposphericModel;
import org.orekit.models.earth.troposphere.ViennaModelCoefficientsLoader;
import org.orekit.models.earth.troposphere.ViennaModelType;
import org.orekit.models.earth.troposphere.ViennaThreeModel;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.OrbitType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.tle.TLE;
import org.orekit.propagation.analytical.tle.TLEPropagator;
import org.orekit.propagation.events.ElevationDetector;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.handlers.ContinueOnEvent;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateComponents;
import org.orekit.time.DatesSelector;
import org.orekit.time.FixedStepSelector;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.TimeStampedPVCoordinates;
/** Measurements generation for checking consistency of observed Doppler conversion.
* @author Luc Maisonobe
*/
public class MeasurementsGeneration {
/** Program entry point.
* @param args program arguments (unused here)
*/
public static void main(String[] args) {
try {
// configures Orekit
File home = new File(System.getProperty("user.home"));
File orekitData = new File(home, "orekit-data");
if (!orekitData.exists()) {
System.err.format(Locale.US, "Failed to find %s folder%n",
orekitData.getAbsolutePath());
System.err.format(Locale.US, "You need to download %s from %s, unzip it in %s and rename it 'orekit-data' for this tutorial to work%n",
"orekit-data-master.zip", "https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip",
home.getAbsolutePath());
System.exit(1);
}
DataProvidersManager manager = DataProvidersManager.getInstance();
manager.addProvider(new DirectoryCrawler(orekitData));
/*
// beware TLE model is not suitable for propagation more than a few days
// we have to select the correct TLE for the dates we will generate measurements
final TLE tle20190611 = new TLE("1 42778U 17036P 19162.69113096 .00000749 00000-0 35781-4 0 9997",
"2 42778 97.3549 220.5119 0010756 245.3704 114.6412 15.21997549109170");
final TLE tle20190612 = new TLE("1 42778U 17036P 19163.74305518 .00000728 00000-0 34875-4 0 9996",
"2 42778 97.3548 221.5408 0010807 241.2329 118.7822 15.21999132109334");
final TLE tle20190617 = new TLE("1 42778U 17036P 19168.73968312 .00000624 00000-0 30313-4 0 9995",
"2 42778 97.3544 226.4276 0011195 221.3881 138.6508 15.22005879110094");
*/
// Imports the International Terrestrial Reference Frame (ITRF) which is defined by the IERS and
// defines the geodesic model of the Earth based on the ITRF files
final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
itrf); //
// Input to define the location of the ground station and it's name
final GeodeticPoint location = new GeodeticPoint(FastMath.toRadians(52.834),
FastMath.toRadians(6.379),
10.0);
final GroundStation station = new GroundStation(new TopocentricFrame(earth, location, "39-CGBSAT-VHF"));
// measurements generation parameters
final double minElevation = FastMath.toRadians(1.0);
final double timeStep = 5.0;
final TimeScale utc = TimeScalesFactory.getUTC();
final AbsoluteDate t0 = new AbsoluteDate("2019-06-11T00:00:00.000", utc);
final double duration = Constants.JULIAN_DAY; // Sets the duration to be the length of a Julian day
//*
// atmosphere model for drag
MarshallSolarActivityFutureEstimation msafe =
new MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
manager.feed(msafe.getSupportedNames(), msafe);
Atmosphere atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
// this orbit is a dummy one, close to the first results from orbit determination from Max Valier satellite
final NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(12, 12);
final AbsoluteDate tOrb = new AbsoluteDate("2019-06-11T16:35:29.149", utc);
final TimeStampedPVCoordinates pvt = new TimeStampedPVCoordinates(tOrb,
new Vector3D(-5253194.0, -4400295.0, 80075.0),
new Vector3D(-559.0, 764.0, 7585.0));
final CartesianOrbit orbit = new CartesianOrbit(pvt, FramesFactory.getEME2000(), gravity.getMu());
final OrbitType type = OrbitType.CARTESIAN;
// set up numerical propagator
final double[][] tol = NumericalPropagator.tolerances(10.0, orbit, type);
final ODEIntegrator integrator = new DormandPrince853Integrator(0.001, 300.0, tol[0], tol[1]);
final NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(type);
propagator.setInitialState(new SpacecraftState(orbit, 16.0)); // the second argument is the spacecraft mass
// add a few realistic force models
final double cd = 2.0;
final double area = 0.25;
propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, gravity));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
propagator.addForceModel(new DragForce(atmosphere, new IsotropicDrag(cd, area)));
//*/
// set up some correction models
// coefficients files for Vienna Model 3 can be found at
// http://vmf.geo.tuwien.ac.at/trop_products/GRID/1x1/VMF3/VMF3_OP/2019/
final ViennaModelCoefficientsLoader loader = new ViennaModelCoefficientsLoader(location.getLatitude(),
location.getLongitude(),
ViennaModelType.VIENNA_THREE);
loader.loadViennaCoefficients(t0.getComponents(utc));
final DiscreteTroposphericModel tropo = new ViennaThreeModel(loader.getA(),
loader.getZenithDelay(),
location.getLatitude(),
location.getLongitude());
final AtmosphericRefractionModel refraction = new EarthITU453AtmosphereRefraction(location.getAltitude());
// set up measurements generation, with realistic models
final Generator generator = new Generator();
//final ObservableSatellite os = generator.addPropagator(TLEPropagator.selectExtrapolator(tle20190611));
final ObservableSatellite os = generator.addPropagator(propagator);
/*
Range rate builder
*/
final double sigma = 15.0; // the value we will tune
final RealMatrix covariance = MatrixUtils.createRealDiagonalMatrix(new double[] { sigma * sigma });
final RandomGenerator random = new Well19937c(0x9e9409e2520c1fc3l); // you can change this seed as you like, it's just a seed
final CorrelatedRandomVectorGenerator crvg =
new CorrelatedRandomVectorGenerator(covariance,
1.0e-10,
new GaussianRandomGenerator(random));
final RangeRateBuilder builder = new RangeRateBuilder(crvg, station, false, sigma, 1.0, os);//*/
final double rangeRateBias = 375.0;
builder.addModifier(new Bias<RangeRate>(new String[] { station.getBaseFrame().getName() + "-range-rate-bias" },
new double[] { rangeRateBias },
new double[] { 1.0 },
new double[] { Double.NEGATIVE_INFINITY },
new double[] { Double.POSITIVE_INFINITY }));
builder.addModifier(new RangeRateTroposphericDelayModifier(tropo, false));
final DatesSelector selector = new FixedStepSelector(timeStep, utc);
final EventDetector detector = new ElevationDetector(station.getBaseFrame()).
withConstantElevation(minElevation).
withRefraction(refraction).
withHandler(new ContinueOnEvent<>());
final Scheduler<RangeRate> scheduler = new EventBasedScheduler<>(builder, selector, generator.getPropagator(os), detector,
SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE);
generator.addScheduler(scheduler);
// generate measurements
final SortedSet<ObservedMeasurement<?>> measurements = generator.generate(t0, t0.shiftedBy(duration));
System.out.println("generated " + measurements.size() + " measurements");
System.out.println("first measurement at " + measurements.first().getDate());
System.out.println("last measurement at " + measurements.last().getDate());
// save measurement in a file
AbsoluteDate mjdRefUTC = new AbsoluteDate(DateComponents.MODIFIED_JULIAN_EPOCH, utc);
try (PrintWriter out = new PrintWriter(new File(home, "generated-doppler-"+station.getBaseFrame().getName()+".dat"))) {
measurements.
stream().
forEach(m -> out.format(Locale.US, "%s %s %s %12.6f%n", // Without MJD coloumn
m.getDate(),
"RANGE_RATE",
station.getBaseFrame().getName(),
(m.getObservedValue()[0]/1000)));
// Uncomment to get mjd column
/*forEach(m -> out.format(Locale.US, "%s %14.8f %s %s %12.6f%n",
m.getDate(),
m.getDate().offsetFrom(mjdRefUTC, utc) / Constants.JULIAN_DAY,
"RANGE-RATE",
station.getBaseFrame().getName(),
m.getObservedValue()[0]));*/
}
} catch (IOException ioe) {
System.err.println(ioe.getLocalizedMessage());
System.exit(1);
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
}
}
}