Line data Source code
1 : /** \file astroDynamics.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Various utilities for astrodynamics.
4 : * \ingroup astrofiles
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2018 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef astroDynamics_hpp
28 : #define astroDynamics_hpp
29 :
30 : #include <cmath>
31 :
32 : #include "sofa.hpp"
33 :
34 : #include "../math/constants.hpp"
35 : #include "../math/geo.hpp"
36 :
37 : namespace mx
38 : {
39 : namespace astro
40 : {
41 :
42 : /// Breaks down a decimal day into hours, minutes and decimal point seconds.
43 : /** Assumes 86400 seconds per day, ignoring leap seconds.
44 : *
45 : * \tparam realT the floating point type for calculations
46 : *
47 : * \retval 0 on success
48 : * \retval -1 on failure.
49 : */
50 : template <typename realT>
51 : int hrMinSecFromDay( int &Dy, ///< [out] the integer day
52 : int &hr, ///< [out] the integer hour
53 : int &min, ///< [out] the integer minute
54 : realT &sec, ///< [out] the integer second
55 : realT day ///< [in] the decimal day to be broken down
56 : )
57 : {
58 : Dy = (int)day;
59 : hr = ( day - (realT)Dy ) * 24.0;
60 : min = ( day - (realT)Dy - ( (realT)hr ) / 24.0 ) * 60. * 24.0;
61 : sec = ( day - (realT)Dy - ( (realT)hr ) / 24.0 - ( (realT)min ) / ( 60.0 * 24.0 ) ) * 3600.0 * 24.0;
62 :
63 : return 0;
64 : }
65 :
66 : /// Returns Greenwich Mean Sidereal Time for a given UTC time.
67 : /** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
68 : *
69 : * \retval 0 on success
70 : */
71 : template <typename realT>
72 : int getGMST( realT &GMST, ///< [out] GMST in hours
73 : int Yr, ///< [in] integer year
74 : int Mo, ///< [in] integer month (1-12)
75 : int Dy, ///< [in] integer day (1-31)
76 : int hr, ///< [in] integer hour (0-23)
77 : int min, ///< [in] integer minute (0-59)
78 : realT sec ///< [in] decimal second (0-59)
79 : )
80 : {
81 : int rv;
82 : double utc0, utc1, tai0, tai1, tt0, tt1, ut10, ut11;
83 :
84 : // Get UTC time in SOFA 2-double format
85 : rv = sofa::iauDtf2d( "UTC", Yr, Mo, Dy, hr, min, static_cast<double>( sec ), &utc0, &utc1 );
86 : if( rv < 0 )
87 : {
88 : /// \todo handle SOFA error.
89 : return -1;
90 : }
91 :
92 : // Convert UTC to TAI
93 : rv = sofa::iauUtctai( utc0, utc1, &tai0, &tai1 );
94 : if( rv < 0 )
95 : {
96 : /// \todo handle SOFA error.
97 : return -1;
98 : }
99 :
100 : // Convert TAI to TT
101 : rv = sofa::iauTaitt( tai0, tai1, &tt0, &tt1 );
102 : if( rv < 0 )
103 : {
104 : /// \todo handle SOFA error.
105 : return -1;
106 : }
107 :
108 : // Convert UTC to UT1
109 : rv = sofa::iauUtcut1( utc0, utc1, -.4, &ut10, &ut11 ); // The current DUT1 - how to keep updated?
110 : if( rv < 0 )
111 : {
112 : /// \todo handle SOFA error.
113 : return -1;
114 : }
115 :
116 : GMST = sofa::iauGmst06( ut10, ut11, tt0, tt1 ) / ( math::two_pi<realT>() ) * static_cast<realT>( 24 );
117 :
118 : return 0;
119 : }
120 :
121 : /// Returns Greenwich Mean Sidereal Time for a given UTC time.
122 : /** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
123 : * \param GMST [output] GMST in hours
124 : * \param Yr [input] integer year
125 : * \param Mo [input] integer month (1-12)
126 : * \param day [input] decimal day (1-31)
127 : * \retval 0 on sucess
128 : */
129 : template <typename realT>
130 : int getGMST( realT &GMST, ///< [out] GMST in hours
131 : int Yr, ///< [in] integer year
132 : int Mo, ///< [in] integer month (1-12)
133 : realT day ///< [in] decimal day (1-31)
134 : )
135 : {
136 : int Dy, hr, min;
137 : double sec;
138 :
139 : hrMinSecFromDay<realT>( Dy, hr, min, sec, day );
140 : return getGMST( GMST, Yr, Mo, Dy, hr, min, sec );
141 : }
142 :
143 : /// Returns Local Mean Sidereal Time for a given UTC time and longitude.
144 : /** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
145 : *
146 : * \retval 0 on success
147 : *
148 : */
149 : template <typename realT>
150 : int getLMST( realT &LMST, ///< [out] LMST in hours
151 : int Yr, ///< [in] integer year
152 : int Mo, ///< [in] integer month (1-12)
153 : int Dy, ///< [in] integer day (1-31)
154 : int hr, ///< [in] integer hour (0-23)
155 : int min, ///< [in] integer minute (0-59)
156 : realT sec, ///< [in] decimal second (0-59)
157 : realT lon ///< [in] longitude in degrees, E+, W-
158 : )
159 : {
160 : int rv;
161 : realT GMST;
162 :
163 : rv = getGMST( GMST, Yr, Mo, Dy, hr, min, sec );
164 :
165 : LMST = GMST + lon / static_cast<realT>( 15 );
166 :
167 : LMST = fmod( LMST, 24 );
168 : if( LMST < 0.0 )
169 : LMST += 24.;
170 :
171 : return rv;
172 : }
173 :
174 : /// Returns Local Mean Sidereal Time for a given UTC time.
175 : /** Uses the SOFA 2006 GMST function, and a (for now) hardcoded value of dUT1
176 : *
177 : * \tparam realT the floating point type for all calculations
178 : *
179 : * \retval 0 on sucess
180 : */
181 : template <typename realT>
182 : int getLMST( realT &LMST, ///< [out] LMST in hours
183 : int Yr, ///< [in] integer year
184 : int Mo, ///< [in] integer month (1-12)
185 : realT day, ///< [in] decimal day (1-31)
186 : realT lon ///< [in] longitude in degrees, E+, W-
187 : )
188 : {
189 : int Dy, hr, min;
190 : realT sec;
191 :
192 : hrMinSecFromDay<realT>( Dy, hr, min, sec, day );
193 :
194 : return getLMST( LMST, Yr, Mo, Dy, hr, min, sec, lon );
195 : }
196 :
197 : /// Calculates the Azimuth and Elevation for a given Hour-Angle, Declination, and Latitude.
198 : /** References: J. Meeus "Astronomical Algorithms", 1991; V. Pisacane "Fundamentals of Space Systems", 2nd ed., 2005.
199 : *
200 : * \tparam realT the real floating point type for calculations
201 : *
202 : * \returns 0 on success
203 : * \returns -1 on error
204 : *
205 : * \ingroup posastro
206 : */
207 : template <typename realT>
208 : int calcAzEl( realT &az, ///< [out] the calculated azimuth [degrees]
209 : realT &el, ///< [out] the calculate elevation [degrees]
210 : realT ha, ///< [in] the hour ange [degrees]
211 : realT dec, ///< [in] the declination [degrees]
212 : realT lat ///< [in] the latitude [degrees]
213 : )
214 : {
215 : realT ha_rad, dec_rad, lat_rad, az_rad, el_rad;
216 :
217 : ha_rad = ha * static_cast<realT>( 15 ) * math::pi<realT>() / static_cast<realT>( 180 );
218 : dec_rad = dec * math::pi<realT>() / static_cast<realT>( 180 );
219 : lat_rad = lat * math::pi<realT>() / static_cast<realT>( 180 );
220 :
221 : az_rad = atan2( sin( ha_rad ), cos( ha_rad ) * sin( lat_rad ) - tan( dec_rad ) * cos( lat_rad ) );
222 :
223 : el_rad = asin( sin( lat_rad ) * sin( dec_rad ) + cos( lat_rad ) * cos( dec_rad ) * cos( ha_rad ) );
224 :
225 : az = az_rad * static_cast<realT>( 180 ) / math::pi<realT>() + static_cast<realT>( 180 );
226 : el = el_rad * static_cast<realT>( 180 ) / math::pi<realT>();
227 :
228 : return 0;
229 : }
230 :
231 : /// Calculate the Parallactic Angle from the pre-calculated trig functions. Result in radians.
232 : /**
233 : * \returns the parallactic angle in radians.
234 : *
235 : * \test Scenario: calculating parallactic angles \ref tests_astrodynamics_parang "[test doc]"
236 : *
237 : * \ingroup posastro
238 : */
239 : template <typename realT>
240 2 : realT parAngRad( realT sinHA, ///< [in] the sine of the target hour angle
241 : realT cosHA, ///< [in] the cosine of the target hour angle
242 : realT sinDec, ///< [in] the sine of the target declination
243 : realT cosDec, ///< [in] the cosine of the target declination
244 : realT tanLat ///< [in] the tangent of the observer latitude
245 : )
246 : {
247 2 : return atan2( sinHA, cosDec * tanLat - sinDec * cosHA );
248 : }
249 :
250 : /// Calculate the Parallactic Angle from the pre-calculated trig functions. Result in degrees.
251 : /**
252 : * \returns the parallactic angle in degrees.
253 : *
254 : *
255 : *
256 : * \ingroup posastro
257 : */
258 : template <typename realT>
259 : realT parAngDeg( realT sinHA, ///< [in] the sine of the target hour angle
260 : realT cosHA, ///< [in] the cosine of the target hour angle
261 : realT sinDec, ///< [in] the sine of the target declination
262 : realT cosDec, ///< [in] the cosine of the target declination
263 : realT tanLat ///< [in] the tangent of the observer latitude
264 : )
265 : {
266 : return math::rtod( parAngRad( sinHA, cosHA, sinDec, cosDec, tanLat ) );
267 : }
268 :
269 : /// Calculate the Parallactic Angle, with angles in radians
270 : /**
271 : * \return the parallactic angle in radians.
272 : *
273 : * \test Scenario: calculating parallactic angles \ref tests_astrodynamics_parang "[test doc]"
274 : *
275 : * \ingroup posastro
276 : */
277 : template <typename realT>
278 2 : realT parAngRad( realT ha, ///< [in] the hour angle, in radians, negative to the east
279 : realT dec, ///< [in] the object declination, in radians.
280 : realT lat ///< [in] the observer latitude, in radians.
281 : )
282 : {
283 2 : return parAngRad( sin( ha ), cos( ha ), sin( dec ), cos( dec ), tan( lat ) );
284 : }
285 :
286 : /// Calculate the Parallactic Angle, with angles in degrees
287 : /**
288 : * \return the parallactic angle in degrees.
289 : *
290 : * \test Scenario: calculating parallactic angles \ref tests_astrodynamics_parang "[test doc]"
291 : *
292 : * \ingroup posastro
293 : *
294 : * \note 2020-Jan-19: re-reordered arguments and added newflag to prevent compilation of current programs so the new
295 : * order is implemented.
296 : */
297 : template <typename realT>
298 2 : realT parAngDeg( realT ha, ///< [in] the hour angle, in degrees, negative to the east
299 : realT dec, ///< [in] the object declination, in degrees.
300 : realT lat, ///< [in] the observer latitude, in degrees.
301 : bool newflag ///< [in] [temporary] to force adaptation of new argument order. Unused, so can be true
302 : ///< or false. Added 2020-Jan-19.
303 : )
304 : {
305 : static_cast<void>( newflag ); // be unused
306 :
307 2 : return math::rtod( parAngRad( math::dtor( ha ), math::dtor( dec ), math::dtor( lat ) ) );
308 : }
309 :
310 : } // namespace astro
311 : } // namespace mx
312 :
313 : #endif // astroDynamics_hpp
|
