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 (Eq. 153 to 155)
91 // with added protection against numerical noise near zenith observation
92 final double cosD = FastMath.min(1.0,
93 scLatRec.sin() * scLatSat.sin() +
94 scLatRec.cos() * scLatSat.cos() * scLon21.cos());
95 final double sinD = FastMath.sqrt(1.0 - cosD * cosD);
96 final double z = FastMath.atan2(sinD, cosD - (r1 / r2));
97 final SinCos scZ = FastMath.sinCos(z);
98
99 // Ray-perigee computation in meters (Eq. 156)
100 this.rp = r1 * scZ.sin();
101
102 // Ray-perigee latitude and longitude
103 if (FastMath.abs(FastMath.abs(lat1) - 0.5 * FastMath.PI) < THRESHOLD) {
104 // receiver is almost at North or South pole
105
106 // Ray-perigee latitude (Eq. 157)
107 this.latP = FastMath.copySign(z, lat1);
108
109 // Ray-perigee longitude (Eq. 164)
110 if (z < 0) {
111 this.lonP = lon2;
112 } else {
113 this.lonP = lon2 + FastMath.PI;
114 }
115
116 } else if (FastMath.abs(scZ.sin()) < THRESHOLD) {
117 // satellite is almost on receiver zenith
118
119 this.latP = recP.getLatitude();
120 this.lonP = recP.getLongitude();
121
122 } else {
123
124 // Ray-perigee latitude (Eq. 158 to 163)
125 final double sinAz = scLon21.sin() * scLatSat.cos() / sinD;
126 final double cosAz = (scLatSat.sin() - cosD * scLatRec.sin()) / (sinD * scLatRec.cos());
127 final double sinLatP = scLatRec.sin() * scZ.sin() - scLatRec.cos() * scZ.cos() * cosAz;
128 final double cosLatP = FastMath.sqrt(1.0 - sinLatP * sinLatP);
129 this.latP = FastMath.atan2(sinLatP, cosLatP);
130
131 // Ray-perigee longitude (Eq. 165 to 167)
132 final double sinLonP = -sinAz * scZ.cos() / cosLatP;
133 final double cosLonP = (scZ.sin() - scLatRec.sin() * sinLatP) / (scLatRec.cos() * cosLatP);
134 this.lonP = FastMath.atan2(sinLonP, cosLonP) + lon1;
135
136 }
137
138 // Sine and cosine of ray-perigee latitude
139 this.scLatP = FastMath.sinCos(latP);
140
141 if (FastMath.abs(FastMath.abs(latP) - 0.5 * FastMath.PI) < THRESHOLD || FastMath.abs(scZ.sin()) < THRESHOLD) {
142 // Eq. 172 and 173
143 this.sinAzP = 0.0;
144 this.cosAzP = -FastMath.copySign(1, latP);
145 } else {
146 final SinCos scLon = FastMath.sinCos(lon2 - lonP);
147 // Sine and cosine of azimuth of satellite as seen from ray-perigee
148 final SinCos scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon));
149 // Eq. 174 and 175
150 this.sinAzP = scLatSat.cos() * scLon.sin() / scPsi.sin();
151 this.cosAzP = (scLatSat.sin() - scLatP.sin() * scPsi.cos()) / (scLatP.cos() * scPsi.sin());
152 }
153
154 // Integration end points s1 and s2 in meters (Eq. 176 and 177)
155 this.s1 = FastMath.sqrt(r1 * r1 - rp * rp);
156 this.s2 = FastMath.sqrt(r2 * r2 - rp * rp);
157 }
158
159 /**
160 * Get receiver altitude.
161 * @return receiver altitude
162 * @since 13.0
163 */
164 public double getRecH() {
165 return recH;
166 }
167
168 /**
169 * Get satellite altitude.
170 * @return satellite altitude
171 * @since 13.0
172 */
173 public double getSatH() {
174 return satH;
175 }
176
177 /**
178 * Get the distance of the first point from the ray perigee.
179 *
180 * @return s1 in meters
181 */
182 public double getS1() {
183 return s1;
184 }
185
186 /**
187 * Get the distance of the second point from the ray perigee.
188 *
189 * @return s2 in meters
190 */
191 public double getS2() {
192 return s2;
193 }
194
195 /**
196 * Get the ray-perigee radius.
197 *
198 * @return the ray-perigee radius in meters
199 */
200 public double getRadius() {
201 return rp;
202 }
203
204 /**
205 * Get the ray-perigee latitude.
206 *
207 * @return the ray-perigee latitude in radians
208 */
209 public double getLatitude() {
210 return latP;
211 }
212
213 /**
214 * Get the ray-perigee latitude sin/cos.
215 *
216 * @return the ray-perigee latitude sin/cos
217 * @since 13.0
218 */
219 public SinCos getScLat() {
220 return scLatP;
221 }
222
223 /**
224 * Get the ray-perigee longitude.
225 *
226 * @return the ray-perigee longitude in radians
227 */
228 public double getLongitude() {
229 return lonP;
230 }
231
232 /**
233 * Get the sine of azimuth of satellite as seen from ray-perigee.
234 *
235 * @return the sine of azimuth
236 */
237 public double getSineAz() {
238 return sinAzP;
239 }
240
241 /**
242 * Get the cosine of azimuth of satellite as seen from ray-perigee.
243 *
244 * @return the cosine of azimuth
245 */
246 public double getCosineAz() {
247 return cosAzP;
248 }
249
250 /**
251 * Compute the great circle angle from ray-perigee to satellite.
252 * <p>
253 * This method used the equations 168 to 171 of the reference document.
254 * </p>
255 *
256 * @param scLat sine and cosine of satellite latitude
257 * @param scLon sine and cosine of satellite longitude minus receiver longitude
258 * @return the great circle angle in radians
259 */
260 private double greatCircleAngle(final SinCos scLat, final SinCos scLon) {
261 if (FastMath.abs(FastMath.abs(latP) - 0.5 * FastMath.PI) < THRESHOLD) {
262 return FastMath.abs(FastMath.asin(scLat.sin()) - latP);
263 } else {
264 final double cosPhi = scLatP.sin() * scLat.sin() + scLatP.cos() * scLat.cos() * scLon.cos();
265 final double sinPhi = FastMath.sqrt(1.0 - cosPhi * cosPhi);
266 return FastMath.atan2(sinPhi, cosPhi);
267 }
268 }
269
270 }