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