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	kusanagi	1.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			}