TidalDisplacement.java

  1. /* Copyright 2002-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.models.earth.displacement;

  18. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  19. import org.hipparchus.util.FastMath;
  20. import org.orekit.data.BodiesElements;
  21. import org.orekit.data.PoissonSeries.CompiledSeries;
  22. import org.orekit.frames.Frame;
  23. import org.orekit.time.AbsoluteDate;
  24. import org.orekit.utils.IERSConventions;
  25. import org.orekit.utils.PVCoordinatesProvider;

  26. /**
  27.  * Modeling of displacement of reference points due to tidal effects.
  28.  * <p>
  29.  * This class implements displacement of reference point (i.e.
  30.  * {@link org.orekit.estimation.measurements.GroundStation ground stations})
  31.  * due to tidal effects, as per IERS conventions.
  32.  * </p>
  33.  * <p>
  34.  * Displacement can be computed with respect to either <em>conventional tide free</em>
  35.  * or <em>mean tide</em> coordinates. The difference between the two systems is
  36.  * about -12cm at poles and +6cm at equator. Selecting one system or the other
  37.  * depends on how the station coordinates have been computed (i.e. it depends
  38.  * whether the coordinates already include the permanent deformation or not).
  39.  * </p>
  40.  * <p>
  41.  * Instances of this class are guaranteed to be immutable
  42.  * </p>
  43.  * @see org.orekit.estimation.measurements.GroundStation
  44.  * @since 9.1
  45.  * @author Luc Maisonobe
  46.  */
  47. public class TidalDisplacement implements StationDisplacement {

  48.     /** Sun motion model. */
  49.     private final PVCoordinatesProvider sun;

  50.     /** Moon motion model. */
  51.     private final PVCoordinatesProvider moon;

  52.     /** Flag for permanent deformation. */
  53.     private final boolean removePermanentDeformation;

  54.     /** Ratio for degree 2 tide generated by Sun. */
  55.     private final double ratio2S;

  56.     /** Ratio for degree 3 tide generated by Sun. */
  57.     private final double ratio3S;

  58.     /** Ratio for degree 2 tide generated by Moon. */
  59.     private final double ratio2M;

  60.     /** Ratio for degree 3 tide generated by Moon. */
  61.     private final double ratio3M;

  62.     /** Displacement Shida number h⁽⁰⁾. */
  63.     private final double hSup0;

  64.     /** Displacement Shida number h⁽²⁾. */
  65.     private final double hSup2;

  66.     /** Displacement Shida number h₃. */
  67.     private final double h3;

  68.     /** Displacement Shida number hI diurnal. */
  69.     private final double hIDiurnal;

  70.     /** Displacement Shida number hI semi-diurnal. */
  71.     private final double hISemiDiurnal;

  72.     /** Displacement Love number l⁽⁰⁾. */
  73.     private final double lSup0;

  74.     /** Displacement Love number l⁽¹⁾ diurnal. */
  75.     private final double lSup1Diurnal;

  76.     /** Displacement Love number l⁽¹⁾ semi-diurnal. */
  77.     private final double lSup1SemiDiurnal;

  78.     /** Displacement Love number l⁽²⁾. */
  79.     private final double lSup2;

  80.     /** Displacement Love number l₃. */
  81.     private final double l3;

  82.     /** Displacement Love number lI diurnal. */
  83.     private final double lIDiurnal;

  84.     /** Displacement Love number lI semi-diurnal. */
  85.     private final double lISemiDiurnal;

  86.     /** Permanent deformation amplitude. */
  87.     private final double h0Permanent;

  88.     /** Function computing corrections in the frequency domain for diurnal tides.
  89.      * <ul>
  90.      *  <li>f[0]: radial correction, longitude cosine part</li>
  91.      *  <li>f[1]: radial correction, longitude sine part</li>
  92.      *  <li>f[2]: North correction, longitude cosine part</li>
  93.      *  <li>f[3]: North correction, longitude sine part</li>
  94.      *  <li>f[4]: East correction, longitude cosine part</li>
  95.      *  <li>f[5]: East correction, longitude sine part</li>
  96.      * </ul>
  97.      */
  98.     private final CompiledSeries frequencyCorrectionDiurnal;

  99.     /** Function computing corrections in the frequency domain for zonal tides.
  100.      * <ul>
  101.      *  <li>f[0]: radial correction</li>
  102.      *  <li>f[1]: North correction, longitude cosine part</li>
  103.      * </ul>
  104.      */
  105.     private final CompiledSeries frequencyCorrectionZonal;

  106.     /** Simple constructor.
  107.      * @param rEarth Earth equatorial radius (from gravity field model)
  108.      * @param sunEarthSystemMassRatio Sun/(Earth + Moon) mass ratio
  109.      * (typically {@link org.orekit.utils.Constants#JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO Constants.JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO})
  110.      * @param earthMoonMassRatio Earth/Moon mass ratio
  111.      * (typically {@link org.orekit.utils.Constants#JPL_SSD_EARTH_MOON_MASS_RATIO Constants.JPL_SSD_EARTH_MOON_MASS_RATIO})
  112.      * @param sun Sun model
  113.      * @param moon Moon model
  114.      * @param conventions IERS conventions to use
  115.      * @param removePermanentDeformation if true, the station coordinates are
  116.      * considered <em>mean tide</em> and already include the permanent deformation, hence
  117.      * it should be removed from the displacement to avoid considering it twice;
  118.      * if false, the station coordinates are considered <em>conventional tide free</em>
  119.      * so the permanent deformation must be included in the displacement
  120.      * @see org.orekit.frames.FramesFactory#getITRF(IERSConventions, boolean)
  121.      * @see org.orekit.frames.FramesFactory#getEOPHistory(IERSConventions, boolean)
  122.      * @see org.orekit.utils.Constants#JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO
  123.      * @see org.orekit.utils.Constants#JPL_SSD_EARTH_MOON_MASS_RATIO
  124.      */
  125.     public TidalDisplacement(final double rEarth,
  126.                              final double sunEarthSystemMassRatio,
  127.                              final double earthMoonMassRatio,
  128.                              final PVCoordinatesProvider sun, final PVCoordinatesProvider moon,
  129.                              final IERSConventions conventions,
  130.                              final boolean removePermanentDeformation) {

  131.         final double sunEarthMassRatio = sunEarthSystemMassRatio * (1 + 1 / earthMoonMassRatio);
  132.         final double moonEarthMassRatio = 1.0 / earthMoonMassRatio;

  133.         this.sun                         = sun;
  134.         this.moon                        = moon;
  135.         this.removePermanentDeformation = removePermanentDeformation;

  136.         final double r2 = rEarth * rEarth;
  137.         final double r4 = r2 * r2;
  138.         this.ratio2S    = r4 * sunEarthMassRatio;
  139.         this.ratio3S    = ratio2S * rEarth;
  140.         this.ratio2M    = r4 * moonEarthMassRatio;
  141.         this.ratio3M    = ratio2M * rEarth;

  142.         // Get the nominal values for the Love and Shiva numbers
  143.         final double[] hl = conventions.getNominalTidalDisplacement();
  144.         hSup0            = hl[ 0];
  145.         hSup2            = hl[ 1];
  146.         h3               = hl[ 2];
  147.         hIDiurnal        = hl[ 3];
  148.         hISemiDiurnal    = hl[ 4];
  149.         lSup0            = hl[ 5];
  150.         lSup1Diurnal     = hl[ 6];
  151.         lSup1SemiDiurnal = hl[ 7];
  152.         lSup2            = hl[ 8];
  153.         l3               = hl[ 9];
  154.         lIDiurnal        = hl[10];
  155.         lISemiDiurnal    = hl[11];
  156.         h0Permanent      = hl[12];

  157.         this.frequencyCorrectionDiurnal = conventions.getTidalDisplacementFrequencyCorrectionDiurnal();
  158.         this.frequencyCorrectionZonal   = conventions.getTidalDisplacementFrequencyCorrectionZonal();

  159.     }

  160.     /** {@inheritDoc} */
  161.     @Override
  162.     public Vector3D displacement(final BodiesElements elements, final Frame earthFrame, final Vector3D referencePoint) {

  163.         final AbsoluteDate date = elements.getDate();

  164.         // preliminary computation (we hold everything in local variables so method is thread-safe)
  165.         final PointData      pointData    = new PointData(referencePoint);
  166.         final Vector3D       sunPosition  = sun.getPVCoordinates(date, earthFrame).getPosition();
  167.         final BodyData       sunData      = new BodyData(sunPosition, ratio2S, ratio3S, pointData);
  168.         final Vector3D       moonPosition = moon.getPVCoordinates(date, earthFrame).getPosition();
  169.         final BodyData       moonData     = new BodyData(moonPosition, ratio2M, ratio3M, pointData);

  170.         // step 1 in IERS procedure: corrections in the time domain
  171.         Vector3D displacement = timeDomainCorrection(pointData, sunData, moonData);

  172.         // step 2 in IERS procedure: corrections in the frequency domain
  173.         displacement = displacement.add(frequencyDomainCorrection(elements, pointData));

  174.         if (removePermanentDeformation) {
  175.             // the station coordinates already include permanent deformation,
  176.             // so it should not be included in the displacement that will be
  177.             // added to these coordinates to avoid considering this deformation twice
  178.             // as step 1 did include permanent deformation, we remove it here
  179.             displacement = displacement.subtract(permanentDeformation(pointData));
  180.         }

  181.         return displacement;

  182.     }

  183.     /** Compute the corrections in the time domain (step 1 in IERS procedure).
  184.      * @param pointData reference point data
  185.      * @param sunData Sun data
  186.      * @param moonData Moon data
  187.      * @return displacement of the reference point
  188.      */
  189.     private Vector3D timeDomainCorrection(final PointData pointData,
  190.                                           final BodyData sunData, final BodyData moonData) {

  191.         final double h2  = hSup0 + hSup2 * pointData.f;
  192.         final double l2  = lSup0 + lSup2 * pointData.f;

  193.         // in-phase, degree 2 (equation 7.5 in IERS conventions 2010)
  194.         final double s2R = sunData.factor2  * 3.0 * l2 * sunData.dot;
  195.         final double s2r = sunData.factor2  * 0.5 * h2 * (3 * sunData.dot2 - 1) - s2R * sunData.dot;
  196.         final double m2R = moonData.factor2 * 3.0 * l2 * moonData.dot;
  197.         final double m2r = moonData.factor2 * 0.5 * h2 * (3 * moonData.dot2 - 1) - m2R * moonData.dot;

  198.         // in-phase, degree 3 (equation 7.6 in IERS conventions 2010)
  199.         final double s3R = sunData.factor3  * l3 * (7.5 * sunData.dot2 - 1.5);
  200.         final double s3r = sunData.factor3  * h3 * sunData.dot * (2.5 * sunData.dot2 - 1.5) - s3R * sunData.dot;
  201.         final double m3R = moonData.factor3 * l3 * (7.5 * moonData.dot2 - 1.5);
  202.         final double m3r = moonData.factor3 * h3 * moonData.dot * (2.5 * moonData.dot2 - 1.5) - m3R * moonData.dot;

  203.         // combine contributions along radial, Sun and Moon directions
  204.         final Vector3D inPhaseDisplacement = new Vector3D(s2r + m2r + s3r + m3r,    pointData.radial,
  205.                                                           (s2R + s3R) / sunData.r,  sunData.position,
  206.                                                           (m2R + m3R) / moonData.r, moonData.position);

  207.         // out-of-phase, degree 2, diurnal tides (equations 7.10a and 7.10b in IERS conventions 2010)
  208.         final double drOd = -0.75 * hIDiurnal * pointData.sin2Phi *
  209.                             (sunData.factor2  * sunData.sin2Phi  * sunData.sinDeltaLambda +
  210.                              moonData.factor2 * moonData.sin2Phi * moonData.sinDeltaLambda);
  211.         final double dnOd = -1.5 * lIDiurnal * pointData.cos2Phi *
  212.                             (sunData.factor2  * sunData.sin2Phi  * sunData.sinDeltaLambda +
  213.                              moonData.factor2 * moonData.sin2Phi * moonData.sinDeltaLambda);
  214.         final double deOd = -1.5 * lIDiurnal * pointData.sinPhi *
  215.                             (sunData.factor2  * sunData.sin2Phi  * sunData.cosDeltaLambda +
  216.                              moonData.factor2 * moonData.sin2Phi * moonData.cosDeltaLambda);

  217.         // out-of-phase, degree 2, semi-diurnal tides (equation 7.11 in IERS conventions 2010)
  218.         final double drOsd = -0.75 * hISemiDiurnal * pointData.cosPhi2 *
  219.                              (sunData.factor2  * sunData.cosPhi2  * sunData.sin2DeltaLambda +
  220.                               moonData.factor2 * moonData.cosPhi2 * moonData.sin2DeltaLambda);
  221.         final double dnOsd = 0.75 * lISemiDiurnal * pointData.sin2Phi *
  222.                              (sunData.factor2  * sunData.cosPhi2  * sunData.sin2DeltaLambda +
  223.                               moonData.factor2 * moonData.cosPhi2 * moonData.sin2DeltaLambda);
  224.         final double deOsd = -1.5 * lISemiDiurnal * pointData.cosPhi *
  225.                              (sunData.factor2  * sunData.cosPhi2  * sunData.cos2DeltaLambda +
  226.                               moonData.factor2 * moonData.cosPhi2 * moonData.cos2DeltaLambda);

  227.         // latitude dependency, diurnal tides (equation 7.8 in IERS conventions 2010)
  228.         final double dnLd = -lSup1Diurnal * pointData.sinPhi2 *
  229.                             (sunData.factor2  * sunData.p21  * sunData.cosDeltaLambda +
  230.                              moonData.factor2 * moonData.p21 * moonData.cosDeltaLambda);
  231.         final double deLd =  lSup1Diurnal * pointData.sinPhi * pointData.cos2Phi *
  232.                             (sunData.factor2  * sunData.p21  * sunData.sinDeltaLambda +
  233.                              moonData.factor2 * moonData.p21 * moonData.sinDeltaLambda);

  234.         // latitude dependency, semi-diurnal tides (equation 7.9 in IERS conventions 2010)
  235.         final double dnLsd = -0.25 * lSup1SemiDiurnal * pointData.sin2Phi *
  236.                              (sunData.factor2  * sunData.p22  * sunData.cos2DeltaLambda +
  237.                               moonData.factor2 * moonData.p22 * moonData.cos2DeltaLambda);
  238.         final double deLsd = -0.25 * lSup1SemiDiurnal * pointData.sin2Phi * pointData.sinPhi *
  239.                              (sunData.factor2  * sunData.p22  * sunData.sin2DeltaLambda +
  240.                               moonData.factor2 * moonData.p22 * moonData.sin2DeltaLambda);

  241.         // combine diurnal and semi-diurnal tides
  242.         final Vector3D outOfPhaseDisplacement = new Vector3D(drOd + drOsd,                pointData.radial,
  243.                                                              dnOd + dnOsd + dnLd + dnLsd, pointData.north,
  244.                                                              deOd + deOsd + deLd + deLsd, pointData.east);

  245.         return inPhaseDisplacement.add(outOfPhaseDisplacement);

  246.     }

  247.     /** Compute the corrections in the frequency domain (step 2 in IERS procedure).
  248.      * @param elements elements affecting Earth orientation
  249.      * @param pointData reference point data
  250.      * @return displacement of the reference point
  251.      */
  252.     private Vector3D frequencyDomainCorrection(final BodiesElements elements, final PointData pointData) {

  253.         // corrections due to diurnal tides
  254.         final double[] cD  = frequencyCorrectionDiurnal.value(elements);
  255.         final double   drD = pointData.sin2Phi * (cD[0] * pointData.cosLambda + cD[1] * pointData.sinLambda);
  256.         final double   dnD = pointData.cos2Phi * (cD[2] * pointData.cosLambda + cD[3] * pointData.sinLambda);
  257.         final double   deD = pointData.sinPhi  * (cD[4] * pointData.cosLambda + cD[5] * pointData.sinLambda);

  258.         // corrections due to zonal long period tides
  259.         final double[] cZ  = frequencyCorrectionZonal.value(elements);
  260.         final double   drZ = (1.5 * pointData.sinPhi2 - 0.5) * cZ[0];
  261.         final double   dnZ = pointData.sin2Phi               * cZ[1];

  262.         return new Vector3D(drD + drZ, pointData.radial,
  263.                             dnD + dnZ, pointData.north,
  264.                             deD,       pointData.east);

  265.     }

  266.     /** Compute the permanent part of the deformation.
  267.      * @param pointData reference point data
  268.      * @return displacement of the reference point
  269.      */
  270.     private Vector3D permanentDeformation(final PointData pointData) {

  271.         final double h2  = hSup0 + hSup2 * pointData.f;
  272.         final double l2  = lSup0 + lSup2 * pointData.f;

  273.         // permanent deformation, which depend only on latitude
  274.         final double factor = FastMath.sqrt(1.25 / FastMath.PI);
  275.         final double dr = factor *       h2 * h0Permanent * pointData.f;
  276.         final double dn = factor * 1.5 * l2 * h0Permanent * pointData.sin2Phi;

  277.         return new Vector3D(dr, pointData.radial,
  278.                             dn, pointData.north);

  279.     }

  280.     /** Holder for various intermediate data related to reference point. */
  281.     private static class PointData {

  282.         /** Reference point position in {@link #getEarthFrame() Earth frame}. */
  283.         private final Vector3D position;

  284.         /** Distance to geocenter. */
  285.         private final double   r;

  286.         /** Sine of geocentric latitude (NOT geodetic latitude). */
  287.         private final double   sinPhi;

  288.         /** Cosine of geocentric latitude (NOT geodetic latitude). */
  289.         private final double   cosPhi;

  290.         /** Square of the sine of the geocentric latitude (NOT geodetic latitude). */
  291.         private final double   sinPhi2;

  292.         /** Square of the cosine of the geocentric latitude (NOT geodetic latitude). */
  293.         private final double   cosPhi2;

  294.         /** Sine of twice the geocentric latitude (NOT geodetic latitude). */
  295.         private final double   sin2Phi;

  296.         /** Cosine of twice the geocentric latitude (NOT geodetic latitude). */
  297.         private final double   cos2Phi;

  298.         /** Sine of longitude. */
  299.         private final double   sinLambda;

  300.         /** Cosine of longitude. */
  301.         private final double   cosLambda;

  302.         /** Unit radial vector (NOT zenith as it starts from geocenter). */
  303.         private final Vector3D radial;

  304.         /** Unit vector in North direction. */
  305.         private final Vector3D north;

  306.         /** Unit vector in East direction. */
  307.         private final Vector3D east;

  308.         /** (3 sin²φ - 1) / 2 where φ is geoCENTRIC latitude of the reference point (NOT geodetic latitude). */
  309.         private final double f;

  310.         /** Simple constructor.
  311.          * @param position reference point position in {@link #getEarthFrame() Earth frame}
  312.          */
  313.         PointData(final Vector3D position) {

  314.             this.position   = position;
  315.             final double x  = position.getX();
  316.             final double y  = position.getY();
  317.             final double z  = position.getZ();
  318.             final double x2 = x * x;
  319.             final double y2 = y * y;
  320.             final double z2 = z * z;

  321.             // preliminary computation related to station position
  322.             final double rho2 = x2 + y2;
  323.             final double rho  = FastMath.sqrt(rho2);
  324.             final double r2   = rho2 + z2;
  325.             r                 = FastMath.sqrt(r2);
  326.             sinPhi            = z / r;
  327.             cosPhi            = rho / r;
  328.             sinPhi2           = sinPhi * sinPhi;
  329.             cosPhi2           = cosPhi * cosPhi;
  330.             sin2Phi           = 2 * sinPhi * cosPhi;
  331.             cos2Phi           = cosPhi2 - sinPhi2;
  332.             if (rho == 0.0) {
  333.                 // at pole
  334.                 sinLambda = 0.0;
  335.                 cosLambda = 1.0;
  336.             } else {
  337.                 // regular point
  338.                 sinLambda = y / rho;
  339.                 cosLambda = x / rho;
  340.             }
  341.             radial = new Vector3D(x / r, y / r, sinPhi);
  342.             north  = new Vector3D(-cosLambda * sinPhi, -sinLambda * sinPhi, cosPhi);
  343.             east   = new Vector3D(-sinLambda, cosLambda, 0);

  344.             // (3 sin²φ - 1) / 2 where φ is geoCENTRIC latitude of the reference point (NOT geodetic latitude)
  345.             f = (z2 - 0.5 * rho2) / r2;

  346.         }

  347.     }

  348.     /** Holder for various intermediate data related to tide generating body. */
  349.     private static class BodyData {

  350.         /** Body position in Earth frame. */
  351.         private final Vector3D position;

  352.         /** Distance to geocenter. */
  353.         private final double r;

  354.         /** Dot product with reference point direction. */
  355.         private final double dot;

  356.         /** Squared dot product with reference point direction. */
  357.         private final double dot2;

  358.         /** Factor for degree 2 tide. */
  359.         private final double factor2;

  360.         /** Factor for degree 3 tide. */
  361.         private final double factor3;

  362.         /** Square of the cosine of the geocentric latitude (NOT geodetic latitude). */
  363.         private final double   cosPhi2;

  364.         /** Sine of twice the geocentric latitude (NOT geodetic latitude). */
  365.         private final double   sin2Phi;

  366.         /** Legendre function P₂¹. */
  367.         private final double p21;

  368.         /** Legendre function P₂². */
  369.         private final double p22;

  370.         /** Sine of the longitude difference with reference point. */
  371.         private final double sinDeltaLambda;

  372.         /** Cosine of the longitude difference with reference point. */
  373.         private final double cosDeltaLambda;

  374.         /** Sine of twice the longitude difference with reference point. */
  375.         private final double sin2DeltaLambda;

  376.         /** Cosine of twice the longitude difference with reference point. */
  377.         private final double cos2DeltaLambda;

  378.         /** Simple constructor.
  379.          * @param position body position in Earth frame
  380.          * @param ratio2 ratio for the degree 2 tide generated by this body
  381.          * @param ratio3 ratio for the degree 3 tide generated by this body
  382.          * @param pointData reference point data
  383.          */
  384.         BodyData(final Vector3D position, final double ratio2, final double ratio3,
  385.                  final PointData pointData) {

  386.             final double x  = position.getX();
  387.             final double y  = position.getY();
  388.             final double z  = position.getZ();
  389.             final double x2 = x * x;
  390.             final double y2 = y * y;
  391.             final double z2 = z * z;

  392.             this.position    = position;
  393.             final double r2  = x2 + y2 + z2;
  394.             r                = FastMath.sqrt(r2);
  395.             dot              = Vector3D.dotProduct(position, pointData.position) / (r * pointData.r);
  396.             dot2             = dot * dot;

  397.             factor2          = ratio2 / (r2 * r);
  398.             factor3          = ratio3 / (r2 * r2);

  399.             final double rho       = FastMath.sqrt(x2 + y2);
  400.             final double sinPhi    = z / r;
  401.             final double cosPhi    = rho / r;
  402.             final double sinCos    = sinPhi * cosPhi;
  403.             cosPhi2                = cosPhi * cosPhi;
  404.             sin2Phi                = 2 * sinCos;
  405.             p21                    = 3 * sinCos;
  406.             p22                    = 3 * cosPhi2;

  407.             final double sinLambda = y / rho;
  408.             final double cosLambda = x / rho;
  409.             sinDeltaLambda  = pointData.sinLambda * cosLambda  - pointData.cosLambda * sinLambda;
  410.             cosDeltaLambda  = pointData.cosLambda * cosLambda  + pointData.sinLambda * sinLambda;
  411.             sin2DeltaLambda = 2 * sinDeltaLambda * cosDeltaLambda;
  412.             cos2DeltaLambda = cosDeltaLambda * cosDeltaLambda - sinDeltaLambda * sinDeltaLambda;

  413.         }

  414.     }

  415. }