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