C*********************************************************************
C
C  PROGRAM SPIKE TO ANALYSE AUTO QUANTITATIVE ANALYSIS
C  FILES FOR THE SPIKING METHOD
C
C          WRITTEN BY
C
C  R. L. SNYDER   NYS COLLEGE OF CERAMICS
C                 ALFRED UNIVERSITY
C                 ALFRED, NEW YORK 14802
C
C and
C
C  C.R. HUBBARD   NATIONAL BUREAU OF STANDARDS
C                 A221 MATL
C                 WASHINGTON D.C. 20234
C
C
C  READS FILES CORRESPONDING TO THE ANALYTE PATTERN PLUS TWO OR MORE
C  PATTERNS OF THE ANALYTE SAMPLE SPIKED WITH KNOWN QUANTITIES OF
C  THE ANALYTE PHASE.
C
C
C
C  DATA FILE STRUCTURE IS DEFINED IN "AUTO - A REAL TIME DIFFRACTOMETER
C  CONTROL SYSTEM", R.L. SNYDER, C.R. HUBBARD AND N. PANAGIOTOPOLOUS
C  NBS PUBLICATION NBSIR 81-2229 (1981)
C***********************************************************************
C
C  THE FOLLOWING PARAMETER STATMENT DEFINES THE DIMENSIONS FOR
C  THE QUANTATATIVE ANALYSIS SUBROUTINES.  TO REDIMENSION
C  SIMPLY CHANGE THIS STATMENT HERE AND IN EACH SUBROUTINE.
C
C      NF = THE NUMBER OF DATA FILES ALLOWED FOR INPUT
C           (EACH REPEAT AND MOUNTING IS CONSIDERED A SEPARATE FILE)
C      NR = THE MAXIMUM NUMBER OF REFLECTIONS WHICH CAN BE PROCESSED
C      NP = THE MAXIMUM NUMBER OF PHASES ALLOWED
C      NT = THE MAXIMUM NUMBER OF DATA POINTS FOR ANY REFLECTION
C      NS = THE MAXIMUM NUMBER OF LINES FOR THE STANDARD PHASE
C      NU = THE MAXIMUM NUMBER OF LINES FOR ALL UNKNOWN PHASES
C
	PARAMETER  (NF=30, NR=100, NP=10, NT=1000, NS=5, NU=20)
	COMMON/IO/ IRD,ITYP,IPRT,ISCRT,ICOBRG,LIMPRT,IDBUG
	COMMON /FILE/ LU,ITYPE,IEND
	COMMON /FILEC/ NAME
	COMMON /RECRD1/ IFLAG(10)
	COMMON /RECR1C/ TITLE,DATE,TIME,RNFLNM
	COMMON /RECRD2/ XLAMBD,TAKEOF,DIVANG,RADIUS,SOLLIN,
     1 SOLLOT,RECSLT,DEADTM,STABIL,TBASE
	COMMON /RECR2C/ UNIT,TARGET,LINE,FOCAL,DIVSLT,MONOCH,DETECT
	COMMON /RECRD6/ XWT(10),SIG(10)
	COMMON /RECR6C/ ELEMNT(10)
	COMMON /RECRD8/ FRELOV(2),SRELOV(2),WTPCPH(8)
	COMMON /RECR8C/ NAMEOV(2),JKLOV(2),ELMTPH(8)
	COMMON/RECRD7/BANG,PANG,EANG,DANG,RATE,MVCNT,TIMLIM,FIXERR,
     * LIMCNT,KODE,IOVLP,IPK,NPK,NBG,FREL,STDREL,STD,SIGSTD,
     * WTPRCT,PKTIME,BKGTIM,IOPTIM,ILOC,CONSTD
	COMMON /RECR7C/ NAM,JKL
	COMMON/DATA/ PEAK,VPEAK,AREA,VAREA,XBKG,VBKG
	COMMON/SVDAT1/ AINT(NR,NF),VINT(NR,NF),NREFL,IFILE
	COMMON/LABELS/ PHASE(NR),HKL(NR)
	COMMON/SVDAT2/ FRELS(NR),SRELS(NR),WT(NF),KODES(NR),STDS(NR),
     * SSTDS(NR),IOVS(NR),IPKS(NR),RELOV(2,NR),
     * RELOVS(2,NR), WTP(8,NR), CONC(NF), XMURHO(NR)
	COMMON /SVDA2C/ MAN(2,NR),KHL(2,NR),ELM(8,NR)
	DIMENSION COMP(NS,NU),VCOMP(NS,NU),SUM(NU),VSUM(NU),
     1 WSUM(NU),L(NU),Y(NF),X(NF),V(NF)
	CHARACTER*6 PHASE,HKL,ALL
	CHARACTER*100 TITLE
	CHARACTER*6 DATE,TIME,MAN,JKL,NAMEOV,JKLOV,NAM,KHL,RNFLNM
	CHARACTER*20 NAME
	CHARACTER*2 ELMTPH, ELM, ELEMNT
	CHARACTER*4 UNIT
	CHARACTER*2 TARGET,LINE,FOCAL,DIVSLT,MONOCH,DETECT
	CHARACTER*8 DAT, TIM
	CHARACTER*10 ESDHDG(2)
       CHARACTER FILNAM*40
	DATA ALL/'  all '/
	ESDHDG(1) = ' Weighted '
	ESDHDG(2) = 'Unweighted'
