ymd2jd <- function(yr=2000, mo=1, dy=1) {

# Year, month, date to Julian day conversion
# Rounds to nearest Julian date at input
# Reference time is 12:00; to get to 00:00 UT on the date, subtract 0.5
#
# 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, 4/18/2008

dy <- round(dy)
mo <- round(mo)
yr <- round(yr)

jd <- dy - 32075 + trunc(1461*(yr + 4800 + trunc((mo - 14)/12))/4) +
          trunc(367*(mo - 2 - 12*trunc((mo-14)/12))/12) -   
          trunc(3*(trunc((yr + 4900 + trunc((mo - 14)/12))/100))/4)

return(jd)

}

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

ut2dmjd <- function(yr=2000, mo=1, dy=1, hr=0, mi=0, se=0, eps=1e-9) {
# Convert from UT day and time to decimal modified Julian date
#
# A. Harris, U. Maryland Astronomy, 4/18/2008
# 
# Uses function ymd2jd to compute Julian date
# eps solves roundoff problem for exactly 0h, 1e-9 adds 1e-4 sec

mjd <- ymd2jd(yr, mo, dy) + (hr + (mi + se/60)/60)/24 - 2400001 + eps
  
return(mjd)
 
}