yodaUtility/sgp4/cJulian.cpp

//
// cJulian.cpp
//
// This class encapsulates Julian dates with the epoch of 12:00 noon (12:00 UT)
// on January 1, 4713 B.C. Some epoch dates:
//    01/01/1990 00:00 UTC - 2447892.5
//    01/01/1990 12:00 UTC - 2447893.0
//    01/01/2000 00:00 UTC - 2451544.5
//    01/01/2001 00:00 UTC - 2451910.5
//
// Note the Julian day begins at noon, which allows astronomers to have all
// the dates in a single observing session the same.
//
// References:
// "Astronomical Formulae for Calculators", Jean Meeus
// "Satellite Communications", Dennis Roddy, 2nd Edition, 1995.
//
// Copyright (c) 2003 Michael F. Henry
//
// mfh 12/24/2003
//
#include "stdafx.h"
#include "globals.h"
#include "cJulian.h"
#include "math.h"
#include "time.h"

//////////////////////////////////////////////////////////////////////////////
// Create a Julian date object from a time_t object. time_t objects store the
// number of seconds since midnight UTC January 1, 1970.
cJulian::cJulian(time_t time)
{
   struct tm *ptm = gmtime(&time);
   assert(ptm);

   int    year = ptm->tm_year + 1900;
   double day  = ptm->tm_yday + 1 +
                 (ptm->tm_hour + 
                  ((ptm->tm_min + 
                   (ptm->tm_sec / 60.0)) / 60.0)) / 24.0;

   Initialize(year, day);
}

//////////////////////////////////////////////////////////////////////////////
// Create a Julian date object from a year and day of year.
// Example parameters: year = 2001, day = 1.5 (Jan 1 12h)
cJulian::cJulian(int year, double day)
{
   Initialize(year, day);
}

//////////////////////////////////////////////////////////////////////////////
// Create a Julian date object.
cJulian::cJulian(int year,               // i.e., 2004
                 int mon,                // 1..12
                 int day,                // 1..31
                 int hour,               // 0..23
                 int min,                // 0..59
                 double sec /* = 0.0 */) // 0..(59.999999...)

{
   // Calculate N, the day of the year (1..366)
   int N;
   int F1 = (int)((275.0 * mon) / 9.0);
   int F2 = (int)((mon + 9.0) / 12.0);

   if (IsLeapYear(year))
   {
      // Leap year
      N = F1 - F2 + day - 30;
   }
   else
   {
      // Common year
      N = F1 - (2 * F2) + day - 30;
   }
   
   double dblDay = N + (hour + (min + (sec / 60.0)) / 60.0) / 24.0;

   Initialize(year, dblDay);
}

//////////////////////////////////////////////////////////////////////////////
void cJulian::Initialize(int year, double day)
{
   // 1582 A.D.: 10 days removed from calendar
   // 3000 A.D.: Arbitrary error checking limit
   assert((year > 1582) && (year < 3000));
   assert((day >= 0.0) && (day <= 366.5));

   // Now calculate Julian date

   year--;

   // Centuries are not leap years unless they divide by 400
   int A = (year / 100);
   int B = 2 - A + (A / 4);

   double NewYears = (int)(365.25 * year) +
                     (int)(30.6001 * 14)  + 
                     1720994.5 + B;  // 1720994.5 = Oct 30, year -1

   m_Date = NewYears + day;
}

//////////////////////////////////////////////////////////////////////////////
// getComponent()
// Return requested components of date.
// Year : Includes the century.
// Month: 1..12
// Day  : 1..31 including fractional part
void cJulian::getComponent(int    *pYear, 
                           int    *pMon  /* = NULL */,
                           double *pDOM  /* = NULL */) const
{
   assert(pYear != NULL);

   double jdAdj = getDate() + 0.5;
   int    Z     = (int)jdAdj;  // integer part
   double F     = jdAdj - Z;   // fractional part
   double alpha = (int)((Z - 1867216.25) / 36524.25);
   double A     = Z + 1 + alpha - (int)(alpha / 4.0);
   double B     = A + 1524.0;
   int    C     = (int)((B - 122.1) / 365.25);
   int    D     = (int)(C * 365.25);
   int    E     = (int)((B - D) / 30.6001);

   double DOM   = B - D - (int)(E * 30.6001) + F;
   int    month = (E < 13.5) ? (E - 1) : (E - 13);
   int    year  = (month > 2.5) ? (C - 4716) : (C - 4715); 
   
   *pYear = year;

   if (pMon != NULL)
      *pMon = month;

   if (pDOM != NULL)
      *pDOM = DOM;
} 

