Index: trunk/SRC/Calendar/caldat.pro
===================================================================
--- trunk/SRC/Calendar/caldat.pro	(revision 133)
+++ trunk/SRC/Calendar/caldat.pro	(revision 133)
@@ -0,0 +1,198 @@
+;+
+;
+; @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
Index: trunk/SRC/Calendar/julday.pro
===================================================================
--- trunk/SRC/Calendar/julday.pro	(revision 133)
+++ trunk/SRC/Calendar/julday.pro	(revision 133)
@@ -0,0 +1,265 @@
+;+
+;
+; @file_comments
+;       Calculate the Julian Day Number for a given month, day, and year.
+;	This is the inverse of the library function CALDAT.
+;	See also caldat, the inverse of this function.
+;
+; @categories Calendar
+;
+; @param MONTH {in}{required} Number of the desired month (1 = January, ..., 12 = December). Can be scalar or array
+;
+; @param DAY {in}{required} Number of day of the month.Can be scalar or array
+;
+; @param YEAR {in}{required} Number of the desired year.Year parameters must be valid
+;               values from the civil calendar.  Years B.C.E. are represented
+;               as negative integers.  Years in the common era are represented
+;               as positive integers.  In particular, note that there is no
+;               year 0 in the civil calendar.  1 B.C.E. (-1) is followed by
+;               1 C.E. (1). Can be scalar or array
+;
+; @param HOUR {in}{optional} Number of the hour of the day. Can be scalar or array
+;
+; @param MINUTE {in}{optional} Number of the minute of the hour. Can be scalar or array
+;
+; @param SECOND {in}{optional} Number of the second of the minute. Can be scalar or array
+;
+; @restrictions The Result will have the same dimensions as the smallest array, or
+;         will be a scalar if all arguments are scalars.
+; 
+;
+; @keywords NDAYSPM {default=30} for using a calendar with fixed number of days per
+;                months.
+;
+; @ returns JULDAY: the Julian Day Number (which begins at noon) of the
+;	specified calendar date.  If Hour, Minute, and Second are not specified,
+;	then the result will be a long integer, otherwise the result is a
+;	double precision floating point number.
+;
+; @uses cm_4cal
+;
+; @restrictions Accuracy using IEEE double precision numbers is approximately
+;   1/10000th of a second, with higher accuracy for smaller (earlier)
+;   Julian dates.
+;
+; @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).
+;
+;	AB, September, 1988
+;	DMS, April, 1995, Added time of day.
+;
+;       Eric Guilyardi, June 1999
+;       Added key_work ndayspm for fixed number of days per months
+;       + change to accept year 0
+;
+;       Sebastien Masson, Aug. 2003
+;       fix bug for negative and large values of month values
+;       eg. julday(349,1,1970)
+;
+;   CT, April 2000, Now accepts vectors or scalars.
+;
+; @version $Id$
+;-
+;
+function JULDAY, MONTH, DAY, YEARin, 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'
+;
+
+  YEAR = long(yearin)
+  zero = where(year EQ 0, cnt)
+  IF cnt NE 0 THEN YEAR[zero] = 654321L
+;
+  CASE key_caltype OF
+    'greg':BEGIN
+
+
+; Gregorian Calander was adopted on Oct. 15, 1582
+; skipping from Oct. 4, 1582 to Oct. 15, 1582
+      GREG = 2299171L           ; incorrect Julian day for Oct. 25, 1582
+
+; Process the input, if all are missing, use todays date.
+      NP = n_params()
+      IF (np EQ 0) THEN RETURN, SYSTIME(/JULIAN)
+      IF (np LT 3) THEN MESSAGE, 'Incorrect number of arguments.'
+
+; Find the dimensions of the Result:
+;  1. Find all of the input arguments that are arrays (ignore scalars)
+;  2. Out of the arrays, find the smallest number of elements
+;  3. Find the dimensions of the smallest array
+
+; Step 1: find all array arguments
+      nDims = [SIZE(month, /N_DIMENSIONS), SIZE(day, /N_DIMENSIONS), $
+               SIZE(year, /N_DIMENSIONS), SIZE(hour, /N_DIMENSIONS), $
+               SIZE(minute, /N_DIMENSIONS), SIZE(second, /N_DIMENSIONS)]
+      arrays = WHERE(nDims GE 1)
+
+      nJulian = 1L              ; assume everything is a scalar
+      IF (arrays[0] GE 0) THEN BEGIN
+                                ; Step 2: find the smallest number of elements
+        nElement = [N_ELEMENTS(month), N_ELEMENTS(day), $
+                    N_ELEMENTS(year), N_ELEMENTS(hour), $
+                    N_ELEMENTS(minute), N_ELEMENTS(second)]
+        nJulian = MIN(nElement[arrays], whichVar)
+                                ; step 3: find dimensions of the smallest array
+        CASE arrays[whichVar] OF
+          0: julianDims = SIZE(month, /DIMENSIONS)
+          1: julianDims = SIZE(day, /DIMENSIONS)
+          2: julianDims = SIZE(year, /DIMENSIONS)
+          3: julianDims = SIZE(hour, /DIMENSIONS)
+          4: julianDims = SIZE(minute, /DIMENSIONS)
+          5: julianDims = SIZE(second, /DIMENSIONS)
+        ENDCASE
+      ENDIF
+
+      d_Second = 0d             ; defaults
+      d_Minute = 0d
+      d_Hour = 0d
+; convert all Arguments to appropriate array size & type
+      SWITCH np OF              ; use switch so we fall thru all arguments...
+        6: d_Second = (SIZE(second, /N_DIMENSIONS) GT 0) ? $
+                      second[0:nJulian-1] : second
+        5: d_Minute = (SIZE(minute, /N_DIMENSIONS) GT 0) ? $
+                      minute[0:nJulian-1] : minute
+        4: d_Hour = (SIZE(hour, /N_DIMENSIONS) GT 0) ? $
+                    hour[0:nJulian-1] : hour
+        3: BEGIN                ; convert m,d,y to type LONG
+          L_MONTH = (SIZE(month, /N_DIMENSIONS) GT 0) ? $
+                    LONG(month[0:nJulian-1]) : LONG(month)
+          L_DAY = (SIZE(day, /N_DIMENSIONS) GT 0) ? $
+                  LONG(day[0:nJulian-1]) : LONG(day)
+          L_YEAR = (SIZE(year, /N_DIMENSIONS) GT 0) ? $
+                   LONG(year[0:nJulian-1]) : LONG(year)
+        END
+      ENDSWITCH
+
+
+      min_calendar = -4716
+      max_calendar = 5000000
+      minn = MIN(l_year, MAX = maxx)
+      IF (minn LT min_calendar) OR (maxx GT max_calendar) THEN MESSAGE, $
+        'Value of Julian date is out of allowed range.'
+; change to accept year 0
+; if (MAX(L_YEAR eq 0) NE 0) then message, $
+; 	'There is no year zero in the civil calendar.'
+;
+; by seb Aug 2003
+
+      tochange = where(L_MONTH LT 0)
+      IF tochange[0] NE -1 THEN BEGIN
+        L_YEAR[tochange] = L_YEAR[tochange]+L_MONTH[tochange]/12-1
+        L_MONTH[tochange] =  12 + L_MONTH[tochange] MOD 12
+      ENDIF
+
+      tochange = where(L_MONTH GT 12)
+      IF tochange[0] NE -1 THEN BEGIN
+        L_YEAR[tochange] = L_YEAR[tochange]+L_MONTH[tochange]/12
+        L_MONTH[tochange] =  L_MONTH[tochange] MOD 12
+      ENDIF
+; by seb Aug 2003 - end
+;
+;
+      bc = (L_YEAR LT 0)
+      L_YEAR = TEMPORARY(L_YEAR) + TEMPORARY(bc)
+      inJanFeb = (L_MONTH LE 2)
+      JY = L_YEAR - inJanFeb
+      JM = L_MONTH + (1b + 12b*TEMPORARY(inJanFeb))
+
+      JUL = floor(365.25d * JY) + floor(30.6001d*TEMPORARY(JM)) + L_DAY + 1720995L
+
+; Test whether to change to Gregorian Calandar.
+      IF (MIN(JUL) GE GREG) THEN BEGIN ; change all dates
+        JA = long(0.01d * TEMPORARY(JY))
+        JUL = TEMPORARY(JUL) + 2L - JA + long(0.25d * JA)
+      ENDIF ELSE BEGIN
+        gregChange = WHERE(JUL ge GREG, ngreg)
+        IF (ngreg GT 0) THEN BEGIN
+          JA = long(0.01d * JY[gregChange])
+          JUL[gregChange] = JUL[gregChange] + 2L - JA + long(0.25d * JA)
+        ENDIF
+      ENDELSE
+
+
+; hour,minute,second?
+      IF (np GT 3) THEN BEGIN   ; yes, compute the fractional Julian date
+; Add a small offset so we get the hours, minutes, & seconds back correctly
+; if we convert the Julian dates back. This offset is proportional to the
+; Julian date, so small dates (a long, long time ago) will be "more" accurate.
+        eps = (MACHAR(/DOUBLE)).eps
+        eps = eps*ABS(jul) > eps
+; For Hours, divide by 24, then subtract 0.5, in case we have unsigned ints.
+        jul = TEMPORARY(JUL) + ( (TEMPORARY(d_Hour)/24d - 0.5d) + $
+                                 TEMPORARY(d_Minute)/1440d + TEMPORARY(d_Second)/86400d + eps )
+      ENDIF
+
+; check to see if we need to reform vector to array of correct dimensions
+      IF (N_ELEMENTS(julianDims) GT 1) THEN $
+        JUL = REFORM(TEMPORARY(JUL), julianDims)
+
+      RETURN, jul
+
+    END 
+    '360d':BEGIN
+;
+; Fixed number of days per month (default=30) :
+;
+      IF keyword_set(ndayspm) THEN BEGIN
+        IF ndayspm EQ 1 THEN ndayspm = 30
+      ENDIF ELSE ndayspm = 30
+
+      L_MONTH = LONG(MONTH)
+      L_DAY = LONG(DAY)
+      L_YEAR = LONG(YEAR)
+
+      neg = where(L_YEAR LT 0)
+      IF neg[0] NE -1 THEN L_YEAR[neg] =  L_YEAR[neg]+1
+
+      JUL = ((L_YEAR-1)*12 + (L_MONTH-1))* ndayspm + L_DAY 
+      if n_elements(Hour) + n_elements(Minute) + n_elements(Second) eq 0 then $
+        return, JUL
+      if n_elements(Hour) eq 0 then Hour = 0
+      if n_elements(Minute) eq 0 then Minute = 0
+      if n_elements(Second) eq 0 then Second = 0
+      
+      IF Hour+Minute+Second EQ 0 THEN return, JUL ELSE $
+        return, JUL + (Hour / 24.0d0) + (Minute/1440.0d0) + (Second / 86400.0d0)
+
+    END 
+    'noleap':BEGIN
+
+      L_MONTH = LONG(MONTH)
+      L_DAY = LONG(DAY)
+      L_YEAR = LONG(YEAR)
+;
+      tochange = where(L_MONTH LT 0)
+      IF tochange[0] NE -1 THEN BEGIN
+        L_YEAR[tochange] = L_YEAR[tochange]+L_MONTH[tochange]/12-1
+        L_MONTH[tochange] =  12 + L_MONTH[tochange] MOD 12
+      ENDIF
+;
+      tochange = where(L_MONTH GT 12)
+      IF tochange[0] NE -1 THEN BEGIN
+        L_YEAR[tochange] = L_YEAR[tochange]+L_MONTH[tochange]/12
+        L_MONTH[tochange] =  L_MONTH[tochange] MOD 12
+      ENDIF
+;
+      L_YEAR =  L_YEAR - 1
+;
+      daysyear = long(total([0, 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30], /cumulative))
+
+      return, 365*L_YEAR + daysyear[L_MONTH] + L_DAY
+
+    END 
+    ELSE:return, report('only 3 types of calendar are accepted: greg, 360d and noleap')
+  ENDCASE
+
+END
