 	SUBROUTINE T68_TO_GMT(T68,GMT)
C
C	CONVERT T68 REAL*8 SECONDS SINCE 5/24/68 (JD 2440000.5)
C	INTO GMT ASCII FORMAT 'YYDDDHHMMSSFFF'
C
	REAL*8 T68
	CHARACTER*14 GMT
	INTEGER IY,IM,ID,IH,MINUTE,I367,JD,JD11,IS,I
	REAL*4 SECOND
	DATA I367/367/
	CALL T68_TO_UTC(T68,IY,IM,ID,IH,MINUTE,SECOND)
	IF (IH.GE.24) THEN
	  IH = IH-24
	  ID = ID+1
	ENDIF
	JD=-7*(IY+(IM+9)/12)/4-3*((IY+(IM-9)/7)/100+1)/4+275*IM/9+ID
	JD11=-7*IY/4-3*((IY-1)/100+1)/4+30
	ID = JD-JD11
	IY = MOD(IY,100)
	IS = 1000.0*SECOND
	WRITE (GMT,FMT='(I2,I3,2I2,I5)') IY,ID,IH,MINUTE,IS
	DO I=1,14
	  IF (GMT(I:I).EQ.' ') GMT(I:I) = '0'
	ENDDO
	RETURN
	END

	REAL*8 FUNCTION GMT2T68(GMT)
C
C	CONVERT THE STRING IN GMT 'YYDDDHHMMSSFFF' FORMAT INTO
C	T68 REAL*8 SECONDS SINCE 5/24/68 00:00:00 UTC FORMAT
C	IF GMT = 'NOW' GET CURRENT SYSTEM TIME AND CONVERT THAT
C
	CHARACTER*14 GMT
	CHARACTER*24 FDATEOUT
	INTEGER IY,ID,MON,IH,MINUTE
	REAL*4 SECOND
 	REAL*8 T68
	IF (INDEX(GMT,'NOW').GT.0) THEN
C	version for Unix library
	  CALL FDATE(FDATEOUT)
	  MON = (INDEX(
	1	' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec '
	2	,FDATEOUT(4:8))-1)/4+1
	  DO WHILE (INDEX(FDATEOUT,':').GT.0)
	    I = INDEX(FDATEOUT,':')
	    FDATEOUT(I:I) = ' '
	  ENDDO
	  READ (FDATEOUT(8:24),*) 
	1	ID, IH, MINUTE, IS, IY 
	  SECOND = IS
	ELSE
 	  READ (GMT,FMT='(I2,I3,2I2,F5.3)') IY,ID,IH,MINUTE,SECOND
	  IY = 1900+IY
	  IF (IY.LT.1968) IY=IY+100
	  MON = 1
	ENDIF
	CALL TIME_TAG(IY,MON,ID,IH,MINUTE,SECOND,T68)
	GMT2T68 = T68
	RETURN
	END

	SUBROUTINE T68_TO_UTC(T,IY,MON,ID,IH,MINUTE,SECOND)
C
C	CONVERTS TIME TAG T INTO UT TIME IH:MINUTE:SECOND ON IM/ID/IY
C	T is seconds since 5/24/68 allowing for leap seconds
C
	REAL*8 T,TT,DT,DTT
	INTEGER*4 MAXDAY(12)
	DATA MAXDAY/31,28,31,30,31,30,31,31,30,31,30,31/
	DAY = T/864.D2+140
	IY = DAY/365.25+1968
	DAY = MOD(DAY,365.25)
	MON = DAY/30.5+1.
	ID = MOD(DAY,30.5)+1.
10	MAXD = MAXDAY(MON)
C	check for Leap year Februaries:
	IF (MON.EQ.2 .AND. ((MOD(IY,4).EQ.0 .AND. MOD(IY,100).NE.0) .OR.
	1	MOD(IY,400).EQ.0) ) MAXD=29
	IF (ID.GT.MAXD) THEN
	  ID = ID-MAXD
	  MON = MON+1
	  IF (MON.GT.12) THEN
	    MON = 1
	    IY = IY+1
	  ENDIF
	  GO TO 10
	ENDIF
	CALL TIME_TAG(IY,MON,ID,0,0,1.,TT)
	DT = T-TT+1.D0
C	MAKE SURE THE DAY IS RIGHT
	IF (DT.LT.86400.D0) THEN
	  IH = DT/36.D2
	  DT = DT-36.D2*IH
	  MINUTE = DT/60.D0
	  SECOND = DT-60.D0*MINUTE
	  RETURN
	ENDIF
C	CHECK FOR A LEAP SECOND
	IF (DT.LT.86401.D0) THEN
	  CALL TIME_TAG(IY,MON,ID+1,0,0,1.,TT)
	  DTT = T-TT+1.D0
	  IF (DTT.GE.0.) THEN
	    ID = ID+1
	    GO TO 10
	  ENDIF
	  IH = 23
	  MINUTE = 59
	  SECOND = DT - 86340.D0
	  RETURN
	ENDIF
	ID = ID+DT/86401.D0
	GO TO 10
	END

	SUBROUTINE TIME_TAG(IY,IM,ID,IH,MINUTE,SECOND,T)
C
C	CONVERTS YEAR, MONTH, DAY, HOUR, MINUTE, SECOND INTO
C	REAL*8 TIMETAG T, SECONDS SINCE 5/24/68 00:00:00 UT
C	The time tag T is REAL*8 seconds since 24 May 1968, UT 00:00:00
C	which was JD 2440000.5
C
C	The seconds to be used are ephemris seconds, or IAT seconds
C	Therefore the conversion from UT date and time must allow for the
C	rates adjustments used before 19?? and the leap seconds used later
C	This routine first calculates UTC ignoring these corrections, then
C	looks in a table to interpolate for DT = ET - UT
C	Linear interpolation is used, so two entries are made for each
C	leap second.  Entry for JDATE should be the JD at noon the day 
C	after the leap second.
C
	REAL*8 T,DT
	REAL*4 SECOND,DELTA_T
	INTEGER*4 JD,I367,JULIAN,JDZERO
	DATA JDZERO/2440001/			!JD at 5/24/68 UT noon
	DATA I367/367/
	JD=-7*(IY+(IM+9)/12)/4-3*((IY+(IM-9)/7)/100+1)/4+275*IM/9+ID
	JULIAN=JD+I367*IY+1721029
C	lookup DT at UT noon to avoid leap seconds
	T=(JULIAN-JDZERO)*864.D2+432.D2
	DT = DELTA_T(T)
	T=(JULIAN-JDZERO)*864.D2+IH*36.D2+MINUTE*6.D1+SECOND+DT
	RETURN
	END

	REAL*4 FUNCTION DELTA_T(T68)
C
C	T68 is REAL*8 seconds since 24 May 1968, UT 00:00:00
C	which was JD 2440000.5
C
C	DELTA_T is the difference between IAT (or ET or TDT) and UTC
C	since 5/24/68
C
C	The seconds to be used are ephemeris seconds, or IAT seconds
C	Therefore the conversion from UT date and time must allow for the
C	rates adjustments used before 1972 and the leap seconds used later
C	This routine first calculates UTC ignoring these corrections, then
C	looks in a table to interpolate for DT = ET - UT
C	Linear interpolation is used, so two entries are made for each
C	leap second.  Entry for JDATE should be the JD at noon the day 
C	after the leap second.
C
	IMPLICIT NONE
	INTEGER*4 NPT
	PARAMETER (NPT=24)
	INTEGER*4 JDATE(NPT),I,NL,NH,NTRY
	REAL*8 DT(NPT),DTNOW,SLOPE
	LOGICAL SAVED/.FALSE./
	INTEGER*4 JD68/2440001/
	REAL*8 DT68/38.66/
	REAL*8 T68,TAB68(NPT)

	SAVE SAVED,TAB68

	DATA (JDATE(I),DT(I),I=1,NPT)/
	1		2440001,	38.66,	      	!	5/24/1968
	1		2445517,	53.18,		!	7/1/1983
	1		2445517,	54.18,		!	7/1/1983
	1		2446248,	54.18,		!	7/1/1985
	1		2446248,	55.18,		!	7/1/1985
	1		2447162,	55.18,		!       1/1/1988
	1		2447162,	56.18,		!	1/1/1988
	1		2447893,	56.18,		!	1/1/1990
	1		2447893,	57.18,		!	1/1/1990
	1		2448258,	57.18,		!	1/1/1991
	1		2448258,	58.18,		!	1/1/1991
	1		2448805,	58.18,		!	7/1/1992
	1		2448805,	59.18,		!	7/1/1992
	1		2449170,	59.18,		!     	7/1/1993
	1		2449170,	60.18,		!     	7/1/1993
	1		2449535,	60.18,		!	7/1/1994
	1		2449535,	61.18,		!	7/1/1994
	1		2450084,	61.18,		!	1/1/1996
	1		2450084,	62.18,		!	1/1/1996
	1		2450631,	62.18,		!	7/1/1997
	1		2450631,	63.18,		!	7/1/1997
	1		2451180,	63.18,		!	1/1/1999
	1		2451180,	64.18,		!	1/1/1999
	1		2455198,	64.18/		!	1/1/2010
C
C---------------------------------------------------
	IF (.NOT. SAVED) THEN
