MultiLayerModel.java

  1. /* Copyright 2013-2022 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;

  19. import java.util.ArrayList;
  20. import java.util.Collections;
  21. import java.util.List;

  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;

  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 {

  41.     /** Observed body ellipsoid. */
  42.     private final ExtendedEllipsoid ellipsoid;

  44.     /** Constant refraction layers. */
  45.     private final List<ConstantRefractionLayer> refractionLayers;

  47.     /** Atmosphere lowest altitude (m). */
  48.     private final double atmosphereLowestAltitude;

  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) {

  58.         super();

  60.         this.ellipsoid = ellipsoid;

  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));

  79.         // get the lowest altitude of the atmospheric model
  80.         this.atmosphereLowestAltitude = refractionLayers.get(refractionLayers.size() - 1).getLowestAltitude();
  81.     }

  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) {

  89.         super();
  90.         // at this stage no optimization is set up: no optimization grid is defined

  92.         this.ellipsoid = ellipsoid;
  93.         this.refractionLayers = new ArrayList<>(refractionLayers);

  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()));

  99.         // get the lowest altitude of the model
  100.         atmosphereLowestAltitude = this.refractionLayers.get(this.refractionLayers.size() - 1).getLowestAltitude();
  101.     }

  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) {

  114.         if (rawIntersection.getAltitude() < atmosphereLowestAltitude) {
  115.             throw new RuggedException(RuggedMessages.NO_LAYER_DATA, rawIntersection.getAltitude(),
  116.                                       atmosphereLowestAltitude);
  117.         }

  119.         Vector3D pos = satPos;
  120.         Vector3D los = satLos.normalize();

  122.         // Compute the intersection point with the lowest layer of atmosphere
  123.         double n1 = -1;
  124.         GeodeticPoint gp = ellipsoid.transform(satPos, ellipsoid.getBodyFrame(), null);

  126.         // Perform the exact computation (no optimization)
  127.         // TBN: the refractionLayers is ordered from the highest to the lowest
  128.         for (ConstantRefractionLayer refractionLayer : refractionLayers) {

  130.             if (refractionLayer.getLowestAltitude() > gp.getAltitude()) {
  131.                 continue;
  132.             }

  134.             final double n2 = refractionLayer.getRefractiveIndex();

  136.             if (n1 > 0) {

  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)

  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

  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

  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.             }

  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.             }

  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);

  228.             n1 = n2;
  229.         }
  230.         return new IntersectionLOS(pos, los);
  231.     }

  233.     /** {@inheritDoc} */
  234.     @Override
  235.     public NormalizedGeodeticPoint applyCorrection(final Vector3D satPos, final Vector3D satLos,
  236.                                                    final NormalizedGeodeticPoint rawIntersection,
  237.                                                    final IntersectionAlgorithm algorithm) {

  239.         final IntersectionLOS intersectionLOS = computeToLowestAtmosphereLayer(satPos, satLos, rawIntersection);
  240.         final Vector3D pos = intersectionLOS.getIntersectionPos();
  241.         final Vector3D los = intersectionLOS.getIntersectionLos();

  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.     }

  249. } // end of class MultiLayerModel

  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 {

  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;

  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.     }

  274.     public Vector3D getIntersectionLos() {
  275.         return intersectionLos;
  276.     }
  277. }
Created with JaCoCo 0.8.7.202105040129