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