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 17 * Copyright (C) 2015 Jeroen Vreeken (jeroen@vreeken.net) 18 */ 19 20 module nova.aberration; 21 22 import std.math; 23 import nova.aberration; 24 import nova.solar; 25 import nova.utility; 26 import nova.ln_types; 27 28 enum TERMS = 36; 29 30 /* data structures to hold arguments and coefficients of Ron-Vondrak theory */ 31 struct arg 32 { 33 double a_L2; 34 double a_L3; 35 double a_L4; 36 double a_L5; 37 double a_L6; 38 double a_L7; 39 double a_L8; 40 double a_LL; 41 double a_D; 42 double a_MM; 43 double a_F; 44 } 45 46 struct XYZ 47 { 48 double sin1; 49 double sin2; 50 double cos1; 51 double cos2; 52 } 53 54 const static arg[TERMS] arguments = [ 55 /* L2 3 4 5 6 7 8 LL D MM F */ 56 {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 57 {0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 58 {0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, 59 {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, 60 {0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 61 {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, 62 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, 63 {0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0}, 64 {0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0}, 65 {0, 2, 0, -1, 0, 0, 0, 0, 0, 0, 0}, 66 {0, 3, -8, 3, 0, 0, 0, 0, 0, 0, 0}, 67 {0, 5, -8, 3, 0, 0, 0, 0, 0, 0, 0}, 68 {2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 69 {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 70 {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, 71 {0, 1, 0, -2, 0, 0, 0, 0, 0, 0, 0}, 72 {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, 73 {0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0}, 74 {2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 75 {0, 1, 0, -1, 0, 0, 0, 0, 0, 0, 0}, 76 {0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 77 {0, 3, 0, -2, 0, 0, 0, 0, 0, 0, 0}, 78 {1, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 79 {2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 80 {0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0}, 81 {2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 82 {0, 3, -2, 0, 0, 0, 0, 0, 0, 0, 0}, 83 {0, 0, 0, 0, 0, 0, 0, 1, 2, -1, 0}, 84 {8, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 85 {8, 14, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 86 {0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0}, 87 {3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 88 {0, 2, 0, -2, 0, 0, 0, 0, 0, 0, 0}, 89 {3, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 90 {0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0}, 91 {0, 0, 0, 0, 0, 0, 0, 1, -2, 0, 0} 92 ]; 93 94 const static XYZ[TERMS] x_coefficients = [ 95 {-1719914, -2, -25, 0}, 96 {6434, 141, 28007, -107}, 97 {715, 0, 0, 0}, 98 {715, 0, 0, 0}, 99 {486, -5, -236, -4}, 100 {159, 0, 0, 0}, 101 {0, 0, 0, 0}, 102 {39, 0, 0, 0}, 103 {33, 0, -10, 0}, 104 {31, 0, 1, 0}, 105 {8, 0, -28, 0}, 106 {8, 0, -28, 0}, 107 {21, 0, 0, 0}, 108 {-19, 0, 0, 0}, 109 {17, 0, 0, 0}, 110 {16, 0, 0, 0}, 111 {16, 0, 0, 0}, 112 {11, 0, -1, 0}, 113 {0, 0, -11, 0}, 114 {-11, 0, -2, 0}, 115 {-7, 0, -8, 0}, 116 {-10, 0, 0, 0}, 117 {-9, 0, 0, 0}, 118 {-9, 0, 0, 0}, 119 {0, 0, -9, 0}, 120 {0, 0, -9, 0}, 121 {8, 0, 0, 0}, 122 {8, 0, 0, 0}, 123 {-4, 0, -7, 0}, 124 {-4, 0, -7, 0}, 125 {-6, 0, -5, 0}, 126 {-1, 0, -1, 0}, 127 {4, 0, -6, 0}, 128 {0, 0, -7, 0}, 129 {5, 0, -5, 0}, 130 {5, 0, 0, 0} 131 ]; 132 133 const static XYZ[TERMS] y_coefficients = [ 134 {25, -13, 1578089, 156}, 135 {25697, -95, -5904, -130}, 136 {6, 0, -657, 0}, 137 {0, 0, -656, 0}, 138 {-216, -4, -446, 5}, 139 {2, 0, -147, 0}, 140 {0, 0, 26, 0}, 141 {0, 0, -36, 0}, 142 {-9, 0, -30, 0}, 143 {1, 0, -28, 0}, 144 {25, 0, 8, 0}, 145 {-25, 0, -8, 0}, 146 {0, 0, -19, 0}, 147 {0, 0, 17, 0}, 148 {0, 0, -16, 0}, 149 {0, 0, 15, 0}, 150 {1, 0, -15, 0}, 151 {-1, 0, -10, 0}, 152 {-10, 0, 0, 0}, 153 {-2, 0, 9, 0}, 154 {-8, 0, 6, 0}, 155 {0, 0, 9, 0}, 156 {0, 0, -9, 0}, 157 {0, 0, -8, 0}, 158 {-8, 0, 0, 0}, 159 {8, 0, 0, 0}, 160 {0, 0, -8, 0}, 161 {0, 0, -7, 0}, 162 {-6, 0, -4, 0}, 163 {6, 0, -4, 0}, 164 {-4, 0, 5, 0}, 165 {-2, 0, -7, 0}, 166 {-5, 0, -4, 0}, 167 {-6, 0, 0, 0}, 168 {-4, 0, -5, 0}, 169 {0, 0, -5, 0} 170 ]; 171 172 const static XYZ[TERMS] z_coefficients = [ 173 {10, 32, 684185, -358}, 174 {11141, -48, -2559, -55}, 175 {-15, 0, -282, 0}, 176 {0, 0, -285, 0}, 177 {-94, 0, -193, 0}, 178 {-6, 0, -61, 0}, 179 {0, 0, 59, 0}, 180 {0, 0, 16, 0}, 181 {-5, 0, -13, 0}, 182 {0, 0, -12, 0}, 183 {11, 0, 3, 0}, 184 {-11, 0, -3, 0}, 185 {0, 0, -8, 0}, 186 {0, 0, 8, 0}, 187 {0, 0, -7, 0}, 188 {1, 0, 7, 0}, 189 {-3, 0, -6, 0}, 190 {-1, 0, 5, 0}, 191 {-4, 0, 0, 0}, 192 {-1, 0, 4, 0}, 193 {-3, 0, 3, 0}, 194 {0, 0, 4, 0}, 195 {0, 0, -4, 0}, 196 {0, 0, -4, 0}, 197 {-3, 0, 0, 0}, 198 {3, 0, 0, 0}, 199 {0, 0, -3, 0}, 200 {0, 0, -3, 0}, 201 {-3, 0, 2, 0}, 202 {3, 0, -2, 0}, 203 {-2, 0, 2, 0}, 204 {1, 0, -4, 0}, 205 {-2, 0, -2, 0}, 206 {-3, 0, 0, 0}, 207 {-2, 0, -2, 0}, 208 {0, 0, -2, 0} 209 ]; 210 211 extern (C) { 212 213 /*! \fn void ln_get_equ_aber(struct ln_equ_posn *mean_position, double JD, struct ln_equ_posn *position) 214 * \param mean_position Mean position of object 215 * \param JD Julian Day 216 * \param position Pointer to store new object position. 217 * 218 * Calculate a stars equatorial coordinates from it's mean equatorial coordinates 219 * with the effects of aberration for a given Julian Day. 220 */ 221 /* Equ 22.3, 22.4 222 */ 223 @nogc void ln_get_equ_aber(const ref ln_equ_posn mean_position, double JD, 224 ref ln_equ_posn position) nothrow 225 { 226 real mean_ra, mean_dec, delta_ra, delta_dec; 227 real L2, L3, L4, L5, L6, L7, L8, LL, D, MM , F, T, X, Y, Z, A; 228 real c; 229 int i; 230 231 /* speed of light in 10-8 au per day */ 232 c = 17314463350.0; 233 234 /* calc T */ 235 T = (JD - 2451545.0) / 36525.0; 236 237 /* calc planetary perturbutions */ 238 L2 = 3.1761467 + 1021.3285546 * T; 239 L3 = 1.7534703 + 628.3075849 * T; 240 L4 = 6.2034809 + 334.0612431 * T; 241 L5 = 0.5995464 + 52.9690965 * T; 242 L6 = 0.8740168 + 21.329909095 * T; 243 L7 = 5.4812939 + 7.4781599 * T; 244 L8 = 5.3118863 + 3.8133036 * T; 245 LL = 3.8103444 + 8399.6847337 * T; 246 D = 5.1984667 + 7771.3771486 * T; 247 MM = 2.3555559 + 8328.6914289 * T; 248 F = 1.6279052 + 8433.4661601 * T; 249 250 X = 0; 251 Y = 0; 252 Z = 0; 253 254 /* sum the terms */ 255 for (i = 0; i < TERMS; i++) { 256 A = arguments[i].a_L2 * L2 + arguments[i].a_L3 * L3 + 257 arguments[i].a_L4 * L4 + arguments[i].a_L5 * L5 + 258 arguments[i].a_L6 * L6 + arguments[i].a_L7 * L7 + 259 arguments[i].a_L8 * L8 + arguments[i].a_LL * LL + 260 arguments[i].a_D * D + arguments[i].a_MM * MM + 261 arguments[i].a_F * F; 262 263 X += (x_coefficients[i].sin1 + x_coefficients[i].sin2 * T) * 264 sin(A) + (x_coefficients[i].cos1 + x_coefficients[i].cos2 * T) * 265 cos(A); 266 Y += (y_coefficients[i].sin1 + y_coefficients[i].sin2 * T) * 267 sin(A) + (y_coefficients[i].cos1 + y_coefficients[i].cos2 * T) * 268 cos(A); 269 Z += (z_coefficients[i].sin1 + z_coefficients[i].sin2 * T) * 270 sin(A) + (z_coefficients[i].cos1 + z_coefficients[i].cos2 * T) * 271 cos(A); 272 } 273 274 /* Equ 22.4 */ 275 mean_ra = ln_deg_to_rad(mean_position.ra); 276 mean_dec = ln_deg_to_rad(mean_position.dec); 277 278 if (mean_dec < PI * 0.4999 ) { 279 delta_ra = (Y * cos(mean_ra) - X * sin(mean_ra)) / cos(mean_dec); 280 delta_ra /= c; 281 delta_dec = (X * cos(mean_ra) + Y * sin(mean_ra)) * sin(mean_dec) - Z * cos(mean_dec); 282 delta_dec /= -c; 283 284 position.ra = ln_rad_to_deg(mean_ra + delta_ra); 285 position.dec = ln_rad_to_deg(mean_dec + delta_dec); 286 } else { 287 /* cos(mean_dec) gets to small when approaching (or at) 90.0 degrees 288 Use an alternative method: 289 - First transform mean position to x,y 290 - Apply offset 291 - Transform back to ra/dec 292 */ 293 real px, py, ra, dec; 294 double cos_dec; 295 296 X /= c; 297 Y /= c; 298 Z /= c; 299 cos_dec = cos(mean_dec); 300 301 px = cos_dec * cos(mean_ra); 302 py = cos_dec * sin(mean_ra); 303 304 px += X; 305 py += Y; 306 307 ra = atan2(py, px); 308 dec = acos(sqrt(px * px + py * py)); 309 310 dec += cos_dec * Z; 311 312 position.ra = ln_rad_to_deg(ra); 313 position.dec = ln_rad_to_deg(dec); 314 } 315 } 316 317 /*! \fn void ln_get_ecl_aber(struct ln_lnlat_posn *mean_position, double JD, struct ln_lnlat_posn *position) 318 * \param mean_position Mean position of object 319 * \param JD Julian Day 320 * \param position Pointer to store new object position. 321 * 322 * Calculate a stars ecliptical coordinates from it's mean ecliptical coordinates 323 * with the effects of aberration for a given Julian Day. 324 */ 325 /* Equ 22.2 pg 139 326 */ 327 @nogc void ln_get_ecl_aber(const ref ln_lnlat_posn mean_position, double JD, 328 ref ln_lnlat_posn position) nothrow 329 { 330 double delta_lng, delta_lat, mean_lng, mean_lat, e, t, k; 331 double true_longitude, T, T2; 332 ln_helio_posn sol_position; 333 334 /* constant of aberration */ 335 k = ln_deg_to_rad(20.49552 * (1.0 / 3600.0)); 336 337 /* Equ 21.1 */ 338 T =(JD - 2451545) / 36525; 339 T2 = T * T; 340 341 /* suns longitude in radians */ 342 ln_get_solar_geom_coords(JD, sol_position); 343 true_longitude = ln_deg_to_rad(sol_position.B); 344 345 /* Earth orbit ecentricity */ 346 e = 0.016708617 - 0.000042037 * T - 0.0000001236 * T2; 347 e = ln_deg_to_rad(e); 348 349 /* longitude of perihelion Earths orbit */ 350 t = 102.93735 + 1.71953 * T + 0.000046 * T2; 351 t = ln_deg_to_rad(t); 352 353 /* change object long/lat to radians */ 354 mean_lng = ln_deg_to_rad(mean_position.lng); 355 mean_lat = ln_deg_to_rad(mean_position.lat); 356 357 /* equ 22.2 */ 358 delta_lng = (-k * cos(true_longitude - mean_lng) 359 + e * k * cos(t - mean_lng)) / cos(mean_lat); 360 delta_lat = -k * sin(mean_lat) * (sin(true_longitude - mean_lng) 361 - e * sin(t - mean_lng)); 362 363 mean_lng += delta_lng; 364 mean_lat += delta_lat; 365 366 position.lng = ln_rad_to_deg(mean_lng); 367 position.lat = ln_rad_to_deg(mean_lat); 368 } 369 370 }