ITURP834AtmosphericRefraction.java

/* Copyright 2013 Applied Defense Solutions, Inc.
 * 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.optim.MaxEval;
import org.hipparchus.optim.nonlinear.scalar.GoalType;
import org.hipparchus.optim.univariate.BrentOptimizer;
import org.hipparchus.optim.univariate.SearchInterval;
import org.hipparchus.optim.univariate.UnivariateObjectiveFunction;
import org.hipparchus.util.FastMath;
import org.orekit.models.AtmosphericRefractionModel;
import org.orekit.models.earth.troposphere.iturp834.ITURP834PathDelay;

/** Implementation of refraction model for Earth exponential atmosphere based on ITU-R P.834 recommendation.
 * <p>
 * This class implements the ray bending part, i.e. section 1 of the recommendation.
 * The excess radio path length part of the model, i.e. section 6 of the recommendation,
 * is implemented in the {@link ITURP834PathDelay} class.
 * </p>
 *
 * @author Thierry Ceolin
 * @since 7.1
 * @see <a href="https://www.itu.int/rec/R-REC-P.834/en">P.834 : Effects of tropospheric refraction on radiowave propagation</a>
 */

public class ITURP834AtmosphericRefraction implements AtmosphericRefractionModel {

    /** Altitude conversion factor. */
    private static final double KM_TO_M = 1000.0;

    /** Coefficients conversion factor. */
    private static final double INV_DEG_TO_INV_RAD = 180.0 / FastMath.PI;

    /** Default a coefficients to compute refractive index for a typical atmosphere. */
    private static final double DEFAULT_CORRECTION_ACOEF = 0.000315;

    /** Default b coefficients to compute refractive index for a typical atmosphere. */
    private static final double DEFAULT_CORRECTION_BCOEF = 0.1361 / KM_TO_M;

    /** Earth ray as defined in ITU-R P.834-9 (m). */
    private static final double EARTH_RAY = 6370.0 * KM_TO_M;

    /** Default coefficients array for Tau function (formula number 9).
     * The coefficients have been converted to SI units
     */
    private static final double[] CCOEF = {
        INV_DEG_TO_INV_RAD * 1.314,  INV_DEG_TO_INV_RAD * 0.6437,  INV_DEG_TO_INV_RAD * 0.02869,
        INV_DEG_TO_INV_RAD * 0.2305 / KM_TO_M, INV_DEG_TO_INV_RAD * 0.09428 / KM_TO_M, INV_DEG_TO_INV_RAD * 0.01096 / KM_TO_M,
        INV_DEG_TO_INV_RAD * 0.008583 / (KM_TO_M * KM_TO_M)
    };

    /** Default coefficients array for TauZero function (formula number 14).
     * The coefficients have been converted to SI units
     */
    private static final double[] CCOEF0 = {
        INV_DEG_TO_INV_RAD * 1.728, INV_DEG_TO_INV_RAD * 0.5411, INV_DEG_TO_INV_RAD * 0.03723,
        INV_DEG_TO_INV_RAD * 0.1815 / KM_TO_M, INV_DEG_TO_INV_RAD * 0.06272 / KM_TO_M, INV_DEG_TO_INV_RAD * 0.011380 / KM_TO_M,
        INV_DEG_TO_INV_RAD * 0.01727 / (KM_TO_M * KM_TO_M), INV_DEG_TO_INV_RAD * 0.008288 / (KM_TO_M * KM_TO_M)
    };

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

    /** station altitude (m). */
    private final double altitude;

    /** minimal elevation angle for the station (rad). */
    private final double thetamin;

    /** minimal elevation angle under free-space propagation (rad). */
    private final double theta0;

    /** elevation where elevation+refraction correction is minimal (near inequality formula number 11 validity domain). */
    private final double elev_star;

    /** refraction correction value where elevation+refraction correction is minimal (near inequality 11 validity domain). */
    private final double refrac_star;

