      REAL CONC(200), INT(200)
      DIMENSION XT(20), YT(20), GT(20)
      DATA XT/4*' ','C','O','N','C','E','N','T','R','A',
     1'T','I','O','N',3*' '/
      DATA YT/3*' ','I','N','T','E','N','S','I','T','Y',' ',
     1'R','A','T','I','O',2*' '/
      DATA GT/3*' ','S','T','A','N','D','A','R','D',' ',
     1'C','U','R','V','E',3*' '/
      DATA NO/'N'/
C
C PROGRAM FOR QUANTITATIVE ANALYSIS BY INTERNAL STANDARD
C  METHOD ON X-RAY POWDER DIFFRACTION DATA
C
      ITYP = 108
      IREAD = 105
      WRITE(ITYP,2)
    2 FORMAT(//' PROGRAM FOR QUANTITATIVE ANALYSIS BY X-RAY POWDER'/
     1' DIFFRACTION INTERNAL STANDARD METHOD'//)
      WRITE(ITYP,3)
    3 FORMAT(' IS A NEW STANDARD CURVE TO BE DETERMINED ?'/
     1' IF NO PROGRAM WILL EXPECT A SLOPE AND Y-INTERCEPT'/
     2' FROM A PREVIOUS RUN'/)
      READ(IREAD,26)NOPE
      IF(NOPE.EQ.NO) GO TO 145
      WRITE(ITYP,4)
   4  FORMAT(' IN THIS METHOD A FIXED CONSTANT AMOUNT OF A'/
     1' SUITABLE STANDARD IS ADDED TO A SERIES OF MIXTURES OF'/
     2' KNOWN CONCENTRATION. A FULLY RESOLVED PEAK FOR THE MATERIAL'/
     3' TO BE DETERMINED AND THE INTERNAL  STANDARD ARE SELECTED.'/
     4' THE BACKGROUNDS IN COUNTS PER SEC ARE THEN DETERMINED'/
     5' ON EACH SIDE OF EACH PEAK AND THE INTEGRATED INTENSITY'/
     6' OF EACH PEAK IS DETERMINED BY SCANNING WITH THE SCALER ON.'/
     7' AND RECORDING TOTAL COUNTS AND TIME.  THIS DATA WILL BE USED'/
     8' TO DETERMINE A STANDARD CURVE FROM WHICH THE CONCENTRATION'/
     9' OF UNKNOWN SAMPLES MAY BE DETERMINED.'//)
      WRITE(ITYP,6)
    6 FORMAT(' TYPE THE VALUES FOR B1,B2,COUNTS,TIME,CONC FOR EACH '/
     1' STANDARD( INTERNAL STANDARD PEAKS COME LATER) IN 5G10.4 WHERE:'/
     2' B1 AND B2 ARE BACKGROUNS ON EACH SIDE OF THE PEAK IN CPS'/
     3' COUNTS IS THE TOTAL # OF COUNTS COLLECTED IN TIME ON THE PEAK'/
     4' CONC IS THE CONCENTRATION OF THE STANDARD'/
     5'   AFTER LAST STANDARD ENTER 0.0 FOR B1'/)
      WRITE(ITYP,8)
    8 FORMAT('          #         #         #         #         #')
      I = 1
  100 READ(IREAD,10)B1,B2,CNT,TIME,CONC(I)
   10 FORMAT(5G10.6)
      IF(B1.EQ.0.0) GOTO 110
      B = (B1 + B2)/2.0
      INT(I) =  CNT - B*TIME
      WRITE(ITYP,12) INT(I)
   12 FORMAT(' INT = ',G12.4)
      I = I + 1
      GO TO 100
  110 NUM = I - 1
      WRITE(ITYP,14)
   14 FORMAT(//' NOW TYPE B1,B2,COUNTS AND TIME FOR THE'/
     1 ' CORRESPONDING INTERNAL STANDARD PEAKS. END WITH B1 = 0.0'/)
      WRITE(ITYP,8)
      I = 1
  120 READ(IREAD,10)B1,B2,CNT,TIME
      IF(B1.EQ.0.0) GO TO130
      B = (B1 + B2)/2.0
      STD = CNT - B*TIME
      WRITE(ITYP,12) STD
      INT(I) = INT(I) / STD
      I = I + 1
      GO TO 120
  130 IF(NUM.EQ.I-1)GO TO 140
      WRITE(ITYP,16)
   16 FORMAT(//' YOU SCEWED UP THE NUMBER OF THE # OF STANDARD PEAKS'/
     1' DOES NOT EQUAL THE NUMBER OF REFRENCE PEAKS RE-TYPE REF. PEAKS')
      I = 1
      GO TO 120
  140 WRITE(ITYP,18)
   18 FORMAT(//' WE NOW WILL DO A LINEAR REGRESSION ON THE DATA TYPED'/
     1' BELOW TO GET THE STANDARD CURVE. IF THERE IS NO SYSTEMATIC'/
     2' ERROR THIS WILL BE A STRAIGHT LINE'//
     3' ************* NOTE  ********************'/
     4' PLOT THIS DATA TO MAKE SURE YOU HAVE A STRAIGHT LINE!'/
     5' ****************************************'//
     6'  I(STD)/I(REF)   CONC'/)
      WRITE(ITYP,20)(INT(I),CONC(I),I=1,NUM)
   20 FORMAT(5XG10.5,2XG10.4)
      CALL LLS(CONC,INT,NUM,SLOPE,YINCP)
      CALL TTPLOT(CONC,INT,XT,YT,GT,NUM)
       GO TO 150
  145 WRITE(ITYP,21)
   21 FORMAT(/'ENTER THE SLOPE AND Y-INTERCEPT FOR YOUR STANDARD'/
     1' CURVE IN 2G10.5'/)
      READ(IREAD,10) SLOPE,YINCP
  150 WRITE(ITYP,22)
   22 FORMAT(//' NOW TYPE B1,B2,COUNT,TIME FOR YOUR UNKNOWN'
     1' FIRST FOR YOUR UNKNOWN')
      READ(IREAD,10)B1,B2,CNT,TIME
       B = (B1 + B2) / 2.0
       STD = CNT - B*TIME
      WRITE(ITYP,12)STD
      WRITE(ITYP,220)
  220 FORMAT(/' NOW FOR THE REFRENCE')
      READ(IREAD,10)B1,B2,CNT,TIME
      B = (B1 + B2 )/ 2.0
      UNKI = CNT - B*TIME
      WRITE(ITYP,12) UNKI
      UNKI = STD/UNKI
      ANS = 100.*(UNKI - YINCP)/SLOPE
      WRITE(ITYP,24)UNKI,ANS
   24 FORMAT(//' INTEGRATED INTENSITY I(UNK)/I(REF) = 'G10.5//
     1' CONCENTRATION OF UNKNOWN = ',G10.5,' PERCENT'//
     2' DO YOU HAVE ANOTHER UNKNOWN FOR SAME CURVE?  (YES OR NO)')
      READ(IREAD,26) NOPE
   26 FORMAT(A1)
      IF ( NOPE.NE.NO)GO TO150
      WRITE(ITYP,27)
   27 FORMAT(/' DO YOU WISH TO ENTER A NEW SLOPE AND Y-INTERCEPT?'/)
      READ(IREAD,26)NOPE
      IF(NOPE.NE.NO) GO TO 145
      CALL EXIT
      END
      SUBROUTINE LLS(X,Y,NUM,SLOPE,YINCP)
      DIMENSION X(200), Y(200)
      AN=NUM
      XB=0.0
      YB=0.0
      SUMX=0.0
      SUMY=0.0
      SUMXY=0.0
      SUMYY=0.0
      SUMXX=0.0
C     XB MEANS X-BAR
C     YB MEANS Y-BAR
      DO 10 J = 1,NUM
      XB=XB+X(J)
      YB=YB+Y(J)
      SUMXX=SUMXX+X(J)*X(J)
      SUMYY=SUMYY+Y(J)*Y(J)
      SUMX=SUMX+X(J)
      SUMY=SUMY+Y(J)
   10 SUMXY=SUMXY+X(J)*Y(J)
      XB=XB/AN
      YB=YB/AN
      SLOPE=(AN*SUMXY-SUMX*SUMY)/(AN*SUMXX-SUMX*SUMX)
C     YINCP MEANS Y-INTERCEPT
      YINCP=(SUMX*SUMXY-SUMXX*SUMY)/(SUMX*SUMX-AN*SUMXX)
C     SQRS MEANS S SQUARED OR THE SIGMA OF A SINGLE MEASURMENT
      SQRS=SUMYY-SUMY*SUMY/AN
      SQRS=SQRS-((SUMXY-SUMX*SUMY/AN)**2)/(SUMXX-SUMX*SUMX/AN)
      SQRS=SQRS/(AN-2.)
      S=SQRTF(SQRS)
C     VARS MEANS THE VARIANCE FOR THE SLOPE
      VARS=SQRS/(SUMXX-AN*XB*XB)
      SS=SQRTF(VARS)
C     SINCP MEANS THE STANDARD DEVIATIONOF THE Y-INTERCEPT
      SINCP=SQRTF(SQRS*SUMXX/(AN*SUMXX-SUMX*SUMX))
      T=YINCP/SINCP
      NF = NUM - 2
      WRITE(108,100) YINCP,SLOPE
  100 FORMAT(//' THE LEAST SQUARES EQUATION IS:',/'  YC =',G14.8,' + ',G
     114.8,' X',//)
      WRITE(108,110) SINCP,SS,S,SQRS,VARS,T,NF
  110 FORMAT(' STATISTICS OF THE STRAIGHT LINE:'/
     1' THE STANDARD DEVIATION OF THE Y INTERCEPT IS =',G14.8/
     2' THE STANDARD DEVIATION OF THE SLOPE IS =',G14.8/
     3' THE STANDARD DEVIATION OF A SINGLE Y VALUE IS =',G14.8/
     4' THE ESTAMITE OF VARIANCE FOR A SINGLE Y VALUE IS =',G14.8/
     5' THE VARIANCE OF THE SLOPE IS =',G14.8/
     6' THE STUDENT ''T'' VALUE IS =',G14.8/
     7' COMPARE THIS VALUE ON A ''T'' TABLE AT ',I5,' DEGREES OF FREEDOM
     8')
      TOTAL = 0.0
      R = 0.0
      WRITE(108,120)
  120 FORMAT(//7X,'X',15X,'Y',14X,'YC',12X,'YC-Y'/)
      DO 20 J=1,NUM
      YCAL = YINCP + SLOPE * X(J)
      DIFF = YCAL - Y(J)
      R = R + DIFF*DIFF
      TOTAL = TOTAL + ABS(DIFF)
      WRITE(108,130) X(J), Y(J), YCAL, DIFF
  130 FORMAT(4(G14.8,2X))
   20 CONTINUE
   40 TOTAL = TOTAL / AN
      WRITE(108,140) TOTAL,R
  140 FORMAT(//' THE AVERAGE OF THE ABSOLUTE VALUES OF YC-Y IS =',G14.8
     1/' THE SUM OF THE SQUARES OF THE RESIDUALS R ='G14.8/1H1)
  150 RETURN
      END
      SUBROUTINE TTPLOT(X,Y,XTITLE,YTITLE,GRAPH,NPTS)
      DIMENSION PRINT(63),X(201),Y(201),XTITLE(20),
     1YTITLE(20),GRAPH(20)
      XMAX=X(1)
      XMIN=X(1)
      YMAX=Y(1)
      YMIN=Y(1)
      DO 601 J=2,NPTS
      IF(X(J).GT.XMAX)XMAX=X(J)
      IF(X(J).LT.XMIN)XMIN=X(J)
      IF(Y(J).GT.YMAX)YMAX=Y(J)
      IF(Y(J).LT.YMIN)YMIN=Y(J)
601   CONTINUE
      RANGEX=ABS(XMAX-XMIN)
      RANGEY=ABS(YMAX-YMIN)
      XMAX=XMAX+(.025*RANGEX)
      XMIN=XMIN-(.025*RANGEX)
      YMAX=YMAX+(.025*RANGEY)
      YMIN=YMIN-(.025*RANGEY)
      XSF=(XMAX-XMIN)/60.0
      YSF=(YMAX-YMIN)/40.0
      DATA DOT/'.'/
      DATA AST/'*'/
      DATA PLS/'+'/
      DATA SPA/' '/
      DO 9000 J=3,53,10
      DO 9001 I=1,9
9001  PRINT(I+J)=DOT
9000  PRINT(J)=PLS
      PRINT(1)=SPA
      ITIMES=0
      IA=4
      PRINT 9
      DO 650 J=40,1,-1
      IP=J+1
      YSCALE=J*YSF+YMIN
      PRINT(3)=DOT
      IF(J.EQ.(IA*10))PRINT(2)=YSCALE
      IF(J.EQ.(IA*10))PRINT(3)=PLS
      IF(J.EQ.(IA*10))PRINT(63)=PLS
      IF(J.LE.30.AND.J.GE.10)ITIMES = ITIMES+1
      IF(ITIMES.GT.0.AND.ITIMES.LE.20)PRINT(1)=YTITLE(ITIMES)
      DO 606 I=1,NPTS
      IF((INT(((Y(I)-YMIN)/YSF)+1.5)).EQ.IP) GO TO 30
      GO TO 606
30    IN=INT(((X(I)-XMIN)/XSF)+3.5)
      PRINT(IN)=AST
606   CONTINUE
      IF(J.EQ.(IA*10))GO TO 31
      PRINT 5,(PRINT(I),I=1,63)
      GO TO 32
31    PRINT 4,(PRINT(I),I=1,63)
32    DO 8000 IT=1,63
8000  PRINT(IT)=SPA
      IF(J.EQ.(IA*10))IA=IA-1
650   CONTINUE
      PRINT(2)=YMIN
      PRINT(3)=PLS
      DO 6000 J=3,53,10
      DO 6001 I=1,9
6001  PRINT(I+J)=DOT
6000  PRINT(J)=PLS
      PRINT(63)=PLS
      DO 7000 I=1,NPTS
      IF((INT(((Y(I)-YMIN)/YSF)+1.5)).EQ.1)  GO TO 300
      GO TO 7000
300   IN=INT(((X(I)-XMIN)/XSF)+3.5)
      PRINT(IN)=AST
7000  CONTINUE
      PRINT 4,(PRINT(I),I=1,63)
      PRINT(1)=XMIN
      DO 7001 J=2,7
7001  PRINT(J)=(J-1)*10*XSF+XMIN
      PRINT 6,(PRINT(I),I=1,7)
      PRINT 7,(XTITLE(I),I=1,20)
      PRINT 8,(GRAPH(I),I=1,20)
      WRITE(108,10)XSF,YSF
      PRINT 9
1     FORMAT(I3)
2     FORMAT(2G)
3     FORMAT(20A1)
4     FORMAT(1X,A1,1X,G8.2,61A1)
5     FORMAT(1X,A1,7X,A1,1X,61A1)
6     FORMAT(3X,E9.3,6(X,E9.3))
7     FORMAT(/,32X,20A1)
8     FORMAT(//,32X,20A1)
9     FORMAT(///,72('-'),///)
10    FORMAT(//,X,'EACH X INCREMENT REPRESENTS',G10.5,
     1'UNITS',/,1X,'EACH Y INCREMENT REPRESENTS',G10.5,
     2'UNITS')
      RETURN
      END
