National Aeronautics and Space Administration

A report in the last EXOSAT Express (vol. 17, p. 3) dealt with the accurate determination of celestial positions using the LE1 CMA. That report mentioned that a calibration error in the Fine Sun Sensor (unit 1, FSS) caused the spacecraft to slowly drift when one star and the sun were used as position references. This drift was not reflected in any data supplied to users. The bulk of this report is made up of a subroutine listing which enables users to correct for the effect of the FSS misalignment (subroutine ATTFSS).

The error originally showed up as a drift in the detector x and y pixel coordinates of a bright source while the time resolved roll, pitch and yaw measurements from the housekeeping data (A511 to 513) indicated no trend in the spacecraft attitude. Figure 1 shows this drift, which has a rate of 0.26 arcseconds per hour.

The following software came originally from the Orbit and Attitude department at ESOC after a report of a new calibration of the FSS by them indicated that it was tilted by 10 arcminutes with respect to the star tracker coordinate system.

The time resolved position of the same source determined using this software is shown in figure 2. The clear drift of the source across the detector is reduced to an immeasurably small drift in celestial position, being less than 0.14 arcseconds per hour. Two hours into this observation the shift in celestial position due to this algorithm is 3.4 arcseconds.

Julian Osborne

	SUBROUTINE ATTFSS(MAT,SRA, SDE, ISTR2, DEBLUR, Y, FMTBUF, DAY)
   
C ON ENTRY:
C	MAT = ATTITUDE MATRIX READ FROM AUXILIARY DATA
C	SRA = R.A. OF GUIDE STAR FROM AUXILIARY DATA
C	SDE = DEC OF GUIDE STAR FROM AUXILIARY DATA
C	ISRT2 = 0 IF 2 STARS WERE USED FOR ATT RECONSTRUCTION (FROM AUX DATA)
C	DEBLUR = .TRUE. IF USER NEEDS CORRECTION FOR CURRENT VALUES OF A511-5
C	Y = ARRAY GIVING SPACECRAFT POSN W.R.T. EARTH (FROM ORBIT DATA)
C	FMTBUF = RECORD OF HOUSEKEEPING DATA AT THE TIME OF REQD ATTITUDE
C	DAY = ESOC MJD (SEE FOT HANDBOOK SECT 2, P5) AT TIME OF REQD ATTITUDE
C
C ON RETURN:
C	MAT = ATTITUDE MATRIX MODIFIED To CORRECT FOR THE FSS INDUCED DRIFT
C
	REAL*8 MAT(9),AMAT(9),AU,DIST
	REAL*8 SRA,SDE,SUN(3),CON,Y(6)
	REAL*8 DAY
	INTEGER*2 FMTBUF(984)
	LOGICAL*2 DEBLUR
	CON=0.01745329252DO	!	2PI/360
	AU=1.495979D8		! 1 AU IN KH
	IF(ISTR2.NE.O)THEN
	  CALL P1950(DAY,AMAT)	!	GET PRECH MATRIX FOR DATE
	  CALL TRANSP(AMAT)
	  CALL SUNDR(DAY,SUN)	! GET POSH OF SUN FROM EARTH
	  CALL MUL(AKAT,SUN)	! PRECESS TO OBS DATE
	  DOLOOP I=1, 3
	    SUN(I)=SUN(I)*AU-Y(I)	! GET POSH OF SUN FROM S/C
	  ENDL00P
	  DIST= OSQRT(SUN(I)*SUN(I)+SUN(2)*SUN(2)+SUN(3)*SUN(3))
	  DOLOOP I=1,3
	    SUN(I)=SUN(I)/DIST
	  ENDLOOP
	  SRA=SRA*CON
	  SDE=SDE*CON
	ENDIF
	CALL FSSCO(ISTR2,FMTBUF,MAT,SUN,SRA,SDE,DEBLUR,MAT,IRET)
	IF(IRET.EQ.1)WRITE(1,12)
12	FORMAT(' FSSCO: NO CONVERGENCE IN 5 ITERATIONS!')
	RETURN
	END