    /** Creates a new default instance.
     * @param altitude altitude of the ground station from which measurement is performed (m)
     */
    public ITURP834AtmosphericRefraction(final double altitude) {
        this.altitude = altitude;
        thetamin = getMinimalElevation(altitude);
        theta0   = thetamin - getTau(thetamin);

        final double rel = 1.e-5;
        final double abs = 1.e-10;
        final BrentOptimizer optimizer = new BrentOptimizer(rel, abs);

        // Call optimizer
        elev_star = optimizer.optimize(new MaxEval(200),
                                       new UnivariateObjectiveFunction(e -> e + getBaseRefraction(e)),
                                       GoalType.MINIMIZE,
                                       new SearchInterval(-FastMath.PI / 30., FastMath.PI / 4)).getPoint();
        refrac_star = getBaseRefraction(elev_star);
    }

    /** Compute the refractive index correction in the case of a typical atmosphere.
     * ITU-R P.834-9, formula number 8, page 3
     * @param alt altitude of the station at the Earth surface (m)
     * @return the refractive index
     */
    private double getRefractiveIndex(final double alt) {

        return 1.0 + DEFAULT_CORRECTION_ACOEF * FastMath.exp(-DEFAULT_CORRECTION_BCOEF * alt);
    }

    /** Compute the minimal elevation angle for a station.
     * ITU-R P.834-9, formula number 10, page 3
     * @param alt altitude of the station at the Earth surface (m)
     * @return the minimal elevation angle (rad)
     */
    private double getMinimalElevation(final double alt) {

        return -FastMath.acos( EARTH_RAY / (EARTH_RAY + alt) * getRefractiveIndex(0.0) / getRefractiveIndex(alt));
    }


    /** Compute the refraction correction in the case of a reference atmosphere.
     * ITU-R P.834-9, formula number 9, page 3
     * @param elevation elevation angle (rad)
     * @return the refraction correction angle (rad)
     */
    private double getTau(final double elevation) {

        final double eld = FastMath.toDegrees(elevation);
        final double tmp0 = CCOEF[0] + (CCOEF[1] + CCOEF[2] * eld) * eld;
        final double tmp1 = altitude * (CCOEF[3] + (CCOEF[4] + CCOEF[5] * eld) * eld);
        final double tmp2 = altitude * altitude * CCOEF[6];
        return 1.0 / (tmp0 + tmp1 + tmp2);
    }


    /** Compute the refraction correction in the case of a reference atmosphere.
     * ITU-R P.834-9, formula number 14, page 4
     * @param elevationZero elevation angle (rad)
     * @return the refraction correction angle (rad)
     */

    private double getTauZero(final double elevationZero) {

        final double eld = FastMath.toDegrees(elevationZero);
        final double tmp0 = CCOEF0[0] + (CCOEF0[1] + CCOEF0[2] * eld) * eld;
        final double tmp1 = altitude * (CCOEF0[3] + (CCOEF0[4] + CCOEF0[5] * eld) * eld);
        final double tmp2 = altitude * altitude * (CCOEF0[6] + CCOEF0[7] * eld);
        return 1.0 / (tmp0 + tmp1 + tmp2);
    }

    /** Compute the refraction correction in the case of a reference atmosphere without validity domain.
     * The computation is done even if the inequality (formula number 11) is not verified
     * ITU-R P.834-9, formula number 14, page 3
     * @param elevation elevation angle (rad)
     * @return the refraction correction angle (rad)
     */
    private double getBaseRefraction(final double elevation) {
        return getTauZero(elevation);
    }

    /** Get the station minimal elevation angle.
     * @return the minimal elevation angle (rad)
     */
    public double getThetaMin() {
        return thetamin;
    }

    /** Get the station elevation angle under free-space propagation .
     * @return the elevation angle under free-space propagation (rad)
     */
    public double getTheta0() {
        return theta0;
    }

    /** {@inheritDoc} */
    @Override
    public double getRefraction(final double elevation) {
        if (elevation < elev_star ) {
            return refrac_star;
        }
        // The validity of the formula is extended for negative elevation,
        // ensuring that the refraction correction angle doesn't make visible a satellite with a too negative elevation
        // elev_star is used instead of thetam (minimal elevation angle).
        return getTauZero(elevation);
    }

}