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.parabolic_motion; 20 21 import std.math; 22 23 import nova.solar; 24 import nova.earth; 25 import nova.transform; 26 import nova.rise_set; 27 import nova.dynamical_time; 28 import nova.sidereal_time; 29 import nova.utility; 30 31 extern (C) { 32 33 /*! \fn double ln_solve_barker (double q, double t); 34 * \param q Perihelion distance in AU 35 * \param t Time since perihelion in days 36 * \return Solution of Barkers equation 37 * 38 * Solve Barkers equation. LIAM add more 39 */ 40 /* Equ 34.3, Barkers Equation */ 41 @nogc double ln_solve_barker(double q, double t) nothrow 42 { 43 double G, Y, W; 44 45 /* equ 34.1 */ 46 W = ((0.03649116245) / (q * sqrt(q))) * t; 47 48 /* equ 34.6 */ 49 G = W / 2.0; 50 Y = cbrt(G + sqrt(G * G + 1.0)); 51 return Y - 1.0 / Y; 52 } 53 54 /*! \fn double ln_get_par_true_anomaly (double q, double t); 55 * \param q Perihelion distance in AU 56 * \param t Time since perihelion 57 * \return True anomaly (degrees) 58 * 59 * Calculate the true anomaly. 60 */ 61 /* equ 30.1 */ 62 @nogc double ln_get_par_true_anomaly(double q, double t) nothrow 63 { 64 double v, s; 65 66 s = ln_solve_barker(q, t); 67 v = 2.0 * atan(s); 68 69 return ln_range_degrees(ln_rad_to_deg(v)); 70 } 71 72 /*! \fn double ln_get_par_radius_vector (double q, double t); 73 * \param q Perihelion distance in AU 74 * \param t Time since perihelion in days 75 * \return Radius vector AU 76 * 77 * Calculate the radius vector. 78 */ 79 /* equ 30.2 */ 80 @nogc double ln_get_par_radius_vector(double q, double t) nothrow 81 { 82 double s; 83 84 s = ln_solve_barker(q, t); 85 return q * (1.0 + s * s); 86 } 87 88 89 /*! \fn void ln_get_par_helio_rect_posn(struct ln_par_orbit *orbit, double JD, struct ln_rect_posn *posn); 90 * \param orbit Orbital parameters of object. 91 * \param JD Julian day 92 * \param posn Position pointer to store objects position 93 * 94 * Calculate the objects rectangular heliocentric position given it's orbital 95 * elements for the given julian day. 96 */ 97 @nogc void ln_get_par_helio_rect_posn(ref ln_par_orbit orbit, double JD, 98 ref ln_rect_posn posn) nothrow 99 { 100 double A, B, C, F, G, H, P, Q, R; 101 double sin_e, cos_e; 102 double a, b, c; 103 double sin_omega, sin_i, cos_omega, cos_i; 104 double r, v, t; 105 106 /* time since perihelion */ 107 t = JD - orbit.JD; 108 109 /* J2000 obliquity of the ecliptic */ 110 sin_e = 0.397777156; 111 cos_e = 0.917482062; 112 113 /* equ 33.7 */ 114 sin_omega = sin(ln_deg_to_rad(orbit.omega)); 115 cos_omega = cos(ln_deg_to_rad(orbit.omega)); 116 sin_i = sin(ln_deg_to_rad(orbit.i)); 117 cos_i = cos(ln_deg_to_rad(orbit.i)); 118 F = cos_omega; 119 G = sin_omega * cos_e; 120 H = sin_omega * sin_e; 121 P = -sin_omega * cos_i; 122 Q = cos_omega * cos_i * cos_e - sin_i * sin_e; 123 R = cos_omega * cos_i * sin_e + sin_i * cos_e; 124 125 /* equ 33.8 */ 126 A = atan2(F, P); 127 B = atan2(G, Q); 128 C = atan2(H, R); 129 a = sqrt(F * F + P * P); 130 b = sqrt(G * G + Q * Q); 131 c = sqrt(H * H + R * R); 132 133 /* get true anomaly */ 134 v = ln_get_par_true_anomaly(orbit.q, t); 135 136 /* get radius vector */ 137 r = ln_get_par_radius_vector(orbit.q, t); 138 139 /* equ 33.9 */ 140 posn.X = r * a * sin(A + ln_deg_to_rad(orbit.w + v)); 141 posn.Y = r * b * sin(B + ln_deg_to_rad(orbit.w + v)); 142 posn.Z = r * c * sin(C + ln_deg_to_rad(orbit.w + v)); 143 } 144 145 146 /*! \fn void ln_get_par_geo_rect_posn(struct ln_par_orbit *orbit, double JD, struct ln_rect_posn *posn); 147 * \param orbit Orbital parameters of object. 148 * \param JD Julian day 149 * \param posn Position pointer to store objects position 150 * 151 * Calculate the objects rectangular geocentric position given it's orbital 152 * elements for the given julian day. 153 */ 154 @nogc void ln_get_par_geo_rect_posn(ref ln_par_orbit orbit, double JD, 155 ref ln_rect_posn posn) nothrow 156 { 157 ln_rect_posn p_posn, e_posn; 158 ln_helio_posn earth; 159 160 /* parabolic helio rect coords */ 161 ln_get_par_helio_rect_posn(orbit, JD, p_posn); 162 163 /* earth rect coords */ 164 ln_get_earth_helio_coords(JD, earth); 165 166 ln_get_rect_from_helio(earth, e_posn); 167 posn.X = p_posn.X - e_posn.X; 168 posn.Y = p_posn.Y - e_posn.Y; 169 posn.Z = p_posn.Z - e_posn.Z; 170 } 171 172 173 /*! 174 * \fn void ln_get_par_body_equ_coords(double JD, struct ln_par_orbit *orbit, struct ln_equ_posn *posn) 175 * \param JD Julian Day. 176 * \param orbit Orbital parameters. 177 * \param posn Pointer to hold asteroid position. 178 * 179 * Calculate a bodies equatorial coordinates for the given julian day. 180 */ 181 @nogc void ln_get_par_body_equ_coords(double JD, ref ln_par_orbit orbit, 182 ref ln_equ_posn posn) nothrow 183 { 184 ln_rect_posn body_rect_posn, sol_rect_posn; 185 double dist, t; 186 double x,y,z; 187 188 /* get solar and body rect coords */ 189 ln_get_par_helio_rect_posn (orbit, JD, body_rect_posn); 190 ln_get_solar_geo_coords(JD, sol_rect_posn); 191 192 /* calc distance and light time */ 193 dist = ln_get_rect_distance (body_rect_posn, sol_rect_posn); 194 t = ln_get_light_time (dist); 195 196 /* repeat calculation with new time (i.e. JD - t) */ 197 ln_get_par_helio_rect_posn (orbit, JD - t, body_rect_posn); 198 199 /* calc equ coords equ 33.10 */ 200 x = sol_rect_posn.X + body_rect_posn.X; 201 y = sol_rect_posn.Y + body_rect_posn.Y; 202 z = sol_rect_posn.Z + body_rect_posn.Z; 203 204 posn.ra = ln_range_degrees(ln_rad_to_deg(atan2(y, x))); 205 posn.dec = ln_rad_to_deg(atan2(z,sqrt(x * x + y * y))); 206 } 207 208 209 /*! 210 * \fn double ln_get_par_body_earth_dist(double JD, struct ln_par_orbit *orbit) 211 * \param JD Julian day. 212 * \param orbit Orbital parameters 213 * \returns Distance in AU 214 * 215 * Calculate the distance between a body and the Earth 216 * for the given julian day. 217 */ 218 @nogc double ln_get_par_body_earth_dist(double JD, ref ln_par_orbit orbit) nothrow 219 { 220 ln_rect_posn body_rect_posn, earth_rect_posn; 221 222 /* get solar and body rect coords */ 223 ln_get_par_geo_rect_posn(orbit, JD, body_rect_posn); 224 earth_rect_posn.X = 0; 225 earth_rect_posn.Y = 0; 226 earth_rect_posn.Z = 0; 227 228 /* calc distance */ 229 return ln_get_rect_distance(body_rect_posn, earth_rect_posn); 230 } 231 232 /*! 233 * \fn double ln_get_par_body_solar_dist(double JD, struct ln_par_orbit *orbit) 234 * \param JD Julian Day. 235 * \param orbit Orbital parameters 236 * \return The distance in AU between the Sun and the body. 237 * 238 * Calculate the distance between a body and the Sun. 239 */ 240 @nogc double ln_get_par_body_solar_dist(double JD, ref ln_par_orbit orbit) nothrow 241 { 242 ln_rect_posn body_rect_posn, sol_rect_posn; 243 244 /* get solar and body rect coords */ 245 ln_get_par_helio_rect_posn(orbit, JD, body_rect_posn); 246 sol_rect_posn.X = 0; 247 sol_rect_posn.Y = 0; 248 sol_rect_posn.Z = 0; 249 250 /* calc distance */ 251 return ln_get_rect_distance(body_rect_posn, sol_rect_posn); 252 } 253 254 /*! \fn double ln_get_par_body_phase_angle(double JD, struct ln_par_orbit *orbit); 255 * \param JD Julian day 256 * \param orbit Orbital parameters 257 * \return Phase angle. 258 * 259 * Calculate the phase angle of the body. The angle Sun - body - Earth. 260 */ 261 @nogc double ln_get_par_body_phase_angle(double JD, ref ln_par_orbit orbit) nothrow 262 { 263 double r,R,d; 264 double t; 265 double phase; 266 267 /* time since perihelion */ 268 t = JD - orbit.JD; 269 270 /* get radius vector */ 271 r = ln_get_par_radius_vector(orbit.q, t); 272 273 /* get solar and Earth-Sun distances */ 274 R = ln_get_earth_solar_dist(JD); 275 d = ln_get_par_body_solar_dist(JD, orbit); 276 277 phase = (r * r + d * d - R * R) / (2.0 * r * d ); 278 return ln_range_degrees(ln_rad_to_deg(acos(phase))); 279 } 280 281 /*! \fn double ln_get_par_body_elong(double JD, struct ln_par_orbit *orbit); 282 * \param JD Julian day 283 * \param orbit Orbital parameters 284 * \return Elongation to the Sun. 285 * 286 * Calculate the bodies elongation to the Sun.. 287 */ 288 @nogc double ln_get_par_body_elong(double JD, ref ln_par_orbit orbit) nothrow 289 { 290 double r,R,d; 291 double t; 292 double elong; 293 294 /* time since perihelion */ 295 t = JD - orbit.JD; 296 297 /* get radius vector */ 298 r = ln_get_par_radius_vector (orbit.q, t); 299 300 /* get solar and Earth-Sun distances */ 301 R = ln_get_earth_solar_dist(JD); 302 d = ln_get_par_body_solar_dist(JD, orbit); 303 304 elong = (R * R + d * d - r * r) / ( 2.0 * R * d ); 305 return ln_range_degrees(ln_rad_to_deg(acos(elong))); 306 } 307 308 /*! \fn double ln_get_par_body_rst(double JD, struct ln_lnlat_posn *observer, struct ln_par_orbit *orbit, struct ln_rst_time *rst); 309 * \param JD Julian day 310 * \param observer Observers position 311 * \param orbit Orbital parameters 312 * \param rst Pointer to store Rise, Set and Transit time in JD 313 * \return 0 for success, 1 for circumpolar (above the horizon), -1 for circumpolar (bellow the horizon) 314 * 315 * Calculate the time the rise, set and transit (crosses the local meridian at upper culmination) 316 * time of a body with a parabolic orbit for the given Julian day. 317 * 318 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 319 * day either above the horizon. Returns -1 when it remains whole day below the horizon. 320 */ 321 @nogc int ln_get_par_body_rst(double JD, const ref ln_lnlat_posn observer, 322 ref ln_par_orbit orbit, ref ln_rst_time rst) nothrow 323 { 324 return ln_get_par_body_rst_horizon(JD, observer, orbit, 325 LN_STAR_STANDART_HORIZON, rst); 326 } 327 328 /*! \fn double ln_get_par_body_rst_horizon(double JD, struct ln_lnlat_posn *observer, struct ln_par_orbit *orbit, double horizon, struct ln_rst_time *rst); 329 * \param JD Julian day 330 * \param observer Observers position 331 * \param orbit Orbital parameters 332 * \param horizon Horizon height 333 * \param rst Pointer to store Rise, Set and Transit time in JD 334 * \return 0 for success, else 1 for circumpolar. 335 * 336 * Calculate the time the rise, set and transit (crosses the local meridian at upper culmination) 337 * time of a body with a parabolic orbit for the given Julian day. 338 * 339 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 340 * day either above the horizon. Returns -1 when it remains whole day below the horizon. 341 */ 342 @nogc int ln_get_par_body_rst_horizon(double JD, const ref ln_lnlat_posn observer, 343 ref ln_par_orbit orbit, double horizon, ref ln_rst_time rst) nothrow 344 { 345 return ln_get_motion_body_rst_horizon(JD, observer, 346 cast(get_motion_body_coords_t)&ln_get_par_body_equ_coords, &orbit, 347 horizon, rst); 348 } 349 350 /*! \fn double ln_get_par_body_next_rst(double JD, struct ln_lnlat_posn *observer, struct ln_par_orbit *orbit, struct ln_rst_time *rst); 351 * \param JD Julian day 352 * \param observer Observers position 353 * \param orbit Orbital parameters 354 * \param rst Pointer to store Rise, Set and Transit time in JD 355 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (bellow the horizon) 356 * 357 * Calculate the time of next rise, set and transit (crosses the local meridian at upper culmination) 358 * time of a body with an parabolic orbit for the given Julian day. 359 * 360 * This function guarantee, that rise, set and transit will be in <JD, JD+1> range. 361 * 362 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 363 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 364 */ 365 @nogc int ln_get_par_body_next_rst(double JD, const ref ln_lnlat_posn observer, 366 ref ln_par_orbit orbit, ref ln_rst_time rst) nothrow 367 { 368 return ln_get_par_body_next_rst_horizon(JD, observer, orbit, 369 LN_STAR_STANDART_HORIZON, rst); 370 } 371 372 /*! \fn double ln_get_par_body_next_rst_horizon(double JD, struct ln_lnlat_posn *observer, struct ln_par_orbit *orbit, double horizon, struct ln_rst_time *rst); 373 * \param JD Julian day 374 * \param observer Observers position 375 * \param orbit Orbital parameters 376 * \param horizon Horizon height 377 * \param rst Pointer to store Rise, Set and Transit time in JD 378 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (bellow the horizon) 379 * 380 * Calculate the time of next rise, set and transit (crosses the local meridian at upper culmination) 381 * time of a body with an parabolic orbit for the given Julian day. 382 * 383 * This function guarantee, that rise, set and transit will be in <JD, JD+1> range. 384 * 385 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 386 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 387 */ 388 @nogc int ln_get_par_body_next_rst_horizon(double JD, const ref ln_lnlat_posn observer, 389 ref ln_par_orbit orbit, double horizon, ref ln_rst_time rst) nothrow 390 { 391 return ln_get_motion_body_next_rst_horizon(JD, observer, 392 cast(get_motion_body_coords_t)&ln_get_par_body_equ_coords, &orbit, 393 horizon, rst); 394 } 395 396 /*! \fn double ln_get_par_body_next_rst_horizon_future(double JD, struct ln_lnlat_posn *observer, struct ln_par_orbit *orbit, double horizon, int day_limit, struct ln_rst_time *rst); 397 * \param JD Julian day 398 * \param observer Observers position 399 * \param orbit Orbital parameters 400 * \param horizon Horizon height 401 * \param day_limit Maximal number of days that will be searched for next rise and set 402 * \param rst Pointer to store Rise, Set and Transit time in JD 403 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (bellow the horizon) 404 * 405 * Calculate the time of next rise, set and transit (crosses the local meridian at upper culmination) 406 * time of a body with an parabolic orbit for the given Julian day. 407 * 408 * This function guarantee, that rise, set and transit will be in <JD, JD + day_limit> range. 409 * 410 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 411 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 412 */ 413 @nogc int ln_get_par_body_next_rst_horizon_future(double JD, 414 const ref ln_lnlat_posn observer, ref ln_par_orbit orbit, 415 double horizon, int day_limit, ref ln_rst_time rst) nothrow 416 { 417 return ln_get_motion_body_next_rst_horizon_future(JD, observer, 418 cast(get_motion_body_coords_t)&ln_get_par_body_equ_coords, &orbit, 419 horizon, day_limit, rst); 420 } 421 422 }