C
CC
C
	SUBROUTINE TRANSP(A), 3-D MATRIX TRANSPOSE
	REAL*8 A(9),S
	S=A(4)
	A(4)=A(2)
	A(2)=S
	S=A(7)
	A(7)=A(3)
	A(3)=S
	S=A(6)
	A(6)=A(S)
	A(8)=S
	RETURN
	END
C
C====== SUNDR AND P1950 SUPPLIED BY A SCHUETZ 14/10/86 FOR SUN POSN AT TIME====
C=========================================================OF OBS===============
C
	SUBROUTINE SUNDR(DAY,X)
C DUBLE PRECISION VERSION
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(1)=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.660400D-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
	REAL*8 XYZ(3,3),W(3,3),XYZN(3,3)
	DOLOOP I=1,3
	  DOLOOP J=1,3
            XYZN(I,J)-O.DO
            DOLOOP K=1, 3
	      XYZN(I,J)=XYZ(I,K)*W(K,J)+XYZN(I,J)
            ENDLOOP
          ENDLOOP
        ENDLOOP
        RETURN
        END
C
C
C
      SUBROUTINE MUL(A.B)
CC CALCULATES THE PRODUCT Of TWO MATRICES (3, 3)*(3,1)=(3, 1)
CC I.E B=A*B
CC J.0. 14/20/86
      REAL*8 A(3,3),B(3),C(3)
      DOLOOP I-1,3
	C(I)=0.D0
    	DOLOOP J=1, 3
          C(l)=A(I,J)*B(J)+C(I)
	ENDLOOP
      ENDLOOP
      DOLOOP I-1,3
        B(I)-C(I)
      ENDLOOP
      RETURN
      END


      SUBROUTINE FSSCO
     +(ISTR2,FMTBUF,ATTI,SUN,SRA,SDE,DEBLUR,ATTC,IRET)
