MariniMurray.java

  1. /* Copyright 2011-2012 Space Applications Services
  2.  * Licensed to CS Communication & Systèmes (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.troposphere;

  19. import java.util.Collections;
  20. import java.util.List;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.util.FastMath;
  24. import org.orekit.bodies.FieldGeodeticPoint;
  25. import org.orekit.bodies.GeodeticPoint;
  26. import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
  27. import org.orekit.models.earth.weather.PressureTemperatureHumidity;
  28. import org.orekit.time.AbsoluteDate;
  29. import org.orekit.time.FieldAbsoluteDate;
  30. import org.orekit.utils.FieldTrackingCoordinates;
  31. import org.orekit.utils.ParameterDriver;
  32. import org.orekit.utils.TrackingCoordinates;
  33. import org.orekit.utils.units.Unit;
  34. import org.orekit.utils.units.UnitsConverter;

  36. /** The Marini-Murray tropospheric delay model for laser ranging.
  37.  *
  38.  * @see "Marini, J.W., and C.W. Murray, correction of Laser Range Tracking Data for
  39.  *      Atmospheric Refraction at Elevations Above 10 degrees, X-591-73-351, NASA GSFC, 1973"
  40.  *
  41.  * @author Joris Olympio
  42.  * @author Luc Maisonobe
  43.  * @since 12.1
  44.  */
  45. public class MariniMurray implements TroposphericModel {

  47.     /** Laser frequency parameter. */
  48.     private final double fLambda;

  50.     /** Create a new Marini-Murray model for the troposphere.
  51.      * @param lambda laser wavelength
  52.      * @param lambdaUnits units in which {@code lambda} is given
  53.      * @see TroposphericModelUtils#MICRO_M
  54.      * @see TroposphericModelUtils#NANO_M
  55.      * @since 12.1
  56.      * */
  57.     public MariniMurray(final double lambda, final Unit lambdaUnits) {

  59.         // compute laser frequency parameter
  60.         final double lambdaMicrometer = new UnitsConverter(lambdaUnits, TroposphericModelUtils.MICRO_M).convert(lambda);
  61.         final double l2 = lambdaMicrometer  * lambdaMicrometer;
  62.         fLambda = 0.9650 + (0.0164 + 0.000228 / l2) / l2;

  64.     }

  66.     /** {@inheritDoc} */
  67.     @Override
  68.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
  69.                                        final PressureTemperatureHumidity weather,
  70.                                        final double[] parameters, final AbsoluteDate date) {

  72.         final double p = weather.getPressure();
  73.         final double t = weather.getTemperature();
  74.         final double e = weather.getWaterVaporPressure();

  76.         // beware since version 12.1 pressures are in Pa and not in hPa, hence the scaling has changed
  77.         final double Ah = 0.00002357 * p;
  78.         final double Aw = 0.00000141 * e;
  79.         final double K = 1.163 - 0.00968 * FastMath.cos(2 * point.getLatitude()) - 0.00104 * t + 0.0000001435 * p;
  80.         final double B = 1.084e-10 * p * t * K + 4.734e-12 * p * (p / t) * (2 * K) / (3 * K - 1);
  81.         final double flambda = getLaserFrequencyParameter();

  83.         final double fsite = getSiteFunctionValue(point);

  85.         final double sinE = FastMath.sin(trackingCoordinates.getElevation());
  86.         final double totalZenith       = (flambda / fsite) * (Ah + Aw + B) / (1.0   + B / ((Ah + Aw + B) * (1.0   + 0.01)));
  87.         final double totalElev         = (flambda / fsite) * (Ah + Aw + B) / (sinE  + B / ((Ah + Aw + B) * (sinE  + 0.01)));
  88.         final double hydrostaticZenith = (flambda / fsite) * (Ah +      B) / (1.0   + B / ((Ah +      B) * (1.0   + 0.01)));
  89.         final double hydrostaticElev   = (flambda / fsite) * (Ah +      B) / (sinE  + B / ((Ah +      B) * (sinE  + 0.01)));
  90.         return new TroposphericDelay(hydrostaticZenith, totalZenith - hydrostaticZenith,
  91.                                      hydrostaticElev,   totalElev   - hydrostaticElev);
  92.     }

  94.     /** {@inheritDoc} */
  95.     @Override
  96.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  97.                                                                                    final FieldGeodeticPoint<T> point,
  98.                                                                                    final FieldPressureTemperatureHumidity<T> weather,
  99.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {

  101.         final T p = weather.getPressure();
  102.         final T t = weather.getTemperature();
  103.         final T e = weather.getWaterVaporPressure();

  105.         // beware since version 12.1 pressures are in Pa and not in hPa, hence the scaling has changed
  106.         final T Ah = p.multiply(0.00002357);
  107.         final T Aw = e.multiply(0.00000141);
  108.         final T K = FastMath.cos(point.getLatitude().multiply(2.)).multiply(0.00968).negate().
  109.                     add(1.163).
  110.                     subtract(t.multiply(0.00104)).
  111.                     add(p.multiply(0.0000001435));
  112.         final T B = K.multiply(t.multiply(p).multiply(1.084e-10 )).
  113.                                add(K.multiply(2.).multiply(p.multiply(p).divide(t).multiply(4.734e-12)).divide(K.multiply(3.).subtract(1.)));
  114.         final double flambda = getLaserFrequencyParameter();

  116.         final T fsite = getSiteFunctionValue(point);

  118.         final T sinE = FastMath.sin(trackingCoordinates.getElevation());
  119.         final T one  = date.getField().getOne();
  120.         final T totalZenith       = fsite.divide(flambda).reciprocal().
  121.                                     multiply(B.add(Ah).add(Aw)).
  122.                                     divide(one.add(one.add(0.01).multiply(B.add(Ah).add(Aw)).divide(B).reciprocal()));
  123.         final T totalElev         = fsite.divide(flambda).reciprocal().
  124.                                     multiply(B.add(Ah).add(Aw)).
  125.                                     divide(sinE.add(sinE.add(0.01).multiply(B.add(Ah).add(Aw)).divide(B).reciprocal()));
  126.         final T hydrostaticZenith = fsite.divide(flambda).reciprocal().
  127.                                     multiply(B.add(Ah)).
  128.                                     divide(one.add(one.add(0.01).multiply(B.add(Ah)).divide(B).reciprocal()));
  129.         final T hydrostaticElev   = fsite.divide(flambda).reciprocal().
  130.                                     multiply(B.add(Ah)).
  131.                                     divide(sinE.add(sinE.add(0.01).multiply(B.add(Ah)).divide(B).reciprocal()));
  132.         return new FieldTroposphericDelay<>(hydrostaticZenith, totalZenith.subtract(hydrostaticZenith),
  133.                                             hydrostaticElev,   totalElev.subtract(hydrostaticElev));
  134.     }

  136.     /** {@inheritDoc} */
  137.     @Override
  138.     public List<ParameterDriver> getParametersDrivers() {
  139.         return Collections.emptyList();
  140.     }

  142.     /** Get the laser frequency parameter f(lambda).
  143.      * It is one for Ruby laser (lambda = 0.6943 micron)
  144.      * For infrared lasers, f(lambda) = 0.97966.
  145.      *
  146.      * @return the laser frequency parameter f(lambda).
  147.      */
  148.     private double getLaserFrequencyParameter() {
  149.         return fLambda;
  150.     }

  152.     /** Get the site parameter.
  153.      *
  154.      * @param point station location
  155.      * @return the site parameter.
  156.      */
  157.     private double getSiteFunctionValue(final GeodeticPoint point) {
  158.         return 1. - 0.0026 * FastMath.cos(2 * point.getLatitude()) - 0.00031 * 0.001 * point.getAltitude();
  159.     }

  161.     /** Get the site parameter.
  162.     *
  163.     * @param <T> type of the elements
  164.     * @param point station location
  165.     * @return the site parameter.
  166.     */
  167.     private <T extends CalculusFieldElement<T>> T getSiteFunctionValue(final FieldGeodeticPoint<T> point) {
  168.         return FastMath.cos(point.getLatitude().multiply(2)).multiply(0.0026).add(point.getAltitude().multiply(0.001).multiply(0.00031)).negate().add(1.);
  169.     }

  171. }
Created with JaCoCo 0.8.12.202403310830