1 /*
2 This program is free software; you can redistribute it and/or modify
3 it under the terms of the GNU General Public License as published by
4 the Free Software Foundation; either version 2 of the License, or
5 (at your option) any later version.
6
7 This program 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
10 GNU 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 2000 Liam Girdwood
17 Copyright 2008-2009 Petr Kubanek*/
18
19 #define _GNU_SOURCE
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <libnova/libnova.h>
24 #ifndef __WIN32__
25 #include <unistd.h>
26 #endif
27
28 #ifdef __WIN32__
29 // the usleep function does not exist in visual studio
usleep(int miliseconds)30 void usleep(int miliseconds)
31 {
32 clock_t endwait;
33 endwait = clock () + miliseconds * CLOCKS_PER_SEC / 1000;
34 while (clock() < endwait) {}
35 }
36
37 #endif
38
39 /*
40 * Define DATE or SYS_DATE for testing VSOP87 theory with other JD
41 */
42 //#define DATE
43 //#define SYS_DATE
44
45 // holds number of tests
46 static int test_number = 0;
47
compare_results(double calc,double expect,double tolerance)48 double compare_results (double calc, double expect, double tolerance)
49 {
50 if (calc - expect > tolerance || calc - expect < (tolerance * -1.0))
51 return (calc - expect);
52 else
53 return (0);
54 }
55
test_result(char * test,double calc,double expect,double tolerance)56 int test_result (char * test, double calc, double expect, double tolerance)
57 {
58 double diff;
59
60 printf ("TEST %s....", test);
61
62 test_number++;
63
64 diff = compare_results (calc, expect, tolerance);
65 if (diff) {
66 printf ("[FAILED]\n");
67 printf (" Expected %8.8f but calculated %8.8f %0.12f error.\n\n", expect, calc, diff);
68 return 1;
69 } else {
70 printf ("[PASSED]\n");
71 printf (" Expected and calculated %8.8f.\n\n", calc);
72 return 0;
73 }
74 }
75
76 /* test julian day calculations */
julian_test(void)77 int julian_test (void)
78 {
79 double JD, JD2;
80 int wday, failed = 0;
81 struct ln_date date, pdate;
82 struct ln_zonedate zonedate;
83
84 time_t now;
85 time_t now_jd;
86
87 /* Get julian day for 04/10/1957 19:00:00 */
88 date.years = 1957;
89 date.months = 10;
90 date.days = 4;
91 date.hours = 19;
92 date.minutes = 0;
93 date.seconds = 0;
94 JD = ln_get_julian_day (&date);
95 failed += test_result ("(Julian Day) JD for 4/10/1957 19:00:00", JD, 2436116.29166667, 0.00001);
96
97 /* Get julian day for 27/01/333 12:00:00 */
98 date.years = 333;
99 date.months = 1;
100 date.days = 27;
101 date.hours = 12;
102 JD = ln_get_julian_day (&date);
103 failed += test_result ("(Julian Day) JD for 27/01/333 12:00:00", JD, 1842713.0, 0.1);
104
105 /* Get julian day for 30/06/1954 00:00:00 */
106 date.years = 1954;
107 date.months = 6;
108 date.days = 30;
109 date.hours = 0;
110 JD = ln_get_julian_day (&date);
111 failed += test_result ("(Julian Day) JD for 30/06/1954 00:00:00", JD, 2434923.5, 0.1);
112
113 wday = ln_get_day_of_week(&date);
114 failed += test_result ("(Julian Day) Weekday No", wday, 3, 0.1);
115
116 /* Test ln_date_to_zonedate and back */
117
118 ln_date_to_zonedate (&date, &zonedate, 7200);
119 ln_zonedate_to_date (&zonedate, &date);
120
121 JD = ln_get_julian_day (&date);
122
123 failed += test_result ("(Julian Day) ln_date_to_zonedate and ln_zonedate_to_date check - JD for 30/06/1954 00:00:00", JD, 2434923.5, 0.1);
124
125 ln_get_date (JD, &pdate);
126 failed += test_result ("(Julian Day) Day from JD for 30/06/1954 00:00:00", pdate.days, 30, 0.1);
127
128 failed += test_result ("(Julian Day) Month from JD for 30/06/1954 00:00:00", pdate.months, 6, 0.1);
129
130 failed += test_result ("(Julian Day) Year from JD for 30/06/1954 00:00:00", pdate.years, 1954, 0.1);
131
132 failed += test_result ("(Julian Day) Hour from JD for 30/06/1954 00:00:00", pdate.hours, 0, 0.1);
133
134 failed += test_result ("(Julian Day) Minute from JD for 30/06/1954 00:00:00", pdate.minutes, 0, 0.1);
135
136 failed += test_result ("(Julian Day) Second from JD for 30/06/1954 00:00:00", pdate.seconds, 0, 0.001);
137
138 JD = ln_get_julian_from_sys ();
139
140 usleep (10000);
141
142 JD2 = ln_get_julian_from_sys ();
143
144 failed += test_result ("(Julian Day) Diferrence between two sucessive ln_get_julian_from_sys () calls (it shall never be zero)", JD2 - JD, 1e-2/86400.0, .99e-1);
145
146 /* Test that we get from JD same value as is in time_t
147 * struct..was problem before introduction of ln_zonedate (on
148 * machines which doesn't run in UTC).
149 */
150
151 time (&now);
152 ln_get_timet_from_julian (ln_get_julian_from_timet (&now), &now_jd);
153
154 failed += test_result ("(Julian Day) Diferrence between time_t from system and from JD", difftime(now, now_jd), 0, 0);
155
156 return failed;
157 }
158
dynamical_test()159 int dynamical_test ()
160 {
161 struct ln_date date;
162 double TD,JD;
163 int failed = 0;
164
165 /* Dynamical Time test for 01/01/2000 00:00:00 */
166 date.years = 2000;
167 date.months = 1;
168 date.days = 1;
169 date.hours = 0;
170 date.minutes = 0;
171 date.seconds = 0;
172
173 JD = ln_get_julian_day (&date);
174 TD = ln_get_jde (JD);
175 failed += test_result ("(Dynamical Time) TD for 01/01/2000 00:00:00", TD, 2451544.50073877, 0.000001);
176 return failed;
177 }
178
heliocentric_test(void)179 int heliocentric_test (void)
180 {
181 struct ln_equ_posn object;
182 struct ln_date date;
183 double JD;
184
185 double diff;
186
187 int failed = 0;
188
189 object.ra = 0;
190
191 date.years = 2000;
192 date.months = 1;
193 date.days = 1;
194 date.hours = 0;
195 date.minutes = 0;
196 date.seconds = 0;
197
198 JD = ln_get_julian_day (&date);
199
200 object.dec = 60;
201
202 diff = ln_get_heliocentric_time_diff (JD, &object);
203
204 failed += test_result ("(Heliocentric time) TD for 01/01, object on 0h +60", diff, 15.0 * 0.0001, 0.0001);
205
206 object.ra = 270;
207 object.dec = 50;
208
209 diff = ln_get_heliocentric_time_diff (JD, &object);
210
211 failed += test_result ("(Heliocentric time) TD for 01/01, object on 18h +50", diff, -16.0 * 0.0001, 0.0001);
212
213 date.months = 8;
214 date.days = 8;
215
216 JD = ln_get_julian_day (&date);
217
218 diff = ln_get_heliocentric_time_diff (JD, &object);
219
220 failed += test_result ("(Heliocentric time) TD for 08/08, object on 18h +50", diff, 12.0 * 0.0001, 0.0001);
221
222 return failed;
223 }
224
nutation_test(void)225 int nutation_test (void)
226 {
227 double JDE, JD;
228 struct ln_nutation nutation;
229 int failed = 0;
230
231 JD = 2446895.5;
232 JDE = ln_get_jde (JD);
233
234 ln_get_nutation (JD, &nutation);
235 failed += test_result ("(Nutation) longitude (deg) for JD 2446895.5", nutation.longitude, -0.00100561, 0.00000001);
236
237 failed += test_result ("(Nutation) obliquity (deg) for JD 2446895.5", nutation.obliquity, 0.00273297, 0.00000001);
238
239 failed += test_result ("(Nutation) ecliptic (deg) for JD 2446895.5", nutation.ecliptic, 23.44367936, 0.00000001);
240 return failed;
241 }
242
transform_test(void)243 int transform_test(void)
244 {
245 struct lnh_equ_posn hobject, hpollux;
246 struct lnh_lnlat_posn hobserver, hecl;
247 struct ln_equ_posn object, pollux, equ;
248 struct ln_hrz_posn hrz;
249 struct ln_lnlat_posn observer, ecl;
250 struct ln_gal_posn gal;
251 double JD;
252 struct ln_date date;
253 int failed = 0;
254
255 /* observers position */
256 hobserver.lng.neg = 0;
257 hobserver.lng.degrees = 282;
258 hobserver.lng.minutes = 56;
259 hobserver.lng.seconds = 4;
260 hobserver.lat.neg = 0;
261 hobserver.lat.degrees = 38;
262 hobserver.lat.minutes = 55;
263 hobserver.lat.seconds = 17;
264
265 /* object position */
266 hobject.ra.hours = 23;
267 hobject.ra.minutes = 9;
268 hobject.ra.seconds = 16.641;
269 hobject.dec.neg = 1;
270 hobject.dec.degrees = 6;
271 hobject.dec.minutes = 43;
272 hobject.dec.seconds = 11.61;
273
274 /* date and time */
275 date.years = 1987;
276 date.months = 4;
277 date.days = 10;
278 date.hours = 19;
279 date.minutes = 21;
280 date.seconds = 0;
281
282 JD = ln_get_julian_day (&date);
283 ln_hequ_to_equ (&hobject, &object);
284 ln_hlnlat_to_lnlat (&hobserver, &observer);
285
286 ln_get_hrz_from_equ (&object, &observer, JD, &hrz);
287 failed += test_result ("(Transforms) Equ to Horiz ALT ", hrz.alt, 15.12426274, 0.00000001);
288 failed += test_result ("(Transforms) Equ to Horiz AZ ", hrz.az, 68.03429264, 0.00000001);
289
290 /* try something close to the pole */
291 object.dec = 90;
292
293 ln_get_hrz_from_equ (&object, &observer, JD, &hrz);
294 failed += test_result ("(Transforms) Equ to Horiz ALT ", hrz.alt, 38.9213888888, 0.00000001);
295 failed += test_result ("(Transforms) Equ to Horiz AZ ", hrz.az, 180.0, 0.00000001);
296
297 object.dec = -90;
298
299 ln_get_hrz_from_equ (&object, &observer, JD, &hrz);
300 failed += test_result ("(Transforms) Equ to Horiz ALT ", hrz.alt, -38.9213888888, 0.00000001);
301 failed += test_result ("(Transforms) Equ to Horiz AZ ", hrz.az, 0.0, 0.00000001);
302
303 observer.lat *= -1;
304
305 ln_get_hrz_from_equ (&object, &observer, JD, &hrz);
306 failed += test_result ("(Transforms) Equ to Horiz ALT ", hrz.alt, 38.9213888888, 0.00000001);
307 failed += test_result ("(Transforms) Equ to Horiz AZ ", hrz.az, 0.0, 0.00000001);
308
309 object.dec = 90;
310
311 ln_get_hrz_from_equ (&object, &observer, JD, &hrz);
312 failed += test_result ("(Transforms) Equ to Horiz ALT ", hrz.alt, -38.9213888888, 0.00000001);
313 failed += test_result ("(Transforms) Equ to Horiz AZ ", hrz.az, 180.0, 0.00000001);
314
315 /* Equ position of Pollux */
316 hpollux.ra.hours = 7;
317 hpollux.ra.minutes = 45;
318 hpollux.ra.seconds = 18.946;
319 hpollux.dec.neg = 0;
320 hpollux.dec.degrees = 28;
321 hpollux.dec.minutes = 1;
322 hpollux.dec.seconds = 34.26;
323
324 ln_hequ_to_equ (&hpollux, &pollux);
325 ln_get_ecl_from_equ(&pollux, JD, &ecl);
326
327 ln_lnlat_to_hlnlat (&ecl, &hecl);
328 failed += test_result ("(Transforms) Equ to Ecl longitude ", ecl.lng, 113.21542105, 0.00000001);
329 failed += test_result ("(Transforms) Equ to Ecl latitude", ecl.lat, 6.68002727, 0.00000001);
330
331 ln_get_equ_from_ecl(&ecl, JD, &equ);
332 failed += test_result ("(Transforms) Ecl to Equ RA ", equ.ra, 116.32894167, 0.00000001);
333 failed += test_result ("(Transforms) Ecl to Equ DEC", equ.dec, 28.02618333, 0.00000001);
334
335 /* Gal pole */
336 gal.l = 0;
337 gal.b = 90;
338
339 ln_get_equ_from_gal (&gal, &equ);
340 failed += test_result ("(Transforms) Gal to Equ RA", equ.ra, 192.25, 0.00000001);
341 failed += test_result ("(Transforms) Gal to Equ DEC", equ.dec, 27.4, 0.00000001);
342
343 ln_get_gal_from_equ (&equ, &gal);
344 failed += test_result ("(Transforms) Equ to Gal b", gal.b, 90, 0.00000001);
345
346 // Swift triger 174738
347
348 equ.ra = 125.2401;
349 equ.dec = +31.9260;
350
351 ln_get_gal_from_equ2000 (&equ, &gal);
352 failed += test_result ("(Transforms) Equ J2000 to Gal l", gal.l, 190.54, 0.005);
353 failed += test_result ("(Transforms) Equ J2000 to Gal b", gal.b, 31.92, 0.005);
354
355 return failed;
356 }
357
sidereal_test()358 int sidereal_test ()
359 {
360 struct ln_date date;
361 double sd;
362 double JD;
363 int failed = 0;
364
365 /* 10/04/1987 19:21:00 */
366 date.years = 1987;
367 date.months = 4;
368 date.days = 10;
369 date.hours = 19;
370 date.minutes = 21;
371 date.seconds = 0;
372
373 JD = ln_get_julian_day (&date);
374 sd = ln_get_mean_sidereal_time (JD);
375
376 failed += test_result ("(Sidereal) mean hours on 10/04/1987 19:21:00 ", sd, 8.58252488, 0.000001);
377 sd = ln_get_apparent_sidereal_time (JD);
378 failed += test_result ("(Sidereal) apparent hours on 10/04/1987 19:21:00 ", sd, 8.58245327, 0.000001);
379 return failed;
380 }
381
solar_coord_test(void)382 int solar_coord_test (void)
383 {
384 struct ln_helio_posn pos;
385 int failed = 0;
386
387 ln_get_solar_geom_coords (2448908.5, &pos);
388 failed += test_result ("(Solar Coords) longitude (deg) on JD 2448908.5 ", pos.L, 200.00810889, 0.00000001);
389 failed += test_result ("(Solar Coords) latitude (deg) on JD 2448908.5 ", pos.B, 0.00018690, 0.00000001);
390 failed += test_result ("(Solar Coords) radius vector (AU) on JD 2448908.5 ", pos.R, 0.99760852, 0.00000001);
391 return failed;
392 }
393
aberration_test(void)394 int aberration_test (void)
395 {
396 struct lnh_equ_posn hobject;
397 struct ln_equ_posn object, pos;
398 struct ln_date date;
399 double JD;
400 int failed = 0;
401
402 /* object position */
403 hobject.ra.hours = 2;
404 hobject.ra.minutes = 44;
405 hobject.ra.seconds = 12.9747;
406 hobject.dec.neg = 0;
407 hobject.dec.degrees = 49;
408 hobject.dec.minutes = 13;
409 hobject.dec.seconds = 39.896;
410
411 /* date */
412 date.years = 2028;
413 date.months = 11;
414 date.days = 13;
415 date.hours = 4;
416 date.minutes = 31;
417 date.seconds = 0;
418
419 JD = ln_get_julian_day (&date);
420
421 ln_hequ_to_equ (&hobject, &object);
422 ln_get_equ_aber (&object, JD, &pos);
423 failed += test_result ("(Aberration) RA ", pos.ra, 41.06238352, 0.00000001);
424 failed += test_result ("(Aberration) DEC ", pos.dec, 49.22962359, 0.00000001);
425 return failed;
426 }
427
precession_test(void)428 int precession_test(void)
429 {
430 double JD;
431 struct ln_equ_posn object, pos, pos2, pm;
432 struct lnh_equ_posn hobject;
433 struct ln_date grb_date;
434 int failed = 0;
435
436 /* object position */
437 hobject.ra.hours = 2;
438 hobject.ra.minutes = 44;
439 hobject.ra.seconds = 11.986;
440 hobject.dec.neg = 0;
441 hobject.dec.degrees = 49;
442 hobject.dec.minutes = 13;
443 hobject.dec.seconds = 42.48;
444
445 JD = 2462088.69;
446 ln_hequ_to_equ (&hobject, &object);
447
448 pm.ra = 0.03425 * (15.0 / 3600.0);
449 pm.dec = -0.0895 / 3600.0;
450
451 ln_get_equ_pm (&object, &pm, JD, &object);
452
453 failed += test_result ("(Proper motion) RA on JD 2462088.69 ", object.ra, 41.054063, 0.00001);
454 failed += test_result ("(Proper motion) DEC on JD 2462088.69 ", object.dec, 49.227750, 0.00001);
455
456 ln_get_equ_prec (&object, JD, &pos);
457 failed += test_result ("(Precession) RA on JD 2462088.69 ", pos.ra, 41.547214, 0.00003);
458 failed += test_result ("(Precession) DEC on JD 2462088.69 ", pos.dec, 49.348483, 0.00001);
459
460 ln_get_equ_prec2 (&object, JD2000, JD, &pos);
461
462 failed += test_result ("(Precession 2) RA on JD 2462088.69 ", pos.ra, 41.547214, 0.00003);
463 failed += test_result ("(Precession 2) DEC on JD 2462088.69 ", pos.dec, 49.348483, 0.00001);
464
465 ln_get_equ_prec2 (&pos, JD, JD2000, &pos2);
466
467 failed += test_result ("(Precession 2) RA on JD 2451545.0 ", pos2.ra, object.ra, 0.00001);
468 failed += test_result ("(Precession 2) DEC on JD 2451545.0 ", pos2.dec, object.dec, 0.00001);
469
470 // INTEGRAL GRB050922A coordinates lead to RA not in <0-360> range
471 pos.ra = 271.2473;
472 pos.dec = -32.0227;
473
474 grb_date.years = 2005;
475 grb_date.months = 9;
476 grb_date.days = 22;
477 grb_date.hours = 13;
478 grb_date.minutes = 43;
479 grb_date.seconds = 18;
480
481 JD = ln_get_julian_day (&grb_date);
482
483 ln_get_equ_prec2 (&pos, JD, JD2000, &pos2);
484
485 failed += test_result ("(Precession 2) RA on JD 2451545.0 ", pos2.ra, 271.1541, 0.0002);
486 failed += test_result ("(Precession 2) DEC on JD 2451545.0 ", pos2.dec, -32.0235, 0.0002);
487
488 // second test from AA, p. 128
489 hobject.ra.hours = 2;
490 hobject.ra.minutes = 31;
491 hobject.ra.seconds = 48.704;
492 hobject.dec.neg = 0;
493 hobject.dec.degrees = 89;
494 hobject.dec.minutes = 15;
495 hobject.dec.seconds = 50.72;
496
497 ln_hequ_to_equ (&hobject, &object);
498
499 // proper motions
500 pm.ra = ((long double) 0.19877) * (15.0 / 3600.0);
501 pm.dec = ((long double) -0.0152) / 3600.0;
502
503 ln_get_equ_pm (&object, &pm, B1900, &pos);
504
505 ln_get_equ_prec2 (&pos, JD2000, B1900, &pos2);
506
507 // the position is so close to pole, that it depends a lot on how precise
508 // functions we will use. So we get such big errors compared to Meeus.
509 // I checked results agains SLAlib on-line calculator and SLAlib performs even worse then we
510
511 failed += test_result ("(Precision 2) RA on B1900 ", pos2.ra, 20.6412499980, 0.002);
512 failed += test_result ("(Precision 2) DEC on B1900 ", pos2.dec, 88.7739388888, 0.0001);
513
514 ln_get_equ_pm (&object, &pm, JD2050, &pos);
515
516 ln_get_equ_prec2 (&pos, JD2000, JD2050, &pos2);
517
518 failed += test_result ("(Precision 2) RA on J2050 ", pos2.ra, 57.0684583320, 0.003);
519 failed += test_result ("(Precision 2) DEC on J2050 ", pos2.dec, 89.4542722222, 0.0001);
520
521 return failed;
522 }
523
apparent_position_test(void)524 int apparent_position_test(void)
525 {
526 double JD;
527 struct lnh_equ_posn hobject, hpm;
528 struct ln_equ_posn object, pm, pos;
529 int failed = 0;
530
531 /* objects position */
532 hobject.ra.hours = 2;
533 hobject.ra.minutes = 44;
534 hobject.ra.seconds = 12.9747;
535 hobject.dec.neg = 0;
536 hobject.dec.degrees = 49;
537 hobject.dec.minutes = 13;
538 hobject.dec.seconds = 39.896;
539
540 /* proper motion of object */
541 hpm.ra.hours = 0;
542 hpm.ra.minutes = 0;
543 hpm.ra.seconds = 0.03425;
544 hpm.dec.neg = 1;
545 hpm.dec.degrees = 0;
546 hpm.dec.minutes = 0;
547 hpm.dec.seconds = 0.0895;
548
549 JD = 2462088.69;
550 ln_hequ_to_equ (&hobject, &object);
551 ln_hequ_to_equ (&hpm, &pm);
552 ln_get_apparent_posn (&object, &pm, JD, &pos);
553
554 failed += test_result ("(Apparent Position) RA on JD 2462088.69 ", pos.ra, 41.55966517, 0.00000001);
555 failed += test_result ("(Apparent Position) DEC on JD 2462088.69 ", pos.dec, 49.34962340, 0.00000001);
556 return failed;
557 }
558
vsop87_test(void)559 int vsop87_test(void)
560 {
561 struct ln_helio_posn pos;
562 struct lnh_equ_posn hequ;
563 struct ln_equ_posn equ;
564 double JD = 2448976.5;
565 double au;
566 int failed = 0;
567
568 #ifdef DATE
569 struct ln_date date;
570 date.years = 2003;
571 date.months = 1;
572 date.days = 29;
573 date.hours = 0;
574 date.minutes = 0;
575 date.seconds = 0;
576
577 JD = ln_get_julian_day (&date);
578 #endif
579 #ifdef SYS_TIME
580 JD = ln_get_julian_from_sys();
581 #endif
582
583 ln_get_solar_equ_coords (JD, &equ);
584 failed += test_result ("(Solar Position) RA on JD 2448976.5 ", equ.ra, 268.32146893, 0.00000001);
585 failed += test_result ("(Solar Position) DEC on JD 2448976.5 ", equ.dec, -23.43026873, 0.00000001);
586
587 ln_get_mercury_helio_coords(JD, &pos);
588 printf("Mercury L %f B %f R %f\n", pos.L, pos.B, pos.R);
589 ln_get_mercury_equ_coords (JD, &equ);
590 ln_equ_to_hequ (&equ, &hequ);
591 printf("Mercury RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
592 au = ln_get_mercury_earth_dist (JD);
593 printf ("mercury -> Earth dist (AU) %f\n",au);
594 au = ln_get_mercury_solar_dist (JD);
595 printf ("mercury -> Sun dist (AU) %f\n",au);
596 au = ln_get_mercury_disk (JD);
597 printf ("mercury -> illuminated disk %f\n",au);
598 au = ln_get_mercury_magnitude (JD);
599 printf ("mercury -> magnitude %f\n",au);
600 au = ln_get_mercury_phase (JD);
601 printf ("mercury -> phase %f\n",au);
602 au = ln_get_mercury_sdiam (JD);
603 printf ("mercury -> sdiam %f\n",au);
604
605 ln_get_venus_helio_coords(JD, &pos);
606 printf("Venus L %f B %f R %f\n", pos.L, pos.B, pos.R);
607 ln_get_venus_equ_coords (JD, &equ);
608 ln_equ_to_hequ (&equ, &hequ);
609 printf("Venus RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
610 au = ln_get_venus_earth_dist (JD);
611 printf ("venus -> Earth dist (AU) %f\n",au);
612 au = ln_get_venus_solar_dist (JD);
613 printf ("venus -> Sun dist (AU) %f\n",au);
614 au = ln_get_venus_disk (JD);
615 printf ("venus -> illuminated disk %f\n",au);
616 au = ln_get_venus_magnitude (JD);
617 printf ("venus -> magnitude %f\n",au);
618 au = ln_get_venus_phase (JD);
619 printf ("venus -> phase %f\n",au);
620 au = ln_get_venus_sdiam (JD);
621 printf ("venus -> sdiam %f\n",au);
622
623 ln_get_earth_helio_coords(JD, &pos);
624 printf("Earth L %f B %f R %f\n", pos.L, pos.B, pos.R);
625 au = ln_get_earth_solar_dist (JD);
626 printf ("earth -> Sun dist (AU) %f\n",au);
627
628 ln_get_mars_helio_coords(JD, &pos);
629 printf("Mars L %f B %f R %f\n", pos.L, pos.B, pos.R);
630 ln_get_mars_equ_coords (JD, &equ);
631 ln_equ_to_hequ (&equ, &hequ);
632 printf("Mars RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
633 au = ln_get_mars_earth_dist (JD);
634 printf ("mars -> Earth dist (AU) %f\n",au);
635 au = ln_get_mars_solar_dist (JD);
636 printf ("mars -> Sun dist (AU) %f\n",au);
637 au = ln_get_mars_disk (JD);
638 printf ("mars -> illuminated disk %f\n",au);
639 au = ln_get_mars_magnitude (JD);
640 printf ("mars -> magnitude %f\n",au);
641 au = ln_get_mars_phase (JD);
642 printf ("mars -> phase %f\n",au);
643 au = ln_get_mars_sdiam (JD);
644 printf ("mars -> sdiam %f\n",au);
645
646 ln_get_jupiter_helio_coords(JD, &pos);
647 printf("Jupiter L %f B %f R %f\n", pos.L, pos.B, pos.R);
648 ln_get_jupiter_equ_coords (JD, &equ);
649 ln_equ_to_hequ (&equ, &hequ);
650 printf("Jupiter RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
651 au = ln_get_jupiter_earth_dist (JD);
652 printf ("jupiter -> Earth dist (AU) %f\n",au);
653 au = ln_get_jupiter_solar_dist (JD);
654 printf ("jupiter -> Sun dist (AU) %f\n",au);
655 au = ln_get_jupiter_disk (JD);
656 printf ("jupiter -> illuminated disk %f\n",au);
657 au = ln_get_jupiter_magnitude (JD);
658 printf ("jupiter -> magnitude %f\n",au);
659 au = ln_get_jupiter_phase (JD);
660 printf ("jupiter -> phase %f\n",au);
661 au = ln_get_jupiter_pol_sdiam (JD);
662 printf ("jupiter -> polar sdiam %f\n",au);
663 au = ln_get_jupiter_equ_sdiam (JD);
664 printf ("jupiter -> equ sdiam %f\n",au);
665
666 ln_get_saturn_helio_coords(JD, &pos);
667 printf("Saturn L %f B %f R %f\n", pos.L, pos.B, pos.R);
668 ln_get_saturn_equ_coords (JD, &equ);
669 ln_equ_to_hequ (&equ, &hequ);
670 printf("Saturn RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
671 au = ln_get_saturn_earth_dist (JD);
672 printf ("saturn -> Earth dist (AU) %f\n",au);
673 au = ln_get_saturn_solar_dist (JD);
674 printf ("saturn -> Sun dist (AU) %f\n",au);
675 au = ln_get_saturn_disk (JD);
676 printf ("saturn -> illuminated disk %f\n",au);
677 au = ln_get_saturn_magnitude (JD);
678 printf ("saturn -> magnitude %f\n",au);
679 au = ln_get_saturn_phase (JD);
680 printf ("saturn -> phase %f\n",au);
681 au = ln_get_saturn_pol_sdiam (JD);
682 printf ("saturn -> polar sdiam %f\n",au);
683 au = ln_get_saturn_equ_sdiam (JD);
684 printf ("saturn -> equ sdiam %f\n",au);
685
686 ln_get_uranus_helio_coords(JD, &pos);
687 printf("Uranus L %f B %f R %f\n", pos.L, pos.B, pos.R);
688 ln_get_uranus_equ_coords (JD, &equ);
689 ln_equ_to_hequ (&equ, &hequ);
690 printf("Uranus RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
691 au = ln_get_uranus_earth_dist (JD);
692 printf ("uranus -> Earth dist (AU) %f\n",au);
693 au = ln_get_uranus_solar_dist (JD);
694 printf ("uranus -> Sun dist (AU) %f\n",au);
695 au = ln_get_uranus_disk (JD);
696 printf ("uranus -> illuminated disk %f\n",au);
697 au = ln_get_uranus_magnitude (JD);
698 printf ("uranus -> magnitude %f\n",au);
699 au = ln_get_uranus_phase (JD);
700 printf ("uranus -> phase %f\n",au);
701 au = ln_get_uranus_sdiam (JD);
702 printf ("uranus -> sdiam %f\n",au);
703
704 ln_get_neptune_helio_coords(JD, &pos);
705 printf("Neptune L %f B %f R %f\n", pos.L, pos.B, pos.R);
706 ln_get_neptune_equ_coords (JD, &equ);
707 ln_equ_to_hequ (&equ, &hequ);
708 printf("Neptune RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
709 au = ln_get_neptune_earth_dist (JD);
710 printf ("neptune -> Earth dist (AU) %f\n",au);
711 au = ln_get_neptune_solar_dist (JD);
712 printf ("neptune -> Sun dist (AU) %f\n",au);
713 au = ln_get_neptune_disk (JD);
714 printf ("neptune -> illuminated disk %f\n",au);
715 au = ln_get_neptune_magnitude (JD);
716 printf ("neptune -> magnitude %f\n",au);
717 au = ln_get_neptune_phase (JD);
718 printf ("neptune -> phase %f\n",au);
719 au = ln_get_neptune_sdiam (JD);
720 printf ("neptune -> sdiam %f\n",au);
721
722 ln_get_pluto_helio_coords(JD, &pos);
723 printf("Pluto L %f B %f R %f\n", pos.L, pos.B, pos.R);
724 ln_get_pluto_equ_coords (JD, &equ);
725 ln_equ_to_hequ (&equ, &hequ);
726 printf("Pluto RA %d:%d:%f Dec %d:%d:%f\n", hequ.ra.hours, hequ.ra.minutes, hequ.ra.seconds, hequ.dec.degrees, hequ.dec.minutes, hequ.dec.seconds);
727 au = ln_get_pluto_earth_dist (JD);
728 printf ("pluto -> Earth dist (AU) %f\n",au);
729 au = ln_get_pluto_solar_dist (JD);
730 printf ("pluto -> Sun dist (AU) %f\n",au);
731 au = ln_get_pluto_disk (JD);
732 printf ("pluto -> illuminated disk %f\n",au);
733 au = ln_get_pluto_magnitude (JD);
734 printf ("pluto -> magnitude %f\n",au);
735 au = ln_get_pluto_phase (JD);
736 printf ("pluto -> phase %f\n",au);
737 au = ln_get_pluto_sdiam (JD);
738 printf ("pluto -> sdiam %f\n",au);
739
740 return failed;
741 }
742
lunar_test()743 int lunar_test ()
744 {
745 double JD = 2448724.5;
746
747 struct ln_rect_posn moon;
748 struct ln_equ_posn equ;
749 struct ln_lnlat_posn ecl;
750 int failed = 0;
751
752 /* JD = get_julian_from_sys();*/
753 /*JD=2448724.5;*/
754 ln_get_lunar_geo_posn (JD, &moon, 0);
755 printf ("lunar x %f y %f z %f\n",moon.X, moon.Y, moon.Z);
756 ln_get_lunar_ecl_coords (JD, &ecl, 0);
757 printf ("lunar long %f lat %f\n",ecl.lng, ecl.lat);
758 ln_get_lunar_equ_coords_prec (JD, &equ, 0);
759 printf ("lunar RA %f Dec %f\n",equ.ra, equ.dec);
760 printf ("lunar distance km %f\n", ln_get_lunar_earth_dist(JD));
761 printf ("lunar disk %f\n", ln_get_lunar_disk(JD));
762 printf ("lunar phase %f\n", ln_get_lunar_phase(JD));
763 printf ("lunar bright limb %f\n", ln_get_lunar_bright_limb(JD));
764 return failed;
765 }
766
elliptic_motion_test()767 int elliptic_motion_test ()
768 {
769 double r,v,l,V,dist;
770 double E, e_JD, o_JD;
771 struct ln_ell_orbit orbit;
772 struct ln_rect_posn posn;
773 struct ln_date epoch_date, obs_date;
774 struct ln_equ_posn equ_posn;
775 int failed = 0;
776
777 obs_date.years = 1990;
778 obs_date.months = 10;
779 obs_date.days = 6;
780 obs_date.hours = 0;
781 obs_date.minutes = 0;
782 obs_date.seconds = 0;
783
784 epoch_date.years = 1990;
785 epoch_date.months = 10;
786 epoch_date.days = 28;
787 epoch_date.hours = 12;
788 epoch_date.minutes = 30;
789 epoch_date.seconds = 0;
790
791 e_JD = ln_get_julian_day (&epoch_date);
792 o_JD = ln_get_julian_day (&obs_date);
793
794 orbit.JD = e_JD;
795 orbit.a = 2.2091404;
796 orbit.e = 0.8502196;
797 orbit.i = 11.94525;
798 orbit.omega = 334.75006;
799 orbit.w = 186.23352;
800 orbit.n = 0;
801
802 E = ln_solve_kepler (0.1, 5.0);
803 failed += test_result ("(Equation of kepler) E when e is 0.1 and M is 5.0 ", E, 5.554589253872320, 0.000000000001);
804
805 v = ln_get_ell_true_anomaly (0.1, E);
806 failed += test_result ("(True Anomaly) v when e is 0.1 and E is 5.5545 ", v, 6.13976152, 0.00000001);
807
808 r = ln_get_ell_radius_vector (0.5, 0.1, E);
809 failed += test_result ("(Radius Vector) r when v is , e is 0.1 and E is 5.5545 ", r, 0.45023478, 0.00000001);
810
811 ln_get_ell_geo_rect_posn (&orbit, o_JD, &posn);
812 failed += test_result ("(Geocentric Rect Coords X) for comet Enckle ", posn.X, 0.72549850, 0.00000001);
813 failed += test_result ("(Geocentric Rect Coords Y) for comet Enckle ", posn.Y, -0.28443537, 0.00000001);
814 failed += test_result ("(Geocentric Rect Coords Z) for comet Enckle ", posn.Z, -0.27031656, 0.00000001);
815
816 ln_get_ell_helio_rect_posn (&orbit, o_JD, &posn);
817 failed += test_result ("(Heliocentric Rect Coords X) for comet Enckle ", posn.X, 0.25017473, 0.00000001);
818 failed += test_result ("(Heliocentric Rect Coords Y) for comet Enckle ", posn.Y, 0.48476422, 0.00000001);
819 failed += test_result ("(Heliocentric Rect Coords Z) for comet Enckle ", posn.Z, 0.35716517, 0.00000001);
820
821 ln_get_ell_body_equ_coords (o_JD, &orbit, &equ_posn);
822 failed += test_result ("(RA) for comet Enckle ", equ_posn.ra, 158.58242653, 0.00000001);
823 failed += test_result ("(Dec) for comet Enckle ", equ_posn.dec, 19.13924815, 0.00000001);
824
825 l = ln_get_ell_orbit_len (&orbit);
826 failed += test_result ("(Orbit Length) for comet Enckle in AU ", l, 10.75710334, 0.00000001);
827
828 V = ln_get_ell_orbit_pvel (&orbit);
829 failed += test_result ("(Orbit Perihelion Vel) for comet Enckle in kms ", V, 70.43130198, 0.00000001);
830
831 V = ln_get_ell_orbit_avel (&orbit);
832 failed += test_result ("(Orbit Aphelion Vel) for comet Enckle in kms ", V, 5.70160892, 0.00000001);
833
834 V = ln_get_ell_orbit_vel (o_JD, &orbit);
835 failed += test_result ("(Orbit Vel JD) for comet Enckle in kms ", V, 48.16148331, 0.00000001);
836
837 dist = ln_get_ell_body_solar_dist (o_JD, &orbit);
838 failed += test_result ("(Body Solar Dist) for comet Enckle in AU ", dist, 0.65203581, 0.00000001);
839
840 dist = ln_get_ell_body_earth_dist (o_JD, &orbit);
841 failed += test_result ("(Body Earth Dist) for comet Enckle in AU ", dist, 0.82481670, 0.00000001);
842
843 // TNO http://www.cfa.harvard.edu/mpec/K05/K05O42.html
844
845 obs_date.years = 2006;
846 obs_date.months = 5;
847 obs_date.days = 5;
848 obs_date.hours = 0;
849 obs_date.minutes = 0;
850 obs_date.seconds = 0;
851
852 epoch_date.years = 2006;
853 epoch_date.months = 3;
854 epoch_date.days = 6;
855 epoch_date.hours = 0;
856 epoch_date.minutes = 0;
857 epoch_date.seconds = 0;
858
859 e_JD = ln_get_julian_day (&epoch_date);
860 o_JD = ln_get_julian_day (&obs_date);
861
862 orbit.JD = e_JD;
863 orbit.a = 45.7082927;
864 orbit.e = 0.1550125;
865 orbit.i = 28.99870;
866 orbit.omega = 79.55499;
867 orbit.w = 296.40937;
868 orbit.n = 0.00318942;
869
870 // MPO refers to Mean anomaly & epoch, we hence need to convert epoch
871 // to perihelion pass
872
873 orbit.JD -= 147.09926 / orbit.n;
874
875 ln_get_ell_body_equ_coords (o_JD, &orbit, &equ_posn);
876 failed += test_result ("(RA) for TNO K05F09Y ", equ_posn.ra, 184.3699999995, 0.001);
877 failed += test_result ("(Dec) for TNO K05F09Y ", equ_posn.dec, 30.3316666666, 0.001);
878
879 return failed;
880 }
881
882 /* need a proper parabolic orbit to properly test */
parabolic_motion_test()883 int parabolic_motion_test ()
884 {
885 double r,v,dist;
886 double e_JD, o_JD;
887 struct ln_par_orbit orbit;
888 struct ln_rect_posn posn;
889 struct ln_date epoch_date, obs_date;
890 struct ln_equ_posn equ_posn;
891 int failed = 0;
892
893 obs_date.years = 2003;
894 obs_date.months = 1;
895 obs_date.days = 11;
896 obs_date.hours = 0;
897 obs_date.minutes = 0;
898 obs_date.seconds = 0;
899
900 epoch_date.years = 2003;
901 epoch_date.months = 1;
902 epoch_date.days = 29;
903 epoch_date.hours = 0;
904 epoch_date.minutes = 6;
905 epoch_date.seconds = 37.44;
906
907 e_JD = ln_get_julian_day (&epoch_date);
908 o_JD = ln_get_julian_day (&obs_date);
909
910 orbit.q = 0.190082;
911 orbit.i = 94.1511;
912 orbit.w = 187.5613;
913 orbit.omega = 119.0676;
914 orbit.JD = e_JD;
915
916 v = ln_get_par_true_anomaly (orbit.q, o_JD - e_JD);
917 failed += test_result ("(True Anomaly) v when e is 0.1 and E is 5.5545 ", v, 247.18968605, 0.00000001);
918
919 r = ln_get_par_radius_vector (orbit.q, o_JD - e_JD);
920 failed += test_result ("(Radius Vector) r when v is , e is 0.1 and E is 5.5545 ", r, 0.62085992, 0.00000001);
921
922 ln_get_par_geo_rect_posn (&orbit, o_JD, &posn);
923 failed += test_result ("(Geocentric Rect Coords X) for comet C/2002 X5 (Kudo-Fujikawa) ", posn.X, 0.29972461, 0.00000001);
924 failed += test_result ("(Geocentric Rect Coords Y) for comet C/2002 X5 (Kudo-Fujikawa) ", posn.Y, -0.93359772, 0.00000001);
925 failed += test_result ("(Geocentric Rect Coords Z) for comet C/2002 X5 (Kudo-Fujikawa) ", posn.Z, 0.24639194, 0.00000001);
926
927 ln_get_par_helio_rect_posn (&orbit, o_JD, &posn);
928 failed += test_result ("(Heliocentric Rect Coords X) for comet C/2002 X5 (Kudo-Fujikawa) ", posn.X, -0.04143700, 0.00000001);
929 failed += test_result ("(Heliocentric Rect Coords Y) for comet C/2002 X5 (Kudo-Fujikawa) ", posn.Y, -0.08736588, 0.00000001);
930 failed += test_result ("(Heliocentric Rect Coords Z) for comet C/2002 X5 (Kudo-Fujikawa) ", posn.Z, 0.61328397, 0.00000001);
931
932 ln_get_par_body_equ_coords (o_JD, &orbit, &equ_posn);
933 failed += test_result ("(RA) for comet C/2002 X5 (Kudo-Fujikawa) ", equ_posn.ra, 287.79617309, 0.00000001);
934 failed += test_result ("(Dec) for comet C/2002 X5 (Kudo-Fujikawa) ", equ_posn.dec, 14.11800859, 0.00000001);
935
936 dist = ln_get_par_body_solar_dist (o_JD, &orbit);
937 failed += test_result ("(Body Solar Dist) for comet C/2002 X5 (Kudo-Fujikawa) in AU ", dist, 0.62085992, 0.00001);
938
939 dist = ln_get_par_body_earth_dist (o_JD, &orbit);
940 failed += test_result ("(Body Earth Dist) for comet C/2002 X5 (Kudo-Fujikawa) in AU ", dist, 1.01101362, 0.00001);
941 return failed;
942 }
943
944 /* data from Meeus, chapter 35 */
hyperbolic_motion_test()945 int hyperbolic_motion_test ()
946 {
947 double r,v,dist;
948 double e_JD, o_JD;
949 struct ln_hyp_orbit orbit;
950 struct ln_date epoch_date, obs_date;
951 struct ln_equ_posn equ_posn;
952 int failed = 0;
953
954 orbit.q = 3.363943;
955 orbit.e = 1.05731;
956
957 // the one from Meeus..
958 v = ln_get_hyp_true_anomaly (orbit.q, orbit.e, 1237.1);
959 failed += test_result ("(True Anomaly) v when q is 3.363943 and e is 1.05731 ", v, 109.40598, 0.00001);
960
961 r = ln_get_hyp_radius_vector (orbit.q, orbit.e, 1237.1);
962 failed += test_result ("(Radius Vector) r when q is 3.363943 and e is 1.05731 ", r, 10.668551, 0.00001);
963
964 // and now something real.. C/2001 Q4 (NEAT)
965 obs_date.years = 2004;
966 obs_date.months = 5;
967 obs_date.days = 15;
968 obs_date.hours = 0;
969 obs_date.minutes = 0;
970 obs_date.seconds = 0;
971
972 epoch_date.years = 2004;
973 epoch_date.months = 5;
974 epoch_date.days = 15;
975 epoch_date.hours = 23;
976 epoch_date.minutes = 12;
977 epoch_date.seconds = 37.44;
978
979 e_JD = ln_get_julian_day (&epoch_date);
980 o_JD = ln_get_julian_day (&obs_date);
981
982 orbit.q = 0.961957;
983 orbit.e = 1.000744;
984 orbit.i = 99.6426;
985 orbit.w = 1.2065;
986 orbit.omega = 210.2785;
987 orbit.JD = e_JD;
988
989 r = ln_get_hyp_radius_vector (orbit.q, orbit.e, o_JD - e_JD);
990 failed += test_result ("(Radius Vector) r for C/2001 Q4 (NEAT) ", r, 0.962, 0.001);
991
992 ln_get_hyp_body_equ_coords (o_JD, &orbit, &equ_posn);
993 failed += test_result ("(RA) for comet C/2001 Q4 (NEAT) ", equ_posn.ra, 128.01, 0.01);
994 failed += test_result ("(Dec) for comet C/2001 Q4 (NEAT) ", equ_posn.dec, 18.3266666666, 0.03);
995
996 dist = ln_get_hyp_body_solar_dist (o_JD, &orbit);
997 failed += test_result ("(Body Solar Dist) for comet C/2001 Q4 (NEAT) in AU ", dist, 0.962, 0.001);
998
999 obs_date.years = 2005;
1000 obs_date.months = 1;
1001
1002 o_JD = ln_get_julian_day (&obs_date);
1003
1004 r = ln_get_hyp_radius_vector (orbit.q, orbit.e, o_JD - e_JD);
1005 failed += test_result ("(Radius Vector) r for C/2001 Q4 (NEAT) ", r, 3.581, 0.001);
1006
1007 ln_get_hyp_body_equ_coords (o_JD, &orbit, &equ_posn);
1008 failed += test_result ("(RA) for comet C/2001 Q4 (NEAT) ", equ_posn.ra, 332.9025, 0.01);
1009 failed += test_result ("(Dec) for comet C/2001 Q4 (NEAT) ", equ_posn.dec, 58.6116666666, 0.001);
1010
1011 return failed;
1012 }
1013
1014
rst_test()1015 int rst_test ()
1016 {
1017 struct ln_lnlat_posn observer;
1018 struct ln_rst_time rst;
1019 struct ln_hms hms;
1020 struct ln_dms dms;
1021 struct ln_date date;
1022 struct ln_equ_posn object;
1023 double JD, JD_next;
1024 int ret;
1025 int failed = 0;
1026
1027 // Arcturus
1028 hms.hours = 14;
1029 hms.minutes = 15;
1030 hms.seconds = 39.67;
1031
1032 dms.neg = 0;
1033 dms.degrees = 19;
1034 dms.minutes = 10;
1035 dms.seconds = 56.7;
1036
1037 object.ra = ln_hms_to_deg (&hms);
1038 object.dec = ln_dms_to_deg (&dms);
1039
1040 date.years = 2006;
1041 date.months = 1;
1042 date.days = 17;
1043
1044 date.hours = 0;
1045 date.minutes = 0;
1046 date.seconds = 0;
1047
1048 JD = ln_get_julian_day (&date);
1049
1050 observer.lng = 15;
1051 observer.lat = 51;
1052
1053 ret = ln_get_object_rst (JD, &observer, &object, &rst);
1054 failed += test_result ("Arcturus sometimes rise at 15 E, 51 N", ret, 0, 0);
1055
1056 if (!ret)
1057 {
1058 ln_get_date (rst.rise, &date);
1059 failed += test_result ("Arcturus rise hour on 2006/01/17 at (15 E,51 N)", date.hours, 21, 0);
1060 failed += test_result ("Arcturus rise minute on 2006/01/17 at (15 E,51 N)", date.minutes, 40, 0);
1061
1062 ln_get_date (rst.transit, &date);
1063 failed += test_result ("Arcturus transit hour on 2006/01/17 at (15 E,51 N)", date.hours, 5, 0);
1064 failed += test_result ("Arcturus transit minute on 2006/01/17 at (15 E,51 N)", date.minutes, 29, 0);
1065
1066 ln_get_date (rst.set, &date);
1067 failed += test_result ("Arcturus set hour on 2006/01/17 at (15 E,51 N)", date.hours, 13, 0);
1068 failed += test_result ("Arcturus set minute on 2006/01/17 at (15 E,51 N)", date.minutes, 14, 0);
1069 }
1070
1071 JD_next = rst.transit - 0.001;
1072 ret = ln_get_object_next_rst (JD_next, &observer, &object, &rst);
1073 failed += test_result ("Arcturus sometimes rise at 15 E, 51 N", ret, 0, 0);
1074
1075 if (!ret)
1076 {
1077 failed += test_result ("Arcturus next date is less then transit time", (JD_next < rst.transit), 1, 0);
1078 failed += test_result ("Arcturus next transit time is less then set time", (rst.transit < rst.set), 1, 0);
1079 failed += test_result ("Arcturus next set time is less then rise time", (rst.set < rst.rise), 1, 0);
1080
1081 ln_get_date (rst.rise, &date);
1082 failed += test_result ("Arcturus next rise hour on 2006/01/17 at (15 E,51 N)", date.hours, 21, 0);
1083 failed += test_result ("Arcturus next rise minute on 2006/01/17 at (15 E,51 N)", date.minutes, 40, 0);
1084
1085 ln_get_date (rst.transit, &date);
1086 failed += test_result ("Arcturus next transit hour on 2006/01/17 at (15 E,51 N)", date.hours, 5, 0);
1087 failed += test_result ("Arcturus next transit minute on 2006/01/17 at (15 E,51 N)", date.minutes, 29, 0);
1088
1089 ln_get_date (rst.set, &date);
1090 failed += test_result ("Arcturus next set hour on 2006/01/17 at (15 E,51 N)", date.hours, 13, 0);
1091 failed += test_result ("Arcturus next set minute on 2006/01/17 at (15 E,51 N)", date.minutes, 14, 0);
1092 }
1093
1094 JD_next = rst.set - 0.001;
1095 ret = ln_get_object_next_rst (JD_next, &observer, &object, &rst);
1096 failed += test_result ("Arcturus sometimes rise at 15 E, 51 N", ret, 0, 0);
1097
1098 if (!ret)
1099 {
1100 failed += test_result ("Arcturus next date is less then set time", (JD_next < rst.set), 1, 0);
1101 failed += test_result ("Arcturus next set time is less then rise time", (rst.set < rst.rise), 1, 0);
1102 failed += test_result ("Arcturus next rise time is less then transit time", (rst.rise < rst.transit), 1, 0);
1103
1104 ln_get_date (rst.rise, &date);
1105 failed += test_result ("Arcturus next rise hour on 2006/01/17 at (15 E,51 N)", date.hours, 21, 0);
1106 failed += test_result ("Arcturus next rise minute on 2006/01/17 at (15 E,51 N)", date.minutes, 40, 0);
1107
1108 ln_get_date (rst.transit, &date);
1109 failed += test_result ("Arcturus next transit hour on 2006/01/18 at (15 E,51 N)", date.hours, 5, 0);
1110 failed += test_result ("Arcturus next transit minute on 2006/01/18 at (15 E,51 N)", date.minutes, 25, 0);
1111
1112 ln_get_date (rst.set, &date);
1113 failed += test_result ("Arcturus next set hour on 2006/01/17 at (15 E,51 N)", date.hours, 13, 0);
1114 failed += test_result ("Arcturus next set minute on 2006/01/17 at (15 E,51 N)", date.minutes, 14, 0);
1115 }
1116
1117 JD_next = rst.rise - 0.001;
1118 ret = ln_get_object_next_rst (JD_next, &observer, &object, &rst);
1119 failed += test_result ("Arcturus sometimes rise at 15 E, 51 N", ret, 0, 0);
1120
1121 if (!ret)
1122 {
1123 failed += test_result ("Arcturus next date is less then rise time", (JD_next < rst.rise), 1, 0);
1124 failed += test_result ("Arcturus next rise time is less then transit time", (rst.rise < rst.transit), 1, 0);
1125 failed += test_result ("Arcturus next transit time is less then set time", (rst.transit < rst.set), 1, 0);
1126
1127 ln_get_date (rst.rise, &date);
1128 failed += test_result ("Arcturus next rise hour on 2006/01/17 at (15 E,51 N)", date.hours, 21, 0);
1129 failed += test_result ("Arcturus next rise minute on 2006/01/17 at (15 E,51 N)", date.minutes, 40, 0);
1130
1131 ln_get_date (rst.transit, &date);
1132 failed += test_result ("Arcturus next transit hour on 2006/01/18 at (15 E,51 N)", date.hours, 5, 0);
1133 failed += test_result ("Arcturus next transit minute on 2006/01/18 at (15 E,51 N)", date.minutes, 25, 0);
1134
1135 ln_get_date (rst.set, &date);
1136 failed += test_result ("Arcturus next set hour on 2006/01/18 at (15 E,51 N)", date.hours, 13, 0);
1137 failed += test_result ("Arcturus next set minute on 2006/01/18 at (15 E,51 N)", date.minutes, 10, 0);
1138 }
1139
1140 JD_next = rst.rise + 0.001;
1141 ret = ln_get_object_next_rst (JD_next, &observer, &object, &rst);
1142 failed += test_result ("Arcturus sometimes rise at 15 E, 51 N", ret, 0, 0);
1143
1144 if (!ret)
1145 {
1146 failed += test_result ("Arcturus next date is less then transit time", (JD_next < rst.transit), 1, 0);
1147 failed += test_result ("Arcturus next transit time is less then set time", (rst.transit < rst.set), 1, 0);
1148 failed += test_result ("Arcturus next set time is less then rise time", (rst.set < rst.rise), 1, 0);
1149
1150 ln_get_date (rst.rise, &date);
1151 failed += test_result ("Arcturus next rise hour on 2006/01/18 at (15 E,51 N)", date.hours, 21, 0);
1152 failed += test_result ("Arcturus next rise minute on 2006/01/18 at (15 E,51 N)", date.minutes, 37, 0);
1153
1154 ln_get_date (rst.transit, &date);
1155 failed += test_result ("Arcturus next transit hour on 2006/01/18 at (15 E,51 N)", date.hours, 5, 0);
1156 failed += test_result ("Arcturus next transit minute on 2006/01/18 at (15 E,51 N)", date.minutes, 25, 0);
1157
1158 ln_get_date (rst.set, &date);
1159 failed += test_result ("Arcturus next set hour on 2006/01/18 at (15 E,51 N)", date.hours, 13, 0);
1160 failed += test_result ("Arcturus next set minute on 2006/01/18 at (15 E,51 N)", date.minutes, 10, 0);
1161 }
1162
1163 ret = ln_get_object_rst_horizon (JD, &observer, &object, 20, &rst);
1164 failed += test_result ("Arcturus sometimes rise above 20 deg at 15 E, 51 N", ret, 0, 0);
1165
1166 if (!ret)
1167 {
1168 ln_get_date (rst.rise, &date);
1169 failed += test_result ("Arcturus rise above 20 deg hour on 2006/01/17 at (15 E,51 N)", date.hours, 0, 0);
1170 failed += test_result ("Arcturus rise above 20 deg minute on 2006/01/17 at (15 E,51 N)", date.minutes, 6, 0);
1171
1172 ln_get_date (rst.transit, &date);
1173 failed += test_result ("Arcturus transit hour on 2006/01/17 at (15 E,51 N)", date.hours, 5, 0);
1174 failed += test_result ("Arcturus transit minute on 2006/01/17 at (15 E,51 N)", date.minutes, 29, 0);
1175
1176 ln_get_date (rst.set, &date);
1177 failed += test_result ("Arcturus set bellow 20 deg hour on 2006/01/17 at (15 E,51 N)", date.hours, 10, 0);
1178 failed += test_result ("Arcturus set bellow 20 deg minute on 2006/01/17 at (15 E,51 N)", date.minutes, 52, 0);
1179 }
1180
1181 JD_next = rst.rise - 0.001;
1182 ret = ln_get_object_next_rst_horizon (JD_next, &observer, &object, 20, &rst);
1183 failed += test_result ("Arcturus sometimes rise above 20 deg at 15 E, 51 N", ret, 0, 0);
1184
1185 if (!ret)
1186 {
1187 failed += test_result ("Arcturus next date is less then rise time", (JD_next < rst.rise), 1, 0);
1188 failed += test_result ("Arcturus next rise time is less then transit time", (rst.rise < rst.transit), 1, 0);
1189 failed += test_result ("Arcturus next transit time is less then set time", (rst.transit < rst.set), 1, 0);
1190
1191 ln_get_date (rst.rise, &date);
1192 failed += test_result ("Arcturus next rise above 20 deg hour on 2006/01/17 at (15 E,51 N)", date.hours, 0, 0);
1193 failed += test_result ("Arcturus next rise above 20 deg minute on 2006/01/17 at (15 E,51 N)", date.minutes, 6, 0);
1194
1195 ln_get_date (rst.transit, &date);
1196 failed += test_result ("Arcturus next transit hour on 2006/01/17 at (15 E,51 N)", date.hours, 5, 0);
1197 failed += test_result ("Arcturus next transit minute on 2006/01/17 at (15 E,51 N)", date.minutes, 29, 0);
1198
1199 ln_get_date (rst.set, &date);
1200 failed += test_result ("Arcturus next set bellow 20 deg hour on 2006/01/17 at (15 E,51 N)", date.hours, 10, 0);
1201 failed += test_result ("Arcturus next set bellow 20 deg minute on 2006/01/17 at (15 E,51 N)", date.minutes, 52, 0);
1202 }
1203
1204 JD_next = rst.transit - 0.001;
1205 ret = ln_get_object_next_rst_horizon (JD_next, &observer, &object, 20, &rst);
1206 failed += test_result ("Arcturus sometimes rise above 20 deg at 15 E, 51 N", ret, 0, 0);
1207
1208 if (!ret)
1209 {
1210 failed += test_result ("Arcturus next date is less then transit time", (JD_next < rst.transit), 1, 0);
1211 failed += test_result ("Arcturus next transit time is less then set time", (rst.transit < rst.set), 1, 0);
1212 failed += test_result ("Arcturus next set time is less then rise time", (rst.set < rst.rise), 1, 0);
1213
1214 ln_get_date (rst.rise, &date);
1215 failed += test_result ("Arcturus next rise above 20 deg hour on 2006/01/18 at (15 E,51 N)", date.hours, 0, 0);
1216 failed += test_result ("Arcturus next rise above 20 deg minute on 2006/01/18 at (15 E,51 N)", date.minutes, 2, 0);
1217
1218 ln_get_date (rst.transit, &date);
1219 failed += test_result ("Arcturus next transit hour on 2006/01/17 at (15 E,51 N)", date.hours, 5, 0);
1220 failed += test_result ("Arcturus next transit minute on 2006/01/17 at (15 E,51 N)", date.minutes, 29, 0);
1221
1222 ln_get_date (rst.set, &date);
1223 failed += test_result ("Arcturus next set bellow 20 deg hour on 2006/01/17 at (15 E,51 N)", date.hours, 10, 0);
1224 failed += test_result ("Arcturus next set bellow 20 deg minute on 2006/01/17 at (15 E,51 N)", date.minutes, 52, 0);
1225 }
1226
1227 JD_next = rst.set - 0.001;
1228 ret = ln_get_object_next_rst_horizon (JD_next, &observer, &object, 20, &rst);
1229 failed += test_result ("Arcturus sometimes rise above 20 deg at 15 E, 51 N", ret, 0, 0);
1230
1231 if (!ret)
1232 {
1233 failed += test_result ("Arcturus next date is less then set time", (JD_next < rst.set), 1, 0);
1234 failed += test_result ("Arcturus next set time is less then rise time", (rst.set < rst.rise), 1, 0);
1235 failed += test_result ("Arcturus next rise time is less then transit time", (rst.rise < rst.transit), 1, 0);
1236
1237 ln_get_date (rst.rise, &date);
1238 failed += test_result ("Arcturus next rise above 20 deg hour on 2006/01/18 at (15 E,51 N)", date.hours, 0, 0);
1239 failed += test_result ("Arcturus next rise above 20 deg minute on 2006/01/18 at (15 E,51 N)", date.minutes, 2, 0);
1240
1241 ln_get_date (rst.transit, &date);
1242 failed += test_result ("Arcturus next transit hour on 2006/01/18 at (15 E,51 N)", date.hours, 5, 0);
1243 failed += test_result ("Arcturus next transit minute on 2006/01/18 at (15 E,51 N)", date.minutes, 25, 0);
1244
1245 ln_get_date (rst.set, &date);
1246 failed += test_result ("Arcturus next set bellow 20 deg hour on 2006/01/17 at (15 E,51 N)", date.hours, 10, 0);
1247 failed += test_result ("Arcturus next set bellow 20 deg minute on 2006/01/17 at (15 E,51 N)", date.minutes, 52, 0);
1248 }
1249
1250 JD_next = rst.set + 0.001;
1251 ret = ln_get_object_next_rst_horizon (JD_next, &observer, &object, 20, &rst);
1252 failed += test_result ("Arcturus sometimes rise above 20 deg at 15 E, 51 N", ret, 0, 0);
1253
1254 if (!ret)
1255 {
1256 failed += test_result ("Arcturus next date is less then rise time", (JD_next < rst.rise), 1, 0);
1257 failed += test_result ("Arcturus next rise time is less then transit time", (rst.rise < rst.transit), 1, 0);
1258 failed += test_result ("Arcturus next transit time is less then set time", (rst.transit < rst.set), 1, 0);
1259
1260 ln_get_date (rst.rise, &date);
1261 failed += test_result ("Arcturus next rise above 20 deg hour on 2006/01/18 at (15 E,51 N)", date.hours, 0, 0);
1262 failed += test_result ("Arcturus next rise above 20 deg minute on 2006/01/18 at (15 E,51 N)", date.minutes, 2, 0);
1263
1264 ln_get_date (rst.transit, &date);
1265 failed += test_result ("Arcturus next transit hour on 2006/01/18 at (15 E,51 N)", date.hours, 5, 0);
1266 failed += test_result ("Arcturus next transit minute on 2006/01/18 at (15 E,51 N)", date.minutes, 25, 0);
1267
1268 ln_get_date (rst.set, &date);
1269 failed += test_result ("Arcturus next set bellow 20 deg hour on 2006/01/18 at (15 E,51 N)", date.hours, 10, 0);
1270 failed += test_result ("Arcturus next set bellow 20 deg minute on 2006/01/18 at (15 E,51 N)", date.minutes, 48, 0);
1271 }
1272
1273 observer.lat = -51;
1274
1275 ret = ln_get_object_rst (JD, &observer, &object, &rst);
1276 failed += test_result ("Arcturus sometimes rise at 15 E, 51 S", ret, 0, 0);
1277
1278 if (!ret)
1279 {
1280 ln_get_date (rst.rise, &date);
1281 failed += test_result ("Arcturus rise hour on 2006/01/17 at (15 E,51 S)", date.hours, 1, 0);
1282 failed += test_result ("Arcturus rise minute on 2006/01/17 at (15 E,51 S)", date.minutes, 7, 0);
1283
1284 ln_get_date (rst.transit, &date);
1285 failed += test_result ("Arcturus transit hour on 2006/01/17 at (15 E,51 S)", date.hours, 5, 0);
1286 failed += test_result ("Arcturus transit minute on 2006/01/17 at (15 E,51 S)", date.minutes, 29, 0);
1287
1288 ln_get_date (rst.set, &date);
1289 failed += test_result ("Arcturus set hour on 2006/01/17 at (15 E,51 S)", date.hours, 9, 0);
1290 failed += test_result ("Arcturus set minute on 2006/01/17 at (15 E,51 S)", date.minutes, 51, 0);
1291 }
1292
1293 ret = ln_get_object_rst_horizon (JD, &observer, &object, -20, &rst);
1294 failed += test_result ("Arcturus sometimes rise above -20 deg at 15 E, 51 S", ret, 0, 0);
1295
1296 if (!ret)
1297 {
1298 ln_get_date (rst.rise, &date);
1299 failed += test_result ("Arcturus rise above -20 deg hour on 2006/01/17 at (15 E,51 S)", date.hours, 22, 0);
1300 failed += test_result ("Arcturus rise above -20 deg minute on 2006/01/17 at (15 E,51 S)", date.minutes, 50, 0);
1301
1302 ln_get_date (rst.transit, &date);
1303 failed += test_result ("Arcturus transit hour on 2006/01/17 at (15 E,51 S)", date.hours, 5, 0);
1304 failed += test_result ("Arcturus transit minute on 2006/01/17 at (15 E,51 S)", date.minutes, 29, 0);
1305
1306 ln_get_date (rst.set, &date);
1307 failed += test_result ("Arcturus set bellow -20 deg hour on 2006/01/17 at (15 E,51 S)", date.hours, 12, 0);
1308 failed += test_result ("Arcturus set bellow -20 deg minute on 2006/01/17 at (15 E,51 S)", date.minutes, 4, 0);
1309 }
1310
1311 ret = ln_get_solar_rst (JD, &observer, &rst);
1312 failed += test_result ("Sun sometimes rise on 2006/01/17 at 15 E, 51 S", ret, 0, 0);
1313
1314 if (!ret)
1315 {
1316 ln_get_date (rst.rise, &date);
1317 failed += test_result ("Sun rise hour on 2006/01/17 at (15 E,51 S)", date.hours, 3, 0);
1318 failed += test_result ("Sun rise minute on 2006/01/17 at (15 E,51 S)", date.minutes, 11, 0);
1319
1320 ln_get_date (rst.transit, &date);
1321 failed += test_result ("Sun transit hour on 2006/01/17 at (15 E,51 S)", date.hours, 11, 0);
1322 failed += test_result ("Sun transit minute on 2006/01/17 at (15 E,51 S)", date.minutes, 9, 0);
1323
1324 ln_get_date (rst.set, &date);
1325 failed += test_result ("Sun set hour on 2006/01/17 at (15 E,51 S)", date.hours, 19, 0);
1326 failed += test_result ("Sun set minute on 2006/01/17 at (15 E,51 S)", date.minutes, 7, 0);
1327 }
1328
1329 observer.lat = 37;
1330
1331 object.dec = -54;
1332 failed += test_result ("Object at dec -54 never rise at 37 N", ln_get_object_rst (JD, &observer, &object, &rst), -1, 0);
1333
1334 object.dec = -52;
1335 failed += test_result ("Object at dec -52 rise at 37 N", ln_get_object_rst (JD, &observer, &object, &rst), 0, 0);
1336
1337 object.dec = 54;
1338 failed += test_result ("Object at dec 54 is always above the horizon at 37 N", ln_get_object_rst (JD, &observer, &object, &rst), 1, 0);
1339
1340 object.dec = 52;
1341 failed += test_result ("Object at dec 52 rise at 37 N", ln_get_object_rst (JD, &observer, &object, &rst), 0, 0);
1342
1343 observer.lat = -37;
1344
1345 object.dec = 54;
1346 failed += test_result ("Object at dec 54 never rise at 37 S", ln_get_object_rst (JD, &observer, &object, &rst), -1, 0);
1347
1348 object.dec = 52;
1349 failed += test_result ("Object at dec 52 rise at 37 S", ln_get_object_rst (JD, &observer, &object, &rst), 0, 0);
1350
1351 object.dec = -54;
1352 failed += test_result ("Object at dec -54 is always above the horizon at 37 S", ln_get_object_rst (JD, &observer, &object, &rst), 1, 0);
1353
1354 object.dec = -52;
1355 failed += test_result ("Object at dec -52 rise at 37 S", ln_get_object_rst (JD, &observer, &object, &rst), 0, 0);
1356
1357 /* Venus on 1988 March 20 at Boston */
1358 date.years = 1988;
1359 date.months = 3;
1360 date.days = 20;
1361
1362 date.hours = 0;
1363 date.minutes = 0;
1364 date.seconds = 0;
1365
1366 JD = ln_get_julian_day (&date);
1367
1368 observer.lng = -71.0833;
1369 observer.lat = 42.3333;
1370
1371 ret = ln_get_venus_rst (JD, &observer, &rst);
1372 failed += test_result ("Venus sometime rise on 1988/03/20 at Boston", ret, 0, 0);
1373
1374 if (!ret)
1375 {
1376 ln_get_date (rst.rise, &date);
1377 failed += test_result ("Venus rise hour on 1988/03/20 at Boston", date.hours, 12, 0);
1378 failed += test_result ("Venus rise minute on 1988/03/20 at Boston", date.minutes, 25, 0);
1379
1380 ln_get_date (rst.transit, &date);
1381 failed += test_result ("Venus transit hour on 1988/03/20 at Boston", date.hours, 19, 0);
1382 failed += test_result ("Venus transit minute on 1988/03/20 at Boston", date.minutes, 41, 0);
1383
1384 ln_get_date (rst.set, &date);
1385 failed += test_result ("Venus set hour on 1988/03/20 at Boston", date.hours, 2, 0);
1386 failed += test_result ("Venus set minute on 1988/03/20 at Boston", date.minutes, 55, 0);
1387 }
1388
1389 return failed;
1390 }
1391
ell_rst_test()1392 int ell_rst_test ()
1393 {
1394 struct ln_lnlat_posn observer;
1395 struct ln_ell_orbit orbit;
1396 struct ln_date date;
1397 struct ln_rst_time rst;
1398 struct ln_equ_posn pos;
1399 double JD;
1400 int ret;
1401 int failed = 0;
1402
1403 /* Comment C/1996 B2 (Hyakutake) somewhere at Japan */
1404
1405 observer.lng = 135;
1406 observer.lat = 35;
1407
1408 date.years = 1996;
1409 date.months = 5;
1410 date.days = 1;
1411
1412 date.hours = 0;
1413 date.minutes = 0;
1414 date.seconds = 0;
1415
1416 orbit.JD = ln_get_julian_day (&date);
1417 orbit.JD += 0.39481;
1418 orbit.a = 1014.2022026431;
1419 orbit.e = 0.9997730;
1420 orbit.i = 124.92379;
1421 orbit.omega = 188.04546;
1422 orbit.w = 130.17654;
1423 orbit.n = 0;
1424
1425 date.years = 1996;
1426 date.months = 3;
1427 date.days = 24;
1428
1429 date.hours = date.minutes = 0;
1430 date.seconds = 0.0;
1431
1432 JD = ln_get_julian_day (&date);
1433
1434 ln_get_ell_body_equ_coords (JD, &orbit, &pos);
1435 failed += test_result ("(RA) for Hyakutake 1996/03/28 00:00", pos.ra, 220.8554, 0.001);
1436 failed += test_result ("(Dec) for Hyakutake 1996/03/28 00:00", pos.dec, 36.5341, 0.001);
1437
1438 date.days = 28;
1439
1440 date.hours = 10;
1441 date.minutes = 42;
1442
1443 JD = ln_get_julian_day (&date);
1444
1445 ln_get_ell_body_equ_coords (JD, &orbit, &pos);
1446 failed += test_result ("(RA) for Hyakutake 1996/03/28 10:42", pos.ra, 56.2140, 0.001);
1447 failed += test_result ("(Dec) for Hyakutake 1996/03/28 10:42", pos.dec, 74.4302, 0.001);
1448
1449 date.days = 23;
1450
1451 date.hours = 17;
1452 date.minutes = 38;
1453 date.seconds = 0;
1454
1455 JD = ln_get_julian_day (&date);
1456
1457 ln_get_ell_body_equ_coords (JD, &orbit, &pos);
1458 failed += test_result ("(RA) for Hyakutake 1996/03/23 17:38", pos.ra, 221.2153, 0.001);
1459 failed += test_result ("(Dec) for Hyakutake 1996/03/23 17:38", pos.dec, 32.4803, 0.001);
1460
1461 JD = ln_get_julian_day (&date);
1462
1463 ret = ln_get_ell_body_rst (JD, &observer, &orbit, &rst);
1464 failed += test_result ("Hyakutake sometime rise on 1996/03/23 at 135 E, 35 N", ret, 0, 0);
1465
1466 if (!ret)
1467 {
1468 ln_get_date (rst.rise, &date);
1469 failed += test_result ("Hyakutake rise hour on 1996/03/23 at 135 E, 35 N", date.hours, 9, 0);
1470 failed += test_result ("Hyakutake rise minute on 1996/03/23 at 135 E, 35 N", date.minutes, 31, 0);
1471
1472 ln_get_date (rst.transit, &date);
1473 failed += test_result ("Hyakutake transit hour on 1996/03/23 at 135 E, 35 N", date.hours, 17, 0);
1474 failed += test_result ("Hyakutake transit minute on 1996/03/23 at 135 E, 35 N", date.minutes, 27, 0);
1475
1476 ln_get_date (rst.set, &date);
1477 failed += test_result ("Hyakutake set hour on 1996/03/23 at 135 E, 35 N", date.hours, 1, 0);
1478 failed += test_result ("Hyakutake set minute on 1996/03/23 at 135 E, 35 N", date.minutes, 49, 0);
1479 }
1480
1481 ret = ln_get_ell_body_next_rst (JD, &observer, &orbit, &rst);
1482 failed += test_result ("Hyakutake sometime rise on 1996/03/23 at 135 E, 35 N", ret, 0, 0);
1483
1484 if (!ret)
1485 {
1486 ln_get_date (rst.rise, &date);
1487 failed += test_result ("Hyakutake next rise hour on 1996/03/23 at 135 E, 35 N", date.hours, 9, 0);
1488 failed += test_result ("Hyakutake next rise minute on 1996/03/23 at 135 E, 35 N", date.minutes, 31, 0);
1489
1490 ln_get_date (rst.transit, &date);
1491 failed += test_result ("Hyakutake next transit hour on 1996/03/24 at 135 E, 35 N", date.hours, 17, 0);
1492 failed += test_result ("Hyakutake next transit minute on 1996/03/24 at 135 E, 35 N", date.minutes, 4, 0);
1493
1494 ln_get_date (rst.set, &date);
1495 failed += test_result ("Hyakutake next set hour on 1996/03/23 at 135 E, 35 N", date.hours, 1, 0);
1496 failed += test_result ("Hyakutake next set minute on 1996/03/23 at 135 E, 35 N", date.minutes, 49, 0);
1497 }
1498
1499 return failed;
1500 }
1501
hyp_future_rst_test()1502 int hyp_future_rst_test ()
1503 {
1504 struct ln_lnlat_posn observer;
1505 struct ln_hyp_orbit orbit;
1506 struct ln_date date;
1507 struct ln_rst_time rst;
1508 double JD;
1509 int ret;
1510 int failed = 0;
1511
1512 observer.lng = 15;
1513 observer.lat = 50;
1514
1515 /* C/2006 P1 (McNaught) */
1516
1517 orbit.q = 0.170742005109787;
1518 orbit.e = 1.00001895427704;
1519 orbit.i = 77.8348999023438;
1520 orbit.w = 155.977096557617;
1521 orbit.omega = 267.414398193359;
1522 orbit.JD = 2454113.251;
1523
1524 date.years = 2007;
1525 date.months = 1;
1526 date.days = 17;
1527
1528 date.hours = 12;
1529 date.minutes = 0;
1530 date.seconds = 0.0;
1531
1532 JD = ln_get_julian_day (&date);
1533
1534 ret = ln_get_hyp_body_next_rst_horizon (JD, &observer, &orbit, 0, &rst);
1535 failed += test_result ("McNaught rise on 2997/01/18 at 15 E, 50 N", ret, 0, 0);
1536
1537 if (!ret)
1538 {
1539 ln_get_date (rst.rise, &date);
1540 failed += test_result ("McNaught rise hour on 2007/01/18 at 15 E, 50 N", date.hours, 9, 0);
1541 failed += test_result ("McNaught rise minute on 2007/01/18 at 15 E, 50 N", date.minutes, 6, 0);
1542
1543 ln_get_date (rst.transit, &date);
1544 failed += test_result ("McNaught transit hour on 2007/01/18 at 15 E, 50 N", date.hours, 11, 0);
1545 failed += test_result ("McNaught transit minute on 2007/01/18 at 15 E, 50 N", date.minutes, 38, 0);
1546
1547 ln_get_date (rst.set, &date);
1548 failed += test_result ("McNaught set hour on 2007/01/17 at 15 E, 50 N", date.hours, 14, 0);
1549 failed += test_result ("McNaught set minute on 2007/01/17 at 15 E, 50 N", date.minutes, 37, 0);
1550 }
1551
1552 ret = ln_get_hyp_body_next_rst_horizon (JD, &observer, &orbit, 15, &rst);
1553 failed += test_result ("McNaught does not rise above 15 degrees on 2007/01/17 at 15 E, 50 N", ret, -1, 0);
1554
1555 return failed;
1556 }
1557
body_future_rst_test()1558 int body_future_rst_test ()
1559 {
1560 struct ln_lnlat_posn observer;
1561 struct ln_date date;
1562 struct ln_rst_time rst;
1563 double JD;
1564 int ret;
1565 int failed = 0;
1566
1567 observer.lng = 0;
1568 observer.lat = 85;
1569
1570 date.years = 2006;
1571 date.months = 1;
1572 date.days = 1;
1573
1574 date.hours = date.minutes = 0;
1575 date.seconds = 0.0;
1576
1577 JD = ln_get_julian_day (&date);
1578
1579 ret = ln_get_body_next_rst_horizon_future (JD, &observer, ln_get_solar_equ_coords, LN_SOLAR_STANDART_HORIZON, 300, &rst);
1580 failed += test_result ("Sun is above horizon sometimes at 0, 85 N", ret, 0, 0);
1581
1582 if (!ret)
1583 {
1584 ln_get_date (rst.rise, &date);
1585 failed += test_result ("Solar next rise years at 0, 85 N", date.years, 2006, 0);
1586 failed += test_result ("Solar next rise months at 0, 85 N", date.months, 3, 0);
1587 failed += test_result ("Solar next rise days at 0, 85 N", date.days, 7, 0);
1588 failed += test_result ("Solar next rise hour at 0, 85 N", date.hours, 10, 0);
1589 failed += test_result ("Solar next rise minute at 0, 85 N", date.minutes, 19, 0);
1590
1591 ln_get_date (rst.transit, &date);
1592 failed += test_result ("Solar next transit years at 0, 85 N", date.years, 2006, 0);
1593 failed += test_result ("Solar next transit months at 0, 85 N", date.months, 3, 0);
1594 failed += test_result ("Solar next transit days at 0, 85 N", date.days, 7, 0);
1595 failed += test_result ("Solar next transit hour at 0 E, 85 N", date.hours, 12, 0);
1596 failed += test_result ("Solar next transit minute at 0 E, 85 N", date.minutes, 10, 0);
1597
1598 ln_get_date (rst.set, &date);
1599 failed += test_result ("Solar next set years at 0, 85 N", date.years, 2006, 0);
1600 failed += test_result ("Solar next set months at 0, 85 N", date.months, 3, 0);
1601 failed += test_result ("Solar next set days at 0, 85 N", date.days, 7, 0);
1602 failed += test_result ("Solar next set hour at 0 E, 85 N", date.hours, 14, 0);
1603 failed += test_result ("Solar next set minute at 0, 85 N", date.minutes, 7, 0);
1604 }
1605
1606 ret = ln_get_body_next_rst_horizon_future (JD, &observer, ln_get_solar_equ_coords, 0, 300, &rst);
1607 failed += test_result ("Sun is above 0 horizon sometimes at 0, 85 N", ret, 0, 0);
1608
1609 if (!ret)
1610 {
1611 ln_get_date (rst.rise, &date);
1612 failed += test_result ("Solar next rise years at 0, 85 N with 0 horizon", date.years, 2006, 0);
1613 failed += test_result ("Solar next rise months at 0, 85 N with 0 horizon", date.months, 3, 0);
1614 failed += test_result ("Solar next rise days at 0, 85 N with 0 horizon", date.days, 6, 0);
1615 failed += test_result ("Solar next rise hour at 0, 85 N with 0 horizon", date.hours, 10, 0);
1616 failed += test_result ("Solar next rise minute at 0, 85 N with 0 horizon", date.minutes, 19, 0);
1617
1618 ln_get_date (rst.transit, &date);
1619 failed += test_result ("Solar next transit years at 0, 85 N with 0 horizon", date.years, 2006, 0);
1620 failed += test_result ("Solar next transit months at 0, 85 N with 0 horizon", date.months, 3, 0);
1621 failed += test_result ("Solar next transit days at 0, 85 N with 0 horizon", date.days, 6, 0);
1622 failed += test_result ("Solar next transit hour at 0 E, 85 N with 0 horizon", date.hours, 12, 0);
1623 failed += test_result ("Solar next transit minute at 0 E, 85 N with 0 horizon", date.minutes, 10, 0);
1624
1625 ln_get_date (rst.set, &date);
1626 failed += test_result ("Solar next set years at 0, 85 N with 0 horizon", date.years, 2006, 0);
1627 failed += test_result ("Solar next set months at 0, 85 N with 0 horizon", date.months, 3, 0);
1628 failed += test_result ("Solar next set days at 0, 85 N with 0 horizon", date.days, 6, 0);
1629 failed += test_result ("Solar next set hour at 0 E, 85 N with 0 horizon", date.hours, 14, 0);
1630 failed += test_result ("Solar next set minute at 0, 85 N with 0 horizon", date.minutes, 7, 0);
1631 }
1632
1633 return failed;
1634 }
1635
parallax_test()1636 int parallax_test ()
1637 {
1638 struct ln_equ_posn mars, parallax;
1639 struct ln_lnlat_posn observer;
1640 struct ln_dms dms;
1641 struct ln_date date;
1642 double jd;
1643 int failed = 0;
1644
1645 dms.neg = 0;
1646 dms.degrees = 33;
1647 dms.minutes = 21;
1648 dms.seconds = 22;
1649
1650 observer.lat = ln_dms_to_deg (&dms);
1651
1652 dms.neg = 1;
1653 dms.degrees = 116;
1654 dms.minutes = 51;
1655 dms.seconds = 47;
1656
1657 observer.lng = ln_dms_to_deg (&dms);
1658
1659 date.years = 2003;
1660 date.months = 8;
1661 date.days = 28;
1662
1663 date.hours = 3;
1664 date.minutes = 17;
1665 date.seconds = 0;
1666
1667 jd = ln_get_julian_day (&date);
1668
1669 ln_get_mars_equ_coords (jd, &mars);
1670
1671 ln_get_parallax (&mars, ln_get_mars_earth_dist (jd), &observer, 1706, jd, ¶llax);
1672
1673 /* parallax is hard to calculate, so we allow relatively big errror */
1674
1675 failed += test_result ("Mars RA parallax for Palomar observatory at 2003/08/28 3:17 UT ", parallax.ra, 0.0053917, 0.00001);
1676 failed += test_result ("Mars DEC parallax for Palomar observatory at 2003/08/28 3:17 UT ", parallax.dec, -14.1 / 3600.0, 0.00002);
1677
1678 return failed;
1679 }
1680
angular_test()1681 int angular_test ()
1682 {
1683 int failed = 0;
1684 double d;
1685 struct ln_equ_posn posn1, posn2;
1686
1687 /* alpha Bootes (Arcturus) */
1688 posn1.ra = 213.9154;
1689 posn1.dec = 19.1825;
1690
1691 /* alpha Virgo (Spica) */
1692 posn2.ra = 201.2983;
1693 posn2.dec = -11.1614;
1694
1695 d = ln_get_angular_separation(&posn1, &posn2);
1696 failed += test_result ("(Angular) Separation of Arcturus and Spica ", d, 32.79302684, 0.00001);
1697
1698 d = ln_get_rel_posn_angle(&posn1, &posn2);
1699 failed += test_result ("(Angular) Position Angle of Arcturus and Spica ", d, 22.39042787, 0.00001);
1700
1701 return failed;
1702 }
1703
utility_test()1704 int utility_test()
1705 {
1706 struct ln_dms dms;
1707 double deg = -1.23, deg2 = 1.23, deg3 = -0.5;
1708
1709 ln_deg_to_dms (deg, &dms);
1710 printf("TEST deg %f ==> deg %c%d min %d sec %f\n", deg, dms.neg ? '-' : '+', dms.degrees, dms.minutes, dms.seconds);
1711 ln_deg_to_dms (deg2, &dms);
1712 printf("TEST deg %f ==> deg %c%d min %d sec %f\n", deg2, dms.neg ? '-' : '+', dms.degrees, dms.minutes, dms.seconds);
1713 ln_deg_to_dms (deg3, &dms);
1714 printf("TEST deg %f ==> deg %c%d min %d sec %f\n", deg3, dms.neg ? '-' : '+', dms.degrees, dms.minutes, dms.seconds);
1715 return 0;
1716 }
1717
airmass_test()1718 int airmass_test()
1719 {
1720 int failed = 0;
1721 double x;
1722
1723 double X = ln_get_airmass (90, 750.0);
1724 failed += test_result ("(Airmass) Airmass at Zenith", X, 1, 0);
1725
1726 X = ln_get_airmass (10, 750.0);
1727 failed += test_result ("(Airmass) Airmass at 10 degrees altitude", X, 5.3, 0.1);
1728
1729 X = ln_get_alt_from_airmass (1, 750.0);
1730 failed += test_result ("(Airmass) Altitude at airmass 1", X, 90, 0);
1731
1732 for (x = -10; x < 90; x += 10.54546456)
1733 {
1734 X = ln_get_alt_from_airmass (ln_get_airmass (x, 750.0), 750.0);
1735 failed += test_result ("(Airmass) Altitude->Airmass->Altitude at 10 degrees", X, x, 0.000000001);
1736 }
1737
1738 return failed;
1739 }
1740
main()1741 int main ()
1742 {
1743 int failed = 0;
1744
1745 failed += julian_test();
1746 failed += dynamical_test();
1747 failed += heliocentric_test ();
1748 failed += sidereal_test();
1749 failed += nutation_test();
1750 failed += transform_test();
1751 failed += solar_coord_test ();
1752 failed += aberration_test();
1753 failed += precession_test();
1754 failed += apparent_position_test ();
1755 failed += vsop87_test();
1756 failed += lunar_test ();
1757 failed += elliptic_motion_test();
1758 failed += parabolic_motion_test ();
1759 failed += hyperbolic_motion_test ();
1760 failed += rst_test ();
1761 failed += ell_rst_test ();
1762 failed += hyp_future_rst_test ();
1763 failed += body_future_rst_test ();
1764 failed += parallax_test ();
1765 failed += angular_test();
1766 failed += utility_test();
1767 failed += airmass_test ();
1768
1769 printf ("Test completed: %d tests, %d errors.\n", test_number, failed);
1770
1771 return (failed > 0);
1772 }
1773