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

  99.         // Ray-perigee computation in meters (Eq. 156)
  100.         this.rp = r1 * scZ.sin();

  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

  106.             // Ray-perigee latitude (Eq. 157)
  107.             this.latP = FastMath.copySign(z, lat1);

  109.             // Ray-perigee longitude (Eq. 164)
  110.             if (z < 0) {
  111.                 this.lonP = lon2;
  112.             } else {
  113.                 this.lonP = lon2 + FastMath.PI;
  114.             }

  116.         } else if (FastMath.abs(scZ.sin()) < THRESHOLD) {
  117.             // satellite is almost on receiver zenith

  119.             this.latP = recP.getLatitude();
  120.             this.lonP = recP.getLongitude();

  122.         } else {

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

  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;

  136.         }

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

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

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

  159.     /**
  160.      * Get receiver altitude.
  161.      * @return receiver altitude
  162.      * @since 13.0
  163.      */
  164.     public double getRecH() {
  165.         return recH;
  166.     }

  168.     /**
  169.      * Get satellite altitude.
  170.      * @return satellite altitude
  171.      * @since 13.0
  172.      */
  173.     public double getSatH() {
  174.         return satH;
  175.     }

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

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

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

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

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

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

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

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

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

  270. }
Created with JaCoCo 0.8.12.202403310830