;+
; NAME:
;	TimeYDoy
; PURPOSE:
;	Converts date in yr, doy to number of days elapsed 
;	since 2000 Jan 1, 0 UT and v.v.
; CATEGORY:
;	gen/idl/toolbox/time
; CALLING SEQUENCE:
	PRO TimeYDoy, yr, doy, day, get_date=get_date, julian=julian, jul2greg=jul2greg
; INPUTS:
;	/get_date NOT set:
;
;	yr			scalar or array; type: integer
;					year
;	doy 		scalar or array; type: any
;					day of year, possibly including a fraction for
;					the time of day.
;
;	/get_date SET:
;
;	day 		scalar or array: type: any
;					number of days elapsed since 2000 Jan 1, UT (Gregorian)
;					possibly including a fraction for the time of day.
;
;					Use TimeShift(d2000=TimeSet(day=day), /to_time) to convert
;					to a time structure (relative to current origin).
;
; OPTIONAL INPUT PARAMETERS:
;	/get_date	by default the number of days elapsed is calculated from
;				the date. If /get_date is set than the date is obtained
;				from the number of days elapsed since 2000 Jan 1, 0 UT.
;				(Gregorian).
;	/julian 	interpret yr, doy as date in Julian calendar
;	/jul2greg
;	jul2greg=jul2greg
;				array[1]; type: time structure
;					time at which the transition from Julian to Gregorian
;					calendar is made.
;				/jul2greg is the same as jul2greg=TimeSet(yr=1582, mon='oct', day=15)
; OUTPUTS:
;	/get_date NOT set:
;
;	day 		scalar or array: type: double
;					number of days elapsed since 200 Jan 1, UT (Gregorian).
;					the fraction of day is copied from the input day of year
;					(if present)
;
;	/get_date SET:
;
;	yr			scalar or array; type: integer
;					year
;	doy 		scalar or array; type: same as input 'day'
;					day of year, the fraction of day is copied from the
;					 input day (if present)
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, IsType, IsTime, TimeSet, TimeYDoy
; SIDE EFFECTS:
;	If /get_date NOT set then yr is converted to integer if it isn't already.
; PROCEDURE:
;	The date yr,doy is a Julian date (if /julian SET) or a Gregorian date
;	(/julian NOT set). 'day' always is # days since 2000 Jan 1.0 (Gregorian)
;	i.e. Julian days minus 2451544.5
;
;
;	JD = 2451557.5  = 2000 Jan 14.0 (Gregorian) = 2000 Jan 1.0 (Julian)
;	JD = 2451544.5  = 2000 Jan  1.0 (Gregorian)
;	JD = 2299160.5  = 1582 Oct 15.0 (Gregorian) = 1582 Oct 5.0 (Julian)
;   		(= -152384.0 days from 2000 Jan 1.0, Gregorian)
;	JD = 2361221.5  = 1752 Sep 14.0 (Gregorian) = 1752 Sep 3.0 (Julian)
;
;	The keyword jul2greg can be used to automatically switch from
;	Julian to Gregorian calendar. In general jul2greg is defined as
;	a UT standard time structure. If specified as keyword /jul2greg
;	then the time 1582 Oct 15.0 (Gregorian is used). The gap between
;	calendars at that time was 10 days.
;
;	Another useful transition would be 1752 Sep 14.0 (Gregorian).
;	The Britisch empire switched from 1752, Sep 2 to 1752 Sep 14
;	(11 day gap; the extra day comes from the year 1700 which is a leap
;	year in the Julian calendar but not in the Gregorian calendar).
; MODIFICATION HISTORY:
;	JAN-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, get_date, /key
InitVar, julian  , /key

IF IsType(jul2greg, /defined) THEN BEGIN

	IF NOT IsTime(jul2greg) THEN jul2greg = TimeSet(yr=1582,month='oct',day=15)

	; Assume that all dates are Gregorian

	IF NOT get_date THEN TimeYDoy, yr, doy, day

	; If the Gregorian interpretation leads to a date
	; preceeding the swith than those dates need to be interpreted as
	; Julian dates.

	TimeYDoy, yr, doy, day, get_date=get_date, julian=day LT TimeGet(/d2000, jul2greg, /scalar)

	RETURN

ENDIF


CASE get_date OF

0: BEGIN							; yr, doy --> day (days since 2000 Jan 1, 0 UT)

	yr = long(yr)

	L = yr-2000						; Whole years from 1 Jan 1.0d 2000 to same of yr
	C = L GT 0
	A = L-C
	B = (1-julian)*A

	day = 365*L+A/4-B/100+B/400+C

	; Day now is the number of days since 2000 Jan 1.0 of the Julian or
	; Gregorian calendar. For the Julian calendar add 13 days
	; so always the number of days since 2000 Jan 1.0 (Gregorian) is returned.

	day += 13*julian

	IF IsType(doy, /generic_int) THEN day += doy-1 ELSE IF (where(doy NE 1))[0] NE -1 THEN day += double(doy-1)

END

1: BEGIN							; day --> yr, doy

	d = day-13*julian

	N = n_elements(d)

	A = floor(d)					; Whole days between previous midnight and midnight of 1 Jan 2000
	f = d-A
	d = A							; Positive day fraction from previous midnight

	yr	= 0*d
	doy = d

	I = lindgen(N)
	B = I
	REPEAT BEGIN
		I		= I[B]
		yr[I]  += doy[I]/366	; Underestimate # full years in IDOY
		A		= yr[I]
		C		= A GT 0
		A	   -= C
		B		= (1-julian)*A
		doy[I]	= d[I]-yr[I]*365L-(A/4-B/100+B/400+C)

		B		= where(abs(doy[I]) GT 365)
	ENDREP UNTIL B[0] EQ -1

	; The repeat loop is left with -365<=Doy<=365. Doy is positive/negative for dates
	; after/before 2000/1/1. The input Julian day corresponds to the day Doy after the
	; start of year Yr (i.e. if Doy lt 0 the date lies in the year Yr-1).

	yr += 2000
	A	= doy LT 0
	yr -= A

	; At this point Yr correctly identifies the year. Where Doy<0, add #(days in year).
	; This makes Doy >=0 everwhere. Note that for Yr=1582, Doy>=0 already, so the
	; correction for the number of days (only 355 in 1582) is not necessary.

	B	 = 365+(yr/4*4 EQ yr AND (julian OR yr/100*100 NE yr OR yr/400*400 EQ yr))
	doy += A*B

	; Check where Doy contains a full year (these are the non-leap years
	; where Doy=+365 after the repeat loop)

	A	 = doy EQ B
	yr	+= A
	doy -= A*B

	; Add 1 to start counting at doy=1 instead of doy=0
	; Restore fraction for time of day

	doy += 1
	IF (where(f NE 0))[0] NE -1 THEN doy += double(f)

END

ENDCASE

RETURN  &  END