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.propagation.numerical.cr3bp;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
24  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25  import org.junit.Assert;
26  import org.junit.Before;
27  import org.junit.Test;
28  import org.orekit.Utils;
29  import org.orekit.bodies.CR3BPFactory;
30  import org.orekit.bodies.CR3BPSystem;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.frames.Frame;
33  import org.orekit.orbits.CR3BPDifferentialCorrection;
34  import org.orekit.orbits.HaloOrbit;
35  import org.orekit.orbits.LibrationOrbitFamily;
36  import org.orekit.orbits.LibrationOrbitType;
37  import org.orekit.orbits.RichardsonExpansion;
38  import org.orekit.propagation.SpacecraftState;
39  import org.orekit.propagation.numerical.NumericalPropagator;
40  import org.orekit.time.AbsoluteDate;
41  import org.orekit.time.TimeScalesFactory;
42  import org.orekit.utils.AbsolutePVCoordinates;
43  import org.orekit.utils.LagrangianPoints;
44  import org.orekit.utils.PVCoordinates;
45  
46  public class CR3BPMultipleShooterTest {
47  
48      @Test
49      public void testHaloOrbit() {
50          
51          final CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
52          final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
53          final HaloOrbit h1 = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L1), 8E6, LibrationOrbitFamily.NORTHERN);
54  
55          // Adaptive stepsize boundaries
56          final double minStep = 1E-12;
57          final double maxstep = 0.001;
58  
59          // Integrator tolerances
60          final double positionTolerance = 1E-5;
61          final double velocityTolerance = 1E-5;
62          final double massTolerance = 1.0e-6;
63          final double[] vecAbsoluteTolerances = {positionTolerance, positionTolerance, positionTolerance, velocityTolerance, velocityTolerance, velocityTolerance, massTolerance };
64          final double[] vecRelativeTolerances =
65                          new double[vecAbsoluteTolerances.length];
66  
67          // Integrator definition
68          final AdaptiveStepsizeIntegrator integrator =
69                          new DormandPrince853Integrator(minStep, maxstep,
70                                                         vecAbsoluteTolerances,
71                                                         vecRelativeTolerances);
72          final int narcs = 1;
73          final List<STMEquations> cr3bpAdditionalEquations = new ArrayList<>(narcs) ;
74          cr3bpAdditionalEquations.add(new STMEquations(syst));
75  
76          // Propagator definition for CR3BP
77          final List<NumericalPropagator> propagatorList = new ArrayList<NumericalPropagator>(narcs);
78          final NumericalPropagator propagator = new NumericalPropagator(integrator);
79          propagator.setOrbitType(null);
80          propagator.setIgnoreCentralAttraction(true);
81          propagator.addForceModel(new CR3BPForceModel(syst));
82  
83          // Add new set of additional equations to the propagator
84          propagator.addAdditionalDerivativesProvider(cr3bpAdditionalEquations.get(0));
85  
86          propagatorList.add(propagator);
87  
88          // First guess trajectory definition
89  
90          final PVCoordinates firstGuess1 = h1.getInitialPV();
91          final PVCoordinates firstGuess2 = new PVCoordinates(firstGuess1.getPosition().add(new Vector3D(0.1,0,0)), firstGuess1.getVelocity().negate());
92          final double arcDuration = h1.getOrbitalPeriod()/2;
93  
94  
95          List<SpacecraftState> firstGuessList = new ArrayList<SpacecraftState>(narcs + 1) ;;
96          firstGuessList.add(new SpacecraftState(new AbsolutePVCoordinates(syst.getRotatingFrame(),
97                                                                           date,
98                                                                           firstGuess1)));
99          firstGuessList.add(new SpacecraftState(new AbsolutePVCoordinates(syst.getRotatingFrame(),
100                                                                          date.shiftedBy(arcDuration),
101                                                                          firstGuess2)));
102 
103         // Multiple Shooting definition
104         final CR3BPMultipleShooter multipleShooting = new CR3BPMultipleShooter(firstGuessList, propagatorList, cr3bpAdditionalEquations, arcDuration, 1E-8, 1);
105         multipleShooting.setPatchPointComponentFreedom(1, 1, false);
106         multipleShooting.setPatchPointComponentFreedom(1, 2, false);
107         multipleShooting.setPatchPointComponentFreedom(1, 3, false);
108         multipleShooting.setPatchPointComponentFreedom(1, 5, false);
109         multipleShooting.setPatchPointComponentFreedom(2, 1, false);
110         multipleShooting.setPatchPointComponentFreedom(2, 3, false);
111         multipleShooting.setPatchPointComponentFreedom(2, 5, false);
112         multipleShooting.setEpochFreedom(1, false);
113         multipleShooting.setEpochFreedom(2, false);
114         multipleShooting.addConstraint(1, 1, 1.0e-5);
115         
116         // Differential correction
117         h1.applyDifferentialCorrection();
118         final PVCoordinates initialPVDC = h1.getInitialPV();
119         final double periodDC = h1.getOrbitalPeriod();
120 
121         // Multiple shooting computation
122         List<SpacecraftState> result = multipleShooting.compute();
123         final AbsolutePVCoordinates initialPVMS = result.get(0).getAbsPVA();
124         final double periodMS = 2 * result.get(1).getDate().durationFrom(result.get(0).getDate()); 
125 
126         Assert.assertEquals(0.0, initialPVDC.getPosition().getY(), 1E-15);
127         Assert.assertEquals(0.0, initialPVDC.getVelocity().getX(), 1E-15);
128         Assert.assertEquals(0.0, initialPVDC.getVelocity().getZ(), 1E-15);
129 
130         Assert.assertEquals(0.0, initialPVMS.getPosition().getY(), 1E-15);
131         Assert.assertEquals(0.0, initialPVMS.getVelocity().getX(), 1E-15);
132         Assert.assertEquals(0.0, initialPVMS.getVelocity().getZ(), 1E-15);
133 
134         Assert.assertEquals(initialPVDC.getPosition().getX(), initialPVMS.getPosition().getX(), 6.6E-4);
135         Assert.assertEquals(initialPVDC.getPosition().getZ(), initialPVMS.getPosition().getZ(), 1.0E-15);
136         Assert.assertEquals(initialPVDC.getVelocity().getY(), initialPVMS.getVelocity().getY(), 7.2E-3);
137 
138         Assert.assertEquals(periodDC, periodMS, 3.0E-2);
139     }
140 
141     @Test(expected=OrekitException.class)
142     public void testLagrangianError() {
143         CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
144         final HaloOrbit h = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L3), 8E6, LibrationOrbitFamily.NORTHERN);
145         h.getClass();
146     }
147 
148     @Test(expected=OrekitException.class)
149     public void testDifferentialCorrectionError() {
150 
151         CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
152 
153         final double orbitalPeriod = 1;
154 
155         final PVCoordinates firstGuess = new PVCoordinates(new Vector3D(0.0, 1.0, 2.0), new Vector3D(3.0, 4.0, 5.0));
156 
157         final PVCoordinates initialConditions =
158                         new CR3BPDifferentialCorrection(firstGuess, syst, orbitalPeriod).compute(LibrationOrbitType.HALO);
159         initialConditions.toString();
160     }
161 
162     @Test(expected=OrekitException.class)
163     public void testSTMError() {
164         // Time settings
165         final AbsoluteDate initialDate =
166                         new AbsoluteDate(1996, 06, 25, 0, 0, 00.000,
167                                          TimeScalesFactory.getUTC());
168         CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
169 
170         final Frame Frame = syst.getRotatingFrame();
171 
172         // Define a Northern Halo orbit around Earth-Moon L1 with a Z-amplitude
173         // of 8 000 km
174         HaloOrbit h = new HaloOrbit(new RichardsonExpansion(syst, LagrangianPoints.L1), 8E6, LibrationOrbitFamily.SOUTHERN);
175 
176         final PVCoordinates pv = new PVCoordinates(new Vector3D(0.0, 1.0, 2.0), new Vector3D(3.0, 4.0, 5.0));
177 
178         final AbsolutePVCoordinates initialAbsPV =
179                         new AbsolutePVCoordinates(Frame, initialDate, pv);
180 
181         // Creating the initial spacecraftstate that will be given to the
182         // propagator
183         final SpacecraftState s = new SpacecraftState(initialAbsPV);
184 
185 
186         final PVCoordinates manifold = h.getManifolds(s, true);
187         manifold.getMomentum();
188     }
189 
190     @Before
191     public void setUp() {
192         Utils.setDataRoot("cr3bp:regular-data");
193     }
194 }