	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 - NiO'/' FL - Fluorophlogopite'/' FS - Fluorophlogopite + Si'
     2/' CU - Copper'/' FC - Fluorophlogopite + Cu'/
     3/'$Default <SI>): ')
	READ(IO5,1010,END=1) ISTD
1030	FORMAT(2A2)
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.'FL') CALL FL(DINT,NUM)
	IF(ISTD.EQ.'FS') CALL FS(DINT,NUM)
	IF(ISTD.EQ.'FC') CALL FC(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
		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
