1 /*****************************************************************************\
2 * AstroOps.cpp
3 *
4 * AstroOps is a 'catch-all' class that performs misc. handy calculations
5 *
6 * author: mark huss (mark@mhuss.com)
7 * Based on Bill Gray's open-source code at projectpluto.com
8 *
9 \*****************************************************************************/
10
11 #include "AstroOps.h"
12
13 #include <math.h>
14
15 // In-class constants cannot be initialized in the class declaration. Sigh.
16 //
17 const double Astro::PI = PI_DEF;
18 const double Astro::TWO_PI = PI_DEF + PI_DEF;
19 const double Astro::PI_OVER_TWO = PI_DEF / 2.;
20 const double Astro::TWO_OVER_PI = 2. / PI_DEF;
21
22 const double Astro::DEG_PER_CIRCLE = 360.;
23 const double Astro::MINUTES_PER_DEGREE = 60.;
24 const double Astro::SECONDS_PER_DEGREE = 3600.;
25
26 const double Astro::TO_RADS = PI_DEF / 180.;
27 const double Astro::TO_CENTURIES = 36525.;
28
29 const double Astro::J2000 = 2451545.0;
30
31 const double Astro::HOURS_PER_DAY = 24.;
32 const int Astro::IHOURS_PER_DAY = 24;
33
34 /*----------------------------------------------------------------------------
35 * The obliquity formula (and all the magic numbers below) come from Meeus,
36 * Astro Algorithms.
37 *
38 * Input t is time in julian centuries from 2000.
39 * Valid range is the years -8000 to +12000 (t = -100 to 100).
40 *
41 * return value is mean obliquity (epsilon sub 0) in radians.
42 */
43
meanObliquity(double t)44 double AstroOps::meanObliquity( double t )
45 {
46 double u, u0;
47 static double t0 = 30000.;
48 static double rval = 0.;
49 static const double rvalStart = 23. * Astro::SECONDS_PER_DEGREE +
50 26. * Astro::MINUTES_PER_DEGREE +
51 21.448;
52 static const int OBLIQ_COEFFS = 10;
53 static const double coeffs[ OBLIQ_COEFFS ] = {
54 -468093., -155., 199925., -5138., -24967.,
55 -3905., 712., 2787., 579., 245.
56 };
57
58 if( t0 != t ) {
59
60 t0 = t;
61 u = u0 = t / 100.; // u is in julian 10000's of years
62 rval = rvalStart;
63
64 for( int i=0; i<OBLIQ_COEFFS; i++ )
65 {
66 rval += u * coeffs[i] / 100.;
67 u *= u0;
68 }
69 // convert from seconds to radians
70 rval = Astro::toRadians( rval / Astro::SECONDS_PER_DEGREE );
71 }
72 return rval;
73 }
74
75 /*----------------------------------------------------------------------------
76 * The nutation formula (and all the magic numbers below) come from
77 * p 132-5, Meeus, Astro Algorithms
78 *
79 * input is time in julian centuries from 2000.
80 *
81 * *pDPhi is nutation (delta phi) in arcseconds
82 * *pDEpsilon is nutation (delta epsilon) in arcseconds
83 *
84 * Either pointer can be NULL, in which case that value is
85 * not computed. (we added this because sometimes, you want
86 * only pDPhi or pDEpsilon; in such cases, computing _both_
87 * is a waste of perfectly good CPU time)
88 */
89
nutation(double t,double * pDPhi,double * pDEpsilon)90 void AstroOps::nutation( double t, double* pDPhi, double* pDEpsilon )
91 {
92 static const double linearPart[5] = {
93 445267.111480, 35999.050340, 477198.867398,
94 483202.017538, -1934.136261 };
95
96 static const double coefficients[5][3] = {
97 { 29785036., -19142., 189474. },
98 { 35752772., - 1603., -300000. },
99 { 13496298., 86972., 56250. },
100 { 9327191., -36825., 327270. },
101 { 12504452., 20708., 450000. }
102 };
103
104 static const int N_NUTATION_COEFFS = 62;
105
106 static const int args[N_NUTATION_COEFFS][3] = {
107 { 324,-13187,5736 },
108 { 1574, -2274, 977 },
109 { 1564, 2062,-895 },
110 { 1687, 1426, 54 },
111 { 1587, 712, -7 },
112 { 449, -517, 224 },
113 { 1573, -386, 200 },
114 { 1599, -301, 129 },
115 { 199, 217, -95 },
116 { 337, -158, 0 },
117 { 323, 129, -70 },
118 { 1549, 123, -53 },
119 { 2812, 63, 0 },
120 { 1588, 63, -33 },
121 { 2799, -59, 26 },
122 { 1538, -58, 32 },
123 { 1598, -51, 27 },
124 { 362, 48, 0 },
125 { 1523, 46, -24 },
126 { 2824, -38, 16 },
127 { 1624, -31, 13 },
128 { 1612, 29, 0 },
129 { 349, 29, -12 },
130 { 1572, 26, 0 },
131 { 322, -22, 0 },
132 { 1548, 21, -10 },
133 { 1812, 17, 0 },
134 { 2788, 16, -8 },
135 { 574, -16, 7 },
136 { 1688, -15, 9 },
137 { 338, -13, 7 },
138 { 1438, -12, 6 },
139 { 1602, 11, 0 },
140 { 2798, -10, 5 },
141 { 2849, -8, 3 },
142 { 1699, 7, -3 },
143 { 462, -7, 0 },
144 { 1449, -7, 3 },
145 { 2823, -7, 3 },
146 { 2837, 6, 0 },
147 { 374, 6, -3 },
148 { 348, 6, -3 },
149 { 2763, -6, 3 },
150 { 2813, -6, 3 },
151 { 1462, 5, 0 },
152 { 198, -5, 3 },
153 { 313, -5, 3 },
154 { 1623, -5, 3 },
155 { 1524, -3, 0 },
156 { 363, 4, 0 },
157 { 448, 4, 0 },
158 { 1474, -3, 0 },
159 { 2674, -3, 0 },
160 { 1649, -3, 0 },
161 { 2699, -3, 0 },
162 { 837, -3, 0 },
163 { 962, -4, 0 },
164 { 437, -4, 0 },
165 { 1577, 4, 0 },
166 { 2187, -4, 0 },
167 { 1712, -3, 0 },
168 { 1597, 3, 0 }
169 };
170
171 static const int timeDependent[16 + 9] = {
172 -16, -2, 2, -34, 1, 12, -4, 0, -5, 0, 1, 0, 0, 1, 0, -1,
173 -31, -5, 5, -1, 0, -6, 0, -1, 3 };
174
175 double terms[5];
176
177 double t2 = t * t;
178 double t3 = t2 * t;
179
180 for( int i=0; i<5; i++ ) {
181 terms[i] = linearPart[i] * t + coefficients[i][0] / 100000.;
182 terms[i] += t2 * coefficients[i][1] * 1.e-7;
183 terms[i] += t3 / coefficients[i][2];
184 terms[i] *= Astro::TO_RADS;
185 }
186
187 // The largest term in pDPhi won't fit into a short int.
188 // Everything else does, luckily.
189 //
190 if( pDPhi )
191 *pDPhi = (-171996. - 174.2 * t) * sin( terms[4] );
192
193 if( pDEpsilon )
194 *pDEpsilon = (92025. + 8.9 * t) * cos( terms[4] );
195
196 for( int i=0; i<N_NUTATION_COEFFS; i++ ) {
197 double totalArg = 0.;
198 int mult = args[i][0];
199
200 for( int j=4; j>=0; j--) {
201 if( mult % 5 != 2)
202 totalArg += (double)( mult % 5 - 2) * terms[j];
203 mult /= 5;
204 }
205
206 double coeff = (double)(args[i][1]);
207 if( i < 16 && 0 != timeDependent[i] )
208 coeff += (double)(timeDependent[i]) * t / 10.;
209 else if( 26 == i || 28 == i )
210 coeff += (double)(27 - i) * t / 10.;
211
212 if( pDPhi )
213 *pDPhi += coeff * sin(totalArg);
214
215 if( 0 != args[i][2] ) {
216 coeff = (double)(args[i][2]);
217
218 if( i < 9 && 0 != timeDependent[i + 16] )
219 coeff += (double)(timeDependent[i + 16]) * t / 10.;
220
221 if( pDEpsilon)
222 *pDEpsilon += coeff * cos( totalArg );
223 }
224 }
225 if( pDPhi )
226 *pDPhi *= .0001;
227
228 if( pDEpsilon )
229 *pDEpsilon *= .0001;
230 }
231
232 /*----------------------------------------------------------------------------
233 * See p 84, in Meeus:
234 * The following returns apparent Greenwich sidereal time for the given jd
235 */
greenwichSiderealTime(double jd)236 double AstroOps::greenwichSiderealTime( double jd )
237 {
238 jd -= Astro::J2000; /* set relative to 2000.0 */
239 double jdm = jd / Astro::TO_CENTURIES; /* convert jd to julian centuries */
240 double intPart = floor( jd );
241 jd -= intPart;
242 double rval = 280.46061837 +
243 360.98564736629 * jd +
244 .98564736629 * intPart +
245 jdm * jdm * ( 3.87933e-4 - jdm / 38710000. );
246
247 return Astro::toRadians( rval );
248 }
249
250 //----------------------------------------------------------------------------
251 // reduce an angle in degrees to (0 <= d < 360)
252 //
normalizeDegrees(double d)253 double AstroOps::normalizeDegrees(double d)
254 {
255 d = fmod( d, Astro::DEG_PER_CIRCLE);
256 if( d < 0.)
257 d += Astro::DEG_PER_CIRCLE;
258
259 return d;
260 }
261
262 //----------------------------------------------------------------------------
263 // reduce an angle in radians to (0 <= r < 2PI)
264 //
normalizeRadians(double r)265 double AstroOps::normalizeRadians(double r)
266 {
267 r = fmod( r, Astro::TWO_PI );
268 if( r < 0. )
269 r += Astro::TWO_PI;
270
271 return r;
272 }
273
274 //----------------------------------------------------------------------------
275