    	PROGRAM INTCAL
C
C **	PROGRAM TO CORRECT FOR INTERNAL CALIBRATION OF ADR FILES
C
C	Subroutines are in SUBS.FOR and INPUT.FOR
C
C **	PROGRAMMER ROBERT L. SNYDER
C		   ALFRED UNIVERSITY
C		   ALFRED, N.Y. 14802
C
C	Modifications by Robert Snyder to:
C		Allow rejecting bad peaks interactively
C		New autoplot file format
C		Allow diffrac 11 or adr output formats
C
C       With modifications by Jim Walker for:
C		Diffrac-11 format file format (subroutine INPUT)
C		Standards in a file (subroutine STDFIL)
C		Changed order of a few questions
C		End of file in dialog takes you back to ENTER FILENAME
C		Loop around whole program to do multiple files
C
	IMPLICIT REAL*8(A-H,O-Z)
	CHARACTER*40 NFILE,IJUNK,LFILE,PLTFIL,ISTD
	CHARACTER*1 IDOT, IBL, NAST, IDN, RUNTYP, IPRIR, IOUT
	BYTE TITLE(80), ITIT(80), IDAT(9)
       	REAL*8 XTHE(200),DTHE(200),YINT,COEF(5),REGR,SD
	DIMENSION IFIND(50)
        DIMENSION T1(200),D1(200),I1(201),I2(201),I3(201), IREJ(5)
	DIMENSION DINT(200),LOC(50),IPEN(20),ISYM(20),ILIN(20),ISET(20)
	REAL*4  PLX(202),PLY(202), F4, CODE, D, AI, FWHM
	REAL*4 THETA1(200), DSP1(200),RELINT(200), INTCPS(200), FWHM1(200)
       INTEGER*2 INT2(2)
	COMMON /INP/ ALO,AHI,TIME,STEP,TITLE
	COMMON /IO/ IO5,IO6,IO7,IO10,IO11,IO20
       EQUIVALENCE (F4, INT2(1))
       DATA CODE/'CALI'/
	DATA IO5,IO6,IO7,IO10,IO11,IO20,IO30/5,6,6,10,11,20,30/
	DATA ISI,IDOT,IBL,IBG/'SI','.',' ','BG'/
C
C  START CODE
C
	DTR = 3.141592654 / 180.
	WRITE(IO7,1000)
1000	FORMAT(//10X'**** Internal Standard Calibration of d-I', 
     1 ' Data ****'//15X'*** Default options are in < > ***'/
     & 14X'*** Press return to take default ***'/)
C
C  Get the file name
C
1	WRITE(IO7,1005)
1005	FORMAT('$Name of the d-I file <return = end session>: ')
	READ(IO5,1010,END=999) NFILE
    	IF(NFILE .EQ. ' ') GO TO 999
1010	FORMAT(A40)
C
C    ****  GET THE INTERNAL STANDARD.
C
	GOTO 1019
C  Error: no such standard.
1011	WRITE(IO7,1012) 
1012	FORMAT(' No such standard.  Try again or ^Z to exit.')
C  Read the file name.
1019	WRITE(IO7,1020)
1020	FORMAT(/' Internal standard used: (PRESS CAPS LOCK)'/
     1' SI - Silicon'/
     1' AG - Silver'/' BI - Bi12GeO20'/' BS - Bi12GeO20 + Si'/
     2' NI - Ni Metal'/' NO - NiO'/' FL - Fluorophlogopite'
     2/' FS - Fluorophlogopite + Si'
     2/' CU - Copper'/' FC - Fluorophlogopite + Cu'/
     3' SB - Silver Behenate'/
     4' SS - Silver Behenate + Si'/
     5/'$Default <SI>): ')
	READ(IO5,1030,END=1) ISTD
1030	FORMAT(A2)
C
C **	FILL DINT ARRAY WITH INTERNAL STANDARD TWO-THETA VALUES
C
	NUM = 0
	IF(ISTD.EQ.'  ') CALL SI(DINT,NUM)
	IF(ISTD.EQ.'^Z') GOTO 998
	IF(ISTD.EQ.'SI') CALL SI(DINT,NUM)
	IF(ISTD.EQ.'AG') CALL AG(DINT,NUM)
	IF(ISTD.EQ.'BI') CALL BI(DINT,NUM)
	IF(ISTD.EQ.'BS') CALL BS(DINT,NUM)
	IF(ISTD.EQ.'NI') CALL NI(DINT,NUM)
	IF(ISTD.EQ.'NO') CALL NO(DINT,NUM)
	IF(ISTD.EQ.'FL') CALL FL(DINT,NUM)
	IF(ISTD.EQ.'FS') CALL FS(DINT,NUM)
	IF(ISTD.EQ.'FC') CALL FC(DINT,NUM)
     	IF(ISTD.EQ.'SB') CALL SB(DINT,NUM)
    	IF(ISTD.EQ.'SS') CALL SS(DINT,NUM)
	IF(ISTD.EQ.'CU') CALL CU(DINT,NUM)
	IF (NUM.EQ.0) CALL STDFIL(DINT,NUM,ISTD,IO7,IO11,200)
	IF (NUM.EQ.0) GOTO 1011
C
C
	WRITE(IO7,1040)
1040	FORMAT(/'$Window in 2-theta to look for int. std. peaks <0.50>: ')
	READ(IO5,1050,END=1) EWIND
1050	FORMAT(F10.4)
	IF(EWIND.EQ.0.0)EWIND = 0.50
C
C **	FIND THE DEFAULT FILE NAME
C
	ITEMP=0
	DO 10 I=1,40
	IF(NFILE(I:I).EQ.IDOT) ITEMP=I
10	CONTINUE
	IF (ITEMP.NE.0) GOTO 18
	DO 12 I=40,1,-1
	IF(NFILE(I:I).EQ.IBL) ITEMP=I
12	CONTINUE
	NFILE(ITEMP:ITEMP+3)='.DIF'
18	LFILE = NFILE
	LFILE(ITEMP+1:ITEMP+3)='CAL'
C
	WRITE(IO7,1140)
1140	FORMAT(/'$Print results on TT, File or Both? (T,F OR <B>) ')
	READ(IO5,1150,END=1) IPRIR
1150	FORMAT(A1)
	IF(IPRIR.EQ.IBL .OR. IPRIR.EQ.'b') IPRIR = 'B'
	IF(IPRIR.EQ.'t') IPRIR = 'T'
	IF(IPRIR.EQ.'f') IPRIR = 'F'
C
	IF(IPRIR.EQ.'T') GOTO 20
	WRITE(IO7,1060) LFILE
1060	FORMAT(/' Name of the corrected data output file '/
    	1'$(Default <',A40,'>) : ')
	READ(IO5,1010,END=1) IJUNK
	IF(IJUNK(1:1).EQ.IBL) GO TO 19
	LFILE = IJUNK
