jd2ymd <- function(jd) {

# Julian day to year, month, date conversion
# Rounds to nearest Julian date at input
#
# From Fliegel & Van Flandern, Comm. ACM 10, 657 (1968)
#   Note: algorithm uses FORTRAN integer mathematics
#   Also: Explanatory Supplement to the Astronomical Almanac
# A. Harris, U. Maryland Astronomy, 3/17/08

L <- round(jd) + 68569
N <- trunc(4 * L/146097)
L <- L - trunc((146097*N + 3)/4)
I <- trunc(4000*(L+1)/1461001)
L <- L - trunc(1461*I/4) + 31
J <- trunc(80*L/2447)
K <- L - trunc(2447*J/80)
L <- trunc(J/11)
J <- J + 2 - 12*L
I <- 100*(N-49) + I + L

list(y=I, m=J, d=K)
}

###############################

dmjd2ut <- function(mjd, tz=0) {

# Decimal Modified Julian Date to UT calculation
# tz is an optional time offset from UT; tz=-4 for EST
# A. Harris, U. Maryland Astronomy, 3/17/08

# Time offset
mjd <- mjd + tz/24

# Convert to Julian date and compute year, month, and day
jd <- mjd + 2400000.5
ymd <- jd2ymd(jd)

# Work out hours, minutes and seconds
dayfrac <- mjd%%1
hr <- dayfrac*24
min <- (hr%%1)*60
sec <- (min%%1)*60

list(yr=ymd$y, mo=ymd$m, dy=ymd$d, hr=trunc(hr), 
     mi=trunc(min), se=sec) 
}