C
C-----LOGICAL UNIT# ASSIGNMENTS
C     IRD IS FOR CRT TERMINAL INPUT OR CARD READER INPUT
C     ITYP IS FOR CRT TERMINAL OUTPUT (REQUESTING USER RESPONSES)
C     IPRT IS FOR 132 COLUMN PRINTER OUTPUT (FINAL RESULTS)
C     LU   IS FOR AUTO DATA FILE INPUT (OPENED BY AUTOLU SUBROUTINE)
C     ICOBRG IS NOT USED IN QUANT82*SPIKE
C     ISCRT  IS NOT USED IN QUANT82*SPIKE
C
	IRD  = 5
	ITYP = 5
	IPRT = 5
C-----INITIALIZE  IDBUG (=1 FOR DEBUG PRINT OUT, =0 OTHERWISE)
	IDBUG=0
	NRFL = 0
	CALL DATIM(DAT,TIM)
C
C-----PRINT SIGN ON MESSAGE
C
	WRITE(ITYP,10)DAT,TIM
10	FORMAT(/17X,'NBS*QUANT84-SPIKE',5X,'Revision Date:',
     & ' 29 FEB 1984'/
     * /' ',18X,'Quantitative Analysis by the Spiking Method'/
     */' ',19X,'Written by R. L. Snyder and C. R. Hubbard '/
     * /' ',24X,'Run on ',A8,5X,'at ',A8//)
C  
C  SELECT OUTPUT DESTINATION
C
        WRITE(ITYP,1)
1	FORMAT(/' Select Output Mode Desired:  (Output is 132 characte',
     1'rs wide)'/'       1 	- causes output to be put on the inter',
     2'active terminal.'/'       newname 	- causes output to go ',
     3'to the named file.'/'       <CR> 	- causes output to go ',
     4'to a file named SPIKE.LIS.'/' Enter the output destination: '$)
       READ(IRD,2) FILNAM
2	FORMAT(A40)
       IF(FILNAM(1:1) .EQ. '1') GO TO 3
       IPRT = 6
       IF(FILNAM(1:1) .EQ. ' ') THEN
		OPEN(UNIT=IPRT,FILE='SPIKE.LIS',STATUS='NEW',
     1 ACCESS='SEQUENTIAL',FORM='FORMATTED')
       		GO TO 3
       ENDIF
	OPEN(UNIT=IPRT,FILE=FILNAM,STATUS='NEW',ACCESS='SEQUENTIAL',
     1 FORM='FORMATTED')
3	LU    = 10
	IF(ITYP.NE.IPRT) WRITE(IPRT,10)DAT,TIM
C-----LIMPRT CONTROL RAW DATA FILE PRINTING AND HEADINGS
C           =0 FOR FULL LISTING
C           =1 FOR HEADING ON FIRST FILE ONLY
C           =2 FOR BRIEFEST RAW DATA LISTING (FILE TITLE ONLY)
	WRITE(ITYP,11)
11	FORMAT(/' Select a print control option (I1):'/
     & ' 0 = Full listing of data files'/
     & ' 1 = Heading on first data file only'/
     & ' 2 = Data file titles only')
	READ(IRD,12) LIMPRT
12	FORMAT(I1)
	IF(LIMPRT.LT.0 .OR. LIMPRT.GT.3) LIMPRT=0
C
C-----ALLOW FOR USE OF UNIT WEIGHTS OR COUNTING STATISTIC WEIGHTS
C-----IN FORMING AVERAGES OF I-REL AND RIR VALUES AS WELL AS RMS
C-----DEVIATIONS.
C
13	WRITE(ITYP,1013)
1013	FORMAT(/' ','Enter weighting choice:'/
     1       ' 0 for counting statistical weights'/
     2       ' 1 for unit weights')
	READ(IRD,*) IWT
	IF( IWT .NE.0 .AND. IWT .NE.1) THEN
	  WRITE(ITYP,1014) IWT
1014	FORMAT(/' ','Input option of ',I3,' not equal to 0 or 1')
	  GOTO 13
	ENDIF
C-------------------------------------------------------------------
C-----READ FILES AND PROCESS RAW DATA INORDER TO DETERMINE THE
C-----INTENSITIES OF LINES IN THE ANALYTE
C-------------------------------------------------------------------
	ITYPE = 3
	CALL QNDATA
	IF(IEND.EQ.1)GO TO 999
C
C-----ALL DATA HAS BEEN READ, PROCESSED AND STORED INTO ARRAYS
C
C********************************************************************
C-----BEGIN CALCULATION OF THE INDIVIDUAL INTENSITY RATIOS AND VARIANCES
C**********************************************************************
C
C-----FIND THE NEXT REFERENCE REFLECTION
C
	WRITE(IPRT,2000) ESDHDG(IWT+1),IFILE
2000	FORMAT(///' Line by Line Composition Analysis ',
     1'from ',A10,   ' Linear Regressions of',i3,
     2' points per analysis')
	WRITE(IPRT,2020)
2020	FORMAT(/' ','Analyte  Line  Reference Line   Weight  ',
     1' Least Sqs'/' ','Phase   h k l   Name    h k l  ',
     1'Percent  est Std Dev',8x,'Slope',12X,'Intercept',6x,
     2 'correlation')
C
C-----PICK OUT THE UNOVERLAPPED REFERENCE LINES (KODES=2)
C
	NSTD=0
	DO 300 I=1,NREFL
	IF(KODES(I).NE.2) GOTO 300
	IF(IOVS(I).EQ.0)GO TO 38
C
C  IF A REFERENCE LINE IS OVERLAPPED PRINT WARNING
C
	WRITE(ITYP,2050)PHASE(I),HKL(I)
	IF(IPRT.NE.ITYP) WRITE(IPRT,2050) PHASE(I),HKL(I)
2050	FORMAT(/' Reference line ',2A6,' is labled as overlapped',
     1' we will assume background has been properly ',
     2'measured and continue')
C
C-----FOUND THE NEXT REFERENCE REFLECTION.  USE IT TO CALCULATE THE
C-----RATIO FOR ALL COMBINATIONS OF THIS REFERENCE REFLECTION WITH UP TO
C-----NU    REFLECTIONS FROM THE "ANALYTE"
C
38	NSTD=NSTD+1
	IF(NSTD.LE.NS) GOTO 50
	WRITE(ITYP,40) NS
	IF(IPRT.NE.ITYP) WRITE(IPRT,40) NS
40	FORMAT(/' # of reference lines exceeds array limit NS=',I3,
     * ' ignoring the extra lines')
	NSTD=NSTD-1
	GOTO 300
C
C-----PICK OUT THE ANALYTE LINES (KODES=1)
C
50	NRFL = 0
	DO 200 J=1,NREFL
	IF(KODES(J).NE.1) GOTO 200
C
C-----FOUND A SAMPLE REFLECTION.  NOW CALCULATE THE COMP FOR ALL OBSERVA
C-----OF THE PAIR (I,J). (IE FOR ALL FILES)
C
	NRFL=NRFL+1
	IF(NRFL.LE.NU) GOTO 70
	WRITE(ITYP,60) NU
	IF(IPRT.NE.ITYP) WRITE(ITYP,60) NU
60	FORMAT(/' # of analyte lines exceeds array limit NU=',I3,
     * ' ignoring the extra line')
	NRFL=NRFL-1
	GOTO 200
C
C  CALCULATE THE ANALYTE COMPOSITION FROM EACH MEASURED LINE
C
70	L(NRFL) = J
	IF(IOVS(J).EQ.0)GO TO 75
	WRITE(ITYP,2100)PHASE(J),HKL(J)
	IF(IPRT.NE.ITYP) WRITE(IPRT,2100) PHASE(J),HKL(J)
2100	FORMAT(/' Reflection ',2A6,' is overlapped and will not be',
     1' used in computation')
	GOTO 200
75	SUMX = 0.0
	SUMY = 0.0
	SUMX2= 0.0
	SUMY2= 0.0
	DO 100 K=1,IFILE
C-----Y IS TO BE THE INTENSITY OF ANALYTE LINE
	Y(K)=AINT(J,K)/AINT(I,K)
C-----V IS TO BE THE VARIANCE (SIGMA**2) OF THE Y VALUES (IWT=0)
	V(K) = 1.0
	IF(IWT.EQ.0) V(K)=Y(K)*Y(K)*(VINT(J,K)/(AINT(J,K)*AINT(J,K))+
     * VINT(I,K)/(AINT(I,K)*AINT(I,K)))
C
C-----X IS TO BE GRAMS OF ANALYTE ADDED PER GRAM OF SAMPLE
C-----WT(J) IS THE WEIGHT FRACTION OF THE ADDED ANALYTE
	X(K) = WT(K) / (1.0 - WT(K))
	SUMX = SUMX + X(K)
	SUMX2= SUMX2+ X(K)*X(K)
	SUMY = SUMY + Y(K)
	SUMY2= SUMY2+ Y(K)*Y(K)
100	CONTINUE
	IF(IDBUG.EQ.1) WRITE(IPRT,102) (X(N),Y(N),V(N),N=1,IFILE)
102	FORMAT(/'   gm added   I/Iref    Weight'/(' ',3F10.3))
C-----CALL SOLVE TO CALCULATE WT% OF THE ANALYTE FOR EACH OF THE N-LINES
	CALL SOLVE(X,Y,V,IFILE,XM,XB,VM,VB)
C-----SINCE AT Y=0 X=-XB/XM AND WT%=-X THEN WT%=+XB/XM
	COMP(NSTD,NRFL) = 100. * XB / XM
	VCOMP(NSTD,NRFL) = COMP(NSTD,NRFL)*COMP(NSTD,NRFL) *
     1 (VB / (XB*XB) + VM / (XM*XM))
	VV = SQRT(VCOMP(NSTD,NRFL))
	SM = SQRT(VM)
	SB = SQRT(VB)
	VARX = SUMX2/IFILE - (SUMX/IFILE)**2
	VARY = SUMY2/IFILE - (SUMY/IFILE)**2
	CORR = XM * SQRT(VARX) / SQRT(VARY)
	WRITE(IPRT,110)PHASE(J),HKL(J),PHASE(I),HKL(I),COMP(NSTD,NRFL),VV,
     * XM,SM,XB,SB,CORR
110	FORMAT(/' ',A6,1X,A6,2X,2(1X,A6),2(2X,F8.2),5X,2(F10.4,'(',F6.4,
     * ')'),F10.3)
C
C-----OUTPUT OBS, CALC, DIFFERENCE INTENSITY RATIOS FOR EXAMINATION
C
	WRITE(IPRT,114)
114	FORMAT(/' ',51X,'gms added      Obs       esd       Calc'
     * '     Delta'/
     *' ',53X, 'per gm ',4(5X,'Ia/Ir'))
	DO 120 K = 1,IFILE
	YCALC = XM*X(K) + XB
	DIFF = YCALC-Y(K)
	SIGV = SQRT(V(K))
	WRITE(IPRT,118) X(K),Y(K),SIGV,YCALC,DIFF
118	FORMAT(' ',50X,5F10.4)
120	CONTINUE
200	CONTINUE
300	CONTINUE
C
C-----CHECK THAT AT LEAST 1 REFERENCE AND AT LEAST 1 SAMPLE LINE WAS INP
C
	IF(NSTD.NE.0 .AND. NRFL.NE.0) GOTO 400
	WRITE(ITYP,320) NSTD,NRFL
	IF(IPRT.NE.ITYP) WRITE(IPRT,320) NSTD,NRFL
320	FORMAT(/' **Input Error** #reference lines=',i3,3x,
     * '# of sample lines=',i3)
	GOTO 999
C
C-----CALCULATE THE AVERAGE COMPOSITION FOR EACH ANALYTE LINE
C
400	WRITE(IPRT,1020)
1020	FORMAT(///' Average and weighted average for each analyte lines:')
	WRITE(IPRT,2070)
2070	FORMAT(/' ','Analyte  Line  Reference Line   Average  Root Mean',
     1 '  Waverage   Weighted  Weighted     # of'/
     1 ' ','Phase   h k l   Name    h k l  ',
     1'   Wt %  Square Dev.    Wt %     e.s.d.   Rms Dev.   Analyses')
	DO 420 I=1,NRFL
	SUM(I) = 0.0
	WSUM(I) = 0.0
	VSUM(I) = 0.0
	DO 410 J=1,NSTD
	SUM(I) = SUM(I) + COMP(J,I)
	WSUM(I)=WSUM(I) + COMP(J,I) / VCOMP(J,I)
	VSUM(I)   = VSUM(I)   + 1.0 / VCOMP(J,I)
410	CONTINUE
C-----CALCULATE THE AVERAGE, WEIGHTED AVERAGE AND ESD
	SUM(I)  = SUM(I) / NSTD
	WSUM(I) = WSUM(I) / VSUM(I)
	VSUM(I) = 1.0 / VSUM(I)
	VV = SQRT(VSUM(I))
	LLL = L(I)
C-----CALCULATED THE WEIGHTED RMS DEVIATION FROM WEIGHTED AVERAGE
	RMSD = 0.0
	WRMSD= 0.0
	IF(NSTD.EQ.1) GOTO 415
	DO 414 J=1,NSTD
	RMSD = RMSD + (SUM(I)-COMP(J,I))**2
414	WRMSD=WRMSD + (SUM(I)-COMP(J,I))**2/VCOMP(J,I)
	RMSD = SQRT(RMSD/(NSTD-1))
	WRMSD= SQRT(WRMSD*VSUM(I)/(NSTD-1))
415	WRITE(IPRT,416)PHASE(LLL),HKL(LLL),ALL,ALL,SUM(I),RMSD,WSUM(I),VV,
     * WRMSD,NSTD
416	FORMAT(2(1X,A6),2X,2(1X,A6),5(2X,F8.2),I10)
420	CONTINUE
C
C  LASTLY AVERAGE OVER THE LINES OF EACH PHASE
C
	WRITE(IPRT,1040)
1040	FORMAT(///' Overall average and weighted average:')
	WRITE(IPRT,2070)
	ICNT = 0
	TSUM = 0.0
	TWSUM= 0.0
	VTSUM = 0.0
	DO 440 J=1,NRFL
	ICNT = ICNT + NSTD
	TSUM = TSUM + SUM(J)
	TWSUM= TWSUM+ WSUM(J) / VSUM(J)
	VTSUM = VTSUM + 1.0 / VSUM(J)
440	CONTINUE
	TSUM = TSUM / NRFL
	TWSUM = TWSUM / VTSUM
	VTSUM = 1.0 / VTSUM
	VTSUM = SQRT(VTSUM)
C-----CALCULATE THE OVERALL UNWEIGHTED AND WEIGHTED RMS DEVIATION
	RMSD = 0.0
	WRMSD = 0.0
	IF (ICNT.EQ.1) GOTO 450
	DEN = 0.0
	DO 448 I=1,NRFL
	DO 448 J=1,NSTD
	RMSD = RMSD + (TSUM-COMP(J,I))**2
	WRMSD = WRMSD + (TSUM-COMP(J,I))**2/VCOMP(J,I)
448	DEN = DEN + 1.0/VCOMP(J,I)
	RMSD = SQRT(RMSD/(ICNT-1))
	WRMSD = SQRT(WRMSD/((ICNT-1)*DEN))
450	WRITE(IPRT,416)PHASE(L(1)),ALL,ALL,ALL,TSUM,RMSD,TWSUM,VTSUM,
     * WRMSD,ICNT
999	CALL EXIT
	END
