C
C     LAST MODIFICATION ON JUN-1998 (on PSI data absorption correction)
C
C23456789012345678901234567890123456789012345678901234567890123456789012
************************************************************************
*       DATA REDUCTION/CORRECTION PACKAGE OF DR. A.CHANDRASEKARAN      *
************************************************************************
	CHARACTER*12 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	CHARACTER*12 OUT1,OUT2,OUT3
	REAL CB(9,100)
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
	COMMON/FILE/ FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OUT/ OUT1,OUT2,OUT3
	COMMON/ABS1/ WAVE,R11,R12,R13,R21,R22,R23,R31,R32,R33,CB,PSIINT
	COMMON/ABS2/ H,K,L,PHI,PHI1,PHI2,ABSS,ANGKI,TETA,MDEG,RTD
	CALL INFORM
	CALL OPTIONS
	IF(OPT1.EQ.0.AND.OPT2.EQ.0.AND.OPT3.EQ.0.AND.OPT4.EQ.0) THEN
	CALL EXTRAS
	ELSE
	ENDIF
	CALL FILES
	CALL INQUIRE
	WRITE(*,*)'  '
	WRITE(*,10)OUT1
10      FORMAT(1X,'SEE THE FILE   ',A12,'    FOR MORE INFORMATION.')
	WRITE(*,*)'  '
	OPEN(UNIT=6,FILE=OUT1,STATUS='UNKNOWN',FORM='FORMATTED')
	CALL INFORM
	CALL CADACS
	CALL DRDACS1
	CALL ABSACS1
	CALL BESTABS
	CLOSE (UNIT=6)
	WRITE(*,*)' '
	STOP '      *****     COMPLETED ALL JOBS.     *****'
	END
************************************************************************
	SUBROUTINE INFORM
	WRITE(6,*)' '
        WRITE(6,*)'              ***   XRAYACS:    ***'
	WRITE(6,*)'   A PROGRAM FOR DATA REDUCTION AND CORRECTIONS.'
	WRITE(6,*)'                  COMPILED BY'
	WRITE(6,*)'              DR. A.CHANDRASEKARAN'
        WRITE(6,*)'                DATED: 01-JUN-98'
	WRITE(6,*)' '
	WRITE(6,*)'   IT HAS 4 PROGRAMS: FIRST 2 WERE WRITTEN BY AND'
	WRITE(6,*)' LAST 2 WERE EXTREMELY MODIFIED/IMPROVED BY DR. ACS.'
	WRITE(6,*)' '
	WRITE(6,*)' THE PROGRAMS ARE: CADACS, DRDACS1, ABSACS1 & BESTABS'
	WRITE(6,*)' ALL THESE ARE COMPATIBLE WITH SHELX PROGRAMS.'
	WRITE(6,*)'         FOR MORE INFORMATION ASK "H".'
	WRITE(6,*)' '
	WRITE(6,*)' ****  IF USEFUL, IT DESERVES TO BE REFERRED!  ****'
        WRITE(6,*)' '
        WRITE(6,*)' A.Chandrasekaran, XRAYACS: program for single crysta
     1l X-ray data corrections,'
        WRITE(6,*)' 1998, Chemistry Department, University of Massachuse
     1tts, Amherst, USA.'
        WRITE(6,*)' '
        WRITE(6,*)'    WRITE PROBLEMS TO: chandru@chemistry.umass.edu'
	WRITE(6,*)' '
	RETURN
	END
************************************************************************
	SUBROUTINE OPTIONS
	CHARACTER*1 OPT
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
10      WRITE(*,*)' DO YOU WANT TO RUN CADACS (Y/N/H)?'
	READ(*,'(A1)')OPT
        IF (OPT.EQ.'Y'.OR.OPT.EQ.'y') OPT1=1
        IF (OPT.EQ.'N'.OR.OPT.EQ.'n') OPT1=0
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'N'.OR.OPT.EQ.'y'.OR.OPT.EQ.'n') GOTO 20
        IF (OPT.EQ.'H'.OR.OPT.EQ.'h') CALL HELPCAD
	GOTO 10
20      WRITE(*,*)' DO YOU WANT TO RUN DRDACS1 (Y/N/H)?'
	READ(*,'(A1)')OPT
        IF (OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	OPT2=1
25      WRITE(*,*)' DO YOU WANT TO OUTPUT F OR F^2? (1=F^1, 2=F^2)'
	READ(*,*,ERR=25)NTY1
	NTY1=ABS(NTY1)
	IF(NTY1.GT.2.OR.NTY1.EQ.0) GOTO 25
	ELSE
	ENDIF
        IF (OPT.EQ.'N'.OR.OPT.EQ.'n') OPT2=0
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'N'.OR.OPT.EQ.'y'.OR.OPT.EQ.'n') GOTO 30
        IF (OPT.EQ.'H'.OR.OPT.EQ.'h') CALL HELPDRD
	GOTO 20
30      WRITE(*,*)' DO YOU WANT TO RUN ABSACS1 (Y/N/H)?'
	READ(*,'(A1)')OPT
        IF (OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	OPT3=1
	IF(NTY1.EQ.1.OR.NTY1.EQ.2) THEN
	NTY2=NTY1
	GOTO 36
	ELSE
35      WRITE(*,*)' IS THE INPUT F OR F^2? (1=F^1, 2=F^2)'
	READ(*,*,ERR=35)NTY2
	NTY2=ABS(NTY2)
	IF(NTY2.GT.2.OR.NTY2.EQ.0) GOTO 35
	ENDIF
36      CONTINUE
	ELSE
	ENDIF
        IF (OPT.EQ.'N'.OR.OPT.EQ.'n') OPT3=0
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'N'.OR.OPT.EQ.'y'.OR.OPT.EQ.'n') GOTO 40
        IF (OPT.EQ.'H'.OR.OPT.EQ.'h') CALL HELPABS
	GOTO 30
40      WRITE(*,*)' DO YOU WANT TO RUN BESTABS (Y/N/H)?'
	READ(*,'(A1)')OPT
        IF (OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	OPT4=1
45      WRITE(*,*)' IS THE INPUT F OR F^2? (1=F^1, 2=F^2)'
	READ(*,*,ERR=45)NTY3
	NTY3=ABS(NTY3)
	IF(NTY3.GT.2.OR.NTY3.EQ.0) GOTO 45
	ELSE
	ENDIF
        IF (OPT.EQ.'N'.OR.OPT.EQ.'n') OPT4=0
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'N'.OR.OPT.EQ.'y'.OR.OPT.EQ.'n') GOTO 50
        IF (OPT.EQ.'H'.OR.OPT.EQ.'h') CALL HELPBES
	GOTO 40
50      RETURN
	END
************************************************************************
	SUBROUTINE HELPCAD
	WRITE(*,*)' '
	WRITE(*,*)' THIS REMOVES THE UNNECESSARY LINES FROM THE'
	WRITE(*,*)' RAW DATA AND MAKES THE INPUT COMPACT.'
	WRITE(*,*)' IT IS NEEDED TO RUN THE NEXT PROGRAM, DRDACS1.'
	WRITE(*,*)' THE INPUT FILE IS  NAME.CAD'
	WRITE(*,*)' THE OUTPUT IS      NAME.DRD'
	WRITE(*,*)' NAME IS ANY NAME.'
	WRITE(*,*)' '
	RETURN
	END
************************************************************************
	SUBROUTINE HELPDRD
	WRITE(*,*)' '
	WRITE(*,*)' THIS PROGRAM APPLIES APPROXIMATE ISOTROPIC DECAY'
	WRITE(*,*)' CORRECTION: LINEAR BETWEEN EACH 2 SETS OF INTENSITY'
	WRITE(*,*)' CONTROLS OR OVERALL-EXPONENTIAL BASED ON THE'
	WRITE(*,*)' NUMBER OF THE REFLN, AND LP CORRECTIONS TO THE DATA.'
	WRITE(*,*)' IT CAN OUTPUT EITHER F OR F^2 BY YOUR CHOICE.'
	WRITE(*,*)' INPUT FILES ARE NAME.CEL AND NAME.DRD'
	WRITE(*,*)' NAME.CEL CONTAINS AF,MA,WAVE,A,B,C,AL,BE,GA,INTCR.'
	WRITE(*,*)' AF= ATTENUATION FACTOR;  MA= MONOCHROMATOR ANGLE;'
	WRITE(*,*)' WAVE= WAVE LENGTH;  INTCR= NUMBER OF INT. CONTROLS.'
        WRITE(*,*)' THE NAME.CEL FILE IS IN FREE FORMAT'
	WRITE(*,*)' NAME.DRD FILE IS PRODUCED BY THE CADACS PART.'
	WRITE(*,*)' NAME IS ANY NAME.'
	WRITE(*,*)' '
	RETURN
	END
************************************************************************
	SUBROUTINE HELPABS
	WRITE(*,*)' '
	WRITE(*,*)' THIS PROGRAM APPLIES ABSORPTION CORRECTION USING'
	WRITE(*,*)' PSI DATA WITH A SEVENTH ORDER SMOOTH CURVE OF'
	WRITE(*,*)' VARIATION OF INTENSITIES IN THE PSI DATA.'
	WRITE(*,*)' IT CAN CORRECT F OR F^2 AND CAN USE SEVERAL'
	WRITE(*,*)' PSI DATA AT ONCE (UPTO 100!).'
	WRITE(*,*)' * IT GIVES A PSIPLOT TOO! *'
        WRITE(*,*)' Please edit the .OUT file and check the columns'
        WRITE(*,*)' INTCAL for each PSI data and if any abnormal'
        WRITE(*,*)' values are there, remove that set of PSI data'
        WRITE(*,*)' and rerun the program.'
	WRITE(*,*)' THE INPUT FILES ARE NAME.PSI & NAME.FF1/NAME.FF2'
	WRITE(*,*)' NAME.PSI SHOULD HAVE:'
	WRITE(*,*)' TITLE                      (1st LINE)'
	WRITE(*,*)' LAMBDA, NO OF DATA SETS    (2nd    ; FREE FORMAT)'
	WRITE(*,*)' ORIENTATION MATRIX         (3rd....; FREE FORMAT!!)'
	WRITE(*,*)' THEN THE PSI DATA SHOULD FOLLOW.'
	WRITE(*,*)' NAME.FF1 IS THE FILE WITH H,K,L,F,DF   AND'
	WRITE(*,*)' NAME.FF2 IS H,K,L,F2,DF2     (3I4,2F8.2)'
	WRITE(*,*)' THESE NAME.FF1/NAME.FF2 COMES FROM DRDACS1.'
	WRITE(*,*)' NAME IS ANY NAME.'
	WRITE(*,*)' '
	RETURN
	END
