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.ionosphere.nequick;
18
19 import java.util.Collections;
20 import java.util.List;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.Field;
24 import org.hipparchus.util.FastMath;
25 import org.hipparchus.util.MathArrays;
26 import org.orekit.bodies.BodyShape;
27 import org.orekit.bodies.FieldGeodeticPoint;
28 import org.orekit.bodies.GeodeticPoint;
29 import org.orekit.frames.TopocentricFrame;
30 import org.orekit.models.earth.ionosphere.IonosphericModel;
31 import org.orekit.propagation.FieldSpacecraftState;
32 import org.orekit.propagation.SpacecraftState;
33 import org.orekit.time.AbsoluteDate;
34 import org.orekit.time.DateComponents;
35 import org.orekit.time.DateTimeComponents;
36 import org.orekit.time.FieldAbsoluteDate;
37 import org.orekit.time.TimeScale;
38 import org.orekit.utils.ParameterDriver;
39 import org.orekit.utils.units.Unit;
40
41 /**
42 * NeQuick ionospheric delay model.
43 *
44 * @author Bryan Cazabonne
45 *
46 * @see "European Union (2016). European GNSS (Galileo) Open Service-Ionospheric Correction
47 * Algorithm for Galileo Single Frequency Users. 1.2."
48 * @see <a href="https://www.itu.int/rec/R-REC-P.531/en">ITU-R P.531</a>
49 *
50 * @since 10.1
51 */
52 public abstract class NeQuickModel implements IonosphericModel {
53
54 /** Mean Earth radius in m (Ref Table 2.5.2). */
55 public static final double RE = 6371200.0;
56
57 /** Factor for the electron density computation. */
58 private static final double DENSITY_FACTOR = 1.0e11;
59
60 /** Factor for the path delay computation. */
61 private static final double DELAY_FACTOR = 40.3e16;
62
63 /** F2 coefficients used by the F2 layer (flatten array for cache efficiency). */
64 private final double[][] flattenF2;
65
66 /** Fm3 coefficients used by the M(3000)F2 layer(flatten array for cache efficiency). */
67 private final double[][] flattenFm3;
68
69 /** UTC time scale. */
70 private final TimeScale utc;
71
72 /** Simple constructor.
73 * @param utc UTC time scale
74 * @since 13.0
75 */
76 protected NeQuickModel(final TimeScale utc) {
77
78 this.utc = utc;
79
80 // F2 layer values
81 this.flattenF2 = new double[12][];
82 this.flattenFm3 = new double[12][];
83
84 }
85
86 /** Get UTC time scale.
87 * @return UTC time scale
88 * @since 13.0
89 */
90 public TimeScale getUtc() {
91 return utc;
92 }
93
94 /** {@inheritDoc} */
95 @Override
96 public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
97 final double frequency, final double[] parameters) {
98 // Point
99 final GeodeticPoint recPoint = baseFrame.getPoint();
100 // Date
101 final AbsoluteDate date = state.getDate();
102
103 // Reference body shape
104 final BodyShape ellipsoid = baseFrame.getParentShape();
105 // Satellite geodetic coordinates
106 final GeodeticPoint satPoint = ellipsoid.transform(state.getPosition(ellipsoid.getBodyFrame()),
107 ellipsoid.getBodyFrame(), state.getDate());
108
109 // Total Electron Content
110 final double tec = stec(date, recPoint, satPoint);
111
112 // Ionospheric delay
113 final double factor = DELAY_FACTOR / (frequency * frequency);
114 return factor * tec;
115 }
116
117 /** {@inheritDoc} */
118 @Override
119 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
120 final TopocentricFrame baseFrame,
121 final double frequency,
122 final T[] parameters) {
123 // Date
124 final FieldAbsoluteDate<T> date = state.getDate();
125 // Point
126 final FieldGeodeticPoint<T> recPoint = baseFrame.getPoint(date.getField());
127
128
129 // Reference body shape
130 final BodyShape ellipsoid = baseFrame.getParentShape();
131 // Satellite geodetic coordinates
132 final FieldGeodeticPoint<T> satPoint = ellipsoid.transform(state.getPosition(ellipsoid.getBodyFrame()), ellipsoid.getBodyFrame(), state.getDate());
133
134 // Total Electron Content
135 final T tec = stec(date, recPoint, satPoint);
136
137 // Ionospheric delay
138 final double factor = DELAY_FACTOR / (frequency * frequency);
139 return tec.multiply(factor);
140 }
141
142 /** {@inheritDoc} */
143 @Override
144 public List<ParameterDriver> getParametersDrivers() {
145 return Collections.emptyList();
146 }
147
148 /**
149 * This method allows the computation of the Slant Total Electron Content (STEC).
150 * @param date current date
151 * @param recP receiver position
152 * @param satP satellite position
153 * @return the STEC in TECUnits
154 */
155 public double stec(final AbsoluteDate date, final GeodeticPoint recP, final GeodeticPoint satP) {
156 return stec(date.getComponents(utc), new Ray(recP, satP));
157 }
158
159 /**
160 * This method allows the computation of the Slant Total Electron Content (STEC).
161 * @param <T> type of the elements
162 * @param date current date
163 * @param recP receiver position
164 * @param satP satellite position
165 * @return the STEC in TECUnits
166 */
167 public <T extends CalculusFieldElement<T>> T stec(final FieldAbsoluteDate<T> date,
168 final FieldGeodeticPoint<T> recP,
169 final FieldGeodeticPoint<T> satP) {
170 return stec(date.getComponents(utc), new FieldRay<>(recP, satP));
171 }
172
173 /** Compute modip for a location.
174 * @param latitude latitude
175 * @param longitude longitude
176 * @return modip at specified location
177 * @since 13.0
178 */
179 protected abstract double computeMODIP(double latitude, double longitude);
180
181 /** Compute modip for a location.
182 * @param <T> type of the field elements
183 * @param latitude latitude
184 * @param longitude longitude
185 * @return modip at specified location
186 * @since 13.0
187 */
188 protected abstract <T extends CalculusFieldElement<T>> T computeMODIP(T latitude, T longitude);
189
190 /**
191 * Compute Fourier time series.
192 * @param dateTime current date time components
193 * @param az effective ionisation level
194 * @return Fourier time series
195 * @since 13.0.1
196 */
197 public FourierTimeSeries computeFourierTimeSeries(final DateTimeComponents dateTime, final double az) {
198
199 // Load the correct CCIR file
200 loadsIfNeeded(dateTime.getDate());
201
202 return new FourierTimeSeries(dateTime, az,
203 flattenF2[dateTime.getDate().getMonth() - 1],
204 flattenFm3[dateTime.getDate().getMonth() - 1]);
205
206 }
207
208 /**
209 * Computes the electron density at a given height.
210 * @param dateTime date
211 * @param az effective ionization level
212 * @param latitude latitude along the integration path
213 * @param longitude longitude along the integration path
214 * @param h height along the integration path in m
215 * @return electron density [m⁻³]
216 * @since 13.0
217 */
218 public double electronDensity(final DateTimeComponents dateTime, final double az,
219 final double latitude, final double longitude, final double h) {
220 return electronDensity(computeFourierTimeSeries(dateTime, az), latitude, longitude, h);
221 }
222
223 /**
224 * Computes the electron density at a given height.
225 * @param fourierTimeSeries Fourier time series for foF2 and M(3000)F2 layer (flatten array)
226 * @param latitude latitude along the integration path
227 * @param longitude longitude along the integration path
228 * @param h height along the integration path in m
229 * @return electron density [m⁻³]
230 * @since 13.0.1
231 */
232 public double electronDensity(final FourierTimeSeries fourierTimeSeries,
233 final double latitude, final double longitude, final double h) {
234
235 final double modip = computeMODIP(latitude, longitude);
236 final NeQuickParameters parameters = new NeQuickParameters(fourierTimeSeries, latitude, longitude, modip);
237
238 // Convert height in kilometers
239 final double hInKm = Unit.KILOMETRE.fromSI(h);
240
241 // Electron density
242 final double n;
243 if (hInKm <= parameters.getHmF2()) {
244 n = bottomElectronDensity(hInKm, parameters);
245 } else {
246 n = topElectronDensity(hInKm, parameters);
247 }
248
249 return n;
250
251 }
252
253 /**
254 * Compute Fourier time series.
255 * @param <T> type of the elements
256 * @param dateTime current date time components
257 * @param az effective ionisation level
258 * @return Fourier time series
259 * @since 13.0.1
260 */
261 public <T extends CalculusFieldElement<T>> FieldFourierTimeSeries<T> computeFourierTimeSeries(final DateTimeComponents dateTime,
262 final T az) {
263
264 // Load the correct CCIR file
265 loadsIfNeeded(dateTime.getDate());
266
267 return new FieldFourierTimeSeries<>(dateTime, az,
268 flattenF2[dateTime.getDate().getMonth() - 1],
269 flattenFm3[dateTime.getDate().getMonth() - 1]);
270
271 }
272
273 /**
274 * Computes the electron density at a given height.
275 * @param <T> type of the elements
276 * @param dateTime date
277 * @param az effective ionization level
278 * @param latitude latitude along the integration path
279 * @param longitude longitude along the integration path
280 * @param h height along the integration path in m
281 * @return electron density [m⁻³]
282 * @since 13.0
283 * CalculusFieldElement, CalculusFieldElement, CalculusFieldElement)}
284 */
285 public <T extends CalculusFieldElement<T>> T electronDensity(final DateTimeComponents dateTime, final T az,
286 final T latitude, final T longitude, final T h) {
287 return electronDensity(computeFourierTimeSeries(dateTime, az), latitude, longitude, h);
288 }
289
290 /**
291 * Computes the electron density at a given height.
292 * @param <T> type of the elements
293 * @param fourierTimeSeries Fourier time series for foF2 and M(3000)F2 layer (flatten array)
294 * @param latitude latitude along the integration path
295 * @param longitude longitude along the integration path
296 * @param h height along the integration path in m
297 * @return electron density [m⁻³]
298 * @since 13.0.1
299 */
300 public <T extends CalculusFieldElement<T>> T electronDensity(final FieldFourierTimeSeries<T> fourierTimeSeries,
301 final T latitude, final T longitude, final T h) {
302
303 final T modip = computeMODIP(latitude, longitude);
304 final FieldNeQuickParameters<T> parameters = new FieldNeQuickParameters<>(fourierTimeSeries, latitude, longitude, modip);
305
306 // Convert height in kilometers
307 final T hInKm = Unit.KILOMETRE.fromSI(h);
308
309 // Electron density
310 final T n;
311 if (hInKm.getReal() <= parameters.getHmF2().getReal()) {
312 n = bottomElectronDensity(hInKm, parameters);
313 } else {
314 n = topElectronDensity(hInKm, parameters);
315 }
316
317 return n;
318
319 }
320
321 /**
322 * Computes the electron density of the bottomside.
323 * @param h height in km
324 * @param parameters NeQuick model parameters
325 * @return the electron density N in m⁻³
326 */
327 private double bottomElectronDensity(final double h, final NeQuickParameters parameters) {
328
329 // Select the relevant B parameter for the current height (Eq. 109 and 110)
330 final double be = parameters.getBE(h);
331 final double bf1 = parameters.getBF1(h);
332 final double bf2 = parameters.getB2Bot();
333
334 // Useful array of constants
335 final double[] ct = new double[] {
336 1.0 / bf2, 1.0 / bf1, 1.0 / be
337 };
338
339 // Compute the exponential argument for each layer (Eq. 111 to 113)
340 final double[] arguments = computeExponentialArguments(h, parameters);
341
342 // S coefficients
343 final double[] s = new double[3];
344 // Array of corrective terms
345 final double[] ds = new double[3];
346
347 // Layer amplitudes
348 final double[] amplitudes = parameters.getLayerAmplitudes();
349
350 // Fill arrays (Eq. 114 to 118)
351 for (int i = 0; i < 3; i++) {
352 if (FastMath.abs(arguments[i]) > 25.0) {
353 s[i] = 0.0;
354 ds[i] = 0.0;
355 } else {
356 final double expA = clipExp(arguments[i]);
357 final double opExpA = 1.0 + expA;
358 s[i] = amplitudes[i] * (expA / (opExpA * opExpA));
359 ds[i] = ct[i] * ((1.0 - expA) / (1.0 + expA));
360 }
361 }
362
363 // Electron density
364 final double aNo = s[0] + s[1] + s[2];
365 if (applyChapmanParameters(h)) {
366 // Chapman parameters (Eq. 119 and 120)
367 final double bc = 1.0 - 10.0 * (MathArrays.linearCombination(s[0], ds[0], s[1], ds[1], s[2], ds[2]) / aNo);
368 final double z = 0.1 * (h - 100.0);
369 // Electron density (Eq. 121)
370 return aNo * clipExp(1.0 - bc * z - clipExp(-z)) * DENSITY_FACTOR;
371 } else {
372 return aNo * DENSITY_FACTOR;
373 }
374
375 }
376
377 /**
378 * Computes the electron density of the bottomside.
379 * @param <T> type of the elements
380 * @param h height in km
381 * @param parameters NeQuick model parameters
382 * @return the electron density N in m⁻³
383 */
384 private <T extends CalculusFieldElement<T>> T bottomElectronDensity(final T h,
385 final FieldNeQuickParameters<T> parameters) {
386
387 // Zero and One
388 final Field<T> field = h.getField();
389 final T zero = field.getZero();
390 final T one = field.getOne();
391
392 // Select the relevant B parameter for the current height (Eq. 109 and 110)
393 final T be = parameters.getBE(h);
394 final T bf1 = parameters.getBF1(h);
395 final T bf2 = parameters.getB2Bot();
396
397 // Useful array of constants
398 final T[] ct = MathArrays.buildArray(field, 3);
399 ct[0] = bf2.reciprocal();
400 ct[1] = bf1.reciprocal();
401 ct[2] = be.reciprocal();
402
403 // Compute the exponential argument for each layer (Eq. 111 to 113)
404 final T[] arguments = computeExponentialArguments(h, parameters);
405
406 // S coefficients
407 final T[] s = MathArrays.buildArray(field, 3);
408 // Array of corrective terms
409 final T[] ds = MathArrays.buildArray(field, 3);
410
411 // Layer amplitudes
412 final T[] amplitudes = parameters.getLayerAmplitudes();
413
414 // Fill arrays (Eq. 114 to 118)
415 for (int i = 0; i < 3; i++) {
416 if (FastMath.abs(arguments[i]).getReal() > 25.0) {
417 s[i] = zero;
418 ds[i] = zero;
419 } else {
420 final T expA = clipExp(arguments[i]);
421 final T opExpA = expA.add(1.0);
422 s[i] = amplitudes[i].multiply(expA.divide(opExpA.multiply(opExpA)));
423 ds[i] = ct[i].multiply(expA.negate().add(1.0).divide(expA.add(1.0)));
424 }
425 }
426
427 // Electron density
428 final T aNo = s[0].add(s[1]).add(s[2]);
429 if (applyChapmanParameters(h.getReal())) {
430 // Chapman parameters (Eq. 119 and 120)
431 final T bc = one.linearCombination(s[0], ds[0], s[1], ds[1], s[2], ds[2]).divide(aNo).multiply(10.0).negate().add(1.0);
432 final T z = h.subtract(100.0).multiply(0.1);
433 // Electron density (Eq. 121)
434 return aNo.multiply(clipExp(bc.multiply(z).add(clipExp(z.negate())).negate().add(1.0))).multiply(DENSITY_FACTOR);
435 } else {
436 return aNo.multiply(DENSITY_FACTOR);
437 }
438
439 }
440
441 /**
442 * Computes the electron density of the topside.
443 * @param h height in km
444 * @param parameters NeQuick model parameters
445 * @return the electron density N in m⁻³
446 */
447 private double topElectronDensity(final double h, final NeQuickParameters parameters) {
448
449 // Constant parameters (Eq. 122 and 123)
450 final double g = 0.125;
451 final double r = 100.0;
452
453 // Arguments deltaH and z (Eq. 124 and 125)
454 final double deltaH = h - parameters.getHmF2();
455 final double h0 = computeH0(parameters);
456 final double z = deltaH / (h0 * (1.0 + (r * g * deltaH) / (r * h0 + g * deltaH)));
457
458 // Exponential (Eq. 126)
459 final double ee = clipExp(z);
460
461 // Electron density (Eq. 127)
462 if (ee > 1.0e11) {
463 return (4.0 * parameters.getNmF2() / ee) * DENSITY_FACTOR;
464 } else {
465 final double opExpZ = 1.0 + ee;
466 return ((4.0 * parameters.getNmF2() * ee) / (opExpZ * opExpZ)) * DENSITY_FACTOR;
467 }
468 }
469
470 /**
471 * Computes the electron density of the topside.
472 * @param <T> type of the elements
473 * @param h height in km
474 * @param parameters NeQuick model parameters
475 * @return the electron density N in m⁻³
476 */
477 private <T extends CalculusFieldElement<T>> T topElectronDensity(final T h,
478 final FieldNeQuickParameters<T> parameters) {
479
480 // Constant parameters (Eq. 122 and 123)
481 final double g = 0.125;
482 final double r = 100.0;
483
484 // Arguments deltaH and z (Eq. 124 and 125)
485 final T deltaH = h.subtract(parameters.getHmF2());
486 final T h0 = computeH0(parameters);
487 final T z = deltaH.divide(h0.multiply(deltaH.multiply(r).multiply(g).divide(h0.multiply(r).add(deltaH.multiply(g))).add(1.0)));
488
489 // Exponential (Eq. 126)
490 final T ee = clipExp(z);
491
492 // Electron density (Eq. 127)
493 if (ee.getReal() > 1.0e11) {
494 return parameters.getNmF2().multiply(4.0).divide(ee).multiply(DENSITY_FACTOR);
495 } else {
496 final T opExpZ = ee.add(1.0);
497 return parameters.getNmF2().multiply(4.0).multiply(ee).divide(opExpZ.multiply(opExpZ)).multiply(DENSITY_FACTOR);
498 }
499 }
500
501 /**
502 * Lazy loading of CCIR data.
503 * @param date current date components
504 */
505 private void loadsIfNeeded(final DateComponents date) {
506
507 // Month index
508 final int monthIndex = date.getMonth() - 1;
509
510 // Check if CCIR has already been loaded for this month
511 if (flattenF2[monthIndex] == null) {
512
513 // Read file
514 final CCIRLoader loader = new CCIRLoader();
515 loader.loadCCIRCoefficients(date);
516
517 // Update arrays
518 this.flattenF2[monthIndex] = flatten(loader.getF2());
519 this.flattenFm3[monthIndex] = flatten(loader.getFm3());
520 }
521 }
522
523 /** Flatten a 3-dimensions array.
524 * <p>
525 * This method convert 3-dimensions arrays into 1-dimension arrays
526 * optimized to avoid cache misses when looping over all elements.
527 * </p>
528 * @param original original array a[i][j][k]
529 * @return flatten array, for embedded loops on j (outer), k (intermediate), i (inner)
530 */
531 private double[] flatten(final double[][][] original) {
532 final double[] flatten = new double[original.length * original[0].length * original[0][0].length];
533 int index = 0;
534 for (int j = 0; j < original[0].length; j++) {
535 for (int k = 0; k < original[0][0].length; k++) {
536 for (final double[][] doubles : original) {
537 flatten[index++] = doubles[j][k];
538 }
539 }
540 }
541 return flatten;
542 }
543
544 /**
545 * A clipped exponential function.
546 * <p>
547 * This function, describe in section F.2.12.2 of the reference document, is
548 * recommended for the computation of exponential values.
549 * </p>
550 * @param power power for exponential function
551 * @return clipped exponential value
552 */
553 protected double clipExp(final double power) {
554 if (power > 80.0) {
555 return 5.5406E34;
556 } else if (power < -80) {
557 return 1.8049E-35;
558 } else {
559 return FastMath.exp(power);
560 }
561 }
562
563 /**
564 * A clipped exponential function.
565 * <p>
566 * This function, describe in section F.2.12.2 of the reference document, is
567 * recommended for the computation of exponential values.
568 * </p>
569 * @param <T> type of the elements
570 * @param power power for exponential function
571 * @return clipped exponential value
572 */
573 protected <T extends CalculusFieldElement<T>> T clipExp(final T power) {
574 if (power.getReal() > 80.0) {
575 return power.newInstance(5.5406E34);
576 } else if (power.getReal() < -80) {
577 return power.newInstance(1.8049E-35);
578 } else {
579 return FastMath.exp(power);
580 }
581 }
582
583 /**
584 * This method allows the computation of the Slant Total Electron Content (STEC).
585 *
586 * @param dateTime current date
587 * @param ray ray-perigee parameters
588 * @return the STEC in TECUnits
589 */
590 abstract double stec(DateTimeComponents dateTime, Ray ray);
591
592 /**
593 * This method allows the computation of the Slant Total Electron Content (STEC).
594 *
595 * @param <T> type of the field elements
596 * @param dateTime current date
597 * @param ray ray-perigee parameters
598 * @return the STEC in TECUnits
599 */
600 abstract <T extends CalculusFieldElement<T>> T stec(DateTimeComponents dateTime, FieldRay<T> ray);
601
602 /**
603 * Check if Chapman parameters should be applied.
604 *
605 * @param hInKm height in kilometers
606 * @return true if Chapman parameters should be applied
607 * @since 13.0
608 */
609 abstract boolean applyChapmanParameters(double hInKm);
610
611 /**
612 * Compute exponential arguments.
613 * @param h height in km
614 * @param parameters NeQuick model parameters
615 * @return exponential arguments
616 * @since 13.0
617 */
618 abstract double[] computeExponentialArguments(double h, NeQuickParameters parameters);
619
620 /**
621 * Compute exponential arguments.
622 * @param <T> type of the field elements
623 * @param h height in km
624 * @param parameters NeQuick model parameters
625 * @return exponential arguments
626 * @since 13.0
627 */
628 abstract <T extends CalculusFieldElement<T>> T[] computeExponentialArguments(T h,
629 FieldNeQuickParameters<T> parameters);
630
631 /**
632 * Compute topside thickness parameter.
633 * @param parameters NeQuick model parameters
634 * @return topside thickness parameter
635 * @since 13.0
636 */
637 abstract double computeH0(NeQuickParameters parameters);
638
639 /**
640 * Compute topside thickness parameter.
641 * @param <T> type of the field elements
642 * @param parameters NeQuick model parameters
643 * @return topside thickness parameter
644 * @since 13.0
645 */
646 abstract <T extends CalculusFieldElement<T>> T computeH0(FieldNeQuickParameters<T> parameters);
647
648 }