API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Nov 15th, 2017
51
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Java 11.95 KB | None | 0 0
raw download clone embed print report
  1. import java.io.File;
  2. import java.util.ArrayList;
  3. import java.util.List;
  4. import java.util.Random;
  5.  
  6. import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer;
  7. import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.Decomposition;
  8. import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
  9. import org.orekit.bodies.CelestialBodyFactory;
  10. import org.orekit.bodies.GeodeticPoint;
  11. import org.orekit.bodies.OneAxisEllipsoid;
  12. import org.orekit.data.DataProvidersManager;
  13. import org.orekit.data.DirectoryCrawler;
  14. import org.orekit.errors.OrekitException;
  15. import org.orekit.estimation.leastsquares.BatchLSEstimator;
  16. import org.orekit.estimation.measurements.AngularAzEl;
  17. import org.orekit.estimation.measurements.GroundStation;
  18. import org.orekit.estimation.measurements.ObservedMeasurement;
  19. import org.orekit.estimation.measurements.Range;
  20. import org.orekit.forces.drag.DragForce;
  21. import org.orekit.forces.drag.DragSensitive;
  22. import org.orekit.forces.drag.IsotropicDrag;
  23. import org.orekit.forces.drag.atmosphere.DTM2000;
  24. import org.orekit.forces.drag.atmosphere.data.MarshallSolarActivityFutureEstimation;
  25. import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
  26. import org.orekit.forces.gravity.Relativity;
  27. import org.orekit.forces.gravity.SolidTides;
  28. import org.orekit.forces.gravity.ThirdBodyAttraction;
  29. import org.orekit.forces.gravity.potential.GravityFieldFactory;
  30. import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
  31. import org.orekit.forces.gravity.potential.TideSystem;
  32. import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
  33. import org.orekit.forces.radiation.SolarRadiationPressure;
  34. import org.orekit.frames.Frame;
  35. import org.orekit.frames.FramesFactory;
  36. import org.orekit.frames.TopocentricFrame;
  37. import org.orekit.orbits.CartesianOrbit;
  38. import org.orekit.orbits.KeplerianOrbit;
  39. import org.orekit.orbits.Orbit;
  40. import org.orekit.orbits.PositionAngle;
  41. import org.orekit.propagation.SpacecraftState;
  42. import org.orekit.propagation.conversion.DormandPrince54IntegratorBuilder;
  43. import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
  44. import org.orekit.propagation.numerical.NumericalPropagator;
  45. import org.orekit.time.AbsoluteDate;
  46. import org.orekit.time.TimeScalesFactory;
  47. import org.orekit.utils.Constants;
  48. import org.orekit.utils.IERSConventions;
  49. import org.orekit.utils.PVCoordinatesProvider;
  50. import org.orekit.utils.ParameterDriversList.DelegatingDriver;
  51. import org.orekit.utils.TimeStampedPVCoordinates;
  52.  
  53. public class main {
  54.     static double mu = 3.986004415e+14;
  55.     static AbsoluteDate mjd2000;
  56.  
  57.     public static void main(String[] args) throws Exception {
  58.         File orekitData = new File("..\\orekit-data");
  59.         DataProvidersManager manager = DataProvidersManager.getInstance();
  60.         manager.addProvider(new DirectoryCrawler(orekitData));
  61.         mjd2000 = new AbsoluteDate(2000, 1, 1, 0, 0, 0.0, TimeScalesFactory.getTT());
  62.  
  63.         orbDeterminationSimulatedObservations();
  64.     }
  65.  
  66.     public static NumericalPropagatorBuilder getPropagator(Orbit initialOrbit) throws OrekitException {
  67.         // Default area and mass.
  68.         return getPropagator(initialOrbit, 50.0, 0.5, 2.2, 0.5, 1.1);
  69.     }
  70.  
  71.     public static NumericalPropagatorBuilder getPropagator(Orbit initialOrbit, double mass, double drgarea,
  72.             double drgcoef, double srparea, double srpcoef) throws OrekitException {
  73.         boolean includeAtmosphere = true;
  74.         boolean includeSolar = true;
  75.         boolean include3rd = true;
  76.         boolean includeRel = true;
  77.         boolean includeTides = true;
  78.         boolean includeGeopotential = true;
  79.  
  80.         // steps limits
  81.         final double minStep = 0.001;
  82.         final double maxStep = 60;
  83.  
  84.         // set up space dynamics propagator
  85.         NumericalPropagatorBuilder propbuilder = new NumericalPropagatorBuilder(initialOrbit,
  86.                 new DormandPrince54IntegratorBuilder(minStep, maxStep, 1), PositionAngle.TRUE, 1);
  87.         propbuilder.setMass(mass);
  88.         // Force model inclusion
  89.         PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
  90.         Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
  91.         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
  92.                 Constants.WGS84_EARTH_FLATTENING, itrf);
  93.         // * Atmosphere DTM2000
  94.         if (includeAtmosphere) {
  95.             MarshallSolarActivityFutureEstimation inMarshall = new MarshallSolarActivityFutureEstimation(
  96.                     "(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)",
  97.                     MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
  98.             earth.setAngularThreshold(1e-10);
  99.             DTM2000 atm = new DTM2000(inMarshall, sun, earth);
  100.             DragSensitive dragsensitive = new IsotropicDrag(drgarea, drgcoef);
  101.             propbuilder.addForceModel(new DragForce(atm, dragsensitive));
  102.         }
  103.  
  104.         // * Relativity
  105.         if (includeRel) {
  106.             propbuilder.addForceModel(new Relativity(mu));
  107.         }
  108.  
  109.         // * Solar radiation pressure
  110.         if (includeSolar) {
  111.             SolarRadiationPressure SRP = new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
  112.                     new IsotropicRadiationSingleCoefficient(srparea, srpcoef));
  113.             propbuilder.addForceModel(SRP);
  114.         }
  115.  
  116.         // * Third body
  117.         if (include3rd) {
  118.             propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
  119.             propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
  120.             propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getJupiter()));
  121.             propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getVenus()));
  122.             propbuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMars()));
  123.         }
  124.         // * Solid Tide
  125.         if (includeTides) {
  126.             propbuilder.addForceModel(new SolidTides(itrf, 6378136.460, mu, TideSystem.TIDE_FREE, true, 60., 2,
  127.                     IERSConventions.IERS_2010, TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true),
  128.                     CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon()));
  129.         }
  130.  
  131.         // * Geopotential
  132.         if (includeGeopotential) {
  133.             NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(16, 16);
  134.             propbuilder.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, gravity));
  135.         }
  136.         return propbuilder;
  137.     }
  138.  
  139.     // orbDeterminationSimulatedObservations propagates an orbit for a time, then
  140.     // takes measurements from a simulated radar and tries to perform OD on the
  141.     // observations.
  142.     public static void orbDeterminationSimulatedObservations() throws Exception {
  143.         double rangeSigma = 100; // m
  144.         double rangeWeight = 1.0;
  145.         double[] azelSigma = { 1.0 * Math.PI / 180., 1.0 * Math.PI / 180. }; // rad
  146.         double[] azelWeight = { 1.0, 1.0 };
  147.         double dopplerSigma = 0.5; // m/s
  148.         double dopplerWeight = 1.0;
  149.         double altitude = 300; // altitude of the circular orbit in km.
  150.  
  151.         double timeSpan = 3; // days to propagate
  152.         double minElev = 10 * Math.PI / 180.; // minimum elevation in rad.
  153.         boolean useRangeMeas = false;
  154.         boolean useAzelMeas = false;
  155.         boolean useDopplerMeas = true;
  156.         double step = 60; // seconds
  157.  
  158.         int simulatedMeasurements = (int) Math.ceil(timeSpan * 86400. / step);
  159.         GroundStation gs = getRadarStation();
  160.         Random rnd = new Random();
  161.  
  162.         // initial state vector.
  163.  
  164.         AbsoluteDate initialDate = mjd2000;
  165.         Orbit initialOrbitCart = new CartesianOrbit(new KeplerianOrbit(6378000 + altitude * 1000, 0.0,
  166.                 50 * Math.PI / 180, 0.0, 0.0, 0.0, PositionAngle.TRUE, FramesFactory.getEME2000(), initialDate, mu));
  167.         NumericalPropagatorBuilder propbuilder = getPropagator(initialOrbitCart);
  168.  
  169.         SpacecraftState currentState = null;
  170.         System.out.println("Initial orbit:");
  171.         printOrbit(initialOrbitCart);
  172.  
  173.         double[] normalizedParameters = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  174.         NumericalPropagator prop = propbuilder.buildPropagator(normalizedParameters);
  175.  
  176.         // Propagation
  177.         System.out.println("Propagating...");
  178.         List<ObservedMeasurement> meas = new ArrayList<>();
  179.         List<TimeStampedPVCoordinates> pv = new ArrayList<>();
  180.         for (int i = 0; i < simulatedMeasurements; i++) {
  181.             AbsoluteDate nextdate = new AbsoluteDate(initialDate, step * i);
  182.             currentState = prop.propagate(nextdate);
  183.             TimeStampedPVCoordinates pos = currentState.getPVCoordinates();
  184.             pv.add(pos);
  185.             // Check if visible from radar.
  186.             double elevation = gs.getBaseFrame().getElevation(pos.getPosition(), currentState.getFrame(),
  187.                     currentState.getDate());
  188.             if (elevation < minElev) {
  189.                 continue;
  190.             }
  191.             if (useAzelMeas) {
  192.                 // Get azel observation
  193.                 double azimuth = gs.getBaseFrame().getAzimuth(pos.getPosition(), currentState.getFrame(),
  194.                         currentState.getDate());
  195.                 double[] azel = { azimuth + rnd.nextGaussian() * azelSigma[0],
  196.                         elevation + rnd.nextGaussian() * azelSigma[1] };
  197.                 // Add noisy AzEl measurement.
  198.                 ObservedMeasurement azelMeas = new AngularAzEl(gs, currentState.getDate(), azel, azelSigma, azelWeight);
  199.                 meas.add(azelMeas);
  200.             }
  201.             if (useRangeMeas) {
  202.                 // Get Range observation.
  203.                 double range = gs.getBaseFrame().getRange(pos.getPosition(), currentState.getFrame(),
  204.                         currentState.getDate());
  205.                 range += rnd.nextGaussian() * rangeSigma;
  206.                 ObservedMeasurement rangeMeas = new Range(gs, currentState.getDate(), range, rangeSigma, rangeWeight);
  207.                 meas.add(rangeMeas);
  208.             }
  209.             if (useDopplerMeas) {
  210.                 // Get Doppler observation.
  211.                 double doppler = gs.getBaseFrame().getRangeRate(pos, currentState.getFrame(), currentState.getDate());
  212.                 doppler += rnd.nextGaussian() * dopplerSigma;
  213.                 ObservedMeasurement dopplerMeas = new Range(gs, currentState.getDate(), doppler, dopplerSigma,
  214.                         dopplerWeight);
  215.                 meas.add(dopplerMeas);
  216.             }
  217.         }
  218.         System.out.println("Number of measurements: " + meas.size());
  219.  
  220.         // Orbit determination.
  221.         System.out.println("Running the orbit determination...");
  222.         GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(Decomposition.QR);
  223.         NumericalPropagatorBuilder propbuilder2 = getPropagator(new CartesianOrbit(initialOrbitCart));
  224.         BatchLSEstimator estimator = new BatchLSEstimator(optimizer, propbuilder2);
  225.         estimator.setObserver(new LeastSquareObserver());
  226.  
  227.         estimator.setMaxEvaluations(20);
  228.         estimator.setMaxIterations(10);
  229.         estimator.setParametersConvergenceThreshold(0.01);
  230.         // Add measurements
  231.         for (ObservedMeasurement oneMeas : meas) {
  232.             estimator.addMeasurement(oneMeas);
  233.         }
  234.         NumericalPropagator[] neworb = estimator.estimate();
  235.         System.out.println("Determined Orbit");
  236.         printOrbit(neworb[0].getInitialState());
  237.         Optimum optimum = estimator.getOptimum();
  238.         System.out.println(
  239.                 "Optimum: " + optimum.getEvaluations() + " " + optimum.getIterations() + " " + optimum.getRMS());
  240.         System.out.println("Residuals: " + optimum.getResiduals());
  241.         System.out.println("Covariance matrix: " + optimum.getCovariances(0.0));
  242.         for (DelegatingDriver param : estimator.getPropagatorParametersDrivers(false).getDrivers()) {
  243.             System.out.println("Param " + param.getName() + " is " + param.getValue());
  244.         }
  245.  
  246.     }
  247.  
  248.     public static GroundStation getRadarStation() throws OrekitException {
  249.         OneAxisEllipsoid refbody2 = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
  250.                 Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, false));
  251.  
  252.         GeodeticPoint stationCoord = new GeodeticPoint(40.41 / 180. * Math.PI, -3.4 / 180. * Math.PI, 700.0);
  253.         TopocentricFrame frame = new TopocentricFrame(refbody2, stationCoord, "Radar station");
  254.  
  255.         GroundStation gs = new GroundStation(frame);
  256.         return gs;
  257.     }
  258.  
  259.     public static void printOrbit(SpacecraftState sc) throws OrekitException {
  260.         printOrbit(sc.getOrbit());
  261.     }
  262.  
  263.     public static void printOrbit(Orbit orbit) throws OrekitException {
  264.         System.out.println("***");
  265.         System.out.println(" Date: " + orbit.getDate().toString(TimeScalesFactory.getTT()));
  266.         TimeStampedPVCoordinates pos = orbit.getPVCoordinates(FramesFactory.getEME2000());
  267.         System.out.println(" Position: " + pos.getPosition());
  268.         System.out.println(" Velocity: " + pos.getVelocity());
  269.         System.out.println(" a: " + orbit.getA() + " / e: " + orbit.getE() + " / i: " + orbit.getI());
  270.     }
  271. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!