C PROGRAM TO DETERMINE PRECISE LATTICE PARAMETERS
C FOR TETRAGONAL AND HEXAGONAL CRYSTAL SYSTEMS
C FROM POWDER DIFFRACTION DATA
      IMPLICIT REAL*8 (A-H)
      IMPLICIT REAL*8 (O-Z)
      REAL*4 XT, YT, GT, XCOS(5), YCOS(5)
      INTEGER TEST,TSET,HEX,HOG,YE,YES
      DIMENSION THET(50),J(50),K(50),L(50),AP(50),CP(50)
      DIMENSION D(50),COST(50),IMP(80),NAM1(9)
      DIMENSION JKL(2,3),T1(2),D1(2),DS1(2),DSS1(2)
      DIMENSION XT(20), YT(20), GT(20)
      DATA XT/4*' ','C','O','S','*','*','2','(','T','H','E','T',
     1'A',')',3*' '/
      DATA YT/3*' ','C','E','L','L',' ','P','A','R',
     1'A','M','E','T','E','R',3*' '/
      DATA GT/4*' ','E','X','T','R','A','P','O','L','A','T',
     1'I','O','N',3*' '/
      DATA TSET/'D'/,HOG/'H'/,YE/'Y'/
      DATA NAM1/'D','=','H','+','K','A','+','L','C'/
      YES = YE
