source: trunk/SRC/Calendar/caldat.pro

;+
;
; @file_comments
; Return the calendar date and time given Julian date.
; This is the inverse of the function <pro>julday</pro>.
; 3 calendars are available according to the value of key_caltype
; (variable of the common file cm_4cal): 'greg', '360d', 'noleap'
;
; @categories
; Calendar
;
; @param julian {in}{required} {type=long integer or double-precision floating-point}
; contains the Julian Day Number (which begins at noon) of the
; specified calendar date.
;
; @param month {out} {type=longword integer}
; Number of the desired month (1 = January, ..., 12 = December).
;
; @param day {out} {type=longword integer}
; Number of day of the month.
;
; @param year {out} {type=longword integer}
; Number of the desired year.
;
; @param hour {out} {type=longword integer}
; Hour of the day
;
; @param minute {out} {type=longword integer}
; Minute of the day
;
; @param second {out} {type=double-precision floating-point}
; Second (and fractions) of the day.
;
; @keyword NDAYSPM {type=integer} {default=30}
; To use a calendar with fixed number of days per month.
; see also the use of key_caltype (variable of the common file cm_4cal)
;
; @uses
; <pro>cm_4cal</pro>
;
; @restrictions
; Accuracy using IEEE double precision numbers is approximately 1/10000th of a
; second.
;
; @history
; Translated from "Numerical Recipes in C", by William H. Press,
; Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling.
; Cambridge University Press, 1988 (second printing).
;
; DMS, July 1992.
; DMS, April 1996, Added HOUR, MINUTE and SECOND keyword
; AB, 7 December 1997, Generalized to handle array input.
;
; Eric Guilyardi, June 1999
; Added key_work ndayspm for fixed number of days per months
;
; AB, 3 January 2000, Make seconds output as DOUBLE in array output.
;
; Sebastien Masson, May 2006, add different calendar with key_caltype
; (variable of the common file cm_4cal)
;
;   CT, Nov 2006: For Hour/Min/Sec, tweak the input to make sure hours
;       and minutes are correct. Restrict hours to 0-23 & min to 0-59.
;
; @version
; $Id$
;
;-
PRO caldat, julian, month, day, year, hour, minute, second, NDAYSPM=ndayspm
;
  compile_opt idl2, strictarrsubs
;
@cm_4cal
;
  ON_ERROR, 2                   ; Return to caller if errors

  IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
  if keyword_set(ndayspm) then key_caltype = '360d'

  nParam = N_PARAMS()
  IF (nParam LT 1) THEN ras = report('Incorrect number of arguments.')

  CASE key_caltype OF
    'greg':BEGIN

      min_julian = -1095
      max_julian = 1827933925
      minn = MIN(julian, MAX = maxx)
      IF (minn LT min_julian) OR (maxx GT max_julian) THEN $
        ras = report('Value of Julian date is out of allowed range.')

      igreg = 2299161L                   ;Beginning of Gregorian calendar
      julLong = FLOOR(julian + 0.5d)     ;Better be long
      minJul = MIN(julLong)

      IF (minJul GE igreg) THEN BEGIN ; all are Gregorian
        jalpha = LONG(((julLong - 1867216L) - 0.25d) / 36524.25d)
        ja = julLong + 1L + jalpha - long(0.25d * jalpha)
      ENDIF ELSE BEGIN
        ja = julLong
        gregChange = WHERE(julLong ge igreg, ngreg)
        IF (ngreg GT 0) THEN BEGIN
          jalpha = long(((julLong[gregChange] - 1867216L) - 0.25d) / 36524.25d)
          ja[gregChange] = julLong[gregChange] + 1L + jalpha - long(0.25d * jalpha)
        ENDIF
      ENDELSE
      jalpha = -1               ; clear memory

      jb = TEMPORARY(ja) + 1524L
      jc = long(6680d + ((jb-2439870L)-122.1d0)/365.25d)
      jd = long(365d * jc + (0.25d * jc))
      je = long((jb - jd) / 30.6001d)

      day = TEMPORARY(jb) - TEMPORARY(jd) - long(30.6001d * je)
      month = TEMPORARY(je) - 1L
      month = ((TEMPORARY(month) - 1L) MOD 12L) + 1L
      year = TEMPORARY(jc) - 4715L
      year = TEMPORARY(year) - (month GT 2)
      year = year - (year LE 0)

    END
    '360d':BEGIN

      IF keyword_set(ndayspm) THEN BEGIN
        IF ndayspm EQ 1 THEN ndayspm = 30
      ENDIF ELSE ndayspm = 30

      ndayspm = long(ndayspm)
      julLong = FLOOR(julian + 0.5d)     ;Better be long
      year = julLong/(12*ndayspm)+1
      month = (julLong-(12*ndayspm)*(year-1))/ndayspm+1
      day = julLong-(12*ndayspm)*(year-1)-ndayspm*(month-1)
      WHILE total(day LT 1) GT 0 DO BEGIN
        tochange = where(day LT 1)
        month[tochange] = month[tochange]-1
        day[tochange] = day[tochange]+ndayspm
      ENDWHILE
      WHILE total(month LT 1) GT 0 DO BEGIN
        tochange = where(month LT 1)
        year[tochange] = year[tochange]-1
        month[tochange] = month[tochange]+12
      ENDWHILE
