/* jdn.c -- Julian Day Number computation
**
** public domain Julian Day Number functions
**
** Based on formulae originally posted by
** Tom Van Flandern / Washington, DC / metares@well.sf.ca.us
** in the UseNet newsgroup sci.astro.
** Reposted 14 May 1991 in FidoNet C Echo conference by
** Paul Schlyter (Stockholm)
** Minor corrections, added JDN to julian, and recast into C by
** Raymond Gardner Englewood, Colorado
**
** Synopsis:
** long ymd_to_jdn(int year, int month, int day, int julian_flag)
** void jdn_to_ymd(long jdn, int *year, int *month, int *day,
** int julian_flag)
** year is negative if BC
** if julian_flag is > 0, use Julian calendar
** if julian_flag is == 0, use Gregorian calendar
** if julian_flag is < 0, routines decide based on date
**
** These routines convert Gregorian and Julian calendar dates to and
** from Julian Day Numbers. Julian Day Numbers (JDN) are used by
** astronomers as a date/time measure independent of calendars and
** convenient for computing the elapsed time between dates. The JDN
** for any date/time is the number of days (including fractional
** days) elapsed since noon, 1 Jan 4713 BC. Julian Day Numbers were
** originated by Joseph Scaliger in 1582 and named after his father
** Julius, not after Julius Caesar. They are not related to the
** Julian calendar.
**
** For dates from 1 Jan 4713 BC thru 12 Dec Feb 32766 AD, ymd_to_jdn()
** will give the JDN for noon on that date. jdn_to_ymd() will compute
** the year, month, and day from the JDN. Years BC are given (and
** returned) as negative numbers. Note that there is no year 0 BC;
** the day before 1 Jan 1 AD is 31 Dec 1 BC. Note also that 1 BC,
** 5 BC, etc. are leap years.
**
** Pope Gregory XIII decreed that the Julian calendar would end on
** 4 Oct 1582 AD and that the next day would be 15 Oct 1582 in the
** Gregorian Calendar. The only other change is that centesimal
** years (years ending in 00) would no longer be leap years
** unless divisible by 400. Britain and its possessions and
** colonies continued to use the Julian calendar up until 2 Sep
** 1752, when the next day became 14 Sep 1752 in the Gregorian
** Calendar. These routines can be compiled to use either
** convention. By default, the British convention will be used.
** Simply #define PAPAL to use Pope Gregory's convention.
**
** Each routine takes, as its last argument, a flag to indicate
** whether to use the Julian or Gregorian calendar convention. If
** this flag is negative, the routines decide based on the date
** itself, using the changeover date described in the preceding
** paragraph. If the flag is zero, Gregorian conventions will be used,
** and if the flag is positive, Julian conventions will be used.
*/
#ifdef PAPAL /* Pope Gregory XIII's decree */
#define LASTJULDATE 15821004L /* last day to use Julian calendar */
#define LASTJULJDN 2299160L /* jdn of same */
#else /* British-American usage */
#define LASTJULDATE 17520902L /* last day to use Julian calendar */
#define LASTJULJDN 2361221L /* jdn of same */
#endif
long ymd_to_jdn(int y, int m, int d, int julian)
{
long jdn;
if (julian < 0) /* set Julian flag if auto set */
julian = (((y * 100L) + m) * 100 + d <= LASTJULDATE);
if (y < 0) /* adjust BC year */
y++;
if (julian)
jdn = 367L * y - 7 * (y + 5001L + (m - 9) / 7) / 4
+ 275 * m / 9 + d + 1729777L;
else
jdn = (long)(d - 32076)
+ 1461L * (y + 4800L + (m - 14) / 12) / 4
+ 367 * (m - 2 - (m - 14) / 12 * 12) / 12
- 3 * ((y + 4900L + (m - 14) / 12) / 100) / 4
+ 1; /* correction by rdg */
return jdn;
}
void jdn_to_ymd(long jdn, int *yy, int *mm, int *dd, int julian)
{
long x, z, m, d, y;
long daysPer400Years = 146097L;
long fudgedDaysPer4000Years = 1460970L + 31;
if (julian < 0) /* set Julian flag if auto set */
julian = (jdn <= LASTJULJDN);
x = jdn + 68569L;
if ( julian )
{
x += 38;
daysPer400Years = 146100L;
fudgedDaysPer4000Years = 1461000L + 1;
}
z = 4 * x / daysPer400Years;
x = x - (daysPer400Years * z + 3) / 4;
y = 4000 * (x + 1) / fudgedDaysPer4000Years;
x = x - 1461 * y / 4 + 31;
m = 80 * x / 2447;
d = x - 2447 * m / 80;
x = m / 11;
m = m + 2 - 12 * x;
y = 100 * (z - 49) + y + x;
*yy = (int)y;
*mm = (int)m;
*dd = (int)d;
if (*yy <= 0) /* adjust BC years */
(*yy)--;
}