Diff: sgp4sdp4/sgp_time.c

Revision:

0:0a841b89d614

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sgp4sdp4/sgp_time.c	Mon Oct 11 10:34:55 2010 +0000
@@ -0,0 +1,428 @@
+/*
+ * Unit SGP_Time
+ *       Author:  Dr TS Kelso
+ * Original Version:  1992 Jun 02
+ * Current Revision:  2000 Jan 22
+ * Modified for Y2K:  1999 Mar 07
+ *          Version:  2.05
+ *        Copyright:  1992-1999, All Rights Reserved
+ * Version 1.50 added Y2K support. Due to limitations in the current
+ * format of the NORAD two-line element sets, however, only dates
+ * through 2056 December 31/2359 UTC are valid.
+ * Version 1.60 modifies Calendar_Date to ensure date matches time
+ * resolution and modifies Time_of_Day to make it more robust.
+ * Version 2.00 adds Julian_Date, Date_Time, and Check_Date to support
+ * checking for valid date/times, permitting the use of Time_to_UTC and
+ * Time_from_UTC for UTC/local time conversions.
+ * Version 2.05 modifies UTC_offset to allow non-integer offsets.
+ *
+ *   Ported to C by: Neoklis Kyriazis  April 9  2001
+ */
+
+#include "sgp4sdp4.h"
+
+/* The function Julian_Date_of_Epoch returns the Julian Date of     */
+/* an epoch specified in the format used in the NORAD two-line      */
+/* element sets. It has been modified to support dates beyond       */
+/* the year 1999 assuming that two-digit years in the range 00-56   */
+/* correspond to 2000-2056. Until the two-line element set format   */
+/* is changed, it is only valid for dates through 2056 December 31. */
+
+double
+Julian_Date_of_Epoch(double epoch)
+{ 
+  double year,day;
+
+  /* Modification to support Y2K */
+  /* Valid 1957 through 2056     */
+  day = modf(epoch*1E-3, &year)*1E3;
+  if( year < 57 )
+    year = year + 2000;
+  else
+    year = year + 1900;
+  /* End modification */
+
+  return( Julian_Date_of_Year(year) + day );
+} /*Function Julian_Date_of_Epoch*/
+
+/*------------------------------------------------------------------*/
+
+/* Converts a Julian epoch to NORAD TLE epoch format */
+double
+Epoch_Time(double jd)
+{  
+  double yr,_time,epoch_time;
+  struct tm edate;
+
+  Calendar_Date(jd, &edate);
+  yr = edate.tm_year - 100*(edate.tm_year/100) ;
+  _time = Frac(jd + 0.5);
+  epoch_time =  yr*1000
+                + DOY(edate.tm_year, edate.tm_mon, edate.tm_mday)
+                + _time;
+
+  return( epoch_time );
+} /*Function Epoch_Time*/
+
+/*------------------------------------------------------------------*/
+
+/* The function DOY calculates the day of the year for the specified */
+/* date. The calculation uses the rules for the Gregorian calendar   */
+/* and is valid from the inception of that calendar system.          */
+int
+DOY(int yr, int mo, int dy)
+{
+  const int days[] = {31,28,31,30,31,30,31,31,30,31,30,31};
+  int i,day;
+
+  day = 0;
+  for( i = 0; i < mo-1; i++ )
+    day += days[i];
+  day = day + dy;
+
+  /* Leap year correction */
+  if( 
+     (yr%4 == 0) && ((yr%100 != 0) || (yr%400 == 0)) && (mo>2)
+     )
+    day++;
+
+  return( day );
+} /*Function DOY*/
+
+/*------------------------------------------------------------------*/
+
+/* Fraction_of_Day calculates the fraction of */
+/* a day passed at the specified input time.  */
+double
+Fraction_of_Day(int hr,int mi,int se)
+{
+  return( (hr + (mi + se/60.0)/60.0)/24.0 );
+} /*Function Fraction_of_Day*/
+
+/*------------------------------------------------------------------*/
+
+/* The function Calendar_Date converts a Julian Date to a struct tm.   */
+/* Only the members tm_year, tm_mon and tm_mday are calculated and set */
+void
+Calendar_Date(double jd, struct tm *cdate)
+{
+  /* Astronomical Formulae for Calculators, Jean Meeus, pages 26-27 */
+  int Z,month;
+  double A,B,C,D,E,F,alpha,day,year,factor;
+
+  factor = 0.5/secday/1000;
+  F = Frac(jd + 0.5);
+  if (F + factor >= 1.0)
+    {
+      jd = jd + factor;
+      F  = 0.0;
+    } /*if*/
+  Z = Round(jd);
+  if( Z < 2299161 )
+    A = Z;
+  else
+    {
+      alpha = Int((Z - 1867216.25)/36524.25);
+      A = Z + 1 + alpha - Int(alpha/4);
+    } /*else*/
+  B = A + 1524;
+  C = Int((B - 122.1)/365.25);
+  D = Int(365.25 * C);
+  E = Int((B - D)/30.6001);
+  day = B - D - Int(30.6001 * E) + F;
+
+  if( E < 13.5 )
+    month = Round(E - 1);
+  else
+    month = Round(E - 13);
+  if( month > 2.5 )
+    year = C - 4716;
+  else
+    year = C - 4715;
+
+  cdate->tm_year = (int) year;
+  cdate->tm_mon = month;
+  cdate->tm_mday = (int) floor(day);
+
+} /*Function Calendar_Date*/
+
+/*------------------------------------------------------------------*/
+
+/* Time_of_Day takes a Julian Date and calculates the clock time */
+/* portion of that date. Only tm_hour, tm_min and tm_sec are set */
+void
+Time_of_Day(double jd, struct tm *cdate)
+{
+  int hr,mn,sc;
+  double _time;
+
+  _time = Frac(jd - 0.5)*secday;
+  _time = Round(_time);
+  hr = floor(_time/3600.0);
+  _time = _time - 3600.0*hr;
+  if( hr == 24 ) hr = 0;
+  mn = floor(_time/60.0);
+  sc = _time - 60.0*mn;
+  cdate->tm_hour = hr;
+  cdate->tm_min = mn;
+  cdate->tm_sec = sc;
+
+} /*Function Time_of_Day*/
+
+/*------------------------------------------------------------------*/
+
+/* The function Julian_Date converts a standard calendar   */
+/* date and time to a Julian Date. The procedure Date_Time */
+/* performs the inverse of this function. */
+double
+Julian_Date(struct tm *cdate)
+{
+  double julian_date;
+
+  julian_date = Julian_Date_of_Year(cdate->tm_year) + 
+                DOY(cdate->tm_year,cdate->tm_mon,cdate->tm_mday) +
+                Fraction_of_Day(cdate->tm_hour,cdate->tm_min,cdate->tm_sec);
+
+  return( julian_date );
+} /*Function Julian_Date */
+
+/*------------------------------------------------------------------*/
+
+
+/*  Date_Time()
+ *
+ *  The function Date_Time() converts a Julian Date to
+ *  standard calendar date and time. The function
+ *  Julian_Date() performs the inverse of this function.
+ */
+
+void
+Date_Time(double julian_date, struct tm *cdate)
+{
+  time_t jtime;
+
+  jtime = (julian_date - 2440587.5)*86400.;
+  *cdate = *gmtime( &jtime );
+
+} /* End of Date_Time() */
+
+
+/*------------------------------------------------------------------*/
+
+/* The procedure Check_Date can be used as a check to see if a calendar    */
+/* date and time are valid. It works by first converting the calendar      */
+/* date and time to a Julian Date (which allows for irregularities, such   */
+/* as a time greater than 24 hours) and then converting back and comparing.*/
+int
+Check_Date(struct tm *cdate)
+{
+  double jt;
+  struct tm chkdate;
+
+  jt = Julian_Date(cdate);
+  Date_Time(jt, &chkdate);
+
+  if( (cdate->tm_year == chkdate.tm_year) &&
+      (cdate->tm_mon  == chkdate.tm_mon ) &&
+      (cdate->tm_mday == chkdate.tm_mday) &&
+      (cdate->tm_hour == chkdate.tm_hour) &&
+      (cdate->tm_min  == chkdate.tm_min ) &&
+      (cdate->tm_sec  == chkdate.tm_sec ) )
+    return ( 1 );
+  else
+    return( 0 );
+
+} /*Procedure Check_Date*/
+
+/*------------------------------------------------------------------*/
+
+/* Procedures Time_to_UTC and Time_from_UTC are used to  */
+/* convert 'struct tm' dates between UTC and local time. */
+/* The procedures JD_to_UTC and JD_from_UTC are used to  */
+/* do the same thing working directly with Julian dates. */
+
+struct tm
+Time_to_UTC(struct tm *cdate)
+{
+  time_t tdate;
+
+  tdate = mktime(cdate);
+  return( *gmtime(&tdate) );
+} /*Procedure Time_to_UTC*/
+
+/*------------------------------------------------------------------*/
+
+struct tm 
+Time_from_UTC(struct tm *cdate)
+{
+  time_t tdate;
+
+  tdate = mktime(cdate);
+  return( *localtime(&tdate) );
+} /*Procedure Time_from_UTC*/
+
+/*------------------------------------------------------------------*/
+
+/*
+ BSD systems don't define the timezone variable, so the following two
+ routines won't work.  They're not used anyway in the example main(),
+ so we might as well comment them out.  
+*/
+
+#if 0
+
+double
+JD_to_UTC(double jt)
+{
+  extern long timezone;
+  struct tm cdate;
+
+  time_t t = 0;
+
+  cdate = *localtime( &t );
+  jt = jt - timezone/secday;
+  if( cdate.tm_isdst )
+    jt= jt - 1.0/24.0;
+
+  return( jt );
+} /*Procedure JD_to_UTC*/
+
+/*------------------------------------------------------------------*/
+
+double
+JD_from_UTC(double jt)
+{
+  extern long timezone;
+  struct tm cdate;
+  time_t t = 0;
+
+  cdate = *localtime( &t );
+  jt = jt + timezone/secday;
+  if( cdate.tm_isdst )
+    jt= jt + 1.0/24.0;
+
+  return( jt );
+} /*Procedure JD_from_UTC*/
+
+#endif
+
+/*------------------------------------------------------------------*/
+
+/* The function Delta_ET has been added to allow calculations on   */
+/* the position of the sun.  It provides the difference between UT */
+/* (approximately the same as UTC) and ET (now referred to as TDT).*/
+/* This function is based on a least squares fit of data from 1950 */
+/* to 1991 and will need to be updated periodically. */
+
+double
+Delta_ET(double year)
+{
+  /* Values determined using data from 1950-1991 in the 1990 
+     Astronomical Almanac.  See DELTA_ET.WQ1 for details. */
+
+  double delta_et;
+
+  delta_et = 26.465 + 0.747622*(year - 1950) +
+             1.886913*sin(twopi*(year - 1975)/33);
+
+  return( delta_et );
+} /*Function Delta_ET*/
+
+/*------------------------------------------------------------------*/
+
+/* The function Julian_Date_of_Year calculates the Julian Date  */
+/* of Day 0.0 of {year}. This function is used to calculate the */
+/* Julian Date of any date by using Julian_Date_of_Year, DOY,   */
+/* and Fraction_of_Day. */
+
+double
+Julian_Date_of_Year(double year)
+{
+  /* Astronomical Formulae for Calculators, Jean Meeus, */
+  /* pages 23-25. Calculate Julian Date of 0.0 Jan year */
+
+  long A,B,i;
+  double jdoy;
+
+  year = year-1;
+  i = year/100;
+  A = i;
+  i = A/4;
+  B = 2-A+i;
+  i = 365.25*year;
+  i += 30.6001*14;
+  jdoy = i+1720994.5+B;
+
+  return (jdoy);
+}  /*Function Julian_Date_of_Year*/
+
+/*------------------------------------------------------------------*/
+
+/* The function ThetaG calculates the Greenwich Mean Sidereal Time */
+/* for an epoch specified in the format used in the NORAD two-line */
+/* element sets. It has now been adapted for dates beyond the year */
+/* 1999, as described above. The function ThetaG_JD provides the   */
+/* same calculation except that it is based on an input in the     */
+/* form of a Julian Date. */
+
+double
+ThetaG(double epoch, deep_arg_t *deep_arg)
+{
+/* Reference:  The 1992 Astronomical Almanac, page B6. */
+
+  double year,day,UT,jd,TU,GMST,_ThetaG;
+
+/* Modification to support Y2K */
+/* Valid 1957 through 2056     */
+  day = modf(epoch*1E-3,&year)*1E3;
+  if(year < 57)
+    year += 2000;
+  else
+    year += 1900;
+  /* End modification */
+
+  UT   = modf(day,&day);
+  jd   = Julian_Date_of_Year(year)+day;
+  TU   = (jd-2451545.0)/36525;
+  GMST = 24110.54841+TU*(8640184.812866+TU*(0.093104-TU* 6.2E-6));
+  GMST = Modulus(GMST+secday*omega_E*UT,secday);
+  _ThetaG = twopi*GMST/secday;
+  deep_arg->ds50 = jd-2433281.5+UT;
+  _ThetaG = FMod2p(6.3003880987*deep_arg->ds50+1.72944494);
+
+  return (_ThetaG);
+} /* Function ThetaG */
+
+/*------------------------------------------------------------------*/
+
+double
+ThetaG_JD(double jd)
+{
+/* Reference:  The 1992 Astronomical Almanac, page B6. */
+
+  double UT,TU,GMST;
+
+  UT   = Frac(jd + 0.5);
+  jd   = jd - UT;
+  TU   = (jd - 2451545.0)/36525;
+  GMST = 24110.54841 + TU * (8640184.812866 + TU * (0.093104 - TU * 6.2E-6));
+  GMST = Modulus(GMST + secday*omega_E*UT,secday);
+
+  return( twopi * GMST/secday );
+} /*Function ThetaG_JD*/
+
+/*------------------------------------------------------------------*/
+
+/* Gets calendar time from time() and produces a UTC calendar date */
+void
+UTC_Calendar_Now( struct tm *cdate )
+{
+  time_t t;
+
+  t = time(0);
+  *cdate = *gmtime(&t);
+  cdate->tm_year += 1900;
+  cdate->tm_mon += 1;
+
+} /* End UTC_Calendar_Now */
+/*------------------------------------------------------------------*/