; year 0 does not exist...
      neg = where(year LT 0)
      IF neg[0] NE -1 THEN year[neg] = year[neg]-1
    END
    'noleap':BEGIN

      julLong = FLOOR(julian + 0.5d)     ;Better be long
      year = julLong/365 + 1
      day = julLong MOD 365L
;
      zero = where(day EQ 0)
;
      month = 1 + (day GT 31) + (day GT 59) + (day GT 90) + (day GT 120) $
              + (day GT 151) + (day GT 181) + (day GT 212) + (day GT 243) $
              + (day GT 273) + (day GT 304) + (day GT 334)
      month = long(month)
;
      day = day - 31L * (day GT 31) - 28L * (day GT 59) - 31L * (day GT 90) $
              - 30L * (day GT 120) - 31L * (day GT 151) - 30L * (day GT 181) $
              - 31L * (day GT 212) - 31L * (day GT 243) - 30L * (day GT 273) $
              - 31L * (day GT 304) - 30L * (day GT 334)
;
      IF zero[0] NE -1 THEN BEGIN
        year[zero] = year[zero]-1
        month[zero] = 12L
        day[zero] = 31L
      ENDIF
;
    END
    ELSE:BEGIN
      ng = report('only 3 types of calendar are accepted: greg, 360d and noleap')
      return
    END
  ENDCASE
;
  zero = where(year ge 600000L, cnt)
  IF cnt NE 0 THEN year[zero] = year[zero]-654321L
;
; see if we need to do hours, minutes, seconds
        IF (nParam GT 4) THEN BEGIN
                fraction = julian + 0.5d - julLong
        eps = 1d-12 > 1d-12*ABS(Temporary(julLong))
                hour = 0 > Floor(fraction * 24d + eps) < 23
                fraction -= hour/24d
                minute = 0 > Floor(fraction*1440d + eps) < 59
                second = 0 > (TEMPORARY(fraction) - minute/1440d)*86400d
        ENDIF

; if julian is an array, reform all output to correct dimensions
  IF (SIZE(julian, /N_DIMENSION) GT 0) THEN BEGIN
    dimensions = SIZE(julian, /DIMENSION)
    month = REFORM(month, dimensions, /overwrite)
    day = REFORM(day, dimensions, /overwrite)
    year = REFORM(year, dimensions, /overwrite)
    IF (nParam GT 4) THEN BEGIN
      hour = REFORM(hour, dimensions, /overwrite)
      minute = REFORM(minute, dimensions, /overwrite)
      second = REFORM(second, dimensions, /overwrite)
    ENDIF
  ENDIF
;
  return

END