Commit d50e49f2 authored by Luc Maisonobe's avatar Luc Maisonobe

Merge remote-tracking branch 'origin/master' into develop

Conflicts:
	src/changes/changes.xml
parents 3bdc5f1a 9853acdf
......@@ -21,6 +21,12 @@
</properties>
<body>
<release version="10.3" date="TBD" description="TBD">
<action dev="luc" type="add" issue="7">
Added tutorial for constant thrust maneuvers.
</action>
<action dev="luc" type="add" issue="7">
Added tutorial for impulsive maneuvers.
</action>
</release>
<release version="10.2" date="2020-07-16"
description="Version 10.2 is the second release of the Orekit Tutorials.
......
......@@ -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,24 @@ 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);
latitude, ascending, mst, propagator);
orbit = improveFrozenEccentricity(tmp2, propagator, nbDays, nbOrbits);
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 +238,18 @@ 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 = new FittedEccentricity(orbit, nbDays, nbOrbits);
fittedEccentricity.fit(propagator);
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 +260,161 @@ 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 {
/** Initial orbit. */
private final CircularOrbit initial;
/** Cycle end date. */
private final AbsoluteDate tEnd;
/** Frozen eccentricity pulsation. */
private final double eta;
/** Sampled points. */
private final List<Observation> observed;
/** 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 initial orbit at start time
* @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);
this.initial = initial;
this.tEnd = initial.getDate().shiftedBy(nbDays * Constants.JULIAN_DAY);
this.observed = new ArrayList<>();
}
/** Perform fitting.
* @param propagator propagator to use
*/
public void fit(final Propagator propagator) {
propagator.resetInitialState(new SpacecraftState(initial));
final AbsoluteDate start = initial.getDate();
// sample orbit for one phasing cycle
propagator.setMasterMode(60, (state, isLast) -> observed.add(new Observation(state)));
propagator.propagate(start, tEnd);
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 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.
* @param date orbit date
* @param frame frame to use for defining orbit
......@@ -428,65 +601,30 @@ 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 orit with improved frozen eccentricity
*/
private CircularOrbit improveFrozenEccentricity(final CircularOrbit previous, final double duration,
final Propagator propagator) {
private CircularOrbit improveFrozenEccentricity(final CircularOrbit previous, final Propagator propagator,
final int nbDays, final int nbOrbits) {
propagator.resetInitialState(new SpacecraftState(previous));
final AbsoluteDate start = previous.getDate();
// fit eccentricity over one phasing cycle
final FittedEccentricity fittedE = new FittedEccentricity(previous, nbDays, nbOrbits);
fittedE.fit(propagator);
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());
}
// 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
// collapse mean evolution circular motion to its center point
return new CircularOrbit(previous.getA(),
previous.getCircularEx() + dex,
previous.getCircularEy() + dey,
previous.getCircularEx() - fittedE.dx0,
previous.getCircularEy() - fittedE.dy0,
previous.getI(), previous.getRightAscensionOfAscendingNode(),
previous.getAlphaV(), PositionAngle.TRUE,
previous.getFrame(), previous.getDate(),
previous.getMu());
}
/** Print ground track grid point.
......@@ -512,16 +650,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),
......
/* Copyright 2002-2020 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.tutorials.maneuvers;
import java.io.File;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.InertialProvider;
import org.orekit.attitudes.LofOffset;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.maneuvers.Maneuver;
import org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel;
import org.orekit.forces.maneuvers.propulsion.PropulsionModel;
import org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers;
import org.orekit.forces.maneuvers.trigger.ManeuverTriggers;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.LOFType;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateComponents;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
/**
* Orekit tutorial for an apogee maneuver.
*
* @author Luc Maisonobe
*/
public class ApogeeManeuver {
/** Private constructor for utility class. */
private ApogeeManeuver() {
// empty
}
/**
* Program entry point.
* @param args program arguments (unused here)
*/
public static void main(final String[] args) {
try {
// configure Orekit
final File home = new File(System.getProperty("user.home"));
final File orekitData = new File(home, "orekit-data");
if (!orekitData.exists()) {
System.err.format(Locale.US, "Failed to find %s folder%n",
orekitData.getAbsolutePath());
System.err.format(Locale.US, "You need to download %s from %s, unzip it in %s and rename it 'orekit-data' for this tutorial to work%n",
"orekit-data-master.zip", "https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip",
home.getAbsolutePath());
System.exit(1);
}
final DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
manager.addProvider(new DirectoryCrawler(orekitData));
// set up initial GTO orbit
final Frame eme2000 = FramesFactory.getEME2000();
final AbsoluteDate date = new AbsoluteDate(new DateComponents(2004, 01, 01),
new TimeComponents(23, 30, 00.000),
TimeScalesFactory.getUTC());
final Orbit orbit =
new KeplerianOrbit(24396159, 0.72831215, FastMath.toRadians(7),
FastMath.toRadians(180), FastMath.toRadians(261),
FastMath.toRadians(0), PositionAngle.TRUE,
eme2000, date,
Constants.EIGEN5C_EARTH_MU);
final SpacecraftState initialState = new SpacecraftState(orbit, 2500.0);
// prepare numerical propagator
final OrbitType orbitType = OrbitType.EQUINOCTIAL;
final double[][] tol = NumericalPropagator.tolerances(1.0, orbit, orbitType);
final AdaptiveStepsizeIntegrator integrator =
new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
integrator.setInitialStepSize(60);
final NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setInitialState(initialState);
propagator.setAttitudeProvider(new LofOffset(eme2000, LOFType.VNC));
// set up an attitude law dedicated to the maneuver
// where the +X axis (direction of acceleration of the thruster)
// points towards a specific direction
final Vector3D direction = new Vector3D(FastMath.toRadians(-7.4978),
FastMath.toRadians(351));
final AttitudeProvider attitudeOverride =
new InertialProvider(new Rotation(direction, Vector3D.PLUS_I),
eme2000);
// maneuver will start at a known date and stop after a known duration
final AbsoluteDate firingDate = new AbsoluteDate(new DateComponents(2004, 1, 2),
new TimeComponents(4, 15, 34.080),
TimeScalesFactory.getUTC());
final double duration = 3653.99;
final ManeuverTriggers triggers = new DateBasedManeuverTriggers(firingDate, duration);
// maneuver has constant thrust
final double thrust = 420;
final double isp = 318;
final PropulsionModel propulsionModel =
new BasicConstantThrustPropulsionModel(thrust, isp,
Vector3D.PLUS_I,
"apogee-engine");
// build maneuver and add it to the propagator as a new force model
propagator.addForceModel(new Maneuver(attitudeOverride, triggers, propulsionModel));
// progress monitoring
propagator.setMasterMode(120.0, (state, isLast) ->
System.out.format(Locale.US, "%s a = %12.3f m, e = %11.9f, m = %8.3f kg%n",
state.getDate(), state.getA(), state.getE(), state.getMass())
);
// propagate orbit, including maneuver
propagator.propagate(firingDate.shiftedBy(-900), firingDate.shiftedBy(duration + 900));
System.exit(0);
} catch (OrekitException e) {
System.err.println(e.getLocalizedMessage());
System.exit(1);
}
}
}
/* Copyright 2002-2020 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.tutorials.maneuvers;
import java.io.File;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.LofOffset;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.maneuvers.ImpulseManeuver;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.LOFType;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.propagation.events.EnablingPredicate;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.EventEnablingPredicateFilter;
import org.orekit.propagation.events.NodeDetector;
import org.orekit.propagation.events.handlers.StopOnIncreasing;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateComponents;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
/**