//////////////////////////////////////////////////////////////////////////////
// toGMST()
// Calculate Greenwich Mean Sidereal Time for the Julian date. The return value
// is the angle, in radians, measuring eastward from the Vernal Equinox to the
// prime meridian. This angle is also referred to as "ThetaG" (Theta GMST).
// 
// References:
//    The 1992 Astronomical Almanac, page B6.
//    Explanatory Supplement to the Astronomical Almanac, page 50.
//    Orbital Coordinate Systems, Part III, Dr. T.S. Kelso, Satellite Times,
//       Nov/Dec 1995
double cJulian::toGMST() const
{
   const double UT = fmod(m_Date + 0.5, 1.0);
   const double TU = (FromJan1_12h_2000() - UT) / 36525.0;

   double GMST = 24110.54841 + TU * 
                 (8640184.812866 + TU * (0.093104 - TU * 6.2e-06));

   GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
   
   if (GMST < 0.0)
      GMST += SEC_PER_DAY;  // "wrap" negative modulo value

   return  (TWOPI * (GMST / SEC_PER_DAY));
 }

//////////////////////////////////////////////////////////////////////////////
// toLMST()
// Calculate Local Mean Sidereal Time for given longitude (for this date).
// The longitude is assumed to be in radians measured west from Greenwich.
// The return value is the angle, in radians, measuring eastward from the
// Vernal Equinox to the given longitude.
double cJulian::toLMST(double lon) const
{
   return fmod(toGMST() + lon, TWOPI);
}

