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 */ 18 19 module nova.elliptic_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 /* number of steps in calculation, 3.32 steps for each significant 32 digit required */ 33 enum KEPLER_STEPS = 53; 34 35 extern (C) { 36 37 /* the BASIC SGN() function for doubles */ 38 static @nogc double sgn(double x) nothrow 39 { 40 if (x == 0.0) 41 return (x); 42 else 43 if (x < 0.0) 44 return (-1.0); 45 else 46 return (1.0); 47 } 48 49 /*! \fn double ln_solve_kepler (double E, double M); 50 * \param E Orbital eccentricity 51 * \param M Mean anomaly 52 * \return Eccentric anomaly 53 * 54 * Calculate the eccentric anomaly. 55 * This method was devised by Roger Sinnott. (Sky and Telescope, Vol 70, pg 159) 56 */ 57 @nogc double ln_solve_kepler(double e, double M) nothrow 58 { 59 double Eo = PI_2; 60 double F, M1; 61 double D = PI_4; 62 int i; 63 64 /* covert to radians */ 65 M = ln_deg_to_rad(M); 66 67 F = sgn(M); 68 M = fabs(M) / (2.0 * PI); 69 M = (M - cast(int)M) * 2.0 * PI * F; 70 71 if (M < 0) 72 M = M + 2.0 * PI; 73 F = 1.0; 74 75 if (M > PI) 76 F = -1.0; 77 78 if (M > PI) 79 M = 2.0 * PI - M; 80 81 for (i = 0; i < KEPLER_STEPS; i++) { 82 M1 = Eo - e * sin(Eo); 83 Eo = Eo + D * sgn(M - M1); 84 D /= 2.0; 85 } 86 Eo *= F; 87 88 /* back to degrees */ 89 Eo = ln_rad_to_deg(Eo); 90 return Eo; 91 } 92 93 /*! \fn double ln_get_ell_mean_anomaly (double n, double delta_JD); 94 * \param n Mean motion (degrees/day) 95 * \param delta_JD Time since perihelion 96 * \return Mean anomaly (degrees) 97 * 98 * Calculate the mean anomaly. 99 */ 100 @nogc double ln_get_ell_mean_anomaly(double n, double delta_JD) nothrow 101 { 102 return delta_JD * n; 103 } 104 105 /*! \fn double ln_get_ell_true_anomaly (double e, double E); 106 * \param e Orbital eccentricity 107 * \param E Eccentric anomaly 108 * \return True anomaly (degrees) 109 * 110 * Calculate the true anomaly. 111 */ 112 /* equ 30.1 */ 113 @nogc double ln_get_ell_true_anomaly(double e, double E) nothrow 114 { 115 double v; 116 117 E = ln_deg_to_rad(E); 118 v = sqrt((1.0 + e) / (1.0 - e)) * tan(E / 2.0); 119 v = 2.0 * atan(v); 120 v = ln_range_degrees(ln_rad_to_deg(v)); 121 return v; 122 } 123 124 /*! \fn double ln_get_ell_radius_vector (double a, double e, double E); 125 * \param a Semi-Major axis in AU 126 * \param e Orbital eccentricity 127 * \param E Eccentric anomaly 128 * \return Radius vector AU 129 * 130 * Calculate the radius vector. 131 */ 132 /* equ 30.2 */ 133 @nogc double ln_get_ell_radius_vector(double a, double e, double E) nothrow 134 { 135 return a * (1.0 - e * cos(ln_deg_to_rad(E))); 136 } 137 138 139 /*! \fn double ln_get_ell_smajor_diam (double e, double q); 140 * \param e Eccentricity 141 * \param q Perihelion distance in AU 142 * \return Semi-major diameter in AU 143 * 144 * Calculate the semi major diameter. 145 */ 146 @nogc double ln_get_ell_smajor_diam(double e, double q) nothrow 147 { 148 return q / (1.0 - e); 149 } 150 151 /*! \fn double ln_get_ell_sminor_diam (double e, double a); 152 * \param e Eccentricity 153 * \param a Semi-Major diameter in AU 154 * \return Semi-minor diameter in AU 155 * 156 * Calculate the semi minor diameter. 157 */ 158 @nogc double ln_get_ell_sminor_diam(double e, double a) nothrow 159 { 160 return a * sqrt(1 - e * e); 161 } 162 163 /*! \fn double ln_get_ell_mean_motion (double a); 164 * \param a Semi major diameter in AU 165 * \return Mean daily motion (degrees/day) 166 * 167 * Calculate the mean daily motion (degrees/day). 168 */ 169 @nogc double ln_get_ell_mean_motion(double a) nothrow 170 { 171 double q = 0.9856076686; /* Gaussian gravitational constant (degrees)*/ 172 return q / (a * sqrt(a)); 173 } 174 175 /*! \fn void ln_get_ell_helio_rect_posn(struct ln_ell_orbit *orbit, double JD, struct ln_rect_posn *posn); 176 * \param orbit Orbital parameters of object. 177 * \param JD Julian day 178 * \param posn Position pointer to store objects position 179 * 180 * Calculate the objects rectangular heliocentric position given it's orbital 181 * elements for the given julian day. 182 */ 183 @nogc void ln_get_ell_helio_rect_posn(ref ln_ell_orbit orbit, double JD, 184 ref ln_rect_posn posn) nothrow 185 { 186 double A,B,C; 187 double F,G,H; 188 double P,Q,R; 189 double sin_e, cos_e; 190 double a,b,c; 191 double sin_omega, sin_i, cos_omega, cos_i; 192 double M,v,E,r; 193 194 /* J2000 obliquity of the ecliptic */ 195 sin_e = 0.397777156; 196 cos_e = 0.917482062; 197 198 /* equ 33.7 */ 199 sin_omega = sin(ln_deg_to_rad(orbit.omega)); 200 cos_omega = cos(ln_deg_to_rad(orbit.omega)); 201 sin_i = sin(ln_deg_to_rad(orbit.i)); 202 cos_i = cos(ln_deg_to_rad(orbit.i)); 203 F = cos_omega; 204 G = sin_omega * cos_e; 205 H = sin_omega * sin_e; 206 P = -sin_omega * cos_i; 207 Q = cos_omega * cos_i * cos_e - sin_i * sin_e; 208 R = cos_omega * cos_i * sin_e + sin_i * cos_e; 209 210 /* equ 33.8 */ 211 A = atan2(F,P); 212 B = atan2(G,Q); 213 C = atan2(H,R); 214 a = sqrt(F * F + P * P); 215 b = sqrt(G * G + Q * Q); 216 c = sqrt(H * H + R * R); 217 218 /* get daily motion */ 219 if (orbit.n == 0.0) 220 orbit.n = ln_get_ell_mean_motion(orbit.a); 221 222 /* get mean anomaly */ 223 M = ln_get_ell_mean_anomaly(orbit.n, JD - orbit.JD); 224 225 /* get eccentric anomaly */ 226 E = ln_solve_kepler(orbit.e, M); 227 228 /* get true anomaly */ 229 v = ln_get_ell_true_anomaly(orbit.e, E); 230 231 /* get radius vector */ 232 r = ln_get_ell_radius_vector(orbit.a, orbit.e, E); 233 234 /* equ 33.9 */ 235 posn.X = r * a * sin(A + ln_deg_to_rad(orbit.w + v)); 236 posn.Y = r * b * sin(B + ln_deg_to_rad(orbit.w + v)); 237 posn.Z = r * c * sin(C + ln_deg_to_rad(orbit.w + v)); 238 } 239 240 /*! \fn void ln_get_ell_geo_rect_posn(struct ln_ell_orbit *orbit, double JD, struct ln_rect_posn *posn); 241 * \param orbit Orbital parameters of object. 242 * \param JD Julian day 243 * \param posn Position pointer to store objects position 244 * 245 * Calculate the objects rectangular geocentric position given it's orbital 246 * elements for the given julian day. 247 */ 248 @nogc void ln_get_ell_geo_rect_posn(ref ln_ell_orbit orbit, double JD, 249 ref ln_rect_posn posn) nothrow 250 { 251 ln_rect_posn p_posn, e_posn; 252 ln_helio_posn earth; 253 254 /* elliptic helio rect coords */ 255 ln_get_ell_helio_rect_posn(orbit, JD, p_posn); 256 257 /* earth rect coords */ 258 ln_get_earth_helio_coords(JD, earth); 259 ln_get_rect_from_helio(earth, e_posn); 260 261 posn.X = e_posn.X - p_posn.X; 262 posn.Y = e_posn.Y - p_posn.Y; 263 posn.Z = e_posn.Z - p_posn.Z; 264 } 265 266 267 /*! 268 * \fn void ln_get_ell_body_equ_coords(double JD, struct ln_ell_orbit *orbit, struct ln_equ_posn *posn) 269 * \param JD Julian Day. 270 * \param orbit Orbital parameters. 271 * \param posn Pointer to hold asteroid position. 272 * 273 * Calculate a body's equatorial coordinates for the given julian day. 274 */ 275 @nogc void ln_get_ell_body_equ_coords(double JD, ref ln_ell_orbit orbit, 276 ref ln_equ_posn posn) nothrow 277 { 278 ln_rect_posn body_rect_posn, sol_rect_posn; 279 double dist, t; 280 double x,y,z; 281 282 /* get solar and body rect coords */ 283 ln_get_ell_helio_rect_posn(orbit, JD, body_rect_posn); 284 ln_get_solar_geo_coords(JD, sol_rect_posn); 285 286 /* calc distance and light time */ 287 dist = ln_get_rect_distance(body_rect_posn, sol_rect_posn); 288 t = ln_get_light_time(dist); 289 290 /* repeat calculation with new time (i.e. JD - t) */ 291 ln_get_ell_helio_rect_posn(orbit, JD - t, body_rect_posn); 292 293 /* calc equ coords equ 33.10 */ 294 x = sol_rect_posn.X + body_rect_posn.X; 295 y = sol_rect_posn.Y + body_rect_posn.Y; 296 z = sol_rect_posn.Z + body_rect_posn.Z; 297 298 posn.ra = ln_range_degrees(ln_rad_to_deg(atan2(y,x))); 299 posn.dec = ln_rad_to_deg(asin(z / sqrt(x * x + y * y + z * z))); 300 } 301 302 303 /*! \fn double ln_get_ell_orbit_len(struct ln_ell_orbit *orbit); 304 * \param orbit Orbital parameters 305 * \return Orbital length in AU 306 * 307 * Calculate the orbital length in AU. 308 * 309 * Accuracy: 310 * - 0.001% for e < 0.88 311 * - 0.01% for e < 0.95 312 * - 1% for e = 0.9997 313 * - 3% for e = 1 314 */ 315 @nogc double ln_get_ell_orbit_len(const ref ln_ell_orbit orbit) nothrow 316 { 317 double A,G,H; 318 double b; 319 320 b = ln_get_ell_sminor_diam(orbit.e, orbit.a); 321 322 A = (orbit.a + b) / 2.0; 323 G = sqrt(orbit.a * b); 324 H = (2.0 * orbit.a * b) / (orbit.a + b); 325 326 /* Meeus, page 239, 2nd edition */ 327 return PI * ((21.0 * A - 2.0 * G - 3.0 * H) / 8.0); 328 } 329 330 /*! \fn double ln_get_ell_orbit_vel(double JD, struct ln_ell_orbit *orbit); 331 * \param JD Julian day. 332 * \param orbit Orbital parameters 333 * \return Orbital velocity in km/s. 334 * 335 * Calculate orbital velocity in km/s for the given julian day. 336 */ 337 @nogc double ln_get_ell_orbit_vel(double JD, ref ln_ell_orbit orbit) nothrow 338 { 339 double V; 340 double r; 341 342 r = ln_get_ell_body_solar_dist(JD, orbit); 343 V = 1.0 / r - 1.0 / (2.0 * orbit.a); 344 V = 42.1219 * sqrt(V); 345 return V; 346 } 347 348 /*! \fn double ln_get_ell_orbit_pvel(struct ln_ell_orbit *orbit); 349 * \param orbit Orbital parameters 350 * \return Orbital velocity in km/s. 351 * 352 * Calculate orbital velocity at perihelion in km/s. 353 */ 354 @nogc double ln_get_ell_orbit_pvel(const ref ln_ell_orbit orbit) nothrow 355 { 356 double V; 357 358 V = 29.7847 / sqrt(orbit.a); 359 V *= sqrt((1.0 + orbit.e) / (1.0 - orbit.e)); 360 return V; 361 } 362 363 /*! \fn double ln_get_ell_orbit_avel(struct ln_ell_orbit *orbit); 364 * \param orbit Orbital parameters 365 * \return Orbital velocity in km/s. 366 * 367 * Calculate the orbital velocity at aphelion in km/s. 368 */ 369 @nogc double ln_get_ell_orbit_avel(const ref ln_ell_orbit orbit) nothrow 370 { 371 double V; 372 373 V = 29.7847 / sqrt(orbit.a); 374 V *= sqrt((1.0 - orbit.e) / (1.0 + orbit.e)); 375 return V; 376 } 377 378 /*! 379 * \fn double ln_get_ell_body_solar_dist(double JD, struct ln_ell_orbit *orbit) 380 * \param JD Julian Day. 381 * \param orbit Orbital parameters 382 * \return The distance in AU between the Sun and the body. 383 * 384 * Calculate the distance between a body and the Sun. 385 */ 386 @nogc double ln_get_ell_body_solar_dist(double JD, ref ln_ell_orbit orbit) nothrow 387 { 388 ln_rect_posn body_rect_posn, sol_rect_posn; 389 390 /* get solar and body rect coords */ 391 ln_get_ell_helio_rect_posn(orbit, JD, body_rect_posn); 392 sol_rect_posn.X = 0; 393 sol_rect_posn.Y = 0; 394 sol_rect_posn.Z = 0; 395 396 /* calc distance */ 397 return ln_get_rect_distance(body_rect_posn, sol_rect_posn); 398 } 399 400 /*! 401 * \fn double ln_get_ell_body_earth_dist(double JD, struct ln_ell_orbit *orbit) 402 * \param JD Julian day. 403 * \param orbit Orbital parameters 404 * \returns Distance in AU 405 * 406 * Calculate the distance between an body and the Earth 407 * for the given julian day. 408 */ 409 @nogc double ln_get_ell_body_earth_dist(double JD, ref ln_ell_orbit orbit) nothrow 410 { 411 ln_rect_posn body_rect_posn, earth_rect_posn; 412 413 /* get solar and body rect coords */ 414 ln_get_ell_geo_rect_posn(orbit, JD, body_rect_posn); 415 earth_rect_posn.X = 0; 416 earth_rect_posn.Y = 0; 417 earth_rect_posn.Z = 0; 418 419 /* calc distance */ 420 return ln_get_rect_distance(body_rect_posn, earth_rect_posn); 421 } 422 423 424 /*! \fn double ln_get_ell_body_phase_angle(double JD, struct ln_ell_orbit *orbit); 425 * \param JD Julian day 426 * \param orbit Orbital parameters 427 * \return Phase angle. 428 * 429 * Calculate the phase angle of the body. The angle Sun - body - Earth. 430 */ 431 @nogc double ln_get_ell_body_phase_angle(double JD, ref ln_ell_orbit orbit) nothrow 432 { 433 double r,R,d; 434 double E,M; 435 double phase; 436 437 /* get mean anomaly */ 438 if (orbit.n == 0.0) 439 orbit.n = ln_get_ell_mean_motion(orbit.a); 440 M = ln_get_ell_mean_anomaly(orbit.n, JD - orbit.JD); 441 442 /* get eccentric anomaly */ 443 E = ln_solve_kepler(orbit.e, M); 444 445 /* get radius vector */ 446 r = ln_get_ell_radius_vector(orbit.a, orbit.e, E); 447 448 /* get solar and Earth distances */ 449 R = ln_get_ell_body_earth_dist(JD, orbit); 450 d = ln_get_ell_body_solar_dist(JD, orbit); 451 452 phase = (r * r + d * d - R * R) / ( 2.0 * r * d); 453 return ln_range_degrees(acos(ln_deg_to_rad(phase))); 454 } 455 456 457 /*! \fn double ln_get_ell_body_elong(double JD, struct ln_ell_orbit *orbit); 458 * \param JD Julian day 459 * \param orbit Orbital parameters 460 * \return Elongation to the Sun. 461 * 462 * Calculate the body's elongation to the Sun. 463 */ 464 @nogc double ln_get_ell_body_elong(double JD, ref ln_ell_orbit orbit) nothrow 465 { 466 double r,R,d; 467 double t; 468 double elong; 469 double E,M; 470 471 /* time since perihelion */ 472 t = JD - orbit.JD; 473 474 /* get mean anomaly */ 475 if (orbit.n == 0.0) 476 orbit.n = ln_get_ell_mean_motion(orbit.a); 477 M = ln_get_ell_mean_anomaly(orbit.n, t); 478 479 /* get eccentric anomaly */ 480 E = ln_solve_kepler(orbit.e, M); 481 482 /* get radius vector */ 483 r = ln_get_ell_radius_vector(orbit.a, orbit.e, E); 484 485 /* get solar and Earth-Sun distances */ 486 R = ln_get_earth_solar_dist(JD); 487 d = ln_get_ell_body_solar_dist(JD, orbit); 488 489 elong = (R * R + d * d - r * r) / ( 2.0 * R * d); 490 return ln_range_degrees(ln_rad_to_deg(acos(elong))); 491 } 492 493 /*! \fn double ln_get_ell_body_rst(double JD, struct ln_lnlat_posn *observer, struct ln_ell_orbit *orbit, struct ln_rst_time *rst); 494 * \param JD Julian day 495 * \param observer Observers position 496 * \param orbit Orbital parameters 497 * \param rst Pointer to store Rise, Set and Transit time in JD 498 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (below the horizon) 499 * 500 * Calculate the time the rise, set and transit (crosses the local meridian at upper culmination) 501 * time of a body with an elliptic orbit for the given Julian day. 502 * 503 * Assumes a "standard" horizon. 504 * 505 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 506 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 507 */ 508 @nogc int ln_get_ell_body_rst(double JD, const ref ln_lnlat_posn observer, 509 ref ln_ell_orbit orbit, ref ln_rst_time rst) nothrow 510 { 511 return ln_get_ell_body_rst_horizon(JD, observer, orbit, 512 LN_STAR_STANDART_HORIZON, rst); 513 } 514 515 /*! \fn double ln_get_ell_body_rst_horizon(double JD, struct ln_lnlat_posn *observer, struct ln_ell_orbit *orbit, double horizon, struct ln_rst_time *rst); 516 * \param JD Julian day 517 * \param observer Observers position 518 * \param orbit Orbital parameters 519 * \param horizon Horizon height 520 * \param rst Pointer to store Rise, Set and Transit time in JD 521 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (below the horizon) 522 * 523 * Calculate the time the rise, set and transit (crosses the local meridian at upper culmination) 524 * time of a body with an elliptic orbit for the given Julian day. 525 * 526 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 527 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 528 */ 529 @nogc int ln_get_ell_body_rst_horizon(double JD, const ref ln_lnlat_posn observer, 530 ref ln_ell_orbit orbit, double horizon, ref ln_rst_time rst) nothrow 531 { 532 return ln_get_motion_body_rst_horizon(JD, observer, 533 cast(get_motion_body_coords_t)&ln_get_ell_body_equ_coords, 534 &orbit, horizon, rst); 535 } 536 537 /*! \fn double ln_get_ell_body_next_rst(double JD, struct ln_lnlat_posn *observer, struct ln_ell_orbit *orbit, struct ln_rst_time *rst); 538 * \param JD Julian day 539 * \param observer Observers position 540 * \param orbit Orbital parameters 541 * \param rst Pointer to store Rise, Set and Transit time in JD 542 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (below the horizon) 543 * 544 * Calculate the time of next rise, set and transit (crosses the local meridian at upper culmination) 545 * time of a body with an elliptic orbit for the given Julian day. 546 * 547 * This function guarantee, that rise, set and transit will be in <JD, JD+1> range. 548 * 549 * Assumes a "standard" horizon. 550 * 551 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 552 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 553 */ 554 @nogc int ln_get_ell_body_next_rst(double JD, const ref ln_lnlat_posn observer, 555 ref ln_ell_orbit orbit, ref ln_rst_time rst) nothrow 556 { 557 return ln_get_ell_body_next_rst_horizon(JD, observer, orbit, 558 LN_STAR_STANDART_HORIZON, rst); 559 } 560 561 /*! \fn double ln_get_ell_body_next_rst_horizon(double JD, struct ln_lnlat_posn *observer, struct ln_ell_orbit *orbit, double horizon, struct ln_rst_time *rst); 562 * \param JD Julian day 563 * \param observer Observers position 564 * \param orbit Orbital parameters 565 * \param horizon Horizon height 566 * \param rst Pointer to store Rise, Set and Transit time in JD 567 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (below the horizon) 568 * 569 * Calculate the time of next rise, set and transit (crosses the local meridian at upper culmination) 570 * time of a body with an elliptic orbit for the given Julian day. 571 * 572 * This function guarantee, that rise, set and transit will be in <JD, JD+1> range. 573 * 574 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 575 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 576 */ 577 @nogc int ln_get_ell_body_next_rst_horizon(double JD, const ref ln_lnlat_posn observer, 578 ref ln_ell_orbit orbit, double horizon, ref ln_rst_time rst) nothrow 579 { 580 return ln_get_motion_body_next_rst_horizon(JD, observer, 581 cast(get_motion_body_coords_t)&ln_get_ell_body_equ_coords, 582 &orbit, horizon, rst); 583 } 584 585 /*! \fn double ln_get_ell_body_next_rst_horizon(double JD, struct ln_lnlat_posn *observer, struct ln_ell_orbit *orbit, double horizon, struct ln_rst_time *rst); 586 * \param JD Julian day 587 * \param observer Observers position 588 * \param orbit Orbital parameters 589 * \param horizon Horizon height 590 * \param day_limit Maximal number of days that will be searched for next rise and set 591 * \param rst Pointer to store Rise, Set and Transit time in JD 592 * \return 0 for success, else 1 for circumpolar (above the horizon), -1 for circumpolar (below the horizon) 593 * 594 * Calculate the time of next rise, set and transit (crosses the local meridian at upper culmination) 595 * time of a body with an elliptic orbit for the given Julian day. 596 * 597 * This function guarantee, that rise, set and transit will be in <JD, JD + day_limit> range. 598 * 599 * Note: this functions returns 1 if the body is circumpolar, that is it remains the whole 600 * day above the horizon. Returns -1 when it remains the whole day below the horizon. 601 */ 602 @nogc int ln_get_ell_body_next_rst_horizon_future(double JD, 603 const ref ln_lnlat_posn observer, ref ln_ell_orbit orbit, 604 double horizon, int day_limit, ref ln_rst_time rst) nothrow 605 { 606 return ln_get_motion_body_next_rst_horizon_future(JD, observer, 607 cast(get_motion_body_coords_t)&ln_get_ell_body_equ_coords, 608 &orbit, horizon, day_limit, rst); 609 } 610 611 /*!\fn double ln_get_ell_last_perihelion (double epoch_JD, double M, double n); 612 * \param epoch_JD Julian day of epoch 613 * \param M Mean anomaly 614 * \param n daily motion in degrees 615 * 616 * Calculate the julian day of the last perihelion. 617 */ 618 @nogc double ln_get_ell_last_perihelion (double epoch_JD, double M, double n) nothrow 619 { 620 return epoch_JD - (M / n); 621 } 622 623 }