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.data;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  import java.util.stream.Collectors;
22  
23  import org.hipparchus.analysis.UnivariateFunction;
24  import org.hipparchus.analysis.interpolation.LinearInterpolator;
25  import org.hipparchus.util.FastMath;
26  import org.orekit.annotation.DefaultDataContext;
27  import org.orekit.data.DataContext;
28  import org.orekit.data.DataProvidersManager;
29  import org.orekit.data.DataSource;
30  import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimationLoader.LineParameters;
31  import org.orekit.time.AbsoluteDate;
32  import org.orekit.time.DateComponents;
33  import org.orekit.time.DateTimeComponents;
34  import org.orekit.time.TimeComponents;
35  import org.orekit.time.TimeScale;
36  import org.orekit.time.TimeStampedDouble;
37  import org.orekit.utils.Constants;
38  import org.orekit.utils.GenericTimeStampedCache;
39  import org.orekit.utils.OrekitConfiguration;
40  import org.orekit.utils.TimeStampedGenerator;
41  
42  /**
43   * This class provides solar activity data needed by atmospheric models: F107 solar flux, Ap and Kp indexes.
44   * <p>
45   * Data comes from the NASA Marshall Solar Activity Future Estimation (MSAFE) as estimates of monthly F10.7
46   * Mean solar flux and Ap geomagnetic parameter
47   * (see <a href="https://www.nasa.gov/solar-cycle-progression-and-forecast"> Marshall Solar Activity website</a>).
48   *
49   * <p>Data can be retrieved at the NASA
50   * <a href="https://www.nasa.gov/solar-cycle-progression-and-forecast/archived-forecast/"> Marshall Solar Activity archived forecast</a>.
51   * Here Kp indices are deduced from Ap indexes, which in turn are tabulated equivalent of retrieved Ap values.
52   *
53   * <p> If several MSAFE files are available, some dates may appear in several files (for example August 2007 is in all files from
54   * the first one published in March 1999 to the February 2008 file). In this case, the data from the most recent file is used
55   * and the older ones are discarded. The date of the file is assumed to be 6 months after its first entry (which explains why
56   * the file having August 2007 as its first entry is the February 2008 file). This implies that MSAFE files must <em>not</em>
57   * be edited to change their time span, otherwise this would break the old entries overriding mechanism.
58   *
59   * <p>With these data, the {@link #getInstantFlux(AbsoluteDate)} and {@link #getMeanFlux(AbsoluteDate)} methods return the same
60   * values and the {@link #get24HoursKp(AbsoluteDate)} and {@link #getThreeHourlyKP(AbsoluteDate)} methods return the same
61   * values.
62   *
63   * <p>Conversion from Ap index values in the MSAFE file to Kp values used by atmosphere models is done using Jacchia's equation
64   * in [1].
65   *
66   * <p>With these data, the {@link #getAp(AbsoluteDate date)} method returns an array of seven times the same daily Ap value,
67   * i.e. it is suited to be used only with the {@link org.orekit.models.earth.atmosphere.NRLMSISE00 NRLMSISE00} atmospheric
68   * model where the switch #9 is set to 1.
69   *
70   * <h2>References</h2>
71   *
72   * <ol> <li> Jacchia, L. G. "CIRA 1972, recent atmospheric models, and improvements in
73   * progress." COSPAR, 21st Plenary Meeting. Vol. 1. 1978. </li> </ol>
74   *
75   * @author Bruno Revelin
76   * @author Luc Maisonobe
77   * @author Evan Ward
78   * @author Pascal Parraud
79   * @author Vincent Cucchietti
80   */
81  public class MarshallSolarActivityFutureEstimation
82          extends AbstractSolarActivityData<LineParameters, MarshallSolarActivityFutureEstimationLoader> {
83  
84      /**
85       * Default regular expression for the supported name that work with all officially published files.
86       *
87       * @since 10.0
88       */
89      public static final String DEFAULT_SUPPORTED_NAMES =
90              "\\p{Alpha}\\p{Lower}\\p{Lower}\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}(?:f|F)10(?:[-_]prd)?\\.(?:txt|TXT)";
91  
92      /** Serializable UID. */
93      private static final long serialVersionUID = -5212198874900835369L;
94  
95      /** Selected strength level of activity. */
96      private final StrengthLevel strengthLevel;
97  
98      /** Cache dedicated to average flux. */
99      private final transient GenericTimeStampedCache<TimeStampedDouble> averageFluxCache;
100 
101     /**
102      * Simple constructor. This constructor uses the {@link DataContext#getDefault() default data context}.
103      * <p>
104      * The original file names used by NASA Marshall space center are of the form: may2019f10_prd.txt or Oct1999F10.TXT. So a
105      * recommended regular expression for the supported name that work with all published files is:
106      * {@link #DEFAULT_SUPPORTED_NAMES}.
107      * <p>
108      * It provides a default configuration for the thread safe cache :
109      * <ul>
110      *     <li>Number of slots : {@code OrekitConfiguration.getCacheSlotsNumber()}</li>
111      *     <li>Max span : {@code Constants.JULIAN_YEAR}</li>
112      *     <li>Max span interval : {@code 0}</li>
113      * </ul>
114      *
115      * @param supportedNames regular expression for supported files names
116      * @param strengthLevel selected strength level of activity
117      *
118      * @see #MarshallSolarActivityFutureEstimation(String, StrengthLevel, DataProvidersManager, TimeScale)
119      */
120     @DefaultDataContext
121     public MarshallSolarActivityFutureEstimation(final String supportedNames,
122                                                  final StrengthLevel strengthLevel) {
123         this(supportedNames, strengthLevel,
124              DataContext.getDefault().getDataProvidersManager(),
125              DataContext.getDefault().getTimeScales().getUTC());
126     }
127 
128     /**
129      * Constructor that allows specifying the source of the MSAFE auxiliary data files.
130      * <p>
131      * It provides a default configuration for the thread safe cache :
132      * <ul>
133      *     <li>Number of slots : {@code OrekitConfiguration.getCacheSlotsNumber()}</li>
134      *     <li>Max span : {@code 31 * Constants.JULIAN_DAY}</li>
135      *     <li>Max interval : {@code 0}</li>
136      *     <li>Minimum step: {@code 27 * Constants.JULIAN_DAY}</li>
137      * </ul>
138      *
139      * @param supportedNames regular expression for supported files names
140      * @param strengthLevel selected strength level of activity
141      * @param dataProvidersManager provides access to auxiliary data files.
142      * @param utc UTC time scale.
143      *
144      * @since 10.1
145      */
146     public MarshallSolarActivityFutureEstimation(final String supportedNames,
147                                                  final StrengthLevel strengthLevel,
148                                                  final DataProvidersManager dataProvidersManager,
149                                                  final TimeScale utc) {
150         this(supportedNames, strengthLevel, dataProvidersManager, utc, OrekitConfiguration.getCacheSlotsNumber(),
151              Constants.JULIAN_DAY * 31, 0, Constants.JULIAN_DAY * 27);
152     }
153 
154     /**
155      * Constructor that allows specifying the source of the MSAFE auxiliary data files and customizable thread safe cache
156      * configuration.
157      *
158      * @param supportedNames regular expression for supported files names
159      * @param strengthLevel selected strength level of activity
160      * @param dataProvidersManager provides access to auxiliary data files.
161      * @param utc UTC time scale.
162      * @param maxSlots maximum number of independent cached time slots in the
163      * {@link GenericTimeStampedCache time-stamped cache}
164      * @param maxSpan maximum duration span in seconds of one slot in the {@link GenericTimeStampedCache time-stamped cache}
165      * @param maxInterval time interval above which a new slot is created in the
166      * {@link GenericTimeStampedCache time-stamped cache}
167      * @param minimumStep overriding minimum step designed for non-homogeneous tabulated values. To be used for example when
168      * caching monthly tabulated values. May be null.
169      *
170      * @since 10.1
171      */
172     public MarshallSolarActivityFutureEstimation(final String supportedNames,
173                                                  final StrengthLevel strengthLevel,
174                                                  final DataProvidersManager dataProvidersManager,
175                                                  final TimeScale utc,
176                                                  final int maxSlots,
177                                                  final double maxSpan,
178                                                  final double maxInterval,
179                                                  final double minimumStep) {
180         super(supportedNames, new MarshallSolarActivityFutureEstimationLoader(strengthLevel, utc),
181               dataProvidersManager, utc, maxSlots, maxSpan, maxInterval, minimumStep);
182 
183         // Initialise fields
184         this.strengthLevel    = strengthLevel;
185         this.averageFluxCache = new GenericTimeStampedCache<>(N_NEIGHBORS, OrekitConfiguration.getCacheSlotsNumber(),
186                                                               Constants.JULIAN_DAY, 0, new AverageFluxGenerator());
187     }
188 
189     /**
190      * Simple constructor which use the {@link DataContext#getDefault() default data context}.
191      * <p>
192      * It provides a default configuration for the thread safe cache :
193      * <ul>
194      *     <li>Number of slots : {@code OrekitConfiguration.getCacheSlotsNumber()}</li>
195      *     <li>Max span : {@code 31 * Constants.JULIAN_DAY}</li>
196      *     <li>Max interval : {@code 0}</li>
197      *     <li>Minimum step: {@code 27 * Constants.JULIAN_DAY}</li>
198      * </ul>
199      *
200      * @param source source for the data
201      * @param strengthLevel selected strength level of activity
202      *
203      * @since 12.0
204      */
205     @DefaultDataContext
206     public MarshallSolarActivityFutureEstimation(final DataSource source,
207                                                  final StrengthLevel strengthLevel) {
208         this(source, strengthLevel, DataContext.getDefault().getTimeScales().getUTC());
209     }
210 
211     /**
212      * Simple constructor.
213      * <p>
214      * It provides a default configuration for the thread safe cache :
215      * <ul>
216      *     <li>Number of slots : {@code OrekitConfiguration.getCacheSlotsNumber()}</li>
217      *     <li>Max span : {@code 31 * Constants.JULIAN_DAY}</li>
218      *     <li>Max interval : {@code 0}</li>
219      *     <li>Minimum step: {@code 27 * Constants.JULIAN_DAY}</li>
220      * </ul>
221      *
222      * @param source source for the data
223      * @param strengthLevel selected strength level of activity
224      * @param utc UTC time scale
225      *
226      * @since 12.0
227      */
228     public MarshallSolarActivityFutureEstimation(final DataSource source,
229                                                  final StrengthLevel strengthLevel,
230                                                  final TimeScale utc) {
231         this(source, strengthLevel, utc, OrekitConfiguration.getCacheSlotsNumber(),
232              Constants.JULIAN_DAY * 31, 0, Constants.JULIAN_DAY * 27);
233     }
234 
235     /**
236      * Constructor with customizable thread safe cache configuration.
237      *
238      * @param source source for the data
239      * @param strengthLevel selected strength level of activity
240      * @param utc UTC time scale
241      * @param maxSlots maximum number of independent cached time slots in the
242      * {@link GenericTimeStampedCache time-stamped cache}
243      * @param maxSpan maximum duration span in seconds of one slot in the {@link GenericTimeStampedCache time-stamped cache}
244      * @param maxInterval time interval above which a new slot is created in the
245      * {@link GenericTimeStampedCache time-stamped cache}
246      * @param minimumStep overriding minimum step designed for non-homogeneous tabulated values. To be used for example when
247      * caching monthly tabulated values. Use {@code Double.NaN} otherwise.
248      *
249      * @since 12.0
250      */
251     public MarshallSolarActivityFutureEstimation(final DataSource source,
252                                                  final StrengthLevel strengthLevel,
253                                                  final TimeScale utc,
254                                                  final int maxSlots,
255                                                  final double maxSpan,
256                                                  final double maxInterval,
257                                                  final double minimumStep) {
258         super(source, new MarshallSolarActivityFutureEstimationLoader(strengthLevel, utc), utc,
259               maxSlots, maxSpan, maxInterval, minimumStep);
260         this.strengthLevel    = strengthLevel;
261         this.averageFluxCache = new GenericTimeStampedCache<>(N_NEIGHBORS, OrekitConfiguration.getCacheSlotsNumber(),
262                                                               Constants.JULIAN_DAY, 0, new AverageFluxGenerator());
263     }
264 
265     /** {@inheritDoc} */
266     public double getInstantFlux(final AbsoluteDate date) {
267         return getMeanFlux(date);
268     }
269 
270     /** {@inheritDoc} */
271     public double getMeanFlux(final AbsoluteDate date) {
272         return getLinearInterpolation(date, LineParameters::getF107);
273     }
274 
275     /** {@inheritDoc} */
276     public double getThreeHourlyKP(final AbsoluteDate date) {
277         return get24HoursKp(date);
278     }
279 
280     /**
281      * Get the date of the file from which data at the specified date comes from.
282      * <p>
283      * If several MSAFE files are available, some dates may appear in several files (for example August 2007 is in all files
284      * from the first one published in March 1999 to the February 2008 file). In this case, the data from the most recent
285      * file is used and the older ones are discarded. The date of the file is assumed to be 6 months after its first entry
286      * (which explains why the file having August 2007 as its first entry is the February 2008 file). This implies that MSAFE
287      * files must <em>not</em> be edited to change their time span, otherwise this would break the old entries overriding
288      * mechanism.
289      * </p>
290      *
291      * @param date date of the solar activity data
292      *
293      * @return date of the file
294      */
295     public DateComponents getFileDate(final AbsoluteDate date) {
296         // Get the neighboring solar activity
297         final LocalSolarActivity localSolarActivity = new LocalSolarActivity(date);
298         final LineParameters     previousParam      = localSolarActivity.getPreviousParam();
299         final LineParameters     currentParam       = localSolarActivity.getNextParam();
300 
301         // Choose which file date to return
302         final double dtP = date.durationFrom(previousParam.getDate());
303         final double dtC = currentParam.getDate().durationFrom(date);
304         return (dtP < dtC) ? previousParam.getFileDate() : currentParam.getFileDate();
305     }
306 
307     /**
308      * The Kp index is derived from the Ap index.
309      * <p>The method used is explained on <a
310      * href="http://www.ngdc.noaa.gov/stp/GEOMAG/kp_ap.html"> NOAA website.</a> as follows:</p>
311      * <p>The scale is 0 to 9 expressed in thirds of a unit, e.g. 5- is 4 2/3,
312      * 5 is 5 and 5+ is 5 1/3. The ap (equivalent range) index is derived from the Kp index as follows:</p>
313      * <table border="1">
314      * <caption>Kp / Ap Conversion Table</caption>
315      * <tbody>
316      * <tr>
317      * <td>Kp</td><td>0o</td><td>0+</td><td>1-</td><td>1o</td><td>1+</td><td>2-</td><td>2o</td><td>2+</td><td>3-</td><td>3o</td><td>3+</td><td>4-</td><td>4o</td><td>4+</td>
318      * </tr>
319      * <tr>
320      * <td>ap</td><td>0</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>9</td><td>12</td><td>15</td><td>18</td><td>22</td><td>27</td><td>32</td>
321      * </tr>
322      * <tr>
323      * <td>Kp</td><td>5-</td><td>5o</td><td>5+</td><td>6-</td><td>6o</td><td>6+</td><td>7-</td><td>7o</td><td>7+</td><td>8-</td><td>8o</td><td>8+</td><td>9-</td><td>9o</td>
324      * </tr>
325      * <tr>
326      * <td>ap</td><td>39</td><td>48</td><td>56</td><td>67</td><td>80</td><td>94</td><td>111</td><td>132</td><td>154</td><td>179</td><td>207</td><td>236</td><td>300</td><td>400</td>
327      * </tr>
328      * </tbody>
329      * </table>
330      *
331      * @param date date of the Kp data
332      *
333      * @return the 24H geomagnetic index
334      */
335     public double get24HoursKp(final AbsoluteDate date) {
336         // get the daily Ap
337         final double ap = getDailyAp(date);
338 
339         // get the corresponding Kp index from
340         // equation 4 in [1] for Ap to Kp conversion
341         return 1.89 * FastMath.asinh(0.154 * ap);
342     }
343 
344     /** {@inheritDoc} */
345     public double getDailyFlux(final AbsoluteDate date) {
346         return getMeanFlux(date.shiftedBy(-Constants.JULIAN_DAY));
347     }
348 
349     public double getAverageFlux(final AbsoluteDate date) {
350         // Extract closest neighbours average
351         final List<TimeStampedDouble> neighbors = averageFluxCache.getNeighbors(date).collect(Collectors.toList());
352 
353         // Create linear interpolating function
354         final double[] x = new double[] { 0, 1 };
355         final double[] y = neighbors.stream().map(TimeStampedDouble::getValue).mapToDouble(Double::doubleValue).toArray();
356 
357         final LinearInterpolator interpolator          = new LinearInterpolator();
358         final UnivariateFunction interpolatingFunction = interpolator.interpolate(x, y);
359 
360         // Interpolate
361         final AbsoluteDate previousDate = neighbors.get(0).getDate();
362         final AbsoluteDate nextDate     = neighbors.get(1).getDate();
363         return interpolatingFunction.value(date.durationFrom(previousDate) / nextDate.durationFrom(previousDate));
364     }
365 
366     /** {@inheritDoc} */
367     public double[] getAp(final AbsoluteDate date) {
368         // Gets the AP for the current date
369         final double ap = getDailyAp(date);
370 
371         // Retuns an array of Ap filled in with the daily Ap only
372         return new double[] { ap, ap, ap, ap, ap, ap, ap };
373     }
374 
375     /**
376      * Gets the daily Ap index.
377      *
378      * @param date the current date
379      *
380      * @return the daily Ap index
381      */
382     private double getDailyAp(final AbsoluteDate date) {
383         return getLinearInterpolation(date, LineParameters::getAp);
384     }
385 
386     /**
387      * Get the strength level for activity.
388      *
389      * @return strength level to set
390      */
391     public StrengthLevel getStrengthLevel() {
392         return strengthLevel;
393     }
394 
395     /** Strength level of activity. */
396     public enum StrengthLevel {
397 
398         /** Strong level of activity. */
399         STRONG,
400 
401         /** Average level of activity. */
402         AVERAGE,
403 
404         /** Weak level of activity. */
405         WEAK
406 
407     }
408 
409     /** Generator generating average flux data between given dates. */
410     private class AverageFluxGenerator implements TimeStampedGenerator<TimeStampedDouble> {
411 
412         /** {@inheritDoc} */
413         @Override
414         public List<TimeStampedDouble> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
415             // No prior data in the cache
416             if (existingDate == null) {
417                 return generateDataFromEarliestToLatestDates(getCurrentDay(date), getNextDay(date));
418             }
419             // Prior data in the cache, fill with data between date and existing date
420             if (date.isBefore(existingDate)) {
421                 return generateDataFromEarliestToLatestDates(date, existingDate);
422             }
423             return generateDataFromEarliestToLatestDates(existingDate, date);
424         }
425 
426         /**
427          * Generate data from earliest to latest date.
428          *
429          * @param earliest earliest date
430          * @param latest latest date
431          *
432          * @return data generated from earliest to latest date
433          */
434         private List<TimeStampedDouble> generateDataFromEarliestToLatestDates(final AbsoluteDate earliest,
435                                                                               final AbsoluteDate latest) {
436             final List<TimeStampedDouble> generated = new ArrayList<>();
437 
438             // Add next computed average until it brackets the latest date
439             AbsoluteDate latestNeighbourDate = getCurrentDay(earliest);
440             while (latestNeighbourDate.isBeforeOrEqualTo(latest)) {
441                 generated.add(computeAverageFlux(latestNeighbourDate));
442                 latestNeighbourDate = getNextDay(latestNeighbourDate);
443             }
444             return generated;
445         }
446 
447         /**
448          * Get the current day at midnight.
449          *
450          * @param date date
451          *
452          * @return current day at midnight.
453          */
454         private AbsoluteDate getCurrentDay(final AbsoluteDate date) {
455             // Find previous day date time components
456             final TimeScale      utc            = getUTC();
457             final DateComponents dateComponents = date.getComponents(utc).getDate();
458 
459             // Create absolute date for previous day
460             return new AbsoluteDate(new DateTimeComponents(dateComponents, TimeComponents.H00), utc);
461         }
462 
463         /**
464          * Get the next day at midnight.
465          *
466          * @param date date
467          *
468          * @return next day at midnight.
469          */
470         private AbsoluteDate getNextDay(final AbsoluteDate date) {
471             // Find previous day date time components
472             final TimeScale      utc               = getUTC();
473             final DateComponents dateComponents    = date.getComponents(utc).getDate();
474             final DateComponents shiftedComponents = new DateComponents(dateComponents, 1);
475 
476             // Create absolute date for previous day
477             return new AbsoluteDate(new DateTimeComponents(shiftedComponents, TimeComponents.H00), utc);
478         }
479 
480         /**
481          * Compute the average flux for given absolute date.
482          *
483          * @param date date at which the average flux is desired
484          *
485          * @return average flux
486          */
487         private TimeStampedDouble computeAverageFlux(final AbsoluteDate date) {
488             // Extract list of neighbors to compute average
489             final TimeStampedGenerator<LineParameters> generator   = getCache().getGenerator();
490             final AbsoluteDate                         initialDate = date.shiftedBy(-40 * Constants.JULIAN_DAY);
491             final AbsoluteDate                         finalDate   = date.shiftedBy(40 * Constants.JULIAN_DAY);
492             final List<LineParameters>                 monthlyData = generator.generate(initialDate, finalDate);
493 
494             // Create interpolator for given data
495             final LinearInterpolator interpolator = new LinearInterpolator();
496 
497             final double[] x = monthlyData.stream().map(param -> param.getDate().durationFrom(initialDate))
498                                           .mapToDouble(Double::doubleValue).toArray();
499             final double[] y = monthlyData.stream().map(LineParameters::getF107).mapToDouble(Double::doubleValue).toArray();
500 
501             final UnivariateFunction interpolatingFunction = interpolator.interpolate(x, y);
502 
503             // Loops over the 81 days centered on the given date
504             double average = 0;
505             for (int i = -40; i < 41; i++) {
506                 average += interpolatingFunction.value(date.shiftedBy(i * Constants.JULIAN_DAY).durationFrom(initialDate));
507             }
508 
509             // Returns the 81 day average flux
510             return new TimeStampedDouble(average / 81, date);
511         }
512     }
513 
514 }