C+
C NAME:
C	Julian
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:
C	subroutine Julian(IDIN,iYr,Doy,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 
c                                       (=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	Doy	real*8		day of year (including fraction for the time of
C                                            day).
C	JD	real*8		Julian day
C	JEpoch	real*8		Julian epoch = time in Julian years
C	BEpoch	real*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				(Doy = 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-
      SUBROUTINE julian(idin,iyr,doy,jdio,jepoch)

      INTEGER	idin
      integer	iyr
      REAL*8	doy
      REAL*8	jdio
      REAL*8	jepoch

      PARAMETER	(i365  =  365)
      PARAMETER	(i1582 = -418)	!1582-2000
      REAL*8		jd
      REAL*8		bepoch
      REAL*8		jd2000
      REAL*8            frac

      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+dble(doy-1.5)	! 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. Doy <= 0 happens for dates earlier
! than 2000 Jan. 1.0; Doy >=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

       doy = idoy+1+frac
       jd = jd-0.5d0

      END IF

      jdio = jd2000+jd

      RETURN
      END