1 /*** File libwcs/dateutil.c
2 *** May 2, 2017
3 *** By Jessica Mink, jmink@cfa.harvard.edu
4 *** Harvard-Smithsonian Center for Astrophysics
5 *** Copyright (C) 1999-2017
6 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
7
8 This library is free software; you can redistribute it and/or
9 modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2 of the License, or (at your option) any later version.
12
13 This library is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with this library; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22 Correspondence concerning WCSTools should be addressed as follows:
23 Internet email: jmink@cfa.harvard.edu
24 Postal address: Jessica Mink
25 Smithsonian Astrophysical Observatory
26 60 Garden St.
27 Cambridge, MA 02138 USA
28 */
29
30 /* Date and time conversion routines using the following conventions:
31 ang = Angle in fractional degrees
32 deg = Angle in degrees as dd:mm:ss.ss
33 doy = 2 floating point numbers: year and day, including fraction, of year
34 *** First day of year is 1, not zero.
35 dt = 2 floating point numbers: yyyy.mmdd, hh.mmssssss
36 ep = fractional year, often epoch of a position including proper motion
37 epb = Besselian epoch = 365.242198781-day years based on 1900.0
38 epj = Julian epoch = 365.25-day years based on 2000.0
39 fd = FITS date string which may be any of the following:
40 yyyy.ffff (fractional year)
41 dd/mm/yy (FITS standard before 2000)
42 dd-mm-yy (nonstandard FITS use before 2000)
43 yyyy-mm-dd (FITS standard after 1999)
44 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999)
45 hr = Sexigesimal hours as hh:mm:dd.ss
46 jd = Julian Date
47 lt = Local time
48 mjd = modified Julian Date = JD - 2400000.5
49 ofd = FITS date string (dd/mm/yy before 2000, else no return)
50 time = use fd2* with no date to convert time as hh:mm:ss.ss to sec, day, year
51 ts = UT seconds since 1950-01-01T00:00 (used for ephemeris computations)
52 tsi = local seconds since 1980-01-01T00:00 (used by IRAF as a time tag)
53 tsu = UT seconds since 1970-01-01T00:00 (used as Unix system time)
54 tsd = UT seconds of current day
55 ut = Universal Time (UTC)
56 et = Ephemeris Time (or TDB or TT) = TAI + 32.184 seconds
57 tai = International Atomic Time (Temps Atomique International) = ET - 32.184 seconds
58 gps = GPS time = TAI - 19 seconds
59 mst = Mean Greenwich Sidereal Time
60 gst = Greenwich Sidereal Time (includes nutation)
61 lst = Local Sidereal Time (includes nutation) (longitude must be set)
62 hjd = Heliocentric Julian Date
63 mhjd = modified Heliocentric Julian Date = HJD - 2400000.5
64
65 * ang2hr (angle)
66 * Convert angle in decimal floating point degrees to hours as hh:mm:ss.ss
67 * ang2deg (angle)
68 * Convert angle in decimal floating point degrees to degrees as dd:mm:ss.ss
69 * deg2ang (angle as dd:mm:ss.ss)
70 * Convert angle in degrees as dd:mm:ss.ss to decimal floating point degrees
71 * ang2hr (angle)
72 * Convert angle in hours as hh:mm:ss.ss to decimal floating point degrees
73 *
74 * doy2dt (year, doy, date, time)
75 * Convert year and day of year to date as yyyy.ddmm and time as hh.mmsss
76 * doy2ep, doy2epb, doy2epj (date, time)
77 * Convert year and day of year to fractional year
78 * doy2fd (year, doy)
79 * Convert year and day of year to FITS date string
80 * doy2mjd (year, doy)
81 * Convert year and day of year to modified Julian date
82 *
83 * dt2doy (date, time, year, doy)
84 * Convert date as yyyy.ddmm and time as hh.mmsss to year and day of year
85 * dt2ep, dt2epb, dt2epj (date, time)
86 * Convert date as yyyy.ddmm and time as hh.mmsss to fractional year
87 * dt2fd (date, time)
88 * Convert date as yyyy.ddmm and time as hh.mmsss to FITS date string
89 * dt2i (date,time,iyr,imon,iday,ihr,imn,sec, ndsec)
90 * Convert yyyy.mmdd hh.mmssss to year month day hours minutes seconds
91 * dt2jd (date,time)
92 * Convert date as yyyy.ddmm and time as hh.mmsss to Julian date
93 * dt2mjd (date,time)
94 * Convert date as yyyy.ddmm and time as hh.mmsss to modified Julian date
95 * dt2ts (date,time)
96 * Convert date (yyyy.ddmm) and time (hh.mmsss) to seconds since 1950-01-01
97 * dt2tsi (date,time)
98 * Convert date (yyyy.ddmm) and time (hh.mmsss) to seconds since 1980-01-01
99 * dt2tsu (date,time)
100 * Convert date (yyyy.ddmm) and time (hh.mmsss) to seconds since 1970-01-01
101 *
102 * ep2dt, epb2dt, epj2dt (epoch,date, time)
103 * Convert fractional year to date as yyyy.ddmm and time as hh.mmsss
104 * ep2fd, epb2fd, epj2fd (epoch)
105 * Convert epoch to FITS ISO date string
106 * ep2i, epb2i, epj2i (epoch,iyr,imon,iday,ihr,imn,sec, ndsec)
107 * Convert fractional year to year month day hours minutes seconds
108 * ep2jd, epb2jd, epj2jd (epoch)
109 * Convert fractional year as used in epoch to Julian date
110 * ep2mjd, epb2mjd, epj2mjd (epoch)
111 * Convert fractional year as used in epoch to modified Julian date
112 * ep2ts, epb2ts, epj2ts (epoch)
113 * Convert fractional year to seconds since 1950.0
114 *
115 * et2fd (string)
116 * Convert from ET (or TDT or TT) in FITS format to UT in FITS format
117 * fd2et (string)
118 * Convert from UT in FITS format to ET (or TDT or TT) in FITS format
119 * jd2jed (dj)
120 * Convert from Julian Date to Julian Ephemeris Date
121 * jed2jd (dj)
122 * Convert from Julian Ephemeris Date to Julian Date
123 * dt2et (date, time)
124 * Convert date (yyyy.ddmm) and time (hh.mmsss) to ephemeris time
125 * edt2dt (date, time)
126 * Convert ephemeris date (yyyy.ddmm) and time (hh.mmsss) to UT
127 * dt2tai (date, time)
128 * Convert date (yyyy.ddmm) and time (hh.mmsss) to TAI date and time
129 * tai2dt (date, time)
130 * Convert TAI date (yyyy.ddmm) and time (hh.mmsss) to UT
131 * ts2ets (tsec)
132 * Convert from UT in seconds since 1950-01-01 to ET in same format
133 * ets2ts (tsec)
134 * Convert from ET in seconds since 1950-01-01 to UT in same format
135 *
136 * fd2ep, fd2epb, fd2epj (string)
137 * Convert FITS date string to fractional year
138 * Convert time alone to fraction of Besselian year
139 * fd2doy (string, year, doy)
140 * Convert FITS standard date string to year and day of year
141 * fd2dt (string, date, time)
142 * Convert FITS date string to date as yyyy.ddmm and time as hh.mmsss
143 * Convert time alone to hh.mmssss with date set to 0.0
144 * fd2i (string,iyr,imon,iday,ihr,imn,sec, ndsec)
145 * Convert FITS standard date string to year month day hours min sec
146 * Convert time alone to hours min sec, year month day are zero
147 * fd2jd (string)
148 * Convert FITS standard date string to Julian date
149 * Convert time alone to fraction of day
150 * fd2mjd (string)
151 * Convert FITS standard date string to modified Julian date
152 * fd2ts (string)
153 * Convert FITS standard date string to seconds since 1950.0
154 * Convert time alone to seconds of day
155 * fd2fd (string)
156 * Convert FITS standard date string to ISO FITS date string
157 * fd2of (string)
158 * Convert FITS standard date string to old-format FITS date and time
159 * fd2ofd (string)
160 * Convert FITS standard date string to old-format FITS date string
161 * fd2oft (string)
162 * Convert time part of FITS standard date string to FITS date string
163 *
164 * jd2doy (dj, year, doy)
165 * Convert Julian date to year and day of year
166 * jd2dt (dj,date,time)
167 * Convert Julian date to date as yyyy.mmdd and time as hh.mmssss
168 * jd2ep, jd2epb, jd2epj (dj)
169 * Convert Julian date to fractional year as used in epoch
170 * jd2fd (dj)
171 * Convert Julian date to FITS ISO date string
172 * jd2i (dj,iyr,imon,iday,ihr,imn,sec, ndsec)
173 * Convert Julian date to year month day hours min sec
174 * jd2mjd (dj)
175 * Convert Julian date to modified Julian date
176 * jd2ts (dj)
177 * Convert Julian day to seconds since 1950.0
178 *
179 * lt2dt()
180 * Return local time as yyyy.mmdd and time as hh.mmssss
181 * lt2fd()
182 * Return local time as FITS ISO date string
183 * lt2tsi()
184 * Return local time as IRAF seconds since 1980-01-01 00:00
185 * lt2tsu()
186 * Return local time as Unix seconds since 1970-01-01 00:00
187 * lt2ts()
188 * Return local time as Unix seconds since 1950-01-01 00:00
189 *
190 * mjd2doy (dj,year,doy)
191 * Convert modified Julian date to date as year and day of year
192 * mjd2dt (dj,date,time)
193 * Convert modified Julian date to date as yyyy.mmdd and time as hh.mmssss
194 * mjd2ep, mjd2epb, mjd2epj (dj)
195 * Convert modified Julian date to fractional year as used in epoch
196 * mjd2fd (dj)
197 * Convert modified Julian date to FITS ISO date string
198 * mjd2i (dj,iyr,imon,iday,ihr,imn,sec, ndsec)
199 * Convert modified Julian date to year month day hours min sec
200 * mjd2jd (dj)
201 * Convert modified Julian date to Julian date
202 * mjd2ts (dj)
203 * Convert modified Julian day to seconds since 1950.0
204 *
205 * ts2dt (tsec,date,time)
206 * Convert seconds since 1950.0 to date as yyyy.ddmm and time as hh.mmsss
207 * ts2ep, ts2epb, ts2epj (tsec)
208 * Convert seconds since 1950.0 to fractional year
209 * ts2fd (tsec)
210 * Convert seconds since 1950.0 to FITS standard date string
211 * ts2i (tsec,iyr,imon,iday,ihr,imn,sec, ndsec)
212 * Convert sec since 1950.0 to year month day hours minutes seconds
213 * ts2jd (tsec)
214 * Convert seconds since 1950.0 to Julian date
215 * ts2mjd (tsec)
216 * Convert seconds since 1950.0 to modified Julian date
217 * tsi2fd (tsec)
218 * Convert seconds since 1980-01-01 to FITS standard date string
219 * tsi2dt (tsec,date,time)
220 * Convert seconds since 1980-01-01 to date as yyyy.ddmm, time as hh.mmsss
221 * tsu2fd (tsec)
222 * Convert seconds since 1970-01-01 to FITS standard date string
223 * tsu2tsi (tsec)
224 * Convert UT seconds since 1970-01-01 to local seconds since 1980-01-01
225 * tsu2dt (tsec,date,time)
226 * Convert seconds since 1970-01-01 to date as yyyy.ddmm, time as hh.mmsss
227 *
228 * tsd2fd (tsec)
229 * Convert seconds since start of day to FITS time, hh:mm:ss.ss
230 * tsd2dt (tsec)
231 * Convert seconds since start of day to hh.mmssss
232 *
233 * fd2gst (string)
234 * convert from FITS date Greenwich Sidereal Time
235 * dt2gst (date, time)
236 * convert from UT as yyyy.mmdd hh.mmssss to Greenwich Sidereal Time
237 * ts2gst (tsec)
238 * Calculate Greenwich Sidereal Time given Universal Time
239 * in seconds since 1951-01-01T0:00:00
240 * fd2mst (string)
241 * convert from FITS UT date to Mean Sidereal Time
242 * dt2gmt (date, time)
243 * convert from UT as yyyy.mmdd hh.mmssss to Mean Sidereal Time
244 * ts2mst (tsec)
245 * Calculate Mean Sidereal Time given Universal Time
246 * in seconds since 1951-01-01T0:00:00
247 * jd2mst (string)
248 * convert from Julian Date to Mean Sidereal Time
249 * mst2fd (string)
250 * convert to current UT in FITS format given Greenwich Mean Sidereal Time
251 * mst2jd (dj)
252 * convert to current UT as Julian Date given Greenwich Mean Sidereal Time
253 * jd2lst (dj)
254 * Calculate Local Sidereal Time from Julian Date
255 * ts2lst (tsec)
256 * Calculate Local Sidereal Time given UT in seconds since 1951-01-01T0:00
257 * fd2lst (string)
258 * Calculate Local Sidereal Time given Universal Time as FITS ISO date
259 * lst2jd (dj, lst)
260 * Calculate Julian Date given current Julian date and Local Sidereal Time
261 * lst2fd (string, lst)
262 * Calculate Julian Date given current UT date and Local Sidereal Time
263 * gst2fd (string)
264 * Calculate current UT given UT date and Greenwich Sidereal Time
265 * gst2jd (dj)
266 * Calculate current UT given UT date and Greenwich Sidereal Time as JD
267 *
268 * compnut (dj, dpsi, deps, eps0)
269 * Compute the longitude and obliquity components of nutation and
270 * mean obliquity from the IAU 1980 theory
271 *
272 * utdt (dj)
273 * Compute difference between UT and dynamical time (ET-UT)
274 * ut2dt (year, doy)
275 * Current Universal Time to year and day of year
276 * ut2dt (date, time)
277 * Current Universal Time to date (yyyy.mmdd) and time (hh.mmsss)
278 * ut2ep(), ut2epb(), ut2epj()
279 * Current Universal Time to fractional year, Besselian, Julian epoch
280 * ut2fd()
281 * Current Universal Time to FITS ISO date string
282 * ut2jd()
283 * Current Universal Time to Julian Date
284 * ut2mjd()
285 * Current Universal Time to Modified Julian Date
286 * ut2tsi()
287 * Current Universal Time to IRAF seconds since 1980-01-01T00:00
288 * ut2tsu()
289 * Current Universal Time to Unix seconds since 1970-01-01T00:00
290 * ut2ts()
291 * Current Universal Time to seconds since 1950-01-01T00:00
292 * isdate (string)
293 * Return 1 if string is a FITS date (old or ISO)
294 *
295 * Internally-used subroutines
296 *
297 * fixdate (iyr, imon, iday, ihr, imn, sec, ndsec)
298 * Round seconds and make sure date and time numbers are within limits
299 * caldays (year, month)
300 * Calculate days in month 1-12 given year (Gregorian calendar only
301 * dint (dnum)
302 * Return integer part of floating point number
303 * dmod (dnum)
304 * Return Mod of floating point number
305 */
306
307 #include <stdlib.h>
308 #include <stdio.h>
309 #include <string.h>
310 #include <math.h>
311 #include <time.h>
312 #include <sys/time.h>
313 #include "wcs.h"
314 #include "fitsfile.h"
315
316 static double suntl();
317 static void fixdate();
318 static int caldays();
319 static double dint();
320 static double dmod();
321
322 static double longitude = 0.0; /* longitude of observatory in degrees (+=west) */
323 void
setlongitude(longitude0)324 setlongitude (longitude0)
325 double longitude0;
326 { longitude = longitude0; return; }
327
328 static int ndec = 3;
329 void
setdatedec(nd)330 setdatedec (nd)
331 int nd;
332 { ndec = nd; return; }
333
334 /* ANG2HR -- Convert angle in fraction degrees to hours as hh:mm:ss.ss */
335
336 void
ang2hr(angle,lstr,string)337 ang2hr (angle, lstr, string)
338
339 double angle; /* Angle in fractional degrees */
340 int lstr; /* Maximum number of characters in string */
341 char *string; /* Character string (hh:mm:ss.ss returned) */
342
343 {
344 angle = angle / 15.0;
345 dec2str (string, lstr, angle, ndec);
346 return;
347 }
348
349
350 /* ANG2DEG -- Convert angle in fraction degrees to degrees as dd:mm:ss.ss */
351
352 void
ang2deg(angle,lstr,string)353 ang2deg (angle, lstr, string)
354
355 double angle; /* Angle in fractional degrees */
356 int lstr; /* Maximum number of characters in string */
357 char *string; /* Character string (dd:mm:ss.ss returned) */
358 {
359 dec2str (string, lstr, angle, ndec);
360 return;
361 }
362
363
364 /* DEG2ANG -- Convert angle in degrees as dd:mm:ss.ss to fractional degrees */
365
366 double
deg2ang(angle)367 deg2ang (angle)
368
369 char *angle; /* Angle as dd:mm:ss.ss */
370 {
371 double deg;
372
373 deg = str2dec (angle);
374 return (deg);
375 }
376
377 /* HR2ANG -- Convert angle in hours as hh:mm:ss.ss to fractional degrees */
378
379 double
hr2ang(angle)380 hr2ang (angle)
381
382 char *angle; /* Angle in sexigesimal hours (hh:mm:ss.sss) */
383
384 {
385 double deg;
386
387 deg = str2dec (angle);
388 deg = deg * 15.0;
389 return (deg);
390 }
391
392
393 /* DT2FD-- convert vigesimal date and time to FITS date, yyyy-mm-ddThh:mm:ss.ss */
394
395 char *
dt2fd(date,time)396 dt2fd (date, time)
397
398 double date; /* Date as yyyy.mmdd
399 yyyy = calendar year (e.g. 1973)
400 mm = calendar month (e.g. 04 = april)
401 dd = calendar day (e.g. 15) */
402 double time; /* Time as hh.mmssxxxx
403 *if time<0, it is time as -(fraction of a day)
404 hh = hour of day (0 .le. hh .le. 23)
405 nn = minutes (0 .le. nn .le. 59)
406 ss = seconds (0 .le. ss .le. 59)
407 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
408 {
409 int iyr,imon,iday,ihr,imn;
410 double sec;
411 int nf;
412 char *string;
413 char tstring[32], dstring[32];
414 char outform[64];
415
416 dt2i (date, time, &iyr,&imon,&iday,&ihr,&imn,&sec, ndec);
417
418 /* Convert to ISO date format */
419 string = (char *) calloc (32, sizeof (char));
420
421 /* Make time string */
422 if (time != 0.0 || ndec > 0) {
423 if (ndec == 0)
424 nf = 2;
425 else
426 nf = 3 + ndec;
427 if (ndec > 0) {
428 sprintf (outform, "%%02d:%%02d:%%0%d.%df", nf, ndec);
429 sprintf (tstring, outform, ihr, imn, sec);
430 }
431 else {
432 sprintf (outform, "%%02d:%%02d:%%0%dd", nf);
433 sprintf (tstring, outform, ihr, imn, (int)(sec+0.5));
434 }
435 }
436
437 /* Make date string */
438 if (date != 0.0)
439 sprintf (dstring, "%4d-%02d-%02d", iyr, imon, iday);
440
441 /* Make FITS (ISO) date string */
442 if (date == 0.0)
443 strcpy (string, tstring);
444 else if (time == 0.0 && ndec < 1)
445 strcpy (string, dstring);
446 else
447 sprintf (string, "%sT%s", dstring, tstring);
448
449 return (string);
450 }
451
452
453 /* DT2JD-- convert from date as yyyy.mmdd and time as hh.mmsss to Julian Date
454 * Return fractional days if date is zero */
455
456 double
dt2jd(date,time)457 dt2jd (date,time)
458
459 double date; /* Date as yyyy.mmdd
460 yyyy = calendar year (e.g. 1973)
461 mm = calendar month (e.g. 04 = april)
462 dd = calendar day (e.g. 15) */
463 double time; /* Time as hh.mmssxxxx
464 *if time<0, it is time as -(fraction of a day)
465 hh = hour of day (0 .le. hh .le. 23)
466 nn = minutes (0 .le. nn .le. 59)
467 ss = seconds (0 .le. ss .le. 59)
468 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
469 {
470 double dj; /* Julian date (returned) */
471 double tsec; /* seconds since 1950.0 */
472
473 tsec = dt2ts (date, time);
474 if (date == 0.0)
475 dj = tsec / 86400.0;
476 else
477 dj = ts2jd (tsec);
478
479 return (dj);
480 }
481
482
483 /* DT2MJD-- convert from date yyyy.mmdd time hh.mmsss to modified Julian Date
484 * Return fractional days if date is zero */
485
486 double
dt2mjd(date,time)487 dt2mjd (date,time)
488
489 double date; /* Date as yyyy.mmdd
490 yyyy = calendar year (e.g. 1973)
491 mm = calendar month (e.g. 04 = april)
492 dd = calendar day (e.g. 15) */
493 double time; /* Time as hh.mmssxxxx
494 *if time<0, it is time as -(fraction of a day)
495 hh = hour of day (0 .le. hh .le. 23)
496 nn = minutes (0 .le. nn .le. 59)
497 ss = seconds (0 .le. ss .le. 59)
498 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
499 {
500 double dj; /* Modified Julian date (returned) */
501 double tsec; /* seconds since 1950.0 */
502
503 tsec = dt2ts (date, time);
504 if (date == 0.0)
505 dj = tsec / 86400.0;
506 else
507 dj = ts2jd (tsec);
508
509 return (dj - 2400000.5);
510 }
511
512
513 /* HJD2JD-- convert Heliocentric Julian Date to (geocentric) Julian date */
514
515 double
hjd2jd(dj,ra,dec,sys)516 hjd2jd (dj, ra, dec, sys)
517
518 double dj; /* Heliocentric Julian date */
519 double ra; /* Right ascension (degrees) */
520 double dec; /* Declination (degrees) */
521 int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */
522 {
523 double lt; /* Light travel difference to the Sun (days) */
524
525 lt = suntl (dj, ra, dec, sys);
526
527 /* Return Heliocentric Julian Date */
528 return (dj - lt);
529 }
530
531
532 /* JD2HJD-- convert (geocentric) Julian date to Heliocentric Julian Date */
533
534 double
jd2hjd(dj,ra,dec,sys)535 jd2hjd (dj, ra, dec, sys)
536
537 double dj; /* Julian date (geocentric) */
538 double ra; /* Right ascension (degrees) */
539 double dec; /* Declination (degrees) */
540 int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */
541 {
542 double lt; /* Light travel difference to the Sun (days) */
543
544 lt = suntl (dj, ra, dec, sys);
545
546 /* Return Heliocentric Julian Date */
547 return (dj + lt);
548 }
549
550
551 /* MHJD2MJD-- convert modified Heliocentric Julian Date to
552 modified geocentric Julian date */
553
554 double
mhjd2mjd(mhjd,ra,dec,sys)555 mhjd2mjd (mhjd, ra, dec, sys)
556
557 double mhjd; /* Modified Heliocentric Julian date */
558 double ra; /* Right ascension (degrees) */
559 double dec; /* Declination (degrees) */
560 int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */
561 {
562 double lt; /* Light travel difference to the Sun (days) */
563 double hjd; /* Heliocentric Julian date */
564
565 hjd = mjd2jd (mhjd);
566
567 lt = suntl (hjd, ra, dec, sys);
568
569 /* Return Heliocentric Julian Date */
570 return (jd2mjd (hjd - lt));
571 }
572
573
574 /* MJD2MHJD-- convert modified geocentric Julian date tp
575 modified Heliocentric Julian Date */
576
577 double
mjd2mhjd(mjd,ra,dec,sys)578 mjd2mhjd (mjd, ra, dec, sys)
579
580 double mjd; /* Julian date (geocentric) */
581 double ra; /* Right ascension (degrees) */
582 double dec; /* Declination (degrees) */
583 int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */
584 {
585 double lt; /* Light travel difference to the Sun (days) */
586 double dj; /* Julian date (geocentric) */
587
588 dj = mjd2jd (mjd);
589
590 lt = suntl (dj, ra, dec, sys);
591
592 /* Return Heliocentric Julian Date */
593 return (jd2mjd (dj + lt));
594 }
595
596
597 /* SUNTL-- compute light travel time to heliocentric correction in days */
598 /* Translated into C from IRAF SPP noao.astutils.asttools.asthjd.x */
599
600 static double
suntl(dj,ra,dec,sys)601 suntl (dj, ra, dec, sys)
602
603 double dj; /* Julian date (geocentric) */
604 double ra; /* Right ascension (degrees) */
605 double dec; /* Declination (degrees) */
606 int sys; /* J2000, B1950, GALACTIC, ECLIPTIC */
607 {
608 double t; /* Number of Julian centuries since J1900 */
609 double manom; /* Mean anomaly of the Earth's orbit (degrees) */
610 double lperi; /* Mean longitude of perihelion (degrees) */
611 double oblq; /* Mean obliquity of the ecliptic (degrees) */
612 double eccen; /* Eccentricity of the Earth's orbit (dimensionless) */
613 double eccen2, eccen3;
614 double tanom; /* True anomaly (approximate formula) (radians) */
615 double slong; /* True longitude of the Sun from the Earth (radians) */
616 double rs; /* Distance to the sun (AU) */
617 double lt; /* Light travel difference to the Sun (days) */
618 double l; /* Longitude of star in orbital plane of Earth (radians) */
619 double b; /* Latitude of star in orbital plane of Earth (radians) */
620 double epoch; /* Epoch of obervation */
621 double rs1,rs2;
622
623 t = (dj - 2415020.0) / 36525.0;
624
625 /* Compute earth orbital parameters */
626 manom = 358.47583 + (t * (35999.04975 - t * (0.000150 + t * 0.000003)));
627 lperi = 101.22083 + (t * (1.7191733 + t * (0.000453 + t * 0.000003)));
628 oblq = 23.452294 - (t * (0.0130125 + t * (0.00000164 - t * 0.000000503)));
629 eccen = 0.01675104 - (t * (0.00004180 + t * 0.000000126));
630 eccen2 = eccen * eccen;
631 eccen3 = eccen * eccen2;
632
633 /* Convert to principle angles */
634 manom = manom - (360.0 * (dint) (manom / 360.0));
635 lperi = lperi - (360.0 * (dint) (lperi / 360.0));
636
637 /* Convert to radians */
638 manom = degrad (manom);
639 lperi = degrad (lperi);
640 oblq = degrad (oblq);
641
642 /* True anomaly */
643 tanom = manom + (2 * eccen - 0.25 * eccen3) * sin (manom) +
644 1.25 * eccen2 * sin (2 * manom) +
645 13./12. * eccen3 * sin (3 * manom);
646
647 /* Distance to the Sun */
648 rs1 = 1.0 - eccen2;
649 rs2 = 1.0 + (eccen * cos (tanom));
650 rs = rs1 / rs2;
651
652 /* True longitude of the Sun seen from the Earth */
653 slong = lperi + tanom + PI;
654
655 /* Longitude and latitude of star in orbital plane of the Earth */
656 epoch = jd2ep (dj);
657 wcscon (sys, WCS_ECLIPTIC, 0.0, 0.0, &ra, &dec, epoch);
658 l = degrad (ra);
659 b = degrad (dec);
660
661 /* Light travel difference to the Sun */
662 lt = -0.005770 * rs * cos (b) * cos (l - slong);
663
664 /* Return light travel difference */
665 return (lt);
666 }
667
668
669 /* JD2DT-- convert Julian date to date as yyyy.mmdd and time as hh.mmssss */
670
671 void
jd2dt(dj,date,time)672 jd2dt (dj,date,time)
673
674 double dj; /* Julian date */
675 double *date; /* Date as yyyy.mmdd (returned) */
676 double *time; /* Time as hh.mmssxxxx (returned) */
677 {
678 int iyr,imon,iday,ihr,imn;
679 double sec;
680
681 /* Convert Julian Date to date and time */
682 jd2i (dj, &iyr, &imon, &iday, &ihr, &imn, &sec, 4);
683
684 /* Convert date to yyyy.mmdd */
685 if (iyr < 0) {
686 *date = (double) (-iyr) + 0.01 * (double) imon + 0.0001 * (double) iday;
687 *date = -(*date);
688 }
689 else
690 *date = (double) iyr + 0.01 * (double) imon + 0.0001 * (double) iday;
691
692 /* Convert time to hh.mmssssss */
693 *time = (double) ihr + 0.01 * (double) imn + 0.0001 * sec;
694
695 return;
696 }
697
698
699 /* JD2I-- convert Julian date to date as year, month, and day, and time hours,
700 minutes, and seconds */
701 /* after Fliegel and Van Flander, CACM 11, 657 (1968) */
702
703
704 void
jd2i(dj,iyr,imon,iday,ihr,imn,sec,ndsec)705 jd2i (dj, iyr, imon, iday, ihr, imn, sec, ndsec)
706
707 double dj; /* Julian date */
708 int *iyr; /* year (returned) */
709 int *imon; /* month (returned) */
710 int *iday; /* day (returned) */
711 int *ihr; /* hours (returned) */
712 int *imn; /* minutes (returned) */
713 double *sec; /* seconds (returned) */
714 int ndsec; /* Number of decimal places in seconds (0=int) */
715
716 {
717 double tsec;
718 double frac, dts, ts, sday;
719 int jd, l, n, i, j;
720
721 tsec = jd2ts (dj);
722 /* ts2i (tsec, iyr, imon, iday, ihr, imn, sec, ndsec); */
723
724 /* Round seconds to 0 - 4 decimal places */
725 if (tsec < 0.0)
726 dts = -0.5;
727 else
728 dts = 0.5;
729 if (ndsec < 1)
730 ts = dint (tsec + dts);
731 else if (ndsec < 2)
732 ts = dint (tsec * 10.0 + dts) / 10.0;
733 else if (ndsec < 3)
734 ts = dint (tsec * 100.0 + dts) / 100.0;
735 else if (ndsec < 4)
736 ts = dint (tsec * 1000.0 + dts) / 1000.0;
737 else
738 ts = dint (tsec * 10000.0 + dts) / 10000.0;
739
740 /* Convert back to Julian Date */
741 dj = ts2jd (ts);
742
743 /* Compute time from fraction of a day */
744 frac = dmod (dj, 1.0);
745 if (frac < 0.5) {
746 jd = (int) (dj - frac);
747 sday = (frac + 0.5) * 86400.0;
748 }
749 else {
750 jd = (int) (dj - frac) + 1;
751 sday = (frac - 0.5) * 86400.0;
752 }
753
754 *ihr = (int) (sday / 3600.0);
755 sday = sday - (double) (*ihr * 3600);
756 *imn = (int) (sday / 60.0);
757 *sec = sday - (double) (*imn * 60);
758
759 /* Compute day, month, year */
760 l = jd + 68569;
761 n = (4 * l) / 146097;
762 l = l - (146097 * n + 3) / 4;
763 i = (4000 * (l + 1)) / 1461001;
764 l = l - (1461 * i) / 4 + 31;
765 j = (80 * l) / 2447;
766 *iday = l - (2447 * j) / 80;
767 l = j / 11;
768 *imon = j + 2 - (12 * l);
769 *iyr = 100 * (n - 49) + i + l;
770
771 return;
772 }
773
774
775 /* JD2MJD-- convert Julian Date to Modified Julian Date */
776
777 double
jd2mjd(dj)778 jd2mjd (dj)
779
780 double dj; /* Julian Date */
781
782 {
783 return (dj - 2400000.5);
784 }
785
786
787 /* JD2EP-- convert Julian date to fractional year as used in epoch */
788
789 double
jd2ep(dj)790 jd2ep (dj)
791
792 double dj; /* Julian date */
793
794 {
795 double date, time;
796 jd2dt (dj, &date, &time);
797 return (dt2ep (date, time));
798 }
799
800
801 /* JD2EPB-- convert Julian date to Besselian epoch */
802
803 double
jd2epb(dj)804 jd2epb (dj)
805
806 double dj; /* Julian date */
807
808 {
809 return (1900.0 + (dj - 2415020.31352) / 365.242198781);
810 }
811
812
813 /* JD2EPJ-- convert Julian date to Julian epoch */
814
815 double
jd2epj(dj)816 jd2epj (dj)
817
818 double dj; /* Julian date */
819
820 {
821 return (2000.0 + (dj - 2451545.0) / 365.25);
822 }
823
824
825 /* LT2DT-- Return local time as yyyy.mmdd and time as hh.mmssss */
826
827 void
lt2dt(date,time)828 lt2dt(date, time)
829
830 double *date; /* Date as yyyy.mmdd (returned) */
831 double *time; /* Time as hh.mmssxxxx (returned) */
832
833 {
834 time_t tsec;
835 struct timeval tp;
836 struct timezone tzp;
837 struct tm *ts;
838
839 gettimeofday (&tp,&tzp);
840
841 tsec = tp.tv_sec;
842 ts = localtime (&tsec);
843
844 if (ts->tm_year < 1000)
845 *date = (double) (ts->tm_year + 1900);
846 else
847 *date = (double) ts->tm_year;
848 *date = *date + (0.01 * (double) (ts->tm_mon + 1));
849 *date = *date + (0.0001 * (double) ts->tm_mday);
850 *time = (double) ts->tm_hour;
851 *time = *time + (0.01 * (double) ts->tm_min);
852 *time = *time + (0.0001 * (double) ts->tm_sec);
853
854 return;
855 }
856
857
858 /* LT2FD-- Return current local time as FITS ISO date string */
859
860 char *
lt2fd()861 lt2fd()
862 {
863 time_t tsec;
864 struct tm *ts;
865 struct timeval tp;
866 struct timezone tzp;
867 int month, day, year, hour, minute, second;
868 char *isotime;
869
870 gettimeofday (&tp,&tzp);
871 tsec = tp.tv_sec;
872
873 ts = localtime (&tsec);
874
875 year = ts->tm_year;
876 if (year < 1000)
877 year = year + 1900;
878 month = ts->tm_mon + 1;
879 day = ts->tm_mday;
880 hour = ts->tm_hour;
881 minute = ts->tm_min;
882 second = ts->tm_sec;
883
884 isotime = (char *) calloc (32, sizeof (char));
885 sprintf (isotime, "%04d-%02d-%02dT%02d:%02d:%02d",
886 year, month, day, hour, minute, second);
887 return (isotime);
888 }
889
890
891 /* LT2TSI-- Return local time as IRAF seconds since 1980-01-01 00:00 */
892
893 int
lt2tsi()894 lt2tsi()
895 {
896 return ((int)(lt2ts() - 946684800.0));
897 }
898
899
900 /* LT2TSU-- Return local time as Unix seconds since 1970-01-01 00:00 */
901
902 time_t
lt2tsu()903 lt2tsu()
904 {
905 return ((time_t)(lt2ts() - 631152000.0));
906 }
907
908 /* LT2TS-- Return local time as Unix seconds since 1950-01-01 00:00 */
909
910 double
lt2ts()911 lt2ts()
912 {
913 double tsec;
914 char *datestring;
915 datestring = lt2fd();
916 tsec = fd2ts (datestring);
917 free (datestring);
918 return (tsec);
919 }
920
921
922 /* MJD2DT-- convert Modified Julian Date to date (yyyy.mmdd) time (hh.mmssss) */
923
924 void
mjd2dt(dj,date,time)925 mjd2dt (dj,date,time)
926
927 double dj; /* Modified Julian Date */
928 double *date; /* Date as yyyy.mmdd (returned)
929 yyyy = calendar year (e.g. 1973)
930 mm = calendar month (e.g. 04 = april)
931 dd = calendar day (e.g. 15) */
932 double *time; /* Time as hh.mmssxxxx (returned)
933 *if time<0, it is time as -(fraction of a day)
934 hh = hour of day (0 .le. hh .le. 23)
935 nn = minutes (0 .le. nn .le. 59)
936 ss = seconds (0 .le. ss .le. 59)
937 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
938 {
939 double tsec;
940
941 tsec = jd2ts (dj + 2400000.5);
942 ts2dt (tsec, date, time);
943
944 return;
945 }
946
947
948 /* MJD2I-- convert Modified Julian Date to date as year, month, day and
949 time as hours, minutes, seconds */
950
951 void
mjd2i(dj,iyr,imon,iday,ihr,imn,sec,ndsec)952 mjd2i (dj, iyr, imon, iday, ihr, imn, sec, ndsec)
953
954 double dj; /* Modified Julian Date */
955 int *iyr; /* year (returned) */
956 int *imon; /* month (returned) */
957 int *iday; /* day (returned) */
958 int *ihr; /* hours (returned) */
959 int *imn; /* minutes (returned) */
960 double *sec; /* seconds (returned) */
961 int ndsec; /* Number of decimal places in seconds (0=int) */
962
963 {
964 double tsec;
965
966 tsec = jd2ts (dj + 2400000.5);
967 ts2i (tsec, iyr, imon, iday, ihr, imn, sec, ndsec);
968 return;
969 }
970
971
972 /* MJD2DOY-- convert Modified Julian Date to Year,Day-of-Year */
973
974 void
mjd2doy(dj,year,doy)975 mjd2doy (dj, year, doy)
976
977 double dj; /* Modified Julian Date */
978 int *year; /* Year (returned) */
979 double *doy; /* Day of year with fraction (returned) */
980
981 {
982 jd2doy (dj + 2400000.5, year, doy);
983 return;
984 }
985
986
987 /* MJD2JD-- convert Modified Julian Date to Julian Date */
988
989 double
mjd2jd(dj)990 mjd2jd (dj)
991
992 double dj; /* Modified Julian Date */
993
994 {
995 return (dj + 2400000.5);
996 }
997
998
999 /* MJD2EP-- convert Modified Julian Date to fractional year */
1000
1001 double
mjd2ep(dj)1002 mjd2ep (dj)
1003
1004 double dj; /* Modified Julian Date */
1005
1006 {
1007 double date, time;
1008 jd2dt (dj + 2400000.5, &date, &time);
1009 return (dt2ep (date, time));
1010 }
1011
1012
1013 /* MJD2EPB-- convert Modified Julian Date to Besselian epoch */
1014
1015 double
mjd2epb(dj)1016 mjd2epb (dj)
1017
1018 double dj; /* Modified Julian Date */
1019
1020 {
1021 return (1900.0 + (dj - 15019.81352) / 365.242198781);
1022 }
1023
1024
1025 /* MJD2EPJ-- convert Modified Julian Date to Julian epoch */
1026
1027 double
mjd2epj(dj)1028 mjd2epj (dj)
1029
1030 double dj; /* Modified Julian Date */
1031
1032 {
1033 return (2000.0 + (dj - 51544.5) / 365.25);
1034 }
1035
1036
1037 /* MJD2FD-- convert modified Julian date to FITS date, yyyy-mm-ddThh:mm:ss.ss */
1038
1039 char *
mjd2fd(dj)1040 mjd2fd (dj)
1041
1042 double dj; /* Modified Julian date */
1043 {
1044 return (jd2fd (dj + 2400000.5));
1045 }
1046
1047
1048 /* MJD2TS-- convert modified Julian date to seconds since 1950.0 */
1049
1050 double
mjd2ts(dj)1051 mjd2ts (dj)
1052
1053 double dj; /* Modified Julian date */
1054 {
1055 return ((dj - 33282.0) * 86400.0);
1056 }
1057
1058
1059 /* EP2FD-- convert fractional year to FITS date, yyyy-mm-ddThh:mm:ss.ss */
1060
1061 char *
ep2fd(epoch)1062 ep2fd (epoch)
1063
1064 double epoch; /* Date as fractional year */
1065 {
1066 double tsec; /* seconds since 1950.0 (returned) */
1067 tsec = ep2ts (epoch);
1068 return (ts2fd (tsec));
1069 }
1070
1071
1072 /* EPB2FD-- convert Besselian epoch to FITS date, yyyy-mm-ddThh:mm:ss.ss */
1073
1074 char *
epb2fd(epoch)1075 epb2fd (epoch)
1076
1077 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1078 {
1079 double dj; /* Julian Date */
1080 dj = epb2jd (epoch);
1081 return (jd2fd (dj));
1082 }
1083
1084
1085 /* EPJ2FD-- convert Julian epoch to FITS date, yyyy-mm-ddThh:mm:ss.ss */
1086
1087 char *
epj2fd(epoch)1088 epj2fd (epoch)
1089
1090 double epoch; /* Julian epoch (fractional 365.25-day years) */
1091 {
1092 double dj; /* Julian Date */
1093 dj = epj2jd (epoch);
1094 return (jd2fd (dj));
1095 }
1096
1097
1098 /* EP2TS-- convert fractional year to seconds since 1950.0 */
1099
1100 double
ep2ts(epoch)1101 ep2ts (epoch)
1102
1103 double epoch; /* Date as fractional year */
1104 {
1105 double dj;
1106 dj = ep2jd (epoch);
1107 return ((dj - 2433282.5) * 86400.0);
1108 }
1109
1110
1111 /* EPB2TS-- convert Besselian epoch to seconds since 1950.0 */
1112
1113 double
epb2ts(epoch)1114 epb2ts (epoch)
1115
1116 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1117 {
1118 double dj;
1119 dj = epb2jd (epoch);
1120 return ((dj - 2433282.5) * 86400.0);
1121 }
1122
1123
1124 /* EPJ2TS-- convert Julian epoch to seconds since 1950.0 */
1125
1126 double
epj2ts(epoch)1127 epj2ts (epoch)
1128
1129 double epoch; /* Julian epoch (fractional 365.25-day years) */
1130 {
1131 double dj;
1132 dj = epj2jd (epoch);
1133 return ((dj - 2433282.5) * 86400.0);
1134 }
1135
1136
1137 /* EPB2EP-- convert Besselian epoch to fractional years */
1138
1139 double
epb2ep(epoch)1140 epb2ep (epoch)
1141
1142 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1143 {
1144 double dj;
1145 dj = epb2jd (epoch);
1146 return (jd2ep (dj));
1147 }
1148
1149
1150 /* EP2EPB-- convert fractional year to Besselian epoch */
1151
1152 double
ep2epb(epoch)1153 ep2epb (epoch)
1154
1155 double epoch; /* Fractional year */
1156 {
1157 double dj;
1158 dj = ep2jd (epoch);
1159 return (jd2epb (dj));
1160 }
1161
1162
1163 /* EPJ2EP-- convert Julian epoch to fractional year */
1164
1165 double
epj2ep(epoch)1166 epj2ep (epoch)
1167
1168 double epoch; /* Julian epoch (fractional 365.25-day years) */
1169 {
1170 double dj;
1171 dj = epj2jd (epoch);
1172 return (jd2ep (dj));
1173 }
1174
1175
1176 /* EP2EPJ-- convert fractional year to Julian epoch */
1177
1178 double
ep2epj(epoch)1179 ep2epj (epoch)
1180
1181 double epoch; /* Fractional year */
1182 {
1183 double dj;
1184 dj = ep2jd (epoch);
1185 return (jd2epj (dj));
1186 }
1187
1188
1189 /* EP2I-- convert fractional year to year month day hours min sec */
1190
1191 void
ep2i(epoch,iyr,imon,iday,ihr,imn,sec,ndsec)1192 ep2i (epoch, iyr, imon, iday, ihr, imn, sec, ndsec)
1193
1194 double epoch; /* Date as fractional year */
1195 int *iyr; /* year (returned) */
1196 int *imon; /* month (returned) */
1197 int *iday; /* day (returned) */
1198 int *ihr; /* hours (returned) */
1199 int *imn; /* minutes (returned) */
1200 double *sec; /* seconds (returned) */
1201 int ndsec; /* Number of decimal places in seconds (0=int) */
1202 {
1203 double date, time;
1204
1205 ep2dt (epoch, &date, &time);
1206 dt2i (date, time, iyr,imon,iday,ihr,imn,sec, ndsec);
1207 return;
1208 }
1209
1210
1211 /* EPB2I-- convert Besselian epoch to year month day hours min sec */
1212
1213 void
epb2i(epoch,iyr,imon,iday,ihr,imn,sec,ndsec)1214 epb2i (epoch, iyr, imon, iday, ihr, imn, sec, ndsec)
1215
1216 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1217 int *iyr; /* year (returned) */
1218 int *imon; /* month (returned) */
1219 int *iday; /* day (returned) */
1220 int *ihr; /* hours (returned) */
1221 int *imn; /* minutes (returned) */
1222 double *sec; /* seconds (returned) */
1223 int ndsec; /* Number of decimal places in seconds (0=int) */
1224 {
1225 double date, time;
1226
1227 epb2dt (epoch, &date, &time);
1228 dt2i (date, time, iyr,imon,iday,ihr,imn,sec, ndsec);
1229 return;
1230 }
1231
1232
1233 /* EPJ2I-- convert Julian epoch to year month day hours min sec */
1234
1235 void
epj2i(epoch,iyr,imon,iday,ihr,imn,sec,ndsec)1236 epj2i (epoch, iyr, imon, iday, ihr, imn, sec, ndsec)
1237
1238 double epoch; /* Julian epoch (fractional 365.25-day years) */
1239 int *iyr; /* year (returned) */
1240 int *imon; /* month (returned) */
1241 int *iday; /* day (returned) */
1242 int *ihr; /* hours (returned) */
1243 int *imn; /* minutes (returned) */
1244 double *sec; /* seconds (returned) */
1245 int ndsec; /* Number of decimal places in seconds (0=int) */
1246 {
1247 double date, time;
1248
1249 epj2dt (epoch, &date, &time);
1250 dt2i (date, time, iyr,imon,iday,ihr,imn,sec, ndsec);
1251 return;
1252 }
1253
1254
1255 /* EP2JD-- convert fractional year as used in epoch to Julian date */
1256
1257 double
ep2jd(epoch)1258 ep2jd (epoch)
1259
1260 double epoch; /* Date as fractional year */
1261
1262 {
1263 double dj; /* Julian date (returned)*/
1264 double date, time;
1265
1266 ep2dt (epoch, &date, &time);
1267 dj = dt2jd (date, time);
1268 return (dj);
1269 }
1270
1271
1272 /* EPB2JD-- convert Besselian epoch to Julian Date */
1273
1274 double
epb2jd(epoch)1275 epb2jd (epoch)
1276
1277 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1278
1279 {
1280 return (2415020.31352 + ((epoch - 1900.0) * 365.242198781));
1281 }
1282
1283
1284 /* EPJ2JD-- convert Julian epoch to Julian Date */
1285
1286 double
epj2jd(epoch)1287 epj2jd (epoch)
1288
1289 double epoch; /* Julian epoch (fractional 365.25-day years) */
1290
1291 {
1292 return (2451545.0 + ((epoch - 2000.0) * 365.25));
1293 }
1294
1295
1296 /* EP2MJD-- convert fractional year as used in epoch to modified Julian date */
1297
1298 double
ep2mjd(epoch)1299 ep2mjd (epoch)
1300
1301 double epoch; /* Date as fractional year */
1302
1303 {
1304 double dj; /* Julian date (returned)*/
1305 double date, time;
1306
1307 ep2dt (epoch, &date, &time);
1308 dj = dt2jd (date, time);
1309 return (dj - 2400000.5);
1310 }
1311
1312
1313 /* EPB2MJD-- convert Besselian epoch to modified Julian Date */
1314
1315 double
epb2mjd(epoch)1316 epb2mjd (epoch)
1317
1318 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1319
1320 {
1321 return (15019.81352 + ((epoch - 1900.0) * 365.242198781));
1322 }
1323
1324
1325 /* EPJ2MJD-- convert Julian epoch to modified Julian Date */
1326
1327 double
epj2mjd(epoch)1328 epj2mjd (epoch)
1329
1330 double epoch; /* Julian epoch (fractional 365.25-day years) */
1331
1332 {
1333 return (51544.5 + ((epoch - 2000.0) * 365.25));
1334 }
1335
1336
1337
1338 /* EPB2EPJ-- convert Besselian epoch to Julian epoch */
1339
1340 double
epb2epj(epoch)1341 epb2epj (epoch)
1342
1343 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1344 {
1345 double dj; /* Julian date */
1346 dj = epb2jd (epoch);
1347 return (jd2epj (dj));
1348 }
1349
1350
1351 /* EPJ2EPB-- convert Julian epoch to Besselian epoch */
1352
1353 double
epj2epb(epoch)1354 epj2epb (epoch)
1355
1356 double epoch; /* Julian epoch (fractional 365.25-day years) */
1357 {
1358 double dj; /* Julian date */
1359 dj = epj2jd (epoch);
1360 return (jd2epb (dj));
1361 }
1362
1363
1364 /* JD2FD-- convert Julian date to FITS date, yyyy-mm-ddThh:mm:ss.ss */
1365
1366 char *
jd2fd(dj)1367 jd2fd (dj)
1368
1369 double dj; /* Julian date */
1370 {
1371 double tsec; /* seconds since 1950.0 (returned) */
1372 tsec = (dj - 2433282.5) * 86400.0;
1373 return (ts2fd (tsec));
1374 }
1375
1376
1377 /* JD2TS-- convert Julian date to seconds since 1950.0 */
1378
1379 double
jd2ts(dj)1380 jd2ts (dj)
1381
1382 double dj; /* Julian date */
1383 {
1384 return ((dj - 2433282.5) * 86400.0);
1385 }
1386
1387
1388 /* JD2TSI-- convert Julian date to IRAF seconds since 1980-01-01T0:00 */
1389
1390 int
jd2tsi(dj)1391 jd2tsi (dj)
1392
1393 double dj; /* Julian date */
1394 {
1395 double ts;
1396 ts = (dj - 2444239.5) * 86400.0;
1397 return ((int) ts);
1398 }
1399
1400
1401 /* JD2TSU-- convert Julian date to Unix seconds since 1970-01-01T0:00 */
1402
1403 time_t
jd2tsu(dj)1404 jd2tsu (dj)
1405
1406 double dj; /* Julian date */
1407 {
1408 return ((time_t)((dj - 2440587.5) * 86400.0));
1409 }
1410
1411
1412 /* DT2DOY-- convert yyyy.mmdd hh.mmss to year and day of year */
1413
1414 void
dt2doy(date,time,year,doy)1415 dt2doy (date, time, year, doy)
1416
1417 double date; /* Date as yyyy.mmdd */
1418 double time; /* Time as hh.mmssxxxx */
1419 int *year; /* Year (returned) */
1420 double *doy; /* Day of year with fraction (returned) */
1421 {
1422 double dj; /* Julian date */
1423 double dj0; /* Julian date on January 1 0:00 */
1424 double date0; /* January first of date's year */
1425 double dyear;
1426
1427 dyear = floor (date);
1428 date0 = dyear + 0.0101;
1429 dj0 = dt2jd (date0, 0.0);
1430 dj = dt2jd (date, time);
1431 *year = (int) (dyear + 0.00000001);
1432 *doy = dj - dj0 + 1.0;
1433 return;
1434 }
1435
1436
1437 /* DOY2DT-- convert year and day of year to yyyy.mmdd hh.mmss */
1438
1439 void
doy2dt(year,doy,date,time)1440 doy2dt (year, doy, date, time)
1441
1442 int year; /* Year */
1443 double doy; /* Day of year with fraction */
1444 double *date; /* Date as yyyy.mmdd (returned) */
1445 double *time; /* Time as hh.mmssxxxx (returned) */
1446 {
1447 double dj; /* Julian date */
1448 double dj0; /* Julian date on January 1 0:00 */
1449 double date0; /* January first of date's year */
1450
1451 date0 = year + 0.0101;
1452 dj0 = dt2jd (date0, 0.0);
1453 dj = dj0 + doy - 1.0;
1454 jd2dt (dj, date, time);
1455 return;
1456 }
1457
1458
1459 /* DOY2EP-- convert year and day of year to fractional year as used in epoch */
1460
1461 double
doy2ep(year,doy)1462 doy2ep (year, doy)
1463
1464 int year; /* Year */
1465 double doy; /* Day of year with fraction */
1466 {
1467 double date, time;
1468 doy2dt (year, doy, &date, &time);
1469 return (dt2ep (date, time));
1470 }
1471
1472
1473
1474 /* DOY2EPB-- convert year and day of year to Besellian epoch */
1475
1476 double
doy2epb(year,doy)1477 doy2epb (year, doy)
1478
1479 int year; /* Year */
1480 double doy; /* Day of year with fraction */
1481 {
1482 double dj;
1483 dj = doy2jd (year, doy);
1484 return (jd2epb (dj));
1485 }
1486
1487
1488 /* DOY2EPJ-- convert year and day of year to Julian epoch */
1489
1490 double
doy2epj(year,doy)1491 doy2epj (year, doy)
1492
1493 int year; /* Year */
1494 double doy; /* Day of year with fraction */
1495 {
1496 double dj;
1497 dj = doy2jd (year, doy);
1498 return (jd2epj (dj));
1499 }
1500
1501
1502 /* DOY2FD-- convert year and day of year to FITS date */
1503
1504 char *
doy2fd(year,doy)1505 doy2fd (year, doy)
1506
1507 int year; /* Year */
1508 double doy; /* Day of year with fraction */
1509 {
1510 double dj; /* Julian date */
1511
1512 dj = doy2jd (year, doy);
1513 return (jd2fd (dj));
1514 }
1515
1516
1517 /* DOY2JD-- convert year and day of year to Julian date */
1518
1519 double
doy2jd(year,doy)1520 doy2jd (year, doy)
1521
1522 int year; /* Year */
1523 double doy; /* Day of year with fraction */
1524 {
1525 double dj0; /* Julian date */
1526 double date; /* Date as yyyy.mmdd (returned) */
1527 double time; /* Time as hh.mmssxxxx (returned) */
1528
1529 date = (double) year + 0.0101;
1530 time = 0.0;
1531 dj0 = dt2jd (date, time);
1532 return (dj0 + doy - 1.0);
1533 }
1534
1535
1536 /* DOY2MJD-- convert year and day of year to Julian date */
1537
1538 double
doy2mjd(year,doy)1539 doy2mjd (year, doy)
1540
1541 int year; /* Year */
1542 double doy; /* Day of year with fraction */
1543 {
1544 double dj0; /* Julian date */
1545 double date; /* Date as yyyy.mmdd (returned) */
1546 double time; /* Time as hh.mmssxxxx (returned) */
1547
1548 date = (double) year + 0.0101;
1549 time = 0.0;
1550 dj0 = dt2jd (date, time);
1551 return (dj0 + doy - 1.0 - 2400000.5);
1552 }
1553
1554
1555 /* DOY2TSU-- convert from FITS date to Unix seconds since 1970-01-01T0:00 */
1556
1557 time_t
doy2tsu(year,doy)1558 doy2tsu (year, doy)
1559
1560 int year; /* Year */
1561 double doy; /* Day of year with fraction */
1562 {
1563 double dj;
1564 dj = doy2jd (year, doy);
1565 return ((time_t)jd2ts (dj));
1566 }
1567
1568
1569 /* DOY2TSI-- convert from FITS date to IRAF seconds since 1980-01-01T0:00 */
1570
1571 int
doy2tsi(year,doy)1572 doy2tsi (year, doy)
1573
1574 int year; /* Year */
1575 double doy; /* Day of year with fraction */
1576 {
1577 double dj;
1578 dj = doy2jd (year, doy);
1579 return ((int)jd2tsi (dj));
1580 }
1581
1582
1583 /* DOY2TS-- convert year, day of year to seconds since 1950 */
1584
1585 double
doy2ts(year,doy)1586 doy2ts (year, doy)
1587
1588 int year; /* Year */
1589 double doy; /* Day of year with fraction */
1590 {
1591 double dj;
1592 dj = doy2jd (year, doy);
1593 return (jd2ts (dj));
1594 }
1595
1596
1597 /* FD2DOY-- convert FITS date to year and day of year */
1598
1599 void
fd2doy(string,year,doy)1600 fd2doy (string, year, doy)
1601
1602 char *string; /* FITS date string, which may be:
1603 fractional year
1604 dd/mm/yy (FITS standard before 2000)
1605 dd-mm-yy (nonstandard use before 2000)
1606 yyyy-mm-dd (FITS standard after 1999)
1607 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1608 int *year; /* Year (returned) */
1609 double *doy; /* Day of year with fraction (returned) */
1610 {
1611 double dj; /* Julian date */
1612
1613 dj = fd2jd (string);
1614 jd2doy (dj, year, doy);
1615 return;
1616 }
1617
1618
1619 /* JD2DOY-- convert Julian date to year and day of year */
1620
1621 void
jd2doy(dj,year,doy)1622 jd2doy (dj, year, doy)
1623
1624 double dj; /* Julian date */
1625 int *year; /* Year (returned) */
1626 double *doy; /* Day of year with fraction (returned) */
1627 {
1628 double date; /* Date as yyyy.mmdd (returned) */
1629 double time; /* Time as hh.mmssxxxx (returned) */
1630 double dj0; /* Julian date at 0:00 on 1/1 */
1631 double dyear;
1632
1633 jd2dt (dj, &date, &time);
1634 *year = (int) date;
1635 dyear = (double) *year;
1636 dj0 = dt2jd (dyear+0.0101, 0.0);
1637 *doy = dj - dj0 + 1.0;
1638 return;
1639 }
1640
1641
1642 /* TS2JD-- convert seconds since 1950.0 to Julian date */
1643
1644 double
ts2jd(tsec)1645 ts2jd (tsec)
1646
1647 double tsec; /* seconds since 1950.0 */
1648 {
1649 return (2433282.5 + (tsec / 86400.0));
1650 }
1651
1652
1653 /* TS2MJD-- convert seconds since 1950.0 to modified Julian date */
1654
1655 double
ts2mjd(tsec)1656 ts2mjd (tsec)
1657
1658 double tsec; /* seconds since 1950.0 */
1659 {
1660 return (33282.0 + (tsec / 86400.0));
1661 }
1662
1663
1664 /* TS2EP-- convert seconds since 1950.0 to fractional year as used in epoch */
1665
1666 double
ts2ep(tsec)1667 ts2ep (tsec)
1668
1669 double tsec; /* Seconds since 1950.0 */
1670
1671 {
1672 double date, time;
1673 ts2dt (tsec, &date, &time);
1674 return (dt2ep (date, time));
1675 }
1676
1677
1678 /* TS2EPB-- convert seconds since 1950.0 to Besselian epoch */
1679
1680 double
ts2epb(tsec)1681 ts2epb (tsec)
1682
1683 double tsec; /* Seconds since 1950.0 */
1684
1685 {
1686 double dj; /* Julian Date */
1687 dj = ts2jd (tsec);
1688 return (jd2epb (dj));
1689 }
1690
1691
1692 /* TS2EPB-- convert seconds since 1950.0 to Julian epoch */
1693
1694 double
ts2epj(tsec)1695 ts2epj (tsec)
1696
1697 double tsec; /* Seconds since 1950.0 */
1698
1699 {
1700 double dj; /* Julian Date */
1701 dj = ts2jd (tsec);
1702 return (jd2epj (dj));
1703 }
1704
1705
1706 /* DT2EP-- convert from date, time as yyyy.mmdd hh.mmsss to fractional year */
1707
1708 double
dt2ep(date,time)1709 dt2ep (date, time)
1710
1711 double date; /* Date as yyyy.mmdd
1712 yyyy = calendar year (e.g. 1973)
1713 mm = calendar month (e.g. 04 = april)
1714 dd = calendar day (e.g. 15) */
1715 double time; /* Time as hh.mmssxxxx
1716 *if time<0, it is time as -(fraction of a day)
1717 hh = hour of day (0 .le. hh .le. 23)
1718 nn = minutes (0 .le. nn .le. 59)
1719 ss = seconds (0 .le. ss .le. 59)
1720 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
1721 {
1722 double epoch; /* Date as fractional year (returned) */
1723 double dj, dj0, dj1, date0, time0, date1;
1724
1725 dj = dt2jd (date, time);
1726 if (date == 0.0)
1727 epoch = dj / 365.2422;
1728 else {
1729 time0 = 0.0;
1730 date0 = dint (date) + 0.0101;
1731 date1 = dint (date) + 1.0101;
1732 dj0 = dt2jd (date0, time0);
1733 dj1 = dt2jd (date1, time0);
1734 epoch = dint (date) + ((dj - dj0) / (dj1 - dj0));
1735 }
1736 return (epoch);
1737 }
1738
1739
1740 /* DT2EPB-- convert from date, time as yyyy.mmdd hh.mmsss to Besselian epoch */
1741
1742 double
dt2epb(date,time)1743 dt2epb (date, time)
1744
1745 double date; /* Date as yyyy.mmdd
1746 yyyy = calendar year (e.g. 1973)
1747 mm = calendar month (e.g. 04 = april)
1748 dd = calendar day (e.g. 15) */
1749 double time; /* Time as hh.mmssxxxx
1750 *if time<0, it is time as -(fraction of a day)
1751 hh = hour of day (0 .le. hh .le. 23)
1752 nn = minutes (0 .le. nn .le. 59)
1753 ss = seconds (0 .le. ss .le. 59)
1754 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
1755 {
1756 double dj; /* Julian date */
1757 double epoch; /* Date as fractional year (returned) */
1758 dj = dt2jd (date, time);
1759 if (date == 0.0)
1760 epoch = dj / 365.242198781;
1761 else
1762 epoch = jd2epb (dj);
1763 return (epoch);
1764 }
1765
1766
1767 /* DT2EPJ-- convert from date, time as yyyy.mmdd hh.mmsss to Julian epoch */
1768
1769 double
dt2epj(date,time)1770 dt2epj (date, time)
1771
1772 double date; /* Date as yyyy.mmdd
1773 yyyy = calendar year (e.g. 1973)
1774 mm = calendar month (e.g. 04 = april)
1775 dd = calendar day (e.g. 15) */
1776 double time; /* Time as hh.mmssxxxx
1777 *if time<0, it is time as -(fraction of a day)
1778 hh = hour of day (0 .le. hh .le. 23)
1779 nn = minutes (0 .le. nn .le. 59)
1780 ss = seconds (0 .le. ss .le. 59)
1781 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
1782 {
1783 double dj; /* Julian date */
1784 double epoch; /* Date as fractional year (returned) */
1785 dj = dt2jd (date, time);
1786 if (date == 0.0)
1787 epoch = dj / 365.25;
1788 else
1789 epoch = jd2epj (dj);
1790 return (epoch);
1791 }
1792
1793
1794 /* EP2DT-- convert from fractional year to date, time as yyyy.mmdd hh.mmsss */
1795
1796 void
ep2dt(epoch,date,time)1797 ep2dt (epoch, date, time)
1798
1799 double epoch; /* Date as fractional year */
1800 double *date; /* Date as yyyy.mmdd (returned)
1801 yyyy = calendar year (e.g. 1973)
1802 mm = calendar month (e.g. 04 = april)
1803 dd = calendar day (e.g. 15) */
1804 double *time; /* Time as hh.mmssxxxx (returned)
1805 *if time<0, it is time as -(fraction of a day)
1806 hh = hour of day (0 .le. hh .le. 23)
1807 nn = minutes (0 .le. nn .le. 59)
1808 ss = seconds (0 .le. ss .le. 59)
1809 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
1810 {
1811 double dj, dj0, dj1, date0, time0, date1, epochi, epochf;
1812
1813 time0 = 0.0;
1814 epochi = dint (epoch);
1815 epochf = epoch - epochi;
1816 date0 = epochi + 0.0101;
1817 date1 = epochi + 1.0101;
1818 dj0 = dt2jd (date0, time0);
1819 dj1 = dt2jd (date1, time0);
1820 dj = dj0 + epochf * (dj1 - dj0);
1821 jd2dt (dj, date, time);
1822 return;
1823 }
1824
1825
1826 /* EPB2DT-- convert from Besselian epoch to date, time as yyyy.mmdd hh.mmsss */
1827
1828 void
epb2dt(epoch,date,time)1829 epb2dt (epoch, date, time)
1830
1831 double epoch; /* Besselian epoch (fractional 365.242198781-day years) */
1832 double *date; /* Date as yyyy.mmdd (returned)
1833 yyyy = calendar year (e.g. 1973)
1834 mm = calendar month (e.g. 04 = april)
1835 dd = calendar day (e.g. 15) */
1836 double *time; /* Time as hh.mmssxxxx (returned)
1837 *if time<0, it is time as -(fraction of a day)
1838 hh = hour of day (0 .le. hh .le. 23)
1839 nn = minutes (0 .le. nn .le. 59)
1840 ss = seconds (0 .le. ss .le. 59)
1841 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
1842 {
1843 double dj; /* Julian date */
1844 dj = epb2jd (epoch);
1845 jd2dt (dj, date, time);
1846 }
1847
1848
1849 /* EPJ2DT-- convert from Julian epoch to date, time as yyyy.mmdd hh.mmsss */
1850
1851 void
epj2dt(epoch,date,time)1852 epj2dt (epoch, date, time)
1853
1854 double epoch; /* Julian epoch (fractional 365.25-day years) */
1855 double *date; /* Date as yyyy.mmdd (returned)
1856 yyyy = calendar year (e.g. 1973)
1857 mm = calendar month (e.g. 04 = april)
1858 dd = calendar day (e.g. 15) */
1859 double *time; /* Time as hh.mmssxxxx (returned)
1860 *if time<0, it is time as -(fraction of a day)
1861 hh = hour of day (0 .le. hh .le. 23)
1862 nn = minutes (0 .le. nn .le. 59)
1863 ss = seconds (0 .le. ss .le. 59)
1864 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
1865 {
1866 double dj; /* Julian date */
1867 dj = epj2jd (epoch);
1868 jd2dt (dj, date, time);
1869 }
1870
1871
1872 /* FD2JD-- convert FITS standard date to Julian date */
1873
1874 double
fd2jd(string)1875 fd2jd (string)
1876
1877 char *string; /* FITS date string, which may be:
1878 fractional year
1879 dd/mm/yy (FITS standard before 2000)
1880 dd-mm-yy (nonstandard use before 2000)
1881 yyyy-mm-dd (FITS standard after 1999)
1882 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1883 {
1884 double date, time;
1885
1886 fd2dt (string, &date, &time);
1887 return (dt2jd (date, time));
1888 }
1889
1890
1891 /* FD2MJD-- convert FITS standard date to modified Julian date */
1892
1893 double
fd2mjd(string)1894 fd2mjd (string)
1895
1896 char *string; /* FITS date string, which may be:
1897 fractional year
1898 dd/mm/yy (FITS standard before 2000)
1899 dd-mm-yy (nonstandard use before 2000)
1900 yyyy-mm-dd (FITS standard after 1999)
1901 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1902 {
1903 return (fd2jd (string) - 2400000.5);
1904 }
1905
1906
1907 /* FD2TSU-- convert from FITS date to Unix seconds since 1970-01-01T0:00 */
1908
1909 time_t
fd2tsu(string)1910 fd2tsu (string)
1911
1912 char *string; /* FITS date string, which may be:
1913 fractional year
1914 dd/mm/yy (FITS standard before 2000)
1915 dd-mm-yy (nonstandard use before 2000)
1916 yyyy-mm-dd (FITS standard after 1999)
1917 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1918 {
1919 double date, time;
1920 fd2dt (string, &date, &time);
1921 return (dt2tsu (date, time));
1922 }
1923
1924
1925 /* FD2TSI-- convert from FITS date to IRAF seconds since 1980-01-01T0:00 */
1926
1927 int
fd2tsi(string)1928 fd2tsi (string)
1929
1930 char *string; /* FITS date string, which may be:
1931 fractional year
1932 dd/mm/yy (FITS standard before 2000)
1933 dd-mm-yy (nonstandard use before 2000)
1934 yyyy-mm-dd (FITS standard after 1999)
1935 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1936 {
1937 double date, time;
1938 fd2dt (string, &date, &time);
1939 return (dt2tsi (date, time));
1940 }
1941
1942
1943 /* FD2TS-- convert FITS standard date to seconds since 1950 */
1944
1945 double
fd2ts(string)1946 fd2ts (string)
1947
1948 char *string; /* FITS date string, which may be:
1949 fractional year
1950 dd/mm/yy (FITS standard before 2000)
1951 dd-mm-yy (nonstandard use before 2000)
1952 yyyy-mm-dd (FITS standard after 1999)
1953 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1954 {
1955 double date, time;
1956 fd2dt (string, &date, &time);
1957 return (dt2ts (date, time));
1958 }
1959
1960
1961 /* FD2FD-- convert any FITS standard date to ISO FITS standard date */
1962
1963 char *
fd2fd(string)1964 fd2fd (string)
1965
1966 char *string; /* FITS date string, which may be:
1967 fractional year
1968 dd/mm/yy (FITS standard before 2000)
1969 dd-mm-yy (nonstandard use before 2000)
1970 yyyy-mm-dd (FITS standard after 1999)
1971 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1972 {
1973 double date, time;
1974 fd2dt (string, &date, &time);
1975 return (dt2fd (date, time));
1976 }
1977
1978
1979 /* FD2OF-- convert any FITS standard date to old FITS standard date time */
1980
1981 char *
fd2of(string)1982 fd2of (string)
1983
1984 char *string; /* FITS date string, which may be:
1985 fractional year
1986 dd/mm/yy (FITS standard before 2000)
1987 dd-mm-yy (nonstandard use before 2000)
1988 yyyy-mm-dd (FITS standard after 1999)
1989 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
1990 {
1991 int iyr,imon,iday,ihr,imn;
1992 double sec;
1993
1994 fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 3);
1995
1996 /* Convert to old FITS date format */
1997 string = (char *) calloc (32, sizeof (char));
1998 if (iyr < 1900)
1999 sprintf (string, "*** date out of range ***");
2000 else if (iyr < 2000)
2001 sprintf (string, "%02d/%02d/%02d %02d:%02d:%06.3f",
2002 iday, imon, iyr-1900, ihr, imn, sec);
2003 else if (iyr < 2900.0)
2004 sprintf (string, "%02d/%02d/%3d %02d:%02d:%6.3f",
2005 iday, imon, iyr-1900, ihr, imn, sec);
2006 else
2007 sprintf (string, "*** date out of range ***");
2008 return (string);
2009 }
2010
2011
2012 /* TAI-UTC from the U.S. Naval Observatory */
2013 /* ftp://maia.usno.navy.mil/ser7/tai-utc.dat */
2014 static double taijd[26]={2441317.5, 2441499.5, 2441683.5, 2442048.5, 2442413.5,
2015 2442778.5, 2443144.5, 2443509.5, 2443874.5, 2444239.5, 2444786.5,
2016 2445151.5, 2445516.5, 2446247.5, 2447161.5, 2447892.5, 2448257.5,
2017 2448804.5, 2449169.5, 2449534.5, 2450083.5, 2450630.5, 2451179.5,
2018 2453736.5, 2454832.5, 2456293.5};
2019 static double taidt[26]={10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,
2020 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,
2021 33.0,34.0,35.0};
2022 static double dttab[173]={13.7,13.4,13.1,12.9,12.7,12.6,12.5,12.5,12.5,12.5,
2023 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,
2024 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,
2025 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,
2026 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,
2027 7.54, 6.97, 6.40, 6.02, 5.41, 4.10, 2.92, 1.82, 1.61, 0.10,-1.02,
2028 -1.28,-2.69,-3.24,-3.64,-4.54,-4.71,-5.11,-5.40,-5.42,-5.20,-5.46,
2029 -5.46,-5.79,-5.63,-5.64,-5.80,-5.66,-5.87,-6.01,-6.19,-6.64,-6.44,
2030 -6.47,-6.09,-5.76,-4.66,-3.74,-2.72,-1.54,-0.02, 1.24, 2.64, 3.86,
2031 5.37, 6.14, 7.75, 9.13,10.46,11.53,13.36,14.65,16.01,17.20,18.24,
2032 19.06,20.25,20.95,21.16,22.25,22.41,23.03,23.49,23.62,23.86,24.49,
2033 24.34,24.08,24.02,24.00,23.87,23.95,23.86,23.93,23.73,23.92,23.96,
2034 24.02,24.33,24.83,25.30,25.70,26.24,26.77,27.28,27.78,28.25,28.71,
2035 29.15,29.57,29.97,30.36,30.72,31.07,31.35,31.68,32.18,32.68,33.15,
2036 33.59,34.00,34.47,35.03,35.73,36.54,37.43,38.29,39.20,40.18,41.17,
2037 42.23};
2038
2039
2040 /* TAI2FD-- convert from TAI in FITS format to UT in FITS format */
2041
2042 char *
tai2fd(string)2043 tai2fd (string)
2044
2045 char *string; /* FITS date string, which may be:
2046 fractional year
2047 dd/mm/yy (FITS standard before 2000)
2048 dd-mm-yy (nonstandard use before 2000)
2049 yyyy-mm-dd (FITS standard after 1999)
2050 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2051 {
2052 double dj0, dj, tsec, dt;
2053
2054 dj0 = fd2jd (string);
2055 dt = utdt (dj0);
2056 dj = dj0 - (dt / 86400.0);
2057 dt = utdt (dj);
2058 tsec = fd2ts (string);
2059 tsec = tsec - dt + 32.184;
2060 return (ts2fd (tsec));
2061 }
2062
2063
2064 /* FD2TAI-- convert from UT in FITS format to TAI in FITS format */
2065
2066 char *
fd2tai(string)2067 fd2tai (string)
2068
2069 char *string; /* FITS date string, which may be:
2070 fractional year
2071 dd/mm/yy (FITS standard before 2000)
2072 dd-mm-yy (nonstandard use before 2000)
2073 yyyy-mm-dd (FITS standard after 1999)
2074 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2075 {
2076 double dj, tsec, dt;
2077
2078 dj = fd2jd (string);
2079 dt = utdt (dj);
2080 tsec = fd2ts (string);
2081 tsec = tsec + dt - 32.184;
2082 return (ts2fd (tsec));
2083 }
2084
2085
2086 /* DT2TAI-- convert from UT as yyyy.mmdd hh.mmssss to TAI in same format */
2087
2088 void
dt2tai(date,time)2089 dt2tai (date, time)
2090 double *date; /* Date as yyyy.mmdd */
2091 double *time; /* Time as hh.mmssxxxx
2092 *if time<0, it is time as -(fraction of a day) */
2093 {
2094 double dj, dt, tsec;
2095
2096 dj = dt2jd (*date, *time);
2097 dt = utdt (dj);
2098 tsec = dt2ts (*date, *time);
2099 tsec = tsec + dt - 32.184;
2100 ts2dt (tsec, date, time);
2101 return;
2102 }
2103
2104
2105 /* TAI2DT-- convert from TAI as yyyy.mmdd hh.mmssss to UT in same format */
2106
2107 void
tai2dt(date,time)2108 tai2dt (date, time)
2109 double *date; /* Date as yyyy.mmdd */
2110 double *time; /* Time as hh.mmssxxxx
2111 *if time<0, it is time as -(fraction of a day) */
2112 {
2113 double dj, dt, tsec, tsec0;
2114
2115 dj = dt2jd (*date, *time);
2116 dt = utdt (dj);
2117 tsec0 = dt2ts (*date, *time);
2118 tsec = tsec0 + dt;
2119 dj = ts2jd (tsec);
2120 dt = utdt (dj);
2121 tsec = tsec0 + dt + 32.184;
2122 ts2dt (tsec, date, time);
2123 return;
2124 }
2125
2126
2127 /* ET2FD-- convert from ET (or TDT or TT) in FITS format to UT in FITS format */
2128
2129 char *
et2fd(string)2130 et2fd (string)
2131
2132 char *string; /* FITS date string, which may be:
2133 fractional year
2134 dd/mm/yy (FITS standard before 2000)
2135 dd-mm-yy (nonstandard use before 2000)
2136 yyyy-mm-dd (FITS standard after 1999)
2137 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2138 {
2139 double dj0, dj, tsec, dt;
2140
2141 dj0 = fd2jd (string);
2142 dt = utdt (dj0);
2143 dj = dj0 - (dt / 86400.0);
2144 dt = utdt (dj);
2145 tsec = fd2ts (string);
2146 tsec = tsec - dt;
2147 return (ts2fd (tsec));
2148 }
2149
2150
2151 /* FD2ET-- convert from UT in FITS format to ET (or TDT or TT) in FITS format */
2152
2153 char *
fd2et(string)2154 fd2et (string)
2155
2156 char *string; /* FITS date string, which may be:
2157 fractional year
2158 dd/mm/yy (FITS standard before 2000)
2159 dd-mm-yy (nonstandard use before 2000)
2160 yyyy-mm-dd (FITS standard after 1999)
2161 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2162 {
2163 double dj, tsec, dt;
2164
2165 dj = fd2jd (string);
2166 dt = utdt (dj);
2167 tsec = fd2ts (string);
2168 tsec = tsec + dt;
2169 return (ts2fd (tsec));
2170 }
2171
2172
2173 /* DT2ET-- convert from UT as yyyy.mmdd hh.mmssss to ET in same format */
2174
2175 void
dt2et(date,time)2176 dt2et (date, time)
2177 double *date; /* Date as yyyy.mmdd */
2178 double *time; /* Time as hh.mmssxxxx
2179 *if time<0, it is time as -(fraction of a day) */
2180 {
2181 double dj, dt, tsec;
2182
2183 dj = dt2jd (*date, *time);
2184 dt = utdt (dj);
2185 tsec = dt2ts (*date, *time);
2186 tsec = tsec + dt;
2187 ts2dt (tsec, date, time);
2188 return;
2189 }
2190
2191
2192 /* EDT2DT-- convert from ET as yyyy.mmdd hh.mmssss to UT in same format */
2193
2194 void
edt2dt(date,time)2195 edt2dt (date, time)
2196 double *date; /* Date as yyyy.mmdd */
2197 double *time; /* Time as hh.mmssxxxx
2198 *if time<0, it is time as -(fraction of a day) */
2199 {
2200 double dj, dt, tsec, tsec0;
2201
2202 dj = dt2jd (*date, *time);
2203 dt = utdt (dj);
2204 tsec0 = dt2ts (*date, *time);
2205 tsec = tsec0 + dt;
2206 dj = ts2jd (tsec);
2207 dt = utdt (dj);
2208 tsec = tsec0 + dt;
2209 ts2dt (tsec, date, time);
2210 return;
2211 }
2212
2213
2214 /* JD2JED-- convert from Julian Date to Julian Ephemeris Date */
2215
2216 double
jd2jed(dj)2217 jd2jed (dj)
2218
2219 double dj; /* Julian Date */
2220 {
2221 double dt;
2222
2223 dt = utdt (dj);
2224 return (dj + (dt / 86400.0));
2225 }
2226
2227
2228 /* JED2JD-- convert from Julian Ephemeris Date to Julian Date */
2229
2230 double
jed2jd(dj)2231 jed2jd (dj)
2232
2233 double dj; /* Julian Ephemeris Date */
2234 {
2235 double dj0, dt;
2236
2237 dj0 = dj;
2238 dt = utdt (dj);
2239 dj = dj0 - (dt / 86400.0);
2240 dt = utdt (dj);
2241 return (dj - (dt / 86400.0));
2242 }
2243
2244
2245 /* TS2ETS-- convert from UT in seconds since 1950-01-01 to ET in same format */
2246
2247 double
ts2ets(tsec)2248 ts2ets (tsec)
2249
2250 double tsec;
2251 {
2252 double dj, dt;
2253
2254 dj = ts2jd (tsec);
2255 dt = utdt (dj);
2256 return (tsec + dt);
2257 }
2258
2259
2260 /* ETS2TS-- convert from ET in seconds since 1950-01-01 to UT in same format */
2261
2262 double
ets2ts(tsec)2263 ets2ts (tsec)
2264
2265 double tsec;
2266 {
2267 double dj, dj0, dt;
2268
2269 dj0 = ts2jd (tsec);
2270 dt = utdt (dj0);
2271 dj = dj0 - (dt / 86400.0);
2272 dt = utdt (dj);
2273 return (tsec - dt);
2274 }
2275
2276
2277 /* UTDT-- Compute difference between UT and dynamical time (ET-UT) */
2278
2279 double
utdt(dj)2280 utdt (dj)
2281
2282 double dj; /* Julian Date (UT) */
2283 {
2284 double dt, date, time, ts, ts1, ts0, date0, yfrac, diff, cj;
2285 int i, iyr, iyear;
2286
2287 /* If after 1972-01-01, use tabulated TAI-UT */
2288 if (dj >= 2441317.5) {
2289 dt = 0.0;
2290 for (i = 22; i > 0; i--) {
2291 if (dj >= taijd[i])
2292 dt = taidt[i];
2293 }
2294 dt = dt + 32.184;
2295 }
2296
2297 /* For 1800-01-01 to 1972-01-01, use table of ET-UT from AE */
2298 else if (dj >= 2378496.5) {
2299 jd2dt (dj, &date, &time);
2300 ts = jd2ts (dj);
2301 iyear = (int) date;
2302 iyr = iyear - 1800;
2303 date0 = (double) iyear + 0.0101;
2304 ts0 = dt2ts (date0, 0.0);
2305 date0 = (double) (iyear + 1) + 0.0101;
2306 ts1 = dt2ts (date0, 0.0);
2307 yfrac = (ts - ts0) / (ts1 - ts0);
2308 diff = dttab[iyr+1] - dttab[iyr];
2309 dt = dttab[iyr] + (diff * yfrac);
2310 }
2311
2312 /* Compute back to 1600 using formula from McCarthy and Babcock (1986) */
2313 else if (dj >= 2305447.5) {
2314 cj = (dj - 2378496.5) / 36525.0;
2315 dt = 5.156 + 13.3066 * (cj - 0.19) * (cj - 0.19);
2316 }
2317
2318 /* Compute back to 948 using formula from Stephenson and Morrison (1984) */
2319 else if (dj >= 2067309.5) {
2320 cj = (dj - 2378496.5) / 36525.0;
2321 dt = 25.5 * cj * cj;
2322 }
2323
2324 /*Compute back to 390 BC using formula from Stephenson and Morrison (1984)*/
2325 else if (dj >= 0.0) {
2326 cj = (dj = 2378496.5) / 36525.0;
2327 dt = 1360.0 + (320.0 * cj) + (44.3 * cj * cj);
2328 }
2329
2330 else
2331 dt = 0.0;
2332 return (dt);
2333 }
2334
2335
2336 /* FD2OFD-- convert any FITS standard datetime string to old FITS standard date */
2337
2338 char *
fd2ofd(string)2339 fd2ofd (string)
2340
2341 char *string; /* FITS date string, which may be:
2342 fractional year
2343 dd/mm/yy (FITS standard before 2000)
2344 dd-mm-yy (nonstandard use before 2000)
2345 yyyy-mm-dd (FITS standard after 1999)
2346 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2347 {
2348 int iyr,imon,iday,ihr,imn;
2349 double sec;
2350 char *dstr;
2351
2352 fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 3);
2353
2354 /* Convert to old FITS date format */
2355 dstr = (char *) calloc (32, sizeof (char));
2356 if (iyr < 1900)
2357 sprintf (dstr, "*** date out of range ***");
2358 else if (iyr < 2000)
2359 sprintf (dstr, "%02d/%02d/%02d", iday, imon, iyr-1900);
2360 else if (iyr < 2900.0)
2361 sprintf (dstr, "%02d/%02d/%3d", iday, imon, iyr-1900);
2362 else
2363 sprintf (dstr, "*** date out of range ***");
2364 return (dstr);
2365 }
2366
2367
2368 /* FD2OFT-- convert any FITS standard date to old FITS standard time */
2369
2370 char *
fd2oft(string)2371 fd2oft (string)
2372
2373 char *string; /* FITS date string, which may be:
2374 fractional year
2375 dd/mm/yy (FITS standard before 2000)
2376 dd-mm-yy (nonstandard use before 2000)
2377 yyyy-mm-dd (FITS standard after 1999)
2378 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2379 {
2380 int iyr,imon,iday,ihr,imn;
2381 double sec;
2382 char *tstr;
2383
2384 fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 3);
2385
2386 /* Convert to old FITS date format */
2387 tstr = (char *) calloc (32, sizeof (char));
2388 sprintf (tstr, "%02d:%02d:%06.3f", ihr, imn, sec);
2389 return (tstr);
2390 }
2391
2392
2393 /* FD2DT-- convert FITS standard date to date, time as yyyy.mmdd hh.mmsss */
2394
2395 void
fd2dt(string,date,time)2396 fd2dt (string, date, time)
2397
2398 char *string; /* FITS date string, which may be:
2399 fractional year
2400 dd/mm/yy (FITS standard before 2000)
2401 dd-mm-yy (nonstandard use before 2000)
2402 yyyy-mm-dd (FITS standard after 1999)
2403 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2404 double *date; /* Date as yyyy.mmdd (returned)
2405 yyyy = calendar year (e.g. 1973)
2406 mm = calendar month (e.g. 04 = april)
2407 dd = calendar day (e.g. 15) */
2408 double *time; /* Time as hh.mmssxxxx (returned)
2409 *if time<0, it is time as -(fraction of a day)
2410 hh = hour of day (0 .le. hh .le. 23)
2411 nn = minutes (0 .le. nn .le. 59)
2412 ss = seconds (0 .le. ss .le. 59)
2413 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
2414 {
2415 int iyr,imon,iday,ihr,imn;
2416 double sec;
2417
2418 fd2i (string,&iyr,&imon,&iday,&ihr,&imn,&sec, 4);
2419
2420 /* Convert date to yyyy.mmdd */
2421 if (iyr < 0) {
2422 *date = (double) (-iyr) + 0.01 * (double) imon + 0.0001 * (double) iday;
2423 *date = -(*date);
2424 }
2425 else
2426 *date = (double) iyr + 0.01 * (double) imon + 0.0001 * (double) iday;
2427
2428 /* Convert time to hh.mmssssss */
2429 *time = (double) ihr + 0.01 * (double) imn + 0.0001 * sec;
2430
2431 return;
2432 }
2433
2434
2435 /* FD2EP-- convert from FITS standard date to fractional year */
2436
2437 double
fd2ep(string)2438 fd2ep (string)
2439
2440 char *string; /* FITS date string, which may be:
2441 yyyy.ffff (fractional year)
2442 dd/mm/yy (FITS standard before 2000)
2443 dd-mm-yy (nonstandard FITS use before 2000)
2444 yyyy-mm-dd (FITS standard after 1999)
2445 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2446
2447 {
2448 double dj; /* Julian date */
2449 dj = fd2jd (string);
2450 if (dj < 1.0)
2451 return (dj / 365.2422);
2452 else
2453 return (jd2ep (dj));
2454 }
2455
2456
2457 /* FD2EPB-- convert from FITS standard date to Besselian epoch */
2458
2459 double
fd2epb(string)2460 fd2epb (string)
2461
2462 char *string; /* FITS date string, which may be:
2463 yyyy.ffff (fractional year)
2464 dd/mm/yy (FITS standard before 2000)
2465 dd-mm-yy (nonstandard FITS use before 2000)
2466 yyyy-mm-dd (FITS standard after 1999)
2467 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2468
2469 {
2470 double dj; /* Julian date */
2471 dj = fd2jd (string);
2472 if (dj < 1.0)
2473 return (dj / 365.242198781);
2474 else
2475 return (jd2epb (dj));
2476 }
2477
2478
2479 /* FD2EPJ-- convert from FITS standard date to Julian epoch */
2480
2481 double
fd2epj(string)2482 fd2epj (string)
2483
2484 char *string; /* FITS date string, which may be:
2485 yyyy.ffff (fractional year)
2486 dd/mm/yy (FITS standard before 2000)
2487 dd-mm-yy (nonstandard FITS use before 2000)
2488 yyyy-mm-dd (FITS standard after 1999)
2489 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2490
2491 {
2492 double dj; /* Julian date */
2493 dj = fd2jd (string);
2494 if (dj < 1.0)
2495 return (dj / 365.25);
2496 else
2497 return (jd2epj (dj));
2498 }
2499
2500
2501 /* DT2TSU-- convert from date and time to Unix seconds since 1970-01-01T0:00 */
2502
2503 time_t
dt2tsu(date,time)2504 dt2tsu (date,time)
2505
2506 double date; /* Date as yyyy.mmdd */
2507 double time; /* Time as hh.mmssxxxx
2508 *if time<0, it is time as -(fraction of a day) */
2509 {
2510 return ((time_t)(dt2ts (date, time) - 631152000.0));
2511 }
2512
2513
2514 /* DT2TSI-- convert from date and time to IRAF seconds since 1980-01-01T0:00 */
2515
2516 int
dt2tsi(date,time)2517 dt2tsi (date,time)
2518
2519 double date; /* Date as yyyy.mmdd */
2520 double time; /* Time as hh.mmssxxxx
2521 *if time<0, it is time as -(fraction of a day) */
2522 {
2523 return ((int)(dt2ts (date, time) - 946684800.0));
2524 }
2525
2526
2527
2528 /* DT2TS-- convert from date, time as yyyy.mmdd hh.mmsss to sec since 1950.0 */
2529
2530 double
dt2ts(date,time)2531 dt2ts (date,time)
2532
2533 double date; /* Date as yyyy.mmdd
2534 yyyy = calendar year (e.g. 1973)
2535 mm = calendar month (e.g. 04 = april)
2536 dd = calendar day (e.g. 15) */
2537 double time; /* Time as hh.mmssxxxx
2538 *if time<0, it is time as -(fraction of a day)
2539 hh = hour of day (0 .le. hh .le. 23)
2540 nn = minutes (0 .le. nn .le. 59)
2541 ss = seconds (0 .le. ss .le. 59)
2542 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
2543 {
2544 double tsec; /* Seconds past 1950.0 (returned) */
2545
2546 double dh,dm,dd;
2547 int iy,im,id;
2548
2549 /* Calculate the number of full years, months, and days already
2550 * elapsed since 0h, March 1, -1 (up to most recent midnight). */
2551
2552 /* convert time of day to elapsed seconds */
2553
2554 /* If time is < 0, it is assumed to be a fractional day */
2555 if (time < 0.0)
2556 tsec = time * -86400.0;
2557 else {
2558 dh = (int) (time + 0.0000000001);
2559 dm = (int) (((time - dh) * 100.0) + 0.0000000001);
2560 tsec = (time * 10000.0) - (dh * 10000.0) - (dm * 100.0);
2561 tsec = (int) (tsec * 100000.0 + 0.0001) / 100000.0;
2562 tsec = tsec + (dm * 60.0) + (dh * 3600.0);
2563 }
2564
2565
2566 /* Calculate the number of full months elapsed since
2567 * the current or most recent March */
2568 if (date >= 0.0301) {
2569 iy = (int) (date + 0.0000000001);
2570 im = (int) (((date - (double) (iy)) * 10000.0) + 0.00000001);
2571 id = im % 100;
2572 im = (im / 100) + 9;
2573 if (im < 12) iy = iy - 1;
2574 im = im % 12;
2575 id = id - 1;
2576
2577 /* starting with March as month 0 and ending with the following
2578 * February as month 11, the calculation of the number of days
2579 * per month reduces to a simple formula. the following statement
2580 * determines the number of whole days elapsed since 3/1/-1 and then
2581 * subtracts the 712163 days between then and 1/1/1950. it converts
2582 * the result to seconds and adds the accumulated seconds above. */
2583 id = id + ((im+1+im/6+im/11)/2 * 31) + ((im-im/6-im/11)/2 * 30) +
2584 (iy / 4) - (iy / 100) + (iy / 400);
2585 dd = (double) id + (365.0 * (double) iy) - 712163.0;
2586 tsec = tsec + (dd * 86400.0);
2587 }
2588
2589 return (tsec);
2590 }
2591
2592
2593 /* TS2DT-- convert seconds since 1950.0 to date, time as yyyy.mmdd hh.mmssss */
2594
2595 void
ts2dt(tsec,date,time)2596 ts2dt (tsec,date,time)
2597
2598 double tsec; /* Seconds past 1950.0 */
2599 double *date; /* Date as yyyy.mmdd (returned)
2600 yyyy = calendar year (e.g. 1973)
2601 mm = calendar month (e.g. 04 = april)
2602 dd = calendar day (e.g. 15) */
2603 double *time; /* Time as hh.mmssxxxx (returned)
2604 *if time<0, it is time as -(fraction of a day)
2605 hh = hour of day (0 .le. hh .le. 23)
2606 nn = minutes (0 .le. nn .le. 59)
2607 ss = seconds (0 .le. ss .le. 59)
2608 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
2609 {
2610 int iyr,imon,iday,ihr,imn;
2611 double sec;
2612
2613 ts2i (tsec,&iyr,&imon,&iday,&ihr,&imn,&sec, 4);
2614
2615 /* Convert date to yyyy.mmdd */
2616 if (iyr < 0) {
2617 *date = (double) (-iyr) + 0.01 * (double) imon + 0.0001 * (double) iday;
2618 *date = -(*date);
2619 }
2620 else
2621 *date = (double) iyr + 0.01 * (double) imon + 0.0001 * (double) iday;
2622
2623 /* Convert time to hh.mmssssss */
2624 *time = (double) ihr + 0.01 * (double) imn + 0.0001 * sec;
2625
2626 return;
2627 }
2628
2629
2630 /* TSI2DT-- Convert seconds since 1980-01-01 to date yyyy.ddmm, time hh.mmsss */
2631
2632 void
tsi2dt(isec,date,time)2633 tsi2dt (isec,date,time)
2634
2635 int isec; /* Seconds past 1980-01-01 */
2636 double *date; /* Date as yyyy.mmdd (returned) */
2637 double *time; /* Time as hh.mmssxxxx (returned) */
2638 {
2639 ts2dt (tsi2ts (isec), date, time);
2640 }
2641
2642
2643 /* TSI2FD-- Convert seconds since 1980-01-01 to FITS standard date string */
2644
2645 char *
tsi2fd(isec)2646 tsi2fd (isec)
2647
2648 int isec; /* Seconds past 1980-01-01 */
2649 {
2650 return (ts2fd (tsi2ts (isec)));
2651 }
2652
2653
2654 /* TSI2TS-- Convert seconds since 1980-01-01 to seconds since 1950-01-01 */
2655
2656 double
tsi2ts(isec)2657 tsi2ts (isec)
2658 int isec; /* Seconds past 1980-01-01 */
2659 {
2660 return ((double) isec + 946684800.0);
2661 }
2662
2663
2664 /* TSU2FD-- Convert seconds since 1970-01-01 to FITS standard date string */
2665
2666 char *
tsu2fd(isec)2667 tsu2fd (isec)
2668 time_t isec; /* Seconds past 1970-01-01 */
2669 {
2670 return (ts2fd (tsu2ts (isec)));
2671 }
2672
2673
2674 /* TSU2DT-- Convert seconds since 1970-01-01 to date yyyy.ddmm, time hh.mmsss */
2675
2676 void
tsu2dt(isec,date,time)2677 tsu2dt (isec,date,time)
2678 time_t isec; /* Seconds past 1970-01-01 */
2679 double *date; /* Date as yyyy.mmdd (returned) */
2680 double *time; /* Time as hh.mmssxxxx (returned) */
2681 {
2682 ts2dt (tsu2ts (isec), date, time);
2683 }
2684
2685
2686 /* TSU2TS-- Convert seconds since 1970-01-01 to seconds since 1950-01-01 */
2687
2688 double
tsu2ts(isec)2689 tsu2ts (isec)
2690 time_t isec; /* Seconds past 1970-01-01 */
2691 {
2692 return ((double) isec + 631152000.0);
2693 }
2694
2695 /* TSU2TSI-- UT seconds since 1970-01-01 to local seconds since 1980-01-01 */
2696
2697 int
tsu2tsi(isec)2698 tsu2tsi (isec)
2699 time_t isec; /* Seconds past 1970-01-01 */
2700 {
2701 double date, time;
2702 struct tm *ts;
2703
2704 /* Get local time from UT seconds */
2705 ts = localtime (&isec);
2706 if (ts->tm_year < 1000)
2707 date = (double) (ts->tm_year + 1900);
2708 else
2709 date = (double) ts->tm_year;
2710 date = date + (0.01 * (double) (ts->tm_mon + 1));
2711 date = date + (0.0001 * (double) ts->tm_mday);
2712 time = (double) ts->tm_hour;
2713 time = time + (0.01 * (double) ts->tm_min);
2714 time = time + (0.0001 * (double) ts->tm_sec);
2715 return ((int)(dt2ts (date, time) - 631152000.0));
2716 }
2717
2718
2719 /* TS2FD-- convert seconds since 1950.0 to FITS date, yyyy-mm-ddThh:mm:ss.ss */
2720
2721 char *
ts2fd(tsec)2722 ts2fd (tsec)
2723
2724 double tsec; /* Seconds past 1950.0 */
2725 {
2726 double date, time;
2727
2728 ts2dt (tsec, &date, &time);
2729 return (dt2fd (date, time));
2730 }
2731
2732
2733 /* TSD2FD-- convert seconds since start of day to FITS time, hh:mm:ss.ss */
2734
2735 char *
tsd2fd(tsec)2736 tsd2fd (tsec)
2737
2738 double tsec; /* Seconds since start of day */
2739 {
2740 double date, time;
2741 char *thms, *fdate;
2742 int lfd, nbc;
2743
2744 ts2dt (tsec, &date, &time);
2745 fdate = dt2fd (date, time);
2746 thms = (char *) calloc (16, 1);
2747 lfd = strlen (fdate);
2748 nbc = lfd - 11;
2749 strncpy (thms, fdate+11, nbc);
2750 return (thms);
2751 }
2752
2753
2754 /* TSD2DT-- convert seconds since start of day to hh.mmssss */
2755
2756 double
tsd2dt(tsec)2757 tsd2dt (tsec)
2758
2759 double tsec; /* Seconds since start of day */
2760 {
2761 double date, time;
2762
2763 ts2dt (tsec, &date, &time);
2764 return (time);
2765 }
2766
2767
2768
2769 /* DT2I-- convert vigesimal date and time to year month day hours min sec */
2770
2771 void
dt2i(date,time,iyr,imon,iday,ihr,imn,sec,ndsec)2772 dt2i (date, time, iyr, imon, iday, ihr, imn, sec, ndsec)
2773
2774 double date; /* Date as yyyy.mmdd (returned)
2775 yyyy = calendar year (e.g. 1973)
2776 mm = calendar month (e.g. 04 = april)
2777 dd = calendar day (e.g. 15) */
2778 double time; /* Time as hh.mmssxxxx (returned)
2779 *if time<0, it is time as -(fraction of a day)
2780 hh = hour of day (0 .le. hh .le. 23)
2781 nn = minutes (0 .le. nn .le. 59)
2782 ss = seconds (0 .le. ss .le. 59)
2783 xxxx = tenths of milliseconds (0 .le. xxxx .le. 9999) */
2784 int *iyr; /* year (returned) */
2785 int *imon; /* month (returned) */
2786 int *iday; /* day (returned) */
2787 int *ihr; /* hours (returned) */
2788 int *imn; /* minutes (returned) */
2789 double *sec; /* seconds (returned) */
2790 int ndsec; /* Number of decimal places in seconds (0=int) */
2791
2792 {
2793 double t,d;
2794
2795 t = time;
2796 if (date < 0.0)
2797 d = -date;
2798 else
2799 d = date;
2800
2801 /* Extract components of time */
2802 *ihr = dint (t + 0.000000001);
2803 t = 100.0 * (t - (double) *ihr);
2804 *imn = dint (t + 0.0000001);
2805 *sec = 100.0 * (t - (double) *imn);
2806
2807 /* Extract components of date */
2808 *iyr = dint (d + 0.00001);
2809 d = 100.0 * (d - (double) *iyr);
2810 if (date < 0.0)
2811 *iyr = - *iyr;
2812 *imon = dint (d + 0.001);
2813 d = 100.0 * (d - (double) *imon);
2814 *iday = dint (d + 0.1);
2815
2816 /* Make sure date and time are legal */
2817 fixdate (iyr, imon, iday, ihr, imn, sec, ndsec);
2818
2819 return;
2820 }
2821
2822
2823 /* FD2I-- convert from FITS standard date to year, mon, day, hours, min, sec */
2824
2825 void
fd2i(string,iyr,imon,iday,ihr,imn,sec,ndsec)2826 fd2i (string, iyr, imon, iday, ihr, imn, sec, ndsec)
2827
2828 char *string; /* FITS date string, which may be:
2829 yyyy.ffff (fractional year)
2830 dd/mm/yy (FITS standard before 2000)
2831 dd-mm-yy (nonstandard FITS use before 2000)
2832 yyyy-mm-dd (FITS standard after 1999)
2833 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
2834 int *iyr; /* year (returned) */
2835 int *imon; /* month (returned) */
2836 int *iday; /* day (returned) */
2837 int *ihr; /* hours (returned) */
2838 int *imn; /* minutes (returned) */
2839 double *sec; /* seconds (returned) */
2840 int ndsec; /* Number of decimal places in seconds (0=int) */
2841
2842 {
2843 double tsec, fday, hr, mn;
2844 int i;
2845 char *sstr, *dstr, *tstr, *cstr, *nval, *fstr;
2846
2847 /* Initialize all returned data to zero */
2848 *iyr = 0;
2849 *imon = 0;
2850 *iday = 0;
2851 *ihr = 0;
2852 *imn = 0;
2853 *sec = 0.0;
2854
2855 /* Return if no input string */
2856 if (string == NULL)
2857 return;
2858
2859 /* Check for various non-numeric characters */
2860 sstr = strchr (string,'/');
2861 dstr = strchr (string,'-');
2862 if (dstr == string)
2863 dstr = strchr (string+1, '-');
2864 fstr = strchr (string, '.');
2865 tstr = strchr (string,'T');
2866 if (tstr == NULL)
2867 tstr = strchr (string, 'Z');
2868 if (tstr == NULL)
2869 tstr = strchr (string, 'S');
2870 if (fstr != NULL && tstr != NULL && fstr > tstr)
2871 fstr = NULL;
2872 cstr = strchr (string,':');
2873
2874 /* Original FITS date format: dd/mm/yy */
2875 if (sstr > string) {
2876 *sstr = '\0';
2877 *iday = (int) atof (string);
2878 if (*iday > 31) {
2879 *iyr = *iday;
2880 if (*iyr >= 0 && *iyr <= 49)
2881 *iyr = *iyr + 2000;
2882 else if (*iyr < 1000)
2883 *iyr = *iyr + 1900;
2884 *sstr = '/';
2885 nval = sstr + 1;
2886 sstr = strchr (nval,'/');
2887 if (sstr > string) {
2888 *sstr = '\0';
2889 *imon = (int) atof (nval);
2890 *sstr = '/';
2891 nval = sstr + 1;
2892 *iday = (int) atof (nval);
2893 }
2894 }
2895 else {
2896 *sstr = '/';
2897 nval = sstr + 1;
2898 sstr = strchr (nval,'/');
2899 if (sstr == NULL)
2900 sstr = strchr (nval,'-');
2901 if (sstr > string) {
2902 *sstr = '\0';
2903 *imon = (int) atof (nval);
2904 *sstr = '/';
2905 nval = sstr + 1;
2906 *iyr = (int) atof (nval);
2907 if (*iyr >= 0 && *iyr <= 49)
2908 *iyr = *iyr + 2000;
2909 else if (*iyr < 1000)
2910 *iyr = *iyr + 1900;
2911 }
2912 }
2913 tstr = strchr (string,'_');
2914 if (tstr == NULL)
2915 return;
2916 }
2917
2918 /* New FITS date format: yyyy-mm-ddThh:mm:ss[.sss] */
2919 else if (dstr > string) {
2920 *dstr = '\0';
2921 *iyr = (int) atof (string);
2922 *dstr = '-';
2923 nval = dstr + 1;
2924 dstr = strchr (nval,'-');
2925 *imon = 1;
2926 *iday = 1;
2927
2928 /* Decode year, month, and day */
2929 if (dstr > string) {
2930 *dstr = '\0';
2931 *imon = (int) atof (nval);
2932 *dstr = '-';
2933 nval = dstr + 1;
2934 if (tstr > string)
2935 *tstr = '\0';
2936 *iday = (int) atof (nval);
2937
2938 /* If fraction of a day is present, turn it into a time */
2939 if (fstr != NULL) {
2940 fday = atof (fstr);
2941 hr = fday * 24.0;
2942 *ihr = (int) hr;
2943 mn = 60.0 * (hr - (double) *ihr);
2944 *imn = (int) mn;
2945 *sec = 60.0 * (mn - (double) *imn);
2946 }
2947
2948 if (tstr > string)
2949 *tstr = 'T';
2950 }
2951
2952 /* If date is > 31, it is really year in old format */
2953 if (*iday > 31) {
2954 i = *iyr;
2955 if (*iday < 100)
2956 *iyr = *iday + 1900;
2957 else
2958 *iyr = *iday;
2959 *iday = i;
2960 }
2961 }
2962
2963 /* In rare cases, a FITS time is entered as an epoch */
2964 else if (tstr == NULL && cstr == NULL && isnum (string)) {
2965 tsec = ep2ts (atof (string));
2966 ts2i (tsec,iyr,imon,iday,ihr,imn,sec, ndsec);
2967 return;
2968 }
2969
2970 /* Extract time, if it is present */
2971 if (tstr > string || cstr > string) {
2972 if (tstr > string)
2973 nval = tstr + 1;
2974 else
2975 nval = string;
2976 cstr = strchr (nval,':');
2977 if (cstr > string) {
2978 *cstr = '\0';
2979 *ihr = (int) atof (nval);
2980 *cstr = ':';
2981 nval = cstr + 1;
2982 cstr = strchr (nval,':');
2983 if (cstr > string) {
2984 *cstr = '\0';
2985 *imn = (int) atof (nval);
2986 *cstr = ':';
2987 nval = cstr + 1;
2988 *sec = atof (nval);
2989 }
2990 else
2991 *imn = (int) atof (nval);
2992 }
2993 else
2994 *ihr = (int) atof (nval);
2995 }
2996 else
2997 ndsec = -1;
2998
2999 /* Make sure date and time are legal */
3000 fixdate (iyr, imon, iday, ihr, imn, sec, ndsec);
3001
3002 return;
3003 }
3004
3005
3006 /* TS2I-- convert sec since 1950.0 to year month day hours minutes seconds */
3007
3008 void
ts2i(tsec,iyr,imon,iday,ihr,imn,sec,ndsec)3009 ts2i (tsec,iyr,imon,iday,ihr,imn,sec, ndsec)
3010
3011 double tsec; /* seconds since 1/1/1950 0:00 */
3012 int *iyr; /* year (returned) */
3013 int *imon; /* month (returned) */
3014 int *iday; /* day (returned) */
3015 int *ihr; /* hours (returned) */
3016 int *imn; /* minutes (returned) */
3017 double *sec; /* seconds (returned) */
3018 int ndsec; /* Number of decimal places in seconds (0=int) */
3019
3020 {
3021 double t,days, ts, dts;
3022 int nc,nc4,nly,ny,m,im;
3023
3024 /* Round seconds to 0 - 4 decimal places */
3025 ts = tsec + 61530883200.0;
3026 if (ts < 0.0)
3027 dts = -0.5;
3028 else
3029 dts = 0.5;
3030 if (ndsec < 1)
3031 t = dint (ts + dts) * 10000.0;
3032 else if (ndsec < 2)
3033 t = dint (ts * 10.0 + dts) * 1000.0;
3034 else if (ndsec < 3)
3035 t = dint (ts * 100.0 + dts) * 100.0;
3036 else if (ndsec < 4)
3037 t = dint (ts * 1000.0 + dts) * 10.0;
3038 else
3039 t = dint (ts * 10000.0 + dts);
3040 ts = t / 10000.0;
3041
3042 /* Time of day (hours, minutes, seconds */
3043 *ihr = (int) (dmod (ts/3600.0, 24.0));
3044 *imn = (int) (dmod (ts/60.0, 60.0));
3045 *sec = dmod (ts, 60.0);
3046
3047 /* Number of days since 0 hr 0/0/0000 */
3048 days = dint ((t / 864000000.0) + 0.000001);
3049
3050 /* Number of leap centuries (400 years) */
3051 nc4 = (int) ((days / 146097.0) + 0.00001);
3052
3053 /* Number of centuries since last /400 */
3054 days = days - (146097.0 * (double) (nc4));
3055 nc = (int) ((days / 36524.0) + 0.000001);
3056 if (nc > 3) nc = 3;
3057
3058 /* Number of leap years since last century */
3059 days = days - (36524.0 * nc);
3060 nly = (int) ((days / 1461.0) + 0.0000000001);
3061
3062 /* Number of years since last leap year */
3063 days = days - (1461.0 * (double) nly);
3064 ny = (int) ((days / 365.0) + 0.00000001);
3065 if (ny > 3) ny = 3;
3066
3067 /* Day of month */
3068 days = days - (365.0 * (double) ny);
3069 if (days < 0) {
3070 m = 0;
3071 *iday = 29;
3072 }
3073 else {
3074 *iday = (int) (days + 0.00000001) + 1;
3075 for (m = 1; m <= 12; m++) {
3076 im = (m + ((m - 1) / 5)) % 2;
3077 /* fprintf (stderr,"%d %d %d %d\n", m, im, *iday, nc); */
3078 if (*iday-1 < im+30) break;
3079 *iday = *iday - im - 30;
3080 }
3081 }
3082
3083 /* Month */
3084 *imon = ((m+1) % 12) + 1;
3085
3086 /* Year */
3087 *iyr = nc4*400 + nc*100 + nly*4 + ny + m/11;
3088
3089 /* Make sure date and time are legal */
3090 fixdate (iyr, imon, iday, ihr, imn, sec, ndsec);
3091
3092 return;
3093 }
3094
3095
3096 /* UT2DOY-- Current Universal Time as year, day of year */
3097
3098 void
ut2doy(year,doy)3099 ut2doy (year, doy)
3100
3101 int *year; /* Year (returned) */
3102 double *doy; /* Day of year (returned) */
3103 {
3104 double date, time;
3105 ut2dt (&date, &time);
3106 dt2doy (date, time, year, doy);
3107 return;
3108 }
3109
3110
3111 /* UT2DT-- Current Universal Time as date (yyyy.mmdd) and time (hh.mmsss) */
3112
3113 void
ut2dt(date,time)3114 ut2dt(date, time)
3115
3116 double *date; /* Date as yyyy.mmdd (returned) */
3117 double *time; /* Time as hh.mmssxxxx (returned) */
3118 {
3119 time_t tsec;
3120 struct timeval tp;
3121 struct timezone tzp;
3122 struct tm *ts;
3123
3124 gettimeofday (&tp,&tzp);
3125
3126 tsec = tp.tv_sec;
3127 ts = gmtime (&tsec);
3128
3129 if (ts->tm_year < 1000)
3130 *date = (double) (ts->tm_year + 1900);
3131 else
3132 *date = (double) ts->tm_year;
3133 *date = *date + (0.01 * (double) (ts->tm_mon + 1));
3134 *date = *date + (0.0001 * (double) ts->tm_mday);
3135 *time = (double) ts->tm_hour;
3136 *time = *time + (0.01 * (double) ts->tm_min);
3137 *time = *time + (0.0001 * (double) ts->tm_sec);
3138
3139 return;
3140 }
3141
3142
3143 /* UT2EP-- Return current Universal Time as fractional year */
3144
3145 double
ut2ep()3146 ut2ep()
3147 {
3148 return (jd2ep (ut2jd()));
3149 }
3150
3151
3152 /* UT2EPB-- Return current Universal Time as Besselian epoch */
3153
3154 double
ut2epb()3155 ut2epb()
3156 {
3157 return (jd2epb (ut2jd()));
3158 }
3159
3160
3161 /* UT2EPJ-- Return current Universal Time as Julian epoch */
3162
3163 double
ut2epj()3164 ut2epj()
3165 {
3166 return (jd2epj (ut2jd()));
3167 }
3168
3169
3170 /* UT2FD-- Return current Universal Time as FITS ISO date string */
3171
3172 char *
ut2fd()3173 ut2fd()
3174 {
3175 int year, month, day, hour, minute, second;
3176 time_t tsec;
3177 struct timeval tp;
3178 struct timezone tzp;
3179 struct tm *ts;
3180 char *isotime;
3181
3182 gettimeofday (&tp,&tzp);
3183 tsec = tp.tv_sec;
3184 ts = gmtime (&tsec);
3185
3186 year = ts->tm_year;
3187 if (year < 1000)
3188 year = year + 1900;
3189 month = ts->tm_mon + 1;
3190 day = ts->tm_mday;
3191 hour = ts->tm_hour;
3192 minute = ts->tm_min;
3193 second = ts->tm_sec;
3194
3195 isotime = (char *) calloc (32, sizeof (char));
3196 sprintf (isotime, "%04d-%02d-%02dT%02d:%02d:%02d",
3197 year, month, day, hour, minute, second);
3198 return (isotime);
3199 }
3200
3201
3202 /* UT2JD-- Return current Universal Time as Julian Date */
3203
3204 double
ut2jd()3205 ut2jd()
3206 {
3207 return (fd2jd (ut2fd()));
3208 }
3209
3210
3211 /* UT2MJD-- convert current UT to Modified Julian Date */
3212
3213 double
ut2mjd()3214 ut2mjd ()
3215
3216 {
3217 return (ut2jd() - 2400000.5);
3218 }
3219
3220 /* UT2TS-- current Universal Time as IRAF seconds since 1950-01-01T00:00 */
3221
3222 double
ut2ts()3223 ut2ts()
3224 {
3225 double tsec;
3226 char *datestring;
3227 datestring = ut2fd();
3228 tsec = fd2ts (datestring);
3229 free (datestring);
3230 return (tsec);
3231 }
3232
3233
3234 /* UT2TSI-- current Universal Time as IRAF seconds since 1980-01-01T00:00 */
3235
3236 int
ut2tsi()3237 ut2tsi()
3238 {
3239 return ((int)(ut2ts() - 946684800.0));
3240 }
3241
3242
3243 /* UT2TSU-- current Universal Time as IRAF seconds since 1970-01-01T00:00 */
3244
3245 time_t
ut2tsu()3246 ut2tsu()
3247 {
3248 return ((time_t)(ut2ts () - 631152000.0));
3249 }
3250
3251
3252 /* FD2GST-- convert from FITS date to Greenwich Sidereal Time */
3253
3254 char *
fd2gst(string)3255 fd2gst (string)
3256
3257 char *string; /* FITS date string, which may be:
3258 fractional year
3259 dd/mm/yy (FITS standard before 2000)
3260 dd-mm-yy (nonstandard use before 2000)
3261 yyyy-mm-dd (FITS standard after 1999)
3262 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
3263 {
3264 double dj, gsec, date, time;
3265
3266 dj = fd2jd (string);
3267 gsec = jd2gst (dj);
3268 ts2dt (gsec, &date, &time);
3269 date = 0.0;
3270 return (dt2fd (date, time));
3271 }
3272
3273
3274 /* DT2GST-- convert from UT as yyyy.mmdd hh.mmssss to Greenwich Sidereal Time*/
3275
3276 void
dt2gst(date,time)3277 dt2gst (date, time)
3278 double *date; /* Date as yyyy.mmdd */
3279 double *time; /* Time as hh.mmssxxxx
3280 *if time<0, it is time as -(fraction of a day) */
3281 {
3282 double dj, gsec;
3283
3284 dj = dt2ts (*date, *time);
3285 gsec = jd2gst (dj);
3286 ts2dt (gsec, date, time);
3287 *date = 0.0;
3288 return;
3289 }
3290
3291
3292 /* JD2LST - Local Sidereal Time in seconds from Julian Date */
3293
3294 double
jd2lst(dj)3295 jd2lst (dj)
3296
3297 double dj; /* Julian Date */
3298 {
3299 double gst, lst;
3300
3301 /* Compute Greenwich Sidereal Time at this epoch */
3302 gst = jd2gst (dj);
3303
3304 /* Subtract longitude (degrees to seconds of time) */
3305 lst = gst - (240.0 * longitude);
3306 if (lst < 0.0)
3307 lst = lst + 86400.0;
3308 else if (lst > 86400.0)
3309 lst = lst - 86400.0;
3310 return (lst);
3311 }
3312
3313
3314 /* FD2LST - Local Sidereal Time as hh:mm:ss.ss
3315 from Universal Time as FITS ISO date */
3316
3317 char *
fd2lst(string)3318 fd2lst (string)
3319
3320 char *string; /* FITS date string, which may be:
3321 fractional year
3322 dd/mm/yy (FITS standard before 2000)
3323 dd-mm-yy (nonstandard use before 2000)
3324 yyyy-mm-dd (FITS standard after 1999) */
3325 {
3326 double dj, date, time, lst;
3327
3328 dj = fd2jd (string);
3329 lst = jd2lst (dj);
3330 ts2dt (lst, &date, &time);
3331 date = 0.0;
3332 return (dt2fd (date, time));
3333 }
3334
3335
3336 /* DT2LST - Local Sidereal Time as hh.mmssss
3337 from Universal Time as yyyy.mmdd hh.mmssss */
3338
3339 void
dt2lst(date,time)3340 dt2lst (date, time)
3341
3342 double *date; /* Date as yyyy.mmdd */
3343 double *time; /* Time as hh.mmssxxxx
3344 *if time<0, it is time as -(fraction of a day) */
3345 {
3346 double dj, lst, date0;
3347
3348 dj = dt2jd (*date, *time);
3349 lst = jd2lst (dj);
3350 date0 = 0.0;
3351 ts2dt (lst, &date0, time);
3352 return;
3353 }
3354
3355
3356 /* TS2LST - Local Sidereal Time in seconds of day
3357 * from Universal Time in seconds since 1951-01-01T0:00:00
3358 */
3359
3360 double
ts2lst(tsec)3361 ts2lst (tsec)
3362
3363 double tsec; /* time since 1950.0 in UT seconds */
3364 {
3365 double gst; /* Greenwich Sidereal Time in seconds since 0:00 */
3366 double lst; /* Local Sidereal Time in seconds since 0:00 */
3367 double gsec, date;
3368
3369 /* Greenwich Sidereal Time */
3370 gsec = ts2gst (tsec);
3371 date = 0.0;
3372 ts2dt (gsec, &date, &gst);
3373
3374 lst = gst - (longitude / 15.0);
3375 if (lst < 0.0)
3376 lst = lst + 86400.0;
3377 else if (lst > 86400.0)
3378 lst = lst - 86400.0;
3379 return (lst);
3380 }
3381
3382
3383 /* LST2FD - calculate current UT given Local Sidereal Time
3384 * plus date in FITS ISO format (yyyy-mm-dd)
3385 * Return UT date and time in FITS ISO format
3386 */
3387
3388 char *
lst2fd(string)3389 lst2fd (string)
3390
3391 char *string; /* UT Date, LST as yyyy-mm-ddShh:mm:ss.ss */
3392 {
3393 double sdj, dj;
3394
3395 sdj = fd2jd (string);
3396
3397 dj = lst2jd (sdj);
3398
3399 return (jd2fd (dj));
3400 }
3401
3402
3403 /* LST2JD - calculate current Julian Date given Local Sidereal Time
3404 * plus current Julian Date (0.5 at 0:00 UT)
3405 * Return UT date and time as Julian Date
3406 */
3407
3408 double
lst2jd(sdj)3409 lst2jd (sdj)
3410
3411 double sdj; /* Julian Date of desired day at 0:00 UT + sidereal time */
3412 {
3413 double gst; /* Greenwich Sidereal Time in seconds since 0:00 */
3414 double lsd; /* Local Sidereal Time in seconds since 0:00 */
3415 double gst0, tsd, dj1, dj0, eqnx;
3416 int idj;
3417
3418 /* Julian date at 0:00 UT */
3419 idj = (int) sdj;
3420 dj0 = (double) idj + 0.5;
3421 if (dj0 > sdj) dj0 = dj0 - 1.0;
3422
3423 /* Greenwich Sidereal Time at 0:00 UT in seconds */
3424 gst0 = jd2gst (dj0);
3425
3426 /* Sidereal seconds since 0:00 */
3427 lsd = (sdj - dj0) * 86400.0;
3428
3429 /* Remove longitude for current Greenwich Sidereal Time in seconds */
3430 /* (convert longitude from degrees to seconds of time) */
3431 gst = lsd + (longitude * 240.0);
3432
3433 /* Time since 0:00 UT */
3434 tsd = (gst - gst0) / 1.0027379093;
3435
3436 /* Julian Date (UT) */
3437 dj1 = dj0 + (tsd / 86400.0);
3438
3439 /* Equation of the equinoxes converted to UT seconds */
3440 eqnx = eqeqnx (dj1) / 1.002739093;
3441
3442 /* Remove equation of equinoxes */
3443 dj1 = dj1 - (eqnx / 86400.0);
3444 if (dj1 < dj0)
3445 dj1 = dj1 + 1.0;
3446
3447 return (dj1);
3448 }
3449
3450
3451 /* MST2FD - calculate current UT given Greenwich Mean Sidereal Time
3452 * plus date in FITS ISO format (yyyy-mm-ddShh:mm:ss.ss)
3453 * Return UT date and time in FITS ISO format
3454 */
3455
3456 char *
mst2fd(string)3457 mst2fd (string)
3458
3459 char *string; /* UT Date, MST as yyyy-mm-ddShh:mm:ss.ss */
3460 {
3461 double sdj, dj;
3462
3463 sdj = fd2jd (string);
3464
3465 dj = mst2jd (sdj);
3466
3467 return (jd2fd (dj));
3468 }
3469
3470
3471 /* MST2JD - calculate current UT given Greenwich Mean Sidereal Time
3472 * plus date in Julian Date (0:00 UT + Mean Sidereal Time)
3473 * Return UT date and time as Julian Date
3474 */
3475
3476 double
mst2jd(sdj)3477 mst2jd (sdj)
3478
3479 double sdj; /* UT Date, MST as Julian Date */
3480 {
3481 double tsd, djd, st0, dj0, dj;
3482
3483 dj0 = (double) ((int) sdj) + 0.5;
3484
3485 /* Greenwich Mean Sidereal Time at 0:00 UT in seconds */
3486 st0 = jd2mst (dj0);
3487
3488 /* Mean Sidereal Time in seconds */
3489 tsd = (sdj - dj0) * 86400.0;
3490 if (tsd < 0.0)
3491 tsd = tsd + 86400.0;
3492
3493 /* Convert to fraction of a day since 0:00 UT */
3494 djd = ((tsd - st0) / 1.0027379093) / 86400.0;
3495
3496 /* Julian Date */
3497 dj = dj0 + djd;
3498 if (dj < dj0)
3499 dj = dj + (1.0 / 1.0027379093);
3500
3501 return (dj);
3502 }
3503
3504
3505
3506 /* GST2FD - calculate current UT given Greenwich Sidereal Time
3507 * plus date in FITS ISO format (yyyy-mm-ddShh:mm:ss.ss)
3508 * Return UT date and time in FITS ISO format
3509 */
3510
3511 char *
gst2fd(string)3512 gst2fd (string)
3513
3514 char *string; /* UT Date, GST as yyyy-mm-ddShh:mm:ss.ss */
3515 {
3516 double sdj, dj;
3517
3518 sdj = fd2jd (string);
3519
3520 dj = gst2jd (sdj);
3521
3522 return (jd2fd (dj));
3523 }
3524
3525
3526 /* GST2JD - calculate current UT given Greenwich Sidereal Time
3527 * plus date as Julian Date (JD at 0:00 UT + sidereal time)
3528 * Return UT date and time as Julian Date
3529 */
3530
3531 double
gst2jd(sdj)3532 gst2jd (sdj)
3533
3534 double sdj; /* UT Date, GST as Julian Date */
3535 {
3536 double dj, tsd, djd, st0, dj0, eqnx;
3537
3538 dj0 = (double) ((int) sdj) + 0.5;
3539
3540 /* Greenwich Mean Sidereal Time at 0:00 UT in seconds */
3541 st0 = jd2mst (dj0);
3542
3543 /* Mean Sidereal Time in seconds */
3544 tsd = (sdj - dj0) * 86400.0;
3545 if (tsd < 0.0)
3546 tsd = tsd + 86400.0;
3547
3548 /* Convert to fraction of a day since 0:00 UT */
3549 djd = ((tsd - st0) / 1.0027379093) / 86400.0;
3550
3551 /* Julian Date */
3552 dj = dj0 + djd;
3553
3554 /* Equation of the equinoxes (converted to UT seconds) */
3555 eqnx = eqeqnx (dj) / 1.002737909;
3556
3557 dj = dj - eqnx / 86400.0;
3558 if (dj < dj0)
3559 dj = dj + 1.0;
3560
3561 return (dj);
3562 }
3563
3564
3565 /* LST2DT - calculate current UT given Local Sidereal Time as hh.mmsss
3566 * plus date as yyyy.mmdd
3567 * Return UT time as hh.mmssss
3568 */
3569
3570 double
lst2dt(date0,time0)3571 lst2dt (date0, time0)
3572
3573 double date0; /* UT date as yyyy.mmdd */
3574 double time0; /* LST as hh.mmssss */
3575 {
3576 double gst; /* Greenwich Sidereal Time in seconds since 0:00 */
3577 double lst; /* Local Sidereal Time in seconds since 0:00 */
3578 double date1; /* UT date as yyyy.mmdd */
3579 double time1; /* UT as hh.mmssss */
3580 double tsec0, gst0, tsd, tsec;
3581
3582 /* Greenwich Sidereal Time at 0:00 UT */
3583 tsec0 = dt2ts (date0, 0.0);
3584 gst0 = ts2gst (tsec0);
3585
3586 /* Current Greenwich Sidereal Time in seconds */
3587 /* (convert longitude from degrees to seconds of time) */
3588 lst = dt2ts (0.0, time0);
3589 gst = lst + (longitude * 240.0);
3590
3591 /* Time since 0:00 UT */
3592 tsd = (gst - gst0) / 1.0027379093;
3593
3594 /* UT date and time */
3595 tsec = tsec0 + tsd;
3596 ts2dt (tsec, &date1, &time1);
3597
3598 return (time1);
3599 }
3600
3601
3602 /* TS2GST - calculate Greenwich Sidereal Time given Universal Time
3603 * in seconds since 1951-01-01T0:00:00
3604 * Return sidereal time of day in seconds
3605 */
3606
3607 double
ts2gst(tsec)3608 ts2gst (tsec)
3609
3610 double tsec; /* time since 1950.0 in UT seconds */
3611 {
3612 double gst; /* Greenwich Sidereal Time in seconds since 0:00 */
3613 double tsd, eqnx, dj;
3614 int its;
3615
3616 /* Elapsed time as of 0:00 UT */
3617 if (tsec >= 0.0) {
3618 its = (int) (tsec + 0.5);
3619 tsd = (double) (its % 86400);
3620 }
3621 else {
3622 its = (int) (-tsec + 0.5);
3623 tsd = (double) (86400 - (its % 86400));
3624 }
3625
3626 /* Mean sidereal time */
3627 gst = ts2mst (tsec);
3628
3629 /* Equation of the equinoxes */
3630 dj = ts2jd (tsec);
3631 eqnx = eqeqnx (dj);
3632
3633 /* Apparent sidereal time at 0:00 ut */
3634 gst = gst + eqnx;
3635
3636 /* Current sidereal time */
3637 gst = gst + (tsd * 1.0027379093);
3638 gst = dmod (gst,86400.0);
3639
3640 return (gst);
3641 }
3642
3643
3644 /* FD2MST-- convert from FITS date Mean Sidereal Time */
3645
3646 char *
fd2mst(string)3647 fd2mst (string)
3648
3649 char *string; /* FITS date string, which may be:
3650 fractional year
3651 dd/mm/yy (FITS standard before 2000)
3652 dd-mm-yy (nonstandard use before 2000)
3653 yyyy-mm-dd (FITS standard after 1999)
3654 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
3655 {
3656 double gsec, date, time, dj;
3657
3658 dj = fd2jd (string);
3659 gsec = jd2mst (dj);
3660 ts2dt (gsec, &date, &time);
3661 date = 0.0;
3662 return (dt2fd (date, time));
3663 }
3664
3665
3666 /* DT2MST-- convert from UT as yyyy.mmdd hh.mmssss to Mean Sidereal Time
3667 in the same format */
3668
3669 void
dt2mst(date,time)3670 dt2mst (date, time)
3671 double *date; /* Date as yyyy.mmdd */
3672 double *time; /* Time as hh.mmssxxxx
3673 *if time<0, it is time as -(fraction of a day) */
3674 {
3675 double date0, gsec, dj;
3676 date0 = *date;
3677 dj = dt2jd (*date, *time);
3678 gsec = jd2mst (dj);
3679 ts2dt (gsec, date, time);
3680 *date = date0;
3681 return;
3682 }
3683
3684
3685 /* TS2MST - calculate Greenwich Mean Sidereal Time given Universal Time
3686 * in seconds since 1951-01-01T0:00:00
3687 */
3688
3689 double
ts2mst(tsec)3690 ts2mst (tsec)
3691
3692 double tsec; /* time since 1950.0 in UT seconds */
3693 {
3694 double dj;
3695
3696 dj = ts2jd (tsec);
3697 return (jd2mst (dj));
3698 }
3699
3700
3701 /* JD2MST - Julian Date to Greenwich Mean Sidereal Time using IAU 2000
3702 * Return sideral time in seconds of time
3703 * (from USNO NOVAS package
3704 * http://aa.usno.navy.mil/software/novas/novas_info.html
3705 */
3706
3707 double
jd2mst2(dj)3708 jd2mst2 (dj)
3709
3710 double dj; /* Julian Date */
3711 {
3712 double dt, t, t2, t3, mst, st;
3713
3714 dt = dj - 2451545.0;
3715 t = dt / 36525.0;
3716 t2 = t * t;
3717 t3 = t2 * t;
3718
3719 /* Compute Greenwich Mean Sidereal Time in seconds */
3720 st = (8640184.812866 * t) + (3155760000.0 * t) - (0.0000062 * t3)
3721 + (0.093104 * t2) + 67310.54841;
3722
3723 mst = dmod (st, 86400.0);
3724 if (mst < 0.0)
3725 mst = mst + 86400.0;
3726 return (mst);
3727 }
3728
3729
3730 /* MJD2MST - Modified Julian Date to Greenwich Mean Sidereal Time using IAU 2000
3731 * Return sideral time in seconds of time
3732 * (from USNO NOVAS package
3733 * http://aa.usno.navy.mil/software/novas/novas_info.html
3734 */
3735
3736 double
mjd2mst(dj)3737 mjd2mst (dj)
3738
3739 double dj; /* Modified Julian Date */
3740 {
3741 double dt, t, t2, t3, mst, st;
3742
3743 dt = dj - 51544.5;
3744 t = dt / 36525.0;
3745 t2 = t * t;
3746 t3 = t2 * t;
3747
3748 /* Compute Greenwich Mean Sidereal Time in seconds */
3749 st = (8640184.812866 * t) + (3155760000.0 * t) - (0.0000062 * t3)
3750 + (0.093104 * t2) + 67310.54841;
3751
3752 mst = dmod (st, 86400.0);
3753 if (mst < 0.0)
3754 mst = mst + 86400.0;
3755 return (mst);
3756 }
3757
3758
3759 /* JD2GST - Julian Date to Greenwich Sideral Time
3760 * Return sideral time in seconds of time
3761 * (Jean Meeus, Astronomical Algorithms, Willmann-Bell, 1991, pp 83-84)
3762 */
3763
3764 double
jd2gst(dj)3765 jd2gst (dj)
3766
3767 double dj; /* Julian Date */
3768 {
3769 double dj0, gmt, gst, tsd, eqnx, ssd, l0;
3770 double ts2ss = 1.00273790935;
3771 int ijd;
3772
3773 /* Julian date at 0:00 UT */
3774 ijd = (int) dj;
3775 dj0 = (double) ijd + 0.5;
3776 if (dj0 > dj) dj0 = dj0 - 1.0;
3777
3778 /* Greenwich mean sidereal time at 0:00 UT in seconds */
3779 l0 = longitude;
3780 longitude = 0.0;
3781 gmt = jd2mst (dj0);
3782 longitude = l0;
3783
3784 /* Equation of the equinoxes */
3785 eqnx = eqeqnx (dj);
3786
3787 /* Apparent sidereal time at 0:00 ut */
3788 gst = gmt + eqnx;
3789
3790 /* UT seconds since 0:00 */
3791 tsd = (dj - dj0) * 86400.0;
3792 ssd = tsd * ts2ss;
3793
3794 /* Current sidereal time */
3795 gst = gst + ssd;
3796 gst = dmod (gst, 86400.0);
3797
3798 return (gst);
3799 }
3800
3801
3802 /* EQEQNX - Compute equation of the equinoxes for apparent sidereal time */
3803
3804 double
eqeqnx(dj)3805 eqeqnx (dj)
3806
3807 double dj; /* Julian Date */
3808
3809 {
3810 double dt, edj, dpsi, deps, obl, eqnx;
3811 double rad2tsec = 13750.98708;
3812
3813 /* Convert UT to Ephemeris Time (TDB or TT)*/
3814 dt = utdt (dj);
3815 edj = dj + dt / 86400.0;
3816
3817 /* Nutation and obliquity */
3818 compnut (edj, &dpsi, &deps, &obl);
3819
3820 /* Correct obliquity for nutation */
3821 obl = obl + deps;
3822
3823 /* Equation of the equinoxes in seconds */
3824 eqnx = (dpsi * cos (obl)) * rad2tsec;
3825
3826 return (eqnx);
3827 }
3828
3829
3830
3831 /* JD2MST - Julian Date to Mean Sideral Time
3832 * Return sideral time in seconds of time
3833 * (Jean Meeus, Astronomical Algorithms, Willmann-Bell, 1991, pp 83-84)
3834 */
3835
3836 double
jd2mst(dj)3837 jd2mst (dj)
3838
3839 double dj; /* Julian Date */
3840 {
3841 double dt, t, mst;
3842
3843 dt = dj - 2451545.0;
3844 t = dt / 36525.0;
3845
3846 /* Compute Greenwich mean sidereal time in degrees (Meeus, page 84) */
3847 mst = 280.46061837 + (360.98564736629 * dt) + (0.000387933 * t * t) -
3848 (t * t * t / 38710000.0);
3849
3850 /* Keep degrees between 0 and 360 */
3851 while (mst > 360.0)
3852 mst = mst - 360.0;
3853 while (mst < 0.0)
3854 mst = mst + 360.0;
3855
3856 /* Convert to time in seconds (3600 / 15) */
3857 mst = mst * 240.0;
3858
3859 /* Subtract longitude (degrees to seconds of time) */
3860 mst = mst - (240.0 * longitude);
3861 if (mst < 0.0)
3862 mst = mst + 86400.0;
3863 else if (mst > 86400.0)
3864 mst = mst - 86400.0;
3865
3866 return (mst);
3867 }
3868
3869
3870 /* COMPNUT - Compute nutation using the IAU 2000b model */
3871 /* Translated from Pat Wallace's Fortran subroutine iau_nut00b (June 26 2007)
3872 into C by Jessica Mink on September 5, 2008 */
3873
3874 #define NLS 77 /* number of terms in the luni-solar nutation model */
3875
3876 void
compnut(dj,dpsi,deps,eps0)3877 compnut (dj, dpsi, deps, eps0)
3878
3879 double dj; /* Julian Date */
3880 double *dpsi; /* Nutation in longitude in radians (returned) */
3881 double *deps; /* Nutation in obliquity in radians (returned) */
3882 double *eps0; /* Mean obliquity in radians (returned) */
3883
3884 /* This routine is translated from the International Astronomical Union's
3885 * Fortran SOFA (Standards Of Fundamental Astronomy) software collection.
3886 *
3887 * notes:
3888 *
3889 * 1) the nutation components in longitude and obliquity are in radians
3890 * and with respect to the equinox and ecliptic of date. the
3891 * obliquity at j2000 is assumed to be the lieske et al. (1977) value
3892 * of 84381.448 arcsec. (the errors that result from using this
3893 * routine with the iau 2006 value of 84381.406 arcsec can be
3894 * neglected.)
3895 *
3896 * the nutation model consists only of luni-solar terms, but includes
3897 * also a fixed offset which compensates for certain long-period
3898 * planetary terms (note 7).
3899 *
3900 * 2) this routine is an implementation of the iau 2000b abridged
3901 * nutation model formally adopted by the iau general assembly in
3902 * 2000. the routine computes the mhb_2000_short luni-solar nutation
3903 * series (luzum 2001), but without the associated corrections for
3904 * the precession rate adjustments and the offset between the gcrs
3905 * and j2000 mean poles.
3906 *
3907 * 3) the full IAU 2000a (mhb2000) nutation model contains nearly 1400
3908 * terms. the IAU 2000b model (mccarthy & luzum 2003) contains only
3909 * 77 terms, plus additional simplifications, yet still delivers
3910 * results of 1 mas accuracy at present epochs. this combination of
3911 * accuracy and size makes the IAU 2000b abridged nutation model
3912 * suitable for most practical applications.
3913 *
3914 * the routine delivers a pole accurate to 1 mas from 1900 to 2100
3915 * (usually better than 1 mas, very occasionally just outside 1 mas).
3916 * the full IAU 2000a model, which is implemented in the routine
3917 * iau_nut00a (q.v.), delivers considerably greater accuracy at
3918 * current epochs; however, to realize this improved accuracy,
3919 * corrections for the essentially unpredictable free-core-nutation
3920 * (fcn) must also be included.
3921 *
3922 * 4) the present routine provides classical nutation. the
3923 * mhb_2000_short algorithm, from which it is adapted, deals also
3924 * with (i) the offsets between the gcrs and mean poles and (ii) the
3925 * adjustments in longitude and obliquity due to the changed
3926 * precession rates. these additional functions, namely frame bias
3927 * and precession adjustments, are supported by the sofa routines
3928 * iau_bi00 and iau_pr00.
3929 *
3930 * 6) the mhb_2000_short algorithm also provides "total" nutations,
3931 * comprising the arithmetic sum of the frame bias, precession
3932 * adjustments, and nutation (luni-solar + planetary). these total
3933 * nutations can be used in combination with an existing IAU 1976
3934 * precession implementation, such as iau_pmat76, to deliver gcrs-to-
3935 * true predictions of mas accuracy at current epochs. however, for
3936 * symmetry with the iau_nut00a routine (q.v. for the reasons), the
3937 * sofa routines do not generate the "total nutations" directly.
3938 * should they be required, they could of course easily be generated
3939 * by calling iau_bi00, iau_pr00 and the present routine and adding
3940 * the results.
3941 *
3942 * 7) the IAU 2000b model includes "planetary bias" terms that are fixed
3943 * in size but compensate for long-period nutations. the amplitudes
3944 * quoted in mccarthy & luzum (2003), namely dpsi = -1.5835 mas and
3945 * depsilon = +1.6339 mas, are optimized for the "total nutations"
3946 * method described in note 6. the luzum (2001) values used in this
3947 * sofa implementation, namely -0.135 mas and +0.388 mas, are
3948 * optimized for the "rigorous" method, where frame bias, precession
3949 * and nutation are applied separately and in that order. during the
3950 * interval 1995-2050, the sofa implementation delivers a maximum
3951 * error of 1.001 mas (not including fcn).
3952 *
3953 * References from original Fortran subroutines:
3954 *
3955 * Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351
3956 *
3957 * Lieske, J.H., Lederle, T., Fricke, W., Morando, B., "Expressions
3958 * for the precession quantities based upon the IAU 1976 system of
3959 * astronomical constants", Astron.Astrophys. 58, 1-2, 1-16. (1977)
3960 *
3961 * Luzum, B., private communication, 2001 (Fortran code
3962 * mhb_2000_short)
3963 *
3964 * McCarthy, D.D. & Luzum, B.J., "An abridged model of the
3965 * precession-nutation of the celestial pole", Cel.Mech.Dyn.Astron.
3966 * 85, 37-49 (2003)
3967 *
3968 * Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
3969 * Francou, G., Laskar, J., Astron.Astrophys. 282, 663-683 (1994)
3970 *
3971 */
3972
3973 {
3974 double as2r = 0.000004848136811095359935899141; /* arcseconds to radians */
3975
3976 double dmas2r = as2r / 1000.0; /* milliarcseconds to radians */
3977
3978 double as2pi = 1296000.0; /* arc seconds in a full circle */
3979
3980 double d2pi = 6.283185307179586476925287; /* 2pi */
3981
3982 double u2r = as2r / 10000000.0; /* units of 0.1 microarcsecond to radians */
3983
3984 double dj0 = 2451545.0; /* reference epoch (j2000), jd */
3985
3986 double djc = 36525.0; /* Days per julian century */
3987
3988 /* Miscellaneous */
3989 double t, el, elp, f, d, om, arg, dp, de, sarg, carg;
3990 double dpsils, depsls, dpsipl, depspl;
3991 int i, j;
3992
3993 int nls = NLS; /* number of terms in the luni-solar nutation model */
3994
3995 /* Fixed offset in lieu of planetary terms (radians) */
3996 double dpplan = - 0.135 * dmas2r;
3997 double deplan = + 0.388 * dmas2r;
3998
3999 /* Tables of argument and term coefficients */
4000
4001 /* Coefficients for fundamental arguments */
4002 /* Luni-solar argument multipliers: */
4003 /* l l' f d om */
4004 static int nals[5*NLS]=
4005 {0, 0, 0, 0, 1,
4006 0, 0, 2, -2, 2,
4007 0, 0, 2, 0, 2,
4008 0, 0, 0, 0, 2,
4009 0, 1, 0, 0, 0,
4010 0, 1, 2, -2, 2,
4011 1, 0, 0, 0, 0,
4012 0, 0, 2, 0, 1,
4013 1, 0, 2, 0, 2,
4014 0, -1, 2, -2, 2,
4015 0, 0, 2, -2, 1,
4016 -1, 0, 2, 0, 2,
4017 -1, 0, 0, 2, 0,
4018 1, 0, 0, 0, 1,
4019 -1, 0, 0, 0, 1,
4020 -1, 0, 2, 2, 2,
4021 1, 0, 2, 0, 1,
4022 -2, 0, 2, 0, 1,
4023 0, 0, 0, 2, 0,
4024 0, 0, 2, 2, 2,
4025 0, -2, 2, -2, 2,
4026 -2, 0, 0, 2, 0,
4027 2, 0, 2, 0, 2,
4028 1, 0, 2, -2, 2,
4029 -1, 0, 2, 0, 1,
4030 2, 0, 0, 0, 0,
4031 0, 0, 2, 0, 0,
4032 0, 1, 0, 0, 1,
4033 -1, 0, 0, 2, 1,
4034 0, 2, 2, -2, 2,
4035 0, 0, -2, 2, 0,
4036 1, 0, 0, -2, 1,
4037 0, -1, 0, 0, 1,
4038 -1, 0, 2, 2, 1,
4039 0, 2, 0, 0, 0,
4040 1, 0, 2, 2, 2,
4041 -2, 0, 2, 0, 0,
4042 0, 1, 2, 0, 2,
4043 0, 0, 2, 2, 1,
4044 0, -1, 2, 0, 2,
4045 0, 0, 0, 2, 1,
4046 1, 0, 2, -2, 1,
4047 2, 0, 2, -2, 2,
4048 -2, 0, 0, 2, 1,
4049 2, 0, 2, 0, 1,
4050 0, -1, 2, -2, 1,
4051 0, 0, 0, -2, 1,
4052 -1, -1, 0, 2, 0,
4053 2, 0, 0, -2, 1,
4054 1, 0, 0, 2, 0,
4055 0, 1, 2, -2, 1,
4056 1, -1, 0, 0, 0,
4057 -2, 0, 2, 0, 2,
4058 3, 0, 2, 0, 2,
4059 0, -1, 0, 2, 0,
4060 1, -1, 2, 0, 2,
4061 0, 0, 0, 1, 0,
4062 -1, -1, 2, 2, 2,
4063 -1, 0, 2, 0, 0,
4064 0, -1, 2, 2, 2,
4065 -2, 0, 0, 0, 1,
4066 1, 1, 2, 0, 2,
4067 2, 0, 0, 0, 1,
4068 -1, 1, 0, 1, 0,
4069 1, 1, 0, 0, 0,
4070 1, 0, 2, 0, 0,
4071 -1, 0, 2, -2, 1,
4072 1, 0, 0, 0, 2,
4073 -1, 0, 0, 1, 0,
4074 0, 0, 2, 1, 2,
4075 -1, 0, 2, 4, 2,
4076 -1, 1, 0, 1, 1,
4077 0, -2, 2, -2, 1,
4078 1, 0, 2, 2, 1,
4079 -2, 0, 2, 2, 2,
4080 -1, 0, 0, 0, 2,
4081 1, 1, 2, -2, 2};
4082
4083 /* Luni-solar nutation coefficients, in 1e-7 arcsec */
4084 /* longitude (sin, t*sin, cos), obliquity (cos, t*cos, sin) */
4085 static double cls[6*NLS]=
4086 {-172064161.0, -174666.0, 33386.0, 92052331.0, 9086.0, 15377.0,
4087 -13170906.0, -1675.0, -13696.0, 5730336.0, -3015.0, -4587.0,
4088 -2276413.0, -234.0, 2796.0, 978459.0, -485.0, 1374.0,
4089 2074554.0, 207.0, -698.0, -897492.0, 470.0, -291.0,
4090 1475877.0, -3633.0, 11817.0, 73871.0, -184.0, -1924.0,
4091 -516821.0, 1226.0, -524.0, 224386.0, -677.0, -174.0,
4092 711159.0, 73.0, -872.0, -6750.0, 0.0, 358.0,
4093 -387298.0, -367.0, 380.0, 200728.0, 18.0, 318.0,
4094 -301461.0, -36.0, 816.0, 129025.0, -63.0, 367.0,
4095 215829.0, -494.0, 111.0, -95929.0, 299.0, 132.0,
4096 128227.0, 137.0, 181.0, -68982.0, -9.0, 39.0,
4097 123457.0, 11.0, 19.0, -53311.0, 32.0, -4.0,
4098 156994.0, 10.0, -168.0, -1235.0, 0.0, 82.0,
4099 63110.0, 63.0, 27.0, -33228.0, 0.0, -9.0,
4100 -57976.0, -63.0, -189.0, 31429.0, 0.0, -75.0,
4101 -59641.0, -11.0, 149.0, 25543.0, -11.0, 66.0,
4102 -51613.0, -42.0, 129.0, 26366.0, 0.0, 78.0,
4103 45893.0, 50.0, 31.0, -24236.0, -10.0, 20.0,
4104 63384.0, 11.0, -150.0, -1220.0, 0.0, 29.0,
4105 -38571.0, -1.0, 158.0, 16452.0, -11.0, 68.0,
4106 32481.0, 0.0, 0.0, -13870.0, 0.0, 0.0,
4107 -47722.0, 0.0, -18.0, 477.0, 0.0, -25.0,
4108 -31046.0, -1.0, 131.0, 13238.0, -11.0, 59.0,
4109 28593.0, 0.0, -1.0, -12338.0, 10.0, -3.0,
4110 20441.0, 21.0, 10.0, -10758.0, 0.0, -3.0,
4111 29243.0, 0.0, -74.0, -609.0, 0.0, 13.0,
4112 25887.0, 0.0, -66.0, -550.0, 0.0, 11.0,
4113 -14053.0, -25.0, 79.0, 8551.0, -2.0, -45.0,
4114 15164.0, 10.0, 11.0, -8001.0, 0.0, -1.0,
4115 -15794.0, 72.0, -16.0, 6850.0, -42.0, -5.0,
4116 21783.0, 0.0, 13.0, -167.0, 0.0, 13.0,
4117 -12873.0, -10.0, -37.0, 6953.0, 0.0, -14.0,
4118 -12654.0, 11.0, 63.0, 6415.0, 0.0, 26.0,
4119 -10204.0, 0.0, 25.0, 5222.0, 0.0, 15.0,
4120 16707.0, -85.0, -10.0, 168.0, -1.0, 10.0,
4121 -7691.0, 0.0, 44.0, 3268.0, 0.0, 19.0,
4122 -11024.0, 0.0, -14.0, 104.0, 0.0, 2.0,
4123 7566.0, -21.0, -11.0, -3250.0, 0.0, -5.0,
4124 -6637.0, -11.0, 25.0, 3353.0, 0.0, 14.0,
4125 -7141.0, 21.0, 8.0, 3070.0, 0.0, 4.0,
4126 -6302.0, -11.0, 2.0, 3272.0, 0.0, 4.0,
4127 5800.0, 10.0, 2.0, -3045.0, 0.0, -1.0,
4128 6443.0, 0.0, -7.0, -2768.0, 0.0, -4.0,
4129 -5774.0, -11.0, -15.0, 3041.0, 0.0, -5.0,
4130 -5350.0, 0.0, 21.0, 2695.0, 0.0, 12.0,
4131 -4752.0, -11.0, -3.0, 2719.0, 0.0, -3.0,
4132 -4940.0, -11.0, -21.0, 2720.0, 0.0, -9.0,
4133 7350.0, 0.0, -8.0, -51.0, 0.0, 4.0,
4134 4065.0, 0.0, 6.0, -2206.0, 0.0, 1.0,
4135 6579.0, 0.0, -24.0, -199.0, 0.0, 2.0,
4136 3579.0, 0.0, 5.0, -1900.0, 0.0, 1.0,
4137 4725.0, 0.0, -6.0, -41.0, 0.0, 3.0,
4138 -3075.0, 0.0, -2.0, 1313.0, 0.0, -1.0,
4139 -2904.0, 0.0, 15.0, 1233.0, 0.0, 7.0,
4140 4348.0, 0.0, -10.0, -81.0, 0.0, 2.0,
4141 -2878.0, 0.0, 8.0, 1232.0, 0.0, 4.0,
4142 -4230.0, 0.0, 5.0, -20.0, 0.0, -2.0,
4143 -2819.0, 0.0, 7.0, 1207.0, 0.0, 3.0,
4144 -4056.0, 0.0, 5.0, 40.0, 0.0, -2.0,
4145 -2647.0, 0.0, 11.0, 1129.0, 0.0, 5.0,
4146 -2294.0, 0.0, -10.0, 1266.0, 0.0, -4.0,
4147 2481.0, 0.0, -7.0, -1062.0, 0.0, -3.0,
4148 2179.0, 0.0, -2.0, -1129.0, 0.0, -2.0,
4149 3276.0, 0.0, 1.0, -9.0, 0.0, 0.0,
4150 -3389.0, 0.0, 5.0, 35.0, 0.0, -2.0,
4151 3339.0, 0.0, -13.0, -107.0, 0.0, 1.0,
4152 -1987.0, 0.0, -6.0, 1073.0, 0.0, -2.0,
4153 -1981.0, 0.0, 0.0, 854.0, 0.0, 0.0,
4154 4026.0, 0.0, -353.0, -553.0, 0.0, -139.0,
4155 1660.0, 0.0, -5.0, -710.0, 0.0, -2.0,
4156 -1521.0, 0.0, 9.0, 647.0, 0.0, 4.0,
4157 1314.0, 0.0, 0.0, -700.0, 0.0, 0.0,
4158 -1283.0, 0.0, 0.0, 672.0, 0.0, 0.0,
4159 -1331.0, 0.0, 8.0, 663.0, 0.0, 4.0,
4160 1383.0, 0.0, -2.0, -594.0, 0.0, -2.0,
4161 1405.0, 0.0, 4.0, -610.0, 0.0, 2.0,
4162 1290.0, 0.0, 0.0, -556.0, 0.0, 0.0};
4163
4164 /* Interval between fundamental epoch J2000.0 and given date (JC) */
4165 t = (dj - dj0) / djc;
4166
4167 /* Luni-solar nutation */
4168
4169 /* Fundamental (delaunay) arguments from Simon et al. (1994) */
4170
4171 /* Mean anomaly of the moon */
4172 el = fmod (485868.249036 + (1717915923.2178 * t), as2pi) * as2r;
4173
4174 /* Mean anomaly of the sun */
4175 elp = fmod (1287104.79305 + (129596581.0481 * t), as2pi) * as2r;
4176
4177 /* Mean argument of the latitude of the moon */
4178 f = fmod (335779.526232 + (1739527262.8478 * t), as2pi) * as2r;
4179
4180 /* Mean elongation of the moon from the sun */
4181 d = fmod (1072260.70369 + (1602961601.2090 * t), as2pi ) * as2r;
4182
4183 /* Mean longitude of the ascending node of the moon */
4184 om = fmod (450160.398036 - (6962890.5431 * t), as2pi ) * as2r;
4185
4186 /* Initialize the nutation values */
4187 dp = 0.0;
4188 de = 0.0;
4189
4190 /* Summation of luni-solar nutation series (in reverse order) */
4191 for (i = nls; i > 0; i=i-1) {
4192 j = i - 1;
4193
4194 /* Argument and functions */
4195 arg = fmod ( (double) (nals[5*j]) * el +
4196 (double) (nals[1+5*j]) * elp +
4197 (double) (nals[2+5*j]) * f +
4198 (double) (nals[3+5*j]) * d +
4199 (double) (nals[4+5*j]) * om, d2pi);
4200 sarg = sin (arg);
4201 carg = cos (arg);
4202
4203 /* Terms */
4204 dp = dp + (cls[6*j] + cls[1+6*j] * t) * sarg + cls[2+6*j] * carg;
4205 de = de + (cls[3+6*j] + cls[4+6*j] * t) * carg + cls[5+6*j] * sarg;
4206 }
4207
4208 /* Convert from 0.1 microarcsec units to radians */
4209 dpsils = dp * u2r;
4210 depsls = de * u2r;
4211
4212 /* In lieu of planetary nutation */
4213
4214 /* Fixed offset to correct for missing terms in truncated series */
4215 dpsipl = dpplan;
4216 depspl = deplan;
4217
4218 /* Results */
4219
4220 /* Add luni-solar and planetary components */
4221 *dpsi = dpsils + dpsipl;
4222 *deps = depsls + depspl;
4223
4224 /* Mean Obliquity in radians (IAU 2006, Hilton, et al.) */
4225 *eps0 = ( 84381.406 +
4226 ( -46.836769 +
4227 ( -0.0001831 +
4228 ( 0.00200340 +
4229 ( -0.000000576 +
4230 ( -0.0000000434 ) * t ) * t ) * t ) * t ) * t ) * as2r;
4231 }
4232
4233
4234 /* ISDATE - Return 1 if string is an old or ISO FITS standard date */
4235
4236 int
isdate(string)4237 isdate (string)
4238
4239 char *string; /* Possible FITS date string, which may be:
4240 dd/mm/yy (FITS standard before 2000)
4241 dd-mm-yy (nonstandard FITS use before 2000)
4242 yyyy-mm-dd (FITS standard after 1999)
4243 yyyy-mm-ddThh:mm:ss.ss (FITS standard after 1999) */
4244
4245 {
4246 int iyr = 0; /* year (returned) */
4247 int imon = 0; /* month (returned) */
4248 int iday = 0; /* day (returned) */
4249 int i;
4250 char *sstr, *dstr, *tstr, *nval;
4251
4252 /* Translate string from ASCII to binary */
4253 if (string == NULL)
4254 return (0);
4255
4256 sstr = strchr (string,'/');
4257 dstr = strchr (string,'-');
4258 if (dstr == string)
4259 dstr = strchr (string+1,'-');
4260 tstr = strchr (string,'T');
4261
4262 /* Original FITS date format: dd/mm/yy */
4263 if (sstr > string) {
4264 *sstr = '\0';
4265 iday = (int) atof (string);
4266 *sstr = '/';
4267 nval = sstr + 1;
4268 sstr = strchr (nval,'/');
4269 if (sstr == NULL)
4270 sstr = strchr (nval,'-');
4271 if (sstr > string) {
4272 *sstr = '\0';
4273 imon = (int) atof (nval);
4274 *sstr = '/';
4275 nval = sstr + 1;
4276 iyr = (int) atof (nval);
4277 if (iyr < 1000)
4278 iyr = iyr + 1900;
4279 }
4280 if (imon > 0 && iday > 0)
4281 return (1);
4282 else
4283 return (0);
4284 }
4285
4286 /* New FITS date format: yyyy-mm-ddThh:mm:ss[.sss] */
4287 else if (dstr > string) {
4288 *dstr = '\0';
4289 iyr = (int) atof (string);
4290 nval = dstr + 1;
4291 *dstr = '-';
4292 dstr = strchr (nval,'-');
4293 imon = 0;
4294 iday = 0;
4295
4296 /* Decode year, month, and day */
4297 if (dstr > string) {
4298 *dstr = '\0';
4299 imon = (int) atof (nval);
4300 *dstr = '-';
4301 nval = dstr + 1;
4302 if (tstr > string)
4303 *tstr = '\0';
4304 iday = (int) atof (nval);
4305 if (tstr > string)
4306 *tstr = 'T';
4307 }
4308
4309 /* If day is > 31, it is really year in old format */
4310 if (iday > 31) {
4311 i = iyr;
4312 if (iday < 100)
4313 iyr = iday + 1900;
4314 else
4315 iyr = iday;
4316 iday = i;
4317 }
4318 if (imon > 0 && iday > 0)
4319 return (1);
4320 else
4321 return (0);
4322 }
4323
4324 /* If FITS date is entered as an epoch, return 0 anyway */
4325 else
4326 return (0);
4327 }
4328
4329
4330 /* Round seconds and make sure date and time numbers are within limits */
4331
4332 static void
fixdate(iyr,imon,iday,ihr,imn,sec,ndsec)4333 fixdate (iyr, imon, iday, ihr, imn, sec, ndsec)
4334
4335 int *iyr; /* year (returned) */
4336 int *imon; /* month (returned) */
4337 int *iday; /* day (returned) */
4338 int *ihr; /* hours (returned) */
4339 int *imn; /* minutes (returned) */
4340 double *sec; /* seconds (returned) */
4341 int ndsec; /* Number of decimal places in seconds (0=int) */
4342 {
4343 double days;
4344
4345 /* Round seconds to 0 - 4 decimal places (no rounding if <0, >4) */
4346 if (ndsec == 0)
4347 *sec = dint (*sec + 0.5);
4348 else if (ndsec < 2)
4349 *sec = dint (*sec * 10.0 + 0.5) / 10.0;
4350 else if (ndsec < 3)
4351 *sec = dint (*sec * 100.0 + 0.5) / 100.0;
4352 else if (ndsec < 4)
4353 *sec = dint (*sec * 1000.0 + 0.5) / 1000.0;
4354 else if (ndsec < 5)
4355 *sec = dint (*sec * 10000.0 + 0.5) / 10000.0;
4356
4357 /* Adjust minutes and hours */
4358 if (*sec > 60.0) {
4359 *sec = *sec - 60.0;
4360 *imn = *imn + 1;
4361 }
4362 if (*imn > 60) {
4363 *imn = *imn - 60;
4364 *ihr = *ihr + 1;
4365 }
4366
4367 /* Return if no date */
4368 if (*iyr == 0 && *imon == 0 && *iday == 0)
4369 return;
4370
4371 /* Adjust date */
4372 if (*ihr > 23) {
4373 *ihr = *ihr - 24;
4374 *iday = *iday + 1;
4375 }
4376 days = caldays (*iyr, *imon);
4377 if (*iday > days) {
4378 *iday = *iday - days;
4379 *imon = *imon + 1;
4380 }
4381 if (*iday < 1) {
4382 *imon = *imon - 1;
4383 if (*imon < 1) {
4384 *imon = *imon + 12;
4385 *iyr = *iyr - 1;
4386 }
4387 days = caldays (*iyr, *imon);
4388 *iday = *iday + days;
4389 }
4390 if (*imon < 1) {
4391 *imon = *imon + 12;
4392 *iyr = *iyr - 1;
4393 days = caldays (*iyr, *imon);
4394 if (*iday > days) {
4395 *iday = *iday - days;
4396 *imon = *imon + 1;
4397 }
4398 }
4399 if (*imon > 12) {
4400 *imon = *imon - 12;
4401 *iyr = *iyr + 1;
4402 }
4403 return;
4404 }
4405
4406
4407 /* Calculate days in month 1-12 given year (Gregorian calendar only) */
4408
4409 static int
caldays(year,month)4410 caldays (year, month)
4411
4412 int year; /* 4-digit year */
4413 int month; /* Month (1=January, 2=February, etc.) */
4414 {
4415 if (month < 1) {
4416 month = month + 12;
4417 year = year + 1;
4418 }
4419 if (month > 12) {
4420 month = month - 12;
4421 year = year + 1;
4422 }
4423 switch (month) {
4424 case 1:
4425 return (31);
4426 case 2:
4427 if (year%400 == 0)
4428 return (29);
4429 else if (year%100 == 0)
4430 return (28);
4431 else if (year%4 == 0)
4432 return (29);
4433 else
4434 return (28);
4435 case 3:
4436 return (31);
4437 case 4:
4438 return (30);
4439 case 5:
4440 return (31);
4441 case 6:
4442 return (30);
4443 case 7:
4444 return (31);
4445 case 8:
4446 return (31);
4447 case 9:
4448 return (30);
4449 case 10:
4450 return (31);
4451 case 11:
4452 return (30);
4453 case 12:
4454 return (31);
4455 default:
4456 return (0);
4457 }
4458 }
4459
4460
4461 static double
dint(dnum)4462 dint (dnum)
4463
4464 double dnum;
4465 {
4466 double dn;
4467
4468 if (dnum < 0.0)
4469 dn = -floor (-dnum);
4470 else
4471 dn = floor (dnum);
4472 return (dn);
4473 }
4474
4475
4476 static double
dmod(dnum,dm)4477 dmod (dnum, dm)
4478
4479 double dnum, dm;
4480 {
4481 double dnumx, dnumi, dnumf;
4482 if (dnum < 0.0)
4483 dnumx = -dnum;
4484 else
4485 dnumx = dnum;
4486 dnumi = dint (dnumx / dm);
4487 if (dnum < 0.0)
4488 dnumf = dnum + (dnumi * dm);
4489 else if (dnum > 0.0)
4490 dnumf = dnum - (dnumi * dm);
4491 else
4492 dnumf = 0.0;
4493 return (dnumf);
4494 }
4495
4496 /* Jul 1 1999 New file, based on iolib/jcon.f and iolib/vcon.f and hgetdate()
4497 * Oct 21 1999 Fix declarations after lint
4498 * Oct 27 1999 Fix bug to return epoch if fractional year input
4499 * Dec 9 1999 Fix bug in ts2jd() found by Pete Ratzlaff (SAO)
4500 * Dec 17 1999 Add all unimplemented conversions
4501 * Dec 20 1999 Add isdate(); leave date, time strings unchanged in fd2i()
4502 * Dec 20 1999 Make all fd2*() subroutines deal with time alone
4503 *
4504 * Jan 3 2000 In old FITS format, year 100 is assumed to be 2000
4505 * Jan 11 2000 Fix epoch to date conversion so .0 is 0:00, not 12:00
4506 * Jan 21 2000 Add separate Besselian and Julian epoch computations
4507 * Jan 28 2000 Add Modified Julian Date conversions
4508 * Mar 2 2000 Implement decimal places for FITS date string
4509 * Mar 14 2000 Fix bug in dealing with 2000-02-29 in ts2i()
4510 * Mar 22 2000 Add lt2* and ut2* to get current time as local and UT
4511 * Mar 24 2000 Fix calloc() calls
4512 * Mar 24 2000 Add tsi2* and tsu2* to convert IRAF and Unix seconds
4513 * May 1 2000 In old FITS format, all years < 1000 get 1900 added to them
4514 * Aug 1 2000 Make ep2jd and jd2ep consistently starting at 1/1 0:00
4515 *
4516 * Jan 11 2001 Print all messages to stderr
4517 * May 21 2001 Add day of year conversions
4518 * May 25 2001 Allow fraction of day in FITS date instead of time
4519 *
4520 * Apr 8 2002 Change all long declaration to time_t
4521 * May 13 2002 Fix bugs found by lint
4522 * Jul 5 2002 Fix bug in fixdate() so fractional seconds come out
4523 * Jul 8 2002 Fix rounding bug in t2i()
4524 * Jul 8 2002 Try Fliegel and Van Flandern's algorithm for JD to UT date
4525 * Jul 8 2002 If first character of string is -, check for other -'s in isdate
4526 * Sep 10 2002 Add ET/TDT/TT conversion from UT subroutines
4527 * Sep 10 2002 Add sidereal time conversions
4528 *
4529 * Jan 30 2003 Fix typo in ts2gst()
4530 * Mar 7 2003 Add conversions for heliocentric julian dates
4531 * May 20 2003 Declare nd in setdatedec()
4532 * Jul 18 2003 Add code to parse Las Campanas dates
4533 *
4534 * Mar 24 2004 If ndec > 0, add UT to FITS date even if it is 0:00:00
4535 *
4536 * Oct 14 2005 Add tsd2fd() and tsd2dt()
4537 *
4538 * May 3 2006 Drop declaration of unused variables
4539 * Jun 20 2006 Initialized uninitialized variables
4540 * Aug 2 2006 Add local sidereal time
4541 * Sep 13 2006 Add more local sidereal time subroutines
4542 * Oct 2 2006 Add UT to old FITS date conversions
4543 * Oct 6 2006 Add eqeqnx() to compute equation of the equinoxes
4544 *
4545 * Jan 8 2007 Remove unused variables
4546 *
4547 * Sep 5 2008 Replace nutation with IAU 2006 model translated from SOFA
4548 * Sep 9 2008 Add ang2hr(), ang2deg(), hr2ang(), deg2ang()
4549 * Sep 10 2008 Add longitude to mean standard time (default = Greenwich)
4550 * Oct 8 2008 Clean up sidereal time computations
4551 *
4552 * Sep 24 2009 Add end to comment "Coefficients for fundamental arguments"
4553 *
4554 * Jan 11 2012 Add TAI, TT, GPS time
4555 * Oct 19 2012 Unused l0 dropped from jd2lst(); ts2ss from jd2mst()
4556 *
4557 * May 2 2017 Allocate new output string for fd2ofd() and fd2oft()
4558 */
4559