RALICE/icepack/iceconvert/rdmc_time.c

   1 /*
   2  * Functions related to the different time conversions in rdmc
   3  */
   4
   5 #include <math.h>
   6 #include <time.h>
   7 #include "rdmc.h"
   8
   9 #define JD2000   (2451545.0) /* Julian day for Jan 1st 2000 12:00 noon */
  10                 /* = standard epoch to be used as reference */
  11 #define JCENTURY  (36525.0) /* number od days per Julian century */
  12
  13 /************************************************************************/
  14 /* rdmc_o_dateconv() converts a YYMMDD-date into the unix time          */
  15 /*              (sec since 1.1.70)                                      */
  16 /************************************************************************/
  17
  18 time_t rdmc_o_dateconv(int date)
  19 {
  20   struct tm d;
  21
  22   if (date <= 0) return RDMC_NA;
  23
  24   d.tm_sec = 0;               /* second (0-61, allows for leap seconds) */
  25   d.tm_min = 0;                                        /* minute (0-59) */
  26   d.tm_hour = 12;                                        /* hour (0-23) */
  27   d.tm_mday = date % 100;                    /* day of the month (1-31) */
  28   d.tm_mon = (date/100) % 100 - 1;                      /* month (0-11) */
  29   d.tm_year = date/10000;                           /* years since 1900 */
  30   d.tm_isdst = 0;        /* non-0 if daylight savings time is in effect */
  31
  32   return mktime(&d);                        /* time in sec since 1.1.70 */
  33
  34 } /* function o_dateconv() */
  35
  36 /************************************************************************/
  37 /* rdmc_o_rdateconv() converts the unix time into the YYMMDD format     */
  38 /************************************************************************/
  39 int rdmc_o_rdateconv(time_t time)
  40 {
  41   struct tm *date;
  42
  43   if (time <= 0) return RDMC_NA;
  44
  45   date = gmtime(&time);
  46   return (date->tm_year * 10000 + (date->tm_mon + 1) * 100 + date->tm_mday);
  47 } /* function o_rdateconv() */
  48
  49 /************************************************************************/
  50 /* rdmc_gps_to_mjd() converts the time in GPS/UT year and day into      */
  51 /*     Modified Julian Day (50084 < mjd < 51545 for 1996-99)            */
  52 /* NOTES: jdyr = (Julian date of 1200 UT 1/1/year) = GPS/UT day one     */
  53 /************************************************************************/
  54 int rdmc_gps_to_mjd(int gpsyear, int gpsday)
  55
  56 /* reference: http://www.starlink.rl.ac.uk/ function cldj.f
  57  * or Zombeck p107 equation for julian date */
  58
  59 {
  60   int mjdyr, mjd;
  61
  62   if ((gpsyear <= 0) || (gpsday <=0 ) )
  63     return RDMC_NA;
  64
  65   mjdyr = (1461*(gpsyear+4799))/4 - (3*((gpsyear+4899)/100))/4 - 2431738.5;
  66   mjd = (mjdyr + gpsday) - 1;
  67   return mjd;
  68
  69 } /* gps_to_mjd() */
  70
  71 /************************************************************************/
  72 /* rdmc_mjd_to_gps() converts the time in Modified Julian Days to       */
  73 /*              GPS/UT year and day within year                         */
  74 /* This uses an approximation which may err high in late December.      */
  75 /* It then checks using the exact calculation of mjdyr and gets gpsday. */
  76 /* NOTES: jdyr = (Julian date of 1200 UT 1/1/year) = GPS/UT day one     */
  77 /************************************************************************/
  78 void rdmc_mjd_to_gps(int mjd, int *gpsyear, int *gpsday)
  79 {
  80   int year,mjdyr;
  81
  82   if ( mjd <=0  ){
  83     *gpsyear=RDMC_NA;
  84     *gpsday=RDMC_NA;
  85     return;
  86   }
  87
  88   year = (4*(mjd + 678960)) / 1461; /* Approx, erring high */
  89   mjdyr = (1461*(year+4799))/4 - (3*((year+4899)/100))/4 - 2431738.5;
  90   if (mjdyr > mjd) {
  91     year -= 1;
  92     mjdyr = (1461*(year+4799))/4 - (3*((year+4899)/100))/4 - 2431738.5;
  93   }
  94   *gpsyear = year;
  95   *gpsday = mjd - mjdyr + 1;
  96   return;
  97 } /* mjd_to_gps() */
  98
  99 /************************************************************************/
 100 /* rdmc_mjd_to_tjd() converts the time in Modified Julian Days to       */
 101 /*              Truncated Julian Days.                                  */
 102 /*              this is a trivial conversion but useful for GRB         */
 103 /*              analysis with the batse catalog                         */
 104 /************************************************************************/
 105 void rdmc_mjd_to_tjd(int mjd, int secs, int ns, int *tjd, double *sec)
 106 {
 107   *tjd=mjd-40000.0;
 108   *sec=(double)secs+(double)ns/1e9;
 109 }
 110
 111 /************************************************************************/
 112 /* rdmc_tjd_to_mjd() converts the time in Truncated Julian Days to      */
 113 /*              Modified Julian Days.                                   */
 114 /*              this is a trivial conversion but useful for GRB         */
 115 /*              analysis with the batse catalog                         */
 116 /************************************************************************/
 117 void rdmc_tjd_to_mjd(int tjd, double sec, int *mjd, int *secs, int *ns)
 118 {
 119   *mjd = tjd+40000.0;
 120   *secs = sec;
 121   *ns = (sec - *secs) * 1e9;
 122 }
 123
 124 /************************************************************************/
 125 /* rdmc_mjd_to_jd() converts the time in Modified Julian Days to       */
 126 /*              Julian Days..                                           */
 127 /*              this is a trivial conversion but stands here for        */
 128 /*              completeness                                            */
 129 /************************************************************************/
 130 void rdmc_mjd_to_jd(int mjd, int secs, int ns, double *jd)
 131 {
 132   *jd=2400000.5+mjd+((secs+ns/1e9)/86400);
 133 }
 134
 135 /************************************************************************/
 136 /* rdmc_jd_to_mjd() converts the time in Julian Days to                 */
 137 /*              Modified Julian Days.                                   */
 138 /*              this is a trivial conversion but stands here for        */
 139 /*              completeness                                            */
 140 /************************************************************************/
 141 void rdmc_jd_to_mjd(double jd, int *mjd, int *secs, int *ns)
 142 {
 143   *mjd=jd-2400000.5;
 144   *secs=(jd-(double)*mjd)*86400;
 145   *ns=((jd-(double)*mjd)*86400-(double)*secs)*1e9;
 146 }
 147 /************************************************************************/
 148 /* rdmc_mjd_to_unixtime() converts the time (unix in sec since 1.1.70)  */
 149 /*                   from mjd and secs                                  */
 150 /************************************************************************/
 151 time_t rdmc_mjd_to_unixtime(int mjd, int secs)
 152 {
 153   time_t unix_time;
 154
 155   if ((mjd <= 0) || (secs <= 0)) return RDMC_NA;  /* if no mjd->no unix time*/
 156   unix_time = (mjd - 40587)* 86400 + secs;     /* 40587 is the 1.1.1970 */
 157   return unix_time;
 158 } /* mjd_to_unixtime() */
 159
 160 /************************************************************************/
 161 /* rdmc_unixtime_to_mjd() converts UNIX time(sec since 1.1.70) into mjd */
 162 /************************************************************************/
 163 int rdmc_unixtime_to_mjd(time_t unix_time)
 164 {
 165   if (unix_time <= 0) return RDMC_NA;
 166   return unix_time/86400 + 40587;
 167 } /* unixtime_to_mjd() */
 168
 169
 170 /************************************************************************/
 171 /* rdmc_unixtime_to_mjdsecs() converts the time (unic in sec since      */
 172 /*                       1.1.70) into secs after mjd                    */
 173 /************************************************************************/
 174 int rdmc_unixtime_to_mjdsecs(time_t unix_time)
 175 {
 176   if (unix_time <= 0) return RDMC_NA;
 177   return unix_time % 86400;
 178 } /* unixtime_to_mjdsecs() */
 179
 180
 181 /************************************************************************/
 182 /* rdmc_jd_to_gmst() converts mjd date and time to Greenwich mean       */
 183 /*                       sidereal time                                  */
 184 /************************************************************************/
 185 void rdmc_jd_to_gmst(double jd, double *gmst)
 186 {
 187 /*
 188  *  GMST in hours
 189  *  GMST from derived NOVAS package sidereal_time() function.
 190  *  NOVAS is online at http://aa.usno.navy.mil/aa/
 191  */
 192
 193   double tdt, gst;
 194   int igst;
 195
 196   tdt = (jd -  JD2000) / JCENTURY;
 197   gst = ((((-6.2e-6 * tdt + 0.093104) * tdt + 184.812866) * tdt
 198            + 67310.54841) + 3.1644e9 * tdt);
 199   igst= (int)(gst/86400.0);
 200   *gmst = gst/3600.0 - (double)(igst*24.0);   /* gmst is smaller than 24 */
 201   if (*gmst < 0.0) *gmst += 24.0;     /* gmst is larger than 0 */
 202
 203 } /*  rdmc_jd_to_gmst() */
 204
 205
 206
 207 /************************************************************************/
 208 /* rdmc_gmst_to_jd() converts Greenwich mean sidereal time  to mjd      */
 209 /*                       sidereal time                                  */
 210 /************************************************************************/
 211 void rdmc_gmst_to_jd(double gmst, int intjd, double *jd)
 212 {
 213 /*
 214  *  GMST in hours, intjd is integer part of jd
 215  *  GMST from derived NOVAS package sidereal_time() function.
 216  *  NOVAS is online at http://aa.usno.navy.mil/aa/
 217  *  Inverted direction according to functions GMST->UT and UT->GMST
 218  *  From "practical Astronomy with your Calculator", Cambridge UP 1988
 219  */
 220
 221
 222   /* Caution: Function untested in any respect */
 223   double tdt, gst, fracjd;
 224
 225   tdt = ((double)intjd -  JD2000) / JCENTURY;
 226   gst = ((((-6.2e-6 * tdt + 0.093104) * tdt + 184.812866) * tdt
 227            + 67310.54841) + 3.1644e9 * tdt);
 228   fracjd = gmst - gst;
 229   fracjd = fmod ((fracjd / 3600.0), 24.0);
 230   if (fracjd < 0.0) fracjd += 24.0;
 231   fracjd=fracjd/24.0*0.9972695663;
 232   *jd=fracjd+(double)intjd;
 233
 234 } /*  rdmc_jd_to_gmst() */
 235