//////////////////////////////////////////////////////////////////////////////
// toTime()
// Convert to type time_t
// Avoid using this function as it discards the fractional seconds of the
// time component.
time_t cJulian::toTime() const
{
   int    nYear;
   int    nMonth;
   double dblDay;

   getComponent(&nYear, &nMonth, &dblDay);

   // dblDay is the fractional Julian Day (i.e., 29.5577).
   // Save the whole number day in nDOM and convert dblDay to
   // the fractional portion of day.
   int nDOM = (int)dblDay;

   dblDay -= nDOM;

   const int SEC_PER_MIN = 60;
   const int SEC_PER_HR  = 60 * SEC_PER_MIN;
   const int SEC_PER_DAY = 24 * SEC_PER_HR;

   int secs = (int)((dblDay * SEC_PER_DAY) + 0.5);

   // Create a "struct tm" type. 
   // NOTE:
   // The "struct tm" type has a 1-second resolution. Any fractional
   // component of the "seconds" time value is discarded.
   struct tm tGMT;
   memset(&tGMT, 0, sizeof(tGMT));

   tGMT.tm_year = nYear - 1900;  // 2001 is 101
   tGMT.tm_mon  = nMonth - 1;    // January is 0
   tGMT.tm_mday = nDOM;          // First day is 1
   tGMT.tm_hour = secs / SEC_PER_HR;
   tGMT.tm_min  = (secs % SEC_PER_HR) / SEC_PER_MIN;
   tGMT.tm_sec  = (secs % SEC_PER_HR) % SEC_PER_MIN;
   tGMT.tm_isdst = 0; // No conversion desired

   time_t tEpoch = mktime(&tGMT);

   if (tEpoch != -1)
   {
      // Valid time_t value returned from mktime().
      // mktime() expects a local time which means that tEpoch now needs 
      // to be adjusted by the difference between this time zone and GMT.
      tEpoch -= timezone;
   }

   return tEpoch;
}
1	//
2	// cJulian.cpp
3	//
4	// This class encapsulates Julian dates with the epoch of 12:00 noon (12:00 UT)
5	// on January 1, 4713 B.C. Some epoch dates:
6	// 01/01/1990 00:00 UTC - 2447892.5
7	// 01/01/1990 12:00 UTC - 2447893.0
8	// 01/01/2000 00:00 UTC - 2451544.5
9	// 01/01/2001 00:00 UTC - 2451910.5
10	//
11	// Note the Julian day begins at noon, which allows astronomers to have all
12	// the dates in a single observing session the same.
13	//
14	// References:
15	// "Astronomical Formulae for Calculators", Jean Meeus
16	// "Satellite Communications", Dennis Roddy, 2nd Edition, 1995.
17	//
18	// Copyright (c) 2003 Michael F. Henry
19	//
20	// mfh 12/24/2003
21	//
22	#include "stdafx.h"
23	#include "globals.h"
24	#include "cJulian.h"
25	#include "math.h"
26	#include "time.h"
27
28	//////////////////////////////////////////////////////////////////////////////
29	// Create a Julian date object from a time_t object. time_t objects store the
30	// number of seconds since midnight UTC January 1, 1970.
31	cJulian::cJulian(time_t time)
32	{
33	struct tm *ptm = gmtime(&time);
34	assert(ptm);
35
36	int year = ptm->tm_year + 1900;
37	double day = ptm->tm_yday + 1 +
38	(ptm->tm_hour +
39	((ptm->tm_min +
40	(ptm->tm_sec / 60.0)) / 60.0)) / 24.0;
41
42	Initialize(year, day);
43	}
44
45	//////////////////////////////////////////////////////////////////////////////
46	// Create a Julian date object from a year and day of year.
47	// Example parameters: year = 2001, day = 1.5 (Jan 1 12h)
48	cJulian::cJulian(int year, double day)
49	{
50	Initialize(year, day);
51	}
52
53	//////////////////////////////////////////////////////////////////////////////
54	// Create a Julian date object.
55	cJulian::cJulian(int year, // i.e., 2004
56	int mon, // 1..12
57	int day, // 1..31
58	int hour, // 0..23
59	int min, // 0..59
60	double sec /* = 0.0 */) // 0..(59.999999...)
61
62	{
63	// Calculate N, the day of the year (1..366)
64	int N;
65	int F1 = (int)((275.0 * mon) / 9.0);
66	int F2 = (int)((mon + 9.0) / 12.0);
67
68	if (IsLeapYear(year))
69	{
70	// Leap year
71	N = F1 - F2 + day - 30;
72	}
73	else
74	{
75	// Common year
76	N = F1 - (2 * F2) + day - 30;
77	}
78
79	double dblDay = N + (hour + (min + (sec / 60.0)) / 60.0) / 24.0;
80
81	Initialize(year, dblDay);
82	}
83
84	//////////////////////////////////////////////////////////////////////////////
85	void cJulian::Initialize(int year, double day)
86	{
87	// 1582 A.D.: 10 days removed from calendar
88	// 3000 A.D.: Arbitrary error checking limit
89	assert((year > 1582) && (year < 3000));
90	assert((day >= 0.0) && (day <= 366.5));
91
92	// Now calculate Julian date
93
94	year--;
95
96	// Centuries are not leap years unless they divide by 400
97	int A = (year / 100);
98	int B = 2 - A + (A / 4);
99
100	double NewYears = (int)(365.25 * year) +
101	(int)(30.6001 * 14) +
102	1720994.5 + B; // 1720994.5 = Oct 30, year -1
103
104	m_Date = NewYears + day;
105	}
106
107	//////////////////////////////////////////////////////////////////////////////
108	// getComponent()
109	// Return requested components of date.
110	// Year : Includes the century.
111	// Month: 1..12
112	// Day : 1..31 including fractional part
113	void cJulian::getComponent(int *pYear,
114	int pMon / = NULL */,
115	double pDOM / = NULL */) const
116	{
117	assert(pYear != NULL);
118
119	double jdAdj = getDate() + 0.5;
120	int Z = (int)jdAdj; // integer part
121	double F = jdAdj - Z; // fractional part
122	double alpha = (int)((Z - 1867216.25) / 36524.25);
123	double A = Z + 1 + alpha - (int)(alpha / 4.0);
124	double B = A + 1524.0;
125	int C = (int)((B - 122.1) / 365.25);
126	int D = (int)(C * 365.25);
127	int E = (int)((B - D) / 30.6001);
128
129	double DOM = B - D - (int)(E * 30.6001) + F;
130	int month = (E < 13.5) ? (E - 1) : (E - 13);
131	int year = (month > 2.5) ? (C - 4716) : (C - 4715);
132
133	*pYear = year;
134
135	if (pMon != NULL)
136	*pMon = month;
137
138	if (pDOM != NULL)
139	*pDOM = DOM;
140	}
141
142	//////////////////////////////////////////////////////////////////////////////
143	// toGMST()
144	// Calculate Greenwich Mean Sidereal Time for the Julian date. The return value
145	// is the angle, in radians, measuring eastward from the Vernal Equinox to the
146	// prime meridian. This angle is also referred to as "ThetaG" (Theta GMST).
147	//
148	// References:
149	// The 1992 Astronomical Almanac, page B6.
150	// Explanatory Supplement to the Astronomical Almanac, page 50.
151	// Orbital Coordinate Systems, Part III, Dr. T.S. Kelso, Satellite Times,
152	// Nov/Dec 1995
153	double cJulian::toGMST() const
154	{
155	const double UT = fmod(m_Date + 0.5, 1.0);
156	const double TU = (FromJan1_12h_2000() - UT) / 36525.0;
157
158	double GMST = 24110.54841 + TU *
159	(8640184.812866 + TU * (0.093104 - TU * 6.2e-06));
160
161	GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
162
163	if (GMST < 0.0)
164	GMST += SEC_PER_DAY; // "wrap" negative modulo value
165
166	return (TWOPI * (GMST / SEC_PER_DAY));
167	}
168
169	//////////////////////////////////////////////////////////////////////////////
170	// toLMST()
171	// Calculate Local Mean Sidereal Time for given longitude (for this date).
172	// The longitude is assumed to be in radians measured west from Greenwich.
173	// The return value is the angle, in radians, measuring eastward from the
174	// Vernal Equinox to the given longitude.
175	double cJulian::toLMST(double lon) const
176	{
177	return fmod(toGMST() + lon, TWOPI);
178	}
179
180	//////////////////////////////////////////////////////////////////////////////
181	// toTime()
182	// Convert to type time_t
183	// Avoid using this function as it discards the fractional seconds of the
184	// time component.
185	time_t cJulian::toTime() const
186	{
187	int nYear;
188	int nMonth;
189	double dblDay;
190
191	getComponent(&nYear, &nMonth, &dblDay);
192
193	// dblDay is the fractional Julian Day (i.e., 29.5577).
194	// Save the whole number day in nDOM and convert dblDay to
195	// the fractional portion of day.
196	int nDOM = (int)dblDay;
197
198	dblDay -= nDOM;
199
200	const int SEC_PER_MIN = 60;
201	const int SEC_PER_HR = 60 * SEC_PER_MIN;
202	const int SEC_PER_DAY = 24 * SEC_PER_HR;
203
204	int secs = (int)((dblDay * SEC_PER_DAY) + 0.5);
205
206	// Create a "struct tm" type.
207	// NOTE:
208	// The "struct tm" type has a 1-second resolution. Any fractional
209	// component of the "seconds" time value is discarded.
210	struct tm tGMT;
211	memset(&tGMT, 0, sizeof(tGMT));
212
213	tGMT.tm_year = nYear - 1900; // 2001 is 101
214	tGMT.tm_mon = nMonth - 1; // January is 0
215	tGMT.tm_mday = nDOM; // First day is 1
216	tGMT.tm_hour = secs / SEC_PER_HR;
217	tGMT.tm_min = (secs % SEC_PER_HR) / SEC_PER_MIN;
218	tGMT.tm_sec = (secs % SEC_PER_HR) % SEC_PER_MIN;
219	tGMT.tm_isdst = 0; // No conversion desired
220
221	time_t tEpoch = mktime(&tGMT);
222
223	if (tEpoch != -1)
224	{
225	// Valid time_t value returned from mktime().
226	// mktime() expects a local time which means that tEpoch now needs
227	// to be adjusted by the difference between this time zone and GMT.
228	tEpoch -= timezone;
229	}
230
231	return tEpoch;
232	}