************************************************************************
	SUBROUTINE HELPBES
	WRITE(*,*)'  '
	WRITE(*,*)' THIS PROGRAM DOES ABSORPTION CORRECTION BASED ON '
	WRITE(*,*)' THE METHOD OF      N.WALKER AND D.STUART,'
	WRITE(*,*)'                 ACTA CRYST., A39, 158-166 (1983).'
	WRITE(*,*)' PROGRAM OF          F.UGOZZOLI'
	WRITE(*,*)'                  COMPUT. CHEM., 11, 109-120 (1987).'
	WRITE(*,*)'  '
	WRITE(*,*)' MODIFIED BY DR. A.CHANDRASEKARAN IN SEP-1994.'
	WRITE(*,*)'  '
	WRITE(*,*)' THE INPUTS FILES ARE NAME.MAT & NAME.FCF'
	WRITE(*,*)' NAME.MAT SHOULD HAVE TITLE(1st LINE),LAMBDA(2nd)'
	WRITE(*,*)' ORIENTATION MATRIX(3rd LINE ONWARDS; FREE FORMAT!!)'
	WRITE(*,*)' THEN NUMBER OF CYCLES OF REFINEMENT (DEF.=2).'
	WRITE(*,*)' NAME.FCF IS THE FILE PRODUCED BY SHELX PROGRAMS.'
	WRITE(*,*)' AFTER HAVING ALL NON-H ATOMS WITH ISOTROPIC '
	WRITE(*,*)' THERMAL PARAMETERS,  REFINE WITH'
	WRITE(*,*)' OMIT 0, LIST 5, AFLS 3, WGHT 0   (SHELX76)   or'
	WRITE(*,*)' LIST 4, L.S. 3, WGHT 0.1 AND NO OMIT (SHELX93)'
	WRITE(*,*)' THE UNNECESSARY LINES SHOULD BE DELETED AND THE'
	WRITE(*,*)' FILE SHOULD BE NAMED AS NAME.FCF'
        WRITE(*,*)' NAME.FFF (H,K,L,FO,DF,FC; 3I4,2F12.2,F10.2) CAN'
        WRITE(*,*)' ALSO BE USED INSTEAD OF NAME.FCF; NAME IS ANY NAME.'
	WRITE(*,*)
	RETURN
	END
************************************************************************
	SUBROUTINE FILES
	CHARACTER*7 FILE
	CHARACTER*12 OUT1,OUT2,OUT3
	CHARACTER*12 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
	COMMON/FILE/ FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OUT/ OUT1,OUT2,OUT3
	WRITE(*,*)' '
	WRITE(*,*)' ALL THE FILES SHOULD HAVE A COMMON FIRST NAME.'
	WRITE(*,*)' PLEASE GIVE THE FILES FIRST NAME:'
	READ(*,'(A7)')FILE
        FILE1=TRIM(FILE)//'.CAD'
        FILE2=TRIM(FILE)//'.CEL'
        FILE3=TRIM(FILE)//'.DRD'
        FILE4=TRIM(FILE)//'.PSI'
        FILE5=TRIM(FILE)//'.FF1'
        IF(NTY1.EQ.2) FILE5=TRIM(FILE)//'.FF2'
        IF(NTY2.EQ.2) FILE5=TRIM(FILE)//'.FF2'
        FILE6=TRIM(FILE)//'.MAT'
        FILE7=TRIM(FILE)//'.FCF'
        FILE8=TRIM(FILE)//'.FFF'
        OUT1=TRIM(FILE)//'.OUT'
        OUT2=TRIM(FILE)//'.ABS'
        OUT3=TRIM(FILE)//'.BES'
	RETURN
	END
************************************************************************
	SUBROUTINE INQUIRE
	CHARACTER*1 OPT
	CHARACTER*12 OUT1,OUT2,OUT3
	CHARACTER*12 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
	COMMON/FILE/ FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OUT/ OUT1,OUT2,OUT3
	WRITE(*,*)' '
	WRITE(*,*)' THE FILES NEEDED ARE:'
	IF(OPT1.EQ.1) WRITE(*,'(23X,A12)')FILE1
	IF(OPT2.EQ.1) WRITE(*,'(23X,A12)')FILE2
	IF(OPT1.EQ.0.AND.OPT2.EQ.1) WRITE(*,'(23X,A12)')FILE3
	IF(OPT3.EQ.1) WRITE(*,'(23X,A12,A7,A12)')FILE4
	IF(OPT2.EQ.0.AND.OPT3.EQ.1) WRITE(*,'(23X,A12)')FILE5
	IF(OPT4.EQ.1) WRITE(*,'(23X,A12)')FILE6
	IF(OPT4.EQ.1) WRITE(*,'(23X,A12,A7,A12)')FILE7,'  or     ',FILE8
	WRITE(*,*)' '
