1 /* Copyright 2002-2026 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 org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.util.FastMath;
21 import org.hipparchus.util.FieldSinCos;
22 import org.orekit.bodies.FieldGeodeticPoint;
23
24 /** Container for ray-periapsis parameters.
25 * <p>By convention, point 1 is at lower height.</p>
26 * @author Bryan Cazabonne
27 * @since 13.0
28 */
29 class FieldRay<T extends CalculusFieldElement<T>> {
30
31 /** Threshold for ray-periapsis parameters computation. */
32 private static final double THRESHOLD = 1.0e-10;
33
34 /** Receiver altitude [m].
35 * @since 13.0
36 */
37 private final T recH;
38
39 /** Satellite altitude [m].
40 * @since 13.0
41 */
42 private final T satH;
43
44 /** Distance of the first point from the ray periapsis [m]. */
45 private final T s1;
46
47 /** Distance of the second point from the ray periapsis [m]. */
48 private final T s2;
49
50 /** Ray-periapsis radius [m]. */
51 private final T rp;
52
53 /** Ray-periapsis latitude [rad]. */
54 private final T latP;
55
56 /** Ray-periapsis longitude [rad]. */
57 private final T lonP;
58
59 /** Sine and cosine of ray-periapsis latitude. */
60 private final FieldSinCos<T> scLatP;
61
62 /** Sine of azimuth of satellite as seen from ray-periapsis. */
63 private final T sinAzP;
64
65 /** Cosine of azimuth of satellite as seen from ray-periapsis. */
66 private final T cosAzP;
67
68 /**
69 * Constructor.
70 *
71 * @param recP receiver position
72 * @param satP satellite position
73 */
74 FieldRay(final FieldGeodeticPoint<T> recP, final FieldGeodeticPoint<T> satP) {
75
76 // Integration limits in meters (Eq. 140 and 141)
77 this.recH = recP.getAltitude();
78 this.satH = satP.getAltitude();
79 final T r1 = recH.add(NeQuickModel.RE);
80 final T r2 = satH.add(NeQuickModel.RE);
81
82 // Useful parameters
83 final T pi = r1.getPi();
84 final T lat1 = recP.getLatitude();
85 final T lat2 = satP.getLatitude();
86 final T lon1 = recP.getLongitude();
87 final T lon2 = satP.getLongitude();
88 final FieldSinCos<T> scLatSat = FastMath.sinCos(lat2);
89 final FieldSinCos<T> scLatRec = FastMath.sinCos(lat1);
90 final FieldSinCos<T> scLon21 = FastMath.sinCos(lon2.subtract(lon1));
91
92 // Zenith angle computation, using:
93 // - Eq. 153 with added protection against numerical noise near zenith observation
94 // - replacing Eq. 154 by a different one more stable near zenith
95 // - replacing Eq. 155 by different ones avoiding trigonometric functions
96 final T cosD = FastMath.min(r1.getField().getOne(),
97 scLatRec.sin().multiply(scLatSat.sin()).
98 add(scLatRec.cos().multiply(scLatSat.cos()).multiply(scLon21.cos())));
99 final T sinDSinM = scLatSat.cos().multiply(scLon21.sin());
100 final T sinDCosM = scLatSat.sin().multiply(scLatRec.cos()).subtract(scLatSat.cos().multiply(scLatRec.sin()).multiply(scLon21.cos()));
101 final T sinD = FastMath.sqrt(sinDSinM.multiply(sinDSinM).add(sinDCosM.multiply(sinDCosM)));
102 final T u = r2.multiply(sinD);
103 final T v = r2.multiply(cosD).subtract(r1);
104 final T inv = FastMath.sqrt(u.multiply(u).add(v.multiply(v))).reciprocal();
105 final T sinZ = u.multiply(inv);
106 final T cosZ = v.multiply(inv);
107
108 // Ray-periapsis computation in meters (Eq. 156)
109 this.rp = r1.multiply(sinZ);
110
111 // Ray-periapsis latitude and longitude
112 if (FastMath.abs(FastMath.abs(lat1).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD) {
113
114 // Ray-periapsis latitude (Eq. 157)
115 this.latP = FastMath.copySign(FastMath.atan2(u, v), lat1);
116
117 // Ray-periapsis longitude (Eq. 164)
118 this.lonP = lon2.add(pi);
119
120 } else if (FastMath.abs(sinZ.getReal()) < THRESHOLD) {
121 // satellite is almost on receiver zenith
122
123 this.latP = recP.getLatitude();
124 this.lonP = recP.getLongitude();
125
126 } else {
127
128 // Ray-periapsis latitude (Eq. 158 to 163)
129 final T sinAz = FastMath.sin(lon2.subtract(lon1)).multiply(scLatSat.cos()).divide(sinD);
130 final T cosAz = scLatSat.sin().subtract(cosD.multiply(scLatRec.sin())).
131 divide(sinD.multiply(scLatRec.cos()));
132 final T sinLatP = scLatRec.sin().multiply(sinZ).
133 subtract(scLatRec.cos().multiply(cosZ).multiply(cosAz));
134 final T cosLatP = FastMath.sqrt(sinLatP.multiply(sinLatP).negate().add(1.0));
135 this.latP = FastMath.atan2(sinLatP, cosLatP);
136
137 // Ray-periapsis longitude (Eq. 165 to 167, plus protection against ray-periapsis along polar axis)
138 if (cosLatP.getReal() < THRESHOLD) {
139 this.lonP = cosLatP.getField().getZero();
140 } else {
141 final T sinLonP = sinAz.negate().multiply(cosZ).divide(cosLatP);
142 final T cosLonP = sinZ.subtract(scLatRec.sin().multiply(sinLatP)).
143 divide(scLatRec.cos().multiply(cosLatP));
144 this.lonP = FastMath.atan2(sinLonP, cosLonP).add(lon1);
145 }
146
147 }
148
149 // Sine and cosine of ray-periapsis latitude
150 this.scLatP = FastMath.sinCos(latP);
151
152 if (FastMath.abs(FastMath.abs(latP).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD || FastMath.abs(
153 sinZ.getReal()) < THRESHOLD) {
154 // Eq. 172 and 173
155 this.sinAzP = pi.getField().getZero();
156 this.cosAzP = FastMath.copySign(pi.getField().getOne(), latP).negate();
157 } else {
158 final FieldSinCos<T> scLon = FastMath.sinCos(lon2.subtract(lonP));
159 // Sine and cosine of azimuth of satellite as seen from ray-periapsis
160 final FieldSinCos<T> scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon));
161 // Eq. 174 and 175
162 this.sinAzP = scLatSat.cos().multiply(scLon.sin()).divide(scPsi.sin());
163 this.cosAzP = scLatSat.sin().subtract(scLatP.sin().multiply(scPsi.cos())).
164 divide(scLatP.cos().multiply(scPsi.sin()));
165 }
166
167 // Integration end points s1 and s2 in meters (Eq. 176 and 177)
168 this.s1 = FastMath.sqrt(r1.multiply(r1).subtract(rp.multiply(rp)));
169 this.s2 = FastMath.sqrt(r2.multiply(r2).subtract(rp.multiply(rp)));
170 }
171
172 /**
173 * Get receiver altitude.
174 * @return receiver altitude
175 * @since 13.0
176 */
177 public T getRecH() {
178 return recH;
179 }
180
181 /**
182 * Get satellite altitude.
183 * @return satellite altitude
184 * @since 13.0
185 */
186 public T getSatH() {
187 return satH;
188 }
189
190 /**
191 * Get the distance of the first point from the ray periapsis.
192 *
193 * @return s1 in meters
194 */
195 public T getS1() {
196 return s1;
197 }
198
199 /**
200 * Get the distance of the second point from the ray periapsis.
201 *
202 * @return s2 in meters
203 */
204 public T getS2() {
205 return s2;
206 }
207
208 /**
209 * Get the ray-periapsis radius.
210 *
211 * @return the ray-periapsis radius in meters
212 */
213 public T getRadius() {
214 return rp;
215 }
216
217 /**
218 * Get the ray-periapsis latitude.
219 *
220 * @return the ray-periapsis latitude in radians
221 */
222 public T getLatitude() {
223 return latP;
224 }
225
226 /**
227 * Get the ray-periapsis latitude sin/cos.
228 *
229 * @return the ray-periapsis latitude sin/cos
230 * @since 13.0
231 */
232 public FieldSinCos<T> getScLat() {
233 return scLatP;
234 }
235
236 /**
237 * Get the ray-periapsis longitude.
238 *
239 * @return the ray-periapsis longitude in radians
240 */
241 public T getLongitude() {
242 return lonP;
243 }
244
245 /**
246 * Get the sine of azimuth of satellite as seen from ray-periapsis.
247 *
248 * @return the sine of azimuth
249 */
250 public T getSineAz() {
251 return sinAzP;
252 }
253
254 /**
255 * Get the cosine of azimuth of satellite as seen from ray-periapsis.
256 *
257 * @return the cosine of azimuth
258 */
259 public T getCosineAz() {
260 return cosAzP;
261 }
262
263 /**
264 * Compute the great circle angle from ray-periapsis to satellite.
265 * <p>
266 * This method used the equations 168 to 171 of the reference document.
267 * </p>
268 *
269 * @param scLat sine and cosine of satellite latitude
270 * @param scLon sine and cosine of satellite longitude minus receiver longitude
271 * @return the great circle angle in radians
272 */
273 private T greatCircleAngle(final FieldSinCos<T> scLat, final FieldSinCos<T> scLon) {
274 if (FastMath.abs(FastMath.abs(latP).getReal() - 0.5 * FastMath.PI) < THRESHOLD) {
275 return FastMath.abs(FastMath.asin(scLat.sin()).subtract(latP));
276 } else {
277 final T cosPhi = scLatP.sin().multiply(scLat.sin()).
278 add(scLatP.cos().multiply(scLat.cos()).multiply(scLon.cos()));
279 final T sinPhi = FastMath.sqrt(cosPhi.multiply(cosPhi).negate().add(1.0));
280 return FastMath.atan2(sinPhi, cosPhi);
281 }
282 }
283
284 }