source: trunk/SRC/Calendar/julday.pro @ 378

;+
;
; @file_comments
; Calculate the Julian Day Number for a given month, day, and year.
; This is the inverse of the library function <pro>caldat</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 month {in}{optional} {type=scalar (integer or double) or array of scalars}
; Number of the desired month (1 = January, ..., 12 = December).
;
; @param day {in}{optional} {type=scalar (integer or double) or array of scalars}
; Number of day of the month.
;
; @param yearin {in}{optional} {type=scalar (integer or double) or array of scalars}
; 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).
; Change: However for climatological year, we do accept the year
; 0 but we change it for year 654321L (the same trick is done in
; <pro>caldat</pro> so caldat, julday(1,1,0) gives you back Jan 1st of year 0)
;
; @param hour {in}{optional} {type=scalar (integer or double) or array of scalars} {default=12}
; Number of the hour of the day.
;
; @param minute {in}{optional} {type=scalar (integer or double) or array of scalars} {default=0}
; Number of the minute of the hour.
;
; @param second {in}{optional} {type=scalar (integer or double) or array of scalars} {default=0}
; Number of the second of the minute.
;
; @restrictions
; The result will have the same dimensions as the smallest array, or
; will be a scalar if all arguments are scalars.
;
; @keyword NDAYSPM {default=30} {type=integer}
; 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)
;
; @returns
; 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
; <pro>cm_4cal</pro>
;
; @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 Recipes 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
;
; CT, April 2000, Now accepts vectors or scalars.
;
; Sebastien Masson, Aug. 2003
; fix bug for negative and large values of month values
; eg. julday(349,1,1970)
;
; Sebastien Masson, May 2006, add different calendat with key_caltype
; (variable of the common file cm_4cal)
;
; @version
; $Id$
;-
FUNCTION julday, month, day, yearin, 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'
  
  NP = n_params()
; Process the input, if all are missing, use today's date.
  IF (np EQ 0) THEN RETURN, SYSTIME(/JULIAN)
  IF (np LT 3) THEN ras = report('Incorrect number of arguments.')
;
  YEAR = long(yearin)
  zero = where(year EQ 0, cnt)
  IF cnt NE 0 THEN YEAR[zero] = 654321L
;
  CASE key_caltype OF
    'greg':BEGIN

; Gregorian Calender was adopted on Oct. 15, 1582
; skipping from Oct. 4, 1582 to Oct. 15, 1582
      GREG = 2299171L           ; incorrect Julian day for Oct. 25, 1582


; 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 $
        ras = report('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 Calendar.
      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 integers.
        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 = 12
      if n_elements(Minute) eq 0 then Minute = 0
      if n_elements(Second) eq 0 then Second = 0

      IF (np LE 3) THEN return, JUL ELSE $
        return, JUL + (Hour / 24.0d0 - 0.5d) + (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))

      JUL = 365*L_YEAR + daysyear[L_MONTH] + 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 = 12
      if n_elements(Minute) eq 0 then Minute = 0
      if n_elements(Second) eq 0 then Second = 0

        IF (np LE 3) THEN return, JUL ELSE $
        return, JUL + (Hour / 24.0d0 - 0.5d) + (Minute/1440.0d0) + (Second / 86400.0d0)

    END
    ELSE:return, report('only 3 types of calendar are accepted: greg, 360d and noleap')
  ENDCASE

END