def JD(year, month, day):
    """
----------------------------------------------------------------------
Purpose:   Calculate Julian day number (Julian date)
Input:     year (nnnn), month (nn), day (nn.n),
Returns:   Julian day number 'jd'
Reference: Meeus, Astronomical formula for Calculators, 2nd ed, 1982
Notes:     Months start at 1. Days start at 1. The Julian day begins at
           Greenwich mean noon, i.e. at 12h. So Jan 1, 1984 at 0h is 
           entered as JD(1984,1,1) and Jan 1, 1984 at 12h is entered 
           as JD(1984,1,1.5)
           There is a jump at JD(1582,10,15) caused by the change of
           calendars. So for dates after 1582-10-15 one enters a date
           from the Julian calendar and before this date you enter a 
           date from the Julian calendar. 
Examples:  Julian date of JD reference:
           print celestial.JD(-4712,1,1.5) ==> 0.0 
           The first day of 1 B.C.:
           print celestial.JD(0,1,1) ==> 1721057.5
           Last day before Gregorian reform:
           print celestial.JD(1582,10,4) ==> 2299159.5
           First day of Gregorian reform:
           print celestial.JD(1582,10,15) ==> 2299170.5
           Half a day later:
           print celestial.JD(1582,10,15.5) ==> 2299161.0
           Unix reference:
           print celestial.JD(1970,1,1) ==> 2440587.5
----------------------------------------------------------------------
   """
    if month > 2:
        y = year
        m = month
    elif month == 1 or month == 2:
        y = year - 1
        m = month + 12
    calday = year + month / 100.0 + day / 10000.0
    if calday > 1582.1015:
        A = int(y / 100.0)
        B = 2 - A + int(A / 4.0)
    else:
        B = 0
    if calday > 0.0229:
        jd = int(365.25 * y) + int(30.6001 * (m + 1)) + day + 1720994.5 + B
    else:
        jd = int(365.25 * y - 0.75) + int(30.6001 * (m + 1)
            ) + day + 1720994.5 + B
    return jd
