1 /* Copyright 2002-2021 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.BufferedReader;
20 import java.io.IOException;
21 import java.io.InputStream;
22 import java.io.InputStreamReader;
23 import java.io.Serializable;
24 import java.nio.charset.StandardCharsets;
25 import java.text.ParseException;
26 import java.util.Iterator;
27 import java.util.SortedSet;
28 import java.util.TreeSet;
29 import java.util.regex.Matcher;
30 import java.util.regex.Pattern;
31
32 import org.hipparchus.util.FastMath;
33 import org.orekit.annotation.DefaultDataContext;
34 import org.orekit.data.AbstractSelfFeedingLoader;
35 import org.orekit.data.DataContext;
36 import org.orekit.data.DataLoader;
37 import org.orekit.data.DataProvidersManager;
38 import org.orekit.errors.OrekitException;
39 import org.orekit.errors.OrekitInternalError;
40 import org.orekit.errors.OrekitMessages;
41 import org.orekit.models.earth.atmosphere.DTM2000InputParameters;
42 import org.orekit.models.earth.atmosphere.NRLMSISE00InputParameters;
43 import org.orekit.time.AbsoluteDate;
44 import org.orekit.time.ChronologicalComparator;
45 import org.orekit.time.DateComponents;
46 import org.orekit.time.Month;
47 import org.orekit.time.TimeScale;
48 import org.orekit.time.TimeStamped;
49 import org.orekit.utils.Constants;
50
51 /**
52 * This class reads and provides solar activity data needed by
53 * atmospheric models: F107 solar flux, Ap and Kp indexes.
54 * <p>
55 * The data are retrieved through the NASA Marshall
56 * Solar Activity Future Estimation (MSAFE) as estimates of monthly
57 * F10.7 Mean solar flux and Ap geomagnetic parameter.
58 * The data can be retrieved at the NASA <a
59 * href="http://sail.msfc.nasa.gov/archive_index.htm">
60 * Marshall Solar Activity website</a>.
61 * Here Kp indices are deduced from Ap indexes, which in turn are tabulated
62 * equivalent of retrieved Ap values.
63 * </p>
64 * <p>
65 * If several MSAFE files are available, some dates may appear in several
66 * files (for example August 2007 is in all files from the first one
67 * published in March 1999 to the February 2008 file). In this case, the
68 * data from the most recent file is used and the older ones are discarded.
69 * The date of the file is assumed to be 6 months after its first entry
70 * (which explains why the file having August 2007 as its first entry is the
71 * February 2008 file). This implies that MSAFE files must <em>not</em> be
72 * edited to change their time span, otherwise this would break the old
73 * entries overriding mechanism.
74 * </p>
75 * <p>
76 * With these data, the {@link #getInstantFlux(AbsoluteDate)} and {@link
77 * #getMeanFlux(AbsoluteDate)} methods return the same values and the {@link
78 * #get24HoursKp(AbsoluteDate)} and {@link #getThreeHourlyKP(AbsoluteDate)}
79 * methods return the same values.
80 * </p>
81 * <p>
82 * Conversion from Ap index values in the MSAFE file to Kp values used by atmosphere
83 * models is done using Jacchia's equation in [1].
84 * </p>
85 * <p>
86 * With these data, the {@link #getAp(AbsoluteDate date)} method returns an array
87 * of seven times the same daily Ap value, i.e. it is suited to be used only with
88 * the {@link org.orekit.models.earth.atmosphere.NRLMSISE00 NRLMSISE00} atmospheric
89 * model where the switch #9 is set to 1.
90 * </p>
91 *
92 * <h2>References</h2>
93 *
94 * <ol> <li> Jacchia, L. G. "CIRA 1972, recent atmospheric models, and improvements in
95 * progress." COSPAR, 21st Plenary Meeting. Vol. 1. 1978. </li> </ol>
96 *
97 * @author Bruno Revelin
98 * @author Luc Maisonobe
99 * @author Evan Ward
100 * @author Pascal Parraud
101 */
102 public class MarshallSolarActivityFutureEstimation extends AbstractSelfFeedingLoader
103 implements DataLoader, DTM2000InputParameters, NRLMSISE00InputParameters {
104
105 /** Default regular expression for the supported name that work with all officially published files.
106 * @since 10.0
107 */
108 public static final String DEFAULT_SUPPORTED_NAMES =
109 "\\p{Alpha}\\p{Lower}\\p{Lower}\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}(?:f|F)10(?:_prd)?\\.(?:txt|TXT)";
110
111 /** Strength level of activity. */
112 public enum StrengthLevel {
113
114 /** Strong level of activity. */
115 STRONG,
116
117 /** Average level of activity. */
118 AVERAGE,
119
120 /** Weak level of activity. */
121 WEAK
122
123 }
124
125 /** Serializable UID. */
126 private static final long serialVersionUID = -5212198874900835369L;
127
128 /** Pattern for the data fields of MSAFE data. */
129 private final Pattern dataPattern;
130
131 /** Data set. */
132 private final SortedSet<TimeStamped> data;
133
134 /** Selected strength level of activity. */
135 private final StrengthLevel strengthLevel;
136
137 /** UTC time scale. */
138 private final TimeScale utc;
139
140 /** First available date. */
141 private AbsoluteDate firstDate;
142
143 /** Last available date. */
144 private AbsoluteDate lastDate;
145
146 /** Previous set of solar activity parameters. */
147 private LineParameters previousParam;
148
149 /** Current set of solar activity parameters. */
150 private LineParameters currentParam;
151
152 /** Simple constructor. This constructor uses the {@link DataContext#getDefault()
153 * default data context}.
154 * <p>
155 * The original file names used by NASA Marshall space center are of the
156 * form: may2019f10_prd.txt or Oct1999F10.TXT. So a recommended regular
157 * expression for the supported name that work with all published files is:
158 * {@link #DEFAULT_SUPPORTED_NAMES}.
159 * </p>
160 * @param supportedNames regular expression for supported files names
161 * @param strengthLevel selected strength level of activity
162 * @see #MarshallSolarActivityFutureEstimation(String, StrengthLevel, DataProvidersManager, TimeScale)
163 */
164 @DefaultDataContext
165 public MarshallSolarActivityFutureEstimation(final String supportedNames,
166 final StrengthLevel strengthLevel) {
167 this(supportedNames, strengthLevel,
168 DataContext.getDefault().getDataProvidersManager(),
169 DataContext.getDefault().getTimeScales().getUTC());
170 }
171
172 /**
173 * Constructor that allows specifying the source of the MSAFE auxiliary data files.
174 *
175 * @param supportedNames regular expression for supported files names
176 * @param strengthLevel selected strength level of activity
177 * @param dataProvidersManager provides access to auxiliary data files.
178 * @param utc UTC time scale.
179 * @since 10.1
180 */
181 public MarshallSolarActivityFutureEstimation(
182 final String supportedNames,
183 final StrengthLevel strengthLevel,
184 final DataProvidersManager dataProvidersManager,
185 final TimeScale utc) {
186 super(supportedNames, dataProvidersManager);
187
188 firstDate = null;
189 lastDate = null;
190 data = new TreeSet<>(new ChronologicalComparator());
191 this.strengthLevel = strengthLevel;
192 this.utc = utc;
193
194 // the data lines have the following form:
195 // 2010.5003 JUL 83.4 81.3 78.7 6.4 5.9 5.2
196 // 2010.5837 AUG 87.3 83.4 78.5 7.0 6.1 4.9
197 // 2010.6670 SEP 90.8 85.5 79.4 7.8 6.2 4.7
198 // 2010.7503 OCT 94.2 87.6 80.4 9.1 6.4 4.9
199 final StringBuilder builder = new StringBuilder("^");
200
201 // first group: year
202 builder.append("\\p{Blank}*(\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit})");
203
204 // month as fraction of year, not stored in a group
205 builder.append("\\.\\p{Digit}+");
206
207 // second group: month as a three upper case letters abbreviation
208 builder.append("\\p{Blank}+(");
209 for (final Month month : Month.values()) {
210 builder.append(month.getUpperCaseAbbreviation());
211 builder.append('|');
212 }
213 builder.delete(builder.length() - 1, builder.length());
214 builder.append(")");
215
216 // third to eighth group: data fields
217 for (int i = 0; i < 6; ++i) {
218 builder.append("\\p{Blank}+([-+]?[0-9]+\\.[0-9]+)");
219 }
220
221 // end of line
222 builder.append("\\p{Blank}*$");
223
224 // compile the pattern
225 dataPattern = Pattern.compile(builder.toString());
226
227 }
228
229 /** Get the strength level for activity.
230 * @return strength level to set
231 */
232 public StrengthLevel getStrengthLevel() {
233 return strengthLevel;
234 }
235
236 /** Find the data bracketing a specified date.
237 * @param date date to bracket
238 */
239 private void bracketDate(final AbsoluteDate date) {
240
241 if (date.durationFrom(firstDate) < 0) {
242 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE,
243 date, firstDate, lastDate, firstDate.durationFrom(date));
244 }
245 if (date.durationFrom(lastDate) > 0) {
246 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER,
247 date, firstDate, lastDate, date.durationFrom(lastDate));
248 }
249
250 // don't search if the cached selection is fine
251 if (previousParam != null &&
252 date.durationFrom(previousParam.getDate()) > 0 &&
253 date.durationFrom(currentParam.getDate()) <= 0 ) {
254 return;
255 }
256
257 if (date.equals(firstDate)) {
258 currentParam = (LineParameters) data.tailSet(date.shiftedBy(1)).first();
259 previousParam = (LineParameters) data.first();
260 } else if (date.equals(lastDate)) {
261 currentParam = (LineParameters) data.last();
262 previousParam = (LineParameters) data.headSet(date.shiftedBy(-1)).last();
263 } else {
264 currentParam = (LineParameters) data.tailSet(date).first();
265 previousParam = (LineParameters) data.headSet(date).last();
266 }
267
268 }
269
270 @Override
271 public String getSupportedNames() {
272 return super.getSupportedNames();
273 }
274
275 /** {@inheritDoc} */
276 public AbsoluteDate getMinDate() {
277 if (firstDate == null) {
278 feed(this);
279 }
280 return firstDate;
281 }
282
283 /** {@inheritDoc} */
284 public AbsoluteDate getMaxDate() {
285 if (lastDate == null) {
286 feed(this);
287 }
288 return lastDate;
289 }
290
291 /** {@inheritDoc} */
292 public double getInstantFlux(final AbsoluteDate date) {
293 return getMeanFlux(date);
294 }
295
296 /** {@inheritDoc} */
297 public double getMeanFlux(final AbsoluteDate date) {
298
299 // get the neighboring dates
300 bracketDate(date);
301
302 // perform a linear interpolation
303 final AbsoluteDate previousDate = previousParam.getDate();
304 final AbsoluteDate currentDate = currentParam.getDate();
305 final double dt = currentDate.durationFrom(previousDate);
306 final double previousF107 = previousParam.getF107();
307 final double currentF107 = currentParam.getF107();
308 final double previousWeight = currentDate.durationFrom(date) / dt;
309 final double currentWeight = date.durationFrom(previousDate) / dt;
310
311 return previousF107 * previousWeight + currentF107 * currentWeight;
312
313 }
314
315 /** {@inheritDoc} */
316 public double getThreeHourlyKP(final AbsoluteDate date) {
317 return get24HoursKp(date);
318 }
319
320 /** Get the date of the file from which data at the specified date comes from.
321 * <p>
322 * If several MSAFE files are available, some dates may appear in several
323 * files (for example August 2007 is in all files from the first one
324 * published in March 1999 to the February 2008 file). In this case, the
325 * data from the most recent file is used and the older ones are discarded.
326 * The date of the file is assumed to be 6 months after its first entry
327 * (which explains why the file having August 2007 as its first entry is the
328 * February 2008 file). This implies that MSAFE files must <em>not</em> be
329 * edited to change their time span, otherwise this would break the old
330 * entries overriding mechanism.
331 * </p>
332 * @param date date of the solar activity data
333 * @return date of the file
334 */
335 public DateComponents getFileDate(final AbsoluteDate date) {
336 bracketDate(date);
337 final double dtP = date.durationFrom(previousParam.getDate());
338 final double dtC = currentParam.getDate().durationFrom(date);
339 return (dtP < dtC) ? previousParam.getFileDate() : currentParam.getFileDate();
340 }
341
342 /** The Kp index is derived from the Ap index.
343 * <p>The method used is explained on <a
344 * href="http://www.ngdc.noaa.gov/stp/GEOMAG/kp_ap.html">
345 * NOAA website.</a> as follows:</p>
346 * <p>The scale is 0 to 9 expressed in thirds of a unit, e.g. 5- is 4 2/3,
347 * 5 is 5 and 5+ is 5 1/3. The ap (equivalent range) index is derived from
348 * the Kp index as follows:</p>
349 * <table border="1">
350 * <caption>Kp / Ap Conversion Table</caption>
351 * <tbody>
352 * <tr>
353 * <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>
354 * </tr>
355 * <tr>
356 * <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>
357 * </tr>
358 * <tr>
359 * <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>
360 * </tr>
361 * <tr>
362 * <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>
363 * </tr>
364 * </tbody>
365 * </table>
366 * @param date date of the Kp data
367 * @return the 24H geomagnetic index
368 */
369 public double get24HoursKp(final AbsoluteDate date) {
370
371 // get the daily Ap
372 final double ap = getDailyAp(date);
373
374 // get the corresponding Kp index from
375 // equation 4 in [1] for Ap to Kp conversion
376 return 1.89 * FastMath.asinh(0.154 * ap);
377 }
378
379 /** {@inheritDoc} */
380 public double getDailyFlux(final AbsoluteDate date) {
381 return getMeanFlux(date.shiftedBy(-Constants.JULIAN_DAY));
382 }
383
384 /** {@inheritDoc} */
385 public double getAverageFlux(final AbsoluteDate date) {
386
387 // Initializes the average flux
388 double average = 0.;
389
390 // Loops over the 81 days centered on the given date
391 for (int i = -40; i < 41; i++) {
392 average += getMeanFlux(date.shiftedBy(i * Constants.JULIAN_DAY));
393 }
394
395 // Returns the 81 day average flux
396 return average / 81;
397 }
398
399 /** {@inheritDoc} */
400 public double[] getAp(final AbsoluteDate date) {
401
402 // Gets the AP for the current date
403 final double ap = getDailyAp(date);
404
405 // Retuns an array of Ap filled in with the daily Ap only
406 return new double[] {ap, ap, ap, ap, ap, ap, ap};
407 }
408
409 /** Gets the daily Ap index.
410 *
411 * @param date the current date
412 * @return the daily Ap index
413 */
414 private double getDailyAp(final AbsoluteDate date) {
415
416 // get the neighboring dates
417 bracketDate(date);
418
419 // perform a linear interpolation
420 final AbsoluteDate previousDate = previousParam.getDate();
421 final AbsoluteDate currentDate = currentParam.getDate();
422 final double dt = currentDate.durationFrom(previousDate);
423 final double previousAp = previousParam.getAp();
424 final double currentAp = currentParam.getAp();
425 final double previousWeight = currentDate.durationFrom(date) / dt;
426 final double currentWeight = date.durationFrom(previousDate) / dt;
427
428 // returns the daily Ap interpolated at the date
429 return previousAp * previousWeight + currentAp * currentWeight;
430 }
431
432 /** Container class for Solar activity indexes. */
433 private static class LineParameters implements TimeStamped, Serializable {
434
435 /** Serializable UID. */
436 private static final long serialVersionUID = 6607862001953526475L;
437
438 /** File date. */
439 private final DateComponents fileDate;
440
441 /** Entry date. */
442 private final AbsoluteDate date;
443
444 /** F10.7 flux at date. */
445 private final double f107;
446
447 /** Ap index at date. */
448 private final double ap;
449
450 /** Simple constructor.
451 * @param fileDate file date
452 * @param date entry date
453 * @param f107 F10.7 flux at date
454 * @param ap Ap index at date
455 */
456 private LineParameters(final DateComponents fileDate, final AbsoluteDate date, final double f107, final double ap) {
457 this.fileDate = fileDate;
458 this.date = date;
459 this.f107 = f107;
460 this.ap = ap;
461 }
462
463 /** Get the file date.
464 * @return file date
465 */
466 public DateComponents getFileDate() {
467 return fileDate;
468 }
469
470 /** Get the current date.
471 * @return current date
472 */
473 public AbsoluteDate getDate() {
474 return date;
475 }
476
477 /** Get the F10.0 flux.
478 * @return f10.7 flux
479 */
480 public double getF107() {
481 return f107;
482 }
483
484 /** Get the Ap index.
485 * @return Ap index
486 */
487 public double getAp() {
488 return ap;
489 }
490
491 }
492
493 /** {@inheritDoc} */
494 public void loadData(final InputStream input, final String name)
495 throws IOException, ParseException, OrekitException {
496
497 // select the groups we want to store
498 final int f107Group;
499 final int apGroup;
500 switch (strengthLevel) {
501 case STRONG :
502 f107Group = 3;
503 apGroup = 6;
504 break;
505 case AVERAGE :
506 f107Group = 4;
507 apGroup = 7;
508 break;
509 default :
510 f107Group = 5;
511 apGroup = 8;
512 break;
513 }
514
515 boolean inData = false;
516 DateComponents fileDate = null;
517 // read the data
518 try (BufferedReader reader = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
519
520 for (String line = reader.readLine(); line != null; line = reader.readLine()) {
521 line = line.trim();
522 if (line.length() > 0) {
523 final Matcher matcher = dataPattern.matcher(line);
524 if (matcher.matches()) {
525
526 // we are in the data section
527 inData = true;
528
529 // extract the data from the line
530 final int year = Integer.parseInt(matcher.group(1));
531 final Month month = Month.parseMonth(matcher.group(2));
532 final AbsoluteDate date = new AbsoluteDate(year, month, 1, this.utc);
533 if (fileDate == null) {
534 // the first entry of each file correspond exactly to 6 months before file publication
535 // so we compute the file date by adding 6 months to its first entry
536 if (month.getNumber() > 6) {
537 fileDate = new DateComponents(year + 1, month.getNumber() - 6, 1);
538 } else {
539 fileDate = new DateComponents(year, month.getNumber() + 6, 1);
540 }
541 }
542
543 // check if there is already an entry for this date or not
544 boolean addEntry = false;
545 final Iterator<TimeStamped> iterator = data.tailSet(date).iterator();
546 if (iterator.hasNext()) {
547 final LineParameters existingEntry = (LineParameters) iterator.next();
548 if (existingEntry.getDate().equals(date)) {
549 // there is an entry for this date
550 if (existingEntry.getFileDate().compareTo(fileDate) < 0) {
551 // the entry was read from an earlier file
552 // we replace it with the new entry as it is fresher
553 iterator.remove();
554 addEntry = true;
555 }
556 } else {
557 // it is the first entry we get for this date
558 addEntry = true;
559 }
560 } else {
561 // it is the first entry we get for this date
562 addEntry = true;
563 }
564 if (addEntry) {
565 // we must add the new entry
566 data.add(new LineParameters(fileDate, date,
567 Double.parseDouble(matcher.group(f107Group)),
568 Double.parseDouble(matcher.group(apGroup))));
569 }
570
571 } else {
572 if (inData) {
573 // we have already read some data, so we are not in the header anymore
574 // however, we don't recognize this non-empty line,
575 // we consider the file is corrupted
576 throw new OrekitException(OrekitMessages.NOT_A_MARSHALL_SOLAR_ACTIVITY_FUTURE_ESTIMATION_FILE,
577 name);
578 }
579 }
580 }
581 }
582
583 }
584
585 if (data.isEmpty()) {
586 throw new OrekitException(OrekitMessages.NOT_A_MARSHALL_SOLAR_ACTIVITY_FUTURE_ESTIMATION_FILE,
587 name);
588 }
589 firstDate = data.first().getDate();
590 lastDate = data.last().getDate();
591
592 }
593
594 /** {@inheritDoc} */
595 public boolean stillAcceptsData() {
596 return true;
597 }
598
599 /** Replace the instance with a data transfer object for serialization.
600 * @return data transfer object that will be serialized
601 */
602 @DefaultDataContext
603 private Object writeReplace() {
604 return new DataTransferObject(getSupportedNames(), strengthLevel);
605 }
606
607 /** Internal class used only for serialization. */
608 @DefaultDataContext
609 private static class DataTransferObject implements Serializable {
610
611 /** Serializable UID. */
612 private static final long serialVersionUID = -5212198874900835369L;
613
614 /** Regular expression that matches the names of the IONEX files. */
615 private final String supportedNames;
616
617 /** Selected strength level of activity. */
618 private final StrengthLevel strengthLevel;
619
620 /** Simple constructor.
621 * @param supportedNames regular expression for supported files names
622 * @param strengthLevel selected strength level of activity
623 */
624 DataTransferObject(final String supportedNames,
625 final StrengthLevel strengthLevel) {
626 this.supportedNames = supportedNames;
627 this.strengthLevel = strengthLevel;
628 }
629
630 /** Replace the deserialized data transfer object with a {@link MarshallSolarActivityFutureEstimation}.
631 * @return replacement {@link MarshallSolarActivityFutureEstimation}
632 */
633 private Object readResolve() {
634 try {
635 return new MarshallSolarActivityFutureEstimation(supportedNames, strengthLevel);
636 } catch (OrekitException oe) {
637 throw new OrekitInternalError(oe);
638 }
639 }
640
641 }
642
643 }