source: trunk/SRC/Calendar/caldat.pro @ 134

;+
;
; @file_comments
;       Return the calendar date and time given julian date.
;       This is the inverse of the function JULDAY.
; 
; @categories Calendar
;
;
; @param JULIAN {in}{required} contains the Julian Day Number (which begins at noon) of the
;       specified calendar date.  It should be a long integer.
;
; @param MONTH {out} Number of the desired month (1 = January, ..., 12 = December).
;
; @param DAY {out} Number of day of the month.
;
; @param YEAR {out} Number of the desired year.
;
; @param HOUR {out} Hour of the day
;
; @param Minute {out} Minute of the day
;
; @param Second {out} Second (and fractions) of the day.
;
;
; @keyword NDAYSPM {default=30} for using a calendar with fixed number of days per
;                months.
;
; @uses cm_4cal
;
;
; @restrictions Accuracy using IEEE double precision numbers is approximately
;       1/10000th of a second.
;
; @history Translated from "Numerical Recipies 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.
;
; @version $Id$
;-
pro CALDAT, julian, month, day, year, hour, minute, second, NDAYSPM = ndayspm
;------------------------------------------------------------
@cm_4cal
;------------------------------------------------------------
  COMPILE_OPT idl2

  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'
  CASE key_caltype OF
    'greg':BEGIN

      nParam = N_PARAMS()
      IF (nParam LT 1) THEN MESSAGE, 'Incorrect number of arguments.'

      min_julian = -1095
      max_julian = 1827933925
      minn = MIN(julian, MAX = maxx)
      IF (minn LT min_julian) OR (maxx GT max_julian) THEN MESSAGE, $
        '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)
; see if we need to do hours, minutes, seconds
      IF (nParam GT 4) THEN BEGIN
        fraction = julian + 0.5d - TEMPORARY(julLong)
        hour = floor(fraction * 24d)
        fraction = TEMPORARY(fraction) - hour/24d
        minute = floor(fraction*1440d)
        second = (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)
        day = REFORM(day, dimensions)
        year = REFORM(year, dimensions)
        IF (nParam GT 4) THEN BEGIN
          hour = REFORM(hour, dimensions)
          minute = REFORM(minute, dimensions)
          second = REFORM(second, dimensions)
        ENDIF
      ENDIF

    END 
    '360d':BEGIN

      jul = long(julian)
      f = (jul - floor(jul))
      IF total(f NE 0.0) GT 0 THEN BEGIN ;Get hours, minutes, seconds.
        hour = floor(f*24.)
        f = f - hour / 24.d
        minute = floor(f*1440)
        second = (f - minute/1440.d0) * 86400.0d0
      ENDIF ELSE BEGIN
        hour = replicate(0L,  n_elements(julian))
        minute = replicate(0L,  n_elements(julian))
        second = replicate(0L,  n_elements(julian))
      ENDELSE

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

      ndayspm = long(ndayspm)
      Z = floor(julian)
      year = z/(12*ndayspm)+1
      month = (z-(12*ndayspm)*(year-1))/ndayspm+1
      day = z-(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

      jul = long(julian)
      year = jul/365 + 1
      day = jul 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
;
  return

END