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.models.earth.weather.PressureTemperatureHumidityProvider;
  29. import org.orekit.time.AbsoluteDate;
  30. import org.orekit.time.FieldAbsoluteDate;
  31. import org.orekit.utils.FieldTrackingCoordinates;
  32. import org.orekit.utils.ParameterDriver;
  33. import org.orekit.utils.TrackingCoordinates;
  34. import org.orekit.utils.units.Unit;
  35. import org.orekit.utils.units.UnitsConverter;

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

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

  51.     /** Provider for pressure, temperature and humidity.
  52.      * @since 13.0
  53.      */
  54.     private final PressureTemperatureHumidityProvider pthProvider;

  56.     /** Create a new Marini-Murray model for the troposphere.
  57.      * @param lambda laser wavelength
  58.      * @param lambdaUnits units in which {@code lambda} is given
  59.      * @param pthProvider provider for pressure, temperature and humidity
  60.      * @see TroposphericModelUtils#MICRO_M
  61.      * @see TroposphericModelUtils#NANO_M
  62.      * @since 12.1
  63.      * */
  64.     public MariniMurray(final double lambda, final Unit lambdaUnits, final PressureTemperatureHumidityProvider pthProvider) {

  66.         this.pthProvider = pthProvider;

  68.         // compute laser frequency parameter
  69.         final double lambdaMicrometer = new UnitsConverter(lambdaUnits, TroposphericModelUtils.MICRO_M).convert(lambda);
  70.         final double l2 = lambdaMicrometer  * lambdaMicrometer;
  71.         this.fLambda = 0.9650 + (0.0164 + 0.000228 / l2) / l2;

  73.     }

  75.     /** {@inheritDoc} */
  76.     @Override
  77.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
  78.                                        final double[] parameters, final AbsoluteDate date) {

  80.         final PressureTemperatureHumidity weather = pthProvider.getWeatherParameters(point, date);
  81.         final double p = weather.getPressure();
  82.         final double t = weather.getTemperature();
  83.         final double e = weather.getWaterVaporPressure();

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

  92.         final double fsite = getSiteFunctionValue(point);

  94.         final double sinE = FastMath.sin(trackingCoordinates.getElevation());
  95.         final double totalZenith       = (flambda / fsite) * (Ah + Aw + B) / (1.0   + B / ((Ah + Aw + B) * (1.0   + 0.01)));
  96.         final double totalElev         = (flambda / fsite) * (Ah + Aw + B) / (sinE  + B / ((Ah + Aw + B) * (sinE  + 0.01)));
  97.         final double hydrostaticZenith = (flambda / fsite) * (Ah +      B) / (1.0   + B / ((Ah +      B) * (1.0   + 0.01)));
  98.         final double hydrostaticElev   = (flambda / fsite) * (Ah +      B) / (sinE  + B / ((Ah +      B) * (sinE  + 0.01)));
  99.         return new TroposphericDelay(hydrostaticZenith, totalZenith - hydrostaticZenith,
  100.                                      hydrostaticElev,   totalElev   - hydrostaticElev);
  101.     }

  103.     /** {@inheritDoc} */
  104.     @Override
  105.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  106.                                                                                    final FieldGeodeticPoint<T> point,
  107.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {

  109.         final FieldPressureTemperatureHumidity<T> weather = pthProvider.getWeatherParameters(point, date);
  110.         final T p = weather.getPressure();
  111.         final T t = weather.getTemperature();
  112.         final T e = weather.getWaterVaporPressure();

  114.         // beware since version 12.1 pressures are in Pa and not in hPa, hence the scaling has changed
  115.         final T Ah = p.multiply(0.00002357);
  116.         final T Aw = e.multiply(0.00000141);
  117.         final T K = FastMath.cos(point.getLatitude().multiply(2.)).multiply(0.00968).negate().
  118.                     add(1.163).
  119.                     subtract(t.multiply(0.00104)).
  120.                     add(p.multiply(0.0000001435));
  121.         final T B = K.multiply(t.multiply(p).multiply(1.084e-10 )).
  122.                                add(K.multiply(2.).multiply(p.multiply(p).divide(t).multiply(4.734e-12)).divide(K.multiply(3.).subtract(1.)));
  123.         final double flambda = getLaserFrequencyParameter();

  125.         final T fsite = getSiteFunctionValue(point);

  127.         final T sinE = FastMath.sin(trackingCoordinates.getElevation());
  128.         final T one  = date.getField().getOne();
  129.         final T totalZenith       = fsite.divide(flambda).reciprocal().
  130.                                     multiply(B.add(Ah).add(Aw)).
  131.                                     divide(one.add(one.add(0.01).multiply(B.add(Ah).add(Aw)).divide(B).reciprocal()));
  132.         final T totalElev         = fsite.divide(flambda).reciprocal().
  133.                                     multiply(B.add(Ah).add(Aw)).
  134.                                     divide(sinE.add(sinE.add(0.01).multiply(B.add(Ah).add(Aw)).divide(B).reciprocal()));
  135.         final T hydrostaticZenith = fsite.divide(flambda).reciprocal().
  136.                                     multiply(B.add(Ah)).
  137.                                     divide(one.add(one.add(0.01).multiply(B.add(Ah)).divide(B).reciprocal()));
  138.         final T hydrostaticElev   = fsite.divide(flambda).reciprocal().
  139.                                     multiply(B.add(Ah)).
  140.                                     divide(sinE.add(sinE.add(0.01).multiply(B.add(Ah)).divide(B).reciprocal()));
  141.         return new FieldTroposphericDelay<>(hydrostaticZenith, totalZenith.subtract(hydrostaticZenith),
  142.                                             hydrostaticElev,   totalElev.subtract(hydrostaticElev));
  143.     }

  145.     /** {@inheritDoc} */
  146.     @Override
  147.     public List<ParameterDriver> getParametersDrivers() {
  148.         return Collections.emptyList();
  149.     }

  151.     /** Get the laser frequency parameter f(lambda).
  152.      * It is one for Ruby laser (lambda = 0.6943 micron)
  153.      * For infrared lasers, f(lambda) = 0.97966.
  154.      *
  155.      * @return the laser frequency parameter f(lambda).
  156.      */
  157.     private double getLaserFrequencyParameter() {
  158.         return fLambda;
  159.     }

  161.     /** Get the site parameter.
  162.      *
  163.      * @param point station location
  164.      * @return the site parameter.
  165.      */
  166.     private double getSiteFunctionValue(final GeodeticPoint point) {
  167.         return 1. - 0.0026 * FastMath.cos(2 * point.getLatitude()) - 0.00031 * 0.001 * point.getAltitude();
  168.     }

  170.     /** Get the site parameter.
  171.     *
  172.     * @param <T> type of the elements
  173.     * @param point station location
  174.     * @return the site parameter.
  175.     */
  176.     private <T extends CalculusFieldElement<T>> T getSiteFunctionValue(final FieldGeodeticPoint<T> point) {
  177.         return FastMath.cos(point.getLatitude().multiply(2)).multiply(0.0026).add(point.getAltitude().multiply(0.001).multiply(0.00031)).negate().add(1.);
  178.     }

  180. }
Created with JaCoCo 0.8.12.202403310830