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