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