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 }