10      WRITE(*,*)' HAVE YOU GOT READY THE ABOVE FILES [(Y)/N/Stop]?'
	READ(*,'(A1)')OPT
        IF (OPT.EQ.'S'.OR.OPT.EQ.'s') THEN
        STOP ' PLEASE RESTART WITH FILES READY.'
        ELSE
        ENDIF
        IF (OPT.NE.'N'.AND.OPT.NE.'n') RETURN
	IF(OPT1.EQ.1) THEN
	WRITE(*,*)' DO YOU WANT TO CHANGE    ',FILE1
	READ(*,'(A1)')OPT
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	WRITE(*,*)' GIVE NEW FILENAME'
	READ(*,'(A12)')FILE1
	ELSE
	ENDIF
	ELSE
	ENDIF
	IF (OPT2.EQ.1) THEN
	WRITE(*,*)' DO YOU WANT TO CHANGE    ',FILE2
	READ(*,'(A1)')OPT
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	WRITE(*,*)' GIVE NEW FILENAME'
	READ(*,'(A12)')FILE2
	ELSE
	ENDIF
	ELSE
	ENDIF
	IF(OPT1.EQ.0.AND.OPT2.EQ.1) THEN
	WRITE(*,*)' DO YOU WANT TO CHANGE    ',FILE3
	READ(*,'(A1)')OPT
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	WRITE(*,*)' GIVE NEW FILENAME'
	READ(*,'(A12)')FILE3
	ELSE
	ENDIF
	ELSE
	ENDIF
	IF(OPT3.EQ.1) THEN
	WRITE(*,*)' DO YOU WANT TO CHANGE    ',FILE4
	READ(*,'(A1)')OPT
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	WRITE(*,*)' GIVE NEW FILENAME'
	READ(*,'(A12)')FILE4
	ELSE
	ENDIF
	ELSE
	ENDIF
	IF(OPT2.EQ.0.AND.OPT3.EQ.1) THEN
	WRITE(*,*)' DO YOU WANT TO CHANGE    ',FILE5
	READ(*,'(A1)')OPT
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	WRITE(*,*)' GIVE NEW FILENAME'
	READ(*,'(A12)')FILE5
	ELSE
	ENDIF
	ELSE
	ENDIF
	IF(OPT4.EQ.1) THEN
	WRITE(*,*)' DO YOU WANT TO CHANGE    ',FILE6
	READ(*,'(A1)')OPT
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	WRITE(*,*)' GIVE NEW FILENAME'
	READ(*,'(A12)')FILE6
	ELSE
	ENDIF
	WRITE(*,*)' DO YOU WANT TO CHANGE    ',FILE7
	READ(*,'(A1)')OPT
        IF(OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	WRITE(*,*)' GIVE NEW FILENAME'
	READ(*,'(A12)')FILE7
	ELSE
	ENDIF
	ELSE
	ENDIF
	RETURN
	END
************************************************************************
************************************************************************
*       -----CADACS.FOR-----                                           *
************************************************************************
	SUBROUTINE CADACS
	CHARACTER*12 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
	COMMON/FILE/ FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	DIMENSION LINE(65)
	IF(OPT1.EQ.0) RETURN
	WRITE(6,*)' '
	WRITE(6,*)' *****    DOING  CADACS.....     *****'
	WRITE(*,*)' '
	WRITE(*,*)' *****    DOING  CADACS.....     *****'
	OPEN(UNIT=1,FILE=FILE1,STATUS='OLD',FORM='FORMATTED')
	OPEN(UNIT=2,FILE=FILE3,STATUS='UNKNOWN',FORM='FORMATTED')
	NN=26
	N=0
 2      READ(1,100,END=5)(LINE(I),I=1,65)
100     FORMAT(65A1)
	IF (LINE(27).EQ.'N'.OR.LINE(27).EQ.'I') THEN
	NN=27
	N=1
	GOTO 1000
	ELSE
	IF (LINE(26).EQ.'N'.OR.LINE(26).EQ.'I') GOTO 1000
	ENDIF
	GOTO 2
1000    CONTINUE
	REWIND 1
 3      READ(1,100,END=4)(LINE(I),I=1,65)
	IF (LINE(NN).NE.'N'.AND.LINE(NN).NE.'I') GOTO 3
	WRITE(2,110)(LINE(I),I=5+N,31+N),(LINE(M),M=40+N,61+N)
110     FORMAT(1X,49A1)
	GO TO 3
  4     CONTINUE
	REWIND 2
	CLOSE (UNIT=1)
	CLOSE (UNIT=2)
	WRITE(6,*)'     *****     CADACS  COMPLETED..     *****'
	WRITE(6,*)' '
	WRITE(*,*)'     *****     CADACS  COMPLETED..     *****'
	WRITE(*,*)' '
5       RETURN
	END
************************************************************************
************************************************************************
*       -----DRDACS1.FOR-----                                          *
************************************************************************
	SUBROUTINE DRDACS1
	CHARACTER*1 CORLOG
	CHARACTER*12 OUT1,OUT2,OUT3
	CHARACTER*12 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
	COMMON/FILE/ FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OUT/ OUT1,OUT2,OUT3
	INTEGER NNN,NN
	PARAMETER (NNN=600,NN=10)
	CHARACTER*27 AMES
	INTEGER NMAX(NNN),HST(NN),KST(NN),LST(NN),NSECT(NNN)
	REAL MF,LB,TC,RB,NF,ESD
	DIMENSION SI(NNN,NN),SCALE1(NNN),SCALE(NNN),TSI(NN),ASCAL(NN)
	IF(OPT2.EQ.0) RETURN
	WRITE(6,*)' *****     DOING  DRDACS1....     *****'
	WRITE(6,*)' '
	WRITE(*,*)' *****     DOING  DRDACS1....     *****'
	OPEN(UNIT=1,FILE=FILE2,STATUS='OLD',FORM='FORMATTED')
	READ(1,*,ERR=10)AF,CMA,WAVE,A,B,C,AL,BE,GA,IMAX
10      CONTINUE
	IF(A.LT.1.0.OR.B.LT.1.0.OR.C.LT.1.0) THEN
	WRITE(*,'(A12)')' PLEASE CHECK THE FILE  ',FILE2
	STOP ' RERUN AGAIN!'
	ELSE
	ENDIF
	WRITE(6,11)AF,CMA,WAVE
 11     FORMAT(///2X,'A.F.=',F7.3,', M.A.=',F6.3,', LAMBDA=',F7.4,','/)
	WRITE(6,12)A,B,C,AL,BE,GA
 12     FORMAT(2X,'A=',F7.3,', B=',F7.3,', C=',F7.3,', AL=',F7.3,
     1  ', BE=',F7.3,', GA=',F7.3,','/)
	WRITE(6,1199)IMAX
 1199   FORMAT(1X,I2,' INTENSITY CONTROL REFLECTIONS.'//)
	RAD=1.0/57.29577951
	COSA=COS(AL*RAD)
	COSB=COS(BE*RAD)
	COSG=COS(GA*RAD)
	SINA=SIN(AL*RAD)
	SINB=SIN(BE*RAD)
	SING=SIN(GA*RAD)
	D=SQRT(1.0+2.0*COSA*COSB*COSG-COSA**2-COSB**2-COSG**2)
	AST=SINA/(D*A)
	BST=SINB/(D*B)
	CST=SING/(D*C)
	DH=(COSB*COSG-COSA)/(SINB*SING)
	DK=(COSA*COSG-COSB)/(SINA*SING)
	DL=(COSA*COSB-COSG)/(SINA*SINB)
	CMA=CMA*RAD*2.0
	CCC=COS(CMA)
	M=1
	N=1
	ALL=0.0
	TN=0.0
	NREF=1
	CLOSE (UNIT=1)
	OPEN(UNIT=2,FILE=FILE3,STATUS='OLD',FORM='FORMATTED')
	IF(IMAX.EQ.0) THEN
	NOPT=2
	SCALED=1.0
	GOTO 146
	ELSE
	ENDIF
100     READ(2,20,END=110,ERR=111)LH,LK,LL,CODE,MF,LB,TC,RB
 20     FORMAT(6X,3I5,1X,A1,5X,F3.0,F6.0,F7.0,F6.0)
	NC=1
	NREF=NREF+1
	IF(LH.EQ.999)GO TO 100
	IF(LH.EQ.0.AND.LK.EQ.0.AND.LL.EQ.0)GO TO 100
	NSECT(M)=NSECT(M)+1
	IF(CODE.NE.'I')GO TO 120
	IF(M.NE.1)GO TO 105
	HST(N)=LH
	KST(N)=LK
	LST(N)=LL
105     CONTINUE
	IF(LH.NE.HST(N).OR.LK.NE.KST(N).OR.LL.NE.LST(N))GO TO 115
	IF(MF.EQ.0.0) GO TO 111
	RINT=(TC-2*(LB+RB))/ABS(MF)
	IF(MF.LT.0) RINT=RINT*AF
	SI(M,N)=RINT
	ALL=ALL+1.0
	N=N+1
	IF(N.LE.IMAX)GO TO 100
	GO TO 130
115     CONTINUE
	SI(M,N)=0.0
	N=N+1
	IF(N.LE.IMAX)GO TO 105
	M=M+1
	N=1
	NC=NC+1
	IF(NC.EQ.2)GO TO 105 
	GO TO 111
130     CONTINUE
	NMAX(M)=ALL
	TN=TN+ALL
	ALL=0.0
	N=1
	M=M+1
120     CONTINUE
	GO TO 100
111     WRITE(*,30)NREF,FILE3
 30     FORMAT(/2X,'ERROR !!!;  CHECK LINE NUMBER ',I5,' IN FILE  ',A12)
	STOP ' RERUN AGAIN!'
110     CONTINUE
	MMAX=M-1
	IF(N.LE.IMAX.AND.N.GT.1) MMAX=M
	ALL=0.0
	DO 135 N=1,IMAX
	TSI(N)=0.0
	DO 136 M=1,MMAX
	IF(SI(M,N).EQ.0)GO TO 136
	TSI(N)=TSI(N)+SI(M,N)
	ALL=ALL+1.0
136     CONTINUE
	ASCAL(N)=TSI(N)/ALL
	ALL=0.0
135     CONTINUE
	WRITE(6,141)
 141    FORMAT(' VARIATION ANALYSIS OF INTENSITY CONTROL REFLECTIONS')
	WRITE(6,142)
 142    FORMAT(4X,'H',3X,'K',3X,'L',6X,'VARIATION',7X,'DECAY(?)',
     1  5X,'NO. OF DATA',6X,'MEAN INTENSITY'/)
	TDECAY=0.0
	TDECAY1=0.0
	DO 137 N=1,IMAX
	ESD=0.0
	I=0
        DO 138 M=2,MMAX
	IF(SI(M,N).EQ.0) GO TO 138
	ESD=ESD+ABS(1.0-SI(M,N)/ASCAL(N))
	I=I+1
138     CONTINUE
	ESD=ESD*100/I
	TDECAY=TDECAY+ESD
        DECAY=(SI(1,N)-SI(MMAX,N))/SI(1,N)*100
	TDECAY1=TDECAY1+DECAY
        IALL=IALL+I+1
        WRITE(6,143)HST(N),KST(N),LST(N),ESD,DECAY,(I+1),ASCAL(N)
 143    FORMAT(1X,3I4,5X,F8.2,' %',4X,F9.2,' %',8X,I4,13X,F8.2)
137     CONTINUE
	NREF=NREF-1
	IF(NREF.LE.0) NREF=1
	NWEAK=0
	NOPT=0
	TDECAY=TDECAY*2.0/IMAX
	TDECAY1=TDECAY1/IMAX
	WRITE(6,144)TDECAY,TDECAY1
	WRITE(*,*)' '
	WRITE(*,144)TDECAY,TDECAY1
 144    FORMAT(' LINEAR MEAN DEACY = ',F6.2,' % ;  FINAL DECAY = ',
     1 F6.2,' %')
	IF(TDECAY.LT.10.0) GO TO 139
	TDECAY=TDECAY/100.0
	TXTAL=1.0-TDECAY1/100.0
	TXTAL=LOG(1.0/TXTAL)
	WRITE(*,*)' '
	WRITE(*,*)' THE DEFAULT DECAY CORRECTION IS LINEAR BETWEEN TWO'
	WRITE(*,*)' CONTROLS. USE THE DEFAULT FOR MULTI-XTAL DATA.'
	WRITE(*,*)' DO YOU WANT EXPONENTIAL DECAY CORRECTION? [Y/(N)]'
	READ(*,'(A1)')CORLOG
	IF(CORLOG.EQ.'Y'.OR.CORLOG.EQ.'y') NOPT=1
139     CONTINUE
	M=1
	N=1
	ALL=0.0
	BSCALE=0.0
140     CONTINUE
	IF(SI(M,N).EQ.0)GO TO 145
	BSCALE=BSCALE+1.0-SI(M,N)/SI(1,N)
	ALL=ALL+1.0
145     CONTINUE
	N=N+1
	IF(N.LE.IMAX)GO TO 140
	SCALE1(M)=BSCALE/ALL
	SCALE(M)=SCALE1(M)-SCALE1(M-1)
	WRONG=(-1.0)*ALL
	IF(SCALE(M).EQ.WRONG) SCALE(M)=0.0
	RDECAY=SCALE(M)*100.0
	WRITE(6,*)' DECAY WITHIN TWO INTENSITY CONTROLS   = ',RDECAY,' %'
	WRITE(6,*)' NUMBER OF REFLNS IN THIS SECTION      = ',NSECT(M)
	M=M+1
	N=1
	ALL=0.0
	BSCALE=0.0
	BSCALE1=0.0
	IF(M.LE.MMAX)GO TO 140
	M=1
	N=1
	REWIND 2
146     CONTINUE        
	I=1
	WRITE(6,*)' '
	WRITE(6,*)' '
	WRITE(6,152)
 152    FORMAT(2X,'BAD REFLECTIONS WITH BACKGROUND RATIO > 5.0')
	WRITE(6,*)' '
	WRITE(*,*)' '
	OPEN(UNIT=3,FILE=FILE5,STATUS='UNKNOWN',FORM='FORMATTED')
	NT=NTY1+2
	WRITE(3,98)NT
98      FORMAT('HKLF -',I1)
150     READ(2,20,END=2000,ERR=111)LH,LK,LL,CODE,MF,LB,TC,RB
	IF(LB.EQ.0) LB=1.0
	IF(RB.EQ.0) RB=1.0
	IF(LH.EQ.999) GO TO 150
	IF(MF.EQ.0.0) GO TO 150
	IF(LH.EQ.0.AND.LK.EQ.0.AND.LL.EQ.0) GO TO 150
C
	NREF1=NREF1+1
	NREF2=NREF2+1
	IF(LB.GT.RB) THEN
	BRR=LB/RB
	ELSE
	BRR=RB/LB
	ENDIF
	IF(BRR.GT.5.0) THEN
	IF(BRR.GE.5.OR.LB.GT.TC.OR.RB.GT.TC) THEN
	AMES=' ** OMIT THIS REFLECTION **'
	WRITE(*,151)LH,LK,LL,LB,TC,RB,AMES
	WRITE(6,151)LH,LK,LL,LB,TC,RB,AMES
 151    FORMAT(1X,3I4,3F10.1,2X,A27)
	ELSE
	WRITE(6,151)LH,LK,LL,LB,TC,RB
	ENDIF
	ELSE
	ENDIF
C
	IF(CODE.NE.'I')GO TO 160
	N=N+1
	IF(N.LE.NMAX(M))GO TO 150
	M=M+1
	N=1
	NREF2=0
	IF(M.LE.MMAX)GO TO 150
	IF(IMAX.EQ.0) GO TO 150
	GO TO 2000
160     CONTINUE
170     CONTINUE
	RINT=(TC-2*(LB+RB))
	IF(MF.LT.0) THEN
	NF=AF
	ELSE
	NF=1.0
	ENDIF
	MF=ABS(MF)
	IF (NOPT-1) 171,172,173
171     CONTINUE
	SCALED=1.0/(1.0-SCALE1(M-1)-SCALE(M)*NREF2/NSECT(M))
	GOTO 173
172     CONTINUE        
	SCALED=TXTAL*NREF1/NREF
	SCALED=EXP(SCALED)
173     CONTINUE
	RINT1=RINT*NF*SCALED/MF
	RINTS=(TC+2*(RB+LB))*(SCALED*NF/MF)**2
	SINTH1=(LH*AST)**2+(LK*BST)**2+(LL*CST)**2
	SINTH2=2*(LH*LK*AST*BST*DL+LH*LL*AST*CST*DK+LK*LL*BST*CST*DH)
	SINTH=SINTH1+SINTH2
	SINTH=SQRT(SINTH)
	SINTH=SINTH*WAVE/2.0
	COSTH=SQRT(1-SINTH**2)
	COS2T=COSTH**2-SINTH**2
	C=COS2T**2
	SIN2T=2*SINTH*COSTH
	SLP=(1+CCC)*SIN2T/(CCC+C)
	FSQ=RINT1*SLP
	SIGF2=SQRT(RINTS+(0.01*RINT1)**2)*SLP
	IF(NTY1.EQ.2) GOTO 180
	IF(FSQ.LE.0.25) THEN
	FF=FSQ
	IF(FF.LT.0.0) FF=0.0
	FF=SQRT(FF)
	SIGF=SIGF2
	ELSE
	FF=SQRT(FSQ)
	SIGF=SQRT((RINTS+(0.01*RINT1)**2)/RINT1*SLP)/2.0
	ENDIF
	IF(FF.LT.4.0*SIGF) NWEAK=NWEAK+1
	WRITE(3,50)LH,LK,LL,FF,SIGF
 50     FORMAT(3I4,2F8.2)
	GOTO 190
180     CONTINUE
	IF(FSQ.LT.2.0*SIGF2) NWEAK=NWEAK+1
        WRITE(3,51)LH,LK,LL,FSQ,SIGF2
 51     FORMAT(3I4,F8.1,F8.2)
190     CONTINUE
	GO TO 150
2000    CONTINUE
        WRITE(3,*)' '
	CLOSE (UNIT=1)
	CLOSE (UNIT=2)
	CLOSE (UNIT=3)
        WRITE(6,60)NREF1
 60     FORMAT(/5X,'TOTAL NUMBER OF MEASUREMENTS =',I5)
        WRITE(6,61)(NREF1-IALL)
 61     FORMAT(/5X,'TOTAL NUMBER OF REFLECTIONS  =',I5/)
	WRITE(6,62)NWEAK
 62     FORMAT(5X,'TOTAL NUMBER OF WEAK MEASUREMENTS =',I5)
3000    CONTINUE
	WRITE(6,*)' '
	WRITE(6,*)'     *****     DRDACS1  COMPLETED.    *****'
	WRITE(*,*)'     *****     DRDACS1  COMPLETED.    *****'
	RETURN
	END
************************************************************************
************************************************************************
*             ABSORPTION CORRECTION USING MANY PSI DATA                *
************************************************************************
	SUBROUTINE ABSACS1
	CHARACTER*12 OUT1,OUT2,OUT3
	CHARACTER*12 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	REAL KAPPA(37),F,DF,F1,POW,START(100),B(9),CB(9,100)
	DIMENSION TITLE(20),FFF(37),PPHI(37),PPSI(37),INT(37),LB(37)
	DIMENSION IB(37),NPI(37),DELTA(37),OMK(37),THETA(37),OME(37)
	INTEGER H,K,L
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
	COMMON/FILE/ FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OUT/ OUT1,OUT2,OUT3
	COMMON/ABS1/ WAVE,R11,R12,R13,R21,R22,R23,R31,R32,R33,CB,PSIINT
	COMMON/ABS2/ H,K,L,PHI,PHI1,PHI2,ABSS,ANGKI,TETA,MDEG,RTD
	IF(OPT3.EQ.0) GOTO 2000
	WRITE(6,*)' '
	WRITE(6,*)' *****     DOING   ABSACS1.....     *****'
	WRITE(*,*)' '
	WRITE(*,*)' *****     DOING   ABSACS1.....     *****'
	WRITE(6,1001)NTY2
 1001   FORMAT(/'  ABSORPTION CORRECTION USING PSI DATA (ON F^',I1,')'/)
	NPTS=37
	MDEG=8
	TMIN=10.0
	TMAX=0.0
	RTD=57.29577951
	DTR=1.0/RTD
	ALFA=50.0
	COSA=COS(ALFA*DTR)
	OPEN(UNIT=1,STATUS='OLD',FILE=FILE4,ERR=10,FORM='FORMATTED')
	GOTO 11
10      WRITE(*,*)'  PSI DATA FILE OPEN ERROR !!!'
	STOP
11      READ(1,1101)(TITLE(I),I=1,15)
 1101   FORMAT(19A4)
	WRITE(6,1002)(TITLE(I),I=1,15)
 1002   FORMAT(2X,' FOR ',19A4)
	READ(1,*)WAVE,NDS
	WRITE(6,1003)WAVE,NDS,NPTS
 1003   FORMAT(/2X,'WAVE LENGTH =',F8.5,', NO. OF DATA SETS = ',I5,
     1/' NO. OF DATA POINTS/SET = ',I3,', FORMAT(3I4,2F8.2)')
	READ(1,*)R11,R12,R13,R21,R22,R23,R31,R32,R33
	WRITE(6,1004)
 1004   FORMAT(/2X,' THE ORIENTATION MATRIX ')
	WRITE(6,1005)R11,R12,R13,R21,R22,R23,R31,R32,R33
 1005   FORMAT(/,5X,'R11= ',F9.6,' R12= ',F9.6,' R13= ',F9.6,
     1  /,5X,'R21= ',F9.6,' R22= ',F9.6,' R23= ',F9.6,
     2  /,5X,'R31= ',F9.6,' R32= ',F9.6,' R33= ',F9.6/)
C-----  PPHI=PSI & PPSI=PHI  OF THE DIFFRACTOMETER ANGLES
	DO 1000 ND=1,NDS
	DO 100 I=1,NPTS
99      READ(1,1102)H,K,L,PPHI(I),NPI(I),LB(I),INT(I),IB(I)
 1102   FORMAT(10X,3I5,7X,F7.2,1X,I3,I6,I7,I6)
	IF(NPI(I).EQ.0.AND.INT(I).EQ.0) GOTO 99
	READ(1,*)NDUM,NDUM1,THETA(I),PPSI(I),OMK(I),KAPPA(I)
100     CONTINUE
	WRITE(6,1006)H,K,L
 1006   FORMAT(/2X,' THE PSI PLOT IS FOR THE REFLECTION ',3I5)
	WRITE(6,1007)
 1007   FORMAT(/7X,'PPHI',2X,'NPI',2X,'LP',2X,'PEAK',2X,'RB',5X,'PPSI',
     1  6X,'OMK',6X,'KAPPA',4X,'INT')
	DO 110 I=1,NPTS
	FFF(I)=(INT(I)-2.0*(LB(I)+IB(I)))/(1.0*NPI(I))
	KAPPA(I)=KAPPA(I)/2.0
	DELTA(I)=ATAN(COSA*TAN(KAPPA(I)*DTR))
	DELTA(I)=DELTA(I)/DTR
	PPSI(I)=PPSI(I)+DELTA(I)
	OME(I)=OMK(I)+DELTA(I)
	WRITE(6,1008)PPSI(I),NPI(I),LB(I),INT(I),IB(I),PPHI(I),OMK(I),
     1  KAPPA(I),FFF(I)
 1008   FORMAT(2X,F10.3,I4,I4,I6,I4,4F10.3)
        IF(PPSI(I))120,120,130
120     PPSI(I)=PPSI(I)+360.0
        IF(PPSI(I).GT.180.0) PPSI(I)=PPSI(I)-180.0
130     CONTINUE
110     CONTINUE
	START(ND)=(FFF(1)+FFF(NPTS))/2.0
	CALL POLY(FFF,PPSI,MDEG,NPTS,B)
	WRITE(6,1009)
 1009   FORMAT(/5X,'THE COEFFICIENTS ARE ...'/)
	WRITE(6,1010)(B(I),I=1,MDEG)
 1010   FORMAT(2X,8(E12.6,1X)/)
	SDIF=0.0
	SUM=0.0
        DO 140 I=1,MDEG+1
140     CB(I,ND)=B(I)
	WRITE(6,1011)
 1011   FORMAT(/,10X,'   PHI',10X,'INTOBS',10X,'INTCAL',10X,'INTDIF',
     1  9X,'CAL/OBS',//)
	DO 150 I=1,NPTS
	ANTCAL=0.0
	DO 160 I1=1,MDEG+1
160     ANTCAL=ANTCAL+B(I1)*(PPSI(I)**(I1-1))
	DIFF=ANTCAL-FFF(I)
	RATIO=ANTCAL/FFF(I)
	SUM=SUM+FFF(I)
	SDIF=SDIF+ABS(DIFF)
        IF(ANTCAL.LE.0.0) THEN
        IF(H.NE.HBAD.AND.K.NE.KBAD.AND.L.NE.LBAD) THEN
        HBAD=H
        KBAD=K
        LBAD=L
        WRITE(*,*)' ERROR IN MODELLING!!!!!!!!!!!!!!!'
        WRITE(*,*)' PLEASE REMOVE ALL THE 37 REFLECTIONS (74 LINES)'
        WRITE(*,*)' BELONGING TO THE REFLECTION WITH THE H K L VALUES'
        WRITE(*,'(3I4)')H,K,L
        WRITE(*,*)' AND RERUN THE ABSORPTION CORRECTION PROGRAM! '
        WRITE(*,*)' '
        ELSE
        ENDIF
        ELSE
        ENDIF
	WRITE(6,1012)PPSI(I),FFF(I),ANTCAL,DIFF,RATIO
 1012   FORMAT(5(6X,F10.3))
150     CONTINUE
	RFAC=SDIF/SUM
	WRITE(6,1013)SUM,SDIF,RFAC
 1013   FORMAT(//,'  SUM =',E17.5,'  SDIF = ',E17.5,' RFAC = ',F10.5/)
C       WRITE(6,1014)
 1014   FORMAT(/,5X,'H',3X,'K',3X,'L',7X,'PHI',6X,'PHI1',6X,'PHI2',
     1  9X,'ABSCOR',6X,' OLD',4X,'  NEW')
1000    CONTINUE
	CALL PSIPLOT(NDS,START,CB,TITLE)
	NUM=0
	SUMABS=0.0
	POW=NTY2/2.0
        WRITE(6,*)' '
        WRITE(6,*)' REFLECTIONS WITH POOR TRANSMISSIONS ARE BELOW (IF PR
     1ESENT):'
        WRITE(6,*)' '
        WRITE(6,*)'  H   K   L     F       SIGMA    TRANS.'
	OPEN(UNIT=2,FILE=FILE5,STATUS='OLD',FORM='FORMATTED')
	OPEN(UNIT=3,FILE=OUT2,STATUS='UNKNOWN',FORM='FORMATTED')
        IF(NTY2.EQ.2) THEN
        WRITE(3,'(A7)')'HKLF -4'
        ELSE
        WRITE(3,'(A7)')'HKLF -3'
        ENDIF
170     READ(2,1103,ERR=170,END=500)H,K,L,F,DF
 1103   FORMAT(3I4,F8.1,F8.2)
	IF(H.EQ.999) GOTO 170
	IF(H.EQ.0.AND.K.EQ.0.AND.L.EQ.0) GOTO 170
	SUMABSS=0.0
	DO 180 ND=1,NDS
	CALL ABSORB(START(ND),ND)
	SUMABSS=SUMABSS+ABSS
180     CONTINUE
	ABSS=SUMABSS/NDS
	TRANS=1.0/ABSS
	ABC=ABS(ABSS)
	ABC=ABC**POW
	F1=F*ABC
	DF=DF*ABC
	IF(TMIN.GT.TRANS) TMIN=TRANS
	IF(TMAX.LT.TRANS) TMAX=TRANS
	SUMABS=SUMABS+TRANS
	NUM=NUM+1
C       WRITE(6,1015)H,K,L,PHI,PHI1,PHI2,ABC,F,F1
 1015   FORMAT(2X,3I4,8(2X,F9.3))
        IF (TRANS.LT.0.1.OR.TRANS.GT.10) THEN
        WRITE(6,'(3I4,2F9.2,F9.4)')H,K,L,F1,DF,TRANS
        ELSE
        ENDIF
        IF(F1.GT.99999) THEN
	WRITE(3,1103)H,K,L,F1,DF
        ELSE
        WRITE(3,1104)H,K,L,F1,DF
1104    FORMAT(3I4,2F8.2)
        ENDIF
	GO TO 170
500     CONTINUE
	CLOSE (UNIT=1)
	CLOSE (UNIT=2)
        WRITE(3,*)' '
	CLOSE (UNIT=3)
	ABSM=SUMABS/NUM
	WRITE(6,*)' '
	WRITE(6,*)' '
        WRITE(6,*)'RESULTS FROM THE FULL HKLF DATA:'
	TMIN=TMIN/TMAX
	ABSM=ABSM/TMAX
	TMAX=1.0
        WRITE(6,1016)TMIN
 1016   FORMAT(/,' TRANSMISSION RATIO; Tmin/Tmax =',F6.3)
	WRITE(6,1017)ABSM,NUM
 1017   FORMAT(/' MEAN OF TRANSMISSION COEFFICIENTS =',F6.3,';
     1  TOTAL REFLECTIONS = ',I5/)
	WRITE(6,*)'     *****     ABSACS1  COMPLETED.     *****'
	WRITE(*,*)'     *****     ABSACS1  COMPLETED.     *****'
2000    RETURN
	END
************************************************************************
C23456789012345678901234567890123456789012345678901234567890123456789012
	SUBROUTINE PSIPLOT(NDS,START,CB,TITLE)
	REAL CB(9,100),START(100),TRANS(60)
	CHARACTER*1 P(18,70)
	DIMENSION TITLE(20)
	INTEGER NY(60)
	TMAX=0.0
	TMIN=10.0
	DO 100 L=1,60
	RATIO=0.0
	DO 110 I=1,NDS
	ANTCAL=0.0
	ANTCAL1=0.0
	DO 120 J=1,9
120     ANTCAL=ANTCAL+CB(J,I)*((3.0*L)**(J-1))
C        RATIO=RATIO+START(I)/ANTCAL
        RATIO=RATIO+ANTCAL/START(I)
110     CONTINUE
        TRANS(L)=RATIO/NDS
        TRANS(L)=1.0/TRANS(L)
	IF(TMAX.LT.TRANS(L)) TMAX=TRANS(L)
	IF(TMIN.GT.TRANS(L)) TMIN=TRANS(L)
100     CONTINUE
	DO 200, L=1,60
	TRANS(L)=TRANS(L)/TMAX
200     CONTINUE
	TMIN=TMIN/TMAX
	WRITE(6,201)
201     FORMAT(1X,'THE BELOW PHI & TRANSM CAN BE USED WITH GOOD GRAPHICS
     1 PROGRAM TO GET GOOD PLOT.')
	WRITE(6,202)
202     FORMAT(3X,'PHI'3X,'TRANSM.(%)')
	DO 300 I=1,60
	WRITE(6,299)3*I,100.0*TRANS(I)
299     FORMAT(1X,I4,',  ',F7.3)
	NY(I)=NINT(50.0-50.0*TRANS(I))+1
	IF(NY(I).LT.1) NY(I)=1
	IF(NY(I).GT.18) NY(I)=18
300     CONTINUE
	DO 401 I=1,18
	DO 400 J=1,70
	P(I,J)=' '
400     CONTINUE
	P(I,5)='.'
	P(I,10)='|'
	P(I,4)='0'
401     CONTINUE
	P(5,2)='T'
	P(6,2)='R'
	P(7,2)='A'
	P(8,2)='N'
	P(9,2)='S'
	P(10,2)='M'
	P(11,2)='I'
	P(12,2)='S'
	P(13,2)='S'
	P(14,2)='I'
	P(15,2)='O'
	P(16,2)='N'
	P(4,9)='='
	P(1,4)='1'
	P(1,6)='0'
	P(1,7)='0'
	P(3,6)='9'
	P(3,7)='6'
	P(5,6)='9'
	P(5,7)='2'
	P(7,6)='8'
	P(7,7)='8'
	P(9,6)='8'
	P(9,7)='4'
	P(11,6)='8'
	P(11,7)='0'
	P(13,6)='7'
	P(13,7)='6'
	P(15,6)='7'
	P(15,7)='2'
	P(17,6)='6'
	P(17,7)='8'
C       START PLOTTING!
	WRITE(6,*)' '
	WRITE(6,298)TMIN
298     FORMAT(2X,'Tmin/Tmax =',F7.4'  FROM PSI DATA ALONE')
        IF(TMIN.LE.0.0) THEN
        WRITE(*,*)' ERROR IN MODELLING!!!!!!!!!!!!!!!'
        WRITE(*,*)' ONE OR MORE PSI DATA HAS PROBLEM IN MODELLING.'
        WRITE(*,*)' check the .OUT file for a column with INTCAL and'
        WRITE(*,*)' if anything is unusual for any data set then'
        WRITE(*,*)' PLEASE REMOVE EACH OF THE 37 REFLECTIONS (74 LINES)'
        WRITE(*,*)' BELONGING TO THOSE REFLECTIONS AND RERUN ABSACS!'
        WRITE(*,*)' '
        ELSE
        ENDIF     
	DO 501 I=1,18
	DO 500 L=1,60
	IF(I.EQ.NY(L)) P(I,L+10)='*'
500     CONTINUE
	IF(I.LE.4) P(I,10)=':'
	NODD=2*ABS(I/2)
	IF(NODD.NE.I) THEN
	P(I,9)='-'
	ELSE
	P(I,4)=' '
	P(I,5)=' '
	ENDIF
501     CONTINUE
	WRITE(6,510)(TITLE(I),I=1,10)
510     FORMAT(/15X,'PSIPLOT BY XRAYACS1 FOR ',10A4/)
	WRITE(6,520)((P(I,J),J=1,70),I=1,18)
520     FORMAT(70A1)
	WRITE(6,530)
530     FORMAT(3X,'0.64',1X,'-',6(':---------'),':')
	WRITE(6,540)
540     FORMAT(9X,'0',8X,'30',8X,'60',8X,'90',8X,'120',7X,'150',
     1  7X,'180')
	WRITE(6,550)
550     FORMAT(36X,'ANGLE')
	RETURN
	END
************************************************************************
	SUBROUTINE ABSORB(RESTAR,ND)
	DIMENSION CB(9,100)
	INTEGER H,K,L
	COMMON/ABS1/ WAVE,R11,R12,R13,R21,R22,R23,R31,R32,R33,CB,PSIINT
	COMMON/ABS2/ H,K,L,PHI,PHI1,PHI2,ABSS,ANGKI,TETA,MDEG,RTD
	DTR=1.0/RTD
	XP=H*R11+K*R12+L*R13
	YP=H*R21+K*R22+L*R23
	ZP=H*R31+K*R32+L*R33
	D=SQRT(XP*XP+YP*YP+ZP*ZP)
	STETA=0.5*WAVE*D
	TETA=ASIN(STETA)*RTD
	RP=SQRT(XP*XP+YP*YP)
	ANGKI=ATAN2(ZP,RP)*RTD
	CP=XP/RP
	SP=YP/RP
	PHI=-ATAN2(CP,SP)*RTD
C       WRITE(6,1000)H,K,L,TETA,PHI,ANGKI
 1000   FORMAT(2X,78X,3I4,3F9.2)
	TGAM2=(SIN(TETA*DTR)*COS(ANGKI*DTR))/COS(TETA*DTR)
	GAM2=ATAN(TGAM2)*RTD
	PHI1=PHI+GAM2
	PHI2=PHI-GAM2
	MARK=0
	SUM=0.0
	PSI=PHI1
	GOTO 100
110     MARK=1
	PSI=PHI2
100     CONTINUE
        IF(PSI)120,120,130
120     PSI=PSI+360.0
130     CONTINUE
        IF(PSI.GT.180.0) PSI=PSI-180.0
	F1=0.0
	DO 140 I2=1,MDEG+1
140     F1=F1+CB(I2,ND)*PSI**(I2-1)
C        F1=RESTAR/F1
        F1=F1/RESTAR
        SUM=SUM+F1
	IF(MARK)110,110,150
C150     ABSS=SUM*0.5
150     ABSS=2.0/SUM
	RETURN
	END
************************************************************************
	SUBROUTINE POLY(W,V,M,N,Y)
	DIMENSION W(37),V(37),X(9,9),Y(9)
        REAL ACS
	M1=M+1
	DO 10 I=1,M1
	DO 10 J=1,M1
	X(I,J)=0.0
10      Y(I)=0.0
	DO 20 I=1,M1
	DO 30 J=1,M1
	DO 40 K=1,N
	ACS=(V(K)/10.0)**(I+J-2)
	ACS=ACS*(10.0**(I+J-12))
	X(I,J)=X(I,J)+ACS
40      CONTINUE
        IF(X(I,J).GT.1.7E28) THEN
        WRITE(6,*)'  OVERFLOW'
        X(I,J)=1.7E38
        ELSE
        X(I,J)=X(I,J)*1.0E10
        ENDIF
30      CONTINUE
	DO 20 K=1,N
20      Y(I)=Y(I)+W(K)*V(K)**(I-1)
	WRITE(6,100)((X(I,J),J=1,M1),Y(I),I=1,M1)
 100    FORMAT(2X,10E11.4)
	CALL SIMQ(X,Y,M1,KS)
	RETURN
	END
************************************************************************
	SUBROUTINE SIMQ(A,B,N,KS)
	REAL A(1),B(1)
C       FORWARD SOLUTION
	TOL=0.0
	KS=0
	JJ=-N
	DO 10 J=1,N
	JY=J+1
	JJ=JJ+N+1
	BIGA=0
	IT=JJ-J
	DO 20 I=J,N
C       SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN
	IJ=IT+I
	IF(ABS(BIGA)-ABS(A(IJ)))30,20,20
30      BIGA=A(IJ)
	IMAX=I
20      CONTINUE
C       TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX)
	IF(ABS(BIGA)-TOL)35,35,40
35      KS=1
	RETURN
C       INTERCHANGE ROWS IF NECESSARY
40      I1=J+N*(J-2)
	IT=IMAX-J
	DO 50 K=J,N
	I1=I1+N
	I2=I1+IT
	SAVE=A(I1)
	A(I1)=A(I2)
	A(I2)=SAVE
C       DIVIDE EQUATION BY LEADING COEFFICIENT
50      A(I1)=A(I1)/BIGA
	SAVE=B(IMAX)
	B(IMAX)=B(J)
	B(J)=SAVE/BIGA
C       ELIMINATE NEXT VARIABLE
	IF(J-N)55,60,55
55      IQS=N*(J-1)
	DO 10 IX=JY,N
	IXJ=IQS+IX
	IT=J-IX
	DO 70 JX=JY,N
	IXJX=N*(JX-1)+IX
	JJX=IXJX+IT
70      A(IXJX)=A(IXJX)-(A(IXJ)*A(JJX))
10      B(IX)=B(IX)-(B(J)*A(IXJ))
C       BACK SOLUTION
60      NY=N-1
	IT=N*N
	DO 80 J=1,NY
	IA=IT-J
	IB=N-J
	IC=N
	DO 80 K=1,J
	B(IB)=B(IB)-A(IA)*B(IC)
	IA=IA-N
80      IC=IC-1
	RETURN
	END
************************************************************************
************************************************************************
*                       BESTABS.FOR                                    *
************************************************************************
* ABSORPTION CORRECTION FOLLOWING THE METHOD OF N.WALKER AND S.STUART, *
*       ACTA CRYST., A39, 158-166 (1983).                              *
************************************************************************
* PROGRAM BY    F.UGOZZOLI, COMPUT. CHEM., 11, 109-120 (1987).         *
************************************************************************
	SUBROUTINE BESTABS
	INTEGER CR,SHX,SHX1,LP,NEW
	CHARACTER*1 TITLE(40)
	REAL*8 UB(3,3)
	DIMENSION MN(2,34)
	REAL*8 PQ(34,34),TN(34),RS(34,34),TN1(34)
	REAL*8 SCALA,OME,FOQ,SCMN(34),SCM(12),SMN,CMN,TET
	REAL*8 SUMFO,SUMFM,SUMFC,DELFOFC,DELFMFC,RINIZ,RFINAL
	CHARACTER*12 OUT1,OUT2,OUT3
	CHARACTER*12 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OPT/ OPT1,OPT2,OPT3,OPT4,NTY1,NTY2,NTY3
	COMMON/FILE/ FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7,FILE8
	COMMON/OUT/ OUT1,OUT2,OUT3
	DATA ((MN(I,J),J=1,34),I=1,2)/
     9 0,0,0,0,0,0,1,2,3,4,1,1,1,1,2,3,4,
     9 0,0,0,0,0,0,1,2,3,4,1,1,1,1,2,3,4,
     9 1,2,3,4,5,6,0,0,0,0,1,2,3,4,2,3,4,
     9 1,2,3,4,5,6,0,0,0,0,1,2,3,4,2,3,4/
C
	IF(OPT4.EQ.0) RETURN
	WRITE(6,*)' '
	WRITE(6,*)' *****     DOING    BESTABS.....     *****'
	WRITE(*,*)' '
	WRITE(*,*)' *****     DOING    BESTABS.....     *****'
	WRITE(6,1)NTY3
	WRITE(*,1)NTY3
 1      FORMAT('                                ON F^',I1)
	CR=1
	SHX=2
	SHX1=3
	LP=6
	NEW=4
	OPEN(UNIT=CR,FILE=FILE6,STATUS='OLD',FORM='FORMATTED')
C
	WRITE(LP,1000)
	READ(CR,1010)TITLE
	WRITE(LP,1020)TITLE
	WRITE(LP,1030)
13      CONTINUE
	READ(CR,*)WAVE
	WRITE(LP,1050)WAVE
	READ(CR,*)((UB(I,J),J=1,3),I=1,3)
	WRITE(LP,1052)((UB(I,J),J=1,3),I=1,3)
	READ(CR,*,ERR=11)NCICLI
11      IF(NCICLI.EQ.0) NCICLI=2
	WRITE(6,*)'  '
	WRITE(6,*)' REQUESTED NO. OF CYCLES= ',NCICLI
	ICICLO=1
	NUM=0
	OPEN(UNIT=SHX,FILE=FILE7,STATUS='OLD',ERR=40,FORM='FORMATTED')
	OPEN(UNIT=SHX1,FILE=FILE8,STATUS='UNKNOWN',FORM='FORMATTED')
	CONTINUE
	IF(NTY3.EQ.1) GOTO 37
35      READ(SHX,102,ERR=40,END=50)IH,IK,IL,FC,FO,SIGFO
	IF(IH.EQ.0.AND.IK.EQ.0.AND.IL.EQ.0) GO TO 35
        WRITE(SHX1,102)IH,IK,IL,FO,SIGFO,FC
	IF(FO.LT.3.0*SIGFO) GO TO 35
	NUM=NUM+1
	SUMFO=SUMFO+FO
	SUMFC=SUMFC+FC
	GO TO 35
37      READ(SHX,103,ERR=40,END=50)IH,IK,IL,FO,FC,SIGFO
	IF(IH.EQ.0.AND.IK.EQ.0.AND.IL.EQ.0) GO TO 37
        WRITE(SHX1,102)IH,IK,IL,FO,SIGFO,FC
	IF(FO.LT.3.0*SIGFO) GO TO 37
	NUM=NUM+1
	SUMFO=SUMFO+FO
	SUMFC=SUMFC+FC
	GO TO 37
40      SHX=SHX1
	OPEN(UNIT=SHX1,FILE=FILE8,STATUS='OLD',FORM='FORMATTED')
45      READ(SHX,102,END=50)IH,IK,IL,FO,SIGFO,FC
	IF(FO.LT.3.0*SIGFO) GO TO 45
	NUM=NUM+1
	SUMFO=SUMFO+FO
	SUMFC=SUMFC+FC
	GO TO 45
50      CONTINUE
	ENDFILE SHX1
C       evaluate scale factor
2000    WRITE(LP,1200)ICICLO
	WRITE(*,*)' CURRENT CYCLE ',ICICLO
	SCALA=SUMFC/SUMFO
	WRITE(LP,1080)NUM,SCALA
C
	APSMIN=1.0E10
	APSMAX=-1.0E10
	ATMIN=1.0E10
	ATMAX=-1.0E10
	TMAX=-1.0E10
	TMIN=1.0E10
C
C       build the 34x34 linear system in polar angles PHI and MU
	DO 100 I=1,34
	TN(I)=0.0
	DO 100 J=1,34
100     PQ(I,J)=0.0
	NRIF=0
	REWIND SHX1
110     READ(SHX1,102,END=350)IH,IK,IL,FO,SIGFO,FC
	IF(FO.LT.3.0*SIGFO) GO TO 110
	FCC=FC
	IF(FCC.LT.3.0*SCALA*SIGFO.OR.FCC.GT.2.0*SCALA*FO) GO TO 110
	NRIF=NRIF+1
C
	CALL ANG(IH,IK,IL,UB,WAVE,TET,CHI,PHI)
	CALL PHIMU(WAVE,TET,CHI,PHI,PHIP,PHIS,U)
	OME=1.0/(SCALA*SIGFO)
	FOQ=FO*FO*OME
	DO 150 I=1,17
	CALL SINMN(MN(1,I),MN(2,I),PHIP,PHIS,U,SCMN(I))
	CALL COSMN(MN(1,I+17),MN(2,I+17),PHIP,PHIS,U,SCMN(I+17))
150     CONTINUE
	DO 200 I=1,17
	DO 200 J=1,34
	PQ(I,J)=PQ(I,J)+FOQ*SCMN(J)*SCMN(I)
	PQ(I+17,J)=PQ(I+17,J)+FOQ*SCMN(J)*SCMN(I+17)
200     CONTINUE
	DO 210 N=1,34
	TN(N)=TN(N)+FO*OME*SCMN(N)*(FC/SCALA-FO)
210     CONTINUE
	GO TO 110
350     CONTINUE
	WRITE(LP,1072)NRIF
C       solve the 34x34 linear system in theta
	CALL GAUSS(PQ,TN,34)
	WRITE(LP,1090)(MN(1,I),MN(2,I),TN(I),I=1,34)
C       build the 12x12 linear system in theta
	DO 360 I=1,12
	TN1(I)=0.0
	DO 360 J=1,12
360     RS(I,J)=0.0
	IRIF=0
	REWIND SHX1
370     READ(SHX1,102,END=600)IH,IK,IL,FO,SIGFO,FC
	IF(FO.LT.3.0*SIGFO) GO TO 370
	FCC=FC
	IF(FCC.LT.3.0*SCALA*SIGFO.OR.FCC.GT.2.0*SCALA*FO) GO TO 370
	IRIF=IRIF+1
	CALL ANG(IH,IK,IL,UB,WAVE,TET,CHI,PHI)
	CALL PHIMU(WAVE,TET,CHI,PHI,PHIP,PHIS,U)
	APS=1.0
	DO 400 I=1,17
	CALL SINMN(MN(1,I),MN(2,I),PHIP,PHIS,U,SMN)
	CALL COSMN(MN(1,I+17),MN(2,I+17),PHIP,PHIS,U,CMN)
	APS=APS+SMN*TN(I)+CMN*TN(I+17)
400     CONTINUE
	IF(APS.LT.0.0) CALL CORR(1,IH,IK,IL,TET,CHI,PHI,FO,APS)
	FM=SCALA*FO*APS
	OME=1.0/(SCALA*SIGFO)
	FMM=FM*FM*OME
	DO 450 M=1,6
	SCM(M)=SIN(M*TET)
	SCM(M+6)=COS(M*TET)
450     CONTINUE
	DO 500 I=1,6
	DO 500 J=1,12
	RS(I,J)=RS(I,J)+FMM*SCM(J)*SCM(I)
	RS(I+6,J)=RS(I+6,J)+FMM*SCM(J)*SCM(I+6)
500     CONTINUE
	DO 510 I=1,12
	TN1(I)=TN1(I)+FM*OME*SCM(I)*(FC-FM)
510     CONTINUE
	GO TO 370
600     CONTINUE
	WRITE(LP,1073)IRIF
C solve the 12x12 linear system
	CALL GAUSS(RS,TN1,12)
	WRITE(LP,1100)
	DO 615 I=1,12
	JJJ=I
	IF(JJJ.GT.6) JJJ=JJJ-6
615     WRITE(LP,1110)JJJ,TN1(I)
C  new cycle of absorption correction
610     CONTINUE
	REWIND SHX1
	OPEN(UNIT=NEW,FILE=OUT3,STATUS='UNKNOWN',FORM='FORMATTED')
611     READ(SHX1,102,END=700)IH,IK,IL,FO,SIGFO,FC
	IF(IH.EQ.0.AND.IK.EQ.0.AND.IL.EQ.0) GO TO 700
	SIGMA=SIGFO
	CALL ANG(IH,IK,IL,UB,WAVE,TET,CHI,PHI)
	CALL PHIMU(WAVE,TET,CHI,PHI,PHIP,PHIS,U)
	APS=1.0
	DO 620 I=1,17
	CALL SINMN(MN(1,I),MN(2,I),PHIP,PHIS,U,SMN)
	CALL COSMN(MN(1,I+17),MN(2,I+17),PHIP,PHIS,U,CMN)
	APS=APS+SMN*TN(I)+CMN*TN(I+17)
620     CONTINUE
	IF(APS.LT.0.0) CALL CORR(1,IH,IK,IL,TET,CHI,PHI,FO,APS)
	FM=SCALA*FO*APS
	ATET=1.0
	DO 650 M=1,6
	SCM(M)=SIN(M*TET)
	SCM(M+6)=COS(M*TET)
	ATET=ATET+SCM(M)*TN1(M)+SCM(M+6)*TN1(M+6)
650     CONTINUE
	IF(ATET.LT.0.0) CALL CORR(2,IH,IK,IL,TET,CHI,PHI,FO,ATET)
	FM=FM*ATET
	SIGFO=SIGFO*SCALA*APS*ATET
        IF(FM.LT.99999) THEN
	WRITE(NEW,101)IH,IK,IL,FM,SIGFO,FC
        ELSE
        WRITE(NEW,1011)IH,IK,IL,FM,SIGFO,FC
        ENDIF
	IF(FO.LT.3.0*SIGMA) GO TO 611
	IF(APS.GT.APSMAX) APSMAX=APS
	IF(APS.LT.APSMIN) APSMIN=APS
	IF(ATET.GT.ATMAX) ATMAX=ATET
	IF(ATET.LT.ATMIN) ATMIN=ATET
	TT=APS*ATET
	TT=TT**(2.0/NTY3)
	IF(TT.GT.TMAX) TMAX=TT
	IF(TT.LT.TMIN) TMIN=TT
C
	GO TO 611
700     CONTINUE
	IF(ICICLO.NE.NCICLI) GO TO 705
705     CONTINUE
	NNN=0
	SUMFC=0.0
	SUMFM=0.0
	DELFOFC=0.0
	DELFMFC=0.0
	ENDFILE NEW
	REWIND NEW
	REWIND SHX1
850     READ(SHX1,102,END=800)IH,IK,IL,FO,SIGFO,FC
851     READ(NEW,1011,END=970)IHN,IKN,ILN,FM,SIGFO,FC
	IF(IH.NE.IHN.OR.IK.NE.IKN.OR.IL.NE.ILN) GO TO 970
	IF(FO.LT.3.0*SIGFO) GO TO 850
	IF(FM.LT.3.0*SCALA*SIGFO) GO TO 850
	NNN=NNN+1
	DELFOFC=DELFOFC+ABS(FO*SCALA-FC)
	SUMFO=SUMFO+FO*SCALA
	DELFMFC=DELFMFC+ABS(FM-FC)
	SUMFM=SUMFM+FM
	SUMFC=SUMFC+FC
	GO TO 850
800     RINIZ=DELFOFC/SUMFO
	RFINAL=DELFMFC/SUMFM
	WRITE(LP,1130)NNN,RINIZ,RFINAL
C
	WRITE(LP,1140)APSMAX,APSMIN,ATMAX,ATMIN
	WRITE(LP,1150)TMAX,TMIN
C   set in new cycle of least squares absorption correction
	IF(ICICLO.EQ.NCICLI) GO TO 2999
	REWIND NEW
	REWIND SHX1
990     READ(NEW,1011,END=980)IH,IK,IL,FM,SIGFO,FC
        WRITE(SHX1,102)IH,IK,IL,FM,SIGFO,FC
	GO TO 990
980     REWIND SHX1
	REWIND NEW
	SUMFC=SUMFC
	SUMFO=SUMFM
	ICICLO=ICICLO+1
	GO TO 2000
970     WRITE(LP,1210)
	STOP 017
C***********************************************************************
C                               FORMATS
C***********************************************************************
101     FORMAT(3I4,F8.2,F8.2,F8.1)
1011    FORMAT(3I4,F8.1,F8.2,F8.1)
102     FORMAT(3I4,2F12.2,F10.2)
103     FORMAT(3I4,2F8.1,32X,F8.2)
104     FORMAT(A1)
105     FORMAT(6A1,1X,6A1)
1000    FORMAT(/1X,20(1H.),'PROGRAM   ABSORB',20(1H.))
1010    FORMAT(40A1)
1020    FORMAT(/1X,'TITLE OF THE JOB: ',40A1)
1030    FORMAT(/1X,'FOUR CIRCLE DIFFRACTION GEOMETRY')
1050    FORMAT(1X,'WAVELENGTH=',F8.6)
1052    FORMAT(//1X,'UB(3,3) ORIENTATION MATRIX',/,(1X,3F15.7))
1053    FORMAT(//1X,'MILLER INDICES TRANSFORMATION MATRIX ROT(3,3) ',
     1  'PREVIOUSLY APPLIED BEFORE SHELX CALCULATIONS.',(/1X,3F8.2))
1054    FORMAT(//1X,'REVERSE MATRIX (ROT)**-1',(/1X,3F8.3))
1055    FORMAT(//1X,'COORECT ORIENTATION MATRIX CALCULATED AS:',
     1  /1X,'(UBP)=(UB)*(ROT)**-1'
     2  /1X,'TO CALCULATE FROM TRANSFORMED MILLER INDICES',
     3  /1X,'THE CORRECT DIFFRACTOMETRIC THETA CHI PHI ANGLES',
     4  (/1X,3F15.7))
1060    FORMAT(3I4,2F8.2,I4,6F8.5)
1072    FORMAT(///,'NUMBER OF REFLECTIONS USED TO CREATE',1X,
     1  'THE LINEAR SYSTEM OF ORDER 34 IN THE POLAR ANGLES',1X,
     2  'PHI AND MU:',I5)
1073    FORMAT(1H1,'NUMBER OF REFLECTIONS USED TO CREATE',1X,
     1  'THE LINEAR SYSTEM OF ORDER 12 IN THE SCATTERING ANGLE',1X,
     2  'THETA:',I5)
1080    FORMAT(/,' SCALE FACTOR (EVALUATED ON',I5,' REFLECTIONS)=',F9.6)
1090    FORMAT(//1X,'COEFFICIENTS FOR FOURIER SERIES EXPANSION',
     1  /1X,'IN PHI AND MU OF THE ABSORPTION CORRECTIONS:',
     2  /5X,'M',5X,'N',11X,'PQ',/,(1X,I5,1X,I5,F15.6))
1100    FORMAT(/1X,'COEFFICIENTS FOR FOURIER SERIES EXPANSION',
     1  /1X,'IN THETA OF THE ABSORPTION CORRECTIONS:',
     2  /5X,'N',8X,'RS')
1110    FORMAT(1X,I5,F15.6)
1130    FORMAT(1H1,'TEST OF FO-FC AGREEMENT ON',I5,' REFLECTIONS:',
     2  //1X,'BEFORE FO''S ABSORPTION CORRECTIONS',
     3  5X,'R=SUM(ABS(K*F-ABS(FC)))/SUM(FO)=',F15.6
     4  //1X,'AFTER FO',1H','S ABSORPTION CORRECTIONS',
     5  5X,'R=SUM(ABS(FM-ABS(FC)))/SUM(FM)=',F15.6
     6  /1X,'WHERE K IS THE SCALE FACTOR GIVEN ABOVE, AND ',
     7  'FM ARE THE FO',1H','S CORRECTED FOR ABSORPTION.')
1140    FORMAT(////1X,'MAX. AND MIN. VALUES FOR ABSORPTION CORRECTIONS',
     1  /1X,'EVALUATED ON THE REFLECTIONS WHICH SATISFY:',
     1  ' FO> 3.0 *SIGMA(FO)',/1X,
     2  'CORRECTION IN PHI AND MU',10X,'MAX=',F12.6,10X,'MIN=',F12.6,
     3  /1X,'CORRECTION IN THETA',13X,'MAX=',F12.6,10X,'MIN=',F12.6)
1150    FORMAT(//2X,'OVERALL TRANSMISSION; MAX= ',F6.3,'  MIN= ',F6.3)
1200    FORMAT(/,1X,'CYCLE',I3)
1210    FORMAT(1X,'SEQUENCE ERROR: PLEASE RERUN WITH NCYCLE=1')
2999    WRITE(6,*)'     *****     BESTABS  COMPLETED.    *****'
	WRITE(*,*)'     *****     BESTABS  COMPLETED.    *****'
	RETURN
	END
C***********************************************************************
	SUBROUTINE ANG(IH,IK,IL,UB,WAVE,TET,CHI,PHI)
	REAL*8 UB(3,3),TET
	REAL NORM
C       IGEOM=1
	LP=13
	PIG=3.14159265
	X=UB(1,1)*IH+UB(1,2)*IK+UB(1,3)*IL
	Y=UB(2,1)*IH+UB(2,2)*IK+UB(2,3)*IL
	Z=UB(3,1)*IH+UB(3,2)*IK+UB(3,3)*IL
	RAD=SQRT(X*X+Y*Y+Z*Z)
	NORM=SQRT(X*X+Y*Y)
C       GO TO (100,200),IGEOM
100     SNCHI=-Z/RAD
	SENT=RAD/2.0
	IF(SENT.LT.1.0) GO TO 110
	SENT=1.0
	WRITE(LP,50)IH,IK,IL
110     CHI=ATAN2(SNCHI,SQRT(1.0-SNCHI*SNCHI))
	IF(CHI.GT.PIG) CHI=CHI-2.0*PIG
	TET=ATAN2(SENT,SQRT(1.0-SENT*SENT))
	SNPHI=X/NORM
	COSPHI=Y/NORM
	PHI=ATAN2(SNPHI,COSPHI)
	RETURN
200     SENT=RAD*WAVE/2.0
	IF(SENT.LT.1.0) GO TO 210
	SENT=1.0
	WRITE(LP,50)IH,IK,IL
210     TET=ATAN2(SENT,SQRT(1.0-SENT*SENT))
	SNCHI=Z/RAD
	COSCHI=NORM/RAD
	CHI=ATAN2(SNCHI,COSCHI)
	SNPHI=Y/NORM
	COSPHI=X/NORM
	PHI=ATAN2(SNPHI,COSPHI)
	IF((CHI*57.29577951).GT.-5.0) RETURN
	CHI=-CHI
	PHI=PHI+PIG
	IF(PHI.GE.2.0*PIG) PHI=PHI-(2.0*PIG)
	RETURN
50      FORMAT(1X,'CONFLICTING BETWEEN MILLER INDICES',3I6,
     1  /1X,'AND UB ORIENTATION MATRIX')
	END
C***********************************************************************
	SUBROUTINE PHIMU(WAVE,TET,CHI,PHI,PHIP,PHIS,U)
	REAL*8 TET
	PIG=3.14159265
	ARG=SIN(TET)*SIN(CHI)
	UU=ATAN2(ARG,(1.0-ARG*ARG))
	U=SQRT(ABS(UU))
	IF(UU.LT.0) U=-U
	U=U*7.0*WAVE
	TG=TAN(TET)
	CS=COS(CHI)
	ARG=1.0/SQRT(1.0+TG*TG*CS*CS)
	DPHI=ATAN2(SQRT(1.0-ARG*ARG),ARG)
	PHIP=PHI+DPHI
	PHIS=PHI-DPHI+PIG
	RETURN
	END
C***********************************************************************
	SUBROUTINE SINMN(M,N,PHIP,PHIS,U,SMN)
	REAL*8 SMN
	SMN=SIN(M*U+N*PHIP)+SIN(M*U+N*PHIS)
	RETURN
	END
C***********************************************************************
	SUBROUTINE COSMN(M,N,PHIP,PHIS,U,CMN)
	REAL*8 CMN
	CMN=COS(M*U+N*PHIP)+COS(M*U+N*PHIS)
	RETURN
	END
C***********************************************************************
	SUBROUTINE GAUSS(A,B,N)
	DOUBLE PRECISION A,B,RP,PV,T,RIS
	DIMENSION A(34,34),B(34),RP(35)
	T=0.0
	NP1=N+1
	NM1=N-1
	DO 10 L=1,N
	PV=0.0
	LM=N-L+1
	DO 1 K=1,LM
	IF(DABS(A(K,1)).LT.PV) GO TO 1
	PV=DABS(A(K,1))
	KK=K
1       CONTINUE
	IF(PV.GT.T) GO TO 2
	RETURN
2       DO 3 J=1,N
	RIS=A(KK,J)
	A(KK,J)=A(1,J)
	A(1,J)=RIS
3       CONTINUE
	RIS=B(KK)
	B(KK)=B(1)
	B(1)=RIS
	DO 4 J=1,NM1
4       RP(J)=A(1,J+1)/A(1,1)
	RP(N)=B(1)/A(1,1)
	RP(NP1)=1.0/A(1,1)
	DO 6 I=2,N
	RIS=A(I,2)-A(I,1)*RP(1)
	DO 5 J=2,NM1
5       A(I,J)=A(I,J+1)-A(I,1)*RP(J)
	A(I,N)=B(I)-A(I,1)*RP(N)
	B(I)=-A(I,1)*RP(NP1)
	A(I,1)=RIS
6       CONTINUE
	DO 8 J=1,N
	DO 7 I=1,NM1
7       A(I,J)=A(I+1,J)
	A(N,J)=RP(J)
8       CONTINUE
	DO 9 I=1,NM1
9       B(I)=B(I+1)
	B(N)=RP(NP1)
10      CONTINUE
	DO 12 I=1,N
	RIS=A(I,1)
	DO 11 J=1,NM1
11      A(I,J)=A(I,J+1)
	A(I,N)=B(I)
	B(I)=RIS
12      CONTINUE
	RETURN
	END
C***********************************************************************
	SUBROUTINE CORR(N,IH,IK,IL,TET,CHI,PHI,FO,A)
	CHARACTER*4 TEXT
      REAL*8 TET
	LP=13
	RAD=57.29577951
	T=TET*RAD
	C=CHI*RAD
	P=PHI*RAD
	TEXT='APS '
	IF(N.EQ.2) TEXT='TET '
	WRITE(LP,100)IH,IK,IL,T,C,P,FO,TEXT,A
100     FORMAT(1X,'H',I3,3X,'K',I3,3X,'L',I3,3X,'THETA',F7.3,3X,
     1  'CHI',F7.3,3X,'PHI ',F7.3,3X,'FO',F7.2,3X,A4,F10.6)
	A=-A
	RETURN
	END
************************************************************************
	SUBROUTINE EXTRAS
	CHARACTER*1 OPT
1       WRITE(*,*)' DO YOU WANT TO CONVERT F^1 TO F^2? (Y/N)?'
	READ(*,'(A1)')OPT
        IF(OPT.NE.'Y'.AND.OPT.NE.'N'.AND.OPT.NE.'y'.AND.OPT.NE.'n')
     1 GOTO 1
        IF (OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	CALL F1TOF2
	ELSE
2       WRITE(*,*)' DO YOU WANT TO CONVERT F^2 TO F^1? (Y/N)?'
	READ(*,'(A1)')OPT
        IF(OPT.NE.'Y'.AND.OPT.NE.'N'.AND.OPT.EQ.'y'.AND.OPT.EQ.'n')
     1 GOTO 2
        IF (OPT.EQ.'Y'.OR.OPT.EQ.'y') THEN
	CALL F2TOF1
	ELSE
	ENDIF
	ENDIF
	STOP  '     *****     NO JOB HAS BEEN REQUESTED!!!!!    *****'
	END
************************************************************************
	SUBROUTINE F1TOF2
	CHARACTER*12 FILE1,FILE2
	INTEGER H,K,L
	WRITE(*,*)' CAUTION: USE ONLY TO CONVERT THE DATA REDUCED BY XRA
     1YACS PROGRAM!'
	WRITE(*,*)' SOME PEOPLE MAY NOT AGREE WITH THIS. BUT YOUR RESULT
     1 WILL BE BETTER ONLY!!'
	WRITE(*,*)' IT IS JUST EQUAL TO MAKE THE -VE F^2s ZERO, WHICH GI
     1VES BETTER RESULTS!!!'
	WRITE(*,*)' LAST DIGIT ACCURACY IS LOST!! BUT IT SHOULD BE OK.'
	WRITE(*,*)' '
	WRITE(*,*)' PLEASE TYPE INPUT FILE NAME (12 LETTERS):'
	READ(*,'(A12)')FILE1
	OPEN(UNIT=1,FILE=FILE1,STATUS='OLD',FORM='FORMATTED')
	WRITE(*,*)' PLEASE TYPE OUTPUT FILE NAME (12 LETTERS):'
	READ(*,'(A12)')FILE2
	OPEN(UNIT=2,FILE=FILE2,STATUS='NEW',FORM='FORMATTED')
1       READ(1,10,END=2)H,K,L,F,S
10      FORMAT(3I4,F8.1,F8.2)
	FSQ=F*F
	IF(F.LT.0.25) THEN
	S=S
	ELSE
	S=2.0*S*F
	ENDIF
	WRITE(2,10)H,K,L,FSQ,S
	GOTO 1
2       STOP ' CONVERSION F^1 TO F^2 DONE!.'
	END
************************************************************************
	SUBROUTINE F2TOF1
	CHARACTER*8 FILE1,FILE2
	INTEGER H,K,L
	WRITE(*,*)' PLEASE TYPE INPUT FILE NAME (12 LETTERS):'
	READ(*,'(A12)')FILE1
	OPEN(UNIT=1,FILE=FILE1,STATUS='OLD',FORM='FORMATTED')
	WRITE(*,*)' PLEASE TYPE OUTPUT FILE NAME (12 LETTERS):'
	READ(*,'(A12)')FILE2
	OPEN(UNIT=2,FILE=FILE2,STATUS='NEW',FORM='FORMATTED')
1       READ(1,10,END=2)H,K,L,F,S
10      FORMAT(3I4,F8.1,F8.2)
	IF(F.LT.0.0) F=0.0
	F=SQRT(F)
	IF(F.LT.0.5) THEN
	S=S
	ELSE
	S=S/(2.0*F)
	ENDIF
	WRITE(2,10)H,K,L,F,S
	GOTO 1
2       STOP ' CONVERSION F^2 TO F^1 DONE!.'
	END
************************************************************************