C
C  OPEN THE OUTPUT FILE
C
*19       WRITE(IO6,1062)
*1062	FORMAT(/' Select the output file format desired:'/' A - ADR',
*     1' ASCII file'/' D - Diffrac 11 binary file'/' <default = A>: '$)
*       READ(IO5,1064)IOUT
*1064	FORMAT(A1)
*       IF(IOUT .EQ. ' ')IOUT = 'A'
*       IF(IOUT .EQ. 'D' .OR. IOUT .EQ. 'd') THEN
19		OPEN ( UNIT=IO20, NAME=LFILE, TYPE='NEW', 
     +	  	FORM='UNFORMATTED', ACCESS='DIRECT', 
     +    	RECORDSIZE=1, ASSOCIATEVARIABLE=MEXT)
       		IOUT = 'D'
*       ELSE  
*		OPEN(UNIT=IO20,NAME=LFILE,TYPE='NEW',
*     1 		ACCESS='SEQUENTIAL',CARRIAGECONTROL='LIST')
*       ENDIF
C
C
C	OPEN D & I FILE FOR READ (AND WRITE DIFFRAC 11 HEADER INFO)
C
20	CALL INPUT(NFILE,TITLE,ALO,AHI,STEP,TIME,WAVEL,
     &     THETA1,DSP1,RELINT,INTCPS,FWHM1,NPEAKS,IOUT)
C
80	INUM=0
C
	WRITE(IO7,1160)
1160	FORMAT(/'$Do you want hkl indices put on output file? <N> ')
	READ(IO5,1150,END=1) IDN
	INDEX = 0
	IF(IDN.EQ.'Y'.OR.IDN.EQ.'y')INDEX = 1
C  JDW put this undocumented feature in to debug STDFIL.
	IF(IDN.EQ.'S') WRITE(IO7,1165) NUM,ISTD,(DINT(I),I=1,NUM)
1165	FORMAT(' There are',I4,' lines for standard ',A40,
     &	   (40(/X,5F10.4)))
C
	WRITE(IO6,1170)
1170	FORMAT(/' Results of peak matching:'/
    	1/5X'Peak No.',4X'Observed 2-T',4X'Calculated 2-T',8X'Delta',
    	1/70('-'))
C
	DO 200 I=1,NPEAKS
	THETA = THETA1(I)
130	THETH = THETA
	IF(I.EQ.1)THETL = THETA
C	IF(NAST.EQ.'*') GO TO 200
	DO 190 II=1,NUM
	DIFF=DINT(II)-THETA
	IF(ABS(DIFF).GT.EWIND) GO TO 190
	INUM=INUM+1
140	XTHE(INUM)=THETA
	DTHE(INUM)=DIFF
	LOC(INUM) = II
        IFIND(INUM) = I
C
C  CHECK IF THE PREVIOUS VALUE WAS ALSO WITHIN THE WINDOW
C
	IF(INUM.EQ.1)GO TO 170
	DIFF1 = XTHE(INUM-1) - THETA
	IF(ABS(DIFF1).GT.EWIND)GO TO 170
C
C  FIND THE SMALLEST
C
C  IF THE DIFFERENCE IS SMALL KEEP IT (FOR BI12GEO20 69 DEG LINE)
	IF(ABS(DTHE(INUM-1))-ABS(DIFF).LT..01)GO TO 150
	IF(ABS(DTHE(INUM-1)).GT.ABS(DIFF))GO TO 160
150	WRITE(IO6,1200)THETA
1200	FORMAT(' Warning second line at ',F9.4,' ignored')
	INUM = INUM - 1
	GO TO 190
C
C  ACCEPT THE NEW LINE
C
160	WRITE(IO6,1210)
1210	FORMAT(' ***** Previous match ignored due to better value***')
	INUM = INUM - 1
	GO TO 140
170	WRITE(IO6,1220) I,THETA,DINT(II),DTHE(INUM)
1220	FORMAT(6X,'#',I3,8X,F9.4,8X,F9.4,8X,F9.4)
190	CONTINUE
200	CONTINUE
C
C  WE NOW GIVE A CHANCE TO REJECT BAD MATCHES
C
201 	WRITE(IO6,1221)
1221	FORMAT(/' If any of the above matches are obviously bad',
     1' you may reject them now,'/' before we try to fit a curve',
     2' through the data.  Enter the numbers next to'/' the #',
     3' signs, of points to be rejected.  Separate numbers by commas.'/
     4' (Hit return to accept all matches): '$)
     	READ(IO5,1320,END=1) (IREJ(I),I=1,5)
     	IF(IREJ(1) .EQ. 0) GO TO 210
C
C  CUT THE BAD PEAKS OUT OF THE LIST
C
     	DO 208 I=5,1,-1
     	IF(IREJ(I) .EQ. 0) GO TO 208
C  FIND THE SERIAL INDEX FROM THE MASTER LIST
       DO 202 III=1,INUM
       IF(IFIND(III) .EQ. IREJ(I)) GO TO 204
202	CONTINUE
       WRITE(IO6,1222) IREJ(I)
1222	FORMAT(/' *** Match ',I3,' is not in the list! Try again.')
       GO TO 201
204	DO 206 II=III,INUM-1
     	XTHE(II) = XTHE(II+1)
     	DTHE(II) = DTHE(II+1)
     	LOC(II) = LOC(II+1)
206	CONTINUE
     	INUM = INUM - 1
208	CONTINUE
C
C 	FIT A POLYNOMIAL BETWEEN DTHE AND XTHE
C
210	IF(INUM.GT.1) GO TO 220
C
C 	PROBLEMS: CAN NOT FIND A SINGLE INT.STD. LINE
C
	IF(INUM.EQ.0) WRITE(IO7,1230)
	IF(INUM.EQ.1) WRITE(IO7,1240)