C				 YYDDDHHMMSSFFF
	  SAVED = .TRUE.
	  DO I=1,NPT
	    TAB68(I)=(JDATE(I)-JD68)*864.D2 + (DT(I)-DT68)
	  ENDDO
	ENDIF
C
C	Find position in Table
C
	NL=1
	NH=NPT-1
	DO WHILE (NL.LT.NH)
	  NTRY=(NL+NH+1)/2
	  IF(TAB68(NTRY).GT.T68) THEN
	    NH=NTRY-1
	  ELSE
	    NL=NTRY
	  ENDIF
	ENDDO
	SLOPE = (DT(NL+1)-DT(NL))/(TAB68(NL+1)-TAB68(NL))
	DTNOW=DT(NL)+ (T68-TAB68(NL)) * SLOPE
	DELTA_T = (DTNOW-DT68)	! time correction ET-UTC in seconds
	RETURN
	END

	SUBROUTINE TTAG(IY,IM,ID,IH,MINUTE,SECOND,T)
C
	INTEGER IY,ID,IM,IH,MINUTE
	REAL*4 SECOND
	REAL*8 T
C
	CALL TIME_TAG(IY,IM,ID,IH,MINUTE,SECOND,T)
	RETURN
	END

	SUBROUTINE UTC(T,IY,MON,ID,IH,MINUTE,SECOND)
C
	REAL*8 T
	INTEGER IY,MON,ID,IH,MINUTE
	REAL*4 SECOND
C
	CALL T68_TO_UTC(T,IY,MON,ID,IH,MINUTE,SECOND)
	END

	REAL*8 FUNCTION JULIAN_DAY(T)
	REAL*8 T,JDZERO,DTZERO
C
C	RETURNS EPHEMERIS JULIAN DAY, GIVEN TIME TAG T
C
	DATA JDZERO/2440001.D0/		!JD at 5/24/68 UT noon	TTTTTT
	DATA DTZERO/38.65D0/		!ET-UT at 5/24/68
	JULIAN_DAY=JDZERO+(T+DTZERO)/864.D2-0.5D0
	RETURN
	END

	SUBROUTINE GET_EPOCH(T,EPOCH)
	REAL*8 T,JULIAN_DAY
	EPOCH=2000.0+(JULIAN_DAY(T)-2451545.D0)/365.25D0
	RETURN
	END

	FUNCTION DOT(A,B)
	DIMENSION A(3),B(3)
	DOT=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
	RETURN
	END

	SUBROUTINE NORM(A)
C
C	NORMALIZES A THREE VECTOR
C
	DIMENSION A(3)
	AM=MAX(ABS(A(1)),ABS(A(2)),ABS(A(3)))
	IF(AM.GT.1.E10) THEN
	  DO I=1,3
	    A(I)=A(I)/AM
	  ENDDO
	ENDIF
	AA=SQRT(A(1)**2+A(2)**2+A(3)**2)
	DO I=1,3
	  A(I)=A(I)/AA
	ENDDO
	RETURN
	END

	SUBROUTINE CROSS(A,B,C)
	DIMENSION A(3),B(3),C(3)
C
C	CROSS PRODUCT C=A X B
C
	C(3)=A(1)*B(2)-A(2)*B(1)
	C(1)=A(2)*B(3)-A(3)*B(2)
	C(2)=A(3)*B(1)-A(1)*B(3)
	RETURN
	END

	SUBROUTINE MAV(A,VSC,VCEL)
C
C	MULTIPLIES MATRIX A TIMES VECTOR VSC GIVING VECTOR VCEL
C
	REAL*4 A(3,3),VSC(3),VCEL(3),SC(3),VC
	DO 5 I=1,3
5	SC(I)=VSC(I)
	DO 15 I=1,3
	VC=0.
	DO 10 J=1,3
10	VC=VC+A(I,J)*SC(J)
15	VCEL(I)=VC
	RETURN
	END

	SUBROUTINE VMAT(A,VCEL,VSC)
C
C	MULTIPLIES TRANSPOSE OF A TIMES VCEL GIVING VSC
C
	REAL*4 A(3,3),VCEL(3),VSC(3),CEL(3),SC
	DO 5 I=1,3
5	CEL(I)=VCEL(I)
	DO 15 I=1,3
	SC=0.
	DO 10 J=1,3
10	SC=SC+A(J,I)*CEL(J)
15	VSC(I)=SC
	RETURN
	END

	FUNCTION GAUSRN(IR,JR)
	EXTERNAL RAN
C       ISW not equal to 0 if we have a precomputed value
	DATA ISW/0/
	SAVE Y,F,ISW
	IF(ISW.EQ.0) GO TO 10
	GAUSRN=Y*F
	ISW=0
	RETURN
C       compute two Gaussian random variables
C       first find X,Y uniform in the unit circle
10	X=2.*RAN(IR,JR)-1.
	Y=2.*RAN(IR,JR)-1.
	RR=X*X+Y*Y
	IF(RR.GE.1.) GO TO 10
C       prob of r < Rg for 2-D Gaussian is 1-exp(-0.5*Rg^2)
C       prob of r < Ru for uniform on unit circle is Ru^2
C       set probabilities equal to transform from unit circle to 2-D gaussian
C       giving Rg = SQRT(-2ln(1-Ru^2)), so F = Rg/Ru as given below
	IF (RR.LT.0.01) THEN
	  F = SQRT(2.*(1+RR*(0.5+RR/3.)))
	ELSE
	  F = SQRT(-2.*ALOG(1.-RR)/RR)
	ENDIF
C       first Gaussian random variable is F*X, with F*Y saved for next call
	GAUSRN=F*X
	ISW=1
	RETURN
	END

	SUBROUTINE ETOC(E,C)
	REAL*4 E(3),C(3),TEMP
C
C	E	INPUT	vector in ecliptic coordinates (J2000.0)
C	C	OUTPUT	vector in celestial coordinates (J2000.0)
C
	PARAMETER (COBL=.917482058,SOBL=.39777716)
	C(1)=E(1)
	TEMP=COBL*E(3)+SOBL*E(2)
	C(2)=COBL*E(2)-SOBL*E(3)
	C(3)=TEMP
	RETURN
	END

	SUBROUTINE CTOE(C,E)
	REAL*4 C(3),E(3),TEMP
C
C	C	INPUT	vector in celestial coordinates
C	E	OUTPUT	vector in ecliptic coordinates
C
	PARAMETER (COBL=.917482058,SOBL=.39777716)
	E(1)=C(1)
	TEMP=COBL*C(3)-SOBL*C(2)
	E(2)=COBL*C(2)+SOBL*C(3)
	E(3)=TEMP
	RETURN
	END

	SUBROUTINE PREMAT(EPOCH,R)
C
C	COMPUTES MATRIX FOR PRECESSION FROM EPOCH TO 2000.
C
C	V(EPOCH)=(I+R)V(2000)
C
C	             T
C	V(2000)=(I+R) V(EPOCH)
C
	REAL*4 R(3,3),RE(3,3),R2(3,3)
	INTEGER IFF
	DATA IFF /0/
	SAVE IFF,R2
	IF (IFF.EQ.0) THEN
	  IFF = 1
	  CALL PREMAT50(2000.0,R2)
	ENDIF
	CALL PREMAT50(EPOCH,RE)
C						      T
C	V(EPOCH)=(I+R)V(2000)=(I+RE)V(50)=(I+RE)(I+R2) V(2000)
C		      T	       T
C	SO R = RE + R2  + RE*R2
C
	DO 30 J=1,3
	  DO 20 I=1,3
	    R(I,J)=RE(I,J)+R2(J,I)
	    DO 10 K=1,3
	      R(I,J) = R(I,J) + RE(I,K)*R2(J,K)
10	    CONTINUE
20	  CONTINUE
30	CONTINUE
	RETURN
	END

	SUBROUTINE PREMAT50(EPOCH,R)
C
C	COMPUTES MATRIX FOR PRECESSION FROM EPOCH TO 1950.
C
C	V(EPOCH)=(I+R)V(50)
C
C	           T
C	V(50)=(I+R) V(EPOCH)
C
	REAL*4 R(3,3),AP(3,3),BP(3,3),CP(3,3),P(3,3)
	DTR=45./ATAN(1.)
	T=(EPOCH-1950.0)/100.
	ZETA=T*(.6402633+T*(.0000839+T*.000005))/DTR
	Z=ZETA+.0002197*T**2/DTR
	THETA=T*(0.5567376-T*(.0001183+T*.0000117))/DTR
	DO 10 I=1,3
	AP(I,3)=0.
	AP(3,I)=0.
	CP(I,3)=0. 
	CP(3,I)=0.
	BP(I,2)=0.
10	BP(2,I)=0.
	AP(1,1)=-2.*SIN(Z/2.)**2
C
C	COS(Z)-1 = -2*SIN(Z/2)**2
C
	AP(2,2)=AP(1,1)
	AP(1,2)=-SIN(Z)
	AP(2,1)=-AP(1,2)
	BP(1,1)=-2.*SIN(THETA/2.)**2
	BP(3,3)=BP(1,1)
	BP(1,3)=-SIN(THETA)
	BP(3,1)=-BP(1,3)
	CP(1,1)=-2.*SIN(ZETA/2.)**2
	CP(2,2)=CP(1,1)
	CP(1,2)=-SIN(ZETA)
	CP(2,1)=-CP(1,2)
