Skip to content

GitLab

Commit d175e4e6 authored by Luc Maisonobe's avatar Luc Maisonobe
1 Browse files
Options

Fixed completely wrong frozen eccentricity handling...

parent b0a55561
Pipeline #668 passed with stages
      in 1 minute and 40 seconds
      Hide whitespace changes
      Inline Side-by-side
      Showing with 209 additions and 70 deletions
      src/main/java/org/orekit/tutorials/bodies/Phasing.java
      View file @ d175e4e6
      ......@@ -23,8 +23,7 @@ import java.io.PrintStream;
      import java.net.URISyntaxException;
      import java.net.URL;
      import java.nio.charset.StandardCharsets;
      import java.text.DecimalFormat;
      import java.text.DecimalFormatSymbols;
      import java.util.ArrayList;
      import java.util.List;
      import java.util.Locale;
      ......@@ -36,8 +35,13 @@ import org.hipparchus.exception.LocalizedCoreFormats;
      import org.hipparchus.exception.MathRuntimeException;
      import org.hipparchus.geometry.euclidean.threed.Vector3D;
      import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
      import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresBuilder;
      import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer;
      import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem;
      import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
      import org.hipparchus.util.FastMath;
      import org.hipparchus.util.MathUtils;
      import org.hipparchus.util.SinCos;
      import org.orekit.bodies.BodyShape;
      import org.orekit.bodies.CelestialBodyFactory;
      import org.orekit.bodies.GeodeticPoint;
      ......@@ -181,9 +185,14 @@ public class Phasing {
      // initial guess for orbit
      CircularOrbit orbit = guessOrbit(date, frame, nbOrbits, nbDays,
      latitude, ascending, mst);
      System.out.println("initial orbit: " + orbit);
      System.out.println("please wait while orbit is adjusted...");
      System.out.println();
      System.out.format(Locale.US, "initial osculating orbit%n");
      System.out.format(Locale.US, " date = %s%n",
      orbit.getDate());
      System.out.format(Locale.US, " a = %14.6f m, ex = %17.10e, ey = %17.10e, i = %12.9f deg, Ω = %12.9f deg%n",
      orbit.getA(), orbit.getCircularEx(), orbit.getCircularEy(),
      FastMath.toDegrees(orbit.getI()),
      FastMath.toDegrees(orbit.getRightAscensionOfAscendingNode()));
      System.out.format(Locale.US, "please wait while orbit is adjusted...%n%n");
      // numerical model for improving orbit
      final double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit, OrbitType.CIRCULAR);
      ......@@ -202,25 +211,33 @@ public class Phasing {
      double deltaV = Double.POSITIVE_INFINITY;
      int counter = 0;
      final DecimalFormat f = new DecimalFormat("0.000E00", new DecimalFormatSymbols(Locale.US));
      while (deltaP > 3.0e-1 || deltaV > 3.0e-4) {
      while (deltaP > 1.0e-3 || deltaV > 1.0e-6) {
      final CircularOrbit previous = orbit;
      final CircularOrbit tmp1 = improveEarthPhasing(previous, nbOrbits, nbDays, propagator);
      final CircularOrbit tmp2 = improveSunSynchronization(tmp1, nbOrbits * tmp1.getKeplerianPeriod(),
      latitude, ascending, mst, propagator);
      orbit = improveFrozenEccentricity(tmp2, nbOrbits * tmp2.getKeplerianPeriod(), propagator);
      FittedEccentricity fittedEccentricity = fitEccentricity(tmp2, propagator, nbDays, nbOrbits);
      // collapse mean evolution circular motion to its center point
      orbit = new CircularOrbit(tmp2.getA(),
      tmp2.getCircularEx() - fittedEccentricity.dx0,
      tmp2.getCircularEy() - fittedEccentricity.dy0,
      tmp2.getI(), tmp2.getRightAscensionOfAscendingNode(),
      tmp2.getAlphaV(), PositionAngle.TRUE,
      tmp2.getFrame(), previous.getDate(),
      tmp2.getMu());
      final double da = orbit.getA() - previous.getA();
      final double dex = orbit.getCircularEx() - previous.getCircularEx();
      final double dey = orbit.getCircularEy() - previous.getCircularEy();
      final double di = FastMath.toDegrees(orbit.getI() - previous.getI());
      final double dr = FastMath.toDegrees(orbit.getRightAscensionOfAscendingNode() -
      previous.getRightAscensionOfAscendingNode());
      System.out.println(" iteration " + (++counter) + ": deltaA = " + f.format(da) +
      " m, deltaEx = " + f.format(dex) + ", deltaEy = " + f.format(dey) +
      ", deltaI = " + f.format(di) + " deg, deltaRAAN = " + f.format(dr) +
      " deg");
      System.out.format(Locale.US,
      " iteration %2d: deltaA = %12.6f m, Δex = %13.6e, Δey = %13.6e, Δi = %12.9f deg, ΔΩ = %12.9f deg%n",
      ++counter, da, dex, dey, di, dr);
      final PVCoordinates delta = new PVCoordinates(previous.getPVCoordinates(),
      orbit.getPVCoordinates());
      ......@@ -230,8 +247,17 @@ public class Phasing {
      }
      // final orbit
      System.out.println();
      System.out.println("final orbit (osculating): " + orbit);
      System.out.format(Locale.US, "%nfinal osculating orbit%n");
      System.out.format(Locale.US, " date = %s%n",
      orbit.getDate());
      System.out.format(Locale.US, " a = %14.6f m, ex = %17.10e, ey = %17.10e, i = %12.9f deg, Ω = %12.9f deg%n",
      orbit.getA(), orbit.getCircularEx(), orbit.getCircularEy(),
      FastMath.toDegrees(orbit.getI()),
      FastMath.toDegrees(orbit.getRightAscensionOfAscendingNode()));
      System.out.format(Locale.US, "%nfinal frozen eccentricity%n");
      final FittedEccentricity fittedEccentricity = fitEccentricity(orbit, propagator, nbDays, nbOrbits);
      System.out.format(Locale.US, " ex_f = %17.10e, ey_f = %17.10e%n",
      fittedEccentricity.cx, fittedEccentricity.cy);
      // generate the ground track grid file
      try (PrintStream output = new PrintStream(new File(input.getParent(), gridOutput), StandardCharsets.UTF_8.name())) {
      ......@@ -242,6 +268,157 @@ public class Phasing {
      }
      /** Fitted eccentricity model.
      * <ul>
      * <li>the mean model is harmonic at frozen eccentricity pulsation</li>
      * <li>the osculating model adds harmonic components with periods T and T/3, where T is orbital period</li>
      * </ul>
      */
      private class FittedEccentricity {
      /** Observation points for the fitting. */
      private final List<Observation> observed;
      /** Frozen eccentricity pulsation. */
      private final double eta;
      /** Model reference date. */
      private AbsoluteDate tR;
      /** Center X component of the mean eccentricity. */
      private double cx;
      /** Center Y component of the mean eccentricity. */
      private double cy;
      /** Initial X offset of mean eccentricity. */
      private double dx0;
      /** Initial Y offset of mean eccentricity. */
      private double dy0;
      /** Simple constructor.
      * @param tR model reference date
      * @param inclination orbit inclination
      * @param nbDays number of days of the phasing cycle
      * @param nbOrbits number of orbits of the phasing cycle
      */
      public FittedEccentricity(final AbsoluteDate tR, final double inclination,
      final int nbDays, final int nbOrbits) {
      // extract gravity field data
      final double referenceRadius = gravityField.getAe();
      final double mu = gravityField.getMu();
      final double[][] unnormalization = GravityFieldFactory.getUnnormalizationFactors(3, 0);
      final double j2 = -unnormalization[2][0] * gravityField.onDate(tR).getNormalizedCnm(2, 0);
      final double period = nbDays * Constants.JULIAN_DAY / nbOrbits;
      final double meanMotion = 2 * FastMath.PI / period;
      final double sinI = FastMath.sin(inclination);
      final double a = FastMath.cbrt(mu / (meanMotion * meanMotion));
      final double rOa = referenceRadius / a;
      this.eta = 3 * meanMotion * j2 * rOa * rOa * (1.25 * sinI * sinI - 1.0);
      this.tR = tR;
      this.observed = new ArrayList<>();
      }
      /** Add a sampled orbit data.
      * @param orbit orbit at sampling date
      */
      public void addSamplingPoint(final Orbit orbit) {
      final CircularOrbit circular = (CircularOrbit) OrbitType.CIRCULAR.convertType(orbit);
      observed.add(new Observation(circular, tR));
      }
      /** Perform fitting.
      */
      public void fit() {
      final LeastSquaresProblem lsp = new LeastSquaresBuilder().
      maxEvaluations(1000).
      maxIterations(1000).
      start(new double[7]).
      target(new double[2 * observed.size()]).
      model(params -> residuals(params),
      params -> jacobian(params)).
      build();
      final LeastSquaresOptimizer.Optimum optimum = new LevenbergMarquardtOptimizer().optimize(lsp);
      // store coefficients (for mean model only)
      cx = optimum.getPoint().getEntry(0);
      cy = optimum.getPoint().getEntry(1);
      dx0 = optimum.getPoint().getEntry(2);
      dy0 = optimum.getPoint().getEntry(3);
      }
      /** Value of the error model.
      * @param params fitting parameters
      * @return model value
      */
      private double[] residuals(final double[] params) {
      final double[] val = new double[2 * observed.size()];
      int i = 0;
      for (final Observation obs : observed) {
      final SinCos sc = FastMath.sinCos(eta * obs.dt);
      final SinCos sc1 = FastMath.sinCos(obs.alphaM);
      final SinCos sc3 = FastMath.sinCos(3 * obs.alphaM);
      val[i++] = params[0] + params[2] * sc.cos() + params[3] * sc.sin() +
      params[4] * sc1.cos() + params[6] * sc3.cos() -
      obs.ex;
      val[i++] = params[1] - params[2] * sc.sin() + params[3] * sc.cos() +
      params[5] * sc1.sin() + params[6] * sc3.sin() -
      obs.ey;
      }
      return val;
      }
      /** Jacobian of the error model.
      * @param params fitting parameters
      * @return model Jacobian
      */
      private double[][] jacobian(final double[] params) {
      final double[][] jac = new double[2 * observed.size()][];
      int i = 0;
      for (final Observation obs : observed) {
      final SinCos sc = FastMath.sinCos(eta * obs.dt);
      final SinCos sc1 = FastMath.sinCos(obs.alphaM);
      final SinCos sc3 = FastMath.sinCos(3 * obs.alphaM);
      jac[i++] = new double[] { 1, 0, sc.cos(), sc.sin(), sc1.cos(), 0, sc3.cos() };
      jac[i++] = new double[] { 0, 1, -sc.sin(), sc.cos(), 0, sc1.sin(), sc3.sin() };
      }
      return jac;
      }
      private class Observation {
      /** Date offset since reference. */
      private double dt;
      /** Mean latitude argument. */
      private double alphaM;
      /** X component of eccentricity. */
      private double ex;
      /** Y component of eccentricity. */
      private double ey;
      /** Simple constructor.
      * @param orbit orbit at observation time
      * @param reference reference time
      */
      Observation(final CircularOrbit orbit, final AbsoluteDate reference) {
      this.dt = orbit.getDate().durationFrom(reference);
      this.alphaM = orbit.getAlphaM();
      this.ex = orbit.getCircularEx();
      this.ey = orbit.getCircularEy();
      }
      }
      }
      /** Guess an initial orbit from theoretical model.
      * @param date orbit date
      * @param frame frame to use for defining orbit
      ......@@ -428,64 +605,26 @@ public class Phasing {
      }
      /** Improve orbit to better match frozen eccentricity property.
      /** Fit eccentricity model.
      * @param previous previous orbit
      * @param duration sampling duration
      * @param propagator propagator to use
      * @return an improved Earth phased, Sun synchronous orbit with frozen eccentricity
      * @param nbDays number of days of the phasing cycle
      * @param nbOrbits number of orbits of the phasing cycle
      * @return fitted eccentricity
      */
      private CircularOrbit improveFrozenEccentricity(final CircularOrbit previous, final double duration,
      final Propagator propagator) {
      private FittedEccentricity fitEccentricity(final CircularOrbit previous, final Propagator propagator,
      final int nbDays, final int nbOrbits) {
      propagator.resetInitialState(new SpacecraftState(previous));
      final AbsoluteDate start = previous.getDate();
      final NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics harmonics =
      gravityField.onDate(previous.getDate());
      final double[][] unnormalization = GravityFieldFactory.getUnnormalizationFactors(2, 0);
      final double a = previous.getA();
      final double sinI = FastMath.sin(previous.getI());
      final double aeOa = gravityField.getAe() / a;
      final double mu = gravityField.getMu();
      final double n = FastMath.sqrt(mu / a) / a;
      final double j2 = -unnormalization[2][0] * harmonics.getNormalizedCnm(2, 0);
      final double frozenPulsation = 3 * n * j2 * aeOa * aeOa * (1 - 1.25 * sinI * sinI);
      // fit the eccentricity to an harmonic model with short and medium periods
      // we will only use the medium periods part for the correction
      final SecularAndHarmonic exModel = new SecularAndHarmonic(0, frozenPulsation, n, 2 * n);
      final SecularAndHarmonic eyModel = new SecularAndHarmonic(0, frozenPulsation, n, 2 * n);
      exModel.resetFitting(start, new double[] {
      previous.getCircularEx(), -1.0e-10, 1.0e-5,
      1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5
      });
      eyModel.resetFitting(start, new double[] {
      previous.getCircularEy(), -1.0e-10, 1.0e-5,
      1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5
      });
      final double step = previous.getKeplerianPeriod() / 5;
      for (double dt = 0; dt < duration; dt += step) {
      final SpacecraftState state = propagator.propagate(start.shiftedBy(dt));
      final CircularOrbit orbit = (CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit());
      exModel.addPoint(state.getDate(), orbit.getCircularEx());
      eyModel.addPoint(state.getDate(), orbit.getCircularEy());
      }
      // sample orbit for one phasing cycle
      final FittedEccentricity fittedE = new FittedEccentricity(start, previous.getI(), nbDays, nbOrbits);
      propagator.setMasterMode(60, (state, isLast) -> fittedE.addSamplingPoint(state.getOrbit()));
      propagator.propagate(start, start.shiftedBy(Constants.JULIAN_DAY * nbDays));
      fittedE.fit();
      // adjust eccentricity
      exModel.fit();
      final double dex = -exModel.getFittedParameters()[1];
      eyModel.fit();
      final double dey = -eyModel.getFittedParameters()[1];
      // put the eccentricity at center of frozen center
      return new CircularOrbit(previous.getA(),
      previous.getCircularEx() + dex,
      previous.getCircularEy() + dey,
      previous.getI(), previous.getRightAscensionOfAscendingNode(),
      previous.getAlphaV(), PositionAngle.TRUE,
      previous.getFrame(), previous.getDate(),
      previous.getMu());
      return fittedE;
      }
      ......@@ -512,16 +651,16 @@ public class Phasing {
      stepSize, propagator);
      // find all other latitude crossings from regular schedule
      final DecimalFormat fTime = new DecimalFormat("0000000.000", new DecimalFormatSymbols(Locale.US));
      final DecimalFormat fAngle = new DecimalFormat("000.00000", new DecimalFormatSymbols(Locale.US));
      for (int i = 0; i < nbOrbits; ++i) {
      final CircularOrbit c = (CircularOrbit) OrbitType.CIRCULAR.convertType(crossing.getOrbit());
      final GeodeticPoint gp = earth.transform(crossing.getPVCoordinates().getPosition(),
      crossing.getFrame(), crossing.getDate());
      out.println(fTime.format(crossing.getDate().durationFrom(start)) +
      " " + fAngle.format(FastMath.toDegrees(gp.getLatitude())) +
      " " + fAngle.format(FastMath.toDegrees(gp.getLongitude())) +
      " " + ascending);
      out.format(Locale.US, "%11.3f %9.5f %9.5f %s %11.5f%n",
      crossing.getDate().durationFrom(start),
      FastMath.toDegrees(gp.getLatitude()), FastMath.toDegrees(gp.getLongitude()),
      ascending,
      FastMath.toDegrees(MathUtils.normalizeAngle(c.getAlphaV(), 0)));
      final AbsoluteDate previousDate = crossing.getDate();
      crossing = findLatitudeCrossing(latitude, previousDate.shiftedBy(period),
      ......
      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