1230	FORMAT(///,'???? Cannot find any int. std. lines :')
1240	FORMAT(///,'???? Found only one int. std. line :')
	GO TO 999
220	WRITE(IO6,1250)INUM
1250	FORMAT(/'$Input the degree of the polynomial to fit to the',I3,
     1 ' points <1>: ')
	READ(IO5,1255,END=1) MCO
1255	FORMAT(I)
	IF(MCO.LE.0) MCO = 1	!JDW
	IF(MCO.LE.INUM-1.AND.MCO.LE.5)GO TO 230
	WRITE(IO6,1260)MCO
1260	FORMAT(/1XI3,' Is too high a degree polynomial')
	GO TO 220
230	SD=REGR(INUM,XTHE,DTHE,MCO,YINT,COEF)
	WRITE(IO6,1270)MCO,YINT,(COEF(I),I=1,MCO)
1270	FORMAT(/' Results of polynomial regression of degree ',I2/
     1 1X,6G16.9)
	WRITE(IO6,1280)
1280	FORMAT(/2X'NUM',2X,'Observed 2T',2X'True 2-T',2X'Delta 2-T',
     1 2X'Calc. 2-T',3X'True-Calc'/70('-'))
	SUM = 0.0
	DO 240 I=1,INUM
	L = LOC(I)
	CALL POLY(TRUE,XTHE(I),MCO,YINT,COEF)
	TRUE = TRUE + XTHE(I)
	DT = DINT(L) - TRUE
	SUM = SUM + ABS(DT)
	WRITE(IO6,1290)I,XTHE(I),DINT(L),DTHE(I),TRUE,DT
1290	FORMAT(2X,I3,2X,5(2X,F9.5))
240	CONTINUE
	SUM = SUM /INUM
	WRITE(IO6,1300)SUM
1300	FORMAT(/' The average delta 2 theta = ',F9.5)
C
C  ASK IF ANY OF THE MATCHES SHOULD BE REJECTED.
C
     	WRITE(IO6,1310)
1310	FORMAT(/' You now have the option to reprocess the data',
     1' rejecting up to 5 of the above'/' data points.  Enter',
     2' the numbers (from the NUM column of the above table) of'/
     3' any peaks you wish to reject separated by commas'/
     4' (Hit return to accept all matches): '$)
     	READ(IO5,1320,END=1) (IREJ(I),I=1,5)
1320	FORMAT(5I)
     	IF(IREJ(1) .EQ. 0) GO TO 250
C
C  CUT THE BAD PEAKS OUT OF THE LIST
C
     	DO 246 I=5,1,-1
     	IF(IREJ(I) .EQ. 0) GO TO 246
     	DO 245 II=IREJ(I),INUM-1
     	XTHE(II) = XTHE(II+1)
     	DTHE(II) = DTHE(II+1)
     	LOC(II) = LOC(II+1)
245	CONTINUE
     	INUM = INUM - 1
246	CONTINUE
     	GO TO 230
C
C **	NOW CORRECT d'S FOR THE INTERNAL STANDARD 
C **    AND WRITE IT ON LFILE
C
C  GET DATE AND TIME
C
250      CALL DATE( IDAT )
C
C  WRITE TITLE AND DATE TO FILE ON LUN 20
C
	IF(IPRIR.EQ.'T')GO TO 255
C
C  SET UP THE ADR FILE HEADER INFO
C
       IF (IOUT .EQ. 'A' .OR. IOUT .EQ. 'a') THEN
		WRITE(IO20,1301) IDAT, TITLE
1301		FORMAT(//'XRAY POWDER DATA ANALYSIS',/3X,
     1		'PATTERN CORRECTED FOR INTERNAL STANDARD BY INTCAL'/
     1  	'DATE: ',9A1,///'SAMPLE ID: ',
     2		80A1,64X,/1X/1X)
C
C  WRITE OUT THE SCAN LIMITS
C
		WRITE(IO20,2001)ALO,TIME,AHI,STEP
2001		FORMAT('STARTING ANGLE : ',F8.4,5X,'COUNT TIME  : ',
     1		F6.1/'ENDING ANGLE  : ',F9.4,/'ANGLE INCREMENT: ',F8.4)
C
C  WRITE OUT THE ADR HEADER RECORD
C
		WRITE(IO20,2002)
2002		FORMAT(//4X'PEAK',3X'2-THETA',6X'd',6X'RELINT',
     1		3X'h  k  l',4X'ABSINT',4X'FWHM'/70('-'))    
       ENDIF
C
C  CALCULATE THE hkl VALUES IF REQUESTED
C
255	NHKL = 0
	IF(INDEX.EQ.1)CALL HKL(I1,I2,I3,T1,D1,NHKL)
	IF(MCO.EQ.1) COEF(2)=0.
	WRITE(IO6,2003)
2003	FORMAT(//4X'PEAK',3X'2-THETA',6X'd',6X'RELINT',3X'h  k  l',
    	14X'ABSINT',4X'FWHM'/70('-'))    
C
C  ADJUST THE DIFFRAC 11 OUTPUT FILE FOR UPDATED INFORMATION
C
       IF(IOUT .EQ. 'D') THEN
C  FIRST SET PROGRAM CODE (WORD 20) TO "CALI"
       		WRITE(IO20'20) CODE
C  		NEXT SUPPLY CODE THAT HKL'S ARE PRESENT
       		IF(INDEX .EQ. 1) THEN
       			READ(IO10'25) F4
       			INT2(2) = 1
       			WRITE(IO20'25) F4
       		ENDIF
       	ENDIF
C
C  NOW BEGIN CORRECTING AND OUTPUTING EACH 2-THETA VALUE
C
       MEXT = 39
	N = 0
       NN = 0
	I1(NHKL+1) = 0
	I2(NHKL+1) = 0
	I3(NHKL+1) = 0
C
	DO 500 I=1,NPEAKS
C	read(io10,1180,end=510) theta,nast,d,ai,bb,fwhm
	THETA = THETA1(I)
	AI = RELINT(I)
	BB = INTCPS(I)
	FWHM = FWHM1(I)
C
C  SKIP ANY ALPHA 2 LINES
C
C430	IF(NAST .EQ. '*') GO TO 500
C
C  SKIP THE INTERNAL STANDARD LINES
C
	DO 440 JJ=1,INUM
	IF(THETA-XTHE(JJ).EQ.0.0)GO TO 500
440	CONTINUE
C
C  CORRECT LINE FOR INTERNAL STANDARD CURVE
C
	CALL POLY(TRUE,THETA,MCO,YINT,COEF)
	THETA = THETA + TRUE
C
C  COMPUTE THE CORRECRED d VALUE
C
	D = 1.540598 / (2.0 * SIN((THETA *DTR) / 2.0))
	JJ = NHKL + 1
	IF(INDEX.EQ.0)GO TO 455
C
C  FIND THE INDICES
C
	DELMIN = 5000.
	DO 450 JJ=1,NHKL
	DELTA = ABS(THETA - T1(JJ))
	IF(DELTA .GT. EWIND) GO TO 450
	IF(DELTA .GE. DELMIN) GO TO 450
	DELMIN = DELTA
	JK = JJ
450	CONTINUE
	IF(DELMIN .GT. EWIND)GO TO 452
	GO TO 460
C
C  NO hkl FOUND PRINT WARNING AND OUTPUT ZERO'S
C
452	WRITE(IO6,1360) THETA
1360	FORMAT(/' **** No indices found for ',F10.5)
455	JK = NHKL + 1
460	N = N + 1
	IF(IPRIR.NE.'T') THEN
C  WRITE THE ADR FILE
       	   IF(IOUT .NE. 'D') THEN
       		WRITE(IO20,1370) N,THETA,NAST,D,AI,I1(JK),
     1  	I2(JK),I3(JK),BB,FWHM
C
C  		WRITE DIFFRAC 11 FILE
C
C		39-40	h,k,l,0 (4I*2)
C		41	D-value
C		42	Intensity
C		43	FWHM
C		44	Residual intensity
C		45-51	Next peak
C
       	   ELSE
C  SKIP TEST AND DO IT       	   	IF(INDEX .EQ.1) THEN
       			INT2(1) = I1(JK)
       			INT2(2) = I2(JK)
       			WRITE(IO20'MEXT) F4
       			INT2(1) = I3(JK)
       			INT2(2) = 0
       			WRITE(IO20'MEXT) F4
C       				ELSE
C		       			NEXT = NEXT + 2
C			       		ENDIF
       	   WRITE(IO20'MEXT) D
       	   WRITE(IO20'MEXT) AI
       	   WRITE(IO20'MEXT) FWHM
       	   F4 = 0.0
       	   WRITE(IO20'MEXT) F4
           NN = NN + 1
       	   ENDIF
       ENDIF
	IF(IPRIR.NE.'F')WRITE(IO6,1370) N,THETA,NAST,D,AI,I1(JK),
     1  I2(JK),I3(JK),BB,FWHM
1370	FORMAT(4X,'#',I3,2X,F9.5,A1,1X,F7.4,2X,F6.2,2X,3I3,1X,F9.1,F9.4)
500	CONTINUE
C
C  PUT OUT THE NUMBER OF PEAKS ONTO THE DIFFRAC 11 FILE
C  WORD 21 = npeak, nresid.peaks	(2I*2)
       INT2(1) = NN
       INT2(2) = 0
       IF(IOUT .EQ. 'D') THEN
       		WRITE(IO20'21) F4
       ENDIF
510    	WRITE(IO6,520)
520	FORMAT(70('-'))
	IF(IPRIR.EQ.'T')GO TO 350
       CLOSE(UNIT = IO10)
	CLOSE(UNIT=IO20)
C
C  LASTLY OUTPUT EACH DATA SET 
C
C  WRITE THE UNCORRECTED DELTA TWO THETAS OUT FOR PLOTTING
C  IN SINGLE PRECISSION
     	PLTFIL = NFILE
	PLTFIL(ITEMP+1:ITEMP+3)='PLT'
	OPEN(UNIT=IO30, NAME=PLTFIL, STATUS='UNKNOWN', 
     1 FORM='FORMATTED', CARRIAGECONTROL='LIST')
    	WRITE(IO30,2400) INUM
	DO 2700 I=1,INUM
	PLX(I) = XTHE(I)
	PLY(I) = DTHE(I)
    	WRITE(IO30,2600) PLX(I), PLY(I)
2700	CONTINUE
	INTVAL = 100
	DX = (THETH - THETL) / INTVAL
	X = THETL - DX
	DO 360 I=1,INTVAL
	X = X + DX
	PLX(I) = X
	CALL POLY(PLY(I),X,MCO,YINT,COEF)
360	CONTINUE
C
C  WRITE THE DATA OUT FOR PLOTTING
	WRITE(IO30,2400) INTVAL
2400	FORMAT(I3)
    	DO 2500 I=1,INTVAL
2500	WRITE(IO30,2600) PLX(I), PLY(I)
2600	FORMAT(2G)
C===============================================================
C  SET UP THE PLOT PARAMETER FILE FOR TKPLOT, GGPLOT ROUTINES
C===============================================================
350     IJUNK = NFILE
       	IJUNK(ITEMP+1:ITEMP+3)='PAR'
C
       	OPEN(UNIT=1, NAME=IJUNK, STATUS='UNKNOWN', 
     1  FORM='FORMATTED', CARRIAGECONTROL='LIST')
C
       	WRITE(1,111)
111	FORMAT('GN3I=  0, 70, 0, 1, 0, 6'/
     1  'XLB1=TWO-THETA (DEGREES)'/'YLB1=DELTA TWO-THETA'/
     1  'TITL=INTERNAL STANDARD CALIBRATION CURVE'/
     1  'SSPI= 1, 2, -3'/'SSPR= 1, 0.1'/'SSPI= 2, 3'/
     1  'LABL= 1,OBSERVED POINTS'/'LABL= 2,CALIBRATED CURVE')
    	CLOSE(UNIT=1)
C
C  PLOT PARAMETER FILE IS COMPLETE
C
	WRITE(IO6,2800) PLTFIL
2800	FORMAT(/10X'*** Plot file written to ',A40,
     1/' To view the plot when finished, run TKPLOT or GGPLOT'//)
    	CLOSE(UNIT=IO30)
C  Go ask for another file.
	GOTO 1
C
C  Explain about the control key.
998	WRITE(IO7,9981)
9981	FORMAT(' That''s not what I meant but I know what you mean.'/
     &  ' In computer jargon ^Z means hold down the control key and',
     &  ' Hit Z')
C
999	WRITE(IO6,9999)
9999	FORMAT(/10X'*** Normal termination of program ***')
    	CALL EXIT
	END

	SUBROUTINE DISP(DT,XTHE,AVDISP,IDN,DTR)
	IMPLICIT REAL*8 (A-H,O-Z)
	CHARACTER*1 IDN
	IF(IDN.EQ.'1')GO TO 222
C
C  HANDLE A REFLECTION SAMPLE DISPLACEMENT FOR DIFFRACTOMETER RADIUS OF R mm
C
	R = 175.0
	IF(IDN.NE.'2')R = 200.0
	DT = (2.0 * AVDISP * COS(XTHE/2.*DTR)) / (R * DTR)
	GO TO 223
C
C  HANDLE A DISPLACEMENT FOR THE MOVING SAMPLE TRANSMISSION GUINIER
C***********NOTE: A SAMPLE ANGLE OF 45 DEG TO THE BEAM IS ASSUMED (IE ALPHA)
C
222	ALPHA = 45.0
	R = 280. * COS( ALPHA * DTR)
	DT = 2.0 * AVDISP / (DTR * R *
     1  (SIN((XTHE+ALPHA)*DTR) - ((COS((XTHE+ALPHA)*DTR) + 1.0)
     1  / TAN(XTHE*DTR))))
223	RETURN
	END

	SUBROUTINE POLY(TRUE,THET,MCO,YINT,COEF)
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION COEF(5)
	TRUE = YINT
	DO 103 J=1,MCO
103	TRUE = TRUE + COEF(J) * THET**J
	RETURN
	END

	SUBROUTINE SI(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR SILICON PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(12)
C
	DATA DSI/ 28.44271,  47.30334,  56.12309,  69.13093,  
     1 76.37729, 88.03148,  94.95399, 106.71030, 114.09426, 
     1 127.54728, 136.89663, 158.63824/
C
	NUM=12
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE SB(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR SILVER BEHENATE PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSB(18)
C
	DATA DSB/1.51203, 3.02432, 4.53715, 6.05076, 7.56543,
     1  9.08142, 10.59901, 12.11847, 13.64007, 15.16409, 16.69081,
     1 18.22052, 19.75351, 21.29008, 22.83054, 24.37519, 25.92435,
     1 27.47834/
C
	NUM=18
C
	DO 10 I=1,NUM
	DINT(I)=DSB(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE SS(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR SILVER BEHENATE + SI PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSS(30)
C
	DATA DSS/1.51203, 3.02432, 4.53715, 6.05076, 7.56543,
     1 9.08142,  10.59901, 12.11847, 13.64007, 15.16409, 16.69081,
     1 18.22052, 19.75351, 21.29008, 22.83054, 24.37519, 25.92435,
     1 27.47834, 28.44271, 47.30334, 56.12309, 69.13093, 76.37729, 
     1 88.03148,  94.95399, 106.71030, 114.09426, 127.54728, 
     1 136.89663, 158.63824/
C/
C
	NUM=30
C
	DO 10 I=1,NUM
	DINT(I)=DSS(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE AG(DINT,NUM)
C
C ** THIS ROUTINE RETURNS SILVER TWO-THETA VALUES IN DINT
C **	NUM IS THE NUMBER OF PEAKS
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DAG(10)
	DATA DAG/38.112,44.295,64.437,77.390,81.533,97.875,
     +	         110.499,114.914,134.871,156.737/
	NUM=10
	DO 10 I=1,NUM
10	DINT(I)=DAG(I)
	RETURN
	END

	SUBROUTINE BI(DINT,NUM)
C
C ** THIS ROUTINE RETURNS BI12GEO20 TWO-THETA VALUES IN DINT
C **	NUM IS THE NUMBER OF PEAKS
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DAG(80)
	DATA DAG/
     1 12.3281, 17.4685, 21.4365, 24.8018, 27.7848, 30.4981, 33.0088,
     1 35.3604, 37.5831, 39.6991, 41.7247, 43.6730, 45.5543, 49.1474,
     1 50.8719, 52.5552, 54.2014, 55.8142, 57.3966, 58.9517, 60.4817,
     1 61.9889, 63.4753, 64.9427, 66.3927, 67.8268, 69.2465, 70.6530,
     1 73.4309, 74.8046, 76.1693, 77.5261, 78.8759, 80.2194, 81.5575,
     1 82.8910, 84.2206, 85.5471, 86.8712, 88.1936, 89.5151, 90.8363,
     1 92.1580, 94.8055, 96.1327, 97.4633, 98.7979,100.1374,101.4824,
     1102.8340,104.1928,105.5598,108.3223,109.7199,111.1297,112.5531,
     1113.9913,116.9181,118.4099,119.9229,121.4594,123.0214,124.6117,
     1126.2330,127.8887,129.5825,131.3188,133.1026,134.9400,136.8382,
     1138.8060,140.8543,145.2514,147.6419,150.2015,152.9787,156.0495,
     1159.5464,163.7428,169.4431/
	NUM=80
	DO 10 I=1,NUM
10	DINT(I)=DAG(I)
	RETURN
	END


	SUBROUTINE BS(DINT,NUM)
C
C ** THIS ROUTINE RETURNS BI12GEO20 + SI TWO-THETA VALUES IN DINT
C **	NUM IS THE NUMBER OF PEAKS
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DAG(92)
	DATA DAG/
     1 12.3281, 17.4685, 21.4365, 24.8018, 27.7848, 30.4981, 33.0088,
     1 35.3604, 37.5831, 39.6991, 41.7247, 43.6730, 45.5543, 49.1474,
     1 50.8719, 52.5552, 54.2014, 55.8142, 57.3966, 58.9517, 60.4817,
     1 61.9889, 63.4753, 64.9427, 66.3927, 67.8268, 69.2465, 70.6530,
     1 73.4309, 74.8046, 76.1693, 77.5261, 78.8759, 80.2194, 81.5575,
     1 82.8910, 84.2206, 85.5471, 86.8712, 88.1936, 89.5151, 90.8363,
     1 92.1580, 94.8055, 96.1327, 97.4633, 98.7979,100.1374,101.4824,
     1102.8340,104.1928,105.5598,108.3223,109.7199,111.1297,112.5531,
     1113.9913,116.9181,118.4099,119.9229,121.4594,123.0214,124.6117,
     1126.2330,127.8887,129.5825,131.3188,133.1026,134.9400,136.8382,
     1138.8060,140.8543,145.2514,147.6419,150.2015,152.9787,156.0495,
     1159.5464,163.7428,169.4431,
     D 28.44271,  47.30334,  56.12309,  69.13093,  
     E 76.37729, 88.03148,  94.95399, 106.71030, 114.09426, 
     F 127.54728, 136.89663, 158.63824/
	NUM=92
C
C  SORT THEM TOGETHER
C
	DO 5 I=1,NUM-1
	DO 5 J=I,NUM
	IF(DAG(I).LE.DAG(J))GO TO 5
	T = DAG(I)
	DAG(I) = DAG(J)
	DAG(J) = T
5	CONTINUE
	DO 10 I=1,NUM
10	DINT(I)=DAG(I)
	RETURN
	END

	SUBROUTINE NI(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR NICKEL PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(8)
	DATA DSI/
     1 44.4970, 51.8505, 76.3829, 92.9395, 98.4446,121.9471,144.6714,
     1155.7024/
C
	NUM=8
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE NO(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR NICKEl OXIDE PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(10)
	DATA DSI/
     1 37.2561, 43.2880, 62.8815, 75.4177, 79.4118, 95.0678, 
     2 107.0016, 111.1259, 129.2350, 146.7792/
C
	NUM=10
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE CU(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR COPPER PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(8)
C
	DATA DSI/ 43.31686, 50.44927, 74.12602, 89.93736, 95.14729,
     2 116.93344, 136.50079, 144.70555/
C
	NUM=8
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE FL(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR FLUOROPHLOGOPITE PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(12)
	DATA DSI/
     1  8.8525, 17.7585, 26.7742, 35.9624, 45.3969, 55.1691, 65.3991,
     1 76.2545, 87.9885,101.0248,116.1928,135.6741/
	NUM = 12
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE FC(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR FLUOROPHLOGOPITE AND CU PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(20)
	DATA DSI/
     1  8.8525, 17.7585, 26.7742, 35.9624, 43.3169, 45.3969, 50.44927, 
     2 55.1691, 65.3991, 74.1260, 76.2545, 87.9885, 89.9374, 95.14729,
     3101.0248,116.1928,116.9334, 135.674,136.5008, 144.70555/
	NUM = 20
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE FS(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR SILICON AND
C **	FLUOROPHLOGOPITE (NBS SRM) PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(24)
C
	DATA DSI/ 28.44271,  47.30334,  56.12309,  69.13093,  
     1 76.37729, 88.03148,  94.95399, 106.71030, 114.09426, 
     1 127.54728, 136.89663, 158.63824,
     1  8.8525, 17.7585, 26.7742, 35.9624, 45.3969, 55.1691, 65.3991,
     1 76.2545, 87.9885,101.0248,116.1928,135.6741/
C
	NUM = 24
C
C  SORT THEM TOGETHER
C
	DO 5 I=1,NUM-1
	DO 5 J=I,NUM
	IF(DSI(I).LE.DSI(J))GO TO 5
	T = DSI(I)
	DSI(I) = DSI(J)
	DSI(J) = T
5	CONTINUE
C
	DO 10 I=1,12
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

      FUNCTION REGR (M, X, Y, N, A, B)
C
C     THIS PROGRAM ORIGINALLY CAME FROM THE BOOK:
C
C     B. CARNAHAN, H. LUTHER, AND J. WILKES, APPLIED NUMERICAL
C     METHODS,  JOHN WILEY AND SONS, NEW YORK, 1969, P580.
C
C
C            FUNCTION REGR, WHICH CARRIES OUT AN NTH-ORDER
C            POLYNOMIAL REGRESSION ON M DATA POINTS CONTAINED
C            IN THE X AND Y ARRAYS.  THE FUNCTION NORMALLY
C            RETURNS THE STANDARD DEVIATION S OF THE POINTS
C            ABOUT THE REGRESSION CURVE.  THE REGRESSION
C            COEFFICIENTS ARE PLACED IN A AND B(1)...B(N).
C            HOWEVER, IF THE SIMULTANEOUS EQUATION SOLVING
C            ROUTINE SIMUL ENCOUNTERS A NEAR-SINGULAR MATRIX,
C            THE FUNCTION RETURNS THE VALUE 0.0 AS REGR.
C
C
C     PARAMETER DEFINITIONS:
C
C     M - THE NUMBER OF DATA POINTS IN X AND Y.
C
C     X,Y - THE ARRAYS CONTAINING THE ABSCISSA AND ORDINATE
C           DATA VALUES RESPECTIVELY
C
C     N - THE ORDER OF THE POLYNOMIAL DESIRED FOR THE FIT.
C
C     A - Y INTERCEPT OF THE POLYNOMIAL.
C
C     B - ARRAY CONTAINING THE COEFFICIENTS OF THE
C         POLYNOMIAL.
C
C     REGR - RETURNS AS THE VALUE OF THE STANDARD
C            DEVIATION OF THE FIT TO THE DATA. IF
C            A NEAR SINGULAR OR IF A SINGULAR
C            MATRIX IS ENCOUNTERED THE VALUE OF
C            0 (ZERO) IS RETURNED AS REGR.
C
C
C     *****  NOTE  *****
C
C     THE POLYNOMIAL HAS THE FORM OF:
C
C     ORDINATE = A + B(I) * (ABSCISSA ** I)
C        FOR I=1 TO N
C
C
      IMPLICIT REAL*8 (A - H, O - Z)
      REAL*8 A, B, REGR, X, Y
      DIMENSION C(11,11), SX(20), SYX(10), CYX(10), X(100),Y(100),
     1B(10)
      DATA EPS/ 1.0D-20/
C
C     .....  COMPUTE SUMS OF POWERS AND PRODUCTS  .....
      IO6 = 6
      NTWO = 2*N
      NP1 = N + 1
      SY = 0.0
      SYY = 0.0
      DO 1 I=1,N
      NPI = N + I
      SX(I) = 0.0
      SX(NPI) = 0.0
      SYX(I) = 0.0
1     CONTINUE
      DO 3 I=1,M
      SY = SY + Y(I)
      SYY = SYY + Y(I)**2
      DUM=1.0
      DO 2 J=1,N
      DUM = DUM*X(I)
      SX(J) = SX(J) + DUM
      SYX(J) = SYX(J) + Y(I)*DUM
2     CONTINUE
      DO 3 J=NP1,NTWO
      DUM = DUM*X(I)
      SX(J) = SX(J) + DUM
3     CONTINUE
C
C     .....  COMPUTE COEFFICIENTS C(I,J)  .....
C
      FM = M
      CYY = SYY - SY*SY/FM
      DO 4 I=1,N
      CYX(I) = SYX(I) - SY*SX(I)/FM
      C(I,NP1) = CYX(I)
      DO 4 J=1,N
      IPJ = I + J
      C(I,J) = SX(IPJ) - SX(I)*SX(J)/FM
4     CONTINUE
C
C     .....  CALL ON SIMUL TO SOLVE SIMULTANEOUS EQUATIONS
C
      DET = SIMUL (N, C, B, EPS, 1, 11)
C
C     I INSERTED THE COMMENTS INFRONT OF THE WRITE
C     STATEMENT BECAUSE I DIDN'T WANT TO PRINT IT OUT.
C
C     WRITE(IO6,200) DET
C200  FORMAT(13H0       DET = ,E14.8)
      IF(DET.NE.0.0) GO TO 6
      REGR = 0.0
      RETURN
C
C     .....  COMPUTE INTERCEPT A AND STANDARD DEVIATION S  .....
C
6     DUM = SY
      TEMP = CYY
      DO 7 I=1,N
      DUM = DUM - B(I)*SX(I)
      TEMP = TEMP - B(I)*CYX(I)
7     CONTINUE
      A = DUM/FM
      DENOM = M - N - 1
C
C     *****  NOTE  *****
C
C     I INSERTED THE FOLLOWING IF STATEMENT TO ELIMINATE
C     TEMP BEING NEAR ZERO BUT SLIGHTLY NEGATIVE.
C     AN EXAMPLE IS IF TEMP EQUALS -.12345E-15.
C
C     THIS STATEMENT IS NOT SUPPLIED IN
C     THE BOOK VERSION, BUT RATHER IS A HOMEBREW PATCH.
C
C     ******************
      IF(TEMP.LE.EPS) TEMP=EPS
      S = DSQRT(TEMP/DENOM)
      REGR = S
      RETURN
      END

      FUNCTION SIMUL( N, A, X, EPS, INDIC, NRC)
C
C     THIS PROGRAM ORIGINALLY CAME FROM THE BOOK:
C
C     B. CARNAHAN, H. LUTHER, AND J. WILKES, APPLIED NUMERICAL
C     METHODS,  JOHN WILEY AND SONS, NEW YORK, 1969, P290
C
C
C        WHEN INDIC IS NEGATIVE, SIMUL COMPUTES THE INVERSE
C        OF THE N BY N MATRIX A IN PLACE.  WHEN INDIC IS ZERO,
C        SIMUL COMPUTES THE N SOLUTIONS X(1)...X(N) CORRESPONDING
C        TO THE SET OF LINEAR EQUATIONS WITH AUGMENTED MATRIX
C        OF COEFFICIENTS IN THE N BY N+1 ARRAY A AND IN
C        ADDITION COMPUTES TH INVERSE OF THE COEFFICIENT
C        MATRIX IN PLACE AS ABOVE.  IF INDIC IS POSITIVE,
C        THE SET OF LINEAR EQUATIONS IS SOLVED BUT THE
C        INVERSE IS NOT COMPUTED IN PLACE.  THE GAUSS-
C        JORDAN COMPLETE ELIMINATION METHOD IS EMPLOYED
C        WITH THE MAXIMUM PIVOT STRATEGY.  ROW AND
C        COLUMN SUBSCRIPTS OF SUCCESIVE PIVOT ELEMENTS ARE
C        SAVED IN ORDER IN THE IROW AND JCOL ARRAYS
C        RESPECTIVELY.  K IS THE PIVOT COUNTER, PIVOT IS
C        THE ALGEBRAIC VALUE OF THE PIVOT ELEMENT. MAX
C        THE NUMBER OF COLUMNS IN A AND DETER THE DETERMINANT
C        OF THE COEFFICIENT MATRIX.  THE SOLUTIONS ARE
C        DETERMINED IN THE N+1TH COLUMN OF A AND THEN
C        UNSCRAMBLED AND PUT IN ORDER IN X(1)...X(N)
C        USING THE PIVOT SUBSCRIPT INFORMATION AVAILABLE
C        IN THE IROW AND JCOL ARRAYS.  THE SIGN OF THE
C        DETERMINANT IS ADJUSTED, IF NECESSARY, BY DETERMINING
C        IF AN EVEN OR ODD NUMBER OF PAIRWISE INTERCHANGES IS
C        REQUIRED TO PUT THE ELEMENTS OF THE JORD ARRAY IN
C        ASCENDING SEQUENCE WHERE JORD(IROW(I))  =  JCOL(I).
C        IF THE INVERSE IS REQUIRED, IT IS UNSCRAMBLED IN PLACE
C        USING Y(1)...Y(N) AS TEMPORARY STORAGE.  THE VALUE OF
C        THE DETERMINANT IS GIVEN AS THE VALUE OF THE FUNCTION.
C        SHOULD THE POTENTIAL PIVOT OF LARGEST MAGNITUDE BE
C        SMALLER IN MAGNITUDE THAN EPS, THE MATRIX IS
C        CONSIDERED TO BE SINGULAR AND A TRUE ZERO IS RETURNED
C        AS THE VALUE OF THE FUNCTION.
C
      IMPLICIT REAL*8(A-H, O-Z)
      REAL*8 A, X, EPS, SIMUL
      DIMENSION IROW(50),JCOL(50),JORD(50),Y(50),A(NRC,NRC),X(N)
C
      MAX = N
      IO6 = 6
      IF(INDIC.GE.0) MAX = N+1
C
C     .....  IS N LARGER THAN 50  .....
      IF(N.LE.50) GO TO 5
      WRITE(IO6,200)
      SIMUL = 0.0
      RETURN
C
C     .....  BEGIN ELIMINATION PROCEDURE  .....
5     DETER = 1.
      DO 18 K = 1,N
      KM1 = K - 1
C
C     .....  SEARCH FOR THE PIVOT ELEMENT  .....
      PIVOT = 0.0
      DO 11 I = 1,N
      DO 11 J = 1,N
C     .....  SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS  .....
      IF(K.EQ.1) GO TO 9
      DO 8 ISCAN = 1,KM1
      DO 8 JSCAN = 1,KM1
      IF(I.EQ.IROW(ISCAN)) GO TO 11
      IF(J.EQ.JCOL(JSCAN)) GO TO 11
8     CONTINUE
9     IF(DABS(A(I,J)).LE.DABS(PIVOT)) GO TO 11
      PIVOT = A(I,J)
      IROW(K) = I
      JCOL(K) = J
11    CONTINUE
C
C     .....  INSURE THAT SELECTED PIVOT IS LARGER THAN EPS  .....
      IF(DABS(PIVOT).GT.EPS) GO TO 13
      SIMUL = 0.0
      RETURN
C
C     .....  UPDATE THE DETERMINANT VALUE  .....
13    IROWK = IROW(K)
      JCOLK = JCOL(K)
      DETER = DETER*PIVOT
C
C     .....  NORMALIZE PIVOT ROW ELEMENTS  .....
      DO 14 J = 1,MAX
      A(IROWK,J) = A(IROWK,J)/PIVOT
14    CONTINUE
C
C     .....  CARRY OUT ELIMINATION AND DEVELOP INVERSE .....
      A(IROWK,JCOLK) = 1./PIVOT
      DO 18 I = 1,N
      AIJCK = A(I,JCOLK)
      IF(I.EQ.IROWK) GO TO 18
      A(I,JCOLK) = - AIJCK/PIVOT
      DO 17 J = 1,MAX
      IF(J.NE.JCOLK) A(I,J) = A(I,J) - AIJCK*A(IROWK,J)
17    CONTINUE
18    CONTINUE
C
C     .....  ORDER SOLUTION VALUES (IF ANY) AND CREATE JORD ARRAY  .....
      DO 20 I = 1,N
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      JORD(IROWI) = JCOLI
      IF(INDIC.GE.0) X(JCOLI) = A(IROWI,MAX)
20    CONTINUE
C
C     .....  ADJUST SIGN OF DETERMINANT  .....
      INTCH = 0
      NM1 = N - 1
      DO 22 I = 1,NM1
      IP1 = I + 1
      DO 22 J = IP1,N
      IF(JORD(J).GE.JORD(I)) GO TO 22
      JTEMP = JORD(J)
      JORD(J) = JORD(I)
      JORD(I) = JTEMP
      INTCH = INTCH + 1
22    CONTINUE
      IF(INTCH.NE.INTCH/2*2) DETER = - DETER
C
C     .....  IF INDIC IS POSITIVE RETURN WITH RESULTS  .....
24    IF(INDIC.LE.0) GO TO 26
      SIMUL = DETER
      RETURN
C
C     .....  IF INDIC IS NEGATIVE OR ZERO, UNSCRAMBLE THE INVERSE
C            FIRST BY ROWS  .....
26    DO 28 J = 1,N
      DO 27 I = 1,N
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      Y(JCOLI) = A(IROWI,J)
27    CONTINUE
      DO 28 I = 1,N
      A(I,J) = Y(I)
28    CONTINUE
C     .....  THEN BY COLUMNS  .....
      DO 30 I = 1,N
      DO 29 J = 1,N
      IROWJ = IROW(J)
      JCOLJ = JCOL(J)
      Y(IROWJ) = A(I,JCOLJ)
29    CONTINUE
      DO 30 J = 1,N
      A(I,J) = Y(J)
30    CONTINUE
C
C     .....  RETURN FOR INDIC NEGATIVE OR ZERO  .....
      SIMUL = DETER
      RETURN
C
C     .....  FORMAT FOR OUTPUT STATEMENT  .....
200   FORMAT(10H0N TOO BIG )
      END

	SUBROUTINE HKL(I1,I2,I3,T1,D1,N)
	IMPLICIT REAL*8(A-H,O-Z)
C
C  2-THETA AND HKL FINDER
C
	INTEGER H
	DIMENSION T1(200),D1(200),I1(200),I2(200),I3(200)
	DATA WAVE / 1.5405981/, IREAD/5/, IWRIT/6/
	DTR = 3.141592654 / 180.
110	WRITE(IWRIT,1000)
1000	FORMAT(/' HKL AND 2 THETA COMPUTATION'//
     1' INPUT ONLY THE NEEDED LATTICE PARAMETERS FOR YOUR SYMMETRY'/
     1'$a, b, c, alpha, beta, gamma: ')
	READ(IREAD,1010)A,B,C,AL,BE,GA
1010	FORMAT(6G10.5)
	WRITE(IWRIT,1005)
1005	FORMAT(/' INPUT THE WAVELENGTH (if other than Cu Ka1)')
	READ(IREAD,1010)W
	IF(W.NE.0.0)WAVE = W
	WRITE(IWRIT,1011)
1011	FORMAT(/' INPUT LATTICE TYPE (P,I OR F)')
	IF(B.EQ.0.0)B = A
	IF(C.EQ.0.0)C = A
	IF(AL.EQ.0.0)AL = 90.
	IF(BE.EQ.0.0)BE = 90.
	IF(GA.EQ.0.0)GA = 90.
	READ(IREAD,1012)LAT
1012	FORMAT(A1)
	ILAT = 1
	IF(LAT.EQ.1HI)ILAT = 2
	IF(LAT.EQ.1HF)ILAT = 3
	WRITE(IWRIT,82)A,B,C,AL,BE,GA
   82 FORMAT(//' REAL CELL DIMENSIONS'/
     1' A = ',F10.5,' B = ',F10.5,' C = ',F10.5/
     2' ALPHA = ',F10.3,' BETA = ',F10.3,' GAMMA = ',F10.3)
	AL = AL * DTR
	BE = BE * DTR
	GA = GA * DTR
  802 S = AL+BE+GA
      S=S/2.0
      V=2.0*A*B*C*SQRT(SIN(S)*SIN(S-AL)*SIN(S-
     1BE)*SIN(S-GA))
      WRITE(IWRIT,83) V, LAT
   83 FORMAT(' VOLUME = ',F10.5, ' CUBIC ANGSTROMS   LATTICE = ',A1)
      AS = B*C*SIN(AL)/V
      BS = A*C*SIN(BE)/V
      CS = A * B* SIN(GA)/V
      V = 1./V
	AL = 180.- AL/DTR
	BE = 180.- BE/DTR
	GA = 180.- GA/DTR
      WRITE(IWRIT,72)AS,BS,CS,AL,BE,GA,V
   72 FORMAT(/' RECIPROCAL CELL DIMENSIONS'/
     1' A* = ',F10.6,' B* = ',F10.6,' C* = ',F10.6/
     2' ALPHA* = ',F10.3,' BETA* = ',F10.3,'GAMMA* = ',F10.6/
     3' RECIPROCAL VOLUME = ',G20.8//)
	AL = AL * DTR
	BE = BE * DTR
	GA = GA * DTR
C      WRITE(IWRIT,92)
C   92 FORMAT(//'  H  K  L     D*        D      2*THETA')
C
C  SET LIMITS FOR DIFFERENT CRYSTAL SYSTEMS
C
C  FIRST TRICLINIC
C
	KMIN = 0
	KMAX = 8
	IF(ABS(AL-90.) .LT. .001)KMIN=-8
C
C  NOW MONOCLINIC
C
	LMIN = 0
	LMAX = 8
	IF(ABS(BE-90.) .LT. .001 .OR. ABS(GA-90.)
     1 .LT. .001)LMIN = -8
	N = 0
      DO 900 J=0,8
      DO 900 K=KMIN,KMAX
      DO900 L=LMIN,LMAX
      IF(J .EQ.0.AND.K.EQ.0.AND.L.EQ.0)GO TO 900
	GO TO(850,830,840),ILAT
C
C  TEST FOR BODY CENTERING
C
830	H = JMOD(J+K+L,2)
	IF(H.EQ.0)GO TO 850
	GO TO 900
C
C  TEST FOR F CENTERING CONDITION
C
840	H = JMOD(J,2) + JMOD(K,2) + JMOD(L,2)
	IF(H.EQ.0.OR.H.EQ.3)GO TO 850
	GO TO 900
850      DS = J*J*AS*AS + K*K*BS*BS + L*L*CS*CS +
     1 2*J*K*COS(GA)*AS*BS +
     2 2*J*L*AS*CS*COS(BE) +
     3 2*K*L*BS*CS*COS(AL)
      DS = SQRT(DS)
      D = 1.0/DS
	T2 = WAVE / (2.0 * D)
	IF(T2.GT.1.0)GO TO 900
C
C  HAS THIS D ALREADY BEEN FOUND?
C
	IF(N.EQ.0)GO TO 160
	DO 150 I=1,N
	IF(D1(I).EQ.D)GO TO 900
150	CONTINUE
160	N = N + 1
	IF(N.GT.200)GO TO 250
	IF(IPRIR.EQ.'T')GO TO 250
	D1(N) = D
	I1(N) = J
	I2(N) = K
	I3(N) = L
	T2 = ASIN( T2 ) / DTR * 2.0
	T1(N) = T2
C      WRITE(IWRIT,93)J,K,L,DS,D,T1(N)
C   93 FORMAT(3I3,3F10.5)
  900 CONTINUE
200	CONTINUE
C
C  SORT ON DECREASING d
250	DO 220 I=1,N-1
	DO 220 J=I+1,N
	IF(D1(I).GT.D1(J))GO TO 220
	T = D1(I)
	D1(I) = D1(J)
	D1(J) = T
	T = T1(I)
	T1(I) = T1(J)
	T1(J) = T
	IT = I1(I)
	I1(I) = I1(J)
	I1(J) = IT
	IT = I2(I)
	I2(I) = I2(J)
	I2(J) = IT
	IT = I3(I)
	I3(I) = I3(J)
	I3(J) = IT
220	CONTINUE
	RETURN
	END

