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.bodies;
18  
19  import org.hipparchus.geometry.euclidean.threed.Rotation;
20  import org.hipparchus.geometry.euclidean.threed.Vector3D;
21  import org.hipparchus.util.FastMath;
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.frames.Frame;
27  import org.orekit.frames.Transform;
28  import org.orekit.time.AbsoluteDate;
29  import org.orekit.time.TimeScale;
30  import org.orekit.time.TimeScalesFactory;
31  import org.orekit.utils.AbsolutePVCoordinates;
32  import org.orekit.utils.LagrangianPoints;
33  import org.orekit.utils.PVCoordinates;
34  
35  public class CR3BPSystemTest {
36  
37      @BeforeEach
38      public void setUp() {
39          Utils.setDataRoot("cr3bp:regular-data");
40      }
41  
42      @Test
43      void testCR3BPSystem() {
44          AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
45          TimeScale timeScale = TimeScalesFactory.getUTC();
46          final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date, timeScale);
47          final double lDim = syst.getDdim();
48          Assertions.assertNotNull(lDim);
49  
50          final double vDim = syst.getVdim();
51          Assertions.assertNotNull(vDim);
52  
53          final double tDim = syst.getTdim();
54          Assertions.assertNotNull(tDim);
55      }
56  
57      @Test
58      void testGetRotatingFrame() {
59          AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
60          TimeScale timeScale = TimeScalesFactory.getUTC();
61          final Frame baryFrame = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getRotatingFrame();
62          Assertions.assertNotNull(baryFrame);
63      }
64  
65      @Test
66      void testGetPrimary() {
67          AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
68          TimeScale timeScale = TimeScalesFactory.getUTC();
69          final CelestialBody primaryBody = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getPrimary();
70          Assertions.assertNotNull(primaryBody);
71      }
72  
73      @Test
74      void testGetSecondary() {
75          AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
76          TimeScale timeScale = TimeScalesFactory.getUTC();
77          final CelestialBody secondaryBody = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getSecondary();
78          Assertions.assertNotNull(secondaryBody);
79      }
80  
81      @Test
82      void testGetMu() {
83          AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
84          TimeScale timeScale = TimeScalesFactory.getUTC();
85          final double mu = CR3BPFactory.getSunJupiterCR3BP(date, timeScale).getMassRatio();
86          Assertions.assertNotNull(mu);
87      }
88  
89      @Test
90      void testGetName() {
91          AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
92          TimeScale timeScale = TimeScalesFactory.getUTC();
93          final String name = CR3BPFactory.getSunEarthCR3BP(date, timeScale).getName();
94          Assertions.assertNotNull(name);
95      }
96  
97      @Test
98      void testGetLPos() {
99          final CR3BPSystem syst = CR3BPFactory.getEarthMoonCR3BP();
100 
101         final Vector3D l1Position = syst.getLPosition(LagrangianPoints.L1);
102         Assertions.assertEquals(3.23E8, l1Position.getX() * syst.getDdim(),3E6);
103         Assertions.assertEquals(0.0, l1Position.getY() * syst.getDdim(),1E3);
104         Assertions.assertEquals(0.0, l1Position.getZ() * syst.getDdim(),1E3);
105 
106         final Vector3D l2Position = syst.getLPosition(LagrangianPoints.L2);
107         Assertions.assertEquals(4.45E8, l2Position.getX() * syst.getDdim(),3E6);
108         Assertions.assertEquals(0.0, l2Position.getY() * syst.getDdim(),1E3);
109         Assertions.assertEquals(0.0, l2Position.getZ() * syst.getDdim(),1E3);
110 
111         final Vector3D l3Position = syst.getLPosition(LagrangianPoints.L3);
112         Assertions.assertEquals(-3.86E8, l3Position.getX() * syst.getDdim(),3E6);
113         Assertions.assertEquals(0.0, l3Position.getY() * syst.getDdim(),1E3);
114         Assertions.assertEquals(0.0, l3Position.getZ() * syst.getDdim(),1E3);
115 
116         final Vector3D l4Position = syst.getLPosition(LagrangianPoints.L4);
117         Assertions.assertEquals(1.87E8, l4Position.getX() * syst.getDdim(),3E6);
118         Assertions.assertEquals(3.32E8, l4Position.getY() * syst.getDdim(),3E6);
119         Assertions.assertEquals(0.0, l4Position.getZ() * syst.getDdim(),1E3);
120 
121         final Vector3D l5Position = syst.getLPosition(LagrangianPoints.L5);
122         Assertions.assertEquals(1.87E8, l5Position.getX() * syst.getDdim(),3E6);
123         Assertions.assertEquals(-3.32E8, l5Position.getY() * syst.getDdim(),3E6);
124         Assertions.assertEquals(0.0, l5Position.getZ() * syst.getDdim(),1E3);
125     }
126 
127     @Test
128     void testGetGamma() {
129         TimeScale timeScale = TimeScalesFactory.getUTC();
130         AbsoluteDate date = new AbsoluteDate(2020, 7, 5, 12, 0, 0.0, timeScale);
131         final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(date, timeScale);
132 
133         final double l1Gamma = syst.getGamma(LagrangianPoints.L1);
134         Assertions.assertEquals(1.497655E9, l1Gamma * syst.getDdim(),1E3);
135 
136 
137         final double l2Gamma = syst.getGamma(LagrangianPoints.L2);
138         Assertions.assertEquals(1.507691E9, l2Gamma * syst.getDdim(),1E3);
139 
140 
141         final double l3Gamma = syst.getGamma(LagrangianPoints.L3);
142         Assertions.assertEquals(1.495981555E11, l3Gamma * syst.getDdim(),9E3);
143     }
144 
145     @Test
146     void testGetRealAPV() {
147 
148         // Time settings
149         final TimeScale timeScale = TimeScalesFactory.getUTC();
150         final AbsoluteDate initialDate =
151             new AbsoluteDate(1996, 06, 26, 0, 0, 00.000,
152                              TimeScalesFactory.getUTC());
153 
154         final CR3BPSystem syst = CR3BPFactory.getSunEarthCR3BP(initialDate, timeScale);
155 
156         final CelestialBody primaryBody = syst.getPrimary();
157         final CelestialBody secondaryBody = syst.getSecondary();
158 
159         final PVCoordinates pv0 = new PVCoordinates(new Vector3D(0,0,1), new Vector3D(0,0,0));
160 
161         final Frame outputFrame = secondaryBody.getInertiallyOrientedFrame();
162 
163         // 1.   Translate the rotating state from the RTBP to a primary-centered rotating state
164         // 2.   Dimensionalize  the  primary-centered  rotating  state  using  the  instantaneously
165         //      defined characteristic quantities
166         // 3.   Apply the transformation matrix
167         // 4.   Apply the transformation to output frame
168 
169         final Frame primaryInertialFrame = primaryBody.getInertiallyOrientedFrame();
170         final PVCoordinates pv21 = secondaryBody.getPVCoordinates(initialDate, primaryInertialFrame);
171 
172         // Distance and Velocity to dimensionalize the state vector
173         final double dist12 = pv21.getPosition().getNorm();
174         final double vCircular  = FastMath.sqrt(secondaryBody.getGM() / dist12);
175 
176         // Dimensionalized state vector centered on primary body
177         final PVCoordinates pvDim = new PVCoordinates((pv0.getPosition().add(new Vector3D(syst.getMassRatio(), 0, 0))).scalarMultiply(dist12),
178                                                       pv0.getVelocity().scalarMultiply(vCircular));
179 
180         // Instantaneous rotation matrix between rotating frame and primary inertial frame
181         final double[][] c = (new Rotation(Vector3D.PLUS_I, Vector3D.PLUS_K,
182                                            pv21.getPosition(), pv21.getMomentum())).getMatrix();
183 
184         // Instantaneous angular velocity of the rotating frame
185         final double theta = pv21.getMomentum().getNorm() / (dist12 * dist12);
186 
187         final double x = pvDim.getPosition().getX();
188         final double y = pvDim.getPosition().getY();
189         final double z = pvDim.getPosition().getZ();
190         final double vx = pvDim.getVelocity().getX();
191         final double vy = pvDim.getVelocity().getY();
192         final double vz = pvDim.getVelocity().getZ();
193 
194         // Position vector in the primary inertial frame
195         final Vector3D newPos = new Vector3D(c[0][0] * x + c[0][1] * y + c[0][2] * z,
196                                              c[1][0] * x + c[1][1] * y + c[1][2] * z,
197                                              c[2][0] * x + c[2][1] * y + c[2][2] * z);
198 
199         final Vector3D vel0 = new Vector3D(c[0][0] * vx + c[0][1] * vy + c[0][2] * vz,
200                                            c[1][0] * vx + c[1][1] * vy + c[1][2] * vz,
201                                            c[2][0] * vx + c[2][1] * vy + c[2][2] * vz);
202         final Vector3D addVel = new Vector3D(c[0][1] * x - c[0][0] * y,
203                                              c[1][1] * x - c[1][0] * y,
204                                              c[2][1] * x - c[2][0] * y);
205         final Vector3D newVel = vel0.add(addVel.scalarMultiply(theta));
206 
207         // State vector in the primary inertial frame
208         final PVCoordinates pv2 = new PVCoordinates(newPos, newVel);
209 
210         // Transformation between primary inertial frame and the output frame
211         final Transform primaryInertialToOutputFrame = primaryInertialFrame.getTransformTo(outputFrame, initialDate);
212 
213         final PVCoordinates pvMat = primaryInertialToOutputFrame.transformPVCoordinates(pv2);
214         final AbsolutePVCoordinates apv0 = new AbsolutePVCoordinates(outputFrame,initialDate,pv0);
215         final AbsolutePVCoordinates apvTrans = syst.getRealAPV(apv0,initialDate,outputFrame);
216 
217         Assertions.assertEquals(pvMat.getPosition().getX(),apvTrans.getPosition().getX(),1E-5);
218         Assertions.assertEquals(pvMat.getPosition().getY(),apvTrans.getPosition().getY(),1E-15);
219         Assertions.assertEquals(pvMat.getPosition().getZ(),apvTrans.getPosition().getZ(),1E-4);
220 
221         Assertions.assertEquals(pvMat.getVelocity().getX(),apvTrans.getVelocity().getX(),1E-2);
222         Assertions.assertEquals(pvMat.getVelocity().getY(),apvTrans.getVelocity().getY(),4E-2);
223         Assertions.assertEquals(pvMat.getVelocity().getZ(),apvTrans.getVelocity().getZ(),2E-2);
224 
225         Assertions.assertEquals(initialDate.durationFrom(AbsoluteDate.J2000_EPOCH),
226                             apvTrans.getDate().durationFrom(AbsoluteDate.J2000_EPOCH),
227                             1E-10);
228     }
229 }