MarshallSolarActivityFutureEstimation.java
/* Copyright 2002-2021 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.models.earth.atmosphere.data;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.io.Serializable;
import java.nio.charset.StandardCharsets;
import java.text.ParseException;
import java.util.Iterator;
import java.util.SortedSet;
import java.util.TreeSet;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import org.hipparchus.util.FastMath;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.data.AbstractSelfFeedingLoader;
import org.orekit.data.DataContext;
import org.orekit.data.DataLoader;
import org.orekit.data.DataProvidersManager;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitInternalError;
import org.orekit.errors.OrekitMessages;
import org.orekit.models.earth.atmosphere.DTM2000InputParameters;
import org.orekit.models.earth.atmosphere.NRLMSISE00InputParameters;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.ChronologicalComparator;
import org.orekit.time.DateComponents;
import org.orekit.time.Month;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeStamped;
import org.orekit.utils.Constants;
/**
* This class reads and provides solar activity data needed by
* atmospheric models: F107 solar flux, Ap and Kp indexes.
* <p>
* The data are retrieved through the NASA Marshall
* Solar Activity Future Estimation (MSAFE) as estimates of monthly
* F10.7 Mean solar flux and Ap geomagnetic parameter.
* The data can be retrieved at the NASA <a
* href="http://sail.msfc.nasa.gov/archive_index.htm">
* Marshall Solar Activity website</a>.
* Here Kp indices are deduced from Ap indexes, which in turn are tabulated
* equivalent of retrieved Ap values.
* </p>
* <p>
* If several MSAFE files are available, some dates may appear in several
* files (for example August 2007 is in all files from the first one
* published in March 1999 to the February 2008 file). In this case, the
* data from the most recent file is used and the older ones are discarded.
* The date of the file is assumed to be 6 months after its first entry
* (which explains why the file having August 2007 as its first entry is the
* February 2008 file). This implies that MSAFE files must <em>not</em> be
* edited to change their time span, otherwise this would break the old
* entries overriding mechanism.
* </p>
* <p>
* With these data, the {@link #getInstantFlux(AbsoluteDate)} and {@link
* #getMeanFlux(AbsoluteDate)} methods return the same values and the {@link
* #get24HoursKp(AbsoluteDate)} and {@link #getThreeHourlyKP(AbsoluteDate)}
* methods return the same values.
* </p>
* <p>
* Conversion from Ap index values in the MSAFE file to Kp values used by atmosphere
* models is done using Jacchia's equation in [1].
* </p>
* <p>
* With these data, the {@link #getAp(AbsoluteDate date)} method returns an array
* of seven times the same daily Ap value, i.e. it is suited to be used only with
* the {@link org.orekit.models.earth.atmosphere.NRLMSISE00 NRLMSISE00} atmospheric
* model where the switch #9 is set to 1.
* </p>
*
* <h2>References</h2>
*
* <ol> <li> Jacchia, L. G. "CIRA 1972, recent atmospheric models, and improvements in
* progress." COSPAR, 21st Plenary Meeting. Vol. 1. 1978. </li> </ol>
*
* @author Bruno Revelin
* @author Luc Maisonobe
* @author Evan Ward
* @author Pascal Parraud
*/
public class MarshallSolarActivityFutureEstimation extends AbstractSelfFeedingLoader
implements DataLoader, DTM2000InputParameters, NRLMSISE00InputParameters {
/** Default regular expression for the supported name that work with all officially published files.
* @since 10.0
*/
public static final String DEFAULT_SUPPORTED_NAMES =
"\\p{Alpha}\\p{Lower}\\p{Lower}\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}(?:f|F)10(?:_prd)?\\.(?:txt|TXT)";
/** Strength level of activity. */
public enum StrengthLevel {
/** Strong level of activity. */
STRONG,
/** Average level of activity. */
AVERAGE,
/** Weak level of activity. */
WEAK
}
/** Serializable UID. */
private static final long serialVersionUID = -5212198874900835369L;
/** Pattern for the data fields of MSAFE data. */
private final Pattern dataPattern;
/** Data set. */
private final SortedSet<TimeStamped> data;
/** Selected strength level of activity. */
private final StrengthLevel strengthLevel;
/** UTC time scale. */
private final TimeScale utc;
/** First available date. */
private AbsoluteDate firstDate;
/** Last available date. */
private AbsoluteDate lastDate;
/** Previous set of solar activity parameters. */
private LineParameters previousParam;
/** Current set of solar activity parameters. */
private LineParameters currentParam;
/** Simple constructor. This constructor uses the {@link DataContext#getDefault()
* default data context}.
* <p>
* The original file names used by NASA Marshall space center are of the
* form: may2019f10_prd.txt or Oct1999F10.TXT. So a recommended regular
* expression for the supported name that work with all published files is:
* {@link #DEFAULT_SUPPORTED_NAMES}.
* </p>
* @param supportedNames regular expression for supported files names
* @param strengthLevel selected strength level of activity
* @see #MarshallSolarActivityFutureEstimation(String, StrengthLevel, DataProvidersManager, TimeScale)
*/
@DefaultDataContext
public MarshallSolarActivityFutureEstimation(final String supportedNames,
final StrengthLevel strengthLevel) {
this(supportedNames, strengthLevel,
DataContext.getDefault().getDataProvidersManager(),
DataContext.getDefault().getTimeScales().getUTC());
}
/**
* Constructor that allows specifying the source of the MSAFE auxiliary data files.
*
* @param supportedNames regular expression for supported files names
* @param strengthLevel selected strength level of activity
* @param dataProvidersManager provides access to auxiliary data files.
* @param utc UTC time scale.
* @since 10.1
*/
public MarshallSolarActivityFutureEstimation(
final String supportedNames,
final StrengthLevel strengthLevel,
final DataProvidersManager dataProvidersManager,
final TimeScale utc) {
super(supportedNames, dataProvidersManager);
firstDate = null;
lastDate = null;
data = new TreeSet<>(new ChronologicalComparator());
this.strengthLevel = strengthLevel;
this.utc = utc;
// the data lines have the following form:
// 2010.5003 JUL 83.4 81.3 78.7 6.4 5.9 5.2
// 2010.5837 AUG 87.3 83.4 78.5 7.0 6.1 4.9
// 2010.6670 SEP 90.8 85.5 79.4 7.8 6.2 4.7
// 2010.7503 OCT 94.2 87.6 80.4 9.1 6.4 4.9
final StringBuilder builder = new StringBuilder("^");
// first group: year
builder.append("\\p{Blank}*(\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit})");
// month as fraction of year, not stored in a group
builder.append("\\.\\p{Digit}+");
// second group: month as a three upper case letters abbreviation
builder.append("\\p{Blank}+(");
for (final Month month : Month.values()) {
builder.append(month.getUpperCaseAbbreviation());
builder.append('|');
}
builder.delete(builder.length() - 1, builder.length());
builder.append(")");
// third to eighth group: data fields
for (int i = 0; i < 6; ++i) {
builder.append("\\p{Blank}+([-+]?[0-9]+\\.[0-9]+)");
}
// end of line
builder.append("\\p{Blank}*$");
// compile the pattern
dataPattern = Pattern.compile(builder.toString());
}
/** Get the strength level for activity.
* @return strength level to set
*/
public StrengthLevel getStrengthLevel() {
return strengthLevel;
}
/** Find the data bracketing a specified date.
* @param date date to bracket
*/
private void bracketDate(final AbsoluteDate date) {
if (date.durationFrom(firstDate) < 0) {
throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE,
date, firstDate, lastDate, firstDate.durationFrom(date));
}
if (date.durationFrom(lastDate) > 0) {
throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER,
date, firstDate, lastDate, date.durationFrom(lastDate));
}
// don't search if the cached selection is fine
if (previousParam != null &&
date.durationFrom(previousParam.getDate()) > 0 &&
date.durationFrom(currentParam.getDate()) <= 0 ) {
return;
}
if (date.equals(firstDate)) {
currentParam = (LineParameters) data.tailSet(date.shiftedBy(1)).first();
previousParam = (LineParameters) data.first();
} else if (date.equals(lastDate)) {
currentParam = (LineParameters) data.last();
previousParam = (LineParameters) data.headSet(date.shiftedBy(-1)).last();
} else {
currentParam = (LineParameters) data.tailSet(date).first();
previousParam = (LineParameters) data.headSet(date).last();
}
}
@Override
public String getSupportedNames() {
return super.getSupportedNames();
}
/** {@inheritDoc} */
public AbsoluteDate getMinDate() {
if (firstDate == null) {
feed(this);
}
return firstDate;
}
/** {@inheritDoc} */
public AbsoluteDate getMaxDate() {
if (lastDate == null) {
feed(this);
}
return lastDate;
}
/** {@inheritDoc} */
public double getInstantFlux(final AbsoluteDate date) {
return getMeanFlux(date);
}
/** {@inheritDoc} */
public double getMeanFlux(final AbsoluteDate date) {
// get the neighboring dates
bracketDate(date);
// perform a linear interpolation
final AbsoluteDate previousDate = previousParam.getDate();
final AbsoluteDate currentDate = currentParam.getDate();
final double dt = currentDate.durationFrom(previousDate);
final double previousF107 = previousParam.getF107();
final double currentF107 = currentParam.getF107();
final double previousWeight = currentDate.durationFrom(date) / dt;
final double currentWeight = date.durationFrom(previousDate) / dt;
return previousF107 * previousWeight + currentF107 * currentWeight;
}
/** {@inheritDoc} */
public double getThreeHourlyKP(final AbsoluteDate date) {
return get24HoursKp(date);
}
/** Get the date of the file from which data at the specified date comes from.
* <p>
* If several MSAFE files are available, some dates may appear in several
* files (for example August 2007 is in all files from the first one
* published in March 1999 to the February 2008 file). In this case, the
* data from the most recent file is used and the older ones are discarded.
* The date of the file is assumed to be 6 months after its first entry
* (which explains why the file having August 2007 as its first entry is the
* February 2008 file). This implies that MSAFE files must <em>not</em> be
* edited to change their time span, otherwise this would break the old
* entries overriding mechanism.
* </p>
* @param date date of the solar activity data
* @return date of the file
*/
public DateComponents getFileDate(final AbsoluteDate date) {
bracketDate(date);
final double dtP = date.durationFrom(previousParam.getDate());
final double dtC = currentParam.getDate().durationFrom(date);
return (dtP < dtC) ? previousParam.getFileDate() : currentParam.getFileDate();
}
/** The Kp index is derived from the Ap index.
* <p>The method used is explained on <a
* href="http://www.ngdc.noaa.gov/stp/GEOMAG/kp_ap.html">
* NOAA website.</a> as follows:</p>
* <p>The scale is 0 to 9 expressed in thirds of a unit, e.g. 5- is 4 2/3,
* 5 is 5 and 5+ is 5 1/3. The ap (equivalent range) index is derived from
* the Kp index as follows:</p>
* <table border="1">
* <caption>Kp / Ap Conversion Table</caption>
* <tbody>
* <tr>
* <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>
* </tr>
* <tr>
* <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>
* </tr>
* <tr>
* <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>
* </tr>
* <tr>
* <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>
* </tr>
* </tbody>
* </table>
* @param date date of the Kp data
* @return the 24H geomagnetic index
*/
public double get24HoursKp(final AbsoluteDate date) {
// get the daily Ap
final double ap = getDailyAp(date);
// get the corresponding Kp index from
// equation 4 in [1] for Ap to Kp conversion
return 1.89 * FastMath.asinh(0.154 * ap);
}
/** {@inheritDoc} */
public double getDailyFlux(final AbsoluteDate date) {
return getMeanFlux(date.shiftedBy(-Constants.JULIAN_DAY));
}
/** {@inheritDoc} */
public double getAverageFlux(final AbsoluteDate date) {
// Initializes the average flux
double average = 0.;
// Loops over the 81 days centered on the given date
for (int i = -40; i < 41; i++) {
average += getMeanFlux(date.shiftedBy(i * Constants.JULIAN_DAY));
}
// Returns the 81 day average flux
return average / 81;
}
/** {@inheritDoc} */
public double[] getAp(final AbsoluteDate date) {
// Gets the AP for the current date
final double ap = getDailyAp(date);
// Retuns an array of Ap filled in with the daily Ap only
return new double[] {ap, ap, ap, ap, ap, ap, ap};
}
/** Gets the daily Ap index.
*
* @param date the current date
* @return the daily Ap index
*/
private double getDailyAp(final AbsoluteDate date) {
// get the neighboring dates
bracketDate(date);
// perform a linear interpolation
final AbsoluteDate previousDate = previousParam.getDate();
final AbsoluteDate currentDate = currentParam.getDate();
final double dt = currentDate.durationFrom(previousDate);
final double previousAp = previousParam.getAp();
final double currentAp = currentParam.getAp();
final double previousWeight = currentDate.durationFrom(date) / dt;
final double currentWeight = date.durationFrom(previousDate) / dt;
// returns the daily Ap interpolated at the date
return previousAp * previousWeight + currentAp * currentWeight;
}
/** Container class for Solar activity indexes. */
private static class LineParameters implements TimeStamped, Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 6607862001953526475L;
/** File date. */
private final DateComponents fileDate;
/** Entry date. */
private final AbsoluteDate date;
/** F10.7 flux at date. */
private final double f107;
/** Ap index at date. */
private final double ap;
/** Simple constructor.
* @param fileDate file date
* @param date entry date
* @param f107 F10.7 flux at date
* @param ap Ap index at date
*/
private LineParameters(final DateComponents fileDate, final AbsoluteDate date, final double f107, final double ap) {
this.fileDate = fileDate;
this.date = date;
this.f107 = f107;
this.ap = ap;
}
/** Get the file date.
* @return file date
*/
public DateComponents getFileDate() {
return fileDate;
}
/** Get the current date.
* @return current date
*/
public AbsoluteDate getDate() {
return date;
}
/** Get the F10.0 flux.
* @return f10.7 flux
*/
public double getF107() {
return f107;
}
/** Get the Ap index.
* @return Ap index
*/
public double getAp() {
return ap;
}
}
/** {@inheritDoc} */
public void loadData(final InputStream input, final String name)
throws IOException, ParseException, OrekitException {
// select the groups we want to store
final int f107Group;
final int apGroup;
switch (strengthLevel) {
case STRONG :
f107Group = 3;
apGroup = 6;
break;
case AVERAGE :
f107Group = 4;
apGroup = 7;
break;
default :
f107Group = 5;
apGroup = 8;
break;
}
boolean inData = false;
DateComponents fileDate = null;
// read the data
try (BufferedReader reader = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
for (String line = reader.readLine(); line != null; line = reader.readLine()) {
line = line.trim();
if (line.length() > 0) {
final Matcher matcher = dataPattern.matcher(line);
if (matcher.matches()) {
// we are in the data section
inData = true;
// extract the data from the line
final int year = Integer.parseInt(matcher.group(1));
final Month month = Month.parseMonth(matcher.group(2));
final AbsoluteDate date = new AbsoluteDate(year, month, 1, this.utc);
if (fileDate == null) {
// the first entry of each file correspond exactly to 6 months before file publication
// so we compute the file date by adding 6 months to its first entry
if (month.getNumber() > 6) {
fileDate = new DateComponents(year + 1, month.getNumber() - 6, 1);
} else {
fileDate = new DateComponents(year, month.getNumber() + 6, 1);
}
}
// check if there is already an entry for this date or not
boolean addEntry = false;
final Iterator<TimeStamped> iterator = data.tailSet(date).iterator();
if (iterator.hasNext()) {
final LineParameters existingEntry = (LineParameters) iterator.next();
if (existingEntry.getDate().equals(date)) {
// there is an entry for this date
if (existingEntry.getFileDate().compareTo(fileDate) < 0) {
// the entry was read from an earlier file
// we replace it with the new entry as it is fresher
iterator.remove();
addEntry = true;
}
} else {
// it is the first entry we get for this date
addEntry = true;
}
} else {
// it is the first entry we get for this date
addEntry = true;
}
if (addEntry) {
// we must add the new entry
data.add(new LineParameters(fileDate, date,
Double.parseDouble(matcher.group(f107Group)),
Double.parseDouble(matcher.group(apGroup))));
}
} else {
if (inData) {
// we have already read some data, so we are not in the header anymore
// however, we don't recognize this non-empty line,
// we consider the file is corrupted
throw new OrekitException(OrekitMessages.NOT_A_MARSHALL_SOLAR_ACTIVITY_FUTURE_ESTIMATION_FILE,
name);
}
}
}
}
}
if (data.isEmpty()) {
throw new OrekitException(OrekitMessages.NOT_A_MARSHALL_SOLAR_ACTIVITY_FUTURE_ESTIMATION_FILE,
name);
}
firstDate = data.first().getDate();
lastDate = data.last().getDate();
}
/** {@inheritDoc} */
public boolean stillAcceptsData() {
return true;
}
/** Replace the instance with a data transfer object for serialization.
* @return data transfer object that will be serialized
*/
@DefaultDataContext
private Object writeReplace() {
return new DataTransferObject(getSupportedNames(), strengthLevel);
}
/** Internal class used only for serialization. */
@DefaultDataContext
private static class DataTransferObject implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = -5212198874900835369L;
/** Regular expression that matches the names of the IONEX files. */
private final String supportedNames;
/** Selected strength level of activity. */
private final StrengthLevel strengthLevel;
/** Simple constructor.
* @param supportedNames regular expression for supported files names
* @param strengthLevel selected strength level of activity
*/
DataTransferObject(final String supportedNames,
final StrengthLevel strengthLevel) {
this.supportedNames = supportedNames;
this.strengthLevel = strengthLevel;
}
/** Replace the deserialized data transfer object with a {@link MarshallSolarActivityFutureEstimation}.
* @return replacement {@link MarshallSolarActivityFutureEstimation}
*/
private Object readResolve() {
try {
return new MarshallSolarActivityFutureEstimation(supportedNames, strengthLevel);
} catch (OrekitException oe) {
throw new OrekitInternalError(oe);
}
}
}
}