LAL  7.1.7.1-2d066e5
XLALCivilTime.c
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2016 Karl Wette
3 * Copyright (C) 2007 Bernd Machenschalk, Jolien Creighton, Kipp Cannon
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with with program; see the file COPYING. If not, write to the
17 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18 * MA 02110-1301 USA
19 */
20 
21 #include <config.h>
22 #include <math.h>
23 #include <time.h>
24 #include <string.h>
25 #include <lal/Date.h>
26 #include <lal/LALAtomicDatatypes.h>
27 #include <lal/XLALError.h>
28 
29 #ifndef HAVE_GMTIME_R
30 #define gmtime_r(timep, result) memcpy((result), gmtime(timep), sizeof(struct tm))
31 #endif
32 
33 #include "XLALLeapSeconds.h" /* contains the leap second table */
34 
35 /**
36  * \defgroup XLALCivilTime_c CivilTime
37  * \ingroup Date_h
38  * \author Chin, D. W. and Creighton, J. D. E.
39  *
40  * \brief XLAL routines for converting civil time structures to GPS times.
41  *
42  * Civil time structures, represented in the C library by the \c struct \c tm
43  * structure, known as a broken down time structure, gives the time by normal
44  * measures (hours, minutes, seconds since midnight, day of the month, month of
45  * the year, and year). LAL uses seconds since the the GPS epoch as its
46  * reference. The GPS epoch is defined to be 0h UTC 6 Jan 1980 in civil time.
47  * The process of converting between a civil time and a GPS time is almost done
48  * by standard C library functions, though the C library functions use a
49  * different reference epoch, which is called the UNIX epoch.
50  *
51  * The tricky bits are:
52  * - What about time zones?
53  * - What about leap seconds?
54  *
55  * This code does not deal with time zones at all. All civil time structures
56  * are taken to be in Coordinated Universal Time or UTC. The user must convert
57  * local times into UTC before using these routines.
58  *
59  * Leap seconds are accounted for. But since the addition and subtraction of
60  * leap seconds is not deterministic, a leap second table needs to be
61  * maintained to account for the number of leap seconds in effect at a
62  * particular time.
63  *
64  * Leap seconds are defined as the difference between UTC and a (yet another)
65  * standard of time called the International Atomic Time or TAI. UTC is
66  * ultimately determined by the position of the stars, so it is occationally
67  * altered by a second to keep the location of fixed points on the celestial
68  * sphere correct to within plus or minus 0.9 seconds. TAI, like the GPS time
69  * used in LAL, is just the number of seconds since some epoch and is not
70  * affected by leap seconds. The difference between TAI and UTC, TAI-UTC,
71  * is the number of leap seconds.
72  *
73  * Note that UTC is fixed to atomic time, though there are an integer number
74  * of leap seconds. The real civil time, as measured in terms of points on
75  * the celestial sphere, is UT1. UTC and UT1 are kept within 0.9 seconds of
76  * each other by the introduction of leap seconds. But if what you want is
77  * UT1 note that UTC can be off by up to about a seconds. For this reason,
78  * we assume that accuracy at the second scale is sufficient, so the routines
79  * have only second precision. If you need more accuracy, you'll need to be
80  * monitoring UT1.
81  *
82  * Another way of representing the civil time in in terms of Julian days.
83  * There is a routine for converting a civil time into Julian days [in the same time system].
84  * The inverse conversion is not attempted.
85  *
86  */
87 /** @{ */
88 
89 /* change in TAI-UTC from previous second:
90  *
91  * return values:
92  * -1: TAI-UTC has decreased by one second (this hasn't happened yet).
93  * In this case, UTC will skip a second going from 23:59:58 at
94  * gpssec-1 to 00:00:00 (of the following day) at gpssec.
95  * 0: This is not a leap second: UTC has just advanced one second
96  * in going from gpssec-1 to gpssec.
97  * +1: TAI-UTC has increased by one second (this hasn't happened yet)
98  * In this case, UTC will add a second going from 23:59:59 at
99  * gpssec-1 to 23:59:60 (of the same day) at gpssec.
100  */
101 static int delta_tai_utc( INT4 gpssec )
102 {
103  int leap;
104 
105  /* assume calling function has already checked this */
106  /*
107  if ( gpssec <= leaps[0].gpssec )
108  {
109  fprintf( stderr, "error: don't know leap seconds before gps time %d\n",
110  leaps[0].gpssec );
111  abort();
112  }
113  */
114 
115  for ( leap = 1; leap < numleaps; ++leap )
116  if ( gpssec == leaps[leap].gpssec )
117  return leaps[leap].taiutc - leaps[leap-1].taiutc;
118 
119  return 0;
120 }
121 
122 /** Returns the leap seconds TAI-UTC at a given GPS second. */
123 int XLALLeapSeconds( INT4 gpssec /**< [In] Seconds relative to GPS epoch.*/ )
124 {
125  int leap;
126 
127  if ( gpssec < leaps[0].gpssec )
128  {
129  XLALPrintError( "XLAL Error - Don't know leap seconds before GPS time %d\n",
130  leaps[0].gpssec );
132  }
133 
134  /* scan leap second table and locate the appropriate interval */
135  for ( leap = 1; leap < numleaps; ++leap )
136  if ( gpssec < leaps[leap].gpssec )
137  break;
138 
139  return leaps[leap-1].taiutc;
140 }
141 
142 
143 /** Returns the leap seconds GPS-UTC at a given GPS second. */
144 int XLALGPSLeapSeconds( INT4 gpssec /**< [In] Seconds relative to GPS epoch.*/ )
145 {
146  int leapTAI;
147  int leapGPS;
148 
149  leapTAI = XLALLeapSeconds ( gpssec );
150  if ( leapTAI < 0 )
152 
153  leapGPS = leapTAI - 19; /* subtract 19 seconds to get leap-seconds wrt to GPS epoch */
154 
155  return leapGPS;
156 }
157 
158 
159 /** Returns the leap seconds TAI-UTC for a given UTC broken down time. */
160 int XLALLeapSecondsUTC( const struct tm *utc /**< [In] UTC as a broken down time.*/ )
161 {
162  REAL8 jd;
163  int leap;
164 
165  jd = XLALConvertCivilTimeToJD( utc );
166  if ( XLAL_IS_REAL8_FAIL_NAN( jd ) )
168 
169  if ( jd < leaps[0].jd )
170  {
171  XLALPrintError( "XLAL Error - Don't know leap seconds before Julian Day %9.1f\n", leaps[0].jd );
173  }
174 
175  /* scan leap second table and locate the appropriate interval */
176  for ( leap = 1; leap < numleaps; ++leap )
177  if ( jd < leaps[leap].jd )
178  break;
179 
180  return leaps[leap-1].taiutc;
181 }
182 
183 
184 /**
185  * Fill in derived fields, e.g. \c tm_wday and \c tm_yday, in a given UTC time structure
186  */
187 struct tm *XLALFillUTC( struct tm *utc /**< [In] UTC time in a broken down time structure. */ )
188 {
189  XLAL_CHECK_NULL(utc != NULL, XLAL_EINVAL);
190 
191  /* Most code in this function comes from time/mktime.c and time/strptime_l.c in
192  glibc v2.23 (commit de51ff8c0516e66554044b27656c6855a9c2ef25), licensed under
193  GNU Lesser General Public License, version 2.1 */
194 
195  /* Copy fields we need from 'utc', then erase it */
196  int sec = utc->tm_sec;
197  int min = utc->tm_min;
198  int hour = utc->tm_hour;
199  int mday = utc->tm_mday;
200  int mon = utc->tm_mon;
201  int year_requested = utc->tm_year;
202  XLAL_INIT_MEM(*utc);
203 
204  /* Ensure that month is in range, and set year accordingly. */
205  int mon_remainder = mon % 12;
206  int negative_mon_remainder = mon_remainder < 0;
207  int mon_years = mon / 12 - negative_mon_remainder;
208  long int lyear_requested = year_requested;
209  long int year = lyear_requested + mon_years;
210 
211  /* How many days come before each month (0-12) */
212  const unsigned short int __mon_yday[2][13] =
213  {
214  /* Normal years. */
215  { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 },
216  /* Leap years. */
217  { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 }
218  };
219 
220  /* Compute day of week */
221  /* We know that January 1st 1970 was a Thursday (= 4). Compute the
222  difference between this data in the one on TM and so determine
223  the weekday. */
224  int corr_year = 1900 + year - (mon < 2);
225  int wday = (-473
226  + (365 * (year - 70))
227  + (corr_year / 4)
228  - ((corr_year / 4) / 25) + ((corr_year / 4) % 25 < 0)
229  + (((corr_year / 4) / 25) / 4)
230  + __mon_yday[0][mon]
231  + mday - 1);
232  wday = ((wday % 7) + 7) % 7;
233 
234  /* Return 1 if 'year' + 1900 is a leap year. */
235  int leapyear =
236  /* Don't add 'year' to 1900, as that might overflow.
237  Also, work even if 'year' is negative. */
238  ((year & 3) == 0
239  && (year % 100 != 0
240  || ((year / 100) & 3) == (- (1900 / 100) & 3)));
241 
242  /* Calculate day of year from year, month, and day of month.
243  The result need not be in range. */
244  int mon_yday = ((__mon_yday[leapyear]
245  [mon_remainder + 12 * negative_mon_remainder])
246  - 1);
247  long int lmday = mday;
248  long int yday = mon_yday + lmday;
249 
250  /* Fill 'utc' */
251  utc->tm_sec = sec;
252  utc->tm_min = min;
253  utc->tm_hour = hour;
254  utc->tm_mday = mday;
255  utc->tm_mon = mon;
256  utc->tm_year = year;
257  utc->tm_wday = wday;
258  utc->tm_yday = yday;
259 
260  return utc;
261 
262 }
263 
264 
265 /**
266  * Returns the GPS seconds since the GPS epoch for a
267  * specified UTC time structure.
268  */
269 INT4 XLALUTCToGPS( const struct tm *utc /**< [In] UTC time in a broken down time structure. */ )
270 {
271  time_t unixsec;
272  INT4 gpssec;
273  int leapsec;
274 
275  /* make sure derived fields in 'utc' are filled correctly */
276  struct tm utc_filled = *utc;
277  utc = XLALFillUTC(&utc_filled);
278  if ( utc == NULL )
280 
281  /* compute leap seconds */
282  leapsec = XLALLeapSecondsUTC( utc );
283  if ( leapsec < 0 )
285  /* compute unix epoch time: seconds since 1970 JAN 1 0h UTC */
286  /* POSIX:2001 definition of seconds since the (UNIX) Epoch */
287  unixsec = utc->tm_sec + utc->tm_min*60 + utc->tm_hour*3600
288  + utc->tm_yday*86400 + (utc->tm_year-70)*31536000
289  + ((utc->tm_year-69)/4)*86400 - ((utc->tm_year-1)/100)*86400
290  + ((utc->tm_year+299)/400)*86400;
291  gpssec = unixsec;
292  gpssec -= XLAL_EPOCH_UNIX_GPS; /* change to gps epoch */
293  gpssec += leapsec - XLAL_EPOCH_GPS_TAI_UTC;
294  /* now check to see if this is an additional leap second */
295  if ( utc->tm_sec == 60 )
296  --gpssec;
297  return gpssec;
298 }
299 
300 
301 /**
302  * Returns a pointer to a tm structure representing the time
303  * specified in seconds since the GPS epoch.
304  */
305 struct tm * XLALGPSToUTC(
306  struct tm *utc, /**< [Out] Pointer to tm struct where result is stored. */
307  INT4 gpssec /**< [In] Seconds since the GPS epoch. */
308  )
309 {
310  time_t unixsec;
311  int leapsec;
312  int delta;
313  leapsec = XLALLeapSeconds( gpssec );
314  if ( leapsec < 0 )
316  unixsec = gpssec - leapsec + XLAL_EPOCH_GPS_TAI_UTC; /* get rid of leap seconds */
317  unixsec += XLAL_EPOCH_UNIX_GPS; /* change to unix epoch */
318  memset( utc, 0, sizeof( *utc ) ); /* blank out utc structure */
319  utc = gmtime_r( &unixsec, utc );
320  /* now check to see if we need to add a 60th second to UTC */
321  if ( ( delta = delta_tai_utc( gpssec ) ) > 0 )
322  utc->tm_sec += 1; /* delta only ever is one, right?? */
323  return utc;
324 }
325 
326 
327 /**
328  * Returns the Julian Day (JD) corresponding to the civil date and time given
329  * in a broken down time structure.
330  *
331  * The time system of the returned JD is the same as that of the input time.
332  *
333  * See \cite esaa1992 and \cite green1985 for details. First, some
334  * definitions:
335  *
336  * Mean Julian Year = 365.25 days
337  * Julian Epoch = 1 Jan 4713BCE, 12:00 GMT (4713 BC Jan 01d.5 GMT)
338  * Fundamental Epoch J2000.0 = 2001-01-01.5 TDB
339  *
340  * Julian Date is the amount of time elapsed since the Julian Epoch,
341  * measured in days and fractions of a day. There are a couple of
342  * complications arising from the length of a year: the Tropical Year is
343  * 365.2422 days. First, the Gregorian correction where 10 days
344  * (1582-10-05 through 1582-10-14) were eliminated. Second, leap years:
345  * years ending with two zeroes (e.g., 1700, 1800) are leap only if
346  * divisible by 400; so, 400 civil years contain 400 * 365.25 - 3 = 146097
347  * days. So, the Julian Date of J2000.0 is JD 2451545.0, and thus the
348  * Julian Epoch = J2000.0 + (JD - 2451545) / 365.25, i.e., number of years
349  * elapsed since J2000.0.
350  *
351  * One algorithm for computing the Julian Day is from \cite vfp1979 based
352  * on a formula in \cite esaa1992 where the algorithm is due to
353  * \cite fvf1968 and ``compactified'' by P. M. Muller and R. N. Wimberly.
354  * The formula is
355  *
356  * \f[
357  * jd = 367 \times y - 7 \times (y + (m + 9)/12)/4 - 3 \times ((y + (m -
358  * 9)/7)/100 + 1)/4 + 275 \times m/9 + d + 1721029
359  * \f]
360  *
361  * where jd is the Julian day number, y is the year, m is the month (1-12),
362  * and d is the day (1-31). This formula is valid only for JD > 0, i.e.,
363  * after -4713 Nov 23 = 4712 BCE Nov 23.
364  *
365  * A shorter formula from the same reference, but which only works for
366  * dates since 1900 March is:
367  *
368  * \f[
369  * jd = 367 \times y - 7 \times (y + (m + 9)/12)/4 + 275 \times m/9 + d +
370  * 1721014
371  * \f]
372  *
373  * We will use this shorter formula since there is unlikely to be any
374  * analyzable data from before 1900 March.
375  *
376  */
377 REAL8 XLALConvertCivilTimeToJD( const struct tm *civil /**< [In] civil time in a broken down time structure. */ )
378 {
379  const int sec_per_day = 60 * 60 * 24; /* seconds in a day */
380  int year, month, day, sec;
381  REAL8 jd;
382 
383  /* this routine only works for dates after 1900 */
384  if ( civil->tm_year <= 0 )
385  {
386  XLALPrintError( "XLAL Error - Year must be after 1900\n" );
388  }
389 
390  year = civil->tm_year + 1900;
391  month = civil->tm_mon + 1; /* month is in range 1-12 */
392  day = civil->tm_mday; /* day is in range 1-31 */
393  sec = civil->tm_sec + 60*(civil->tm_min + 60*(civil->tm_hour)); /* seconds since midnight */
394 
395  jd = 367*year - 7*(year + (month + 9)/12)/4 + 275*month/9 + day + 1721014;
396  /* note: Julian days start at noon: subtract half a day */
397  jd += (REAL8)sec/(REAL8)sec_per_day - 0.5;
398 
399  return jd;
400 } // XLALConvertCivilTimeToJD()
401 
402 /**
403  * Returns the Modified Julian Day MJD corresponding to the civil date and time given
404  * in a broken down time structure (using the same time system as the input).
405  *
406  * Note:
407  * - MJD numbers starts at midnight rather than noon.
408  *
409  */
410 REAL8 XLALConvertCivilTimeToMJD( const struct tm *civil /**< [In] civil time in a broken down time structure. */ )
411 {
412  return XLAL_JD_TO_MJD ( XLALConvertCivilTimeToJD ( civil ) );
413 }
414 
415 /** \deprecated Use XLALConvertCivilTimeToJD() instead, this is only provided for pylal backwards compatibility.
416  * (see #1856)
417  */
418 REAL8
419 XLALJulianDay ( const struct tm *civil )
420 {
421  XLAL_PRINT_DEPRECATION_WARNING ( "XLALConvertCivilTimeToJD" );
422  return XLALConvertCivilTimeToJD ( civil );
423 } // XLALJulianDay()
424 
425 /** \deprecated Use XLALConvertCivilTimeToJD() instead, this is only provided for pylal backwards compatibility.
426  * (see #1856)
427  */
428 INT4
429 XLALModifiedJulianDay ( const struct tm *civil )
430 {
431  XLAL_PRINT_DEPRECATION_WARNING ( "XLALConvertCivilTimeToJD" );
432  return XLALConvertCivilTimeToJD ( civil );
433 } // XLALModifiedJulianDay()
434 
435 
436 /** @} */
struct tm * XLALFillUTC(struct tm *utc)
Fill in derived fields, e.g.
#define XLAL_EPOCH_UNIX_GPS
The UNIX time of the GPS origin epoch.
Definition: Date.h:82
#define XLAL_ERROR(...)
Macro to invoke a failure from a XLAL routine returning an integer.
Definition: XLALError.h:701
#define XLAL_CHECK_NULL(assertion,...)
Macro to test an assertion and invoke a failure if it is not true in a function that returns a pointe...
Definition: XLALError.h:826
#define XLAL_JD_TO_MJD(jd)
Modified Julian Day for specified civil time structure.
Definition: Date.h:90
int XLALLeapSecondsUTC(const struct tm *utc)
Returns the leap seconds TAI-UTC for a given UTC broken down time.
REAL8 XLALJulianDay(const struct tm *civil)
REAL8 XLALConvertCivilTimeToJD(const struct tm *civil)
Returns the Julian Day (JD) corresponding to the civil date and time given in a broken down time stru...
static const struct leaps_table leaps[]
int XLALLeapSeconds(INT4 gpssec)
Returns the leap seconds TAI-UTC at a given GPS second.
#define XLAL_EPOCH_GPS_TAI_UTC
Leap seconds (TAI-UTC) on the GPS epoch (1980 JAN 6 0h UTC)
Definition: Date.h:88
Input domain error.
Definition: XLALError.h:411
#define XLAL_ERROR_NULL(...)
Macro to invoke a failure from a XLAL routine returning a pointer.
Definition: XLALError.h:714
double REAL8
Double precision real floating-point number (8 bytes).
static const int numleaps
#define XLAL_IS_REAL8_FAIL_NAN(val)
Tests if val is a XLAL REAL8 failure NaN.
Definition: XLALError.h:394
static int delta_tai_utc(INT4 gpssec)
REAL8 XLALConvertCivilTimeToMJD(const struct tm *civil)
Returns the Modified Julian Day MJD corresponding to the civil date and time given in a broken down t...
int XLALPrintError(const char *fmt,...)
Definition: XLALError.c:68
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
Prints a deprecation warning at the "warning" verbosity level.
Definition: XLALError.h:229
#define gmtime_r(timep, result)
Definition: XLALCivilTime.c:30
INT4 XLALModifiedJulianDay(const struct tm *civil)
INT4 XLALUTCToGPS(const struct tm *utc)
Returns the GPS seconds since the GPS epoch for a specified UTC time structure.
#define XLAL_ERROR_REAL8(...)
Macro to invoke a failure from a XLAL routine returning a REAL8.
Definition: XLALError.h:753
#define XLAL_INIT_MEM(x)
MACRO to initialize arbitrary variable &#39;x&#39; to zero.
int32_t INT4
Four-byte signed integer.
Internal function call failed bit: "or" this with existing error number.
Definition: XLALError.h:463
Invalid argument.
Definition: XLALError.h:410
struct tm * XLALGPSToUTC(struct tm *utc, INT4 gpssec)
Returns a pointer to a tm structure representing the time specified in seconds since the GPS epoch...
int XLALGPSLeapSeconds(INT4 gpssec)
Returns the leap seconds GPS-UTC at a given GPS second.