sgp_time.c

00001 /*
00002  * Unit SGP_Time
00003  *       Author:  Dr TS Kelso
00004  * Original Version:  1992 Jun 02
00005  * Current Revision:  2000 Jan 22
00006  * Modified for Y2K:  1999 Mar 07
00007  *          Version:  2.05
00008  *        Copyright:  1992-1999, All Rights Reserved
00009  * Version 1.50 added Y2K support. Due to limitations in the current
00010  * format of the NORAD two-line element sets, however, only dates
00011  * through 2056 December 31/2359 UTC are valid.
00012  * Version 1.60 modifies Calendar_Date to ensure date matches time
00013  * resolution and modifies Time_of_Day to make it more robust.
00014  * Version 2.00 adds Julian_Date, Date_Time, and Check_Date to support
00015  * checking for valid date/times, permitting the use of Time_to_UTC and
00016  * Time_from_UTC for UTC/local time conversions.
00017  * Version 2.05 modifies UTC_offset to allow non-integer offsets.
00018  *
00019  *   Ported to C by: Neoklis Kyriazis  April 9  2001
00020  */
00021 
00022 #include "sgp4sdp4.h"
00023 
00024 /* The function Julian_Date_of_Epoch returns the Julian Date of     */
00025 /* an epoch specified in the format used in the NORAD two-line      */
00026 /* element sets. It has been modified to support dates beyond       */
00027 /* the year 1999 assuming that two-digit years in the range 00-56   */
00028 /* correspond to 2000-2056. Until the two-line element set format   */
00029 /* is changed, it is only valid for dates through 2056 December 31. */
00030 
00031 double
00032 Julian_Date_of_Epoch(double epoch)
00033 { 
00034   double year,day;
00035 
00036   /* Modification to support Y2K */
00037   /* Valid 1957 through 2056     */
00038   day = modf(epoch*1E-3, &year)*1E3;
00039   if( year < 57 )
00040     year = year + 2000;
00041   else
00042     year = year + 1900;
00043   /* End modification */
00044 
00045   return( Julian_Date_of_Year(year) + day );
00046 } /*Function Julian_Date_of_Epoch*/
00047 
00048 /*------------------------------------------------------------------*/
00049 
00050 /* Converts a Julian epoch to NORAD TLE epoch format */
00051 double
00052 Epoch_Time(double jd)
00053 {  
00054   double yr,_time,epoch_time;
00055   struct tm edate;
00056 
00057   Calendar_Date(jd, &edate);
00058   yr = edate.tm_year - 100*(edate.tm_year/100) ;
00059   _time = Frac(jd + 0.5);
00060   epoch_time =  yr*1000
00061                 + DOY(edate.tm_year, edate.tm_mon, edate.tm_mday)
00062                 + _time;
00063 
00064   return( epoch_time );
00065 } /*Function Epoch_Time*/
00066 
00067 /*------------------------------------------------------------------*/
00068 
00069 /* The function DOY calculates the day of the year for the specified */
00070 /* date. The calculation uses the rules for the Gregorian calendar   */
00071 /* and is valid from the inception of that calendar system.          */
00072 int
00073 DOY(int yr, int mo, int dy)
00074 {
00075   const int days[] = {31,28,31,30,31,30,31,31,30,31,30,31};
00076   int i,day;
00077 
00078   day = 0;
00079   for( i = 0; i < mo-1; i++ )
00080     day += days[i];
00081   day = day + dy;
00082 
00083   /* Leap year correction */
00084   if( 
00085      (yr%4 == 0) && ((yr%100 != 0) || (yr%400 == 0)) && (mo>2)
00086      )
00087     day++;
00088 
00089   return( day );
00090 } /*Function DOY*/
00091 
00092 /*------------------------------------------------------------------*/
00093 
00094 /* Fraction_of_Day calculates the fraction of */
00095 /* a day passed at the specified input time.  */
00096 double
00097 Fraction_of_Day(int hr,int mi,int se)
00098 {
00099   return( (hr + (mi + se/60.0)/60.0)/24.0 );
00100 } /*Function Fraction_of_Day*/
00101 
00102 /*------------------------------------------------------------------*/
00103 
00104 /* The function Calendar_Date converts a Julian Date to a struct tm.   */
00105 /* Only the members tm_year, tm_mon and tm_mday are calculated and set */
00106 void
00107 Calendar_Date(double jd, struct tm *cdate)
00108 {
00109   /* Astronomical Formulae for Calculators, Jean Meeus, pages 26-27 */
00110   int Z,month;
00111   double A,B,C,D,E,F,alpha,day,year,factor;
00112 
00113   factor = 0.5/secday/1000;
00114   F = Frac(jd + 0.5);
00115   if (F + factor >= 1.0)
00116     {
00117       jd = jd + factor;
00118       F  = 0.0;
00119     } /*if*/
00120   Z = Round(jd);
00121   if( Z < 2299161 )
00122     A = Z;
00123   else
00124     {
00125       alpha = Int((Z - 1867216.25)/36524.25);
00126       A = Z + 1 + alpha - Int(alpha/4);
00127     } /*else*/
00128   B = A + 1524;
00129   C = Int((B - 122.1)/365.25);
00130   D = Int(365.25 * C);
00131   E = Int((B - D)/30.6001);
00132   day = B - D - Int(30.6001 * E) + F;
00133 
00134   if( E < 13.5 )
00135     month = Round(E - 1);
00136   else
00137     month = Round(E - 13);
00138   if( month > 2.5 )
00139     year = C - 4716;
00140   else
00141     year = C - 4715;
00142 
00143   cdate->tm_year = (int) year;
00144   cdate->tm_mon = month;
00145   cdate->tm_mday = (int) floor(day);
00146 
00147 } /*Function Calendar_Date*/
00148 
00149 /*------------------------------------------------------------------*/
00150 
00151 /* Time_of_Day takes a Julian Date and calculates the clock time */
00152 /* portion of that date. Only tm_hour, tm_min and tm_sec are set */
00153 void
00154 Time_of_Day(double jd, struct tm *cdate)
00155 {
00156   int hr,mn,sc;
00157   double _time;
00158 
00159   _time = Frac(jd - 0.5)*secday;
00160   _time = Round(_time);
00161   hr = floor(_time/3600.0);
00162   _time = _time - 3600.0*hr;
00163   if( hr == 24 ) hr = 0;
00164   mn = floor(_time/60.0);
00165   sc = _time - 60.0*mn;
00166   cdate->tm_hour = hr;
00167   cdate->tm_min = mn;
00168   cdate->tm_sec = sc;
00169 
00170 } /*Function Time_of_Day*/
00171 
00172 /*------------------------------------------------------------------*/
00173 
00174 /* The function Julian_Date converts a standard calendar   */
00175 /* date and time to a Julian Date. The procedure Date_Time */
00176 /* performs the inverse of this function. */
00177 double
00178 Julian_Date(struct tm *cdate)
00179 {
00180   double julian_date;
00181 
00182   julian_date = Julian_Date_of_Year(cdate->tm_year) + 
00183                 DOY(cdate->tm_year,cdate->tm_mon,cdate->tm_mday) +
00184                 Fraction_of_Day(cdate->tm_hour,cdate->tm_min,cdate->tm_sec);
00185 
00186   return( julian_date );
00187 } /*Function Julian_Date */
00188 
00189 /*------------------------------------------------------------------*/
00190 
00191 
00192 /*  Date_Time()
00193  *
00194  *  The function Date_Time() converts a Julian Date to
00195  *  standard calendar date and time. The function
00196  *  Julian_Date() performs the inverse of this function.
00197  */
00198 
00199 void
00200 Date_Time(double julian_date, struct tm *cdate)
00201 {
00202   time_t jtime;
00203 
00204   jtime = (julian_date - 2440587.5)*86400.;
00205   *cdate = *gmtime( &jtime );
00206 
00207 } /* End of Date_Time() */
00208 
00209 
00210 /*------------------------------------------------------------------*/
00211 
00212 /* The procedure Check_Date can be used as a check to see if a calendar    */
00213 /* date and time are valid. It works by first converting the calendar      */
00214 /* date and time to a Julian Date (which allows for irregularities, such   */
00215 /* as a time greater than 24 hours) and then converting back and comparing.*/
00216 int
00217 Check_Date(struct tm *cdate)
00218 {
00219   double jt;
00220   struct tm chkdate;
00221 
00222   jt = Julian_Date(cdate);
00223   Date_Time(jt, &chkdate);
00224 
00225   if( (cdate->tm_year == chkdate.tm_year) &&
00226       (cdate->tm_mon  == chkdate.tm_mon ) &&
00227       (cdate->tm_mday == chkdate.tm_mday) &&
00228       (cdate->tm_hour == chkdate.tm_hour) &&
00229       (cdate->tm_min  == chkdate.tm_min ) &&
00230       (cdate->tm_sec  == chkdate.tm_sec ) )
00231     return ( 1 );
00232   else
00233     return( 0 );
00234 
00235 } /*Procedure Check_Date*/
00236 
00237 /*------------------------------------------------------------------*/
00238 
00239 /* Procedures Time_to_UTC and Time_from_UTC are used to  */
00240 /* convert 'struct tm' dates between UTC and local time. */
00241 /* The procedures JD_to_UTC and JD_from_UTC are used to  */
00242 /* do the same thing working directly with Julian dates. */
00243 
00244 struct tm
00245 Time_to_UTC(struct tm *cdate)
00246 {
00247   time_t tdate;
00248 
00249   tdate = mktime(cdate);
00250   return( *gmtime(&tdate) );
00251 } /*Procedure Time_to_UTC*/
00252 
00253 /*------------------------------------------------------------------*/
00254 
00255 struct tm 
00256 Time_from_UTC(struct tm *cdate)
00257 {
00258   time_t tdate;
00259 
00260   tdate = mktime(cdate);
00261   return( *localtime(&tdate) );
00262 } /*Procedure Time_from_UTC*/
00263 
00264 /*------------------------------------------------------------------*/
00265 
00266 /*
00267  BSD systems don't define the timezone variable, so the following two
00268  routines won't work.  They're not used anyway in the example main(),
00269  so we might as well comment them out.  
00270 */
00271 
00272 #if 0
00273 
00274 double
00275 JD_to_UTC(double jt)
00276 {
00277   extern long timezone;
00278   struct tm cdate;
00279 
00280   time_t t = 0;
00281 
00282   cdate = *localtime( &t );
00283   jt = jt - timezone/secday;
00284   if( cdate.tm_isdst )
00285     jt= jt - 1.0/24.0;
00286 
00287   return( jt );
00288 } /*Procedure JD_to_UTC*/
00289 
00290 /*------------------------------------------------------------------*/
00291 
00292 double
00293 JD_from_UTC(double jt)
00294 {
00295   extern long timezone;
00296   struct tm cdate;
00297   time_t t = 0;
00298 
00299   cdate = *localtime( &t );
00300   jt = jt + timezone/secday;
00301   if( cdate.tm_isdst )
00302     jt= jt + 1.0/24.0;
00303 
00304   return( jt );
00305 } /*Procedure JD_from_UTC*/
00306 
00307 #endif
00308 
00309 /*------------------------------------------------------------------*/
00310 
00311 /* The function Delta_ET has been added to allow calculations on   */
00312 /* the position of the sun.  It provides the difference between UT */
00313 /* (approximately the same as UTC) and ET (now referred to as TDT).*/
00314 /* This function is based on a least squares fit of data from 1950 */
00315 /* to 1991 and will need to be updated periodically. */
00316 
00317 double
00318 Delta_ET(double year)
00319 {
00320   /* Values determined using data from 1950-1991 in the 1990 
00321      Astronomical Almanac.  See DELTA_ET.WQ1 for details. */
00322 
00323   double delta_et;
00324 
00325   delta_et = 26.465 + 0.747622*(year - 1950) +
00326              1.886913*sin(twopi*(year - 1975)/33);
00327 
00328   return( delta_et );
00329 } /*Function Delta_ET*/
00330 
00331 /*------------------------------------------------------------------*/
00332 
00333 /* The function Julian_Date_of_Year calculates the Julian Date  */
00334 /* of Day 0.0 of {year}. This function is used to calculate the */
00335 /* Julian Date of any date by using Julian_Date_of_Year, DOY,   */
00336 /* and Fraction_of_Day. */
00337 
00338 double
00339 Julian_Date_of_Year(double year)
00340 {
00341   /* Astronomical Formulae for Calculators, Jean Meeus, */
00342   /* pages 23-25. Calculate Julian Date of 0.0 Jan year */
00343 
00344   long A,B,i;
00345   double jdoy;
00346 
00347   year = year-1;
00348   i = year/100;
00349   A = i;
00350   i = A/4;
00351   B = 2-A+i;
00352   i = 365.25*year;
00353   i += 30.6001*14;
00354   jdoy = i+1720994.5+B;
00355 
00356   return (jdoy);
00357 }  /*Function Julian_Date_of_Year*/
00358 
00359 /*------------------------------------------------------------------*/
00360 
00361 /* The function ThetaG calculates the Greenwich Mean Sidereal Time */
00362 /* for an epoch specified in the format used in the NORAD two-line */
00363 /* element sets. It has now been adapted for dates beyond the year */
00364 /* 1999, as described above. The function ThetaG_JD provides the   */
00365 /* same calculation except that it is based on an input in the     */
00366 /* form of a Julian Date. */
00367 
00368 double
00369 ThetaG(double epoch, deep_arg_t *deep_arg)
00370 {
00371 /* Reference:  The 1992 Astronomical Almanac, page B6. */
00372 
00373   double year,day,UT,jd,TU,GMST,_ThetaG;
00374 
00375 /* Modification to support Y2K */
00376 /* Valid 1957 through 2056     */
00377   day = modf(epoch*1E-3,&year)*1E3;
00378   if(year < 57)
00379     year += 2000;
00380   else
00381     year += 1900;
00382   /* End modification */
00383 
00384   UT   = modf(day,&day);
00385   jd   = Julian_Date_of_Year(year)+day;
00386   TU   = (jd-2451545.0)/36525;
00387   GMST = 24110.54841+TU*(8640184.812866+TU*(0.093104-TU* 6.2E-6));
00388   GMST = Modulus(GMST+secday*omega_E*UT,secday);
00389   _ThetaG = twopi*GMST/secday;
00390   deep_arg->ds50 = jd-2433281.5+UT;
00391   _ThetaG = FMod2p(6.3003880987*deep_arg->ds50+1.72944494);
00392 
00393   return (_ThetaG);
00394 } /* Function ThetaG */
00395 
00396 /*------------------------------------------------------------------*/
00397 
00398 double
00399 ThetaG_JD(double jd)
00400 {
00401 /* Reference:  The 1992 Astronomical Almanac, page B6. */
00402 
00403   double UT,TU,GMST;
00404 
00405   UT   = Frac(jd + 0.5);
00406   jd   = jd - UT;
00407   TU   = (jd - 2451545.0)/36525;
00408   GMST = 24110.54841 + TU * (8640184.812866 + TU * (0.093104 - TU * 6.2E-6));
00409   GMST = Modulus(GMST + secday*omega_E*UT,secday);
00410 
00411   return( twopi * GMST/secday );
00412 } /*Function ThetaG_JD*/
00413 
00414 /*------------------------------------------------------------------*/
00415 
00416 /* Gets calendar time from time() and produces a UTC calendar date */
00417 void
00418 UTC_Calendar_Now( struct tm *cdate )
00419 {
00420   time_t t;
00421 
00422   t = time(0);
00423   *cdate = *gmtime(&t);
00424   cdate->tm_year += 1900;
00425   cdate->tm_mon += 1;
00426 
00427 } /* End UTC_Calendar_Now */
00428 /*------------------------------------------------------------------*/