/* 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); } } }