source: branches/FACT++_lidctrl_usb/pal/palTest.c

Last change on this file was 18712, checked in by tbretz, 8 years ago
Updated to PAL 0.9.7 (adds mainly light deflection to abberation which was previously missing)
File size: 60.4 KB
Line 
1/*
2*+
3* Name:
4* palTest
5
6* Purpose:
7* Test the PAL library
8
9* Language:
10* Starlink ANSI C
11
12* Type of Module:
13* Application
14
15* Description:
16* Test the PAL library is functioning correctly. Uses some of the SLA test code.
17
18* Authors:
19* TIMJ: Tim Jenness (JAC, Hawaii)
20* {enter_new_authors_here}
21
22* History:
23* 2012-02-08 (TIMJ):
24* Initial version
25* Adapted with permission from the Fortran SLALIB library.
26* {enter_further_changes_here}
27
28* Copyright:
29* Copyright (C) 2012 Science and Technology Facilities Council.
30* All Rights Reserved.
31
32* Licence:
33* This program is free software; you can redistribute it and/or
34* modify it under the terms of the GNU General Public License as
35* published by the Free Software Foundation; either version 3 of
36* the License, or (at your option) any later version.
37*
38* This program is distributed in the hope that it will be
39* useful, but WITHOUT ANY WARRANTY; without even the implied
40* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
41* PURPOSE. See the GNU General Public License for more details.
42*
43* You should have received a copy of the GNU General Public License
44* along with this program; if not, write to the Free Software
45* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301
46* USA.
47
48* Bugs:
49* {note_any_bugs_here}
50*-
51*/
52
53#include <stdlib.h>
54#include <stdio.h>
55#include <string.h>
56
57#include "pal.h"
58#include "palmac.h"
59
60static int verbose = 1;
61
62/* Support functions to allow to test results.
63 viv and vvd match the SOFA/ERFA implementations */
64
65static void viv(int ival, int ivalok, const char *func, const char *test,
66 int *status)
67/*
68** - - - -
69** v i v
70** - - - -
71**
72** Validate an integer result.
73**
74** Internal function used by t_sofa_c program.
75**
76** Given:
77** ival int value computed by function under test
78** ivalok int correct value
79** func char[] name of function under test
80** test char[] name of individual test
81**
82** Given and returned:
83** status int set to FALSE if test fails
84**
85** This revision: 2009 November 4
86*/
87{
88 if (ival != ivalok) {
89 *status = 1;
90 printf("%s failed: %s want %d got %d\n",
91 func, test, ivalok, ival);
92 } else if (verbose) {
93 printf("%s passed: %s want %d got %d\n",
94 func, test, ivalok, ival);
95 }
96 return;
97}
98
99static void vvd(double val, double valok, double dval,
100 const char *func, const char *test, int *status)
101/*
102** - - - -
103** v v d
104** - - - -
105**
106** Validate a double result.
107**
108** Internal function used by t_sofa_c program.
109**
110** Given:
111** val double value computed by function under test
112** valok double expected value
113** dval double maximum allowable error
114** func char[] name of function under test
115** test char[] name of individual test
116**
117** Given and returned:
118** status int set to FALSE if test fails
119**
120** This revision: 2008 June 8
121*/
122{
123 double a, f; /* absolute and fractional error */
124
125
126 a = val - valok;
127 if (fabs(a) > dval) {
128 f = fabs(valok / a);
129 *status = 1;
130 printf("%s failed: %s want %.20g got %.20g (1/%.3g)\n",
131 func, test, valok, val, f);
132 } else if (verbose) {
133 printf("%s passed: %s want %.20g got %.20g\n",
134 func, test, valok, val);
135 }
136 return;
137}
138
139/* Verify a string */
140static void vcs( const char * val, const char * valok,
141 const char * func, const char * test,
142 int *status ) {
143
144 if (strcmp(val, valok) != 0) {
145 *status = 1;
146 printf("%s failed: %s want %s got %s\n",
147 func, test, valok, val );
148 } else if (verbose) {
149 printf("%s passed: %s want %s got %s\n",
150 func, test, valok, val );
151 }
152 return;
153
154}
155
156/* Verify the 3x3 rmat matrix */
157static void
158vrmat( double rmat[3][3], double expected[3][3], const char * func,
159 double dval, int * status ) {
160 int i;
161 char buf[10];
162 for( i = 0; i < 3; i++ ) {
163 int j;
164 for( j = 0; j < 3; j++ ) {
165 sprintf( buf, "%d,%d", i, j );
166 vvd( rmat[i][j], expected[i][j], dval, func, buf, status );
167 }
168 }
169}
170
171/* Verify a vector */
172static void
173vvec( int len, double *vec, double *expected, const char *func,
174 int *status ) {
175 int i;
176 char buf[10];
177 for( i = 0; i < len; i++ ) {
178 sprintf( buf, "%d", i );
179 vvd( vec[i], expected[i], 1e-12, func, buf, status );
180 }
181}
182
183/******************************************************************/
184/* TEST FUNCTIONS */
185
186/* Adding E-terms */
187
188static void t_addet( int *status ) {
189 double r1,d1,r2,d2;
190 double rm = 2.;
191 double dm = -1.;
192 double eq = 1975.;
193
194
195 palAddet ( rm, dm, eq, &r1, &d1 );
196 vvd ( r1 - rm, 2.983864874295250e-6, 1e-12, "palAddet",
197 "R", status );
198 vvd ( d1 - dm, 2.379650804185118e-7, 1e-12, "palAddet",
199 "D", status );
200
201 palSubet ( r1, d1, eq, &r2, &d2 );
202 vvd ( r2 - rm, 0, 1e-12, "palSubet", "R", status );
203 vvd ( d2 - dm, 0, 1e-12, "palSubet", "D", status );
204
205}
206
207static void t_afin( int * status ) {
208
209 int j;
210 int i = 1;
211 double d = 0.0;
212 const char * s = "12 34 56.7 |";
213 const char * s2 = "45 00 00.000 ";
214
215 palDafin (s, &i, &d, &j);
216 viv ( i, 12, "palDafin", "I", status );
217 vvd ( d, 0.2196045986911432, 1e-12, "palDafin", "A",
218 status );
219 viv ( j, 0, "palDafin", "J", status );
220
221 i = 1;
222 palDafin (s2, &i, &d, &j);
223 viv ( i, 14, "palDafin", "I", status );
224 vvd ( d, PAL__DPI/4.0, 1e-12, "palDafin", "A",
225 status );
226 viv ( j, 0, "palDafin", "J", status );
227}
228
229/* Altaz */
230
231static void t_altaz( int *status ) {
232 double az, azd, azdd, el, eld, eldd, pa, pad, padd;
233 palAltaz( 0.7, -0.7, -0.65,
234 &az, &azd, &azdd, &el, &eld, &eldd, &pa, &pad, &padd );
235
236 vvd ( az, 4.400560746660174, 1e-12, "palAltaz",
237 "AZ", status );
238 vvd ( azd, -0.2015438937145421, 1e-13, "palAltaz",
239 "AZD", status );
240 vvd ( azdd, -0.4381266949668748, 1e-13, "palAltaz",
241 "AZDD", status );
242 vvd ( el, 1.026646506651396, 1e-12, "palAltaz",
243 "EL", status );
244 vvd ( eld, -0.7576920683826450, 1e-13, "palAltaz",
245 "ELD", status );
246 vvd ( eldd, 0.04922465406857453, 1e-14, "palAltaz",
247 "ELDD", status );
248 vvd ( pa, 1.707639969653937, 1e-12, "palAltaz",
249 "PA", status );
250 vvd ( pad, 0.4717832355365627, 1e-13, "palAltaz",
251 "PAD", status );
252 vvd ( padd, -0.2957914128185515, 1e-13, "palAltaz",
253 "PADD", status );
254}
255
256/* Airmass */
257
258static void t_airmas( int *status ) {
259 vvd ( palAirmas ( 1.2354 ), 3.015698990074724,
260 1e-12, "palAirmas", " ", status );
261}
262
263/* Apparent to mean place */
264
265static void t_amp ( int *status ) {
266 double rm, dm;
267
268 /* Original SLA test is not accurate since palMapqk
269 differs from slaMapqk */
270 palAmp ( 2.345, -1.234, 50100, 1990, &rm, &dm );
271 vvd ( rm, 2.344472180027961, 1e-6, "palAmp", "R",
272 status );
273 vvd ( dm, -1.233573099847705, 1e-7, "palAmp", "D",
274 status );
275
276 /* This is the palMapqk test */
277 palAmp( 1.234, -0.567, 55927.0, 2010.0, &rm, &dm );
278 vvd( rm, 1.233512033578303857, 1.0E-12, "palAmp", "rm", status );
279 vvd( dm, -0.56702909748530827549, 1.0E-12, "palAmp", "dm", status );
280}
281
282/* Apparent to Observed place */
283
284static void t_aop ( int *status ) {
285
286 int i;
287
288 double rap, dap, date, dut, elongm, phim, hm, xp, yp,
289 tdk, pmb, rh, wl, tlr, aob, zob, hob, dob, rob, aoprms[14];
290
291 dap = -0.1234;
292 date = 51000.1;
293 dut = 25.0;
294 elongm = 2.1;
295 phim = 0.5;
296 hm = 3000.0;
297 xp = -0.5e-6;
298 yp = 1.0e-6;
299 tdk = 280.0;
300 pmb = 550.0;
301 rh = 0.6;
302 tlr = 0.006;
303
304 for (i=1; i<=3; i++) {
305
306 if ( i == 1 ) {
307 rap = 2.7;
308 wl = 0.45;
309 } else if ( i == 2 ) {
310 rap = 2.345;
311 } else {
312 wl = 1.0e6;
313 }
314
315 palAop ( rap, dap, date, dut, elongm, phim, hm, xp, yp,
316 tdk, pmb, rh, wl, tlr, &aob, &zob, &hob, &dob, &rob );
317
318 if ( i == 1 ) {
319 vvd( aob, 1.812817787123283034, 1e-10, "palAop",
320 "lo aob", status );
321 vvd( zob, 1.393860816635714034, 1e-8, "palAop",
322 "lo zob", status );
323 vvd( hob, -1.297808009092456683, 1e-8, "palAop",
324 "lo hob", status );
325 vvd( dob, -0.122967060534561, 1e-8, "palAop",
326 "lo dob", status );
327 vvd( rob, 2.699270287872084, 1e-8, "palAop",
328 "lo rob", status );
329 } else if ( i == 2 ) {
330 vvd( aob, 2.019928026670621442, 1e-10, "palAop",
331 "aob/o", status );
332 vvd( zob, 1.101316172427482466, 1e-10, "palAop",
333 "zob/o", status );
334 vvd( hob, -0.9432923558497740862, 1e-10, "palAop",
335 "hob/o", status );
336 vvd( dob, -0.1232144708194224, 1e-10, "palAop",
337 "dob/o", status );
338 vvd( rob, 2.344754634629428, 1e-10, "palAop",
339 "rob/o", status );
340 } else {
341 vvd( aob, 2.019928026670621442, 1e-10, "palAop",
342 "aob/r", status );
343 vvd( zob, 1.101267532198003760, 1e-10, "palAop",
344 "zob/r", status );
345 vvd( hob, -0.9432533138143315937, 1e-10, "palAop",
346 "hob/r", status );
347 vvd( dob, -0.1231850665614878, 1e-10, "palAop",
348 "dob/r", status );
349 vvd( rob, 2.344715592593984, 1e-10, "palAop",
350 "rob/r", status );
351 }
352 }
353
354 date = 48000.3;
355 wl = 0.45;
356
357 palAoppa ( date, dut, elongm, phim, hm, xp, yp, tdk,
358 pmb, rh, wl, tlr, aoprms );
359 vvd( aoprms[0], 0.4999993892136306, 1e-13, "palAoppa",
360 "0", status );
361 vvd( aoprms[1], 0.4794250025886467, 1e-13, "palAoppa",
362 "1", status );
363 vvd( aoprms[2], 0.8775828547167932, 1e-13, "palAoppa",
364 "2", status );
365 vvd( aoprms[3], 1.363180872136126e-6, 1e-13, "palAoppa",
366 "3", status );
367 vvd( aoprms[4], 3000.0, 1e-10, "palAoppa", "4",
368 status );
369 vvd( aoprms[5], 280.0, 1e-11, "palAoppa", "5",
370 status );
371 vvd( aoprms[6], 550.0, 1e-11, "palAoppa", "6",
372 status );
373 vvd( aoprms[7], 0.6, 1e-13, "palAoppa", "7",
374 status );
375 vvd( aoprms[8], 0.45, 1e-13, "palAoppa", "8",
376 status );
377 vvd( aoprms[9], 0.006, 1e-15, "palAoppa", "9",
378 status );
379 vvd( aoprms[10], 0.0001562803328459898, 1e-13,
380 "palAoppa", "10", status );
381 vvd( aoprms[11], -1.792293660141e-7, 1e-13,
382 "palAoppa", "11", status );
383 vvd( aoprms[12], 2.101874231495843, 1e-13,
384 "palAoppa", "12", status );
385 vvd( aoprms[13], 7.601916802079765, 1e-8,
386 "palAoppa", "13", status );
387
388 palOap ( "r", 1.6, -1.01, date, dut, elongm, phim,
389 hm, xp, yp, tdk, pmb, rh, wl, tlr, &rap, &dap );
390 vvd( rap, 1.601197569844787, 1e-10, "palOap",
391 "rr", status );
392 vvd( dap, -1.012528566544262, 1e-10, "palOap",
393 "rd", status );
394 palOap ( "h", -1.234, 2.34, date, dut, elongm, phim,
395 hm, xp, yp, tdk, pmb, rh, wl, tlr, &rap, &dap );
396 vvd( rap, 5.693087688154886463, 1e-10, "palOap",
397 "hr", status );
398 vvd( dap, 0.8010281167405444, 1e-10, "palOap",
399 "hd", status );
400 palOap ( "a", 6.1, 1.1, date, dut, elongm, phim,
401 hm, xp, yp, tdk, pmb, rh, wl, tlr, &rap, &dap );
402 vvd( rap, 5.894305175192448940, 1e-10, "palOap",
403 "ar", status );
404 vvd( dap, 1.406150707974922, 1e-10, "palOap",
405 "ad", status );
406
407 palOapqk ( "r", 2.1, -0.345, aoprms, &rap, &dap );
408 vvd( rap, 2.10023962776202, 1e-10, "palOapqk",
409 "rr", status );
410 vvd( dap, -0.3452428692888919, 1e-10, "palOapqk",
411 "rd", status );
412 palOapqk ( "h", -0.01, 1.03, aoprms, &rap, &dap );
413 vvd( rap, 1.328731933634564995, 1e-10, "palOapqk",
414 "hr", status );
415 vvd( dap, 1.030091538647746, 1e-10, "palOapqk",
416 "hd", status );
417 palOapqk ( "a", 4.321, 0.987, aoprms, &rap, &dap );
418 vvd( rap, 0.4375507112075065923, 1e-10, "palOapqk",
419 "ar", status );
420 vvd( dap, -0.01520898480744436, 1e-10, "palOapqk",
421 "ad", status );
422
423 palAoppat ( date + PAL__DS2R, aoprms );
424 vvd( aoprms[13], 7.602374979243502, 1e-8, "palAoppat",
425 " ", status );
426}
427
428/* Bearings */
429
430static void t_bear( int *status ) {
431 double a1 = 1.234;
432 double b1 = -0.123;
433 double a2 = 2.345;
434 double b2 = 0.789;
435
436 double d1[3];
437 double d2[3];
438
439 vvd ( palDbear ( a1, b1, a2, b2 ), 0.7045970341781791,
440 1e-12, "palDbear", " ", status );
441 palDcs2c ( a1, b1, d1 );
442 palDcs2c ( a2, b2, d2 );
443
444 vvd ( palDpav ( d1, d2 ), 0.7045970341781791,
445 1e-12, "palDpav", " ", status );
446
447}
448
449/* Calendar to MJD */
450
451static void t_caldj( int *status ) {
452 int j;
453 double djm;
454
455 palCaldj ( 1999, 12, 31, &djm, &j );
456 vvd ( djm, 51543, 0, "palCaldj", " ", status );
457 viv ( j, 0, "palCaldj", "J", status );
458}
459
460/* palDaf2r */
461
462static void t_caf2r( int * status ) {
463 int j;
464 double dr;
465
466 palDaf2r ( 76, 54, 32.1, &dr, &j );
467 vvd ( dr, 1.342313819975276, 1e-12, "palDaf2r",
468 "r", status );
469 viv ( j, 0, "palDaf2r", "j", status );
470}
471
472/* Test palDcc2s routines */
473
474static void t_cc2s( int * status ) {
475 double dv[3] = { 100., -50., 25. };
476 double da, db;
477
478 palDcc2s ( dv, &da, &db );
479 vvd ( da, -0.4636476090008061, 1e-12, "palDcc2s",
480 "A", status );
481 vvd ( db, 0.2199879773954594, 1e-12, "palDcc2s",
482 "B", status );
483}
484
485/* palDd2tf */
486
487static void t_cd2tf( int *status ) {
488 int ihmsf[4];
489 char s;
490
491 palDd2tf ( 4, -0.987654321, &s, ihmsf );
492 viv ( s, '-', "palDd2tf", "S", status );
493 viv ( ihmsf[0], 23, "palDd2tf", "(1)", status );
494 viv ( ihmsf[1], 42, "palDd2tf", "(2)", status );
495 viv ( ihmsf[2], 13, "palDd2tf", "(3)", status );
496 viv ( ihmsf[3], 3333, "palDd2tf", "(4)", status );
497}
498
499/* Calendar to MJD */
500
501static void t_cldj( int *status ) {
502 double d;
503 int j;
504
505 palCldj ( 1899, 12, 31, &d, &j );
506 vvd ( d, 15019, 0, "palCldj", "D", status );
507 viv ( j, 0, "palCldj", "J", status );
508}
509
510/* palDr2af */
511
512static void t_cr2af( int *status ) {
513 char s;
514 int idmsf[4];
515 palDr2af ( 4, 2.345, &s, idmsf );
516 viv ( s, '+', "palDr2af", "S", status );
517 viv ( idmsf[0], 134, "palDr2af", "(1)", status );
518 viv ( idmsf[1], 21, "palDr2af", "(2)", status );
519 viv ( idmsf[2], 30, "palDr2af", "(3)", status );
520 viv ( idmsf[3], 9706, "palDr2af", "(4)", status );
521}
522
523/* palDr2tf */
524
525static void t_cr2tf( int *status ) {
526 char s;
527 int ihmsf[4];
528 palDr2tf ( 4, -3.01234, &s, ihmsf );
529 viv ( s, '-', "palDr2tf", "S", status );
530 viv ( ihmsf[0], 11, "palDr2tf", "(1)", status );
531 viv ( ihmsf[1], 30, "palDr2tf", "(2)", status );
532 viv ( ihmsf[2], 22, "palDr2tf", "(3)", status );
533 viv ( ihmsf[3], 6484, "palDr2tf", "(4)", status );
534}
535
536/* palDtf2d */
537
538static void t_ctf2d( int *status ) {
539 double dd;
540 int j;
541
542 palDtf2d (23, 56, 59.1, &dd, &j);
543 vvd ( dd, 0.99790625, 1e-12, "palDtf2d", "D", status );
544 viv ( j, 0, "palDtf2d", "J", status );
545}
546
547/* palDtf2r */
548
549static void t_ctf2r( int *status ) {
550 double dr;
551 int j;
552
553 palDtf2r (23, 56, 59.1, &dr, &j);
554 vvd ( dr, 6.270029887942679, 1e-12, "palDtf2r",
555 "R", status );
556 viv ( j, 0, "palDtf2r", "J", status );
557}
558
559static void t_dat ( int *status ) {
560 vvd ( palDat ( 43900 ), 18, 0, "palDat",
561 " ", status );
562 vvd ( palDtt ( 40404 ), 39.709746, 1e-12, "palDtt",
563 " ", status );
564 vvd ( palDt ( 500 ), 4686.7, 1e-10, "palDt",
565 "500", status );
566 vvd ( palDt ( 1400 ), 408, 1e-11, "palDt",
567 "1400", status );
568 vvd ( palDt ( 1950 ), 27.99145626, 1e-12, "palDt",
569 "1950", status );
570}
571
572/* Dates */
573
574static void t_djcal( int *status ) {
575 const double djm = 50123.9999;
576 int iy, im, id;
577 int iydmf[4];
578 int j;
579 double f;
580
581 palDjcal ( 4, djm, iydmf, &j );
582 viv ( iydmf[0], 1996, "palDjcal", "Y", status );
583 viv ( iydmf[1], 2, "palDjcal", "M", status );
584 viv ( iydmf[2], 10, "palDjcal", "D", status );
585 viv ( iydmf[3], 9999, "palDjcal", "F", status );
586 viv ( j, 0, "palDjcal", "J", status );
587
588 palDjcl ( djm, &iy, &im, &id, &f, &j );
589 viv ( iy, 1996, "palDjcl", "Y", status );
590 viv ( im, 2, "palDjcl", "M", status );
591 viv ( id, 10, "palDjcl", "D", status );
592 vvd ( f, 0.9999, 1e-7, "palDjcl", "F", status );
593 viv ( j, 0, "palDjcl", "J", status );
594
595}
596
597/* Matrix inversion */
598
599static void t_dmat( int *status ) {
600 int j;
601 int iw[3];
602 double dd;
603 double da[9] = {
604 2.22, 1.6578, 1.380522,
605 1.6578, 1.380522, 1.22548578,
606 1.380522, 1.22548578, 1.1356276122
607 };
608 double dv[3] = {
609 2.28625, 1.7128825, 1.429432225
610 };
611
612 palDmat( 3, da, dv, &dd, &j, iw );
613 vvd ( da[0], 18.02550629769198,
614 1e-10, "palDmat", "a[0]", status );
615 vvd ( da[1], -52.16386644917280607,
616 1e-10, "palDmat", "a[1]", status );
617 vvd ( da[2], 34.37875949717850495,
618 1e-10, "palDmat", "a[2]", status );
619 vvd ( da[3], -52.16386644917280607,
620 1e-10, "palDmat", "a[3]", status );
621 vvd ( da[4], 168.1778099099805627,
622 1e-10, "palDmat", "a[4]", status );
623 vvd ( da[5], -118.0722869694232670,
624 1e-10, "palDmat", "a[5]", status );
625 vvd ( da[6], 34.37875949717850495,
626 1e-10, "palDmat", "a[6]", status );
627 vvd ( da[7], -118.0722869694232670,
628 1e-10, "palDmat", "a[7]", status );
629 vvd ( da[8], 86.50307003740151262,
630 1e-10, "palDmat", "a[8]", status );
631 vvd ( dv[0], 1.002346480763383,
632 1e-12, "palDmat", "v[0]", status );
633 vvd ( dv[1], 0.03285594016974583489,
634 1e-12, "palDmat", "v[1]", status );
635 vvd ( dv[2], 0.004760688414885247309,
636 1e-12, "palDmat", "v[2]", status );
637 vvd ( dd, 0.003658344147359863,
638 1e-12, "palDmat", "D", status );
639 viv ( j, 0, "palDmat", "J", status );
640
641}
642
643/* Test palDe2h and palDh2e routines */
644
645static void t_e2h( int *status ) {
646 double dh, dd, dp, da, de;
647
648 dh = -0.3;
649 dd = -1.1;
650 dp = -0.7;
651
652 palDe2h( dh, dd, dp, &da, &de );
653 vvd( da, 2.820087515852369, 1e-12, "palDe2h",
654 "AZ", status);
655 vvd( de, 1.132711866443304, 1e-12, "palDe2h",
656 "EL", status );
657
658 palDh2e( da, de, dp, &dh, &dd );
659 vvd( dh, -0.3, 1e-12, "palDh2e", "HA", status);
660 vvd( dd, -1.1, 1e-12, "palDh2e", "DEC", status );
661
662}
663
664/* Epochs */
665
666static void t_epb( int *status ) {
667 vvd ( palEpb( 45123 ), 1982.419793168669, 1e-8,
668 "palEpb", " ", status );
669}
670
671static void t_epb2d( int *status ) {
672 vvd ( palEpb2d( 1975.5 ), 42595.5995279655, 1e-7,
673 "palEpb2d", " ", status );
674}
675
676static void t_epco( int *status ) {
677 vvd ( palEpco ( 'B', 'J', 2000 ), 2000.001277513665,
678 1e-7, "palEpco", "BJ", status );
679 vvd ( palEpco ( 'J', 'B', 1950 ), 1949.999790442300,
680 1e-7, "palEpco", "JB", status );
681 vvd ( palEpco ( 'J', 'j', 2000 ), 2000,
682 1e-7, "palEpco", "JJ", status );
683}
684
685static void t_epj( int *status ) {
686 vvd ( palEpj( 42999 ), 1976.603696098563,
687 1e-7, "palEpj", " ", status );
688}
689
690static void t_epj2d( int *status ) {
691 vvd ( palEpj2d( 2010.077 ), 55225.124250,
692 1e-6, "palEpj2d", " ", status );
693}
694
695/* Equation of the equinoxes */
696
697/* Use SOFA/ERFA test because of change in precession model */
698static void t_eqeqx (int *status ) {
699 vvd ( palEqeqx( 53736. ), -0.8834195072043790156e-5,
700 1e-15, "palEqeqx", " ", status );
701}
702
703/* E-terms */
704
705static void t_etrms( int * status ) {
706 double ev[3];
707
708 palEtrms ( 1976.9, ev );
709
710 vvd ( ev[0], -1.621617102537041e-6, 1e-18, "palEtrms",
711 "X", status );
712 vvd ( ev[1], -3.310070088507914e-7, 1e-18, "palEtrms",
713 "Y", status );
714 vvd ( ev[2], -1.435296627515719e-7, 1e-18, "palEtrms",
715 "Z", status );
716}
717
718/* J2000 to Galactic */
719
720static void t_eqgal( int *status ) {
721 double dl, db;
722
723 palEqgal ( 5.67, -1.23, &dl, &db );
724
725 vvd ( dl, 5.612270780904526, 1e-12, "palEqgal",
726 "DL", status );
727 vvd ( db, -0.6800521449061520, 1e-12, "palEqgal",
728 "DB", status );
729}
730
731/* Galactic to J2000 equatorial */
732
733static void t_galeq( int *status ) {
734 double dr, dd;
735
736 palGaleq ( 5.67, -1.23, &dr, &dd );
737
738 vvd ( dr, 0.04729270418071426, 1e-12, "palGaleq",
739 "DR", status );
740 vvd ( dd, -0.7834003666745548, 1e-12, "palGaleq",
741 "DD", status );
742}
743
744/* Galactic to supergalactic */
745static void t_galsup(int *status ) {
746 double dsl, dsb;
747
748 palGalsup ( 6.1, -1.4, &dsl, &dsb );
749
750 vvd ( dsl, 4.567933268859171, 1e-12, "palGalsup",
751 "DSL", status );
752 vvd ( dsb, -0.01862369899731829, 1e-12, "palGalsup",
753 "DSB", status );
754}
755
756/* Geocentric coordinates */
757
758/* This is not from sla_test.f */
759
760static void t_geoc( int *status ) {
761 double r;
762 double z;
763 /* JCMT */
764 const double lat = 19.822838905884 * PAL__DD2R;
765 const double alt = 4120.0;
766 palGeoc( lat, alt, &r, &z );
767
768 /* Note the lower tolerance than normal since the models in SLA
769 differ from the more up to date model in SOFA/ERFA */
770 vvd( r, 4.01502667039618e-05, 1e-10, "palGeoc", "R", status );
771 vvd( z, 1.43762411970295e-05, 1e-10, "palGeoc", "Z", status );
772
773}
774
775/* Galactic to Fk4 */
776
777static void t_ge50 ( int *status ) {
778 double dr, dd;
779 palGe50( 6.1, -1.55, &dr, &dd );
780 vvd ( dr, 0.1966825219934508, 1e-12, "palGe50",
781 "DR", status );
782 vvd ( dd, -0.4924752701678960, 1e-12, "palGe50",
783 "DD", status );
784
785
786}
787
788/* GMST */
789
790/* We use the SOFA/ERFA test values rather than the values from SLA
791 because the precession models have changed */
792
793static void t_gmst( int *status ) {
794 vvd ( palGmst( 53736. ), 1.754174971870091203,
795 1e-12, "palGmst", " ", status );
796
797 vvd ( palGmsta( 53736., 0.0 ), 1.754174971870091203,
798 1e-12, "palGmsta", " ", status );
799}
800
801/* FK5 */
802
803static void t_fk52h ( int *status ) {
804 double r5, d5, dr5, dd5, rh, dh;
805
806 palFk5hz ( 1.234, -0.987, 1980, &rh, &dh );
807 vvd ( rh, 1.234000136713611301, 1e-13, "palFk5hz",
808 "R", status );
809 vvd ( dh, -0.9869999702020807601, 1e-13, "palFk5hz",
810 "D", status );
811 palHfk5z ( rh, dh, 1980, &r5, &d5, &dr5, &dd5 );
812 vvd ( r5, 1.234, 1e-13, "palHfk5z", "R", status );
813 vvd ( d5, -0.987, 1e-13, "palHfk5z", "D", status );
814 vvd ( dr5, 0.000000006822074, 1e-13, "palHfk5z",
815 "DR", status );
816 vvd ( dd5, -0.000000002334012, 1e-13, "palHfk5z",
817 "DD", status );
818
819}
820
821static void t_intin( int *status ) {
822 const char s[] = " -12345, , -0 2000 + ";
823 /* 1234567890123456789012345678 */
824 int i = 1;
825 long n = 0;
826 int j;
827
828 palIntin ( s, &i, &n, &j );
829 viv ( i, 10, "palIntin", "I1", status );
830 viv ( n, -12345, "palIntin", "V1", status );
831 viv ( j, -1, "palIntin", "J1", status );
832
833 palIntin ( s, &i, &n, &j );
834 viv ( i, 12, "palIntin", "I2", status );
835 viv ( n, -12345, "palIntin", "V2", status );
836 viv ( j, 1, "palIntin", "J2", status );
837
838 palIntin ( s, &i, &n, &j );
839 viv ( i, 17, "palIntin", "I3", status );
840 viv ( n, 0, "palIntin", "V3", status );
841 viv ( j, -1, "palIntin", "J3", status );
842
843 palIntin ( s, &i, &n, &j );
844 viv ( i, 23, "palIntin", "I4", status );
845 viv ( n, 2000, "palIntin", "V4", status );
846 viv ( j, 0, "palIntin", "J4", status );
847
848 palIntin ( s, &i, &n, &j );
849 viv ( i, 29, "palIntin", "I5", status );
850 viv ( n, 2000, "palIntin", "V5", status );
851 viv ( j, 1, "palIntin", "J5", status ); /* Note that strtol does not care about a + */
852
853}
854
855
856/* Moon */
857
858static void t_moon ( int *status ) {
859 double pv[6];
860
861 double expected1[] = {
862 0.00229161514616454,
863 0.000973912029208393,
864 0.000669931538978146,
865 -3.44709700068209e-09,
866 5.44477533462392e-09,
867 2.11785724844417e-09
868 };
869
870 /* SLA test only include slaMoon so we use the
871 example from SUN/67 */
872 palDmoon( 48634.4687174074, pv );
873 vvec( 6, pv, expected1, "palDmoon", status );
874}
875
876/* Nutation */
877
878static void t_nut( int *status ) {
879 double dpsi, deps, eps0;
880
881 double expected[3][3] = {
882 { 9.999999969492166e-1, 7.166577986249302e-5, 3.107382973077677e-5 },
883 { -7.166503970900504e-5, 9.999999971483732e-1, -2.381965032461830e-5 },
884 { -3.107553669598237e-5, 2.381742334472628e-5, 9.999999992335206818e-1 }
885 };
886
887 double rmatn[3][3];
888
889 /* SLA tests with low precision */
890 palNut( 46012.32, rmatn );
891 vrmat( rmatn, expected, "palNut", 1.0e-3, status );
892
893 /* Use the SOFA/ERFA tests */
894 palNutc( 54388.0, &dpsi, &deps, &eps0 );
895 vvd( eps0, 0.4090749229387258204, 1e-14,
896 "palNutc", "eps0", status);
897
898 palNutc( 53736.0, &dpsi, &deps, &eps0 );
899 vvd(dpsi, -0.9630912025820308797e-5, 1e-13,
900 "palNutc", "dpsi", status);
901 vvd(deps, 0.4063238496887249798e-4, 1e-13,
902 "palNutc", "deps", status);
903}
904
905/* palPrebn */
906
907static void t_prebn( int *status ) {
908 double rmatp[3][3];
909 double prebn_expected[3][3] = {
910 { 9.999257613786738e-1, -1.117444640880939e-2, -4.858341150654265e-3 },
911 { 1.117444639746558e-2, 9.999375635561940e-1, -2.714797892626396e-5 },
912 { 4.858341176745641e-3, -2.714330927085065e-5, 9.999881978224798e-1 },
913 };
914
915 palPrebn ( 1925., 1975., rmatp );
916 vrmat( rmatp, prebn_expected, "palPrebn", 1.0e-12, status );
917}
918
919/* Range */
920
921static void t_range( int *status ) {
922 vvd ( palDrange ( -4 ), 2.283185307179586,
923 1e-12, "palDrange", " ", status );
924}
925
926static void t_ranorm( int *status ) {
927 vvd ( palDranrm ( -0.1 ), 6.183185307179587,
928 1e-12, "palDranrm", "2", status );
929}
930
931/* Separation routines */
932
933static void t_sep( int *status ) {
934 double d1[3] = { 1.0, 0.1, 0.2 };
935 double d2[3] = { -3.0, 1e-3, 0.2 };
936 double ad1, bd1, ad2, bd2;
937
938 palDcc2s( d1, &ad1, &bd1 );
939 palDcc2s( d2, &ad2, &bd2 );
940
941 vvd ( palDsep ( ad1, bd1, ad2, bd2 ),
942 2.8603919190246608, 1e-7, "palDsep", " ", status );
943 vvd ( palDsepv ( d1, d2 ),
944 2.8603919190246608, 1e-7, "palDsepv", " ", status );
945
946}
947
948/* Supergalactic */
949
950static void t_supgal( int *status ) {
951 double dl, db;
952
953 palSupgal ( 6.1, -1.4, &dl, &db );
954
955 vvd ( dl, 3.798775860769474, 1e-12, "palSupgal",
956 "DL", status );
957 vvd ( db, -0.1397070490669407, 1e-12, "palSupgal",
958 "DB", status );
959
960}
961
962/* Test spherical tangent-plane-projection routines */
963static void t_tp( int *status ) {
964
965 int j;
966 double dr0, dd0, dr1, dd1, dx, dy, dr2, dd2, dr01,
967 dd01, dr02, dd02;
968
969 dr0 = 3.1;
970 dd0 = -0.9;
971 dr1 = dr0 + 0.2;
972 dd1 = dd0 - 0.1;
973 palDs2tp( dr1, dd1, dr0, dd0, &dx, &dy, &j );
974 vvd( dx, 0.1086112301590404, 1e-12, "palDs2tp",
975 "x", status );
976 vvd( dy, -0.1095506200711452, 1e-12, "palDs2tp",
977 "y", status );
978 viv( j, 0, "palDs2tp", "j", status );
979
980 palDtp2s( dx, dy, dr0, dd0, &dr2, &dd2 );
981 vvd( dr2 - dr1, 0., 1e-12, "palDtp2s", "r", status );
982 vvd( dd2 - dd1, 0., 1e-12, "palDtp2s", "d", status );
983
984 palDtps2c( dx, dy, dr2, dd2, &dr01, &dd01, &dr02, &dd02, &j );
985 vvd( dr01, 3.1, 1e-12, "palDtps2c", "r1", status);
986 vvd( dd01, -0.9, 1e-12, "palDtps2c", "d1", status);
987 vvd( dr02, 0.3584073464102072, 1e-12, "palDtps2c",
988 "r2", status);
989 vvd( dd02, -2.023361658234722, 1e-12, "palDtps2c",
990 "d2", status );
991 viv( j, 1, "palDtps2c", "n", status );
992
993}
994
995/* Test all the 3-vector and 3x3 matrix routines. */
996
997static void t_vecmat( int * status ) {
998 int i;
999
1000 /* palDav2m */
1001 double drm1[3][3];
1002 double dav[3] = { -0.123, 0.0987, 0.0654 };
1003 double dav2m_expected[3][3] = {
1004 { 0.9930075842721269, 0.05902743090199868, -0.1022335560329612 },
1005 { -0.07113807138648245, 0.9903204657727545, -0.1191836812279541 },
1006 { 0.09420887631983825, 0.1256229973879967, 0.9875948309655174 },
1007 };
1008
1009 /* palDeuler */
1010 double drm2[3][3];
1011 double deuler_expected[3][3] = {
1012 { -0.1681574770810878, 0.1981362273264315, 0.9656423242187410 },
1013 { -0.2285369373983370, 0.9450659587140423, -0.2337117924378156 },
1014 { -0.9589024617479674, -0.2599853247796050, -0.1136384607117296 } };
1015
1016 /* palDmxm */
1017 double drm[3][3];
1018 double dmxm_expected[3][3] = {
1019 { -0.09010460088585805, 0.3075993402463796, 0.9472400998581048 },
1020 { -0.3161868071070688, 0.8930686362478707, -0.3200848543149236 },
1021 { -0.9444083141897035, -0.3283459407855694, 0.01678926022795169 },
1022 };
1023
1024 /* palDcs2c et al */
1025 double dv1[3];
1026 double dv2[3];
1027 double dv3[3];
1028 double dv4[3];
1029 double dv5[3];
1030 double dv6[3];
1031 double dv7[3];
1032 double dvm;
1033
1034 /* palDav2m */
1035 palDav2m( dav, drm1 );
1036 vrmat( drm1, dav2m_expected, "palDav2m", 1.0e-12, status );
1037
1038 /* Test palDeuler */
1039 palDeuler( "YZY", 2.345, -0.333, 2.222, drm2 );
1040 vrmat( drm2, deuler_expected, "palDeuler", 1.0e-12, status );
1041
1042 /* palDmxm */
1043 palDmxm( drm2, drm1, drm );
1044 vrmat( drm, dmxm_expected, "palDmxm", 1.0e-12, status );
1045
1046 /* palDcs2c */
1047 palDcs2c( 3.0123, -0.999, dv1 );
1048 vvd ( dv1[0], -0.5366267667260525, 1e-12,
1049 "palDcs2c", "x", status );
1050 vvd ( dv1[1], 0.06977111097651444, 1e-12,
1051 "palDcs2c", "y", status );
1052 vvd ( dv1[2], -0.8409302618566215, 1e-12,
1053 "palDcs2c", "z", status );
1054
1055 /* palDmxv */
1056 palDmxv( drm1, dv1, dv2 );
1057 palDmxv( drm2, dv2, dv3 );
1058 vvd ( dv3[0], -0.7267487768696160, 1e-12,
1059 "palDmxv", "x", status );
1060 vvd ( dv3[1], 0.5011537352639822, 1e-12,
1061 "palDmxv", "y", status );
1062 vvd ( dv3[2], 0.4697671220397141, 1e-12,
1063 "palDmxv", "z", status );
1064
1065 /* palDimxv */
1066 palDimxv( drm, dv3, dv4 );
1067 vvd ( dv4[0], -0.5366267667260526, 1e-12,
1068 "palDimxv", "X", status );
1069 vvd ( dv4[1], 0.06977111097651445, 1e-12,
1070 "palDimxv", "Y", status );
1071 vvd ( dv4[2], -0.8409302618566215, 1e-12,
1072 "palDimxv", "Z", status );
1073
1074 /* palDm2av */
1075 palDm2av( drm, dv5 );
1076 vvd ( dv5[0], 0.006889040510209034, 1e-12,
1077 "palDm2av", "X", status );
1078 vvd ( dv5[1], -1.577473205461961, 1e-12,
1079 "palDm2av", "Y", status );
1080 vvd ( dv5[2], 0.5201843672856759, 1e-12,
1081 "palDm2av", "Z", status );
1082
1083 for (i=0; i<3; i++) {
1084 dv5[i] *= 1000.0;
1085 }
1086
1087 /* palDvn */
1088 palDvn( dv5, dv6, &dvm );
1089 vvd ( dv6[0], 0.004147420704640065, 1e-12,
1090 "palDvn", "X", status );
1091 vvd ( dv6[1], -0.9496888606842218, 1e-12,
1092 "palDvn", "Y", status );
1093 vvd ( dv6[2], 0.3131674740355448, 1e-12,
1094 "palDvn", "Z", status );
1095 vvd ( dvm, 1661.042127339937, 1e-9, "palDvn",
1096 "M", status );
1097
1098 vvd ( palDvdv ( dv6, dv1 ), -0.3318384698006295,
1099 1e-12, "palDvn", " ", status );
1100
1101 /* palDvxv */
1102 palDvxv( dv6, dv1, dv7 );
1103 vvd ( dv7[0], 0.7767720597123304, 1e-12,
1104 "palDvxv", "X", status );
1105 vvd ( dv7[1], -0.1645663574562769, 1e-12,
1106 "palDvxv", "Y", status );
1107 vvd ( dv7[2], -0.5093390925544726, 1e-12,
1108 "palDvxv", "Z", status );
1109
1110}
1111
1112static void t_ecleq( int *status ) {
1113 double dr;
1114 double dd;
1115 palEcleq( 1.234, -0.123, 43210.0, &dr, &dd );
1116 vvd( dr, 1.229910118208851, 1e-5, "palEcleq",
1117 "RA", status );
1118 vvd( dd, 0.2638461400411088, 1e-5, "palEcleq",
1119 "Dec", status );
1120}
1121
1122static void t_ecmat( int *status ) {
1123 double rmat[3][3];
1124 double expected[3][3] = {
1125 { 1.0, 0.0, 0.0 },
1126 { 0.0, 0.91749307789883549624, 0.3977517467060596168 },
1127 { 0.0, -0.3977517467060596168, 0.91749307789883549624 } };
1128
1129 palEcmat( 55966.46, rmat );
1130 vrmat( rmat, expected, "palEcmat", 1.0e-12, status );
1131}
1132
1133static void t_eqecl ( int *status ) {
1134 double dl, db;
1135 palEqecl ( 0.789, -0.123, 46555, &dl, &db );
1136
1137 /* Slight changes from SLA for 2006 precession/nutation */
1138 vvd ( dl, 0.7036566430349022, 1e-6, "palEqecl",
1139 "L", status );
1140 vvd ( db, -0.4036047164116848, 1e-6, "palEqecl",
1141 "B", status );
1142}
1143
1144static void t_prec( int *status ) {
1145 double rmat[3][3];
1146 double expected[3][3] = {
1147 { 0.9999856154510, -0.0049192906204, -0.0021376320580 },
1148 { 0.0049192906805, 0.9999879002027, -5.2297405698747e-06 },
1149 { 0.0021376319197, -5.2859681191735e-06, 0.9999977152483 } };
1150
1151 palPrec( 1990.0, 2012.0, rmat );
1152 vrmat( rmat, expected, "palPrec", 1.0e-12, status );
1153}
1154
1155static void t_preces( int *status ) {
1156 double ra;
1157 double dc;
1158 ra = 6.28;
1159 dc = -1.123;
1160 palPreces ( "FK4", 1925, 1950, &ra, &dc );
1161 vvd ( ra, 0.002403604864728447, 1e-12, "palPreces",
1162 "R", status );
1163 vvd ( dc, -1.120570643322045, 1e-12, "palPreces",
1164 "D", status );
1165
1166 /* This is the SLA test but PAL now uses the IAU 2006
1167 precession model so we need to loosen the comparison */
1168 ra = 0.0123;
1169 dc = 1.0987;
1170 palPreces ( "FK5", 2050, 1990, &ra, &dc );
1171 vvd ( ra, 6.282003602708382, 1e-6, "palPreces",
1172 "R", status );
1173 vvd ( dc, 1.092870326188383, 1e-6, "palPreces",
1174 "D", status );
1175
1176}
1177
1178static void t_evp( int *status ) {
1179 double dvb[3],dpb[3],dvh[3],dph[3];
1180 double vbex[3] = { 1.6957348127008098514e-07,
1181 -9.1093446116039685966e-08,
1182 -3.9528532243991863036e-08 };
1183 double pbex[3] = {-0.49771075259730546136,
1184 -0.80273812396332311359,
1185 -0.34851593942866060383 };
1186 double vhex[3] = { 1.6964379181455713805e-07,
1187 -9.1147224045727438391e-08,
1188 -3.9553158272334222497e-08 };
1189 double phex[3] = { -0.50169124421419830639,
1190 -0.80650980174901798492,
1191 -0.34997162028527262212 };
1192
1193 double vbex2[3] = {
1194 -0.0109187426811683,
1195 -0.0124652546173285,
1196 -0.0054047731809662
1197 };
1198 double pbex2[3] = {
1199 -0.7714104440491060,
1200 +0.5598412061824225,
1201 +0.2425996277722475
1202 };
1203 double vhex2[3] = {
1204 -0.0109189182414732,
1205 -0.0124718726844084,
1206 -0.0054075694180650
1207 };
1208 double phex2[3] = {
1209 -0.7757238809297653,
1210 +0.5598052241363390,
1211 +0.2426998466481708
1212 };
1213
1214 palEvp( 2010.0, 2012.0, dvb, dpb, dvh, dph );
1215
1216 vvec( 3, dvb, vbex, "palEvp", status );
1217 vvec( 3, dpb, pbex, "palEvp", status );
1218 vvec( 3, dvh, vhex, "palEvp", status );
1219 vvec( 3, dph, phex, "palEvp", status );
1220
1221 palEpv ( 53411.52501161, dph, dvh, dpb, dvb );
1222
1223 vvec( 3, dvb, vbex2, "palEpv", status );
1224 vvec( 3, dpb, pbex2, "palEpv", status );
1225 vvec( 3, dvh, vhex2, "palEpv", status );
1226 vvec( 3, dph, phex2, "palEpv", status );
1227
1228}
1229
1230static void t_map( int *status ) {
1231 double ra, da;
1232 palMap ( 6.123, -0.999, 1.23e-5, -0.987e-5,
1233 0.123, 32.1, 1999, 43210.9, &ra, &da );
1234
1235 /* These are the SLA tests but and they agree to 0.1 arcsec
1236 with PAL/SOFA/ERFA. We expect a slight difference from the change
1237 to nutation models. */
1238 vvd ( ra, 6.117130429775647, 1e-6, "palMap",
1239 "RA", status );
1240 vvd ( da, -1.000880769038632, 1e-8, "palMap",
1241 "DA", status );
1242}
1243
1244static void t_mappa( int *status ) {
1245 double amprms[21];
1246 double expected[21] = {1.9986310746064646082,
1247 -0.1728200754134739392,
1248 0.88745394651412767839,
1249 0.38472374350184274094,
1250 -0.17245634725219796679,
1251 0.90374808622520386159,
1252 0.3917884696321610738,
1253 2.0075929387510784968e-08,
1254 -9.9464149073251757597e-05,
1255 -1.6125306981057062306e-05,
1256 -6.9897255793245634435e-06,
1257 0.99999999489900059935,
1258 0.99999983777998024959,
1259 -0.00052248206600935195865,
1260 -0.00022683144398381763045,
1261 0.00052248547063364874764,
1262 0.99999986339269864022,
1263 1.4950491424992534218e-05,
1264 0.00022682360163333854623,
1265 -1.5069005133483779417e-05,
1266 0.99999997416198904698};
1267
1268 palMappa( 2010.0, 55927.0, amprms );
1269 vvec( 21, amprms, expected, "palMappa", status );
1270}
1271
1272static void t_mapqk( int *status ) {
1273 /* Test mapqk by taking the geocentric apparent positions of Arcturus
1274 as downloaded from aa.usno.mil/data/docs/geocentric.php and trying
1275 to calculate it from Arcturus' mean position, proper motion, parallax,
1276 and radial velocity */
1277
1278 double amprms[21];
1279 double ra_0, dec_0; /* mean position */
1280 double ra_app, dec_app; /* geocentric apparent position */
1281 double ra_test, dec_test;
1282 double px, pm_ra, pm_dec, v_rad;
1283 int j;
1284 int iposn;
1285
1286 /*
1287 The data below represents the position of Arcturus on
1288 JD (UT) 2457000.375 as reported by
1289 http://aa.usno.navy.mil/data/docs/geocentric.php
1290 */
1291
1292 const char radec_0[] = "14 15 39.67207 19 10 56.673";
1293 const char radec_app[] = "14 16 19.59 19 6 19.56";
1294
1295 iposn = 1;
1296 palDafin(radec_0, &iposn, &ra_0, &j);
1297 palDafin(radec_0, &iposn, &dec_0, &j);
1298 ra_0 *= 15.0;
1299
1300 pm_ra = -1.0939*PAL__DAS2R;
1301 pm_ra /= cos(dec_0);
1302 pm_dec = -2.00006*PAL__DAS2R;
1303 v_rad = -5.19;
1304 px = 0.08883*PAL__DAS2R;
1305
1306 palMappa(2000.0, 56999.87537249177, amprms); /* time is the TDB MJD calculated from
1307 a JD of 2457000.375 with astropy.time */
1308
1309 palMapqk(ra_0, dec_0, pm_ra, pm_dec, px, v_rad, amprms, &ra_test, &dec_test);
1310
1311 iposn = 1;
1312 palDafin(radec_app, &iposn, &ra_app, &j);
1313 palDafin(radec_app, &iposn, &dec_app, &j);
1314 ra_app *= 15.0;
1315
1316 /* find the angular distance from the known mean position
1317 to the calculated mean postion */
1318 double dd;
1319 dd = palDsep(ra_test, dec_test, ra_app, dec_app);
1320 dd *= PAL__DR2AS;
1321 vvd( dd, 0.0, 0.1, "palMapqk", "distance", status);
1322}
1323
1324static void t_mapqkz( int *status ) {
1325 double amprms[21], ra, da, ra_c, da_c;
1326
1327 /* Run inputs through mapqk with zero proper motion, parallax
1328 and radial velocity. Then run the same inputs through mapqkz.
1329 Verify that the results are the same */
1330 palMappa( 2010.0, 55927.0, amprms );
1331 palMapqk(1.234, -0.567, 0.0, 0.0, 0.0, 0.0, amprms, &ra_c, &da_c);
1332 palMapqkz( 1.234, -0.567, amprms, &ra, &da );
1333 vvd( ra, ra_c, 1.0E-12, "palMapqkz", "ra", status );
1334 vvd( da, da_c, 1.0E-12, "palMapqkz", "da", status );
1335}
1336
1337static void t_ampqk( int *status ) {
1338 double amprms[21], rm, dm;
1339 palMappa( 2010.0, 55927.0, amprms );
1340 palAmpqk( 1.234, -0.567, amprms, &rm, &dm );
1341 vvd( rm, 1.233512033578303857, 1.0E-12, "palAmpqk", "rm", status );
1342 vvd( dm, -0.56702909748530827549, 1.0E-12, "palAmpqk", "dm", status );
1343}
1344
1345static void t_fk45z( int *status ) {
1346 double r2000, d2000;
1347 palFk45z( 1.2, -0.3, 1960.0, &r2000, &d2000 );
1348 vvd( r2000, 1.2097812228966762227, 1.0E-12, "palFk45z", "r2000", status );
1349 vvd( d2000, -0.29826111711331398935, 1.0E-12, "palFk45z", "d2000", status );
1350}
1351
1352static void t_fk54z( int *status ) {
1353 double r1950, d1950, dr1950, dd1950;
1354 palFk54z( 1.2, -0.3, 1960.0, &r1950, &d1950, &dr1950, &dd1950 );
1355 vvd( r1950, 1.1902221805755279771, 1.0E-12, "palFk54z", "r1950", status );
1356 vvd( d1950, -0.30178317645793828472, 1.0E-12, "palFk54z", "d1950", status );
1357 vvd( dr1950, -1.7830874775952945507e-08, 1.0E-12, "palFk54z", "dr1950", status );
1358 vvd( dd1950, 7.196059425334821089e-09, 1.0E-12, "palFk54z", "dd1950", status );
1359}
1360
1361static void t_fk524( int *status ) {
1362 double r1950, d1950, dr1950, dd1950, p1950, v1950;
1363 palFk524(4.567, -1.23, -3e-5, 8e-6, 0.29,
1364 -35.0, &r1950, &d1950, &dr1950, &dd1950, &p1950, &v1950);
1365 vvd(r1950, 4.543778603272084, 1e-12, "palFk524", "r", status);
1366 vvd(d1950, -1.229642790187574, 1e-12, "palFk524", "d", status);
1367 vvd(dr1950, -2.957873121769244e-5, 1e-17, "palFk524", "dr", status);
1368 vvd(dd1950, 8.117725309659079e-6, 1e-17, "palFk524", "dd", status);
1369 vvd(p1950, 0.2898494999992917, 1e-12, "palFk524", "p", status);
1370 vvd(v1950, -35.026862824252680, 1e-11, "palFk524", "v", status);
1371}
1372
1373static void t_flotin( int * status ) {
1374
1375 int j;
1376 const char * s = " 12.345, , -0 1E3-4 2000 E ";
1377 /* 123456789012345678901234567890123 */
1378 int i = 1;
1379 double dv = 0.0;
1380
1381 palDfltin ( s, &i, &dv, &j );
1382 viv ( i, 10, "palDfltin", "I1", status );
1383 vvd ( dv, 12.345, 1e-12, "palDfltin", "V1", status );
1384 viv ( j, 0, "palDfltin", "J1", status );
1385
1386 palDfltin ( s, &i, &dv, &j );
1387 viv ( i, 12, "palDfltin", "I2", status );
1388 vvd ( dv, 12.345, 1e-12, "palDfltin", "V2", status );
1389 viv ( j, 1, "palDfltin", "J2", status );
1390
1391 palDfltin ( s, &i, &dv, &j );
1392 viv ( i, 16, "palDfltin", "I3", status );
1393 vvd ( dv, 0, 0, "palDfltin", "V3", status );
1394 viv ( j, -1, "palDfltin", "J3", status );
1395
1396 palDfltin ( s, &i, &dv, &j );
1397 viv ( i, 19, "palDfltin", "I4", status );
1398 vvd ( dv, 1000, 0, "palDfltin", "V4", status );
1399 viv ( j, 0, "palDfltin", "J4", status );
1400
1401 palDfltin ( s, &i, &dv, &j );
1402 viv ( i, 22, "palDfltin", "I5", status );
1403 vvd ( dv, -4, 0, "palDfltin", "V5", status );
1404 viv ( j, -1, "palDfltin", "J5", status );
1405
1406 palDfltin ( s, &i, &dv, &j );
1407 viv ( i, 28, "palDfltin", "I6", status );
1408 vvd ( dv, 2000, 0, "palDfltin", "V6", status );
1409 viv ( j, 0, "palDfltin", "J6", status );
1410
1411 palDfltin ( s, &i, &dv, &j );
1412 viv ( i, 34, "palDfltin", "I7", status );
1413 vvd ( dv, 2000, 0, "palDfltin", "V7", status );
1414 viv ( j, 1, "palDfltin", "J7", status ); /* differs from slaDfltin */
1415
1416 /* Now test overflow and underflow */
1417 i = 1;
1418 palDfltin( " 1D600 ", &i, &dv, &j );
1419 viv ( i, 8, "palDfltin", "I8", status );
1420 vvd ( dv, HUGE_VAL, 0, "palDfltin", "V8", status );
1421 viv ( j, 2, "palDfltin", "J8", status );
1422
1423}
1424
1425static void t_obs( int * status ) {
1426
1427 char shortname[11];
1428 char longname[41];
1429 double w, p, h;
1430 int lstat;
1431
1432 lstat = palObs( 0, "MMT", shortname, sizeof(shortname),
1433 longname, sizeof(longname), &w, &p, &h );
1434 vcs ( shortname, "MMT", "palObs", "1/C", status );
1435 vcs ( longname, "MMT 6.5m, Mt Hopkins", "palObs", "1/NAME",
1436 status );
1437 vvd ( w, 1.935300584055477, 1e-8, "palObs",
1438 "1/W", status );
1439 vvd ( p, 0.5530735081550342238, 1e-10, "palObs",
1440 "1/P", status );
1441 vvd ( h, 2608, 1e-10, "palObs",
1442 "1/H", status );
1443 viv( lstat, 0, "palObs", "retval", status );
1444
1445 lstat = palObs ( 61, NULL, shortname, sizeof(shortname),
1446 longname, sizeof(longname), &w, &p, &h );
1447 vcs ( shortname, "KECK1", "palObs", "2/C", status );
1448 vcs ( longname, "Keck 10m Telescope #1", "palObs",
1449 "2/NAME", status );
1450 vvd ( w, 2.713545757918895, 1e-8, "palObs",
1451 "2/W", status );
1452 vvd ( p, 0.3460280563536619, 1e-8, "palObs",
1453 "2/P", status );
1454 vvd ( h, 4160, 1e-10, "palObs",
1455 "2/H", status );
1456 viv( lstat, 0, "palObs", "retval", status );
1457
1458 lstat = palObs ( 83, NULL, shortname, sizeof(shortname),
1459 longname, sizeof(longname), &w, &p, &h );
1460 vcs ( shortname, "MAGELLAN2", "palObs", "3/C", status );
1461 vcs ( longname, "Magellan 2, 6.5m, Las Campanas",
1462 "palObs", "3/NAME", status );
1463 vvd ( w, 1.233819305534497, 1e-8, "palObs",
1464 "3/W", status );
1465 vvd ( p, -0.506389344359954, 1e-8, "palObs",
1466 "3/P", status );
1467 vvd ( h, 2408, 1e-10, "palObs",
1468 "3/H", status );
1469 viv( lstat, 0, "palObs", "retval", status );
1470
1471 /* the first argument here should be 1 greater than the number of items
1472 * in const struct telData defined in palObs.c
1473 */
1474 lstat = palObs ( 86, NULL, shortname, sizeof(shortname),
1475 longname, sizeof(longname), &w, &p, &h );
1476 vcs ( longname, "?", "palObs", "4/NAME", status );
1477 viv( lstat, -1, "palObs", "retval", status );
1478
1479 lstat = palObs ( 0, "MISSING", shortname, sizeof(shortname),
1480 longname, sizeof(longname), &w, &p, &h );
1481 vcs ( longname, "?", "palObs", "5/NAME", status );
1482 viv( lstat, -1, "palObs", "retval", status );
1483
1484 lstat = palObs( 0, "mmt", shortname, sizeof(shortname),
1485 longname, sizeof(longname), &w, &p, &h );
1486 vcs ( shortname, "MMT", "palObs", "6/C", status );
1487 vcs ( longname, "MMT 6.5m, Mt Hopkins", "palObs", "6/NAME",
1488 status );
1489 vvd ( w, 1.935300584055477, 1e-8, "palObs",
1490 "6/W", status );
1491 vvd ( p, 0.5530735081550342238, 1e-10, "palObs",
1492 "6/P", status );
1493 vvd ( h, 2608, 1e-10, "palObs",
1494 "6/H", status );
1495 viv( lstat, 0, "palObs", "retval", status );
1496
1497}
1498
1499static void t_pa( int *status ) {
1500 vvd ( palPa ( -1.567, 1.5123, 0.987 ),
1501 -1.486288540423851, 1e-12, "palPa", " ", status );
1502 vvd ( palPa ( 0, 0.789, 0.789 ),
1503 0, 0, "palPa", "zenith", status );
1504}
1505
1506static void t_pcd( int *status ) {
1507 double disco, x, y;
1508 disco = 178.585;
1509 x = 0.0123;
1510 y = -0.00987;
1511
1512 palPcd ( disco, &x, &y );
1513 vvd ( x, 0.01284630845735895, 1e-14, "palPcd", "x", status );
1514 vvd ( y, -0.01030837922553926, 1e-14, "palPcd", "y", status );
1515
1516 palUnpcd ( disco, &x, &y );
1517 vvd ( x, 0.0123, 1e-14, "palUnpcd", "x", status );
1518 vvd ( y, -0.00987, 1e-14, "palUnpcd", "y,", status );
1519
1520 /* Now negative disco round trip */
1521 disco = -0.3333333;
1522 x = 0.0123;
1523 y = -0.00987;
1524 palPcd ( disco, &x, &y );
1525 palUnpcd ( disco, &x, &y );
1526 vvd ( x, 0.0123, 1e-14, "palUnpcd", "x", status );
1527 vvd ( y, -0.00987, 1e-14, "palUnpcd", "y,", status );
1528}
1529
1530static void t_planet( int * status ) {
1531 int j;
1532 double pv[6];
1533 double u[13];
1534 double expected1[6] = { 0., 0., 0., 0., 0., 0. };
1535 double expectedue1[13] = {
1536 1.000878908362435284, -0.3336263027874777288, 50000.,
1537 2.840425801310305210, 0.1264380368035014224, -0.2287711835229143197,
1538 -0.01301062595106185195, 0.5657102158104651697, 0.2189745287281794885,
1539 2.852427310959998500, -0.01552349065435120900,
1540 50000., 0.0
1541 };
1542 double expectedue2[13] = {
1543 1.00006, -4.856142884511782, 50000., 0.3, -0.2,
1544 0.1, -0.4520378601821727, 0.4018114312730424,
1545 -.3515850023639121, 0.3741657386773941,
1546 -0.2511321445456515, 50000., 0.
1547 };
1548 double expectedue3[13] = {
1549 1.000000000000000,
1550 -0.3329769417028020949,
1551 50100.,
1552 2.638884303608524597,
1553 1.070994304747824305,
1554 0.1544112080167568589,
1555 -0.2188240619161439344,
1556 0.5207557453451906385,
1557 0.2217782439275216936,
1558 2.852118859689216658,
1559 0.01452010174371893229,
1560 50100.,
1561 0.
1562 };
1563 double expectedpv[6] = {
1564 0.07944764084631667011, -0.04118141077419014775,
1565 0.002915180702063625400, -0.6890132370721108608e-6,
1566 0.4326690733487621457e-6, -0.1763249096254134306e-6,
1567 };
1568 double expectedpv2[6] = {
1569 1.947628959288897677,
1570 -1.013736058752235271,
1571 -0.3536409947732733647,
1572 2.742247411571786194e-8,
1573 1.170467244079075911e-7,
1574 3.709878268217564005e-8
1575 };
1576
1577 double ra,dec,diam, r;
1578 int jform;
1579 double epoch, orbinc, anode, perih, aorq, e, aorl,
1580 dm;
1581
1582 /* palEl2ue */
1583 palEl2ue ( 50000, 1, 49000, 0.1, 2, 0.2,
1584 3, 0.05, 3, 0.003312, u, &j );
1585 vvec( 13, u, expectedue1, "palEl2ue", status );
1586 viv ( j, 0, "palEl2ue", "J", status );
1587
1588 /* palPertel */
1589 palPertel ( 2, 43000., 43200., 43000.,
1590 0.2, 3, 4, 5, 0.02, 6,
1591 &epoch, &orbinc, &anode, &perih, &aorq, &e, &aorl, &j );
1592 vvd ( epoch, 43200., 1e-10, "palPertel",
1593 "EPOCH", status );
1594 vvd ( orbinc, 0.1995661466545422381, 1e-7, "palPertel",
1595 "ORBINC", status );
1596 vvd ( anode, 2.998052737821591215, 1e-7, "palPertel",
1597 "ANODE", status );
1598 vvd ( perih, 4.009516448441143636, 1e-6, "palPertel",
1599 "PERIH", status );
1600 vvd ( aorq, 5.014216294790922323, 1e-7, "palPertel",
1601 "AORQ", status );
1602 vvd ( e, 0.02281386258309823607, 1e-7, "palPertel",
1603 "E", status );
1604 vvd ( aorl, 0.01735248648779583748, 1e-6, "palPertel",
1605 "AORL", status );
1606 viv ( j, 0, "palPertel", "J", status );
1607
1608 /* palPertue */
1609 palPertue ( 50100, u, &j );
1610 vvec( 13, u, expectedue3, "palPertue", status );
1611 viv ( j, 0, "palPertue", "J", status );
1612
1613 /* palPlanel */
1614 palPlanel ( 50600, 2, 50500, 0.1, 3, 5,
1615 2, 0.3, 4, 0, pv, &j );
1616 vvec( 6, pv, expectedpv2, "palPlanel", status );
1617 viv ( j, 0, "palPlanel", "J", status );
1618
1619 /* palPlanet */
1620
1621 palPlanet( 1e6, 0, pv, &j );
1622 vvec( 6, pv, expected1, "palPlanet 1", status );
1623 viv ( j, -1, "palPlanet", "J 1", status );
1624
1625 palPlanet( 1e6, 9, pv, &j);
1626 viv ( j, -1, "palPlanet", "J 2", status );
1627
1628 palPlanet ( -320000, 3, pv, &j );
1629 vvd ( pv[0], 0.9308038666827242603, 1e-11, "palPlanet",
1630 "pv[0] 3", status );
1631 vvd ( pv[1], 0.3258319040252137618, 1e-11, "palPlanet",
1632 "pv[1] 3", status );
1633 vvd ( pv[2], 0.1422794544477122021, 1e-11, "palPlanet",
1634 "pv[2] 3", status );
1635 vvd ( pv[3], -7.441503423889371696e-8, 1e-17, "palPlanet",
1636 "pv[3] 3", status );
1637 vvd ( pv[4], 1.699734557528650689e-7, 1e-17, "palPlanet",
1638 "pv[4] 3", status );
1639 vvd ( pv[5], 7.415505123001430864e-8, 1e-17, "palPlanet",
1640 "pv[5] 3", status );
1641 viv ( j, 1, "palPlanet", "J 3", status );
1642
1643 palPlanet ( 43999.9, 1, pv, &j );
1644 vvd ( pv[0], 0.2945293959257422246, 1e-11, "palPlanet",
1645 "pv[0] 4", status );
1646 vvd ( pv[1], -0.2452204176601052181, 1e-11, "palPlanet",
1647 "pv[1] 4", status );
1648 vvd ( pv[2], -0.1615427700571978643, 1e-11, "palPlanet",
1649 "pv[2] 4", status );
1650 vvd ( pv[3], 1.636421147459047057e-7, 1e-18, "palPlanet",
1651 "pv[3] 4", status );
1652 vvd ( pv[4], 2.252949422574889753e-7, 1e-18, "palPlanet",
1653 "pv[4] 4", status );
1654 vvd ( pv[5], 1.033542799062371839e-7, 1e-18, "palPlanet",
1655 "pv[5] 4", status );
1656 viv ( j, 0, "palPlanet", "J 4", status );
1657
1658 /* palPlante test would go here */
1659
1660 palPlante ( 50600., -1.23, 0.456, 2, 50500.,
1661 0.1, 3., 5., 2., 0.3, 4.,
1662 0., &ra, &dec, &r, &j );
1663 vvd ( ra, 6.222958101333794007, 1e-6, "palPlante",
1664 "RA", status );
1665 vvd ( dec, 0.01142220305739771601, 1e-6, "palPlante",
1666 "DEC", status );
1667 vvd ( r, 2.288902494080167624, 1e-8, "palPlante",
1668 "R", status );
1669 viv ( j, 0, "palPlante", "J", status );
1670
1671 u[0] = 1.0005;
1672 u[1] = -0.3;
1673 u[2] = 55000.;
1674 u[3] = 2.8;
1675 u[4] = 0.1;
1676 u[5] = -0.2;
1677 u[6] = -0.01;
1678 u[7] = 0.5;
1679 u[8] = 0.22;
1680 u[9] = 2.8;
1681 u[10] = -0.015;
1682 u[11] = 55001.;
1683 u[12] = 0;
1684
1685 /* palPlantu */
1686
1687 palPlantu ( 55001., -1.23, 0.456, u, &ra, &dec, &r, &j );
1688 vvd ( ra, 0.3531814831241686647, 1e-6, "palPlantu",
1689 "RA", status );
1690 vvd ( dec, 0.06940344580567131328, 1e-6, "palPlantu",
1691 "DEC", status );
1692 vvd ( r, 3.031687170873274464, 1e-8, "palPlantu",
1693 "R", status );
1694 viv ( j, 0, "palPlantu", "J", status );
1695
1696 /* palPv2el */
1697
1698 pv[0] = 0.3;
1699 pv[1] = -0.2;
1700 pv[2] = 0.1;
1701 pv[3] = -0.9e-7;
1702 pv[4] = 0.8e-7;
1703 pv[5] = -0.7e-7;
1704
1705 palPv2el ( pv, 50000, 0.00006, 1,
1706 &jform, &epoch, &orbinc, &anode, &perih,
1707 &aorq, &e, &aorl, &dm, &j );
1708 viv ( jform, 1, "palPv2el", "JFORM", status );
1709 vvd ( epoch, 50000, 1e-10, "palPv2el",
1710 "EPOCH", status );
1711 vvd ( orbinc, 1.52099895268912, 1e-12, "palPv2el",
1712 "ORBINC", status );
1713 vvd ( anode, 2.720503180538650, 1e-12, "palPv2el",
1714 "ANODE", status );
1715 vvd ( perih, 2.194081512031836, 1e-12, "palPv2el",
1716 "PERIH", status );
1717 vvd ( aorq, 0.2059371035373771, 1e-12, "palPv2el",
1718 "AORQ", status );
1719 vvd ( e, 0.9866822985810528, 1e-12, "palPv2el",
1720 "E", status );
1721 vvd ( aorl, 0.2012758344836794, 1e-12, "palPv2el",
1722 "AORL", status );
1723 vvd ( dm, 0.1840740507951820, 1e-12, "palPv2el",
1724 "DM", status );
1725 viv ( j, 0, "palPv2el", "J", status );
1726
1727 /* palPv2ue */
1728 palPv2ue ( pv, 50000., 0.00006, u, &j );
1729 vvec( 13, u, expectedue2, "palPv2ue", status );
1730 viv ( j, 0, "palPv2ue", "J", status );
1731
1732 /* Planets */
1733 palRdplan ( 40999.9, 0, 0.1, -0.9, &ra, &dec, &diam );
1734 vvd ( ra, 5.772270359389275837, 1e-6, "palRdplan",
1735 "ra 0", status );
1736 vvd ( dec, -0.2089207338795416192, 1e-7, "palRdplan",
1737 "dec 0", status );
1738 vvd ( diam, 9.415338935229717875e-3, 1e-10, "palRdplan",
1739 "diam 0", status );
1740 palRdplan ( 41999.9, 1, 1.1, -0.9, &ra, &dec, &diam );
1741 vvd ( ra, 3.866363420052936653, 1e-6, "palRdplan",
1742 "ra 1", status );
1743 vvd ( dec, -0.2594430577550113130, 1e-7, "palRdplan",
1744 "dec 1", status );
1745 vvd ( diam, 4.638468996795023071e-5, 1e-14, "palRdplan",
1746 "diam 1", status );
1747 palRdplan ( 42999.9, 2, 2.1, 0.9, &ra, &dec, &diam );
1748 vvd ( ra, 2.695383203184077378, 1e-6, "palRdplan",
1749 "ra 2", status );
1750 vvd ( dec, 0.2124044506294805126, 1e-7, "palRdplan",
1751 "dec 2", status );
1752 vvd ( diam, 4.892222838681000389e-5, 1e-14, "palRdplan",
1753 "diam 2", status );
1754 palRdplan ( 43999.9, 3, 3.1, 0.9, &ra, &dec, &diam );
1755 vvd ( ra, 2.908326678461540165, 1e-7, "palRdplan",
1756 "ra 3", status );
1757 vvd ( dec, 0.08729783126905579385, 1e-7, "palRdplan",
1758 "dec 3", status );
1759 vvd ( diam, 8.581305866034962476e-3, 1e-7, "palRdplan",
1760 "diam 3", status );
1761 palRdplan ( 44999.9, 4, -0.1, 1.1, &ra, &dec, &diam );
1762 vvd ( ra, 3.429840787472851721, 1e-6, "palRdplan",
1763 "ra 4", status );
1764 vvd ( dec, -0.06979851055261161013, 1e-7, "palRdplan",
1765 "dec 4", status );
1766 vvd ( diam, 4.540536678439300199e-5, 1e-14, "palRdplan",
1767 "diam 4", status );
1768 palRdplan ( 45999.9, 5, -1.1, 0.1, &ra, &dec, &diam );
1769 vvd ( ra, 4.864669466449422548, 1e-6, "palRdplan",
1770 "ra 5", status );
1771 vvd ( dec, -0.4077714497908953354, 1e-7, "palRdplan",
1772 "dec 5", status );
1773 vvd ( diam, 1.727945579027815576e-4, 1e-14, "palRdplan",
1774 "diam 5", status );
1775 palRdplan ( 46999.9, 6, -2.1, -0.1, &ra, &dec, &diam );
1776 vvd ( ra, 4.432929829176388766, 1e-6, "palRdplan",
1777 "ra 6", status );
1778 vvd ( dec, -0.3682820877854730530, 1e-7, "palRdplan",
1779 "dec 6", status );
1780 vvd ( diam, 8.670829016099083311e-5, 1e-14, "palRdplan",
1781 "diam 6", status );
1782 palRdplan ( 47999.9, 7, -3.1, -1.1, &ra, &dec, &diam );
1783 vvd ( ra, 4.894972492286818487, 1e-6, "palRdplan",
1784 "ra 7", status );
1785 vvd ( dec, -0.4084068901053653125, 1e-7, "palRdplan",
1786 "dec 7", status );
1787 vvd ( diam, 1.793916783975974163e-5, 1e-14, "palRdplan",
1788 "diam 7", status );
1789 palRdplan ( 48999.9, 8, 0, 0, &ra, &dec, &diam );
1790 vvd ( ra, 5.066050284760144000, 1e-6, "palRdplan",
1791 "ra 8", status );
1792 vvd ( dec, -0.3744690779683850609, 1e-7, "palRdplan",
1793 "dec 8", status );
1794 vvd ( diam, 1.062210086082700563e-5, 1e-14, "palRdplan",
1795 "diam 8", status );
1796
1797 /* palUe2el */
1798 palUe2el ( u, 1, &jform, &epoch, &orbinc, &anode, &perih,
1799 &aorq, &e, &aorl, &dm, &j );
1800 viv ( jform, 1, "palUe2el", "JFORM", status );
1801 vvd ( epoch, 50000.00000000000, 1e-10, "palUe2el",
1802 "EPOCH", status );
1803 vvd ( orbinc, 1.520998952689120, 1e-12, "palUe2el",
1804 "ORBINC", status );
1805 vvd ( anode, 2.720503180538650, 1e-12, "palUe2el",
1806 "ANODE", status );
1807 vvd ( perih, 2.194081512031836, 1e-12, "palUe2el",
1808 "PERIH", status );
1809 vvd ( aorq, 0.2059371035373771, 1e-12, "palUe2el",
1810 "AORQ", status );
1811 vvd ( e, 0.9866822985810528, 1e-12, "palUe2el",
1812 "E", status );
1813 vvd ( aorl, 0.2012758344836794, 1e-12, "palUe2el",
1814 "AORL", status );
1815 viv ( j, 0, "palUe2el", "J", status );
1816
1817 /* palUe2pv */
1818 palUe2pv( 50010., u, pv, &j );
1819
1820 /* Update the final two elements of the expecte UE array */
1821 expectedue2[11] = 50010.;
1822 expectedue2[12] = 0.7194308220038886856;
1823
1824 vvec( 13, u, expectedue2, "palUe2pv", status );
1825 vvec( 6, pv, expectedpv, "palUe2pv", status );
1826 viv ( j, 0, "palUe2pv", "J", status );
1827
1828}
1829
1830static void t_pm( int * status ) {
1831 double ra2, dec2;
1832 double ra1, dec1, pmr1, pmd1, px1, rv1;
1833
1834 ra1 = 5.43;
1835 dec1 = -0.87;
1836 pmr1 = -0.33e-5;
1837 pmd1 = 0.77e-5;
1838 px1 = 0.7;
1839 rv1 = 50.3*365.2422/365.25;
1840
1841 palPm ( ra1, dec1, pmr1, pmd1, px1, rv1,
1842 1899, 1943,
1843 &ra2, &dec2 );
1844 vvd ( ra2, 5.429855087793875, 1e-10, "palPm",
1845 "R", status );
1846 vvd ( dec2, -0.8696617307805072, 1e-10, "palPm",
1847 "D", status );
1848
1849 /* SOFA/ERFA test */
1850 ra1 = 0.01686756;
1851 dec1 = -1.093989828;
1852 pmr1 = -1.78323516e-5;
1853 pmd1 = 2.336024047e-6;
1854 px1 = 0.74723;
1855 rv1 = -21.6;
1856
1857 palPm(ra1, dec1, pmr1, pmd1, px1, rv1,
1858 palEpj(50083.0), palEpj(53736.0),
1859 &ra2, &dec2);
1860 vvd(ra2, 0.01668919069414242368, 1e-13,
1861 "palPm", "ra", status);
1862 vvd(dec2, -1.093966454217127879, 1e-13,
1863 "palPm", "dec", status);
1864
1865
1866}
1867
1868static void t_polmo( int *status ) {
1869 double elong, phi, daz;
1870
1871 palPolmo( 0.7, -0.5, 1.0e-6, -2.0e-6, &elong, &phi, &daz );
1872 vvd(elong, 0.7000004837322044, 1.0e-12, "palPolmo", "elong", status );
1873 vvd(phi, -0.4999979467222241, 1.0e-12, "palPolmo", "phi", status );
1874 vvd(daz, 1.008982781275728e-6, 1.0e-12, "palPolmo", "daz", status );
1875}
1876
1877static void t_pvobs( int *status ) {
1878 double pv[6];
1879 double expected[6] = { -4.7683600138836167813e-06,
1880 1.0419056712717953176e-05,
1881 4.099831053320363277e-05,
1882 -7.5976959740661272483e-10,
1883 -3.4771429582640930371e-10,
1884 0.0};
1885 palPvobs( 1.3, 10000.0, 2.0, pv );
1886 vvec( 6, pv, expected, "palPvobs", status );
1887}
1888
1889static void t_rv( int *status ) {
1890 vvd ( palRverot ( -0.777, 5.67, -0.3, 3.19 ),
1891 -0.1948098355075913, 1e-6,
1892 "palRverot", " ", status );
1893 vvd ( palRvgalc ( 1.11E0, -0.99E0 ),
1894 158.9630759840254, 1e-3, "palRvgalc", " ", status );
1895 vvd ( palRvlg ( 3.97E0, 1.09E0 ),
1896 -197.818762175363, 1e-3, "palRvlg", " ", status );
1897 vvd ( palRvlsrd ( 6.01E0, 0.1E0 ),
1898 -4.082811335150567, 1e-4, "palRvlsrd", " ", status );
1899 vvd ( palRvlsrk ( 6.01E0, 0.1E0 ),
1900 -5.925180579830265, 1e-4, "palRvlsrk", " ", status );
1901}
1902
1903static void t_rvgalc( int *status ) {
1904 double rv;
1905 rv = palRvgalc( 2.7, -1.0 );
1906 vvd( rv, 213.98084425751144977, 1.0E-12, "palRvgalc", "rv", status );
1907}
1908
1909static void t_rvlg( int *status ) {
1910 double rv;
1911 rv = palRvlg( 2.7, -1.0 );
1912 vvd( rv, 291.79205281252404802, 1.0E-12, "palRvlg", "rv", status );
1913}
1914
1915static void t_rvlsrd( int *status ) {
1916 double rv;
1917 rv = palRvlsrd( 2.7, -1.0 );
1918 vvd( rv, 9.620674692097630043, 1.0E-12, "palRvlsrd", "rv", status );
1919}
1920
1921static void t_rvlsrk( int *status ) {
1922 double rv;
1923 rv = palRvlsrk( 2.7, -1.0 );
1924 vvd( rv, 12.556356851411955233, 1.0E-12, "palRvlsrk", "rv", status );
1925}
1926
1927static void t_refco( int *status ) {
1928 double phpa, tc, rh, wl, refa, refb;
1929 phpa = 800.0;
1930 tc = 10.0 + 273.15; /* SLA uses kelvin */
1931 rh = 0.9;
1932 wl = 0.4;
1933 palRefcoq(tc, phpa, rh, wl, &refa, &refb);
1934 vvd(refa, 0.2264949956241415009e-3, 1e-15,
1935 "palRefcoq", "refa", status);
1936 vvd(refb, -0.2598658261729343970e-6, 1e-18,
1937 "palRefcoq", "refb", status);
1938}
1939
1940static void t_ref( int *status ) {
1941 double ref, refa, refb, refa2, refb2, vu[3], vr[3], zr;
1942
1943 palRefro( 1.4, 3456.7, 280, 678.9, 0.9, 0.55,
1944 -0.3, 0.006, 1e-9, &ref );
1945 vvd( ref, 0.00106715763018568, 1e-12, "palRefro",
1946 "o", status );
1947
1948 palRefro( 1.4, 3456.7, 280, 678.9, 0.9, 1000,
1949 -0.3, 0.006, 1e-9, &ref );
1950 vvd( ref, 0.001296416185295403, 1e-12, "palRefro",
1951 "r", status );
1952
1953 palRefcoq( 275.9, 709.3, 0.9, 101, &refa, &refb );
1954 vvd( refa, 2.324736903790639e-4, 1e-12, "palRefcoq",
1955 "a/r", status );
1956 vvd( refb, -2.442884551059e-7, 1e-15, "palRefcoq",
1957 "b/r", status );
1958
1959 palRefco( 2111.1, 275.9, 709.3, 0.9, 101,
1960 -1.03, 0.0067, 1e-12, &refa, &refb );
1961 vvd( refa, 2.324673985217244e-4, 1e-12, "palRefco",
1962 "a/r", status );
1963 vvd( refb, -2.265040682496e-7, 1e-15, "palRefco",
1964 "b/r", status );
1965
1966 palRefcoq( 275.9, 709.3, 0.9, 0.77, &refa, &refb );
1967 vvd( refa, 2.007406521596588e-4, 1e-12, "palRefcoq",
1968 "a", status );
1969 vvd( refb, -2.264210092590e-7, 1e-15, "palRefcoq",
1970 "b", status );
1971
1972 palRefco( 2111.1, 275.9, 709.3, 0.9, 0.77,
1973 -1.03, 0.0067, 1e-12, &refa, &refb );
1974 vvd( refa, 2.007202720084551e-4, 1e-12, "palRefco",
1975 "a", status );
1976 vvd( refb, -2.223037748876e-7, 1e-15, "palRefco",
1977 "b", status );
1978
1979 palAtmdsp ( 275.9, 709.3, 0.9, 0.77,
1980 refa, refb, 0.5, &refa2, &refb2 );
1981 vvd ( refa2, 2.034523658888048e-4, 1e-12, "palAtmdsp",
1982 "a", status );
1983 vvd ( refb2, -2.250855362179e-7, 1e-15, "palAtmdsp",
1984 "b", status );
1985
1986 palDcs2c ( 0.345, 0.456, vu );
1987 palRefv ( vu, refa, refb, vr );
1988 vvd ( vr[0], 0.8447487047790478, 1e-12, "palRefv",
1989 "x1", status );
1990 vvd ( vr[1], 0.3035794890562339, 1e-12, "palRefv",
1991 "y1", status );
1992 vvd ( vr[2], 0.4407256738589851, 1e-12, "palRefv",
1993 "z1", status );
1994
1995 palDcs2c ( 3.7, 0.03, vu );
1996 palRefv ( vu, refa, refb, vr );
1997 vvd ( vr[0], -0.8476187691681673, 1e-12, "palRefv",
1998 "x2", status );
1999 vvd ( vr[1], -0.5295354802804889, 1e-12, "palRefv",
2000 "y2", status );
2001 vvd ( vr[2], 0.0322914582168426, 1e-12, "palRefv",
2002 "z2", status );
2003
2004 palRefz ( 0.567, refa, refb, &zr );
2005 vvd ( zr, 0.566872285910534, 1e-12, "palRefz",
2006 "hi el", status );
2007
2008 palRefz ( 1.55, refa, refb, &zr );
2009 vvd ( zr, 1.545697350690958, 1e-12, "palRefz",
2010 "lo el", status );
2011
2012
2013
2014}
2015
2016static void t_vers( int *status ) {
2017 char verstring[32];
2018
2019 int ver = palVers( verstring, sizeof(verstring));
2020 printf("PAL Version %s (%d)\n", verstring, ver);
2021 if ( ver < 6000 ) {
2022 *status = 1; /* palVers introduced at v0.6.0 */
2023 }
2024}
2025
2026/**********************************************************************/
2027
2028int main (void) {
2029
2030 /* Use the SLA and SOFA/ERFA conventions */
2031 int status = 0; /* Unix and SAE convention */
2032
2033 t_addet(&status);
2034 t_afin(&status);
2035 t_altaz(&status);
2036 t_ampqk(&status);
2037 t_aop(&status);
2038 t_airmas(&status);
2039 t_amp(&status);
2040 t_bear(&status);
2041 t_caf2r(&status);
2042 t_caldj(&status);
2043 t_cc2s(&status);
2044 t_cd2tf(&status);
2045 t_cldj(&status);
2046 t_cr2af(&status);
2047 t_cr2tf(&status);
2048 t_ctf2d(&status);
2049 t_ctf2r(&status);
2050 t_dat(&status);
2051 t_djcal(&status);
2052 t_dmat(&status);
2053 t_epb(&status);
2054 t_epb2d(&status);
2055 t_epco(&status);
2056 t_epj(&status);
2057 t_epj2d(&status);
2058 t_eqecl(&status);
2059 t_eqeqx(&status);
2060 t_etrms(&status);
2061 t_eqgal(&status);
2062 t_evp(&status);
2063 t_fk45z(&status);
2064 t_fk54z(&status);
2065 t_fk524(&status);
2066 t_flotin(&status);
2067 t_galeq(&status);
2068 t_galsup(&status);
2069 t_ge50(&status);
2070 t_geoc(&status);
2071 t_gmst(&status);
2072 t_fk52h(&status);
2073 t_intin(&status);
2074 t_prec(&status);
2075 t_preces(&status);
2076 t_ecleq(&status);
2077 t_ecmat(&status);
2078 t_e2h(&status);
2079 t_map(&status);
2080 t_mappa(&status);
2081 t_mapqk(&status);
2082 t_mapqkz(&status);
2083 t_moon(&status);
2084 t_nut(&status);
2085 t_obs(&status);
2086 t_pa(&status);
2087 t_pcd(&status);
2088 t_planet(&status);
2089 t_pm(&status);
2090 t_polmo(&status);
2091 t_prebn(&status);
2092 t_pvobs(&status);
2093 t_range(&status);
2094 t_ranorm(&status);
2095 t_ref(&status);
2096 t_refco(&status);
2097 t_rv(&status);
2098 t_rvgalc(&status);
2099 t_rvlg(&status);
2100 t_rvlsrd(&status);
2101 t_rvlsrk(&status);
2102 t_sep(&status);
2103 t_supgal(&status);
2104 t_tp(&status);
2105 t_vecmat(&status);
2106 t_vers(&status);
2107 return status;
2108}
Note: See TracBrowser for help on using the repository browser.