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.models.earth.atmosphere;
18  
19  
20  import java.util.TimeZone;
21  
22  import org.hipparchus.geometry.euclidean.threed.Rotation;
23  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.util.Decimal64;
26  import org.hipparchus.util.FastMath;
27  import org.junit.Assert;
28  import org.junit.Before;
29  import org.junit.Test;
30  import org.orekit.SolarInputs97to05;
31  import org.orekit.Utils;
32  import org.orekit.bodies.CelestialBodyFactory;
33  import org.orekit.bodies.OneAxisEllipsoid;
34  import org.orekit.frames.Frame;
35  import org.orekit.frames.FramesFactory;
36  import org.orekit.frames.Transform;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.time.TimeScale;
39  import org.orekit.time.TimeScalesFactory;
40  import org.orekit.utils.Constants;
41  import org.orekit.utils.IERSConventions;
42  import org.orekit.utils.PVCoordinatesProvider;
43  
44  public class DTM2000Test {
45  
46      @Test
47      public void testWithOriginalTestsCases() {
48  
49          Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
50          PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
51          OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf);
52          SolarInputs97to05 in = SolarInputs97to05.getInstance();
53          earth.setAngularThreshold(1e-10);
54          DTM2000 atm = new DTM2000(in, sun, earth);
55          double roTestCase;
56          double myRo;
57  
58          // Inputs :
59  //      alt=800.
60  //      lat=40.
61  //      day=185.
62  //      hl=16.
63  //      xlon=0.
64  //      fm(1)=150.
65  //      f(1) =fm(1)
66  //      fm(2)=0.
67  //      f(2)=0.
68  //      akp(1)=0.
69  //      akp(2)=0.
70  //      akp(3)=0.
71  //      akp(4)=0.
72  
73          // Outputs :
74          roTestCase = 1.8710001353820e-17 * 1000;
75  
76          // Computation and results
77          myRo = atm.getDensity(185, 800*1000, 0, FastMath.toRadians(40), 16*FastMath.PI/12, 150, 150, 0, 0);
78          Assert.assertEquals(roTestCase, myRo , roTestCase * 1e-14);
79  
80  //      IDEM., day=275
81  
82          roTestCase=    2.8524195214905e-17* 1000;
83  
84          myRo = atm.getDensity(275, 800*1000, 0, FastMath.toRadians(40), 16*FastMath.PI/12, 150, 150, 0, 0);
85          Assert.assertEquals(roTestCase, myRo , roTestCase * 1e-14);
86  
87  //      IDEM., day=355
88  
89          roTestCase=    1.7343324462212e-17* 1000;
90  
91          myRo = atm.getDensity(355, 800*1000, 0, FastMath.toRadians(40), 16*FastMath.PI/12, 150, 150, 0, 0);
92          Assert.assertEquals(roTestCase, myRo , roTestCase * 2e-14);
93  //      IDEM., day=85
94  
95          roTestCase=    2.9983740796297e-17* 1000;
96  
97          myRo = atm.getDensity(85, 800*1000, 0, FastMath.toRadians(40), 16*FastMath.PI/12, 150, 150, 0, 0);
98          Assert.assertEquals(roTestCase, myRo , roTestCase * 1e-14);
99  
100 
101 //      alt=500.
102 //      lat=-70.      NB: the subroutine requires latitude in rad
103 //      day=15.
104 //      hl=16.        NB: the subroutine requires local time in rad (0hr=0 rad)
105 //      xlon=0.
106 //      fm(1)=70.
107 //      f(1) =fm(1)
108 //      fm(2)=0.
109 //      f(2)=0.
110 //      akp(1)=0.
111 //      akp(2)=0.
112 //      akp(3)=0.
113 //      akp(4)=0.
114 //      ro=    1.3150282384722D-16
115 //      tz=    793.65487014559
116 //      tinf=    793.65549802348
117         // note that the values above are the ones present in the original fortran source comments
118         // however, running this original source (converted to double precision) does
119         // not yield the results in the comments, but instead gives the following results
120         // as all the other tests cases do behave correctly, we assume the comments are wrong
121         // and prefer to ensure we get the same result as the original CODE
122         // as we don't have access to any other tests cases, we can't really decide if this is
123         // the best approach. Indeed, we are able to get the same results as original fortran
124         roTestCase =    1.5699108952425600E-016* 1000;
125 
126         myRo = atm.getDensity(15, 500*1000, 0, FastMath.toRadians(-70), 16*FastMath.PI/12, 70, 70, 0, 0);
127         Assert.assertEquals(roTestCase, myRo , roTestCase * 1e-14);
128 
129 //      IDEM., alt=800.
130 //      ro=    1.9556768571305D-18
131 //      tz=    793.65549797919
132 //      tinf=    793.65549802348
133         // note that the values above are the ones present in the original fortran source comments
134         // however, running this original source (converted to double precision) does
135         // not yield the results in the comments, but instead gives the following results
136         // as all the other tests cases do behave correctly, we assume the comments are wrong
137         // and prefer to ensure we get the same result as the original CODE
138         // as we don't have access to any other tests cases, we can't really decide if this is
139         // the best approach. Indeed, we are able to get the same results as original fortran
140         roTestCase =    2.4123751406975562E-018* 1000;
141         myRo = atm.getDensity(15, 800*1000, 0, FastMath.toRadians(-70), 16*FastMath.PI/12, 70, 70, 0, 0);
142         Assert.assertEquals(roTestCase, myRo , roTestCase * 1e-14);
143 
144     }
145 
146     @Test
147     public void testNonEarthRotationAxisAlignedFrame() {
148         //setup
149         AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
150         Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
151         Rotation rotation = new Rotation(Vector3D.PLUS_I, FastMath.PI / 2, RotationConvention.VECTOR_OPERATOR);
152         Frame frame = new Frame(ecef, new Transform(date, rotation), "other");
153         Vector3D pEcef = new Vector3D(6378137 + 300e3, 0, 0);
154         Vector3D pFrame = ecef.getTransformTo(frame, date)
155                 .transformPosition(pEcef);
156         PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
157         OneAxisEllipsoid earth = new OneAxisEllipsoid(
158                 6378136.460, 1.0 / 298.257222101, ecef);
159         SolarInputs97to05 in = SolarInputs97to05.getInstance();
160         earth.setAngularThreshold(1e-10);
161         DTM2000 atm = new DTM2000(in, sun, earth);
162 
163         //action
164         final double actual = atm.getDensity(date, pFrame, frame);
165 
166         //verify
167         Assert.assertEquals(atm.getDensity(date, pEcef, ecef), actual, 0.0);
168     }
169 
170     @Test
171     public void testField() {
172         Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
173         PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
174         OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf);
175         SolarInputs97to05 in = SolarInputs97to05.getInstance();
176         earth.setAngularThreshold(1e-10);
177         DTM2000 atm = new DTM2000(in, sun, earth);
178 
179         // Computation and results
180         for (double alti = 400; alti < 1000; alti += 50) {
181             for (double lon = 0; lon < 6; lon += 0.5) {
182                 for (double lat = -1.5; lat < 1.5; lat += 0.5) {
183                     for (double hl = 0; hl < 6; hl += 0.5) {
184                         double rhoD = atm.getDensity(185, alti*1000, lon, lat, hl, 50, 150, 0, 0);
185                         Decimal64 rho64 = atm.getDensity(185, new Decimal64(alti*1000),
186                                                          new Decimal64(lon), new Decimal64(lat),
187                                                          new Decimal64(hl), 50, 150, 0, 0);
188                         Assert.assertEquals(rhoD, rho64.getReal(), rhoD * 1e-14);
189                     }
190                 }
191             }
192         }
193 
194     }
195     
196     /** Test issue 539. Density computation should be independent of user's default time zone.
197      * See <a href="https://gitlab.orekit.org/orekit/orekit/issues/539"> issue 539 on Orekit forge.</a>
198      */
199     @Test
200     public void testTimeZoneIndependantIssue539() {
201         
202         // Prepare input: Choose a date in summer time for "GMT+1" time zone.
203         // So that after 22h in GMT we are in the next day in local time
204         TimeScale utc     = TimeScalesFactory.getUTC();
205         AbsoluteDate date = new AbsoluteDate("2000-04-01T22:30:00.000", utc);
206         
207         // LEO random position, in GCRF
208         Vector3D position = new Vector3D(-1038893.194, -4654348.144, 5021579.14);
209         Frame gcrf = FramesFactory.getGCRF();
210         
211         // Get ITRF and Earth ellipsoid
212         Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
213         PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
214         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.GRIM5C1_EARTH_EQUATORIAL_RADIUS,
215                                                       Constants.GRIM5C1_EARTH_FLATTENING, itrf);
216         SolarInputs97to05 in = SolarInputs97to05.getInstance();
217         earth.setAngularThreshold(1e-10);
218         DTM2000 atm = new DTM2000(in, sun, earth);
219         
220         // Store user time zone
221         TimeZone defaultTZ = TimeZone.getDefault();
222         
223         // Set default time zone to UTC & get density
224         TimeZone.setDefault(TimeZone.getTimeZone("Etc/UTC"));        
225         double rhoUtc = atm.getDensity(date, position, gcrf);
226         
227         // Set default time zone to GMT+1 & get density
228         TimeZone.setDefault(TimeZone.getTimeZone("Europe/Paris"));
229         double rhoParis = atm.getDensity(date, position, gcrf);
230         
231         // Check that the 2 densities are equal
232         Assert.assertEquals(0., rhoUtc - rhoParis, 0.);
233         
234         // Set back default time zone to what it was before the test, to avoid any interference with another routine
235         TimeZone.setDefault(defaultTZ);
236     }
237 
238     @Before
239     public void setUp() {
240         Utils.setDataRoot("regular-data");
241     }
242 
243 }