RichardsonExpansion.java

  1. /* Copyright 2002-2025 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.orbits;

  19. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.hipparchus.util.FastMath;
  21. import org.hipparchus.util.SinCos;
  22. import org.orekit.bodies.CR3BPSystem;
  23. import org.orekit.errors.OrekitException;
  24. import org.orekit.errors.OrekitMessages;
  25. import org.orekit.utils.LagrangianPoints;
  26. import org.orekit.utils.PVCoordinates;

  28. /** Class implementing the Third-Order Richardson Expansion.
  29.  * @see "Dynamical systems, the three-body problem, and space mission design, Koon, Lo, Marsden, Ross"
  30.  * @author Vincent Mouraux
  31.  * @since 10.2
  32.  */
  33. public class RichardsonExpansion {

  35.     /** Distance between a Lagrangian Point and its closest primary body, meters. */
  36.     private final double gamma;

  38.     /** Mass ratio of the CR3BP System. */
  39.     private final double mu;

  41.     /** Distance between the two primary bodies, meters. */
  42.     private final double dDim;

  44.     /** Halo Orbit frequency. */
  45.     private final double wp;

  47.     /** Simple Halo Orbit coefficient. */
  48.     private final double k;

  50.     /** Simple Richardson coefficient. */
  51.     private final double a21;

  53.     /** Simple Richardson coefficient. */
  54.     private final double a22;

  56.     /** Simple Richardson coefficient. */
  57.     private final double a23;

  59.     /** Simple Richardson coefficient. */
  60.     private final double a24;

  62.     /** Simple Richardson coefficient. */
  63.     private final double b21;

  65.     /** Simple Richardson coefficient. */
  66.     private final double b22;

  68.     /** Simple Richardson coefficient. */
  69.     private final double d21;

  71.     /** Simple Richardson coefficient. */
  72.     private final double a31;

  74.     /** Simple Richardson coefficient. */
  75.     private final double a32;

  77.     /** Simple Richardson coefficient. */
  78.     private final double b31;

  80.     /** Simple Richardson coefficient. */
  81.     private final double b32;

  83.     /** Simple Richardson coefficient. */
  84.     private final double d31;

  86.     /** Simple Richardson coefficient. */
  87.     private final double d32;

  89.     /** Simple Richardson coefficient. */
  90.     private final double s1;

  92.     /** Simple Richardson coefficient. */
  93.     private final double s2;

  95.     /** Simple Richardson coefficient. */
  96.     private final double l1;

  98.     /** Simple Richardson coefficient. */
  99.     private final double l2;

  101.     /** Simple Halo Orbit coefficient. */
  102.     private final double delta;

  104.     /** Orbital Period of the Halo Orbit, seconds. */
  105.     private double orbitalPeriod;

  107.     /** Lagrangian Point considered. */
  108.     private LagrangianPoints point;

  110.     /** CR3BP System considered. */
  111.     private CR3BPSystem cr3bpSystem;

  113.     /** Simple Constructor.
  114.      * @param cr3bpSystem CR3BP System considered
  115.      * @param point Lagrangian Point considered
  116.     */
  117.     public RichardsonExpansion(final CR3BPSystem cr3bpSystem,
  118.                                       final LagrangianPoints point) {

  120.         this.point       = point;
  121.         this.cr3bpSystem = cr3bpSystem;
  122.         this.mu          = cr3bpSystem.getMassRatio();
  123.         this.dDim        = cr3bpSystem.getDdim();
  124.         this.gamma       = cr3bpSystem.getGamma(point);

  126.         final double c2 = getCn(2.0);
  127.         final double c3 = getCn(3.0);
  128.         final double c4 = getCn(4.0);

  130.         wp = FastMath.sqrt(0.5 * (2.0 - c2 + FastMath.sqrt(9.0 * c2 * c2 - 8.0 * c2)));

  132.         final double ld = FastMath.sqrt(0.5 * (2.0 - c2 + FastMath.sqrt(9.0 * c2 * c2 - 8.0 * c2)));

  134.         k = (wp * wp + 1.0 + 2.0 * c2) / (2.0 * wp);

  136.         final double d1 = 3.0 * ld * ld * (k * (6.0 * ld * ld - 1.0) - 2.0 * ld) / k;
  137.         final double d2 = 8.0 * ld * ld * (k * (11.0 * ld * ld - 1.0) - 2.0 * ld) / k;

  139.         a21 = 3.0 * c3 * (k * k - 2.0) / (4.0 * (1.0 + 2.0 * c2));

  141.         a22 = 3.0 * c3 / (4.0 * (1.0 + 2.0 * c2));

  143.         a23 = -3.0 * c3 * ld * (3.0 * k * k * k * ld - 6.0 * k * (k - ld) + 4.0) / (4.0 * k * d1);

  145.         a24 = -3.0 * c3 * ld * (2.0 + 3.0 * k * ld) / (4.0 * k * d1);

  147.         b21 = -3.0 * c3 * ld * (3.0 * k * ld - 4.0) / (2.0 * d1);

  149.         b22 = 3.0 * c3 * ld / d1;

  151.         d21 = -c3 / (2.0 * ld * ld);

  153.         a31 =
  154.             -9.0 * ld * (4.0 * c3 * (k * a23 - b21) + k * c4 * (4.0 + k * k)) / (4.0 * d2) +
  155.             (9.0 * ld * ld + 1.0 - c2) / (2.0 * d2) * (2.0 * c3 * (2.0 * a23 - k * b21) + c4 * (2.0 + 3.0 * k * k));

  157.         a32 =
  158.             -9.0 * ld * (4.0 * c3 * (k * a24 - b22) + k * c4) / (4.0 * d2) -
  159.             3 * (9.0 * ld * ld + 1.0 - c2) * (c3 * (k * b22 + d21 - 2.0 * a24) - c4) / (2.0 * d2);

  161.         b31 =
  162.             3.0 * 8.0 * ld * (3.0 * c3 * (k * b21 - 2.0 * a23) - c4 * (2.0 + 3.0 * k * k)) / (8.0 * d2) +
  163.             3.0 * ((9.0 * ld * ld + 1.0 + 2.0 * c2) * (4.0 * c3 * (k * a23 - b21) + k * c4 * (4.0 + k * k))) / (8.0 * d2);

  165.         b32 =
  166.             9.0 * ld * (c3 * (k * b22 + d21 - 2.0 * a24) - c4) / d2 +
  167.             3.0 * ((9.0 * ld * ld + 1.0 + 2.0 * c2) * (4.0 * c3 * (k * a24 - b22) +  k * c4)) / (8.0 * d2);

  169.         d31 = 3.0 / (64.0 * ld * ld) * (4.0 * c3 * a24 + c4);

  171.         d32 = 3.0 / (64.0 * ld * ld) * (4.0 * c3 * (a23 - d21) + c4 * (4.0 + k * k));

  173.         s1 =
  174.             (3.0 * c3 * (2.0 * a21 * (k * k - 2.0) - a23 * (k * k + 2.0) - 2.0 * k * b21) / 2.0 -
  175.             3.0 * c4 * (3.0 * k * k * k * k - 8.0 * k * k + 8.0) / 8.0) / (2.0 * ld * (ld * (1.0 + k * k) - 2.0 * k));

  177.         s2 =
  178.             (3.0 * c3 * (2.0 * a22 * (k * k - 2.0) + a24 * (k * k + 2.0) + 2.0 * k * b22 + 5.0 * d21) / 2.0 +
  179.             3.0 * c4 * (12.0 - k * k) / 8.0) / (2.0 * ld * (ld * (1.0 + k * k) - 2.0 * k));

  181.         l1 = -3.0 * c3 * (2.0 * a21 + a23 + 5.0 * d21) / 2.0 - 3.0 * c4 * (12.0 - k * k) / 8.0 + 2.0 * ld * ld * s1;

  183.         l2 = 3.0 * c3 * (a24 - 2.0 * a22) / 2.0 + (9.0 * c4 / 8.0) + (2.0 * ld * ld * s2);

  185.         delta = wp * wp - c2;
  186.     }


  189.     /** Calculate Cn Richardson Coefficient.
  190.      * @param order Order 'n' of the coefficient needed.
  191.      * @return cn Cn Richardson Coefficient
  192.     */
  193.     private double getCn(final double order) {
  194.         final double gamma3 = gamma * gamma * gamma;
  195.         final double cn;
  196.         switch (point) {
  197.             case L1:
  198.                 cn = (1.0 / gamma3) * (mu + FastMath.pow(-1, order) * (1 - mu) * FastMath.pow(gamma, order + 1) / FastMath.pow(1 - gamma, order + 1));
  199.                 break;
  200.             case L2:
  201.                 cn = (1.0 / gamma3) * (FastMath.pow(-1, order) * mu + FastMath.pow(-1, order) * (1 - mu) * FastMath.pow(gamma, order + 1) / FastMath.pow(1 + gamma, order + 1));
  202.                 break;
  203.             default:
  204.                 throw new OrekitException(OrekitMessages.CANNOT_COMPUTE_LAGRANGIAN, point);
  205.         }
  206.         return cn;
  207.     }

  209.     /** Calculate first Guess.
  210.      * @param azr z-axis Amplitude of the required Halo Orbit, meters
  211.      * @param type type of the Halo Orbit ("Northern" or "Southern")
  212.      * @param t Orbit time, seconds (must be greater than 0)
  213.      * @param phi Orbit phase, rad
  214.      * @return firstGuess PVCoordinates of the first guess
  215.     */
  216.     public PVCoordinates computeHaloFirstGuess(final double azr, final LibrationOrbitFamily type,
  217.                                            final double t, final double phi) {

  219.         // Z-Axis Halo Orbit Amplitude
  220.         final double az = azr / (gamma * dDim);

  222.         // X-Axis Halo Orbit Amplitude
  223.         final double ax      = FastMath.sqrt((delta + l2 * az * az) / -l1);
  224.         final double nu      = 1.0 + s1 * ax * ax + s2 * az * az;
  225.         final double tau     = nu * t;
  226.         final double tau1    = wp * tau + phi;
  227.         final double m       = (type == LibrationOrbitFamily.NORTHERN) ? 1 : 3;
  228.         final double dm      = 2 - m;
  229.         final SinCos scTau1  = FastMath.sinCos(tau1);
  230.         final SinCos sc2Tau1 = SinCos.sum(scTau1, scTau1);
  231.         final SinCos sc3Tau1 = SinCos.sum(scTau1, sc2Tau1);

  233.         // First guess position relative to its Lagrangian point
  234.         final double firstx =
  235.                         a21 * ax * ax +
  236.                         a22 * az * az - ax * scTau1.cos() +
  237.                         (a23 * ax * ax - a24 * az * az) * sc2Tau1.cos() +
  238.                         (a31 * ax * ax * ax - a32 * ax * az * az) * sc3Tau1.cos();

  240.         final double firsty =
  241.                         k * ax * scTau1.sin() +
  242.                         (b21 * ax * ax - b22 * az * az) * sc2Tau1.sin() +
  243.                         (b31 * ax * ax * ax - b32 * ax * az * az) * sc3Tau1.sin();

  245.         final double firstz =
  246.                         dm * az * scTau1.cos() +
  247.                         dm * d21 * ax * az * (sc2Tau1.cos() - 3.0) +
  248.                         dm * (d32 * az * ax * ax - d31 * az * az * az) * sc3Tau1.cos();

  250.         // First guess Velocity relative to its Lagrangian point
  251.         final double vx =
  252.                         ax * wp * nu * scTau1.sin() -
  253.                         2.0 * (a23 * ax * ax - a24 * az * az) * wp * nu * sc2Tau1.sin() -
  254.                         3.0 * (a31 * ax * ax * ax - a32 * ax * az * az) * wp * nu * sc3Tau1.sin();

  256.         final double vy =
  257.                         k * ax * wp * nu * scTau1.cos() +
  258.                         2.0 * wp * nu * (b21 * ax * ax - b22 * az * az) * sc2Tau1.cos() +
  259.                         3.0 * wp * nu * (b31 * ax * ax * ax - b32 * ax * az * az) * sc3Tau1.cos();

  261.         final double vz =
  262.                         -dm * az * wp * nu * scTau1.sin() -
  263.                         2.0 * dm * d21 * ax * az * wp * nu * sc2Tau1.sin() -
  264.                         3.0 * dm * (d32 * az * ax * ax - d31 * az * az * az) * wp * nu * sc3Tau1.sin();

  266.         // Return PV Coordinates
  267.         return point == LagrangianPoints.L1 ?
  268.                 new PVCoordinates(new Vector3D(firstx * gamma + 1.0 - mu - gamma, firsty * gamma, firstz * gamma), new Vector3D(vx * gamma, vy * gamma, vz * gamma)) :
  269.                 new PVCoordinates(new Vector3D(firstx * gamma + 1.0 - mu + gamma, firsty * gamma, firstz * gamma), new Vector3D(vx * gamma, vy * gamma, vz * gamma));
  270.     }

  272.     /** Calculate first Guess.
  273.      * @param ayr x-axis Amplitude of the required Lyapunov Orbit, meters
  274.      * @param t time
  275.      * @param phi Orbit phase, rad
  276.      * @return firstGuess PVCoordinates of the first guess
  277.     */
  278.     public PVCoordinates computeLyapunovFirstGuess(final double ayr, final double t, final double phi) {

  280.         // Z-Axis Lyapunov Orbit Amplitude
  281.         final double az = 0;

  283.         // y-Axis Lyapunov Orbit Amplitude
  284.         final double ay   = ayr / (gamma * dDim);
  285.         final double ax   = ay / k;
  286.         final double nu   = 1.0 + s1 * ax * ax + s2 * az * az;
  287.         final double tau  = nu * t;
  288.         final double tau1 = wp * tau + phi;

  290.         // First guess position relative to its Lagrangian point
  291.         final double firstx = -ax * FastMath.cos(tau1);
  292.         final double firsty = 0.0;
  293.         final double firstz = 0.0;

  295.         // First guess Velocity relative to its Lagrangian point
  296.         final double vx = 0.0;
  297.         final double vy = k * ax * wp * nu * FastMath.cos(tau1);
  298.         final double vz = 0.0;

  300.         // Return PV Coordinates
  301.         return point == LagrangianPoints.L1 ?
  302.                  new PVCoordinates(new Vector3D(firstx * gamma + 1.0 - mu - gamma, firsty * gamma, firstz * gamma),
  303.                                    new Vector3D(vx * gamma, vy * gamma, vz * gamma)) :
  304.                  new PVCoordinates(new Vector3D(firstx * gamma + 1.0 - mu + gamma, firsty * gamma, firstz * gamma),
  305.                                    new Vector3D(vx * gamma, vy * gamma, vz * gamma));
  306.     }

  308.     /** Return the orbital period of the Halo Orbit.
  309.      * @param azr z-axis Amplitude of the required Halo Orbit, meters
  310.      * @return the orbitalPeriod
  311.      */
  312.     public double getHaloOrbitalPeriod(final double azr) {
  313.         final double az = azr / (gamma * dDim);
  314.         final double ax = FastMath.sqrt((delta + l2 * az * az) / -l1);
  315.         final double nu = 1.0 + s1 * ax * ax + s2 * az * az;
  316.         orbitalPeriod   = FastMath.abs(2.0 * FastMath.PI / (wp * nu));
  317.         return orbitalPeriod;
  318.     }

  320.     /** Return the orbital period of the Halo Orbit.
  321.      * @param axr x-axis Amplitude of the required Lyapunov Orbit, meters
  322.      * @return the orbitalPeriod
  323.      */
  324.     public double getLyapunovOrbitalPeriod(final double axr) {
  325.         final double az = 0;
  326.         final double ax = axr / (gamma * dDim);
  327.         final double nu = 1.0 + s1 * ax * ax + s2 * az * az;
  328.         orbitalPeriod   = FastMath.abs(2.0 * FastMath.PI / (wp * nu));
  329.         return orbitalPeriod;
  330.     }

  332.     /** Get the considered CR3BP system.
  333.      * @return CRR3BP system
  334.      */
  335.     public CR3BPSystem getCr3bpSystem() {
  336.         return cr3bpSystem;
  337.     }

  339.     /** Get the considered lagrangian point.
  340.      * @return lagrangian point
  341.      */
  342.     public LagrangianPoints getLagrangianPoint() {
  343.         return point;
  344.     }

  346. }
Created with JaCoCo 0.8.12.202403310830