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.Field;
23 import org.hipparchus.CalculusFieldElement;
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.annotation.DefaultDataContext;
31 import org.orekit.bodies.FieldGeodeticPoint;
32 import org.orekit.bodies.GeodeticPoint;
33 import org.orekit.data.DataContext;
34 import org.orekit.frames.TopocentricFrame;
35 import org.orekit.propagation.FieldSpacecraftState;
36 import org.orekit.propagation.SpacecraftState;
37 import org.orekit.time.AbsoluteDate;
38 import org.orekit.time.DateTimeComponents;
39 import org.orekit.time.FieldAbsoluteDate;
40 import org.orekit.time.TimeScale;
41 import org.orekit.utils.Constants;
42 import org.orekit.utils.ParameterDriver;
43
44 /**
45 * Klobuchar ionospheric delay model.
46 * Klobuchar ionospheric delay model is designed as a GNSS correction model.
47 * The parameters for the model are provided by the GPS satellites in their broadcast
48 * messsage.
49 * This model is based on the assumption the electron content is concentrated
50 * in 350 km layer.
51 *
52 * The delay refers to L1 (1575.42 MHz).
53 * If the delay is sought for L2 (1227.60 MHz), multiply the result by 1.65 (Klobuchar, 1996).
54 * More generally, since ionospheric delay is inversely proportional to the square of the signal
55 * frequency f, to adapt this model to other GNSS frequencies f, multiply by (L1 / f)^2.
56 *
57 * References:
58 * ICD-GPS-200, Rev. C, (1997), pp. 125-128
59 * Klobuchar, J.A., Ionospheric time-delay algorithm for single-frequency GPS users,
60 * IEEE Transactions on Aerospace and Electronic Systems, Vol. 23, No. 3, May 1987
61 * Klobuchar, J.A., "Ionospheric Effects on GPS", Global Positioning System: Theory and
62 * Applications, 1996, pp.513-514, Parkinson, Spilker.
63 *
64 * @author Joris Olympio
65 * @since 7.1
66 *
67 */
68 public class KlobucharIonoModel implements IonosphericModel {
69
70 /** The 4 coefficients of a cubic equation representing the amplitude of the vertical delay. Units are sec/semi-circle^(i-1) for the i-th coefficient, i=1, 2, 3, 4. */
71 private final double[] alpha;
72
73 /** The 4 coefficients of a cubic equation representing the period of the model. Units are sec/semi-circle^(i-1) for the i-th coefficient, i=1, 2, 3, 4. */
74 private final double[] beta;
75
76 /** GPS time scale. */
77 private final TimeScale gps;
78
79 /** Create a new Klobuchar ionospheric delay model, when a single frequency system is used.
80 * This model accounts for at least 50 percent of RMS error due to ionospheric propagation effect (ICD-GPS-200)
81 *
82 * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
83 *
84 * @param alpha coefficients of a cubic equation representing the amplitude of the vertical delay.
85 * @param beta coefficients of a cubic equation representing the period of the model.
86 * @see #KlobucharIonoModel(double[], double[], TimeScale)
87 */
88 @DefaultDataContext
89 public KlobucharIonoModel(final double[] alpha, final double[] beta) {
90 this(alpha, beta, DataContext.getDefault().getTimeScales().getGPS());
91 }
92
93 /**
94 * Create a new Klobuchar ionospheric delay model, when a single frequency system is
95 * used. This model accounts for at least 50 percent of RMS error due to ionospheric
96 * propagation effect (ICD-GPS-200)
97 *
98 * @param alpha coefficients of a cubic equation representing the amplitude of the
99 * vertical delay.
100 * @param beta coefficients of a cubic equation representing the period of the
101 * model.
102 * @param gps GPS time scale.
103 * @since 10.1
104 */
105 public KlobucharIonoModel(final double[] alpha,
106 final double[] beta,
107 final TimeScale gps) {
108 this.alpha = alpha.clone();
109 this.beta = beta.clone();
110 this.gps = gps;
111 }
112
113 /**
114 * Calculates the ionospheric path delay for the signal path from a ground
115 * station to a satellite.
116 * <p>
117 * The path delay is computed for any elevation angle.
118 * </p>
119 * @param date current date
120 * @param geo geodetic point of receiver/station
121 * @param elevation elevation of the satellite in radians
122 * @param azimuth azimuth of the satellite in radians
123 * @param frequency frequency of the signal in Hz
124 * @param parameters ionospheric model parameters
125 * @return the path delay due to the ionosphere in m
126 */
127 public double pathDelay(final AbsoluteDate date, final GeodeticPoint geo,
128 final double elevation, final double azimuth, final double frequency,
129 final double[] parameters) {
130
131 // Sine and cosine of the azimuth
132 final SinCos sc = FastMath.sinCos(azimuth);
133
134 // degrees to semicircles
135 final double rad2semi = 1. / FastMath.PI;
136 final double semi2rad = FastMath.PI;
137
138 // Earth Centered angle
139 final double psi = 0.0137 / (elevation / FastMath.PI + 0.11) - 0.022;
140
141 // Subionospheric latitude: the latitude of the IPP (Ionospheric Pierce Point)
142 // in [-0.416, 0.416], semicircle
143 final double latIono = FastMath.min(
144 FastMath.max(geo.getLatitude() * rad2semi + psi * sc.cos(), -0.416),
145 0.416);
146
147 // Subionospheric longitude: the longitude of the IPP
148 // in semicircle
149 final double lonIono = geo.getLongitude() * rad2semi + (psi * sc.sin() / FastMath.cos(latIono * semi2rad));
150
151 // Geomagnetic latitude, semicircle
152 final double latGeom = latIono + 0.064 * FastMath.cos((lonIono - 1.617) * semi2rad);
153
154 // day of week and tow (sec)
155 // Note: Sunday=0, Monday=1, Tuesday=2, Wednesday=3, Thursday=4, Friday=5, Saturday=6
156 final DateTimeComponents dtc = date.getComponents(gps);
157 final int dofweek = dtc.getDate().getDayOfWeek();
158 final double secday = dtc.getTime().getSecondsInLocalDay();
159 final double tow = dofweek * 86400. + secday;
160
161 final double t = 43200. * lonIono + tow;
162 final double tsec = t - FastMath.floor(t / 86400.) * 86400; // Seconds of day
163
164 // Slant factor, semicircle
165 final double slantFactor = 1.0 + 16.0 * FastMath.pow(0.53 - elevation / FastMath.PI, 3);
166
167 // Period of model, seconds
168 final double period = FastMath.max(72000., beta[0] + (beta[1] + (beta[2] + beta[3] * latGeom) * latGeom) * latGeom);
169
170 // Phase of the model, radians
171 // (Max at 14.00 = 50400 sec local time)
172 final double x = 2.0 * FastMath.PI * (tsec - 50400.0) / period;
173
174 // Amplitude of the model, seconds
175 final double amplitude = FastMath.max(0, alpha[0] + (alpha[1] + (alpha[2] + alpha[3] * latGeom) * latGeom) * latGeom);
176
177 // Ionospheric correction (L1)
178 double ionoTimeDelayL1 = slantFactor * (5. * 1e-9);
179 if (FastMath.abs(x) < 1.570) {
180 ionoTimeDelayL1 += slantFactor * (amplitude * (1.0 - FastMath.pow(x, 2) / 2.0 + FastMath.pow(x, 4) / 24.0));
181 }
182
183 // Ionospheric delay for the L1 frequency, in meters, with slant correction.
184 final double ratio = FastMath.pow(1575.42e6 / frequency, 2);
185 return ratio * Constants.SPEED_OF_LIGHT * ionoTimeDelayL1;
186 }
187
188 /** {@inheritDoc} */
189 @Override
190 public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
191 final double frequency, final double[] parameters) {
192
193 // Elevation in radians
194 final Vector3D position = state.getPosition(baseFrame);
195 final double elevation = position.getDelta();
196
197 // Only consider measures above the horizon
198 if (elevation > 0.0) {
199 // Date
200 final AbsoluteDate date = state.getDate();
201 // Geodetic point
202 final GeodeticPoint geo = baseFrame.getPoint();
203 // Azimuth angle in radians
204 double azimuth = FastMath.atan2(position.getX(), position.getY());
205 if (azimuth < 0.) {
206 azimuth += MathUtils.TWO_PI;
207 }
208 // Delay
209 return pathDelay(date, geo, elevation, azimuth, frequency, parameters);
210 }
211
212 return 0.0;
213 }
214
215 /**
216 * Calculates the ionospheric path delay for the signal path from a ground
217 * station to a satellite.
218 * <p>
219 * The path delay is computed for any elevation angle.
220 * </p>
221 * @param <T> type of the elements
222 * @param date current date
223 * @param geo geodetic point of receiver/station
224 * @param elevation elevation of the satellite in radians
225 * @param azimuth azimuth of the satellite in radians
226 * @param frequency frequency of the signal in Hz
227 * @param parameters ionospheric model parameters
228 * @return the path delay due to the ionosphere in m
229 */
230 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldAbsoluteDate<T> date, final FieldGeodeticPoint<T> geo,
231 final T elevation, final T azimuth, final double frequency,
232 final T[] parameters) {
233
234 // Sine and cosine of the azimuth
235 final FieldSinCos<T> sc = FastMath.sinCos(azimuth);
236
237 // Field
238 final Field<T> field = date.getField();
239 final T zero = field.getZero();
240 final T one = field.getOne();
241
242 // degrees to semicircles
243 final T pi = one.getPi();
244 final T rad2semi = pi.reciprocal();
245
246 // Earth Centered angle
247 final T psi = elevation.divide(pi).add(0.11).divide(0.0137).reciprocal().subtract(0.022);
248
249 // Subionospheric latitude: the latitude of the IPP (Ionospheric Pierce Point)
250 // in [-0.416, 0.416], semicircle
251 final T latIono = FastMath.min(
252 FastMath.max(geo.getLatitude().multiply(rad2semi).add(psi.multiply(sc.cos())), zero.subtract(0.416)),
253 zero.newInstance(0.416));
254
255 // Subionospheric longitude: the longitude of the IPP
256 // in semicircle
257 final T lonIono = geo.getLongitude().multiply(rad2semi).add(psi.multiply(sc.sin()).divide(FastMath.cos(latIono.multiply(pi))));
258
259 // Geomagnetic latitude, semicircle
260 final T latGeom = latIono.add(FastMath.cos(lonIono.subtract(1.617).multiply(pi)).multiply(0.064));
261
262 // day of week and tow (sec)
263 // Note: Sunday=0, Monday=1, Tuesday=2, Wednesday=3, Thursday=4, Friday=5, Saturday=6
264 final DateTimeComponents dtc = date.getComponents(gps);
265 final int dofweek = dtc.getDate().getDayOfWeek();
266 final double secday = dtc.getTime().getSecondsInLocalDay();
267 final double tow = dofweek * 86400. + secday;
268
269 final T t = lonIono.multiply(43200.).add(tow);
270 final T tsec = t.subtract(FastMath.floor(t.divide(86400.)).multiply(86400.)); // Seconds of day
271
272 // Slant factor, semicircle
273 final T slantFactor = FastMath.pow(elevation.divide(pi).negate().add(0.53), 3).multiply(16.0).add(one);
274
275 // Period of model, seconds
276 final T period = FastMath.max(zero.newInstance(72000.), latGeom.multiply(latGeom.multiply(latGeom.multiply(beta[3]).add(beta[2])).add(beta[1])).add(beta[0]));
277
278 // Phase of the model, radians
279 // (Max at 14.00 = 50400 sec local time)
280 final T x = tsec.subtract(50400.0).multiply(pi.multiply(2.0)).divide(period);
281
282 // Amplitude of the model, seconds
283 final T amplitude = FastMath.max(zero, latGeom.multiply(latGeom.multiply(latGeom.multiply(alpha[3]).add(alpha[2])).add(alpha[1])).add(alpha[0]));
284
285 // Ionospheric correction (L1)
286 T ionoTimeDelayL1 = slantFactor.multiply(5. * 1e-9);
287 if (FastMath.abs(x.getReal()) < 1.570) {
288 ionoTimeDelayL1 = ionoTimeDelayL1.add(slantFactor.multiply(amplitude.multiply(one.subtract(FastMath.pow(x, 2).multiply(0.5)).add(FastMath.pow(x, 4).divide(24.0)))));
289 }
290
291 // Ionospheric delay for the L1 frequency, in meters, with slant correction.
292 final double ratio = FastMath.pow(1575.42e6 / frequency, 2);
293 return ionoTimeDelayL1.multiply(Constants.SPEED_OF_LIGHT).multiply(ratio);
294 }
295
296 /** {@inheritDoc} */
297 @Override
298 public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
299 final double frequency, final T[] parameters) {
300
301 // Elevation and azimuth in radians
302 final FieldVector3D<T> position = state.getPosition(baseFrame);
303 final T elevation = position.getDelta();
304
305 if (elevation.getReal() > 0.0) {
306 // Date
307 final FieldAbsoluteDate<T> date = state.getDate();
308 // Geodetic point
309 final FieldGeodeticPoint<T> geo = baseFrame.getPoint(date.getField());
310 // Azimuth angle in radians
311 T azimuth = FastMath.atan2(position.getX(), position.getY());
312 if (azimuth.getReal() < 0.) {
313 azimuth = azimuth.add(MathUtils.TWO_PI);
314 }
315 // Delay
316 return pathDelay(date, geo, elevation, azimuth, frequency, parameters);
317 }
318
319 return elevation.getField().getZero();
320 }
321
322 @Override
323 public List<ParameterDriver> getParametersDrivers() {
324 return Collections.emptyList();
325 }
326 }