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;
18  
19  import org.hipparchus.exception.DummyLocalizable;
20  import org.orekit.errors.OrekitException;
21  import org.orekit.errors.OrekitMessages;
22  import org.orekit.models.earth.atmosphere.DTM2000InputParameters;
23  import org.orekit.time.AbsoluteDate;
24  import org.orekit.time.ChronologicalComparator;
25  import org.orekit.time.TimeScalesFactory;
26  import org.orekit.time.TimeStamped;
27  import org.orekit.utils.Constants;
28  
29  import java.io.BufferedReader;
30  import java.io.IOException;
31  import java.io.InputStream;
32  import java.io.InputStreamReader;
33  import java.io.Serializable;
34  import java.util.Calendar;
35  import java.util.GregorianCalendar;
36  import java.util.SortedSet;
37  import java.util.TimeZone;
38  import java.util.TreeSet;
39  
40  
41  /** This class reads and provides solar activity data needed by the
42   * two atmospheric models. The data are furnished at the <a
43   * href="http://sol.spacenvironment.net/~JB2006/">
44   * official JB2006 website.</a>
45   *
46   * @author Fabien Maussion
47   */
48  public class SolarInputs97to05 implements DTM2000InputParameters {
49  
50      /** Serializable UID. */
51      private static final long serialVersionUID = -3687601846334870069L;
52  
53      private static final double third = 1.0/3.0;
54  
55      private static final double[] kpTab = new double[] {
56          0, 0+third, 1-third, 1, 1+third, 2-third, 2, 2+third,
57          3-third, 3, 3+third, 4-third, 4, 4+third, 5-third, 5,
58          5+third, 6-third, 6, 6+third, 7-third, 7, 7+third,
59          8-third, 8, 8+third, 9-third, 9
60      };
61  
62      private static final double[] apTab = new double[] {
63          0, 2, 3, 4, 5, 6, 7, 9, 12, 15, 18, 22, 27, 32,
64          39, 48, 56, 67, 80, 94, 111, 132, 154 , 179, 207, 236, 300, 400
65      };
66  
67      /** All entries. */
68      private SortedSet<TimeStamped> data;
69  
70      private LineParameters currentParam;
71      private AbsoluteDate firstDate;
72      private AbsoluteDate lastDate;
73  
74      /** Simple constructor.
75       * Data file address is set internally, nothing to be done here.
76       */
77      private SolarInputs97to05() {
78  
79          data = new TreeSet<TimeStamped>(new ChronologicalComparator());
80          InputStream in = SolarInputs97to05.class.getResourceAsStream("/atmosphere/JB_All_97-05.txt");
81          BufferedReader rFlux = new BufferedReader(new InputStreamReader(in));
82  
83  
84          in = SolarInputs97to05.class.getResourceAsStream("/atmosphere/NOAA_ap_97-05.dat.txt");
85          BufferedReader rAp = new BufferedReader(new InputStreamReader(in));
86  
87          try {
88              read(rFlux, rAp);
89          } catch (IOException e) {
90              throw new RuntimeException(e);
91          }
92      }
93  
94      /** Singleton getter.
95       * @return the unique instance of this class.
96       */
97      public static SolarInputs97to05 getInstance() {
98          if (LazyHolder.instance == null) {
99              throw LazyHolder.orekitException;
100         }
101         return LazyHolder.instance;
102     }
103 
104     private void read(BufferedReader rFlux, BufferedReader rAp) throws IOException, OrekitException {
105 
106         rFlux.readLine();
107         rFlux.readLine();
108         rFlux.readLine();
109         rFlux.readLine();
110         rAp.readLine();
111         String lineAp;
112         String[] flux;
113         String[] ap;
114         Calendar cal = new GregorianCalendar();
115         cal.setTimeZone(TimeZone.getTimeZone("UTC"));
116         cal.set(0, 0, 0, 0, 0, 0);
117         cal.set(Calendar.MILLISECOND, 0);
118 
119         AbsoluteDate date = null;
120         boolean first = true;
121 
122         for (String lineFlux = rFlux.readLine(); lineFlux != null; lineFlux = rFlux.readLine()) {
123 
124             flux = lineFlux.trim().split("\\s+");
125 
126             lineAp = rAp.readLine();
127             if (lineAp == null) {
128                 throw new OrekitException(new DummyLocalizable("inconsistent JB2006 and geomagnetic indices files"));
129             }
130             ap = lineAp.trim().split("\\s+");
131 
132             int fluxYear = Integer.parseInt(flux[0]);
133             int fluxDay = Integer.parseInt(flux[1]);
134             int apYear  = Integer.parseInt(ap[11]);
135 
136             if (fluxDay != Integer.parseInt(ap[0])) {
137                 throw new OrekitException(new DummyLocalizable("inconsistent JB2006 and geomagnetic indices files"));
138             }
139             if (((fluxYear <  2000) && ((fluxYear - 1900) != apYear)) ||
140                 ((fluxYear >= 2000) && ((fluxYear - 2000) != apYear))) {
141                 throw new OrekitException(new DummyLocalizable("inconsistent JB2006 and geomagnetic indices files"));
142             }
143 
144             cal.set(Calendar.YEAR, fluxYear);
145             cal.set(Calendar.DAY_OF_YEAR, fluxDay);
146 
147             date = new AbsoluteDate(cal.getTime(), TimeScalesFactory.getUTC());
148 
149             if (first) {
150                 first = false;
151                 firstDate = date;
152             }
153 
154             data.add(new LineParameters(date,
155                                         new double[] {
156                     Double.parseDouble(ap[3]),
157                     Double.parseDouble(ap[4]),
158                     Double.parseDouble(ap[5]),
159                     Double.parseDouble(ap[6]),
160                     Double.parseDouble(ap[7]),
161                     Double.parseDouble(ap[8]),
162                     Double.parseDouble(ap[9]),
163                     Double.parseDouble(ap[10]),
164 
165             },
166             Double.parseDouble(flux[3]),
167             Double.parseDouble(flux[4]),
168             Double.parseDouble(flux[5]),
169             Double.parseDouble(flux[6]),
170             Double.parseDouble(flux[7]),
171             Double.parseDouble(flux[8])));
172 
173         }
174         lastDate = date;
175 
176     }
177 
178     private void findClosestLine(AbsoluteDate date) {
179 
180         if ((date.durationFrom(firstDate) < 0) || (date.durationFrom(lastDate) > Constants.JULIAN_DAY)) {
181             throw new OrekitException(OrekitMessages.OUT_OF_RANGE_DATE, date, firstDate, lastDate);
182         }
183 
184         // don't search if the cached selection is fine
185         if ((currentParam != null) && (date.durationFrom(currentParam.date) >= 0) &&
186                 (date.durationFrom(currentParam.date) < Constants.JULIAN_DAY )) {
187             return;
188         }
189         LineParameters before = new LineParameters(date.shiftedBy(-Constants.JULIAN_DAY), null, 0, 0, 0, 0, 0, 0);
190 
191         // search starting from entries a few steps before the target date
192         SortedSet<TimeStamped> tailSet = data.tailSet(before);
193         if (tailSet != null) {
194             currentParam = (LineParameters) tailSet.first();
195             if (currentParam.date.durationFrom(date) == -Constants.JULIAN_DAY) {
196                 currentParam = (LineParameters) data.tailSet(date).first();
197             }
198         } else {
199             throw new OrekitException(new DummyLocalizable("unable to find data for date {0}"), date);
200         }
201     }
202 
203     /** Container class for Solar activity indexes.  */
204     private static class LineParameters implements TimeStamped, Serializable {
205 
206         /** Serializable UID. */
207         private static final long serialVersionUID = -1127762834954768272L;
208 
209         /** Entries */
210         private  final AbsoluteDate date;
211         private final double[] ap;
212         private final double f10;
213         private final double f10B;
214         private final double s10;
215         private final double s10B;
216         private final double xm10;
217         private final double xm10B;
218 
219 
220         /** Simple constructor. */
221         private LineParameters(AbsoluteDate date, double[]  ap, double f10,
222                                double f10B, double s10, double s10B,
223                                double xm10, double xm10B) {
224             this.date = date;
225             this.ap = ap;
226             this.f10 = f10;
227             this.f10B = f10B;
228             this.s10 = s10;
229             this.s10B = s10B;
230             this.xm10 = xm10;
231             this.xm10B = xm10B;
232 
233         }
234 
235         /** Get the current date */
236         public AbsoluteDate getDate() {
237             return date;
238         }
239 
240     }
241 
242     public double getAp(AbsoluteDate date) {
243         double result = Double.NaN;
244         try {
245             findClosestLine(date);
246             Calendar cal = new GregorianCalendar();
247             cal.setTimeZone(TimeZone.getTimeZone("UTC"));
248             cal.setTime(date.toDate(TimeScalesFactory.getUTC()));
249             int hour = cal.get(Calendar.HOUR_OF_DAY);
250             for (int i= 0; i<8; i++) {
251                 if ((hour >= (i * 3)) && (hour < ((i + 1) * 3))) {
252                     result = currentParam.ap[i];
253                 }
254             }
255         }
256         catch (OrekitException e) {
257             // nothing
258         }
259         return result;
260     }
261 
262     public double getF10(AbsoluteDate date) {
263         double result = Double.NaN;
264         try {
265             findClosestLine(date);
266             result=currentParam.f10;
267         }
268         catch (OrekitException e) {
269             // nothing
270         }
271         return result;
272     }
273 
274     public double getF10B(AbsoluteDate date) {
275         double result = Double.NaN;
276         try {
277             findClosestLine(date);
278             result=currentParam.f10B;
279         }
280         catch (OrekitException e) {
281             // nothing
282         }
283         return result;
284     }
285 
286     public AbsoluteDate getMaxDate() {
287         return lastDate.shiftedBy(Constants.JULIAN_DAY);
288     }
289 
290     public AbsoluteDate getMinDate() {
291         return firstDate;
292     }
293 
294     public double getS10(AbsoluteDate date) {
295         double result = Double.NaN;
296         try {
297             findClosestLine(date);
298             result=currentParam.s10;
299         }
300         catch (OrekitException e) {
301             // nothing
302         }
303         return result;
304     }
305 
306     public double getS10B(AbsoluteDate date) {
307         double result = Double.NaN;
308         try {
309             findClosestLine(date);
310             result=currentParam.s10B;
311         }
312         catch (OrekitException e) {
313             // nothing
314         }
315         return result;
316     }
317 
318     public double getXM10(AbsoluteDate date) {
319         double result = Double.NaN;
320         try {
321             findClosestLine(date);
322             result=currentParam.xm10;
323         }
324         catch (OrekitException e) {
325             // nothing
326         }
327         return result;
328     }
329 
330     public double getXM10B(AbsoluteDate date) {
331         double result = Double.NaN;
332         try {
333             findClosestLine(date);
334             result=currentParam.xm10B;
335         }
336         catch (OrekitException e) {
337             // nothing
338         }
339         return result;
340     }
341 
342     public double get24HoursKp(AbsoluteDate date) {
343         double result = 0;
344         AbsoluteDate myDate = date;
345 
346         for (int i=0; i<8; i++) {
347             result += getThreeHourlyKP(date);
348             myDate = myDate.shiftedBy(3 * 3600);
349         }
350 
351         return result/8;
352     }
353 
354     public double getInstantFlux(AbsoluteDate date) {
355         return getF10(date);
356     }
357 
358     public double getMeanFlux(AbsoluteDate date) {
359         return getF10B(date);
360     }
361 
362     /** The 3-H Kp is derived from the Ap index.
363      * The used method is explained on <a
364      * href="http://www.ngdc.noaa.gov/stp/GEOMAG/kp_ap.shtml">
365      * NOAA website.</a>. Here is the corresponding tab :
366      * <pre>
367      * The scale is O to 9 expressed in thirds of a unit, e.g. 5- is 4 2/3,
368      * 5 is 5 and 5+ is 5 1/3.
369      *
370      * The 3-hourly ap (equivalent range) index is derived from the Kp index as follows:
371      *
372      * Kp = 0o   0+   1-   1o   1+   2-   2o   2+   3-   3o   3+   4-   4o   4+
373      * ap =  0    2    3    4    5    6    7    9   12   15   18   22   27   32
374      * Kp = 5-   5o   5+   6-   6o   6+   7-   7o   7+   8-   8o   8+   9-   9o
375      * ap = 39   48   56   67   80   94  111  132  154  179  207  236  300  400
376      *
377      * </pre>
378      */
379     public double getThreeHourlyKP(AbsoluteDate date) {
380 
381         double ap = getAp(date);
382         int i = 0;
383         for ( i= 0; ap>=apTab[i]; i++) {
384             if (i==apTab.length-1) {
385                 i++;
386                 break;
387             }
388         }
389         return kpTab[i-1];
390     }
391 
392     /** Holder for the singleton.
393      * <p>We use the Initialization on demand holder idiom to store
394      * the singleton, as it is both thread-safe, efficient (no
395      * synchronization) and works with all versions of java.</p>
396      */
397     private static class LazyHolder {
398         private static final SolarInputs97to05 instance;
399         private static final OrekitException orekitException;
400         static {
401             SolarInputs97to05 tmpInstance = null;
402             OrekitException tmpException = null;
403             try {
404                 tmpInstance = new SolarInputs97to05();
405             } catch (OrekitException oe) {
406                 tmpException = oe;
407             }
408             instance        = tmpInstance;
409             orekitException = tmpException;
410         }
411     }
412 
413 }