diff options
Diffstat (limited to 'wcs/dateutil.c')
| -rw-r--r-- | wcs/dateutil.c | 4554 |
1 files changed, 4554 insertions, 0 deletions
diff --git a/wcs/dateutil.c b/wcs/dateutil.c new file mode 100644 index 0000000..ada0c95 --- /dev/null +++ b/wcs/dateutil.c @@ -0,0 +1,4554 @@ +/*** File libwcs/dateutil.c + *** October 19, 2012 + *** By Jessica Mink, jmink@cfa.harvard.edu + *** Harvard-Smithsonian Center for Astrophysics + *** Copyright (C) 1999-2012 + *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA + + This library 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 2 of the License, or (at your option) any later version. + + This library 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 Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + Correspondence concerning WCSTools should be addressed as follows: + Internet email: jmink@cfa.harvard.edu + Postal address: Jessica Mink + Smithsonian Astrophysical Observatory + 60 Garden St. + Cambridge, MA 02138 USA + */ + +/* Date and time conversion routines using the following conventions: + ang = Angle in fractional degrees + deg = Angle in degrees as dd:mm:ss.ss + doy = 2 floating point numbers: year and day, including fraction, of year + *** First day of year is 1, not zero. + dt = 2 floating point numbers: yyyy.mmdd, hh.mmssssss + ep = fractional year, often epoch of a position including proper motion + epb = Besselian epoch = 365.242198781-day years based on 1900.0 + epj = Julian epoch = 365.25-day years based on 2000.0 + fd = FITS date string which may be any of the following: + yyyy.ffff (fractional year) + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard FITS use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) + hr = Sexigesimal hours as hh:mm:dd.ss + jd = Julian Date + lt = Local time + mjd = modified Julian Date = JD - 2400000.5 + ofd = FITS date string (dd/mm/yy before 2000, else no return) + time = use fd2* with no date to convert time as hh:mm:ss.ss to sec, day, year + ts = UT seconds since 1950-01-01T00:00 (used for ephemeris computations) + tsi = local seconds since 1980-01-01T00:00 (used by IRAF as a time tag) + tsu = UT seconds since 1970-01-01T00:00 (used as Unix system time) + tsd = UT seconds of current day + ut = Universal Time (UTC) + et = Ephemeris Time (or TDB or TT) = TAI + 32.184 seconds + tai = International Atomic Time (Temps Atomique International) = ET - 32.184 seconds + gps = GPS time = TAI - 19 seconds + mst = Mean Greenwich Sidereal Time + gst = Greenwich Sidereal Time (includes nutation) + lst = Local Sidereal Time (includes nutation) (longitude must be set) + hjd = Heliocentric Julian Date + mhjd = modified Heliocentric Julian Date = HJD - 2400000.5 + + * ang2hr (angle) + * Convert angle in decimal floating point degrees to hours as hh:mm:ss.ss + * ang2deg (angle) + * Convert angle in decimal floating point degrees to degrees as dd:mm:ss.ss + * deg2ang (angle as dd:mm:ss.ss) + * Convert angle in degrees as dd:mm:ss.ss to decimal floating point degrees + * ang2hr (angle) + * Convert angle in hours as hh:mm:ss.ss to decimal floating point degrees + * + * doy2dt (year, doy, date, time) + * Convert year and day of year to date as yyyy.ddmm and time as hh.mmsss + * doy2ep, doy2epb, doy2epj (date, time) + * Convert year and day of year to fractional year + * doy2fd (year, doy) + * Convert year and day of year to FITS date string + * doy2mjd (year, doy) + * Convert year and day of year to modified Julian date + * + * dt2doy (date, time, year, doy) + * Convert date as yyyy.ddmm and time as hh.mmsss to year and day of year + * dt2ep, dt2epb, dt2epj (date, time) + * Convert date as yyyy.ddmm and time as hh.mmsss to fractional year + * dt2fd (date, time) + * Convert date as yyyy.ddmm and time as hh.mmsss to FITS date string + * dt2i (date,time,iyr,imon,iday,ihr,imn,sec, ndsec) + * Convert yyyy.mmdd hh.mmssss to year month day hours minutes seconds + * dt2jd (date,time) + * Convert date as yyyy.ddmm and time as hh.mmsss to Julian date + * dt2mjd (date,time) + * Convert date as yyyy.ddmm and time as hh.mmsss to modified Julian date + * dt2ts (date,time) + * Convert date (yyyy.ddmm) and time (hh.mmsss) to seconds since 1950-01-01 + * dt2tsi (date,time) + * Convert date (yyyy.ddmm) and time (hh.mmsss) to seconds since 1980-01-01 + * dt2tsu (date,time) + * Convert date (yyyy.ddmm) and time (hh.mmsss) to seconds since 1970-01-01 + * + * ep2dt, epb2dt, epj2dt (epoch,date, time) + * Convert fractional year to date as yyyy.ddmm and time as hh.mmsss + * ep2fd, epb2fd, epj2fd (epoch) + * Convert epoch to FITS ISO date string + * ep2i, epb2i, epj2i (epoch,iyr,imon,iday,ihr,imn,sec, ndsec) + * Convert fractional year to year month day hours minutes seconds + * ep2jd, epb2jd, epj2jd (epoch) + * Convert fractional year as used in epoch to Julian date + * ep2mjd, epb2mjd, epj2mjd (epoch) + * Convert fractional year as used in epoch to modified Julian date + * ep2ts, epb2ts, epj2ts (epoch) + * Convert fractional year to seconds since 1950.0 + * + * et2fd (string) + * Convert from ET (or TDT or TT) in FITS format to UT in FITS format + * fd2et (string) + * Convert from UT in FITS format to ET (or TDT or TT) in FITS format + * jd2jed (dj) + * Convert from Julian Date to Julian Ephemeris Date + * jed2jd (dj) + * Convert from Julian Ephemeris Date to Julian Date + * dt2et (date, time) + * Convert date (yyyy.ddmm) and time (hh.mmsss) to ephemeris time + * edt2dt (date, time) + * Convert ephemeris date (yyyy.ddmm) and time (hh.mmsss) to UT + * dt2tai (date, time) + * Convert date (yyyy.ddmm) and time (hh.mmsss) to TAI date and time + * tai2dt (date, time) + * Convert TAI date (yyyy.ddmm) and time (hh.mmsss) to UT + * ts2ets (tsec) + * Convert from UT in seconds since 1950-01-01 to ET in same format + * ets2ts (tsec) + * Convert from ET in seconds since 1950-01-01 to UT in same format + * + * fd2ep, fd2epb, fd2epj (string) + * Convert FITS date string to fractional year + * Convert time alone to fraction of Besselian year + * fd2doy (string, year, doy) + * Convert FITS standard date string to year and day of year + * fd2dt (string, date, time) + * Convert FITS date string to date as yyyy.ddmm and time as hh.mmsss + * Convert time alone to hh.mmssss with date set to 0.0 + * fd2i (string,iyr,imon,iday,ihr,imn,sec, ndsec) + * Convert FITS standard date string to year month day hours min sec + * Convert time alone to hours min sec, year month day are zero + * fd2jd (string) + * Convert FITS standard date string to Julian date + * Convert time alone to fraction of day + * fd2mjd (string) + * Convert FITS standard date string to modified Julian date + * fd2ts (string) + * Convert FITS standard date string to seconds since 1950.0 + * Convert time alone to seconds of day + * fd2fd (string) + * Convert FITS standard date string to ISO FITS date string + * fd2of (string) + * Convert FITS standard date string to old-format FITS date and time + * fd2ofd (string) + * Convert FITS standard date string to old-format FITS date string + * fd2oft (string) + * Convert time part of FITS standard date string to FITS date string + * + * jd2doy (dj, year, doy) + * Convert Julian date to year and day of year + * jd2dt (dj,date,time) + * Convert Julian date to date as yyyy.mmdd and time as hh.mmssss + * jd2ep, jd2epb, jd2epj (dj) + * Convert Julian date to fractional year as used in epoch + * jd2fd (dj) + * Convert Julian date to FITS ISO date string + * jd2i (dj,iyr,imon,iday,ihr,imn,sec, ndsec) + * Convert Julian date to year month day hours min sec + * jd2mjd (dj) + * Convert Julian date to modified Julian date + * jd2ts (dj) + * Convert Julian day to seconds since 1950.0 + * + * lt2dt() + * Return local time as yyyy.mmdd and time as hh.mmssss + * lt2fd() + * Return local time as FITS ISO date string + * lt2tsi() + * Return local time as IRAF seconds since 1980-01-01 00:00 + * lt2tsu() + * Return local time as Unix seconds since 1970-01-01 00:00 + * lt2ts() + * Return local time as Unix seconds since 1950-01-01 00:00 + * + * mjd2doy (dj,year,doy) + * Convert modified Julian date to date as year and day of year + * mjd2dt (dj,date,time) + * Convert modified Julian date to date as yyyy.mmdd and time as hh.mmssss + * mjd2ep, mjd2epb, mjd2epj (dj) + * Convert modified Julian date to fractional year as used in epoch + * mjd2fd (dj) + * Convert modified Julian date to FITS ISO date string + * mjd2i (dj,iyr,imon,iday,ihr,imn,sec, ndsec) + * Convert modified Julian date to year month day hours min sec + * mjd2jd (dj) + * Convert modified Julian date to Julian date + * mjd2ts (dj) + * Convert modified Julian day to seconds since 1950.0 + * + * ts2dt (tsec,date,time) + * Convert seconds since 1950.0 to date as yyyy.ddmm and time as hh.mmsss + * ts2ep, ts2epb, ts2epj (tsec) + * Convert seconds since 1950.0 to fractional year + * ts2fd (tsec) + * Convert seconds since 1950.0 to FITS standard date string + * ts2i (tsec,iyr,imon,iday,ihr,imn,sec, ndsec) + * Convert sec since 1950.0 to year month day hours minutes seconds + * ts2jd (tsec) + * Convert seconds since 1950.0 to Julian date + * ts2mjd (tsec) + * Convert seconds since 1950.0 to modified Julian date + * tsi2fd (tsec) + * Convert seconds since 1980-01-01 to FITS standard date string + * tsi2dt (tsec,date,time) + * Convert seconds since 1980-01-01 to date as yyyy.ddmm, time as hh.mmsss + * tsu2fd (tsec) + * Convert seconds since 1970-01-01 to FITS standard date string + * tsu2tsi (tsec) + * Convert UT seconds since 1970-01-01 to local seconds since 1980-01-01 + * tsu2dt (tsec,date,time) + * Convert seconds since 1970-01-01 to date as yyyy.ddmm, time as hh.mmsss + * + * tsd2fd (tsec) + * Convert seconds since start of day to FITS time, hh:mm:ss.ss + * tsd2dt (tsec) + * Convert seconds since start of day to hh.mmssss + * + * fd2gst (string) + * convert from FITS date Greenwich Sidereal Time + * dt2gst (date, time) + * convert from UT as yyyy.mmdd hh.mmssss to Greenwich Sidereal Time + * ts2gst (tsec) + * Calculate Greenwich Sidereal Time given Universal Time + * in seconds since 1951-01-01T0:00:00 + * fd2mst (string) + * convert from FITS UT date to Mean Sidereal Time + * dt2gmt (date, time) + * convert from UT as yyyy.mmdd hh.mmssss to Mean Sidereal Time + * ts2mst (tsec) + * Calculate Mean Sidereal Time given Universal Time + * in seconds since 1951-01-01T0:00:00 + * jd2mst (string) + * convert from Julian Date to Mean Sidereal Time + * mst2fd (string) + * convert to current UT in FITS format given Greenwich Mean Sidereal Time + * mst2jd (dj) + * convert to current UT as Julian Date given Greenwich Mean Sidereal Time + * jd2lst (dj) + * Calculate Local Sidereal Time from Julian Date + * ts2lst (tsec) + * Calculate Local Sidereal Time given UT in seconds since 1951-01-01T0:00 + * fd2lst (string) + * Calculate Local Sidereal Time given Universal Time as FITS ISO date + * lst2jd (dj, lst) + * Calculate Julian Date given current Julian date and Local Sidereal Time + * lst2fd (string, lst) + * Calculate Julian Date given current UT date and Local Sidereal Time + * gst2fd (string) + * Calculate current UT given UT date and Greenwich Sidereal Time + * gst2jd (dj) + * Calculate current UT given UT date and Greenwich Sidereal Time as JD + * + * compnut (dj, dpsi, deps, eps0) + * Compute the longitude and obliquity components of nutation and + * mean obliquity from the IAU 1980 theory + * + * utdt (dj) + * Compute difference between UT and dynamical time (ET-UT) + * ut2dt (year, doy) + * Current Universal Time to year and day of year + * ut2dt (date, time) + * Current Universal Time to date (yyyy.mmdd) and time (hh.mmsss) + * ut2ep(), ut2epb(), ut2epj() + * Current Universal Time to fractional year, Besselian, Julian epoch + * ut2fd() + * Current Universal Time to FITS ISO date string + * ut2jd() + * Current Universal Time to Julian Date + * ut2mjd() + * Current Universal Time to Modified Julian Date + * ut2tsi() + * Current Universal Time to IRAF seconds since 1980-01-01T00:00 + * ut2tsu() + * Current Universal Time to Unix seconds since 1970-01-01T00:00 + * ut2ts() + * Current Universal Time to seconds since 1950-01-01T00:00 + * isdate (string) + * Return 1 if string is a FITS date (old or ISO) + * + * Internally-used subroutines + * + * fixdate (iyr, imon, iday, ihr, imn, sec, ndsec) + * Round seconds and make sure date and time numbers are within limits + * caldays (year, month) + * Calculate days in month 1-12 given year (Gregorian calendar only + * dint (dnum) + * Return integer part of floating point number + * dmod (dnum) + * Return Mod of floating point number + */ + +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <math.h> +#include <time.h> +#include <sys/time.h> +#include "wcs.h" +#include "fitsfile.h" + +static double suntl(); +static void fixdate(); +static int caldays(); +static double dint(); +static double dmod(); + +static double longitude = 0.0; /* longitude of observatory in degrees (+=west) */ +void +setlongitude (longitude0) +double longitude0; +{ longitude = longitude0; return; } + +static int ndec = 3; +void +setdatedec (nd) +int nd; +{ ndec = nd; return; } + +/* ANG2HR -- Convert angle in fraction degrees to hours as hh:mm:ss.ss */ + +void +ang2hr (angle, lstr, string) + +double angle; /* Angle in fractional degrees */ +int lstr; /* Maximum number of characters in string */ +char *string; /* Character string (hh:mm:ss.ss returned) */ + +{ + angle = angle / 15.0; + dec2str (string, lstr, angle, ndec); + return; +} + + +/* ANG2DEG -- Convert angle in fraction degrees to degrees as dd:mm:ss.ss */ + +void +ang2deg (angle, lstr, string) + +double angle; /* Angle in fractional degrees */ +int lstr; /* Maximum number of characters in string */ +char *string; /* Character string (dd:mm:ss.ss returned) */ +{ + dec2str (string, lstr, angle, ndec); + return; +} + + +/* DEG2ANG -- Convert angle in degrees as dd:mm:ss.ss to fractional degrees */ + +double +deg2ang (angle) + +char *angle; /* Angle as dd:mm:ss.ss */ +{ + double deg; + + deg = str2dec (angle); + return (deg); +} + +/* HR2ANG -- Convert angle in hours as hh:mm:ss.ss to fractional degrees */ + +double +hr2ang (angle) + +char *angle; /* Angle in sexigesimal hours (hh:mm:ss.sss) */ + +{ + double deg; + + deg = str2dec (angle); + deg = deg * 15.0; + return (deg); +} + + +/* DT2FD-- convert vigesimal date and time to FITS date, yyyy-mm-ddThh:mm:ss.ss */ + +char * +dt2fd (date, time) + +double date; /* Date as yyyy.mmdd + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + int iyr,imon,iday,ihr,imn; + double sec; + int nf; + char *string; + char tstring[32], dstring[32]; + char outform[64]; + + dt2i (date, time, &iyr,&imon,&iday,&ihr,&imn,&sec, ndec); + + /* Convert to ISO date format */ + string = (char *) calloc (32, sizeof (char)); + + /* Make time string */ + if (time != 0.0 || ndec > 0) { + if (ndec == 0) + nf = 2; + else + nf = 3 + ndec; + if (ndec > 0) { + sprintf (outform, "%%02d:%%02d:%%0%d.%df", nf, ndec); + sprintf (tstring, outform, ihr, imn, sec); + } + else { + sprintf (outform, "%%02d:%%02d:%%0%dd", nf); + sprintf (tstring, outform, ihr, imn, (int)(sec+0.5)); + } + } + + /* Make date string */ + if (date != 0.0) + sprintf (dstring, "%4d-%02d-%02d", iyr, imon, iday); + + /* Make FITS (ISO) date string */ + if (date == 0.0) + strcpy (string, tstring); + else if (time == 0.0 && ndec < 1) + strcpy (string, dstring); + else + sprintf (string, "%sT%s", dstring, tstring); + + return (string); +} + + +/* DT2JD-- convert from date as yyyy.mmdd and time as hh.mmsss to Julian Date + * Return fractional days if date is zero */ + +double +dt2jd (date,time) + +double date; /* Date as yyyy.mmdd + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double dj; /* Julian date (returned) */ + double tsec; /* seconds since 1950.0 */ + + tsec = dt2ts (date, time); + if (date == 0.0) + dj = tsec / 86400.0; + else + dj = ts2jd (tsec); + + return (dj); +} + + +/* DT2MJD-- convert from date yyyy.mmdd time hh.mmsss to modified Julian Date + * Return fractional days if date is zero */ + +double +dt2mjd (date,time) + +double date; /* Date as yyyy.mmdd + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double dj; /* Modified Julian date (returned) */ + double tsec; /* seconds since 1950.0 */ + + tsec = dt2ts (date, time); + if (date == 0.0) + dj = tsec / 86400.0; + else + dj = ts2jd (tsec); + + return (dj - 2400000.5); +} + + +/* HJD2JD-- convert Heliocentric Julian Date to (geocentric) Julian date */ + +double +hjd2jd (dj, ra, dec, sys) + +double dj; /* Heliocentric Julian date */ +double ra; /* Right ascension (degrees) */ +double dec; /* Declination (degrees) */ +int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */ +{ + double lt; /* Light travel difference to the Sun (days) */ + + lt = suntl (dj, ra, dec, sys); + + /* Return Heliocentric Julian Date */ + return (dj - lt); +} + + +/* JD2HJD-- convert (geocentric) Julian date to Heliocentric Julian Date */ + +double +jd2hjd (dj, ra, dec, sys) + +double dj; /* Julian date (geocentric) */ +double ra; /* Right ascension (degrees) */ +double dec; /* Declination (degrees) */ +int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */ +{ + double lt; /* Light travel difference to the Sun (days) */ + + lt = suntl (dj, ra, dec, sys); + + /* Return Heliocentric Julian Date */ + return (dj + lt); +} + + +/* MHJD2MJD-- convert modified Heliocentric Julian Date to + modified geocentric Julian date */ + +double +mhjd2mjd (mhjd, ra, dec, sys) + +double mhjd; /* Modified Heliocentric Julian date */ +double ra; /* Right ascension (degrees) */ +double dec; /* Declination (degrees) */ +int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */ +{ + double lt; /* Light travel difference to the Sun (days) */ + double hjd; /* Heliocentric Julian date */ + + hjd = mjd2jd (mhjd); + + lt = suntl (hjd, ra, dec, sys); + + /* Return Heliocentric Julian Date */ + return (jd2mjd (hjd - lt)); +} + + +/* MJD2MHJD-- convert modified geocentric Julian date tp + modified Heliocentric Julian Date */ + +double +mjd2mhjd (mjd, ra, dec, sys) + +double mjd; /* Julian date (geocentric) */ +double ra; /* Right ascension (degrees) */ +double dec; /* Declination (degrees) */ +int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */ +{ + double lt; /* Light travel difference to the Sun (days) */ + double dj; /* Julian date (geocentric) */ + + dj = mjd2jd (mjd); + + lt = suntl (dj, ra, dec, sys); + + /* Return Heliocentric Julian Date */ + return (jd2mjd (dj + lt)); +} + + +/* SUNTL-- compute light travel time to heliocentric correction in days */ +/* Translated into C from IRAF SPP noao.astutils.asttools.asthjd.x */ + +static double +suntl (dj, ra, dec, sys) + +double dj; /* Julian date (geocentric) */ +double ra; /* Right ascension (degrees) */ +double dec; /* Declination (degrees) */ +int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */ +{ + double t; /* Number of Julian centuries since J1900 */ + double manom; /* Mean anomaly of the Earth's orbit (degrees) */ + double lperi; /* Mean longitude of perihelion (degrees) */ + double oblq; /* Mean obliquity of the ecliptic (degrees) */ + double eccen; /* Eccentricity of the Earth's orbit (dimensionless) */ + double eccen2, eccen3; + double tanom; /* True anomaly (approximate formula) (radians) */ + double slong; /* True longitude of the Sun from the Earth (radians) */ + double rs; /* Distance to the sun (AU) */ + double lt; /* Light travel difference to the Sun (days) */ + double l; /* Longitude of star in orbital plane of Earth (radians) */ + double b; /* Latitude of star in orbital plane of Earth (radians) */ + double epoch; /* Epoch of obervation */ + double rs1,rs2; + + t = (dj - 2415020.0) / 36525.0; + + /* Compute earth orbital parameters */ + manom = 358.47583 + (t * (35999.04975 - t * (0.000150 + t * 0.000003))); + lperi = 101.22083 + (t * (1.7191733 + t * (0.000453 + t * 0.000003))); + oblq = 23.452294 - (t * (0.0130125 + t * (0.00000164 - t * 0.000000503))); + eccen = 0.01675104 - (t * (0.00004180 + t * 0.000000126)); + eccen2 = eccen * eccen; + eccen3 = eccen * eccen2; + + /* Convert to principle angles */ + manom = manom - (360.0 * (dint) (manom / 360.0)); + lperi = lperi - (360.0 * (dint) (lperi / 360.0)); + + /* Convert to radians */ + manom = degrad (manom); + lperi = degrad (lperi); + oblq = degrad (oblq); + + /* True anomaly */ + tanom = manom + (2 * eccen - 0.25 * eccen3) * sin (manom) + + 1.25 * eccen2 * sin (2 * manom) + + 13./12. * eccen3 * sin (3 * manom); + + /* Distance to the Sun */ + rs1 = 1.0 - eccen2; + rs2 = 1.0 + (eccen * cos (tanom)); + rs = rs1 / rs2; + + /* True longitude of the Sun seen from the Earth */ + slong = lperi + tanom + PI; + + /* Longitude and latitude of star in orbital plane of the Earth */ + epoch = jd2ep (dj); + wcscon (sys, WCS_ECLIPTIC, 0.0, 0.0, &ra, &dec, epoch); + l = degrad (ra); + b = degrad (dec); + + /* Light travel difference to the Sun */ + lt = -0.005770 * rs * cos (b) * cos (l - slong); + + /* Return light travel difference */ + return (lt); +} + + +/* JD2DT-- convert Julian date to date as yyyy.mmdd and time as hh.mmssss */ + +void +jd2dt (dj,date,time) + +double dj; /* Julian date */ +double *date; /* Date as yyyy.mmdd (returned) */ +double *time; /* Time as hh.mmssxxxx (returned) */ +{ + int iyr,imon,iday,ihr,imn; + double sec; + + /* Convert Julian Date to date and time */ + jd2i (dj, &iyr, &imon, &iday, &ihr, &imn, &sec, 4); + + /* Convert date to yyyy.mmdd */ + if (iyr < 0) { + *date = (double) (-iyr) + 0.01 * (double) imon + 0.0001 * (double) iday; + *date = -(*date); + } + else + *date = (double) iyr + 0.01 * (double) imon + 0.0001 * (double) iday; + + /* Convert time to hh.mmssssss */ + *time = (double) ihr + 0.01 * (double) imn + 0.0001 * sec; + + return; +} + + +/* JD2I-- convert Julian date to date as year, month, and day, and time hours, + minutes, and seconds */ +/* after Fliegel and Van Flander, CACM 11, 657 (1968) */ + + +void +jd2i (dj, iyr, imon, iday, ihr, imn, sec, ndsec) + +double dj; /* Julian date */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ + +{ + double tsec; + double frac, dts, ts, sday; + int jd, l, n, i, j; + + tsec = jd2ts (dj); + /* ts2i (tsec, iyr, imon, iday, ihr, imn, sec, ndsec); */ + + /* Round seconds to 0 - 4 decimal places */ + if (tsec < 0.0) + dts = -0.5; + else + dts = 0.5; + if (ndsec < 1) + ts = dint (tsec + dts); + else if (ndsec < 2) + ts = dint (tsec * 10.0 + dts) / 10.0; + else if (ndsec < 3) + ts = dint (tsec * 100.0 + dts) / 100.0; + else if (ndsec < 4) + ts = dint (tsec * 1000.0 + dts) / 1000.0; + else + ts = dint (tsec * 10000.0 + dts) / 10000.0; + + /* Convert back to Julian Date */ + dj = ts2jd (ts); + + /* Compute time from fraction of a day */ + frac = dmod (dj, 1.0); + if (frac < 0.5) { + jd = (int) (dj - frac); + sday = (frac + 0.5) * 86400.0; + } + else { + jd = (int) (dj - frac) + 1; + sday = (frac - 0.5) * 86400.0; + } + + *ihr = (int) (sday / 3600.0); + sday = sday - (double) (*ihr * 3600); + *imn = (int) (sday / 60.0); + *sec = sday - (double) (*imn * 60); + + /* Compute day, month, year */ + l = jd + 68569; + n = (4 * l) / 146097; + l = l - (146097 * n + 3) / 4; + i = (4000 * (l + 1)) / 1461001; + l = l - (1461 * i) / 4 + 31; + j = (80 * l) / 2447; + *iday = l - (2447 * j) / 80; + l = j / 11; + *imon = j + 2 - (12 * l); + *iyr = 100 * (n - 49) + i + l; + + return; +} + + +/* JD2MJD-- convert Julian Date to Modified Julian Date */ + +double +jd2mjd (dj) + +double dj; /* Julian Date */ + +{ + return (dj - 2400000.5); +} + + +/* JD2EP-- convert Julian date to fractional year as used in epoch */ + +double +jd2ep (dj) + +double dj; /* Julian date */ + +{ + double date, time; + jd2dt (dj, &date, &time); + return (dt2ep (date, time)); +} + + +/* JD2EPB-- convert Julian date to Besselian epoch */ + +double +jd2epb (dj) + +double dj; /* Julian date */ + +{ + return (1900.0 + (dj - 2415020.31352) / 365.242198781); +} + + +/* JD2EPJ-- convert Julian date to Julian epoch */ + +double +jd2epj (dj) + +double dj; /* Julian date */ + +{ + return (2000.0 + (dj - 2451545.0) / 365.25); +} + + +/* LT2DT-- Return local time as yyyy.mmdd and time as hh.mmssss */ + +void +lt2dt(date, time) + +double *date; /* Date as yyyy.mmdd (returned) */ +double *time; /* Time as hh.mmssxxxx (returned) */ + +{ + time_t tsec; + struct timeval tp; + struct timezone tzp; + struct tm *ts; + + gettimeofday (&tp,&tzp); + + tsec = tp.tv_sec; + ts = localtime (&tsec); + + if (ts->tm_year < 1000) + *date = (double) (ts->tm_year + 1900); + else + *date = (double) ts->tm_year; + *date = *date + (0.01 * (double) (ts->tm_mon + 1)); + *date = *date + (0.0001 * (double) ts->tm_mday); + *time = (double) ts->tm_hour; + *time = *time + (0.01 * (double) ts->tm_min); + *time = *time + (0.0001 * (double) ts->tm_sec); + + return; +} + + +/* LT2FD-- Return current local time as FITS ISO date string */ + +char * +lt2fd() +{ + time_t tsec; + struct tm *ts; + struct timeval tp; + struct timezone tzp; + int month, day, year, hour, minute, second; + char *isotime; + + gettimeofday (&tp,&tzp); + tsec = tp.tv_sec; + + ts = localtime (&tsec); + + year = ts->tm_year; + if (year < 1000) + year = year + 1900; + month = ts->tm_mon + 1; + day = ts->tm_mday; + hour = ts->tm_hour; + minute = ts->tm_min; + second = ts->tm_sec; + + isotime = (char *) calloc (32, sizeof (char)); + sprintf (isotime, "%04d-%02d-%02dT%02d:%02d:%02d", + year, month, day, hour, minute, second); + return (isotime); +} + + +/* LT2TSI-- Return local time as IRAF seconds since 1980-01-01 00:00 */ + +int +lt2tsi() +{ + return ((int)(lt2ts() - 946684800.0)); +} + + +/* LT2TSU-- Return local time as Unix seconds since 1970-01-01 00:00 */ + +time_t +lt2tsu() +{ + return ((time_t)(lt2ts() - 631152000.0)); +} + +/* LT2TS-- Return local time as Unix seconds since 1950-01-01 00:00 */ + +double +lt2ts() +{ + double tsec; + char *datestring; + datestring = lt2fd(); + tsec = fd2ts (datestring); + free (datestring); + return (tsec); +} + + +/* MJD2DT-- convert Modified Julian Date to date (yyyy.mmdd) time (hh.mmssss) */ + +void +mjd2dt (dj,date,time) + +double dj; /* Modified Julian Date */ +double *date; /* Date as yyyy.mmdd (returned) + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double *time; /* Time as hh.mmssxxxx (returned) + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double tsec; + + tsec = jd2ts (dj + 2400000.5); + ts2dt (tsec, date, time); + + return; +} + + +/* MJD2I-- convert Modified Julian Date to date as year, month, day and + time as hours, minutes, seconds */ + +void +mjd2i (dj, iyr, imon, iday, ihr, imn, sec, ndsec) + +double dj; /* Modified Julian Date */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ + +{ + double tsec; + + tsec = jd2ts (dj + 2400000.5); + ts2i (tsec, iyr, imon, iday, ihr, imn, sec, ndsec); + return; +} + + +/* MJD2DOY-- convert Modified Julian Date to Year,Day-of-Year */ + +void +mjd2doy (dj, year, doy) + +double dj; /* Modified Julian Date */ +int *year; /* Year (returned) */ +double *doy; /* Day of year with fraction (returned) */ + +{ + jd2doy (dj + 2400000.5, year, doy); + return; +} + + +/* MJD2JD-- convert Modified Julian Date to Julian Date */ + +double +mjd2jd (dj) + +double dj; /* Modified Julian Date */ + +{ + return (dj + 2400000.5); +} + + +/* MJD2EP-- convert Modified Julian Date to fractional year */ + +double +mjd2ep (dj) + +double dj; /* Modified Julian Date */ + +{ + double date, time; + jd2dt (dj + 2400000.5, &date, &time); + return (dt2ep (date, time)); +} + + +/* MJD2EPB-- convert Modified Julian Date to Besselian epoch */ + +double +mjd2epb (dj) + +double dj; /* Modified Julian Date */ + +{ + return (1900.0 + (dj - 15019.81352) / 365.242198781); +} + + +/* MJD2EPJ-- convert Modified Julian Date to Julian epoch */ + +double +mjd2epj (dj) + +double dj; /* Modified Julian Date */ + +{ + return (2000.0 + (dj - 51544.5) / 365.25); +} + + +/* MJD2FD-- convert modified Julian date to FITS date, yyyy-mm-ddThh:mm:ss.ss */ + +char * +mjd2fd (dj) + +double dj; /* Modified Julian date */ +{ + return (jd2fd (dj + 2400000.5)); +} + + +/* MJD2TS-- convert modified Julian date to seconds since 1950.0 */ + +double +mjd2ts (dj) + +double dj; /* Modified Julian date */ +{ + return ((dj - 33282.0) * 86400.0); +} + + +/* EP2FD-- convert fractional year to FITS date, yyyy-mm-ddThh:mm:ss.ss */ + +char * +ep2fd (epoch) + +double epoch; /* Date as fractional year */ +{ + double tsec; /* seconds since 1950.0 (returned) */ + tsec = ep2ts (epoch); + return (ts2fd (tsec)); +} + + +/* EPB2FD-- convert Besselian epoch to FITS date, yyyy-mm-ddThh:mm:ss.ss */ + +char * +epb2fd (epoch) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ +{ + double dj; /* Julian Date */ + dj = epb2jd (epoch); + return (jd2fd (dj)); +} + + +/* EPJ2FD-- convert Julian epoch to FITS date, yyyy-mm-ddThh:mm:ss.ss */ + +char * +epj2fd (epoch) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ +{ + double dj; /* Julian Date */ + dj = epj2jd (epoch); + return (jd2fd (dj)); +} + + +/* EP2TS-- convert fractional year to seconds since 1950.0 */ + +double +ep2ts (epoch) + +double epoch; /* Date as fractional year */ +{ + double dj; + dj = ep2jd (epoch); + return ((dj - 2433282.5) * 86400.0); +} + + +/* EPB2TS-- convert Besselian epoch to seconds since 1950.0 */ + +double +epb2ts (epoch) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ +{ + double dj; + dj = epb2jd (epoch); + return ((dj - 2433282.5) * 86400.0); +} + + +/* EPJ2TS-- convert Julian epoch to seconds since 1950.0 */ + +double +epj2ts (epoch) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ +{ + double dj; + dj = epj2jd (epoch); + return ((dj - 2433282.5) * 86400.0); +} + + +/* EPB2EP-- convert Besselian epoch to fractional years */ + +double +epb2ep (epoch) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ +{ + double dj; + dj = epb2jd (epoch); + return (jd2ep (dj)); +} + + +/* EP2EPB-- convert fractional year to Besselian epoch */ + +double +ep2epb (epoch) + +double epoch; /* Fractional year */ +{ + double dj; + dj = ep2jd (epoch); + return (jd2epb (dj)); +} + + +/* EPJ2EP-- convert Julian epoch to fractional year */ + +double +epj2ep (epoch) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ +{ + double dj; + dj = epj2jd (epoch); + return (jd2ep (dj)); +} + + +/* EP2EPJ-- convert fractional year to Julian epoch */ + +double +ep2epj (epoch) + +double epoch; /* Fractional year */ +{ + double dj; + dj = ep2jd (epoch); + return (jd2epj (dj)); +} + + +/* EP2I-- convert fractional year to year month day hours min sec */ + +void +ep2i (epoch, iyr, imon, iday, ihr, imn, sec, ndsec) + +double epoch; /* Date as fractional year */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ +{ + double date, time; + + ep2dt (epoch, &date, &time); + dt2i (date, time, iyr,imon,iday,ihr,imn,sec, ndsec); + return; +} + + +/* EPB2I-- convert Besselian epoch to year month day hours min sec */ + +void +epb2i (epoch, iyr, imon, iday, ihr, imn, sec, ndsec) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ +{ + double date, time; + + epb2dt (epoch, &date, &time); + dt2i (date, time, iyr,imon,iday,ihr,imn,sec, ndsec); + return; +} + + +/* EPJ2I-- convert Julian epoch to year month day hours min sec */ + +void +epj2i (epoch, iyr, imon, iday, ihr, imn, sec, ndsec) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ +{ + double date, time; + + epj2dt (epoch, &date, &time); + dt2i (date, time, iyr,imon,iday,ihr,imn,sec, ndsec); + return; +} + + +/* EP2JD-- convert fractional year as used in epoch to Julian date */ + +double +ep2jd (epoch) + +double epoch; /* Date as fractional year */ + +{ + double dj; /* Julian date (returned)*/ + double date, time; + + ep2dt (epoch, &date, &time); + dj = dt2jd (date, time); + return (dj); +} + + +/* EPB2JD-- convert Besselian epoch to Julian Date */ + +double +epb2jd (epoch) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ + +{ + return (2415020.31352 + ((epoch - 1900.0) * 365.242198781)); +} + + +/* EPJ2JD-- convert Julian epoch to Julian Date */ + +double +epj2jd (epoch) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ + +{ + return (2451545.0 + ((epoch - 2000.0) * 365.25)); +} + + +/* EP2MJD-- convert fractional year as used in epoch to modified Julian date */ + +double +ep2mjd (epoch) + +double epoch; /* Date as fractional year */ + +{ + double dj; /* Julian date (returned)*/ + double date, time; + + ep2dt (epoch, &date, &time); + dj = dt2jd (date, time); + return (dj - 2400000.5); +} + + +/* EPB2MJD-- convert Besselian epoch to modified Julian Date */ + +double +epb2mjd (epoch) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ + +{ + return (15019.81352 + ((epoch - 1900.0) * 365.242198781)); +} + + +/* EPJ2MJD-- convert Julian epoch to modified Julian Date */ + +double +epj2mjd (epoch) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ + +{ + return (51544.5 + ((epoch - 2000.0) * 365.25)); +} + + + +/* EPB2EPJ-- convert Besselian epoch to Julian epoch */ + +double +epb2epj (epoch) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ +{ + double dj; /* Julian date */ + dj = epb2jd (epoch); + return (jd2epj (dj)); +} + + +/* EPJ2EPB-- convert Julian epoch to Besselian epoch */ + +double +epj2epb (epoch) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ +{ + double dj; /* Julian date */ + dj = epj2jd (epoch); + return (jd2epb (dj)); +} + + +/* JD2FD-- convert Julian date to FITS date, yyyy-mm-ddThh:mm:ss.ss */ + +char * +jd2fd (dj) + +double dj; /* Julian date */ +{ + double tsec; /* seconds since 1950.0 (returned) */ + tsec = (dj - 2433282.5) * 86400.0; + return (ts2fd (tsec)); +} + + +/* JD2TS-- convert Julian date to seconds since 1950.0 */ + +double +jd2ts (dj) + +double dj; /* Julian date */ +{ + return ((dj - 2433282.5) * 86400.0); +} + + +/* JD2TSI-- convert Julian date to IRAF seconds since 1980-01-01T0:00 */ + +int +jd2tsi (dj) + +double dj; /* Julian date */ +{ + double ts; + ts = (dj - 2444239.5) * 86400.0; + return ((int) ts); +} + + +/* JD2TSU-- convert Julian date to Unix seconds since 1970-01-01T0:00 */ + +time_t +jd2tsu (dj) + +double dj; /* Julian date */ +{ + return ((time_t)((dj - 2440587.5) * 86400.0)); +} + + +/* DT2DOY-- convert yyyy.mmdd hh.mmss to year and day of year */ + +void +dt2doy (date, time, year, doy) + +double date; /* Date as yyyy.mmdd */ +double time; /* Time as hh.mmssxxxx */ +int *year; /* Year (returned) */ +double *doy; /* Day of year with fraction (returned) */ +{ + double dj; /* Julian date */ + double dj0; /* Julian date on January 1 0:00 */ + double date0; /* January first of date's year */ + double dyear; + + dyear = floor (date); + date0 = dyear + 0.0101; + dj0 = dt2jd (date0, 0.0); + dj = dt2jd (date, time); + *year = (int) (dyear + 0.00000001); + *doy = dj - dj0 + 1.0; + return; +} + + +/* DOY2DT-- convert year and day of year to yyyy.mmdd hh.mmss */ + +void +doy2dt (year, doy, date, time) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +double *date; /* Date as yyyy.mmdd (returned) */ +double *time; /* Time as hh.mmssxxxx (returned) */ +{ + double dj; /* Julian date */ + double dj0; /* Julian date on January 1 0:00 */ + double date0; /* January first of date's year */ + + date0 = year + 0.0101; + dj0 = dt2jd (date0, 0.0); + dj = dj0 + doy - 1.0; + jd2dt (dj, date, time); + return; +} + + +/* DOY2EP-- convert year and day of year to fractional year as used in epoch */ + +double +doy2ep (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double date, time; + doy2dt (year, doy, &date, &time); + return (dt2ep (date, time)); +} + + + +/* DOY2EPB-- convert year and day of year to Besellian epoch */ + +double +doy2epb (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj; + dj = doy2jd (year, doy); + return (jd2epb (dj)); +} + + +/* DOY2EPJ-- convert year and day of year to Julian epoch */ + +double +doy2epj (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj; + dj = doy2jd (year, doy); + return (jd2epj (dj)); +} + + +/* DOY2FD-- convert year and day of year to FITS date */ + +char * +doy2fd (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj; /* Julian date */ + + dj = doy2jd (year, doy); + return (jd2fd (dj)); +} + + +/* DOY2JD-- convert year and day of year to Julian date */ + +double +doy2jd (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj0; /* Julian date */ + double date; /* Date as yyyy.mmdd (returned) */ + double time; /* Time as hh.mmssxxxx (returned) */ + + date = (double) year + 0.0101; + time = 0.0; + dj0 = dt2jd (date, time); + return (dj0 + doy - 1.0); +} + + +/* DOY2MJD-- convert year and day of year to Julian date */ + +double +doy2mjd (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj0; /* Julian date */ + double date; /* Date as yyyy.mmdd (returned) */ + double time; /* Time as hh.mmssxxxx (returned) */ + + date = (double) year + 0.0101; + time = 0.0; + dj0 = dt2jd (date, time); + return (dj0 + doy - 1.0 - 2400000.5); +} + + +/* DOY2TSU-- convert from FITS date to Unix seconds since 1970-01-01T0:00 */ + +time_t +doy2tsu (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj; + dj = doy2jd (year, doy); + return ((time_t)jd2ts (dj)); +} + + +/* DOY2TSI-- convert from FITS date to IRAF seconds since 1980-01-01T0:00 */ + +int +doy2tsi (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj; + dj = doy2jd (year, doy); + return ((int)jd2tsi (dj)); +} + + +/* DOY2TS-- convert year, day of year to seconds since 1950 */ + +double +doy2ts (year, doy) + +int year; /* Year */ +double doy; /* Day of year with fraction */ +{ + double dj; + dj = doy2jd (year, doy); + return (jd2ts (dj)); +} + + +/* FD2DOY-- convert FITS date to year and day of year */ + +void +fd2doy (string, year, doy) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +int *year; /* Year (returned) */ +double *doy; /* Day of year with fraction (returned) */ +{ + double dj; /* Julian date */ + + dj = fd2jd (string); + jd2doy (dj, year, doy); + return; +} + + +/* JD2DOY-- convert Julian date to year and day of year */ + +void +jd2doy (dj, year, doy) + +double dj; /* Julian date */ +int *year; /* Year (returned) */ +double *doy; /* Day of year with fraction (returned) */ +{ + double date; /* Date as yyyy.mmdd (returned) */ + double time; /* Time as hh.mmssxxxx (returned) */ + double dj0; /* Julian date at 0:00 on 1/1 */ + double dyear; + + jd2dt (dj, &date, &time); + *year = (int) date; + dyear = (double) *year; + dj0 = dt2jd (dyear+0.0101, 0.0); + *doy = dj - dj0 + 1.0; + return; +} + + +/* TS2JD-- convert seconds since 1950.0 to Julian date */ + +double +ts2jd (tsec) + +double tsec; /* seconds since 1950.0 */ +{ + return (2433282.5 + (tsec / 86400.0)); +} + + +/* TS2MJD-- convert seconds since 1950.0 to modified Julian date */ + +double +ts2mjd (tsec) + +double tsec; /* seconds since 1950.0 */ +{ + return (33282.0 + (tsec / 86400.0)); +} + + +/* TS2EP-- convert seconds since 1950.0 to fractional year as used in epoch */ + +double +ts2ep (tsec) + +double tsec; /* Seconds since 1950.0 */ + +{ + double date, time; + ts2dt (tsec, &date, &time); + return (dt2ep (date, time)); +} + + +/* TS2EPB-- convert seconds since 1950.0 to Besselian epoch */ + +double +ts2epb (tsec) + +double tsec; /* Seconds since 1950.0 */ + +{ + double dj; /* Julian Date */ + dj = ts2jd (tsec); + return (jd2epb (dj)); +} + + +/* TS2EPB-- convert seconds since 1950.0 to Julian epoch */ + +double +ts2epj (tsec) + +double tsec; /* Seconds since 1950.0 */ + +{ + double dj; /* Julian Date */ + dj = ts2jd (tsec); + return (jd2epj (dj)); +} + + +/* DT2EP-- convert from date, time as yyyy.mmdd hh.mmsss to fractional year */ + +double +dt2ep (date, time) + +double date; /* Date as yyyy.mmdd + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double epoch; /* Date as fractional year (returned) */ + double dj, dj0, dj1, date0, time0, date1; + + dj = dt2jd (date, time); + if (date == 0.0) + epoch = dj / 365.2422; + else { + time0 = 0.0; + date0 = dint (date) + 0.0101; + date1 = dint (date) + 1.0101; + dj0 = dt2jd (date0, time0); + dj1 = dt2jd (date1, time0); + epoch = dint (date) + ((dj - dj0) / (dj1 - dj0)); + } + return (epoch); +} + + +/* DT2EPB-- convert from date, time as yyyy.mmdd hh.mmsss to Besselian epoch */ + +double +dt2epb (date, time) + +double date; /* Date as yyyy.mmdd + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double dj; /* Julian date */ + double epoch; /* Date as fractional year (returned) */ + dj = dt2jd (date, time); + if (date == 0.0) + epoch = dj / 365.242198781; + else + epoch = jd2epb (dj); + return (epoch); +} + + +/* DT2EPJ-- convert from date, time as yyyy.mmdd hh.mmsss to Julian epoch */ + +double +dt2epj (date, time) + +double date; /* Date as yyyy.mmdd + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double dj; /* Julian date */ + double epoch; /* Date as fractional year (returned) */ + dj = dt2jd (date, time); + if (date == 0.0) + epoch = dj / 365.25; + else + epoch = jd2epj (dj); + return (epoch); +} + + +/* EP2DT-- convert from fractional year to date, time as yyyy.mmdd hh.mmsss */ + +void +ep2dt (epoch, date, time) + +double epoch; /* Date as fractional year */ +double *date; /* Date as yyyy.mmdd (returned) + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double *time; /* Time as hh.mmssxxxx (returned) + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double dj, dj0, dj1, date0, time0, date1, epochi, epochf; + + time0 = 0.0; + epochi = dint (epoch); + epochf = epoch - epochi; + date0 = epochi + 0.0101; + date1 = epochi + 1.0101; + dj0 = dt2jd (date0, time0); + dj1 = dt2jd (date1, time0); + dj = dj0 + epochf * (dj1 - dj0); + jd2dt (dj, date, time); + return; +} + + +/* EPB2DT-- convert from Besselian epoch to date, time as yyyy.mmdd hh.mmsss */ + +void +epb2dt (epoch, date, time) + +double epoch; /* Besselian epoch (fractional 365.242198781-day years) */ +double *date; /* Date as yyyy.mmdd (returned) + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double *time; /* Time as hh.mmssxxxx (returned) + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double dj; /* Julian date */ + dj = epb2jd (epoch); + jd2dt (dj, date, time); +} + + +/* EPJ2DT-- convert from Julian epoch to date, time as yyyy.mmdd hh.mmsss */ + +void +epj2dt (epoch, date, time) + +double epoch; /* Julian epoch (fractional 365.25-day years) */ +double *date; /* Date as yyyy.mmdd (returned) + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double *time; /* Time as hh.mmssxxxx (returned) + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double dj; /* Julian date */ + dj = epj2jd (epoch); + jd2dt (dj, date, time); +} + + +/* FD2JD-- convert FITS standard date to Julian date */ + +double +fd2jd (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double date, time; + + fd2dt (string, &date, &time); + return (dt2jd (date, time)); +} + + +/* FD2MJD-- convert FITS standard date to modified Julian date */ + +double +fd2mjd (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + return (fd2jd (string) - 2400000.5); +} + + +/* FD2TSU-- convert from FITS date to Unix seconds since 1970-01-01T0:00 */ + +time_t +fd2tsu (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double date, time; + fd2dt (string, &date, &time); + return (dt2tsu (date, time)); +} + + +/* FD2TSI-- convert from FITS date to IRAF seconds since 1980-01-01T0:00 */ + +int +fd2tsi (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double date, time; + fd2dt (string, &date, &time); + return (dt2tsi (date, time)); +} + + +/* FD2TS-- convert FITS standard date to seconds since 1950 */ + +double +fd2ts (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double date, time; + fd2dt (string, &date, &time); + return (dt2ts (date, time)); +} + + +/* FD2FD-- convert any FITS standard date to ISO FITS standard date */ + +char * +fd2fd (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double date, time; + fd2dt (string, &date, &time); + return (dt2fd (date, time)); +} + + +/* FD2OF-- convert any FITS standard date to old FITS standard date time */ + +char * +fd2of (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + int iyr,imon,iday,ihr,imn; + double sec; + + fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 3); + + /* Convert to old FITS date format */ + string = (char *) calloc (32, sizeof (char)); + if (iyr < 1900) + sprintf (string, "*** date out of range ***"); + else if (iyr < 2000) + sprintf (string, "%02d/%02d/%02d %02d:%02d:%06.3f", + iday, imon, iyr-1900, ihr, imn, sec); + else if (iyr < 2900.0) + sprintf (string, "%02d/%02d/%3d %02d:%02d:%6.3f", + iday, imon, iyr-1900, ihr, imn, sec); + else + sprintf (string, "*** date out of range ***"); + return (string); +} + + +/* TAI-UTC from the U.S. Naval Observatory */ +/* ftp://maia.usno.navy.mil/ser7/tai-utc.dat */ +static double taijd[26]={2441317.5, 2441499.5, 2441683.5, 2442048.5, 2442413.5, + 2442778.5, 2443144.5, 2443509.5, 2443874.5, 2444239.5, 2444786.5, + 2445151.5, 2445516.5, 2446247.5, 2447161.5, 2447892.5, 2448257.5, + 2448804.5, 2449169.5, 2449534.5, 2450083.5, 2450630.5, 2451179.5, + 2453736.5, 2454832.5, 2456293.5}; +static double taidt[26]={10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0, + 20.0,21.0,22.0,23.0,24.0,25.0,26.0,27.0,28.0,29.0,30.0,31.0,32.0, + 33.0,34.0,35.0}; +static double dttab[173]={13.7,13.4,13.1,12.9,12.7,12.6,12.5,12.5,12.5,12.5, + 12.5,12.5,12.5,12.5,12.5,12.5,12.5,12.4,12.3,12.2,12.0,11.7,11.4, + 11.1,10.6,10.2, 9.6, 9.1, 8.6, 8.0, 7.5, 7.0, 6.6, 6.3, 6.0, 5.8, + 5.7, 5.6, 5.6, 5.6, 5.7, 5.8, 5.9, 6.1, 6.2, 6.3, 6.5, 6.6, 6.8, + 6.9, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.7, 7.8, 7.8,7.88,7.82, + 7.54, 6.97, 6.40, 6.02, 5.41, 4.10, 2.92, 1.82, 1.61, 0.10,-1.02, + -1.28,-2.69,-3.24,-3.64,-4.54,-4.71,-5.11,-5.40,-5.42,-5.20,-5.46, + -5.46,-5.79,-5.63,-5.64,-5.80,-5.66,-5.87,-6.01,-6.19,-6.64,-6.44, + -6.47,-6.09,-5.76,-4.66,-3.74,-2.72,-1.54,-0.02, 1.24, 2.64, 3.86, + 5.37, 6.14, 7.75, 9.13,10.46,11.53,13.36,14.65,16.01,17.20,18.24, + 19.06,20.25,20.95,21.16,22.25,22.41,23.03,23.49,23.62,23.86,24.49, + 24.34,24.08,24.02,24.00,23.87,23.95,23.86,23.93,23.73,23.92,23.96, + 24.02,24.33,24.83,25.30,25.70,26.24,26.77,27.28,27.78,28.25,28.71, + 29.15,29.57,29.97,30.36,30.72,31.07,31.35,31.68,32.18,32.68,33.15, + 33.59,34.00,34.47,35.03,35.73,36.54,37.43,38.29,39.20,40.18,41.17, + 42.23}; + + +/* TAI2FD-- convert from TAI in FITS format to UT in FITS format */ + +char * +tai2fd (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double dj0, dj, tsec, dt; + + dj0 = fd2jd (string); + dt = utdt (dj0); + dj = dj0 - (dt / 86400.0); + dt = utdt (dj); + tsec = fd2ts (string); + tsec = tsec - dt + 32.184; + return (ts2fd (tsec)); +} + + +/* FD2TAI-- convert from UT in FITS format to TAI in FITS format */ + +char * +fd2tai (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double dj, tsec, dt; + + dj = fd2jd (string); + dt = utdt (dj); + tsec = fd2ts (string); + tsec = tsec + dt - 32.184; + return (ts2fd (tsec)); +} + + +/* DT2TAI-- convert from UT as yyyy.mmdd hh.mmssss to TAI in same format */ + +void +dt2tai (date, time) +double *date; /* Date as yyyy.mmdd */ +double *time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + double dj, dt, tsec; + + dj = dt2jd (*date, *time); + dt = utdt (dj); + tsec = dt2ts (*date, *time); + tsec = tsec + dt - 32.184; + ts2dt (tsec, date, time); + return; +} + + +/* TAI2DT-- convert from TAI as yyyy.mmdd hh.mmssss to UT in same format */ + +void +tai2dt (date, time) +double *date; /* Date as yyyy.mmdd */ +double *time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + double dj, dt, tsec, tsec0; + + dj = dt2jd (*date, *time); + dt = utdt (dj); + tsec0 = dt2ts (*date, *time); + tsec = tsec0 + dt; + dj = ts2jd (tsec); + dt = utdt (dj); + tsec = tsec0 + dt + 32.184; + ts2dt (tsec, date, time); + return; +} + + +/* ET2FD-- convert from ET (or TDT or TT) in FITS format to UT in FITS format */ + +char * +et2fd (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double dj0, dj, tsec, dt; + + dj0 = fd2jd (string); + dt = utdt (dj0); + dj = dj0 - (dt / 86400.0); + dt = utdt (dj); + tsec = fd2ts (string); + tsec = tsec - dt; + return (ts2fd (tsec)); +} + + +/* FD2ET-- convert from UT in FITS format to ET (or TDT or TT) in FITS format */ + +char * +fd2et (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double dj, tsec, dt; + + dj = fd2jd (string); + dt = utdt (dj); + tsec = fd2ts (string); + tsec = tsec + dt; + return (ts2fd (tsec)); +} + + +/* DT2ET-- convert from UT as yyyy.mmdd hh.mmssss to ET in same format */ + +void +dt2et (date, time) +double *date; /* Date as yyyy.mmdd */ +double *time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + double dj, dt, tsec; + + dj = dt2jd (*date, *time); + dt = utdt (dj); + tsec = dt2ts (*date, *time); + tsec = tsec + dt; + ts2dt (tsec, date, time); + return; +} + + +/* EDT2DT-- convert from ET as yyyy.mmdd hh.mmssss to UT in same format */ + +void +edt2dt (date, time) +double *date; /* Date as yyyy.mmdd */ +double *time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + double dj, dt, tsec, tsec0; + + dj = dt2jd (*date, *time); + dt = utdt (dj); + tsec0 = dt2ts (*date, *time); + tsec = tsec0 + dt; + dj = ts2jd (tsec); + dt = utdt (dj); + tsec = tsec0 + dt; + ts2dt (tsec, date, time); + return; +} + + +/* JD2JED-- convert from Julian Date to Julian Ephemeris Date */ + +double +jd2jed (dj) + +double dj; /* Julian Date */ +{ + double dt; + + dt = utdt (dj); + return (dj + (dt / 86400.0)); +} + + +/* JED2JD-- convert from Julian Ephemeris Date to Julian Date */ + +double +jed2jd (dj) + +double dj; /* Julian Ephemeris Date */ +{ + double dj0, dt; + + dj0 = dj; + dt = utdt (dj); + dj = dj0 - (dt / 86400.0); + dt = utdt (dj); + return (dj - (dt / 86400.0)); +} + + +/* TS2ETS-- convert from UT in seconds since 1950-01-01 to ET in same format */ + +double +ts2ets (tsec) + +double tsec; +{ + double dj, dt; + + dj = ts2jd (tsec); + dt = utdt (dj); + return (tsec + dt); +} + + +/* ETS2TS-- convert from ET in seconds since 1950-01-01 to UT in same format */ + +double +ets2ts (tsec) + +double tsec; +{ + double dj, dj0, dt; + + dj0 = ts2jd (tsec); + dt = utdt (dj0); + dj = dj0 - (dt / 86400.0); + dt = utdt (dj); + return (tsec - dt); +} + + +/* UTDT-- Compute difference between UT and dynamical time (ET-UT) */ + +double +utdt (dj) + +double dj; /* Julian Date (UT) */ +{ + double dt, date, time, ts, ts1, ts0, date0, yfrac, diff, cj; + int i, iyr, iyear; + + /* If after 1972-01-01, use tabulated TAI-UT */ + if (dj >= 2441317.5) { + dt = 0.0; + for (i = 22; i > 0; i--) { + if (dj >= taijd[i]) + dt = taidt[i]; + } + dt = dt + 32.184; + } + + /* For 1800-01-01 to 1972-01-01, use table of ET-UT from AE */ + else if (dj >= 2378496.5) { + jd2dt (dj, &date, &time); + ts = jd2ts (dj); + iyear = (int) date; + iyr = iyear - 1800; + date0 = (double) iyear + 0.0101; + ts0 = dt2ts (date0, 0.0); + date0 = (double) (iyear + 1) + 0.0101; + ts1 = dt2ts (date0, 0.0); + yfrac = (ts - ts0) / (ts1 - ts0); + diff = dttab[iyr+1] - dttab[iyr]; + dt = dttab[iyr] + (diff * yfrac); + } + + /* Compute back to 1600 using formula from McCarthy and Babcock (1986) */ + else if (dj >= 2305447.5) { + cj = (dj - 2378496.5) / 36525.0; + dt = 5.156 + 13.3066 * (cj - 0.19) * (cj - 0.19); + } + + /* Compute back to 948 using formula from Stephenson and Morrison (1984) */ + else if (dj >= 2067309.5) { + cj = (dj - 2378496.5) / 36525.0; + dt = 25.5 * cj * cj; + } + + /*Compute back to 390 BC using formula from Stephenson and Morrison (1984)*/ + else if (dj >= 0.0) { + cj = (dj = 2378496.5) / 36525.0; + dt = 1360.0 + (320.0 * cj) + (44.3 * cj * cj); + } + + else + dt = 0.0; + return (dt); +} + + +/* FD2OFD-- convert any FITS standard date to old FITS standard date */ + +char * +fd2ofd (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + int iyr,imon,iday,ihr,imn; + double sec; + + fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 3); + + /* Convert to old FITS date format */ + string = (char *) calloc (32, sizeof (char)); + if (iyr < 1900) + sprintf (string, "*** date out of range ***"); + else if (iyr < 2000) + sprintf (string, "%02d/%02d/%02d", iday, imon, iyr-1900); + else if (iyr < 2900.0) + sprintf (string, "%02d/%02d/%3d", iday, imon, iyr-1900); + else + sprintf (string, "*** date out of range ***"); + return (string); +} + + +/* FD2OFT-- convert any FITS standard date to old FITS standard time */ + +char * +fd2oft (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + int iyr,imon,iday,ihr,imn; + double sec; + + fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 3); + + /* Convert to old FITS date format */ + string = (char *) calloc (32, sizeof (char)); + sprintf (string, "%02d:%02d:%06.3f", ihr, imn, sec); + return (string); +} + + +/* FD2DT-- convert FITS standard date to date, time as yyyy.mmdd hh.mmsss */ + +void +fd2dt (string, date, time) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +double *date; /* Date as yyyy.mmdd (returned) + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double *time; /* Time as hh.mmssxxxx (returned) + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + int iyr,imon,iday,ihr,imn; + double sec; + + fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 4); + + /* Convert date to yyyy.mmdd */ + if (iyr < 0) { + *date = (double) (-iyr) + 0.01 * (double) imon + 0.0001 * (double) iday; + *date = -(*date); + } + else + *date = (double) iyr + 0.01 * (double) imon + 0.0001 * (double) iday; + + /* Convert time to hh.mmssssss */ + *time = (double) ihr + 0.01 * (double) imn + 0.0001 * sec; + + return; +} + + +/* FD2EP-- convert from FITS standard date to fractional year */ + +double +fd2ep (string) + +char *string; /* FITS date string, which may be: + yyyy.ffff (fractional year) + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard FITS use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ + +{ + double dj; /* Julian date */ + dj = fd2jd (string); + if (dj < 1.0) + return (dj / 365.2422); + else + return (jd2ep (dj)); +} + + +/* FD2EPB-- convert from FITS standard date to Besselian epoch */ + +double +fd2epb (string) + +char *string; /* FITS date string, which may be: + yyyy.ffff (fractional year) + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard FITS use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ + +{ + double dj; /* Julian date */ + dj = fd2jd (string); + if (dj < 1.0) + return (dj / 365.242198781); + else + return (jd2epb (dj)); +} + + +/* FD2EPJ-- convert from FITS standard date to Julian epoch */ + +double +fd2epj (string) + +char *string; /* FITS date string, which may be: + yyyy.ffff (fractional year) + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard FITS use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ + +{ + double dj; /* Julian date */ + dj = fd2jd (string); + if (dj < 1.0) + return (dj / 365.25); + else + return (jd2epj (dj)); +} + + +/* DT2TSU-- convert from date and time to Unix seconds since 1970-01-01T0:00 */ + +time_t +dt2tsu (date,time) + +double date; /* Date as yyyy.mmdd */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + return ((time_t)(dt2ts (date, time) - 631152000.0)); +} + + +/* DT2TSI-- convert from date and time to IRAF seconds since 1980-01-01T0:00 */ + +int +dt2tsi (date,time) + +double date; /* Date as yyyy.mmdd */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + return ((int)(dt2ts (date, time) - 946684800.0)); +} + + + +/* DT2TS-- convert from date, time as yyyy.mmdd hh.mmsss to sec since 1950.0 */ + +double +dt2ts (date,time) + +double date; /* Date as yyyy.mmdd + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + double tsec; /* Seconds past 1950.0 (returned) */ + + double dh,dm,dd; + int iy,im,id; + +/* Calculate the number of full years, months, and days already + * elapsed since 0h, March 1, -1 (up to most recent midnight). */ + + /* convert time of day to elapsed seconds */ + + /* If time is < 0, it is assumed to be a fractional day */ + if (time < 0.0) + tsec = time * -86400.0; + else { + dh = (int) (time + 0.0000000001); + dm = (int) (((time - dh) * 100.0) + 0.0000000001); + tsec = (time * 10000.0) - (dh * 10000.0) - (dm * 100.0); + tsec = (int) (tsec * 100000.0 + 0.0001) / 100000.0; + tsec = tsec + (dm * 60.0) + (dh * 3600.0); + } + + + /* Calculate the number of full months elapsed since + * the current or most recent March */ + if (date >= 0.0301) { + iy = (int) (date + 0.0000000001); + im = (int) (((date - (double) (iy)) * 10000.0) + 0.00000001); + id = im % 100; + im = (im / 100) + 9; + if (im < 12) iy = iy - 1; + im = im % 12; + id = id - 1; + + /* starting with March as month 0 and ending with the following + * February as month 11, the calculation of the number of days + * per month reduces to a simple formula. the following statement + * determines the number of whole days elapsed since 3/1/-1 and then + * subtracts the 712163 days between then and 1/1/1950. it converts + * the result to seconds and adds the accumulated seconds above. */ + id = id + ((im+1+im/6+im/11)/2 * 31) + ((im-im/6-im/11)/2 * 30) + + (iy / 4) - (iy / 100) + (iy / 400); + dd = (double) id + (365.0 * (double) iy) - 712163.0; + tsec = tsec + (dd * 86400.0); + } + + return (tsec); +} + + +/* TS2DT-- convert seconds since 1950.0 to date, time as yyyy.mmdd hh.mmssss */ + +void +ts2dt (tsec,date,time) + +double tsec; /* Seconds past 1950.0 */ +double *date; /* Date as yyyy.mmdd (returned) + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double *time; /* Time as hh.mmssxxxx (returned) + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +{ + int iyr,imon,iday,ihr,imn; + double sec; + + ts2i (tsec,&iyr,&imon,&iday,&ihr,&imn,&sec, 4); + + /* Convert date to yyyy.mmdd */ + if (iyr < 0) { + *date = (double) (-iyr) + 0.01 * (double) imon + 0.0001 * (double) iday; + *date = -(*date); + } + else + *date = (double) iyr + 0.01 * (double) imon + 0.0001 * (double) iday; + + /* Convert time to hh.mmssssss */ + *time = (double) ihr + 0.01 * (double) imn + 0.0001 * sec; + + return; +} + + +/* TSI2DT-- Convert seconds since 1980-01-01 to date yyyy.ddmm, time hh.mmsss */ + +void +tsi2dt (isec,date,time) + +int isec; /* Seconds past 1980-01-01 */ +double *date; /* Date as yyyy.mmdd (returned) */ +double *time; /* Time as hh.mmssxxxx (returned) */ +{ + ts2dt (tsi2ts (isec), date, time); +} + + +/* TSI2FD-- Convert seconds since 1980-01-01 to FITS standard date string */ + +char * +tsi2fd (isec) + +int isec; /* Seconds past 1980-01-01 */ +{ + return (ts2fd (tsi2ts (isec))); +} + + +/* TSI2TS-- Convert seconds since 1980-01-01 to seconds since 1950-01-01 */ + +double +tsi2ts (isec) +int isec; /* Seconds past 1980-01-01 */ +{ + return ((double) isec + 946684800.0); +} + + +/* TSU2FD-- Convert seconds since 1970-01-01 to FITS standard date string */ + +char * +tsu2fd (isec) +time_t isec; /* Seconds past 1970-01-01 */ +{ + return (ts2fd (tsu2ts (isec))); +} + + +/* TSU2DT-- Convert seconds since 1970-01-01 to date yyyy.ddmm, time hh.mmsss */ + +void +tsu2dt (isec,date,time) +time_t isec; /* Seconds past 1970-01-01 */ +double *date; /* Date as yyyy.mmdd (returned) */ +double *time; /* Time as hh.mmssxxxx (returned) */ +{ + ts2dt (tsu2ts (isec), date, time); +} + + +/* TSU2TS-- Convert seconds since 1970-01-01 to seconds since 1950-01-01 */ + +double +tsu2ts (isec) +time_t isec; /* Seconds past 1970-01-01 */ +{ + return ((double) isec + 631152000.0); +} + +/* TSU2TSI-- UT seconds since 1970-01-01 to local seconds since 1980-01-01 */ + +int +tsu2tsi (isec) +time_t isec; /* Seconds past 1970-01-01 */ +{ + double date, time; + struct tm *ts; + + /* Get local time from UT seconds */ + ts = localtime (&isec); + if (ts->tm_year < 1000) + date = (double) (ts->tm_year + 1900); + else + date = (double) ts->tm_year; + date = date + (0.01 * (double) (ts->tm_mon + 1)); + date = date + (0.0001 * (double) ts->tm_mday); + time = (double) ts->tm_hour; + time = time + (0.01 * (double) ts->tm_min); + time = time + (0.0001 * (double) ts->tm_sec); + return ((int)(dt2ts (date, time) - 631152000.0)); +} + + +/* TS2FD-- convert seconds since 1950.0 to FITS date, yyyy-mm-ddThh:mm:ss.ss */ + +char * +ts2fd (tsec) + +double tsec; /* Seconds past 1950.0 */ +{ + double date, time; + + ts2dt (tsec, &date, &time); + return (dt2fd (date, time)); +} + + +/* TSD2FD-- convert seconds since start of day to FITS time, hh:mm:ss.ss */ + +char * +tsd2fd (tsec) + +double tsec; /* Seconds since start of day */ +{ + double date, time; + char *thms, *fdate; + int lfd, nbc; + + ts2dt (tsec, &date, &time); + fdate = dt2fd (date, time); + thms = (char *) calloc (16, 1); + lfd = strlen (fdate); + nbc = lfd - 11; + strncpy (thms, fdate+11, nbc); + return (thms); +} + + +/* TSD2DT-- convert seconds since start of day to hh.mmssss */ + +double +tsd2dt (tsec) + +double tsec; /* Seconds since start of day */ +{ + double date, time; + + ts2dt (tsec, &date, &time); + return (time); +} + + + +/* DT2I-- convert vigesimal date and time to year month day hours min sec */ + +void +dt2i (date, time, iyr, imon, iday, ihr, imn, sec, ndsec) + +double date; /* Date as yyyy.mmdd (returned) + yyyy = calendar year (e.g. 1973) + mm = calendar month (e.g. 04 = april) + dd = calendar day (e.g. 15) */ +double time; /* Time as hh.mmssxxxx (returned) + *if time<0, it is time as -(fraction of a day) + hh = hour of day (0 .le. hh .le. 23) + nn = minutes (0 .le. nn .le. 59) + ss = seconds (0 .le. ss .le. 59) + xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ + +{ + double t,d; + + t = time; + if (date < 0.0) + d = -date; + else + d = date; + + /* Extract components of time */ + *ihr = dint (t + 0.000000001); + t = 100.0 * (t - (double) *ihr); + *imn = dint (t + 0.0000001); + *sec = 100.0 * (t - (double) *imn); + + /* Extract components of date */ + *iyr = dint (d + 0.00001); + d = 100.0 * (d - (double) *iyr); + if (date < 0.0) + *iyr = - *iyr; + *imon = dint (d + 0.001); + d = 100.0 * (d - (double) *imon); + *iday = dint (d + 0.1); + + /* Make sure date and time are legal */ + fixdate (iyr, imon, iday, ihr, imn, sec, ndsec); + + return; +} + + +/* FD2I-- convert from FITS standard date to year, mon, day, hours, min, sec */ + +void +fd2i (string, iyr, imon, iday, ihr, imn, sec, ndsec) + +char *string; /* FITS date string, which may be: + yyyy.ffff (fractional year) + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard FITS use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ + +{ + double tsec, fday, hr, mn; + int i; + char *sstr, *dstr, *tstr, *cstr, *nval, *fstr; + + /* Initialize all returned data to zero */ + *iyr = 0; + *imon = 0; + *iday = 0; + *ihr = 0; + *imn = 0; + *sec = 0.0; + + /* Return if no input string */ + if (string == NULL) + return; + + /* Check for various non-numeric characters */ + sstr = strchr (string,'/'); + dstr = strchr (string,'-'); + if (dstr == string) + dstr = strchr (string+1, '-'); + fstr = strchr (string, '.'); + tstr = strchr (string,'T'); + if (tstr == NULL) + tstr = strchr (string, 'Z'); + if (tstr == NULL) + tstr = strchr (string, 'S'); + if (fstr != NULL && tstr != NULL && fstr > tstr) + fstr = NULL; + cstr = strchr (string,':'); + + /* Original FITS date format: dd/mm/yy */ + if (sstr > string) { + *sstr = '\0'; + *iday = (int) atof (string); + if (*iday > 31) { + *iyr = *iday; + if (*iyr >= 0 && *iyr <= 49) + *iyr = *iyr + 2000; + else if (*iyr < 1000) + *iyr = *iyr + 1900; + *sstr = '/'; + nval = sstr + 1; + sstr = strchr (nval,'/'); + if (sstr > string) { + *sstr = '\0'; + *imon = (int) atof (nval); + *sstr = '/'; + nval = sstr + 1; + *iday = (int) atof (nval); + } + } + else { + *sstr = '/'; + nval = sstr + 1; + sstr = strchr (nval,'/'); + if (sstr == NULL) + sstr = strchr (nval,'-'); + if (sstr > string) { + *sstr = '\0'; + *imon = (int) atof (nval); + *sstr = '/'; + nval = sstr + 1; + *iyr = (int) atof (nval); + if (*iyr >= 0 && *iyr <= 49) + *iyr = *iyr + 2000; + else if (*iyr < 1000) + *iyr = *iyr + 1900; + } + } + tstr = strchr (string,'_'); + if (tstr == NULL) + return; + } + + /* New FITS date format: yyyy-mm-ddThh:mm:ss[.sss] */ + else if (dstr > string) { + *dstr = '\0'; + *iyr = (int) atof (string); + *dstr = '-'; + nval = dstr + 1; + dstr = strchr (nval,'-'); + *imon = 1; + *iday = 1; + + /* Decode year, month, and day */ + if (dstr > string) { + *dstr = '\0'; + *imon = (int) atof (nval); + *dstr = '-'; + nval = dstr + 1; + if (tstr > string) + *tstr = '\0'; + *iday = (int) atof (nval); + + /* If fraction of a day is present, turn it into a time */ + if (fstr != NULL) { + fday = atof (fstr); + hr = fday * 24.0; + *ihr = (int) hr; + mn = 60.0 * (hr - (double) *ihr); + *imn = (int) mn; + *sec = 60.0 * (mn - (double) *imn); + } + + if (tstr > string) + *tstr = 'T'; + } + + /* If date is > 31, it is really year in old format */ + if (*iday > 31) { + i = *iyr; + if (*iday < 100) + *iyr = *iday + 1900; + else + *iyr = *iday; + *iday = i; + } + } + + /* In rare cases, a FITS time is entered as an epoch */ + else if (tstr == NULL && cstr == NULL && isnum (string)) { + tsec = ep2ts (atof (string)); + ts2i (tsec,iyr,imon,iday,ihr,imn,sec, ndsec); + return; + } + + /* Extract time, if it is present */ + if (tstr > string || cstr > string) { + if (tstr > string) + nval = tstr + 1; + else + nval = string; + cstr = strchr (nval,':'); + if (cstr > string) { + *cstr = '\0'; + *ihr = (int) atof (nval); + *cstr = ':'; + nval = cstr + 1; + cstr = strchr (nval,':'); + if (cstr > string) { + *cstr = '\0'; + *imn = (int) atof (nval); + *cstr = ':'; + nval = cstr + 1; + *sec = atof (nval); + } + else + *imn = (int) atof (nval); + } + else + *ihr = (int) atof (nval); + } + else + ndsec = -1; + + /* Make sure date and time are legal */ + fixdate (iyr, imon, iday, ihr, imn, sec, ndsec); + + return; +} + + +/* TS2I-- convert sec since 1950.0 to year month day hours minutes seconds */ + +void +ts2i (tsec,iyr,imon,iday,ihr,imn,sec, ndsec) + +double tsec; /* seconds since 1/1/1950 0:00 */ +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ + +{ + double t,days, ts, dts; + int nc,nc4,nly,ny,m,im; + + /* Round seconds to 0 - 4 decimal places */ + ts = tsec + 61530883200.0; + if (ts < 0.0) + dts = -0.5; + else + dts = 0.5; + if (ndsec < 1) + t = dint (ts + dts) * 10000.0; + else if (ndsec < 2) + t = dint (ts * 10.0 + dts) * 1000.0; + else if (ndsec < 3) + t = dint (ts * 100.0 + dts) * 100.0; + else if (ndsec < 4) + t = dint (ts * 1000.0 + dts) * 10.0; + else + t = dint (ts * 10000.0 + dts); + ts = t / 10000.0; + + /* Time of day (hours, minutes, seconds */ + *ihr = (int) (dmod (ts/3600.0, 24.0)); + *imn = (int) (dmod (ts/60.0, 60.0)); + *sec = dmod (ts, 60.0); + + /* Number of days since 0 hr 0/0/0000 */ + days = dint ((t / 864000000.0) + 0.000001); + + /* Number of leap centuries (400 years) */ + nc4 = (int) ((days / 146097.0) + 0.00001); + + /* Number of centuries since last /400 */ + days = days - (146097.0 * (double) (nc4)); + nc = (int) ((days / 36524.0) + 0.000001); + if (nc > 3) nc = 3; + + /* Number of leap years since last century */ + days = days - (36524.0 * nc); + nly = (int) ((days / 1461.0) + 0.0000000001); + + /* Number of years since last leap year */ + days = days - (1461.0 * (double) nly); + ny = (int) ((days / 365.0) + 0.00000001); + if (ny > 3) ny = 3; + + /* Day of month */ + days = days - (365.0 * (double) ny); + if (days < 0) { + m = 0; + *iday = 29; + } + else { + *iday = (int) (days + 0.00000001) + 1; + for (m = 1; m <= 12; m++) { + im = (m + ((m - 1) / 5)) % 2; + /* fprintf (stderr,"%d %d %d %d\n", m, im, *iday, nc); */ + if (*iday-1 < im+30) break; + *iday = *iday - im - 30; + } + } + + /* Month */ + *imon = ((m+1) % 12) + 1; + + /* Year */ + *iyr = nc4*400 + nc*100 + nly*4 + ny + m/11; + + /* Make sure date and time are legal */ + fixdate (iyr, imon, iday, ihr, imn, sec, ndsec); + + return; +} + + +/* UT2DOY-- Current Universal Time as year, day of year */ + +void +ut2doy (year, doy) + +int *year; /* Year (returned) */ +double *doy; /* Day of year (returned) */ +{ + double date, time; + ut2dt (&date, &time); + dt2doy (date, time, year, doy); + return; +} + + +/* UT2DT-- Current Universal Time as date (yyyy.mmdd) and time (hh.mmsss) */ + +void +ut2dt(date, time) + +double *date; /* Date as yyyy.mmdd (returned) */ +double *time; /* Time as hh.mmssxxxx (returned) */ +{ + time_t tsec; + struct timeval tp; + struct timezone tzp; + struct tm *ts; + + gettimeofday (&tp,&tzp); + + tsec = tp.tv_sec; + ts = gmtime (&tsec); + + if (ts->tm_year < 1000) + *date = (double) (ts->tm_year + 1900); + else + *date = (double) ts->tm_year; + *date = *date + (0.01 * (double) (ts->tm_mon + 1)); + *date = *date + (0.0001 * (double) ts->tm_mday); + *time = (double) ts->tm_hour; + *time = *time + (0.01 * (double) ts->tm_min); + *time = *time + (0.0001 * (double) ts->tm_sec); + + return; +} + + +/* UT2EP-- Return current Universal Time as fractional year */ + +double +ut2ep() +{ + return (jd2ep (ut2jd())); +} + + +/* UT2EPB-- Return current Universal Time as Besselian epoch */ + +double +ut2epb() +{ + return (jd2epb (ut2jd())); +} + + +/* UT2EPJ-- Return current Universal Time as Julian epoch */ + +double +ut2epj() +{ + return (jd2epj (ut2jd())); +} + + +/* UT2FD-- Return current Universal Time as FITS ISO date string */ + +char * +ut2fd() +{ + int year, month, day, hour, minute, second; + time_t tsec; + struct timeval tp; + struct timezone tzp; + struct tm *ts; + char *isotime; + + gettimeofday (&tp,&tzp); + tsec = tp.tv_sec; + ts = gmtime (&tsec); + + year = ts->tm_year; + if (year < 1000) + year = year + 1900; + month = ts->tm_mon + 1; + day = ts->tm_mday; + hour = ts->tm_hour; + minute = ts->tm_min; + second = ts->tm_sec; + + isotime = (char *) calloc (32, sizeof (char)); + sprintf (isotime, "%04d-%02d-%02dT%02d:%02d:%02d", + year, month, day, hour, minute, second); + return (isotime); +} + + +/* UT2JD-- Return current Universal Time as Julian Date */ + +double +ut2jd() +{ + return (fd2jd (ut2fd())); +} + + +/* UT2MJD-- convert current UT to Modified Julian Date */ + +double +ut2mjd () + +{ + return (ut2jd() - 2400000.5); +} + +/* UT2TS-- current Universal Time as IRAF seconds since 1950-01-01T00:00 */ + +double +ut2ts() +{ + double tsec; + char *datestring; + datestring = ut2fd(); + tsec = fd2ts (datestring); + free (datestring); + return (tsec); +} + + +/* UT2TSI-- current Universal Time as IRAF seconds since 1980-01-01T00:00 */ + +int +ut2tsi() +{ + return ((int)(ut2ts() - 946684800.0)); +} + + +/* UT2TSU-- current Universal Time as IRAF seconds since 1970-01-01T00:00 */ + +time_t +ut2tsu() +{ + return ((time_t)(ut2ts () - 631152000.0)); +} + + +/* FD2GST-- convert from FITS date to Greenwich Sidereal Time */ + +char * +fd2gst (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double dj, gsec, date, time; + + dj = fd2jd (string); + gsec = jd2gst (dj); + ts2dt (gsec, &date, &time); + date = 0.0; + return (dt2fd (date, time)); +} + + +/* DT2GST-- convert from UT as yyyy.mmdd hh.mmssss to Greenwich Sidereal Time*/ + +void +dt2gst (date, time) +double *date; /* Date as yyyy.mmdd */ +double *time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + double dj, gsec; + + dj = dt2ts (*date, *time); + gsec = jd2gst (dj); + ts2dt (gsec, date, time); + *date = 0.0; + return; +} + + +/* JD2LST - Local Sidereal Time in seconds from Julian Date */ + +double +jd2lst (dj) + +double dj; /* Julian Date */ +{ + double gst, lst; + + /* Compute Greenwich Sidereal Time at this epoch */ + gst = jd2gst (dj); + + /* Subtract longitude (degrees to seconds of time) */ + lst = gst - (240.0 * longitude); + if (lst < 0.0) + lst = lst + 86400.0; + else if (lst > 86400.0) + lst = lst - 86400.0; + return (lst); +} + + +/* FD2LST - Local Sidereal Time as hh:mm:ss.ss + from Universal Time as FITS ISO date */ + +char * +fd2lst (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) */ +{ + double dj, date, time, lst; + + dj = fd2jd (string); + lst = jd2lst (dj); + ts2dt (lst, &date, &time); + date = 0.0; + return (dt2fd (date, time)); +} + + +/* DT2LST - Local Sidereal Time as hh.mmssss + from Universal Time as yyyy.mmdd hh.mmssss */ + +void +dt2lst (date, time) + +double *date; /* Date as yyyy.mmdd */ +double *time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + double dj, lst, date0; + + dj = dt2jd (*date, *time); + lst = jd2lst (dj); + date0 = 0.0; + ts2dt (lst, &date0, time); + return; +} + + +/* TS2LST - Local Sidereal Time in seconds of day + * from Universal Time in seconds since 1951-01-01T0:00:00 + */ + +double +ts2lst (tsec) + +double tsec; /* time since 1950.0 in UT seconds */ +{ + double gst; /* Greenwich Sidereal Time in seconds since 0:00 */ + double lst; /* Local Sidereal Time in seconds since 0:00 */ + double gsec, date; + + /* Greenwich Sidereal Time */ + gsec = ts2gst (tsec); + date = 0.0; + ts2dt (gsec, &date, &gst); + + lst = gst - (longitude / 15.0); + if (lst < 0.0) + lst = lst + 86400.0; + else if (lst > 86400.0) + lst = lst - 86400.0; + return (lst); +} + + +/* LST2FD - calculate current UT given Local Sidereal Time + * plus date in FITS ISO format (yyyy-mm-dd) + * Return UT date and time in FITS ISO format + */ + +char * +lst2fd (string) + +char *string; /* UT Date, LST as yyyy-mm-ddShh:mm:ss.ss */ +{ + double sdj, dj; + + sdj = fd2jd (string); + + dj = lst2jd (sdj); + + return (jd2fd (dj)); +} + + +/* LST2JD - calculate current Julian Date given Local Sidereal Time + * plus current Julian Date (0.5 at 0:00 UT) + * Return UT date and time as Julian Date + */ + +double +lst2jd (sdj) + +double sdj; /* Julian Date of desired day at 0:00 UT + sidereal time */ +{ + double gst; /* Greenwich Sidereal Time in seconds since 0:00 */ + double lsd; /* Local Sidereal Time in seconds since 0:00 */ + double gst0, tsd, dj1, dj0, eqnx; + int idj; + + /* Julian date at 0:00 UT */ + idj = (int) sdj; + dj0 = (double) idj + 0.5; + if (dj0 > sdj) dj0 = dj0 - 1.0; + + /* Greenwich Sidereal Time at 0:00 UT in seconds */ + gst0 = jd2gst (dj0); + + /* Sidereal seconds since 0:00 */ + lsd = (sdj - dj0) * 86400.0; + + /* Remove longitude for current Greenwich Sidereal Time in seconds */ + /* (convert longitude from degrees to seconds of time) */ + gst = lsd + (longitude * 240.0); + + /* Time since 0:00 UT */ + tsd = (gst - gst0) / 1.0027379093; + + /* Julian Date (UT) */ + dj1 = dj0 + (tsd / 86400.0); + + /* Equation of the equinoxes converted to UT seconds */ + eqnx = eqeqnx (dj1) / 1.002739093; + + /* Remove equation of equinoxes */ + dj1 = dj1 - (eqnx / 86400.0); + if (dj1 < dj0) + dj1 = dj1 + 1.0; + + return (dj1); +} + + +/* MST2FD - calculate current UT given Greenwich Mean Sidereal Time + * plus date in FITS ISO format (yyyy-mm-ddShh:mm:ss.ss) + * Return UT date and time in FITS ISO format + */ + +char * +mst2fd (string) + +char *string; /* UT Date, MST as yyyy-mm-ddShh:mm:ss.ss */ +{ + double sdj, dj; + + sdj = fd2jd (string); + + dj = mst2jd (sdj); + + return (jd2fd (dj)); +} + + +/* MST2JD - calculate current UT given Greenwich Mean Sidereal Time + * plus date in Julian Date (0:00 UT + Mean Sidereal Time) + * Return UT date and time as Julian Date + */ + +double +mst2jd (sdj) + +double sdj; /* UT Date, MST as Julian Date */ +{ + double tsd, djd, st0, dj0, dj; + + dj0 = (double) ((int) sdj) + 0.5; + + /* Greenwich Mean Sidereal Time at 0:00 UT in seconds */ + st0 = jd2mst (dj0); + + /* Mean Sidereal Time in seconds */ + tsd = (sdj - dj0) * 86400.0; + if (tsd < 0.0) + tsd = tsd + 86400.0; + + /* Convert to fraction of a day since 0:00 UT */ + djd = ((tsd - st0) / 1.0027379093) / 86400.0; + + /* Julian Date */ + dj = dj0 + djd; + if (dj < dj0) + dj = dj + (1.0 / 1.0027379093); + + return (dj); +} + + + +/* GST2FD - calculate current UT given Greenwich Sidereal Time + * plus date in FITS ISO format (yyyy-mm-ddShh:mm:ss.ss) + * Return UT date and time in FITS ISO format + */ + +char * +gst2fd (string) + +char *string; /* UT Date, GST as yyyy-mm-ddShh:mm:ss.ss */ +{ + double sdj, dj; + + sdj = fd2jd (string); + + dj = gst2jd (sdj); + + return (jd2fd (dj)); +} + + +/* GST2JD - calculate current UT given Greenwich Sidereal Time + * plus date as Julian Date (JD at 0:00 UT + sidereal time) + * Return UT date and time as Julian Date + */ + +double +gst2jd (sdj) + +double sdj; /* UT Date, GST as Julian Date */ +{ + double dj, tsd, djd, st0, dj0, eqnx; + + dj0 = (double) ((int) sdj) + 0.5; + + /* Greenwich Mean Sidereal Time at 0:00 UT in seconds */ + st0 = jd2mst (dj0); + + /* Mean Sidereal Time in seconds */ + tsd = (sdj - dj0) * 86400.0; + if (tsd < 0.0) + tsd = tsd + 86400.0; + + /* Convert to fraction of a day since 0:00 UT */ + djd = ((tsd - st0) / 1.0027379093) / 86400.0; + + /* Julian Date */ + dj = dj0 + djd; + + /* Equation of the equinoxes (converted to UT seconds) */ + eqnx = eqeqnx (dj) / 1.002737909; + + dj = dj - eqnx / 86400.0; + if (dj < dj0) + dj = dj + 1.0; + + return (dj); +} + + +/* LST2DT - calculate current UT given Local Sidereal Time as hh.mmsss + * plus date as yyyy.mmdd + * Return UT time as hh.mmssss + */ + +double +lst2dt (date0, time0) + +double date0; /* UT date as yyyy.mmdd */ +double time0; /* LST as hh.mmssss */ +{ + double gst; /* Greenwich Sidereal Time in seconds since 0:00 */ + double lst; /* Local Sidereal Time in seconds since 0:00 */ + double date1; /* UT date as yyyy.mmdd */ + double time1; /* UT as hh.mmssss */ + double tsec0, gst0, tsd, tsec; + + /* Greenwich Sidereal Time at 0:00 UT */ + tsec0 = dt2ts (date0, 0.0); + gst0 = ts2gst (tsec0); + + /* Current Greenwich Sidereal Time in seconds */ + /* (convert longitude from degrees to seconds of time) */ + lst = dt2ts (0.0, time0); + gst = lst + (longitude * 240.0); + + /* Time since 0:00 UT */ + tsd = (gst - gst0) / 1.0027379093; + + /* UT date and time */ + tsec = tsec0 + tsd; + ts2dt (tsec, &date1, &time1); + + return (time1); +} + + +/* TS2GST - calculate Greenwich Sidereal Time given Universal Time + * in seconds since 1951-01-01T0:00:00 + * Return sidereal time of day in seconds + */ + +double +ts2gst (tsec) + +double tsec; /* time since 1950.0 in UT seconds */ +{ + double gst; /* Greenwich Sidereal Time in seconds since 0:00 */ + double tsd, eqnx, dj; + int its; + + /* Elapsed time as of 0:00 UT */ + if (tsec >= 0.0) { + its = (int) (tsec + 0.5); + tsd = (double) (its % 86400); + } + else { + its = (int) (-tsec + 0.5); + tsd = (double) (86400 - (its % 86400)); + } + + /* Mean sidereal time */ + gst = ts2mst (tsec); + + /* Equation of the equinoxes */ + dj = ts2jd (tsec); + eqnx = eqeqnx (dj); + + /* Apparent sidereal time at 0:00 ut */ + gst = gst + eqnx; + + /* Current sidereal time */ + gst = gst + (tsd * 1.0027379093); + gst = dmod (gst,86400.0); + + return (gst); +} + + +/* FD2MST-- convert from FITS date Mean Sidereal Time */ + +char * +fd2mst (string) + +char *string; /* FITS date string, which may be: + fractional year + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ +{ + double gsec, date, time, dj; + + dj = fd2jd (string); + gsec = jd2mst (dj); + ts2dt (gsec, &date, &time); + date = 0.0; + return (dt2fd (date, time)); +} + + +/* DT2MST-- convert from UT as yyyy.mmdd hh.mmssss to Mean Sidereal Time + in the same format */ + +void +dt2mst (date, time) +double *date; /* Date as yyyy.mmdd */ +double *time; /* Time as hh.mmssxxxx + *if time<0, it is time as -(fraction of a day) */ +{ + double date0, gsec, dj; + date0 = *date; + dj = dt2jd (*date, *time); + gsec = jd2mst (dj); + ts2dt (gsec, date, time); + *date = date0; + return; +} + + +/* TS2MST - calculate Greenwich Mean Sidereal Time given Universal Time + * in seconds since 1951-01-01T0:00:00 + */ + +double +ts2mst (tsec) + +double tsec; /* time since 1950.0 in UT seconds */ +{ + double dj; + + dj = ts2jd (tsec); + return (jd2mst (dj)); +} + + +/* JD2MST - Julian Date to Greenwich Mean Sidereal Time using IAU 2000 + * Return sideral time in seconds of time + * (from USNO NOVAS package + * http://aa.usno.navy.mil/software/novas/novas_info.html + */ + +double +jd2mst2 (dj) + +double dj; /* Julian Date */ +{ + double dt, t, t2, t3, mst, st; + + dt = dj - 2451545.0; + t = dt / 36525.0; + t2 = t * t; + t3 = t2 * t; + + /* Compute Greenwich Mean Sidereal Time in seconds */ + st = (8640184.812866 * t) + (3155760000.0 * t) - (0.0000062 * t3) + + (0.093104 * t2) + 67310.54841; + + mst = dmod (st, 86400.0); + if (mst < 0.0) + mst = mst + 86400.0; + return (mst); +} + + +/* MJD2MST - Modified Julian Date to Greenwich Mean Sidereal Time using IAU 2000 + * Return sideral time in seconds of time + * (from USNO NOVAS package + * http://aa.usno.navy.mil/software/novas/novas_info.html + */ + +double +mjd2mst (dj) + +double dj; /* Modified Julian Date */ +{ + double dt, t, t2, t3, mst, st; + + dt = dj - 51544.5; + t = dt / 36525.0; + t2 = t * t; + t3 = t2 * t; + + /* Compute Greenwich Mean Sidereal Time in seconds */ + st = (8640184.812866 * t) + (3155760000.0 * t) - (0.0000062 * t3) + + (0.093104 * t2) + 67310.54841; + + mst = dmod (st, 86400.0); + if (mst < 0.0) + mst = mst + 86400.0; + return (mst); +} + + +/* JD2GST - Julian Date to Greenwich Sideral Time + * Return sideral time in seconds of time + * (Jean Meeus, Astronomical Algorithms, Willmann-Bell, 1991, pp 83-84) + */ + +double +jd2gst (dj) + +double dj; /* Julian Date */ +{ + double dj0, gmt, gst, tsd, eqnx, ssd, l0; + double ts2ss = 1.00273790935; + int ijd; + + /* Julian date at 0:00 UT */ + ijd = (int) dj; + dj0 = (double) ijd + 0.5; + if (dj0 > dj) dj0 = dj0 - 1.0; + + /* Greenwich mean sidereal time at 0:00 UT in seconds */ + l0 = longitude; + longitude = 0.0; + gmt = jd2mst (dj0); + longitude = l0; + + /* Equation of the equinoxes */ + eqnx = eqeqnx (dj); + + /* Apparent sidereal time at 0:00 ut */ + gst = gmt + eqnx; + + /* UT seconds since 0:00 */ + tsd = (dj - dj0) * 86400.0; + ssd = tsd * ts2ss; + + /* Current sidereal time */ + gst = gst + ssd; + gst = dmod (gst, 86400.0); + + return (gst); +} + + +/* EQEQNX - Compute equation of the equinoxes for apparent sidereal time */ + +double +eqeqnx (dj) + +double dj; /* Julian Date */ + +{ + double dt, edj, dpsi, deps, obl, eqnx; + double rad2tsec = 13750.98708; + + /* Convert UT to Ephemeris Time (TDB or TT)*/ + dt = utdt (dj); + edj = dj + dt / 86400.0; + + /* Nutation and obliquity */ + compnut (edj, &dpsi, &deps, &obl); + + /* Correct obliquity for nutation */ + obl = obl + deps; + + /* Equation of the equinoxes in seconds */ + eqnx = (dpsi * cos (obl)) * rad2tsec; + + return (eqnx); +} + + + +/* JD2MST - Julian Date to Mean Sideral Time + * Return sideral time in seconds of time + * (Jean Meeus, Astronomical Algorithms, Willmann-Bell, 1991, pp 83-84) + */ + +double +jd2mst (dj) + +double dj; /* Julian Date */ +{ + double dt, t, mst; + + dt = dj - 2451545.0; + t = dt / 36525.0; + + /* Compute Greenwich mean sidereal time in degrees (Meeus, page 84) */ + mst = 280.46061837 + (360.98564736629 * dt) + (0.000387933 * t * t) - + (t * t * t / 38710000.0); + + /* Keep degrees between 0 and 360 */ + while (mst > 360.0) + mst = mst - 360.0; + while (mst < 0.0) + mst = mst + 360.0; + + /* Convert to time in seconds (3600 / 15) */ + mst = mst * 240.0; + + /* Subtract longitude (degrees to seconds of time) */ + mst = mst - (240.0 * longitude); + if (mst < 0.0) + mst = mst + 86400.0; + else if (mst > 86400.0) + mst = mst - 86400.0; + + return (mst); +} + + +/* COMPNUT - Compute nutation using the IAU 2000b model */ +/* Translated from Pat Wallace's Fortran subroutine iau_nut00b (June 26 2007) + into C by Jessica Mink on September 5, 2008 */ + +#define NLS 77 /* number of terms in the luni-solar nutation model */ + +void +compnut (dj, dpsi, deps, eps0) + +double dj; /* Julian Date */ +double *dpsi; /* Nutation in longitude in radians (returned) */ +double *deps; /* Nutation in obliquity in radians (returned) */ +double *eps0; /* Mean obliquity in radians (returned) */ + +/* This routine is translated from the International Astronomical Union's + * Fortran SOFA (Standards Of Fundamental Astronomy) software collection. + * + * notes: + * + * 1) the nutation components in longitude and obliquity are in radians + * and with respect to the equinox and ecliptic of date. the + * obliquity at j2000 is assumed to be the lieske et al. (1977) value + * of 84381.448 arcsec. (the errors that result from using this + * routine with the iau 2006 value of 84381.406 arcsec can be + * neglected.) + * + * the nutation model consists only of luni-solar terms, but includes + * also a fixed offset which compensates for certain long-period + * planetary terms (note 7). + * + * 2) this routine is an implementation of the iau 2000b abridged + * nutation model formally adopted by the iau general assembly in + * 2000. the routine computes the mhb_2000_short luni-solar nutation + * series (luzum 2001), but without the associated corrections for + * the precession rate adjustments and the offset between the gcrs + * and j2000 mean poles. + * + * 3) the full IAU 2000a (mhb2000) nutation model contains nearly 1400 + * terms. the IAU 2000b model (mccarthy & luzum 2003) contains only + * 77 terms, plus additional simplifications, yet still delivers + * results of 1 mas accuracy at present epochs. this combination of + * accuracy and size makes the IAU 2000b abridged nutation model + * suitable for most practical applications. + * + * the routine delivers a pole accurate to 1 mas from 1900 to 2100 + * (usually better than 1 mas, very occasionally just outside 1 mas). + * the full IAU 2000a model, which is implemented in the routine + * iau_nut00a (q.v.), delivers considerably greater accuracy at + * current epochs; however, to realize this improved accuracy, + * corrections for the essentially unpredictable free-core-nutation + * (fcn) must also be included. + * + * 4) the present routine provides classical nutation. the + * mhb_2000_short algorithm, from which it is adapted, deals also + * with (i) the offsets between the gcrs and mean poles and (ii) the + * adjustments in longitude and obliquity due to the changed + * precession rates. these additional functions, namely frame bias + * and precession adjustments, are supported by the sofa routines + * iau_bi00 and iau_pr00. + * + * 6) the mhb_2000_short algorithm also provides "total" nutations, + * comprising the arithmetic sum of the frame bias, precession + * adjustments, and nutation (luni-solar + planetary). these total + * nutations can be used in combination with an existing IAU 1976 + * precession implementation, such as iau_pmat76, to deliver gcrs-to- + * true predictions of mas accuracy at current epochs. however, for + * symmetry with the iau_nut00a routine (q.v. for the reasons), the + * sofa routines do not generate the "total nutations" directly. + * should they be required, they could of course easily be generated + * by calling iau_bi00, iau_pr00 and the present routine and adding + * the results. + * + * 7) the IAU 2000b model includes "planetary bias" terms that are fixed + * in size but compensate for long-period nutations. the amplitudes + * quoted in mccarthy & luzum (2003), namely dpsi = -1.5835 mas and + * depsilon = +1.6339 mas, are optimized for the "total nutations" + * method described in note 6. the luzum (2001) values used in this + * sofa implementation, namely -0.135 mas and +0.388 mas, are + * optimized for the "rigorous" method, where frame bias, precession + * and nutation are applied separately and in that order. during the + * interval 1995-2050, the sofa implementation delivers a maximum + * error of 1.001 mas (not including fcn). + * + * References from original Fortran subroutines: + * + * Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351 + * + * Lieske, J.H., Lederle, T., Fricke, W., Morando, B., "Expressions + * for the precession quantities based upon the IAU 1976 system of + * astronomical constants", Astron.Astrophys. 58, 1-2, 1-16. (1977) + * + * Luzum, B., private communication, 2001 (Fortran code + * mhb_2000_short) + * + * McCarthy, D.D. & Luzum, B.J., "An abridged model of the + * precession-nutation of the celestial pole", Cel.Mech.Dyn.Astron. + * 85, 37-49 (2003) + * + * Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M., + * Francou, G., Laskar, J., Astron.Astrophys. 282, 663-683 (1994) + * + */ + +{ + double as2r = 0.000004848136811095359935899141; /* arcseconds to radians */ + + double dmas2r = as2r / 1000.0; /* milliarcseconds to radians */ + + double as2pi = 1296000.0; /* arc seconds in a full circle */ + + double d2pi = 6.283185307179586476925287; /* 2pi */ + + double u2r = as2r / 10000000.0; /* units of 0.1 microarcsecond to radians */ + + double dj0 = 2451545.0; /* reference epoch (j2000), jd */ + + double djc = 36525.0; /* Days per julian century */ + + /* Miscellaneous */ + double t, el, elp, f, d, om, arg, dp, de, sarg, carg; + double dpsils, depsls, dpsipl, depspl; + int i, j; + + int nls = NLS; /* number of terms in the luni-solar nutation model */ + + /* Fixed offset in lieu of planetary terms (radians) */ + double dpplan = - 0.135 * dmas2r; + double deplan = + 0.388 * dmas2r; + +/* Tables of argument and term coefficients */ + + /* Coefficients for fundamental arguments */ + /* Luni-solar argument multipliers: */ + /* l l' f d om */ +static int nals[5*NLS]= + {0, 0, 0, 0, 1, + 0, 0, 2, -2, 2, + 0, 0, 2, 0, 2, + 0, 0, 0, 0, 2, + 0, 1, 0, 0, 0, + 0, 1, 2, -2, 2, + 1, 0, 0, 0, 0, + 0, 0, 2, 0, 1, + 1, 0, 2, 0, 2, + 0, -1, 2, -2, 2, + 0, 0, 2, -2, 1, + -1, 0, 2, 0, 2, + -1, 0, 0, 2, 0, + 1, 0, 0, 0, 1, + -1, 0, 0, 0, 1, + -1, 0, 2, 2, 2, + 1, 0, 2, 0, 1, + -2, 0, 2, 0, 1, + 0, 0, 0, 2, 0, + 0, 0, 2, 2, 2, + 0, -2, 2, -2, 2, + -2, 0, 0, 2, 0, + 2, 0, 2, 0, 2, + 1, 0, 2, -2, 2, + -1, 0, 2, 0, 1, + 2, 0, 0, 0, 0, + 0, 0, 2, 0, 0, + 0, 1, 0, 0, 1, + -1, 0, 0, 2, 1, + 0, 2, 2, -2, 2, + 0, 0, -2, 2, 0, + 1, 0, 0, -2, 1, + 0, -1, 0, 0, 1, + -1, 0, 2, 2, 1, + 0, 2, 0, 0, 0, + 1, 0, 2, 2, 2, + -2, 0, 2, 0, 0, + 0, 1, 2, 0, 2, + 0, 0, 2, 2, 1, + 0, -1, 2, 0, 2, + 0, 0, 0, 2, 1, + 1, 0, 2, -2, 1, + 2, 0, 2, -2, 2, + -2, 0, 0, 2, 1, + 2, 0, 2, 0, 1, + 0, -1, 2, -2, 1, + 0, 0, 0, -2, 1, + -1, -1, 0, 2, 0, + 2, 0, 0, -2, 1, + 1, 0, 0, 2, 0, + 0, 1, 2, -2, 1, + 1, -1, 0, 0, 0, + -2, 0, 2, 0, 2, + 3, 0, 2, 0, 2, + 0, -1, 0, 2, 0, + 1, -1, 2, 0, 2, + 0, 0, 0, 1, 0, + -1, -1, 2, 2, 2, + -1, 0, 2, 0, 0, + 0, -1, 2, 2, 2, + -2, 0, 0, 0, 1, + 1, 1, 2, 0, 2, + 2, 0, 0, 0, 1, + -1, 1, 0, 1, 0, + 1, 1, 0, 0, 0, + 1, 0, 2, 0, 0, + -1, 0, 2, -2, 1, + 1, 0, 0, 0, 2, + -1, 0, 0, 1, 0, + 0, 0, 2, 1, 2, + -1, 0, 2, 4, 2, + -1, 1, 0, 1, 1, + 0, -2, 2, -2, 1, + 1, 0, 2, 2, 1, + -2, 0, 2, 2, 2, + -1, 0, 0, 0, 2, + 1, 1, 2, -2, 2}; + + /* Luni-solar nutation coefficients, in 1e-7 arcsec */ + /* longitude (sin, t*sin, cos), obliquity (cos, t*cos, sin) */ +static double cls[6*NLS]= + {-172064161.0, -174666.0, 33386.0, 92052331.0, 9086.0, 15377.0, + -13170906.0, -1675.0, -13696.0, 5730336.0, -3015.0, -4587.0, + -2276413.0, -234.0, 2796.0, 978459.0, -485.0, 1374.0, + 2074554.0, 207.0, -698.0, -897492.0, 470.0, -291.0, + 1475877.0, -3633.0, 11817.0, 73871.0, -184.0, -1924.0, + -516821.0, 1226.0, -524.0, 224386.0, -677.0, -174.0, + 711159.0, 73.0, -872.0, -6750.0, 0.0, 358.0, + -387298.0, -367.0, 380.0, 200728.0, 18.0, 318.0, + -301461.0, -36.0, 816.0, 129025.0, -63.0, 367.0, + 215829.0, -494.0, 111.0, -95929.0, 299.0, 132.0, + 128227.0, 137.0, 181.0, -68982.0, -9.0, 39.0, + 123457.0, 11.0, 19.0, -53311.0, 32.0, -4.0, + 156994.0, 10.0, -168.0, -1235.0, 0.0, 82.0, + 63110.0, 63.0, 27.0, -33228.0, 0.0, -9.0, + -57976.0, -63.0, -189.0, 31429.0, 0.0, -75.0, + -59641.0, -11.0, 149.0, 25543.0, -11.0, 66.0, + -51613.0, -42.0, 129.0, 26366.0, 0.0, 78.0, + 45893.0, 50.0, 31.0, -24236.0, -10.0, 20.0, + 63384.0, 11.0, -150.0, -1220.0, 0.0, 29.0, + -38571.0, -1.0, 158.0, 16452.0, -11.0, 68.0, + 32481.0, 0.0, 0.0, -13870.0, 0.0, 0.0, + -47722.0, 0.0, -18.0, 477.0, 0.0, -25.0, + -31046.0, -1.0, 131.0, 13238.0, -11.0, 59.0, + 28593.0, 0.0, -1.0, -12338.0, 10.0, -3.0, + 20441.0, 21.0, 10.0, -10758.0, 0.0, -3.0, + 29243.0, 0.0, -74.0, -609.0, 0.0, 13.0, + 25887.0, 0.0, -66.0, -550.0, 0.0, 11.0, + -14053.0, -25.0, 79.0, 8551.0, -2.0, -45.0, + 15164.0, 10.0, 11.0, -8001.0, 0.0, -1.0, + -15794.0, 72.0, -16.0, 6850.0, -42.0, -5.0, + 21783.0, 0.0, 13.0, -167.0, 0.0, 13.0, + -12873.0, -10.0, -37.0, 6953.0, 0.0, -14.0, + -12654.0, 11.0, 63.0, 6415.0, 0.0, 26.0, + -10204.0, 0.0, 25.0, 5222.0, 0.0, 15.0, + 16707.0, -85.0, -10.0, 168.0, -1.0, 10.0, + -7691.0, 0.0, 44.0, 3268.0, 0.0, 19.0, + -11024.0, 0.0, -14.0, 104.0, 0.0, 2.0, + 7566.0, -21.0, -11.0, -3250.0, 0.0, -5.0, + -6637.0, -11.0, 25.0, 3353.0, 0.0, 14.0, + -7141.0, 21.0, 8.0, 3070.0, 0.0, 4.0, + -6302.0, -11.0, 2.0, 3272.0, 0.0, 4.0, + 5800.0, 10.0, 2.0, -3045.0, 0.0, -1.0, + 6443.0, 0.0, -7.0, -2768.0, 0.0, -4.0, + -5774.0, -11.0, -15.0, 3041.0, 0.0, -5.0, + -5350.0, 0.0, 21.0, 2695.0, 0.0, 12.0, + -4752.0, -11.0, -3.0, 2719.0, 0.0, -3.0, + -4940.0, -11.0, -21.0, 2720.0, 0.0, -9.0, + 7350.0, 0.0, -8.0, -51.0, 0.0, 4.0, + 4065.0, 0.0, 6.0, -2206.0, 0.0, 1.0, + 6579.0, 0.0, -24.0, -199.0, 0.0, 2.0, + 3579.0, 0.0, 5.0, -1900.0, 0.0, 1.0, + 4725.0, 0.0, -6.0, -41.0, 0.0, 3.0, + -3075.0, 0.0, -2.0, 1313.0, 0.0, -1.0, + -2904.0, 0.0, 15.0, 1233.0, 0.0, 7.0, + 4348.0, 0.0, -10.0, -81.0, 0.0, 2.0, + -2878.0, 0.0, 8.0, 1232.0, 0.0, 4.0, + -4230.0, 0.0, 5.0, -20.0, 0.0, -2.0, + -2819.0, 0.0, 7.0, 1207.0, 0.0, 3.0, + -4056.0, 0.0, 5.0, 40.0, 0.0, -2.0, + -2647.0, 0.0, 11.0, 1129.0, 0.0, 5.0, + -2294.0, 0.0, -10.0, 1266.0, 0.0, -4.0, + 2481.0, 0.0, -7.0, -1062.0, 0.0, -3.0, + 2179.0, 0.0, -2.0, -1129.0, 0.0, -2.0, + 3276.0, 0.0, 1.0, -9.0, 0.0, 0.0, + -3389.0, 0.0, 5.0, 35.0, 0.0, -2.0, + 3339.0, 0.0, -13.0, -107.0, 0.0, 1.0, + -1987.0, 0.0, -6.0, 1073.0, 0.0, -2.0, + -1981.0, 0.0, 0.0, 854.0, 0.0, 0.0, + 4026.0, 0.0, -353.0, -553.0, 0.0, -139.0, + 1660.0, 0.0, -5.0, -710.0, 0.0, -2.0, + -1521.0, 0.0, 9.0, 647.0, 0.0, 4.0, + 1314.0, 0.0, 0.0, -700.0, 0.0, 0.0, + -1283.0, 0.0, 0.0, 672.0, 0.0, 0.0, + -1331.0, 0.0, 8.0, 663.0, 0.0, 4.0, + 1383.0, 0.0, -2.0, -594.0, 0.0, -2.0, + 1405.0, 0.0, 4.0, -610.0, 0.0, 2.0, + 1290.0, 0.0, 0.0, -556.0, 0.0, 0.0}; + + /* Interval between fundamental epoch J2000.0 and given date (JC) */ + t = (dj - dj0) / djc; + +/* Luni-solar nutation */ + +/* Fundamental (delaunay) arguments from Simon et al. (1994) */ + + /* Mean anomaly of the moon */ + el = fmod (485868.249036 + (1717915923.2178 * t), as2pi) * as2r; + + /* Mean anomaly of the sun */ + elp = fmod (1287104.79305 + (129596581.0481 * t), as2pi) * as2r; + + /* Mean argument of the latitude of the moon */ + f = fmod (335779.526232 + (1739527262.8478 * t), as2pi) * as2r; + + /* Mean elongation of the moon from the sun */ + d = fmod (1072260.70369 + (1602961601.2090 * t), as2pi ) * as2r; + + /* Mean longitude of the ascending node of the moon */ + om = fmod (450160.398036 - (6962890.5431 * t), as2pi ) * as2r; + + /* Initialize the nutation values */ + dp = 0.0; + de = 0.0; + + /* Summation of luni-solar nutation series (in reverse order) */ + for (i = nls; i > 0; i=i-1) { + j = i - 1; + + /* Argument and functions */ + arg = fmod ( (double) (nals[5*j]) * el + + (double) (nals[1+5*j]) * elp + + (double) (nals[2+5*j]) * f + + (double) (nals[3+5*j]) * d + + (double) (nals[4+5*j]) * om, d2pi); + sarg = sin (arg); + carg = cos (arg); + + /* Terms */ + dp = dp + (cls[6*j] + cls[1+6*j] * t) * sarg + cls[2+6*j] * carg; + de = de + (cls[3+6*j] + cls[4+6*j] * t) * carg + cls[5+6*j] * sarg; + } + + /* Convert from 0.1 microarcsec units to radians */ + dpsils = dp * u2r; + depsls = de * u2r; + +/* In lieu of planetary nutation */ + + /* Fixed offset to correct for missing terms in truncated series */ + dpsipl = dpplan; + depspl = deplan; + +/* Results */ + + /* Add luni-solar and planetary components */ + *dpsi = dpsils + dpsipl; + *deps = depsls + depspl; + + /* Mean Obliquity in radians (IAU 2006, Hilton, et al.) */ + *eps0 = ( 84381.406 + + ( -46.836769 + + ( -0.0001831 + + ( 0.00200340 + + ( -0.000000576 + + ( -0.0000000434 ) * t ) * t ) * t ) * t ) * t ) * as2r; +} + + +/* ISDATE - Return 1 if string is an old or ISO FITS standard date */ + +int +isdate (string) + +char *string; /* Possible FITS date string, which may be: + dd/mm/yy (FITS standard before 2000) + dd-mm-yy (nonstandard FITS use before 2000) + yyyy-mm-dd (FITS standard after 1999) + yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */ + +{ + int iyr = 0; /* year (returned) */ + int imon = 0; /* month (returned) */ + int iday = 0; /* day (returned) */ + int i; + char *sstr, *dstr, *tstr, *nval; + + /* Translate string from ASCII to binary */ + if (string == NULL) + return (0); + + sstr = strchr (string,'/'); + dstr = strchr (string,'-'); + if (dstr == string) + dstr = strchr (string+1,'-'); + tstr = strchr (string,'T'); + + /* Original FITS date format: dd/mm/yy */ + if (sstr > string) { + *sstr = '\0'; + iday = (int) atof (string); + *sstr = '/'; + nval = sstr + 1; + sstr = strchr (nval,'/'); + if (sstr == NULL) + sstr = strchr (nval,'-'); + if (sstr > string) { + *sstr = '\0'; + imon = (int) atof (nval); + *sstr = '/'; + nval = sstr + 1; + iyr = (int) atof (nval); + if (iyr < 1000) + iyr = iyr + 1900; + } + if (imon > 0 && iday > 0) + return (1); + else + return (0); + } + + /* New FITS date format: yyyy-mm-ddThh:mm:ss[.sss] */ + else if (dstr > string) { + *dstr = '\0'; + iyr = (int) atof (string); + nval = dstr + 1; + *dstr = '-'; + dstr = strchr (nval,'-'); + imon = 0; + iday = 0; + + /* Decode year, month, and day */ + if (dstr > string) { + *dstr = '\0'; + imon = (int) atof (nval); + *dstr = '-'; + nval = dstr + 1; + if (tstr > string) + *tstr = '\0'; + iday = (int) atof (nval); + if (tstr > string) + *tstr = 'T'; + } + + /* If day is > 31, it is really year in old format */ + if (iday > 31) { + i = iyr; + if (iday < 100) + iyr = iday + 1900; + else + iyr = iday; + iday = i; + } + if (imon > 0 && iday > 0) + return (1); + else + return (0); + } + + /* If FITS date is entered as an epoch, return 0 anyway */ + else + return (0); +} + + +/* Round seconds and make sure date and time numbers are within limits */ + +static void +fixdate (iyr, imon, iday, ihr, imn, sec, ndsec) + +int *iyr; /* year (returned) */ +int *imon; /* month (returned) */ +int *iday; /* day (returned) */ +int *ihr; /* hours (returned) */ +int *imn; /* minutes (returned) */ +double *sec; /* seconds (returned) */ +int ndsec; /* Number of decimal places in seconds (0=int) */ +{ + double days; + + /* Round seconds to 0 - 4 decimal places (no rounding if <0, >4) */ + if (ndsec == 0) + *sec = dint (*sec + 0.5); + else if (ndsec < 2) + *sec = dint (*sec * 10.0 + 0.5) / 10.0; + else if (ndsec < 3) + *sec = dint (*sec * 100.0 + 0.5) / 100.0; + else if (ndsec < 4) + *sec = dint (*sec * 1000.0 + 0.5) / 1000.0; + else if (ndsec < 5) + *sec = dint (*sec * 10000.0 + 0.5) / 10000.0; + + /* Adjust minutes and hours */ + if (*sec > 60.0) { + *sec = *sec - 60.0; + *imn = *imn + 1; + } + if (*imn > 60) { + *imn = *imn - 60; + *ihr = *ihr + 1; + } + + /* Return if no date */ + if (*iyr == 0 && *imon == 0 && *iday == 0) + return; + + /* Adjust date */ + if (*ihr > 23) { + *ihr = *ihr - 24; + *iday = *iday + 1; + } + days = caldays (*iyr, *imon); + if (*iday > days) { + *iday = *iday - days; + *imon = *imon + 1; + } + if (*iday < 1) { + *imon = *imon - 1; + if (*imon < 1) { + *imon = *imon + 12; + *iyr = *iyr - 1; + } + days = caldays (*iyr, *imon); + *iday = *iday + days; + } + if (*imon < 1) { + *imon = *imon + 12; + *iyr = *iyr - 1; + days = caldays (*iyr, *imon); + if (*iday > days) { + *iday = *iday - days; + *imon = *imon + 1; + } + } + if (*imon > 12) { + *imon = *imon - 12; + *iyr = *iyr + 1; + } + return; +} + + +/* Calculate days in month 1-12 given year (Gregorian calendar only) */ + +static int +caldays (year, month) + +int year; /* 4-digit year */ +int month; /* Month (1=January, 2=February, etc.) */ +{ + if (month < 1) { + month = month + 12; + year = year + 1; + } + if (month > 12) { + month = month - 12; + year = year + 1; + } + switch (month) { + case 1: + return (31); + case 2: + if (year%400 == 0) + return (29); + else if (year%100 == 0) + return (28); + else if (year%4 == 0) + return (29); + else + return (28); + case 3: + return (31); + case 4: + return (30); + case 5: + return (31); + case 6: + return (30); + case 7: + return (31); + case 8: + return (31); + case 9: + return (30); + case 10: + return (31); + case 11: + return (30); + case 12: + return (31); + default: + return (0); + } +} + + +static double +dint (dnum) + +double dnum; +{ + double dn; + + if (dnum < 0.0) + dn = -floor (-dnum); + else + dn = floor (dnum); + return (dn); +} + + +static double +dmod (dnum, dm) + +double dnum, dm; +{ + double dnumx, dnumi, dnumf; + if (dnum < 0.0) + dnumx = -dnum; + else + dnumx = dnum; + dnumi = dint (dnumx / dm); + if (dnum < 0.0) + dnumf = dnum + (dnumi * dm); + else if (dnum > 0.0) + dnumf = dnum - (dnumi * dm); + else + dnumf = 0.0; + return (dnumf); +} + +/* Jul 1 1999 New file, based on iolib/jcon.f and iolib/vcon.f and hgetdate() + * Oct 21 1999 Fix declarations after lint + * Oct 27 1999 Fix bug to return epoch if fractional year input + * Dec 9 1999 Fix bug in ts2jd() found by Pete Ratzlaff (SAO) + * Dec 17 1999 Add all unimplemented conversions + * Dec 20 1999 Add isdate(); leave date, time strings unchanged in fd2i() + * Dec 20 1999 Make all fd2*() subroutines deal with time alone + * + * Jan 3 2000 In old FITS format, year 100 is assumed to be 2000 + * Jan 11 2000 Fix epoch to date conversion so .0 is 0:00, not 12:00 + * Jan 21 2000 Add separate Besselian and Julian epoch computations + * Jan 28 2000 Add Modified Julian Date conversions + * Mar 2 2000 Implement decimal places for FITS date string + * Mar 14 2000 Fix bug in dealing with 2000-02-29 in ts2i() + * Mar 22 2000 Add lt2* and ut2* to get current time as local and UT + * Mar 24 2000 Fix calloc() calls + * Mar 24 2000 Add tsi2* and tsu2* to convert IRAF and Unix seconds + * May 1 2000 In old FITS format, all years < 1000 get 1900 added to them + * Aug 1 2000 Make ep2jd and jd2ep consistently starting at 1/1 0:00 + * + * Jan 11 2001 Print all messages to stderr + * May 21 2001 Add day of year conversions + * May 25 2001 Allow fraction of day in FITS date instead of time + * + * Apr 8 2002 Change all long declaration to time_t + * May 13 2002 Fix bugs found by lint + * Jul 5 2002 Fix bug in fixdate() so fractional seconds come out + * Jul 8 2002 Fix rounding bug in t2i() + * Jul 8 2002 Try Fliegel and Van Flandern's algorithm for JD to UT date + * Jul 8 2002 If first character of string is -, check for other -'s in isdate + * Sep 10 2002 Add ET/TDT/TT conversion from UT subroutines + * Sep 10 2002 Add sidereal time conversions + * + * Jan 30 2003 Fix typo in ts2gst() + * Mar 7 2003 Add conversions for heliocentric julian dates + * May 20 2003 Declare nd in setdatedec() + * Jul 18 2003 Add code to parse Las Campanas dates + * + * Mar 24 2004 If ndec > 0, add UT to FITS date even if it is 0:00:00 + * + * Oct 14 2005 Add tsd2fd() and tsd2dt() + * + * May 3 2006 Drop declaration of unused variables + * Jun 20 2006 Initialized uninitialized variables + * Aug 2 2006 Add local sidereal time + * Sep 13 2006 Add more local sidereal time subroutines + * Oct 2 2006 Add UT to old FITS date conversions + * Oct 6 2006 Add eqeqnx() to compute equation of the equinoxes + * + * Jan 8 2007 Remove unused variables + * + * Sep 5 2008 Replace nutation with IAU 2006 model translated from SOFA + * Sep 9 2008 Add ang2hr(), ang2deg(), hr2ang(), deg2ang() + * Sep 10 2008 Add longitude to mean standard time (default = Greenwich) + * Oct 8 2008 Clean up sidereal time computations + * + * Sep 24 2009 Add end to comment "Coefficients for fundamental arguments" + * + * Jan 11 2012 Add TAI, TT, GPS time + * Oct 19 2012 Unused l0 dropped from jd2lst(); ts2ss from jd2mst() + */ |