Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import java.io.File;
- import java.util.ArrayList;
- import java.util.List;
- import java.util.Random;
- import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer;
- import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.Decomposition;
- import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
- 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.leastsquares.BatchLSEstimator;
- import org.orekit.estimation.measurements.AngularAzEl;
- import org.orekit.estimation.measurements.GroundStation;
- import org.orekit.estimation.measurements.ObservedMeasurement;
- import org.orekit.estimation.measurements.Range;
- import org.orekit.forces.drag.DragForce;
- import org.orekit.forces.drag.DragSensitive;
- import org.orekit.forces.drag.IsotropicDrag;
- import org.orekit.forces.drag.atmosphere.DTM2000;
- import org.orekit.forces.drag.atmosphere.data.MarshallSolarActivityFutureEstimation;
- import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
- 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.gravity.potential.TideSystem;
- import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
- import org.orekit.forces.radiation.SolarRadiationPressure;
- import org.orekit.frames.Frame;
- import org.orekit.frames.FramesFactory;
- import org.orekit.frames.TopocentricFrame;
- import org.orekit.orbits.CartesianOrbit;
- import org.orekit.orbits.KeplerianOrbit;
- import org.orekit.orbits.Orbit;
- import org.orekit.orbits.PositionAngle;
- import org.orekit.propagation.SpacecraftState;
- import org.orekit.propagation.conversion.DormandPrince54IntegratorBuilder;
- import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
- import org.orekit.propagation.numerical.NumericalPropagator;
- import org.orekit.time.AbsoluteDate;
- import org.orekit.time.TimeScalesFactory;
- import org.orekit.utils.Constants;
- import org.orekit.utils.IERSConventions;
- import org.orekit.utils.PVCoordinatesProvider;
- import org.orekit.utils.ParameterDriversList.DelegatingDriver;
- import org.orekit.utils.TimeStampedPVCoordinates;
- public class main {
- static double mu = 3.986004415e+14;
- static AbsoluteDate mjd2000;
- public static void main(String[] args) throws Exception {
- File orekitData = new File("..\\orekit-data");
- DataProvidersManager manager = DataProvidersManager.getInstance();
- manager.addProvider(new DirectoryCrawler(orekitData));
- mjd2000 = new AbsoluteDate(2000, 1, 1, 0, 0, 0.0, TimeScalesFactory.getTT());
- orbDeterminationSimulatedObservations();
- }
- public static NumericalPropagatorBuilder getPropagator(Orbit initialOrbit) throws OrekitException {
- // Default area and mass.
- return getPropagator(initialOrbit, 50.0, 0.5, 2.2, 0.5, 1.1);
- }
- public static NumericalPropagatorBuilder getPropagator(Orbit initialOrbit, double mass, double drgarea,
- double drgcoef, double srparea, double srpcoef) throws OrekitException {
- boolean includeAtmosphere = true;
- boolean includeSolar = true;
- boolean include3rd = true;
- boolean includeRel = true;
- boolean includeTides = true;
- boolean includeGeopotential = true;
- // steps limits
- final double minStep = 0.001;
- final double maxStep = 60;
- // set up space dynamics propagator
- NumericalPropagatorBuilder propbuilder = new NumericalPropagatorBuilder(initialOrbit,
- new DormandPrince54IntegratorBuilder(minStep, maxStep, 1), PositionAngle.TRUE, 1);
- propbuilder.setMass(mass);
- // Force model inclusion
- PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
- Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
- OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
- Constants.WGS84_EARTH_FLATTENING, itrf);
- // * Atmosphere DTM2000
- if (includeAtmosphere) {
- MarshallSolarActivityFutureEstimation inMarshall = new MarshallSolarActivityFutureEstimation(
- "(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)",
- MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
- earth.setAngularThreshold(1e-10);
- DTM2000 atm = new DTM2000(inMarshall, sun, earth);
- DragSensitive dragsensitive = new IsotropicDrag(drgarea, drgcoef);
- propbuilder.addForceModel(new DragForce(atm, dragsensitive));
- }
- // * Relativity
- if (includeRel) {
- propbuilder.addForceModel(new Relativity(mu));
- }
- // * Solar radiation pressure
- if (includeSolar) {
- SolarRadiationPressure SRP = new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
- new IsotropicRadiationSingleCoefficient(srparea, srpcoef));
- propbuilder.addForceModel(SRP);
- }
- // * Third body
- if (include3rd) {
- propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
- propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
- propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getJupiter()));
- propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getVenus()));
- propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMars()));
- }
- // * Solid Tide
- if (includeTides) {
- propbuilder.addForceModel(new SolidTides(itrf, 6378136.460, mu, TideSystem.TIDE_FREE, true, 60., 2,
- IERSConventions.IERS_2010, TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true),
- CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon()));
- }
- // * Geopotential
- if (includeGeopotential) {
- NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(16, 16);
- propbuilder.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, gravity));
- }
- return propbuilder;
- }
- // orbDeterminationSimulatedObservations propagates an orbit for a time, then
- // takes measurements from a simulated radar and tries to perform OD on the
- // observations.
- public static void orbDeterminationSimulatedObservations() throws Exception {
- double rangeSigma = 100; // m
- double rangeWeight = 1.0;
- double[] azelSigma = { 1.0 * Math.PI / 180., 1.0 * Math.PI / 180. }; // rad
- double[] azelWeight = { 1.0, 1.0 };
- double dopplerSigma = 0.5; // m/s
- double dopplerWeight = 1.0;
- double altitude = 250; // altitude of the circular orbit in km.
- double timeSpan = 3; // days to propagate
- double minElev = 10 * Math.PI / 180.; // minimum elevation in rad.
- boolean useRangeMeas = false;
- boolean useAzelMeas = true;
- boolean useDopplerMeas = false;
- double step = 60; // seconds
- int simulatedMeasurements = (int) Math.ceil(timeSpan * 86400. / step);
- GroundStation gs = getRadarStation();
- Random rnd = new Random();
- // initial state vector.
- AbsoluteDate initialDate = mjd2000;
- Orbit initialOrbitCart = new CartesianOrbit(new KeplerianOrbit(6378000 + altitude * 1000, 0.0,
- 50 * Math.PI / 180, 0.0, 0.0, 0.0, PositionAngle.TRUE, FramesFactory.getEME2000(), initialDate, mu));
- NumericalPropagatorBuilder propbuilder = getPropagator(initialOrbitCart);
- SpacecraftState currentState = null;
- System.out.println("Initial orbit:");
- printOrbit(initialOrbitCart);
- double[] normalizedParameters = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
- NumericalPropagator prop = propbuilder.buildPropagator(normalizedParameters);
- // Propagation
- System.out.println("Propagating...");
- List<ObservedMeasurement> meas = new ArrayList<>();
- List<TimeStampedPVCoordinates> pv = new ArrayList<>();
- for (int i = 0; i < simulatedMeasurements; i++) {
- AbsoluteDate nextdate = new AbsoluteDate(initialDate, step * i);
- currentState = prop.propagate(nextdate);
- TimeStampedPVCoordinates pos = currentState.getPVCoordinates();
- pv.add(pos);
- // Check if visible from radar.
- double elevation = gs.getBaseFrame().getElevation(pos.getPosition(), currentState.getFrame(),
- currentState.getDate());
- if (elevation < minElev) {
- continue;
- }
- if (useAzelMeas) {
- // Get azel observation
- double azimuth = gs.getBaseFrame().getAzimuth(pos.getPosition(), currentState.getFrame(),
- currentState.getDate());
- double[] azel = { azimuth + rnd.nextGaussian() * azelSigma[0],
- elevation + rnd.nextGaussian() * azelSigma[1] };
- // Add noisy AzEl measurement.
- ObservedMeasurement azelMeas = new AngularAzEl(gs, currentState.getDate(), azel, azelSigma, azelWeight);
- meas.add(azelMeas);
- }
- if (useRangeMeas) {
- // Get Range observation.
- double range = gs.getBaseFrame().getRange(pos.getPosition(), currentState.getFrame(),
- currentState.getDate());
- range += rnd.nextGaussian() * rangeSigma;
- ObservedMeasurement rangeMeas = new Range(gs, currentState.getDate(), range, rangeSigma, rangeWeight);
- meas.add(rangeMeas);
- }
- if (useDopplerMeas) {
- // Get Doppler observation.
- double doppler = gs.getBaseFrame().getRangeRate(pos, currentState.getFrame(), currentState.getDate());
- doppler += rnd.nextGaussian() * dopplerSigma;
- ObservedMeasurement dopplerMeas = new Range(gs, currentState.getDate(), doppler, dopplerSigma,
- dopplerWeight);
- meas.add(dopplerMeas);
- }
- }
- System.out.println("Number of measurements: " + meas.size());
- // Orbit determination.
- System.out.println("Running the orbit determination...");
- GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(Decomposition.QR);
- NumericalPropagatorBuilder propbuilder2 = getPropagator(new CartesianOrbit(initialOrbitCart));
- BatchLSEstimator estimator = new BatchLSEstimator(optimizer, propbuilder2);
- estimator.setObserver(new LeastSquareObserver());
- estimator.setMaxEvaluations(20);
- estimator.setMaxIterations(10);
- estimator.setParametersConvergenceThreshold(0.01);
- // Add measurements
- for (ObservedMeasurement oneMeas : meas) {
- estimator.addMeasurement(oneMeas);
- }
- NumericalPropagator[] neworb = estimator.estimate();
- System.out.println("Determined Orbit");
- printOrbit(neworb[0].getInitialState());
- Optimum optimum = estimator.getOptimum();
- System.out.println(
- "Optimum: " + optimum.getEvaluations() + " " + optimum.getIterations() + " " + optimum.getRMS());
- System.out.println("Residuals: " + optimum.getResiduals());
- System.out.println("Covariance matrix: " + optimum.getCovariances(0.0));
- for (DelegatingDriver param : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
- System.out.println("Param " + param.getName() + " is " + param.getValue());
- }
- }
- public static GroundStation getRadarStation() throws OrekitException {
- OneAxisEllipsoid refbody2 = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
- Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, false));
- GeodeticPoint stationCoord = new GeodeticPoint(40.41 / 180. * Math.PI, -3.4 / 180. * Math.PI, 700.0);
- TopocentricFrame frame = new TopocentricFrame(refbody2, stationCoord, "Radar station");
- GroundStation gs = new GroundStation(frame);
- return gs;
- }
- public static void printOrbit(SpacecraftState sc) throws OrekitException {
- printOrbit(sc.getOrbit());
- }
- public static void printOrbit(Orbit orbit) throws OrekitException {
- System.out.println("***");
- System.out.println(" Date: " + orbit.getDate().toString(TimeScalesFactory.getTT()));
- TimeStampedPVCoordinates pos = orbit.getPVCoordinates(FramesFactory.getEME2000());
- System.out.println(" Position: " + pos.getPosition());
- System.out.println(" Velocity: " + pos.getVelocity());
- System.out.println(" a: " + orbit.getA() + " / e: " + orbit.getE() + " / i: " + orbit.getI());
- }
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement