C+
C NAME:
C	Julian8
C PURPOSE:
C	Of the four quantities, date, Julian day number, Julian and Besselian
C	epoch, one is supplied. The other three are then calculated.
C CATEGORY:
C	Time keeping
C CALLING SEQUENCE:
	subroutine Julian8(IDIN,iYr,Doy8,JDio,JEpoch)
C INPUTS:
C	IDIN	integer   ID=0  date    --> JD, Julian epoch, Besselian epoch
C				=1  JD      --> Date, Julian epoch, Besselian epoch
C				=2  Julian epoch    --> Date, JD, Besselian epoch
C	(not active)	=3  Besselian epoch --> Date, JD, Julian epoch
C				Add 10: use modified Julian days (=JD-2400000.5)
C				Add 20: use days relative to Jan. 1, 2000 (=JD-2451544.5)
C  INPUTS/OUTPUTS: (depending on value of ID)
C	iYr	integer		year; the year xxxBC should be entered as -xxx+1.
C	Doy8	real*8		day of year (including fraction for the time of day).
C	JD	real*8		Julian day
C	JEpochreal*8		Julian epoch = time in Julian years
C	BEpochreal*8		Besselian epoch
C PROCEDURE:
C	Dates before 5 october 1582 are interpreted as Julian dates;
C	after 15 october as Gregorian dates.
C	5 october 1582 (Julian) = 15 october 1582 (Gregorian).
C
C	RJD	'relative' Julian day = (Julian day - 2451544.5)      
C	J2000.0 = 2000 January 1.5d TDB = JD 2451545.0
C	JD 2451544.5 = 2000 January 1.0d (i.e. 1 January, midnight)
C	Check:	JD 2398220.0 =  1854 Jan.  1.5d
C		JD 2299161.0 =  1582 Oct. 15.5d (Gregorian calender)
C			     =  1582 Oct.  5.5d (Julian calender)
C				(Doy8 = 278.5)
C		JD       0.0 = -4712 Jan.  1.5d
C
C	JEpoch 1900 = JD 2415020.0 = Jan 0.5 1900
C	JEpoch 2000 = JD 2451545.0 = Jan 1.5 2000
C
C	BEpoch = 1900+(JD -2415020.31352)/365.242198781
C	       = 1900+(RJD+  36524.18648)/365.242198781
C	Check:  BEpoch 1950.0 = JD 2433282.423 = Jan 0.923 1950
C		BEpoch 1982.0 = JD 2444970.174 = Jan 0.674 1982
C
C MODIFICATION HISTORY:
C	1989-1990, Paul Hick (MPAE,UCSD/CASS; pphick@ucsd.edu)
C-
	    integer	IDIN
	    integer	iYr
	    real*8	Doy8
	    real*8	JDio
	    real*8	JEpoch
	    real*8	Frac

	parameter	(I365  =  365)
	parameter	(I1582 = -418)	!1582-2000
	real*8		JD
	real*8		BEpoch
	real*8		JD2000

	JD2000 = 2451545	! Jan 1.5, 2000 (noon)
	if (IDIN .ge. 10) JD2000 = 51544.5d0
	if (IDIN .ge. 20) JD2000 = 0.5d0

	ID = mod(IDIN,10)	! ID=0,1,2,3

	if (ID .eq. 0) then	! Date ---> Julian day, etc.
	    Laps = iYr-2000	! Whole years from 1 Jan 1.0d 2000 to idem of iYr
	    I = Laps
	    if (Laps .gt. 0) I = Laps-1
	    J = max(I1582,I)
	    Leap = I/4-J/100+J/400
	    if (Laps .gt. 0) Leap = Leap+1	! # leap years
	    JD = I365*Laps+Leap+(Doy8-1.5D0)	! Rel. to 2000 Jan. 1.5 (noon)
	    if (iYr .le. 1582) JD = JD+10

	    JEpoch = 2000+JD/365.25
	    BEpoch = 1900+(JD+36524.68648)/365.242198781

	else				! Julian day, etc. ---> Date
	    if (ID .eq. 2) then		! Julian epoch
		JD = JEpoch-2000
		BEpoch = 2000.0012775+1.000021359*JD
		JD = JD*365.25
	    else if (ID .eq. 3) then
		JD = BEpoch-1900
		JEpoch = 1900.000858+.999978642*JD
		JD = -36524.68648+JD*365.242198781
	    else			! Assume ID=1
		JD = JDio-JD2000
		JEpoch = 2000+JD/365.25
		BEpoch = 1900+(JD+36524.68648)/365.242198781
	    end if

	    JD = JD+0.5D0		! Shift to previous midnight
	    JD0 = int(JD)		! Whole days between previous midnight ...
	    if (JD .lt. JD0) JD0 = JD0-1! .. and midnight of 1 Jan 2000
	    Frac = JD-JD0		! Day fraction from previous midnight
	    Laps = 0
	    iDoy = JD0
	    do while (iabs(iDoy) .gt. 365)
		Laps = Laps+iDoy/366	! Underestimate number of full years in iDoy
		I = Laps
		if (Laps .gt. 0) I = Laps-1
		J = max(I1582,I)
		Leap = I/4-J/100+J/400
		if (Laps .gt. 0) Leap = Leap+1	! # leap years
		iDoy = JD0-Laps*I365-Leap
		if (Laps .le. I1582) iDoy = iDoy-10
	    end do

	    !-----
	    ! Leave do while loop with -365<=iDoy<=365. Doy8 <= 0 happens for dates earlier
	    ! than 2000 Jan. 1.0; Doy8 >=0 for dates later than 2000 Jan 1.0.

	    iYr = 2000+Laps
	    if (iDoy .lt. 0) then
		iYr = iYr-1
		if (iYr .eq. 1582) then
		    iDoy = iDoy+355	! 355: # days in 1582
		    if (iDoy .lt. 0) iYr = iYr-1
		end if
	    end if

	    Leap = 0			! Suppose it's not a leap year
	    if (iYr/4*4 .eq. iYr .and. (iYr .lt. 1582 .or. !It's a leap year
     &		iYr/100*100 .ne. iYr .or. iYr/400*400 .eq. iYr)) Leap = 1
	    if (iDoy .lt. 0) iDoy = 365+Leap+iDoy	! Now iDoy>=0

	    if (iDoy .eq. 365+Leap) then
		iDoy = 0
		iYr = iYr+1
	    end if

	    Doy8 = iDoy+1+Frac
	    JD = JD-0.5D0

	end if

	JDio = JD2000+JD

	return
	end