Commit 424bbc36 authored by Luc Maisonobe's avatar Luc Maisonobe

Added m-dailies to eccentricity fitting.

parent 67efd526
Pipeline #722 passed with stages
in 2 minutes and 3 seconds
......@@ -163,6 +163,7 @@ public class Phasing {
final Frame frame = inputData.getOrbit().getInertialFrame();
final int nbOrbits = inputData.getPhasingOrbitsNumber();
final int nbDays = inputData.getPhasingDaysNumer();
final int maxMDaily = inputData.getMaxMDaily();
final double latitude = FastMath.toRadians(inputData.getReferenceLatitude());
final boolean ascending = inputData.isReferenceAscending();
final double mst = TimeComponents.parseTime(inputData.getMeanSolarTime()).getSecondsInUTCDay() / 3600;
......@@ -197,9 +198,9 @@ public class Phasing {
// numerical model for improving orbit
final OrbitType orbitType = OrbitType.CIRCULAR;
final double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit, orbitType);
final double[][] tolerances = NumericalPropagator.tolerances(0.001, orbit, orbitType);
final DormandPrince853Integrator integrator =
new DormandPrince853Integrator(1.0e-4 * orbit.getKeplerianPeriod(),
new DormandPrince853Integrator(1.0e-5 * orbit.getKeplerianPeriod(),
1.0e-1 * orbit.getKeplerianPeriod(),
tolerances[0], tolerances[1]);
integrator.setInitialStepSize(1.0e-2 * orbit.getKeplerianPeriod());
......@@ -221,7 +222,7 @@ public class Phasing {
final CircularOrbit tmp1 = improveEarthPhasing(previous, nbOrbits, nbDays, propagator);
final CircularOrbit tmp2 = improveSunSynchronization(tmp1, nbOrbits * tmp1.getKeplerianPeriod(),
latitude, ascending, mst, propagator);
orbit = improveFrozenEccentricity(tmp2, propagator, nbDays, nbOrbits);
orbit = improveFrozenEccentricity(tmp2, propagator, maxMDaily, nbDays, nbOrbits);
final double da = orbit.getA() - previous.getA();
final double dex = orbit.getCircularEx() - previous.getCircularEx();
......@@ -229,15 +230,15 @@ public class Phasing {
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: Δa = %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());
deltaP = delta.getPosition().getNorm();
deltaV = delta.getVelocity().getNorm();
System.out.format(Locale.US,
" iteration %2d: Δa = %12.6f m, Δex = %13.6e, Δey = %13.6e, Δi = %12.9f deg, ΔΩ = %12.9f deg, ΔP = %11.6f m, ΔV = %11.9f m/s%n",
++counter, da, dex, dey, di, dr, deltaP, deltaV);
}
// final orbit
......@@ -249,7 +250,7 @@ public class Phasing {
FastMath.toDegrees(orbit.getI()),
FastMath.toDegrees(orbit.getRightAscensionOfAscendingNode()));
System.out.format(Locale.US, "%nfinal frozen eccentricity%n");
final MeanEccentricityFitter eccentricityFitter = new MeanEccentricityFitter(orbit, nbDays, nbOrbits);
MeanEccentricityFitter eccentricityFitter = new MeanEccentricityFitter(orbit, maxMDaily, nbDays, nbOrbits);
eccentricityFitter.fit(propagator);
System.out.format(Locale.US, " ex_f = %17.10e, ey_f = %17.10e%n",
eccentricityFitter.cx, eccentricityFitter.cy);
......@@ -277,6 +278,9 @@ public class Phasing {
/** Cycle end date. */
private final AbsoluteDate tEnd;
/** Maximum index of m-daily terms. */
private final int maxMDaily;
/** Observation points for the fitting. */
private final List<CircularOrbit> observed;
......@@ -306,14 +310,16 @@ public class Phasing {
/** Simple constructor.
* @param initial orbit at start time
* @param maxMDaily maximum index of m-daily terms
* @param nbDays number of days of the phasing cycle
* @param nbOrbits number of orbits of the phasing cycle
*/
public MeanEccentricityFitter(final CircularOrbit initial,
public MeanEccentricityFitter(final CircularOrbit initial, final int maxMDaily,
final int nbDays, final int nbOrbits) {
this.initial = initial;
this.tEnd = initial.getDate().shiftedBy(nbDays * Constants.JULIAN_DAY);
this.observed = new ArrayList<>();
this.initial = initial;
this.tEnd = initial.getDate().shiftedBy(nbDays * Constants.JULIAN_DAY);
this.maxMDaily = maxMDaily;
this.observed = new ArrayList<>();
}
/** Perform fitting.
......@@ -359,7 +365,7 @@ public class Phasing {
final LeastSquaresProblem lsp = new LeastSquaresBuilder().
maxEvaluations(1000).
maxIterations(1000).
start(new double[4]).
start(new double[4 + 4 * maxMDaily]).
target(new double[2 * observed.size()]).
model(params -> residuals(params),
params -> jacobian(params)).
......@@ -367,10 +373,10 @@ public class Phasing {
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);
cx = optimum.getPoint().getEntry(0);
cy = optimum.getPoint().getEntry(1);
dx0 = optimum.getPoint().getEntry(2);
dy0 = optimum.getPoint().getEntry(3);
}
......@@ -387,10 +393,22 @@ public class Phasing {
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() +
x1 * sc1.cos() + xy3 * sc3.cos() - c.getCircularEx();
val[i++] = params[1] - params[2] * sc.sin() + params[3] * sc.cos() +
y1 * sc1.sin() + xy3 * sc3.sin() - c.getCircularEy();
final SinCos[] scM = new SinCos[maxMDaily];
for (int k = 0; k < scM.length; ++k) {
scM[k] = FastMath.sinCos(2 * (k + 1) * FastMath.PI * dt / Constants.JULIAN_DAY);
};
final double exM = params[0] + params[2] * sc.cos() + params[3] * sc.sin();
final double eyM = params[1] - params[2] * sc.sin() + params[3] * sc.cos();
final double exJ2 = x1 * sc1.cos() + xy3 * sc3.cos();
final double eyJ2 = y1 * sc1.sin() + xy3 * sc3.sin();
double exS = 0;
double eyS = 0;
for (int k = 0; k < scM.length; ++k) {
exS += params[4 * k + 4] * scM[k].cos() + params[4 * k + 5] * scM[k].sin();
eyS += params[4 * k + 6] * scM[k].cos() + params[4 * k + 7] * scM[k].sin();
}
val[i++] = exM + exJ2 + exS - c.getCircularEx();
val[i++] = eyM + eyJ2 + eyS - c.getCircularEy();
}
return val;
}
......@@ -405,8 +423,28 @@ public class Phasing {
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() };
final SinCos[] scM = new SinCos[maxMDaily];
for (int k = 0; k < scM.length; ++k) {
scM[k] = FastMath.sinCos(2 * (k + 1) * FastMath.PI * dt / Constants.JULIAN_DAY);
};
final double[] jacX = new double[4 + 4 * scM.length];
final double[] jacY = new double[4 + 4 * scM.length];
jacX[0] = 1;
jacX[1] = 0;
jacX[2] = sc.cos();
jacX[3] = sc.sin();
jacY[0] = 0;
jacY[1] = 1;
jacY[2] = -sc.sin();
jacY[3] = sc.cos();
for (int k = 0; k < scM.length; ++k) {
jacX[4 + 4 * k] = scM[k].cos();
jacX[5 + 4 * k] = scM[k].sin();
jacY[6 + 4 * k] = scM[k].cos();
jacY[7 + 4 * k] = scM[k].sin();
}
jac[i++] = jacX;
jac[i++] = jacY;
}
return jac;
}
......@@ -602,15 +640,16 @@ public class Phasing {
/** Fit eccentricity model.
* @param previous previous orbit
* @param propagator propagator to use
* @param maxMDaily maximum index of m-daily terms
* @param nbDays number of days of the phasing cycle
* @param nbOrbits number of orbits of the phasing cycle
* @return orit with improved frozen eccentricity
*/
private CircularOrbit improveFrozenEccentricity(final CircularOrbit previous, final Propagator propagator,
final int nbDays, final int nbOrbits) {
final int maxMDaily, final int nbDays, final int nbOrbits) {
// fit mean eccentricity over one phasing cycle
final MeanEccentricityFitter eccentricityFitter = new MeanEccentricityFitter(previous, nbDays, nbOrbits);
final MeanEccentricityFitter eccentricityFitter = new MeanEccentricityFitter(previous, maxMDaily, nbDays, nbOrbits);
eccentricityFitter.fit(propagator);
// collapse mean evolution circular motion to its center point
......@@ -813,6 +852,9 @@ public class Phasing {
/** Gravity data. */
private TutorialGravity gravity;
/** Maximum order of m-daily terms for eccentricity fitting. */
private int maxMDaily;
/** List of grid data. */
private List<TutorialGrid> grids;
......@@ -931,7 +973,23 @@ public class Phasing {
this.gravity = gravity;
}
/**)
/**
* Set the maximum order of m-daily terms for eccentricity fitting.
* @param maximum order of m-daily terms for eccentricity fitting
*/
public void setMaxMDaily(final int maxMDaily) {
this.maxMDaily = maxMDaily;
}
/**
* Get the maximum order of m-daily terms for eccentricity fitting.
* @return maximum order of m-daily terms for eccentricity fitting
*/
public int getMaxMDaily() {
return maxMDaily;
}
/**
* Get the list of grid data.
* @return the list of grid data
*/
......
......@@ -20,6 +20,9 @@ gravity:
degree: 10
order: 10
# Eccentricity fitting parameter
maxMDaily: 6
# ground track grid (angles in degrees)
grids:
- latitude: -60.0
......
......@@ -20,6 +20,9 @@ gravity:
degree: 10
order: 10
# Eccentricity fitting parameter
maxMDaily: 6
# ground track grid (angles in degrees)
grids:
- latitude: -60.0
......
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