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