MariniMurrayModel.java

/* Copyright 2011-2012 Space Applications Services
 * Licensed to CS Communication & Systèmes (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.models.earth;

import org.hipparchus.util.FastMath;

/** The Marini-Murray tropospheric delay model for laser ranging.
 *
 * @see "Marini, J.W., and C.W. Murray, correction of Laser Range Tracking Data for
 *      Atmospheric Refraction at Elevations Above 10 degrees, X-591-73-351, NASA GSFC, 1973"
 *
 * @author Joris Olympio
 */
public class MariniMurrayModel implements TroposphericModel {

    /** Serializable UID. */
    private static final long serialVersionUID = 8442906721207317886L;

    /** The temperature at the station, K. */
    private double T0;

    /** The atmospheric pressure, mbar. */
    private double P0;

    /** water vapor pressure at the laser site, mbar. */
    private double e0;

    /** Geodetic site latitude, radians. */
    private double latitude;

    /** Laser wavelength, micrometers. */
    private double lambda;

    /** Create a new Marini-Murray model for the troposphere using the given
     * environmental conditions.
     * @param t0 the temperature at the station, K
     * @param p0 the atmospheric pressure at the station, mbar
     * @param rh the humidity at the station, percent (50% -> 0.5)
     * @param latitude site latitude
     * @param lambda laser wavelength (c/f), nm
     */
    public MariniMurrayModel(final double t0, final double p0, final double rh, final double latitude, final double lambda) {
        this.T0 = t0;
        this.P0 = p0;

        this.e0 = getWaterVapor(rh);

        this.latitude = latitude;

        this.lambda = lambda * 1e-3;
    }

    /** Create a new Marini-Murray model using a standard atmosphere model.
     *
     * <ul>
     * <li>temperature: 20 degree Celsius
     * <li>pressure: 1013.25 mbar
     * <li>humidity: 50%
     * </ul>
     *
     * @param latitude site latitude
     * @param frequency laser frequency, Hz
     *
     * @return a Saastamoinen model with standard environmental values
     */
    public static MariniMurrayModel getStandardModel(final double latitude, final double frequency) {
        return new MariniMurrayModel(273.15 + 20, 1013.25, 0.5, latitude, frequency);
    }

    @Override
    public double pathDelay(final double elevation, final double height) {
        final double A = 0.002357 * P0 + 0.000141 * e0;
        final double K = 1.163 * 0.00968 * FastMath.cos(2 * latitude) - 0.00104 * T0 + 0.00001435 * P0;
        final double B = (1.084 * 1e-8) * P0 * T0 * K + (4.734 * 1e-8) * P0 * (P0 / T0) * (2 * K) / (3 * K - 1);
        final double flambda = getLaserFrequencyParameter();

        final double fsite = getSiteFunctionValue(height / 1000.);

        final double sinE = FastMath.sin(elevation);
        final double dR = (flambda / fsite) * (A + B) / (sinE + B / ((A + B) * (sinE + 0.01)) );
        return dR;
    }

    /** Get the laser frequency parameter f(lambda).
     * It is one for Ruby laser (lambda = 0.6943 micron)
     * For infrared lasers, f(lambda) = 0.97966.
     *
     * @return the laser frequency parameter f(lambda).
     */
    private double getLaserFrequencyParameter() {
        return 0.9650 + 0.0164 * FastMath.pow(lambda, -2) + 0.000228 * FastMath.pow(lambda, -4);
    }

    /** Get the laser frequency parameter f(lambda).
     *
     * @param height height above the geoid, km
     * @return the laser frequency parameter f(lambda).
     */
    private double getSiteFunctionValue(final double height) {
        return 1. - 0.0026 * FastMath.cos(2 * latitude) - 0.00031 * height;
    }

    /** Get the water vapor.
     * The water vapor model is the one of Giacomo and Davis as indicated in IERS TN 32, chap. 9.
     *
     * See: Giacomo, P., Equation for the dertermination of the density of moist air, Metrologia, V. 18, 1982
     *
     * @param rh relative humidity, in percent.
     * @return the water vapor, in mbar (1 mbar = 100 Pa).
     */
    private double getWaterVapor(final double rh) {

        // saturation water vapor, equation (3) of reference paper, in mbar
        // with amended 1991 values (see reference paper)
        final double es = 0.01 * FastMath.exp((1.2378847 * 1e-5) * T0 * T0 -
                                              (1.9121316 * 1e-2) * T0 +
                                              33.93711047 -
                                              (6.3431645 * 1e3) * 1. / T0);

        // enhancement factor, equation (4) of reference paper
        final double fw = 1.00062 + (3.14 * 1e-6) * P0 + (5.6 * 1e-7) * FastMath.pow(T0 - 273.15, 2);

        final double e = rh * fw * es;
        return e;
    }
}