C
C	AP=A-I, BP=B-I, CP=C-I
C	NOW COMPUTE BC-I=AP*BP+AP+BP
C
	DO 30 I=1,3
	DO 30 J=1,3
	S=0.
	DO 20 K=1,3
20	S=S+BP(I,K)*CP(K,J)
30	P(I,J)=S+BP(I,J)+CP(I,J)
C
C	NOW COMPUTE R=ABC-I=AP*(BC-I)+AP+(BC-I)
C
	DO 50 I=1,3
	DO 50 J=1,3
	S=0.
	DO 40 K=1,3
40	S=S+AP(I,K)*P(K,J)
50	R(I,J)=S+AP(I,J)+P(I,J)
	RETURN
	END

	SUBROUTINE CTOG(C2000,G)
C
C	CONVERTS VECTOR C FROM CELESTIAL COORDINATES (2000.0)
C	TO GALACTIC COORDINATE VECTOR G
C
	DIMENSION G(3),C2000(3),R(3,3)
	DATA R/ -0.0548809,  0.4941080, -0.8676667,
	1       -0.8734369, -0.4448323, -0.1980717,
 	2       -0.4838349,  0.7469817,  0.4559848/
	CALL MAV(R,C2000,G)
	RETURN
	ENTRY GTOC(G,C2000)
C
C	CONVERTS GALACTIC IN G TO CELESTIAL IN C
C
	CALL VMAT(R,G,C2000)
	RETURN
	END

	SUBROUTINE C50TOG(C50,G)
C
C	CONVERTS VECTOR C FROM CELESTIAL COORDINATES (1950.0)
C	TO GALACTIC COORDINATE VECTOR G
C
	DIMENSION G(3),C50(3),R(3,3)
	DATA R/-.0669888,.4927285,-.8676008,
	1       -.8727558,-.450347,-.1883746,
	2       -.4835389,.7445847,.4601998/
	CALL MAV(R,C50,G)
	RETURN
	ENTRY GTOC50(G,C50)
C
C	CONVERTS GALACTIC IN G TO CELESTIAL IN C
C
	CALL VMAT(R,G,C50)
	RETURN
	END

	FUNCTION DVS(A,B)
C
C	COMPUTES DISTANCE BETWEEN VECTORS SQUARED
C	DVS=(A-B)**2
C
	DIMENSION A(3),B(3)
	DVS=(A(1)-B(1))**2+(A(2)-B(2))**2+(A(3)-B(3))**2
	RETURN
	END

	FUNCTION DMR_BEAM(SKY,CENTER)
C
C	APPROXIMATES DMR BEAM SHAPE.  VALUE=1 IN THE CENTER
C	FWHM = 7 DEGREES, VALUE = 10^-4 AT 25 DEGREES
C
	PARAMETER (XE = 0.005381, XS = 0.003798805)
	X = DVS(SKY,CENTER)
C	change from quadratic to linear at XS
	IF (X.GT.XS) X = XS*(2.0*SQRT(X/XS)-1.0)
	DMR_BEAM = EXP(-X/XE)
	RETURN
	END
	
	FUNCTION BEAM(SKY,CENTER,WIDTH)
C
C	CALULATES A TEST FUNCTION BEAM SHAPE
C	WIDTH IS THE HWZI IN RADIANS
C
	DIMENSION SKY(3),CENTER(3)
	X=DVS(SKY,CENTER)/WIDTH**2
	BEAM=0.
	IF(X.LT.0.98) BEAM=EXP(-X/(1.-X))
	RETURN
	END

C	QUATRNION.FTN	ROUTINES TO MANIPULATE QUATERNIONS

	SUBROUTINE QXQ(P,Q,R)
C
C	PRODUCT OF TWO QUATERNIONS P*Q -> R
	REAL*4 P(4),Q(4),R(4)
C
C	ORDER OF COMPONENTS:	i,j,k,1
C	Reminder: 1*1=1; i*i=j*j=k*k=-1
C	1*i=i;  1*j=j;  1*k=k;
C	i*j=-j*i=k;  j*k=-k*j=i; k*i=-i*k=j
C
	R(1)=P(4)*Q(1)+Q(4)*P(1)+P(2)*Q(3)-Q(2)*P(3)
	R(2)=P(4)*Q(2)+Q(4)*P(2)+P(3)*Q(1)-Q(3)*P(1)
	R(3)=P(4)*Q(3)+Q(4)*P(3)+P(1)*Q(2)-Q(1)*P(2)
	R(4)=P(4)*Q(4)-P(1)*Q(1)-P(2)*Q(2)-P(3)*Q(3)
	RETURN
	END

	SUBROUTINE QINV(P,PINV)
C
C	COMPUTES 1/P -> PINV FOR QUATERNIONS
C
	REAL*4 P(4),PINV(4)
	PP=P(1)**2+P(2)**2+P(3)**2+P(4)**2
	PINV(1)=-P(1)/PP
	PINV(2)=-P(2)/PP
	PINV(3)=-P(3)/PP
	PINV(4)=P(4)/PP
	RETURN
	END

	SUBROUTINE QPQ(P,Q,R)
C
C	QUATERNION SUM P+Q -> R
	REAL*4 P(4),Q(4),R(4)
	DO 10 I=1,4
10	R(I)=P(I)+Q(I)
	RETURN
	END

	SUBROUTINE RXQ(SCALAR,Q,R)
C
C	REAL SCALAR*Q -> R
	REAL*4 SCALAR,Q(4),R(4)
	DO 10 I=1,4
10	R(I)=SCALAR*Q(I)
	RETURN
	END

	SUBROUTINE Q2M(Q,A)
C
C	CONVERTS QUATERNION Q INTO ROTATION MATRIX A
C
	REAL*4 Q(4),A(3,3)
C
C	CONVERT TO CSC STYLE QUATERNION: I USE V(CEL)=A*V(SC)
C		BUT CSC USES V(SC)=A*V(CEL), SO USE -SPACE PARTS
C
	Q1=-Q(1)
	Q2=-Q(2)
	Q3=-Q(3)
	Q4=Q(4)
C
10	Q11=Q1*Q1
	Q22=Q2*Q2
	Q33=Q3*Q3
	Q44=Q4*Q4
	S=Q11+Q22+Q33+Q44
	IF(ABS(S-1.).GT.1.E-4) THEN
	  TYPE *,' NORMALIZATION ERROR IN NEW_Q2M'
	  S=SQRT(S)
	  Q1=Q1/S
	  Q2=Q2/S
	  Q3=Q3/S
	  Q4=Q4/S
	  GO TO 10
	ENDIF
	Q12=Q1*Q2
	Q13=Q1*Q3
	Q14=Q1*Q4
	Q23=Q2*Q3
	Q24=Q2*Q4
	Q34=Q3*Q4
C
	A(1,1) = Q11 - Q22 - Q33 + Q44
	A(1,2) = 2. * ( Q12 + Q34 )
	A(1,3) = 2. * ( Q13 - Q24 )
	A(2,1) = 2. * ( Q12 - Q34 )
	A(2,2) = -Q11 + Q22 - Q33 + Q44
	A(2,3) = 2. * ( Q23 + Q14 )
	A(3,1) = 2. * ( Q13 + Q24 )
	A(3,2) = 2. * ( Q23 - Q14 )
	A(3,3) = -Q11 - Q22 + Q33 + Q44
C
	RETURN
	END

	SUBROUTINE M2Q(A,Q)
C
C	CONVERTS ROTATION MATRIX A INTO QUATERNION Q
C
	REAL*4 A(3,3),Q(4)
C
	TRACE=A(1,1)+A(2,2)+A(3,3)
	TPMM=A(1,1)-A(2,2)-A(3,3)
	TMPM=-A(1,1)+A(2,2)-A(3,3)
	TMMP=A(3,3)-A(2,2)-A(1,1)
	IF(TMMP.GE.MAX(TRACE,TPMM,TMPM)) THEN
	  Q3=.5*SQRT(1.+TMMP)
	  FOURQ3=4.*Q3
	  Q1=(A(1,3)+A(3,1))/FOURQ3
	  Q2=(A(2,3)+A(3,2))/FOURQ3
	  Q4=(A(1,2)-A(2,1))/FOURQ3
	ELSE IF(TMPM.GT.MAX(TRACE,TPMM)) THEN
	  Q2=.5*SQRT(1.+TMPM)
	  FOURQ2=4.*Q2
	  Q1=(A(1,2)+A(2,1))/FOURQ2
	  Q3=(A(3,2)+A(2,3))/FOURQ2
	  Q4=(A(3,1)-A(1,3))/FOURQ2
	ELSE IF(TRACE.GT.TPMM) THEN
	  Q4=.5*SQRT(1.+TRACE)
	  FOURQ4=4.*Q4
	  Q1=(A(2,3)-A(3,2))/FOURQ4
	  Q2=(A(3,1)-A(1,3))/FOURQ4
	  Q3=(A(1,2)-A(2,1))/FOURQ4
	ELSE
	  Q1=.5*SQRT(1.+TPMM)
	  FOURQ1=4.*Q1
	  Q2=(A(1,2)+A(2,1))/FOURQ1
	  Q3=(A(1,3)+A(3,1))/FOURQ1
	  Q4=(A(2,3)-A(3,2))/FOURQ1
	ENDIF
