MarshallSolarActivityFutureEstimation.java

  1. /* Copyright 2002-2016 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.forces.drag;

  18. import java.io.BufferedReader;
  19. import java.io.IOException;
  20. import java.io.InputStream;
  21. import java.io.InputStreamReader;
  22. import java.io.Serializable;
  23. import java.text.ParseException;
  24. import java.util.Arrays;
  25. import java.util.Iterator;
  26. import java.util.SortedSet;
  27. import java.util.TreeSet;
  28. import java.util.regex.Matcher;
  29. import java.util.regex.Pattern;

  30. import org.orekit.data.DataLoader;
  31. import org.orekit.data.DataProvidersManager;
  32. import org.orekit.errors.OrekitException;
  33. import org.orekit.errors.OrekitMessages;
  34. import org.orekit.time.AbsoluteDate;
  35. import org.orekit.time.ChronologicalComparator;
  36. import org.orekit.time.DateComponents;
  37. import org.orekit.time.Month;
  38. import org.orekit.time.TimeScale;
  39. import org.orekit.time.TimeScalesFactory;
  40. import org.orekit.time.TimeStamped;

  41. /**
  42.  * This class reads and provides solar activity data needed by
  43.  * atmospheric models: F107 solar flux and Kp indexes.
  44.  * <p>
  45.  * The data are retrieved through the NASA Marshall
  46.  * Solar Activity Future Estimation (MSAFE) as estimates of monthly
  47.  * F10.7 Mean solar flux and Ap geomagnetic parameter.
  48.  * The data can be retrieved at the NASA <a
  49.  * href="http://sail.msfc.nasa.gov/archive_index.htm">
  50.  * Marshall Solar Activity website</a>.
  51.  * Here Kp indices are deduced from Ap indexes, which in turn are tabulated
  52.  * equivalent of retrieved Ap values.
  53.  * </p>
  54.  * <p>
  55.  * If several MSAFE files are available, some dates may appear in several
  56.  * files (for example August 2007 is in all files from the first one
  57.  * published in March 1999 to the February 2008 file). In this case, the
  58.  * data from the most recent file is used and the older ones are discarded.
  59.  * The date of the file is assumed to be 6 months after its first entry
  60.  * (which explains why the file having August 2007 as its first entry is the
  61.  * February 2008 file). This implies that MSAFE files must <em>not</em> be
  62.  * edited to change their time span, otherwise this would break the old
  63.  * entries overriding mechanism.
  64.  * </p>
  65.  * <p>
  66.  * With these data, the {@link #getInstantFlux(AbsoluteDate)} and {@link
  67.  * #getMeanFlux(AbsoluteDate)} methods return the same values and the {@link
  68.  * #get24HoursKp(AbsoluteDate)} and {@link #getThreeHourlyKP(AbsoluteDate)}
  69.  * methods return the same values.
  70.  * </p>
  71.  * @author Bruno Revelin
  72.  * @author Luc Maisonobe
  73.  */
  74. public class MarshallSolarActivityFutureEstimation implements DTM2000InputParameters, DataLoader {

  75.     /** Strength level of activity. */
  76.     public enum StrengthLevel {

  77.         /** Strong level of activity. */
  78.         STRONG,

  79.         /** Average level of activity. */
  80.         AVERAGE,

  81.         /** Weak level of activity. */
  82.         WEAK

  83.     };

  84.     /** Serializable UID. */
  85.     private static final long serialVersionUID = -5212198874900835369L;

  86.     /** 1/3. */
  87.     private static final double ONE_THIRD = 1.0 / 3.0;

  88.     /** 3 hours geomagnetic indices array. */
  89.     private static final double[] KP_ARRAY = new double[] {
  90.         0, 0 + ONE_THIRD, 1 - ONE_THIRD, 1, 1 + ONE_THIRD, 2 - ONE_THIRD,
  91.         2, 2 + ONE_THIRD, 3 - ONE_THIRD, 3, 3 + ONE_THIRD, 4 - ONE_THIRD,
  92.         4, 4 + ONE_THIRD, 5 - ONE_THIRD, 5, 5 + ONE_THIRD, 6 - ONE_THIRD,
  93.         6, 6 + ONE_THIRD, 7 - ONE_THIRD, 7, 7 + ONE_THIRD, 8 - ONE_THIRD,
  94.         8, 8 + ONE_THIRD, 9 - ONE_THIRD, 9
  95.     };

  96.     /** Mean geomagnetic indices array. */
  97.     private static final double[] AP_ARRAY = new double[] {
  98.         0, 2, 3, 4, 5, 6, 7, 9, 12, 15, 18, 22, 27, 32,
  99.         39, 48, 56, 67, 80, 94, 111, 132, 154, 179, 207, 236, 300, 400
  100.     };

  101.     /** Pattern for the data fields of MSAFE data. */
  102.     private final Pattern dataPattern;

  103.     /** Data set. */
  104.     private final SortedSet<TimeStamped> data;

  105.     /** Selected strength level of activity. */
  106.     private final StrengthLevel strengthLevel;

  107.     /** First available date. */
  108.     private AbsoluteDate firstDate;

  109.     /** Last available date. */
  110.     private AbsoluteDate lastDate;

  111.     /** Previous set of solar activity parameters. */
  112.     private LineParameters previousParam;

  113.     /** Current set of solar activity parameters. */
  114.     private LineParameters currentParam;

  115.     /** Regular expression for supported files names. */
  116.     private final String supportedNames;

  117.     /** Simple constructor.
  118.      * <p>
  119.      * The original file names used by NASA Marshall space center are of the
  120.      * form: Dec2010F10.txt or Oct1999F10.TXT. So a recommended regular
  121.      * expression for the supported name that work with all published files is:
  122.      * "(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)"
  123.      * </p>
  124.      * @param supportedNames regular expression for supported files names
  125.      * @param strengthLevel selected strength level of activity
  126.      */
  127.     public MarshallSolarActivityFutureEstimation(final String supportedNames,
  128.                                                  final StrengthLevel strengthLevel) {

  129.         firstDate           = null;
  130.         lastDate            = null;
  131.         data                = new TreeSet<TimeStamped>(new ChronologicalComparator());
  132.         this.supportedNames = supportedNames;
  133.         this.strengthLevel  = strengthLevel;

  134.         // the data lines have the following form:
  135.         // 2010.5003   JUL    83.4      81.3      78.7       6.4       5.9       5.2
  136.         // 2010.5837   AUG    87.3      83.4      78.5       7.0       6.1       4.9
  137.         // 2010.6670   SEP    90.8      85.5      79.4       7.8       6.2       4.7
  138.         // 2010.7503   OCT    94.2      87.6      80.4       9.1       6.4       4.9
  139.         final StringBuilder builder = new StringBuilder("^");

  140.         // first group: year
  141.         builder.append("\\p{Blank}*(\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit})");

  142.         // month as fraction of year, not stored in a group
  143.         builder.append("\\.\\p{Digit}+");

  144.         // second group: month as a three upper case letters abbreviation
  145.         builder.append("\\p{Blank}+(");
  146.         for (final Month month : Month.values()) {
  147.             builder.append(month.getUpperCaseAbbreviation());
  148.             builder.append('|');
  149.         }
  150.         builder.delete(builder.length() - 1, builder.length());
  151.         builder.append(")");

  152.         // third to eighth group: data fields
  153.         for (int i = 0; i < 6; ++i) {
  154.             builder.append("\\p{Blank}+([-+]?[0-9]+\\.[0-9]+)");
  155.         }

  156.         // end of line
  157.         builder.append("\\p{Blank}*$");

  158.         // compile the pattern
  159.         dataPattern = Pattern.compile(builder.toString());

  160.     }

  161.     /** Get the strength level for activity.
  162.      * @return strength level to set
  163.      */
  164.     public StrengthLevel getStrengthLevel() {
  165.         return strengthLevel;
  166.     }

  167.     /** Find the data bracketing a specified date.
  168.      * @param date date to bracket
  169.      * @throws OrekitException if specified date is out of range
  170.      */
  171.     private void bracketDate(final AbsoluteDate date) throws OrekitException {

  172.         if ((date.durationFrom(firstDate) < 0) || (date.durationFrom(lastDate) > 0)) {
  173.             throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE,
  174.                                       date, firstDate, lastDate);
  175.         }

  176.         // don't search if the cached selection is fine
  177.         if ((previousParam != null) &&
  178.             (date.durationFrom(previousParam.getDate()) > 0) &&
  179.             (date.durationFrom(currentParam.getDate()) <= 0 )) {
  180.             return;
  181.         }

  182.         if (date.equals(firstDate)) {
  183.             currentParam  = (LineParameters) data.tailSet(date.shiftedBy(1)).first();
  184.             previousParam = (LineParameters) data.first();
  185.         } else if (date.equals(lastDate)) {
  186.             currentParam  = (LineParameters) data.last();
  187.             previousParam = (LineParameters) data.headSet(date.shiftedBy(-1)).last();
  188.         } else {
  189.             currentParam  = (LineParameters) data.tailSet(date).first();
  190.             previousParam = (LineParameters) data.headSet(date).last();
  191.         }

  192.     }

  193.     /** Get the supported names for data files.
  194.      * @return regular expression for the supported names for data files
  195.      */
  196.     public String getSupportedNames() {
  197.         return supportedNames;
  198.     }

  199.     /** {@inheritDoc} */
  200.     public AbsoluteDate getMinDate() throws OrekitException {
  201.         if (firstDate == null) {
  202.             DataProvidersManager.getInstance().feed(getSupportedNames(), this);
  203.         }
  204.         return firstDate;
  205.     }

  206.     /** {@inheritDoc} */
  207.     public AbsoluteDate getMaxDate() throws OrekitException {
  208.         if (firstDate == null) {
  209.             DataProvidersManager.getInstance().feed(getSupportedNames(), this);
  210.         }
  211.         return lastDate;
  212.     }

  213.     /** {@inheritDoc} */
  214.     public double getInstantFlux(final AbsoluteDate date) throws OrekitException {
  215.         return getMeanFlux(date);
  216.     }

  217.     /** {@inheritDoc} */
  218.     public double getMeanFlux(final AbsoluteDate date) throws OrekitException {

  219.         // get the neighboring dates
  220.         bracketDate(date);

  221.         // perform a linear interpolation
  222.         final AbsoluteDate previousDate = previousParam.getDate();
  223.         final AbsoluteDate currentDate  = currentParam.getDate();
  224.         final double dt                 = currentDate.durationFrom(previousDate);
  225.         final double previousF107       = previousParam.getF107();
  226.         final double currentF107        = currentParam.getF107();
  227.         final double previousWeight     = currentDate.durationFrom(date)  / dt;
  228.         final double currentWeight      = date.durationFrom(previousDate) / dt;

  229.         return previousF107 * previousWeight + currentF107 * currentWeight;

  230.     }

  231.     /** {@inheritDoc} */
  232.     public double getThreeHourlyKP(final AbsoluteDate date) throws OrekitException {
  233.         return get24HoursKp(date);
  234.     }

  235.     /** Get the date of the file from which data at the specified date comes from.
  236.      * <p>
  237.      * If several MSAFE files are available, some dates may appear in several
  238.      * files (for example August 2007 is in all files from the first one
  239.      * published in March 1999 to the February 2008 file). In this case, the
  240.      * data from the most recent file is used and the older ones are discarded.
  241.      * The date of the file is assumed to be 6 months after its first entry
  242.      * (which explains why the file having August 2007 as its first entry is the
  243.      * February 2008 file). This implies that MSAFE files must <em>not</em> be
  244.      * edited to change their time span, otherwise this would break the old
  245.      * entries overriding mechanism.
  246.      * </p>
  247.      * @param date date of the solar activity data
  248.      * @return date of the file
  249.      * @exception OrekitException if specified date is out of range
  250.      */
  251.     public DateComponents getFileDate(final AbsoluteDate date) throws OrekitException {
  252.         bracketDate(date);
  253.         final double dtP = date.durationFrom(previousParam.getDate());
  254.         final double dtC = currentParam.getDate().durationFrom(date);
  255.         return (dtP < dtC) ? previousParam.getFileDate() : currentParam.getFileDate();
  256.     }

  257.     /** The Kp index is derived from the Ap index.
  258.      * <p>The method used is explained on <a
  259.      * href="http://www.ngdc.noaa.gov/stp/GEOMAG/kp_ap.html">
  260.      * NOAA website.</a> as follows:</p>
  261.      * <p>The scale is 0 to 9 expressed in thirds of a unit, e.g. 5- is 4 2/3,
  262.      * 5 is 5 and 5+ is 5 1/3. The ap (equivalent range) index is derived from
  263.      * the Kp index as follows:</p>
  264.      * <table border="1">
  265.      * <tbody>
  266.      * <tr>
  267.      * <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>
  268.      * </tr>
  269.      * <tr>
  270.      * <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>
  271.      * </tr>
  272.      * <tr>
  273.      * <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>
  274.      * </tr>
  275.      * <tr>
  276.      * <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>
  277.      * </tr>
  278.      * </tbody>
  279.      * </table>
  280.      * @param date date of the Kp data
  281.      * @return the 24H geomagnetic index
  282.      * @exception OrekitException if the date is out of range of available data
  283.      */
  284.     public double get24HoursKp(final AbsoluteDate date) throws OrekitException {

  285.         // get the neighboring dates
  286.         bracketDate(date);

  287.         // perform a linear interpolation
  288.         final AbsoluteDate previousDate = previousParam.getDate();
  289.         final AbsoluteDate currentDate  = currentParam.getDate();
  290.         final double dt                 = currentDate.durationFrom(previousDate);
  291.         final double previousAp         = previousParam.getAp();
  292.         final double currentAp          = currentParam.getAp();
  293.         final double previousWeight     = currentDate.durationFrom(date)  / dt;
  294.         final double currentWeight      = date.durationFrom(previousDate) / dt;
  295.         final double ap                 = previousAp * previousWeight + currentAp * currentWeight;

  296.         // calculating Ap index, then corresponding Kp index
  297.         final int i = Arrays.binarySearch(AP_ARRAY, ap);
  298.         if (i >= 0) {
  299.             // the exact value for ap has been found, return the corresponding Kp
  300.             return KP_ARRAY[i];
  301.         } else {
  302.             // the exact value has not been found, we have an insertion point
  303.             final int jSup = -(i + 1);
  304.             final int jInf = jSup - 1;
  305.             if ((ap - AP_ARRAY[jInf]) < (AP_ARRAY[jSup] - ap)) {
  306.                 return KP_ARRAY[jInf];
  307.             } else {
  308.                 return KP_ARRAY[jSup];
  309.             }
  310.         }

  311.     }

  312.     /** Container class for Solar activity indexes.  */
  313.     private static class LineParameters implements TimeStamped, Serializable {

  314.         /** Serializable UID. */
  315.         private static final long serialVersionUID = 6607862001953526475L;

  316.         /** File date. */
  317.         private final DateComponents fileDate;

  318.         /** Entry date. */
  319.         private  final AbsoluteDate date;

  320.         /** F10.7 flux at date. */
  321.         private final double f107;

  322.         /** Ap index at date. */
  323.         private final double ap;

  324.         /** Simple constructor.
  325.          * @param fileDate file date
  326.          * @param date entry date
  327.          * @param f107 F10.7 flux at date
  328.          * @param ap Ap index at date
  329.          */
  330.         private LineParameters(final DateComponents fileDate, final AbsoluteDate date, final double f107, final double ap) {
  331.             this.fileDate = fileDate;
  332.             this.date     = date;
  333.             this.f107     = f107;
  334.             this.ap       = ap;
  335.         }

  336.         /** Get the file date.
  337.          * @return file date
  338.          */
  339.         public DateComponents getFileDate() {
  340.             return fileDate;
  341.         }

  342.         /** Get the current date.
  343.          * @return current date
  344.          */
  345.         public AbsoluteDate getDate() {
  346.             return date;
  347.         }

  348.         /** Get the F10.0 flux.
  349.          * @return f10.7 flux
  350.          */
  351.         public double getF107() {
  352.             return f107;
  353.         }

  354.         /** Get the Ap index.
  355.          * @return Ap index
  356.          */
  357.         public double getAp() {
  358.             return ap;
  359.         }

  360.     }

  361.     /** {@inheritDoc} */
  362.     public void loadData(final InputStream input, final String name)
  363.         throws IOException, ParseException, OrekitException {

  364.         // select the groups we want to store
  365.         final int f107Group;
  366.         final int apGroup;
  367.         switch (strengthLevel) {
  368.             case STRONG :
  369.                 f107Group = 3;
  370.                 apGroup   = 6;
  371.                 break;
  372.             case AVERAGE :
  373.                 f107Group = 4;
  374.                 apGroup   = 7;
  375.                 break;
  376.             default :
  377.                 f107Group = 5;
  378.                 apGroup   = 8;
  379.                 break;
  380.         }

  381.         // read the data
  382.         final BufferedReader reader = new BufferedReader(new InputStreamReader(input, "UTF-8"));
  383.         boolean inData = false;
  384.         final TimeScale utc = TimeScalesFactory.getUTC();
  385.         DateComponents fileDate = null;
  386.         for (String line = reader.readLine(); line != null; line = reader.readLine()) {
  387.             line = line.trim();
  388.             if (line.length() > 0) {
  389.                 final Matcher matcher = dataPattern.matcher(line);
  390.                 if (matcher.matches()) {

  391.                     // we are in the data section
  392.                     inData = true;

  393.                     // extract the data from the line
  394.                     final int year = Integer.parseInt(matcher.group(1));
  395.                     final Month month = Month.parseMonth(matcher.group(2));
  396.                     final AbsoluteDate date = new AbsoluteDate(year, month, 1, utc);
  397.                     if (fileDate == null) {
  398.                         // the first entry of each file correspond exactly to 6 months before file publication
  399.                         // so we compute the file date by adding 6 months to its first entry
  400.                         if (month.getNumber() > 6) {
  401.                             fileDate = new DateComponents(year + 1, month.getNumber() - 6, 1);
  402.                         } else {
  403.                             fileDate = new DateComponents(year, month.getNumber() + 6, 1);
  404.                         }
  405.                     }

  406.                     // check if there is already an entry for this date or not
  407.                     boolean addEntry = false;
  408.                     final Iterator<TimeStamped> iterator = data.tailSet(date).iterator();
  409.                     if (iterator.hasNext()) {
  410.                         final LineParameters existingEntry = (LineParameters) iterator.next();
  411.                         if (existingEntry.getDate().equals(date)) {
  412.                             // there is an entry for this date
  413.                             if (existingEntry.getFileDate().compareTo(fileDate) < 0) {
  414.                                 // the entry was read from an earlier file
  415.                                 // we replace it with the new entry as it is fresher
  416.                                 iterator.remove();
  417.                                 addEntry = true;
  418.                             }
  419.                         } else {
  420.                             // it is the first entry we get for this date
  421.                             addEntry = true;
  422.                         }
  423.                     } else {
  424.                         // it is the first entry we get for this date
  425.                         addEntry = true;
  426.                     }
  427.                     if (addEntry) {
  428.                         // we must add the new entry
  429.                         data.add(new LineParameters(fileDate, date,
  430.                                                     Double.parseDouble(matcher.group(f107Group)),
  431.                                                     Double.parseDouble(matcher.group(apGroup))));
  432.                     }

  433.                 } else {
  434.                     if (inData) {
  435.                         // we have already read some data, so we are not in the header anymore
  436.                         // however, we don't recognize this non-empty line,
  437.                         // we consider the file is corrupted
  438.                         throw new OrekitException(OrekitMessages.NOT_A_MARSHALL_SOLAR_ACTIVITY_FUTURE_ESTIMATION_FILE,
  439.                                                   name);
  440.                     }
  441.                 }
  442.             }
  443.         }

  444.         if (data.isEmpty()) {
  445.             throw new OrekitException(OrekitMessages.NOT_A_MARSHALL_SOLAR_ACTIVITY_FUTURE_ESTIMATION_FILE,
  446.                                       name);
  447.         }
  448.         firstDate = data.first().getDate();
  449.         lastDate  = data.last().getDate();

  450.     }

  451.     /** {@inheritDoc} */
  452.     public boolean stillAcceptsData() {
  453.         return true;
  454.     }

  455. }