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;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
21  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
22  import org.junit.jupiter.api.Assertions;
23  import org.junit.jupiter.api.BeforeEach;
24  import org.junit.jupiter.api.Test;
25  import org.orekit.Utils;
26  import org.orekit.bodies.CR3BPFactory;
27  import org.orekit.bodies.CR3BPSystem;
28  import org.orekit.errors.OrekitException;
29  import org.orekit.frames.Frame;
30  import org.orekit.propagation.SpacecraftState;
31  import org.orekit.propagation.numerical.NumericalPropagator;
32  import org.orekit.propagation.numerical.cr3bp.CR3BPForceModel;
33  import org.orekit.propagation.numerical.cr3bp.STMEquations;
34  import org.orekit.time.AbsoluteDate;
35  import org.orekit.time.TimeScalesFactory;
36  import org.orekit.utils.AbsolutePVCoordinates;
37  import org.orekit.utils.LagrangianPoints;
38  import org.orekit.utils.PVCoordinates;
39  
40  public class HaloOrbitTest {
41  
42  
43      @Test
44      public void testHaloOrbit() {
45          CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
46          final PVCoordinates firstGuess = new PVCoordinates(new Vector3D(0.0, 1.0, 2.0), new Vector3D(3.0, 4.0, 5.0));
47          final HaloOrbit h1 = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L1), 8E6, LibrationOrbitFamily.NORTHERN);
48          final HaloOrbit h2 = new HaloOrbit(syst, firstGuess, 2.0);
49          final HaloOrbit h3 = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L2), 8E6, LibrationOrbitFamily.SOUTHERN);
50  
51          final double orbitalPeriod1 = h1.getOrbitalPeriod();
52          final double orbitalPeriod2 = h2.getOrbitalPeriod();
53          final double orbitalPeriod3 = h3.getOrbitalPeriod();
54  
55          final PVCoordinates firstGuess1 = h1.getInitialPV();
56          final PVCoordinates firstGuess2 = h2.getInitialPV();
57          final PVCoordinates firstGuess3 = h3.getInitialPV();
58  
59          Assertions.assertNotEquals(0.0, orbitalPeriod1, 0.5);
60          Assertions.assertNotEquals(0.0, orbitalPeriod3, 0.5);
61          Assertions.assertEquals(2.0, orbitalPeriod2, 1E-15);
62  
63          Assertions.assertNotEquals(0.0, firstGuess1.getPosition().getX(), 0.6);
64          Assertions.assertEquals(0.0, firstGuess1.getPosition().getY(), 1E-15);
65          Assertions.assertEquals(0.0, firstGuess1.getVelocity().getX(), 1E-15);
66          Assertions.assertNotEquals(0.0, firstGuess1.getVelocity().getY(), 0.01);
67          Assertions.assertEquals(0.0, firstGuess1.getVelocity().getZ(), 1E-15);
68  
69          Assertions.assertNotEquals(0.0, firstGuess3.getPosition().getX(), 1);
70          Assertions.assertEquals(0.0, firstGuess3.getPosition().getY(), 1E-15);
71          Assertions.assertEquals(0.0, firstGuess3.getVelocity().getX(), 1E-15);
72          Assertions.assertNotEquals(0.0, firstGuess3.getVelocity().getY(), 0.01);
73          Assertions.assertEquals(0.0, firstGuess3.getVelocity().getZ(), 1E-15);
74  
75          Assertions.assertEquals(firstGuess.getPosition().getX(), firstGuess2.getPosition().getX(), 1E-15);
76          Assertions.assertEquals(firstGuess.getPosition().getY(), firstGuess2.getPosition().getY(), 1E-15);
77          Assertions.assertEquals(firstGuess.getPosition().getZ(), firstGuess2.getPosition().getZ(), 1E-15);
78          Assertions.assertEquals(firstGuess.getVelocity().getX(), firstGuess2.getVelocity().getX(), 1E-15);
79          Assertions.assertEquals(firstGuess.getVelocity().getY(), firstGuess2.getVelocity().getY(), 1E-15);
80          Assertions.assertEquals(firstGuess.getVelocity().getZ(), firstGuess2.getVelocity().getZ(), 1E-15);
81      }
82  
83      @Test
84          public void testLagrangianError() {
85          Assertions.assertThrows(OrekitException.class, () -> {
86              CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
87              final HaloOrbit h = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L3), 8E6, LibrationOrbitFamily.NORTHERN);
88              h.getClass();
89          });
90      }
91  
92      @Test
93      public void testManifolds() {
94  
95          // Time settings
96          final AbsoluteDate initialDate =
97              new AbsoluteDate(1996, 06, 25, 0, 0, 00.000,
98                               TimeScalesFactory.getUTC());
99          CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
100 
101         syst.getPrimary().getPVCoordinates(initialDate, syst.getSecondary().getInertiallyOrientedFrame());
102 
103         final Frame Frame = syst.getRotatingFrame();
104 
105         // Define a Northern Halo orbit around Earth-Moon L1 with a Z-amplitude
106         // of 8 000 km
107         HaloOrbit h = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L1), 8E6, LibrationOrbitFamily.SOUTHERN);
108 
109         final double orbitalPeriod = h.getOrbitalPeriod();
110 
111         double integrationTime = orbitalPeriod * 0.9;
112 
113         final PVCoordinates firstGuess = h.getInitialPV();
114 
115         final PVCoordinates initialConditions =
116             new CR3BPDifferentialCorrection(firstGuess, syst, orbitalPeriod).compute(LibrationOrbitType.HALO);
117 
118         final AbsolutePVCoordinates initialAbsPV =
119             new AbsolutePVCoordinates(Frame, initialDate, initialConditions);
120 
121         // Creating the initial spacecraftstate that will be given to the
122         // propagator
123         final SpacecraftState initialState = new SpacecraftState(initialAbsPV);
124 
125         // Integration parameters
126         // These parameters are used for the Dormand-Prince integrator, a
127         // variable step integrator,
128         // these limits prevent the integrator to spend too much time when the
129         // equations are too stiff,
130         // as well as the reverse situation.
131         final double minStep = 1E-10;
132         final double maxstep = 1E-2;
133 
134         // tolerances for integrators
135         // Used by the integrator to estimate its variable integration step
136         final double positionTolerance = 0.0001;
137         final double velocityTolerance = 0.0001;
138         final double massTolerance = 1.0e-6;
139         final double[] vecAbsoluteTolerances = { positionTolerance, positionTolerance, positionTolerance,
140                 velocityTolerance, velocityTolerance, velocityTolerance, massTolerance };
141         final double[] vecRelativeTolerances = new double[vecAbsoluteTolerances.length];
142 
143         // Defining the numerical integrator that will be used by the propagator
144         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, vecAbsoluteTolerances,
145                 vecRelativeTolerances);
146 
147         final STMEquations stm = new STMEquations(syst);
148         final SpacecraftState augmentedInitialState =
149                         stm.setInitialPhi(initialState);
150         NumericalPropagator propagator = new NumericalPropagator(integrator);
151         propagator.setOrbitType(null);
152         propagator.setIgnoreCentralAttraction(true);
153         propagator.addForceModel(new CR3BPForceModel(syst));
154         propagator.addAdditionalDerivativesProvider(stm);
155         propagator.setInitialState(augmentedInitialState);
156         final SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(integrationTime));
157 
158         final PVCoordinates initialUnstableManifold = h.getManifolds(finalState, false);
159         final PVCoordinates initialStableManifold = h.getManifolds(finalState, true);
160 
161         Assertions.assertNotEquals(finalState.getPosition().getX(), initialUnstableManifold.getPosition().getX(), 1E-7);
162         Assertions.assertNotEquals(finalState.getPosition().getY(), initialUnstableManifold.getPosition().getY(), 1E-7);
163         Assertions.assertNotEquals(finalState.getPosition().getZ(), initialUnstableManifold.getPosition().getZ(), 1E-7);
164 
165         Assertions.assertNotEquals(finalState.getPosition().getX(), initialStableManifold.getPosition().getX(), 1E-7);
166         Assertions.assertNotEquals(finalState.getPosition().getY(), initialStableManifold.getPosition().getY(), 1E-7);
167         Assertions.assertNotEquals(finalState.getPosition().getZ(), initialStableManifold.getPosition().getZ(), 1E-7);
168     }
169 
170     @Test
171     public void testDifferentialCorrectionError() {
172         Assertions.assertThrows(OrekitException.class, () -> {
173             CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
174 
175             final double orbitalPeriod = 1;
176 
177             final PVCoordinates firstGuess = new PVCoordinates(new Vector3D(0.0, 1.0, 2.0), new Vector3D(3.0, 4.0, 5.0));
178 
179             final PVCoordinates initialConditions =
180                     new CR3BPDifferentialCorrection(firstGuess, syst, orbitalPeriod).compute(LibrationOrbitType.HALO);
181             initialConditions.toString();
182         });
183     }
184 
185     @Test
186     public void testSTMError() {
187         Assertions.assertThrows(OrekitException.class, () -> {
188             // Time settings
189             final AbsoluteDate initialDate =
190                     new AbsoluteDate(1996, 06, 25, 0, 0, 00.000,
191                             TimeScalesFactory.getUTC());
192             CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
193 
194             final Frame Frame = syst.getRotatingFrame();
195 
196             // Define a Northern Halo orbit around Earth-Moon L1 with a Z-amplitude
197             // of 8 000 km
198             HaloOrbit h = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L1), 8E6, LibrationOrbitFamily.SOUTHERN);
199 
200             final PVCoordinates pv = new PVCoordinates(new Vector3D(0.0, 1.0, 2.0), new Vector3D(3.0, 4.0, 5.0));
201 
202             final AbsolutePVCoordinates initialAbsPV =
203                     new AbsolutePVCoordinates(Frame, initialDate, pv);
204 
205             // Creating the initial spacecraftstate that will be given to the
206             // propagator
207             final SpacecraftState s = new SpacecraftState(initialAbsPV);
208 
209 
210             final PVCoordinates manifold = h.getManifolds(s, true);
211             manifold.getMomentum();
212         });
213     }
214 
215     @BeforeEach
216     public void setUp() {
217         Utils.setDataRoot("cr3bp:regular-data");
218     }
219 }