Commit b0c64152 authored by Luc Maisonobe's avatar Luc Maisonobe

Remove eccentricity short periods analytically.

parent e719b9e0
Pipeline #704 passed with stages
in 2 minutes and 23 seconds
......@@ -26,6 +26,7 @@ import java.nio.charset.StandardCharsets;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import java.util.stream.Collectors;
import org.hipparchus.analysis.UnivariateFunction;
import org.hipparchus.analysis.solvers.BaseUnivariateSolver;
......@@ -248,10 +249,10 @@ public class Phasing {
FastMath.toDegrees(orbit.getI()),
FastMath.toDegrees(orbit.getRightAscensionOfAscendingNode()));
System.out.format(Locale.US, "%nfinal frozen eccentricity%n");
final FittedEccentricity fittedEccentricity = new FittedEccentricity(orbit, nbDays, nbOrbits);
fittedEccentricity.fit(propagator);
final MeanEccentricityFitter eccentricityFitter = new MeanEccentricityFitter(orbit, nbDays, nbOrbits);
eccentricityFitter.fit(propagator);
System.out.format(Locale.US, " ex_f = %17.10e, ey_f = %17.10e%n",
fittedEccentricity.cx, fittedEccentricity.cy);
eccentricityFitter.cx, eccentricityFitter.cy);
// generate the ground track grid file
try (PrintStream output = new PrintStream(new File(input.getParent(), gridOutput), StandardCharsets.UTF_8.name())) {
......@@ -262,13 +263,13 @@ public class Phasing {
}
/** Fitted eccentricity model.
/** Fitter for mean 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>
* <li>the short period terms with periods T and T/3 are removed analytically during fitting</li>
* </ul>
*/
private class FittedEccentricity {
private class MeanEccentricityFitter {
/** Initial orbit. */
private final CircularOrbit initial;
......@@ -276,11 +277,20 @@ public class Phasing {
/** Cycle end date. */
private final AbsoluteDate tEnd;
/** Observation points for the fitting. */
private final List<CircularOrbit> observed;
/** Frozen eccentricity pulsation. */
private final double eta;
private double eta;
/** X component of period T harmonic. */
private double x1;
/** Sampled points. */
private final List<Observation> observed;
/** Y component of period T harmonic. */
private double y1;
/** Common X/Y component of period T/3 harmonic. */
private double xy3;
/** Center X component of the mean eccentricity. */
private double cx;
......@@ -299,25 +309,11 @@ public class Phasing {
* @param nbDays number of days of the phasing cycle
* @param nbOrbits number of orbits of the phasing cycle
*/
public FittedEccentricity(final CircularOrbit initial,
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(initial.getDate()).getNormalizedCnm(2, 0);
final double period = nbDays * Constants.JULIAN_DAY / nbOrbits;
final double meanMotion = 2 * FastMath.PI / period;
final double sinI = FastMath.sin(initial.getI());
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);
public MeanEccentricityFitter(final CircularOrbit initial,
final int nbDays, final int nbOrbits) {
this.initial = initial;
this.tEnd = initial.getDate().shiftedBy(nbDays * Constants.JULIAN_DAY);
this.observed = new ArrayList<>();
}
/** Perform fitting.
......@@ -329,13 +325,41 @@ public class Phasing {
final AbsoluteDate start = initial.getDate();
// sample orbit for one phasing cycle
propagator.setMasterMode(60, (state, isLast) -> observed.add(new Observation(state)));
propagator.setMasterMode(60,
(state, isLast) ->
observed.add((CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit())));
propagator.propagate(start, tEnd);
// compute mean semi-major axis and mean inclination
final double meanA = observed.stream().collect(Collectors.averagingDouble(c -> c.getA()));
final double meanI = observed.stream().collect(Collectors.averagingDouble(c -> c.getI()));
// extract gravity field data
final double referenceRadius = gravityField.getAe();
final double mu = gravityField.getMu();
final double[][] unnormalization = GravityFieldFactory.getUnnormalizationFactors(2, 0);
final double j2 = -unnormalization[2][0] * gravityField.onDate(initial.getDate()).getNormalizedCnm(2, 0);
// compute coefficient of harmonic terms
final double meanMotion = FastMath.sqrt(mu / meanA) / meanA;
final double sinI = FastMath.sin(meanI);
final double sinI2 = sinI * sinI;
final double rOa = referenceRadius / meanA;
// mean frozen eccentricity pulsation (long period)
this.eta = 3 * meanMotion * j2 * rOa * rOa * (1.25 * sinI2 - 1.0);
// short periods
final double kappa1 = 1.5 * j2 * rOa * rOa;
final double kappa2 = kappa1 / (1.0 - kappa1 * (3.0 - 4.0 * sinI2));
x1 = kappa2 * (1.0 - 1.25 * sinI2);
y1 = kappa2 * (1.0 - 1.75 * sinI2);
xy3 = kappa2 * 7.0 / 12.0 * sinI2;
final LeastSquaresProblem lsp = new LeastSquaresBuilder().
maxEvaluations(1000).
maxIterations(1000).
start(new double[7]).
start(new double[4]).
target(new double[2 * observed.size()]).
model(params -> residuals(params),
params -> jacobian(params)).
......@@ -357,16 +381,16 @@ public class Phasing {
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);
for (final CircularOrbit c : observed) {
final double dt = c.getDate().durationFrom(initial.getDate());
final double alphaM = c.getAlphaM();
final SinCos sc = FastMath.sinCos(eta * dt);
final SinCos sc1 = FastMath.sinCos(alphaM);
final SinCos sc3 = FastMath.sinCos(3 * alphaM);
val[i++] = params[0] + params[2] * sc.cos() + params[3] * sc.sin() +
params[4] * sc1.cos() + params[6] * sc3.cos() -
obs.ex;
x1 * sc1.cos() + xy3 * sc3.cos() - c.getCircularEx();
val[i++] = params[1] - params[2] * sc.sin() + params[3] * sc.cos() +
params[5] * sc1.sin() + params[6] * sc3.sin() -
obs.ey;
y1 * sc1.sin() + xy3 * sc3.sin() - c.getCircularEy();
}
return val;
}
......@@ -378,43 +402,15 @@ public class Phasing {
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() };
for (final CircularOrbit c : observed) {
final double dt = c.getDate().durationFrom(initial.getDate());
final SinCos sc = FastMath.sinCos(eta * dt);
jac[i++] = new double[] { 1, 0, sc.cos(), sc.sin() };
jac[i++] = new double[] { 0, 1, -sc.sin(), sc.cos() };
}
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 state spacecraft state at observation time
*/
Observation(final SpacecraftState state) {
final CircularOrbit orbit = (CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit());
this.dt = orbit.getDate().durationFrom(initial.getDate());
this.alphaM = orbit.getAlphaM();
this.ex = orbit.getCircularEx();
this.ey = orbit.getCircularEy();
}
}
}
/** Guess an initial orbit from theoretical model.
......@@ -613,14 +609,14 @@ public class Phasing {
private CircularOrbit improveFrozenEccentricity(final CircularOrbit previous, final Propagator propagator,
final int nbDays, final int nbOrbits) {
// fit eccentricity over one phasing cycle
final FittedEccentricity fittedE = new FittedEccentricity(previous, nbDays, nbOrbits);
fittedE.fit(propagator);
// fit mean eccentricity over one phasing cycle
final MeanEccentricityFitter eccentricityFitter = new MeanEccentricityFitter(previous, nbDays, nbOrbits);
eccentricityFitter.fit(propagator);
// collapse mean evolution circular motion to its center point
return new CircularOrbit(previous.getA(),
previous.getCircularEx() - fittedE.dx0,
previous.getCircularEy() - fittedE.dy0,
previous.getCircularEx() - eccentricityFitter.dx0,
previous.getCircularEy() - eccentricityFitter.dy0,
previous.getI(), previous.getRightAscensionOfAscendingNode(),
previous.getAlphaV(), PositionAngle.TRUE,
previous.getFrame(), previous.getDate(),
......@@ -935,7 +931,7 @@ public class Phasing {
this.gravity = gravity;
}
/**
/**)
* Get the list of grid data.
* @return the list of grid data
*/
......
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