J2DifferentialEffect.java

  1. /* Copyright 2002-2024 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.propagation.analytical;

  18. import org.hipparchus.util.FastMath;
  19. import org.hipparchus.util.SinCos;
  20. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  21. import org.orekit.orbits.EquinoctialOrbit;
  22. import org.orekit.orbits.Orbit;
  23. import org.orekit.orbits.OrbitType;
  24. import org.orekit.orbits.PositionAngleType;
  25. import org.orekit.propagation.SpacecraftState;
  26. import org.orekit.time.AbsoluteDate;

  27. /** Analytical model for J2 effect.
  28.  * <p>
  29.  * This class computes the differential effect of J2 due to an initial orbit
  30.  * offset. A typical case is when an inclination maneuver changes an orbit
  31.  * inclination at time t₀. As ascending node drift rate depends on
  32.  * inclination, the change induces a time-dependent change in ascending node
  33.  * for later dates.
  34.  * </p>
  35.  * @see org.orekit.forces.maneuvers.SmallManeuverAnalyticalModel
  36.  * @author Luc Maisonobe
  37.  */
  38. public class J2DifferentialEffect
  39.     implements AdapterPropagator.DifferentialEffect {

  40.     /** Reference date. */
  41.     private final AbsoluteDate referenceDate;

  42.     /** Differential drift on perigee argument. */
  43.     private final double dPaDot;

  44.     /** Differential drift on ascending node. */
  45.     private final double dRaanDot;

  46.     /** Indicator for applying effect before reference date. */
  47.     private final boolean applyBefore;

  48.     /** Simple constructor.
  49.      * <p>
  50.      * The {@code applyBefore} parameter is mainly used when the differential
  51.      * effect is associated with a maneuver. In this case, the parameter must be
  52.      * set to {@code false}.
  53.      * </p>
  54.      * @param original original state at reference date
  55.      * @param directEffect direct effect changing the orbit
  56.      * @param applyBefore if true, effect is applied both before and after
  57.      * reference date, if false it is only applied after reference date
  58.      * @param gravityField gravity field to use
  59.      */
  60.     public J2DifferentialEffect(final SpacecraftState original,
  61.                                 final AdapterPropagator.DifferentialEffect directEffect,
  62.                                 final boolean applyBefore,
  63.                                 final UnnormalizedSphericalHarmonicsProvider gravityField) {
  64.         this(original, directEffect, applyBefore,
  65.              gravityField.getAe(), gravityField.getMu(),
  66.              -gravityField.onDate(original.getDate()).getUnnormalizedCnm(2, 0));
  67.     }

  68.     /** Simple constructor.
  69.          * <p>
  70.          * The {@code applyBefore} parameter is mainly used when the differential
  71.          * effect is associated with a maneuver. In this case, the parameter must be
  72.          * set to {@code false}.
  73.          * </p>
  74.          * @param orbit0 original orbit at reference date
  75.          * @param orbit1 shifted orbit at reference date
  76.          * @param applyBefore if true, effect is applied both before and after
  77.          * reference date, if false it is only applied after reference date
  78.          * @param gravityField gravity field to use
  79.          */
  80.     public J2DifferentialEffect(final Orbit orbit0, final Orbit orbit1, final boolean applyBefore,
  81.                                 final UnnormalizedSphericalHarmonicsProvider gravityField) {
  82.         this(orbit0, orbit1, applyBefore,
  83.              gravityField.getAe(), gravityField.getMu(),
  84.              -gravityField.onDate(orbit0.getDate()).getUnnormalizedCnm(2, 0));
  85.     }

  86.     /** Simple constructor.
  87.      * <p>
  88.      * The {@code applyBefore} parameter is mainly used when the differential
  89.      * effect is associated with a maneuver. In this case, the parameter must be
  90.      * set to {@code false}.
  91.      * </p>
  92.      * @param original original state at reference date
  93.      * @param directEffect direct effect changing the orbit
  94.      * @param applyBefore if true, effect is applied both before and after
  95.      * reference date, if false it is only applied after reference date
  96.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  97.      * @param mu central attraction coefficient (m³/s²)
  98.      * @param j2 un-normalized zonal coefficient (about +1.08e-3 for Earth)
  99.      */
  100.     public J2DifferentialEffect(final SpacecraftState original,
  101.                                 final AdapterPropagator.DifferentialEffect directEffect,
  102.                                 final boolean applyBefore,
  103.                                 final double referenceRadius, final double mu, final double j2) {
  104.         this(original.getOrbit(),
  105.              directEffect.apply(original.shiftedBy(0.001)).getOrbit().shiftedBy(-0.001),
  106.              applyBefore, referenceRadius, mu, j2);
  107.     }

  108.     /** Simple constructor.
  109.      * <p>
  110.      * The {@code applyBefore} parameter is mainly used when the differential
  111.      * effect is associated with a maneuver. In this case, the parameter must be
  112.      * set to {@code false}.
  113.      * </p>
  114.      * @param orbit0 original orbit at reference date
  115.      * @param orbit1 shifted orbit at reference date
  116.      * @param applyBefore if true, effect is applied both before and after
  117.      * reference date, if false it is only applied after reference date
  118.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  119.      * @param mu central attraction coefficient (m³/s²)
  120.      * @param j2 un-normalized zonal coefficient (about +1.08e-3 for Earth)
  121.      */
  122.     public J2DifferentialEffect(final Orbit orbit0, final Orbit orbit1, final boolean applyBefore,
  123.                                 final double referenceRadius, final double mu, final double j2) {

  124.         this.referenceDate = orbit0.getDate();
  125.         this.applyBefore   = applyBefore;

  126.         // extract useful parameters
  127.         final double a0 = orbit0.getA();
  128.         final double e0 = orbit0.getE();
  129.         final double i0 = orbit0.getI();
  130.         final double a1 = orbit1.getA();
  131.         final double e1 = orbit1.getE();
  132.         final double i1 = orbit1.getI();

  133.         // compute reference drifts
  134.         final double oMe2       = 1 - e0 * e0;
  135.         final double ratio      = referenceRadius / (a0 * oMe2);
  136.         final SinCos scI        = FastMath.sinCos(i0);
  137.         final double n          = FastMath.sqrt(mu / a0) / a0;
  138.         final double c          = ratio * ratio * n * j2;
  139.         final double refPaDot   =  0.75 * c * (4 - 5 * scI.sin() * scI.sin());
  140.         final double refRaanDot = -1.5  * c * scI.cos();

  141.         // differential model on perigee argument drift
  142.         final double dPaDotDa = -3.5 * refPaDot / a0;
  143.         final double dPaDotDe = 4 * refPaDot * e0 / oMe2;
  144.         final double dPaDotDi = -7.5 * c * scI.sin() * scI.cos();
  145.         dPaDot = dPaDotDa * (a1 - a0) + dPaDotDe * (e1 - e0) + dPaDotDi * (i1 - i0);

  146.         // differential model on ascending node drift
  147.         final double dRaanDotDa = -3.5 * refRaanDot / a0;
  148.         final double dRaanDotDe = 4 * refRaanDot * e0 / oMe2;
  149.         final double dRaanDotDi = -refRaanDot * FastMath.tan(i0);
  150.         dRaanDot = dRaanDotDa * (a1 - a0) + dRaanDotDe * (e1 - e0) + dRaanDotDi * (i1 - i0);

  151.     }

  152.     /** Compute the effect of the maneuver on an orbit.
  153.      * @param orbit1 original orbit at t₁, without maneuver
  154.      * @return orbit at t₁, taking the maneuver
  155.      * into account if t₁ &gt; t₀
  156.      * @see #apply(SpacecraftState)
  157.      */
  158.     public Orbit apply(final Orbit orbit1) {

  159.         if (orbit1.getDate().compareTo(referenceDate) <= 0 && !applyBefore) {
  160.             // the orbit change has not occurred yet, don't change anything
  161.             return orbit1;
  162.         }

  163.         return updateOrbit(orbit1);

  164.     }

  165.     /** {@inheritDoc} */
  166.     public SpacecraftState apply(final SpacecraftState state1) {

  167.         if (state1.getDate().compareTo(referenceDate) <= 0 && !applyBefore) {
  168.             // the orbit change has not occurred yet, don't change anything
  169.             return state1;
  170.         }

  171.         return new SpacecraftState(updateOrbit(state1.getOrbit()),
  172.                                    state1.getAttitude(), state1.getMass());

  173.     }

  174.     /** Compute the differential effect of J2 on an orbit.
  175.      * @param orbit1 original orbit at t₁, without differential J2
  176.      * @return orbit at t₁, always taking the effect into account
  177.      */
  178.     private Orbit updateOrbit(final Orbit orbit1) {

  179.         // convert current orbital state to equinoctial elements
  180.         final EquinoctialOrbit original =
  181.                 (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(orbit1);

  182.         // compute differential effect
  183.         final AbsoluteDate date = original.getDate();
  184.         final double dt         = date.durationFrom(referenceDate);
  185.         final double dPaRaan    = (dPaDot + dRaanDot) * dt;
  186.         final SinCos scPaRaan   = FastMath.sinCos(dPaRaan);
  187.         final double dRaan      = dRaanDot * dt;
  188.         final SinCos scRaan     = FastMath.sinCos(dRaan);

  189.         final double ex         = original.getEquinoctialEx() * scPaRaan.cos() -
  190.                                   original.getEquinoctialEy() * scPaRaan.sin();
  191.         final double ey         = original.getEquinoctialEx() * scPaRaan.sin() +
  192.                                   original.getEquinoctialEy() * scPaRaan.cos();
  193.         final double hx         = original.getHx() * scRaan.cos() - original.getHy() * scRaan.sin();
  194.         final double hy         = original.getHx() * scRaan.sin() + original.getHy() * scRaan.cos();
  195.         final double lambda     = original.getLv() + dPaRaan;

  196.         // build updated orbit
  197.         final EquinoctialOrbit updated =
  198.                 new EquinoctialOrbit(original.getA(), ex, ey, hx, hy, lambda, PositionAngleType.TRUE,
  199.                                      original.getFrame(), date, original.getMu());

  200.         // convert to required type
  201.         return orbit1.getType().convertType(updated);

  202.     }

  203. }