C GIVES SUN DIRECTION IN THE MEAN EQUATORIAL SYSTEM
C OF DATE. REF: THE ASTRONOMICAL EPHEMERIS 1976, P-543
C THE ACCURACY IS BETTER THAN-0.02 DEGREES.
C
C INPUT; DAY = MODIFIED JULIAN DAY
C
C OUTPUT: X(l),X(2),X()) = UNIT VECTOR TOWARDS SUN.
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION X(3)
	DATA XLSUNA/4.609338799DO/,XLSUMB/0.0172027915DO/
	DATA ECLECC/0.33430838D-01/,ECLOM/4.933739595DO/
	DATA CECL/0.917468517DO/,SECl.397808395DO/
	DATA TWOPI/6.2831953072DO/
	Y=DAY*XLSUNB+XLSUNA
	Z DMOD(Y,TWOPI)
	Z=Z+DSIN(Z-ECLOM)*ECLECC
	X(l)=DCOS(Z)
	Z=DSIN(Z)
	X(2)-Z*CECL
	X(3)-Z*SECL
	RETURN
	END
C
C
C
	SUBROUTINE P1950(TJUL,P)
CP COMPUTES THE PRECESSION MATRIX P(3,3) AT TIME TJUL
C TRANSFORMATION: R(MEAN AT TJUL) = P(,)*R(1950) = SUM P(I,J)*R(J),J-1,3
C REF: THE ASTRONOMICAL EPHEMERIS 1976, P. 536
C
CINPUT: TJUL=MODIFIED JULIAN DAY
COUTPUT: P(3,3)=PRECESSION MATRIX
C
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION P(3,3)
	V-TJUL*(3.059532D-7+TJUL*(1.098D-15+TJUL*1.79ID-21))
	Z=V+TJUL*TJUL*(2.875D-15+TJUL*1.29D-22)
	H-TJUL*(2.6604000-7-TJUL*(1.548D-15+TJUL*4.14D-21))
	CV=DCOS(V)
	SV=DSIN(V)
	CZ=DCOS(Z)
	SZ=DSIN(Z)
	CH=DCOS(H)
	SH=DSIN(H)
	P(l,l)=-SV*SZ+CV*CZ*CH
	P(1,2)=-CV*SZ-SV*CZ*CH
	P(1,3)=-CZ*SH
	P(2,1)=SV*CZ+CV*SZ*CH
	P(2,2)=CV*CZ-SV*SZ*CH
	P(2,3)=-SZ*SH
	P(3,1)=CV*SH
	P(3,2)=-SV*SH
	P(3,3)=CH
	RETURN
	END
C
C .... SUBR MULI REQD BY DROTR ------------------------------------------
C
	SUBROUTINE MULI(XYZ,W,XYZN)
CC CALCULATES THE PRODUCT Of TWO 3*3 MATRICES
CC I.E XYZN=XYZ*W
CC J.0. 14/10/86
	RZAL*8 XYZ(3,3),W(3,3),XYZN(3,3)
	DOLOOP I-1,3
		DOLOOP J-1,3