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;
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.geometry.euclidean.threed.FieldVector3D;
25 import org.hipparchus.geometry.euclidean.threed.Vector3D;
26 import org.hipparchus.util.FastMath;
27 import org.hipparchus.util.FieldSinCos;
28 import org.hipparchus.util.MathUtils;
29 import org.hipparchus.util.SinCos;
30 import org.orekit.bodies.FieldGeodeticPoint;
31 import org.orekit.bodies.GeodeticPoint;
32 import org.orekit.frames.FieldStaticTransform;
33 import org.orekit.frames.StaticTransform;
34 import org.orekit.frames.TopocentricFrame;
35 import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201;
36 import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Data;
37 import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Header;
38 import org.orekit.propagation.FieldSpacecraftState;
39 import org.orekit.propagation.SpacecraftState;
40 import org.orekit.time.AbsoluteDate;
41 import org.orekit.time.FieldAbsoluteDate;
42 import org.orekit.utils.Constants;
43 import org.orekit.utils.FieldLegendrePolynomials;
44 import org.orekit.utils.LegendrePolynomials;
45 import org.orekit.utils.ParameterDriver;
46
47 /**
48 * Ionospheric model based on SSR IM201 message.
49 * <p>
50 * Within this message, the ionospheric VTEC is provided
51 * using spherical harmonic expansions. For a given ionospheric
52 * layer, the slant TEC value is calculated using the satellite
53 * elevation and the height of the corresponding layer. The
54 * total slant TEC is computed by the sum of the individual slant
55 * TEC for each layer.
56 * </p>
57 * @author Bryan Cazabonne
58 * @since 11.0
59 * @see "IGS State Space Representation (SSR) Format, Version 1.00, October 2020."
60 */
61 public class SsrVtecIonosphericModel implements IonosphericModel, IonosphericDelayModel {
62
63 /** Earth radius in meters (see reference). */
64 private static final double EARTH_RADIUS = 6370000.0;
65
66 /** Multiplication factor for path delay computation. */
67 private static final double FACTOR = 40.3e16;
68
69 /** SSR Ionosphere VTEC Spherical Harmonics Message.. */
70 private final transient SsrIm201 vtecMessage;
71
72 /**
73 * Constructor.
74 * @param vtecMessage SSR Ionosphere VTEC Spherical Harmonics Message.
75 */
76 public SsrVtecIonosphericModel(final SsrIm201 vtecMessage) {
77 this.vtecMessage = vtecMessage;
78 }
79
80 /** {@inheritDoc} */
81 @Override
82 public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
83 final double frequency, final double[] parameters) {
84 return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
85 }
86
87 /** {@inheritDoc} */
88 @Override
89 public double pathDelay(final SpacecraftState state,
90 final TopocentricFrame baseFrame, final AbsoluteDate receptionDate,
91 final double frequency, final double[] parameters) {
92
93 // we use transform from base frame to inert frame and invert it
94 // because base frame could be peered with inertial frame (hence improving performances with caching)
95 // but the reverse is almost never used
96 final StaticTransform base2Inert = baseFrame.getStaticTransformTo(state.getFrame(), receptionDate);
97 final Vector3D position = base2Inert.getInverse().transformPosition(state.getPosition());
98
99 // Elevation in radians
100 final double elevation = position.getDelta();
101
102 // Only consider measures above the horizon
103 if (elevation > 0.0) {
104
105 // Azimuth angle in radians
106 double azimuth = FastMath.atan2(position.getX(), position.getY());
107 if (azimuth < 0.) {
108 azimuth += MathUtils.TWO_PI;
109 }
110
111 // Initialize slant TEC
112 double stec = 0.0;
113
114 // Message header
115 final SsrIm201Header header = vtecMessage.getHeader();
116
117 // Loop on ionospheric layers
118 for (final SsrIm201Data data : vtecMessage.getData()) {
119 stec += stecIonosphericLayer(data, header, elevation, azimuth, baseFrame.getPoint());
120 }
121
122 // Return the path delay
123 return FACTOR * stec / (frequency * frequency);
124
125 }
126
127 // Delay is equal to 0.0
128 return 0.0;
129
130 }
131
132 /** {@inheritDoc} */
133 @Override
134 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
135 final double frequency, final T[] parameters) {
136 return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
137 }
138
139 /** {@inheritDoc} */
140 @Override
141 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
142 final TopocentricFrame baseFrame, final FieldAbsoluteDate<T> receptionDate,
143 final double frequency, final T[] parameters) {
144
145 // Field
146 final Field<T> field = state.getDate().getField();
147
148 // we use transform from base frame to inert frame and invert it
149 // because base frame could be peered with inertial frame (hence improving performances with caching)
150 // but the reverse is almost never used
151 final FieldStaticTransform<T> base2Inert = baseFrame.getStaticTransformTo(state.getFrame(), receptionDate);
152 final FieldVector3D<T> position = base2Inert.getInverse().transformPosition(state.getPosition());
153
154 // Elevation in radians
155 final T elevation = position.getDelta();
156
157 // Only consider measures above the horizon
158 if (elevation.getReal() > 0.0) {
159
160 // Azimuth angle in radians
161 T azimuth = FastMath.atan2(position.getX(), position.getY());
162 if (azimuth.getReal() < 0.) {
163 azimuth = azimuth.add(MathUtils.TWO_PI);
164 }
165
166 // Initialize slant TEC
167 T stec = field.getZero();
168
169 // Message header
170 final SsrIm201Header header = vtecMessage.getHeader();
171
172 // Loop on ionospheric layers
173 for (SsrIm201Data data : vtecMessage.getData()) {
174 stec = stec.add(stecIonosphericLayer(data, header, elevation, azimuth, baseFrame.getPoint(field)));
175 }
176
177 // Return the path delay
178 return stec.multiply(FACTOR).divide(frequency * frequency);
179
180 }
181
182 // Delay is equal to 0.0
183 return field.getZero();
184
185 }
186
187 /** {@inheritDoc} */
188 @Override
189 public List<ParameterDriver> getParametersDrivers() {
190 return Collections.emptyList();
191 }
192
193 /**
194 * Calculates the slant TEC for a given ionospheric layer.
195 * @param im201Data ionospheric data for the current layer
196 * @param im201Header container for data contained in the header
197 * @param elevation satellite elevation angle [rad]
198 * @param azimuth satellite azimuth angle [rad]
199 * @param point geodetic point
200 * @return the slant TEC for the current ionospheric layer
201 */
202 private static double stecIonosphericLayer(final SsrIm201Data im201Data, final SsrIm201Header im201Header,
203 final double elevation, final double azimuth,
204 final GeodeticPoint point) {
205
206 // Geodetic point data
207 final double phiR = point.getLatitude();
208 final double lambdaR = point.getLongitude();
209 final double hR = point.getAltitude();
210
211 // Data contained in the message
212 final double hI = im201Data.getHeightIonosphericLayer();
213 final int degree = im201Data.getSphericalHarmonicsDegree();
214 final int order = im201Data.getSphericalHarmonicsOrder();
215 final double[][] cnm = im201Data.getCnm();
216 final double[][] snm = im201Data.getSnm();
217
218 // Spherical Earth's central angle
219 final double psiPP = calculatePsi(hR, hI, elevation);
220
221 // Sine and cosine of useful angles
222 final SinCos scA = FastMath.sinCos(azimuth);
223 final SinCos scPhiR = FastMath.sinCos(phiR);
224 final SinCos scPsi = FastMath.sinCos(psiPP);
225
226 // Pierce point latitude and longitude
227 final double phiPP = calculatePiercePointLatitude(scPhiR, scPsi, scA);
228 final double lambdaPP = calculatePiercePointLongitude(scA, phiPP, psiPP, phiR, lambdaR);
229
230 // Mean sun fixed longitude (modulo 2pi)
231 final double lambdaS = calculateSunLongitude(im201Header, lambdaPP);
232
233 // VTEC
234 // According to the documentation, negative VTEC values must be ignored and shall be replaced by 0.0
235 final double vtec = FastMath.max(0.0, calculateVTEC(degree, order, cnm, snm, phiPP, lambdaS));
236
237 // Return STEC for the current ionospheric layer
238 return vtec / FastMath.sin(elevation + psiPP);
239
240 }
241
242 /**
243 * Calculates the slant TEC for a given ionospheric layer.
244 * @param im201Data ionospheric data for the current layer
245 * @param im201Header container for data contained in the header
246 * @param elevation satellite elevation angle [rad]
247 * @param azimuth satellite azimuth angle [rad]
248 * @param point geodetic point
249 * @param <T> type of the elements
250 * @return the slant TEC for the current ionospheric layer
251 */
252 private static <T extends CalculusFieldElement<T>> T stecIonosphericLayer(final SsrIm201Data im201Data, final SsrIm201Header im201Header,
253 final T elevation, final T azimuth,
254 final FieldGeodeticPoint<T> point) {
255
256 // Geodetic point data
257 final T phiR = point.getLatitude();
258 final T lambdaR = point.getLongitude();
259 final T hR = point.getAltitude();
260
261 // Data contained in the message
262 final double hI = im201Data.getHeightIonosphericLayer();
263 final int degree = im201Data.getSphericalHarmonicsDegree();
264 final int order = im201Data.getSphericalHarmonicsOrder();
265 final double[][] cnm = im201Data.getCnm();
266 final double[][] snm = im201Data.getSnm();
267
268 // Spherical Earth's central angle
269 final T psiPP = calculatePsi(hR, hI, elevation);
270
271 // Sine and cosine of useful angles
272 final FieldSinCos<T> scA = FastMath.sinCos(azimuth);
273 final FieldSinCos<T> scPhiR = FastMath.sinCos(phiR);
274 final FieldSinCos<T> scPsi = FastMath.sinCos(psiPP);
275
276 // Pierce point latitude and longitude
277 final T phiPP = calculatePiercePointLatitude(scPhiR, scPsi, scA);
278 final T lambdaPP = calculatePiercePointLongitude(scA, phiPP, psiPP, phiR, lambdaR);
279
280 // Mean sun fixed longitude (modulo 2pi)
281 final T lambdaS = calculateSunLongitude(im201Header, lambdaPP);
282
283 // VTEC
284 // According to the documentation, negative VTEC values must be ignored and shall be replaced by 0.0
285 final T vtec = FastMath.max(phiR.getField().getZero(), calculateVTEC(degree, order, cnm, snm, phiPP, lambdaS));
286
287 // Return STEC for the current ionospheric layer
288 return vtec.divide(FastMath.sin(elevation.add(psiPP)));
289
290 }
291
292 /**
293 * Calculates the spherical Earth’s central angle between station position and
294 * the projection of the pierce point to the spherical Earth’s surface.
295 * @param hR height of station position in meters
296 * @param hI height of ionospheric layer in meters
297 * @param elevation satellite elevation angle in radians
298 * @return the spherical Earth’s central angle in radians
299 */
300 private static double calculatePsi(final double hR, final double hI,
301 final double elevation) {
302 final double ratio = (EARTH_RADIUS + hR) / (EARTH_RADIUS + hI);
303 return MathUtils.SEMI_PI - elevation - FastMath.asin(ratio * FastMath.cos(elevation));
304 }
305
306 /**
307 * Calculates the spherical Earth’s central angle between station position and
308 * the projection of the pierce point to the spherical Earth’s surface.
309 * @param hR height of station position in meters
310 * @param hI height of ionospheric layer in meters
311 * @param elevation satellite elevation angle in radians
312 * @param <T> type of the elements
313 * @return the spherical Earth’s central angle in radians
314 */
315 private static <T extends CalculusFieldElement<T>> T calculatePsi(final T hR, final double hI,
316 final T elevation) {
317 final T ratio = hR.add(EARTH_RADIUS).divide(EARTH_RADIUS + hI);
318 return hR.getPi().multiply(0.5).subtract(elevation).subtract(FastMath.asin(ratio.multiply(FastMath.cos(elevation))));
319 }
320
321 /**
322 * Calculates the latitude of the pierce point in the spherical Earth model.
323 * @param scPhiR sine and cosine of the geocentric latitude of the station
324 * @param scPsi sine and cosine of the spherical Earth's central angle
325 * @param scA sine and cosine of the azimuth angle
326 * @return the latitude of the pierce point in the spherical Earth model in radians
327 */
328 private static double calculatePiercePointLatitude(final SinCos scPhiR, final SinCos scPsi, final SinCos scA) {
329 return FastMath.asin(scPhiR.sin() * scPsi.cos() + scPhiR.cos() * scPsi.sin() * scA.cos());
330 }
331
332 /**
333 * Calculates the latitude of the pierce point in the spherical Earth model.
334 * @param scPhiR sine and cosine of the geocentric latitude of the station
335 * @param scPsi sine and cosine of the spherical Earth's central angle
336 * @param scA sine and cosine of the azimuth angle
337 * @param <T> type of the elements
338 * @return the latitude of the pierce point in the spherical Earth model in radians
339 */
340 private static <T extends CalculusFieldElement<T>> T calculatePiercePointLatitude(final FieldSinCos<T> scPhiR,
341 final FieldSinCos<T> scPsi,
342 final FieldSinCos<T> scA) {
343 return FastMath.asin(scPhiR.sin().multiply(scPsi.cos()).add(scPhiR.cos().multiply(scPsi.sin()).multiply(scA.cos())));
344 }
345
346 /**
347 * Calculates the longitude of the pierce point in the spherical Earth model.
348 * @param scA sine and cosine of the azimuth angle
349 * @param phiPP the latitude of the pierce point in the spherical Earth model in radians
350 * @param psiPP the spherical Earth’s central angle in radians
351 * @param phiR the geocentric latitude of the station in radians
352 * @param lambdaR the geocentric longitude of the station
353 * @return the longitude of the pierce point in the spherical Earth model in radians
354 */
355 private static double calculatePiercePointLongitude(final SinCos scA,
356 final double phiPP, final double psiPP,
357 final double phiR, final double lambdaR) {
358
359 // arcSin(sin(PsiPP) * sin(Azimuth) / cos(PhiPP))
360 final double arcSin = FastMath.asin(FastMath.sin(psiPP) * scA.sin() / FastMath.cos(phiPP));
361
362 // Return
363 return verifyCondition(scA.cos(), psiPP, phiR) ? lambdaR + FastMath.PI - arcSin : lambdaR + arcSin;
364
365 }
366
367 /**
368 * Calculates the longitude of the pierce point in the spherical Earth model.
369 * @param scA sine and cosine of the azimuth angle
370 * @param phiPP the latitude of the pierce point in the spherical Earth model in radians
371 * @param psiPP the spherical Earth’s central angle in radians
372 * @param phiR the geocentric latitude of the station in radians
373 * @param lambdaR the geocentric longitude of the station
374 * @param <T> type of the elements
375 * @return the longitude of the pierce point in the spherical Earth model in radians
376 */
377 private static <T extends CalculusFieldElement<T>> T calculatePiercePointLongitude(final FieldSinCos<T> scA,
378 final T phiPP, final T psiPP,
379 final T phiR, final T lambdaR) {
380
381 // arcSin(sin(PsiPP) * sin(Azimuth) / cos(PhiPP))
382 final T arcSin = FastMath.asin(FastMath.sin(psiPP).multiply(scA.sin()).divide(FastMath.cos(phiPP)));
383
384 // Return
385 return verifyCondition(scA.cos().getReal(), psiPP.getReal(), phiR.getReal()) ?
386 lambdaR.add(arcSin.getPi()).subtract(arcSin) : lambdaR.add(arcSin);
387
388 }
389
390 /**
391 * Calculate the mean sun fixed longitude phase.
392 * @param im201Header header of the IM201 message
393 * @param lambdaPP the longitude of the pierce point in the spherical Earth model in radians
394 * @return the mean sun fixed longitude phase in radians
395 */
396 private static double calculateSunLongitude(final SsrIm201Header im201Header, final double lambdaPP) {
397 final double t = getTime(im201Header);
398 return MathUtils.normalizeAngle(lambdaPP + (t - 50400.0) * FastMath.PI / 43200.0, FastMath.PI);
399 }
400
401 /**
402 * Calculate the mean sun fixed longitude phase.
403 * @param im201Header header of the IM201 message
404 * @param lambdaPP the longitude of the pierce point in the spherical Earth model in radians
405 * @param <T> type of the elements
406 * @return the mean sun fixed longitude phase in radians
407 */
408 private static <T extends CalculusFieldElement<T>> T calculateSunLongitude(final SsrIm201Header im201Header, final T lambdaPP) {
409 final double t = getTime(im201Header);
410 return MathUtils.normalizeAngle(lambdaPP.add(lambdaPP.getPi().multiply(t - 50400.0).divide(43200.0)), lambdaPP.getPi());
411 }
412
413 /**
414 * Calculate the VTEC contribution for a given ionospheric layer.
415 * @param degree degree of spherical expansion
416 * @param order order of spherical expansion
417 * @param cnm cosine coefficients for the layer in TECU
418 * @param snm sine coefficients for the layer in TECU
419 * @param phiPP geocentric latitude of ionospheric pierce point for the layer in radians
420 * @param lambdaS mean sun fixed and phase shifted longitude of ionospheric pierce point
421 * @return the VTEC contribution for the current ionospheric layer in TECU
422 */
423 private static double calculateVTEC(final int degree, final int order,
424 final double[][] cnm, final double[][] snm,
425 final double phiPP, final double lambdaS) {
426
427 // Initialize VTEC value
428 double vtec = 0.0;
429
430 // Compute Legendre Polynomials Pnm(sin(phiPP))
431 final LegendrePolynomials p = new LegendrePolynomials(degree, order, FastMath.sin(phiPP));
432
433 // Calculate VTEC
434 for (int n = 0; n <= degree; n++) {
435
436 for (int m = 0; m <= FastMath.min(n, order); m++) {
437
438 // Legendre coefficients
439 final SinCos sc = FastMath.sinCos(m * lambdaS);
440 final double pCosmLambda = p.getPnm(n, m) * sc.cos();
441 final double pSinmLambda = p.getPnm(n, m) * sc.sin();
442
443 // Update VTEC value
444 vtec += cnm[n][m] * pCosmLambda + snm[n][m] * pSinmLambda;
445
446 }
447
448 }
449
450 // Return the VTEC
451 return vtec;
452
453 }
454
455 /**
456 * Calculate the VTEC contribution for a given ionospheric layer.
457 * @param degree degree of spherical expansion
458 * @param order order of spherical expansion
459 * @param cnm cosine coefficients for the layer in TECU
460 * @param snm sine coefficients for the layer in TECU
461 * @param phiPP geocentric latitude of ionospheric pierce point for the layer in radians
462 * @param lambdaS mean sun fixed and phase shifted longitude of ionospheric pierce point
463 * @param <T> type of the elements
464 * @return the VTEC contribution for the current ionospheric layer in TECU
465 */
466 private static <T extends CalculusFieldElement<T>> T calculateVTEC(final int degree, final int order,
467 final double[][] cnm, final double[][] snm,
468 final T phiPP, final T lambdaS) {
469
470 // Initialize VTEC value
471 T vtec = phiPP.getField().getZero();
472
473 // Compute Legendre Polynomials Pnm(sin(phiPP))
474 final FieldLegendrePolynomials<T> p = new FieldLegendrePolynomials<>(degree, order, FastMath.sin(phiPP));
475
476 // Calculate VTEC
477 for (int n = 0; n <= degree; n++) {
478
479 for (int m = 0; m <= FastMath.min(n, order); m++) {
480
481 // Legendre coefficients
482 final FieldSinCos<T> sc = FastMath.sinCos(lambdaS.multiply(m));
483 final T pCosmLambda = p.getPnm(n, m).multiply(sc.cos());
484 final T pSinmLambda = p.getPnm(n, m).multiply(sc.sin());
485
486 // Update VTEC value
487 vtec = vtec.add(pCosmLambda.multiply(cnm[n][m]).add(pSinmLambda.multiply(snm[n][m])));
488
489 }
490
491 }
492
493 // Return the VTEC
494 return vtec;
495
496 }
497
498 /**
499 * Get the SSR epoch time of computation modulo 86400 seconds.
500 * @param im201Header header data
501 * @return the SSR epoch time of computation modulo 86400 seconds
502 */
503 private static double getTime(final SsrIm201Header im201Header) {
504 final double ssrEpochTime = im201Header.getSsrEpoch1s();
505 return ssrEpochTime - FastMath.floor(ssrEpochTime / Constants.JULIAN_DAY) * Constants.JULIAN_DAY;
506 }
507
508 /**
509 * Verify the condition for the calculation of the pierce point longitude.
510 * @param scACos cosine of the azimuth angle
511 * @param psiPP the spherical Earth’s central angle in radians
512 * @param phiR the geocentric latitude of the station in radians
513 * @return true if the condition is respected
514 */
515 private static boolean verifyCondition(final double scACos, final double psiPP,
516 final double phiR) {
517
518 // tan(PsiPP) * cos(Azimuth)
519 final double tanPsiCosA = FastMath.tan(psiPP) * scACos;
520
521 // Verify condition
522 return phiR >= 0 && tanPsiCosA > FastMath.tan(MathUtils.SEMI_PI - phiR) ||
523 phiR < 0 && -tanPsiCosA > FastMath.tan(MathUtils.SEMI_PI + phiR);
524
525 }
526
527 }