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 (Eq. 153 to 155)
  93.         final T cosD = scLatRec.sin().multiply(scLatSat.sin()).
  94.                        add(scLatRec.cos().multiply(scLatSat.cos()).multiply(scLon21.cos()));
  95.         final T sinD = FastMath.sqrt(cosD.multiply(cosD).negate().add(1.0));
  96.         final T z    = FastMath.atan2(sinD, cosD.subtract(r1.divide(r2)));
  97.         final FieldSinCos<T> scZ = FastMath.sinCos(z);

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

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

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

  108.             // Ray-perigee longitude (Eq. 164)
  109.             if (z.getReal() < 0) {
  110.                 this.lonP = lon2;
  111.             } else {
  112.                 this.lonP = lon2.add(pi);
  113.             }

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

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

  121.         } else {

  123.             // Ray-perigee latitude (Eq. 158 to 163)
  124.             final T sinAz   = FastMath.sin(lon2.subtract(lon1)).multiply(scLatSat.cos()).divide(sinD);
  125.             final T cosAz   = scLatSat.sin().subtract(cosD.multiply(scLatRec.sin())).
  126.                               divide(sinD.multiply(scLatRec.cos()));
  127.             final T sinLatP = scLatRec.sin().multiply(scZ.sin()).
  128.                               subtract(scLatRec.cos().multiply(scZ.cos()).multiply(cosAz));
  129.             final T cosLatP = FastMath.sqrt(sinLatP.multiply(sinLatP).negate().add(1.0));
  130.             this.latP       = FastMath.atan2(sinLatP, cosLatP);

  132.             // Ray-perigee longitude (Eq. 165 to 167)
  133.             final T sinLonP = sinAz.negate().multiply(scZ.cos()).divide(cosLatP);
  134.             final T cosLonP = scZ.sin().subtract(scLatRec.sin().multiply(sinLatP)).
  135.                               divide(scLatRec.cos().multiply(cosLatP));
  136.             this.lonP       = FastMath.atan2(sinLonP, cosLonP).add(lon1);

  138.         }

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

  143.         if (FastMath.abs(FastMath.abs(latP).subtract(pi.multiply(0.5)).getReal()) < THRESHOLD || FastMath.abs(
  144.             scZ.sin().getReal()) < THRESHOLD) {
  145.             // Eq. 172 and 173
  146.             this.sinAzP = pi.getField().getZero();
  147.             this.cosAzP = FastMath.copySign(pi.getField().getOne(), latP).negate();
  148.         } else {
  149.             final FieldSinCos<T> scLon = FastMath.sinCos(lon2.subtract(lonP));
  150.             // Sine and cosine of azimuth of satellite as seen from ray-perigee
  151.             final FieldSinCos<T> scPsi = FastMath.sinCos(greatCircleAngle(scLatSat, scLon));
  152.             // Eq. 174 and 175
  153.             this.sinAzP = scLatSat.cos().multiply(scLon.sin()).divide(scPsi.sin());
  154.             this.cosAzP = scLatSat.sin().subtract(scLatP.sin().multiply(scPsi.cos())).
  155.                           divide(scLatP.cos().multiply(scPsi.sin()));
  156.         }

  158.         // Integration end points s1 and s2 in meters (Eq. 176 and 177)
  159.         this.s1 = FastMath.sqrt(r1.multiply(r1).subtract(rp.multiply(rp)));
  160.         this.s2 = FastMath.sqrt(r2.multiply(r2).subtract(rp.multiply(rp)));
  161.     }

  163.     /**
  164.      * Get receiver altitude.
  165.      * @return receiver altitude
  166.      * @since 13.0
  167.      */
  168.     public T getRecH() {
  169.         return recH;
  170.     }

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

  181.     /**
  182.      * Get the distance of the first point from the ray perigee.
  183.      *
  184.      * @return s1 in meters
  185.      */
  186.     public T getS1() {
  187.         return s1;
  188.     }

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

  199.     /**
  200.      * Get the ray-perigee radius.
  201.      *
  202.      * @return the ray-perigee radius in meters
  203.      */
  204.     public T getRadius() {
  205.         return rp;
  206.     }

  208.     /**
  209.      * Get the ray-perigee latitude.
  210.      *
  211.      * @return the ray-perigee latitude in radians
  212.      */
  213.     public T getLatitude() {
  214.         return latP;
  215.     }

  217.     /**
  218.      * Get the ray-perigee latitude sin/cos.
  219.      *
  220.      * @return the ray-perigee latitude sin/cos
  221.      * @since 13.0
  222.      */
  223.     public FieldSinCos<T> getScLat() {
  224.         return scLatP;
  225.     }

  227.     /**
  228.      * Get the ray-perigee longitude.
  229.      *
  230.      * @return the ray-perigee longitude in radians
  231.      */
  232.     public T getLongitude() {
  233.         return lonP;
  234.     }

  236.     /**
  237.      * Get the sine of azimuth of satellite as seen from ray-perigee.
  238.      *
  239.      * @return the sine of azimuth
  240.      */
  241.     public T getSineAz() {
  242.         return sinAzP;
  243.     }

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

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

  275. }
Created with JaCoCo 0.8.12.202403310830