sgp_time.c
1 /* 2 * Unit SGP_Time 3 * Author: Dr TS Kelso 4 * Original Version: 1992 Jun 02 5 * Current Revision: 2000 Jan 22 6 * Modified for Y2K: 1999 Mar 07 7 * Version: 2.05 8 * Copyright: 1992-1999, All Rights Reserved 9 * Version 1.50 added Y2K support. Due to limitations in the current 10 * format of the NORAD two-line element sets, however, only dates 11 * through 2056 December 31/2359 UTC are valid. 12 * Version 1.60 modifies Calendar_Date to ensure date matches time 13 * resolution and modifies Time_of_Day to make it more robust. 14 * Version 2.00 adds Julian_Date, Date_Time, and Check_Date to support 15 * checking for valid date/times, permitting the use of Time_to_UTC and 16 * Time_from_UTC for UTC/local time conversions. 17 * Version 2.05 modifies UTC_offset to allow non-integer offsets. 18 * 19 * Ported to C by: Neoklis Kyriazis April 9 2001 20 */ 21 #define _POSIX_C_SOURCE 1 // gmtime_r() 22 #include <time.h> 23 #include "sgp4sdp4.h" 24 25 /* The function Julian_Date_of_Epoch returns the Julian Date of */ 26 /* an epoch specified in the format used in the NORAD two-line */ 27 /* element sets. It has been modified to support dates beyond */ 28 /* the year 1999 assuming that two-digit years in the range 00-56 */ 29 /* correspond to 2000-2056. Until the two-line element set format */ 30 /* is changed, it is only valid for dates through 2056 December 31. */ 31 double Julian_Date_of_Epoch(double epoch) 32 { 33 double year, day; 34 35 /* Modification to support Y2K */ 36 /* Valid 1957 through 2056 */ 37 day = modf(epoch * 1E-3, &year) * 1E3; 38 if (year < 57) 39 year = year + 2000; 40 else 41 year = year + 1900; 42 /* End modification */ 43 44 return (Julian_Date_of_Year(year) + day); 45 } 46 47 /* Converts a Julian epoch to NORAD TLE epoch format */ 48 double Epoch_Time(double jd) 49 { 50 double yr, _time, epoch_time; 51 struct tm edate; 52 53 Calendar_Date(jd, &edate); 54 yr = edate.tm_year - 100 * (edate.tm_year / 100); 55 _time = Frac(jd + 0.5); 56 epoch_time = yr * 1000 57 + DOY(edate.tm_year, edate.tm_mon, edate.tm_mday) + _time; 58 59 return (epoch_time); 60 } 61 62 /* The function DOY calculates the day of the year for the specified */ 63 /* date. The calculation uses the rules for the Gregorian calendar */ 64 /* and is valid from the inception of that calendar system. */ 65 int DOY(int yr, int mo, int dy) 66 { 67 const int days[] = 68 { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 }; 69 int i, day; 70 71 day = 0; 72 for (i = 0; i < mo - 1; i++) 73 day += days[i]; 74 day = day + dy; 75 76 /* Leap year correction */ 77 if ((yr % 4 == 0) && ((yr % 100 != 0) || (yr % 400 == 0)) && (mo > 2)) 78 day++; 79 80 return (day); 81 } 82 83 /* Fraction_of_Day calculates the fraction of */ 84 /* a day passed at the specified input time. */ 85 double Fraction_of_Day(int hr, int mi, int se) 86 { 87 return ((hr + (mi + se / 60.0) / 60.0) / 24.0); 88 } 89 90 /* The function Calendar_Date converts a Julian Date to a struct tm. */ 91 /* Only the members tm_year, tm_mon and tm_mday are calculated and set */ 92 void Calendar_Date(double jd, struct tm *cdate) 93 { 94 /* Astronomical Formulae for Calculators, Jean Meeus, pages 26-27 */ 95 int Z, month; 96 double A, B, C, D, E, F, alpha, day, year, factor; 97 98 factor = 0.5 / secday / 1000; 99 F = Frac(jd + 0.5); 100 if (F + factor >= 1.0) 101 { 102 jd = jd + factor; 103 F = 0.0; 104 } /*if */ 105 Z = Round(jd); 106 if (Z < 2299161) 107 A = Z; 108 else 109 { 110 alpha = Int((Z - 1867216.25) / 36524.25); 111 A = Z + 1 + alpha - Int(alpha / 4); 112 } /*else */ 113 B = A + 1524; 114 C = Int((B - 122.1) / 365.25); 115 D = Int(365.25 * C); 116 E = Int((B - D) / 30.6001); 117 day = B - D - Int(30.6001 * E) + F; 118 119 if (E < 13.5) 120 month = Round(E - 1); 121 else 122 month = Round(E - 13); 123 if (month > 2.5) 124 year = C - 4716; 125 else 126 year = C - 4715; 127 128 cdate->tm_year = (int)year; 129 cdate->tm_mon = month; 130 cdate->tm_mday = (int)floor(day); 131 132 } 133 134 /* Time_of_Day takes a Julian Date and calculates the clock time */ 135 /* portion of that date. Only tm_hour, tm_min and tm_sec are set */ 136 void Time_of_Day(double jd, struct tm *cdate) 137 { 138 int hr, mn, sc; 139 double _time; 140 141 _time = Frac(jd - 0.5) * secday; 142 _time = Round(_time); 143 hr = floor(_time / 3600.0); 144 _time = _time - 3600.0 * hr; 145 if (hr == 24) 146 hr = 0; 147 mn = floor(_time / 60.0); 148 sc = _time - 60.0 * mn; 149 cdate->tm_hour = hr; 150 cdate->tm_min = mn; 151 cdate->tm_sec = sc; 152 153 } 154 155 /* The function Julian_Date converts a standard calendar */ 156 /* date and time (FIXME: GMT?) to a Julian Date. The procedure Date_Time */ 157 /* performs the inverse of this function. */ 158 double Julian_Date(struct tm *cdate) 159 { 160 double julian_date; 161 162 julian_date = Julian_Date_of_Year(cdate->tm_year) + 163 DOY(cdate->tm_year, cdate->tm_mon, cdate->tm_mday) + 164 Fraction_of_Day(cdate->tm_hour, cdate->tm_min, cdate->tm_sec); 165 166 return (julian_date); 167 } 168 169 /* Date_Time() 170 * 171 * The function Date_Time() converts a Julian Date to 172 * standard calendar date and time (GMT). The function 173 * Julian_Date() performs the inverse of this function. 174 */ 175 176 void Date_Time(double julian_date, struct tm *cdate) 177 { 178 time_t jtime; 179 180 jtime = (julian_date - 2440587.5) * 86400.; 181 *cdate = *gmtime(&jtime); 182 183 } 184 185 /* The procedure Check_Date can be used as a check to see if a calendar */ 186 /* date and time are valid. It works by first converting the calendar */ 187 /* date and time to a Julian Date (which allows for irregularities, such */ 188 /* as a time greater than 24 hours) and then converting back and comparing.*/ 189 int Check_Date(struct tm *cdate) 190 { 191 double jt; 192 struct tm chkdate; 193 194 jt = Julian_Date(cdate); 195 Date_Time(jt, &chkdate); 196 197 if ((cdate->tm_year == chkdate.tm_year) && 198 (cdate->tm_mon == chkdate.tm_mon) && 199 (cdate->tm_mday == chkdate.tm_mday) && 200 (cdate->tm_hour == chkdate.tm_hour) && 201 (cdate->tm_min == chkdate.tm_min) && (cdate->tm_sec == chkdate.tm_sec)) 202 return (1); 203 else 204 return (0); 205 206 } 207 208 /* Procedure gmtime_r is used to convert the calendar */ 209 /* time time_t to a broken-down time representation, */ 210 /* expressed in Coordinated Universal Time (UTC). */ 211 /* This function is not present in WIN32 */ 212 213 #ifdef WIN32 214 struct tm *gmtime_r(const time_t * timer, struct tm *result) 215 { 216 struct tm *local_result; 217 218 local_result = gmtime(timer); 219 220 if (local_result == NULL || result == NULL) 221 return NULL; 222 223 return memcpy(result, local_result, sizeof(*result)); 224 225 } 226 #endif 227 228 /* Procedures Time_to_UTC and Time_from_UTC are used to */ 229 /* convert 'struct tm' dates between UTC and local time. */ 230 /* The procedures JD_to_UTC and JD_from_UTC are used to */ 231 /* do the same thing working directly with Julian dates. */ 232 233 void Time_to_UTC(struct tm *cdate, struct tm *odate) 234 { 235 time_t tdate; 236 237 /* 238 functions such as Julian_Date work with tm_year in being 239 the year since 0 AD and mktime uses it as the year since 1900 240 241 tm_isdst = -1 forces the mktime call to lookup the daylight 242 savings for the current timezone instead of using whatever 243 was in the structure when it was created. 244 */ 245 cdate->tm_year -= 1900; 246 cdate->tm_isdst = -1; 247 tdate = mktime(cdate); 248 gmtime_r(&tdate, odate); 249 odate->tm_year += 1900; 250 } 251 252 struct tm Time_from_UTC(struct tm *cdate) 253 { 254 time_t tdate; 255 256 tdate = mktime(cdate); 257 return (*localtime(&tdate)); 258 } 259 260 /* The function Delta_ET has been added to allow calculations on */ 261 /* the position of the sun. It provides the difference between UT */ 262 /* (approximately the same as UTC) and ET (now referred to as TDT).*/ 263 /* This function is based on a least squares fit of data from 1950 */ 264 /* to 1991 and will need to be updated periodically. */ 265 double Delta_ET(double year) 266 { 267 /* Values determined using data from 1950-1991 in the 1990 268 Astronomical Almanac. See DELTA_ET.WQ1 for details. */ 269 270 double delta_et; 271 272 delta_et = 26.465 + 0.747622 * (year - 1950) + 273 1.886913 * sin(twopi * (year - 1975) / 33); 274 275 return (delta_et); 276 } 277 278 /* The function Julian_Date_of_Year calculates the Julian Date */ 279 /* of Day 0.0 of {year}. This function is used to calculate the */ 280 /* Julian Date of any date by using Julian_Date_of_Year, DOY, */ 281 /* and Fraction_of_Day. */ 282 double Julian_Date_of_Year(double year) 283 { 284 /* Astronomical Formulae for Calculators, Jean Meeus, */ 285 /* pages 23-25. Calculate Julian Date of 0.0 Jan year */ 286 287 long A, B, i; 288 double jdoy; 289 290 year = year - 1; 291 i = year / 100; 292 A = i; 293 i = A / 4; 294 B = 2 - A + i; 295 i = 365.25 * year; 296 i += 30.6001 * 14; 297 jdoy = i + 1720994.5 + B; 298 299 return (jdoy); 300 } 301 302 /* The function ThetaG calculates the Greenwich Mean Sidereal Time */ 303 /* for an epoch specified in the format used in the NORAD two-line */ 304 /* element sets. It has now been adapted for dates beyond the year */ 305 /* 1999, as described above. The function ThetaG_JD provides the */ 306 /* same calculation except that it is based on an input in the */ 307 /* form of a Julian Date. */ 308 double ThetaG(double epoch, deep_arg_t * deep_arg) 309 { 310 /* Reference: The 1992 Astronomical Almanac, page B6. */ 311 312 double year, day, UT, jd, TU, GMST, _ThetaG; 313 314 /* Modification to support Y2K */ 315 /* Valid 1957 through 2056 */ 316 day = modf(epoch * 1E-3, &year) * 1E3; 317 if (year < 57) 318 year += 2000; 319 else 320 year += 1900; 321 /* End modification */ 322 323 UT = modf(day, &day); 324 jd = Julian_Date_of_Year(year) + day; 325 TU = (jd - 2451545.0) / 36525; 326 GMST = 24110.54841 + TU * (8640184.812866 + TU * (0.093104 - TU * 6.2E-6)); 327 GMST = Modulus(GMST + secday * omega_E * UT, secday); 328 _ThetaG = twopi * GMST / secday; 329 deep_arg->ds50 = jd - 2433281.5 + UT; 330 _ThetaG = FMod2p(6.3003880987 * deep_arg->ds50 + 1.72944494); 331 332 return (_ThetaG); 333 } 334 335 double ThetaG_JD(double jd) 336 { 337 /* Reference: The 1992 Astronomical Almanac, page B6. */ 338 339 double UT, TU, GMST; 340 341 UT = Frac(jd + 0.5); 342 jd = jd - UT; 343 TU = (jd - 2451545.0) / 36525; 344 GMST = 24110.54841 + TU * (8640184.812866 + TU * (0.093104 - TU * 6.2E-6)); 345 GMST = Modulus(GMST + secday * omega_E * UT, secday); 346 347 return (twopi * GMST / secday); 348 } 349 350 /* Gets calendar time from time() and produces a UTC calendar date */ 351 void UTC_Calendar_Now(struct tm *cdate) 352 { 353 time_t t; 354 355 t = time(0); 356 *cdate = *gmtime(&t); 357 cdate->tm_year += 1900; 358 cdate->tm_mon += 1; 359 360 }