Ray.java

  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;

  19. import org.hipparchus.util.FastMath;
  20. import org.hipparchus.util.SinCos;
  21. import org.orekit.bodies.GeodeticPoint;

  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 {

  30.     /** Threshold for ray-perigee parameters computation. */
  31.     private static final double THRESHOLD = 1.0e-10;

  33.     /** Receiver altitude [m].
  34.      * @since 13.0
  35.      */
  36.     private final double recH;

  38.     /** Satellite altitude [m].
  39.      * @since 13.0
  40.      */
  41.     private final double satH;

  43.     /** Distance of the first point from the ray perigee [m]. */
  44.     private final double s1;

  46.     /** Distance of the second point from the ray perigee [m]. */
  47.     private final double s2;

  49.     /** Ray-perigee radius [m]. */
  50.     private final double rp;

  52.     /** Ray-perigee latitude [rad]. */
  53.     private final double latP;

  55.     /** Ray-perigee longitude [rad]. */
  56.     private final double lonP;

  58.     /** Sine and cosine of ray-perigee latitude. */
  59.     private final SinCos scLatP;

  61.     /** Sine of azimuth of satellite as seen from ray-perigee. */
  62.     private final double sinAzP;

  64.     /** Cosine of azimuth of satellite as seen from ray-perigee. */
  65.     private final double cosAzP;

  67.     /**
  68.      * Constructor.
  69.      *
  70.      * @param recP receiver position
  71.      * @param satP satellite position
  72.      */
  73.     public Ray(final GeodeticPoint recP, final GeodeticPoint satP) {

  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;

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

  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;

  106.         // Ray-perigee computation in meters (Eq. 156)
  107.         this.rp = r1 * sinZ;

  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

  113.             // Ray-perigee latitude (Eq. 157)
  114.             this.latP = FastMath.copySign(FastMath.atan2(u, v), lat1);

  116.             // Ray-perigee longitude (Eq. 164)
  117.             this.lonP = lon2 + FastMath.PI;

  119.         } else if (FastMath.abs(sinZ) < THRESHOLD) {
  120.             // satellite is almost on receiver zenith

  122.             this.latP = recP.getLatitude();
  123.             this.lonP = recP.getLongitude();

  125.         } else {

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

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

  143.         }

  145.         // Sine and cosine of ray-perigee latitude
  146.         this.scLatP = FastMath.sinCos(latP);

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

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

  166.     /**
  167.      * Get receiver altitude.
  168.      * @return receiver altitude
  169.      * @since 13.0
  170.      */
  171.     public double getRecH() {
  172.         return recH;
  173.     }

  175.     /**
  176.      * Get satellite altitude.
  177.      * @return satellite altitude
  178.      * @since 13.0
  179.      */
  180.     public double getSatH() {
  181.         return satH;
  182.     }

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

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

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

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

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

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

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

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

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

  277. }
Created with JaCoCo 0.8.12.202403310830