1   /* Copyright 2013 Applied Defense Solutions, Inc.
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;
18  
19  import org.hipparchus.optim.MaxEval;
20  import org.hipparchus.optim.nonlinear.scalar.GoalType;
21  import org.hipparchus.optim.univariate.BrentOptimizer;
22  import org.hipparchus.optim.univariate.SearchInterval;
23  import org.hipparchus.optim.univariate.UnivariateObjectiveFunction;
24  import org.hipparchus.util.FastMath;
25  import org.orekit.models.AtmosphericRefractionModel;
26  import org.orekit.models.earth.troposphere.iturp834.ITURP834PathDelay;
27  
28  /** Implementation of refraction model for Earth exponential atmosphere based on ITU-R P.834 recommendation.
29   * <p>
30   * This class implements the ray bending part, i.e. section 1 of the recommendation.
31   * The excess radio path length part of the model, i.e. section 6 of the recommendation,
32   * is implemented in the {@link ITURP834PathDelay} class.
33   * </p>
34   *
35   * @author Thierry Ceolin
36   * @since 7.1
37   * @see <a href="https://www.itu.int/rec/R-REC-P.834/en">P.834 : Effects of tropospheric refraction on radiowave propagation</a>
38   */
39  
40  public class ITURP834AtmosphericRefraction implements AtmosphericRefractionModel {
41  
42      /** Altitude conversion factor. */
43      private static final double KM_TO_M = 1000.0;
44  
45      /** Coefficients conversion factor. */
46      private static final double INV_DEG_TO_INV_RAD = 180.0 / FastMath.PI;
47  
48      /** Default a coefficients to compute refractive index for a typical atmosphere. */
49      private static final double DEFAULT_CORRECTION_ACOEF = 0.000315;
50  
51      /** Default b coefficients to compute refractive index for a typical atmosphere. */
52      private static final double DEFAULT_CORRECTION_BCOEF = 0.1361 / KM_TO_M;
53  
54      /** Earth ray as defined in ITU-R P.834-9 (m). */
55      private static final double EARTH_RAY = 6370.0 * KM_TO_M;
56  
57      /** Default coefficients array for Tau function (formula number 9).
58       * The coefficients have been converted to SI units
59       */
60      private static final double[] CCOEF = {
61          INV_DEG_TO_INV_RAD * 1.314,  INV_DEG_TO_INV_RAD * 0.6437,  INV_DEG_TO_INV_RAD * 0.02869,
62          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,
63          INV_DEG_TO_INV_RAD * 0.008583 / (KM_TO_M * KM_TO_M)
64      };
65  
66      /** Default coefficients array for TauZero function (formula number 14).
67       * The coefficients have been converted to SI units
68       */
69      private static final double[] CCOEF0 = {
70          INV_DEG_TO_INV_RAD * 1.728, INV_DEG_TO_INV_RAD * 0.5411, INV_DEG_TO_INV_RAD * 0.03723,
71          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,
72          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)
73      };
74  
75      /** Serializable UID. */
76      private static final long serialVersionUID = 20160118L;
77  
78      /** station altitude (m). */
79      private final double altitude;
80  
81      /** minimal elevation angle for the station (rad). */
82      private final double thetamin;
83  
84      /** minimal elevation angle under free-space propagation (rad). */
85      private final double theta0;
86  
87      /** elevation where elevation+refraction correction is minimal (near inequality formula number 11 validity domain). */
88      private final double elev_star;
89  
90      /** refraction correction value where elevation+refraction correction is minimal (near inequality 11 validity domain). */
91      private final double refrac_star;
92  
93      /** Creates a new default instance.
94       * @param altitude altitude of the ground station from which measurement is performed (m)
95       */
96      public ITURP834AtmosphericRefraction(final double altitude) {
97          this.altitude = altitude;
98          thetamin = getMinimalElevation(altitude);
99          theta0   = thetamin - getTau(thetamin);
100 
101         final double rel = 1.e-5;
102         final double abs = 1.e-10;
103         final BrentOptimizer optimizer = new BrentOptimizer(rel, abs);
104 
105         // Call optimizer
106         elev_star = optimizer.optimize(new MaxEval(200),
107                                        new UnivariateObjectiveFunction(e -> e + getBaseRefraction(e)),
108                                        GoalType.MINIMIZE,
109                                        new SearchInterval(-FastMath.PI / 30., FastMath.PI / 4)).getPoint();
110         refrac_star = getBaseRefraction(elev_star);
111     }
112 
113     /** Compute the refractive index correction in the case of a typical atmosphere.
114      * ITU-R P.834-9, formula number 8, page 3
115      * @param alt altitude of the station at the Earth surface (m)
116      * @return the refractive index
117      */
118     private double getRefractiveIndex(final double alt) {
119 
120         return 1.0 + DEFAULT_CORRECTION_ACOEF * FastMath.exp(-DEFAULT_CORRECTION_BCOEF * alt);
121     }
122 
123     /** Compute the minimal elevation angle for a station.
124      * ITU-R P.834-9, formula number 10, page 3
125      * @param alt altitude of the station at the Earth surface (m)
126      * @return the minimal elevation angle (rad)
127      */
128     private double getMinimalElevation(final double alt) {
129 
130         return -FastMath.acos( EARTH_RAY / (EARTH_RAY + alt) * getRefractiveIndex(0.0) / getRefractiveIndex(alt));
131     }
132 
133 
134     /** Compute the refraction correction in the case of a reference atmosphere.
135      * ITU-R P.834-9, formula number 9, page 3
136      * @param elevation elevation angle (rad)
137      * @return the refraction correction angle (rad)
138      */
139     private double getTau(final double elevation) {
140 
141         final double eld = FastMath.toDegrees(elevation);
142         final double tmp0 = CCOEF[0] + (CCOEF[1] + CCOEF[2] * eld) * eld;
143         final double tmp1 = altitude * (CCOEF[3] + (CCOEF[4] + CCOEF[5] * eld) * eld);
144         final double tmp2 = altitude * altitude * CCOEF[6];
145         return 1.0 / (tmp0 + tmp1 + tmp2);
146     }
147 
148 
149     /** Compute the refraction correction in the case of a reference atmosphere.
150      * ITU-R P.834-9, formula number 14, page 4
151      * @param elevationZero elevation angle (rad)
152      * @return the refraction correction angle (rad)
153      */
154 
155     private double getTauZero(final double elevationZero) {
156 
157         final double eld = FastMath.toDegrees(elevationZero);
158         final double tmp0 = CCOEF0[0] + (CCOEF0[1] + CCOEF0[2] * eld) * eld;
159         final double tmp1 = altitude * (CCOEF0[3] + (CCOEF0[4] + CCOEF0[5] * eld) * eld);
160         final double tmp2 = altitude * altitude * (CCOEF0[6] + CCOEF0[7] * eld);
161         return 1.0 / (tmp0 + tmp1 + tmp2);
162     }
163 
164     /** Compute the refraction correction in the case of a reference atmosphere without validity domain.
165      * The computation is done even if the inequality (formula number 11) is not verified
166      * ITU-R P.834-9, formula number 14, page 3
167      * @param elevation elevation angle (rad)
168      * @return the refraction correction angle (rad)
169      */
170     private double getBaseRefraction(final double elevation) {
171         return getTauZero(elevation);
172     }
173 
174     /** Get the station minimal elevation angle.
175      * @return the minimal elevation angle (rad)
176      */
177     public double getThetaMin() {
178         return thetamin;
179     }
180 
181     /** Get the station elevation angle under free-space propagation .
182      * @return the elevation angle under free-space propagation (rad)
183      */
184     public double getTheta0() {
185         return theta0;
186     }
187 
188     /** {@inheritDoc} */
189     @Override
190     public double getRefraction(final double elevation) {
191         if (elevation < elev_star ) {
192             return refrac_star;
193         }
194         // The validity of the formula is extended for negative elevation,
195         // ensuring that the refraction correction angle doesn't make visible a satellite with a too negative elevation
196         // elev_star is used instead of thetam (minimal elevation angle).
197         return getTauZero(elevation);
198     }
199 
200 }