1 /* 2 * This library is free software; you can redistribute it and/or 3 * modify it under the terms of the GNU Lesser General Public 4 * License as published by the Free Software Foundation; either 5 * version 2 of the License, or (at your option) any later version. 6 * 7 * This library is distributed in the hope that it will be useful, 8 * but WITHOUT ANY WARRANTY; without even the implied warranty of 9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 10 * Lesser General Public License for more details. 11 * 12 * You should have received a copy of the GNU General Public License 13 * along with this program; if not, write to the Free Software 14 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 15 * 16 * Copyright (C) 2000 - 2005 Liam Girdwood <lgirdwood@gmail.com> 17 */ 18 19 module nova.nutation; 20 21 import std.math; 22 import nova.nutation; 23 import nova.dynamical_time; 24 import nova.utility; 25 import nova.ln_types; 26 27 enum TERMS = 63; 28 enum LN_NUTATION_EPOCH_THRESHOLD = 0.1; 29 30 struct nutation_arguments { 31 double D; 32 double M; 33 double MM; 34 double F; 35 double O; 36 } 37 38 struct nutation_coefficients { 39 double longitude1; 40 double longitude2; 41 double obliquity1; 42 double obliquity2; 43 } 44 45 /* arguments and coefficients taken from table 21A on page 133 */ 46 47 const static nutation_arguments[TERMS] arguments = [ 48 {0.0, 0.0, 0.0, 0.0, 1.0}, 49 {-2.0, 0.0, 0.0, 2.0, 2.0}, 50 {0.0, 0.0, 0.0, 2.0, 2.0}, 51 {0.0, 0.0, 0.0, 0.0, 2.0}, 52 {0.0, 1.0, 0.0, 0.0, 0.0}, 53 {0.0, 0.0, 1.0, 0.0, 0.0}, 54 {-2.0, 1.0, 0.0, 2.0, 2.0}, 55 {0.0, 0.0, 0.0, 2.0, 1.0}, 56 {0.0, 0.0, 1.0, 2.0, 2.0}, 57 {-2.0, -1.0, 0.0, 2.0, 2.0}, 58 {-2.0, 0.0, 1.0, 0.0, 0.0}, 59 {-2.0, 0.0, 0.0, 2.0, 1.0}, 60 {0.0, 0.0, -1.0, 2.0, 2.0}, 61 {2.0, 0.0, 0.0, 0.0, 0.0}, 62 {0.0, 0.0, 1.0, 0.0, 1.0}, 63 {2.0, 0.0, -1.0, 2.0, 2.0}, 64 {0.0, 0.0, -1.0, 0.0, 1.0}, 65 {0.0, 0.0, 1.0, 2.0, 1.0}, 66 {-2.0, 0.0, 2.0, 0.0, 0.0}, 67 {0.0, 0.0, -2.0, 2.0, 1.0}, 68 {2.0, 0.0, 0.0, 2.0, 2.0}, 69 {0.0, 0.0, 2.0, 2.0, 2.0}, 70 {0.0, 0.0, 2.0, 0.0, 0.0}, 71 {-2.0, 0.0, 1.0, 2.0, 2.0}, 72 {0.0, 0.0, 0.0, 2.0, 0.0}, 73 {-2.0, 0.0, 0.0, 2.0, 0.0}, 74 {0.0, 0.0, -1.0, 2.0, 1.0}, 75 {0.0, 2.0, 0.0, 0.0, 0.0}, 76 {2.0, 0.0, -1.0, 0.0, 1.0}, 77 {-2.0, 2.0, 0.0, 2.0, 2.0}, 78 {0.0, 1.0, 0.0, 0.0, 1.0}, 79 {-2.0, 0.0, 1.0, 0.0, 1.0}, 80 {0.0, -1.0, 0.0, 0.0, 1.0}, 81 {0.0, 0.0, 2.0, -2.0, 0.0}, 82 {2.0, 0.0, -1.0, 2.0, 1.0}, 83 {2.0, 0.0, 1.0, 2.0, 2.0}, 84 {0.0, 1.0, 0.0, 2.0, 2.0}, 85 {-2.0, 1.0, 1.0, 0.0, 0.0}, 86 {0.0, -1.0, 0.0, 2.0, 2.0}, 87 {2.0, 0.0, 0.0, 2.0, 1.0}, 88 {2.0, 0.0, 1.0, 0.0, 0.0}, 89 {-2.0, 0.0, 2.0, 2.0, 2.0}, 90 {-2.0, 0.0, 1.0, 2.0, 1.0}, 91 {2.0, 0.0, -2.0, 0.0, 1.0}, 92 {2.0, 0.0, 0.0, 0.0, 1.0}, 93 {0.0, -1.0, 1.0, 0.0, 0.0}, 94 {-2.0, -1.0, 0.0, 2.0, 1.0}, 95 {-2.0, 0.0, 0.0, 0.0, 1.0}, 96 {0.0, 0.0, 2.0, 2.0, 1.0}, 97 {-2.0, 0.0, 2.0, 0.0, 1.0}, 98 {-2.0, 1.0, 0.0, 2.0, 1.0}, 99 {0.0, 0.0, 1.0, -2.0, 0.0}, 100 {-1.0, 0.0, 1.0, 0.0, 0.0}, 101 {-2.0, 1.0, 0.0, 0.0, 0.0}, 102 {1.0, 0.0, 0.0, 0.0, 0.0}, 103 {0.0, 0.0, 1.0, 2.0, 0.0}, 104 {0.0, 0.0, -2.0, 2.0, 2.0}, 105 {-1.0, -1.0, 1.0, 0.0, 0.0}, 106 {0.0, 1.0, 1.0, 0.0, 0.0}, 107 {0.0, -1.0, 1.0, 2.0, 2.0}, 108 {2.0, -1.0, -1.0, 2.0, 2.0}, 109 {0.0, 0.0, 3.0, 2.0, 2.0}, 110 {2.0, -1.0, 0.0, 2.0, 2.0} 111 ]; 112 113 const static nutation_coefficients[TERMS] coefficients = [ 114 {-171996.0, -174.2, 92025.0,8.9}, 115 {-13187.0, -1.6, 5736.0, -3.1}, 116 {-2274.0, -0.2, 977.0, -0.5}, 117 {2062.0, 0.2, -895.0, 0.5}, 118 {1426.0, -3.4, 54.0, -0.1}, 119 {712.0, 0.1, -7.0, 0.0}, 120 {-517.0, 1.2, 224.0, -0.6}, 121 {-386.0, -0.4, 200.0, 0.0}, 122 {-301.0, 0.0, 129.0, -0.1}, 123 {217.0, -0.5, -95.0, 0.3}, 124 {-158.0, 0.0, 0.0, 0.0}, 125 {129.0, 0.1, -70.0, 0.0}, 126 {123.0, 0.0, -53.0, 0.0}, 127 {63.0, 0.0, 0.0, 0.0}, 128 {63.0, 0.1, -33.0, 0.0}, 129 {-59.0, 0.0, 26.0, 0.0}, 130 {-58.0, -0.1, 32.0, 0.0}, 131 {-51.0, 0.0, 27.0, 0.0}, 132 {48.0, 0.0, 0.0, 0.0}, 133 {46.0, 0.0, -24.0, 0.0}, 134 {-38.0, 0.0, 16.0, 0.0}, 135 {-31.0, 0.0, 13.0, 0.0}, 136 {29.0, 0.0, 0.0, 0.0}, 137 {29.0, 0.0, -12.0, 0.0}, 138 {26.0, 0.0, 0.0, 0.0}, 139 {-22.0, 0.0, 0.0, 0.0}, 140 {21.0, 0.0, -10.0, 0.0}, 141 {17.0, -0.1, 0.0, 0.0}, 142 {16.0, 0.0, -8.0, 0.0}, 143 {-16.0, 0.1, 7.0, 0.0}, 144 {-15.0, 0.0, 9.0, 0.0}, 145 {-13.0, 0.0, 7.0, 0.0}, 146 {-12.0, 0.0, 6.0, 0.0}, 147 {11.0, 0.0, 0.0, 0.0}, 148 {-10.0, 0.0, 5.0, 0.0}, 149 {-8.0, 0.0, 3.0, 0.0}, 150 {7.0, 0.0, -3.0, 0.0}, 151 {-7.0, 0.0, 0.0, 0.0}, 152 {-7.0, 0.0, 3.0, 0.0}, 153 {-7.0, 0.0, 3.0, 0.0}, 154 {6.0, 0.0, 0.0, 0.0}, 155 {6.0, 0.0, -3.0, 0.0}, 156 {6.0, 0.0, -3.0, 0.0}, 157 {-6.0, 0.0, 3.0, 0.0}, 158 {-6.0, 0.0, 3.0, 0.0}, 159 {5.0, 0.0, 0.0, 0.0}, 160 {-5.0, 0.0, 3.0, 0.0}, 161 {-5.0, 0.0, 3.0, 0.0}, 162 {-5.0, 0.0, 3.0, 0.0}, 163 {4.0, 0.0, 0.0, 0.0}, 164 {4.0, 0.0, 0.0, 0.0}, 165 {4.0, 0.0, 0.0, 0.0}, 166 {-4.0, 0.0, 0.0, 0.0}, 167 {-4.0, 0.0, 0.0, 0.0}, 168 {-4.0, 0.0, 0.0, 0.0}, 169 {3.0, 0.0, 0.0, 0.0}, 170 {-3.0, 0.0, 0.0, 0.0}, 171 {-3.0, 0.0, 0.0, 0.0}, 172 {-3.0, 0.0, 0.0, 0.0}, 173 {-3.0, 0.0, 0.0, 0.0}, 174 {-3.0, 0.0, 0.0, 0.0}, 175 {-3.0, 0.0, 0.0, 0.0}, 176 {-3.0, 0.0, 0.0, 0.0} 177 ]; 178 179 /* cache values */ 180 static real c_JD = 0.0, c_longitude = 0.0, c_obliquity = 0.0, c_ecliptic = 0.0; 181 182 extern (C) { 183 184 /*! \fn void ln_get_nutation(double JD, struct ln_nutation *nutation) 185 * \param JD Julian Day. 186 * \param nutation Pointer to store nutation 187 * 188 * Calculate nutation of longitude and obliquity in degrees from Julian Ephemeris Day 189 */ 190 /* Chapter 21 pg 131-134 Using Table 21A 191 */ 192 /* TODO: add argument to specify this */ 193 /* TODO: use JD or JDE. confirm */ 194 @nogc void ln_get_nutation(double JD, ref ln_nutation nutation) nothrow 195 { 196 real D, M, MM, F, O, T, T2, T3, JDE; 197 real coeff_sine, coeff_cos; 198 real argument; 199 int i; 200 201 /* should we bother recalculating nutation */ 202 if (fabs(JD - c_JD) > LN_NUTATION_EPOCH_THRESHOLD) { 203 /* set the new epoch */ 204 c_JD = JD; 205 c_longitude = 0; 206 c_obliquity = 0; 207 208 /* get julian ephemeris day */ 209 JDE = cast(real)ln_get_jde(JD); 210 211 /* calc T */ 212 T = (JDE - 2451545.0) / 36525.0; 213 T2 = T * T; 214 T3 = T2 * T; 215 216 /* calculate D,M,M',F and Omega */ 217 D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474.0; 218 M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300000.0; 219 MM = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250.0; 220 F = 93.2719100 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270.0; 221 O = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0; 222 223 /* convert to radians */ 224 D = ln_deg_to_rad(D); 225 M = ln_deg_to_rad(M); 226 MM = ln_deg_to_rad(MM); 227 F = ln_deg_to_rad(F); 228 O = ln_deg_to_rad(O); 229 230 /* calc sum of terms in table 21A */ 231 for (i = 0; i < TERMS; i++) { 232 /* calc coefficients of sine and cosine */ 233 coeff_sine = (coefficients[i].longitude1 + 234 (coefficients[i].longitude2 * T)); 235 coeff_cos = (coefficients[i].obliquity1 + 236 (coefficients[i].obliquity2 * T)); 237 238 argument = arguments[i].D * D 239 + arguments[i].M * M 240 + arguments[i].MM * MM 241 + arguments[i].F * F 242 + arguments[i].O * O; 243 244 c_longitude += coeff_sine * sin(argument); 245 c_obliquity += coeff_cos * cos(argument); 246 } 247 248 /* change to arcsecs */ 249 c_longitude /= 10000.0; 250 c_obliquity /= 10000.0; 251 252 /* change to degrees */ 253 c_longitude /= (60.0 * 60.0); 254 c_obliquity /= (60.0 * 60.0); 255 256 /* calculate mean ecliptic - Meeus 2nd edition, eq. 22.2 */ 257 c_ecliptic = 23.0 + 26.0 / 60.0 + 21.448 / 3600.0 258 - 46.8150 / 3600.0 * T 259 - 0.00059 / 3600.0 * T2 260 + 0.001813 / 3600.0 * T3; 261 262 /* c_ecliptic += c_obliquity; * Uncomment this if function should 263 return true obliquity rather than 264 mean obliquity */ 265 } 266 267 /* return results */ 268 nutation.longitude = c_longitude; 269 nutation.obliquity = c_obliquity; 270 nutation.ecliptic = c_ecliptic; 271 } 272 273 /*! \fn void ln_get_equ_nut(struct ln_equ_posn *mean_position, double JD, struct ln_equ_posn *position) 274 * \param mean_position Mean position of object 275 * \param JD Julian Day. 276 * \param position Pointer to store new object position. 277 * 278 * Calculate a stars equatorial coordinates from it's mean equatorial coordinates 279 * with the effects of nutation for a given Julian Day. 280 */ 281 /* Equ 22.1 282 */ 283 @nogc void ln_get_equ_nut(const ref ln_equ_posn mean_position, double JD, 284 ref ln_equ_posn position) nothrow 285 { 286 ln_nutation nut; 287 ln_get_nutation (JD, nut); 288 289 real mean_ra, mean_dec, delta_ra, delta_dec; 290 291 mean_ra = ln_deg_to_rad(mean_position.ra); 292 mean_dec = ln_deg_to_rad(mean_position.dec); 293 294 // Equ 22.1 295 296 real nut_ecliptic = ln_deg_to_rad(nut.ecliptic + nut.obliquity); 297 real sin_ecliptic = sin(nut_ecliptic); 298 299 real sin_ra = sin(mean_ra); 300 real cos_ra = cos(mean_ra); 301 302 real tan_dec = tan(mean_dec); 303 304 delta_ra = (cos (nut_ecliptic) + sin_ecliptic * sin_ra * tan_dec) * nut.longitude - cos_ra * tan_dec * nut.obliquity; 305 delta_dec = (sin_ecliptic * cos_ra) * nut.longitude + sin_ra * nut.obliquity; 306 307 position.ra = mean_position.ra + delta_ra; 308 position.dec = mean_position.dec + delta_dec; 309 } 310 311 }