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