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.propagation.semianalytical.dsst.forces;
18  
19  import org.hipparchus.geometry.euclidean.threed.Rotation;
20  import org.hipparchus.geometry.euclidean.threed.RotationOrder;
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
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.attitudes.Attitude;
27  import org.orekit.attitudes.AttitudeProvider;
28  import org.orekit.attitudes.FrameAlignedProvider;
29  import org.orekit.attitudes.LofOffset;
30  import org.orekit.bodies.CelestialBody;
31  import org.orekit.bodies.CelestialBodyFactory;
32  import org.orekit.bodies.OneAxisEllipsoid;
33  import org.orekit.forces.BoxAndSolarArraySpacecraft;
34  import org.orekit.frames.Frame;
35  import org.orekit.frames.FramesFactory;
36  import org.orekit.frames.LOFType;
37  import org.orekit.orbits.EquinoctialOrbit;
38  import org.orekit.orbits.Orbit;
39  import org.orekit.orbits.PositionAngleType;
40  import org.orekit.propagation.PropagationType;
41  import org.orekit.propagation.SpacecraftState;
42  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
43  import org.orekit.time.AbsoluteDate;
44  import org.orekit.time.DateComponents;
45  import org.orekit.time.TimeComponents;
46  import org.orekit.time.TimeScalesFactory;
47  import org.orekit.utils.Constants;
48  import org.orekit.utils.IERSConventions;
49  import org.orekit.utils.TimeStampedAngularCoordinates;
50  
51  import java.io.IOException;
52  import java.text.ParseException;
53  import java.util.ArrayList;
54  import java.util.Arrays;
55  import java.util.List;
56  
57  class DSSTSolarRadiationPressureTest {
58  
59      @Test
60      void testGetMeanElementRate() throws IllegalArgumentException {
61  
62          final Frame earthFrame = FramesFactory.getGCRF();
63          final AbsoluteDate initDate = new AbsoluteDate(2003, 9, 16, 0, 0, 0, TimeScalesFactory.getUTC());
64          final double mu = 3.986004415E14;
65          // a  = 42166258 m
66          // ex = 6.532127416888538E-6
67          // ey = 9.978642849310487E-5
68          // hx = -5.69711879850274E-6
69          // hy = 6.61038518895005E-6
70          // lM = 8.56084687583949 rad
71          final Orbit orbit = new EquinoctialOrbit(4.2166258E7,
72                                                   6.532127416888538E-6,
73                                                   9.978642849310487E-5,
74                                                   -5.69711879850274E-6,
75                                                   6.61038518895005E-6,
76                                                   8.56084687583949,
77                                                   PositionAngleType.TRUE,
78                                                   earthFrame,
79                                                   initDate,
80                                                   mu);
81  
82          // SRP Force Model
83          DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
84                                                              new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
85                                                                                   Constants.WGS84_EARTH_FLATTENING,
86                                                                                   FramesFactory.getITRF(IERSConventions.IERS_2010, false)),
87                                                              mu);
88          // Attitude of the satellite
89          Rotation rotation =  new Rotation(0.9999999999999984,
90                                            1.653020584550675E-8,
91                                            -4.028108631990782E-8,
92                                            -3.539139805514139E-8,
93                                            false);
94          Vector3D rotationRate = new Vector3D(0., 0., 0.);
95          Vector3D rotationAcceleration = new Vector3D(0., 0., 0.);
96          TimeStampedAngularCoordinates orientation = new TimeStampedAngularCoordinates(initDate,
97                                                                                        rotation,
98                                                                                        rotationRate,
99                                                                                        rotationAcceleration);
100         final Attitude att = new Attitude(earthFrame, orientation);
101 
102         // Spacecraft state
103         final SpacecraftState state = new SpacecraftState(orbit, att, 1000.0);
104         final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
105 
106         // Force model parameters
107         final double[] parameters = srp.getParameters(orbit.getDate());
108         // Initialize force model
109         srp.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, parameters);
110 
111         // Register the attitude provider to the force model
112         AttitudeProvider attitudeProvider = new FrameAlignedProvider(rotation);
113         srp.registerAttitudeProvider(attitudeProvider );
114 
115         // Compute the mean element rate
116         final double[] elements = new double[7];
117         Arrays.fill(elements, 0.0);
118         final double[] daidt = srp.getMeanElementRate(state, auxiliaryElements, parameters);
119         for (int i = 0; i < daidt.length; i++) {
120             elements[i] = daidt[i];
121         }
122 
123         Assertions.assertEquals( 6.840790448823038E-8,    elements[0], 1.e-23);
124         Assertions.assertEquals(-2.990943627915497E-11,   elements[1], 1.e-26);
125         Assertions.assertEquals(-2.538400074176317E-10,   elements[2], 1.e-25);
126         Assertions.assertEquals( 2.037839945151859E-13,   elements[3], 1.e-28);
127         Assertions.assertEquals(-2.3338909771295392E-14,  elements[4], 1.e-29);
128         Assertions.assertEquals( 1.6082478750869883E-11,  elements[5], 1.e-26);
129 
130     }
131 
132     @Test
133     void testShortPeriodTerms() throws IllegalArgumentException {
134 
135         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
136 
137         final Orbit orbit = new EquinoctialOrbit(7069219.9806427825,
138                                                  -4.5941811292223825E-4,
139                                                  1.309932339472599E-4,
140                                                  -1.002996107003202,
141                                                  0.570979900577994,
142                                                  2.62038786211518,
143                                                  PositionAngleType.TRUE,
144                                                  FramesFactory.getEME2000(),
145                                                  initDate,
146                                                  3.986004415E14);
147 
148         final SpacecraftState meanState = new SpacecraftState(orbit);
149 
150         final CelestialBody    sun   = CelestialBodyFactory.getSun();
151 
152         final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0,
153                                                                                      sun,
154                                                                                      50.0, Vector3D.PLUS_J,
155                                                                                      2.0, 0.1,
156                                                                                      0.2, 0.6);
157 
158         final AttitudeProvider attitudeProvider = new LofOffset(meanState.getFrame(),
159                                                                 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
160                                                                 0.0, 0.0, 0.0);
161 
162         final DSSTForceModel srp = new DSSTSolarRadiationPressure(sun,
163                                                                   new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
164                                                                                        Constants.WGS84_EARTH_FLATTENING,
165                                                                                        FramesFactory.getITRF(IERSConventions.IERS_2010, false)),
166                                                                   boxAndWing,
167                                                                   meanState.getOrbit().getMu());
168 
169         //Create the auxiliary object
170         final AuxiliaryElements aux = new AuxiliaryElements(meanState.getOrbit(), 1);
171 
172         // Set the force models
173         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
174 
175         srp.registerAttitudeProvider(attitudeProvider);
176         shortPeriodTerms.addAll(srp.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, srp.getParameters(meanState.getDate())));
177         srp.updateShortPeriodTerms(srp.getParametersAllValues(), meanState);
178 
179         double[] y = new double[6];
180         for (final ShortPeriodTerms spt : shortPeriodTerms) {
181             final double[] shortPeriodic = spt.value(meanState.getOrbit());
182             for (int i = 0; i < shortPeriodic.length; i++) {
183                 y[i] += shortPeriodic[i];
184             }
185         }
186         Assertions.assertEquals( 0.3668654523023707,   y[0], 1.e-15);
187         Assertions.assertEquals(-2.5673332283107E-10,  y[1], 1.e-23);
188         Assertions.assertEquals(-3.84959877691969E-9,  y[2], 1.e-23);
189         Assertions.assertEquals(-3.069285299519558E-9, y[3], 1.e-24);
190         Assertions.assertEquals(-4.908870542277221E-9, y[4], 1.e-24);
191         Assertions.assertEquals(-2.38549338428359E-9,  y[5], 1.e-23);
192     }
193 
194     @BeforeEach
195     public void setUp() throws IOException, ParseException {
196         Utils.setDataRoot("regular-data:potential/shm-format");
197     }
198 }