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