1 /*
2 *+
3 *  Name:
4 *     palOne2One
5 
6 *  Purpose:
7 *     File containing simple PAL wrappers for SLA routines that are identical in SOFA
8 
9 *  Invocation:
10 *     Matches SLA API
11 
12 *  Language:
13 *     Starlink ANSI C
14 
15 *  Type of Module:
16 *     Library routine
17 
18 *  Description:
19 *     Some SOFA/ERFA routines are identical to their SLA counterparts. PAL provides
20 *     direct counterparts to these although it is generally a better idea to
21 *     use the SOFA/ERFA routine directly in new code.
22 *
23 *     The PAL routines with direct equivalents in SOFA/ERFA are:
24 *     - palCldj
25 *     - palDbear
26 *     - palDaf2r
27 *     - palDav2m
28 *     - palDcc2s
29 *     - palDcs2c
30 *     - palDd2tf
31 *     - palDimxv
32 *     - palDm2av
33 *     - palDjcl
34 *     - palDmxm
35 *     - palDmxv
36 *     - palDpav
37 *     - palDr2af
38 *     - palDr2tf
39 *     - palDranrm
40 *     - palDsep
41 *     - palDsepv
42 *     - palDtf2d
43 *     - palDtf2r
44 *     - palDvdv
45 *     - palDvn
46 *     - palDvxv
47 *     - palEpb
48 *     - palEpb2d
49 *     - palEpj
50 *     - palEpj2d
51 *     - palEqeqx
52 *     - palFk5hz
53 *     - palGmst
54 *     - palGmsta
55 *     - palHfk5z
56 *     - palRefcoq
57 
58 *  Authors:
59 *     TIMJ: Tim Jenness (JAC, Hawaii)
60 *     DSB: David S Berry (JAC, Hawaii)
61 *     {enter_new_authors_here}
62 
63 *  Notes:
64 *     - Do not call these functions from other PAL functions. Always use
65 *       the SOFA/ERFA routines directly in new code.
66 *     - These are implemented as real functions rather than C preprocessor
67 *       macros so there may be a performance penalty in using the PAL
68 *       version instead of the SOFA/ERFA version.
69 *     - Routines that take MJDs have SOFA/ERFA equivalents that have an explicit
70 *       MJD offset included.
71 *     - palEqeqx, palGmst and palGmsta use the IAU 2006 precession model.
72 
73 *  History:
74 *     2012-02-10 (TIMJ):
75 *        Initial version
76 *        Adapted with permission from the Fortran SLALIB library.
77 *     2012-03-23 (TIMJ):
78 *        Update prologue.
79 *     2012-05-09 (DSBJ):
80 *        Move palDrange into a separate file.
81 *     2014-07-15 (TIMJ):
82 *        SOFA now has palRefcoq equivalent.
83 *     {enter_further_changes_here}
84 
85 *  Copyright:
86 *     Copyright (C) 2014 Tim Jenness
87 *     Copyright (C) 2012 Science and Technology Facilities Council.
88 *     All Rights Reserved.
89 
90 *  Licence:
91 *     This program is free software: you can redistribute it and/or
92 *     modify it under the terms of the GNU Lesser General Public
93 *     License as published by the Free Software Foundation, either
94 *     version 3 of the License, or (at your option) any later
95 *     version.
96 *
97 *     This program is distributed in the hope that it will be useful,
98 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
99 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
100 *     GNU Lesser General Public License for more details.
101 *
102 *     You should have received a copy of the GNU Lesser General
103 *     License along with this program.  If not, see
104 *     <http://www.gnu.org/licenses/>.
105 
106 *  Bugs:
107 *     {note_any_bugs_here}
108 *-
109 */
110 
111 #include "pal.h"
112 #include "palmac.h"
113 #include "pal1sofa.h"
114 
115 /*
116 *+
117 *  Name:
118 *     palCldj
119 
120 *  Purpose:
121 *     Gregorian Calendar to Modified Julian Date
122 
123 *  Language:
124 *     Starlink ANSI C
125 
126 *  Type of Module:
127 *     Library routine
128 
129 *  Invocation:
130 *     palCldj( int iy, int im, int id, double *djm, int *j );
131 
132 *  Arguments:
133 *     iy = int (Given)
134 *        Year in Gregorian calendar
135 *     im = int (Given)
136 *        Month in Gregorian calendar
137 *     id = int (Given)
138 *        Day in Gregorian calendar
139 *     djm = double * (Returned)
140 *        Modified Julian Date (JD-2400000.5) for 0 hrs
141 *     j = int * (Returned)
142 *        status: 0 = OK, 1 = bad year (MJD not computed),
143 *        2 = bad month (MJD not computed), 3 = bad day (MJD computed).
144 
145 *  Description:
146 *     Gregorian calendar to Modified Julian Date.
147 
148 *  Notes:
149 *     - Uses eraCal2jd(). See SOFA/ERFA documentation for details.
150 
151 *-
152 */
153 
palCldj(int iy,int im,int id,double * djm,int * j)154 void palCldj ( int iy, int im, int id, double *djm, int *j ) {
155   double djm0;
156   *j = eraCal2jd( iy, im, id, &djm0, djm );
157 }
158 
159 /*
160 *+
161 *  Name:
162 *     palDbear
163 
164 *  Purpose:
165 *     Bearing (position angle) of one point on a sphere relative to another
166 
167 *  Language:
168 *     Starlink ANSI C
169 
170 *  Type of Module:
171 *     Library routine
172 
173 *  Invocation:
174 *     pa = palDbear( double a1, double b1, double a2, double b2 );
175 
176 *  Arguments:
177 *     a1 = double (Given)
178 *        Longitude of point 1 (e.g. RA) in radians.
179 *     b1 = double (Given)
180 *        Latitude of point 1 (e.g. Dec) in radians.
181 *     a2 = double (Given)
182 *        Longitude of point 2 in radians.
183 *     b2 = double (Given)
184 *        Latitude of point 2 in radians.
185 
186 *  Returned Value:
187 *     The result is the bearing (position angle), in radians, of point
188 *     A2,B2 as seen from point A1,B1.  It is in the range +/- pi.  If
189 *     A2,B2 is due east of A1,B1 the bearing is +pi/2.  Zero is returned
190 *     if the two points are coincident.
191 
192 *  Description:
193 *     Bearing (position angle) of one point in a sphere relative to another.
194 
195 *  Notes:
196 *     - Uses eraPas(). See SOFA/ERFA documentation for details.
197 
198 *-
199 */
200 
palDbear(double a1,double b1,double a2,double b2)201 double palDbear ( double a1, double b1, double a2, double b2 ) {
202   return eraPas( a1, b1, a2, b2 );
203 }
204 
205 /*
206 *+
207 *  Name:
208 *     palDaf2r
209 
210 *  Purpose:
211 *     Convert degrees, arcminutes, arcseconds to radians
212 
213 *  Language:
214 *     Starlink ANSI C
215 
216 *  Type of Module:
217 *     Library routine
218 
219 *  Invocation:
220 *     palDaf2r( int ideg, int iamin, double asec, double *rad, int *j );
221 
222 *  Arguments:
223 *     ideg = int (Given)
224 *        Degrees.
225 *     iamin = int (Given)
226 *        Arcminutes.
227 *     iasec = double (Given)
228 *        Arcseconds.
229 *     rad = double * (Returned)
230 *        Angle in radians.
231 *     j = int * (Returned)
232 *        Status: 0 = OK, 1 = "ideg" out of range 0-359,
233 *                2 = "iamin" outside of range 0-59,
234 *                2 = "asec" outside range 0-59.99999
235 
236 *  Description:
237 *     Convert degrees, arcminutes, arcseconds to radians.
238 
239 *  Notes:
240 *     - Uses eraAf2a(). See SOFA/ERFA documentation for details.
241 
242 *-
243 */
244 
245 /* Arguments differ slightly. Assumes that the sign is always positive
246    and dealt with externally. */
palDaf2r(int ideg,int iamin,double asec,double * rad,int * j)247 void palDaf2r ( int ideg, int iamin, double asec, double *rad, int *j ) {
248   *j = eraAf2a( ' ', ideg, iamin, asec, rad );
249 }
250 
251 /*
252 *+
253 *  Name:
254 *     palDav2m
255 
256 *  Purpose:
257 *     Form the rotation matrix corresponding to a given axial vector.
258 
259 *  Language:
260 *     Starlink ANSI C
261 
262 *  Type of Module:
263 *     Library routine
264 
265 *  Invocation:
266 *     palDav2m( double axvec[3], double rmat[3][3] );
267 
268 *  Arguments:
269 *     axvec = double [3] (Given)
270 *       Axial vector (radians)
271 *     rmat = double [3][3] (Returned)
272 *       Rotation matrix.
273 
274 *  Description:
275 *     A rotation matrix describes a rotation about some arbitrary axis,
276 *     called the Euler axis.  The "axial vector" supplied to this routine
277 *     has the same direction as the Euler axis, and its magnitude is the
278 *     amount of rotation in radians.
279 
280 *  Notes:
281 *     - Uses eraRv2m(). See SOFA/ERFA documentation for details.
282 
283 *-
284 */
285 
palDav2m(double axvec[3],double rmat[3][3])286 void palDav2m ( double axvec[3], double rmat[3][3] ) {
287   eraRv2m( axvec, rmat );
288 }
289 
290 /*
291 *+
292 *  Name:
293 *     palDcc2s
294 
295 *  Purpose:
296 *     Cartesian to spherical coordinates
297 
298 *  Language:
299 *     Starlink ANSI C
300 
301 *  Type of Module:
302 *     Library routine
303 
304 *  Invocation:
305 *     palDcc2s( double v[3], double *a, double *b );
306 
307 *  Arguments:
308 *     v = double [3] (Given)
309 *        x, y, z vector.
310 *     a = double * (Returned)
311 *        Spherical coordinate (radians)
312 *     b = double * (Returned)
313 *        Spherical coordinate (radians)
314 
315 *  Description:
316 *     The spherical coordinates are longitude (+ve anticlockwise looking
317 *     from the +ve latitude pole) and latitude.  The Cartesian coordinates
318 *     are right handed, with the x axis at zero longitude and latitude, and
319 *     the z axis at the +ve latitude pole.
320 
321 *  Notes:
322 *     - Uses eraC2s(). See SOFA/ERFA documentation for details.
323 
324 *-
325 */
326 
palDcc2s(double v[3],double * a,double * b)327 void palDcc2s ( double v[3], double *a, double *b ) {
328   eraC2s( v, a, b );
329 }
330 
331 /*
332 *+
333 *  Name:
334 *     palDcs2c
335 
336 *  Purpose:
337 *     Spherical coordinates to direction cosines
338 
339 *  Language:
340 *     Starlink ANSI C
341 
342 *  Type of Module:
343 *     Library routine
344 
345 *  Invocation:
346 *     palDcs2c( double a, double b, double v[3] );
347 
348 *  Arguments:
349 *     a = double (Given)
350 *        Spherical coordinate in radians (ra, long etc).
351 *     b = double (Given)
352 *        Spherical coordinate in radians (dec, lat etc).
353 *     v = double [3] (Returned)
354 *        x, y, z vector
355 
356 *  Description:
357 *     The spherical coordinates are longitude (+ve anticlockwise looking
358 *     from the +ve latitude pole) and latitude.  The Cartesian coordinates
359 *     are right handed, with the x axis at zero longitude and latitude, and
360 *     the z axis at the +ve latitude pole.
361 
362 *  Notes:
363 *     - Uses eraS2c(). See SOFA/ERFA documentation for details.
364 
365 *-
366 */
367 
palDcs2c(double a,double b,double v[3])368 void palDcs2c ( double a, double b, double v[3] ) {
369   eraS2c( a, b, v );
370 }
371 
372 /*
373 *+
374 *  Name:
375 *     palDd2tf
376 
377 *  Purpose:
378 *     Convert an interval in days into hours, minutes, seconds
379 
380 *  Language:
381 *     Starlink ANSI C
382 
383 *  Type of Module:
384 *     Library routine
385 
386 *  Invocation:
387 *     palDd2tf( int ndp, double days, char *sign, int ihmsf[4] );
388 
389 *  Arguments:
390 *     ndp = int (Given)
391 *        Number of decimal places of seconds
392 *     days = double (Given)
393 *        Interval in days
394 *     sign = char * (Returned)
395 *        '+' or '-' (single character, not string)
396 *     ihmsf = int [4] (Returned)
397 *        Hours, minutes, seconds, fraction
398 
399 *  Description:
400 *     Convert and interval in days into hours, minutes, seconds.
401 
402 *  Notes:
403 *     - Uses eraD2tf(). See SOFA/ERFA documentation for details.
404 
405 *-
406 */
407 
palDd2tf(int ndp,double days,char * sign,int ihmsf[4])408 void palDd2tf ( int ndp, double days, char *sign, int ihmsf[4] ) {
409   eraD2tf( ndp, days, sign, ihmsf );
410 }
411 
412 /*
413 *+
414 *  Name:
415 *     palDimxv
416 
417 *  Purpose:
418 *     Perform the 3-D backward unitary transformation
419 
420 *  Language:
421 *     Starlink ANSI C
422 
423 *  Type of Module:
424 *     Library routine
425 
426 *  Invocation:
427 *     palDimxv( double dm[3][3], double va[3], double vb[3] );
428 
429 *  Arguments:
430 *     dm = double [3][3] (Given)
431 *        Matrix
432 *     va = double [3] (Given)
433 *        vector
434 *     vb = double [3] (Returned)
435 *        Result vector
436 
437 *  Description:
438 *     Perform the 3-D backward unitary transformation.
439 
440 *  Notes:
441 *     - Uses eraTrxp(). See SOFA/ERFA documentation for details.
442 
443 *-
444 */
445 
palDimxv(double dm[3][3],double va[3],double vb[3])446 void palDimxv ( double dm[3][3], double va[3], double vb[3] ) {
447   eraTrxp( dm, va, vb );
448 }
449 
450 /*
451 *+
452 *  Name:
453 *     palDm2av
454 
455 *  Purpose:
456 *     From a rotation matrix, determine the corresponding axial vector
457 
458 *  Language:
459 *     Starlink ANSI C
460 
461 *  Type of Module:
462 *     Library routine
463 
464 *  Invocation:
465 *     palDm2av( double rmat[3][3], double axvec[3] );
466 
467 *  Arguments:
468 *     rmat = double [3][3] (Given)
469 *        Rotation matrix
470 *     axvec = double [3] (Returned)
471 *        Axial vector (radians)
472 
473 *  Description:
474 *     A rotation matrix describes a rotation about some arbitrary axis,
475 *     called the Euler axis.  The "axial vector" returned by this routine
476 *     has the same direction as the Euler axis, and its magnitude is the
477 *     amount of rotation in radians.  (The magnitude and direction can be
478 *     separated by means of the routine palDvn.)
479 
480 *  Notes:
481 *     - Uses eraRm2v(). See SOFA/ERFA documentation for details.
482 
483 *-
484 */
485 
palDm2av(double rmat[3][3],double axvec[3])486 void palDm2av ( double rmat[3][3], double axvec[3] ) {
487   eraRm2v( rmat, axvec );
488 }
489 
490 /*
491 *+
492 *  Name:
493 *     palDjcl
494 
495 *  Purpose:
496 *     Modified Julian Date to Gregorian year, month, day and fraction of day
497 
498 *  Language:
499 *     Starlink ANSI C
500 
501 *  Type of Module:
502 *     Library routine
503 
504 *  Invocation:
505 *     palDjcl( double djm, int *iy, int *im, int *id, double *fd, int *j );
506 
507 *  Arguments:
508 *     djm = double (Given)
509 *        modified Julian Date (JD-2400000.5)
510 *     iy = int * (Returned)
511 *        year
512 *     im = int * (Returned)
513 *        month
514 *     id = int * (Returned)
515 *        day
516 *     fd = double * (Returned)
517 *        Fraction of day.
518 
519 *  Description:
520 *     Modified Julian Date to Gregorian year, month, day and fraction of day.
521 
522 *  Notes:
523 *     - Uses eraJd2cal(). See SOFA/ERFA documentation for details.
524 
525 *-
526 */
527 
528 /* Requires additional SLA MJD reference date */
palDjcl(double djm,int * iy,int * im,int * id,double * fd,int * j)529 void palDjcl ( double djm, int *iy, int *im, int *id, double *fd, int *j ) {
530   *j = eraJd2cal( PAL__MJD0, djm, iy, im, id, fd );
531 }
532 
533 /*
534 *+
535 *  Name:
536 *     palDmxm
537 
538 *  Purpose:
539 *     Product of two 3x3 matrices
540 
541 *  Language:
542 *     Starlink ANSI C
543 
544 *  Type of Module:
545 *     Library routine
546 
547 *  Invocation:
548 *     palDmxm( double a[3][3], double b[3][3], double c[3][3] );
549 
550 *  Arguments:
551 *     a = double [3][3] (Given)
552 *        Matrix
553 *     b = double [3][3] (Given)
554 *        Matrix
555 *     c = double [3][3] (Returned)
556 *        Matrix result
557 
558 *  Description:
559 *     Product of two 3x3 matrices.
560 
561 *  Notes:
562 *     - Uses eraRxr(). See SOFA/ERFA documentation for details.
563 
564 *-
565 */
566 
palDmxm(double a[3][3],double b[3][3],double c[3][3])567 void palDmxm ( double a[3][3], double b[3][3], double c[3][3] ) {
568   eraRxr( a, b, c );
569 }
570 
571 /*
572 *+
573 *  Name:
574 *     palDmxv
575 
576 *  Purpose:
577 *     Performs the 3-D forward unitary transformation
578 
579 *  Language:
580 *     Starlink ANSI C
581 
582 *  Type of Module:
583 *     Library routine
584 
585 *  Invocation:
586 *     palDmxv( double dm[3][3], double va[3], double vb[3] );
587 
588 *  Arguments:
589 *     dm = double [3][3] (Given)
590 *        matrix
591 *     va = double [3] (Given)
592 *        vector
593 *     dp = double [3] (Returned)
594 *        result vector
595 
596 *  Description:
597 *     Performs the 3-D forward unitary transformation.
598 
599 *  Notes:
600 *     - Uses eraRxp(). See SOFA/ERFA documentation for details.
601 
602 *-
603 */
604 
palDmxv(double dm[3][3],double va[3],double vb[3])605 void palDmxv ( double dm[3][3], double va[3], double vb[3] ) {
606   eraRxp( dm, va, vb );
607 }
608 
609 /*
610 *+
611 *  Name:
612 *     palDpav
613 
614 *  Purpose:
615 *     Position angle of one celestial direction with respect to another
616 
617 *  Language:
618 *     Starlink ANSI C
619 
620 *  Type of Module:
621 *     Library routine
622 
623 *  Invocation:
624 *     pa = palDpav( double v1[3], double v2[3] );
625 
626 *  Arguments:
627 *     v1 = double [3] (Given)
628 *        direction cosines of one point.
629 *     v2 = double [3] (Given)
630 *        direction cosines of the other point.
631 
632 *  Returned Value:
633 *     The result is the bearing (position angle), in radians, of point
634 *     V2 with respect to point V1.  It is in the range +/- pi.  The
635 *     sense is such that if V2 is a small distance east of V1, the
636 *     bearing is about +pi/2.  Zero is returned if the two points
637 *     are coincident.
638 
639 *  Description:
640 *     Position angle of one celestial direction with respect to another.
641 
642 *  Notes:
643 *     - The coordinate frames correspond to RA,Dec, Long,Lat etc.
644 *     - Uses eraPap(). See SOFA/ERFA documentation for details.
645 
646 *-
647 */
648 
palDpav(double v1[3],double v2[3])649 double palDpav ( double v1[3], double v2[3] ) {
650   return eraPap( v1, v2 );
651 }
652 
653 /*
654 *+
655 *  Name:
656 *     palDr2af
657 
658 *  Purpose:
659 *     Convert an angle in radians to degrees, arcminutes, arcseconds
660 
661 *  Language:
662 *     Starlink ANSI C
663 
664 *  Type of Module:
665 *     Library routine
666 
667 *  Invocation:
668 *     palDr2af( int ndp, double angle, char *sign, int idmsf[4] );
669 
670 *  Arguments:
671 *     ndp = int (Given)
672 *        number of decimal places of arcseconds
673 *     angle = double (Given)
674 *        angle in radians
675 *     sign = char * (Returned)
676 *        '+' or '-' (single character)
677 *     idmsf = int [4] (Returned)
678 *        Degrees, arcminutes, arcseconds, fraction
679 
680 *  Description:
681 *     Convert an angle in radians to degrees, arcminutes, arcseconds.
682 
683 *  Notes:
684 *     - Uses eraA2af(). See SOFA/ERFA documentation for details.
685 
686 *-
687 */
688 
palDr2af(int ndp,double angle,char * sign,int idmsf[4])689 void palDr2af ( int ndp, double angle, char *sign, int idmsf[4] ) {
690   eraA2af( ndp, angle, sign, idmsf );
691 }
692 
693 /*
694 *+
695 *  Name:
696 *     palDr2tf
697 
698 *  Purpose:
699 *     Convert an angle in radians to hours, minutes, seconds
700 
701 *  Language:
702 *     Starlink ANSI C
703 
704 *  Type of Module:
705 *     Library routine
706 
707 *  Invocation:
708 *     palDr2tf ( int ndp, double angle, char *sign, int ihmsf[4] );
709 
710 *  Arguments:
711 *     ndp = int (Given)
712 *        number of decimal places of arcseconds
713 *     angle = double (Given)
714 *        angle in radians
715 *     sign = char * (Returned)
716 *        '+' or '-' (single character)
717 *     idmsf = int [4] (Returned)
718 *        Hours, minutes, seconds, fraction
719 
720 *  Description:
721 *     Convert an angle in radians to hours, minutes, seconds.
722 
723 *  Notes:
724 *     - Uses eraA2tf(). See SOFA/ERFA documentation for details.
725 
726 *-
727 */
728 
palDr2tf(int ndp,double angle,char * sign,int ihmsf[4])729 void palDr2tf( int ndp, double angle, char *sign, int ihmsf[4] ) {
730   eraA2tf( ndp, angle, sign, ihmsf );
731 }
732 
733 /*
734 *+
735 *  Name:
736 *     palDranrm
737 
738 *  Purpose:
739 *     Normalize angle into range 0-2 pi
740 
741 *  Language:
742 *     Starlink ANSI C
743 
744 *  Type of Module:
745 *     Library routine
746 
747 *  Invocation:
748 *     norm = palDranrm( double angle );
749 
750 *  Arguments:
751 *     angle = double (Given)
752 *        angle in radians
753 
754 *  Returned Value:
755 *     Angle expressed in the range 0-2 pi
756 
757 *  Description:
758 *     Normalize angle into range 0-2 pi.
759 
760 *  Notes:
761 *     - Uses eraAnp(). See SOFA/ERFA documentation for details.
762 
763 *-
764 */
765 
palDranrm(double angle)766 double palDranrm ( double angle ) {
767   return eraAnp( angle );
768 }
769 
770 /*
771 *+
772 *  Name:
773 *     palDsep
774 
775 *  Purpose:
776 *     Angle between two points on a sphere
777 
778 *  Language:
779 *     Starlink ANSI C
780 
781 *  Type of Module:
782 *     Library routine
783 
784 *  Invocation:
785 *     ang = palDsep( double a1, double b1, double a2, double b2 );
786 
787 *  Arguments:
788 *     a1 = double (Given)
789 *        Spherical coordinate of one point (radians)
790 *     b1 = double (Given)
791 *        Spherical coordinate of one point (radians)
792 *     a2 = double (Given)
793 *        Spherical coordinate of other point (radians)
794 *     b2 = double (Given)
795 *        Spherical coordinate of other point (radians)
796 
797 *  Returned Value:
798 *     Angle, in radians, between the two points. Always positive.
799 
800 *  Description:
801 *     Angle between two points on a sphere.
802 
803 *  Notes:
804 *     - The spherical coordinates are [RA,Dec], [Long,Lat] etc, in radians.
805 *     - Uses eraSeps(). See SOFA/ERFA documentation for details.
806 
807 *-
808 */
809 
palDsep(double a1,double b1,double a2,double b2)810 double palDsep ( double a1, double b1, double a2, double b2 ) {
811   return eraSeps( a1, b1, a2, b2 );
812 }
813 
814 /*
815 *+
816 *  Name:
817 *     palDsepv
818 
819 *  Purpose:
820 *     Angle between two vectors
821 
822 *  Language:
823 *     Starlink ANSI C
824 
825 *  Type of Module:
826 *     Library routine
827 
828 *  Invocation:
829 *     ang = palDsepv( double v1[3], double v2[3] );
830 
831 *  Arguments:
832 *     v1 = double [3] (Given)
833 *        First vector
834 *     v2 = double [3] (Given)
835 *        Second vector
836 
837 *  Returned Value:
838 *     Angle, in radians, between the two points. Always positive.
839 
840 *  Description:
841 *     Angle between two vectors.
842 
843 *  Notes:
844 *     - Uses eraSepp(). See SOFA/ERFA documentation for details.
845 
846 *-
847 */
848 
palDsepv(double v1[3],double v2[3])849 double palDsepv ( double v1[3], double v2[3] ) {
850   return eraSepp( v1, v2 );
851 }
852 
853 /*
854 *+
855 *  Name:
856 *     palDtf2d
857 
858 *  Purpose:
859 *     Convert hours, minutes, seconds to days
860 
861 *  Language:
862 *     Starlink ANSI C
863 
864 *  Type of Module:
865 *     Library routine
866 
867 *  Invocation:
868 *     palDtf2d( int ihour, int imin, double sec, double *days, int *j );
869 
870 *  Arguments:
871 *     ihour = int (Given)
872 *        Hours
873 *     imin = int (Given)
874 *        Minutes
875 *     sec = double (Given)
876 *        Seconds
877 *     days = double * (Returned)
878 *        Interval in days
879 *     j = int * (Returned)
880 *        status: 0 = ok, 1 = ihour outside range 0-23,
881 *        2 = imin outside range 0-59, 3 = sec outside range 0-59.999...
882 
883 *  Description:
884 *     Convert hours, minutes, seconds to days.
885 
886 *  Notes:
887 *     - Uses eraTf2d(). See SOFA/ERFA documentation for details.
888 
889 *-
890 */
891 
892 /* Assumes that the sign is always positive and is dealt with externally */
palDtf2d(int ihour,int imin,double sec,double * days,int * j)893 void palDtf2d ( int ihour, int imin, double sec, double *days, int *j ) {
894   *j = eraTf2d( ' ', ihour, imin, sec, days );
895 }
896 
897 /*
898 *+
899 *  Name:
900 *     palDtf2r
901 
902 *  Purpose:
903 *     Convert hours, minutes, seconds to radians
904 
905 *  Language:
906 *     Starlink ANSI C
907 
908 *  Type of Module:
909 *     Library routine
910 
911 *  Invocation:
912 *     palDtf2r( int ihour, int imin, double sec, double *rad, int *j );
913 
914 *  Arguments:
915 *     ihour = int (Given)
916 *        Hours
917 *     imin = int (Given)
918 *        Minutes
919 *     sec = double (Given)
920 *        Seconds
921 *     days = double * (Returned)
922 *        Angle in radians
923 *     j = int * (Returned)
924 *        status: 0 = ok, 1 = ihour outside range 0-23,
925 *        2 = imin outside range 0-59, 3 = sec outside range 0-59.999...
926 
927 *  Description:
928 *     Convert hours, minutes, seconds to radians.
929 
930 *  Notes:
931 *     - Uses eraTf2a(). See SOFA/ERFA documentation for details.
932 
933 *-
934 */
935 
936 /* Assumes that the sign is dealt with outside this routine */
palDtf2r(int ihour,int imin,double sec,double * rad,int * j)937 void palDtf2r ( int ihour, int imin, double sec, double *rad, int *j ) {
938   *j = eraTf2a( ' ', ihour, imin, sec, rad );
939 }
940 
941 /*
942 *+
943 *  Name:
944 *     palDvdv
945 
946 *  Purpose:
947 *     Scalar product of two 3-vectors
948 
949 *  Language:
950 *     Starlink ANSI C
951 
952 *  Type of Module:
953 *     Library routine
954 
955 *  Invocation:
956 *     prod = palDvdv ( double va[3], double vb[3] );
957 
958 *  Arguments:
959 *     va = double [3] (Given)
960 *        First vector
961 *     vb = double [3] (Given)
962 *        Second vector
963 
964 *  Returned Value:
965 *     Scalar product va.vb
966 
967 *  Description:
968 *     Scalar product of two 3-vectors.
969 
970 *  Notes:
971 *     - Uses eraPdp(). See SOFA/ERFA documentation for details.
972 
973 *-
974 */
975 
palDvdv(double va[3],double vb[3])976 double palDvdv ( double va[3], double vb[3] ) {
977   return eraPdp( va, vb );
978 }
979 
980 /*
981 *+
982 *  Name:
983 *     palDvn
984 
985 *  Purpose:
986 *     Normalizes a 3-vector also giving the modulus
987 
988 *  Language:
989 *     Starlink ANSI C
990 
991 *  Type of Module:
992 *     Library routine
993 
994 *  Invocation:
995 *     palDvn( double v[3], double uv[3], double *vm );
996 
997 *  Arguments:
998 *     v = double [3] (Given)
999 *        vector
1000 *     uv = double [3] (Returned)
1001 *        unit vector in direction of "v"
1002 *     vm = double * (Returned)
1003 *        modulus of "v"
1004 
1005 *  Description:
1006 *     Normalizes a 3-vector also giving the modulus.
1007 
1008 *  Notes:
1009 *     - Uses eraPn(). See SOFA/ERFA documentation for details.
1010 
1011 *-
1012 */
1013 
1014 /* Note that the arguments are flipped */
palDvn(double v[3],double uv[3],double * vm)1015 void palDvn ( double v[3], double uv[3], double *vm ) {
1016   eraPn( v, vm, uv );
1017 }
1018 
1019 /*
1020 *+
1021 *  Name:
1022 *     palDvxv
1023 
1024 *  Purpose:
1025 *     Vector product of two 3-vectors
1026 
1027 *  Language:
1028 *     Starlink ANSI C
1029 
1030 *  Type of Module:
1031 *     Library routine
1032 
1033 *  Invocation:
1034 *     palDvxv( double va[3], double vb[3], double vc[3] );
1035 
1036 *  Arguments:
1037 *     va = double [3] (Given)
1038 *        First vector
1039 *     vb = double [3] (Given)
1040 *        Second vector
1041 *     vc = double [3] (Returned)
1042 *        Result vector
1043 
1044 *  Description:
1045 *     Vector product of two 3-vectors.
1046 
1047 *  Notes:
1048 *     - Uses eraPxp(). See SOFA/ERFA documentation for details.
1049 
1050 *-
1051 */
1052 
palDvxv(double va[3],double vb[3],double vc[3])1053 void palDvxv ( double va[3], double vb[3], double vc[3] ) {
1054   eraPxp( va, vb, vc );
1055 }
1056 
1057 /*
1058 *+
1059 *  Name:
1060 *     palEpb
1061 
1062 *  Purpose:
1063 *     Conversion of modified Julian Data to Besselian Epoch
1064 
1065 *  Language:
1066 *     Starlink ANSI C
1067 
1068 *  Type of Module:
1069 *     Library routine
1070 
1071 *  Invocation:
1072 *     epb = palEpb ( double date );
1073 
1074 *  Arguments:
1075 *     date = double (Given)
1076 *        Modified Julian Date (JD - 2400000.5)
1077 
1078 *  Returned Value:
1079 *      Besselian epoch.
1080 
1081 *  Description:
1082 *     Conversion of modified Julian Data to Besselian Epoch.
1083 
1084 *  Notes:
1085 *     - Uses eraEpb(). See SOFA/ERFA documentation for details.
1086 
1087 *-
1088 */
1089 
1090 /* Requires additional SLA MJD reference date */
palEpb(double date)1091 double palEpb ( double date ) {
1092   return eraEpb( PAL__MJD0, date );
1093 }
1094 
1095 /*
1096 *+
1097 *  Name:
1098 *     palEpb2d
1099 
1100 *  Purpose:
1101 *     Conversion of Besselian Epoch to Modified Julian Date
1102 
1103 *  Language:
1104 *     Starlink ANSI C
1105 
1106 *  Type of Module:
1107 *     Library routine
1108 
1109 *  Invocation:
1110 *     mjd = palEpb2d ( double epb );
1111 
1112 *  Arguments:
1113 *     epb = double (Given)
1114 *        Besselian Epoch
1115 
1116 *  Returned Value:
1117 *     Modified Julian Date (JD - 2400000.5)
1118 
1119 *  Description:
1120 *     Conversion of Besselian Epoch to Modified Julian Date.
1121 
1122 *  Notes:
1123 *     - Uses eraEpb2jd(). See SOFA/ERFA documentation for details.
1124 
1125 *-
1126 */
1127 
1128 
palEpb2d(double epb)1129 double palEpb2d ( double epb ) {
1130   double djm0, djm;
1131   eraEpb2jd( epb, &djm0, &djm );
1132   return djm;
1133 }
1134 
1135 /*
1136 *+
1137 *  Name:
1138 *     palEpj
1139 
1140 *  Purpose:
1141 *     Conversion of Modified Julian Date to Julian Epoch
1142 
1143 *  Language:
1144 *     Starlink ANSI C
1145 
1146 *  Type of Module:
1147 *     Library routine
1148 
1149 *  Invocation:
1150 *     epj = palEpj ( double date );
1151 
1152 *  Arguments:
1153 *     date = double (Given)
1154 *        Modified Julian Date (JD - 2400000.5)
1155 
1156 *  Returned Value:
1157 *     The Julian Epoch.
1158 
1159 *  Description:
1160 *     Conversion of Modified Julian Date to Julian Epoch.
1161 
1162 *  Notes:
1163 *     - Uses eraEpj(). See SOFA/ERFA documentation for details.
1164 
1165 *-
1166 */
1167 
1168 /* Requires additional SLA MJD reference date */
palEpj(double date)1169 double palEpj ( double date ) {
1170   return eraEpj( PAL__MJD0, date );
1171 }
1172 
1173 /*
1174 *+
1175 *  Name:
1176 *     palEpj2d
1177 
1178 *  Purpose:
1179 *     Conversion of Julian Epoch to Modified Julian Date
1180 
1181 *  Language:
1182 *     Starlink ANSI C
1183 
1184 *  Type of Module:
1185 *     Library routine
1186 
1187 *  Invocation:
1188 *     mjd = palEpj2d ( double epj );
1189 
1190 *  Arguments:
1191 *     epj = double (Given)
1192 *        Julian Epoch.
1193 
1194 *  Returned Value:
1195 *     Modified Julian Date (JD - 2400000.5)
1196 
1197 *  Description:
1198 *     Conversion of Julian Epoch to Modified Julian Date.
1199 
1200 *  Notes:
1201 *     - Uses eraEpj2d(). See SOFA/ERFA documentation for details.
1202 
1203 *-
1204 */
palEpj2d(double epj)1205 double palEpj2d ( double epj ) {
1206   double djm0, djm;
1207   eraEpj2jd( epj, &djm0, &djm );
1208   return djm;
1209 }
1210 
1211 /*
1212 *+
1213 *  Name:
1214 *     palEqeqx
1215 
1216 *  Purpose:
1217 *     Equation of the equinoxes (IAU 2000/2006)
1218 
1219 *  Language:
1220 *     Starlink ANSI C
1221 
1222 *  Type of Module:
1223 *     Library routine
1224 
1225 *  Invocation:
1226 *     palEqeqx( double date );
1227 
1228 *  Arguments:
1229 *     date = double (Given)
1230 *        TT as Modified Julian Date (JD-400000.5)
1231 
1232 *  Description:
1233 *     Equation of the equinoxes (IAU 2000/2006).
1234 
1235 *  Notes:
1236 *     - Uses eraEe06a(). See SOFA/ERFA documentation for details.
1237 
1238 *-
1239 */
1240 
1241 /* Requires additional SLA MJD reference date */
palEqeqx(double date)1242 double palEqeqx ( double date ) {
1243   return eraEe06a( PAL__MJD0, date );
1244 }
1245 
1246 /*
1247 *+
1248 *  Name:
1249 *     palFk5hz
1250 
1251 *  Purpose:
1252 *     Transform an FK5 (J2000) star position into the frame of the
1253 *     Hipparcos catalogue.
1254 
1255 *  Language:
1256 *     Starlink ANSI C
1257 
1258 *  Type of Module:
1259 *     Library routine
1260 
1261 *  Invocation:
1262 *     palFk5hz ( double r5, double d5, double epoch,
1263 *                double *rh, double *dh );
1264 
1265 *  Arguments:
1266 *     r5 = double (Given)
1267 *        FK5 RA (radians), equinox J2000, epoch "epoch"
1268 *     d5 = double (Given)
1269 *        FK5 dec (radians), equinox J2000, epoch "epoch"
1270 *     epoch = double (Given)
1271 *        Julian epoch
1272 *     rh = double * (Returned)
1273 *        RA (radians)
1274 *     dh = double * (Returned)
1275 *        Dec (radians)
1276 
1277 *  Description:
1278 *     Transform an FK5 (J2000) star position into the frame of the
1279 *     Hipparcos catalogue.
1280 
1281 *  Notes:
1282 *     - Assumes zero Hipparcos proper motion.
1283 *     - Uses eraEpj2jd() and eraFk5hz.
1284 *       See SOFA/ERFA documentation for details.
1285 
1286 *-
1287 */
1288 
palFk5hz(double r5,double d5,double epoch,double * rh,double * dh)1289 void palFk5hz ( double r5, double d5, double epoch,
1290                 double *rh, double *dh ) {
1291   /* Need to convert epoch to Julian date first */
1292   double date1, date2;
1293   eraEpj2jd( epoch, &date1, &date2 );
1294   eraFk5hz( r5, d5, date1, date2, rh, dh );
1295 }
1296 
1297 /*
1298 *+
1299 *  Name:
1300 *     palGmst
1301 
1302 *  Purpose:
1303 *     Greenwich mean sidereal time (consistent with IAU 2006 precession).
1304 
1305 *  Language:
1306 *     Starlink ANSI C
1307 
1308 *  Type of Module:
1309 *     Library routine
1310 
1311 *  Invocation:
1312 *     mst = palGmst ( double ut1 );
1313 
1314 *  Arguments:
1315 *     ut1 = double (Given)
1316 *        Universal time (UT1) expressed as modified Julian Date (JD-2400000.5)
1317 
1318 *  Returned Value:
1319 *     Greenwich mean sidereal time
1320 
1321 *  Description:
1322 *     Greenwich mean sidereal time (consistent with IAU 2006 precession).
1323 
1324 *  Notes:
1325 *     - Uses eraGmst06(). See SOFA/ERFA documentation for details.
1326 
1327 *-
1328 */
1329 
1330 /* Note that SOFA/ERFA has more accurate time arguments
1331    and we use the 2006 precession model */
palGmst(double ut1)1332 double palGmst ( double ut1 ) {
1333   return eraGmst06( PAL__MJD0, ut1, PAL__MJD0, ut1 );
1334 }
1335 
1336 /*
1337 *+
1338 *  Name:
1339 *     palGmsta
1340 
1341 *  Purpose:
1342 *     Greenwich mean sidereal time (consistent with IAU 2006 precession).
1343 
1344 *  Language:
1345 *     Starlink ANSI C
1346 
1347 *  Type of Module:
1348 *     Library routine
1349 
1350 *  Invocation:
1351 *     mst = palGmsta ( double date, double ut1 );
1352 
1353 *  Arguments:
1354 *     date = double (Given)
1355 *        UT1 date (MJD: integer part of JD-2400000.5)
1356 *     ut1 = double (Given)
1357 *        UT1 time (fraction of a day)
1358 
1359 *  Returned Value:
1360 *     Greenwich mean sidereal time (in range 0 to 2 pi)
1361 
1362 *  Description:
1363 *     Greenwich mean sidereal time (consistent with IAU 2006 precession).
1364 
1365 *  Notes:
1366 *     - For best accuracy use eraGmst06() directly.
1367 *     - Uses eraGmst06(). See SOFA/ERFA documentation for details.
1368 
1369 *-
1370 */
1371 
1372 /* Slightly better but still not as accurate as SOFA/ERFA */
1373 
palGmsta(double date,double ut)1374 double palGmsta( double date, double ut ) {
1375   date += PAL__MJD0;
1376   return eraGmst06( date, ut, date, ut );
1377 }
1378 
1379 /*
1380 *+
1381 *  Name:
1382 *     palHfk5z
1383 
1384 *  Purpose:
1385 *     Hipparcos star position to FK5 J2000
1386 
1387 *  Language:
1388 *     Starlink ANSI C
1389 
1390 *  Type of Module:
1391 *     Library routine
1392 
1393 *  Invocation:
1394 *     palHfk5z( double rh, double dh, double epoch,
1395 *               double *r5, double *d5, double *dr5, double *dd5 );
1396 
1397 *  Arguments:
1398 *     rh = double (Given)
1399 *        Hipparcos RA (radians)
1400 *     dh = double (Given)
1401 *        Hipparcos Dec (radians)
1402 *     epoch = double (Given)
1403 *        Julian epoch (TDB)
1404 *     r5 = double * (Returned)
1405 *        RA (radians, FK5, equinox J2000, epoch "epoch")
1406 *     d5 = double * (Returned)
1407 *        Dec (radians, FK5, equinox J2000, epoch "epoch")
1408 
1409 *  Description:
1410 *     Transform a Hipparcos star position into FK5 J2000, assuming
1411 *     zero Hipparcos proper motion.
1412 
1413 *  Notes:
1414 *     - Uses eraEpj2jd and eraHfk5z(). See SOFA/ERFA documentation for details.
1415 
1416 *-
1417 */
1418 
palHfk5z(double rh,double dh,double epoch,double * r5,double * d5,double * dr5,double * dd5)1419 void palHfk5z ( double rh, double dh, double epoch,
1420                 double *r5, double *d5, double *dr5, double *dd5 ) {
1421   /* Need to convert epoch to Julian date first */
1422   double date1, date2;
1423   eraEpj2jd( epoch, &date1, &date2 );
1424   eraHfk5z( rh, dh, date1, date2, r5, d5, dr5, dd5 );
1425 }
1426 
1427 /*
1428 *+
1429 *  Name:
1430 *     palRefcoq
1431 
1432 *  Purpose:
1433 *     Determine the constants A and B in the atmospheric refraction model
1434 
1435 *  Language:
1436 *     Starlink ANSI C
1437 
1438 *  Type of Module:
1439 *     Library routine
1440 
1441 *  Invocation:
1442 *     palRefcoq( double tdk, double pmb, double rh, double wl,
1443 *                double *refa, double *refb );
1444 
1445 *  Arguments:
1446 *     tdk = double (Given)
1447 *        Ambient temperature at the observer (K)
1448 *     pmb = double (Given)
1449 *        Pressure at the observer (millibar)
1450 *     rh =  double (Given)
1451 *        Relative humidity at the observer (range 0-1)
1452 *     wl =  double (Given)
1453 *        Effective wavelength of the source (micrometre).
1454 *        Radio refraction is chosen by specifying wl > 100 micrometres.
1455 *     refa = double * (Returned)
1456 *        tan Z coefficient (radian)
1457 *     refb = double * (Returned)
1458 *        tan**3 Z coefficient (radian)
1459 
1460 *  Description:
1461 *     Determine the constants A and B in the atmospheric refraction
1462 *     model dZ = A tan Z + B tan**3 Z.  This is a fast alternative
1463 *     to the palRefco routine.
1464 *
1465 *     Z is the "observed" zenith distance (i.e. affected by refraction)
1466 *     and dZ is what to add to Z to give the "topocentric" (i.e. in vacuo)
1467 *     zenith distance.
1468 
1469 *  Notes:
1470 *     - Uses eraRefco(). See SOFA/ERFA documentation for details.
1471 *     - Note that the SOFA/ERFA routine uses different order of
1472 *       of arguments and uses deg C rather than K.
1473 
1474 *-
1475 */
1476 
palRefcoq(double tdk,double pmb,double rh,double wl,double * refa,double * refb)1477 void palRefcoq ( double tdk, double pmb, double rh, double wl,
1478                  double *refa, double *refb ) {
1479   /* Note that SLA (and therefore PAL) uses units of kelvin
1480      but SOFA/ERFA uses deg C */
1481   eraRefco( pmb, tdk - 273.15, rh, wl, refa, refb );
1482 }
1483