1 /* Copyright 2020 Clément Jonglez
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 * Clément Jonglez 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
18 package org.orekit.models.earth.atmosphere.data;
19
20 import java.util.List;
21 import java.util.stream.Collectors;
22
23 import org.orekit.annotation.DefaultDataContext;
24 import org.orekit.data.AbstractSelfFeedingLoader;
25 import org.orekit.data.DataContext;
26 import org.orekit.data.DataProvidersManager;
27 import org.orekit.errors.OrekitException;
28 import org.orekit.errors.OrekitMessages;
29 import org.orekit.models.earth.atmosphere.DTM2000InputParameters;
30 import org.orekit.models.earth.atmosphere.NRLMSISE00InputParameters;
31 import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherDataLoader.LineParameters;
32 import org.orekit.time.AbsoluteDate;
33 import org.orekit.time.TimeScale;
34 import org.orekit.time.TimeStamped;
35 import org.orekit.utils.Constants;
36 import org.orekit.utils.ImmutableTimeStampedCache;
37
38 /**
39 * This class provides three-hourly and daily solar activity data needed by atmospheric
40 * models: F107 solar flux, Ap and Kp indexes.
41 * The {@link org.orekit.data.DataLoader} implementation and the parsing is handled by the class {@link CssiSpaceWeatherDataLoader}.
42 * <p>
43 * The data are retrieved through space weather files offered by AGI/CSSI on the AGI
44 * <a href="ftp://ftp.agi.com/pub/DynamicEarthData/SpaceWeather-All-v1.2.txt">FTP</a> as
45 * well as on the CelesTrack <a href="http://celestrak.com/SpaceData/">website</a>.
46 * These files are updated several times a day by using several sources mentioned in the
47 * <a href="http://celestrak.com/SpaceData/SpaceWx-format.php">Celestrak space
48 * weather data documentation</a>.
49 * </p>
50 *
51 * @author Clément Jonglez
52 * @since 10.2
53 */
54 public class CssiSpaceWeatherData extends AbstractSelfFeedingLoader
55 implements DTM2000InputParameters, NRLMSISE00InputParameters {
56
57 /** Default regular expression for supported names that works with all officially published files. */
58 public static final String DEFAULT_SUPPORTED_NAMES = "^S(?:pace)?W(?:eather)?-(?:All)?.*\\.txt$";
59
60 /** Serializable UID. */
61 private static final long serialVersionUID = 4249411710645968978L;
62
63 /** Size of the list. */
64 private static final int N_NEIGHBORS = 2;
65
66 /** Data set. */
67 private final transient ImmutableTimeStampedCache<TimeStamped> data;
68
69 /** UTC time scale. */
70 private final TimeScale utc;
71
72 /** First available date. */
73 private final AbsoluteDate firstDate;
74
75 /** Date of last data before the prediction starts. */
76 private final AbsoluteDate lastObservedDate;
77
78 /** Date of last daily prediction before the monthly prediction starts. */
79 private final AbsoluteDate lastDailyPredictedDate;
80
81 /** Last available date. */
82 private final AbsoluteDate lastDate;
83
84 /** Previous set of solar activity parameters. */
85 private LineParameters previousParam;
86
87 /** Current set of solar activity parameters. */
88 private LineParameters nextParam;
89
90 /**
91 * Simple constructor. This constructor uses the default data context.
92 * <p>
93 * The original file names provided by AGI/CSSI are of the form:
94 * SpaceWeather-All-v1.2.txt (AGI's ftp) or SW-Last5Years.txt (CelesTrak's website).
95 * So a recommended regular expression for the supported names that works
96 * with all published files is: {@link #DEFAULT_SUPPORTED_NAMES}.
97 * </p>
98 *
99 * @param supportedNames regular expression for supported AGI/CSSI space weather files names
100 */
101 @DefaultDataContext
102 public CssiSpaceWeatherData(final String supportedNames) {
103 this(supportedNames, DataContext.getDefault().getDataProvidersManager(),
104 DataContext.getDefault().getTimeScales().getUTC());
105 }
106
107 /**
108 * Constructor that allows specifying the source of the CSSI space weather
109 * file.
110 *
111 * @param supportedNames regular expression for supported AGI/CSSI space weather files names
112 * @param dataProvidersManager provides access to auxiliary data files.
113 * @param utc UTC time scale.
114 */
115 public CssiSpaceWeatherData(final String supportedNames,
116 final DataProvidersManager dataProvidersManager,
117 final TimeScale utc) {
118 super(supportedNames, dataProvidersManager);
119
120 this.utc = utc;
121 final CssiSpaceWeatherDataLoader loader =
122 new CssiSpaceWeatherDataLoader(utc);
123 this.feed(loader);
124 data =
125 new ImmutableTimeStampedCache<>(N_NEIGHBORS, loader.getDataSet());
126 firstDate = loader.getMinDate();
127 lastDate = loader.getMaxDate();
128 lastObservedDate = loader.getLastObservedDate();
129 lastDailyPredictedDate = loader.getLastDailyPredictedDate();
130 }
131
132 /** {@inheritDoc} */
133 public AbsoluteDate getMinDate() {
134 return firstDate;
135 }
136
137 /** {@inheritDoc} */
138 public AbsoluteDate getMaxDate() {
139 return lastDate;
140 }
141
142 /**
143 * Find the data bracketing a specified date.
144 *
145 * @param date date to bracket
146 */
147 private void bracketDate(final AbsoluteDate date) {
148 if (date.durationFrom(firstDate) < 0) {
149 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE,
150 date, firstDate, lastDate, firstDate.durationFrom(date));
151 }
152 if (date.durationFrom(lastDate) > 0) {
153 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER,
154 date, firstDate, lastDate, date.durationFrom(lastDate));
155 }
156
157 // don't search if the cached selection is fine
158 if (previousParam != null && date.durationFrom(previousParam.getDate()) > 0 &&
159 date.durationFrom(nextParam.getDate()) <= 0) {
160 return;
161 }
162
163 final List<TimeStamped> neigbors = data.getNeighbors(date).collect(Collectors.toList());
164 previousParam = (LineParameters) neigbors.get(0);
165 nextParam = (LineParameters) neigbors.get(1);
166 if (previousParam.getDate().compareTo(date) > 0) { // TODO delete dead code
167 /**
168 * Throwing exception if neighbors are unbalanced because we are at the
169 * beginning of the data set
170 */
171 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, date, firstDate, lastDate);
172 }
173 }
174
175 /**
176 * Performs a linear interpolation between two values The weights are computed
177 * from the time delta between previous date, current date, next date.
178 *
179 * @param date the current date
180 * @param previousValue the value at previous date
181 * @param nextValue the value at next date
182 * @return the value interpolated for the current date
183 */
184 private double getLinearInterpolation(final AbsoluteDate date, final double previousValue, final double nextValue) {
185 // perform a linear interpolation
186 final AbsoluteDate previousDate = previousParam.getDate();
187 final AbsoluteDate currentDate = nextParam.getDate();
188 final double dt = currentDate.durationFrom(previousDate);
189 final double previousWeight = currentDate.durationFrom(date) / dt;
190 final double nextWeight = date.durationFrom(previousDate) / dt;
191
192 // returns the data interpolated at the date
193 return previousValue * previousWeight + nextValue * nextWeight;
194 }
195
196 /** {@inheritDoc} */
197 public double getInstantFlux(final AbsoluteDate date) {
198 // Interpolating two neighboring daily fluxes
199 // get the neighboring dates
200 bracketDate(date);
201 return getLinearInterpolation(date, previousParam.getF107Obs(), nextParam.getF107Obs());
202 }
203
204 /** {@inheritDoc} */
205 public double getMeanFlux(final AbsoluteDate date) {
206 return getAverageFlux(date);
207 }
208
209 /** {@inheritDoc} */
210 public double getThreeHourlyKP(final AbsoluteDate date) {
211 if (date.compareTo(lastObservedDate) <= 0) {
212 /**
213 * If observation data is available, it contains three-hourly data
214 */
215 bracketDate(date);
216 final double hourOfDay = date.offsetFrom(previousParam.getDate(), utc) / 3600;
217 int i_kp = (int) (hourOfDay / 3);
218 if (i_kp >= 8) {
219 /**
220 * hourOfDay can take the value 24.0 at midnight due to floating point precision
221 * when bracketing the dates or during a leap second because the hour of day is
222 * computed in UTC view
223 */
224 i_kp = 7;
225 }
226 return previousParam.getThreeHourlyKp(i_kp);
227 } else {
228 /**
229 * Only predictions are available, there are no three-hourly data
230 */
231 return get24HoursKp(date);
232 }
233 }
234
235 /** {@inheritDoc} */
236 public double get24HoursKp(final AbsoluteDate date) {
237 if (date.compareTo(lastDailyPredictedDate) <= 0) {
238 // Daily data is available, just taking the daily average
239 bracketDate(date);
240 return previousParam.getKpSum() / 8;
241 } else {
242 // Only monthly data is available, better interpolate between two months
243 // get the neighboring dates
244 bracketDate(date);
245 return getLinearInterpolation(date, previousParam.getKpSum() / 8, nextParam.getKpSum() / 8);
246 }
247 }
248
249 /** {@inheritDoc} */
250 public double getDailyFlux(final AbsoluteDate date) {
251 // Getting the value for the previous day
252 return getDailyFluxOnDay(date.shiftedBy(-Constants.JULIAN_DAY));
253 }
254
255 /**
256 * Gets the daily flux on the current day.
257 *
258 * @param date the current date
259 * @return the daily F10.7 flux (observed)
260 */
261 private double getDailyFluxOnDay(final AbsoluteDate date) {
262 if (date.compareTo(lastDailyPredictedDate) <= 0) {
263 // Getting the value for the previous day
264 bracketDate(date);
265 return previousParam.getF107Obs();
266 } else {
267 // Only monthly data is available, better interpolate between two months
268 // get the neighboring dates
269 bracketDate(date);
270 return getLinearInterpolation(date, previousParam.getF107Obs(), nextParam.getF107Obs());
271 }
272 }
273
274 /** {@inheritDoc} */
275 public double getAverageFlux(final AbsoluteDate date) {
276 if (date.compareTo(lastDailyPredictedDate) <= 0) {
277 bracketDate(date);
278 return previousParam.getCtr81Obs();
279 } else {
280 // Only monthly data is available, better interpolate between two months
281 // get the neighboring dates
282 bracketDate(date);
283 return getLinearInterpolation(date, previousParam.getCtr81Obs(), nextParam.getCtr81Obs());
284 }
285 }
286
287 /** {@inheritDoc} */
288 public double[] getAp(final AbsoluteDate date) {
289 final double[] apArray = new double[7];
290 apArray[0] = getDailyAp(date);
291 apArray[1] = getThreeHourlyAp(date);
292 apArray[2] = getThreeHourlyAp(date.shiftedBy(-3.0 * 3600.0));
293 apArray[3] = getThreeHourlyAp(date.shiftedBy(-6.0 * 3600.0));
294 apArray[4] = getThreeHourlyAp(date.shiftedBy(-9.0 * 3600.0));
295 apArray[5] = get24HoursAverageAp(date.shiftedBy(-12.0 * 3600.0));
296 apArray[6] = get24HoursAverageAp(date.shiftedBy(-36.0 * 3600.0));
297 return apArray;
298 }
299
300 /**
301 * Gets the value of the three-hourly Ap index for the given date.
302 *
303 * @param date the current date
304 * @return the current three-hourly Ap index
305 */
306 private double getThreeHourlyAp(final AbsoluteDate date) {
307 if (date.compareTo(lastObservedDate.shiftedBy(Constants.JULIAN_DAY)) < 0) {
308 /**
309 * If observation data is available, it contains three-hourly data.
310 */
311 bracketDate(date);
312 final double hourOfDay = date.offsetFrom(previousParam.getDate(), utc) / 3600;
313 int i_ap = (int) (hourOfDay / 3);
314 if (i_ap >= 8) {
315 /**
316 * hourOfDay can take the value 24.0 at midnight due to floating point precision
317 * when bracketing the dates or during a leap second because the hour of day is
318 * computed in UTC view
319 */
320 i_ap = 7;
321 }
322 return previousParam.getThreeHourlyAp(i_ap);
323 } else {
324 /**
325 * Only predictions are available, there are no three-hourly data
326 */
327 return getDailyAp(date);
328 }
329 }
330
331 /**
332 * Gets the running average of the 8 three-hourly Ap indices prior to current
333 * time If three-hourly data is available, the result is different than
334 * getDailyAp.
335 *
336 * @param date the current date
337 * @return the 24 hours running average of the Ap index
338 */
339 private double get24HoursAverageAp(final AbsoluteDate date) {
340 if (date.compareTo(lastDailyPredictedDate) <= 0) {
341 // Computing running mean
342 double apSum = 0.0;
343 for (int i = 0; i < 8; i++) {
344 apSum += getThreeHourlyAp(date.shiftedBy(-3.0 * 3600 * i));
345 }
346 return apSum / 8;
347 } else {
348 /**
349 * Only monthly predictions are available, no need to compute the average from
350 * three hourly data
351 */
352 return getDailyAp(date);
353 }
354 }
355
356 /**
357 * Get the daily Ap index for the given date.
358 *
359 * @param date the current date
360 * @return the daily Ap index
361 */
362 private double getDailyAp(final AbsoluteDate date) {
363 if (date.compareTo(lastDailyPredictedDate) <= 0) {
364 // Daily data is available, just taking the daily average
365 bracketDate(date);
366 return previousParam.getApAvg();
367 } else {
368 // Only monthly data is available, better interpolate between two months
369 // get the neighboring dates
370 bracketDate(date);
371 return getLinearInterpolation(date, previousParam.getApAvg(), nextParam.getApAvg());
372 }
373 }
374
375 public String getSupportedNames() {
376 return super.getSupportedNames();
377 }
378 }