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