1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (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.util.FastMath;
20  import org.junit.jupiter.api.Assertions;
21  import org.junit.jupiter.api.BeforeEach;
22  import org.junit.jupiter.api.Test;
23  import org.orekit.Utils;
24  import org.orekit.bodies.GeodeticPoint;
25  import org.orekit.bodies.OneAxisEllipsoid;
26  import org.orekit.frames.FramesFactory;
27  import org.orekit.frames.TopocentricFrame;
28  import org.orekit.utils.Constants;
29  import org.orekit.utils.IERSConventions;
30  
31  public class EarthITURP834AtmosphericRefractionTest {
32  
33      private final double onehundredth = 1e-2;
34      private final double twohundredth = 2e-2;
35      private final double onethousandth = 1e-3;
36      private final double epsilon = 1e-15;
37  
38      // Table (ref. Astronomical Refraction, Michael E. Thomas and Richard I. Joseph)
39      //       (JOHNS HOPKINS APL TECHNICAL DIGEST, VOLUME 17, NUMBER 3 (1996))
40      // elevation (deg)
41      private final double[] ref_elevation  = new double[] {0.00, 0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 4.50, 5.00,
42                                                            6.00, 7.00, 8.00, 9.00, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0,
43                                                            25.0, 30.0, 35.0, 50.0, 55.0, 60.0, 65.0, 70.0, 80.0, 90.0};
44      // refraction correction angle (in arcminutes)
45      private final double[] ref_refraction = new double[] {34.5, 31.4, 28.7, 26.4, 24.3, 22.5, 20.9, 19.5, 18.3, 17.2, 16.1, 15.2, 14.4, 10.7, 9.90,
46                                                            8.50, 7.40, 6.60, 5.90, 5.30, 4.90, 4.50, 4.10, 3.80, 3.60, 3.30, 3.10, 2.90, 2.80, 2.60,
47                                                            2.10, 1.70, 1.40, 0.80, 0.70, 0.60, 0.50, 0.40, 0.20, 0.00};
48  
49  
50      // Kiruna-2 ESTRACK Station (Sweden)
51      private TopocentricFrame stationk;
52      private final String namek = "Kiruna-2";
53      // Hartebeesthoek IGS Station (South Africa)
54      // lowest elevation angle that verify inequality number 11 : theta0 = -1.039 degree;
55      private TopocentricFrame stationh;
56      private final String nameh = "Hartebeesthoek";
57      // Everest Fake Station (China/Nepal)
58      private TopocentricFrame statione;
59      private final String namee = "Everest";
60      // Dead Sea Fake Station (Israel)
61      private TopocentricFrame stationd;
62      private final String named = "Dead Sea";
63      // Altitude0 Fake Station ()
64      private TopocentricFrame stationa;
65      private final String namea = "Alt0";
66  
67      @BeforeEach
68      public void setUp() throws Exception {
69          Utils.setDataRoot("regular-data:potential:tides");
70          IERSConventions  conventions = IERSConventions.IERS_2010;
71          OneAxisEllipsoid earth       = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
72                                                              Constants.WGS84_EARTH_FLATTENING,
73                                                              FramesFactory.getITRF(conventions, true));
74  
75          // Kiruna-2 (Sweden)
76          final GeodeticPoint kir = new GeodeticPoint(FastMath.toRadians(67.858428),
77                                                     FastMath.toRadians(20.966880),
78                                                     385.8);
79  
80          // Hartebeesthoek (South Africa)
81          final GeodeticPoint har = new GeodeticPoint(FastMath.toRadians(-24.110243),
82                                                      FastMath.toRadians(27.685308),
83                                                      1415.821);
84  
85          // Everest (fake station)
86          final GeodeticPoint eve = new GeodeticPoint(FastMath.toRadians(27.988333),
87                                                      FastMath.toRadians(86.991944),
88                                                      8848.0);
89  
90          // Dead Sea (fake station)
91          final GeodeticPoint des = new GeodeticPoint(FastMath.toRadians(31.500000),
92                                                      FastMath.toRadians(35.500000),
93                                                      -422.0);
94  
95          // Alt0 (fake station)
96          final GeodeticPoint alt = new GeodeticPoint(FastMath.toRadians(31.500000),
97                                                      FastMath.toRadians(35.500000),
98                                                      0.0);
99          stationk = new TopocentricFrame(earth, kir, namek);
100         stationh = new TopocentricFrame(earth, har, nameh);
101         statione = new TopocentricFrame(earth, eve, namee);
102         stationd = new TopocentricFrame(earth, des, named);
103         stationa = new TopocentricFrame(earth, alt, namea);
104     }
105 
106 
107     @Test
108     public void testEarthITU453AtmosphereRefractionHighest() {
109 
110         // elevation angle of the space station under free-space propagation conditions
111         final double elevation = FastMath.toRadians(2.0);
112 
113         // Station altitude
114         final double altitude = statione.getPoint().getAltitude();
115         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
116 
117         // refraction correction in degrees
118         double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
119         Assertions.assertEquals(0.11458177523385392, refraction, epsilon);
120     }
121 
122     @Test
123     public void testEarthITU453AtmosphereRefractionLowest() {
124 
125         // elevation angle of the space station under free-space propagation conditions
126         final double elevation = FastMath.toRadians(2.0);
127 
128         // Station altitude
129         final double altitude = stationd.getPoint().getAltitude();
130         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
131 
132         // refraction correction in degrees
133         double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
134         Assertions.assertEquals(0.3550620274090111, refraction, epsilon);
135     }
136 
137     @Test
138     public void testEarthITU453AtmosphereRefraction2degree() {
139 
140         // elevation angle of the space station under free-space propagation conditions
141         final double elevation = FastMath.toRadians(2.0);
142 
143         // Station altitude
144         final double altitude = stationk.getPoint().getAltitude();
145         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
146 
147         // refraction correction in degrees
148         final double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
149         Assertions.assertEquals(refraction, 0.32, onehundredth);
150 
151         final double thetamin = FastMath.toDegrees(modelTropo.getThetaMin());
152         Assertions.assertEquals(-0.5402509318003884, thetamin, epsilon);
153         final double theta0 = FastMath.toDegrees(modelTropo.getTheta0());
154         Assertions.assertEquals(-1.4959064751203384, theta0, epsilon);
155 
156     }
157 
158     @Test
159     public void testEarthITU453AtmosphereRefraction4degree() {
160 
161         // elevation angle of the space station under free-space propagation conditions
162         final double elevation = FastMath.toRadians(4.0);
163 
164         // Station altitude
165         final double altitude = stationk.getPoint().getAltitude();
166         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
167 
168         // refraction correction in degrees
169         double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
170         Assertions.assertEquals(0.21, refraction, onehundredth);
171     }
172 
173     @Test
174     public void testEarthITU453AtmosphereRefraction10degree() {
175 
176         // elevation angle of the space station under free-space propagation conditions
177         final double elevation = FastMath.toRadians(10.0);
178 
179         // Station altitude
180         final double altitude = stationk.getPoint().getAltitude();
181         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
182 
183         // refraction correction in degrees
184         double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
185         Assertions.assertEquals(0.10, refraction, twohundredth);
186     }
187 
188     @Test
189     public void testEarthITU453AtmosphereRefraction30degree() {
190 
191         // elevation angle of the space station under free-space propagation conditions
192         final double elevation = FastMath.toRadians(30.0);
193 
194         // Station altitude
195         final double altitude = stationk.getPoint().getAltitude();
196         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
197 
198         // refraction correction in degrees
199         double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
200         Assertions.assertEquals(0.02, refraction, onehundredth);
201     }
202 
203     @Test
204     public void testEarthITU453AtmosphereRefraction90degree() {
205 
206         // elevation angle of the space station under free-space propagation conditions
207         final double elevation = FastMath.toRadians(90.0);
208 
209         // Station altitude
210         final double altitude = stationk.getPoint().getAltitude();
211         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
212 
213         // refraction correction in degrees
214         double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
215         Assertions.assertEquals(0.002, refraction, onethousandth);
216 
217     }
218     @Test
219     public void testEarthITU453AtmosphereRefractionminusdegree() {
220 
221         // elevation angle of the space station under free-space propagation conditions
222         final double elevation = FastMath.toRadians(-10.);
223 
224         // Station altitude
225         final double altitude = stationh.getPoint().getAltitude();
226         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
227 
228         // refraction correction in degrees
229         double refraction = FastMath.toDegrees(modelTropo.getRefraction(elevation));
230         Assertions.assertEquals(1.7367073234643113, refraction, onethousandth);
231     }
232 
233     @Test
234     public void testEarthITU453AtmosphereRefractiontable() {
235 
236         // Station altitude
237         final double altitude = stationa.getPoint().getAltitude();
238         ITURP834AtmosphericRefraction modelTropo = new ITURP834AtmosphericRefraction(altitude);
239 
240         for (int itab=0; itab<40; itab++) {
241             // elevation angle of the space station under free-space propagation conditions
242             final double elevation = FastMath.toRadians(ref_elevation[itab]);
243 
244             // refraction correction in arcminutes
245             final double refraction = 60.0 * FastMath.toDegrees(modelTropo.getRefraction(elevation));
246             Assertions.assertEquals(ref_refraction[itab], refraction, 2.1);
247         }
248     }
249 }