FieldRay.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.CalculusFieldElement;
  20. import org.hipparchus.util.FastMath;
  21. import org.hipparchus.util.FieldSinCos;
  22. import org.orekit.bodies.FieldGeodeticPoint;

  24. /** Container for ray-perigee parameters.
  25.  * <p>By convention, point 1 is at lower height.</p>
  26.  * @author Bryan Cazabonne
  27.  * @since 13.0
  28.  */
  29. class FieldRay<T extends CalculusFieldElement<T>> {

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

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

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

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

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

  50.     /** Ray-perigee radius [m]. */
  51.     private final T rp;

  53.     /** Ray-perigee latitude [rad]. */
  54.     private final T latP;

  56.     /** Ray-perigee longitude [rad]. */
  57.     private final T lonP;

  59.     /** Sine and cosine of ray-perigee latitude. */
  60.     private final FieldSinCos<T> scLatP;

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

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

  68.     /**
  69.      * Constructor.
  70.      *
  71.      * @param recP  receiver position
  72.      * @param satP  satellite position
  73.      */
  74.     FieldRay(final FieldGeodeticPoint<T> recP, final FieldGeodeticPoint<T> satP) {

  76.         // Integration limits in meters (Eq. 140 and 141)
  77.         this.recH  = recP.getAltitude();
  78.         this.satH  = satP.getAltitude();
  79.         final T r1 = recH.add(NeQuickModel.RE);
  80.         final T r2 = satH.add(NeQuickModel.RE);

  82.         // Useful parameters
  83.         final T pi = r1.getPi();
  84.         final T lat1 = recP.getLatitude();
  85.         final T lat2 = satP.getLatitude();
  86.         final T lon1 = recP.getLongitude();
  87.         final T lon2 = satP.getLongitude();
  88.         final FieldSinCos<T> scLatSat = FastMath.sinCos(lat2);
  89.         final FieldSinCos<T> scLatRec = FastMath.sinCos(lat1);
  90.         final FieldSinCos<T> scLon21 = FastMath.sinCos(lon2.subtract(lon1));

  92.         // Zenith angle computation, using:
  93.         //   - Eq. 153 with added protection against numerical noise near zenith observation
  94.         //   - replacing Eq. 154 by a different one more stable near zenith
  95.         //   - replacing Eq. 155 by different ones avoiding trigonometric functions
  96.         final T cosD     = FastMath.min(r1.getField().getOne(),
  97.                                         scLatRec.sin().multiply(scLatSat.sin()).
  98.                                         add(scLatRec.cos().multiply(scLatSat.cos()).multiply(scLon21.cos())));
  99.         final T sinDSinM = scLatSat.cos().multiply(scLon21.sin());
  100.         final T sinDCosM = scLatSat.sin().multiply(scLatRec.cos()).subtract(scLatSat.cos().multiply(scLatRec.sin()).multiply(scLon21.cos()));
  101.         final T sinD     = FastMath.sqrt(sinDSinM.multiply(sinDSinM).add(sinDCosM.multiply(sinDCosM)));
  102.         final T u        = r2.multiply(sinD);
  103.         final T v        = r2.multiply(cosD).subtract(r1);
  104.         final T inv      = FastMath.sqrt(u.multiply(u).add(v.multiply(v))).reciprocal();
  105.         final T sinZ     = u.multiply(inv);
  106.         final T cosZ     = v.multiply(inv);

  108.         // Ray-perigee computation in meters (Eq. 156)
  109.         this.rp = r1.multiply(sinZ);

  111.         // Ray-perigee latitude and longitude
  112.         if (FastMath.abs(FastMath.abs(lat1).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD) {

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

  117.             // Ray-perigee longitude (Eq. 164)
  118.             this.lonP = lon2.add(pi);

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

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

  126.         } else {

  128.             // Ray-perigee latitude (Eq. 158 to 163)
  129.             final T sinAz   = FastMath.sin(lon2.subtract(lon1)).multiply(scLatSat.cos()).divide(sinD);
  130.             final T cosAz   = scLatSat.sin().subtract(cosD.multiply(scLatRec.sin())).
  131.                               divide(sinD.multiply(scLatRec.cos()));
  132.             final T sinLatP = scLatRec.sin().multiply(sinZ).
  133.                               subtract(scLatRec.cos().multiply(cosZ).multiply(cosAz));
  134.             final T cosLatP = FastMath.sqrt(sinLatP.multiply(sinLatP).negate().add(1.0));
  135.             this.latP       = FastMath.atan2(sinLatP, cosLatP);

  137.             // Ray-perigee longitude (Eq. 165 to 167, plus protection against ray-perigee along polar axis)
  138.             if (cosLatP.getReal() < THRESHOLD) {
  139.                 this.lonP = cosLatP.getField().getZero();
  140.             } else {
  141.                 final T sinLonP = sinAz.negate().multiply(cosZ).divide(cosLatP);
  142.                 final T cosLonP = sinZ.subtract(scLatRec.sin().multiply(sinLatP)).
  143.                         divide(scLatRec.cos().multiply(cosLatP));
  144.                 this.lonP = FastMath.atan2(sinLonP, cosLonP).add(lon1);
  145.             }

  147.         }

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

  152.         if (FastMath.abs(FastMath.abs(latP).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD || FastMath.abs(
  153.             sinZ.getReal()) < THRESHOLD) {
  154.             // Eq. 172 and 173
  155.             this.sinAzP = pi.getField().getZero();
  156.             this.cosAzP = FastMath.copySign(pi.getField().getOne(), latP).negate();
  157.         } else {
  158.             final FieldSinCos<T> scLon = FastMath.sinCos(lon2.subtract(lonP));
  159.             // Sine and cosine of azimuth of satellite as seen from ray-perigee
  160.             final FieldSinCos<T> scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon));
  161.             // Eq. 174 and 175
  162.             this.sinAzP = scLatSat.cos().multiply(scLon.sin()).divide(scPsi.sin());
  163.             this.cosAzP = scLatSat.sin().subtract(scLatP.sin().multiply(scPsi.cos())).
  164.                           divide(scLatP.cos().multiply(scPsi.sin()));
  165.         }

  167.         // Integration end points s1 and s2 in meters (Eq. 176 and 177)
  168.         this.s1 = FastMath.sqrt(r1.multiply(r1).subtract(rp.multiply(rp)));
  169.         this.s2 = FastMath.sqrt(r2.multiply(r2).subtract(rp.multiply(rp)));
  170.     }

  172.     /**
  173.      * Get receiver altitude.
  174.      * @return receiver altitude
  175.      * @since 13.0
  176.      */
  177.     public T getRecH() {
  178.         return recH;
  179.     }

  181.     /**
  182.      * Get satellite altitude.
  183.      * @return satellite altitude
  184.      * @since 13.0
  185.      */
  186.     public T getSatH() {
  187.         return satH;
  188.     }

  190.     /**
  191.      * Get the distance of the first point from the ray perigee.
  192.      *
  193.      * @return s1 in meters
  194.      */
  195.     public T getS1() {
  196.         return s1;
  197.     }

  199.     /**
  200.      * Get the distance of the second point from the ray perigee.
  201.      *
  202.      * @return s2 in meters
  203.      */
  204.     public T getS2() {
  205.         return s2;
  206.     }

  208.     /**
  209.      * Get the ray-perigee radius.
  210.      *
  211.      * @return the ray-perigee radius in meters
  212.      */
  213.     public T getRadius() {
  214.         return rp;
  215.     }

  217.     /**
  218.      * Get the ray-perigee latitude.
  219.      *
  220.      * @return the ray-perigee latitude in radians
  221.      */
  222.     public T getLatitude() {
  223.         return latP;
  224.     }

  226.     /**
  227.      * Get the ray-perigee latitude sin/cos.
  228.      *
  229.      * @return the ray-perigee latitude sin/cos
  230.      * @since 13.0
  231.      */
  232.     public FieldSinCos<T> getScLat() {
  233.         return scLatP;
  234.     }

  236.     /**
  237.      * Get the ray-perigee longitude.
  238.      *
  239.      * @return the ray-perigee longitude in radians
  240.      */
  241.     public T getLongitude() {
  242.         return lonP;
  243.     }

  245.     /**
  246.      * Get the sine of azimuth of satellite as seen from ray-perigee.
  247.      *
  248.      * @return the sine of azimuth
  249.      */
  250.     public T getSineAz() {
  251.         return sinAzP;
  252.     }

  254.     /**
  255.      * Get the cosine of azimuth of satellite as seen from ray-perigee.
  256.      *
  257.      * @return the cosine of azimuth
  258.      */
  259.     public T getCosineAz() {
  260.         return cosAzP;
  261.     }

  263.     /**
  264.      * Compute the great circle angle from ray-perigee to satellite.
  265.      * <p>
  266.      * This method used the equations 168 to 171 of the reference document.
  267.      * </p>
  268.      *
  269.      * @param scLat sine and cosine of satellite latitude
  270.      * @param scLon sine and cosine of satellite longitude minus receiver longitude
  271.      * @return the great circle angle in radians
  272.      */
  273.     private T greatCircleAngle(final FieldSinCos<T> scLat, final FieldSinCos<T> scLon) {
  274.         if (FastMath.abs(FastMath.abs(latP).getReal() - 0.5 * FastMath.PI) < THRESHOLD) {
  275.             return FastMath.abs(FastMath.asin(scLat.sin()).subtract(latP));
  276.         } else {
  277.             final T cosPhi = scLatP.sin().multiply(scLat.sin()).
  278.                              add(scLatP.cos().multiply(scLat.cos()).multiply(scLon.cos()));
  279.             final T sinPhi = FastMath.sqrt(cosPhi.multiply(cosPhi).negate().add(1.0));
  280.             return FastMath.atan2(sinPhi, cosPhi);
  281.         }
  282.     }

  284. }
Created with JaCoCo 0.8.12.202403310830