J2DifferentialEffect.java
/* Copyright 2002-2023 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.propagation.analytical;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.SinCos;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.orbits.EquinoctialOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.time.AbsoluteDate;
/** Analytical model for J2 effect.
* <p>
* This class computes the differential effect of J2 due to an initial orbit
* offset. A typical case is when an inclination maneuver changes an orbit
* inclination at time t₀. As ascending node drift rate depends on
* inclination, the change induces a time-dependent change in ascending node
* for later dates.
* </p>
* @see org.orekit.forces.maneuvers.SmallManeuverAnalyticalModel
* @author Luc Maisonobe
*/
public class J2DifferentialEffect
implements AdapterPropagator.DifferentialEffect {
/** Reference date. */
private final AbsoluteDate referenceDate;
/** Differential drift on perigee argument. */
private final double dPaDot;
/** Differential drift on ascending node. */
private final double dRaanDot;
/** Indicator for applying effect before reference date. */
private final boolean applyBefore;
/** Simple constructor.
* <p>
* The {@code applyBefore} parameter is mainly used when the differential
* effect is associated with a maneuver. In this case, the parameter must be
* set to {@code false}.
* </p>
* @param original original state at reference date
* @param directEffect direct effect changing the orbit
* @param applyBefore if true, effect is applied both before and after
* reference date, if false it is only applied after reference date
* @param gravityField gravity field to use
*/
public J2DifferentialEffect(final SpacecraftState original,
final AdapterPropagator.DifferentialEffect directEffect,
final boolean applyBefore,
final UnnormalizedSphericalHarmonicsProvider gravityField) {
this(original, directEffect, applyBefore,
gravityField.getAe(), gravityField.getMu(),
-gravityField.onDate(original.getDate()).getUnnormalizedCnm(2, 0));
}
/** Simple constructor.
* <p>
* The {@code applyBefore} parameter is mainly used when the differential
* effect is associated with a maneuver. In this case, the parameter must be
* set to {@code false}.
* </p>
* @param orbit0 original orbit at reference date
* @param orbit1 shifted orbit at reference date
* @param applyBefore if true, effect is applied both before and after
* reference date, if false it is only applied after reference date
* @param gravityField gravity field to use
*/
public J2DifferentialEffect(final Orbit orbit0, final Orbit orbit1, final boolean applyBefore,
final UnnormalizedSphericalHarmonicsProvider gravityField) {
this(orbit0, orbit1, applyBefore,
gravityField.getAe(), gravityField.getMu(),
-gravityField.onDate(orbit0.getDate()).getUnnormalizedCnm(2, 0));
}
/** Simple constructor.
* <p>
* The {@code applyBefore} parameter is mainly used when the differential
* effect is associated with a maneuver. In this case, the parameter must be
* set to {@code false}.
* </p>
* @param original original state at reference date
* @param directEffect direct effect changing the orbit
* @param applyBefore if true, effect is applied both before and after
* reference date, if false it is only applied after reference date
* @param referenceRadius reference radius of the Earth for the potential model (m)
* @param mu central attraction coefficient (m³/s²)
* @param j2 un-normalized zonal coefficient (about +1.08e-3 for Earth)
*/
public J2DifferentialEffect(final SpacecraftState original,
final AdapterPropagator.DifferentialEffect directEffect,
final boolean applyBefore,
final double referenceRadius, final double mu, final double j2) {
this(original.getOrbit(),
directEffect.apply(original.shiftedBy(0.001)).getOrbit().shiftedBy(-0.001),
applyBefore, referenceRadius, mu, j2);
}
/** Simple constructor.
* <p>
* The {@code applyBefore} parameter is mainly used when the differential
* effect is associated with a maneuver. In this case, the parameter must be
* set to {@code false}.
* </p>
* @param orbit0 original orbit at reference date
* @param orbit1 shifted orbit at reference date
* @param applyBefore if true, effect is applied both before and after
* reference date, if false it is only applied after reference date
* @param referenceRadius reference radius of the Earth for the potential model (m)
* @param mu central attraction coefficient (m³/s²)
* @param j2 un-normalized zonal coefficient (about +1.08e-3 for Earth)
*/
public J2DifferentialEffect(final Orbit orbit0, final Orbit orbit1, final boolean applyBefore,
final double referenceRadius, final double mu, final double j2) {
this.referenceDate = orbit0.getDate();
this.applyBefore = applyBefore;
// extract useful parameters
final double a0 = orbit0.getA();
final double e0 = orbit0.getE();
final double i0 = orbit0.getI();
final double a1 = orbit1.getA();
final double e1 = orbit1.getE();
final double i1 = orbit1.getI();
// compute reference drifts
final double oMe2 = 1 - e0 * e0;
final double ratio = referenceRadius / (a0 * oMe2);
final SinCos scI = FastMath.sinCos(i0);
final double n = FastMath.sqrt(mu / a0) / a0;
final double c = ratio * ratio * n * j2;
final double refPaDot = 0.75 * c * (4 - 5 * scI.sin() * scI.sin());
final double refRaanDot = -1.5 * c * scI.cos();
// differential model on perigee argument drift
final double dPaDotDa = -3.5 * refPaDot / a0;
final double dPaDotDe = 4 * refPaDot * e0 / oMe2;
final double dPaDotDi = -7.5 * c * scI.sin() * scI.cos();
dPaDot = dPaDotDa * (a1 - a0) + dPaDotDe * (e1 - e0) + dPaDotDi * (i1 - i0);
// differential model on ascending node drift
final double dRaanDotDa = -3.5 * refRaanDot / a0;
final double dRaanDotDe = 4 * refRaanDot * e0 / oMe2;
final double dRaanDotDi = -refRaanDot * FastMath.tan(i0);
dRaanDot = dRaanDotDa * (a1 - a0) + dRaanDotDe * (e1 - e0) + dRaanDotDi * (i1 - i0);
}
/** Compute the effect of the maneuver on an orbit.
* @param orbit1 original orbit at t₁, without maneuver
* @return orbit at t₁, taking the maneuver
* into account if t₁ > t₀
* @see #apply(SpacecraftState)
*/
public Orbit apply(final Orbit orbit1) {
if (orbit1.getDate().compareTo(referenceDate) <= 0 && !applyBefore) {
// the orbit change has not occurred yet, don't change anything
return orbit1;
}
return updateOrbit(orbit1);
}
/** {@inheritDoc} */
public SpacecraftState apply(final SpacecraftState state1) {
if (state1.getDate().compareTo(referenceDate) <= 0 && !applyBefore) {
// the orbit change has not occurred yet, don't change anything
return state1;
}
return new SpacecraftState(updateOrbit(state1.getOrbit()),
state1.getAttitude(), state1.getMass());
}
/** Compute the differential effect of J2 on an orbit.
* @param orbit1 original orbit at t₁, without differential J2
* @return orbit at t₁, always taking the effect into account
*/
private Orbit updateOrbit(final Orbit orbit1) {
// convert current orbital state to equinoctial elements
final EquinoctialOrbit original =
(EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(orbit1);
// compute differential effect
final AbsoluteDate date = original.getDate();
final double dt = date.durationFrom(referenceDate);
final double dPaRaan = (dPaDot + dRaanDot) * dt;
final SinCos scPaRaan = FastMath.sinCos(dPaRaan);
final double dRaan = dRaanDot * dt;
final SinCos scRaan = FastMath.sinCos(dRaan);
final double ex = original.getEquinoctialEx() * scPaRaan.cos() -
original.getEquinoctialEy() * scPaRaan.sin();
final double ey = original.getEquinoctialEx() * scPaRaan.sin() +
original.getEquinoctialEy() * scPaRaan.cos();
final double hx = original.getHx() * scRaan.cos() - original.getHy() * scRaan.sin();
final double hy = original.getHx() * scRaan.sin() + original.getHy() * scRaan.cos();
final double lambda = original.getLv() + dPaRaan;
// build updated orbit
final EquinoctialOrbit updated =
new EquinoctialOrbit(original.getA(), ex, ey, hx, hy, lambda, PositionAngleType.TRUE,
original.getFrame(), date, original.getMu());
// convert to required type
return orbit1.getType().convertType(updated);
}
}