;+
;
; @file_comments
; Calculate the Julian Day Number for a given month, day, and year.
; This is the inverse of the library function
; caldat.
; 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}{required} {type=scalar (integer or double) or array of scalars}
; Number of the desired month (1 = January, ..., 12 = December).
;
; @param DAY {in}{required} {type=scalar (integer or double) or array of scalars}
; Number of day of the month.
;
; @param YEARin {in}{required} {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
; caldat 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
; 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
;
; 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'
;
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
; 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 ras= report('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 $
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 total([hour NE 12, minute NE 0, second NE 0]) EQ 0 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 total([hour NE 12, minute NE 0, second NE 0]) EQ 0 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