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.models.earth.atmosphere;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.util.Binary64;
22  import org.hipparchus.util.Binary64Field;
23  import org.hipparchus.util.FastMath;
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.bodies.CelestialBody;
29  import org.orekit.bodies.CelestialBodyFactory;
30  import org.orekit.bodies.OneAxisEllipsoid;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.errors.OrekitMessages;
33  import org.orekit.frames.Frame;
34  import org.orekit.frames.FramesFactory;
35  import org.orekit.time.AbsoluteDate;
36  import org.orekit.time.FieldAbsoluteDate;
37  import org.orekit.time.TimeScalesFactory;
38  import org.orekit.utils.IERSConventions;
39  
40  class JB2006Test {
41  
42      private static double EPSILON = 1.0e-12;
43  
44      @BeforeEach
45      public void setUp() {
46          Utils.setDataRoot("regular-data:atmosphere");
47      }
48  
49      @Test
50      public void testWithOriginalTestsCases()  {
51  
52          Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
53          CelestialBody sun = CelestialBodyFactory.getSun();
54          OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf);
55          earth.setAngularThreshold(1e-10);
56  
57          // Static inputs
58          final AbsoluteDate date = new AbsoluteDate("2001-07-19T00:00:00.000Z", TimeScalesFactory.getUTC());
59          final double sunRa = FastMath.toRadians(90.);
60          final double sunDecli = FastMath.toRadians(20.);
61          final double subLon = FastMath.toRadians(90.);
62          final double subLat = FastMath.toRadians(45.);
63  
64          InputParameters inputParameters = new InputParameters();
65          JB2006 atm = new JB2006(inputParameters, sun, earth);
66  
67          // alt = 400
68          double computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 400000);
69          double referenceDensity = 4.066e-12;
70          Assertions.assertEquals(referenceDensity * 1e12, FastMath.round(computedDensity * 1e15) / 1e3, EPSILON);
71  
72          // alt = 89.999km
73          try {
74              atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 89999.0);
75              Assertions.fail("an exception should have been thrown");
76          } catch (OrekitException oe) {
77              Assertions.assertEquals(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, oe.getSpecifier());
78              Assertions.assertEquals(89999.0, (Double) oe.getParts()[0], 1.0e-15);
79              Assertions.assertEquals(90000.0, (Double) oe.getParts()[1], 1.0e-15);
80          }
81  
82          // alt = 90
83          computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat,90000);
84          referenceDensity = 0.3285e-05;
85          Assertions.assertEquals(referenceDensity * 1e05, FastMath.round(computedDensity * 1e09) / 1e4, EPSILON);
86  
87          // alt = 110
88          computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 110000);
89          referenceDensity = 0.7587e-07;
90          Assertions.assertEquals(referenceDensity * 1e07, FastMath.round(computedDensity * 1e11) / 1e4, EPSILON);
91  
92          // alt = 180
93          computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 180000);
94          referenceDensity = 0.5439;
95          Assertions.assertEquals(referenceDensity, FastMath.round(computedDensity * 1e13) / 1e4, EPSILON);
96  
97          // alt = 230
98          computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 230000);
99          referenceDensity = 0.1250e-09;
100         Assertions.assertEquals(referenceDensity * 1e09, FastMath.round(computedDensity * 1e13) / 1e4, EPSILON);
101 
102         // alt = 270
103         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 270000);
104         referenceDensity = 0.4818e-10;
105         Assertions.assertEquals(referenceDensity * 1e10, FastMath.round(computedDensity * 1e14) / 1e4, EPSILON);
106 
107         // alt = 660
108         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 660000);
109         referenceDensity = 0.9451e-13;
110         Assertions.assertEquals(referenceDensity * 1e13, FastMath.round(computedDensity*1e17) / 1e4, EPSILON);
111 
112         //  alt = 890
113         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 890000);
114         referenceDensity = 0.8305e-14;
115         Assertions.assertEquals(referenceDensity * 1e14, FastMath.round(computedDensity * 1e18) / 1e4, EPSILON);
116 
117         //  alt = 1320
118         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 1320000);
119         referenceDensity = 0.2004e-14;
120         Assertions.assertEquals(referenceDensity * 1e14, FastMath.round(computedDensity * 1e18) / 1e4, EPSILON);
121 
122         //  alt = 1600
123         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, 1600000);
124         referenceDensity = 0.1159e-14;
125         Assertions.assertEquals(referenceDensity * 1e14, FastMath.round(computedDensity * 1e18) / 1e4, EPSILON);
126     }
127 
128     @Test
129     public void testWithOriginalTestsCasesField() {
130         doTestWithOriginalTestsCasesField(Binary64Field.getInstance());
131     }
132 
133     private <T extends CalculusFieldElement<T>> void doTestWithOriginalTestsCasesField(Field<T> field)  {
134 
135         final T zero = field.getZero();
136         Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
137         CelestialBody sun = CelestialBodyFactory.getSun();
138         OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1.0 / 298.257222101, itrf);
139         earth.setAngularThreshold(1e-10);
140 
141         InputParameters inputParameters = new InputParameters();
142         JB2006 atm = new JB2006(inputParameters, sun, earth);
143 
144         // Static inputs
145         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2001-07-19T00:00:00.000Z", TimeScalesFactory.getUTC());
146         T sunRa = zero.add(FastMath.toRadians(90.));
147         T sunDecli = zero.add(FastMath.toRadians(20.));
148         T subLon = zero.add(FastMath.toRadians(90.));
149         T subLat = zero.add(FastMath.toRadians(45.));
150 
151         // alt = 400
152         T computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(400000));
153         double referenceDensity = 4.066e-12;
154         Assertions.assertEquals(referenceDensity * 1e12, FastMath.round(computedDensity.getReal() * 1e15) / 1e3, EPSILON);
155 
156         // alt = 89.999km
157         try {
158             atm.computeDensity(date, sunRa, sunDecli, subLon, subLat,zero.add(89999.0));
159             Assertions.fail("an exception should have been thrown");
160         } catch (OrekitException oe) {
161             Assertions.assertEquals(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, oe.getSpecifier());
162             Assertions.assertEquals(89999.0, (Double) oe.getParts()[0], 1.0e-15);
163             Assertions.assertEquals(90000.0, (Double) oe.getParts()[1], 1.0e-15);
164         }
165 
166         // alt = 90
167         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(90000));
168         referenceDensity = 0.3285e-05;
169         Assertions.assertEquals(referenceDensity * 1e05, FastMath.round(computedDensity.getReal() * 1e09) / 1e4, EPSILON);
170 
171         // alt = 110
172         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(110000));
173         referenceDensity = 0.7587e-07;
174         Assertions.assertEquals(referenceDensity * 1e07, FastMath.round(computedDensity.getReal() * 1e11) / 1e4, EPSILON);
175 
176         // alt = 180
177         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(180000));
178         referenceDensity = 0.5439;
179         Assertions.assertEquals(referenceDensity, FastMath.round(computedDensity.getReal() * 1e13) / 1e4, EPSILON);
180 
181         // alt = 230
182         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(230000));
183         referenceDensity = 0.1250e-09;
184         Assertions.assertEquals(referenceDensity * 1e09, FastMath.round(computedDensity.getReal() * 1e13) / 1e4, EPSILON);
185 
186         // alt = 270
187         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(270000));
188         referenceDensity = 0.4818e-10;
189         Assertions.assertEquals(referenceDensity * 1e10, FastMath.round(computedDensity.getReal() * 1e14) / 1e4, EPSILON);
190 
191         // alt = 660
192         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(660000));
193         referenceDensity = 0.9451e-13;
194         Assertions.assertEquals(referenceDensity * 1e13, FastMath.round(computedDensity.getReal() * 1e17) / 1e4, EPSILON);
195 
196         //  alt = 890
197         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat, zero.add(890000));
198         referenceDensity = 0.8305e-14;
199         Assertions.assertEquals(referenceDensity * 1e14, FastMath.round(computedDensity.getReal() * 1e18) / 1e4, EPSILON);
200 
201         //  alt = 1320
202         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat,zero.add(1320000));
203         referenceDensity = 0.2004e-14;
204         Assertions.assertEquals(referenceDensity*1e14, FastMath.round(computedDensity.getReal() * 1e18) / 1e4, EPSILON);
205 
206         //  alt = 1600
207         computedDensity = atm.computeDensity(date, sunRa, sunDecli, subLon, subLat,zero.add(1600000));
208         referenceDensity = 0.1159e-14;
209         Assertions.assertEquals(referenceDensity * 1e14, FastMath.round(computedDensity.getReal() * 1e18) / 1e4, EPSILON);
210     }
211 
212     private class InputParameters implements JB2006InputParameters {
213 
214         @Override
215         public double getF10(AbsoluteDate date) {
216             return 135;
217         }
218 
219         @Override
220         public double getF10B(AbsoluteDate date) {
221             return 95;
222         }
223 
224         @Override
225         public double getS10(AbsoluteDate date) {
226             return 140;
227         }
228 
229         @Override
230         public double getS10B(AbsoluteDate date) {
231             return 100;
232         }
233 
234         @Override
235         public double getXM10(AbsoluteDate date) {
236             return 130;
237         }
238 
239         @Override
240         public double getXM10B(AbsoluteDate date) {
241             return 95;
242         }
243 
244         @Override
245         public double getAp(AbsoluteDate date) {
246             return 30;
247         }
248 
249         @Override
250         public AbsoluteDate getMinDate() {
251             return AbsoluteDate.PAST_INFINITY;
252         }
253 
254         @Override
255         public AbsoluteDate getMaxDate() {
256             return AbsoluteDate.FUTURE_INFINITY;
257         }
258     }
259 
260 }