sidereal.py: Internal maintenance specification

Converts a JulianDate instance to a Python datetime value. The algorithm here is rather inscrutable; refer to Duffett-Smith for the details.

# - - -   J u l i a n D a t e . d a t e t i m e

    def datetime ( self ):
        """Convert to a standard Python datetime object in UT.
        """

Step 1 of the algorithm reads: “Add 0.5 to JD. Set I=integer part and F=fractional part.” Once we have extracted the fractional days part, we can safely add the JULIAN_BIAS back in with no loss of significance.

        #-- 1 --
        # [ i  :=  int(self.j + 0.5)
        #   f  :=  (self.j + 0.5) % 1.0 ]
        i, f  =  divmod ( self.j + 0.5, 1.0 )
        i  +=  JULIAN_BIAS

Step 2: “If I is larger than 2,299,160, calculate A = integer part of ((I- 1,867,216.25) / 36,524.25); calculate B = I + 1 + a - integer part of (A/4); otherwise set A = I.” The last part of this sentence should probably set “...otherwise set B = I,” since A is not referenced after this step, but B certainly is.

        #-- 2 --
        if  i > 2299160:
            a  =  int((i-1867216.25)/36524.25)
            b  =  i + 1 + a - int ( a / 4.0 )
        else:
            b  =  i

Step 3: “Calculate C = B + 1524.”

        #-- 3 --
        c = b + 1524

Step 4: “Calculate D = integer part of (C-122.1)/365.25.”

        #-- 4 --
        d = int((c-122.1)/365.25)

Step 5: “Calculate E = integer part of (365.25 × D ).”

        #-- 5 --
        e = int(365.25*d)

Step 6: “Calculate G = integer part of ((C-E)/30.6001).”

        #-- 6 --
        g = int((c-e)/30.6001)

Step 7: “Calculate d=C-E+F-integer part of (30.6001 × G). This is the day of the month, including the decimal fraction of the day.” We will remove the fractional day and convert it to hours, minutes, and seconds.

        #-- 7 --
        dayFrac = c - e + f - int ( 30.6001 * g )
        day, frac = divmod ( dayFrac, 1.0 )
        dd = int(day)
        hr, mn, sc = dmsUnits.singleToMix ( 24.0*frac )

Step 8: “Calculate m=G-1 if G is less than 13.5, or m=G-13 if G is more than 13.5. This is the month number.”

        #-- 8 --
        if  g < 13.5:  mm = int(g - 1)
        else:             mm = int(g - 13)

Step 9: “Calculate y = D-4716 if m is more than 2.5, or y = D - 4715 if m is less than 2.5. This is the year.”

        #-- 9 --
        if  mm > 2.5:  yyyy = int(d-4716)
        else:          yyyy = int(d-4715)

Now that we have all the pieces, assemble them into a datetime and return it. One further elaboration: the value of sc is a float, and the datetime.datetime() constructor wants fractional seconds passed as a “microseconds” argument.

        #-- 10 --
        sec, fracSec = divmod(sc, 1.0)
        usec = int(fracSec * 1e6)
        return datetime.datetime ( yyyy, mm, dd, hr, mn, int(sec),
                                   usec )