C DEFINE I/O UNITS
      ITYP = 6
      IREAD = 5
      RAD = 6.283184/360.
      WRITE(ITYP,1)
    1 FORMAT(/' PROGRAM TO DETERMINE PRECISION  LATTICE'/
     1'  PARAMETERS FOR TETRAGONAL AND HEXAGONAL SYSTEMS'/
     2'  USING AN ITERATIVE EXTRAPOLATION METHOD'///
     3' TO BEGIN WE MUST OBTAIN CRUDE VALUES OF A AND C TO'/
     4' DETERMINE THE C/A RATIO.  TO DO THIS YOU MUST'/
     5' CHOOSE TWO HIGH ORDER REFLECTIONS WITH RELATIVELY'/
     6' HIGH H(OR K) AND L INDICES.  WHAT EQUATION RELATES'/
     7' THE UNKNOWN A AND C VALUES TO KNOWN PARAMETERS?')
      IC = 0
  100 READ(IREAD,3)IMP
      IC = IC + 1
    3 FORMAT(80A1)
      M = 1
      DO120 I = 1,80
      IF(IMP(I) - NAM1(M))120,110,120
  110 M = M + 1
      IF (M.EQ.10) GO TO 130
  120 CONTINUE
      WRITE(ITYP,5)
    5 FORMAT(//' NO DUMMY!! THE EQ. STARTS WITH D*(SQ) = ...'/
     1'  TRY AGAIN'/)
      IF(IC.GT.3) GO TO 121
      GO TO 100
  121 WRITE(ITYP,6)
    6 FORMAT(/' THE EQUATION ENDS WITH :'/
     1' = (H(SQ) + K(SQ))*A*(SQ) + L(SQ)*C*(SQ)'/
     2' NOW SEE IF YOU CAN AT LEAST TYPE IT ')
      GO TO 100
  130 WRITE(ITYP,7)
    7 FORMAT(//' CORRECT! WE WILL USE TWO EQUATIONS OF THIS'/
     1' TYPE AND SOLVE SIMULTAENOUSLY FOR THE CRUDE A AND C'/
     2' VALUES.  BUT FIRST I MUST KNOW IF YOU WISH TO ENTER'/
     3' D OR TWO THETA VALUES?'/)
      READ(IREAD,3)TEST
      WRITE(ITYP,9)
    9 FORMAT(/' OK BUT I  ALSO NEED LAMBDA (G10.6)'/)
      READ(IREAD,11)WAVE
   11 FORMAT(G10.6)
      WRITE(ITYP,13)
   13 FORMAT(/' THAT LOOKS LIKE A REASONABLE LAMBDA SO I NOW'/
     1' NEED H,K,L AND TWO THETA( OR D) FOR TWO REFLECTIONS'/
     2' TYPE VALUES IN 3I3,G10.3 ONE REFLECTION AT A TIME')
      READ(IREAD,15)((JKL(I,1),JKL(I,2),JKL(I,3),T1(I)),I=1,2)
   15 FORMAT(3I3,G10.4)
      WRITE(ITYP,14)
   14 FORMAT(15X,'  H  K  L',2X,'2-THETA',7X,'D',7X,'D*',6X,'D*(SQ)')
      DO 170 I = 1,2
      IF(TEST-TSET)150,140,150
  140 D1(I) = T1(I)
      T1(I)    = WAVE/(2.0*D1(I))
      T1(I)    = 2.0*DATAN(T1(I)/DSQRT(1-T1(I)**2))/RAD
      GO TO  160
  150 D1(I) = WAVE/(2.*SIN(T1(I)*RAD/2))
  160 DS1(I) = 1.0/D1(I)
      DSS1(I) = DS1(I)*DS1(I)
  170 WRITE(ITYP,17)I,(JKL(I,JJ),JJ=1,3),T1(I),D1(I),DS1(I),DSS1(I)
   17 FORMAT(' REFLECTION',I2,2X3I3,2G10.4,2G13.7)
      WRITE(ITYP,19)
   19 FORMAT(//' BEFORE I CAN CALCULATE A AND C FROM A* AND C* I MUST'/
     1' KNOW WEATHER THE SYSTEM IS TETRAGONAL OR HEXAGONAL BECAUSE'/
     2' IN HEX A=2/(SQRT(3)*A*)  WHILE IN TETRAG. A=1/A*'/
     3' WHICH IS IT?')
      READ(IREAD,3)HEX
      B1 = JKL(1,3)**2
      B2 = JKL(2,3)**2
      IF(HEX-HOG)173,175,173
  173 A1 = JKL(1,1)**2 + JKL(1,2)**2
      A2 = JKL(2,1)**2 + JKL(2,2)**2
      GO TO 178
  175 A1 = 4.0*(JKL(1,1)**2 + JKL(1,1)*JKL(2,1) +
     1 JKL(2,1)**2)/3.0
      A2 = 4.0*(JKL(2,1)**2 + JKL(2,1)*JKL(2,2) +
     1 JKL(2,2)**2)/3.0
  178 CALL SIMULT(A1,B1,DSS1(1),A2,B2,DSS1(2),ASS,CSS)
      ICYCLE = 0
      IF(HEX-HOG)180,190,180
C CALCULATE A AND C FOR TETRAGONAL
  180 A = 1.0/DSQRT(ASS)
      C = 1.0/DSQRT(CSS)
      GO TO 200
  190 A = 2.0/(DSQRT(ASS*3.0))
      C = 1.0/DSQRT(CSS)
