update ast 8.6.2

40 files changed, 6669 insertions, 0 deletions

diff --git a/ast/pal/pal.h b/ast/pal/pal.h
new file mode 100644
index 0000000..a588168
--- /dev/null
+++ b/ast/pal/pal.h
@@ -0,0 +1,551 @@
+#ifndef PALHDEF
+#define PALHDEF
+
+/*
+*+
+*  Name:
+*     pal.h
+
+*  Purpose:
+*     Function prototypes for PAL routines.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Include file
+
+*  Description:
+*     Function prototypes for PAL routines.
+
+*  Authors:
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version. Define all SLA prototypes in PAL form even
+*        though none are implemented.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+
+void palAddet ( double rm, double dm, double eq, double *rc, double *dc );
+
+void palAfin ( const char *string, int *iptr, float *a, int *j );
+
+double palAirmas ( double zd );
+
+void palAltaz ( double ha, double dec, double phi,
+                double *az, double *azd, double *azdd,
+                double *el, double *eld, double *eldd,
+                double *pa, double *pad, double *padd );
+
+void palAmp ( double ra, double da, double date, double eq,
+              double *rm, double *dm );
+
+void palAmpqk ( double ra, double da, double amprms[21],
+                double *rm, double *dm );
+
+void palAop ( double rap, double dap, double date, double dut,
+              double elongm, double phim, double hm, double xp,
+              double yp, double tdk, double pmb, double rh,
+              double wl, double tlr,
+              double *aob, double *zob, double *hob,
+              double *dob, double *rob );
+
+void palAoppa ( double date, double dut, double elongm, double phim,
+                double hm, double xp, double yp, double tdk, double pmb,
+                double rh, double wl, double tlr, double aoprms[14] );
+
+void palAoppat ( double date, double aoprms[14] );
+
+void palAopqk ( double rap, double dap, const double aoprms[14],
+                double *aob, double *zob, double *hob,
+                double *dob, double *rob );
+
+void palAtmdsp ( double tdk, double pmb, double rh, double wl1,
+                 double a1, double b1, double wl2, double *a2, double *b2 );
+
+void palAv2m ( float axvec[3], float rmat[3][3] );
+
+float palBear ( float a1, float b1, float a2, float b2 );
+
+void palCaf2r ( int ideg, int iamin, float asec, float *rad, int *j );
+
+void palCaldj ( int iy, int im, int id, double *djm, int *j );
+
+void palCalyd ( int iy, int im, int id, int *ny, int *nd, int *j );
+
+void palCc2s ( float v[3], float *a, float *b );
+
+void palCc62s ( float v[6], float *a, float *b, float *r,
+                float *ad, float *bd, float *rd );
+
+void palCd2tf ( int ndp, float days, char *sign, int ihmsf[4] );
+
+void palCldj ( int iy, int im, int id, double *djm, int *j );
+
+void palClyd ( int iy, int im, int id, int *ny, int *nd, int *jstat );
+
+void palCombn ( int nsel, int ncand, int list[], int *j );
+
+void palCr2af ( int ndp, float angle, char *sign, int idmsf[4] );
+
+void palCr2tf ( int ndp, float angle, char *sign, int ihmsf[4] );
+
+void palCs2c ( float a, float b, float v[3] );
+
+void palCs2c6 ( float a, float b, float r, float ad,
+                float bd, float rd, float v[6] );
+
+void palCtf2d ( int ihour, int imin, float sec, float *days, int *j );
+
+void palCtf2r ( int ihour, int imin, float sec, float *rad, int *j );
+
+void palDaf2r ( int ideg, int iamin, double asec, double *rad, int *j );
+
+void palDafin ( const char *string, int *iptr, double *a, int *j );
+
+double palDat ( double dju );
+
+void palDav2m ( double axvec[3], double rmat[3][3] );
+
+double palDbear ( double a1, double b1, double a2, double b2 );
+
+void palDbjin ( const char *string, int *nstrt,
+                double *dreslt, int *jf1, int *jf2 );
+
+void palDc62s ( double v[6], double *a, double *b, double *r,
+                double *ad, double *bd, double *rd );
+
+void palDcc2s ( double v[3], double *a, double *b );
+
+void palDcmpf ( double coeffs[6], double *xz, double *yz, double *xs,
+                double *ys, double *perp, double *orient );
+
+void palDcs2c ( double a, double b, double v[3] );
+
+void palDd2tf ( int ndp, double days, char *sign, int ihmsf[4] );
+
+void palDe2h ( double ha, double dec, double phi,
+               double *az, double *el );
+
+void palDeuler ( const char *order, double phi, double theta, double psi,
+                 double rmat[3][3] );
+
+void palDfltin ( const char *string, int *nstrt, double *dreslt, int *jflag );
+
+void palDh2e ( double az, double el, double phi, double *ha, double *dec);
+
+void palDimxv ( double dm[3][3], double va[3], double vb[3] );
+
+void palDjcal ( int ndp, double djm, int iymdf[4], int *j );
+
+void palDjcl ( double djm, int *iy, int *im, int *id, double *fd, int *j );
+
+void palDm2av ( double rmat[3][3], double axvec[3] );
+
+void palDmat ( int n, double *a, double *y, double *d, int *jf, int *iw );
+
+void palDmoon ( double date, double pv[6] );
+
+void palDmxm ( double a[3][3], double b[3][3], double c[3][3] );
+
+void palDmxv ( double dm[3][3], double va[3], double vb[3] );
+
+double palDpav ( double v1[3], double v2[3] );
+
+void palDr2af ( int ndp, double angle, char *sign, int idmsf[4] );
+
+void palDr2tf ( int ndp, double angle, char *sign, int ihmsf[4] );
+
+double palDrange ( double angle );
+
+double palDranrm ( double angle );
+
+void palDs2c6 ( double a, double b, double r, double ad, double bd,
+                double rd, double v[6] );
+
+void palDs2tp ( double ra, double dec, double raz, double decz,
+                double *xi, double *eta, int *j );
+
+double palDsep ( double a1, double b1, double a2, double b2 );
+
+double palDsepv ( double v1[3], double v2[3] );
+
+double palDt ( double epoch );
+
+void palDtf2d ( int ihour, int imin, double sec, double *days, int *j );
+
+void palDtf2r ( int ihour, int imin, double sec, double *rad, int *j );
+
+void palDtp2s ( double xi, double eta, double raz, double decz,
+                double *ra, double *dec );
+
+void palDtp2v ( double xi, double eta, double v0[3], double v[3] );
+
+void palDtps2c ( double xi, double eta, double ra, double dec,
+                 double *raz1, double *decz1,
+                 double *raz2, double *decz2, int *n );
+
+void palDtpv2c ( double xi, double eta, double v[3],
+                 double v01[3], double v02[3], int *n );
+
+double palDtt ( double dju );
+
+void palDv2tp ( double v[3], double v0[3], double *xi, double *eta, int *j );
+
+double palDvdv ( double va[3], double vb[3] );
+
+void palDvn ( double v[3], double uv[3], double *vm );
+
+void palDvxv ( double va[3], double vb[3], double vc[3] );
+
+void palE2h ( float ha, float dec, float phi, float *az, float *el );
+
+void palEarth ( int iy, int id, float fd, float posvel[6] );
+
+void palEcleq ( double dl, double db, double date, double *dr, double *dd );
+
+void palEcmat ( double date, double rmat[3][3] );
+
+void palEcor ( float rm, float dm, int iy, int id, float fd,
+               float *rv, float *tl );
+
+void palEg50 ( double dr, double dd, double *dl, double *db );
+
+void palEl2ue ( double date, int jform, double epoch, double orbinc,
+                double anode, double perih, double aorq, double e,
+                double aorl, double dm, double u[13], int *jstat );
+
+double palEpb ( double date );
+
+double palEpb2d ( double epb );
+
+double palEpco ( char k0, char k, double e );
+
+double palEpj ( double date );
+
+double palEpj2d ( double epj );
+
+void palEpv( double date, double ph[3], double vh[3],
+             double pb[3], double vb[3] );
+
+void palEqecl ( double dr, double dd, double date, double *dl, double *db );
+
+double palEqeqx ( double date );
+
+void palEqgal ( double dr, double dd, double *dl, double *db );
+
+void palEtrms ( double ep, double ev[3] );
+
+void palEuler ( const char *order, float phi, float theta, float psi,
+                float rmat[3][3] );
+
+void palEvp ( double date, double deqx,
+              double dvb[3], double dpb[3],
+              double dvh[3], double dph[3] );
+
+void palFitxy ( int itype, int np, double xye[][2], double xym[][2],
+                double coeffs[6], int *j );
+
+void palFk425 ( double r1950, double d1950, double dr1950,
+                double dd1950, double p1950, double v1950,
+                double *r2000, double *d2000, double *dr2000,
+                double *dd2000, double *p2000, double *v2000 );
+
+void palFk45z ( double r1950, double d1950, double bepoch,
+                double *r2000, double *d2000 );
+
+void palFk524 ( double r2000, double d2000, double dr2000,
+                double dd2000, double p2000, double v2000,
+                double *r1950, double *d1950, double *dr1950,
+                double *dd1950, double *p1950, double *v1950 );
+
+void palFk52h ( double r5, double d5, double dr5, double dd5,
+                double *dr, double *dh, double *drh, double *ddh );
+
+void palFk54z ( double r2000, double d2000, double bepoch,
+                double *r1950, double *d1950,
+                double *dr1950, double *dd1950 );
+
+void palFk5hz ( double r5, double d5, double epoch,
+                double *rh, double *dh );
+
+void palFlotin ( const char *string, int *nstrt, float *reslt, int *jflag );
+
+void palGaleq ( double dl, double db, double *dr, double *dd );
+
+void palGalsup ( double dl, double db, double *dsl, double *dsb );
+
+void palGe50 ( double dl, double db, double *dr, double *dd );
+
+void palGeoc ( double p, double h, double *r, double *z );
+
+double palGmst ( double ut1 );
+
+double palGmsta ( double date, double ut1 );
+
+void palH2e ( float az, float el, float phi, float *ha, float *dec );
+
+void palH2fk5 ( double dr, double dh, double drh, double ddh,
+                double *r5, double *d5, double *dr5, double *dd5 );
+
+void palHfk5z ( double rh, double dh, double epoch,
+                double *r5, double *d5, double *dr5, double *dd5 );
+
+void palImxv ( float rm[3][3], float va[3], float vb[3] );
+
+void palInt2in ( const char *string, int *nstrt, int *ireslt, int *jflag );
+
+void palIntin ( const char *string, int *nstrt, long *ireslt, int *jflag );
+
+void palInvf ( double fwds[6], double bkwds[6], int *j );
+
+void palKbj ( int jb, double e, char *k, int *j );
+
+void palM2av ( float rmat[3][3], float axvec[3] );
+
+void palMap ( double rm, double dm, double pr, double pd,
+              double px, double rv, double eq, double date,
+              double *ra, double *da );
+
+void palMappa ( double eq, double date, double amprms[21] );
+
+void palMapqk ( double rm, double dm, double pr, double pd,
+                double px, double rv, double amprms[21],
+                double *ra, double *da );
+
+void palMapqkz ( double rm, double dm, double amprms[21],
+                 double *ra, double *da );
+
+void palMoon ( int iy, int id, float fd, float posvel[6] );
+
+void palMxm ( float a[3][3], float b[3][3], float c[3][3] );
+
+void palMxv ( float rm[3][3], float va[3], float vb[3] );
+
+void palNut ( double date, double rmatn[3][3] );
+
+void palNutc ( double date, double *dpsi, double *deps, double *eps0 );
+
+void palNutc80 ( double date, double *dpsi, double *deps, double *eps0 );
+
+void palOap ( const char *type, double ob1, double ob2, double date,
+              double dut, double elongm, double phim, double hm,
+              double xp, double yp, double tdk, double pmb,
+              double rh, double wl, double tlr,
+              double *rap, double *dap );
+
+void palOapqk ( const char *type, double ob1, double ob2, const double aoprms[14],
+                double *rap, double *dap );
+
+int palObs( size_t n, const char * c,
+            char * ident, size_t identlen,
+            char * name, size_t namelen,
+            double * w, double * p, double * h );
+
+double palPa ( double ha, double dec, double phi );
+
+double palPav ( float v1[3], float v2[3] );
+
+void palPcd ( double disco, double *x, double *y );
+
+void palPda2h ( double p, double d, double a,
+                double *h1, int *j1, double *h2, int *j2 );
+
+void palPdq2h ( double p, double d, double q,
+                double *h1, int *j1, double *h2, int *j2 );
+
+void palPermut ( int n, int istate[], int iorder[], int *j );
+
+void palPertel (int jform, double date0, double date1,
+                double epoch0, double orbi0, double anode0,
+                double perih0, double aorq0, double e0, double am0,
+                double *epoch1, double *orbi1, double *anode1,
+                double *perih1, double *aorq1, double *e1, double *am1,
+                int *jstat );
+
+void palPertue ( double date, double u[13], int *jstat );
+
+void palPlanel ( double date, int jform, double epoch, double orbinc,
+                 double anode, double perih, double aorq,  double e,
+                 double aorl, double dm, double pv[6], int *jstat );
+
+void palPlanet ( double date, int np, double pv[6], int *j );
+
+void palPlante ( double date, double elong, double phi, int jform,
+                 double epoch, double orbinc, double anode, double perih,
+                 double aorq, double e, double aorl, double dm,
+                 double *ra, double *dec, double *r, int *jstat );
+
+void palPlantu ( double date, double elong, double phi, const double u[13],
+                 double *ra, double *dec, double *r, int *jstat );
+
+void palPm ( double r0, double d0, double pr, double pd,
+             double px, double rv, double ep0, double ep1,
+             double *r1, double *d1 );
+
+void palPolmo ( double elongm, double phim, double xp, double yp,
+                double *elong, double *phi, double *daz );
+
+void palPrebn ( double bep0, double bep1, double rmatp[3][3] );
+
+void palPrec ( double ep0, double ep1, double rmatp[3][3] );
+
+void palPrecl ( double ep0, double ep1, double rmatp[3][3] );
+
+void palPreces ( const char sys[3], double ep0, double ep1,
+                 double *ra, double *dc );
+
+void palPrenut ( double epoch, double date, double rmatpn[3][3] );
+
+void palPv2el ( const double pv[6], double date, double pmass, int jformr,
+                int *jform, double *epoch, double *orbinc,
+                double *anode, double *perih, double *aorq, double *e,
+                double *aorl, double *dm, int *jstat );
+
+void palPv2ue ( const double pv[6], double date, double pmass,
+                double u[13], int *jstat );
+
+void palPvobs ( double p, double h, double stl, double pv[6] );
+
+void palPxy ( int np, double xye[][2], double xym[][2],
+              double coeffs[6],
+              double xyp[][2], double *xrms, double *yrms, double *rrms );
+
+float palRange ( float angle );
+
+float palRanorm ( float angle );
+
+double palRcc ( double tdb, double ut1, double wl, double u, double v );
+
+void palRdplan ( double date, int np, double elong, double phi,
+                 double *ra, double *dec, double *diam );
+
+void palRefco ( double hm, double tdk, double pmb, double rh,
+                double wl, double phi, double tlr, double eps,
+                double *refa, double *refb );
+
+void palRefcoq ( double tdk, double pmb, double rh, double wl,
+                double *refa, double *refb );
+
+void palRefro ( double zobs, double hm, double tdk, double pmb,
+                double rh, double wl, double phi, double tlr, double eps,
+                double *ref );
+
+void palRefv ( double vu[3], double refa, double refb, double vr[3] );
+
+void palRefz ( double zu, double refa, double refb, double *zr );
+
+double palRverot ( double phi, double ra, double da, double st );
+
+double palRvgalc ( double r2000, double d2000 );
+
+double palRvlg ( double r2000, double d2000 );
+
+double palRvlsrd ( double r2000, double d2000 );
+
+double palRvlsrk ( double r2000, double d2000 );
+
+void palS2tp ( float ra, float dec, float raz, float decz,
+               float *xi, float *eta, int *j );
+
+float palSep ( float a1, float b1, float a2, float b2 );
+
+float palSepv ( float v1[3], float v2[3] );
+
+void palSmat ( int n, float *a, float *y, float *d, int *jf, int *iw );
+
+void palSubet ( double rc, double dc, double eq,
+                double *rm, double *dm );
+
+void palSupgal ( double dsl, double dsb, double *dl, double *db );
+
+void palSvd ( int m, int n, int mp, int np,
+              double *a, double *w, double *v, double *work,
+              int *jstat );
+
+void palSvdcov ( int n, int np, int nc,
+                 double *w, double *v, double *work, double *cvm );
+
+void palSvdsol ( int m, int n, int mp, int np,
+                 double *b, double *u, double *w, double *v,
+                 double *work, double *x );
+
+void palTp2s ( float xi, float eta, float raz, float decz,
+               float *ra, float *dec );
+
+void palTp2v ( float xi, float eta, float v0[3], float v[3] );
+
+void palTps2c ( float xi, float eta, float ra, float dec,
+                float *raz1, float *decz1,
+                float *raz2, float *decz2, int *n );
+
+void palTpv2c ( float xi, float eta, float v[3],
+                float v01[3], float v02[3], int *n );
+
+void palUe2el ( const double u[13], int jformr,
+                int *jform, double *epoch, double *orbinc,
+                double *anode, double *perih, double *aorq, double *e,
+                double *aorl, double *dm, int *jstat );
+
+void palUe2pv ( double date, double u[13], double pv[], int *jstat );
+
+void palUnpcd ( double disco, double *x, double *y );
+
+void palV2tp ( float v[3], float v0[3], float *xi, float *eta, int *j );
+
+float palVdv ( float va[3], float vb[3] );
+
+int palVers ( char * verstring, size_t verlen );
+
+void palVn ( float v[3], float uv[3], float *vm );
+
+void palVxv ( float va[3], float vb[3], float vc[3] );
+
+void palXy2xy ( double x1, double y1, double coeffs[6],
+                double *x2, double *y2 );
+
+double palZd ( double ha, double dec, double phi );
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/ast/pal/pal1sofa.h b/ast/pal/pal1sofa.h
new file mode 100644
index 0000000..11f7446
--- /dev/null
+++ b/ast/pal/pal1sofa.h
@@ -0,0 +1,142 @@
+/*
+*+
+*  Name:
+*     pal1sofa.h
+
+*  Purpose:
+*     Mappings of ERFA names to SOFA names
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Include file
+
+*  Invocation:
+*     #include "pal1sofa.h"
+
+*  Description:
+*     PAL will work with both SOFA and ERFA libraries and the
+*     difference is generally a change in prefix. This include
+*     file maps the ERFA form of functions to the SOFA form
+*     and includes the relevant sofa.h vs erfa.h file.
+
+*  Authors:
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - PAL uses the ERFA form by default.
+
+*  History:
+*     2014-07-29 (TIMJ):
+*        Initial version
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2014 Tim Jenness
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#ifndef PAL1SOFAHDEF
+#define PAL1SOFAHDEF
+
+#if HAVE_CONFIG_H
+#  include <config.h>
+#endif
+
+# if HAVE_SOFA_H
+
+#  include "sofa.h"
+#  include "sofam.h"
+
+   /* Must replace ERFA with SOFA */
+
+#  define eraA2af iauA2af
+#  define eraA2tf iauA2tf
+#  define eraAf2a iauAf2a
+#  define eraAnp iauAnp
+#  define eraAnpm iauAnpm
+#  define eraC2s iauC2s
+#  define eraCal2jd iauCal2jd
+#  define eraD2tf iauD2tf
+#  define eraDat iauDat
+#  define eraEe06a iauEe06a
+#  define eraEpb iauEpb
+#  define eraEpb2jd iauEpb2jd
+#  define eraEpj iauEpj
+#  define eraEpj2jd iauEpj2jd
+#  define eraEpv00 iauEpv00
+#  define eraFk5hz iauFk5hz
+#  define eraGd2gc iauGd2gc
+#  define eraGmst06 iauGmst06
+#  define eraHfk5z iauHfk5z
+#  define eraIr iauIr
+#  define eraJd2cal iauJd2cal
+#  define eraNut06a iauNut06a
+#  define eraObl06 iauObl06
+#  define eraP06e iauP06e
+#  define eraPap iauPap
+#  define eraPas iauPas
+#  define eraPdp iauPdp
+#  define eraPlan94 iauPlan94
+#  define eraPmat06 iauPmat06
+#  define eraPn iauPn
+#  define eraPnm06a iauPnm06a
+#  define eraPxp iauPxp
+#  define eraRefco iauRefco
+#  define eraRm2v iauRm2v
+#  define eraRv2m iauRv2m
+#  define eraRx iauRx
+#  define eraRxp iauRxp
+#  define eraRxpv iauRxpv
+#  define eraRxr iauRxr
+#  define eraRy iauRy
+#  define eraRz iauRz
+#  define eraS2c iauS2c
+#  define eraSepp iauSepp
+#  define eraSeps iauSeps
+#  define eraStarpm iauStarpm
+#  define eraTf2a iauTf2a
+#  define eraTf2d iauTf2d
+#  define eraTr iauTr
+#  define eraTrxp iauTrxp
+
+/* These are from sofam.h */
+
+#  define ERFA_WGS84 WGS84
+
+#  define ERFA_DJ00 DJ00
+#  define ERFA_DJY DJY
+#  define ERFA_DAU DAU
+
+# else
+
+#  include "erfa.h"
+#  include "erfam.h"
+
+/* No further action required */
+
+# endif
+
+#endif
diff --git a/ast/pal/palAddet.c b/ast/pal/palAddet.c
new file mode 100644
index 0000000..ce3ddab
--- /dev/null
+++ b/ast/pal/palAddet.c
@@ -0,0 +1,112 @@
+/*
+*+
+*  Name:
+*     palAddet
+
+*  Purpose:
+*     Add the E-terms to a pre IAU 1976 mean place
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palAddet ( double rm, double dm, double eq,
+*                     double *rc, double *dc );
+
+*  Arguments:
+*     rm = double (Given)
+*        RA without E-terms (radians)
+*     dm = double (Given)
+*        Dec without E-terms (radians)
+*     eq = double (Given)
+*        Besselian epoch of mean equator and equinox
+*     rc = double * (Returned)
+*        RA with E-terms included (radians)
+*     dc = double * (Returned)
+*        Dec with E-terms included (radians)
+
+*  Description:
+*     Add the E-terms (elliptic component of annual aberration)
+*     to a pre IAU 1976 mean place to conform to the old
+*     catalogue convention.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     Most star positions from pre-1984 optical catalogues (or
+*     derived from astrometry using such stars) embody the
+*     E-terms.  If it is necessary to convert a formal mean
+*     place (for example a pulsar timing position) to one
+*     consistent with such a star catalogue, then the RA,Dec
+*     should be adjusted using this routine.
+
+*  See Also:
+*     Explanatory Supplement to the Astronomical Ephemeris,
+*     section 2D, page 48.
+
+*  History:
+*     2012-02-12(TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1999 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palAddet ( double rm, double dm, double eq, double *rc, double *dc ) {
+  double a[3];   /* The E-terms */
+  double v[3];
+  int i;
+
+  /* Note the preference for IAU routines */
+
+  /* Retrieve the E-terms */
+  palEtrms( eq, a );
+
+  /* Spherical to Cartesian */
+  eraS2c( rm, dm, v );
+
+  /* Include the E-terms */
+  for (i=0; i<3; i++) {
+    v[i] += a[i];
+  }
+
+  /* Cartesian to spherical */
+  eraC2s( v, rc, dc );
+
+  /* Bring RA into conventional range */
+  *rc = eraAnp( *rc );
+
+}
diff --git a/ast/pal/palAmpqk.c b/ast/pal/palAmpqk.c
new file mode 100644
index 0000000..ce698b8
--- /dev/null
+++ b/ast/pal/palAmpqk.c
@@ -0,0 +1,159 @@
+/*
+*+
+*  Name:
+*     palAmpqk
+
+*  Purpose:
+*     Convert star RA,Dec from geocentric apparent to mean place.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palAmpqk ( double ra, double da, double amprms[21],
+*                     double *rm, double *dm )
+
+*  Arguments:
+*     ra = double (Given)
+*        Apparent RA (radians).
+*     da = double (Given)
+*        Apparent Dec (radians).
+*     amprms = double[21] (Given)
+*        Star-independent mean-to-apparent parameters (see palMappa):
+*        (0)      time interval for proper motion (Julian years)
+*        (1-3)    barycentric position of the Earth (AU)
+*        (4-6)    heliocentric direction of the Earth (unit vector)
+*        (7)      (grav rad Sun)*2/(Sun-Earth distance)
+*        (8-10)   abv: barycentric Earth velocity in units of c
+*        (11)     sqrt(1-v*v) where v=modulus(abv)
+*        (12-20)  precession/nutation (3,3) matrix
+*     rm = double (Returned)
+*        Mean RA (radians).
+*     dm = double (Returned)
+*        Mean Dec (radians).
+
+*  Description:
+*     Convert star RA,Dec from geocentric apparent to mean place. The "mean"
+*     coordinate system is in fact close to ICRS. Use of this function
+*     is appropriate when efficiency is important and where many star
+*     positions are all to be transformed for one epoch and equinox.  The
+*     star-independent parameters can be obtained by calling the palMappa
+*     function.
+
+*  Note:
+*     Iterative techniques are used for the aberration and
+*     light deflection corrections so that the routines
+*     palAmp (or palAmpqk) and palMap (or palMapqk) are
+*     accurate inverses;  even at the edge of the Sun's disc
+*     the discrepancy is only about 1 nanoarcsecond.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-13 (PTW):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     2016-12-19 (TIMJ):
+*        Add in light deflection (was missed in the initial port).
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2000 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     Copyright (C) 2016 Tim Jenness
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palAmpqk ( double ra, double da, double amprms[21], double *rm,
+                double *dm ){
+
+/* Local Variables: */
+   double ab1;                 /* sqrt(1-v*v) where v=modulus of Earth vel */
+   double abv[3];              /* Earth velocity wrt SSB (c, FK5) */
+   double p1[3], p2[3], p3[3]; /* work vectors */
+   double ab1p1, p1dv, p1dvp1, w;
+   double gr2e, pde, pdep1, ehn[3], p[3];
+   int i, j;
+
+/* Unpack some of the parameters */
+   gr2e = amprms[7];
+   ab1  = amprms[11];
+   for( i = 0; i < 3; i++ ) {
+      ehn[i] = amprms[i + 4];
+      abv[i] = amprms[i + 8];
+   }
+
+/* Apparent RA,Dec to Cartesian */
+   eraS2c( ra, da, p3 );
+
+/* Precession and nutation */
+   eraTrxp( (double(*)[3]) &amprms[12], p3, p2 );
+
+/* Aberration */
+   ab1p1 = ab1 + 1.0;
+   for( i = 0; i < 3; i++ ) {
+      p1[i] = p2[i];
+   }
+   for( j = 0; j < 2; j++ ) {
+      p1dv = eraPdp( p1, abv );
+      p1dvp1 = 1.0 + p1dv;
+      w = 1.0 + p1dv / ab1p1;
+      for( i = 0; i < 3; i++ ) {
+         p1[i] = ( p1dvp1 * p2[i] - w * abv[i] ) / ab1;
+      }
+      eraPn( p1, &w, p3 );
+      for( i = 0; i < 3; i++ ) {
+         p1[i] = p3[i];
+      }
+   }
+
+/* Light deflection */
+   for( i = 0; i < 3; i++ ) {
+       p[i] = p1[i];
+   }
+   for( j = 0; j < 5; j++ ) {
+      pde = eraPdp( p, ehn );
+      pdep1 = 1.0 + pde;
+      w = pdep1 - gr2e*pde;
+      for( i = 0; i < 3; i++ ) {
+         p[i] = (pdep1*p1[i] - gr2e*ehn[i])/w;
+      }
+      eraPn( p, &w, p2 );
+      for( i = 0; i < 3; i++ ) {
+         p[i] = p2[i];
+      }
+   }
+
+/* Mean RA,Dec */
+   eraC2s( p, rm, dm );
+   *rm = eraAnp( *rm );
+}
diff --git a/ast/pal/palCaldj.c b/ast/pal/palCaldj.c
new file mode 100644
index 0000000..03c6b97
--- /dev/null
+++ b/ast/pal/palCaldj.c
@@ -0,0 +1,99 @@
+/*
+*+
+*  Name:
+*     palCaldj
+
+*  Purpose:
+*     Gregorian Calendar to Modified Julian Date
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palCaldj ( int iy, int im, int id, double *djm, int *j );
+
+*  Arguments:
+*     iy = int (Given)
+*        Year in the Gregorian calendar
+*     im = int (Given)
+*        Month in the Gergorian calendar
+*     id = int (Given)
+*        Day in the Gregorian calendar
+*     djm = double * (Returned)
+*        Modified Julian Date (JD-2400000.5) for 0 hrs
+*     j = status (Returned)
+*       0 = OK. See eraCal2jd for other values.
+
+*  Description:
+*     Modified Julian Date to Gregorian Calendar with special
+*     behaviour for 2-digit years relating to 1950 to 2049.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-11 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Notes:
+*     - Uses eraCal2jd
+*     - Unlike eraCal2jd this routine treats the years 0-100 as
+*       referring to the end of the 20th Century and beginning of
+*       the 21st Century. If this behaviour is not acceptable
+*       use the SOFA/ERFA routine directly or palCldj.
+*       Acceptable years are 00-49, interpreted as 2000-2049,
+*                            50-99,     "       "  1950-1999,
+*                            all others, interpreted literally.
+*     - Unlike SLA this routine will work with negative years.
+
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+void palCaldj ( int iy, int im, int id, double *djm, int *j ) {
+  int adj = 0;   /* Year adjustment */
+  double djm0;
+
+  if (iy >= 0 && iy <= 49) {
+    adj = 2000;
+  } else if (iy >= 50 && iy <= 99) {
+    adj = 1900;
+  }
+  iy += adj;
+
+  *j = eraCal2jd( iy, im, id, &djm0, djm );
+}
diff --git a/ast/pal/palDat.c b/ast/pal/palDat.c
new file mode 100644
index 0000000..f869681
--- /dev/null
+++ b/ast/pal/palDat.c
@@ -0,0 +1,95 @@
+/*
+*+
+*  Name:
+*     palDat
+
+*  Purpose:
+*     Return offset between UTC and TAI
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     dat = palDat( double utc );
+
+*  Arguments:
+*     utc = double (Given)
+*        UTC date as a modified JD (JD-2400000.5)
+
+*  Returned Value:
+*     dat = double
+*        TAI-UTC in seconds
+
+*  Description:
+*     Increment to be applied to Coordinated Universal Time UTC to give
+*     International Atomic Time (TAI).
+
+*  Authors:
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - This routine converts the MJD argument to calendar date before calling
+*       the SOFA/ERFA eraDat function.
+*     - This routine matches the slaDat interface which differs from the eraDat
+*       interface. Consider coding directly to the SOFA/ERFA interface.
+*     - See eraDat for a description of error conditions when calling this function
+*       with a time outside of the UTC range.
+*     - The status argument from eraDat is ignored. This is reasonable since the
+*       error codes are mainly related to incorrect calendar dates when calculating
+*       the JD internally.
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version
+*        Adapted with permission from the Fortran SLALIB library
+*        although the core algorithm is now from SOFA.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+
+#include "pal1sofa.h"
+
+double palDat ( double dju ) {
+  int iy;
+  int im;
+  int id;
+  int status;
+  double fd;
+  double deltat;
+
+  eraJd2cal( PAL__MJD0, dju,
+             &iy, &im, &id, &fd );
+
+  status = eraDat( iy, im, id, fd, &deltat );
+  return deltat;
+}
diff --git a/ast/pal/palDe2h.c b/ast/pal/palDe2h.c
new file mode 100644
index 0000000..a250e9e
--- /dev/null
+++ b/ast/pal/palDe2h.c
@@ -0,0 +1,142 @@
+/*
+*+
+*  Name:
+*     palDe2h
+
+*  Purpose:
+*     Equatorial to horizon coordinates: HA,Dec to Az,E
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDe2h( double ha, double dec, double phi, double * az, double * el );
+
+*  Arguments:
+*     ha = double * (Given)
+*        Hour angle (radians)
+*     dec = double * (Given)
+*        Declination (radians)
+*     phi = double (Given)
+*        Observatory latitude (radians)
+*     az = double * (Returned)
+*        Azimuth (radians)
+*     el = double * (Returned)
+*        Elevation (radians)
+
+*  Description:
+*     Convert equatorial to horizon coordinates.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - All the arguments are angles in radians.
+*     - Azimuth is returned in the range 0-2pi;  north is zero,
+*       and east is +pi/2.  Elevation is returned in the range
+*       +/-pi/2.
+*     - The latitude must be geodetic.  In critical applications,
+*       corrections for polar motion should be applied.
+*     - In some applications it will be important to specify the
+*       correct type of hour angle and declination in order to
+*       produce the required type of azimuth and elevation.  In
+*       particular, it may be important to distinguish between
+*       elevation as affected by refraction, which would
+*       require the "observed" HA,Dec, and the elevation
+*       in vacuo, which would require the "topocentric" HA,Dec.
+*       If the effects of diurnal aberration can be neglected, the
+*       "apparent" HA,Dec may be used instead of the topocentric
+*       HA,Dec.
+*     - No range checking of arguments is carried out.
+*     - In applications which involve many such calculations, rather
+*       than calling the present routine it will be more efficient to
+*       use inline code, having previously computed fixed terms such
+*       as sine and cosine of latitude, and (for tracking a star)
+*       sine and cosine of declination.
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include <math.h>
+
+void
+palDe2h ( double ha, double dec, double phi, double *az, double *el) {
+
+  double sh;
+  double ch;
+  double sd;
+  double cd;
+  double sp;
+  double cp;
+
+  double a;
+
+  double x;
+  double y;
+  double z;
+  double r;
+
+  /*  Useful trig functions */
+  sh = sin(ha);
+  ch = cos(ha);
+  sd = sin(dec);
+  cd = cos(dec);
+  sp = sin(phi);
+  cp = cos(phi);
+
+  /*  Az,El as x,y,z */
+  x = -ch * cd * sp + sd * cp;
+  y = -sh * cd;
+  z = ch * cd * cp + sd * sp;
+
+  /*  To spherical */
+  r = sqrt(x * x + y * y);
+  if (r == 0.) {
+    a = 0.;
+  } else {
+    a = atan2(y, x);
+  }
+  if (a < 0.) {
+    a += PAL__D2PI;
+  }
+  *az = a;
+  *el = atan2(z, r);
+
+  return;
+}
diff --git a/ast/pal/palDeuler.c b/ast/pal/palDeuler.c
new file mode 100644
index 0000000..17e536b
--- /dev/null
+++ b/ast/pal/palDeuler.c
@@ -0,0 +1,141 @@
+/*
+*+
+*  Name:
+*     palDeuler
+
+*  Purpose:
+*     Form a rotation matrix from the Euler angles
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palDeuler ( const char *order, double phi, double theta, double psi,
+*                      double rmat[3][3] );
+
+*  Arguments:
+*     order = const char[] (Given)
+*        Specifies about which axes the rotation occurs
+*     phi = double (Given)
+*        1st rotation (radians)
+*     theta = double (Given)
+*        2nd rotation (radians)
+*     psi = double (Given)
+*        3rd rotation (radians)
+*     rmat = double[3][3] (Given & Returned)
+*        Rotation matrix
+
+*  Description:
+*     A rotation is positive when the reference frame rotates
+*     anticlockwise as seen looking towards the origin from the
+*     positive region of the specified axis.
+*
+*     The characters of ORDER define which axes the three successive
+*     rotations are about.  A typical value is 'ZXZ', indicating that
+*     RMAT is to become the direction cosine matrix corresponding to
+*     rotations of the reference frame through PHI radians about the
+*     old Z-axis, followed by THETA radians about the resulting X-axis,
+*     then PSI radians about the resulting Z-axis.
+*
+*     The axis names can be any of the following, in any order or
+*     combination:  X, Y, Z, uppercase or lowercase, 1, 2, 3.  Normal
+*     axis labelling/numbering conventions apply;  the xyz (=123)
+*     triad is right-handed.  Thus, the 'ZXZ' example given above
+*     could be written 'zxz' or '313' (or even 'ZxZ' or '3xZ').  ORDER
+*     is terminated by length or by the first unrecognized character.
+*
+*     Fewer than three rotations are acceptable, in which case the later
+*     angle arguments are ignored.  If all rotations are zero, the
+*     identity matrix is produced.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1997 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void
+palDeuler( const char *order, double phi, double theta, double psi,
+                 double rmat[3][3] ) {
+  int i = 0;
+  double rotations[3];
+
+  /* Initialise rmat */
+  eraIr( rmat );
+
+  /* copy the rotations into an array */
+  rotations[0] = phi;
+  rotations[1] = theta;
+  rotations[2] = psi;
+
+  /* maximum three rotations */
+  while (i < 3 && order[i] != '\0') {
+
+    switch (order[i]) {
+    case 'X':
+    case 'x':
+    case '1':
+      eraRx( rotations[i], rmat );
+      break;
+
+    case 'Y':
+    case 'y':
+    case '2':
+      eraRy( rotations[i], rmat );
+      break;
+
+    case 'Z':
+    case 'z':
+    case '3':
+      eraRz( rotations[i], rmat );
+      break;
+
+    default:
+      /* break out the loop if we do not recognize something */
+      i = 3;
+
+    }
+
+    /* Go to the next position */
+    i++;
+  }
+
+  return;
+}
diff --git a/ast/pal/palDh2e.c b/ast/pal/palDh2e.c
new file mode 100644
index 0000000..65434b7
--- /dev/null
+++ b/ast/pal/palDh2e.c
@@ -0,0 +1,133 @@
+/*
+*+
+*  Name:
+*     palDh2e
+
+*  Purpose:
+*     Horizon to equatorial coordinates: Az,El to HA,Dec
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDh2e( double az, double el, double phi, double * ha, double * dec );
+
+*  Arguments:
+*     az = double (Given)
+*        Azimuth (radians)
+*     el = double (Given)
+*        Elevation (radians)
+*     phi = double (Given)
+*        Observatory latitude (radians)
+*     ha = double * (Returned)
+*        Hour angle (radians)
+*     dec = double * (Returned)
+*        Declination (radians)
+
+*  Description:
+*     Convert horizon to equatorial coordinates.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - All the arguments are angles in radians.
+*     - The sign convention for azimuth is north zero, east +pi/2.
+*     - HA is returned in the range +/-pi.  Declination is returned
+*       in the range +/-pi/2.
+*     - The latitude is (in principle) geodetic.  In critical
+*       applications, corrections for polar motion should be applied.
+*     - In some applications it will be important to specify the
+*       correct type of elevation in order to produce the required
+*       type of HA,Dec.  In particular, it may be important to
+*       distinguish between the elevation as affected by refraction,
+*       which will yield the "observed" HA,Dec, and the elevation
+*       in vacuo, which will yield the "topocentric" HA,Dec.  If the
+*       effects of diurnal aberration can be neglected, the
+*       topocentric HA,Dec may be used as an approximation to the
+*       "apparent" HA,Dec.
+*     - No range checking of arguments is done.
+*     - In applications which involve many such calculations, rather
+*       than calling the present routine it will be more efficient to
+*       use inline code, having previously computed fixed terms such
+*       as sine and cosine of latitude.
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1996 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include <math.h>
+
+void
+palDh2e ( double az, double el, double phi, double *ha, double *dec) {
+
+  double sa;
+  double ca;
+  double se;
+  double ce;
+  double sp;
+  double cp;
+
+  double x;
+  double y;
+  double z;
+  double r;
+
+  /*  Useful trig functions */
+  sa = sin(az);
+  ca = cos(az);
+  se = sin(el);
+  ce = cos(el);
+  sp = sin(phi);
+  cp = cos(phi);
+
+  /*  HA,Dec as x,y,z */
+  x = -ca * ce * sp + se * cp;
+  y = -sa * ce;
+  z = ca * ce * cp + se * sp;
+
+  /*  To HA,Dec */
+  r = sqrt(x * x + y * y);
+  if (r == 0.) {
+    *ha = 0.;
+  } else {
+    *ha = atan2(y, x);
+  }
+  *dec = atan2(z, r);
+
+  return;
+}
diff --git a/ast/pal/palDjcal.c b/ast/pal/palDjcal.c
new file mode 100644
index 0000000..b6aac9e
--- /dev/null
+++ b/ast/pal/palDjcal.c
@@ -0,0 +1,97 @@
+/*
+*+
+*  Name:
+*     palDjcal
+
+*  Purpose:
+*     Modified Julian Date to Gregorian Calendar
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palDjcal ( int ndp, double djm, int iymdf[4], int *j );
+
+*  Arguments:
+*     ndp = int (Given)
+*        Number of decimal places of days in fraction.
+*     djm = double (Given)
+*        Modified Julian Date (JD-2400000.5)
+*     iymdf[4] = int[] (Returned)
+*       Year, month, day, fraction in Gregorian calendar.
+*     j = status (Returned)
+*       0 = OK. See eraJd2cal for other values.
+
+*  Description:
+*     Modified Julian Date to Gregorian Calendar, expressed
+*     in a form convenient for formatting messages (namely
+*     rounded to a specified precision, and with the fields
+*     stored in a single array)
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-10 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Notes:
+*     - Uses eraJd2cal
+
+*  Copyright:
+*     Copyright (C) 2004 Patrick T. Wallace
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include <math.h>
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+void palDjcal ( int ndp, double djm, int iymdf[4], int *j ) {
+  double frac = 0.0;
+  double nfd;
+
+  *j = eraJd2cal( PAL__MJD0, djm, &(iymdf[0]),
+		  &(iymdf[1]), &(iymdf[2]),
+		  &frac);
+
+  /* Convert ndp to a power of 10 */
+  nfd = pow( 10., (double)ndp );
+
+  /* Multiply the fraction */
+  frac *= nfd;
+
+  /* and now we want to round to the nearest integer */
+  iymdf[3] = (int)DNINT(frac);
+
+}
diff --git a/ast/pal/palDmat.c b/ast/pal/palDmat.c
new file mode 100644
index 0000000..2948f92
--- /dev/null
+++ b/ast/pal/palDmat.c
@@ -0,0 +1,182 @@
+/*
+*+
+*  Name:
+*     palDmat
+
+*  Purpose:
+*     Matrix inversion & solution of simultaneous equations
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palDmat( int n, double *a, double *y, double *d, int *jf,
+*                    int *iw );
+
+*  Arguments:
+*     n = int (Given)
+*        Number of simultaneous equations and number of unknowns.
+*     a = double[] (Given & Returned)
+*        A non-singular NxN matrix (implemented as a contiguous block
+*        of memory). After calling this routine "a" contains the
+*        inverse of the matrix.
+*     y = double[] (Given & Returned)
+*        On input the vector of N knowns. On exit this vector contains the
+*        N solutions.
+*     d = double * (Returned)
+*        The determinant.
+*     jf = int * (Returned)
+*        The singularity flag.  If the matrix is non-singular, jf=0
+*        is returned.  If the matrix is singular, jf=-1 & d=0.0 are
+*        returned.  In the latter case, the contents of array "a" on
+*        return are undefined.
+*     iw = int[] (Given)
+*        Integer workspace of size N.
+
+*  Description:
+*     Matrix inversion & solution of simultaneous equations
+*     For the set of n simultaneous equations in n unknowns:
+*          A.Y = X
+*     this routine calculates the inverse of A, the determinant
+*     of matrix A and the vector of N unknowns.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-11 (TIMJ):
+*        Combination of a port of the Fortran and a comparison
+*        with the obfuscated GPL C routine.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Notes:
+*     - Implemented using Gaussian elimination with partial pivoting.
+*     - Optimized for speed rather than accuracy with errors 1 to 4
+*       times those of routines optimized for accuracy.
+
+*  Copyright:
+*     Copyright (C) 2001 Rutherford Appleton Laboratory.
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+
+void palDmat ( int n, double *a, double *y, double *d, int *jf, int *iw ) {
+
+  const double SFA = 1e-20;
+
+  int k;
+  double*aoff;
+
+  *jf=0;
+  *d=1.0;
+  for(k=0,aoff=a; k<n; k++, aoff+=n){
+    int imx;
+    double * aoff2 = aoff;
+    double amx=fabs(aoff[k]);
+    imx=k;
+    if(k!=n){
+      int i;
+      double *apos2;
+      for(i=k+1,apos2=aoff+n;i<n;i++,apos2+=n){
+	double t=fabs(apos2[k]);
+	if(t>amx){
+	  amx=t;
+	  imx=i;
+	  aoff2=apos2;
+	}
+      }
+    }
+    if(amx<SFA){
+      *jf=-1;
+    } else {
+      if(imx!=k){
+	double t;
+	int j;
+	for(j=0;j<n;j++){
+	  t=aoff[j];
+	  aoff[j]=aoff2[j];
+	  aoff2[j]=t;
+	}
+	t=y[k];
+	y[k]=y[imx];
+	y[imx]=t;*d=-*d;
+      }
+      iw[k]=imx;
+      *d*=aoff[k];
+      if(fabs(*d)<SFA){
+	*jf=-1;
+      } else {
+	double yk;
+	double * apos2;
+	int i, j;
+	aoff[k]=1.0/aoff[k];
+	for(j=0;j<n;j++){
+	  if(j!=k){
+	    aoff[j]*=aoff[k];
+	  }
+	}
+	yk=y[k]*aoff[k];
+	y[k]=yk;
+	for(i=0,apos2=a;i<n;i++,apos2+=n){
+	  if(i!=k){
+	    for(j=0;j<n;j++){
+	      if(j!=k){
+		apos2[j]-=apos2[k]*aoff[j];
+	      }
+	    }
+	    y[i]-=apos2[k]*yk;
+	  }
+	}
+	for(i=0,apos2=a;i<n;i++,apos2+=n){
+	  if(i!=k){
+	    apos2[k]*=-aoff[k];
+	  }
+	}
+      }
+    }
+  }
+  if(*jf!=0){
+    *d=0.0;
+  } else {
+    for(k=n;k-->0;){
+      int ki=iw[k];
+      if(k!=ki){
+	int i;
+	double *apos = a;
+	for(i=0;i<n;i++,apos+=n){
+	  double t=apos[k];
+	  apos[k]=apos[ki];
+	  apos[ki]=t;
+	}
+      }
+    }
+  }
+}
diff --git a/ast/pal/palDrange.c b/ast/pal/palDrange.c
new file mode 100644
index 0000000..804ed55
--- /dev/null
+++ b/ast/pal/palDrange.c
@@ -0,0 +1,77 @@
+/*
+*+
+*  Name:
+*     palDrange
+
+*  Purpose:
+*     Normalize angle into range +/- pi
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDrange( double angle )
+
+*  Arguments:
+*     angle = double (Given)
+*        The angle in radians.
+
+*  Description:
+*     The result is "angle" expressed in the range +/- pi. If the
+*     supplied value for "angle" is equal to +/- pi, it is returned
+*     unchanged.
+
+*  Authors:
+*     DSB: David S Berry (JAC, Hawaii)
+*     PTW: Patrick T. Wallace
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-05-09 (DSB):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include <math.h>
+
+double palDrange( double angle ){
+   double result = fmod( angle, PAL__D2PI );
+   if( result > PAL__DPI ) {
+      result -= PAL__D2PI;
+   } else if( result < -PAL__DPI ) {
+      result += PAL__D2PI;
+   }
+   return result;
+}
+
diff --git a/ast/pal/palDs2tp.c b/ast/pal/palDs2tp.c
new file mode 100644
index 0000000..cbf582a
--- /dev/null
+++ b/ast/pal/palDs2tp.c
@@ -0,0 +1,127 @@
+/*
+*+
+*  Name:
+*     palDs2tp
+
+*  Purpose:
+*     Spherical to tangent plane projection
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDs2tp( double ra, double dec, double raz, double decz,
+*               double *xi, double *eta, int *j );
+
+*  Arguments:
+*     ra = double (Given)
+*        RA spherical coordinate of point to be projected (radians)
+*     dec = double (Given)
+*        Dec spherical coordinate of point to be projected (radians)
+*     raz = double (Given)
+*        RA spherical coordinate of tangent point (radians)
+*     decz = double (Given)
+*        Dec spherical coordinate of tangent point (radians)
+*     xi = double * (Returned)
+*        First rectangular coordinate on tangent plane (radians)
+*     eta = double * (Returned)
+*        Second rectangular coordinate on tangent plane (radians)
+*     j = int * (Returned)
+*        status: 0 = OK, star on tangent plane
+*                1 = error, star too far from axis
+*                2 = error, antistar on tangent plane
+*                3 = error, antistar too far from axis
+
+*  Description:
+*     Projection of spherical coordinates onto tangent plane:
+*     "gnomonic" projection - "standard coordinates"
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1996 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include <math.h>
+
+void
+palDs2tp ( double ra, double dec, double raz, double decz,
+           double *xi, double *eta, int *j ) {
+
+  const double TINY = 1.0e-6;
+
+  double cdec;
+  double sdec;
+  double radif;
+  double cdecz;
+  double denom;
+  double sdecz;
+  double cradif;
+  double sradif;
+
+  /*  Trig functions */
+  sdecz = sin(decz);
+  sdec = sin(dec);
+  cdecz = cos(decz);
+  cdec = cos(dec);
+  radif = ra - raz;
+  sradif = sin(radif);
+  cradif = cos(radif);
+
+  /*  Reciprocal of star vector length to tangent plane */
+  denom = sdec * sdecz + cdec * cdecz * cradif;
+
+  /*  Handle vectors too far from axis */
+  if (denom > TINY) {
+    *j = 0;
+  } else if (denom >= 0.) {
+    *j = 1;
+    denom = TINY;
+  } else if (denom > -TINY) {
+    *j = 2;
+    denom = -TINY;
+  } else {
+    *j = 3;
+  }
+
+  /*  Compute tangent plane coordinates (even in dubious cases) */
+  *xi = cdec * sradif / denom;
+  *eta = (sdec * cdecz - cdec * sdecz * cradif) / denom;
+
+  return;
+}
diff --git a/ast/pal/palDtp2s.c b/ast/pal/palDtp2s.c
new file mode 100644
index 0000000..583a56d
--- /dev/null
+++ b/ast/pal/palDtp2s.c
@@ -0,0 +1,95 @@
+/*
+*+
+*  Name:
+*     palDtp2s
+
+*  Purpose:
+*     Tangent plane to spherical coordinates
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDtp2s( double xi, double eta, double raz, double decz,
+*               double *ra, double *dec);
+
+*  Arguments:
+*     xi = double (Given)
+*        First rectangular coordinate on tangent plane (radians)
+*     eta = double (Given)
+*        Second rectangular coordinate on tangent plane (radians)
+*     raz = double (Given)
+*        RA spherical coordinate of tangent point (radians)
+*     decz = double (Given)
+*        Dec spherical coordinate of tangent point (radians)
+*     ra = double * (Returned)
+*        RA spherical coordinate of point to be projected (radians)
+*     dec = double * (Returned)
+*        Dec spherical coordinate of point to be projected (radians)
+
+*  Description:
+*     Transform tangent plane coordinates into spherical.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+#include <math.h>
+
+void
+palDtp2s ( double xi, double eta, double raz, double decz,
+           double *ra, double *dec ) {
+
+  double cdecz;
+  double denom;
+  double sdecz;
+  double d;
+
+  sdecz = sin(decz);
+  cdecz = cos(decz);
+  denom = cdecz - eta * sdecz;
+  d = atan2(xi, denom) + raz;
+  *ra = eraAnp(d);
+  *dec = atan2(sdecz + eta * cdecz, sqrt(xi * xi + denom * denom));
+
+  return;
+}
diff --git a/ast/pal/palDtps2c.c b/ast/pal/palDtps2c.c
new file mode 100644
index 0000000..ecb3090
--- /dev/null
+++ b/ast/pal/palDtps2c.c
@@ -0,0 +1,151 @@
+/*
+*+
+*  Name:
+*     palDtps2c
+
+*  Purpose:
+*     Determine RA,Dec of tangent point from coordinates
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDtps2c( double xi, double eta, double ra, double dec,
+*                double * raz1, double decz1,
+*                double * raz2, double decz2, int *n);
+
+*  Arguments:
+*     xi = double (Given)
+*        First rectangular coordinate on tangent plane (radians)
+*     eta = double (Given)
+*        Second rectangular coordinate on tangent plane (radians)
+*     ra = double (Given)
+*        RA spherical coordinate of star (radians)
+*     dec = double (Given)
+*        Dec spherical coordinate of star (radians)
+*     raz1 = double * (Returned)
+*        RA spherical coordinate of tangent point, solution 1 (radians)
+*     decz1 = double * (Returned)
+*        Dec spherical coordinate of tangent point, solution 1 (radians)
+*     raz2 = double * (Returned)
+*        RA spherical coordinate of tangent point, solution 2 (radians)
+*     decz2 = double * (Returned)
+*        Dec spherical coordinate of tangent point, solution 2 (radians)
+*     n = int * (Returned)
+*        number of solutions: 0 = no solutions returned (note 2)
+*                             1 = only the first solution is useful (note 3)
+*                             2 = both solutions are useful (note 3)
+
+
+*  Description:
+*     From the tangent plane coordinates of a star of known RA,Dec,
+*     determine the RA,Dec of the tangent point.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - The RAZ1 and RAZ2 values are returned in the range 0-2pi.
+*     - Cases where there is no solution can only arise near the poles.
+*       For example, it is clearly impossible for a star at the pole
+*       itself to have a non-zero XI value, and hence it is
+*       meaningless to ask where the tangent point would have to be
+*       to bring about this combination of XI and DEC.
+*     - Also near the poles, cases can arise where there are two useful
+*       solutions.  The argument N indicates whether the second of the
+*       two solutions returned is useful.  N=1 indicates only one useful
+*       solution, the usual case;  under these circumstances, the second
+*       solution corresponds to the "over-the-pole" case, and this is
+*       reflected in the values of RAZ2 and DECZ2 which are returned.
+*     - The DECZ1 and DECZ2 values are returned in the range +/-pi, but
+*       in the usual, non-pole-crossing, case, the range is +/-pi/2.
+*     - This routine is the spherical equivalent of the routine sla_DTPV2C.
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+#include <math.h>
+
+void
+palDtps2c( double xi, double eta, double ra, double dec,
+           double * raz1, double * decz1,
+           double * raz2, double * decz2, int *n) {
+
+  double x2;
+  double y2;
+  double sd;
+  double cd;
+  double sdf;
+  double r2;
+
+  x2 = xi * xi;
+  y2 = eta * eta;
+  sd = sin(dec);
+  cd = cos(dec);
+  sdf = sd * sqrt(x2 + 1. + y2);
+  r2 = cd * cd * (y2 + 1.) - sd * sd * x2;
+  if (r2 >= 0.) {
+    double r;
+    double s;
+    double c;
+
+    r = sqrt(r2);
+    s = sdf - eta * r;
+    c = sdf * eta + r;
+    if (xi == 0. && r == 0.) {
+      r = 1.;
+    }
+    *raz1 = eraAnp(ra - atan2(xi, r));
+    *decz1 = atan2(s, c);
+    r = -r;
+    s = sdf - eta * r;
+    c = sdf * eta + r;
+    *raz2 = eraAnp(ra - atan2(xi, r));
+    *decz2 = atan2(s, c);
+    if (fabs(sdf) < 1.) {
+      *n = 1;
+    } else {
+      *n = 2;
+    }
+  } else {
+    *n = 0;
+  }
+  return;
+}
diff --git a/ast/pal/palDtt.c b/ast/pal/palDtt.c
new file mode 100644
index 0000000..f6a3714
--- /dev/null
+++ b/ast/pal/palDtt.c
@@ -0,0 +1,77 @@
+/*
+*+
+*  Name:
+*     palDtt
+
+*  Purpose:
+*     Return offset between UTC and TT
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     dtt = palDtt( double utc );
+
+*  Arguments:
+*     utc = double (Given)
+*        UTC date as a modified JD (JD-2400000.5)
+
+*  Returned Value:
+*     dtt = double
+*        TT-UTC in seconds
+
+*  Description:
+*     Increment to be applied to Coordinated Universal Time UTC to give
+*     Terrestrial Time TT (formerly Ephemeris Time ET)
+
+*  Authors:
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     PTW: Patrick T. Wallace
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - Consider a comprehensive upgrade to use the time transformations in SOFA's time
+*       cookbook:  http://www.iausofa.org/sofa_ts_c.pdf.
+*     - See eraDat for a description of error conditions when calling this function
+*       with a time outside of the UTC range. This behaviour differs from slaDtt.
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+
+double palDtt( double utc ) {
+  return 32.184 + palDat( utc );
+}
diff --git a/ast/pal/palEcmat.c b/ast/pal/palEcmat.c
new file mode 100644
index 0000000..e9d9aeb
--- /dev/null
+++ b/ast/pal/palEcmat.c
@@ -0,0 +1,82 @@
+/*
+*+
+*  Name:
+*     palEcmat
+
+*  Purpose:
+*     Form the equatorial to ecliptic rotation matrix - IAU 2006
+*     precession model.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palEcmat( double date, double rmat[3][3] )
+
+*  Arguments:
+*     date = double (Given)
+*        TT as Modified Julian Date (JD-2400000.5). The difference
+*        between TT and TDB is of the order of a millisecond or two
+*        (i.e. about 0.02 arc-seconds).
+*     rmat = double[3][3] (Returned)
+*        Rotation matrix
+
+*  Description:
+*     The equatorial to ecliptic rotation matrix is found and returned.
+*     The matrix is in the sense   V(ecl)  =  RMAT * V(equ);  the
+*     equator, equinox and ecliptic are mean of date.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-10 (DSB):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1996 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+void palEcmat( double date, double rmat[3][3] ) {
+
+/* Mean obliquity (the angle between the ecliptic and mean equator of
+   date). */
+   double eps0 = eraObl06( PAL__MJD0, date );
+
+/* Matrix */
+   palDeuler( "X", eps0, 0.0, 0.0, rmat );
+
+}
diff --git a/ast/pal/palEqgal.c b/ast/pal/palEqgal.c
new file mode 100644
index 0000000..9df3d09
--- /dev/null
+++ b/ast/pal/palEqgal.c
@@ -0,0 +1,118 @@
+/*
+*+
+*  Name:
+*     palEqgal
+
+*  Purpose:
+*     Convert from J2000.0 equatorial coordinates to Galactic
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palEqgal ( double dr, double dd, double *dl, double *db );
+
+*  Arguments:
+*     dr = double (Given)
+*       J2000.0 RA (radians)
+*     dd = double (Given)
+*       J2000.0 Dec (radians
+*     dl = double * (Returned)
+*       Galactic longitude (radians).
+*     db = double * (Returned)
+*       Galactic latitude (radians).
+
+*  Description:
+*     Transformation from J2000.0 equatorial coordinates
+*     to IAU 1958 galactic coordinates.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     The equatorial coordinates are J2000.0.  Use the routine
+*     palGe50 if conversion to B1950.0 'FK4' coordinates is
+*     required.
+
+*  See Also:
+*     Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960)
+
+*  History:
+*     2012-02-12(TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1998 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palEqgal ( double dr, double dd, double *dl, double *db ) {
+
+  double v1[3];
+  double v2[3];
+
+/*
+*  L2,B2 system of galactic coordinates
+*
+*  P = 192.25       RA of galactic north pole (mean B1950.0)
+*  Q =  62.6        inclination of galactic to mean B1950.0 equator
+*  R =  33          longitude of ascending node
+*
+*  P,Q,R are degrees
+*
+*  Equatorial to galactic rotation matrix (J2000.0), obtained by
+*  applying the standard FK4 to FK5 transformation, for zero proper
+*  motion in FK5, to the columns of the B1950 equatorial to
+*  galactic rotation matrix:
+*/
+  double rmat[3][3] = {
+    { -0.054875539726,-0.873437108010,-0.483834985808 },
+    { +0.494109453312,-0.444829589425,+0.746982251810 },
+    { -0.867666135858,-0.198076386122,+0.455983795705 }
+  };
+
+  /* Spherical to Cartesian */
+  eraS2c( dr, dd, v1 );
+
+  /* Equatorial to Galactic */
+  eraRxp( rmat, v1, v2 );
+
+  /* Cartesian to spherical */
+  eraC2s( v2, dl, db );
+
+  /* Express in conventional ranges */
+  *dl = eraAnp( *dl );
+  *db = eraAnpm( *db );
+
+}
diff --git a/ast/pal/palEtrms.c b/ast/pal/palEtrms.c
new file mode 100644
index 0000000..4484682
--- /dev/null
+++ b/ast/pal/palEtrms.c
@@ -0,0 +1,106 @@
+/*
+*+
+*  Name:
+*     palEtrms
+
+*  Purpose:
+*     Compute the E-terms vector
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palEtrms ( double ep, double ev[3] );
+
+*  Arguments:
+*     ep = double (Given)
+*        Besselian epoch
+*     ev = double [3] (Returned)
+*        E-terms as (dx,dy,dz)
+
+*  Description:
+*     Computes the E-terms (elliptic component of annual aberration)
+*     vector.
+*
+*     Note the use of the J2000 aberration constant (20.49552 arcsec).
+*     This is a reflection of the fact that the E-terms embodied in
+*     existing star catalogues were computed from a variety of
+*     aberration constants.  Rather than adopting one of the old
+*     constants the latest value is used here.
+*
+*  See also:
+*     - Smith, C.A. et al., 1989.  Astr.J. 97, 265.
+*     - Yallop, B.D. et al., 1989.  Astr.J. 97, 274.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-12 (TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1996 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+
+void palEtrms ( double ep, double ev[3] ) {
+
+  /* Use the J2000 aberration constant */
+  const double ABCONST = 20.49552;
+
+  double t, e, e0, p, ek, cp;
+
+  /*  Julian centuries since B1950 */
+  t = (ep - 1950.) * .0100002135903;
+
+  /*  Eccentricity */
+  e = .01673011 - (t * 1.26e-7 + 4.193e-5) * t;
+
+  /*  Mean obliquity */
+  e0 = (84404.836 - ((t * .00181 + .00319) * t + 46.8495) * t) *
+    PAL__DAS2R;
+
+  /*  Mean longitude of perihelion */
+  p = (((t * .012 + 1.65) * t + 6190.67) * t + 1015489.951) *
+    PAL__DAS2R;
+
+  /*  E-terms */
+  ek = e * ABCONST * PAL__DAS2R;
+  cp = cos(p);
+  ev[0] = ek * sin(p);
+  ev[1] = -ek * cp * cos(e0);
+  ev[2] = -ek * cp * sin(e0);
+
+}
diff --git a/ast/pal/palEvp.c b/ast/pal/palEvp.c
new file mode 100644
index 0000000..b30a3a9
--- /dev/null
+++ b/ast/pal/palEvp.c
@@ -0,0 +1,110 @@
+/*
+*+
+*  Name:
+*     palEvp
+
+*  Purpose:
+*     Returns the barycentric and heliocentric velocity and position of the
+*     Earth.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palEvp( double date, double deqx, double dvb[3], double dpb[3],
+*                  double dvh[3], double dph[3] )
+
+*  Arguments:
+*     date = double (Given)
+*        TDB (loosely ET) as a Modified Julian Date (JD-2400000.5)
+*     deqx = double (Given)
+*        Julian epoch (e.g. 2000.0) of mean equator and equinox of the
+*        vectors returned.  If deqx <= 0.0, all vectors are referred to the
+*        mean equator and equinox (FK5) of epoch date.
+*     dvb = double[3] (Returned)
+*        Barycentric velocity (AU/s, AU)
+*     dpb = double[3] (Returned)
+*        Barycentric position (AU/s, AU)
+*     dvh = double[3] (Returned)
+*        heliocentric velocity (AU/s, AU)
+*     dph = double[3] (Returned)
+*        Heliocentric position (AU/s, AU)
+
+*  Description:
+*     Returns the barycentric and heliocentric velocity and position of the
+*     Earth at a given epoch, given with respect to a specified equinox.
+*     For information about accuracy, see the function eraEpv00.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-13 (PTW):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2005 Patrick T. Wallace
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+void palEvp( double date, double deqx, double dvb[3], double dpb[3],
+             double dvh[3], double dph[3] ){
+
+/* Local Variables; */
+   int i;
+   double pvh[2][3], pvb[2][3], d1, d2, r[3][3];
+
+/* BCRS PV-vectors. */
+   eraEpv00 ( 2400000.5, date, pvh, pvb );
+
+/* Was precession to another equinox requested? */
+   if ( deqx > 0.0 ) {
+
+/* Yes: compute precession matrix from J2000.0 to deqx. */
+      eraEpj2jd ( deqx, &d1, &d2 );
+      eraPmat06 ( d1, d2, r );
+
+/* Rotate the PV-vectors. */
+      eraRxpv ( r, pvh, pvh );
+      eraRxpv ( r, pvb, pvb );
+   }
+
+/* Return the required vectors. */
+   for ( i = 0; i < 3; i++ ) {
+      dvh[i] = pvh[1][i] / PAL__SPD;
+      dvb[i] = pvb[1][i] / PAL__SPD;
+      dph[i] = pvh[0][i];
+      dpb[i] = pvb[0][i];
+   }
+}
diff --git a/ast/pal/palFk45z.c b/ast/pal/palFk45z.c
new file mode 100644
index 0000000..643fd74
--- /dev/null
+++ b/ast/pal/palFk45z.c
@@ -0,0 +1,186 @@
+/*
+*+
+*  Name:
+*     palFk45z
+
+*  Purpose:
+*     Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero
+*     proper motion in the FK5 frame
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palFk45z( double r1950, double d1950, double bepoch, double *r2000,
+*               double *d2000 )
+
+*  Arguments:
+*     r1950 = double (Given)
+*        B1950.0 FK4 RA at epoch (radians).
+*     d1950 = double (Given)
+*        B1950.0 FK4 Dec at epoch (radians).
+*     bepoch = double (Given)
+*        Besselian epoch (e.g. 1979.3)
+*     r2000 = double (Returned)
+*        J2000.0 FK5 RA (Radians).
+*     d2000 = double (Returned)
+*        J2000.0 FK5 Dec(Radians).
+
+*  Description:
+*     Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero
+*     proper motion in the FK5 frame (double precision)
+*
+*     This function converts stars from the Bessel-Newcomb, FK4
+*     system to the IAU 1976, FK5, Fricke system, in such a
+*     way that the FK5 proper motion is zero.  Because such a star
+*     has, in general, a non-zero proper motion in the FK4 system,
+*     the routine requires the epoch at which the position in the
+*     FK4 system was determined.
+*
+*     The method is from Appendix 2 of Ref 1, but using the constants
+*     of Ref 4.
+
+*  Notes:
+*     - The epoch BEPOCH is strictly speaking Besselian, but if a
+*     Julian epoch is supplied the result will be affected only to
+*     a negligible extent.
+*
+*     - Conversion from Besselian epoch 1950.0 to Julian epoch 2000.0
+*     only is provided for.  Conversions involving other epochs will
+*     require use of the appropriate precession, proper motion, and
+*     E-terms routines before and/or after palFk45z is called.
+*
+*     - In the FK4 catalogue the proper motions of stars within 10
+*     degrees of the poles do not embody the differential E-term effect
+*     and should, strictly speaking, be handled in a different manner
+*     from stars outside these regions. However, given the general lack
+*     of homogeneity of the star data available for routine astrometry,
+*     the difficulties of handling positions that may have been
+*     determined from astrometric fields spanning the polar and non-polar
+*     regions, the likelihood that the differential E-terms effect was not
+*     taken into account when allowing for proper motion in past
+*     astrometry, and the undesirability of a discontinuity in the
+*     algorithm, the decision has been made in this routine to include the
+*     effect of differential E-terms on the proper motions for all stars,
+*     whether polar or not.  At epoch 2000, and measuring on the sky rather
+*     than in terms of dRA, the errors resulting from this simplification
+*     are less than 1 milliarcsecond in position and 1 milliarcsecond per
+*     century in proper motion.
+*
+*  References:
+*     - Aoki,S., et al, 1983.  Astron.Astrophys., 128, 263.
+*     - Smith, C.A. et al, 1989.  "The transformation of astrometric
+*       catalog systems to the equinox J2000.0".  Astron.J. 97, 265.
+*     - Yallop, B.D. et al, 1989.  "Transformation of mean star places
+*       from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
+*       Astron.J. 97, 274.
+*     - Seidelmann, P.K. (ed), 1992.  "Explanatory Supplement to
+*       the Astronomical Almanac", ISBN 0-935702-68-7.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-10 (DSB):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1998 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+void palFk45z( double r1950, double d1950, double bepoch, double *r2000,
+               double *d2000 ){
+
+/* Local Variables: */
+   double w;
+   int i;
+   int j;
+   double r0[3], a1[3], v1[3], v2[6]; /* Position and position+velocity vectors */
+
+
+/* CANONICAL CONSTANTS  (see references) */
+
+/* Vector A. */
+   double a[3] = { -1.62557E-6, -0.31919E-6, -0.13843E-6 };
+
+/* Vectors Adot. */
+   double ad[3] = { 1.245E-3, -1.580E-3, -0.659E-3 };
+
+/* Matrix M (only half of which is needed here). */
+   double em[6][3] = { {0.9999256782, -0.0111820611, -0.0048579477},
+                       {0.0111820610, 0.9999374784, -0.0000271765},
+                       {0.0048579479, -0.0000271474, 0.9999881997},
+                       {-0.000551, -0.238565, 0.435739},
+                       {0.238514, -0.002667, -0.008541},
+                       {-0.435623, 0.012254, 0.002117} };
+
+
+/* Spherical to Cartesian. */
+   eraS2c( r1950, d1950, r0 );
+
+/* Adjust vector A to give zero proper motion in FK5. */
+   w = ( bepoch - 1950.0 )/PAL__PMF;
+   for( i = 0; i < 3; i++ ) {
+      a1[ i ] = a[ i ] + w*ad[ i ];
+   }
+
+/* Remove e-terms. */
+   w = r0[ 0 ]*a1[ 0 ] + r0[ 1 ]*a1[ 1 ] + r0[ 2 ]*a1[ 2 ];
+   for( i = 0; i < 3; i++ ) {
+      v1[ i ] = r0[ i ] - a1[ i ] + w*r0[ i ];
+   }
+
+/* Convert position vector to Fricke system. */
+   for( i = 0; i < 6; i++ ) {
+      w = 0.0;
+      for( j = 0; j < 3; j++ ) {
+         w += em[ i ][ j ]*v1[ j ];
+      }
+      v2[ i ] = w;
+   }
+
+/* Allow for fictitious proper motion in FK4. */
+   w = ( palEpj( palEpb2d( bepoch ) ) - 2000.0 )/PAL__PMF;
+   for( i = 0; i < 3; i++ ) {
+      v2[ i ] += w*v2[ i + 3 ];
+   }
+
+/* Revert to spherical coordinates. */
+   eraC2s( v2, &w, d2000 );
+   *r2000 = eraAnp( w );
+}
+
+
diff --git a/ast/pal/palFk524.c b/ast/pal/palFk524.c
new file mode 100644
index 0000000..47c3abb
--- /dev/null
+++ b/ast/pal/palFk524.c
@@ -0,0 +1,259 @@
+/*
+*+
+*  Name:
+*     palFk524
+
+*  Purpose:
+*     Convert J2000.0 FK5 star data to B1950.0 FK4.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palFk524( double r2000, double d2000, double dr2000, double dd2000,
+*               double p2000, double v2000, double *r1950, double *d1950,
+*               double *dr1950, double *dd1950, double *p1950, double *v1950 )
+
+*  Arguments:
+*     r2000 = double (Given)
+*        J2000.0 FK5 RA (radians).
+*     d2000 = double (Given)
+*        J2000.0 FK5 Dec (radians).
+*     dr2000 = double (Given)
+*        J2000.0 FK5 RA proper motion (rad/Jul.yr)
+*     dd2000 = double (Given)
+*        J2000.0 FK5 Dec proper motion (rad/Jul.yr)
+*     p2000 = double (Given)
+*        J2000.0 FK5 parallax (arcsec)
+*     v2000 = double (Given)
+*         J2000.0 FK5 radial velocity (km/s, +ve = moving away)
+*     r1950 = double * (Returned)
+*        B1950.0 FK4 RA (radians).
+*     d1950 = double * (Returned)
+*        B1950.0 FK4 Dec (radians).
+*     dr1950 = double * (Returned)
+*        B1950.0 FK4 RA proper motion (rad/Jul.yr)
+*     dd1950 = double * (Returned)
+*        B1950.0 FK4 Dec proper motion (rad/Jul.yr)
+*     p1950 = double * (Returned)
+*        B1950.0 FK4 parallax (arcsec)
+*     v1950 = double * (Returned)
+*         B1950.0 FK4 radial velocity (km/s, +ve = moving away)
+
+*  Description:
+*     This function converts stars from the IAU 1976, FK5, Fricke
+*     system, to the Bessel-Newcomb, FK4 system.  The precepts
+*     of Smith et al (Ref 1) are followed, using the implementation
+*     by Yallop et al (Ref 2) of a matrix method due to Standish.
+*     Kinoshita's development of Andoyer's post-Newcomb precession is
+*     used.  The numerical constants from Seidelmann et al (Ref 3) are
+*     used canonically.
+
+*  Notes:
+*     - The proper motions in RA are dRA/dt rather than
+*     cos(Dec)*dRA/dt, and are per year rather than per century.
+*     - Note that conversion from Julian epoch 2000.0 to Besselian
+*     epoch 1950.0 only is provided for.  Conversions involving
+*     other epochs will require use of the appropriate precession,
+*     proper motion, and E-terms routines before and/or after
+*     FK524 is called.
+*     - In the FK4 catalogue the proper motions of stars within
+*     10 degrees of the poles do not embody the differential
+*     E-term effect and should, strictly speaking, be handled
+*     in a different manner from stars outside these regions.
+*     However, given the general lack of homogeneity of the star
+*     data available for routine astrometry, the difficulties of
+*     handling positions that may have been determined from
+*     astrometric fields spanning the polar and non-polar regions,
+*     the likelihood that the differential E-terms effect was not
+*     taken into account when allowing for proper motion in past
+*     astrometry, and the undesirability of a discontinuity in
+*     the algorithm, the decision has been made in this routine to
+*     include the effect of differential E-terms on the proper
+*     motions for all stars, whether polar or not.  At epoch 2000,
+*     and measuring on the sky rather than in terms of dRA, the
+*     errors resulting from this simplification are less than
+*     1 milliarcsecond in position and 1 milliarcsecond per
+*     century in proper motion.
+*
+*  References:
+*     - Smith, C.A. et al, 1989.  "The transformation of astrometric
+*       catalog systems to the equinox J2000.0".  Astron.J. 97, 265.
+*     - Yallop, B.D. et al, 1989.  "Transformation of mean star places
+*       from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
+*       Astron.J. 97, 274.
+*     - Seidelmann, P.K. (ed), 1992.  "Explanatory Supplement to
+*       the Astronomical Almanac", ISBN 0-935702-68-7.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-13 (DSB):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "math.h"
+
+void palFk524( double r2000, double d2000, double dr2000, double dd2000,
+               double p2000, double v2000, double *r1950, double *d1950,
+               double *dr1950, double *dd1950, double *p1950, double *v1950 ){
+
+/* Local Variables; */
+   double r, d, ur, ud, px, rv;
+   double sr, cr, sd, cd, x, y, z, w;
+   double v1[ 6 ], v2[ 6 ];
+   double xd, yd, zd;
+   double rxyz, wd, rxysq, rxy;
+   int i, j;
+
+/* Small number to avoid arithmetic problems. */
+   static const double tiny = 1.0E-30;
+
+/* Canonical constants (see references). Constant vector and matrix. */
+   double a[ 6 ] = { -1.62557E-6,  -0.31919E-6, -0.13843E-6,
+                     +1.245E-3,    -1.580E-3,   -0.659E-3 };
+   double emi[ 6 ][ 6 ] = {
+                 { 0.9999256795,      0.0111814828,      0.0048590039,
+                  -0.00000242389840, -0.00000002710544, -0.00000001177742},
+                 {-0.0111814828,      0.9999374849,     -0.0000271771,
+                   0.00000002710544, -0.00000242392702,  0.00000000006585 },
+                 {-0.0048590040,     -0.0000271557,      0.9999881946,
+                   0.00000001177742,  0.00000000006585, -0.00000242404995 },
+                 {-0.000551,          0.238509,         -0.435614,
+                   0.99990432,        0.01118145,        0.00485852 },
+                 {-0.238560,         -0.002667,          0.012254,
+                  -0.01118145,        0.99991613,       -0.00002717},
+                 { 0.435730,         -0.008541,          0.002117,
+                  -0.00485852,       -0.00002716,        0.99996684 } };
+
+/* Pick up J2000 data (units radians and arcsec/JC). */
+   r = r2000;
+   d = d2000;
+   ur = dr2000*PAL__PMF;
+   ud = dd2000*PAL__PMF;
+   px = p2000;
+   rv = v2000;
+
+/* Spherical to Cartesian. */
+   sr = sin( r );
+   cr = cos( r );
+   sd = sin( d );
+   cd = cos( d );
+
+   x = cr*cd;
+   y = sr*cd;
+   z =    sd;
+
+   w = PAL__VF*rv*px;
+
+   v1[ 0 ] = x;
+   v1[ 1 ] = y;
+   v1[ 2 ] = z;
+
+   v1[ 3 ] = -ur*y - cr*sd*ud + w*x;
+   v1[ 4 ] =  ur*x - sr*sd*ud + w*y;
+   v1[ 5 ] =            cd*ud + w*z;
+
+/* Convert position+velocity vector to BN system. */
+   for( i = 0; i < 6; i++ ) {
+      w = 0.0;
+      for( j = 0; j < 6; j++ ) {
+         w += emi[ i ][ j ]*v1[ j ];
+      }
+      v2[ i ] = w;
+   }
+
+/* Position vector components and magnitude. */
+   x = v2[ 0 ];
+   y = v2[ 1 ];
+   z = v2[ 2 ];
+   rxyz = sqrt( x*x + y*y + z*z );
+
+/* Apply E-terms to position. */
+   w = x*a[ 0 ] + y*a[ 1 ] + z*a[ 2 ];
+   x += a[ 0 ]*rxyz - w*x;
+   y += a[ 1 ]*rxyz - w*y;
+   z += a[ 2 ]*rxyz - w*z;
+
+/* Recompute magnitude. */
+   rxyz = sqrt( x*x + y*y + z*z );
+
+/* Apply E-terms to both position and velocity. */
+   x = v2[ 0 ];
+   y = v2[ 1 ];
+   z = v2[ 2 ];
+   w = x*a[ 0 ] + y*a[ 1 ] + z*a[ 2 ];
+   wd = x*a[ 3 ] + y*a[ 4 ] + z*a[ 5 ];
+   x += a[ 0 ]*rxyz - w*x;
+   y += a[ 1 ]*rxyz - w*y;
+   z += a[ 2 ]*rxyz - w*z;
+   xd = v2[ 3 ] + a[ 3 ]*rxyz - wd*x;
+   yd = v2[ 4 ] + a[ 4 ]*rxyz - wd*y;
+   zd = v2[ 5 ] + a[ 5 ]*rxyz - wd*z;
+
+/* Convert to spherical. */
+   rxysq = x*x + y*y;
+   rxy = sqrt( rxysq );
+
+   if( x == 0.0 && y == 0.0 ) {
+      r = 0.0;
+   } else {
+      r = atan2( y, x );
+      if( r <  0.0 ) r += PAL__D2PI;
+   }
+   d = atan2( z, rxy );
+
+   if( rxy > tiny ) {
+      ur = ( x*yd - y*xd )/rxysq;
+      ud = ( zd*rxysq - z*( x*xd + y*yd ) )/( ( rxysq + z*z )*rxy );
+   }
+
+/* Radial velocity and parallax. */
+   if( px > tiny ) {
+      rv = ( x*xd + y*yd + z*zd )/( px*PAL__VF*rxyz );
+      px /= rxyz;
+   }
+
+/* Return results. */
+   *r1950 = r;
+   *d1950 = d;
+   *dr1950 = ur/PAL__PMF;
+   *dd1950 = ud/PAL__PMF;
+   *p1950 = px;
+   *v1950 = rv;
+}
diff --git a/ast/pal/palFk54z.c b/ast/pal/palFk54z.c
new file mode 100644
index 0000000..b902b0d
--- /dev/null
+++ b/ast/pal/palFk54z.c
@@ -0,0 +1,113 @@
+/*
+*+
+*  Name:
+*     palFk54z
+
+*  Purpose:
+*     Convert a J2000.0 FK5 star position to B1950.0 FK4 assuming
+*     zero proper motion and parallax.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palFk54z( double r2000, double d2000, double bepoch, double *r1950,
+*               double *d1950, double *dr1950, double *dd1950 )
+
+*  Arguments:
+*     r2000 = double (Given)
+*        J2000.0 FK5 RA (radians).
+*     d2000 = double (Given)
+*        J2000.0 FK5 Dec (radians).
+*     bepoch = double (Given)
+*         Besselian epoch (e.g. 1950.0).
+*     r1950 = double * (Returned)
+*        B1950 FK4 RA (radians) at epoch "bepoch".
+*     d1950 = double * (Returned)
+*        B1950 FK4 Dec (radians) at epoch "bepoch".
+*     dr1950 = double * (Returned)
+*        B1950 FK4 proper motion (RA) (radians/trop.yr)).
+*     dr1950 = double * (Returned)
+*        B1950 FK4 proper motion (Dec) (radians/trop.yr)).
+
+*  Description:
+*     This function converts star positions from the IAU 1976,
+*     FK5, Fricke system to the Bessel-Newcomb, FK4 system.
+
+*  Notes:
+*     - The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.
+*     - Conversion from Julian epoch 2000.0 to Besselian epoch 1950.0
+*     only is provided for.  Conversions involving other epochs will
+*     require use of the appropriate precession functions before and
+*     after this function is called.
+*     - The FK5 proper motions, the parallax and the radial velocity
+*      are presumed zero.
+*     - It is the intention that FK5 should be a close approximation
+*     to an inertial frame, so that distant objects have zero proper
+*     motion;  such objects have (in general) non-zero proper motion
+*     in FK4, and this function returns those fictitious proper
+*     motions.
+*     - The position returned by this function is in the B1950
+*     reference frame but at Besselian epoch BEPOCH.  For comparison
+*     with catalogues the "bepoch" argument will frequently be 1950.0.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-13 (DSB):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palFk54z( double r2000, double d2000, double bepoch, double *r1950,
+               double *d1950, double *dr1950, double *dd1950 ){
+
+/* Local Variables: */
+   double r, d, px, rv, y;
+
+/* FK5 equinox J2000 (any epoch) to FK4 equinox B1950 epoch B1950. */
+   palFk524( r2000, d2000, 0.0, 0.0, 0.0, 0.0, &r, &d, dr1950, dd1950,
+             &px, &rv );
+
+/* Fictitious proper motion to epoch "bepoch". */
+   y = bepoch - 1950.0;
+   *r1950 = r + *dr1950*y;
+   *d1950 = d + *dd1950*y;
+
+}
+
diff --git a/ast/pal/palGaleq.c b/ast/pal/palGaleq.c
new file mode 100644
index 0000000..d0da1ee
--- /dev/null
+++ b/ast/pal/palGaleq.c
@@ -0,0 +1,118 @@
+/*
+*+
+*  Name:
+*     palGaleq
+
+*  Purpose:
+*     Convert from galactic to J2000.0 equatorial coordinates
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palGaleq ( double dl, double db, double *dr, double *dd );
+
+*  Arguments:
+*     dl = double (Given)
+*       Galactic longitude (radians).
+*     db = double (Given)
+*       Galactic latitude (radians).
+*     dr = double * (Returned)
+*       J2000.0 RA (radians)
+*     dd = double * (Returned)
+*       J2000.0 Dec (radians)
+
+*  Description:
+*     Transformation from IAU 1958 galactic coordinates to
+*     J2000.0 equatorial coordinates.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     The equatorial coordinates are J2000.0.  Use the routine
+*     palGe50 if conversion to B1950.0 'FK4' coordinates is
+*     required.
+
+*  See Also:
+*     Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960)
+
+*  History:
+*     2012-02-12(TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1998 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palGaleq ( double dl, double db, double *dr, double *dd ) {
+
+  double v1[3];
+  double v2[3];
+
+/*
+*  L2,B2 system of galactic coordinates
+*
+*  P = 192.25       RA of galactic north pole (mean B1950.0)
+*  Q =  62.6        inclination of galactic to mean B1950.0 equator
+*  R =  33          longitude of ascending node
+*
+*  P,Q,R are degrees
+*
+*  Equatorial to galactic rotation matrix (J2000.0), obtained by
+*  applying the standard FK4 to FK5 transformation, for zero proper
+*  motion in FK5, to the columns of the B1950 equatorial to
+*  galactic rotation matrix:
+*/
+  double rmat[3][3] = {
+    { -0.054875539726,-0.873437108010,-0.483834985808 },
+    { +0.494109453312,-0.444829589425,+0.746982251810 },
+    { -0.867666135858,-0.198076386122,+0.455983795705 }
+  };
+
+  /* Spherical to Cartesian */
+  eraS2c( dl, db, v1 );
+
+  /* Galactic to equatorial */
+  eraTrxp( rmat, v1, v2 );
+
+  /* Cartesian to spherical */
+  eraC2s( v2, dr, dd );
+
+  /* Express in conventional ranges */
+  *dr = eraAnp( *dr );
+  *dd = eraAnpm( *dd );
+
+}
diff --git a/ast/pal/palGalsup.c b/ast/pal/palGalsup.c
new file mode 100644
index 0000000..e469fb7
--- /dev/null
+++ b/ast/pal/palGalsup.c
@@ -0,0 +1,116 @@
+/*
+*+
+*  Name:
+*     palGalsup
+
+*  Purpose:
+*     Convert from galactic to supergalactic coordinates
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palGalsup ( double dl, double db, double *dsl, double *dsb );
+
+*  Arguments:
+*     dl = double (Given)
+*       Galactic longitude.
+*     db = double (Given)
+*       Galactic latitude.
+*     dsl = double * (Returned)
+*       Supergalactic longitude.
+*     dsb = double * (Returned)
+*       Supergalactic latitude.
+
+*  Description:
+*     Transformation from IAU 1958 galactic coordinates to
+*     de Vaucouleurs supergalactic coordinates.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  See Also:
+*     - de Vaucouleurs, de Vaucouleurs, & Corwin, Second Reference
+*       Catalogue of Bright Galaxies, U. Texas, page 8.
+*     - Systems & Applied Sciences Corp., Documentation for the
+*       machine-readable version of the above catalogue,
+*       Contract NAS 5-26490.
+*
+*     (These two references give different values for the galactic
+*     longitude of the supergalactic origin.  Both are wrong;  the
+*     correct value is L2=137.37.)
+
+*  History:
+*     2012-02-12(TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1999 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palGalsup ( double dl, double db, double *dsl, double *dsb ) {
+
+  double v1[3];
+  double v2[3];
+
+/*
+*  System of supergalactic coordinates:
+*
+*    SGL   SGB        L2     B2      (deg)
+*     -    +90      47.37  +6.32
+*     0     0         -      0
+*
+*  Galactic to supergalactic rotation matrix:
+*/
+  double rmat[3][3] = {
+    { -0.735742574804,+0.677261296414,+0.000000000000 },
+    { -0.074553778365,-0.080991471307,+0.993922590400 },
+    { +0.673145302109,+0.731271165817,+0.110081262225 }
+  };
+
+  /* Spherical to Cartesian */
+  eraS2c( dl, db, v1 );
+
+  /* Galactic to Supergalactic */
+  eraRxp( rmat, v1, v2 );
+
+  /* Cartesian to spherical */
+  eraC2s( v2, dsl, dsb );
+
+  /* Express in conventional ranges */
+  *dsl = eraAnp( *dsl );
+  *dsb = eraAnpm( *dsb );
+
+}
diff --git a/ast/pal/palGeoc.c b/ast/pal/palGeoc.c
new file mode 100644
index 0000000..9954d92
--- /dev/null
+++ b/ast/pal/palGeoc.c
@@ -0,0 +1,83 @@
+/*
+*+
+*  Name:
+*     palGeoc
+
+*  Purpose:
+*     Convert geodetic position to geocentric
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palGeoc( double p, double h, double * r, double *z );
+
+*  Arguments:
+*     p = double (Given)
+*       latitude (radians)
+*     h = double (Given)
+*       height above reference spheroid (geodetic, metres)
+*     r = double * (Returned)
+*       distance from Earth axis (AU)
+*     z = double * (Returned)
+*       distance from plane of Earth equator (AU)
+
+*  Description:
+*     Convert geodetic position to geocentric.
+
+*  Authors:
+*     PTW: Patrick T. Wallace
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - Geocentric latitude can be obtained by evaluating atan2(z,r)
+*     - Uses WGS84 reference ellipsoid and calls eraGd2gc
+
+*  History:
+*     2012-03-01 (TIMJ):
+*        Initial version moved from palOne2One
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2004 Patrick T. Wallace
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palGeoc ( double p, double h, double *r, double *z ) {
+  double xyz[3];
+  const double elong = 0.0;   /* Use zero longitude */
+  const double AU = 1.49597870E11;
+  /* WGS84 looks to be the closest match */
+  eraGd2gc( ERFA_WGS84, elong, p, h, xyz );
+  *r = xyz[0] / (AU * cos(elong) );
+  *z = xyz[2] / AU;
+}
diff --git a/ast/pal/palMappa.c b/ast/pal/palMappa.c
new file mode 100644
index 0000000..4e7ee64
--- /dev/null
+++ b/ast/pal/palMappa.c
@@ -0,0 +1,129 @@
+/*
+*+
+*  Name:
+*     palMappa
+
+*  Purpose:
+*     Compute parameters needed by palAmpqk and palMapqk.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palMappa( double eq, double date, double amprms[21] )
+
+*  Arguments:
+*     eq = double (Given)
+*        epoch of mean equinox to be used (Julian)
+*     date = double (Given)
+*        TDB (JD-2400000.5)
+*     amprms =   double[21]  (Returned)
+*        star-independent mean-to-apparent parameters:
+*        - (0)      time interval for proper motion (Julian years)
+*        - (1-3)    barycentric position of the Earth (AU)
+*        - (4-6)    heliocentric direction of the Earth (unit vector)
+*        - (7)      (Schwarzschild radius of Sun)/(Sun-Earth distance)
+*        - (8-10)   abv: barycentric Earth velocity in units of c
+*        - (11)     sqrt(1-v^2) where v=modulus(abv)
+*        - (12-20)  precession/nutation (3,3) matrix
+
+*  Description:
+*     Compute star-independent parameters in preparation for
+*     transformations between mean place and geocentric apparent place.
+*
+*     The parameters produced by this function are required in the
+*     parallax, aberration, and nutation/bias/precession parts of the
+*     mean/apparent transformations.
+*
+*     The reference systems and timescales used are IAU 2006.
+
+*  Notes:
+*     - For date, the distinction between the required TDB and TT
+*     is always negligible.  Moreover, for all but the most
+*     critical applications UTC is adequate.
+*     - The vector amprms(1-3) is referred to the mean equinox and
+*     equator of epoch eq.
+*     - The parameters amprms produced by this function are used by
+*     palAmpqk, palMapqk and palMapqkz.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-13 (PTW):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2003 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+#include <string.h>
+
+void palMappa( double eq, double date, double amprms[21] ){
+
+/* Local constants */
+
+/*  Gravitational radius of the Sun x 2 (2*mu/c**2, AU) */
+  const double GR2 = 2.0 * 9.87063e-9;
+
+/* Local Variables; */
+   int i;
+   double ebd[ 3 ], ehd[ 3 ], eh[ 3 ], e, vn[ 3 ], vm;
+
+/* Initialise so that unsused values are returned holding zero */
+   memset( amprms, 0, 21*sizeof( *amprms ) );
+
+/* Time interval for proper motion correction. */
+   amprms[ 0 ] = eraEpj( PAL__MJD0, date ) - eq;
+
+/* Get Earth barycentric and heliocentric position and velocity. */
+   palEvp( date, eq, ebd, &amprms[ 1 ], ehd, eh );
+
+/* Heliocentric direction of Earth (normalized) and modulus. */
+   eraPn( eh, &e, &amprms[ 4 ] );
+
+/* Light deflection parameter */
+   amprms[7] = GR2 / e;
+
+/* Aberration parameters. */
+   for( i = 0; i < 3; i++ ) {
+      amprms[ i + 8 ] = ebd[ i ]*PAL__CR;
+   }
+   eraPn( &amprms[8], &vm, vn );
+   amprms[ 11 ] = sqrt( 1.0 - vm*vm );
+
+/* NPB matrix. */
+   palPrenut( eq, date, (double(*)[ 3 ]) &amprms[ 12 ] );
+}
diff --git a/ast/pal/palMapqkz.c b/ast/pal/palMapqkz.c
new file mode 100644
index 0000000..44a2c79
--- /dev/null
+++ b/ast/pal/palMapqkz.c
@@ -0,0 +1,150 @@
+/*
+*+
+*  Name:
+*     palMapqkz
+
+*  Purpose:
+*     Quick mean to apparent place (no proper motion or parallax).
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palMapqkz( double rm, double dm, double amprms[21],
+*                     double *ra, double *da )
+
+*  Arguments:
+*     rm = double (Given)
+*        Mean RA (radians).
+*     dm = double (Given)
+*        Mean Dec (radians).
+*     amprms = double[21] (Given)
+*        Star-independent mean-to-apparent parameters (see palMappa):
+*        (0-3)    not used
+*        (4-6)    heliocentric direction of the Earth (unit vector)
+*        (7)      not used
+*        (8-10)   abv: barycentric Earth velocity in units of c
+*        (11)     sqrt(1-v^2) where v=modulus(abv)
+*        (12-20)  precession/nutation (3,3) matrix
+*     ra = double * (Returned)
+*        Apparent RA (radians).
+*     da = double * (Returned)
+*        Apparent Dec (radians).
+
+*  Description:
+*     Quick mean to apparent place:  transform a star RA,dec from
+*     mean place to geocentric apparent place, given the
+*     star-independent parameters, and assuming zero parallax
+*     and proper motion.
+*
+*     Use of this function is appropriate when efficiency is important
+*     and where many star positions, all with parallax and proper
+*     motion either zero or already allowed for, and all referred to
+*     the same equator and equinox, are to be transformed for one
+*     epoch.  The star-independent parameters can be obtained by
+*     calling the palMappa function.
+*
+*     The corresponding function for the case of non-zero parallax
+*     and proper motion is palMapqk.
+
+*  Notes:
+*     - The reference systems and timescales used are IAU 2006.
+*     - The mean place rm, dm and the vectors amprms[1-3] and amprms[4-6]
+*       are referred to the mean equinox and equator of the epoch
+*       specified when generating the precession/nutation matrix
+*       amprms[12-20].  In the call to palMappa (q.v.) normally used
+*       to populate amprms, this epoch is the first argument (eq).
+*     - The vector amprms(4-6) is referred to the mean equinox and
+*       equator of epoch eq.
+*     - Strictly speaking, the routine is not valid for solar-system
+*       sources, though the error will usually be extremely small.
+*       However, to prevent gross errors in the case where the
+*       position of the Sun is specified, the gravitational
+*       deflection term is restrained within about 920 arcsec of the
+*       centre of the Sun's disc.  The term has a maximum value of
+*       about 1.85 arcsec at this radius, and decreases to zero as
+*       the centre of the disc is approached.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-13 (PTW):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1999 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palMapqkz ( double rm, double dm, double amprms[21], double *ra,
+                 double *da ){
+
+/* Local Variables: */
+   int i;
+   double ab1, abv[3], p[3], w, p1dv, p2[3], p3[3];
+   double gr2e, pde, pdep1, ehn[3], p1[3];
+
+/* Unpack scalar and vector parameters. */
+   ab1 = amprms[11];
+   gr2e = amprms[7];
+   for( i = 0; i < 3; i++ ) {
+      abv[i] = amprms[i+8];
+      ehn[i] = amprms[i+4];
+   }
+
+/* Spherical to x,y,z. */
+   eraS2c( rm, dm, p );
+
+/* Light deflection (restrained within the Sun's disc) */
+   pde = eraPdp( p, ehn );
+   pdep1 = pde + 1.0;
+   w = gr2e / ( pdep1 > 1.0e-5 ? pdep1 : 1.0e-5 );
+   for( i = 0; i < 3; i++) {
+      p1[i] = p[i] + w * ( ehn[i] - pde * p[i] );
+   }
+
+/* Aberration. */
+   p1dv = eraPdp( p1, abv );
+   w = 1.0 + p1dv / ( ab1 + 1.0 );
+   for( i = 0; i < 3; i++ ) {
+      p2[i] = ( ( ab1 * p1[i] ) + ( w * abv[i] ) );
+   }
+
+/* Precession and nutation. */
+   eraRxp( (double(*)[3]) &amprms[12], p2, p3 );
+
+/* Geocentric apparent RA,dec. */
+   eraC2s( p3, ra, da );
+   *ra = eraAnp( *ra );
+}
diff --git a/ast/pal/palOne2One.c b/ast/pal/palOne2One.c
new file mode 100644
index 0000000..a8d8072
--- /dev/null
+++ b/ast/pal/palOne2One.c
@@ -0,0 +1,1482 @@
+/*
+*+
+*  Name:
+*     palOne2One
+
+*  Purpose:
+*     File containing simple PAL wrappers for SLA routines that are identical in SOFA
+
+*  Invocation:
+*     Matches SLA API
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Description:
+*     Some SOFA/ERFA routines are identical to their SLA counterparts. PAL provides
+*     direct counterparts to these although it is generally a better idea to
+*     use the SOFA/ERFA routine directly in new code.
+*
+*     The PAL routines with direct equivalents in SOFA/ERFA are:
+*     - palCldj
+*     - palDbear
+*     - palDaf2r
+*     - palDav2m
+*     - palDcc2s
+*     - palDcs2c
+*     - palDd2tf
+*     - palDimxv
+*     - palDm2av
+*     - palDjcl
+*     - palDmxm
+*     - palDmxv
+*     - palDpav
+*     - palDr2af
+*     - palDr2tf
+*     - palDranrm
+*     - palDsep
+*     - palDsepv
+*     - palDtf2d
+*     - palDtf2r
+*     - palDvdv
+*     - palDvn
+*     - palDvxv
+*     - palEpb
+*     - palEpb2d
+*     - palEpj
+*     - palEpj2d
+*     - palEqeqx
+*     - palFk5hz
+*     - palGmst
+*     - palGmsta
+*     - palHfk5z
+*     - palRefcoq
+
+*  Authors:
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     DSB: David S Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     - Do not call these functions from other PAL functions. Always use
+*       the SOFA/ERFA routines directly in new code.
+*     - These are implemented as real functions rather than C preprocessor
+*       macros so there may be a performance penalty in using the PAL
+*       version instead of the SOFA/ERFA version.
+*     - Routines that take MJDs have SOFA/ERFA equivalents that have an explicit
+*       MJD offset included.
+*     - palEqeqx, palGmst and palGmsta use the IAU 2006 precession model.
+
+*  History:
+*     2012-02-10 (TIMJ):
+*        Initial version
+*        Adapted with permission from the Fortran SLALIB library.
+*     2012-03-23 (TIMJ):
+*        Update prologue.
+*     2012-05-09 (DSBJ):
+*        Move palDrange into a separate file.
+*     2014-07-15 (TIMJ):
+*        SOFA now has palRefcoq equivalent.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2014 Tim Jenness
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+/*
+*+
+*  Name:
+*     palCldj
+
+*  Purpose:
+*     Gregorian Calendar to Modified Julian Date
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palCldj( int iy, int im, int id, double *djm, int *j );
+
+*  Arguments:
+*     iy = int (Given)
+*        Year in Gregorian calendar
+*     im = int (Given)
+*        Month in Gregorian calendar
+*     id = int (Given)
+*        Day in Gregorian calendar
+*     djm = double * (Returned)
+*        Modified Julian Date (JD-2400000.5) for 0 hrs
+*     j = int * (Returned)
+*        status: 0 = OK, 1 = bad year (MJD not computed),
+*        2 = bad month (MJD not computed), 3 = bad day (MJD computed).
+
+*  Description:
+*     Gregorian calendar to Modified Julian Date.
+
+*  Notes:
+*     - Uses eraCal2jd(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palCldj ( int iy, int im, int id, double *djm, int *j ) {
+  double djm0;
+  *j = eraCal2jd( iy, im, id, &djm0, djm );
+}
+
+/*
+*+
+*  Name:
+*     palDbear
+
+*  Purpose:
+*     Bearing (position angle) of one point on a sphere relative to another
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     pa = palDbear( double a1, double b1, double a2, double b2 );
+
+*  Arguments:
+*     a1 = double (Given)
+*        Longitude of point A (e.g. RA) in radians.
+*     a2 = double (Given)
+*        Latitude of point A (e.g. Dec) in radians.
+*     b1 = double (Given)
+*        Longitude of point B in radians.
+*     b2 = double (Given)
+*        Latitude of point B in radians.
+
+*  Returned Value:
+*     The result is the bearing (position angle), in radians, of point
+*     A2,B2 as seen from point A1,B1.  It is in the range +/- pi.  If
+*     A2,B2 is due east of A1,B1 the bearing is +pi/2.  Zero is returned
+*     if the two points are coincident.
+
+*  Description:
+*     Bearing (position angle) of one point in a sphere relative to another.
+
+*  Notes:
+*     - Uses eraPas(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+double palDbear ( double a1, double b1, double a2, double b2 ) {
+  return eraPas( a1, b1, a2, b2 );
+}
+
+/*
+*+
+*  Name:
+*     palDaf2r
+
+*  Purpose:
+*     Convert degrees, arcminutes, arcseconds to radians
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDaf2r( int ideg, int iamin, double asec, double *rad, int *j );
+
+*  Arguments:
+*     ideg = int (Given)
+*        Degrees.
+*     iamin = int (Given)
+*        Arcminutes.
+*     iasec = double (Given)
+*        Arcseconds.
+*     rad = double * (Returned)
+*        Angle in radians.
+*     j = int * (Returned)
+*        Status: 0 = OK, 1 = "ideg" out of range 0-359,
+*                2 = "iamin" outside of range 0-59,
+*                2 = "asec" outside range 0-59.99999
+
+*  Description:
+*     Convert degrees, arcminutes, arcseconds to radians.
+
+*  Notes:
+*     - Uses eraAf2a(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Arguments differ slightly. Assumes that the sign is always positive
+   and dealt with externally. */
+void palDaf2r ( int ideg, int iamin, double asec, double *rad, int *j ) {
+  *j = eraAf2a( ' ', ideg, iamin, asec, rad );
+}
+
+/*
+*+
+*  Name:
+*     palDav2m
+
+*  Purpose:
+*     Form the rotation matrix corresponding to a given axial vector.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDav2m( double axvec[3], double rmat[3][3] );
+
+*  Arguments:
+*     axvec = double [3] (Given)
+*       Axial vector (radians)
+*     rmat = double [3][3] (Returned)
+*       Rotation matrix.
+
+*  Description:
+*     A rotation matrix describes a rotation about some arbitrary axis,
+*     called the Euler axis.  The "axial vector" supplied to this routine
+*     has the same direction as the Euler axis, and its magnitude is the
+*     amount of rotation in radians.
+
+*  Notes:
+*     - Uses eraRv2m(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDav2m ( double axvec[3], double rmat[3][3] ) {
+  eraRv2m( axvec, rmat );
+}
+
+/*
+*+
+*  Name:
+*     palDcc2s
+
+*  Purpose:
+*     Cartesian to spherical coordinates
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDcc2s( double v[3], double *a, double *b );
+
+*  Arguments:
+*     v = double [3] (Given)
+*        x, y, z vector.
+*     a = double * (Returned)
+*        Spherical coordinate (radians)
+*     b = double * (Returned)
+*        Spherical coordinate (radians)
+
+*  Description:
+*     The spherical coordinates are longitude (+ve anticlockwise looking
+*     from the +ve latitude pole) and latitude.  The Cartesian coordinates
+*     are right handed, with the x axis at zero longitude and latitude, and
+*     the z axis at the +ve latitude pole.
+
+*  Notes:
+*     - Uses eraC2s(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDcc2s ( double v[3], double *a, double *b ) {
+  eraC2s( v, a, b );
+}
+
+/*
+*+
+*  Name:
+*     palDcs2c
+
+*  Purpose:
+*     Spherical coordinates to direction cosines
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDcs2c( double a, double b, double v[3] );
+
+*  Arguments:
+*     a = double (Given)
+*        Spherical coordinate in radians (ra, long etc).
+*     b = double (Given)
+*        Spherical coordinate in radians (dec, lat etc).
+*     v = double [3] (Returned)
+*        x, y, z vector
+
+*  Description:
+*     The spherical coordinates are longitude (+ve anticlockwise looking
+*     from the +ve latitude pole) and latitude.  The Cartesian coordinates
+*     are right handed, with the x axis at zero longitude and latitude, and
+*     the z axis at the +ve latitude pole.
+
+*  Notes:
+*     - Uses eraS2c(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDcs2c ( double a, double b, double v[3] ) {
+  eraS2c( a, b, v );
+}
+
+/*
+*+
+*  Name:
+*     palDd2tf
+
+*  Purpose:
+*     Convert an interval in days into hours, minutes, seconds
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDd2tf( int ndp, double days, char *sign, int ihmsf[4] );
+
+*  Arguments:
+*     ndp = int (Given)
+*        Number of decimal places of seconds
+*     days = double (Given)
+*        Interval in days
+*     sign = char * (Returned)
+*        '+' or '-' (single character, not string)
+*     ihmsf = int [4] (Returned)
+*        Hours, minutes, seconds, fraction
+
+*  Description:
+*     Convert and interval in days into hours, minutes, seconds.
+
+*  Notes:
+*     - Uses eraD2tf(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDd2tf ( int ndp, double days, char *sign, int ihmsf[4] ) {
+  eraD2tf( ndp, days, sign, ihmsf );
+}
+
+/*
+*+
+*  Name:
+*     palDimxv
+
+*  Purpose:
+*     Perform the 3-D backward unitary transformation
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDimxv( double dm[3][3], double va[3], double vb[3] );
+
+*  Arguments:
+*     dm = double [3][3] (Given)
+*        Matrix
+*     va = double [3] (Given)
+*        vector
+*     vb = double [3] (Returned)
+*        Result vector
+
+*  Description:
+*     Perform the 3-D backward unitary transformation.
+
+*  Notes:
+*     - Uses eraTrxp(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDimxv ( double dm[3][3], double va[3], double vb[3] ) {
+  eraTrxp( dm, va, vb );
+}
+
+/*
+*+
+*  Name:
+*     palDm2av
+
+*  Purpose:
+*     From a rotation matrix, determine the corresponding axial vector
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDm2av( double rmat[3][3], double axvec[3] );
+
+*  Arguments:
+*     rmat = double [3][3] (Given)
+*        Rotation matrix
+*     axvec = double [3] (Returned)
+*        Axial vector (radians)
+
+*  Description:
+*     A rotation matrix describes a rotation about some arbitrary axis,
+*     called the Euler axis.  The "axial vector" returned by this routine
+*     has the same direction as the Euler axis, and its magnitude is the
+*     amount of rotation in radians.  (The magnitude and direction can be
+*     separated by means of the routine palDvn.)
+
+*  Notes:
+*     - Uses eraRm2v(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDm2av ( double rmat[3][3], double axvec[3] ) {
+  eraRm2v( rmat, axvec );
+}
+
+/*
+*+
+*  Name:
+*     palDjcl
+
+*  Purpose:
+*     Modified Julian Date to Gregorian year, month, day and fraction of day
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDjcl( double djm, int *iy, int *im, int *id, double *fd, int *j );
+
+*  Arguments:
+*     djm = double (Given)
+*        modified Julian Date (JD-2400000.5)
+*     iy = int * (Returned)
+*        year
+*     im = int * (Returned)
+*        month
+*     id = int * (Returned)
+*        day
+*     fd = double * (Returned)
+*        Fraction of day.
+
+*  Description:
+*     Modified Julian Date to Gregorian year, month, day and fraction of day.
+
+*  Notes:
+*     - Uses eraJd2cal(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Requires additional SLA MJD reference date */
+void palDjcl ( double djm, int *iy, int *im, int *id, double *fd, int *j ) {
+  *j = eraJd2cal( PAL__MJD0, djm, iy, im, id, fd );
+}
+
+/*
+*+
+*  Name:
+*     palDmxm
+
+*  Purpose:
+*     Product of two 3x3 matrices
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDmxm( double a[3][3], double b[3][3], double c[3][3] );
+
+*  Arguments:
+*     a = double [3][3] (Given)
+*        Matrix
+*     b = double [3][3] (Given)
+*        Matrix
+*     c = double [3][3] (Returned)
+*        Matrix result
+
+*  Description:
+*     Product of two 3x3 matrices.
+
+*  Notes:
+*     - Uses eraRxr(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDmxm ( double a[3][3], double b[3][3], double c[3][3] ) {
+  eraRxr( a, b, c );
+}
+
+/*
+*+
+*  Name:
+*     palDmxv
+
+*  Purpose:
+*     Performs the 3-D forward unitary transformation
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDmxv( double dm[3][3], double va[3], double vb[3] );
+
+*  Arguments:
+*     dm = double [3][3] (Given)
+*        matrix
+*     va = double [3] (Given)
+*        vector
+*     dp = double [3] (Returned)
+*        result vector
+
+*  Description:
+*     Performs the 3-D forward unitary transformation.
+
+*  Notes:
+*     - Uses eraRxp(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDmxv ( double dm[3][3], double va[3], double vb[3] ) {
+  eraRxp( dm, va, vb );
+}
+
+/*
+*+
+*  Name:
+*     palDpav
+
+*  Purpose:
+*     Position angle of one celestial direction with respect to another
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     pa = palDpav( double v1[3], double v2[3] );
+
+*  Arguments:
+*     v1 = double [3] (Given)
+*        direction cosines of one point.
+*     v2 = double [3] (Given)
+*        direction cosines of the other point.
+
+*  Returned Value:
+*     The result is the bearing (position angle), in radians, of point
+*     V2 with respect to point V1.  It is in the range +/- pi.  The
+*     sense is such that if V2 is a small distance east of V1, the
+*     bearing is about +pi/2.  Zero is returned if the two points
+*     are coincident.
+
+*  Description:
+*     Position angle of one celestial direction with respect to another.
+
+*  Notes:
+*     - The coordinate frames correspond to RA,Dec, Long,Lat etc.
+*     - Uses eraPap(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+double palDpav ( double v1[3], double v2[3] ) {
+  return eraPap( v1, v2 );
+}
+
+/*
+*+
+*  Name:
+*     palDr2af
+
+*  Purpose:
+*     Convert an angle in radians to degrees, arcminutes, arcseconds
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDr2af( int ndp, double angle, char *sign, int idmsf[4] );
+
+*  Arguments:
+*     ndp = int (Given)
+*        number of decimal places of arcseconds
+*     angle = double (Given)
+*        angle in radians
+*     sign = char * (Returned)
+*        '+' or '-' (single character)
+*     idmsf = int [4] (Returned)
+*        Degrees, arcminutes, arcseconds, fraction
+
+*  Description:
+*     Convert an angle in radians to degrees, arcminutes, arcseconds.
+
+*  Notes:
+*     - Uses eraA2af(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDr2af ( int ndp, double angle, char *sign, int idmsf[4] ) {
+  eraA2af( ndp, angle, sign, idmsf );
+}
+
+/*
+*+
+*  Name:
+*     palDr2tf
+
+*  Purpose:
+*     Convert an angle in radians to hours, minutes, seconds
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDr2tf ( int ndp, double angle, char *sign, int ihmsf[4] );
+
+*  Arguments:
+*     ndp = int (Given)
+*        number of decimal places of arcseconds
+*     angle = double (Given)
+*        angle in radians
+*     sign = char * (Returned)
+*        '+' or '-' (single character)
+*     idmsf = int [4] (Returned)
+*        Hours, minutes, seconds, fraction
+
+*  Description:
+*     Convert an angle in radians to hours, minutes, seconds.
+
+*  Notes:
+*     - Uses eraA2tf(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDr2tf( int ndp, double angle, char *sign, int ihmsf[4] ) {
+  eraA2tf( ndp, angle, sign, ihmsf );
+}
+
+/*
+*+
+*  Name:
+*     palDranrm
+
+*  Purpose:
+*     Normalize angle into range 0-2 pi
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     norm = palDranrm( double angle );
+
+*  Arguments:
+*     angle = double (Given)
+*        angle in radians
+
+*  Returned Value:
+*     Angle expressed in the range 0-2 pi
+
+*  Description:
+*     Normalize angle into range 0-2 pi.
+
+*  Notes:
+*     - Uses eraAnp(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+double palDranrm ( double angle ) {
+  return eraAnp( angle );
+}
+
+/*
+*+
+*  Name:
+*     palDsep
+
+*  Purpose:
+*     Angle between two points on a sphere
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     ang = palDsep( double a1, double b1, double a2, double b2 );
+
+*  Arguments:
+*     a1 = double (Given)
+*        Spherical coordinate of one point (radians)
+*     b1 = double (Given)
+*        Spherical coordinate of one point (radians)
+*     a2 = double (Given)
+*        Spherical coordinate of other point (radians)
+*     b2 = double (Given)
+*        Spherical coordinate of other point (radians)
+
+*  Returned Value:
+*     Angle, in radians, between the two points. Always positive.
+
+*  Description:
+*     Angle between two points on a sphere.
+
+*  Notes:
+*     - The spherical coordinates are [RA,Dec], [Long,Lat] etc, in radians.
+*     - Uses eraSeps(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+double palDsep ( double a1, double b1, double a2, double b2 ) {
+  return eraSeps( a1, b1, a2, b2 );
+}
+
+/*
+*+
+*  Name:
+*     palDsepv
+
+*  Purpose:
+*     Angle between two vectors
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     ang = palDsepv( double v1[3], double v2[3] );
+
+*  Arguments:
+*     v1 = double [3] (Given)
+*        First vector
+*     v2 = double [3] (Given)
+*        Second vector
+
+*  Returned Value:
+*     Angle, in radians, between the two points. Always positive.
+
+*  Description:
+*     Angle between two vectors.
+
+*  Notes:
+*     - Uses eraSepp(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+double palDsepv ( double v1[3], double v2[3] ) {
+  return eraSepp( v1, v2 );
+}
+
+/*
+*+
+*  Name:
+*     palDtf2d
+
+*  Purpose:
+*     Convert hours, minutes, seconds to days
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDtf2d( int ihour, int imin, double sec, double *days, int *j );
+
+*  Arguments:
+*     ihour = int (Given)
+*        Hours
+*     imin = int (Given)
+*        Minutes
+*     sec = double (Given)
+*        Seconds
+*     days = double * (Returned)
+*        Interval in days
+*     j = int * (Returned)
+*        status: 0 = ok, 1 = ihour outside range 0-23,
+*        2 = imin outside range 0-59, 3 = sec outside range 0-59.999...
+
+*  Description:
+*     Convert hours, minutes, seconds to days.
+
+*  Notes:
+*     - Uses eraTf2d(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Assumes that the sign is always positive and is dealt with externally */
+void palDtf2d ( int ihour, int imin, double sec, double *days, int *j ) {
+  *j = eraTf2d( ' ', ihour, imin, sec, days );
+}
+
+/*
+*+
+*  Name:
+*     palDtf2r
+
+*  Purpose:
+*     Convert hours, minutes, seconds to radians
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDtf2r( int ihour, int imin, double sec, double *rad, int *j );
+
+*  Arguments:
+*     ihour = int (Given)
+*        Hours
+*     imin = int (Given)
+*        Minutes
+*     sec = double (Given)
+*        Seconds
+*     days = double * (Returned)
+*        Angle in radians
+*     j = int * (Returned)
+*        status: 0 = ok, 1 = ihour outside range 0-23,
+*        2 = imin outside range 0-59, 3 = sec outside range 0-59.999...
+
+*  Description:
+*     Convert hours, minutes, seconds to radians.
+
+*  Notes:
+*     - Uses eraTf2a(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Assumes that the sign is dealt with outside this routine */
+void palDtf2r ( int ihour, int imin, double sec, double *rad, int *j ) {
+  *j = eraTf2a( ' ', ihour, imin, sec, rad );
+}
+
+/*
+*+
+*  Name:
+*     palDvdv
+
+*  Purpose:
+*     Scalar product of two 3-vectors
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     prod = palDvdv ( double va[3], double vb[3] );
+
+*  Arguments:
+*     va = double [3] (Given)
+*        First vector
+*     vb = double [3] (Given)
+*        Second vector
+
+*  Returned Value:
+*     Scalar product va.vb
+
+*  Description:
+*     Scalar product of two 3-vectors.
+
+*  Notes:
+*     - Uses eraPdp(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+double palDvdv ( double va[3], double vb[3] ) {
+  return eraPdp( va, vb );
+}
+
+/*
+*+
+*  Name:
+*     palDvn
+
+*  Purpose:
+*     Normalizes a 3-vector also giving the modulus
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDvn( double v[3], double uv[3], double *vm );
+
+*  Arguments:
+*     v = double [3] (Given)
+*        vector
+*     uv = double [3] (Returned)
+*        unit vector in direction of "v"
+*     vm = double * (Returned)
+*        modulus of "v"
+
+*  Description:
+*     Normalizes a 3-vector also giving the modulus.
+
+*  Notes:
+*     - Uses eraPn(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Note that the arguments are flipped */
+void palDvn ( double v[3], double uv[3], double *vm ) {
+  eraPn( v, vm, uv );
+}
+
+/*
+*+
+*  Name:
+*     palDvxv
+
+*  Purpose:
+*     Vector product of two 3-vectors
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palDvxv( double va[3], double vb[3], double vc[3] );
+
+*  Arguments:
+*     va = double [3] (Given)
+*        First vector
+*     vb = double [3] (Given)
+*        Second vector
+*     vc = double [3] (Returned)
+*        Result vector
+
+*  Description:
+*     Vector product of two 3-vectors.
+
+*  Notes:
+*     - Uses eraPxp(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palDvxv ( double va[3], double vb[3], double vc[3] ) {
+  eraPxp( va, vb, vc );
+}
+
+/*
+*+
+*  Name:
+*     palEpb
+
+*  Purpose:
+*     Conversion of modified Julian Data to Besselian Epoch
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     epb = palEpb ( double date );
+
+*  Arguments:
+*     date = double (Given)
+*        Modified Julian Date (JD - 2400000.5)
+
+*  Returned Value:
+*      Besselian epoch.
+
+*  Description:
+*     Conversion of modified Julian Data to Besselian Epoch.
+
+*  Notes:
+*     - Uses eraEpb(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Requires additional SLA MJD reference date */
+double palEpb ( double date ) {
+  return eraEpb( PAL__MJD0, date );
+}
+
+/*
+*+
+*  Name:
+*     palEpb2d
+
+*  Purpose:
+*     Conversion of Besselian Epoch to Modified Julian Date
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     mjd = palEpb2d ( double epb );
+
+*  Arguments:
+*     epb = double (Given)
+*        Besselian Epoch
+
+*  Returned Value:
+*     Modified Julian Date (JD - 2400000.5)
+
+*  Description:
+*     Conversion of Besselian Epoch to Modified Julian Date.
+
+*  Notes:
+*     - Uses eraEpb2jd(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+
+double palEpb2d ( double epb ) {
+  double djm0, djm;
+  eraEpb2jd( epb, &djm0, &djm );
+  return djm;
+}
+
+/*
+*+
+*  Name:
+*     palEpj
+
+*  Purpose:
+*     Conversion of Modified Julian Date to Julian Epoch
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     epj = palEpj ( double date );
+
+*  Arguments:
+*     date = double (Given)
+*        Modified Julian Date (JD - 2400000.5)
+
+*  Returned Value:
+*     The Julian Epoch.
+
+*  Description:
+*     Conversion of Modified Julian Date to Julian Epoch.
+
+*  Notes:
+*     - Uses eraEpj(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Requires additional SLA MJD reference date */
+double palEpj ( double date ) {
+  return eraEpj( PAL__MJD0, date );
+}
+
+/*
+*+
+*  Name:
+*     palEpj2d
+
+*  Purpose:
+*     Conversion of Julian Epoch to Modified Julian Date
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     mjd = palEpj2d ( double epj );
+
+*  Arguments:
+*     epj = double (Given)
+*        Julian Epoch.
+
+*  Returned Value:
+*     Modified Julian Date (JD - 2400000.5)
+
+*  Description:
+*     Conversion of Julian Epoch to Modified Julian Date.
+
+*  Notes:
+*     - Uses eraEpj2d(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+double palEpj2d ( double epj ) {
+  double djm0, djm;
+  eraEpj2jd( epj, &djm0, &djm );
+  return djm;
+}
+
+/*
+*+
+*  Name:
+*     palEqeqx
+
+*  Purpose:
+*     Equation of the equinoxes (IAU 2000/2006)
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palEqeqx( double date );
+
+*  Arguments:
+*     date = double (Given)
+*        TT as Modified Julian Date (JD-400000.5)
+
+*  Description:
+*     Equation of the equinoxes (IAU 2000/2006).
+
+*  Notes:
+*     - Uses eraEe06a(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Requires additional SLA MJD reference date */
+double palEqeqx ( double date ) {
+  return eraEe06a( PAL__MJD0, date );
+}
+
+/*
+*+
+*  Name:
+*     palFk5hz
+
+*  Purpose:
+*     Transform an FK5 (J2000) star position into the frame of the
+*     Hipparcos catalogue.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palFk5hz ( double r5, double d5, double epoch,
+*                double *rh, double *dh );
+
+*  Arguments:
+*     r5 = double (Given)
+*        FK5 RA (radians), equinox J2000, epoch "epoch"
+*     d5 = double (Given)
+*        FK5 dec (radians), equinox J2000, epoch "epoch"
+*     epoch = double (Given)
+*        Julian epoch
+*     rh = double * (Returned)
+*        RA (radians)
+*     dh = double * (Returned)
+*        Dec (radians)
+
+*  Description:
+*     Transform an FK5 (J2000) star position into the frame of the
+*     Hipparcos catalogue.
+
+*  Notes:
+*     - Assumes zero Hipparcos proper motion.
+*     - Uses eraEpj2jd() and eraFk5hz.
+*       See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palFk5hz ( double r5, double d5, double epoch,
+                double *rh, double *dh ) {
+  /* Need to convert epoch to Julian date first */
+  double date1, date2;
+  eraEpj2jd( epoch, &date1, &date2 );
+  eraFk5hz( r5, d5, date1, date2, rh, dh );
+}
+
+/*
+*+
+*  Name:
+*     palGmst
+
+*  Purpose:
+*     Greenwich mean sidereal time (consistent with IAU 2006 precession).
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     mst = palGmst ( double ut1 );
+
+*  Arguments:
+*     ut1 = double (Given)
+*        Universal time (UT1) expressed as modified Julian Date (JD-2400000.5)
+
+*  Returned Value:
+*     Greenwich mean sidereal time
+
+*  Description:
+*     Greenwich mean sidereal time (consistent with IAU 2006 precession).
+
+*  Notes:
+*     - Uses eraGmst06(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Note that SOFA/ERFA has more accurate time arguments
+   and we use the 2006 precession model */
+double palGmst ( double ut1 ) {
+  return eraGmst06( PAL__MJD0, ut1, PAL__MJD0, ut1 );
+}
+
+/*
+*+
+*  Name:
+*     palGmsta
+
+*  Purpose:
+*     Greenwich mean sidereal time (consistent with IAU 2006 precession).
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     mst = palGmsta ( double date, double ut1 );
+
+*  Arguments:
+*     date = double (Given)
+*        UT1 date (MJD: integer part of JD-2400000.5)
+*     ut1 = double (Given)
+*        UT1 time (fraction of a day)
+
+*  Returned Value:
+*     Greenwich mean sidereal time (in range 0 to 2 pi)
+
+*  Description:
+*     Greenwich mean sidereal time (consistent with IAU 2006 precession).
+
+*  Notes:
+*     - For best accuracy use eraGmst06() directly.
+*     - Uses eraGmst06(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+/* Slightly better but still not as accurate as SOFA/ERFA */
+
+double palGmsta( double date, double ut ) {
+  date += PAL__MJD0;
+  return eraGmst06( date, ut, date, ut );
+}
+
+/*
+*+
+*  Name:
+*     palHfk5z
+
+*  Purpose:
+*     Hipparcos star position to FK5 J2000
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palHfk5z( double rh, double dh, double epoch,
+*               double *r5, double *d5, double *dr5, double *dd5 );
+
+*  Arguments:
+*     rh = double (Given)
+*        Hipparcos RA (radians)
+*     dh = double (Given)
+*        Hipparcos Dec (radians)
+*     epoch = double (Given)
+*        Julian epoch (TDB)
+*     r5 = double * (Returned)
+*        RA (radians, FK5, equinox J2000, epoch "epoch")
+*     d5 = double * (Returned)
+*        Dec (radians, FK5, equinox J2000, epoch "epoch")
+
+*  Description:
+*     Transform a Hipparcos star position into FK5 J2000, assuming
+*     zero Hipparcos proper motion.
+
+*  Notes:
+*     - Uses eraEpj2jd and eraHfk5z(). See SOFA/ERFA documentation for details.
+
+*-
+*/
+
+void palHfk5z ( double rh, double dh, double epoch,
+                double *r5, double *d5, double *dr5, double *dd5 ) {
+  /* Need to convert epoch to Julian date first */
+  double date1, date2;
+  eraEpj2jd( epoch, &date1, &date2 );
+  eraHfk5z( rh, dh, date1, date2, r5, d5, dr5, dd5 );
+}
+
+/*
+*+
+*  Name:
+*     palRefcoq
+
+*  Purpose:
+*     Determine the constants A and B in the atmospheric refraction model
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palRefcoq( double tdk, double pmb, double rh, double wl,
+*                double *refa, double *refb );
+
+*  Arguments:
+*     tdk = double (Given)
+*        Ambient temperature at the observer (K)
+*     pmb = double (Given)
+*        Pressure at the observer (millibar)
+*     rh =  double (Given)
+*        Relative humidity at the observer (range 0-1)
+*     wl =  double (Given)
+*        Effective wavelength of the source (micrometre).
+*        Radio refraction is chosen by specifying wl > 100 micrometres.
+*     refa = double * (Returned)
+*        tan Z coefficient (radian)
+*     refb = double * (Returned)
+*        tan**3 Z coefficient (radian)
+
+*  Description:
+*     Determine the constants A and B in the atmospheric refraction
+*     model dZ = A tan Z + B tan**3 Z.  This is a fast alternative
+*     to the palRefco routine.
+*
+*     Z is the "observed" zenith distance (i.e. affected by refraction)
+*     and dZ is what to add to Z to give the "topocentric" (i.e. in vacuo)
+*     zenith distance.
+
+*  Notes:
+*     - Uses eraRefco(). See SOFA/ERFA documentation for details.
+*     - Note that the SOFA/ERFA routine uses different order of
+*       of arguments and uses deg C rather than K.
+
+*-
+*/
+
+void palRefcoq ( double tdk, double pmb, double rh, double wl,
+                 double *refa, double *refb ) {
+  /* Note that SLA (and therefore PAL) uses units of kelvin
+     but SOFA/ERFA uses deg C */
+  eraRefco( pmb, tdk - 273.15, rh, wl, refa, refb );
+}
diff --git a/ast/pal/palPrebn.c b/ast/pal/palPrebn.c
new file mode 100644
index 0000000..989ce59
--- /dev/null
+++ b/ast/pal/palPrebn.c
@@ -0,0 +1,98 @@
+/*
+*+
+*  Name:
+*     palPrebn
+
+*  Purpose:
+*     Generate the matrix of precession between two objects (old)
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palPrebn ( double bep0, double bep1, double rmatp[3][3] );
+
+*  Arguments:
+*     bep0 = double (Given)
+*        Beginning Besselian epoch.
+*     bep1 = double (Given)
+*        Ending Besselian epoch
+*     rmatp = double[3][3] (Returned)
+*        precession matrix in the sense V(BEP1) = RMATP * V(BEP0)
+
+*  Description:
+*     Generate the matrix of precession between two epochs,
+*     using the old, pre-IAU1976, Bessel-Newcomb model, using
+*     Kinoshita's formulation
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  See Also:
+*     Kinoshita, H. (1975) 'Formulas for precession', SAO Special
+*     Report No. 364, Smithsonian Institution Astrophysical
+*     Observatory, Cambridge, Massachusetts.
+
+*  History:
+*     2012-02-12(TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1996 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+
+void palPrebn ( double bep0, double bep1, double rmatp[3][3] ) {
+
+  double t,bigt, zeta, theta, z, tas2r, w;
+
+  /* Interval between basic epoch B1850.0 and beginning epoch in TC */
+  bigt = (bep0-1850)/100.;
+
+  /*  Interval over which precession required, in tropical centuries */
+  t = (bep1-bep0)/100.;
+
+  /* Euler angles */
+  tas2r = t * PAL__DAS2R;
+  w = 2303.5548 + ( 1.39720 + 0.000059 * bigt) * bigt;
+
+  zeta = ( w + ( 0.30242 - 0.000269 * bigt + 0.017996 * t ) * t ) * tas2r;
+  z = ( w + ( 1.09478 + 0.000387 * bigt + 0.018324 * t ) * t ) * tas2r;
+  theta = ( 2005.1125 + ( -0.85294 - 0.000365 * bigt ) * bigt +
+	    (-0.42647 - 0.000365 * bigt - 0.041802 * t ) * t ) * tas2r;
+
+  /*  Rotation matrix */
+  palDeuler("ZYZ", -zeta, theta, -z, rmatp);
+
+}
diff --git a/ast/pal/palPrec.c b/ast/pal/palPrec.c
new file mode 100644
index 0000000..678770d
--- /dev/null
+++ b/ast/pal/palPrec.c
@@ -0,0 +1,107 @@
+/*
+*+
+*  Name:
+*     palPrec
+
+*  Purpose:
+*     Form the matrix of precession between two epochs (IAU 2006)
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palPrec( double ep0, double ep1, double rmatp[3][3] )
+
+*  Arguments:
+*     ep0 = double (Given)
+*        Beginning epoch
+*     ep1 = double (Given)
+*        Ending epoch
+*     rmatp = double[3][3] (Returned)
+*        Precession matrix
+
+*  Description:
+*     The IAU 2006 precession matrix from ep0 to ep1 is found and
+*     returned. The matrix is in the sense  V(EP1)  =  RMATP * V(EP0).
+*     The epochs are TDB (loosely TT) Julian epochs.
+*
+*     Though the matrix method itself is rigorous, the precession
+*     angles are expressed through canonical polynomials which are
+*     valid only for a limited time span of a few hundred years around
+*     the current epoch.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-10 (DSB):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1996 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palPrec( double ep0, double ep1, double rmatp[3][3] ){
+
+/* Local Variables: */
+   double rmatq[3][3];
+   double ep0_days;
+   double ep1_days;
+
+/* Convert supplied dates to days since J2000 */
+   ep0_days = ( ep0 - 2000.0 )*ERFA_DJY;
+   ep1_days = ( ep1 - 2000.0 )*ERFA_DJY;
+
+/* If beginning epoch is J2000, just return the rotation matrix from
+   J2000 to EP1. */
+   if( ep0 == 2000.0 ) {
+      eraPmat06( ERFA_DJ00, ep1_days, rmatp );
+
+/* If end epoch is J2000, get the rotation matrix from J2000 to EP0 and
+   then transpose it to get the rotation matrix from EP0 to J2000. */
+   } else if( ep1 == 2000.0 ) {
+      eraPmat06( ERFA_DJ00, ep0_days, rmatp );
+      eraTr( rmatp, rmatp );
+
+/* Otherwise. get the two matrices used above and multiply them
+   together. */
+   } else {
+      eraPmat06( ERFA_DJ00, ep0_days, rmatp );
+      eraTr( rmatp, rmatp );
+      eraPmat06( ERFA_DJ00, ep1_days, rmatq );
+      eraRxr( rmatp, rmatq, rmatp );
+   }
+
+}
diff --git a/ast/pal/palPrenut.c b/ast/pal/palPrenut.c
new file mode 100644
index 0000000..04e36f3
--- /dev/null
+++ b/ast/pal/palPrenut.c
@@ -0,0 +1,111 @@
+/*
+*+
+*  Name:
+*     palPrenut
+
+*  Purpose:
+*     Form the matrix of bias-precession-nutation (IAU 2006/2000A)
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palPrenut( double epoch, double date, double rmatpn[3][3] )
+
+*  Arguments:
+*     epoch = double (Returned)
+*        Julian epoch for mean coordinates.
+*     date = double (Returned)
+*        Modified Julian Date (JD-2400000.5) for true coordinates.
+*     rmatpn = double[3][3] (Returned)
+*        combined NPB matrix
+
+*  Description:
+*     Form the matrix of bias-precession-nutation (IAU 2006/2000A).
+*     The epoch and date are TT (but TDB is usually close enough).
+*     The matrix is in the sense   v(true)  =  rmatpn * v(mean).
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-10 (PTW):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+void palPrenut ( double epoch, double date, double rmatpn[3][3] ){
+
+/* Local Variables: */
+   double bpa;
+   double bpia;
+   double bqa;
+   double chia;
+   double d1;
+   double d2;
+   double eps0;
+   double epsa;
+   double gam;
+   double oma;
+   double pa;
+   double phi;
+   double pia;
+   double psi;
+   double psia;
+   double r1[3][3];
+   double r2[3][3];
+   double thetaa;
+   double za;
+   double zetaa;
+
+/* Specified Julian epoch as a 2-part JD. */
+   eraEpj2jd( epoch, &d1, &d2 );
+
+/* P matrix, from specified epoch to J2000.0. */
+   eraP06e( d1, d2, &eps0, &psia, &oma, &bpa, &bqa, &pia, &bpia, &epsa,
+            &chia, &za, &zetaa, &thetaa, &pa, &gam, &phi, &psi );
+   eraIr( r1 );
+   eraRz( -chia, r1 );
+   eraRx( oma, r1 );
+   eraRz( psia, r1 );
+   eraRx( -eps0, r1 );
+
+/* NPB matrix, from J2000.0 to date. */
+   eraPnm06a( PAL__MJD0, date, r2 );
+
+/* NPB matrix, from specified epoch to date. */
+   eraRxr( r2, r1, rmatpn );
+}
diff --git a/ast/pal/palPvobs.c b/ast/pal/palPvobs.c
new file mode 100644
index 0000000..763d202
--- /dev/null
+++ b/ast/pal/palPvobs.c
@@ -0,0 +1,108 @@
+/*
+*+
+*  Name:
+*     palPvobs
+
+*  Purpose:
+*     Position and velocity of an observing station.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     palPvobs( double p, double h, double stl, double pv[6] )
+
+*  Arguments:
+*     p = double (Given)
+*        Latitude (geodetic, radians).
+*     h = double (Given)
+*        Height above reference spheroid (geodetic, metres).
+*     stl = double (Given)
+*        Local apparent sidereal time (radians).
+*     pv = double[ 6 ] (Returned)
+*        position/velocity 6-vector (AU, AU/s, true equator
+*                                    and equinox of date).
+
+*  Description:
+*     Returns the position and velocity of an observing station.
+
+*  Notes:
+*     - The WGS84 reference ellipsoid is used.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-16 (DSB):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "palmac.h"
+#include "palmac.h"
+#include "pal1sofa.h"
+
+void palPvobs( double p, double h, double stl, double pv[6] ){
+
+/* Local Variables: */
+   double xyz[3], z, r, s, c, v;
+
+/* Geodetic to geocentric conversion (WGS84 reference ellipsoid). */
+   eraGd2gc( ERFA_WGS84, 0.0, p, h, xyz );
+
+/* Convert from metres to AU */
+   r = xyz[ 0 ]/ERFA_DAU;
+   z = xyz[ 2 ]/ERFA_DAU;
+
+/* Functions of ST. */
+   s = sin( stl );
+   c = cos( stl );
+
+/* Speed. */
+   v = PAL__SR*r;
+
+/* Position. */
+   pv[ 0 ] = r*c;
+   pv[ 1 ] = r*s;
+   pv[ 2 ] = z;
+
+/* Velocity. */
+   pv[ 3 ] = -v*s;
+   pv[ 4 ] = v*c;
+   pv[ 5 ] = 0.0;
+
+}
+
+
diff --git a/ast/pal/palRvgalc.c b/ast/pal/palRvgalc.c
new file mode 100644
index 0000000..2f01576
--- /dev/null
+++ b/ast/pal/palRvgalc.c
@@ -0,0 +1,111 @@
+/*
+*+
+*  Name:
+*     palRvgalc
+
+*  Purpose:
+*     Velocity component in a given direction due to the rotation
+*     of the Galaxy.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     double palRvgalc( double r2000, double d2000 )
+
+*  Arguments:
+*     r2000 = double (Given)
+*        J2000.0 mean RA (radians)
+*     d2000 = double (Given)
+*        J2000.0 mean Dec (radians)
+
+*  Returned Value:
+*     Component of dynamical LSR motion in direction R2000,D2000 (km/s).
+
+*  Description:
+*     This function returns the Component of dynamical LSR motion in
+*     the direction of R2000,D2000. The result is +ve when the dynamical
+*     LSR is receding from the given point on the sky.
+*
+*  Notes:
+*     - The Local Standard of Rest used here is a point in the
+*     vicinity of the Sun which is in a circular orbit around
+*     the Galactic centre.  Sometimes called the "dynamical" LSR,
+*     it is not to be confused with a "kinematical" LSR, which
+*     is the mean standard of rest of star catalogues or stellar
+*     populations.
+*
+*  Reference:
+*     - The orbital speed of 220 km/s used here comes from Kerr &
+*     Lynden-Bell (1986), MNRAS, 221, p1023.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-16 (DSB):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+double palRvgalc( double r2000, double d2000 ){
+
+/* Local Variables: */
+   double vb[ 3 ];
+
+/*
+*  LSR velocity due to Galactic rotation
+*
+*  Speed = 220 km/s
+*  Apex  = L2,B2  90deg, 0deg
+*        = RA,Dec  21 12 01.1  +48 19 47  J2000.0
+*
+*  This is expressed in the form of a J2000.0 x,y,z vector:
+*
+*      VA(1) = X = -SPEED*COS(RA)*COS(DEC)
+*      VA(2) = Y = -SPEED*SIN(RA)*COS(DEC)
+*      VA(3) = Z = -SPEED*SIN(DEC)
+*/
+
+   double va[ 3 ] = { -108.70408, +97.86251, -164.33610 };
+
+/* Convert given J2000 RA,Dec to x,y,z. */
+   eraS2c( r2000, d2000, vb );
+
+/* Compute dot product with LSR motion vector. */
+   return eraPdp( va, vb );
+}
diff --git a/ast/pal/palRvlg.c b/ast/pal/palRvlg.c
new file mode 100644
index 0000000..255801e
--- /dev/null
+++ b/ast/pal/palRvlg.c
@@ -0,0 +1,106 @@
+/*
+*+
+*  Name:
+*     palRvlg
+
+*  Purpose:
+*     Velocity component in a given direction due to Galactic rotation
+*     and motion of the local group.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     double palRvlg( double r2000, double d2000 )
+
+*  Arguments:
+*     r2000 = double (Given)
+*        J2000.0 mean RA (radians)
+*     d2000 = double (Given)
+*        J2000.0 mean Dec (radians)
+
+*  Returned Value:
+*     Component of SOLAR motion in direction R2000,D2000 (km/s).
+
+*  Description:
+*     This function returns the velocity component in a given
+*     direction due to the combination of the rotation of the
+*     Galaxy and the motion of the Galaxy relative to the mean
+*     motion of the local group. The result is +ve when the Sun
+*     is receding from the given point on the sky.
+*
+*  Reference:
+*     - IAU Trans 1976, 168, p201.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-16 (DSB):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+double palRvlg( double r2000, double d2000 ){
+
+/* Local Variables: */
+   double vb[ 3 ];
+
+/*
+*
+*  Solar velocity due to Galactic rotation and translation
+*
+*  Speed = 300 km/s
+*
+*  Apex  = L2,B2  90deg, 0deg
+*        = RA,Dec  21 12 01.1  +48 19 47  J2000.0
+*
+*  This is expressed in the form of a J2000.0 x,y,z vector:
+*
+*      VA(1) = X = -SPEED*COS(RA)*COS(DEC)
+*      VA(2) = Y = -SPEED*SIN(RA)*COS(DEC)
+*      VA(3) = Z = -SPEED*SIN(DEC)
+*/
+
+   double va[ 3 ] = { -148.23284, +133.44888, -224.09467 };
+
+/* Convert given J2000 RA,Dec to x,y,z. */
+   eraS2c( r2000, d2000, vb );
+
+/* Compute dot product with Solar motion vector. */
+   return eraPdp( va, vb );
+}
diff --git a/ast/pal/palRvlsrd.c b/ast/pal/palRvlsrd.c
new file mode 100644
index 0000000..6302c4f
--- /dev/null
+++ b/ast/pal/palRvlsrd.c
@@ -0,0 +1,116 @@
+/*
+*+
+*  Name:
+*     palRvlsrd
+
+*  Purpose:
+*     Velocity component in a given direction due to the Sun's motion
+*     with respect to the dynamical Local Standard of Rest.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     double palRvlsrd( double r2000, double d2000 )
+
+*  Arguments:
+*     r2000 = double (Given)
+*        J2000.0 mean RA (radians)
+*     d2000 = double (Given)
+*        J2000.0 mean Dec (radians)
+
+*  Returned Value:
+*     Component of "peculiar" solar motion in direction R2000,D2000 (km/s).
+
+*  Description:
+*     This function returns the velocity component in a given direction
+*     due to the Sun's motion with respect to the dynamical Local Standard
+*     of Rest. The result is +ve when the Sun is receding from the given
+*     point on the sky.
+
+*  Notes:
+*     - The Local Standard of Rest used here is the "dynamical" LSR,
+*     a point in the vicinity of the Sun which is in a circular orbit
+*     around the Galactic centre.  The Sun's motion with respect to the
+*     dynamical LSR is called the "peculiar" solar motion.
+*     - There is another type of LSR, called a "kinematical" LSR.  A
+*     kinematical LSR is the mean standard of rest of specified star
+*     catalogues or stellar populations, and several slightly different
+*     kinematical LSRs are in use.  The Sun's motion with respect to an
+*     agreed kinematical LSR is known as the "standard" solar motion.
+*     To obtain a radial velocity correction with respect to an adopted
+*     kinematical LSR use the routine palRvlsrk.
+
+*  Reference:
+*     - Delhaye (1965), in "Stars and Stellar Systems", vol 5, p73.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-16 (DSB):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+double palRvlsrd( double r2000, double d2000 ){
+
+/* Local Variables: */
+   double vb[ 3 ];
+
+/*
+*  Peculiar solar motion from Delhaye 1965: in Galactic Cartesian
+*  coordinates (+9,+12,+7) km/s.  This corresponds to about 16.6 km/s
+*  towards Galactic coordinates L2 = 53 deg, B2 = +25 deg, or RA,Dec
+*  17 49 58.7 +28 07 04 J2000.
+*
+*  The solar motion is expressed here in the form of a J2000.0
+*  equatorial Cartesian vector:
+*
+*      VA(1) = X = -SPEED*COS(RA)*COS(DEC)
+*      VA(2) = Y = -SPEED*SIN(RA)*COS(DEC)
+*      VA(3) = Z = -SPEED*SIN(DEC)
+*/
+
+   double va[ 3 ] = { +0.63823, +14.58542, -7.80116 };
+
+/* Convert given J2000 RA,Dec to x,y,z. */
+   eraS2c( r2000, d2000, vb );
+
+/* Compute dot product with Solar motion vector. */
+   return eraPdp( va, vb );
+}
diff --git a/ast/pal/palRvlsrk.c b/ast/pal/palRvlsrk.c
new file mode 100644
index 0000000..968188a
--- /dev/null
+++ b/ast/pal/palRvlsrk.c
@@ -0,0 +1,116 @@
+/*
+*+
+*  Name:
+*     palRvlsrk
+
+*  Purpose:
+*     Velocity component in a given direction due to the Sun's motion
+*     with respect to an adopted kinematic Local Standard of Rest.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     double palRvlsrk( double r2000, double d2000 )
+
+*  Arguments:
+*     r2000 = double (Given)
+*        J2000.0 mean RA (radians)
+*     d2000 = double (Given)
+*        J2000.0 mean Dec (radians)
+
+*  Returned Value:
+*     Component of "standard" solar motion in direction R2000,D2000 (km/s).
+
+*  Description:
+*     This function returns the velocity component in a given direction
+*     due to the Sun's motion with respect to an adopted kinematic
+*     Local Standard of Rest. The result is +ve when the Sun is receding
+*     from the given point on the sky.
+
+*  Notes:
+*     - The Local Standard of Rest used here is one of several
+*     "kinematical" LSRs in common use.  A kinematical LSR is the mean
+*     standard of rest of specified star catalogues or stellar
+*     populations.  The Sun's motion with respect to a kinematical LSR
+*     is known as the "standard" solar motion.
+*     - There is another sort of LSR, the "dynamical" LSR, which is a
+*     point in the vicinity of the Sun which is in a circular orbit
+*     around the Galactic centre.  The Sun's motion with respect to
+*     the dynamical LSR is called the "peculiar" solar motion.  To
+*     obtain a radial velocity correction with respect to the
+*     dynamical LSR use the routine palRvlsrd.
+
+*  Reference:
+*     - Delhaye (1965), in "Stars and Stellar Systems", vol 5, p73.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     2012-02-16 (DSB):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+double palRvlsrk( double r2000, double d2000 ){
+
+/* Local Variables: */
+   double vb[ 3 ];
+
+/*
+*  Standard solar motion (from Methods of Experimental Physics, ed Meeks,
+*  vol 12, part C, sec 6.1.5.2, p281):
+*
+*  20 km/s towards RA 18h Dec +30d (1900).
+*
+*  The solar motion is expressed here in the form of a J2000.0
+*  equatorial Cartesian vector:
+*
+*      VA(1) = X = -SPEED*COS(RA)*COS(DEC)
+*      VA(2) = Y = -SPEED*SIN(RA)*COS(DEC)
+*      VA(3) = Z = -SPEED*SIN(DEC)
+*/
+
+   double va[ 3 ] = { -0.29000, +17.31726, -10.00141 };
+
+/* Convert given J2000 RA,Dec to x,y,z. */
+   eraS2c( r2000, d2000, vb );
+
+/* Compute dot product with Solar motion vector. */
+   return eraPdp( va, vb );
+}
diff --git a/ast/pal/palSubet.c b/ast/pal/palSubet.c
new file mode 100644
index 0000000..ada122a
--- /dev/null
+++ b/ast/pal/palSubet.c
@@ -0,0 +1,112 @@
+/*
+*+
+*  Name:
+*     palSubet
+
+*  Purpose:
+*     Remove the E-terms from a pre IAU 1976 catalogue RA,Dec
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palSubet ( double rc, double dc, double eq,
+*                     double *rm, double *dm );
+
+*  Arguments:
+*     rc = double (Given)
+*        RA with E-terms included (radians)
+*     dc = double (Given)
+*        Dec with E-terms included (radians)
+*     eq = double (Given)
+*        Besselian epoch of mean equator and equinox
+*     rm = double * (Returned)
+*        RA without E-terms (radians)
+*     dm = double * (Returned)
+*        Dec without E-terms (radians)
+
+*  Description:
+*     Remove the E-terms (elliptic component of annual aberration)
+*     from a pre IAU 1976 catalogue RA,Dec to give a mean place.
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*     Most star positions from pre-1984 optical catalogues (or
+*     derived from astrometry using such stars) embody the
+*     E-terms.  This routine converts such a position to a
+*     formal mean place (allowing, for example, comparison with a
+*     pulsar timing position).
+
+*  See Also:
+*     Explanatory Supplement to the Astronomical Ephemeris,
+*     section 2D, page 48.
+
+*  History:
+*     2012-02-12(TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palSubet ( double rc, double dc, double eq, double *rm, double *dm ) {
+  double a[3];   /* The E-terms */
+  double v[3];
+  double f;
+  int i;
+
+  /* Note the preference for IAU routines */
+
+  /* Retrieve the E-terms */
+  palEtrms( eq, a );
+
+  /* Spherical to Cartesian */
+  eraS2c( rc, dc, v );
+
+  /* Include the E-terms */
+  f = 1.0 + eraPdp( v, a );
+  for (i=0; i<3; i++) {
+    v[i] = f*v[i] - a[i];
+  }
+
+  /* Cartesian to spherical */
+  eraC2s( v, rm, dm );
+
+  /* Bring RA into conventional range */
+  *rm = eraAnp( *rm );
+
+}
diff --git a/ast/pal/palSupgal.c b/ast/pal/palSupgal.c
new file mode 100644
index 0000000..1f4aeb2
--- /dev/null
+++ b/ast/pal/palSupgal.c
@@ -0,0 +1,116 @@
+/*
+*+
+*  Name:
+*     palSupgal
+
+*  Purpose:
+*     Convert from supergalactic to galactic coordinates
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     void palSupgal ( double dsl, double dsb, double *dl, double *db );
+
+*  Arguments:
+*     dsl = double (Given)
+*       Supergalactic longitude.
+*     dsb = double (Given)
+*       Supergalactic latitude.
+*     dl = double * (Returned)
+*       Galactic longitude.
+*     db = double * (Returned)
+*       Galactic latitude.
+
+*  Description:
+*     Transformation from de Vaucouleurs supergalactic coordinates
+*     to IAU 1958 galactic coordinates
+
+*  Authors:
+*     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  See Also:
+*     - de Vaucouleurs, de Vaucouleurs, & Corwin, Second Reference
+*       Catalogue of Bright Galaxies, U. Texas, page 8.
+*     - Systems & Applied Sciences Corp., Documentation for the
+*       machine-readable version of the above catalogue,
+*       Contract NAS 5-26490.
+*
+*     (These two references give different values for the galactic
+*     longitude of the supergalactic origin.  Both are wrong;  the
+*     correct value is L2=137.37.)
+
+*  History:
+*     2012-02-12(TIMJ):
+*        Initial version with documentation taken from Fortran SLA
+*        Adapted with permission from the Fortran SLALIB library.
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 1995 Rutherford Appleton Laboratory
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+#include "pal.h"
+#include "pal1sofa.h"
+
+void palSupgal ( double dsl, double dsb, double *dl, double *db ) {
+
+  double v1[3];
+  double v2[3];
+
+/*
+*  System of supergalactic coordinates:
+*
+*    SGL   SGB        L2     B2      (deg)
+*     -    +90      47.37  +6.32
+*     0     0         -      0
+*
+*  Galactic to supergalactic rotation matrix:
+*/
+  double rmat[3][3] = {
+    { -0.735742574804,+0.677261296414,+0.000000000000 },
+    { -0.074553778365,-0.080991471307,+0.993922590400 },
+    { +0.673145302109,+0.731271165817,+0.110081262225 }
+  };
+
+  /* Spherical to Cartesian */
+  eraS2c( dsl, dsb, v1 );
+
+  /* Supergalactic to galactic */
+  eraTrxp( rmat, v1, v2 );
+
+  /* Cartesian to spherical */
+  eraC2s( v2, dl, db );
+
+  /* Express in conventional ranges */
+  *dl = eraAnp( *dl );
+  *db = eraAnpm( *db );
+
+}
diff --git a/ast/pal/palmac.h b/ast/pal/palmac.h
new file mode 100644
index 0000000..207b843
--- /dev/null
+++ b/ast/pal/palmac.h
@@ -0,0 +1,136 @@
+#ifndef PALMACDEF
+#define PALMACDEF
+
+/*
+*+
+*  Name:
+*     palmac.h
+
+*  Purpose:
+*     Macros used by the PAL library
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Include file
+
+*  Description:
+*     A collection of useful macros provided and used by the PAL library
+
+*  Authors:
+*     TIMJ: Tim Jenness (JAC, Hawaii)
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  Notes:
+*
+
+*  History:
+*     2012-02-08 (TIMJ):
+*        Initial version.
+*        Adapted with permission from the Fortran SLALIB library.
+*     2012-04-13 (DSB):
+*        Added PAL__DR2H and PAL__DR2S
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2012 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software: you can redistribute it and/or
+*     modify it under the terms of the GNU Lesser General Public
+*     License as published by the Free Software Foundation, either
+*     version 3 of the License, or (at your option) any later
+*     version.
+*
+*     This program is distributed in the hope that it will be useful,
+*     but WITHOUT ANY WARRANTY; without even the implied warranty of
+*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*     GNU Lesser General Public License for more details.
+*
+*     You should have received a copy of the GNU Lesser General
+*     License along with this program.  If not, see
+*     <http://www.gnu.org/licenses/>.
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+/* Pi */
+static const double PAL__DPI = 3.1415926535897932384626433832795028841971693993751;
+
+/* 2Pi */
+static const double PAL__D2PI = 6.2831853071795864769252867665590057683943387987502;
+
+/* pi/2:  90 degrees in radians */
+static const double PAL__DPIBY2 = 1.5707963267948966192313216916397514420985846996876;
+
+/* pi/180:  degrees to radians */
+static const double PAL__DD2R = 0.017453292519943295769236907684886127134428718885417;
+
+/* Radians to arcseconds */
+static const double PAL__DR2AS = 2.0626480624709635515647335733077861319665970087963e5;
+
+/* Arcseconds to radians */
+static const double PAL__DAS2R = 4.8481368110953599358991410235794797595635330237270e-6;
+
+/* Radians to degrees */
+static const double PAL__DR2D = 57.295779513082320876798154814105170332405472466564;
+
+/* Hours to radians */
+static const double PAL__DH2R = 0.26179938779914943653855361527329190701643078328126;
+
+/* Radians to hours */
+static const double PAL__DR2H = 3.8197186342054880584532103209403446888270314977709;
+
+/* Radians to seconds of time */
+static const double PAL__DR2S = 1.3750987083139757010431557155385240879777313391975e4;
+
+/* Seconds of time to radians */
+static const double PAL__DS2R = 7.272205216643039903848712e-5;
+
+/* Start of SLA modified Julian date epoch */
+static const double PAL__MJD0 = 2400000.5;
+
+/* Light time for 1 AU (sec) */
+static const double PAL__CR = 499.004782;
+
+/* Seconds per day */
+static const double PAL__SPD = 86400.0;
+
+/* Km per sec to AU per tropical century
+   = 86400 * 36524.2198782 / 149597870 */
+static const double PAL__VF = 21.095;
+
+/*  Radians per year to arcsec per century. This needs to be a macro since it
+    is an expression including other constants. */
+#define PAL__PMF (100.0*60.0*60.0*360.0/PAL__D2PI);
+
+/* Mean sidereal rate - the rotational angular velocity of Earth
+   in radians/sec from IERS Conventions (2003). */
+static const double PAL__SR = 7.2921150e-5;
+
+/*  Gaussian gravitational constant (exact) */
+static const double PAL__GCON = 0.01720209895;
+
+/* DINT(A) - truncate to nearest whole number towards zero (double) */
+#define DINT(A) ((A)<0.0?ceil(A):floor(A))
+
+/* DNINT(A) - round to nearest whole number (double) */
+#define DNINT(A) ((A)<0.0?ceil((A)-0.5):floor((A)+0.5))
+
+/* DMAX(A,B) - return maximum value - evaluates arguments multiple times */
+#define DMAX(A,B) ((A) > (B) ? (A) : (B) )
+
+/* DMIN(A,B) - return minimum value - evaluates arguments multiple times */
+#define DMIN(A,B) ((A) < (B) ? (A) : (B) )
+
+/* We actually prefer to use C99 copysign() but we define this here as a backup
+   but it will not detect -0.0 so is not useful for palDfltin. */
+/* DSIGN(A,B) - magnitude of A with sign of B (double) */
+#define DSIGN(A,B) ((B)<0.0?-fabs(A):fabs(A))
+
+#endif