source: trunk/FACT++/pal/palTest.c@ 19506

Last change on this file since 19506 was 18712, checked in by tbretz, 10 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.