C
C BEGIN EXTRAPOLATION LOOP
C
  200 ICYCLE = ICYCLE + 1
      CONA = C/A
      IF(ICYCLE-1)210,210,220
  210 WRITE(ITYP,21)
   21 FORMAT(//' THE CRUDE VALUES FOR A,C AND C/A OBTAINED FROM'/
     1' SIMULTAENOUS SOLUTION OF TWO EQUATIONS ARE:'/)
      GO TO 230
  220 WRITE(ITYP,23)ICYCLE
   23 FORMAT(//' THE REFINED VALUES FOR A AND C BEFORE CYCLE ',I2/
     1' ARE LISTED BELOW IF ANALYSIS OF THE STANDARD DEVIATIONS SHOW NO'
     2/' SIGNIFICANT IMPROVEMENT THEN YOU ARE FINISHED IF NOT DO'/
     3' ANOTHER CYCLE'/)
  230  WRITE (ITYP,25)A,C,CONA
   25 FORMAT(' A = ',G12.5,5X,'C = ',G12.5,5X,'C/A = ',G13.6/)
      IF(ICYCLE-1)240,240,290
  240 WRITE(ITYP,27)
   27 FORMAT(' NOW BASED ON THE PSUDO-INDEPENDANT EQUATIONS RELATING'/
     2' FIRST A AND THEN C TO H,K,L,D AND C/A FOR YOUR PARTICULAR'/
     3' CRYSTAL SYSTEM I WILL CALCULATE AN A AND C VALUE FOR EACH '/
     4' REFLECTION YOU INPUT.  TYPE H,K,L AND 2-THETA(OR D)'/
     5' IN 3I3,G10.4 ONE ON A LINE. AFTER LAST VALUES ENTER ALL 0''S'/)
      I = 1
  250 READ(IREAD,15)J(I),K(I),L(I),THET(I)
      IF(THET(I).EQ.0.0) GO TO 300
      IF(TEST-TSET)270,260,270
  260 D(I) = THET(I)
      THET(I)  = WAVE/(2.0*D(I))
      THET(I)  = 2.0*DATAN(THET(I)/DSQRT(1-THET(I)**2))/RAD
      GO TO 280
  270 D(I) = WAVE/(2.*SIN(THET(I)*RAD/2))
  280 COST(I) = DCOS(THET(I)*RAD/2.0)**2
      NUM = I
      I = I + 1
      GO TO 250
  290 WRITE(ITYP,29)
   29 FORMAT(/' DO YOU WISH ANOTHER CYCLE? YES OR NO')
      READ(IREAD,3)YES
      IF(YES.NE.YE) GO TO 900
      WRITE(ITYP,30)
   30 FORMAT(/' DO YOU WANT FULL RESULTS OF LEAST SQUARES'/
     1' LISTED? YES OR NO...BE SURE TO LIST LAST CYCLE')
      READ(IREAD,3)YES
  300 CONA = CONA**2
      WRITE(ITYP,31)
   31 FORMAT(//'CALCULATED VALUES FOR A AND C AND COS**2(THETA):'//
     1'  H  K  L',2X,'2-THETA',6X,'D',7X,'A CALC',4X,'C CALC',2X,'COS**2(THETA)
     2(THETA)')
      DO 400 I = 1,NUM
      LS = L(I) * L(I)
      IF(HEX-HOG)310,320,310
  310 JS = J(I)**2 + K(I)**2
      AP(I) = D(I) * DSQRT(JS +LS/CONA)
      CP(I) = D(I) * DSQRT(JS*CONA + LS)
      GO TO 400
  320 JS = 4./3.*(J(I)**2 + J(I)*K(I) +K(I)**2)
      AP(I) = D(I) * DSQRT(JS + LS/CONA)
      CP(I) = D(I) * DSQRT(JS*CONA + LS)
  400 WRITE(ITYP,33)J(I),K(I),L(I),THET(I),D(I),AP(I),CP(I),COST(I)
   33 FORMAT(3I3,4G10.4,G13.7)
      IF(ICYCLE-1)410,410,420
  410 WRITE(ITYP,34)
  34  FORMAT(//' WE ARE NOW READY TO GET BETTER VALUES FOR A AND C'/
     1 ' BY CARRING OUT A LEAST SQUARES REGRESSION ON FIRST A AND THEN'/
     2 ' C VS. COS**2(THETA). THE EXTRAPOLATED Y INTERCEPT WILL BE'/
     3 ' THE IMPROVED VALUE')
  420 WRITE(ITYP,35)
   35 FORMAT(//' THE FOLLOWING APPLY TO THE EXTRAPOLATION OF A'/
     1' USING ONLY THE 5 REFLECTIONS WITH HIGHEST H INDICES'/)
      CALL SORT(J,K,L,THET,D,AP,CP,COST,NUM)
      CALL LLS(COST,AP,5,A,YES,YE)
      DO 425 I = 1,5
      XCOS(I) = COST(I)
  425 YCOS(I) = AP(I)
      CALL TTPLOT(XCOS,YCOS,XT,YT,GT,5)
      WRITE(ITYP,36)
   36 FORMAT(//' THE FOLLOWING APPLY TO THE EXTRAPOLATION OF C'/
     1' USING ONLY THE 5 REFLECTIONS WITH HIGHEST L INDICES'/)
      CALL SORT(L,K,J,THET,D,AP,CP,COST,NUM)
      CALL LLS(COST,CP,5,C,YES,YE)
      DO 430 I = 1,5
      XCOS(I) = COST(I)
  430 YCOS(I) = CP(I)
      CALL TTPLOT(XCOS,YCOS,XT,YT,GT,5)
      GO TO 200
  900 WRITE(ITYP,40)
   40 FORMAT(//' OK THEN WE''LL CALL IT QUITS'/
     1' I DON''T KNOW WHY YOUR TIRED I DID ALL THE WORK!!')
      CALL EXIT
      END
      SUBROUTINE LLS(X,Y,NUM,YINCP,YES,YE)
      IMPLICIT REAL*8 (A-H)
      IMPLICIT REAL*8 (O-Z)
      INTEGER YES,YE
      DIMENSION X(50),  Y(50)
      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=DSQRT(SQRS)
C     VARS MEANS THE VARIANCE FOR THE SLOPE
      VARS=SQRS/(SUMXX-AN*XB*XB)
      SS=DSQRT(VARS)
C     SINCP MEANS THE STANDARD DEVIATIONOF THE Y-INTERCEPT
      SINCP=DSQRT(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',//)
      IF( YES.NE.YE ) GO TO 150
      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 SIMULT(A1,B1,C1,A2,B2,C2,X,Y)
      IMPLICIT REAL*8 (A-H)
      IMPLICIT REAL*8 (O-Z)
C
C SUBROUTINE TO SOLVE TWO EQUATIONS WITH TWO UNKNOWNS
C
C MEANING OF PARAMETERS:
CEQUATION1:   A1*X  +  B1*Y  =  C1
CEQUATION2:   A2*X  +  B2*Y  =  C2
C
C  USER SUPPLIES A,B, AND C VALUES ROUTINE CALCULATES X AND Y
C
      Y= (C1/A1 - C2/A2)/(B1/A1 - B2/A2)
      X = C2/A2 - Y*B2/A2
      RETURN
      END
      SUBROUTINE SORT(J1,J2,J3,A1,A2,A3,A4,A5,NUM)
      IMPLICIT REAL*8 (A-H)
      IMPLICIT REAL*8 (O-Z)
      DIMENSION J1(2),J2(2),J3(2),A1(2),A2(2),A3(2),A4(2),A5(2)
      NN = NUM - 1
      DO 100 I = 1,NN
      J = I + 1
      DO 100 II = J,NUM
      IF(J1(I).GE.J1(II)) GO TO 100
      JT=J1(I)
      J1(I) = J1(II)
      J1(II) = JT
      JT=J2(I)
      J2(I) = J2(II)
      J2(II) = JT
      JT=J3(I)
      J3(I) = J3(II)
      J3(II) = JT
      AT=A1(I)
      A1(I) = A1(II)
      A1(II) = AT
      AT=A2(I)
      A2(I) = A2(II)
      A2(II) = AT
      AT=A3(I)
      A3(I) = A3(II)
      A3(II) = AT
      AT=A4(I)
      A4(I) = A4(II)
      A4(II) = AT
      AT=A5(I)
      A5(I) = A5(II)
      A5(II) = AT
  100 CONTINUE
      RETURN
      END
      SUBROUTINE TTPLOT(X,Y,XTITLE,YTITLE,GRAPH,NPTS)
      DIMENSION PRINT(63),X(101),Y(101),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
      PRINT(63)=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
