1   /* Copyright 2002-2016 CS Systèmes d'Information
2    * Licensed to CS Systèmes d'Information (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.forces.drag;
18  
19  
20  import java.io.FileNotFoundException;
21  import java.text.ParseException;
22  
23  import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
24  import org.apache.commons.math3.util.FastMath;
25  import org.junit.Assert;
26  import org.junit.Before;
27  import org.junit.Test;
28  import org.orekit.SolarInputs97to05;
29  import org.orekit.Utils;
30  import org.orekit.bodies.CelestialBodyFactory;
31  import org.orekit.bodies.GeodeticPoint;
32  import org.orekit.bodies.OneAxisEllipsoid;
33  import org.orekit.errors.OrekitException;
34  import org.orekit.errors.OrekitMessages;
35  import org.orekit.frames.Frame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.time.DateComponents;
39  import org.orekit.time.TimeComponents;
40  import org.orekit.time.TimeScalesFactory;
41  import org.orekit.utils.Constants;
42  import org.orekit.utils.IERSConventions;
43  import org.orekit.utils.PVCoordinatesProvider;
44  
45  public class JB2006Test {
46  
47      @Test
48      public void testWithOriginalTestsCases() throws OrekitException, ParseException {
49  
50          Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
51          PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
52          OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf);
53  
54          SolarInputs97to05 in = SolarInputs97to05.getInstance();
55          earth.setAngularThreshold(1e-10);
56          JB2006 atm = new JB2006(in, sun, earth);
57          double myRo;
58          double PI = 3.1415927;
59  
60          // SET SOLAR INDICES USE 1 DAY LAG FOR EUV AND F10 INFLUENCE
61          double S10  = 140;
62          double S10B = 100;
63          double F10  = 135;
64          double F10B = 95;
65          // USE 5 DAY LAG FOR MG FUV INFLUENCE
66          double XM10  = 130;
67          double XM10B = 95;
68  
69          // USE 6.7 HR LAG FOR ap INFLUENCE
70          double AP = 30;
71          // SET TIME OF INTEREST
72          double IYR  = 01;
73          double IDAY = 200;
74          if (IYR<50) IYR = IYR + 100;
75          double IYY  = (IYR-50)*365 + ((IYR-1)/4-12);
76          double ID1950 = IYY + IDAY;
77          double D1950  = ID1950;
78          double AMJD   = D1950 + 33281;
79  
80          // COMPUTE DENSITY KG/M3 RHO
81  
82          // alt = 400
83          myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*400,F10, F10B, AP,S10,S10B,XM10,XM10B);
84          double roTestCase = 0.4066e-11;
85          double tzTestCase=1137.7;
86          double tinfTestCase=1145.8;
87          Assert.assertEquals(roTestCase*1e12, FastMath.round(myRo*1e15)/1e3,0);
88          Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
89          Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
90  
91          // alt = 89.999km
92          try {
93              atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,89999.0,F10, F10B, AP,S10,S10B,XM10,XM10B);
94              Assert.fail("an exception should have been thrown");
95          } catch (OrekitException oe) {
96              Assert.assertEquals(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, oe.getSpecifier());
97              Assert.assertEquals(89999.0, (Double) oe.getParts()[0], 1.0e-15);
98              Assert.assertEquals(90000.0, (Double) oe.getParts()[1], 1.0e-15);
99          }
100 
101         // alt = 90
102         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*90,F10, F10B, AP,S10,S10B,XM10,XM10B);
103         roTestCase = 0.3285e-05;
104         tzTestCase=183.0;
105         tinfTestCase=1142.8;
106         Assert.assertEquals(roTestCase*1e05, FastMath.round(myRo*1e09)/1e4,0);
107         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
108         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
109 
110         // alt = 110
111         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*110,F10, F10B, AP,S10,S10B,XM10,XM10B);
112         roTestCase = 0.7587e-07;
113         tzTestCase=257.4;
114         tinfTestCase=1142.8;
115         Assert.assertEquals(roTestCase*1e07, FastMath.round(myRo*1e11)/1e4,0);
116         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
117         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
118 
119         // alt = 180
120         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*180,F10, F10B, AP,S10,S10B,XM10,XM10B);
121         roTestCase = 0.5439; // *1e-9
122         tzTestCase=915.0;
123         tinfTestCase=1130.9;
124         Assert.assertEquals(roTestCase, FastMath.round(myRo*1e13)/1e4,0);
125         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
126         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
127 
128         // alt = 230
129         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*230,F10, F10B, AP,S10,S10B,XM10,XM10B);
130         roTestCase =0.1250e-09;
131         tzTestCase=1047.5;
132         tinfTestCase=1137.4;
133         Assert.assertEquals(roTestCase*1e09, FastMath.round(myRo*1e13)/1e4,0);
134         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
135         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
136 
137         // alt = 270
138         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*270,F10, F10B, AP,S10,S10B,XM10,XM10B);
139         roTestCase =0.4818e-10;
140         tzTestCase=1095.6;
141         tinfTestCase=1142.5;
142         Assert.assertEquals(roTestCase*1e10, FastMath.round(myRo*1e14)/1e4,0);
143         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
144         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
145 
146         // alt = 660
147         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*660,F10, F10B, AP,S10,S10B,XM10,XM10B);
148         roTestCase =0.9451e-13;
149         tzTestCase=1149.0;
150         tinfTestCase=1149.9 ;
151         Assert.assertEquals(roTestCase*1e13, FastMath.round(myRo*1e17)/1e4,0);
152         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
153         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
154 
155         //  alt = 890
156         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*890,F10, F10B, AP,S10,S10B,XM10,XM10B);
157         roTestCase =0.8305e-14;
158         tzTestCase=1142.5;
159         tinfTestCase=1142.8 ;
160         Assert.assertEquals(roTestCase*1e14, FastMath.round(myRo*1e18)/1e4,0);
161         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
162         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
163 
164         //  alt = 1320
165         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*1320,F10, F10B, AP,S10,S10B,XM10,XM10B);
166         roTestCase =0.2004e-14;
167         tzTestCase=1142.7;
168         tinfTestCase=1142.8 ;
169         Assert.assertEquals(roTestCase*1e14, FastMath.round(myRo*1e18)/1e4,0);
170         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
171         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
172 
173         //  alt = 1600
174         myRo = atm.getDensity(AMJD, 90.*PI/180., 20.*PI/180.,90.*PI/180.,45.*PI/180.,1000*1600,F10, F10B, AP,S10,S10B,XM10,XM10B);
175         roTestCase = 0.1159e-14;
176         tzTestCase=1142.8;
177         tinfTestCase=1142.8 ;
178         Assert.assertEquals(roTestCase*1e14, FastMath.round(myRo*1e18)/1e4,0);
179         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
180         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
181 
182 
183         //  OTHER entries
184         AMJD +=50;
185         myRo = atm.getDensity(AMJD, 45.*PI/180., 10.*PI/180.,45.*PI/180.,-10.*PI/180.,400*1000,F10, F10B, AP,S10,S10B,XM10,XM10B);
186         roTestCase = 0.4838e-11;
187         tzTestCase=1137.4 ;
188         tinfTestCase= 1145.4 ;
189         Assert.assertEquals(roTestCase*1e11, FastMath.round(myRo*1e15)/1e4,0);
190         Assert.assertEquals(tzTestCase, FastMath.round(atm.getLocalTemp()*10)/10.0,0);
191         Assert.assertEquals(tinfTestCase, FastMath.round(atm.getExosphericTemp()*10)/10.0,0);
192 
193     }
194 
195     public void testComparisonWithDTM2000() throws OrekitException, ParseException, FileNotFoundException {
196 
197         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 01, 01),
198                                              TimeComponents.H00,
199                                              TimeScalesFactory.getUTC());
200         Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
201         PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
202         OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf);
203 
204         SolarInputs97to05 in = SolarInputs97to05.getInstance();
205         earth.setAngularThreshold(1e-10);
206         JB2006 jb = new JB2006(in, sun, earth);
207         DTM2000 dtm = new DTM2000(in, sun, earth);
208         // Positions
209 
210         Vector3D pos = new Vector3D(6500000.0,
211                                     -1234567.0,
212                                     4000000.0);
213 
214         // COMPUTE DENSITY KG/M3 RHO
215 
216         // alt = 400
217         double roJb = jb.getDensity(date, pos, FramesFactory.getEME2000());
218         double roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000());
219 
220         pos = new Vector3D(3011109.360780633,
221                            -5889822.669411588,
222                            4002170.0385907636);
223 
224         // COMPUTE DENSITY KG/M3 RHO
225 
226         // alt = 400
227         roJb = jb.getDensity(date, pos, FramesFactory.getEME2000());
228         roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000());
229 
230         pos =new Vector3D(-1033.4793830*1000,
231                           7901.2952754*1000,
232                           6380.3565958 *1000);
233 
234         // COMPUTE DENSITY KG/M3 RHO
235 
236         // alt = 400
237         roJb = jb.getDensity(date, pos, FramesFactory.getEME2000());
238         roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000());
239 
240         GeodeticPoint point;
241         for (int i = 0; i<367; i++) {
242             date = date.shiftedBy(Constants.JULIAN_DAY);
243             point = new GeodeticPoint(FastMath.toRadians(40), 0, 300*1000);
244             pos = earth.transform(point);
245             roJb = jb.getDensity(date, pos, FramesFactory.getEME2000());
246             roDtm = dtm.getDensity(date, pos, FramesFactory.getEME2000());
247             Assert.assertEquals(roDtm, roJb, roJb);
248 
249         }
250 
251     }
252 
253     public void testSolarInputs() throws OrekitException, ParseException {
254 
255         AbsoluteDate date = new AbsoluteDate(new DateComponents(2001, 01, 14),
256                                              TimeComponents.H00,
257                                              TimeScalesFactory.getUTC());
258 
259         SolarInputs97to05 in = SolarInputs97to05.getInstance();
260 
261 //      2001  14   2451924.0 176.3 164.4 180.0 180.4 163.4 169.2
262 //      14 176 164   9  12   9   6   4   4   9   7
263         Assert.assertEquals(176.3, in.getF10(date),0);
264         Assert.assertEquals(164.4, in.getF10B(date),0);
265         Assert.assertEquals(180.0, in.getS10(date),0);
266         Assert.assertEquals(180.4, in.getS10B(date),0);
267         Assert.assertEquals(163.4, in.getXM10(date),0);
268         Assert.assertEquals(169.2, in.getXM10B(date),0);
269         Assert.assertEquals(9 , in.getAp(date),0);
270 
271 
272         date = date.shiftedBy(11 * 3600);
273         Assert.assertEquals(6 , in.getAp(date),0);
274 
275         date = new AbsoluteDate(new DateComponents(1998, 02, 02),
276                                 new TimeComponents(18, 00, 00),
277                                 TimeScalesFactory.getUTC());
278 //      1998  33   2450847.0  89.1  95.1  95.8  97.9  81.3  92.0  1
279 //      33  89  95   4   5   4   4   2   0   0   3                          98
280         Assert.assertEquals(89.1, in.getF10(date),0);
281         Assert.assertEquals(95.1, in.getF10B(date),0);
282         Assert.assertEquals(95.8, in.getS10(date),0);
283         Assert.assertEquals(97.9, in.getS10B(date),0);
284         Assert.assertEquals(81.3, in.getXM10(date),0);
285         Assert.assertEquals(92.0, in.getXM10B(date),0);
286         Assert.assertEquals(0 , in.getAp(date),0);
287         date = date.shiftedBy(6 * 3600 - 1);
288         Assert.assertEquals(3 , in.getAp(date),0);
289     }
290 
291     @Before
292     public void setUp() {
293         Utils.setDataRoot("regular-data");
294     }
295 
296 }