;; ************************************************************************
;;   -DOC
;;   -NAME:  Mtim_ymd_to_jday 
;;   (\p{jd}, \p{yr}, \p{mo}, \p{da}, \p{offset}, \p{mjd}, \p{tjd}, \p{invert})
;;
;;   -DESCRIPTION:
;;   
;;   Convert dates in the year/month/day notation to and from  Julian Day
;;   notation.
;;
;;   -PARAMETERS:
;;   \p{julday}        & points to returned array, which will contain
;;                        the Julian Day equivalents. \\
;;   \p{yr}        & points to array of year numbers.  Valid years 
;;                        must be greater than -4700 (i.e., BCE).\\
;;   \p{mo}        & points to array of month numbers, ranging between
;;                        one (January) and twelve (December). \\
;;   \p{da}     & points to array of day numbers, starting with 1.0 at
;;                midnight on the beginning of the month (not at
;;                noon!).  Example: the midnight between 31 Dec 1997
;;                and 01 Jan 1998 would be day number 1.0 of Jan.  ;;
;;                Days may be fractional in order to indicate
;;                partial days.\\
;;   \p{offset}  & numerical offset, in days, to subtract from
;;                 returned day values.  To return Julian days, use an
;;                 offset of ;; 0.0D.  If \p{mjd} or \p{tjd} are set,
;;                 then \p{offset} is set to 2400000.5 or 2440000.5,
;;                 respectively. With the proper choice of \p{offset},
;;                 better than nanosecond precision can be maintained.
;;   \p{invert}  & if set, then the calendar date is found from the
;;                 passed values of \p{jd}.\\
;;
;;   -RETURN:
;;   
;;   Return value is the array of julian days.
;;
;;   -REFERTO:
;;
;;   -INDEX: date conversion, time conversion, Julian Day, 
;;           conversion to Julian Day
;;
;;   -AUTHOR:   Craig Markwardt
;;   -HISTORY:
;;   \Mitem[08 Sep 1997] - Written, cm
;;   \Mitem[09 Feb 1998] - Added offset parameter, convert to IDL
;;   \Mitem[14 May 2002] - Added JDAY_XMOD to make standalone
;;
;;   -REFERENCE:
;;   \Mcite{Hatcher, D.A. 1984, {\it Royal
;;   Astronomical Society Quarterly Journal}, 25, 53}
;;   -IMPLEMENTATION_NOTES:
;;   -ENDDOC
;;

forward_function jday_xmod
function jday_xmod, x, m
  return, (((x MOD m) + m) MOD m)
end

function jday, yr, mo, da, jday=jd, $
               mjd=mjd, tjd=tjd, offset=offset, invert=invert, $
               juliancal=juliancal, now=now

  if keyword_set(mjd) then offset = 2400000.5D
  if keyword_set(tjd) then offset = 2440000.5D
  if n_elements(offset) EQ 0 then offset=0.D
  if NOT keyword_set(now) then begin
      if NOT keyword_set(invert) AND n_params() EQ 0 then begin
          print, 'USAGE:'
          print, '  JD  = JDAY(YEAR, MONTH, DAY, [/TJD, /MJD]) - convert YMD to JD'
          print, '  YMD = JDAY(JDAY=JD, /INVERT, [/TJD, /MJD]) - convert JD to YMD'
          print, '  JD  = JDAY(/NOW, [/INVERT, /TJD, /MJD])    - use the present time'
          return, -1D
      endif
      if keyword_set(invert) AND n_elements(jd) EQ 0 then begin
          message, 'ERROR: You must pass JDAY when /INVERTing.'
          return, -1D
      endif
      if NOT keyword_set(invert) AND $
        (n_elements(yr) EQ 0 OR n_elements(mo) EQ 0 $
         OR n_elements(da) EQ 0) then begin
          message, 'ERROR: You must pass YEAR, MONTH and DAY.'
          return, -1D
      endif
  endif

  offset_int = floor(offset)         ;; Integer part of offset
  offset_fra = offset - offset_int   ;; Fractional part of offset

  ;; WARNING: does not account for leap seconds
  if keyword_set(now) then begin
      j0 = 2440587.5D
      jd = (j0-offset) + systime(1)/86400D
  endif

  ;; This is the conversion TO Julian Days
  if NOT keyword_set(invert) then begin

      if keyword_set(now) then return, jd

      d0 = 365.25D * 4712.D - offset_int

      ;; Convert to "offset years" (see Hatcher), which start on 01 Mar
      yr2 = yr - floor((12.-mo)/10)
      mo2 = jday_xmod(mo - 3.D, 12.D)

      ;; Julian day conversion
      yy = floor(d0 + 365.25D*yr2)
      dd = floor(30.6D * mo2 + 0.5D)
      jd = yy + dd + 59.D + da - 0.5D
      
      ;; Correction if dates are in Gregorian calendar
      if NOT keyword_set(juliancal) then $
        jd = jd - (floor(floor((yr2/100.D) + 49.D) * 0.75D) - 38.D)

      jd = jd - offset_fra

      return, jd
  endif else begin
      ;; This is the conversion FROM Julian days to calendar dates
      nn = offset_fra + jd
      ;; Strangeness here is to account for beginnings and endings of months
      jd_fra = jday_xmod(nn+0.5D, 1D) - 0.5D
      nn = nn + offset_int - jd_fra

      if NOT keyword_set(juliancal) then $
        nn = nn + (floor(floor((nn - 4479.5D)/36524.25D) * 0.75D + 0.5D)-37.D)

      yr = long(   floor(nn/365.25D) - 4712.D   )
      dd = floor(  jday_xmod(nn-59.25D, 365.25D)  )

      mo = floor(  jday_xmod( floor((dd+0.5D)/30.6D) + 2.D, 12.D ) + 1.D  )
      da = floor( jday_xmod(dd+0.5D, 30.6D) + 1.D ) + 0.5D + jd_fra

      if n_params() EQ 0 then begin
          zz1 = transpose([yr])
          zz2 = transpose([mo])
          zz3 = transpose([da])
          return, [ double(zz1), double(zz2), double(zz3) ]
      endif else return, [double(yr(0)), mo(0), da(0)]
  endelse

  return, jd
end