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.analysis.differentiation.DSFactory;
22  import org.hipparchus.analysis.differentiation.DerivativeStructure;
23  import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.util.Binary64;
27  import org.hipparchus.util.Binary64Field;
28  import org.hipparchus.util.FastMath;
29  import org.hipparchus.util.MathUtils;
30  import org.junit.jupiter.api.Assertions;
31  import org.junit.jupiter.api.BeforeEach;
32  import org.junit.jupiter.api.Test;
33  import org.orekit.Utils;
34  import org.orekit.bodies.CelestialBody;
35  import org.orekit.bodies.CelestialBodyFactory;
36  import org.orekit.bodies.GeodeticPoint;
37  import org.orekit.bodies.OneAxisEllipsoid;
38  import org.orekit.errors.OrekitException;
39  import org.orekit.errors.OrekitMessages;
40  import org.orekit.frames.Frame;
41  import org.orekit.frames.FramesFactory;
42  import org.orekit.models.earth.atmosphere.data.JB2008SpaceEnvironmentData;
43  import org.orekit.time.AbsoluteDate;
44  import org.orekit.time.DateComponents;
45  import org.orekit.time.FieldAbsoluteDate;
46  import org.orekit.time.TimeComponents;
47  import org.orekit.time.TimeScale;
48  import org.orekit.time.TimeScalesFactory;
49  import org.orekit.utils.Constants;
50  import org.orekit.utils.IERSConventions;
51  
52  import java.text.ParseException;
53  
54  class JB2008Test {
55  
56      private static final TimeScale TT = TimeScalesFactory.getTT();
57  
58      @Test
59      void testLegacy() {
60          final boolean print = false;
61          // Build the model
62          final CelestialBody sun = CelestialBodyFactory.getSun();
63          final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
64          final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
65                                                              Constants.WGS84_EARTH_FLATTENING, itrf);
66          final JB2008 atm = new JB2008(new LegacyInputParams(), sun, earth);
67  
68          // Loop over cases
69          for (int i = 0; i < 12; i++) {
70              final double rho = atm.computeDensity(LegacyInputParams.TC[i], LegacyInputParams.SUNRA[i], LegacyInputParams.SUNDEC[i],
71                                                    RAP(LegacyInputParams.TC[i], LegacyInputParams.SATLON[i]),
72                                                    FastMath.toRadians(LegacyInputParams.SATLAT[i]),
73                                                    LegacyInputParams.SATALT[i] * 1000.);
74              checkLegacy(i, rho, print);
75          }
76  
77      }
78  
79      @Test
80      void testDensityWithLocalSolarActivityData() {
81          // First case of "testAltitude"
82          final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
83          final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
84                                                              Constants.WGS84_EARTH_FLATTENING, itrf);
85          final CelestialBody sun = CelestialBodyFactory.getSun();
86          final JB2008 model = new JB2008(null, sun, earth);
87          double referenceDensity = 0.27945654e-05;
88          double computedDensity = model.getDensity(52951.003805740744, 3.046653643566772, -0.285987757544287,
89                                                    1.28211886851503, -1.4877186543999,
90                                                    91.0e+3, 91.00, 137.10,
91                                                    108.80, 123.80, 116.70, 128.50,
92                                                    168.00, 138.60, 43.);
93          Assertions.assertEquals(referenceDensity, computedDensity, referenceDensity * 2.e-5);
94      }
95  
96      @Test
97      void testDensityWithLocalSolarActivityDataField() {
98          doTestDensityWithLocalSolarActivityData(Binary64Field.getInstance());
99      }
100 
101     <T extends CalculusFieldElement<T>> void doTestDensityWithLocalSolarActivityData(Field<T> field) {
102         // First case of "testAltitude"
103         T zero = field.getZero();
104         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
105         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
106                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
107         final CelestialBody sun = CelestialBodyFactory.getSun();
108         final JB2008 model = new JB2008(null, sun, earth);
109         double referenceDensity = 0.27945654e-05;
110         double computedDensity = model.getDensity(zero.add(52951.003805740744), zero.add(3.046653643566772), zero.add(-0.285987757544287),
111                                                   zero.add(1.28211886851503), zero.add(-1.4877186543999),
112                                                   zero.add(91.0e+3), 91.00, 137.10,
113                                                   108.80, 123.80, 116.70, 128.50,
114                                                   168.00, 138.60, 43.).getReal();
115         Assertions.assertEquals(referenceDensity, computedDensity, referenceDensity * 2.e-5);
116     }
117 
118     @Test
119     void testAltitude() {
120         final boolean print = false;
121         // Build the iput params provider
122         final InputParams ip = new InputParams();
123         // Get Sun
124         final CelestialBody sun = CelestialBodyFactory.getSun();
125         // Get Earth body shape
126         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
127         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
128                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
129         earth.setAngularThreshold(1e-10);
130         // Build the model
131         final JB2008 atm = new JB2008(ip, sun, earth);
132 
133         // Reference locations as {lat, lon, alt}
134         final double[][] loc = {{-85.24,  73.46,   91.0e+3},
135                                 {-18.65, 218.68,  110.0e+3},
136                                 {-68.05, 145.28,  122.0e+3},
137                                 { 37.71,  34.55,  150.0e+3},
138                                 { 74.36,  46.25,  220.0e+3},
139                                 { -8.85, 216.56,  270.0e+3},
140                                 {-39.64,  32.00,  400.0e+3},
141                                 {-51.93,  38.83,  550.0e+3},
142                                 {-21.25, 213.67,  700.0e+3},
143                                 { 46.43,  29.37,  900.0e+3},
144                                 { 65.97,  38.12, 1200.0e+3},
145                                 {-21.31, 211.78, 1700.0e+3},
146                                 {-51.87,  26.64, 2300.0e+3}};
147 
148         // Loop over cases
149         for (int i = 0; i < 13; i++) {
150             // Build the point
151             final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(loc[i][0]),
152                                                           FastMath.toRadians(loc[i][1]),
153                                                           loc[i][2]);
154             // Run
155             final double rho = atm.getDensity(InputParams.TC[i], earth.transform(point), atm.getFrame());
156 
157             // Check
158             checkAltitude(i, rho, print);
159         }
160    }
161 
162     @Test
163     void testException() throws ParseException {
164 
165         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
166         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
167                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
168 
169         final JB2008 atm = new JB2008(new InputParams(), CelestialBodyFactory.getSun(), earth);
170 
171         // alt = 89.999km
172         try {
173             atm.computeDensity(AbsoluteDate.ARBITRARY_EPOCH, 0., 0., 0., 0., 89999.0);
174             Assertions.fail("an exception should have been thrown");
175         } catch (OrekitException oe) {
176             Assertions.assertEquals(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, oe.getSpecifier());
177             Assertions.assertEquals(89999.0, (Double) oe.getParts()[0], 1.0e-15);
178             Assertions.assertEquals(90000.0, (Double) oe.getParts()[1], 1.0e-15);
179         }
180 
181     }
182 
183     @Test
184     void testDensityField() {
185 
186         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
187         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
188                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
189 
190         final JB2008 atm = new JB2008(new InputParams(), CelestialBodyFactory.getSun(), earth);
191 
192         final AbsoluteDate date = InputParams.TC[4];
193 
194         for (double alt = 100; alt < 1000; alt += 50) {
195             for (double lat = -1.2; lat < 1.2; lat += 0.4) {
196                 for (double lon = 0; lon < 6.28; lon += 0.8) {
197 
198                     final GeodeticPoint point = new GeodeticPoint(lat, lon, alt * 1000.);
199                     final Vector3D pos = earth.transform(point);
200                     Field<Binary64> field = Binary64Field.getInstance();
201 
202                     // Run
203                     final double    rho = atm.getDensity(date, pos, itrf);
204                     final Binary64 rho64 = atm.getDensity(new FieldAbsoluteDate<>(field, date),
205                                                            new FieldVector3D<>(field, pos),
206                                                            itrf);
207 
208                     Assertions.assertEquals(rho, rho64.getReal(), rho * 4.0e-10);
209 
210                 }
211             }
212         }
213 
214     }
215 
216     @Test
217     void testDensityGradient() {
218 
219         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
220         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
221                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
222 
223         final JB2008 atm = new JB2008(new InputParams(), CelestialBodyFactory.getSun(), earth);
224 
225         final AbsoluteDate date = InputParams.TC[6];
226 
227         // Build the position
228         final double alt = 400.;
229         final double lat =  60.;
230         final double lon = -70.;
231         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(lat),
232                                                       FastMath.toRadians(lon),
233                                                       alt * 1000.);
234         final Vector3D pos = earth.transform(point);
235 
236         // Run
237         DerivativeStructure zero = new DSFactory(1, 1).variable(0, 0.0);
238         FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 10.0);
239         DerivativeStructure  rhoX = differentiator.
240                         differentiate((double x) -> {
241                             try {
242                                 return atm.getDensity(date, new Vector3D(1, pos, x, Vector3D.PLUS_I), itrf);
243                             } catch (OrekitException oe) {
244                                 return Double.NaN;
245                             }
246                         }). value(zero);
247         DerivativeStructure  rhoY = differentiator.
248                         differentiate((double y) -> {
249                             try {
250                                 return atm.getDensity(date, new Vector3D(1, pos, y, Vector3D.PLUS_J), itrf);
251                             } catch (OrekitException oe) {
252                                 return Double.NaN;
253                             }
254                         }). value(zero);
255         DerivativeStructure  rhoZ = differentiator.
256                         differentiate((double z) -> {
257                             try {
258                                 return atm.getDensity(date, new Vector3D(1, pos, z, Vector3D.PLUS_K), itrf);
259                             } catch (OrekitException oe) {
260                                 return Double.NaN;
261                             }
262                         }). value(zero);
263 
264         DSFactory factory3 = new DSFactory(3, 1);
265         Field<DerivativeStructure> field = factory3.getDerivativeField();
266         final DerivativeStructure rhoDS = atm.getDensity(new FieldAbsoluteDate<>(field, date),
267                                                          new FieldVector3D<>(factory3.variable(0, pos.getX()),
268                                                                              factory3.variable(1, pos.getY()),
269                                                                              factory3.variable(2, pos.getZ())),
270                                                          itrf);
271 
272         Assertions.assertEquals(rhoX.getValue(), rhoDS.getReal(), rhoX.getValue() * 2.0e-14);
273         Assertions.assertEquals(rhoY.getValue(), rhoDS.getReal(), rhoY.getValue() * 2.0e-14);
274         Assertions.assertEquals(rhoZ.getValue(), rhoDS.getReal(), rhoZ.getValue() * 2.0e-14);
275         Assertions.assertEquals(rhoX.getPartialDerivative(1),
276                             rhoDS.getPartialDerivative(1, 0, 0),
277                             FastMath.abs(6.0e-10 * rhoX.getPartialDerivative(1)));
278         Assertions.assertEquals(rhoY.getPartialDerivative(1),
279                             rhoDS.getPartialDerivative(0, 1, 0),
280                             FastMath.abs(6.0e-10 * rhoY.getPartialDerivative(1)));
281         Assertions.assertEquals(rhoZ.getPartialDerivative(1),
282                             rhoDS.getPartialDerivative(0, 0, 1),
283                             FastMath.abs(6.0e-10 * rhoY.getPartialDerivative(1)));
284 
285     }
286 
287     @Test
288     void testComparisonWithReference() {
289 
290         // The objective of this test is to compare Orekit results with
291         // the reference JB2008 Code Files provided by Space Environment
292         // (i.e., JB2008.for and JB08DRVY2K.for)
293 
294         // Earth
295         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
296         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
297                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
298 
299         // Input
300         final int    year  = 2004;
301         final int    month = 1;
302         final int    day   = 2;
303         final int    hour  = 12;
304         final int    min   = 0;
305         final double sec   = 0.0;
306         final double lat   = FastMath.toRadians(45.0);
307         final double lon   = FastMath.toRadians(45.0);
308         final double alt   = 250.0e3;
309 
310         // Initialize JB2008 model
311         final JB2008SpaceEnvironmentData JBData = new JB2008SpaceEnvironmentData("SOLFSMY_trunc.txt", "DTCFILE_trunc.TXT");
312         final JB2008 atm = new JB2008(JBData, CelestialBodyFactory.getSun(), earth);
313 
314         // Compute density
315         final GeodeticPoint point = new GeodeticPoint(lat, lon, alt);
316         final Vector3D pos = earth.transform(point);
317         final AbsoluteDate date = new AbsoluteDate(year, month, day, hour, min, sec, TimeScalesFactory.getUTC());
318         final double density = atm.getDensity(date, pos, itrf);
319 
320         // Verify
321         final double ref = 6.6862e-11;
322         Assertions.assertEquals(ref, density, 1.0e-15);
323 
324     }
325 
326     /** Convert longitude of position to right ascension of position.
327      * @param d1950 duration from fifties epoch
328      * @param satLon longitude of position (°)
329      * @return right ascension of position (rad)
330      */
331     private double RAP(final AbsoluteDate date, final double satLon) {
332         double d1950 = date.getMJD() - 33281.;
333         double theta;
334         final double nbday = FastMath.floor(d1950);
335         if (nbday < 7305.) {
336             theta = 1.7294446614 + 1.72027915246e-2 * nbday + 6.3003880926 * (d1950 - nbday);
337         } else {
338             final double ts70 = d1950 - 7305.;
339             final double ids70 = FastMath.floor(ts70);
340             final double tfrac = ts70 - ids70;
341             theta = 1.73213438565 + 1.720279169407e-2 * ids70 +
342                     (1.720279169407e-2 + MathUtils.TWO_PI) * tfrac +
343                     5.0755141943e-15 * ts70 * ts70;
344         }
345         theta = MathUtils.normalizeAngle(theta, FastMath.PI);
346 
347         return MathUtils.normalizeAngle(theta + FastMath.toRadians(satLon), 0.);
348     }
349 
350     /** Check legacy results.
351      * @param id legacy test number
352      * @param rho computed density
353      * @param print if true print else check
354      */
355     private void checkLegacy(final int id, final double rho, final boolean print) {
356         final double[] rhoRef  = {0.18730056e-11, 0.25650339e-11, 0.57428913e-11,
357                                   0.83266893e-11, 0.82238726e-11, 0.48686457e-11,
358                                   0.67210914e-11, 0.74215571e-11, 0.31821075e-11,
359                                   0.29553578e-11, 0.64122627e-11, 0.79559727e-11};
360         final double dRho = 2.e-5;
361         final int nb = id + 1;
362         if (print) {
363             System.out.printf("Case #%d\n", nb);
364             System.out.printf("Rho:  %12.5e  %12.5e\n", rhoRef[id], rho);
365         } else {
366             Assertions.assertEquals(rhoRef[id],  rho,  rhoRef[id]  * dRho);
367         }
368 
369     }
370 
371     /** Check altiude results.
372      * @param id test number
373      * @param rho computed density
374      * @param print if true print else check
375      */
376     private void checkAltitude(final int id, final double rho, final boolean print) {
377         final double[] rhoRef  = {0.27945654e-05, 0.94115202e-07, 0.15025977e-07, 0.21128330e-08,
378                                   0.15227435e-09, 0.54609767e-10, 0.45899746e-11, 0.14922800e-12,
379                                   0.17392987e-13, 0.35250121e-14, 0.13482414e-14, 0.77684879e-15,
380                                   0.19900569e-15};
381         final double dRho = 3.e-4;
382         final int nb = id + 1;
383         if (print) {
384           System.out.printf("Case #%d\n", nb);
385           System.out.printf("Rho:  %12.5e  %12.5e\n", rhoRef[id], rho);
386         } else {
387             Assertions.assertEquals(rhoRef[id],  rho,  rhoRef[id]  * dRho);
388         }
389 
390     }
391 
392     @BeforeEach
393     public void setUp() {
394         Utils.setDataRoot("regular-data:atmosphere");
395     }
396 
397     private static class LegacyInputParams implements JB2008InputParameters {
398 
399         private static final AbsoluteDate[] TC  = {AbsoluteDate.createMJDDate(MJD(20035), 0.00454861111 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
400                                                    AbsoluteDate.createMJDDate(MJD(20035), 0.50362 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
401                                                    AbsoluteDate.createMJDDate(MJD(20036), 0.00468 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
402                                                    AbsoluteDate.createMJDDate(MJD(20036), 0.50375 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
403                                                    AbsoluteDate.createMJDDate(MJD(20037), 0. * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
404                                                    AbsoluteDate.createMJDDate(MJD(20037), 0.50556 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
405                                                    AbsoluteDate.createMJDDate(MJD(20038), 0.00088 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
406                                                    AbsoluteDate.createMJDDate(MJD(20038), 0.50644 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
407                                                    AbsoluteDate.createMJDDate(MJD(20039), 0.00543 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
408                                                    AbsoluteDate.createMJDDate(MJD(20039), 0.50450* Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
409                                                    AbsoluteDate.createMJDDate(MJD(20040), 0. * Constants.JULIAN_DAY, TimeScalesFactory.getUTC()),
410                                                    AbsoluteDate.createMJDDate(MJD(20040), 0.50556 * Constants.JULIAN_DAY, TimeScalesFactory.getUTC())};
411 
412         private static final double[] SUNRA  = {3.8826, 3.8914, 3.9001, 3.9089, 3.9176, 3.9265,
413                                                 3.9352, 3.9441, 3.9529, 3.9618, 3.9706, 3.9795};
414 
415         private static final double[] SUNDEC = {-0.2847, -0.2873, -0.2898, -0.2923, -0.2948, -0.2973,
416                                                 -0.2998, -0.3022, -0.3046, -0.3070, -0.3094, -0.3117};
417 
418         private static final double[] SATLON = {73.46, 218.68, 34.55, 46.25, 216.56, 32.00,
419                                                 38.83, 213.67, 29.37, 38.12, 211.78, 26.64};
420 
421         private static final double[] SATLAT = {-85.24, -18.65, 37.71, 74.36,  -8.85, -39.64,
422                                                 -51.93, -21.25, 46.43, 65.97, -21.31, -51.87};
423 
424         private static final double[] SATALT = {398.91, 376.75, 373.45, 380.61, 374.03, 385.05,
425                                                 389.83, 376.98, 374.56, 378.97, 377.76, 390.09};
426 
427         private static final double[] F10    = {128.80, 128.80, 129.60, 129.60, 124.10, 124.10,
428                                                 140.90, 140.90, 104.60, 104.60,  94.90,  94.90};
429 
430         private static final double[] F10B   = {105.60, 105.60, 105.60, 105.60, 105.60, 105.60,
431                                                 105.60, 105.60, 105.70, 105.70, 105.90, 105.90};
432 
433         private static final double[] S10    = {103.50, 103.50, 110.20, 110.20, 109.40, 109.40,
434                                                 108.60, 108.60, 107.40, 107.40, 110.90, 110.90};
435 
436         private static final double[] S10B   = {103.80, 103.80, 103.80, 103.80, 103.80, 103.80,
437                                                 103.70, 103.70, 103.70, 103.70, 103.70, 103.70};
438 
439         private static final double[] XM10   = {110.90, 110.90, 115.60, 115.60, 110.00, 110.00,
440                                                 110.00, 110.00, 106.90, 106.90, 102.20, 102.20};
441 
442         private static final double[] XM10B  = {106.90, 106.90, 106.90, 106.90, 106.90, 106.90,
443                                                 107.00, 107.00, 107.10, 107.10, 107.10, 107.10};
444 
445         private static final double[] Y10    = {127.90, 127.90, 125.90, 125.90, 127.70, 127.70,
446                                                 125.60, 125.60, 126.60, 126.60, 126.80, 126.80};
447 
448         private static final double[] Y10B   = {112.90, 112.90, 112.90, 112.90, 113.00, 113.00,
449                                                 113.20, 113.20, 113.20, 113.20, 113.30, 113.30};
450 
451         private static final double[] DSTDTC = {  3.,  80., 240., 307., 132.,  40.,
452                                                 327., 327., 118.,  25.,  85., 251.};
453 
454         @Override
455         public AbsoluteDate getMinDate() {
456             return TC[0];
457         }
458 
459         @Override
460         public AbsoluteDate getMaxDate() {
461             return TC[TC.length - 1];
462         }
463 
464         @Override
465         public double getF10(AbsoluteDate date)
466         {
467             for (int i = 0; i < TC.length; i++) {
468                 if (date.equals(TC[i])) {
469                     return F10[i];
470                 }
471             }
472             return Double.NaN;
473         }
474 
475         @Override
476         public double getF10B(AbsoluteDate date)
477         {
478             for (int i = 0; i < TC.length; i++) {
479                 if (date.equals(TC[i])) {
480                     return F10B[i];
481                 }
482             }
483             return Double.NaN;
484         }
485 
486         @Override
487         public double getS10(AbsoluteDate date)
488         {
489             for (int i = 0; i < TC.length; i++) {
490                 if (date.equals(TC[i])) {
491                     return S10[i];
492                 }
493             }
494             return Double.NaN;
495         }
496 
497         @Override
498         public double getS10B(AbsoluteDate date)
499         {
500             for (int i = 0; i < TC.length; i++) {
501                 if (date.equals(TC[i])) {
502                     return S10B[i];
503                 }
504             }
505             return Double.NaN;
506         }
507 
508         @Override
509         public double getXM10(AbsoluteDate date)
510         {
511             for (int i = 0; i < TC.length; i++) {
512                 if (date.equals(TC[i])) {
513                     return XM10[i];
514                 }
515             }
516             return Double.NaN;
517         }
518 
519         @Override
520         public double getXM10B(AbsoluteDate date)
521         {
522             for (int i = 0; i < TC.length; i++) {
523                 if (date.equals(TC[i])) {
524                     return XM10B[i];
525                 }
526             }
527             return Double.NaN;
528         }
529 
530         @Override
531         public double getY10(AbsoluteDate date)
532         {
533             for (int i = 0; i < TC.length; i++) {
534                 if (date.equals(TC[i])) {
535                     return Y10[i];
536                 }
537             }
538             return Double.NaN;
539         }
540 
541         @Override
542         public double getY10B(AbsoluteDate date)
543         {
544             for (int i = 0; i < TC.length; i++) {
545                 if (date.equals(TC[i])) {
546                     return Y10B[i];
547                 }
548             }
549             return Double.NaN;
550         }
551 
552         @Override
553         public double getDSTDTC(AbsoluteDate date)
554         {
555             for (int i = 0; i < TC.length; i++) {
556                 if (date.equals(TC[i])) {
557                     return DSTDTC[i];
558                 }
559             }
560             return Double.NaN;
561         }
562 
563 
564         private static int MJD(final int d1950) {
565             return d1950 + 33281;
566         }
567 
568     }
569 
570     private static class InputParams implements JB2008InputParameters {
571 
572         /** Serializable UID. */
573         private static final long serialVersionUID = 5091441542522297257L;
574 
575         /** Dates. */
576         public static final AbsoluteDate[] TC = new AbsoluteDate[] {
577             new AbsoluteDate(new DateComponents(2003, 312), new TimeComponents( 0,  6, 33.0), TT),
578             new AbsoluteDate(new DateComponents(2003, 312), new TimeComponents(12,  5, 13.0), TT),
579             new AbsoluteDate(new DateComponents(2003, 312), new TimeComponents(18, 45,  3.0), TT),
580             new AbsoluteDate(new DateComponents(2003, 313), new TimeComponents( 0,  6, 44.0), TT),
581             new AbsoluteDate(new DateComponents(2003, 313), new TimeComponents(12,  5, 24.0), TT),
582             new AbsoluteDate(new DateComponents(2003, 314), new TimeComponents( 0,  0,  0.0), TT),
583             new AbsoluteDate(new DateComponents(2003, 314), new TimeComponents(12,  8,  0.0), TT),
584             new AbsoluteDate(new DateComponents(2003, 315), new TimeComponents( 0,  1, 16.0), TT),
585             new AbsoluteDate(new DateComponents(2003, 315), new TimeComponents(12,  9, 16.0), TT),
586             new AbsoluteDate(new DateComponents(2003, 316), new TimeComponents( 0,  7, 49.0), TT),
587             new AbsoluteDate(new DateComponents(2003, 316), new TimeComponents(12,  6, 29.0), TT),
588             new AbsoluteDate(new DateComponents(2003, 317), new TimeComponents( 0,  0,  0.0), TT),
589             new AbsoluteDate(new DateComponents(2003, 317), new TimeComponents(12,  8,  0.0), TT)
590         };
591 
592         /** F10 data. */
593         private static final double[] F10 = new double[] {
594             91.00, 91.00, 91.00, 92.70, 92.70, 93.00,
595             93.00, 94.60, 94.60, 95.60, 95.60, 98.70, 98.70
596         };
597 
598         /** F10B data. */
599         private static final double[] F10B = new double[] {
600             137.10, 137.10, 137.10, 136.90, 136.90, 136.80,
601             136.80, 136.70, 136.70, 136.70, 136.70, 136.90, 136.90
602         };
603 
604         /** S10 data. */
605         private static final double[] S10 = new double[] {
606             108.80, 108.80, 108.80, 104.20, 104.20, 102.60,
607             102.60, 100.30, 100.30,  99.50,  99.50, 101.20, 101.20
608         };
609 
610         /** S10B data. */
611         private static final double[] S10B = new double[] {
612             123.80, 123.80, 123.80, 123.70, 123.70, 123.60,
613             123.60, 123.50, 123.50, 123.50, 123.50, 123.60, 123.60
614         };
615 
616         /** XM10 data. */
617         private static final double[] XM10 = new double[] {
618             116.70, 116.70, 116.70, 109.60, 109.60, 100.20,
619             100.20,  97.00,  97.00,  95.40,  95.40,  95.40,  95.40
620         };
621 
622         /** XM10B data. */
623         private static final double[] XM10B = new double[] {
624             128.50, 128.50, 128.50, 128.00, 128.00, 127.70,
625             127.70, 127.60, 127.60, 127.60, 127.60, 127.70, 127.70
626         };
627 
628         /** Y10 data. */
629         private static final double[] Y10 = new double[] {
630             168.00, 168.00, 168.00, 147.90, 147.90, 131.60,
631             131.60, 122.60, 122.60, 114.30, 114.30, 112.70, 112.70
632         };
633 
634         /** Y10B data. */
635         private static final double[] Y10B = new double[] {
636             138.60, 138.60, 138.60, 138.40, 138.40, 138.10,
637             138.10, 137.90, 137.90, 137.90, 137.90, 137.80, 137.80
638         };
639 
640         /** DSTDTC data. */
641         private static final double[] DSTDTC = new double[] {
642              43.00,  30.00,  67.00,  90.00,  87.00, 115.00,
643             114.00, 106.00, 148.00, 148.00, 105.00, 146.00, 119.00
644         };
645 
646         /** Constructor. */
647         public InputParams() {
648         }
649 
650         @Override
651         public AbsoluteDate getMinDate() {
652             return TC[0];
653         }
654 
655         @Override
656         public AbsoluteDate getMaxDate() {
657             return TC[TC.length - 1];
658         }
659 
660         @Override
661         public double getF10(AbsoluteDate date)
662             {
663             for (int i = 0; i < TC.length; i++) {
664                 if (date.equals(TC[i])) {
665                     return F10[i];
666                 }
667             }
668             return Double.NaN;
669         }
670 
671         @Override
672         public double getF10B(AbsoluteDate date)
673             {
674             for (int i = 0; i < TC.length; i++) {
675                 if (date.equals(TC[i])) {
676                     return F10B[i];
677                 }
678             }
679             return Double.NaN;
680         }
681 
682         @Override
683         public double getS10(AbsoluteDate date)
684             {
685             for (int i = 0; i < TC.length; i++) {
686                 if (date.equals(TC[i])) {
687                     return S10[i];
688                 }
689             }
690             return Double.NaN;
691         }
692 
693         @Override
694         public double getS10B(AbsoluteDate date)
695             {
696             for (int i = 0; i < TC.length; i++) {
697                 if (date.equals(TC[i])) {
698                     return S10B[i];
699                 }
700             }
701             return Double.NaN;
702         }
703 
704         @Override
705         public double getXM10(AbsoluteDate date)
706             {
707             for (int i = 0; i < TC.length; i++) {
708                 if (date.equals(TC[i])) {
709                     return XM10[i];
710                 }
711             }
712             return Double.NaN;
713         }
714 
715         @Override
716         public double getXM10B(AbsoluteDate date)
717             {
718             for (int i = 0; i < TC.length; i++) {
719                 if (date.equals(TC[i])) {
720                     return XM10B[i];
721                 }
722             }
723             return Double.NaN;
724         }
725 
726         @Override
727         public double getY10(AbsoluteDate date)
728             {
729             for (int i = 0; i < TC.length; i++) {
730                 if (date.equals(TC[i])) {
731                     return Y10[i];
732                 }
733             }
734             return Double.NaN;
735         }
736 
737         @Override
738         public double getY10B(AbsoluteDate date)
739             {
740             for (int i = 0; i < TC.length; i++) {
741                 if (date.equals(TC[i])) {
742                     return Y10B[i];
743                 }
744             }
745             return Double.NaN;
746         }
747 
748         @Override
749         public double getDSTDTC(AbsoluteDate date)
750             {
751             for (int i = 0; i < TC.length; i++) {
752                 if (date.equals(TC[i])) {
753                     return DSTDTC[i];
754                 }
755             }
756             return Double.NaN;
757         }
758 
759     }
760 
761 }