1 /* Copyright 2002-2025 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;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.util.FastMath;
21 import org.hipparchus.util.FieldSinCos;
22 import org.hipparchus.util.MathUtils;
23 import org.hipparchus.util.SinCos;
24 import org.orekit.annotation.DefaultDataContext;
25 import org.orekit.bodies.BodyShape;
26 import org.orekit.data.DataContext;
27 import org.orekit.time.AbsoluteDate;
28 import org.orekit.time.FieldAbsoluteDate;
29 import org.orekit.time.TimeScale;
30
31 import org.orekit.utils.Constants;
32 import org.orekit.utils.ExtendedPositionProvider;
33
34 /** This is the realization of the Jacchia-Bowman 2008 atmospheric model.
35 * <p>
36 * It is described in the paper:<br>
37 * <a href="https://www.researchgate.net/publication/228621668_A_New_Empirical_Thermospheric_Density_Model_JB2008_Using_New_Solar_and_Geomagnetic_Indices">A
38 * New Empirical Thermospheric Density Model JB2008 Using New Solar Indices</a><br>
39 * <i>Bruce R. Bowman & al.</i><br>
40 * AIAA 2008-6438<br>
41 * </p>
42 * <p>
43 * Two computation methods are proposed to the user:
44 * <ul>
45 * <li>one OREKIT independent and compliant with initial FORTRAN routine entry values:
46 * {@link #getDensity(double, double, double, double, double, double, double, double, double, double, double, double, double, double, double)}. </li>
47 * <li>one compliant with OREKIT Atmosphere interface, necessary to the
48 * {@link org.orekit.forces.drag.DragForce drag force model} computation.</li>
49 * </ul>
50 * <p>
51 * This model provides dense output for all altitudes and positions. Output data are :
52 * <ul>
53 * <li>Exospheric Temperature above Input Position (deg K)</li>
54 * <li>Temperature at Input Position (deg K)</li>
55 * <li>Total Mass-Density at Input Position (kg/m³)</li>
56 * </ul>
57 * <p>
58 * The model needs geographical and time information to compute general values,
59 * but also needs space weather data : mean and daily solar flux, retrieved through
60 * different indices, and planetary geomagnetic indices.<br>
61 * More information on these indices can be found on the <a
62 * href="http://sol.spacenvironment.net/~JB2008/indices.html">
63 * official JB2008 website.</a>
64 * </p>
65 *
66 * @author Bruce R Bowman (HQ AFSPC, Space Analysis Division), 2008: FORTRAN routine
67 * @author Pascal Parraud (java translation)
68 */
69 public class JB2008 extends AbstractJacchiaBowmanModel {
70
71 /** FZ global model values (1997-2006 fit). */
72 private static final double[] FZM = {
73 0.2689e+00, -0.1176e-01, 0.2782e-01, -0.2782e-01, 0.3470e-03
74 };
75
76 /** GT global model values (1997-2006 fit). */
77 private static final double[] GTM = {
78 -0.3633e+00, 0.8506e-01, 0.2401e+00, -0.1897e+00, -0.2554e+00,
79 -0.1790e-01, 0.5650e-03, -0.6407e-03, -0.3418e-02, -0.1252e-02
80 };
81
82 /** External data container. */
83 private final JB2008InputParameters inputParams;
84
85 /** Constructor with space environment information for internal computation.
86 *
87 * <p>This method uses the {@link DataContext#getDefault() default data context}.
88 *
89 * @param parameters the solar and magnetic activity data
90 * @param sun the sun position
91 * @param earth the earth body shape
92 * @see #JB2008(JB2008InputParameters, ExtendedPositionProvider, BodyShape, TimeScale)
93 */
94 @DefaultDataContext
95 public JB2008(final JB2008InputParameters parameters,
96 final ExtendedPositionProvider sun, final BodyShape earth) {
97 this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
98 }
99
100 /**
101 * Constructor with space environment information for internal computation.
102 *
103 * @param parameters the solar and magnetic activity data
104 * @param sun the sun position
105 * @param earth the earth body shape
106 * @param utc UTC time scale. Used to computed the day fraction.
107 * @since 10.1
108 */
109 public JB2008(final JB2008InputParameters parameters,
110 final ExtendedPositionProvider sun,
111 final BodyShape earth,
112 final TimeScale utc) {
113 super(sun, utc, earth,
114 parameters == null ? AbsoluteDate.PAST_INFINITY : parameters.getMinDate(),
115 parameters == null ? AbsoluteDate.FUTURE_INFINITY : parameters.getMaxDate());
116 this.inputParams = parameters;
117 }
118
119 /** Get the local density with initial entries.
120 * <p>
121 * The method creates a new instance of {@link JB2008} model and set the input
122 * solar activity data equal to the provided one. These data are then
123 * available form {@link AbsoluteDate#PAST_INFINITY} to {@link AbsoluteDate#FUTURE_INFINITY}.
124 * </p>
125 * @param dateMJD date and time, in modified julian days and fraction
126 * @param sunRA Right Ascension of Sun (radians).
127 * @param sunDecli Declination of Sun (radians).
128 * @param satLon Right Ascension of position (radians)
129 * @param satLat Geocentric latitude of position (radians)
130 * @param satAlt Height of position (m)
131 * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
132 * (Tabular time 1.0 day earlier)
133 * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
134 * (Tabular time 1.0 day earlier)
135 * @param s10 EUV index (26-34 nm) scaled to F10<br>
136 * (Tabular time 1 day earlier)
137 * @param s10B UV 81-day averaged centered index
138 * (Tabular time 1 day earlier)
139 * @param xm10 MG2 index scaled to F10<br>
140 * (Tabular time 2.0 days earlier)
141 * @param xm10B MG2 81-day ave. centered index<br>
142 * (Tabular time 2.0 days earlier)
143 * @param y10 Solar X-Ray & Lya index scaled to F10<br>
144 * (Tabular time 5.0 days earlier)
145 * @param y10B Solar X-Ray & Lya 81-day ave. centered index<br>
146 * (Tabular time 5.0 days earlier)
147 * @param dstdtc Temperature change computed from Dst index
148 * @return total mass-Density at input position (kg/m³)
149 */
150 public double getDensity(final double dateMJD, final double sunRA, final double sunDecli,
151 final double satLon, final double satLat, final double satAlt,
152 final double f10, final double f10B, final double s10,
153 final double s10B, final double xm10, final double xm10B,
154 final double y10, final double y10B, final double dstdtc) {
155 final LocalProvider provider = new LocalProvider(f10, f10B, s10, s10B, xm10, xm10B, y10, y10B, dstdtc);
156 final JB2008 modelWithLocalData = new JB2008(provider, getSun(), getEarth(), getUtc());
157 final double mjd = FastMath.floor(dateMJD);
158 final AbsoluteDate computationEpoch = AbsoluteDate.createMJDDate((int) mjd, (dateMJD - mjd) * Constants.JULIAN_DAY, getUtc());
159 return modelWithLocalData.computeDensity(computationEpoch, sunRA, sunDecli, satLon, satLat, satAlt);
160 }
161
162 /** Get the local density with initial entries.
163 * <p>
164 * The method creates a new instance of {@link JB2008} model and set the input
165 * solar activity data equal to the provided one. These data are then
166 * available form {@link AbsoluteDate#PAST_INFINITY} to {@link AbsoluteDate#FUTURE_INFINITY}.
167 * </p>
168 * @param dateMJD date and time, in modified julian days and fraction
169 * @param sunRA Right Ascension of Sun (radians).
170 * @param sunDecli Declination of Sun (radians).
171 * @param satLon Right Ascension of position (radians)
172 * @param satLat Geocentric latitude of position (radians)
173 * @param satAlt Height of position (m)
174 * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
175 * (Tabular time 1.0 day earlier)
176 * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
177 * (Tabular time 1.0 day earlier)
178 * @param s10 EUV index (26-34 nm) scaled to F10<br>
179 * (Tabular time 1 day earlier)
180 * @param s10B UV 81-day averaged centered index
181 * (Tabular time 1 day earlier)
182 * @param xm10 MG2 index scaled to F10<br>
183 * (Tabular time 2.0 days earlier)
184 * @param xm10B MG2 81-day ave. centered index<br>
185 * (Tabular time 2.0 days earlier)
186 * @param y10 Solar X-Ray & Lya index scaled to F10<br>
187 * (Tabular time 5.0 days earlier)
188 * @param y10B Solar X-Ray & Lya 81-day ave. centered index<br>
189 * (Tabular time 5.0 days earlier)
190 * @param dstdtc Temperature change computed from Dst index
191 * @param <T> type of the elements
192 * @return total mass-Density at input position (kg/m³)
193 */
194 public <T extends CalculusFieldElement<T>> T getDensity(final T dateMJD, final T sunRA, final T sunDecli,
195 final T satLon, final T satLat, final T satAlt,
196 final double f10, final double f10B, final double s10,
197 final double s10B, final double xm10, final double xm10B,
198 final double y10, final double y10B, final double dstdtc) {
199 final LocalProvider provider = new LocalProvider(f10, f10B, s10, s10B, xm10, xm10B, y10, y10B, dstdtc);
200 final JB2008 modelWithLocalData = new JB2008(provider, getSun(), getEarth(), getUtc());
201 final T mjd = FastMath.floor(dateMJD);
202 final FieldAbsoluteDate<T> computationEpoch = FieldAbsoluteDate.createMJDDate((int) mjd.getReal(), dateMJD.subtract(mjd).multiply(Constants.JULIAN_DAY), getUtc());
203 return modelWithLocalData.computeDensity(computationEpoch, sunRA, sunDecli, satLon, satLat, satAlt);
204 }
205
206 /** {@inheritDoc} */
207 @Override
208 protected double computeTInf(final AbsoluteDate date, final double tsubl, final double dtclst) {
209 return tsubl + inputParams.getDSTDTC(date) + dtclst;
210 }
211
212 /** {@inheritDoc} */
213 @Override
214 protected double computeTc(final AbsoluteDate date) {
215 final double f10 = inputParams.getF10(date);
216 final double f10B = inputParams.getF10B(date);
217 final double s10 = inputParams.getS10(date);
218 final double s10B = inputParams.getS10B(date);
219 final double xm10 = inputParams.getXM10(date);
220 final double xm10B = inputParams.getXM10B(date);
221 final double y10 = inputParams.getY10(date);
222 final double y10B = inputParams.getY10B(date);
223 final double fn = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
224 final double fsb = f10B * fn + s10B * (1. - fn);
225 return 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) + 0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);
226 }
227
228 /** {@inheritDoc} */
229 @Override
230 protected double getF10(final AbsoluteDate date) {
231 return inputParams.getF10(date);
232 }
233
234 /** {@inheritDoc} */
235 @Override
236 protected double getF10B(final AbsoluteDate date) {
237 return inputParams.getF10B(date);
238 }
239
240 /** {@inheritDoc} */
241 @Override
242 protected <T extends CalculusFieldElement<T>> T computeTInf(final AbsoluteDate date, final T tsubl, final T dtclst) {
243 return tsubl.add(inputParams.getDSTDTC(date)).add(dtclst);
244 }
245
246 /** {@inheritDoc} */
247 @Override
248 protected double semian(final AbsoluteDate date, final double day, final double altKm) {
249
250 final double f10B = inputParams.getF10B(date);
251 final double s10B = inputParams.getS10B(date);
252 final double xm10B = inputParams.getXM10B(date);
253
254 final double htz = altKm / 1000.0;
255
256 // COMPUTE NEW 81-DAY CENTERED SOLAR INDEX FOR FZ
257 double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;
258
259 // SEMIANNUAL AMPLITUDE
260 final double fzz = FZM[0] + fsmb * (FZM[1] + htz * (FZM[2] + FZM[3] * htz + FZM[4] * fsmb));
261
262 // COMPUTE DAILY 81-DAY CENTERED SOLAR INDEX FOR GT
263 fsmb = f10B - 0.75 * s10B - 0.37 * xm10B;
264
265 // SEMIANNUAL PHASE FUNCTION
266 final double tau = MathUtils.TWO_PI * (day - 1.0) / 365;
267 final SinCos sc1P = FastMath.sinCos(tau);
268 final SinCos sc2P = SinCos.sum(sc1P, sc1P);
269 final double gtz = GTM[0] + GTM[1] * sc1P.sin() + GTM[2] * sc1P.cos() + GTM[3] * sc2P.sin() + GTM[4] * sc2P.cos() +
270 fsmb * (GTM[5] + GTM[6] * sc1P.sin() + GTM[7] * sc1P.cos() + GTM[8] * sc2P.sin() + GTM[9] * sc2P.cos());
271
272 return FastMath.max(1.0e-6, fzz) * gtz;
273
274 }
275
276 /** {@inheritDoc} */
277 @Override
278 protected <T extends CalculusFieldElement<T>> T semian(final AbsoluteDate date, final T day, final T altKm) {
279
280 final double f10B = inputParams.getF10B(date);
281 final double s10B = inputParams.getS10B(date);
282 final double xm10B = inputParams.getXM10B(date);
283
284 final T htz = altKm.divide(1000.0);
285
286 // COMPUTE NEW 81-DAY CENTERED SOLAR INDEX FOR FZ
287 double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;
288
289 // SEMIANNUAL AMPLITUDE
290 final T fzz = htz.multiply(FZM[3]).add(FZM[2] + FZM[4] * fsmb).multiply(htz).add(FZM[1]).multiply(fsmb).add(FZM[0]);
291
292 // COMPUTE DAILY 81-DAY CENTERED SOLAR INDEX FOR GT
293 fsmb = f10B - 0.75 * s10B - 0.37 * xm10B;
294
295 // SEMIANNUAL PHASE FUNCTION
296 final T tau = day.subtract(1).divide(365).multiply(MathUtils.TWO_PI);
297 final FieldSinCos<T> sc1P = FastMath.sinCos(tau);
298 final FieldSinCos<T> sc2P = FieldSinCos.sum(sc1P, sc1P);
299 final T gtz = sc2P.cos().multiply(GTM[9]).
300 add(sc2P.sin().multiply(GTM[8])).
301 add(sc1P.cos().multiply(GTM[7])).
302 add(sc1P.sin().multiply(GTM[6])).
303 add(GTM[5]).multiply(fsmb).
304 add( sc2P.cos().multiply(GTM[4]).
305 add(sc2P.sin().multiply(GTM[3])).
306 add(sc1P.cos().multiply(GTM[2])).
307 add(sc1P.sin().multiply(GTM[1])).
308 add(GTM[0]));
309
310 return fzz.getReal() > 1.0e-6 ? gtz.multiply(fzz) : gtz.multiply(1.0e-6);
311
312 }
313
314 /** Local provider for solar activity data. */
315 private static class LocalProvider implements JB2008InputParameters {
316
317 /** 10.7-cm Solar flux. */
318 private double f10;
319
320 /** 10.7-cm Solar Flux, averaged 81-day centered on the input time. */
321 private double f10B;
322
323 /** EUV index (26-34 nm) scaled to F10. */
324 private double s10;
325
326 /** UV 81-day averaged centered index. */
327 private double s10B;
328
329 /** MG2 index scaled to F10. */
330 private double xm10;
331
332 /** MG2 81-day ave. centered index. */
333 private double xm10B;
334
335 /** Solar X-Ray & Lya index scaled to F10. */
336 private double y10;
337
338 /** Solar X-Ray & Lya 81-day ave. centered index. */
339 private double y10B;
340
341 /** Temperature change computed from Dst index. */
342 private double dstdtc;
343
344 /** Constructor.
345 * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
346 * (Tabular time 1.0 day earlier)
347 * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
348 * (Tabular time 1.0 day earlier)
349 * @param s10 EUV index (26-34 nm) scaled to F10<br>
350 * (Tabular time 1 day earlier)
351 * @param s10B UV 81-day averaged centered index
352 * (Tabular time 1 day earlier)
353 * @param xm10 MG2 index scaled to F10<br>
354 * (Tabular time 2.0 days earlier)
355 * @param xm10B MG2 81-day ave. centered index<br>
356 * (Tabular time 2.0 days earlier)
357 * @param y10 Solar X-Ray & Lya index scaled to F10<br>
358 * (Tabular time 5.0 days earlier)
359 * @param y10B Solar X-Ray & Lya 81-day ave. centered index<br>
360 * (Tabular time 5.0 days earlier)
361 * @param dstdtc Temperature change computed from Dst index
362 */
363 LocalProvider(final double f10, final double f10B, final double s10,
364 final double s10B, final double xm10, final double xm10B,
365 final double y10, final double y10B, final double dstdtc) {
366 this.f10 = f10;
367 this.f10B = f10B;
368 this.s10 = s10;
369 this.s10B = s10B;
370 this.xm10 = xm10;
371 this.xm10B = xm10B;
372 this.y10 = y10;
373 this.y10B = y10B;
374 this.dstdtc = dstdtc;
375 }
376
377 /** {@inheritDoc} */
378 @Override
379 public AbsoluteDate getMinDate() {
380 return AbsoluteDate.PAST_INFINITY;
381 }
382
383 /** {@inheritDoc} */
384 @Override
385 public AbsoluteDate getMaxDate() {
386 return AbsoluteDate.FUTURE_INFINITY;
387 }
388
389 /** {@inheritDoc} */
390 @Override
391 public double getF10(final AbsoluteDate date) {
392 return f10;
393 }
394
395 /** {@inheritDoc} */
396 @Override
397 public double getF10B(final AbsoluteDate date) {
398 return f10B;
399 }
400
401 /** {@inheritDoc} */
402 @Override
403 public double getS10(final AbsoluteDate date) {
404 return s10;
405 }
406
407 /** {@inheritDoc} */
408 @Override
409 public double getS10B(final AbsoluteDate date) {
410 return s10B;
411 }
412
413 /** {@inheritDoc} */
414 @Override
415 public double getXM10(final AbsoluteDate date) {
416 return xm10;
417 }
418
419 /** {@inheritDoc} */
420 @Override
421 public double getXM10B(final AbsoluteDate date) {
422 return xm10B;
423 }
424
425 /** {@inheritDoc} */
426 @Override
427 public double getY10(final AbsoluteDate date) {
428 return y10;
429 }
430
431 /** {@inheritDoc} */
432 @Override
433 public double getY10B(final AbsoluteDate date) {
434 return y10B;
435 }
436
437 /** {@inheritDoc} */
438 @Override
439 public double getDSTDTC(final AbsoluteDate date) {
440 return dstdtc;
441 }
442 }
443 }