C
	S=SQRT(Q1*Q1+Q2*Q2+Q3*Q3+Q4*Q4)
	IF(ABS(S-1.).GT.0.00001) TYPE *,'NORM ERROR IN NEW_M2Q'
	Q(1)=-Q1/S
	Q(2)=-Q2/S
	Q(3)=-Q3/S
	Q(4)=Q4/S
C
	RETURN
	END

	SUBROUTINE ORTH(V,W,X)
	REAL*4 V(3),W(3),X(3)
C
C	CONSTRUCT AN ORTHONORMAL TRIPLET BASED ON V
C
	B=0.
	DO 10 I=1,3
	AV=ABS(V(I))
	IF(AV.LE.B) GO TO 10
	   B=AV
	   M=I+1
10	W(I)=0.
	IF(M.EQ.4) M=1
	W(M)=1.
	CALL ORTH2(V,W,X)
	RETURN
	END

	SUBROUTINE ORTH2(V,W,X)
C
C	CONSTRUCT AN ORTHONORMAL TRIPLET WITH X-AXIS ALONG V, STORED IN V
C	Y-AXIS IN THE VW PLANE, STORED IN W
C	AND Z-AXIS FORMING RIGHT-HANDED SYSTEM, STORED IN X
C
	REAL*4 V(3),W(3),X(3)
	CALL NORM(V)
	VW=DOT(V,W)
	DO 10 I=1,3
10	W(I)=W(I)-VW*V(I)
	CALL NORM(W)
	CALL CROSS(V,W,X)
	RETURN
	END

	FUNCTION PLANCK(ANU,T)
C
C	GIVES ERGS/CM**2/SR/SEC/WAVENUMBER
C	T IN KELVIN
C	ANU IN 1/CM
C
	REAL*4 THCC,HCK
	PARAMETER (THCC=1.1910439E-5,HCK=1.438769)
	IF(T.GT.0.) GO TO 10
5	PLANCK=0.
	RETURN
10	X=HCK*ANU/T
	IF(X.GT.50.) GO TO 5
	IF(X.GT..01) GO TO 15
	D=X*(1.+X*(.5+X/6.))
	GO TO 20
15	D=EXP(X)-1.
20	PLANCK=THCC*ANU**3/D
	RETURN
	END

	FUNCTION TBRITE(ANU,PLANCK)
	REAL*4 THCC,HCK
	PARAMETER (THCC=1.1910439E-5,HCK=1.438769)
	IF(PLANCK.GT.0.) GO TO 10
5	TBRITE=0.
	RETURN
10	IF(ANU.LE.0.) GO TO 5
	D=THCC*ANU**3/PLANCK
	IF(D.GT..01) GO TO 15
	X=D*(1.-D*(.5-D/3.))
	GO TO 20
15	X=ALOG(1.+D)
20	TBRITE=HCK*ANU/X
	RETURN
	END

	SUBROUTINE EZODI(E,R,NNU,NU,INU)
	INTEGER NNU
	REAL*4 E(3),R(3),NU(NNU),INU(NNU)