C	ORIGINAL VERSION DATE	<870824.1904>
C
CP	FROM THE NOMINAL EXOSAT SPACECRAFT ATTITUDE, IT CALCULATES THE
CP	INSTANTANEOUS ATTITUDE CORRECTED (IF NECESSARY) FOR MOTION
CP	DUE TO USE OF FSS,FOR ATTITUDE CONTROL AND OPTIONALLY FOR MOTION
CP	DUE TO LIMIT CYCLING. IT SHOULD BE USED ONLY FOR STABLE POINTING
CP	PERIODS. (USE SLWAT FOR PERIODS WHEN OUTER LOOP IS OPEN).
CP	ONE MAY CALL FSSCO EVEN IF NOT USING THE FSS - ATTC WILL THEN
CP	JUST BE SET TO ATTI.
C
CC	PROJ-EXO,SUBJ=ATT,UTIL-SPEC,AUTH-A.SCKUETZ, JUNE 84
CC	MODIFIED BY J. STERNBERG TO SPLIT FSS STUFF OFF FROM SLEW STUFF(09/86)
C	ARGUMENTS :
C
CI      INPUT	FMTBUF	I2   (984) TELEMETRY FORMAT - HK RECORD
CI			           AT THE TIME FOR WHICH ATTITUDE IS REQUIRED
CI				   HALF-WORDS
Cl				   0- 39	IDENTIFICATION
CI				   40- 71	MODE ANALYSIS RESULTS
CI				   72- 583	HX TELEMETRY FORMAT
CI				   584- 711	HK DERIVED PARAMETERS
CI				   712- 726	TIME/QUALITY INFORMATION
CI				   727- 983	NOT RELEVANT
C
CI             ATTI	R8  (3,3)  3-AXIS ATTITUDE FROM AUX DATA
CI	       	                   (I,1) STAR TRACKER X-AXIS
CI	                           (1,2) STAR TRACKER Y-AXIS
CI	                           (1,3) STAR TRACKER Z-AXIS
C
C
C
CI            SUN	R8  (3)    S/C TO SUN UNIT VECTOR
C
CI	      SRA	R8  (1)	   RIGHT ASCENSION OF GUIDE STAR I (RADIANS)
CI	      SDE	R8  (1)    DECLINATION OF GUIDE STAR I (RADIANS)
CI
C
CI	      DEBLUR	L   (1)	   A FLAG TO INDICATE WHETHER ATTITUDE
Cl				   MATRIX SHOULD IN ADDITION BE CORRECTED
CI				   FOR MOTION DUE TO LIMIT CYCLE:
CI					.TRUE.=CORRECT IT
CI					.FALSE.= DON'T CORRECT IT, BECAUSE
CI					THE BLURRING IS EITHER TO BE IGNORED
CI					OR TO BE HANDLED OUTSIDE THIS S/W.
C
CO     OUTPUT ATTC	R8 (3,3)   3-AXIS ATTITUDE AT TIME OF "FMTBUF"
CO				   CORRECTED IF NECESSARY FOR FSS/DEBLURRING
C
CO	      IRET	I  (1)	   RETURN CODE	-
CO				   = 0 NO ERROR
CO				   BIT 15 (=LSB) SET IT ITERATIONS DON'T
CO				   CONVERGE IN "SUNCO" ROUTINE
CO				   BIT 14 SET IF UNRECOGNISABLE FORMAT
C
CS	DROTX	ROTATES	AN ATTITUDE MATRIX AROUND X,Y OR Z AXIS
CS	SUNCO	RECOMPUTE THE ATTITUDE MATRIX AS FUNCTION OF
CS	 		STAR REFERENCE POSITION AND FSS MEASUREMENTS
CS	MVW	MOVE n 16-BIT WORDS FROM 2ND ARGUMENT TO 3RD.
C
C
	IMPLICIT REAL*8 (A-H,O-Z)
	LOGICAL DEBLUR,SUMPRE
	INTEGER*2 FMTBUF(l),A5289,A103,AO59,LI,L2,L3
	INTEGER*4 I4,D32768
	REAL*8 ATTI(3,3),ATTC(3,3),SUN(3),ERR(3)
	INTEGER*2 I2(2)
	EQUIVALENCE (I4,I2)
	DATA D32768/32768/
	DATA SECA/.206264806247096D+06/	!RADS TO ARC SECS
	DATA SECFSH/2417.1656982082/	!RADS TO FSS HIGH-GAIN UNITS (256/3")
	DATA SECFSL/38674.651171331/	!RADS TO FSS LOW-GAIN UNITS (16/3")
	DATA A5289/709/
	DATA A103/182/
	DATA A059/175/
	IRET=0
C COPY ATTI INTO ATTC
	CALL MVW(36,ATTI,ATTC)
C CHECK FOR 2 STAR MODE
	IF(ISTR2.NE.0)THEN
C EXTRACT ICRS, WHICH IS PARAMETER A528. IT IS A 4-BIT PARAMETER, BUT
C HERE ONLY ONE BIT IS USED (THE LSB) CALLED "ICRS1"
C EXTRACT ICFL TOO (A529)
		I=ISHFT(FMTBUF(A5289),-7)
		ICFL=IAND(I,1)
		ICRSI=IAND(I/2,1)
C IF THE MONITORING IS ON ONE STAR AND SUN THE INERTIAL X (ROLL)
C ERROR IS NOT RELIABLE AND NEEDS TO BE RECOMPUTED FROM FSS MEAS.
C ALSO, CHECK SUN PRESENCE INDICATOR FOR FSS IX (A103, FR 13, WD 10 BT 6)
		SUNPRE=IAND(FMTBUF(A103),OO1OOO8).NE.0
		IF(ICRSI.EQ.0.AND.ICFL.EQ.O.AND.SUNPRE)THEN
C RETRIEVE FSS GAIN SELECTION (PARAM.AO59, FRAME 12,WORD 13, BIT 0)
			IGAIN=IAND(ISHFT(FMTBUF(A059),-7),1)
C RETRIEVE AND SUM ALL FOUR FSS MEASUREMENTS FROM TLM FORMAT
C PARAM.A163, FRAME 1, WORD 15, BITS 0-7
			JCOUNT=0
			NFSS=0
			L1=88
			L2=472
			L3=128
			DOLOOP IAD=LI, L2, L3
			   IFSS-IAND(FMTBUF(D),255)
			   IF(IFSS.GT.127) IFSS=IFSS-256
			   JCOUNT=JCOUNT+l
			   NFSS-NFSS + IFSS
			ENDLOOP
C CONVERT INTO RADIANS
    		       IF(IGAIN.NE.O)THEN
C    IGAIN=1=> COARSE RESOLUTION => HIGH GAIN
		         ROLL=NFSS/SECFSH
		       ELSE
C    IGAIN=O => FINE RESOLUTION => L0W GAIN
		         ROLL=NFSS/SECFSL
 		       ENDIF
		       IF(JCOUNT.NE.0) ROLL=ROLL/JCOUNT
C	CALL SUNCO ROUTINE TO GET ATTITUDE MATRIX UPDATED FOR ROLL
		       CALL SUNCO(ATTC,SUN,SRA,SDE,ROLL,ICOUNT)
C CHECK	# OF ITERATIONS, JUST OUT OF INTEREST.
		       IF(ICOUNT.GT.5)IRET=1
		ENDIF
	ENDIF
C
	IF(DEBLUR)THEN
C GET PARAMETERS A511, 512, 513 (TRUE XYZ ERRORS)
           DOLOOP I=1, 3
		I2(1)=0
		I2(2)=FMTBUF(697+I)
		ERR(I)=(I4-D32768)/SECA*0.33DO
  	   ENDLOOP
C APPLY ROTATIONS ERRORS AROUND Y AND Z
	   CALL DROTX(ATTC,ERR(2),2)
	   CALL DROTX(ATTC,ERR(3),3)
C	APPLY ROTATION ABOUT X, IF ROLL ERROR AVAILABLE
C	(IF ROLL ERROR CORRECTION ALREADY PERFORMED THEN SKIP)
	  IF(ICRS1.NE.O.AND.ICFL.EQ.O.AND.SUNPRE)CALL DROTX(ATTC,ERR(1),1)
	ENDIF
	RETURN
	END
C
C
C
	SUBROUTINE SUNCO(ATT,SUN,SRA,SDE,ROLL,ICOUNT)
C
CP SUBROUTINE TO RECOMPUTE THE EXOSAT ATTITUDE MATRIX AS FUNCTION OF
CP STAR REFERENCE POSITION AND SUN EFFECTIVE POSITION
C
CI ATT (3, 3)	R*8	OLD ATTITUDE MATRIX
CI			ATT(l,l),ATT(2,l),ATT(3,I) =	STAR TRACKER X-AXIS
CI			ATT(I,2),ATT(2,2),ATT(3,2) =	STAR TRACKER Y-AXIS
CI			ATT(1,3),ATT(2,3),ATT(3,3) =	STAR TRACKER Z-AXIS
CI SUN(3)	R*8	S/C TO SUN UNIT VECTOR
CI SRA		R*8	RIGHT ASCENSION OF REFERENCE STAR (RADIANS)
CI SDE		R*8	DECLINATION	OF REFERENCE STAR (RADIANS)
CI ROLL		R*8	ACTUAL ROLL ERROR MEASURED BY FSS (RADIANS)
CI			(POSITIVE WHEN SUN IS BELOW X-Y PLANE)
CO ATT(3,3)	R*8	UPDATED ATTITUDE MATRIX (TWO-WAY PARAMETER!)
C
CO ICOUNT	I	ITERATION COUNT (FOR TESTING)
C
C
CC   PROJ=EXO, SUBJ=ATT, UTIL=SPEC, AUTH=A.SCHUETZ/ECD/OAD,AUG'86
C
CS SUBROUTINES CALLED
CS DRCTR ROTATES A SET OF REFERENCE AXES AROUND AN ARBITRARY AXIS
CS DROTX ROTATES A SET OF REFERENCE AXES AROUND ONE OF THEM
CS MVW	MOVE n 16-BIT WORDS FROM 2ND ARGUMENT TO 3RD ARGUMENT
C
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION ATT(3,3),SUN(3).FSSDAT(3)
C
CV LOCAL VARIABLES
CV FSSDAT(3) R*8	FINE SUN SENSOR DATA
CV			(1) = MISALIGNEMENT IN ROLL	(RADIANS)
CV			(2) = MISALIGNEMENT IN PITCH	(RADIANS)
CV			(3) = SCALE FACTOR
CV RLIM	   MINIMUM ROLL ERROR CORRECTION ON RADIANS	(=10 ARCSEC)
CV STAR    STAR UNIT VECTOR
CV FSSR    FINE SUN SENSOR REFERENCE AXES IN INERTIAL REFERENCE SYSTEM
CV SUNS    SUN UNIT VECTOR IN FSS REFERENCE SYSTEM
C
	DIMENSION STAR(3),SUNB(3),FSSR(3,3)
	DATA FSSDAT/.00285DO,.00049DO,I.ODO/	! = 9.8	1.7 ARC MINS
	DATA RLIM/4.848D-5/
C
C  IF PRESENT CONVERT STAR POSITION INTO A UNIT VECTOR (USE SUNB)
C  IF NOT ASSUME STAR IS ON X-AXIS
	IF(SRA.EQ.O.DO.AND.SDE.EQ.O.DO)THEN
C  COMMENT THIS OCCURS DURING CERTAIN CONTINGENCIES, BUT SUNCO
CI WILL STILL BEHAVE CORRECTLY.
		CALL MVW(12,ATT(1,1),SUNB)
	ELSE
C  COMMENT NORMAL CASE - VALID STAR COORDINATES - CONVERT
		SUNB(1)=DCOS(SDE)*DCOS(SRA)
		SUNB(2)=DCOS(SDE)*DSIN(SRA)
		SUNB(3)=DSIN(SDE)
	ENDIF
C
C  EXPRESS STAR POSITION IN STAR TRACKER REFERENCE MATRIX
	STAR(1) =ATT(1,1)*SUNB(1)+ATT(2,1)*SUNB(2)+ATT(3,1)*SUNB(3)
	STAR(2) =ATT(1,2)*SUNB(1)+ATT(2,2)*SUNB(2)+ATT(3,2)*SUNB(3)
	STAR(3) =ATT(1,3)*SUNS(l)+ATT(2,3)*SUNB(2)+AIT(3,3)*SUNB(3)
C
	ICOUNT=0
C
    10 CONTINUE
	ICOUNT=ICOUNT+1
C
C INITIALISE FSS REFERENCE MATRIX TO STAR TRACKER REFERENCE MATRIX
	CALL MVW(36,ATT,FSSR)
C
C APPLY MISALlGNMEHTS IN ROLL AND PITCH
	CALL DROTX(FSSR,FSSDAT(1),1)
	CALL DROTX(FSSR,FSSDAT(2),2)
C
C COMPUTE PREDICTED SUN POSITION IN FSS REFERENCE
	SUNB(1)=FSSR(1,1)*SUN(1) + FSSR(2,1)*SUN(2) + FSSR(3,1)*SUN(3)
	SUNB(2)=FSSR(1,2)*SUN(1) + FSSR(2,2)*SUN(2) + FSSR(3,2)*SUN(3)
	SUNB(3)=FSSR(1,3)*SUN(1) + FSSR(2,3)*SUN(2) + FSSR(3,3)*SUN(3)
C
C DERIVE PREDICTED ROLL ANGLE ERROR	(= MINUS ROLL ANGLE)
	PROLL=-DATAN(SUNB(3)/SUNB(2))
C
C RESIDUAL ROLL ERROR = ACTUAL - PREDICTED
	DROLL=ROLL*FSSDAT(3) - PROLL
C
C IF RESIDUAL ROLL ERROR IS SMALL ENOUGH RETURN
	IF(DABS(DROLL).LT.RLIM.OR.ICOUNT.GE.10)RETURN
C
C APPLY THAT CORRECTION TO THE INITIAL ATTITUDE MATRIX WITH THE
C AXE OF ROTATION BEING ALONG THE STAR VECTOR AND ITERATE
	CALL DROTR(ATT, DROLL, STAR)
	GO TO 10
	END
C
C---DROTR & DROTX SUPPLIED BY A SCHUETZ 14/10/86 -----------------------
C
	SUBROUTINE DROTR(XYZ,PHI,A)
C
CP	ROTATE AROUND ARBITRARY AXIS BY ANGLE PHI
C
CP	PROJ-EXO,SUBJ=ATT,UTIL=MULT,AUTH=J.FERTIG,3/10/60
C
C	CALLING SEQUENCE
C
Cl	INPUT	XYZ	R8	(3,3)	XYZ(1,1),XYZ(2,I),XYZ(3,1) FOR X-AXIS
CI					XYZ(1,2),XYZ(2,2),XYZ(3,2) FOR Y-AXIS
CI					XYZ(1,3),XYZ(2,3),XYZ(3,3) FOR Z-AXIS
CI		PHI	R8	(1)	ROTATION ANGLE (RAD)
CI		A	R8	(3)	AXIS OF ROTATION IN S/C SYSTEM
CI					OF COORDINATES
CO	OUTPUT  XYZ	R8	(3,3)	X,Y,Z-AXES AFTER ROTATION
C
C
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION XYZ(3,3),A(3),W(3,3),XYZN(3,3)
	C=DCOS(PHI)
	S=DSIH(PHI)
	W(1,1)=C+(1.-C)*A(1)*A(1)
	W(2,1)=  (1.-C)*A(2)*A(1)+S*A(3)
	W(3,1)=  (1.-C)*A(3)*A(l)-S*A(2)
	W(1,2)=  (1.-C)*A(1)*A(2)-S*A(3)
	W(2,2)=C+(1.-C)*A(2)*A(2)
	W(3,2)=  (1.-C)*A(3)-A(2)+S*A(1)
	W(1,3)=  (1.-C)*A(1)*A(3)+S*A(2)
	W(2,3)=  (1.-C)*A(2)*A(3)-S*A(1)
	W(3,3)=C+(1.-C)*A(3)*A(3)
	CALL MUL1(XYZ,W,XYZN)	! REPLACES DMAMA3 CALL
	DO 20 I=1,3
	 DO 10 J=1,3
		XYZ(I,J)=XYZN(l,J)
 10	 CONTINUE
 20	CONTINUE
	RETURN
	END
C
C
C
	SUBROUTINE DROTX(A,FI,I)
CP ROTATES A SET OF AXES X,Y,Z AROUND ONE OF THEN BY AN ANGLE FI
CC PROJ=GEN,SUBJ=NUM,UTIL=MULT,AUTH=FRAITURE/SAB
CI A	R8	3*3	A(1,1),A(2,1),A(3,1) FOR X-AXIS
CI			A(1,2).A(2,2),A(3,2) FOR Y-AXIS
CI			A(1,3),A(2.3).A(3,3) FOR Z-AXIS
CI FI	R8	3	ROTATION ANGLE IN RADIANS
CI I	I4		ROTATION AXIS NUMBER (1 FOR AXIS OX)
CO A	R8	3*3	X,Y,Z-AXES AFTER ROTATION
C----------------------------------------------------------------
	IMPLICIT  REAL*8  (A-H,O-Z)
	DIMENSION A(3,J),J(4)
	SAVE J
	DATA J/2,3,1,2/
C ---------------------------------------------------------------
	K=J(I)
	L=J(I+1)
	CSF=DCOS(FI)
	SNF=DSIN(FI)
	DO 10 N=1,3
	  U=A(N,K)*CSF+A(M,L)*SNF
	  A(N,L)=-A(N,K)*SNF+A(N,L)*CSF
	  A(N,K)=U
10	CONTINUE
	RETURN
	END