1   /*
2    * Copyright 2020-2025 Clément Jonglez
3    * Licensed to CS GROUP (CS) under one or more
4    * contributor license agreements.  See the NOTICE file distributed with
5    * this work for additional information regarding copyright ownership.
6    * Clément Jonglez licenses this file to You under the Apache License, Version 2.0
7    * (the "License"); you may not use this file except in compliance with
8    * the License.  You may obtain a copy of the License at
9    *
10   *   http://www.apache.org/licenses/LICENSE-2.0
11   *
12   * Unless required by applicable law or agreed to in writing, software
13   * distributed under the License is distributed on an "AS IS" BASIS,
14   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15   * See the License for the specific language governing permissions and
16   * limitations under the License.
17   */
18  
19  package org.orekit.models.earth.atmosphere.data;
20  
21  
22  import static org.hamcrest.MatcherAssert.assertThat;
23  import static org.junit.jupiter.api.Timeout.ThreadMode.SEPARATE_THREAD;
24  import static org.orekit.OrekitMatchers.closeTo;
25  import static org.orekit.OrekitMatchers.pvCloseTo;
26  
27  import java.net.URISyntaxException;
28  import java.net.URL;
29  import java.util.ArrayList;
30  import java.util.Comparator;
31  import java.util.List;
32  import java.util.SortedSet;
33  import java.util.concurrent.Callable;
34  import java.util.concurrent.ExecutorService;
35  import java.util.concurrent.Executors;
36  import java.util.concurrent.TimeUnit;
37  import java.util.concurrent.atomic.AtomicReference;
38  import java.util.stream.Collectors;
39  
40  import org.hipparchus.ode.ODEIntegrator;
41  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
42  import org.hipparchus.util.FastMath;
43  import org.junit.jupiter.api.Assertions;
44  import org.junit.jupiter.api.BeforeEach;
45  import org.junit.jupiter.api.DisplayName;
46  import org.junit.jupiter.api.RepeatedTest;
47  import org.junit.jupiter.api.Test;
48  import org.junit.jupiter.api.Timeout;
49  import org.orekit.Utils;
50  import org.orekit.bodies.CelestialBody;
51  import org.orekit.bodies.CelestialBodyFactory;
52  import org.orekit.bodies.OneAxisEllipsoid;
53  import org.orekit.data.DataContext;
54  import org.orekit.data.DataProvidersManager;
55  import org.orekit.data.DataSource;
56  import org.orekit.errors.OrekitException;
57  import org.orekit.errors.OrekitMessages;
58  import org.orekit.forces.drag.DragForce;
59  import org.orekit.forces.drag.IsotropicDrag;
60  import org.orekit.frames.Frame;
61  import org.orekit.frames.FramesFactory;
62  import org.orekit.models.earth.atmosphere.Atmosphere;
63  import org.orekit.models.earth.atmosphere.DTM2000;
64  import org.orekit.models.earth.atmosphere.DTM2000InputParameters;
65  import org.orekit.models.earth.atmosphere.NRLMSISE00;
66  import org.orekit.models.earth.atmosphere.NRLMSISE00InputParameters;
67  import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherDataLoader.LineParameters;
68  import org.orekit.orbits.KeplerianOrbit;
69  import org.orekit.orbits.Orbit;
70  import org.orekit.orbits.OrbitType;
71  import org.orekit.orbits.PositionAngleType;
72  import org.orekit.propagation.SpacecraftState;
73  import org.orekit.propagation.numerical.NumericalPropagator;
74  import org.orekit.propagation.sampling.OrekitStepHandler;
75  import org.orekit.time.AbsoluteDate;
76  import org.orekit.time.TimeScale;
77  import org.orekit.time.TimeScalesFactory;
78  import org.orekit.time.TimeStampedDouble;
79  import org.orekit.utils.Constants;
80  import org.orekit.utils.GenericTimeStampedCache;
81  import org.orekit.utils.IERSConventions;
82  
83  /**
84   *
85   * @author Clément Jonglez
86   */
87  public class CssiSpaceWeatherLoaderTest {
88      private TimeScale utc;
89  
90      @BeforeEach
91      public void setUp() {
92          Utils.setDataRoot("regular-data:atmosphere");
93          utc = TimeScalesFactory.getUTC();
94      }
95  
96      private CssiSpaceWeatherData loadCswl() {
97          return new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES);
98      }
99  
100     @Test
101     public void testIssue1117() throws URISyntaxException {
102         final URL url = CssiSpaceWeatherLoaderTest.class.getClassLoader().getResource("atmosphere/SpaceWeather-All-v1.2_reduced.txt");
103         CssiSpaceWeatherData cswl = new CssiSpaceWeatherData(new DataSource(url.toURI()));
104         Assertions.assertEquals(new AbsoluteDate("2020-02-19", utc), cswl.getMinDate());
105         Assertions.assertEquals(new AbsoluteDate("2020-02-22", utc), cswl.getMaxDate());
106     }
107 
108     @Test
109     public void testMinDate() {
110         CssiSpaceWeatherData cswl = loadCswl();
111         Assertions.assertEquals(new AbsoluteDate("1957-10-01", utc), cswl.getMinDate());
112     }
113 
114     @Test
115     public void testMaxDate() {
116         CssiSpaceWeatherData cswl = loadCswl();
117         Assertions.assertEquals(new AbsoluteDate("2044-06-01", utc), cswl.getMaxDate());
118     }
119 
120     @Test
121     public void testThreeHourlyKpObserved() {
122         CssiSpaceWeatherData cswl = loadCswl();
123         AbsoluteDate date = new AbsoluteDate(1957, 10, 10, 0, 0, 0.0, utc);
124         final double kp = cswl.getThreeHourlyKP(date);
125         assertThat(kp, closeTo(3.0, 1e-10));
126     }
127 
128     /** Requests a date between two months, requiring interpolation */
129     @Test
130     public void testThreeHourlyKpMonthlyPredicted() {
131         CssiSpaceWeatherData cswl = loadCswl();
132         AbsoluteDate date = new AbsoluteDate(2038, 6, 16, 0, 0, 0.0, utc);
133         final double kp = cswl.getThreeHourlyKP(date);
134         assertThat(kp, closeTo((2.7 + 4.1) / 2, 1e-3));
135     }
136 
137     /**
138      * Testing first day of data
139      * Because the Ap up to 57 hours prior to the date are asked,
140      * this should return an exception
141      */
142     @Test
143     public void testThreeHourlyApObservedFirstDay() {
144         CssiSpaceWeatherData cswl = loadCswl();
145         AbsoluteDate date = new AbsoluteDate(1957, 10, 1, 5, 17, 0.0, utc);
146         try {
147             cswl.getAp(date);
148             Assertions.fail("an exception should have been thrown");
149         } catch (OrekitException oe) {
150             Assertions.assertEquals(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE, oe.getSpecifier());
151         }
152     }
153 
154     /**
155      * Testing second day of data
156      * Because the Ap up to 57 hours prior to the date are asked,
157      * this should return an exception
158      */
159     @Test
160     public void testThreeHourlyApObservedSecondDay() {
161         CssiSpaceWeatherData cswl = loadCswl();
162         AbsoluteDate date = new AbsoluteDate(1957, 10, 2, 3, 14, 0.0, utc);
163         try {
164             cswl.getAp(date);
165             Assertions.fail("an exception should have been thrown");
166         } catch (OrekitException oe) {
167             Assertions.assertEquals(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE, oe.getSpecifier());
168         }
169     }
170 
171     /**
172      * Testing third day of data
173      * Because the Ap up to 57 hours prior to the date are asked,
174      * this should return an exception
175      */
176     @Test
177     public void testThreeHourlyApObservedThirdDay() {
178         CssiSpaceWeatherData cswl = loadCswl();
179         AbsoluteDate date = new AbsoluteDate(1957, 10, 3, 3, 14, 0.0, utc);
180         try {
181             cswl.getAp(date);
182             Assertions.fail("an exception should have been thrown");
183         } catch (OrekitException oe) {
184             Assertions.assertEquals(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE, oe.getSpecifier());
185         }
186     }
187 
188     /** Here, no more side effects are expected */
189     @Test
190     public void testThreeHourlyApObserved() {
191         CssiSpaceWeatherData cswl = loadCswl();
192         AbsoluteDate date = new AbsoluteDate(1957, 10, 10, 3, 14, 0.0, utc);
193         final double[] apExpected = new double[] { 18, 27, 15, 9, 7, 5.625, 3.625 };
194         final double[] ap = cswl.getAp(date);
195         for (int i = 0; i < 7; i++) {
196             assertThat(ap[i], closeTo(apExpected[i], 1e-10));
197         }
198     }
199 
200     /**
201      * This test is very approximate, at least to check that the two proper months were used for the interpolation
202      * But the manual interpolation of all 7 coefficients would have been a pain
203      */
204     @Test
205     public void testThreeHourlyApMonthlyPredicted() {
206         CssiSpaceWeatherData cswl = loadCswl();
207         AbsoluteDate date = new AbsoluteDate(2038, 6, 16, 0, 0, 0.0, utc);
208         final double[] ap = cswl.getAp(date);
209         for (int i = 0; i < 7; i++) {
210             assertThat(ap[i], closeTo((12 + 29) / 2, 1.0));
211         }
212     }
213 
214     @Test
215     public void testDailyFlux() {
216         CssiSpaceWeatherData cswl = loadCswl();
217         AbsoluteDate date = new AbsoluteDate(2000, 1, 1, 0, 0, 0.0, utc);
218         final double dailyFlux = cswl.getDailyFlux(date);
219         // The daily flux is supposed to be the one from 1999-12-31
220         assertThat(dailyFlux, closeTo(130.1, 1e-10));
221     }
222 
223     @Test
224     public void testDailyFluxMonthlyPredicted() {
225         CssiSpaceWeatherData cswl = loadCswl();
226         AbsoluteDate date = new AbsoluteDate(2034, 6, 17, 0, 0, 0.0, utc);
227         final double dailyFlux = cswl.getDailyFlux(date);
228         assertThat(dailyFlux, closeTo((132.7 + 137.3) / 2, 1e-3));
229     }
230 
231     @Test
232     public void testMeanFlux() {
233         CssiSpaceWeatherData cswl = loadCswl();
234         AbsoluteDate date = new AbsoluteDate(2000, 1, 1, 0, 0, 0.0, utc);
235         final double meanFlux = cswl.getMeanFlux(date);
236         assertThat(meanFlux, closeTo(165.6, 1e-10));
237     }
238 
239     @Test
240     public void testMeanFluxMonthlyPredicted() {
241         CssiSpaceWeatherData cswl = loadCswl();
242         AbsoluteDate date = new AbsoluteDate(2034, 6, 16, 0, 0, 0.0, utc);
243         final double meanFlux = cswl.getMeanFlux(date);
244         assertThat(meanFlux, closeTo((132.1 + 134.9) / 2, 1e-3));
245     }
246 
247     @Test
248     public void testAverageFlux() {
249         CssiSpaceWeatherData cswl = loadCswl();
250         AbsoluteDate date = new AbsoluteDate(2000, 1, 1, 0, 0, 0.0, utc);
251         final double averageFlux = cswl.getAverageFlux(date);
252         assertThat(averageFlux, closeTo(165.6, 1e-10));
253     }
254 
255     @Test
256     public void testAverageFluxMonthlyPredicted() {
257         CssiSpaceWeatherData cswl = loadCswl();
258         AbsoluteDate date = new AbsoluteDate(2034, 6, 16, 0, 0, 0.0, utc);
259         final double averageFlux = cswl.getAverageFlux(date);
260         assertThat(averageFlux, closeTo((132.1 + 134.9) / 2, 1e-3));
261     }
262 
263     @Test
264     public void testInstantFlux() {
265         CssiSpaceWeatherData cswl = loadCswl();
266         AbsoluteDate date = new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, utc);
267         final double instantFlux = cswl.getInstantFlux(date);
268         assertThat(instantFlux, closeTo((129.9 + 132.9) / 2, 1e-10));
269     }
270 
271     /**
272      * This is to test that daily Kp values with three figures are correctly parsed
273      */
274     @Test
275     public void testDailyKp2003SolarStorm() {
276         CssiSpaceWeatherData cswl = loadCswl();
277         AbsoluteDate date = new AbsoluteDate(2003, 10, 29, 23, 0, 0.0, utc);
278         final double kp = cswl.get24HoursKp(date);
279         assertThat(kp, closeTo(583.0 / 10.0 / 8, 1e-10));
280     }
281 
282     /**
283      * This is to test that Ap values with three figures are correctly parsed
284      */
285     @Test
286     public void testAp2003SolarStorm() {
287         CssiSpaceWeatherData cswl = loadCswl();
288         AbsoluteDate date = new AbsoluteDate(2003, 10, 29, 23, 0, 0.0, utc);
289         final double[] ap = cswl.getAp(date);
290         assertThat(ap[0], closeTo(204, 1e-10));
291         assertThat(ap[1], closeTo(300, 1e-10));
292         assertThat(ap[2], closeTo(300, 1e-10));
293         assertThat(ap[3], closeTo(179, 1e-10));
294         assertThat(ap[4], closeTo(179, 1e-10));
295     }
296 
297     /**
298      * Check integration error is small when integrating the same equations over the same
299      * interval.
300      */
301     @Test
302     public void testWithPropagator() {
303         CelestialBody sun = CelestialBodyFactory.getSun();
304         final Frame eci = FramesFactory.getGCRF();
305         final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
306         AbsoluteDate date = new AbsoluteDate(2004, 1, 1, utc);
307         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
308                 Constants.WGS84_EARTH_FLATTENING, ecef);
309         Orbit orbit = new KeplerianOrbit(6378137 + 400e3, 1e-3, FastMath.toRadians(50), 0, 0, 0, PositionAngleType.TRUE,
310                 eci, date, Constants.EIGEN5C_EARTH_MU);
311         final SpacecraftState ic = new SpacecraftState(orbit);
312 
313         final AbsoluteDate end = date.shiftedBy(5 * Constants.JULIAN_DAY);
314         final AbsoluteDate resetDate = date.shiftedBy(0.8 * Constants.JULIAN_DAY + 0.1);
315 
316         final SpacecraftState[] lastState = new SpacecraftState[1];
317         final OrekitStepHandler stepSaver = interpolator -> {
318             final AbsoluteDate start = interpolator.getPreviousState().getDate();
319             if (start.compareTo(resetDate) < 0) {
320                 lastState[0] = interpolator.getPreviousState();
321             }
322         };
323 
324         // propagate with state rest to take slightly different path
325         NumericalPropagator propagator = getNumericalPropagatorWithDTM(sun, earth, ic);
326         propagator.setStepHandler(stepSaver);
327         propagator.propagate(resetDate);
328         propagator.resetInitialState(lastState[0]);
329         propagator.clearStepHandlers();
330         SpacecraftState actual = propagator.propagate(end);
331 
332         // propagate straight through
333         propagator = getNumericalPropagatorWithDTM(sun, earth, ic);
334         propagator.resetInitialState(ic);
335         propagator.clearStepHandlers();
336         SpacecraftState expected = propagator.propagate(end);
337 
338         assertThat(actual.getPVCoordinates(), pvCloseTo(expected.getPVCoordinates(), 1.0));
339 
340         // propagate with state rest to take slightly different path
341         propagator = getNumericalPropagatorWithMSIS(sun, earth, ic);
342         propagator.setStepHandler(stepSaver);
343         propagator.propagate(resetDate);
344         propagator.resetInitialState(lastState[0]);
345         propagator.clearStepHandlers();
346         actual = propagator.propagate(end);
347 
348         // propagate straight through
349         propagator = getNumericalPropagatorWithMSIS(sun, earth, ic);
350         propagator.resetInitialState(ic);
351         propagator.clearStepHandlers();
352         expected = propagator.propagate(end);
353 
354         assertThat(actual.getPVCoordinates(), pvCloseTo(expected.getPVCoordinates(), 1.0));
355     }
356 
357     @Test
358     public void testIssue841() throws OrekitException {
359         final CssiSpaceWeatherDataLoader loader = new CssiSpaceWeatherDataLoader(utc);
360         DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
361         manager.feed("SpaceWeather-All-v1.2_reduced.txt", loader);
362         final SortedSet<LineParameters> set = loader.getDataSet();
363         Assertions.assertEquals(4, set.size());
364 
365         CssiSpaceWeatherData cswl = new CssiSpaceWeatherData("SpaceWeather-All-v1.2_reduced.txt");
366         Assertions.assertEquals(71.6, cswl.getInstantFlux(new AbsoluteDate("2020-02-20T00:00:00.000", utc)), 0.01);
367     }
368 
369     /**
370      * This test in a multi-threaded environment would not necessarily fail without the fix (even though it will very likely
371      * fail).
372      * <p>
373      * However, it cannot fail with the fix.
374      */
375     @RepeatedTest(10)
376     @DisplayName("Test in a multi-threaded environment")
377     void testIssue1072() {
378         // GIVEN
379         final CssiSpaceWeatherData weatherData = loadCswl();
380 
381         // Create date sample at which flux will be evaluated
382         final AbsoluteDate       initialDate = new AbsoluteDate();
383         final List<AbsoluteDate> dates       = new ArrayList<>();
384         final int                sampleSize  = 100;
385         for (int i = 0; i < sampleSize + 1; i++) {
386             dates.add(initialDate.shiftedBy(i * Constants.JULIAN_DAY * 30));
387         }
388 
389         // Create list of tasks to run in parallel
390         final AtomicReference<List<TimeStampedDouble>> results = new AtomicReference<>(new ArrayList<>());
391         final List<Callable<List<TimeStampedDouble>>>  tasks   = new ArrayList<>();
392         for (int i = 0; i < sampleSize + 1; i++) {
393             final AbsoluteDate currentDate = dates.get(i);
394             // Each task will evaluate value at specific date and store this value and associated date in a shared list
395             tasks.add(() -> (results.getAndUpdate((listToUpdate) -> {
396                 final List<TimeStampedDouble> newList = new ArrayList<>(listToUpdate);
397                 newList.add(new TimeStampedDouble(weatherData.get24HoursKp(currentDate), currentDate));
398                 return newList;
399             })));
400         }
401 
402         // Create multithreading environment
403         ExecutorService service = Executors.newFixedThreadPool(sampleSize);
404 
405         // WHEN & THEN
406         try {
407             service.invokeAll(tasks);
408             results.get().sort(Comparator.comparing(TimeStampedDouble::getDate));
409             final List<Double> sortedComputedResults = results.get().stream().map(TimeStampedDouble::getValue).collect(
410                     Collectors.toList());
411 
412             // Compare to expected result
413             for (int i = 0; i < sampleSize + 1; i++) {
414                 final AbsoluteDate currentDate = dates.get(i);
415                 Assertions.assertEquals(weatherData.get24HoursKp(currentDate), sortedComputedResults.get(i));
416             }
417 
418             try {
419                 // wait for proper ending
420                 service.shutdown();
421                 service.awaitTermination(5, TimeUnit.SECONDS);
422             } catch (InterruptedException ie) {
423                 // Restore interrupted state...
424                 Thread.currentThread().interrupt();
425             }
426         }
427         catch (Exception e) {
428             // Should not fail
429             Assertions.fail();
430         }
431     }
432 
433     @Test
434     @Timeout(value = 5, threadMode = SEPARATE_THREAD)
435     public void testIssue1269() {
436         // GIVEN
437         final NRLMSISE00InputParameters solarActivity =
438                 new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES);
439 
440         // WHEN & THEN
441         solarActivity.getAverageFlux(new AbsoluteDate("2025-02-01T00:00:00.000", TimeScalesFactory.getUTC()));
442         solarActivity.getDailyFlux(new AbsoluteDate("2025-02-01T00:00:00.000", TimeScalesFactory.getUTC()));
443         solarActivity.getAp(new AbsoluteDate("2025-02-01T00:00:00.000", TimeScalesFactory.getUTC()));
444 
445     }
446 
447     @Test
448     void testExpectedCacheConfigurationAndCalls() {
449         // GIVEN
450         final AbsoluteDate date = new AbsoluteDate(2020, 2, 25, 2, 0, 0, TimeScalesFactory.getUTC());
451 
452         final CssiSpaceWeatherData atm = new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES);
453 
454         // WHEN
455         // Call flux at instants that shall generate slots
456         atm.getInstantFlux(date.shiftedBy(-1*Constants.JULIAN_DAY));
457         atm.getInstantFlux(date);
458         atm.getInstantFlux(date.shiftedBy(1*Constants.JULIAN_DAY));
459         atm.getInstantFlux(date.shiftedBy(3*Constants.JULIAN_DAY));
460         atm.getInstantFlux(date.shiftedBy(2*Constants.JULIAN_DAY));
461 
462         // Call flux at instants that shall not generate slots
463         atm.getInstantFlux(date.shiftedBy(-0.6 * Constants.JULIAN_DAY));
464         atm.getInstantFlux(date.shiftedBy(1.8 * Constants.JULIAN_DAY));
465 
466         // THEN
467         final GenericTimeStampedCache<LineParameters> cache = atm.getCache();
468         Assertions.assertEquals(5, cache.getGenerateCalls());
469         Assertions.assertEquals(5, cache.getSlots());
470         Assertions.assertEquals(10, cache.getEntries());
471         Assertions.assertEquals(7, cache.getGetNeighborsCalls());
472     }
473 
474     /**
475      * Configure a numerical propagator with DTM2000 atmosphere.
476      *
477      * @param sun   Sun.
478      * @param earth Earth.
479      * @param ic    initial condition.
480      * @return a propagator with DTM2000 atmosphere.
481      */
482     private NumericalPropagator getNumericalPropagatorWithDTM(CelestialBody sun, OneAxisEllipsoid earth,
483             SpacecraftState ic) {
484         // some non-integer step size to induce truncation error in flux interpolation
485         final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(120 + 0.1);
486         NumericalPropagator propagator = new NumericalPropagator(integrator);
487         DTM2000InputParameters flux = loadCswl();
488         final Atmosphere atmosphere = new DTM2000(flux, sun, earth);
489         final IsotropicDrag satellite = new IsotropicDrag(1, 3.2);
490         propagator.addForceModel(new DragForce(atmosphere, satellite));
491 
492         propagator.setInitialState(ic);
493         propagator.setOrbitType(OrbitType.CARTESIAN);
494 
495         return propagator;
496     }
497 
498     /**
499      * Configure a numerical propagator with NRLMSISE00 atmosphere.
500      *
501      * @param sun   Sun.
502      * @param earth Earth.
503      * @param ic    initial condition.
504      * @return a propagator with NRLMSISE00 atmosphere.
505      */
506     private NumericalPropagator getNumericalPropagatorWithMSIS(CelestialBody sun, OneAxisEllipsoid earth,
507             SpacecraftState ic) {
508         // some non-integer step size to induce truncation error in flux interpolation
509         final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(120 + 0.1);
510         NumericalPropagator propagator = new NumericalPropagator(integrator);
511         NRLMSISE00InputParameters flux = loadCswl();
512         final Atmosphere atmosphere = new NRLMSISE00(flux, sun, earth);
513         final IsotropicDrag satellite = new IsotropicDrag(1, 3.2);
514         propagator.addForceModel(new DragForce(atmosphere, satellite));
515 
516         propagator.setInitialState(ic);
517         propagator.setOrbitType(OrbitType.CARTESIAN);
518 
519         return propagator;
520     }
521 }