Commit 2e89c7d8 authored by Luc Maisonobe's avatar Luc Maisonobe

Revert "Fixed completely wrong frozen eccentricity handling..."

This reverts commit d175e4e6.
parent d5a7babd
Pipeline #696 passed with stages
in 1 minute and 58 seconds
......@@ -23,7 +23,8 @@ import java.io.PrintStream;
import java.net.URISyntaxException;
import java.net.URL;
import java.nio.charset.StandardCharsets;
import java.util.ArrayList;
import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.util.List;
import java.util.Locale;
......@@ -35,13 +36,8 @@ 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;
......@@ -185,14 +181,9 @@ public class Phasing {
// initial guess for orbit
CircularOrbit orbit = guessOrbit(date, frame, nbOrbits, nbDays,
latitude, ascending, mst);
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");
System.out.println("initial orbit: " + orbit);
System.out.println("please wait while orbit is adjusted...");
System.out.println();
// numerical model for improving orbit
final double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit, OrbitType.CIRCULAR);
......@@ -211,33 +202,25 @@ public class Phasing {
double deltaV = Double.POSITIVE_INFINITY;
int counter = 0;
while (deltaP > 1.0e-3 || deltaV > 1.0e-6) {
final DecimalFormat f = new DecimalFormat("0.000E00", new DecimalFormatSymbols(Locale.US));
while (deltaP > 3.0e-1 || deltaV > 3.0e-4) {
final CircularOrbit previous = orbit;
final CircularOrbit tmp1 = improveEarthPhasing(previous, nbOrbits, nbDays, propagator);
final CircularOrbit tmp2 = improveSunSynchronization(tmp1, nbOrbits * tmp1.getKeplerianPeriod(),
latitude, ascending, mst, 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());
orbit = improveFrozenEccentricity(tmp2, nbOrbits * tmp2.getKeplerianPeriod(), propagator);
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.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);
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");
final PVCoordinates delta = new PVCoordinates(previous.getPVCoordinates(),
orbit.getPVCoordinates());
......@@ -247,17 +230,8 @@ public class Phasing {
}
// final 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);
System.out.println();
System.out.println("final orbit (osculating): " + orbit);
// generate the ground track grid file
try (PrintStream output = new PrintStream(new File(input.getParent(), gridOutput), StandardCharsets.UTF_8.name())) {
......@@ -268,157 +242,6 @@ 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
......@@ -605,26 +428,64 @@ public class Phasing {
}
/** Fit eccentricity model.
/** Improve orbit to better match frozen eccentricity property.
* @param previous previous orbit
* @param duration sampling duration
* @param propagator propagator to use
* @param nbDays number of days of the phasing cycle
* @param nbOrbits number of orbits of the phasing cycle
* @return fitted eccentricity
* @return an improved Earth phased, Sun synchronous orbit with frozen eccentricity
*/
private FittedEccentricity fitEccentricity(final CircularOrbit previous, final Propagator propagator,
final int nbDays, final int nbOrbits) {
private CircularOrbit improveFrozenEccentricity(final CircularOrbit previous, final double duration,
final Propagator propagator) {
propagator.resetInitialState(new SpacecraftState(previous));
final AbsoluteDate start = previous.getDate();
// 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();
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());
}
return fittedE;
// 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());
}
......@@ -651,16 +512,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.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)));
out.println(fTime.format(crossing.getDate().durationFrom(start)) +
" " + fAngle.format(FastMath.toDegrees(gp.getLatitude())) +
" " + fAngle.format(FastMath.toDegrees(gp.getLongitude())) +
" " + ascending);
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