1 /* Copyright 2013-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.rugged.refraction;
18
19 import java.util.ArrayList;
20 import java.util.Collections;
21 import java.util.List;
22
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.util.FastMath;
25 import org.orekit.bodies.GeodeticPoint;
26 import org.orekit.rugged.errors.RuggedException;
27 import org.orekit.rugged.errors.RuggedMessages;
28 import org.orekit.rugged.intersection.IntersectionAlgorithm;
29 import org.orekit.rugged.utils.ExtendedEllipsoid;
30 import org.orekit.rugged.utils.NormalizedGeodeticPoint;
31
32 /**
33 * Atmospheric refraction model based on multiple layers with associated refractive index.
34 * @author Sergio Esteves
35 * @author Luc Maisonobe
36 * @author Guylaine Prat
37 * @since 2.0
38 */
39 public class MultiLayerModel extends AtmosphericRefraction {
40
41 /** Observed body ellipsoid. */
42 private final ExtendedEllipsoid ellipsoid;
43
44 /** Constant refraction layers. */
45 private final List<ConstantRefractionLayer> refractionLayers;
46
47 /** Atmosphere lowest altitude (m). */
48 private final double atmosphereLowestAltitude;
49
50 /** Simple constructor.
51 * <p>
52 * This model uses a built-in set of layers.
53 * </p>
54 * @param ellipsoid the ellipsoid to be used.
55 */
56 public MultiLayerModel(final ExtendedEllipsoid ellipsoid) {
57
58 super();
59
60 this.ellipsoid = ellipsoid;
61
62 this.refractionLayers = new ArrayList<>(15);
63 this.refractionLayers.add(new ConstantRefractionLayer(100000.00, 1.000000));
64 this.refractionLayers.add(new ConstantRefractionLayer( 50000.00, 1.000000));
65 this.refractionLayers.add(new ConstantRefractionLayer( 40000.00, 1.000001));
66 this.refractionLayers.add(new ConstantRefractionLayer( 30000.00, 1.000004));
67 this.refractionLayers.add(new ConstantRefractionLayer( 23000.00, 1.000012));
68 this.refractionLayers.add(new ConstantRefractionLayer( 18000.00, 1.000028));
69 this.refractionLayers.add(new ConstantRefractionLayer( 14000.00, 1.000052));
70 this.refractionLayers.add(new ConstantRefractionLayer( 11000.00, 1.000083));
71 this.refractionLayers.add(new ConstantRefractionLayer( 9000.00, 1.000106));
72 this.refractionLayers.add(new ConstantRefractionLayer( 7000.00, 1.000134));
73 this.refractionLayers.add(new ConstantRefractionLayer( 5000.00, 1.000167));
74 this.refractionLayers.add(new ConstantRefractionLayer( 3000.00, 1.000206));
75 this.refractionLayers.add(new ConstantRefractionLayer( 1000.00, 1.000252));
76 this.refractionLayers.add(new ConstantRefractionLayer( 0.00, 1.000278));
77 this.refractionLayers.add(new ConstantRefractionLayer( -1000.00, 1.000306));
78
79 // get the lowest altitude of the atmospheric model
80 this.atmosphereLowestAltitude = refractionLayers.get(refractionLayers.size() - 1).getLowestAltitude();
81 }
82
83 /** Simple constructor.
84 * @param ellipsoid the ellipsoid to be used.
85 * @param refractionLayers the refraction layers to be used with this model (layers can be in any order).
86 */
87 public MultiLayerModel(final ExtendedEllipsoid ellipsoid, final List<ConstantRefractionLayer> refractionLayers) {
88
89 super();
90 // at this stage no optimization is set up: no optimization grid is defined
91
92 this.ellipsoid = ellipsoid;
93 this.refractionLayers = new ArrayList<>(refractionLayers);
94
95 // sort the layers from the highest (index = 0) to the lowest (index = size - 1)
96 Collections.sort(this.refractionLayers,
97 (l1, l2) -> Double.compare(l2.getLowestAltitude(), l1.getLowestAltitude()));
98
99 // get the lowest altitude of the model
100 atmosphereLowestAltitude = this.refractionLayers.get(this.refractionLayers.size() - 1).getLowestAltitude();
101 }
102
103 /** Compute the (position, LOS) of the intersection with the lowest atmospheric layer.
104 * @param satPos satellite position, in body frame
105 * @param satLos satellite line of sight, in body frame
106 * @param rawIntersection intersection point without refraction correction
107 * @return the intersection position and LOS with the lowest atmospheric layer
108 * @since 2.1
109 */
110 private IntersectionLOS computeToLowestAtmosphereLayer(
111 final Vector3D satPos, final Vector3D satLos,
112 final NormalizedGeodeticPoint rawIntersection) {
113
114 if (rawIntersection.getAltitude() < atmosphereLowestAltitude) {
115 throw new RuggedException(RuggedMessages.NO_LAYER_DATA, rawIntersection.getAltitude(),
116 atmosphereLowestAltitude);
117 }
118
119 Vector3D pos = satPos;
120 Vector3D los = satLos.normalize();
121
122 // Compute the intersection point with the lowest layer of atmosphere
123 double n1 = -1;
124 GeodeticPoint gp = ellipsoid.transform(satPos, ellipsoid.getBodyFrame(), null);
125
126 // Perform the exact computation (no optimization)
127 // TBN: the refractionLayers is ordered from the highest to the lowest
128 for (ConstantRefractionLayer refractionLayer : refractionLayers) {
129
130 if (refractionLayer.getLowestAltitude() > gp.getAltitude()) {
131 continue;
132 }
133
134 final double n2 = refractionLayer.getRefractiveIndex();
135
136 if (n1 > 0) {
137
138 // when we get here, we have already performed one iteration in the loop
139 // so gp is the los intersection with the layers interface (it was a
140 // point on ground at loop initialization, but is overridden at each iteration)
141
142 // get new los by applying Snell's law at atmosphere layers interfaces
143 // we avoid computing sequences of inverse-trigo/trigo/inverse-trigo functions
144 // we just use linear algebra and square roots, it is faster and more accurate
145
146 // at interface crossing, the interface normal is z, the local zenith direction
147 // the ray direction (i.e. los) is u in the upper layer and v in the lower layer
148 // v is in the (u, zenith) plane, so we can say
149 // (1) v = α u + β z
150 // with α>0 as u and v are roughly in the same direction as the ray is slightly bent
151
152 // let θ₁ be the los incidence angle at interface crossing
153 // θ₁ = π - angle(u, zenith) is between 0 and π/2 for a downwards observation
154 // let θ₂ be the exit angle at interface crossing
155 // from Snell's law, we have n₁ sin θ₁ = n₂ sin θ₂ and θ₂ is also between 0 and π/2
156 // we have:
157 // (2) u·z = -cos θ₁
158 // (3) v·z = -cos θ₂
159 // combining equations (1), (2) and (3) and remembering z·z = 1 as z is normalized , we get
160 // (4) β = α cos θ₁ - cos θ₂
161 // with all the expressions above, we can rewrite the fact v is normalized:
162 // 1 = v·v
163 // = α² u·u + 2αβ u·z + β² z·z
164 // = α² - 2αβ cos θ₁ + β²
165 // = α² - 2α² cos² θ₁ + 2 α cos θ₁ cos θ₂ + α² cos² θ₁ - 2 α cos θ₁ cos θ₂ + cos² θ₂
166 // = α²(1 - cos² θ₁) + cos² θ₂
167 // hence α² = (1 - cos² θ₂)/(1 - cos² θ₁)
168 // = sin² θ₂ / sin² θ₁
169 // as α is positive, and both θ₁ and θ₂ are between 0 and π/2, we finally get
170 // α = sin θ₂ / sin θ₁
171 // (5) α = n₁/n₂
172 // the α coefficient is independent from the incidence angle,
173 // it depends only on the ratio of refractive indices!
174 //
175 // back to equation (4) and using again the fact θ₂ is between 0 and π/2, we can now write
176 // β = α cos θ₁ - cos θ₂
177 // = n₁/n₂ cos θ₁ - cos θ₂
178 // = n₁/n₂ cos θ₁ - √(1 - sin² θ₂)
179 // = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² sin² θ₁)
180 // = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² (1 - cos² θ₁))
181 // = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² + (n₁/n₂)² cos² θ₁)
182 // (6) β = -k - √(k² - ζ)
183 // where ζ = (n₁/n₂)² - 1 and k = n₁/n₂ u·z, which is negative, and close to -1 for
184 // nadir observations. As we expect atmosphere models to have small transitions between
185 // layers, we have to handle accurately the case where n₁/n₂ ≈ 1 so ζ ≈ 0. In this case,
186 // a cancellation occurs inside the square root: √(k² - ζ) ≈ √k² ≈ -k (because k is negative).
187 // So β ≈ -k + k ≈ 0 and another cancellation occurs, outside of the square root.
188 // This means that despite equation (6) is mathematically correct, it is prone to numerical
189 // errors when consecutive layers have close refractive indices. A better equivalent
190 // expression is needed. The fact β is close to 0 in this case was expected because
191 // equation (1) reads v = α u + β z, and α = n₁/n₂, so when n₁/n₂ ≈ 1, we have
192 // α ≈ 1 and β ≈ 0, so v ≈ u: when two layers have similar refractive indices, the
193 // propagation direction is almost unchanged.
194 //
195 // The first step for the accurate computation of β is to compute ζ = (n₁/n₂)² - 1
196 // accurately and avoid a cancellation just after a division (which is the least accurate
197 // of the four operations) and a squaring. We will simply use:
198 // ζ = (n₁/n₂)² - 1
199 // = (n₁ - n₂) (n₁ + n₂) / n₂²
200 // The cancellation is still there, but it occurs in the subtraction n₁ - n₂, which are
201 // the most accurate values we can get.
202 // The second step for the accurate computation of β is to rewrite equation (6)
203 // by both multiplying and dividing by the dual factor -k + √(k² - ζ):
204 // β = -k - √(k² - ζ)
205 // = (-k - √(k² - ζ)) * (-k + √(k² - ζ)) / (-k + √(k² - ζ))
206 // = (k² - (k² - ζ)) / (-k + √(k² - ζ))
207 // (7) β = ζ / (-k + √(k² - ζ))
208 // expression (7) is more stable numerically than expression (6), because when ζ ≈ 0
209 // its denominator is about -2k, there are no cancellation anymore after the square root.
210 // β is computed with the same accuracy as ζ
211 final double alpha = n1 / n2;
212 final double k = alpha * Vector3D.dotProduct(los, gp.getZenith());
213 final double zeta = (n1 - n2) * (n1 + n2) / (n2 * n2);
214 final double beta = zeta / (FastMath.sqrt(k * k - zeta) - k);
215 los = new Vector3D(alpha, los, beta, gp.getZenith());
216 }
217
218 // In case the altitude of the intersection without atmospheric refraction
219 // is above the lowest altitude of the atmosphere: stop the search
220 if (rawIntersection.getAltitude() > refractionLayer.getLowestAltitude()) {
221 break;
222 }
223
224 // Get for the intersection point: the position and the LOS
225 pos = ellipsoid.pointAtAltitude(pos, los, refractionLayer.getLowestAltitude());
226 gp = ellipsoid.transform(pos, ellipsoid.getBodyFrame(), null);
227
228 n1 = n2;
229 }
230 return new IntersectionLOS(pos, los);
231 }
232
233 /** {@inheritDoc} */
234 @Override
235 public NormalizedGeodeticPoint applyCorrection(final Vector3D satPos, final Vector3D satLos,
236 final NormalizedGeodeticPoint rawIntersection,
237 final IntersectionAlgorithm algorithm) {
238
239 final IntersectionLOS intersectionLOS = computeToLowestAtmosphereLayer(satPos, satLos, rawIntersection);
240 final Vector3D pos = intersectionLOS.getIntersectionPos();
241 final Vector3D los = intersectionLOS.getIntersectionLos();
242
243 // at this stage the pos belongs to the lowest atmospheric layer.
244 // We can compute the intersection of line of sight (los) with Digital Elevation Model
245 // as usual (without atmospheric refraction).
246 return algorithm.refineIntersection(ellipsoid, pos, los, rawIntersection);
247 }
248
249 } // end of class MultiLayerModel
250
251 /** Container for the (position, LOS) of the intersection with the lowest atmospheric layer.
252 * @author Guylaine Prat
253 * @since 2.1
254 */
255 class IntersectionLOS {
256
257 /** Position of the intersection with the lowest atmospheric layer. */
258 private Vector3D intersectionPos;
259 /** LOS of the intersection with the lowest atmospheric layer. */
260 private Vector3D intersectionLos;
261
262 /** Default constructor.
263 * @param intersectionPos position of the intersection
264 * @param intersectionLos los of the intersection
265 */
266 IntersectionLOS(final Vector3D intersectionPos, final Vector3D intersectionLos) {
267 this.intersectionPos = intersectionPos;
268 this.intersectionLos = intersectionLos;
269 }
270 public Vector3D getIntersectionPos() {
271 return intersectionPos;
272 }
273
274 public Vector3D getIntersectionLos() {
275 return intersectionLos;
276 }
277 }