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