2.1 Gregorian or Julian date to Julian Day Number - Algorithm Development

Algorithm development begins by choosing a base year. Converting from a Gregorian calendar date implies that the base year is relative to the year zero. If the given date is relative to some particular year in the Gregorian calendar, then that date becomes our base year.

The month number and year number will be converted so that January and February are considered months 13 and 14 rather than months 1 and 2. January and February of year Y are converted so that they occur in the year Y - 1. In other words, a system is used where the first month of the year is March (designated as month 3) and the last month of the year is February of the following year and designated as month 14. The purpose of this scheme is to cause the extra day in a leap year to occur at the end of the "year." As will be seen, this simplifies the calculation of the number of leap years from year zero to a particular date and it simplifies the function which relates month number to the number of days in the preceding months.

Since the base year is March 1 of year 0, the Julian day number is determined relative to that date. From table 2, it can be seen that March 1.0 of year 0 is JD 1721119.5, assuming the date is in the Gregorian calendar. Thus March 1.5 of year 0 is JD 1721120.0 and is the Julian Day Number if a date of March 1 of year 0 is specified. Thus 1721118.5 must be added to the day number (of one) to get the correct JD. Note, however, that dates in the Gregorian calendar are assumed even though the base date is before the Gregorian calendar was invented. This works fine as long as the dates input to the algorithm are Gregorian calendar dates or dates in the Gregorian proleptic calendar.

Next, the days in the year are accounted for. This is easily done using the following function:

Y*365 +[Y/4] -[Y/100] + [Y/400]

This expression directly encodes the rule that there are 365 days in a year plus an extra day for leap years, where leap years occur every 4 years except for century years unless they are also divisible by 400. Here [x] is the greatest integer function, that is, [x] is the greatest integer that is not greater than x. Thus [2]=2, [2.5]=2, [0]=0, [-1.5]=-2, and [-3]=-3. For positive numbers, the greatest integer is the same as a function which truncates any fractional part of a number.

Note that although 0 is a leap year, the number of days relative to March 1 of the year zero are being considered. Relative to March 1 of the year zero, the first leap year is year 4 and February 29, 4, is part of the year beginning March 1, 3. The expression - [Y/100] eliminates all century years as leap years while [Y/400] adds back century years that are divisible by 400.

Finally, the number of days in the months prior to the specified month are accounted for. The following gives the range and domain of the function explicitly:

monthdays in previous months

Since, all the months except February are either 30 or 31 days long, one might suspect that a linear approximation can be used and the fractional part of the value removed to generate a convenient form of the function. One can not be certain that this will work on a particular function before it is examined in more detail, but, as will be shown, the function given in the table can be successfully generated in this way.

The function will take the form:
X = month number,
M = the slope of the line,
B= the Y intercept of the line, and
FIX() is the truncation function.

Now since FIX() removes the fraction, the values of MX+B that produce the same FIX function value are values of MX + B from some integer, say Y, up to, but not including, the value Y+1. The values for each Y can be graphed, but because of the scale, it is difficult to see precisely how the line might fit through all the function values or whether or not it is even possible to do so. However, it isn't difficult to write a computer program that determines what equations can pass through these intervals. Doing so shows that the line with maximum slope goes from point (4,31) to (12,276), has slope= 30.625 and Y intercept=-91.5. The line with minimum slope goes from point (14,337) to (7,123), has a slope= 30.57142857142857 and Y intercept=-91.0000000000000. These two lines cross at (9.333333333333067,194.3333333333252) and I call this the pivot point, since there will be lines within the specified intervals that go through this point and take on any of the slopes from 30.57142857142857 to 30.625. Given slope M, the values 194.3333333333252 - M * 9.333333333333067 can be used as the Y intercept. However, there will usually be many other Y intercept values that work as well.

From a computational point of view, one probably notices that a slope of 30.6 contains a minimum number of digits and shouldn't be too difficult to produce using integer arithmetic. If for slope = 30.6, the line goes through the pivot, the Y intercept is -91.26666666666669 but there are lots of other Y intercept values that work as well. These Y intercept values can range from -91.4 to -91.2, although the value of -91.2 can not actually be taken since the line would go through a couple of points that are integers and the truncation function would not produce the desired result. Thus the function


can be used to produce our desired function values. The function can also be converted to other forms for convenience. For example

FIX(30.6*X-91.4) is equivalent to (153*m-457)\5 since 153/5=30.6 and -457/5=-91.4.

Other coefficients for our generating function are also possible, and they can be searched for using the bounds described above. Suppose, for example, that the function is to have the form (A*X+B)\32 so that shifts can be used in the binary representation to implement integer division. Since the slope must range from 30.57142857142857 to 30.625, A/32 must also be in this range. This implies that A must range from 30.57142857142857 * 32 = 978.28571428571424 to 30.625*32= 980.0. The integer 979 is an obvious choice. This means the slope is 979.0/32=30.59375. This slope implies that the Y intercept is in the range -91.15625 to -91.3125. Thus B must be in the range of -91.15625*32= -2,917.0 to -91.3125*32=-2,922, where the larger of the intercepts (-2,917) is excluded as before. Thus the range choice for B is -2,922 through -2918 inclusive. Thus the alternative formula: (979*m-2918)\32 can be used to generate the required function. Here "\" is a shift which implements integer division with truncation rather than the greatest integer function. This is possible because variable "m" has a value of 3 or greater so that the numerator is always positive.