Changeset 9 for trunk/ToBeReviewed/CALENDRIER/julday.pro

Timestamp:

04/26/06 16:29:38 (18 years ago)

Author:

pinsard

Message:

upgrade of CALENDRIER/Calendar according to cerbere.lodyc.jussieu.fr:/usr/home/smasson/SAXO_RD/ : files

File:

• trunk/ToBeReviewed/CALENDRIER/julday.pro (copied) (copied from trunk/CALENDRIER/julday.pro) (7 diffs)

trunk/ToBeReviewed/CALENDRIER/julday.pro

@@ -1 +1 @@
 ; $Id$
 ;
-; Copyright (c) 1988-1998, Research Systems, Inc.  All rights reserved.
+; Copyright (c) 1988-2003, Research Systems, Inc.  All rights reserved.
 ;       Unauthorized reproduction prohibited.
 
-function JULDAY, MONTH, DAY, YEAR, Hour, Minute, Second, NDAYSPM = ndayspm, _extra=ex
 ;+
 ; NAME:
@@ -13 +12 @@
 ;       This is the inverse of the library function CALDAT.
 ;       See also caldat, the inverse of this function.
+;
 ; CATEGORY:
 ;       Misc.
 ;
 ; CALLING SEQUENCE:
-;       Result = JULDAY(Month, Day, Year)
+;       Result = JULDAY([[[[Month, Day, Year], Hour], Minute], Second])
 ;
 ; INPUTS:
@@ -24 +24 @@
 ;       DAY:    Number of day of the month.
 ;
-;       YEAR:   Number of the desired year.
+;       YEAR:   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).
 ;
 ;       HOUR:   Number of the hour of the day.
@@ -30 +35 @@
 ;       MINUTE: Number of the minute of the hour.
 ;
-;       SECOND: 
+;       SECOND: Number of the second of the minute.
+;
+;   Note: Month, Day, Year, Hour, Minute, and Second can all be arrays.
+;         The Result will have the same dimensions as the smallest array, or
+;         will be a scalar if all arguments are scalars.
 ;
 ; OPTIONAL INPUT PARAMETERS:
@@ -37 +46 @@
 ; KEYWORD PARAMETERS:
 ;
-;       NDAYSPM: developpe par eric, ca veut dire: nombre de jours par mois! 
-;                par defaut c''est 30, sinon le specifier en donnant
-;                une valeur a ndayspm
-;                pour passer a un calendrier avec un nombre de jours constant par
-;                mois. Utilise en particulier ds julday et caldat
+;       NDAYSPM: for using a calendar with fixed number of days per
+;                months. defaut value of NDAYSPM=30
 ;
 ; OUTPUTS:
-;       JULDAY returns the Julian Day Number (which begins at noon) of the 
-;       specified calendar date.  If the time of day (Hr, Min, Sec), is 0,
-;       the result will be a long integer, otherwise the result is a 
+;       JULDAY 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.
 ;
-; COMMON BLOCKS:
-;       None.
+; COMMON BLOCKS: cm_4cal
 ;
 ; SIDE EFFECTS:
@@ -57 +62 @@
 ; RESTRICTIONS:
 ;       Accuracy using IEEE double precision numbers is approximately
-;       1/10000th of a second.
+;   1/10000th of a second, with higher accuracy for smaller (earlier)
+;   Julian dates.
 ;
 ; MODIFICATION HISTORY:
@@ -69 +75 @@
 ;       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.
 ;-
 ;
-   ON_ERROR, 2                  ; Return to caller if errors
-
-   MONTHS = ['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG', $
-             'SEP','OCT','NOV','DEC']
-
-   IF NOT keyword_set(ndayspm)  THEN BEGIN 
+function JULDAY, 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
+
 
 ; Gregorian Calander was adopted on Oct. 15, 1582
-      GREG = 15L + 31L * (10L + 12L * 1582L)
+; 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 begin
-         DATE = systime()
-         L_MONTH = long(where(strupcase(strmid(DATE, 4, 3)) eq MONTHS))
-         L_MONTH = L_MONTH[0] + 1 ; Scalarize it...
-         L_DAY = long(strmid(DATE, 8, 2))
-         L_YEAR = long(strmid(DATE, 20, 4))
-      endif else if np ge 3 then begin
-         L_MONTH = LONG(MONTH)
-         L_DAY = LONG(DAY)
-         L_YEAR=LONG(YEAR)
-;
-;***************************************************
-; Change test to allow year 0 for climatological data
-;********************************************************
-;         if (L_YEAR eq 0) then message, 'There is no year zero.'
-      endif else message, 'Wrong number of parameters.'
-;
-;***************************************************
-; Change test to allow year 0 for climatological data
-;********************************************************
-;      if (L_YEAR lt 0) then L_YEAR = L_YEAR + 1
-      if (L_YEAR le 0) then L_YEAR = L_YEAR + 1
-      if (L_MONTH gt 2) then begin
-         JY = L_YEAR
-         JM = L_MONTH + 1
-      endif else begin
-         JY = L_YEAR - 1
-         JM = L_MONTH + 13
-      endelse
-
-      JUL = floor(365.25 * JY) + floor(30.6001 * JM) + L_DAY + 1720995
+      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 ((L_DAY + 31L * (L_MONTH + 12L * L_YEAR)) ge GREG) then begin
-         JA = long(0.01 * JY)
-         JUL = JUL + 2 - JA + long(0.25 * JA)
-      endif
-
-      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
-
-      return, JUL + (Hour / 24.0d0 - .5d0) + (Minute/1440.0d0) + (Second / 86400.0d0)
-
-   ENDIF ELSE BEGIN 
+      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 ndayspm EQ 1 THEN ndayspm = 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)
-
-      JUL = ((L_YEAR-1)*12. + (L_MONTH-1))* ndayspm + L_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
+        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
       
-      return, JUL + (Hour / 24.0d0) + (Minute/1440.0d0) + (Second / 86400.0d0)
-
-   ENDELSE 
-end
+      IF Hour+Minute+Second EQ 0 THEN return, JUL ELSE $
+        return, JUL + (Hour / 24.0d0) + (Minute/1440.0d0) + (Second / 86400.0d0)
+
+    END 
+    'noleap':BEGIN
+    END 
+    ELSE:return, report('only 3 types of calendar are accepted: greg, 360d and noleap')
+  ENDCASE
+
+END