MultiLayerModel.java
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.rugged.refraction;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.errors.RuggedMessages;
import org.orekit.rugged.intersection.IntersectionAlgorithm;
import org.orekit.rugged.utils.ExtendedEllipsoid;
import org.orekit.rugged.utils.NormalizedGeodeticPoint;
/**
* Atmospheric refraction model based on multiple layers with associated refractive index.
* @author Sergio Esteves
* @author Luc Maisonobe
* @author Guylaine Prat
* @since 2.0
*/
public class MultiLayerModel extends AtmosphericRefraction {
/** Observed body ellipsoid. */
private final ExtendedEllipsoid ellipsoid;
/** Constant refraction layers. */
private final List<ConstantRefractionLayer> refractionLayers;
/** Atmosphere lowest altitude (m). */
private final double atmosphereLowestAltitude;
/** Simple constructor.
* <p>
* This model uses a built-in set of layers.
* </p>
* @param ellipsoid the ellipsoid to be used.
*/
public MultiLayerModel(final ExtendedEllipsoid ellipsoid) {
super();
this.ellipsoid = ellipsoid;
this.refractionLayers = new ArrayList<>(15);
this.refractionLayers.add(new ConstantRefractionLayer(100000.00, 1.000000));
this.refractionLayers.add(new ConstantRefractionLayer( 50000.00, 1.000000));
this.refractionLayers.add(new ConstantRefractionLayer( 40000.00, 1.000001));
this.refractionLayers.add(new ConstantRefractionLayer( 30000.00, 1.000004));
this.refractionLayers.add(new ConstantRefractionLayer( 23000.00, 1.000012));
this.refractionLayers.add(new ConstantRefractionLayer( 18000.00, 1.000028));
this.refractionLayers.add(new ConstantRefractionLayer( 14000.00, 1.000052));
this.refractionLayers.add(new ConstantRefractionLayer( 11000.00, 1.000083));
this.refractionLayers.add(new ConstantRefractionLayer( 9000.00, 1.000106));
this.refractionLayers.add(new ConstantRefractionLayer( 7000.00, 1.000134));
this.refractionLayers.add(new ConstantRefractionLayer( 5000.00, 1.000167));
this.refractionLayers.add(new ConstantRefractionLayer( 3000.00, 1.000206));
this.refractionLayers.add(new ConstantRefractionLayer( 1000.00, 1.000252));
this.refractionLayers.add(new ConstantRefractionLayer( 0.00, 1.000278));
this.refractionLayers.add(new ConstantRefractionLayer( -1000.00, 1.000306));
// get the lowest altitude of the atmospheric model
this.atmosphereLowestAltitude = refractionLayers.get(refractionLayers.size() - 1).getLowestAltitude();
}
/** Simple constructor.
* @param ellipsoid the ellipsoid to be used.
* @param refractionLayers the refraction layers to be used with this model (layers can be in any order).
*/
public MultiLayerModel(final ExtendedEllipsoid ellipsoid, final List<ConstantRefractionLayer> refractionLayers) {
super();
// at this stage no optimization is set up: no optimization grid is defined
this.ellipsoid = ellipsoid;
this.refractionLayers = new ArrayList<>(refractionLayers);
// sort the layers from the highest (index = 0) to the lowest (index = size - 1)
Collections.sort(this.refractionLayers,
(l1, l2) -> Double.compare(l2.getLowestAltitude(), l1.getLowestAltitude()));
// get the lowest altitude of the model
atmosphereLowestAltitude = this.refractionLayers.get(this.refractionLayers.size() - 1).getLowestAltitude();
}
/** Compute the (position, LOS) of the intersection with the lowest atmospheric layer.
* @param satPos satellite position, in body frame
* @param satLos satellite line of sight, in body frame
* @param rawIntersection intersection point without refraction correction
* @return the intersection position and LOS with the lowest atmospheric layer
* @since 2.1
*/
private IntersectionLOS computeToLowestAtmosphereLayer(
final Vector3D satPos, final Vector3D satLos,
final NormalizedGeodeticPoint rawIntersection) {
if (rawIntersection.getAltitude() < atmosphereLowestAltitude) {
throw new RuggedException(RuggedMessages.NO_LAYER_DATA, rawIntersection.getAltitude(),
atmosphereLowestAltitude);
}
Vector3D pos = satPos;
Vector3D los = satLos.normalize();
// Compute the intersection point with the lowest layer of atmosphere
double n1 = -1;
GeodeticPoint gp = ellipsoid.transform(satPos, ellipsoid.getBodyFrame(), null);
// Perform the exact computation (no optimization)
// TBN: the refractionLayers is ordered from the highest to the lowest
for (ConstantRefractionLayer refractionLayer : refractionLayers) {
if (refractionLayer.getLowestAltitude() > gp.getAltitude()) {
continue;
}
final double n2 = refractionLayer.getRefractiveIndex();
if (n1 > 0) {
// when we get here, we have already performed one iteration in the loop
// so gp is the los intersection with the layers interface (it was a
// point on ground at loop initialization, but is overridden at each iteration)
// get new los by applying Snell's law at atmosphere layers interfaces
// we avoid computing sequences of inverse-trigo/trigo/inverse-trigo functions
// we just use linear algebra and square roots, it is faster and more accurate
// at interface crossing, the interface normal is z, the local zenith direction
// the ray direction (i.e. los) is u in the upper layer and v in the lower layer
// v is in the (u, zenith) plane, so we can say
// (1) v = α u + β z
// with α>0 as u and v are roughly in the same direction as the ray is slightly bent
// let θ₁ be the los incidence angle at interface crossing
// θ₁ = π - angle(u, zenith) is between 0 and π/2 for a downwards observation
// let θ₂ be the exit angle at interface crossing
// from Snell's law, we have n₁ sin θ₁ = n₂ sin θ₂ and θ₂ is also between 0 and π/2
// we have:
// (2) u·z = -cos θ₁
// (3) v·z = -cos θ₂
// combining equations (1), (2) and (3) and remembering z·z = 1 as z is normalized , we get
// (4) β = α cos θ₁ - cos θ₂
// with all the expressions above, we can rewrite the fact v is normalized:
// 1 = v·v
// = α² u·u + 2αβ u·z + β² z·z
// = α² - 2αβ cos θ₁ + β²
// = α² - 2α² cos² θ₁ + 2 α cos θ₁ cos θ₂ + α² cos² θ₁ - 2 α cos θ₁ cos θ₂ + cos² θ₂
// = α²(1 - cos² θ₁) + cos² θ₂
// hence α² = (1 - cos² θ₂)/(1 - cos² θ₁)
// = sin² θ₂ / sin² θ₁
// as α is positive, and both θ₁ and θ₂ are between 0 and π/2, we finally get
// α = sin θ₂ / sin θ₁
// (5) α = n₁/n₂
// the α coefficient is independent from the incidence angle,
// it depends only on the ratio of refractive indices!
//
// back to equation (4) and using again the fact θ₂ is between 0 and π/2, we can now write
// β = α cos θ₁ - cos θ₂
// = n₁/n₂ cos θ₁ - cos θ₂
// = n₁/n₂ cos θ₁ - √(1 - sin² θ₂)
// = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² sin² θ₁)
// = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² (1 - cos² θ₁))
// = n₁/n₂ cos θ₁ - √(1 - (n₁/n₂)² + (n₁/n₂)² cos² θ₁)
// (6) β = -k - √(k² - ζ)
// where ζ = (n₁/n₂)² - 1 and k = n₁/n₂ u·z, which is negative, and close to -1 for
// nadir observations. As we expect atmosphere models to have small transitions between
// layers, we have to handle accurately the case where n₁/n₂ ≈ 1 so ζ ≈ 0. In this case,
// a cancellation occurs inside the square root: √(k² - ζ) ≈ √k² ≈ -k (because k is negative).
// So β ≈ -k + k ≈ 0 and another cancellation occurs, outside of the square root.
// This means that despite equation (6) is mathematically correct, it is prone to numerical
// errors when consecutive layers have close refractive indices. A better equivalent
// expression is needed. The fact β is close to 0 in this case was expected because
// equation (1) reads v = α u + β z, and α = n₁/n₂, so when n₁/n₂ ≈ 1, we have
// α ≈ 1 and β ≈ 0, so v ≈ u: when two layers have similar refractive indices, the
// propagation direction is almost unchanged.
//
// The first step for the accurate computation of β is to compute ζ = (n₁/n₂)² - 1
// accurately and avoid a cancellation just after a division (which is the least accurate
// of the four operations) and a squaring. We will simply use:
// ζ = (n₁/n₂)² - 1
// = (n₁ - n₂) (n₁ + n₂) / n₂²
// The cancellation is still there, but it occurs in the subtraction n₁ - n₂, which are
// the most accurate values we can get.
// The second step for the accurate computation of β is to rewrite equation (6)
// by both multiplying and dividing by the dual factor -k + √(k² - ζ):
// β = -k - √(k² - ζ)
// = (-k - √(k² - ζ)) * (-k + √(k² - ζ)) / (-k + √(k² - ζ))
// = (k² - (k² - ζ)) / (-k + √(k² - ζ))
// (7) β = ζ / (-k + √(k² - ζ))
// expression (7) is more stable numerically than expression (6), because when ζ ≈ 0
// its denominator is about -2k, there are no cancellation anymore after the square root.
// β is computed with the same accuracy as ζ
final double alpha = n1 / n2;
final double k = alpha * Vector3D.dotProduct(los, gp.getZenith());
final double zeta = (n1 - n2) * (n1 + n2) / (n2 * n2);
final double beta = zeta / (FastMath.sqrt(k * k - zeta) - k);
los = new Vector3D(alpha, los, beta, gp.getZenith());
}
// In case the altitude of the intersection without atmospheric refraction
// is above the lowest altitude of the atmosphere: stop the search
if (rawIntersection.getAltitude() > refractionLayer.getLowestAltitude()) {
break;
}
// Get for the intersection point: the position and the LOS
pos = ellipsoid.pointAtAltitude(pos, los, refractionLayer.getLowestAltitude());
gp = ellipsoid.transform(pos, ellipsoid.getBodyFrame(), null);
n1 = n2;
}
return new IntersectionLOS(pos, los);
}
/** {@inheritDoc} */
@Override
public NormalizedGeodeticPoint applyCorrection(final Vector3D satPos, final Vector3D satLos,
final NormalizedGeodeticPoint rawIntersection,
final IntersectionAlgorithm algorithm) {
final IntersectionLOS intersectionLOS = computeToLowestAtmosphereLayer(satPos, satLos, rawIntersection);
final Vector3D pos = intersectionLOS.getIntersectionPos();
final Vector3D los = intersectionLOS.getIntersectionLos();
// at this stage the pos belongs to the lowest atmospheric layer.
// We can compute the intersection of line of sight (los) with Digital Elevation Model
// as usual (without atmospheric refraction).
return algorithm.refineIntersection(ellipsoid, pos, los, rawIntersection);
}
} // end of class MultiLayerModel
/** Container for the (position, LOS) of the intersection with the lowest atmospheric layer.
* @author Guylaine Prat
* @since 2.1
*/
class IntersectionLOS {
/** Position of the intersection with the lowest atmospheric layer. */
private Vector3D intersectionPos;
/** LOS of the intersection with the lowest atmospheric layer. */
private Vector3D intersectionLos;
/** Default constructor.
* @param intersectionPos position of the intersection
* @param intersectionLos los of the intersection
*/
IntersectionLOS(final Vector3D intersectionPos, final Vector3D intersectionLos) {
this.intersectionPos = intersectionPos;
this.intersectionLos = intersectionLos;
}
public Vector3D getIntersectionPos() {
return intersectionPos;
}
public Vector3D getIntersectionLos() {
return intersectionLos;
}
}