C
C	zodiacal light model for arbitrary direction and point of view
C	by Ned Wright, UCLA
C
C	E is a unit vector direction of view
C	R is the viewing location, in AU.  Both E and R are in the 
C	ecliptic J2000 frame.  Using 1950 ecliptic is close enough...
C
C	NNU is the number of wavenumbers
C	NU is the list of wavenumbers for which the ZODI light should be found
C	INU is the OUTPUT list of intensites in erg/cm**2/sec/sr/cm**-1
C
C	EXTERNALS: 	DOT	dot product of two 3 vectors
C			NORM	normalized a 3 vector
C			MAV	3x3 matrix times 3 vector -> 3 vector
C			PLANCK	blackbody of NU (/cm) and T (K) in
C				erg/cm^2/sec/sr/(/cm)	(Eplee's)
C
	REAL*4 TSUN,FSCAT,FDUST(2),TDUST(2)
C
C	ZODI LIGHT PARAMETERS: DENSITY LAW
	DATA ALPHA,BETA/1.3,2.1/
C	TEMPERATURES AND EMISSIVITIES
	DATA TSUN,FSCAT/5500.,12.E-14/
C	DUST
	DATA FDUST,TDUST/1.5E-10,9.E-10,265.2,148.1/
C
	REAL*4 A(3,3),X(3),Y(3),B(3)
C
C	MATRIX TO TRANSFORM TO ZODIACAL CLOUD FRAME
	DATA A/.9975,0.,+.05,0.,1.,0.,-.05,0.,.9975/
	CALL MAV(A,E,Y)
	CALL MAV(A,R,X)
	CALL NORM(Y)
C	GET IMPACT PARAMETER
	XY=DOT(X,Y)
	DO 100 I=1,3
100	B(I)=X(I)-XY*Y(I)
	BB=SQRT(DOT(B,B))
	T1=XY/BB
	S=SQRT(1.+T1**2)
	W1=1./(S*(T1+S))
C
C	INTEGRATION VARIABLE IS w = 1 - sin(theta)
C		where X = B + b tan(theta) Y
C		Range of w is from W1 to 0.
C
	NW=4
	DW=W1/NW
	W=DW/2.
C	ZERO THE SUMS
	DO 110 I=1,NNU
110	  INU(I)=0.
C	DO THE INTEGRATION OVER LINE OF SIGHT
	DO 150 IW=1,NW
	  C=SQRT(W*(2.-W))		!cos(theta)
	  S=1.-W
	  T=BB*S/C
C	FIND THE VECTOR POSITION
	  DO 120 I=1,3
120	    X(I)=B(I)+T*Y(I)
	  RR=SQRT(DOT(X,X))
	  Z=X(3)
C	SOLAR FLUX VARIES LIKE R**-2
	  FSUN=1./RR**2
C	OVER 1/LAMBDA DUST T GOES LIKE FSUN**1/5
	  TF=FSUN**0.2
C	DENSITY LAW FROM HELIOS MODEL
	  DENSE=EXP(-BETA*ABS(Z)/RR)/RR**ALPHA
	  DS=BB*DW/C**3
	  PHASE=DOT(Y,X)/RR
C	TRY RAYLEIGH PHASE LAW ON SCATTERED LIGHT, BUT ALBEDO IS GRAY
	  SCAT=FSCAT*DENSE*(1.+PHASE**2)*FSUN*DS
C	ADD SCATTERED LIGHT TO INU
	  DO 130 I=1,NNU
130	    INU(I)=INU(I)+SCAT*PLANCK(NU(I),TSUN)
C	TWO TYPES OF DUST
	  DO 140 ID=1,2
	  THERM=FDUST(ID)*DENSE*DS
	  TD=TDUST(ID)*TF
C	ADD DUST TO INU
	    DO 140 I=1,NNU
140	      INU(I)=INU(I)+THERM*PLANCK(NU(I),TD)*NU(I)
150	  W=W+DW
	RETURN
	END

	SUBROUTINE SUNPOS(T,SCEL)
C
C	NEW VERSION OF SUNPOS THAT USES EPHEM ROUTINES
C	30-AUG-83
C
C	INPUT T REAL*8 SECONDS SINCE 5/24/68
C
C	OUTPUT SCEL, REAL*4 SCEL(3), UNIT VECTOR POINTING AT SUN
C
	REAL*8 T,T1,T2
	REAL*4 SCEL(3),S1(3),S2(3)
	DATA T1,T2/1.D20,-1.D20/
	IF(T1.LE.T .AND. T.LE.T2) GO TO 100
	   T1=T-100.D0
	   CALL EPHEM(-1,T1,S1)
	   CALL NORM(S1)
	   T2=T1+3800.D0
	   CALL EPHEM(-1,T2,S2)
	   CALL NORM(S2)
100	X=(T-T1)/(T2-T1)
	Y=1.-X
	DO 150 I=1,3
150	SCEL(I)=X*S2(I)+Y*S1(I)
	RETURN
	END

	SUBROUTINE GCMOON(T,PMOON)
C
C	NEW VERSION USING EPHEM, 30-AUG-83
C
C	INPUT T REAL*8 SECONDS SINCE 5/24/68
C
C	OUTPUT REAL*4 PMOON(3) GEOCENTRIC VECTOR POS'N OF MOON
C	IN KILOMETERS
C
	PARAMETER (AEARTH=6378.16)	! equatorial radius of Earth
	REAL*8 T
	REAL*4 PMOON(3)
	CALL EPHEM(0,T,PMOON)
	DO 100 I=1,3
100	PMOON(I)=AEARTH*PMOON(I)
	RETURN
	END

	SUBROUTINE EPHEM(N,T,P)
C
C	EPHERMIS ROUTINE FOR COBE
C	INPUT N INTEGER*4 PLANETS NUMBER
C		-1 FOR SUN
C		0 FOR MOON
C		1 FOR MARS
C		2 FOR JUPITER
C		3 FOR SATURN
C		4 FOR URANUS
C		5 FOR NEPTUNE
C		6 FOR 1st asteroid IN FILE 'ASTORBITS', 7 FOR 2nd etc.
C
C	INPUT T REAL*8 TIME SINCE 5/24/68 0 HOURS UTC
C
C	OUTPUT VECTOR P, REAL*4 P(3)
C		DIRECTION OF P IS GEOCENTRIC POSITION OF PLANET
C		IN CELESTIAL COORDINATE SYSTEM WITH EPOCH 2000
C		MAGNITUDE IS DISTANCE IN AU EXCEPT FOR MOON IN EARTH RADII
C
C	Based on Van Flandern & Pulkinnen, 1979, Ap.J.Supp. v. 41, pp. 391-411
C	Only the heliocentric series for PLON, BETA and RP are used.
C	The U, V, W series have been deleted for compactness.
C
	INTEGER*4 N
	REAL*4 P(3),PE(3),SE(3),CONST(33),DBAR(10),STR,PRECES(3,3)
	REAL*8 T,TT,TTT,BETA,RP,PLON,ARGS(33),U,V,W,SLOPES(33)
	REAL*8 JD,JULIAN_DAY,CYCLES,TOLD
	DATA TOLD/-1.D20/
	COMMON /FMCMMN/ PLON,U,V,W,ARGS
	PARAMETER (PI=3.1415927,TWOPI=6.2831853,DTR=57.29578)
	DATA CONST/.606434,.374897,.259091,.827362,.347343,0.,
	1   .779072,.993126,.700695,.485541,.566441,.505498,
	2   .140023,.292498,.987353,.053856,.849694,.089608,
	3   .056531,.814794,.133295,.882987,.821218,.870169,
	4   .400589,.664614,.846912,.725368,.480856,0.,
	5   .663854,.041020,.357335/
	DATA SLOPES/.3660110129D-1,.3629164709D-1,.367481952D-1,
	1   .3386319198D-1,-.14709391D-3,0.D0,		!Moon
	2   .273790931D-2,.27377785D-2,			!Sun
	3   .11367714D-1,.1136759566D-1,.1136762384D-1,	!Mercury
	4   .445046867D-2,.445036173D-2,.445040017D-2,	!Venus
	5   .145575328D-2,.145561327D-2,.145569465D-2,	!Mars
	6   3*.23080893D-3,				!Jupiter???
	7   3*.9294371D-4,				!Saturn???
	8   2*.3269438D-4,.3265562D-4,			!Uranus
	9   2*.1672092D-4,.1663715D-4,	!Neptune fixed 9/21/92
	9   0.D0,.1115482D-4,2*.1104864D-4/		!Pluto
	DATA DBAR/1.00021,60.40974,1.83094,	!Sun, Moon, Mars
	1   5.30693,9.61711,19.24877,30.08900,	!Jupiter,Saturn,Uranus,Neptune
	2   41.3268,1.07693,1.23437/		!Pluto,Mercury,Venus
C
C
C	COMPUTE ARGUMENTS
C
	IF(DABS(T-TOLD).LE.1.E-2) GO TO 150
	STR=3600.*DTR
	CCC=86400.*2.9975E10/1.49E13		!Speed of light in AU/day
	JD=JULIAN_DAY(T)
	TT=JD-2451545.D0			!Days since J2000.0
	TTT=TT/36525.D0+1.D0			!Centuries since 1900.0
	EPOCH=1900.0+100.0*TTT
	CALL PREMAT(EPOCH,PRECES)
	DO 100 I=1,33
	CYCLES=CONST(I)+SLOPES(I)*TT
100	ARGS(I)=TWOPI*DMOD(CYCLES,1.D0)
	TOLD=T
C	
C	CALL APPROPRIATE ROUTINE
C
150	DO 152 I=1,3
	   PE(I)=0.				
152	   SE(I)=0.
	IF(N.NE.0) THEN
	   CALL SUN(BETASUN,RP,TTT)		
	   SUNLON=ARGS(7)+PLON/STR		
	   SE(1)=COS(SUNLON)*RP
	   SE(2)=SIN(SUNLON)*RP
	   ENDIF
	IF(N.EQ.-1) GO TO 200			
	IF(N.EQ.0) THEN
	   CALL MOONVF(BETA,RP,TTT)
	   NL=1
	   ENDIF
	IF(N.EQ.1) THEN
	   CALL MARS(BETA,RP,TTT)
	   NL=15
	   ENDIF
	IF(N.EQ.2) THEN
	   CALL JUPITER(BETA,RP,TTT)
	   NL=18
	   ENDIF
	IF(N.EQ.3) THEN
	   CALL SATURN(BETA,RP,TTT)
	   NL=21
	   ENDIF
	IF(N.EQ.4) THEN
	   CALL URANUS(BETA,RP,TTT)
	   NL=24
	   ENDIF
	IF(N.EQ.5) THEN
	   CALL NEPTUNE(BETA,RP,TTT)
	   NL=27
	   ENDIF
	IF(N.GE.6) THEN
	   DJD=2./CCC
	   DO ITER=1,2
	      CALL ASTEROID_ORBIT(N-5,JD-DJD,PE)
	      DO I=1,3
	         PE(I)=PE(I)+SE(I)
	      ENDDO
	      DJD=SQRT(DOT(PE,PE))/CCC
	   ENDDO
	   GO TO 250
	ENDIF
	IF(N.LT.-1) STOP 'N TOO LOW IN EPHEM'
C
	PLONG=ARGS(NL)+PLON/STR
	BRAD=BETA/STR
C
C	Compute heliocentric vector in ecliptic coordinates
C
	PE(3)=RP*SIN(BRAD)			
	PE(2)=RP*COS(BRAD)			
	PE(1)=PE(2)*COS(PLONG)
	PE(2)=PE(2)*SIN(PLONG)
200	DO 220 I=1,3
220	   PE(I)=PE(I)+SE(I)			!add in Sun to get geocentric
250	EPS=(84428.-47.*TTT+9.*COS(ARGS(5)))/STR
	CEPS=COS(EPS)
	SEPS=SIN(EPS)
	DLAM=-17.*SIN(ARGS(5))/STR	!accurate conversion to celestial
	P(1)=PE(1)-DLAM*PE(2)
	PE(2)=PE(2)+DLAM*PE(1)
	P(2)=PE(2)*CEPS-PE(3)*SEPS
	P(3)=PE(3)*CEPS+PE(2)*SEPS
C
C	CONVERT TO 2000
C
	CALL VMAT(PRECES,P,PE)
	DO I=1,3
	  P(I)=P(I)+PE(I)
	ENDDO
	RETURN
	END

	SUBROUTINE SUN (BETA,RP,T)
C
C	SUN CODE
C
C
	IMPLICIT REAL*8 (A-H,L,O-Z)
	IMPLICIT INTEGER*4 (I-K,M-N)
	COMMON /FMCMMN/ PLON,U,V,W,A1,A2,A3,A4,A5,A6,A7,A8,A9,
	1   A10,A11,A12,A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23,
	1   A24,A25,A26,A27,A28,A29,A30,A31,A32,A33
C
C
	PLON=6910*DSIN(A8)+72*DSIN(2*A8)-17*T*DSIN(A8)-7*DCOS(A8-A19)
	1   +6*DSIN(A1-A7)+5*DSIN(4*A8-8*A16+3*A19)-5*DCOS(2*A8-2*A13)
	2   -4*DSIN(A8-A13)+4*DCOS(4*A8-8*A16+3*A19)+3*DSIN(2*A8-2*A13)
	3   -3*DSIN(A19)-3*DSIN(2*A8-2*A19)
C	DONE
	BETA=0.0
C	DONE
	RP=1.00014D0-0.01675*DCOS(A8)-0.00014*DCOS(2*A8)
C	DONE
	RETURN
	END

	SUBROUTINE MOONVF (BETA,RP,T)
C
C	MOON CODE
C
	IMPLICIT REAL*8 (A-H,L,O-Z)
	IMPLICIT INTEGER*4 (I-K,M-N)
	COMMON /FMCMMN/ PLON,U,V,W,A1,A2,A3,A4,A5,A6,A7,A8,A9,
	1   A10,A11,A12,A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23,
	1   A24,A25,A26,A27,A28,A29,A30,A31,A32,A33
C
	PL1=22640.D0*DSIN(A2)-4586.D0*DSIN(A2-2*A4)+2370*DSIN(2*A4)
	1   +769*DSIN(2*A2)-668*DSIN(A8)-412*DSIN(2*A3)-212*DSIN(2*A2-2*A4)
	2   -206*DSIN(A2-2*A4+A8)+192*DSIN(A2+2*A4)+165*DSIN(2*A4-A8)
	3   +148*DSIN(A2-A8)-125*DSIN(A4)-110*DSIN(A2+A8)-55*DSIN(2*A3-2*A4)
	4   -45*DSIN(A2+2*A3)+40*DSIN(A2-2*A3)-38*DSIN(A2-4*A4)
	5   +36*DSIN(3*A2)-31*DSIN(2*A2-4*A4)+28*DSIN(A2-2*A4-A8)
	6   -24*DSIN(2*A4+A8)+19*DSIN(A2-A4)+18*DSIN(A4+A8)
	7   +15*DSIN(A2+2*A4-A8)+14*DSIN(2*A2+2*A4)+14*DSIN(4*A4)
	PL2=-13*DSIN(3*A2-2*A4)-11*DSIN(A2+16*A7-18*A12)+10*DSIN(2*A2-A8)
	9   +9*DSIN(A2-2*A3-2*A4)+9*DCOS(A2+16*A7-18*A12)
	9   -9*DSIN(2*A2-2*A4+A8)-8*DSIN(A2+A4)+8*DSIN(2*A4-2*A8)
	1   -8*DSIN(2*A2+A8)-7*DSIN(2*A8)-7*DSIN(A2-2*A4+2*A8)
	2   +7*DSIN(A5)-6*DSIN(A2-2*A3+2*A4)-6*DSIN(2*A3+2*A4)
	3    -4*DSIN(A2-4*A4+A8)+4*T*DCOS(A2+16*A7-18*A12)-4*DSIN(2*A2+2*A3)
	PL3=4*T*DSIN(A2+16*A7-18*A12)+3*DSIN(A2-3*A4)-3*DSIN(A2+2*A4+A8)
	1   -3*DSIN(2*A2-4*A4+A8)+3*DSIN(A2-2*A8)+3*DSIN(A2-2*A4-2*A8)
	2   -2*DSIN(2*A2-2*A4-A8)-2*DSIN(2*A3-2*A4+A8)+2*DSIN(A2+4*A4)
	3   +2*DSIN(4*A2)+2*DSIN(4*A4-A8)+2*DSIN(2*A2-A4)
	PLON=PL1+PL2+PL3
C	DONE
	BETA=18461.D0*DSIN(A3)+1010*DSIN(A2+A3)+1000*DSIN(A2-A3)
	1   -624*DSIN(A3-2*A4)-199*DSIN(A2-A3-2*A4)-167*DSIN(A2+A3-2*A4)
	2   +117*DSIN(A3+2*A4)+62*DSIN(2*A2+A3)+33*DSIN(A2-A3+2*A4)
	3   +32*DSIN(2*A2-A3)-30*DSIN(A3-2*A4+A8)-16*DSIN(2*A2+A3-2*A4)
	4   +15*DSIN(A2+A3+2*A4)+12*DSIN(A3-2*A4-A8)-9*DSIN(A2-A3-2*A4+A8)
	5   -8*DSIN(A3+A5)+8*DSIN(A3+2*A4-A8)-7*DSIN(A2+A3-2*A4+A8)
	6   +7*DSIN(A2+A3-A8)-7*DSIN(A2+A3-4*A4)-6*DSIN(A3+A8)
	7   -6*DSIN(3*A3)+6*DSIN(A2-A3-A8)-5*DSIN(A3+A4)-5*DSIN(A2+A3+A8)
	8   -5*DSIN(A2-A3+A8)+5*DSIN(A3-A8)+5*DSIN(A3-A4)+4*DSIN(3*A2+A3)
	9   -4*DSIN(A3-4*A4)-3*DSIN(A2-A3-4*A4)+3*DSIN(A2-3*A3)
	9   -2*DSIN(2*A2-A3-4*A4)-2*DSIN(3*A3-2*A4)+2*DSIN(2*A2-A3+2*A4)
	1   +2*DSIN(A2-A3+2*A4-A8)+2*DSIN(2*A2-A3-2*A4)+2*DSIN(3*A2-A3)
C	DONE
	RP=60.36298D0-3.27746D0*DCOS(A2)-0.57994D0*DCOS(A2-2*A4)
	1   -0.46357*DCOS(2*A4)-0.08904*DCOS(2*A2)+0.03865*DCOS(2*A2-2*A4)
	2   -0.03237*DCOS(2*A4-A8)-0.02688*DCOS(A2+2*A4)
	3   -0.02358*DCOS(A2-2*A4+A8)-0.02030*DCOS(A2-A8)+0.01719*DCOS(A4)
	4   +0.01671*DCOS(A2+A8)+0.01247*DCOS(A2-2*A3)+0.00704*DCOS(A8)
	5   +0.00529*DCOS(2*A4+A8)-0.00524*DCOS(A2-4*A4)
	6   +0.00398*DCOS(A2-2*A4-A8)-0.00366*DCOS(3*A2)
	7   -0.00295*DCOS(2*A2-4*A4)-0.00263*DCOS(A4+A8)
	8   +0.00249*DCOS(3*A2-2*A4)-0.00221*DCOS(A2+2*A4-A8)
	9   +0.00185*DCOS(2*A3-2*A4)-0.00161*DCOS(2*A4-2*A8)
	9   +0.00147*DCOS(A2+2*A3-2*A4)-0.00142*DCOS(4*A4)
	1   +0.00139*DCOS(2*A2-2*A4+A8)-0.00118*DCOS(A2-4*A4+A8)
	2   -0.00116*DCOS(2*A2+2*A4)-0.00110*DCOS(2*A2-A8)
C	DONE
	RETURN
	END

	SUBROUTINE MARS (BETA,RP,T)
C
C	MARS CODE
C
	IMPLICIT REAL*8 (A-H,L,O-Z)
	IMPLICIT INTEGER*4 (I-K,M-N)
	COMMON /FMCMMN/ PLON,U,V,W,A1,A2,A3,A4,A5,A6,A7,A8,A9,
	1   A10,A11,A12,A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23,
	1   A24,A25,A26,A27,A28,A29,A30,A31,A32,A33
C
C
	PLON=38451.D0*DSIN(A16)+2238.D0*DSIN(2*A16)+181.D0*DSIN(3*A16)
	1   -52*DSIN(2*A17)+37*T*DSIN(A16)-22*DCOS(A16-2*A19)
	2   -19*DSIN(A16-A19)+17*DCOS(A16-A19)+17*DSIN(4*A16)
	3   -16*DCOS(2*A16-2*A19)+13*DCOS(A8-2*A16)-10*DSIN(A16-2*A17)
	4   -10*DSIN(A16+2*A17)+7*DCOS(A8-A16)-7*DCOS(2*A8-3*A16)
	5   -5*DSIN(A13-3*A16)-5*DSIN(A8-A16)-5*DSIN(A8-2*A16)
	6   -4*DCOS(2*A8-4*A16)+4*T*DSIN(2*A16)+4*DCOS(A19)
	7   +3*DCOS(A13-3*A16)+3*DSIN(2*A16-2*A19)
C	DONE
	BETA=6603.D0*DSIN(A17)+622*DSIN(A16-A17)+615*DSIN(A16+A17)
	1   +64*DSIN(2*A16+A17)
C	DONE
	RP=1.53031D0-0.1417D0*DCOS(A16)-0.0066*DCOS(2*A16)
	1   -0.00047*DCOS(3*A16)
C	DONE
	RETURN
	END

	SUBROUTINE JUPITER (BETA,RP,T)
C
C	JUPITER CODE
C
	IMPLICIT REAL*8 (A-H,L,O-Z)
	IMPLICIT INTEGER*4 (I-K,M-N)
	COMMON /FMCMMN/ PLON,U,V,W,A1,A2,A3,A4,A5,A6,A7,A8,A9,
	1   A10,A11,A12,A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23,
	1   A24,A25,A26,A27,A28,A29,A30,A31,A32,A33
C
C
	PLON1=19934.D0*DSIN(A19)+5023.D0*T+2511+1093.D0*DCOS(2*A19-5*A22)
	1   +601*DSIN(2*A19)-479*DSIN(2*A19-5*A22)-185*DSIN(2*A19-2*A22)
	2   +137*DSIN(3*A19-5*A22)-131*DSIN(A19-2*A22)+79*DCOS(A19-A22)
	3   -76*DCOS(2*A19-2*A22)-74*T*DCOS(A19)+68*T*DSIN(A19)
	4   +66*DCOS(2*A19-3*A22)+63*DCOS(3*A19-5*A22)+53*DCOS(A19-5*A22)
	5   +49*DSIN(2*A19-3*A22)-43*T*DSIN(2*A19-5*A22)-37*DCOS(A19)
	6   +25*DSIN(2*A18)+25*DSIN(3*A19)-23*DSIN(A19-5*A22)
	7   -19*T*DCOS(2*A19-5*A22)+17*DCOS(2*A19-4*A22)
	8   +17*DCOS(3*A19-3*A22)-14*DSIN(A19-A22)-13*DSIN(3*A19-4*A22)
	9   -9*DCOS(2*A18)+9*DCOS(A22)-9*DSIN(A22)-9*DSIN(3*A19-2*A22)
	9   +9*DSIN(4*A19-5*A22)+9*DSIN(2*A19-6*A22+3*A25)-8*DCOS(4*A19-10*A22
	1   )
	2   +7*DCOS(3*A19-4*A22)-7*DCOS(A19-3*A22)-7*DSIN(4*A19-10*A22)
	3   -7*DSIN(A19-3*A22)+6*DCOS(4*A19-5*A22)-6*DSIN(3*A19-3*A22)
	PLON2=+5*DCOS(2*A22)-4*DSIN(4*A19-4*A22)-4*DCOS(3*A22)
	1   +4*DCOS(2*A19-A22)-4*DCOS(3*A19-2*A22)-4*T*DCOS(2*A19)
	2   +3*T*DSIN(2*A19)+3*DCOS(5*A22)+3*DCOS(5*A19-10*A22)
	3   +3*DSIN(2*A22)-2*DSIN(2*A18-A19)+2*DSIN(2*A18+A19)
	4   -2*T*DSIN(3*A19-5*A22)-2*T*DSIN(A19-5*A22)
	PLON=PLON1+PLON2
C	DONE
	BETA=-4692*DCOS(A19)+259*DSIN(A19)+227-227*DCOS(2*A19)
	1   +30*T*DSIN(A19)+21*T*DCOS(A19)+16*DSIN(3*A19-5*A22)
	2   -13*DSIN(A19-5*A22)-12*DCOS(3*A19)+12*DSIN(2*A19)
	3   +7*DCOS(3*A19-5*A22)-5*DCOS(A19-5*A22)
C	DONE
	RP=5.20883D0-0.25122D0*DCOS(A19)-0.00604*DCOS(2*A19)
	1   +0.00260*DCOS(2*A19-2*A22)-0.00170*DCOS(3*A19-5*A22)
	2   -0.00106*DSIN(2*A19-2*A22)-0.00091*T*DSIN(A19)
	3   -0.00084*T*DCOS(A19)+0.00069*DSIN(2*A19-3*A22)
	4   -0.00067*DSIN(A19-5*A22)+0.00066*DSIN(3*A19-5*A22)
	5   +0.00063*DSIN(A19-A22)-0.00051*DCOS(2*A19-3*A22)
	6   -0.00046*DSIN(A19)-0.00029*DCOS(A19-5*A22)
	7   +0.00027*DCOS(A19-2*A22)-0.00022*DCOS(3*A19)
	8   -0.00021*DSIN(2*A19-5*A22)
C	DONE
	RETURN
	END

	SUBROUTINE SATURN (BETA,RP,T)
C
C	SATURN CODE
C
	IMPLICIT REAL*8 (A-H,L,O-Z)
	IMPLICIT INTEGER*4 (I-K,M-N)
	COMMON /FMCMMN/ PLON,U,V,W,A1,A2,A3,A4,A5,A6,A7,A8,A9,
	1   A10,A11,A12,A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23,
	1   A24,A25,A26,A27,A28,A29,A30,A31,A32,A33
C
C
	PLON1=23045.D0*DSIN(A22)+5014.D0*T-2689.D0*DCOS(2*A19-5*A22)+2507
	1   +1177*DSIN(2*A19-5*A22)-826*DCOS(2*A19-4*A22)
	2   +802*DSIN(2*A22)+425*DSIN(A19-2*A22)-229*T*DCOS(A22)
	3   -153*DCOS(2*A19-6*A22)-142*T*DSIN(A22)-114*DCOS(A22)
	4   +101*T*DSIN(2*A19-5*A22)-70*DCOS(2*A21)+67*DSIN(2*A21)
	5   +66*DSIN(2*A19-6*A22)+60*T*DCOS(2*A19-5*A22)
	6   +41*DSIN(A19-3*A22)+39*DSIN(3*A22)+31*DSIN(A19-A22)
	7   +31*DSIN(2*A19-2*A22)-29*DCOS(2*A19-3*A22)
	8   -28*DSIN(2*A19-6*A22+3*A25)+28*DCOS(A19-3*A22)
	9   +22*T*DSIN(2*A19-4*A22)-22*DSIN(A22-3*A25)+20*DSIN(2*A19-3*A22)
	9   +20*DCOS(4*A19-10*A22)+19*DCOS(2*A22-3*A25)
	1   +19*DSIN(4*A19-10*A22)-17*T*DCOS(2*A22)-16*DCOS(A22-3*A25)
	2   -12*DSIN(2*A19-4*A22)+12*DCOS(A19)-12*DSIN(2*A22-2*A25)
	3   -11*T*DSIN(2*A22)-11*DCOS(2*A19-7*A22)+10*DSIN(2*A22-3*A25)
	4   +10*DCOS(2*A19-2*A22)+9*DSIN(4*A19-9*A22)
	PLON2=-8*DSIN(A22-2*A25)-8*DCOS(2*A21+A22)+8*DCOS(2*A21-A22)
	1   +8*DCOS(A22-A25)-8*DSIN(2*A21-A22)+7*DSIN(2*A21+A22)
	2   -7*DCOS(A19-2*A22)-7*DCOS(2*A22)-6*T*DSIN(4*A19-10*A22)
	3   +6*T*DCOS(4*A19-10*A22)+6*T*DSIN(2*A19-6*A22)
	4   -5*DSIN(3*A19-7*A22)-5*DCOS(3*A19-3*A22)-5*DCOS(2*A22-2*A25)
	5   +5*DSIN(3*A19-4*A22)+5*DSIN(2*A19-7*A22)+4*DSIN(3*A19-3*A22)
	6   +4*DSIN(3*A19-5*A22)+4*T*DCOS(A19-2*A22)+3*T*DCOS(2*A19-4*A22)
	7   +3*DCOS(2*A19-6*A22+3*A25)-3*T*DSIN(2*A21)+3*T*DCOS(2*A19-6*A22)
	8   -3*T*DCOS(2*A21)+3*DCOS(3*A19-7*A22)
	PLON3=+3*DCOS(4*A19-9*A22)+3*DSIN(3*A19-6*A22)+3*DSIN(2*A19-A22)
	1   +3*DSIN(A19-4*A22)+2*DCOS(3*A22-3*A25)+2*T*DSIN(A19-2*A22)
	2   +2*DSIN(4*A22)-2*DCOS(3*A19-4*A22)-2*DCOS(2*A19-A22)
	3   -2*DSIN(2*A19-7*A22+3*A25)+2*DCOS(A19-4*A22)+2*DCOS(4*A19-11*A22)
	4   -2*DSIN(A22-A25)
	PLON=PLON1+PLON2+PLON3
C	DONE
	BETA=8297*DSIN(A22)-3346*DCOS(A22)+462*DSIN(2*A22)
	1   -189*DCOS(2*A22)+185+79*T*DCOS(A22)-71*DCOS(2*A19-4*A22)
	2   +46*DSIN(2*A19-6*A22)-45*DCOS(2*A19-6*A22)+29*DSIN(3*A22)
	3   -20*DCOS(2*A19-3*A22)+18*T*DSIN(A22)-14*DCOS(2*A19-5*A22)
	4   -11*DCOS(3*A22)-10*T+9*DSIN(A19-3*A22)+8*DSIN(A19-A22)
	5   -6*DSIN(2*A19-3*A22)+5*DSIN(2*A19-7*A22)-5*DCOS(2*A19-7*A22)
	6   +4*DSIN(2*A19-5*A22)-4*T*DSIN(2*A22)-3*DCOS(A19-A22)
	7   +3*DCOS(A19-3*A22)+3*T*DSIN(2*A19-4*A22)+3*DSIN(A19-2*A22)
	8   +2*DSIN(4*A22)-2*DCOS(2*A19-2*A22)
C	DONE
	RP=9.55774D0-0.53252D0*DCOS(A22)-0.01878D0*DSIN(2*A19-4*A22)
	1   -0.01482*DCOS(2*A22)+0.00817*DSIN(A19-A22)
	2   -0.00539*DCOS(A19-2*A22)-0.00524*T*DSIN(A22)
	3   +0.00349*DSIN(2*A19-5*A22)+0.00347*DSIN(2*A19-6*A22)
	4   +0.00328*T*DCOS(A22)-0.00225*DSIN(A22)
	5   +0.00149*DCOS(2*A19-6*A22)-0.00126*DCOS(2*A19-2*A22)
	6   +0.00104*DCOS(A19-A22)+0.00101*DCOS(2*A19-5*A22)
	7   +0.00098*DCOS(A19-3*A22)-0.00073*DCOS(2*A19-3*A22)
	8   -0.00062*DCOS(3*A22)+0.00042*DSIN(2*A22-3*A25)
	9   +0.00041*DSIN(2*A19-2*A22)-0.00040*DSIN(A19-3*A22)
	9   +0.00040*DCOS(2*A19-4*A22)-0.00028*T-0.00023*DSIN(A19)
	1   +0.00020*DSIN(2*A19-7*A22)
C	DONE
	RETURN
	END

	SUBROUTINE URANUS (BETA,RP,T)
C
C	URANUS CODE
C
	IMPLICIT REAL*8 (A-H,L,O-Z)
	IMPLICIT INTEGER*4 (I-K,M-N)
	COMMON /FMCMMN/ PLON,U,V,W,A1,A2,A3,A4,A5,A6,A7,A8,A9,
	1   A10,A11,A12,A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23,
	1   A24,A25,A26,A27,A28,A29,A30,A31,A32,A33
C
C
	PLON=19397.D0*DSIN(A25)+570*DSIN(2*A25)-536*T*DCOS(A25)
	1   +143*DSIN(A22-2*A25)+110*T*DSIN(A25)+102*DSIN(A22-3*A25)
	2   +76*DCOS(A22-3*A25)-49*DSIN(A19-A25)+32*T*T-30*T*DCOS(2*A25)
	3   +29*DSIN(2*A19-6*A22+3*A25)+29*DCOS(2*A25-2*A28)
	4   -28*DCOS(A25-A28)+23*DSIN(3*A25)-21*DCOS(A19-A25)
	5   +20*DSIN(A25-A28)+20*DCOS(A22-2*A25)-19*DCOS(A22-A25)
	6   +17*DSIN(2*A25-3*A28)+14*DSIN(3*A25-3*A28)+13*DSIN(A22-A25)
	7   -12*T*T*DCOS(A25)-12*DCOS(A25)+10*DSIN(2*A25-2*A28)
	8   -9*DSIN(2*A26)-9*T*T*DSIN(A25)+9*DCOS(2*A25-3*A28)
	9   +8*T*DCOS(A22-2*A25)+7*T*DCOS(A22-3*A25)-7*T*DSIN(A22-3*A25)
	9   +7*T*DSIN(2*A25)+6*DSIN(2*A19-6*A22+2*A25)
	1   +6*DCOS(2*A19-6*A22+2*A25)+5*DSIN(A22-4*A25)-4*DSIN(3*A25-4*A28)
	2   +4*DCOS(3*A25-3*A28)-3*DCOS(A28)-2*DSIN(A28)
C	DONE
	BETA=2775*DSIN(A26)+131*DSIN(A25-A26)+130*DSIN(A25+A26)
C	DONE
	RP=19.21216D0-0.90154D0*DCOS(A25)-0.02488*T*DSIN(A25)
	1   -0.02121*DCOS(2*A25)-0.00585*DCOS(A22-2*A25)
	2   -0.00508*T*DCOS(A25)-0.00451*DCOS(A19-A25)
	3   +0.00336*DSIN(A22-A25)+0.00198*DSIN(A19-A25)
	4   +0.00118*DCOS(A22-3*A25)+0.00107*DSIN(A22-2*A25)
	5   -0.00103*T*DSIN(2*A25)-0.00081*DCOS(3*A25-3*A28)
C	DONE
	RETURN
	END

	SUBROUTINE NEPTUNE (BETA,RP,T)
C
C	NEPTUNE CODE
C
	IMPLICIT REAL*8 (A-H,L,O-Z)
	IMPLICIT INTEGER*4 (I-K,M-N)
	COMMON /FMCMMN/ PLON,U,V,W,A1,A2,A3,A4,A5,A6,A7,A8,A9,
	1   A10,A11,A12,A13,A14,A15,A16,A17,A18,A19,A20,A21,A22,A23,
	1   A24,A25,A26,A27,A28,A29,A30,A31,A32,A33
C
C
	PLON=3523*DSIN(A28)-50*DSIN(2*A29)-43*T*DCOS(A28)
	1   +29*DSIN(A19-A28)+19*DSIN(2*A28)-18*DCOS(A19-A28)
	2   +13*DCOS(A22-A28)+13*DSIN(A22-A28)-9*DSIN(2*A25-3*A28)
	3   +9*DCOS(2*A25-2*A28)-5*DCOS(2*A25-3*A28)-4*T*DSIN(A28)
	4   +4*DCOS(A25-2*A28)+4*T*T*DSIN(A28)
C	DONE
	BETA=6404*DSIN(A29)+55*DSIN(A28+A29)+55*DSIN(A28-A29)
	1   -33*T*DSIN(A29)
C	DONE
	RP=30.07175D0-0.25701D0*DCOS(A28)-0.00787*DCOS(2*A24-A25-2*A27)
	1   +0.00409*DCOS(A19-A28)-0.00314*T*DSIN(A28)+0.0025*DSIN(A19-A28)
	2   -0.00194*DSIN(A22-A28)+0.00185*DCOS(A22-A28)
C	DONE
	RETURN
	END

	SUBROUTINE ASTEROID_ORBIT(NA,JDIN,RR)
	INTEGER NA
	REAL*8 JDIN
	REAL*4 RR(3)
C
C	CONVERTS ELLIPTICAL ELEMENTS IN ORB INTO RR AT JDIN
C	NA: asteroid number in 'ASTORBITS' file
C	JDIN: Julian date for position
C	RR: heliocentric current ecliptic vector position in AU
C
	PARAMETER (NAST=26)
	PARAMETER (PI=3.1415927,TWOPI=6.2831853,DTR=57.29578)
	COMMON /ASTEROID_ORBITS/ JD(2),ND,ORB(10,2,NAST)
	REAL*8 JD
	INTEGER*4 ND
	REAL*4 ORB
C
	REAL*4 R(3,2)
	DO ID=1,ND
	  A=ORB(2,ID,NA)
	  T=JDIN-JD(1)
	  AN=(T-ORB(1,ID,NA))*ORB(10,ID,NA)
	  AN = MOD(AN,TWOPI)
	  E=ORB(3,ID,NA)
	  NIT=0
	  EE=AN+E*SIN(AN)
10	    ERR=AN-EE+E*SIN(EE)
  	    IF(ABS(ERR).LT.1.E-6) GO TO 20
	    D=E*COS(EE)-1.
	    DE=-ERR/D
	    EE=EE+DE
	    NIT=NIT+1
	    IF(NIT.LT.10) GO TO 10
	    TYPE 900,ERR
900	FORMAT(' NOT CONVERGED IN ORBRV, ERR=',F12.7)
20	  F=1.-E*E
	  CE=COS(EE)
	  SE=SIN(EE)
	  CP=(CE-E)/(1.-E*CE)
	  SP=SE*(1.+E*CP)/SQRT(F)
	  RHO=A*F/(1.+E*CP)
	  DO 30 I=1,3
30	    R(I,ID)=RHO*(CP*ORB(I+3,ID,NA)+SP*ORB(I+6,ID,NA))
	ENDDO
	IF(ND.EQ.1) THEN
	  DO I=1,3
	    RR(I)=R(I,1)
	  ENDDO
	  GO TO 100
	ENDIF
	X=(JDIN-JD(1))/(JD(2)-JD(1))
	IF(X.LE.0.) THEN
	  V=1.
	ELSE IF (X.GE.1.) THEN
	  V=0.
	ELSE
	  V=1.-3.*X**2+2.*X**3
	ENDIF
	W=1.-V
	DO I=1,3
	  RR(I)=V*R(I,1)+W*R(I,2)
	ENDDO
	EPS=SQRT(DVS(R(1,1),R(1,2))/DOT(R,R))
	TYPE 9666,JDIN,NA,EPS
9666	FORMAT(F12.1,I8,F10.6)
C
C	CONVERT FROM J2000.0 TO CURRENT EPOCH, NEGLECTING THE "b" TERM
C	WHICH IS  < 7" FOR 1985-2015
C
100	T=(JDIN-2451545.D0)/36525.D0
	A=T*(1.396971+0.0003086*T)/DTR
	SA=SIN(A)
	CA=COS(A)
	T=RR(1)*CA-RR(2)*SA
	RR(2)=RR(2)*CA+RR(1)*SA
	RR(1)=T
	RETURN
	END

	SUBROUTINE READ_ORBITS
C
C	READS ASTEROID INFORMATION INTO COMMON BLOCK
C
C	ASTEROID COMMON BLOCK
C
	PARAMETER (NAST=26)
	COMMON /ASTEROID_ORBITS/ JD(2),ND,ORB(10,2,NAST)
C
	REAL*8 JD
	INTEGER*4 ND
	REAL*4 ORB
C
	REAL*4 INC,NODE,MM,ME
	OPEN (UNIT=18,FILE='ASTORBITS',READONLY,STATUS='OLD')
	READ (18,900) JD
900	FORMAT(2D10.0)
901	FORMAT(10X,7F10.4)
	ND=2
	IF(JD(2).LE.0.) ND=1
	DO N=1,NAST
	  DO I=1,ND
	    READ(18,901,END=100) INC,NODE,ARG,A,MM,E,ME
	    TOSC=JD(I)-JD(1)
	    CALL ORBIT_CONVERT(ORB(1,I,N),TOSC,INC,NODE,ARG,A,MM,E,ME)
	  ENDDO
	ENDDO
100	CLOSE(UNIT=18)
	RETURN
	END

	SUBROUTINE ORBIT_CONVERT(ORB,TOSC,INC,NODE,ARG,A,RATE,E,ME)
C	ROUTINE TO CONVERT ORBIT IN ECLIPTIC COORDINATES TO INTERNAL FORM
	DIMENSION ORB(10)
	REAL INC,NODE,ME
	PARAMETER (PI=3.1415927,TWOPI=6.2831853,DTR=57.29578)
	ORB(2)=A
	ORB(3)=E
	OMEGA=RATE/DTR
	ORB(10)=OMEGA
	DT=ME/DTR/OMEGA
	ORB(1)=TOSC-DT
	DO 10 I=4,9
10	ORB(I)=0.
	ORB(4)=1.
	ORB(8)=1.
	SA=SIN(ARG/DTR)
	CA=COS(ARG/DTR)
	DO 20 I=4,7,3
	T=CA*ORB(I)-SA*ORB(I+1)
	ORB(I+1)=CA*ORB(I+1)+SA*ORB(I)
20	ORB(I)=T
	SI=SIN(INC/DTR)
	CI=COS(INC/DTR)
	DO 30 I=5,8,3
	T=CI*ORB(I)-SI*ORB(I+1)
	ORB(I+1)=CI*ORB(I+1)+SI*ORB(I)
30	ORB(I)=T
	SN=SIN(NODE/DTR)
	CN=COS(NODE/DTR)
	DO 40 I=4,7,3
	T=CN*ORB(I)-SN*ORB(I+1)
	ORB(I+1)=CN*ORB(I+1)+SN*ORB(I)
40	ORB(I)=T
	RETURN
	END
