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