diff options
Diffstat (limited to 'ast/pal')
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]) &rms[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, &rms[ 1 ], ehd, eh ); + +/* Heliocentric direction of Earth (normalized) and modulus. */ + eraPn( eh, &e, &rms[ 4 ] ); + +/* Light deflection parameter */ + amprms[7] = GR2 / e; + +/* Aberration parameters. */ + for( i = 0; i < 3; i++ ) { + amprms[ i + 8 ] = ebd[ i ]*PAL__CR; + } + eraPn( &rms[8], &vm, vn ); + amprms[ 11 ] = sqrt( 1.0 - vm*vm ); + +/* NPB matrix. */ + palPrenut( eq, date, (double(*)[ 3 ]) &rms[ 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]) &rms[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 |