PROGRAM REFPK C C ROBERT H. BLESSING C MEDICAL FOUNDATION OF BUFFALO C 73 HIGH STREET C BUFFALO, NEW YORK 14203, USA C TELEPHONE: (716) 856-9600 C ELECTRONIC MAIL: Blessing@MFB.Buffalo.Edu C C JULY 1993 C PARAMETER (NMAX=100) CHARACTER ATIME*8,ADATE*9,TITLE*80 COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2 COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2 COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2, & Q1(3,3),Q2(3,3) COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL COMMON /BLOCK7/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT COMMON /BLOCKC/ NCOND,ICOND(38) CALL INPUT CALL DATA CALL TFIT CALL UQFIT CALL OUTPUT STOP 'PROGRAM REFPK FINIS!' END C----------------------------------------------------------------------- SUBROUTINE ACCUM1 (N,X,Y,A,B) C C ACCUMULATE LEAST-SQUARES NORMAL MATRIX AND VECTOR ELEMENTS. C DIMENSION X(N),A(N,N),B(N) DOUBLE PRECISION A DO 1 I=1,N DO 2 J=I,N A(I,J)=A(I,J)+X(I)*X(J) A(J,I)=A(I,J) 2 CONTINUE B(I)=B(I)+X(I)*Y 1 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE ACCUM2 (N,X,Y1,Y2,A,B1,B2) C C ACCUMULATE LEAST-SQUARES NORMAL MATRIX AND VECTOR ELEMENTS. C DIMENSION X(N),A(N,N),B1(N),B2(N) DOUBLE PRECISION A DO 1 I=1,N DO 2 J=I,N A(I,J)=A(I,J)+X(I)*X(J) A(J,I)=A(I,J) 2 CONTINUE B1(I)=B1(I)+X(I)*Y1 B2(I)=B2(I)+X(I)*Y2 1 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE AGREE0 (IFILE,U,NOBS,NPAR,CHISQ1,CHISQ2,CHISQ3,RMSD0, & RMSD1,RMSDX) C C STATISTICS OF FIT FOR REFLECTION PEAK POSITIONS C DIMENSION U(3,3),H(3),Y(3),Z(3) PI=ACOS(-1.0) RAD=PI/180 DEG=180/PI SUMDX=0 SUMD0=0 SUMD1=0 CHISQ1=0 CHISQ2=0 CHISQ3=0 REWIND IFILE DO 5 N=1,NOBS READ (IFILE) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z, & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH Y1=Y(1) Y2=Y(2) Y3=Y(3) DO 1 J=1,3 Y(J)=0 DO 1 K=1,3 Y(J)=Y(J)+H(K)*U(K,J) 1 CONTINUE DO 2 J=1,3 Y(J)=Y(J)/DSTAR 2 CONTINUE IF (TWOTH.LT.0) THEN DO 3 J=1,3 Y(J)=-Y(J) 3 CONTINUE END IF CHISQ1=CHISQ1+(Y1-Y(1))**2 CHISQ2=CHISQ2+(Y2-Y(2))**2 CHISQ3=CHISQ3+(Y3-Y(3))**2 COSPH=COS(PHI*RAD) SINPH=SIN(PHI*RAD) COSCH=COS(CHI*RAD) SINCH=SIN(CHI*RAD) SINOM=Y(1)*COSPH-Y(2)*SINPH COSOM=(Y(1)*SINPH+Y(2)*COSPH)*COSCH-Y(3)*SINCH DELTA=ARCTAN(SINOM,COSOM)*DEG-OMEGA DELTA=A180(DELTA) IF (TWOTH.LT.0) DELTA=-DELTA SUMD0=SUMD0+DELTA**2 X0CALC=XM-DELTA SUMD1=SUMD1+(X0-X0CALC)**2 DX=(W1+X2-X1+W2)/(L2-L1+1) SUMDX=SUMDX+DX**2 5 CONTINUE NFREE=NOBS-NPAR CHISQ1=CHISQ1/NFREE CHISQ2=CHISQ2/NFREE CHISQ3=CHISQ3/NFREE RMSDX=SQRT(SUMDX/NOBS) RMSD0=SQRT(SUMD0/NOBS) RMSD1=SQRT(SUMD1/(NOBS-3*NPAR)) RETURN END C----------------------------------------------------------------------- SUBROUTINE AGREE1 (MODEL,IFILE,Q1,Q2,T1,T2,TANTHM,F,WL0, & NOBS,NPAR,CHISQ1,CHISQ2,RMSD1,RMSD2) CHARACTER MODEL*4 C C STATISTICS OF FIT FOR REFLECTION WIDTHS C DIMENSION Q1(3,3),Q2(3,3),H(3),Y(3),Z(3) IF (NPAR.LT.6) THEN U1=Q1(1,1) U2=Q2(1,1) END IF CHISQ1=0 CHISQ2=0 SUMD1=0 SUMD2=0 REWIND IFILE DO 5 N=1,NOBS READ (IFILE) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z, & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH IF (NPAR.GE.6) THEN U1=0 U2=0 DO 1 I=1,3 DO 1 J=1,3 U1=U1+Z(J)*Q1(I,J)*Z(I) U2=U2+Z(J)*Q2(I,J)*Z(I) 1 CONTINUE END IF IF (MODEL.EQ.'LRNZ') THEN W1CALC=SQRT(U1)+T1*TANTH W2CALC=SQRT(U2)+T2*TANTH ELSE IF (MODEL.EQ.'GAUS') THEN W1CALC=SQRT(U1+(T1*TANTH)**2) W2CALC=SQRT(U2+(T2*TANTH)**2) ELSE IF (MODEL.EQ.'PLMC') THEN W1CALC=SQRT(U1-2*F*T1*TANTH/TANTHM+T1*(TANTH/TANTHM)**2) W2CALC=SQRT(U2-2*F*T2*TANTH/TANTHM+T2*(TANTH/TANTHM)**2) END IF IF (MODEL.EQ.'LRNZ') THEN CHISQ1=CHISQ1+(W1-W1CALC)**2 CHISQ2=CHISQ2+(W2-W2CALC)**2 ELSE CHISQ1=CHISQ1+(W1**2-W1CALC**2)**2 CHISQ2=CHISQ2+(W2**2-W2CALC**2)**2 END IF SUMD1=SUMD1+(W1-W1CALC)**2 SUMD2=SUMD2+(W2-W2CALC)**2 5 CONTINUE NFREE=NOBS-NPAR RMSD1=SQRT(SUMD1/NFREE) RMSD2=SQRT(SUMD2/NFREE) CHISQ1=CHISQ1/NFREE CHISQ2=CHISQ2/NFREE RETURN END C----------------------------------------------------------------------- SUBROUTINE CENTROID (WINDOW,CUTOFF,XX,YY,VX,VY,I1,I2,X0,SIGX0) C C INTENSITY-WEIGHTED CENTROID OF PEAK ABOVE LINEAR BACKGROUND C DIMENSION XX(1),YY(1),VX(1),VY(1) PARAMETER (NMAX=100) DIMENSION X(NMAX),Y(NMAX),V(NMAX) EQUIVALENCE (V(1),X(1)) C C FIND THE SMALLEST X AND LARGEST X FOR WHICH THE LOCAL AVERAGE Y(X) C INCREASES SIGNIFICANTLY ABOVE THE LOCAL AVERAGE BACKGROUND. C N=0 DO 11 I=I1,I2,+1 N=N+1 Y(N)=YY(I) V(N)=VY(I) 11 CONTINUE CALL FINDL (N,Y,V,WINDOW,CUTOFF,L) IF (L.EQ.0) GO TO 9 L1=I1+L N=0 DO 12 I=I2,I1,-1 N=N+1 Y(N)=YY(I) V(N)=VY(I) 12 CONTINUE CALL FINDL (N,Y,V,WINDOW,CUTOFF,L) IF (L.EQ.0) GO TO 9 L2=I2-L IF (L1.GE.L2) GO TO 9 C C FIT A STRAIGHT LINE TO THE BACKGROUND NEAR THE PEAK LIMITS. C L=NINT(WINDOW*N) N=0 DO 1 I=L1-L,L2+L IF (I.LT.L1.OR.I.GT.L2) THEN N=N+1 X(N)=XX(I) Y(N)=YY(I) END IF 1 CONTINUE CALL FITLINE (N,X,Y,C0,C1,V0,V1,COV) C C FIND THE INTENSITY-WEIGHTED CENTROID BETWEEN THE PEAK LIMITS. C SUMY=0 SUMYX=0 SUMX2V=0 SUMXV=0 SUMV=0 DO 2 I=L1,L2 XI=XX(I) YI=YY(I)-C0-C1*XI VI=VY(I)+V0+V1*XI**2+2*COV*XI DX=0.50*(XX(I+1)-XX(I-1)) VD=0.25*(VX(I+1)+VX(I-1)) VI=YI**2*VD+DX**2*VI YI=YI*DX SUMY=SUMY+YI SUMYX=SUMYX+YI*XI SUMX2V=SUMX2V+XI**2*VI SUMXV=SUMXV+XI*VI SUMV=SUMV+VI 2 CONTINUE C C TEST FOR A SIGNIFICANTLY POSITIVE PEAK ABOVE BACKGROUND BETWEEN C THE TRUNCATED PEAK LIMITS. C IF (SUMY.LT.CUTOFF*SQRT(SUMV)) GO TO 9 X0=SUMYX/SUMY SIGX0=SQRT(ABS(SUMX2V-2*X0*SUMXV+X0**2*SUMV)/SUMY**2) RETURN 9 CONTINUE C C RETURN X0 = 90 DEGREES TO FLAG AN ERROR CONDITION. C X0=90 RETURN END C----------------------------------------------------------------------- SUBROUTINE DATA C C SELECT REFLECTION PROFILES FOR THE LEAST-SQUARES CALCULATIONS OF C PEAK POSITION AND PEAK WIDTH PARAMETERS. C PARAMETER (NMAX=100) CHARACTER ATIME*8,ADATE*9,TITLE*80 COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2 COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2 COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2, & Q1(3,3),Q2(3,3) COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL COMMON /BLOCKC/ NCOND,ICOND(38) DIMENSION H(3),ANGLES(4),Y(3),Z(3) PI=ACOS(-1.0) DEG=180/PI RAD=PI/180 C C INITIALIZE COUNTERS FOR REFLECTIONS THAT FAIL THE ACCEPTANCE TESTS. C N0=0 N1=0 N2=0 C C INITIALIZE MINIMUM AND MAXIMUM SERIAL NUMBERS AND EXPOSURE TIMES. C IIMIN=+1E6 IIMAX=-1E6 XTMIN=+1E5 XTMAX=-1E5 C C OPEN A SCRATCH FILE FOR LEAST-SQUARES DATA. C OPEN (UNIT=IO3,FILE='SCRATCH',STATUS='SCRATCH',FORM='UNFORMATTED') C C LOOP THROUGH ANALYSIS OF REFLECTION PROFILE DATA. C IEND=0 1 CALL READR (IO1,IEND,II,IH,IK,IL,ANGLES,NSTEP,XX,YY,VX,VY,XTIME) IF (IEND.NE.0) GO TO 9 C C SKIP REFLECTIONS THAT ARE NOT WITHIN THE SERIAL NUMBER OR EXPOSURE C TIME LIMITS. C IF (ABS(II).LT.II1.OR.ABS(II).GT.II2) GO TO 1 IF (XTIME.LT.XT1.OR.XTIME.GT.XT2) GO TO 1 C C TEST AGAINST CONDITIONS LIMITING POSSIBLE REFLECTIONS. C CALL HKLTEST (IH,IK,IL,NCOND,ICOND,IABSENT) IF (IABSENT.NE.0) GO TO 1 N0=N0+1 C C STORE MINIMUM AND MAXIMUM SERIAL NUMBERS AND EXPOSURE TIMES. C IIMIN=MIN(IIMIN,ABS(II)) IIMAX=MAX(IIMAX,ABS(II)) XTMIN=MIN(XTMIN,XTIME) XTMAX=MAX(XTMAX,XTIME) C C SMOOTH THE STATISTICAL FLUCTUATIONS IN THE SCAN PROFILE. C CALL SMOOTH (NSMOOTH,NSTEP,YY) CALL SMOOTH (NSMOOTH,NSTEP,VY) C C SET THE SCAN LIMITS (DEFAULT VALUES ARE SC1 = 0, SC2 = 1). C DX=XX(NSTEP)-XX(1) I1=NINDEX(XX(1)+SC1*DX) I2=NINDEX(XX(1)+SC2*DX) C C FIND THE INTENSITY-WEIGHTED CENTROID OF THE PEAK. C CALL CENTROID (WINDOW,CUTOFF,XX,YY,VX,VY,I1,I2,X0,SIGX0) IF (X0.EQ.90) THEN N1=N1+1 GO TO 1 END IF C C STORE THE INDICES IN A REAL ARRAY. C H(1)=IH H(2)=IK H(3)=IL C C CALCULATE SIN(THETA)/LAMBDA. C S=SINTHL(H,GINV) C C CALCULATE SPECTRAL THETA VALUES. C THE=ASIN(S*WLE)*DEG TH1=ASIN(S*WL1)*DEG TH0=ASIN(S*WL0)*DEG TH2=ASIN(S*WL2)*DEG C C CALCULATE THE POSITIONS OF THE THE ALPHA-ONE AND ALPHA-TWO PEAKS, C ASSUMING THE CENTROID TO BE THE ALPHA-BAR POSITION. C X1=X0+(TH1-TH0) X2=X0+(TH2-TH0) CALL XCHECK (X1,X0,X2,XX,YY,I1,I2) C C ALLOW FOR THE BETA FILTER ABSORPTION EDGE. C J1=I1 IF (WLE.NE.0) THEN XE=X0+(THE-TH0) I1=MAX(I1,NINDEX(XE)) END IF C C FIND THE PEAK LIMITS. C CALL LEHLAR (L) IF (L.EQ.0) THEN N2=N2+1 GO TO 1 END IF I1=J1 C C CALCULATE ANGULAR BASE-WIDTHS OF HALF-PEAKS BELOW THETA(ALPHA-ONE) C AND ABOVE THETA(ALPHA-TWO). (L1 AND L2 ARE THE PEAK LIMITS, I.E., C STEPS L1 AND L2 ARE IN THE PEAK.) C W1=X1-XX(L1) W2=XX(L2)-X2 C C TRANSFORM TO SETTING ANGLES AS DEFINED BY WALTER C. HAMILTON C (1974). INTERNATIONAL TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV, C PP. 273-284. C CALL GEOM (IDIFF,ANGLES,TWOTH,OMEGA,CHI,PHI) C C CALCULATE THE DISPLACEMENT OF THE PEAK CENTROID FROM THE MIDPOINT C OF THE SCAN RANGE. C XM=0.5*(XX(I1)+XX(I2)) DELTA=X0-XM C C GET UNIT DIRECTION VECTORS FOR THE DIFFRACTOMETER Y- AND Z-AXES. C CALL YZ (DELTA,TWOTH,OMEGA,CHI,PHI,Y,Z) C C CALCULATE DSTAR AND TAN(THETA). C DSTAR=2*S TANTH=TAN(ASIN(S*WL0)) C C WRITE DATA TO A FILE OF LEAST-SQUARES OBSERVATIONS. C WRITE (IO3) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z, & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH C C LOOP BACK FOR NEXT REFLECTION. C GO TO 1 C C END LOOP THROUGH REFLECTION DATA. C 9 ENDFILE IO3 CLOSE (UNIT=IO1) C C WRITE DATA SET LIMITS TO BGLP.DAT FILE. C WRITE (IO2,601) IIMIN,IIMAX,XTMIN-0.005,XTMAX+0.005 601 FORMAT (1X,2I10,2F10.3) C C PRINT SOME DATA SET STATISTICS. C WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,611) IIMIN,IIMAX,XTMIN-0.005,XTMAX+0.005 WRITE (ILP,612) N0,N1,N2,N0-(N1+N2) 610 FORMAT ('1'/'0PROGRAM REFPK. ',A,1X,A,'. ',A) 611 FORMAT ('0BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY & EXPOSURE TIMES (HOURS)'/'0 IIMIN IIMAX XTMI &N XTMAX'/' '2I10,10X,2F10.2) 612 FORMAT ('0STATISTICS ON REFLECTIONS PROCESSED:'/'0N0 = ',I5,' TOT &AL REFLECTIONS (EXCLUDING SYMMETRY FORBIDDEN REFLECTIONS)'/' N1 ', &'= ',I5,' REFLECTIONS DID NOT HAVE A SIGNIFICANT PEAK ABOVE LINEA &R BACKGROUND'/' N2 = ',I5,' ADDITIONAL REFLECTIONS DID NOT GIVE W &ELL-DEFINED PEAK LIMITS'/' N0 - (N1 + N2) = ',I5,' REFLECTIONS US &ED FOR LEAST-SQUARES CALCULATIONS OF PEAK POSITION AND WIDTH PARAM &ETERS') RETURN END C----------------------------------------------------------------------- SUBROUTINE DPRINT (II,IH,IK,IL,IPRINT) C C IPRINT = 1 FOR H00, 0K0, 00L, HH0, 0KK, H0H, AND HHH. C IPRINT = 0 OTHERWISE. C IF (II.LT.0) GO TO 10 J=ABS(IH) K=ABS(IK) L=ABS(IL) IF (J.EQ.0.AND.K.EQ.0) GO TO 11 IF (J.EQ.0.AND.L.EQ.0) GO TO 11 IF (K.EQ.0.AND.L.EQ.0) GO TO 11 IF (J.EQ.K.AND.L.EQ.0) GO TO 11 IF (J.EQ.L.AND.K.EQ.0) GO TO 11 IF (K.EQ.L.AND.J.EQ.0) GO TO 11 IF (J.EQ.K.AND.K.EQ.L) GO TO 11 10 IPRINT=0 RETURN 11 IPRINT=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE EIGEN (N,M,A,R) C C COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX C C A - ORIGINAL MATRIX, DESTROYED IN COMPUTATION. C RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF MATRIX A IN C DESCENDING ORDER. C R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, IN SAME C SEQUENCE AS EIGENVALUES) C N - ORDER OF MATRICES A AND R C - LOGICAL SIZE OF ARRAYS A AND R C M - PHYSICAL SIZE OF ARRAYS A AND R, M .GE. N C C E.G., FOR N = 3 SUBROUTINE EIGEN RETURNS C C (W1 0 0 ) C A = (0 W2 0 ) , WITH W1 .GE. W2 .GE. W3, AND C (0 0 W3) C C (X1 X2 X3) C R = (Y1 Y2 Y3) . C (Z1 Z2 Z3) C C THE DIAGONAL ELEMENTS OF A ARE THE EIGENVALUES, AND THE COLUMNS OF C R ARE THE CORRESPONDING EIGENVECTORS. C C DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED BY VON C NEUMANN FOR LARGE COMPUTERS AS FOUND IN 'MATHEMATICAL METHODS FOR C DIGITAL COMPUTERS', VOL. I, EDITED BY A. RALSTON AND H. S. WILF, C JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7. C DIMENSION A(M,M),R(M,M) C C GENERATE IDENTITY MATRIX AS FIRST APPROXIMATION TO THE MATRIX OF C EIGENVECTORS. C DO 10 I=1,N R(I,I)=1 10 CONTINUE DO 15 I=1,N-1 DO 15 J=I+1,N R(I,J)=0 R(J,I)=0 15 CONTINUE C C CALCULATE INITIAL AND FINAL OFF-DIAGONAL NORMS FOR SETTING AND C TESTING THRESHOLD. C Q=0 DO 20 I=1,N-1 DO 20 J=I+1,N Q=Q+A(I,J)**2 20 CONTINUE Q=SQRT(2*Q) E=1E-10 IF (Q.LT.E) GO TO 41 QF=Q*E/N C C LOWER THRESHOLD AT EACH ITERATION. C DO 40 ICYCLE=1,50 Q=Q/N C C LOOP THROUGH OFF-DIAGONAL PIVOT ELEMENTS. C DO 30 I=1,N-1 DO 30 J=I+1,N IF (ABS(A(I,J)).GT.Q) THEN C C CALCULATE ROTATION TO ANNHILATE OFF-DIAGONAL ELEMENT A(I,J). C U=-A(I,J) V=(A(I,I)-A(J,J))/2 W=SIGN(1.0,V)*U/SQRT(U**2+V**2) C C SINE AND COSINE OF ROTATION ANGLE C S=W/SQRT(2*(1+SQRT(1-W**2))) C=SQRT(1-S**2) C2=C**2 S2=S**2 CS=C*S C C TRANSFORM AND REPLACE THE ELEMENTS IN THE ITH AND JTH COLUMNS OF C THE A MATRIX AND OF THE R MATRIX. C DO 35 K=1,N C C SKIP ELEMENTS OF THE PIVOT SET A(I,J), (A(J,I)), A(I,I), A(J,J). C IF (K.NE.I.AND.K.NE.J) THEN BKI=A(K,I)*C-A(K,J)*S BKJ=A(K,I)*S+A(K,J)*C A(K,I)=BKI A(K,J)=BKJ A(I,K)=A(K,I) A(J,K)=A(K,J) END IF BKI=R(K,I)*C-R(K,J)*S BKJ=R(K,I)*S+R(K,J)*C R(K,I)=BKI R(K,J)=BKJ 35 CONTINUE C C COMPUTE THE NEW OFF-DIAGONAL PIVOT ELEMENT AND ITS ASSOCIATED C DIAGONAL ELEMENTS. C BII=A(I,I)*C2+A(J,J)*S2-2*A(I,J)*CS BJJ=A(I,I)*S2+A(J,J)*C2+2*A(I,J)*CS BIJ=(A(I,I)-A(J,J))*CS+A(I,J)*(C2-S2) A(I,I)=BII A(J,J)=BJJ A(I,J)=BIJ A(J,I)=A(I,J) END IF 30 CONTINUE IF (Q.LT.QF) GO TO 41 40 CONTINUE 41 CONTINUE C C SORT ON DECREASING EIGENVALUE. C DO 45 I=1,N-1 DO 45 J=I+1,N IF (A(I,I).LT.A(J,J)) THEN T=A(I,I) A(I,I)=A(J,J) A(J,J)=T DO 50 K=1,N T=R(K,I) R(K,I)=R(K,J) R(K,J)=T 50 CONTINUE END IF 45 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE EXTREM (IXMIN,IXMAX,IX) PARAMETER (N1=4,N2=16,N3=50) DIMENSION IXMIN(N1,N2,N3),IXMAX(N1,N2,N3),ITYPE(N1),IX(N2) DATA ITYPE /7, 9, 15, 16/ C C TABULATE IN RANKED ORDER THE N3 SMALLEST AND N3 LARGEST C VALUES OF THE VARIABLE OF TYPE I1. C C I1 TYPE OF VARIABLE C I2 VALUE OF VARIABLE C I3 RANK IN LIST C DO 10 I1=1,N1 I2=ITYPE(I1) DO 11 I3=1,N3 IF (IX(I2).GE.IXMIN(I1,I2,I3)) GO TO 11 DO 12 J3=N3,I3+1,-1 DO 12 I2=1,N2 12 IXMIN(I1,I2,J3)=IXMIN(I1,I2,J3-1) DO 13 I2=1,N2 13 IXMIN(I1,I2,I3)=IX(I2) GO TO 20 11 CONTINUE 20 CONTINUE I2=ITYPE(I1) DO 21 I3=1,N3 IF (IX(I2).LE.IXMAX(I1,I2,I3)) GO TO 21 DO 22 J3=N3,I3+1,-1 DO 22 I2=1,N2 22 IXMAX(I1,I2,J3)=IXMAX(I1,I2,J3-1) DO 23 I2=1,N2 23 IXMAX(I1,I2,I3)=IX(I2) GO TO 10 21 CONTINUE 10 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE FINDL (N,Y,V,WINDOW,CUTOFF,L) C C FIND THE X AT WHICH THE LOCAL AVERAGE Y(X) FIRST INCREASES C SIGNIFICANTLY ABOVE THE LOCAL AVERAGE BACKGROUND. C DIMENSION Y(N),V(N) L=NINT(WINDOW*N) DO 1 I=L+1,N-L YB=0 VB=0 YI=0 VI=0 DO 2 J=0,L YB=YB+Y(I-J) VB=VB+V(I-J) YI=YI+Y(I+J) VI=VI+V(I+J) 2 CONTINUE IF (YI-YB.GT.CUTOFF*SQRT(VI+VB)) THEN L=I RETURN END IF 1 CONTINUE L=0 RETURN END C----------------------------------------------------------------------- SUBROUTINE FITLINE (N,X,Y,C0,C1,V0,V1,COV) C C LEAST-SQUARES STRAIGHT LINE BACKGROUND C C GIVE EQUAL WEIGHT TO ALL BACKGROUND POINTS BECAUSE: (1) IT IS C EXPECTED THAT THE BACKGROUND SHOULD BE APPROXIMATELY CONSTANT, AND C THUS CONSTANT WEIGHTS SHOULD BE APPROPRIATE; AND (2) WEIGHTS C PROPORTIONAL TO THE RECIPROCALS OF THE VARIANCES OF THE C EXPERIMENTAL MEASUREMENTS, ALTHOUGH APPROPRIATE FOR NORMALLY C DISTRIBUTED DATA, ARE NOT APPROPRIATE FOR POISSON DISTRIBUTED C DATA. WEIGHT = 1/VARIANCE WOULD BIAS THE FIT TOWARD THE C BACKGROUND POINTS WITH THE LOWEST COUNT RATES. (SEE P. F. PRICE C (1979). ACTA CRYST. A35, 57-60.) C DIMENSION X(N),Y(N) SUMX=0 SUMY=0 SUMX2=0 SUMXY=0 DO 1 I=1,N XI=X(I) YI=Y(I) SUMX=SUMX+XI SUMY=SUMY+YI SUMX2=SUMX2+XI*XI SUMXY=SUMXY+XI*YI 1 CONTINUE DET=N*SUMX2-SUMX*SUMX C0=(SUMX2*SUMY-SUMX*SUMXY)/DET C1=(-SUMX*SUMY+N*SUMXY)/DET V0=SUMX2/DET V1=N/DET COV=-SUMX/DET CHISQ=0 DO 2 I=1,N CHISQ=CHISQ+(Y(I)-C0-C1*X(I))**2 2 CONTINUE CHISQ=CHISQ/(N-2) V0=CHISQ*V0 V1=CHISQ*V1 COV=CHISQ*COV RETURN END C----------------------------------------------------------------------- SUBROUTINE GEOM (IDIFF,ANGLES,TWOTH,OMEGA,CHI,PHI) C C DIFFRACTOMETER GEOMETRIES C DIMENSION ANGLES(4) REAL KAPPA C C KAPPA AXIS INCLINATION ANGLE, ALPHA = PI/3.6 = 50 DEGREES C DATA SINA, COSA /0.766044, 0.642788/ DATA DEG, RAD /57.2957795, 0.017453293/ IF (IDIFF.EQ.1) THEN C C HAMILTON'S DIFFRACTOMETER AXES C C WALTER C. HAMILTON (1974). INTERNATIONAL TABLES FOR X-RAY C CRYSTALLOGRAPHY, VOL. IV, PP. 273-284. C TWOTH = ANGLES(1) OMEGA = ANGLES(2) CHI = ANGLES(3) PHI = ANGLES(4) ELSE IF (IDIFF.EQ.2) THEN C C BUSING'S AND LEVY'S DIFFRACTOMETER AXES C C W. R. BUSING AND H. A. LEVY (1967). ACTA CRYST. 22, 457-464. C TWOTH = ANGLES(1) OMEGA = ANGLES(2) CHI = ANGLES(3) PHI = ANGLES(4) C C TRANSFORM TO SETTING ANGLES AS DEFINED BY HAMILTON. C OMEGA = -OMEGA CHI = -CHI PHI = -PHI ELSE IF (IDIFF.EQ.3) THEN C C SIEMENS (NEE NICOLET, NEE SYNTEX) P3 DIFFRACTOMETER C TWOTH = ANGLES(1) OMEGA = ANGLES(2) PHI = ANGLES(3) CHI = ANGLES(4) C C TRANSFORM TO SETTING ANGLES AS DEFINED BY HAMILTON. C THETA = TWOTH/2 OMEGA = OMEGA - THETA OMEGA = -OMEGA ELSE IF (IDIFF.EQ.4) THEN C C ENRAF-NONIUS CAD4 DIFFRACTOMETER KAPPA ANGLES C THETA = ANGLES(1) PHI = ANGLES(2) OMEGA = ANGLES(3) KAPPA = ANGLES(4) C C TRANSFORM FROM KAPPA TO EULERIAN SETTING ANGLES. C KAPPA = KAPPA*RAD SINX = SINA*SIN(KAPPA/2) COSX = SQRT(COSA**2 + SINA**2*COS(KAPPA/2)**2) CHI = 2*ARCTAN(SINX,COSX) C C WITH ALPHA = 50 DEGREES, KAPPA IS MECHANICALLY LIMITED TO THE RANGE 0 C TO 180 DEGREES. THIS LIMITS CHI TO THE RANGE -100 TO +100 DEGREES. C SINX = COSA*SIN(KAPPA/2)/COS(CHI/2) COSX = COS(KAPPA/2)/COS(CHI/2) DELTA = ARCTAN(SINX,COSX)*DEG OMEGA = OMEGA + DELTA PHI = PHI + DELTA CHI = CHI*DEG C C TRANSFORM TO SETTING ANGLES AS DEFINED BY HAMILTON. C TWOTH = 2*THETA OMEGA = OMEGA - THETA OMEGA = -OMEGA CHI = -CHI PHI = -PHI END IF C C RETURN ANGLES -180 .LE. A .LE. +180. C TWOTH = A180(TWOTH) OMEGA = A180(OMEGA) PHI = A180(PHI) CHI = A180(CHI) RETURN END C----------------------------------------------------------------------- FUNCTION ARCTAN(SINX,COSX) IF (SINX.LT.-1) SINX=-1 IF (SINX.GT.+1) SINX=+1 IF (COSX.LT.-1) COSX=-1 IF (COSX.GT.+1) COSX=+1 ARCTAN=ATAN2(SINX,COSX) RETURN END C----------------------------------------------------------------------- FUNCTION A180(A) A=MOD(A,360.0) IF (A.LT.-180) A=A+360 IF (A.GT.+180) A=A-360 A180=A RETURN END C----------------------------------------------------------------------- FUNCTION A360(A) A=MOD(A,360.0) IF (A.LT.0) A=A+360 A360=A RETURN END C----------------------------------------------------------------------- SUBROUTINE HKLCOND (H,I) C C CONDITIONS LIMITING POSSIBLE REFLECTIONS C CHARACTER*18 H,A DIMENSION A(38) DATA A &/'HKL,H+K=2N ','HKL,K+L=2N ','HKL,L+H=2N ', & 'HKL,H+K,K+L,L+H=2N','HKL,H+K+L=2N ','HKL,-H+K+L=3N ', & 'HKL,H-K+L=3N ','HKL,H-K=3N ','HK0,H=2N ', & 'HK0,K=2N ','HK0,H+K=2N ','HK0,H+K=4N ', & '0KL,K=2N ','0KL,L=2N ','0KL,K+L=2N ', & '0KL,K+L=4N ','H0L,L=2N ','H0L,H=2N ', & 'H0L,L+H=2N ','H0L,L+H=4N ','HH(-2H)L,L=2N ', & 'H(-H)0L,L=2N ','HHL,L=2N(R-AXES) ','HHL,L=2N ', & 'HKH,K=2N ','HKK,H=2N ','HHL,2H+L=4N ', & 'HKH,2H+K=4N ','HKK,2K+H=4N ','H00,H=2N ', & 'H00,H=4N ','0K0,K=2N ','0K0,K=4N ', & '00L,L=2N ','00L,L=4N ','000L,L=2N ', & '000L,L=3N ','000L,L=6N '/ IF (I.EQ.0) THEN DO 1 I=1,38 IF (H.EQ.A(I)) RETURN 1 CONTINUE STOP 'INPUT ERROR. HKL CONDITION IS NOT IN THE TABLE.' ELSE H=A(I) RETURN END IF END C----------------------------------------------------------------------- SUBROUTINE HKLTEST (H,K,L,NCOND,ICOND,IABSENT) C C CONDITIONS LIMITING POSSIBLE REFLECTIONS C C ( 1) HKL,H+K=2N C ( 2) HKL,K+L=2N C ( 3) HKL,L+H=2N C ( 4) HKL,H+K,K+L=2N C ( 5) HKL,H+K+L=2N C ( 6) HKL,-H+K+L=3N C ( 7) HKL,H-K+L=3N C ( 8) HKL,H-K=3N C ( 9) HK0,H=2N C (10) HK0,K=2N C (11) HK0,H+K=2N C (12) HK0,H+K=4N C (13) 0KL,K=2N C (14) 0KL,L=2N C (15) 0KL,K+L=2N C (16) 0KL,K+L=4N C (17) H0L,L=2N C (18) H0L,H=2N C (19) H0L,L+H=2N C (20) H0L,L+H=4N C (21) HH(-2H)L,L=2N C (22) H(-H)0L,L=2N C (23) HHL,L=2N (RHOMBOHEDRAL AXES) C (24) HHL,L=2N C (25) HKH,K=2N C (26) HKK,H=2N C (27) HHL,2H+L=4N C (28) HKH,2H+K=4N C (29) HKK,2K+H=4N C (30) H00,H=2N C (31) H00,H=4N C (32) 0K0,K=2N C (33) 0K0,K=4N C (34) 00L,L=2N C (35) 00L,L=4N C (36) 000L,L=2N C (37) 000L,L=3N C (38) 000L,L=6N C C SET IABSENT = 0 FOR AN ALLOWED REFLECTION, C = 1 FOR A FORBIDDEN REFLECTION. C INTEGER H DIMENSION ICOND(38) IF (NCOND.EQ.0) THEN IABSENT=0 RETURN END IF DO 90 I=1,NCOND GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22, & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38) ICOND(I) 1 IF (MOD(H+K,2).NE.0) GO TO 91 GO TO 90 2 IF (MOD(K+L,2).NE.0) GO TO 91 GO TO 90 3 IF (MOD(L+H,2).NE.0) GO TO 91 GO TO 90 4 IF (MOD(H+K,2).NE.0.OR.MOD(K+L,2).NE.0) GO TO 91 GO TO 90 5 IF (MOD(H+K+L,2).NE.0) GO TO 91 GO TO 90 6 IF (MOD(-H+K+L,3).NE.0) GO TO 91 GO TO 90 7 IF (MOD(H-K+L,3).NE.0) GO TO 91 GO TO 90 8 IF (MOD(H-K,3).NE.0) GO TO 91 GO TO 90 9 IF (L.EQ.0.AND.MOD(H,2).NE.0) GO TO 91 GO TO 90 10 IF (L.EQ.0.AND.MOD(K,2).NE.0) GO TO 91 GO TO 90 11 IF (L.EQ.0.AND.MOD(H+K,2).NE.0) GO TO 91 GO TO 90 12 IF (L.EQ.0.AND.MOD(H+K,4).NE.0) GO TO 91 GO TO 90 13 IF (H.EQ.0.AND.MOD(K,2).NE.0) GO TO 91 GO TO 90 14 IF (H.EQ.0.AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 15 IF (H.EQ.0.AND.MOD(K+L,2).NE.0) GO TO 91 GO TO 90 16 IF (H.EQ.0.AND.MOD(K+L,4).NE.0) GO TO 91 GO TO 90 17 IF (K.EQ.0.AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 18 IF (K.EQ.0.AND.MOD(H,2).NE.0) GO TO 91 GO TO 90 19 IF (K.EQ.0.AND.MOD(L+H,2).NE.0) GO TO 91 GO TO 90 20 IF (K.EQ.0.AND.MOD(L+H,4).NE.0) GO TO 91 GO TO 90 21 IF (K.EQ.H.AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 22 IF (K.EQ.-H.AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 23 IF (K.EQ.H.AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 24 IF (ABS(K).EQ.ABS(H).AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 25 IF (ABS(L).EQ.ABS(H).AND.MOD(K,2).NE.0) GO TO 91 GO TO 90 26 IF (ABS(L).EQ.ABS(K).AND.MOD(H,2).NE.0) GO TO 91 GO TO 90 27 IF (ABS(K).EQ.ABS(H).AND.MOD(2*H+L,4).NE.0) GO TO 91 GO TO 90 28 IF (ABS(L).EQ.ABS(H).AND.MOD(2*H+K,4).NE.0) GO TO 91 GO TO 90 29 IF (ABS(L).EQ.ABS(K).AND.MOD(2*K+H,4).NE.0) GO TO 91 GO TO 90 30 IF (K.EQ.0.AND.L.EQ.0.AND.MOD(H,2).NE.0) GO TO 91 GO TO 90 31 IF (K.EQ.0.AND.L.EQ.0.AND.MOD(H,4).NE.0) GO TO 91 GO TO 90 32 IF (H.EQ.0.AND.L.EQ.0.AND.MOD(K,2).NE.0) GO TO 91 GO TO 90 33 IF (H.EQ.0.AND.L.EQ.0.AND.MOD(K,4).NE.0) GO TO 91 GO TO 90 34 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 35 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,4).NE.0) GO TO 91 GO TO 90 36 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,2).NE.0) GO TO 91 GO TO 90 37 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,3).NE.0) GO TO 91 GO TO 90 38 IF (H.EQ.0.AND.K.EQ.0.AND.MOD(L,6).NE.0) GO TO 91 90 CONTINUE IABSENT=0 RETURN 91 IABSENT=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE INPUT C C READ CONTROL DATA, AND SET DEFAULT VALUES. C PARAMETER (NMAX=100) CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE1*80,FILE2*80,HCOND*18 CHARACTER DIFF*4,SOURCE*2,TARGET*2,BETFIL*2,FILTER*2,MONOCR*2 COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2 COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2 COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2, & Q1(3,3),Q2(3,3) COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL COMMON /BLOCK7/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT COMMON /BLOCKC/ NCOND,ICOND(38) DIMENSION CELL(6) DIMENSION TARGET(3),WLA1(3),WLA2(3),DWLA1(3),DWLA2(3) DIMENSION BETFIL(5),ABSEDG(5) DIMENSION THMC(2,3),RHOMC(4) DIMENSION ER(2,4),CC(4,3,4) PI=ACOS(-1.0) C C K-ALPHA X-RAY WAVELENGTHS C DATA TARGET / 'CU', 'MO', 'AG'/ DATA WLA1 /1.54056, 0.70930, 0.55941/ DATA WLA2 /1.54439, 0.71359, 0.56380/ C C K-ALPHA SPECTRAL LINE-WIDTHS (FULL-WIDTHS AT HALF-HEIGHT) C DATA DWLA1 /0.00058, 0.00029, 0.00028/ DATA DWLA2 /0.00077, 0.00032, 0.00029/ C C K-BETA FILTER ABSORPTION EDGE WAVELENGTHS C DATA BETFIL / 'NI', 'NB', 'ZR', 'PD', 'RH'/ DATA ABSEDG /1.48807, 0.65298, 0.68883, 0.50920, 0.53395/ C C MONOCHROMATOR BRAGG ANGLES C C GRAPHITE (0 0 0 2) REFLECTION C QUARTZ (1 0-1 1) REFLECTION C C GRAPHITE, P 6(3)/M M C, Z = 4, A = 2.455(2), C = 6.69 A, C D(0 0 0 2) = 3.345 A C C QUARTZ, P 3(1,2) 2 1, Z = 3, A = 4.910(10), C = 5.394(10) A, C D(1 0-1 1) = 3.339 A C C GRAPHITE QUARTZ C C CU TH11 TH21 C MO TH12 TH22 C AG TH13 TH23 C DATA THMC /13.3, 13.3, & 6.1, 6.1, & 4.8, 4.8/ C C MONOCHROMATOR-DIFFRACTOMETER GEOMETRY ANGLE C C PARALLEL RHO = 0 OR C PERPENDICULAR RHO = 90 DEGREES C C INT.TAB. BUS.LEV. P3 CAD4 C DATA RHOMC / 0, 0, 0, 90/ C C ERROR ESTIMATES FOR DIFFRACTOMETER SCAN ANGLE (DEGREES THETA) AND C SCAN-ANGLE DRIVE-PULSE TIMING (MICORSECONDS) C C SIGTHE SIGTIM C C INT.TAB. ER11 ER21 C BUS.LEV. ER12 ER22 C P3 ER13 ER23 C CAD4 ER14 ER24 C DATA ER / 0.001, 1.0, & 0.001, 1.0, & 0.005, 10.0, & 0.005, 10.0/ C C COINCIDENCE CORRECTION DATA (DEAD TIME (MICROSECONDS) AND C ATTENUATOR FACTOR) C C TAU SIGTAU ATT SIGATT C C INT.TAB. CU CC111 CC211 CC311 CC411 C MO CC121 CC221 CC321 CC421 C AG CC131 CC231 CC331 CC431 C BUS.LEV. CU CC112 CC212 CC312 CC412 C MO CC122 CC222 CC322 CC422 C AG CC132 CC232 CC332 CC432 C P3 CU CC113 CC213 CC313 CC413 C MO CC123 CC223 CC323 CC423 C AG CC133 CC233 CC333 CC433 C CAD4 CU CC114 CC214 CC314 CC414 C MO CC124 CC224 CC324 CC424 C AG CC134 CC234 CC334 CC434 C DATA CC / 0.0, 0.0, 1.0, 0.0, & 0.0, 0.0, 1.0, 0.0, & 0.0, 0.0, 1.0, 0.0, & 0.0, 0.0, 1.0, 0.0, & 0.0, 0.0, 1.0, 0.0, & 0.0, 0.0, 1.0, 0.0, & 2.5, 0.2, 31.4, 0.5, & 2.48, 0.12, 27.01, 0.41, & 0.0, 0.0, 1.0, 0.0, & 0.0, 0.0, 25.6, 0.3, & 0.0, 0.0, 21.5, 0.3, & 0.0, 0.0, 1.0, 0.0/ C C DUMMY CONTROL VARIABLES C FILTER=' ' MONOCR=' ' THM=0 RHOM=0 FRACTD=0 SIGFD=0 F=0 WL1=0 WL2=0 DWL1=0 DWL2=0 WLE=0 SIGTHE=0 SIGTIM=0 TAU=0 SIGTAU=0 ATT=0 SIGATT=0 II1=0 II2=0 XT1=0 XT2=0 SC1=0 SC2=0 WINDOW=0 CUTOFF=0 FLL=0 WT2=0 NSMOOTH=0 JMODEL=0 C C START. C CALL TIME (ATIME) CALL DATE (ADATE) C C I/O DEVICES. LOGICAL UNIT NUMBERS C IO1=1 IO2=2 IO3=3 ILP=6 C C READ CONTROL DATA. C OPEN (UNIT=IO1,FILE='refpk.dat',STATUS='OLD') READ (IO1,500) TITLE READ (IO1,500) FILE1 READ (IO1,500) FILE2 READ (IO1,501) CELL READ (IO1,505) DIFF READ (IO1,502) SOURCE READ (IO1,502,END=9) FILTER READ (IO1,502,END=9) MONOCR READ (IO1,501,END=9) THM,RHOM,FRACTD,SIGFD READ (IO1,501,END=9) WL1,WL2,DWL1,DWL2,WLE READ (IO1,501,END=9) SIGTHE,SIGTIM READ (IO1,501,END=9) TAU,SIGTAU,ATT,SIGATT READ (IO1,503,END=9) II1,II2,XT1,XT2 READ (IO1,501,END=9) SC1,SC2 READ (IO1,501,END=9) WINDOW READ (IO1,501,END=9) CUTOFF READ (IO1,501,END=9) FLL READ (IO1,501,END=9) WT2 READ (IO1,504,END=9) NSMOOTH READ (IO1,504,END=9) JMODEL 500 FORMAT (A) 501 FORMAT (8F10.0) 505 FORMAT (A4) 502 FORMAT (A2) 503 FORMAT (2I10,2F10.0) 504 FORMAT (8I10) 9 CLOSE (UNIT=IO1) C C RECIPROCAL LATTICE METRIC TENSOR C CALL METRIC (CELL,GINV) C C DIFFRACTOMETER TYPE C IDIFF=0 IF (DIFF.EQ.'H ') IDIFF=1 IF (DIFF.EQ.'BL ') IDIFF=2 IF (DIFF.EQ.'P3 ') IDIFF=3 IF (DIFF.EQ.'CAD4') IDIFF=4 IF (IDIFF.EQ.0) STOP 'UNKNOWN DIFFRACTOMETER TYPE' C C WAVELENGTHS AND SPECTRAL WIDTHS C IXRAY=0 INEUTRON=0 DO 10 I=1,3 IF (SOURCE.EQ.TARGET(I)) IXRAY=I 10 CONTINUE IF (SOURCE.EQ.'N '.OR.SOURCE.EQ.'NE') INEUTRON=1 IF (IXRAY.NE.0) THEN I=IXRAY IF (WL1.EQ.0) WL1=WLA1(I) IF (WL2.EQ.0) WL2=WLA2(I) IF (DWL1.EQ.0) DWL1=DWLA1(I) IF (DWL2.EQ.0) DWL2=DWLA2(I) END IF IF (WT2.EQ.0) WT2=0.5 WL0=(WL1+WL2*WT2)/(1+WT2) C C BETA FILTER ABSORPTION EDGE C I=0 DO 20 J=1,5 IF (FILTER.EQ.BETFIL(J)) I=J 20 CONTINUE IF (I.NE.0.AND.WLE.EQ.0) WLE=ABSEDG(I) C C MONOCHROMATOR VARIABLES C IF (MONOCR.NE.' '.OR.THM.NE.0) THEN C C BRAGG ANGLE C IF (MONOCR.EQ.'GR') ICRYST=1 IF (MONOCR.EQ.'QU') ICRYST=2 IF (THM.EQ.0) THM=THMC(ICRYST,IXRAY) TANTHM=TAN(THM*PI/180) C C RHO ANGLE (PARALLEL OR PERPENDICULAR GEOMETRY) C IF (RHOM.EQ.0) RHOM=RHOMC(IDIFF) C C DYNAMIC DIFFRACTION (PERFECT CRYSTAL) FRACTION C IF (FRACTD.NE.0) SIGFD=0.05*FRACTD END IF C C ESTIMATED ERRORS FOR DIFFRACTOMETER SCAN ANGLE AND DRIVE PULSE C TIMING C J=IDIFF IF (SIGTHE.EQ.0) SIGTHE=ER(1,J) IF (SIGTIM.EQ.0) SIGTIM=ER(2,J) C C CONSTANTS FOR COINCIDENCE CORRECTIONS C J=IXRAY K=IDIFF IF (TAU.EQ.0) TAU=CC(1,J,K) IF (SIGTAU.EQ.0) SIGTAU=CC(2,J,K) IF (ATT.EQ.0) ATT=CC(3,J,K) IF (SIGATT.EQ.0) SIGATT=CC(4,J,K) C C BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY C EXPOSURE TIMES C IF (II1.EQ.0) II1=-1E6 IF (II2.EQ.0) II2=+1E6 IF (XT1.EQ.0) XT1=-1E5 IF (XT2.EQ.0) XT2=+1E5 C C INITIAL SCAN LIMITS (DECIMAL FRACTIONS OF SCAN WIDTH) C IF (SC1.LT.0) SC1=0 IF (SC2.LE.0.OR.SC2.GT.1) SC2=1 C C WINDOW WIDTH FOR INITIAL SCAN FOR PEAK LIMITS C IF (WINDOW.EQ.0) WINDOW=1/12.0 C C E.S.D. MULTIPLIER FOR STATISTICAL SIGNIFICANCE TESTS C IF (CUTOFF.EQ.0) CUTOFF=3 C C NUMBER OF DATA-SMOOTHING PASSES FOR PEAK LOCATION C IF (NSMOOTH.EQ.0) NSMOOTH=2 IF (NSMOOTH.LT.0) NSMOOTH=0 C C MULTIPLIER FOR LEHMANN-LARSEN BIAS CORRECTION C IF (FLL.EQ.0) FLL=2 C C CALCULATE SPECTRAL DISPERSION COEFFICIENTS FOR ALPHA-ONE AND C ALPHA-TWO LINES. C C DELTA(THETA) = (180/PI)*(DELTA(LAMBDA)/LAMBDA)*TAN(THETA) C C LORENTZIAN LINE SHAPE C C Y/YMAX = 1/(1 + (X - X0)**2/(W/2)**2) C C AT (X - X0) = 5*(W/2), C WHERE W/2 = HALF-WIDTH AT HALF-HEIGHT, C (X - X0) = 2.5*W C Y/YMAX = 1/(1 + 5**2) = 1/26 = 0.03846 C C GAUSSIAN LINE SHAPE C C Y/YMAX = EXP(-(X - X0)**2/(2*SIGMA**2)) C C AT Y/YMAX = 1/2, C 1/2 = EXP(-(W/2)**2/(2*SIGMA**2)) C LN(1/2) = -(W/2)**2/(2*SIGMA**2) C (W/2)/SIGMA = SQRT(-2*LN(1/2)) = SQRT(2*LN(2)) = 1.1774 C C W/2 = 1.1774*SIGMA C SIGMA = 0.8493*(W/2) C C AT Y/YMAX = 1/26, C (X - X0) = SIGMA*SQRT(-2*LN(1/26)) = 2.5527*SIGMA C (X - X0) = (W/2)*SQRT(LN(1/26)/LN(1/2)) = 2.1680*(W/2) C (X - X0) = 1.084*W C CL=2.500*180/PI CG=1.084*180/PI T1=DWL1/WL1 T2=DWL2/WL2 T1L=CL*T1 T2L=CL*T2 T1G=CG*T1 T2G=CG*T2 C C THE CHARACTERISTIC EMISSION LINE SHAPES SHOULD BE APPROXIMATELY C LORENTZIAN. C T1=T1L T2=T2L SIGT1=0 SIGT2=0 C C SET RATIO OF DISPERSION COEFFICIENTS FOR PARALLEL-GEOMETRY C MONOCHROMATOR. C C THE WIDTH OF THE TWICE-REFLECTED BEAM IS GIVEN BY C C W = SQRT(A - 2*B*TANTH/TANTHM + C*(TANTH/TANTHM)**2), C C WHERE C C B = ALPHA**2 + 2*ETA**2, C C = ALPHA**2 + ETA**2, C C ALPHA IS THE DIVERGENCE OF THE BEAM TO THE MONOCHROMATOR CRYSTAL, C AND ETA IS THE MOSAIC SPREAD OF THE MONOCHROMATOR CRYSTAL. C C TH AND THM HAVE THE SAME SIGN FOR PARALLEL GEOMETRY, OPPOSITE C SIGNS FOR FOR ANTI-PARALLEL GEOMETRY. C C W = SQRT(A - 2*F*C*TANTH/TANTHM + C*(TANTH/TANTHM)**2), C C WHERE F = B/C = (ALPHA**2 + 2*ETA**2)/(ALPHA**2 + ETA**2), AND C 1 .LE. F .LE. 2. C C FOR FIXED- OR ROTATING-ANODE DIVERGENT X-RAY SOURCES, EXPECT C EFFECTIVELY ALPHA = 2*ETA AND THEREFORE F = 6/5. C IF (RHOM.EQ.0.AND.F.EQ.0) F=1.2 C C WRITE CONTROL DATA TO BGLP.DAT FILE. C OPEN (UNIT=IO2,FILE='bglp.dat',STATUS='NEW') WRITE (IO2,600) ATIME WRITE (IO2,600) ADATE WRITE (IO2,600) TITLE WRITE (IO2,600) FILE1 WRITE (IO2,600) FILE2 WRITE (IO2,601) (CELL(I),I=1,6) WRITE (IO2,602) IDIFF,INEUTRON WRITE (IO2,601) THM,RHOM,FRACTD,SIGFD,F WRITE (IO2,601) WLE,WL1,WL0,WL2 WRITE (IO2,601) SIGTHE,SIGTIM WRITE (IO2,601) TAU,SIGTAU,ATT,SIGATT WRITE (IO2,601) WINDOW WRITE (IO2,601) CUTOFF WRITE (IO2,602) NSMOOTH WRITE (IO2,601) SC1,SC2 600 FORMAT (1X,A) 601 FORMAT (1X,6F10.5) 602 FORMAT (1X,2I10) C C PRINT CONTROL DATA. C OPEN (UNIT=ILP,FILE='refpk.lp',STATUS='NEW',CARRIAGECONTROL= & 'FORTRAN') WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,611) FILE1 WRITE (ILP,612) FILE2 WRITE (ILP,613) CELL IF (IDIFF.EQ.1) WRITE (ILP,614) IF (IDIFF.EQ.2) WRITE (ILP,619) IF (IDIFF.EQ.3) WRITE (ILP,615) IF (IDIFF.EQ.4) WRITE (ILP,616) WRITE (ILP,629) SOURCE IF (FILTER.NE.' ') WRITE (ILP,630) FILTER IF (THM.NE.0) WRITE (ILP,618) MONOCR,THM,RHOM,FRACTD,SIGFD WRITE (ILP,631) WLE,WL1,WL0,WL2,DWL1,DWL2 WRITE (ILP,632) T1,T2 IF (WT2.NE.0.5) WRITE (ILP,622) WT2 WRITE (ILP,617) SIGTHE,SIGTIM WRITE (ILP,641) TAU,SIGTAU,ATT,SIGATT IF (II1.GT.-1E6.OR.II2.LT.+1E6.OR.XT1.GT.-1E5.OR.XT2.LT.+1E5) & WRITE (ILP,640) II1,II2,XT1,XT2 IF (SC1.GT.0.OR.SC2.LT.1) WRITE (ILP,642) SC1,SC2 WRITE (ILP,648) WINDOW WRITE (ILP,646) CUTOFF IF (FLL.NE.2) WRITE (ILP,621) FLL IF (NSMOOTH.NE.2) WRITE (ILP,647) NSMOOTH C C READ AND PRINT CONDITIONS LIMITING POSSIBLE REFLECTIONS. C OPEN (UNIT=IO1,FILE='hklcond.dat',STATUS='UNKNOWN') WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,643) NCOND=0 7 READ (IO1,500,END=8) HCOND WRITE (ILP,644) HCOND I=0 CALL HKLCOND (HCOND,I) NCOND=NCOND+1 ICOND(NCOND)=I GO TO 7 8 IF (NCOND.EQ.0) THEN WRITE (ILP,645) CLOSE (UNIT=IO1,DISPOSE='DELETE') ELSE CLOSE (UNIT=IO1) END IF C C CONVERT FROM MICROSECONDS TO SECONDS. C TAU=TAU*1E-6 SIGTAU=SIGTAU*1E-6 SIGTIM=SIGTIM*1E-6 C C CALCULATE VARIANCE VALUES. C VARTAU=SIGTAU**2 VARATT=SIGATT**2 VARTHE=SIGTHE**2 VARTIM=SIGTIM**2 C C OPEN INPUT REFLECTION DATA FILE. C OPEN (UNIT=IO1,FILE=FILE1,STATUS='OLD',FORM='UNFORMATTED', & READONLY) 610 FORMAT ('1'/'0PROGRAM REFPK. ',A,1X,A,'. ',A) 611 FORMAT ('0INPUT FILE: ',A) 612 FORMAT ('0OUTPUT FILE: ',A) 613 FORMAT ('0LATTICE PARAMETERS'/'0',9X,'A',9X,'B',9X,'C',5X,'ALPHA', &6X,'BETA',5X,'GAMMA'/' ',3F10.4,3F10.3) 614 FORMAT ('0INTERNATIONAL TABLES VOL. IV DIFFRACTOMETER AXES') 619 FORMAT ('0BUSING AND LEVY DIFFRACTOMETER AXES') 615 FORMAT ('0SIEMENS (NEE NICOLET, NEE SYNTEX) P3 DIFFRACTOMETER') 616 FORMAT ('0ENRAF-NONIUS CAD4 DIFFRACTOMETER') 618 FORMAT ('0CRYSTAL MONOCHROMATOR, BRAGG ANGLE, GEOMETRY ANGLE, AN', &'D DYNAMIC DIFFRACTION (PERFECT CRYSTAL) FRACTION'/'0CRYSTAL T &HETA(M) RHO(M) FRAC(DYN) SIG(FRAC)'/' ',A2,8X,2F10.2,2E10.3) 629 FORMAT (' '/'$RADIATION: ',A) 630 FORMAT ('+ FILTER: ',A) 631 FORMAT ('0WAVELENGTHS'/'0FILT. EDGE ALPHA-ONE ALPHA-BAR ALPHA-TWO' &/' ',4F10.5/'0SPECTRAL LINE WIDTHS'/'0 ALPHA-ONE ', &' ALPHA-TWO'/' ',2(10X,F10.5)) 632 FORMAT ('0EXPECTED, LOWER LIMIT VALUES FOR T*TAN(THETA) SPECTRAL L &INE BROADENING COEFFICIENTS'/'0 T1 ', &' T2'/' ',2(10X,F10.4)) 622 FORMAT ('0INTENSITY OF ALPHA-TWO RELATIVE TO ALPHA-ONE'/'0 W &T2'/' ',F10.2) 617 FORMAT ('0ERROR ESTIMATES FOR DIFFRACTOMETER SCAN ANGLE (DEGREES T &HETA) AND SCAN-ANGLE DRIVE-PULSE TIMING (MICROSECONDS)'/'0 SIGT &HE SIGTIM'/' ',F10.4,F10.1) 641 FORMAT ('0COUNTER DEAD TIME (MICROSECONDS) AND BEAM ATTENUATOR FAC &TOR'/'0 TAU SIGTAU ATT SIGATT'/' ',4F10.2) 640 FORMAT ('0BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY & EXPOSURE TIMES (HOURS)'/'0 IIMIN IIMAX XTMI &N XTMAX'/' ',2I10,10X,2F10.2) 642 FORMAT ('0SCAN LIMITS (DECIMAL FRACTION OF SCAN WIDTH)'/'0 ', &' I1 I2'/' ',2F10.4) 648 FORMAT ('0WINDOW WIDTH FOR INITIAL SCAN FOR PEAK LIMITS (DECIMAL F &RACTION OF SCAN WIDTH)'/'0 WINDOW'/' ',F10.4) 646 FORMAT ('0E.S.D. MULTIPLIER FOR STATISTICAL SIGNIFICANCE TESTS'/'0 & CUTOFF'/' ',F10.2) 621 FORMAT ('0PROPORTIONALITY FACTOR FOR THE LEHMANN-LARSEN BIAS CORRE &CTION'/'0 FLL'/' ',F10.2) 647 FORMAT ('0NUMBER OF DATA-SMOOTHING PASSES FOR PEAK LOCATION'/'0 ', &' NSMOOTH'/' ',I10) 643 FORMAT ('0CONDITIONS LIMITING POSSIBLE REFLECTIONS'/) 644 FORMAT (' ',2X,A) 645 FORMAT (' NONE') RETURN END C----------------------------------------------------------------------- SUBROUTINE LEHLAR (L) C C LOCATE THE PEAK LIMITS BY THE METHOD OF LEHMANN AND LARSEN, WHICH C SEEKS THE TWO POINTS IN THE PROFILE BETWEEN WHICH INTEGRATION OF C THE NET INTENSITY GIVES MAXIMUM I/SIGMA(I). C C M. S. LEHMANN AND F. K. LARSEN (1974). ACTA CRYST. A30, 580-584. C C R. H. BLESSING, P. COPPENS, AND P. J. BECKER (1974). J. APPL. C CRYST. 7, 448-492. C C RETURN L = 0 TO FLAG AN ERROR CONDITION. C PARAMETER (NMAX=100) COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK6/ I1,I2,L1,L2,X1,X2,FLL DIMENSION Y(NMAX),V(NMAX) C C ANALYZE LOWER ANGLE ALPHA-ONE HALF-PEAK. C N=0 DO 1 I=I1,NINDEX(X1),+1 N=N+1 Y(N)=YY(I) V(N)=VY(I) 1 CONTINUE CALL LLIMIT (N,Y,V,FLL,L) IF (L.EQ.0) RETURN L1=I1+L-1 C C ANALYZE HIGHER ANGLE ALPHA-TWO HALF-PEAK. C N=0 DO 2 I=I2,NINDEX(X2),-1 N=N+1 Y(N)=YY(I) V(N)=VY(I) 2 CONTINUE CALL LLIMIT (N,Y,V,FLL,L) IF (L.EQ.0) RETURN L2=I2-L+1 RETURN END C----------------------------------------------------------------------- SUBROUTINE LLCORR (N,Y,F,L) C C CORRECT FOR THE MAXIMUM I/SIGMA(I) PEAK LIMIT BIAS. C C THE METHOD FINDS PEAK LIMITS THAT ARE SLIGHTLY WITHIN THE PEAKS. C LEHMANN AND LARSEN HAVE SHOWN THAT THE BIAS OF THE PEAK LIMITS IS C PROPORTIONAL TO 1/(2 + R), WHERE R IS THE PEAK-TO-BACKGROUND C RATIO. C C RETURN L = 0 TO FLAG AN ERROR CONDITION. C DIMENSION Y(N) P=0 B=0 DO 1 I=1,L-1 B=B+Y(I) 1 CONTINUE DO 2 I=L,N P=P+Y(I) 2 CONTINUE NP=N-L+1 NB=L-1 B=B*NP/NB R=(P-B)/B IF (R.LE.0) THEN L=0 RETURN END IF Q=F/(2+R) C C RMIN = 0, QMAX = 0.5*F. DEFAULT: F = 2.0, QMAX = 1. C L=L-NINT(Q*NP) IF (L.LT.0) L=0 RETURN END C----------------------------------------------------------------------- SUBROUTINE LLIMIT (N,Y,V,F,L) C C FIND MAXIMUM I/SIGMA(I) PEAK LIMIT. C DIMENSION Y(N),V(N) IF (N.LT.4) GO TO 9 P=0 B=0 VP=0 VB=0 XP=N XB=0 QMAX=0 DO 1 I=1,N P=P+Y(I) VP=VP+V(I) 1 CONTINUE DO 2 I=1,N-1 P=P-Y(I) B=B+Y(I) VP=VP-V(I) VB=VB+V(I) XP=XP-1 XB=XB+1 Q=XP/XB Q=(P-Q*B)/SQRT(VP+Q**2*VB) IF (Q.GT.QMAX) THEN QMAX=Q L=I END IF C C STORE THE Q = I/SIGMA(I) DATA IN THE NO-LONGER NEEDED V-ARRAY. C V(I)=Q 2 CONTINUE C C FIT A THREE- TO NINE-POINT PARABOLA TO THE I/SIGMA(I) DATA, AND C FIND ITS MAXIMUM. C M=MIN(L-1,N-L-1,4) IF (M.LT.1) GO TO 9 J=0 DO 11 I=L-M,L+M J=J+1 V(J)=V(I) 11 CONTINUE K=J DO 12 I=L-M,L+M K=K+1 V(K)=I 12 CONTINUE CALL PFIT (J,V(J+1),V(1)) C1=V(J+2) C2=V(J+3) IF (C2.GE.0) GO TO 9 L=NINT(-0.5*C1/C2) IF (L.LE.1.OR.L.GE.N-1) GO TO 9 CALL LLCORR (N,Y,F,L) RETURN 9 CONTINUE C C RETURN L = 0 TO FLAG AN ERROR CONDITION. C L=0 RETURN END C----------------------------------------------------------------------- SUBROUTINE MATINV (N,M,A,D) C C INVERT A SQUARE MATRIX BY GAUSS-JORDAN ELIMINATION, AND CALCULATE C ITS DETERMINANT. C C A - INPUT MATRIX, WHICH IS REPLACED BY ITS INVERSE C N - ORDER OF MATRIX C - LOGICAL SIZE OF ARRAY A C M - PHYSICAL SIZE OF ARRAY A, M .GE. N C D - DETERMINANT OF THE INPUT MATRIX C C RETURNS D = 0 TO FLAG A SINGULAR MATRIX. C C MATRIX IS SCALED TO AVOID ARITHMETIC OVERFLOW OR UNDERFLOW DURING C INVERSION. C C MATRIX SCALING USES A DIAGONAL MATRIX Q WITH C C Q(I,I) = 1/SQRT(ABS(A(I,I))). C C FORM B = Q A Q, C C B(I,J) = Q(I,I)*A(I,J)*Q(J,J), C C AND INVERT B TO OBTAIN B-1. THEN, C C B-1 = (Q A Q)-1 C = (A Q)-1 Q-1 C B-1 Q = Q-1 A-1 C Q B-1 Q = A-1 C C IN ADDITION, C C DET(A) = DET(B)/PROD(I=1,N) Q(I,I)**2 C C IF N EXCEEDS MAXN, THE WORK VECTORS Q, IK, AND JK MUST BE RE- C DIMENSIONED TO A LENGTH OF N. C PARAMETER (MAXN=25) DIMENSION A(M,M),Q(MAXN),IK(MAXN),JK(MAXN) DOUBLE PRECISION A,Q,AMAX,T C C STORE RECIPROCAL SQUARE-ROOT DIAGONAL MAGNITUDES FOR SCALING. C DO 10 I=1,N T=SQRT(ABS(A(I,I))) IF (T.EQ.0) THEN D=0 RETURN ELSE Q(I)=1/T END IF 10 CONTINUE C C SCALE MATRIX. C DO 11 I=1,N DO 11 J=1,N A(I,J)=A(I,J)*Q(I)*Q(J) 11 CONTINUE C C INVERT SCALED MATRIX AND CALCULATE ITS DETERMINANT. C D=1 DO 100 K=1,N C C FIND LARGEST ELEMENT A(I,J) IN REST OF MATRIX. C AMAX=0 DO 30 I=K,N DO 30 J=K,N IF (ABS(A(I,J)).GT.ABS(AMAX)) THEN AMAX=A(I,J) IK(K)=I JK(K)=J END IF 30 CONTINUE D=D*AMAX IF (D.EQ.0) RETURN C C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN A(K,K). C I=IK(K) IF (I.GT.K) THEN DO 50 J=1,N T=A(K,J) A(K,J)=A(I,J) A(I,J)=-T 50 CONTINUE END IF J=JK(K) IF (J.GT.K) THEN DO 60 I=1,N T=A(I,K) A(I,K)=A(I,J) A(I,J)=-T 60 CONTINUE END IF C C ACCUMULATE ELEMENTS OF INVERSE MATRIX. C DO 70 I=1,N IF (I.NE.K) A(I,K)=-A(I,K)/AMAX 70 CONTINUE DO 80 I=1,N DO 80 J=1,N IF (I.NE.K.AND.J.NE.K) A(I,J)=A(I,J)+A(I,K)*A(K,J) 80 CONTINUE DO 90 J=1,N IF (J.NE.K) A(K,J)=+A(K,J)/AMAX 90 CONTINUE A(K,K)=1/AMAX 100 CONTINUE C C RESTORE ORDERING OF MATRIX. C DO 130 K=N,1,-1 J=IK(K) IF (J.GT.K) THEN DO 110 I=1,N T=A(I,K) A(I,K)=-A(I,J) A(I,J)=T 110 CONTINUE END IF I=JK(K) IF (I.GT.K) THEN DO 120 J=1,N T=A(K,J) A(K,J)=-A(I,J) A(I,J)=T 120 CONTINUE END IF 130 CONTINUE C C SCALE INVERSE MATRIX. C DO 140 I=1,N DO 140 J=1,N A(I,J)=A(I,J)*Q(I)*Q(J) 140 CONTINUE C C SCALE DETERMINANT, IF POSSIBLE. C T=0 DO 150 I=1,N T=T+LOG(Q(I)) 150 CONTINUE T=LOG(ABS(D))-2*T C C TEST AGAINST RANGE OF MACHINE ALLOWED MAGNITUDES. C C EIGHT-BIT EXPONENT HAS MAXIMUM MAGNITUDE 2**7 = 128. C XMIN = 2**(-128) = 0.29E-38 LOG(XMIN) = -88.7 C XMAX = (2**(+128))/2 = 1.70E+38 LOG(XMAX) = +88.0 C IF (-87.LT.T.AND.T.LT.+87) D=(D/ABS(D))*EXP(T) RETURN END C----------------------------------------------------------------------- SUBROUTINE MATINV3 (A,B,D) C C INVERT A SYMMETRIC 3 X 3 MATRIX AND CALCULATE ITS DETERMINANT. C C A IS MATRIX, B IS INVERSE, D IS DETERMINANT OF A. C DIMENSION A(3,3),B(3,3) DO 1 I=1,3 K=1+MOD(I+1,3) M=1+MOD(I+3,3) DO 1 J=1,3 L=1+MOD(J+1,3) N=1+MOD(J+3,3) B(I,J)=A(L,K)*A(N,M)-A(L,M)*A(N,K) 1 CONTINUE D=0 DO 2 I=1,3 D=D+A(I,1)*B(1,I) 2 CONTINUE DO 3 I=1,3 DO 3 J=1,3 B(I,J)=B(I,J)/D 3 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE METRIC (A,GB) C C CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND RECIPROCAL C LATTICE METRIC TENSORS FROM DIRECT LATTICE PARAMETERS. C DIMENSION A(6),B(6),GA(3,3),GB(3,3) PI=ACOS(-1.0) DEG=180/PI RAD=PI/180 CA=COS(A(4)*RAD) CB=COS(A(5)*RAD) CG=COS(A(6)*RAD) SA=SIN(A(4)*RAD) SB=SIN(A(5)*RAD) SG=SIN(A(6)*RAD) V=A(1)*A(2)*A(3)*SQRT(1-CA**2-CB**2-CG**2+2*CA*CB*CG) B(1)=A(2)*A(3)*SA/V B(2)=A(1)*A(3)*SB/V B(3)=A(1)*A(2)*SG/V B(4)=(CB*CG-CA)/(SB*SG) B(5)=(CA*CG-CB)/(SA*SG) B(6)=(CA*CB-CG)/(SA*SB) GA(1,1)=A(1)*A(1) GA(1,2)=A(1)*A(2)*CG GA(1,3)=A(1)*A(3)*CB GA(2,1)=GA(1,2) GA(2,2)=A(2)*A(2) GA(2,3)=A(2)*A(3)*CA GA(3,1)=GA(1,3) GA(3,2)=GA(2,3) GA(3,3)=A(3)*A(3) GB(1,1)=B(1)*B(1) GB(1,2)=B(1)*B(2)*B(6) GB(1,3)=B(1)*B(3)*B(5) GB(2,1)=GB(1,2) GB(2,2)=B(2)*B(2) GB(2,3)=B(2)*B(3)*B(4) GB(3,1)=GB(1,3) GB(3,2)=GB(2,3) GB(3,3)=B(3)*B(3) B(4)=ACOS(B(4))*DEG B(5)=ACOS(B(5))*DEG B(6)=ACOS(B(6))*DEG RETURN END C----------------------------------------------------------------------- SUBROUTINE MV (N,A,B,X) C C MATRIX A TIMES COLUMN VECTOR B C DIMENSION A(N,N),B(N),X(N) DOUBLE PRECISION A DO 1 I=1,N X(I)=0 DO 1 J=1,N X(I)=X(I)+A(I,J)*B(J) 1 CONTINUE RETURN END C----------------------------------------------------------------------- FUNCTION NINDEX(X) C C FIND THE INDEX I OF THE ABSCISSA XX(I) NEAREST THE VALUE X. C C THE ABSCISSAE IN THE ARRAY XX(NSTEP) MUST BE ORDERED SUCH THAT C XX(I) INCREASES WITH I. C PARAMETER (NMAX=100) COMMON /BLOCK5/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) IF (X.LT.XX(1)) THEN NINDEX=0 RETURN END IF IF (X.GT.XX(NSTEP)) THEN NINDEX=NSTEP+1 RETURN END IF DO 1 I=2,NSTEP IF (XX(I).GE.X) THEN IF (X-XX(I-1).LT.XX(I)-X) THEN NINDEX=I-1 ELSE NINDEX=I END IF RETURN END IF 1 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE OUTPUT C C PRINTED OUTPUT COMPARING OBSERVED AND CALCULATED PEAK C POSITIONS AND PEAK LIMITS C PARAMETER (NMAX=100,M1=4,M2=16,M3=50) CHARACTER ATIME*8,ADATE*9,TITLE*80,HEADING*120 DIMENSION HEADING(M1),JXMIN(M1,M2,M3),JXMAX(M1,M2,M3),JX(M2) COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2 COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2, & Q1(3,3),Q2(3,3) DIMENSION H(3),Y(3),Z(3) DATA DEG, RAD /57.2957795, 0.017453293/ CHARACTER TYPE1*4,TYPE2*4 IF (MODEL1.EQ.1) TYPE1='LRNZ' IF (MODEL1.EQ.2) TYPE1='GAUS' IF (MODEL1.EQ.3) TYPE1='PLMC' IF (MODEL2.EQ.1) TYPE2='LRNZ' IF (MODEL2.EQ.2) TYPE2='GAUS' IF (MODEL2.EQ.3) TYPE2='PLMC' C C INITIALIZE ARRAYS FOR CATALOGING PROFILES WITH EXTREME C FEATURES. C DO 10 J1=1,M1 DO 10 J2=1,M2 DO 10 J3=1,M3 JXMIN(J1,J2,J3)=+1E6 10 JXMAX(J1,J2,J3)=-1E6 C C LOOP THROUGH DATA FILE OF LEAST-SQUARES OBSERVATIONS. C WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,635) NOBS=0 REWIND IO3 1 READ (IO3,END=9) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z, & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH IF (II.LT.0) GO TO 1 NOBS=NOBS+1 C C CALCULATE PEAK CENTROID. C DO 11 I=1,3 Y(I)=0 DO 11 J=1,3 11 Y(I)=Y(I)+H(J)*U(J,I) DO 12 I=1,3 12 Y(I)=Y(I)/DSTAR IF (TWOTH.LT.0) THEN DO 13 I=1,3 13 Y(I)=-Y(I) END IF COSPH=COS(PHI*RAD) SINPH=SIN(PHI*RAD) COSCH=COS(CHI*RAD) SINCH=SIN(CHI*RAD) SINOM=Y(1)*COSPH-Y(2)*SINPH COSOM=(Y(1)*SINPH+Y(2)*COSPH)*COSCH-Y(3)*SINCH DELTA=ARCTAN(SINOM,COSOM)*DEG-OMEGA DELTA=A180(DELTA) IF (TWOTH.LT.0) DELTA=-DELTA X0CALC=XM-DELTA C C CALCULATE PEAK LIMITS. C U1=0 U2=0 Q=0 DO 15 I=1,3 DO 15 J=1,3 U1=U1+Z(J)*Q1(I,J)*Z(I) 15 U2=U2+Z(J)*Q2(I,J)*Z(I) IF (TYPE1.EQ.'LRNZ') W1CALC=SQRT(U1)+T1*TANTH IF (TYPE1.EQ.'GAUS') W1CALC=SQRT(U1+(T1*TANTH)**2) IF (TYPE1.EQ.'PLMC') W1CALC=SQRT(U1-2*F*T1*TANTH/TANTHM & +T1*(TANTH/TANTHM)**2) IF (TYPE2.EQ.'LRNZ') W2CALC=SQRT(U2)+T2*TANTH IF (TYPE2.EQ.'GAUS') W2CALC=SQRT(U2+(T2*TANTH)**2) IF (TYPE2.EQ.'PLMC') W2CALC=SQRT(U2-2*F*T2*TANTH/TANTHM & +T2*(TANTH/TANTHM)**2) L1CALC=NINDEX(X1-W1CALC) L2CALC=NINDEX(X2+W2CALC) C C PRINT SELECTED REFLECTION RECORDS. C I=II J=H(1) K=H(2) L=H(3) CALL DPRINT (I,J,K,L,IPRINT) IF (IPRINT.EQ.1) THEN I0=NINDEX(X0) I0CALC=NINDEX(X0CALC) WRITE (ILP,636) I,J,K,L,I0,I0CALC,L1,L2,L1CALC,L2CALC END IF C C CATALOG PROFILES WITH EXTREME FEATURES. C JX(1)=I JX(2)=J JX(3)=K JX(4)=L JX(5)=NINDEX(X0) JX(6)=NINDEX(X0CALC) JX(7)=1000*(X0-XM) JX(8)=1000*(X0CALC-XM) JX(9)=1000*(X0-X0CALC) JX(10)=L1 JX(11)=L2 JX(12)=L1CALC JX(13)=L2CALC W=W1+X2-X1+W2 WCALC=W1CALC+X2-X1+W2CALC JX(14)=1000*W JX(15)=1000*WCALC JX(16)=1000*(W-WCALC) CALL EXTREM (JXMIN,JXMAX,JX) C C LOOP BACK FOR NEXT REFLECTION. C GO TO 1 C C END LOOP THROUGH REFLECTION DATA. C 9 CLOSE (UNIT=IO3) C C PRINT CATALOGS OF EXTREME CASES. C WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,649) DO 29 J1=1,M1 WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,645) HEADING(J1) DO 39 J3=1,MIN(M3,NOBS) I1=JXMIN(J1,1,J3) I2=JXMIN(J1,2,J3) I3=JXMIN(J1,3,J3) I4=JXMIN(J1,4,J3) I5=JXMIN(J1,5,J3) I6=JXMIN(J1,6,J3) F7=JXMIN(J1,7,J3)*0.001 F8=JXMIN(J1,8,J3)*0.001 F9=JXMIN(J1,9,J3)*0.001 I10=JXMIN(J1,10,J3) I11=JXMIN(J1,11,J3) I12=JXMIN(J1,12,J3) I13=JXMIN(J1,13,J3) F14=JXMIN(J1,14,J3)*0.001 F15=JXMIN(J1,15,J3)*0.001 F16=JXMIN(J1,16,J3)*0.001 IF (F8.LT.-999.99) F8=-999.99 IF (F8.GT.+999.99) F8=+999.99 IF (F9.LT.-999.99) F9=-999.99 IF (F9.GT.+999.99) F9=+999.99 39 WRITE (ILP,646) I1,I2,I3,I4,I5,I6,F7,F8,F9,I10,I11,I12,I13, & F14,F15,F16 WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,645) HEADING(J1) DO 49 J3=1,MIN(M3,NOBS) I1=JXMAX(J1,1,J3) I2=JXMAX(J1,2,J3) I3=JXMAX(J1,3,J3) I4=JXMAX(J1,4,J3) I5=JXMAX(J1,5,J3) I6=JXMAX(J1,6,J3) F7=JXMAX(J1,7,J3)*0.001 F8=JXMAX(J1,8,J3)*0.001 F9=JXMAX(J1,9,J3)*0.001 I10=JXMAX(J1,10,J3) I11=JXMAX(J1,11,J3) I12=JXMAX(J1,12,J3) I13=JXMAX(J1,13,J3) F14=JXMAX(J1,14,J3)*0.001 F15=JXMAX(J1,15,J3)*0.001 F16=JXMAX(J1,16,J3)*0.001 IF (F8.LT.-999.99) F8=-999.99 IF (F8.GT.+999.99) F8=+999.99 IF (F9.LT.-999.99) F9=-999.99 IF (F9.GT.+999.99) F9=+999.99 49 WRITE (ILP,646) I1,I2,I3,I4,I5,I6,F7,F8,F9,I10,I11,I12,I13, & F14,F15,F16 29 CONTINUE 610 FORMAT ('1'/'0PROGRAM REFPK. ',A,1X,A,'. ',A) 635 FORMAT ('0 I H K L X0O X0C L1O L2O L1C L2C &'/) 636 FORMAT (' ',I6,2X,3I4,3(4X,2I4)) 649 FORMAT ( &'0THE FOLLOWING TABLES LIST REFLECTIONS WITH EXTREME VALUES OF'/ &' SEVERAL CHARACTERISTIC PROPERTIES. THE TABLES ARE PRESENTED'/ &' IN PAIRS IN WHICH THE FIRST TABLE LISTS THE REFLECTIONS WITH'/ &' THE SMALLEST VALUES OF A GIVEN PROPERTY AND THE SECOND THOSE'/ &' WITH THE LARGEST VALUES. THE PROPERTY ON WHICH THE TABLE IS'/ &' BASED IS STATED IN THE HEADING ABOVE THE TABLE COLUMN'/ &' HEADINGS, AND THE ENTRIES IN THE TABLE ARE SORTED IN ORDER OF'/ &' INCREASING VALUE OF THE PROPERTY IN THE TABLES OF SMALLEST'/ &' VALUES, AND DECREASING VALUE IN THE TABLES OF LARGEST VALUES.') DATA HEADING / &'DX0O, DISPLACEMENT OF OBSERVED PEAK CENTROID FROM MIDPOINT OF SCA &N (DEGREES THETA)', &'D0=X0O-X0C, DIFFERENCE BETWEEN OBSERVED AND CALCULATED CENTROIDS &(DEGREES THETA)', &'WC, CALCULATED BASE WIDTH (DEGREES THETA)', &'WO-WC, DIFFERENCE BETWEEN OBSERVED AND CALCULATED FULL-PEAK BASE &WIDTH (DEGREES THETA)'/ 645 FORMAT ('0',A/ &'0 I H K L X0O X0C DX0O DX0C X0O-X0C', &' L1O L2O L1C L2C L2O-L1O L2C-L1C WO-WC' &/ &' D0', &' WO WC'/) 646 FORMAT (' ',I6,2X,3I4,4X,2I4,2X,3F8.2,4X,2I4,2X,2I4,2X,3F8.2) RETURN END C----------------------------------------------------------------------- SUBROUTINE PFIT (N,X,Y) C C FIT A PARABOLA, Y = C0 + C1*X + C2*X**2. C C RETURN C0, C1, AND C2 IN X(1), X(2), AND X(3), RESPECTIVELY. C DIMENSION X(N),Y(N),A(3,3),AINV(3,3),B(3) IF (N.LT.3) GO TO 9 X0=0.5*(X(1)+X(N)) IF (N.GT.3) GO TO 5 C C FOR THREE POINTS, FIND THE PARABOLA THROUGH THE POINTS. C X1=X(1)-X0 X2=X(2)-X0 X3=X(3)-X0 Y1=Y(1) Y2=Y(2) Y3=Y(3) A(1,1)=1 A(1,2)=X1 A(1,3)=X1*X1 A(2,1)=1 A(2,2)=X2 A(2,3)=X2*X2 A(3,1)=1 A(3,2)=X3 A(3,3)=X3*X3 B(1)=Y1 B(2)=Y2 B(3)=Y3 GO TO 2 5 CONTINUE C C FOR MORE THAN THREE POINTS, FIND THE LEAST-SQUARES PARABOLA. C SUMX=0 SUMX2=0 SUMX3=0 SUMX4=0 SUMY=0 SUMXY=0 SUMX2Y=0 DO 1 I=1,N XI=X(I)-X0 YI=Y(I) SUMX=SUMX+XI SUMX2=SUMX2+XI*XI SUMX3=SUMX3+XI*XI*XI SUMX4=SUMX4+XI*XI*XI*XI SUMY=SUMY+YI SUMXY=SUMXY+XI*YI SUMX2Y=SUMX2Y+XI*XI*YI 1 CONTINUE A(1,1)=N A(1,2)=SUMX A(1,3)=SUMX2 A(2,1)=SUMX A(2,2)=SUMX2 A(2,3)=SUMX3 A(3,1)=SUMX2 A(3,2)=SUMX3 A(3,3)=SUMX4 B(1)=SUMY B(2)=SUMXY B(3)=SUMX2Y 2 CALL MATINV3 (A,AINV,DETA) IF (DETA.EQ.0) GO TO 9 DO 3 I=1,3 X(I)=0 DO 3 J=1,3 X(I)=X(I)+AINV(I,J)*B(J) 3 CONTINUE C C ADJUST FOR SHIFTED ORIGIN. C X(1)=X(1)-X(2)*X0+X(3)*X0**2 X(2)=X(2)-2*X(3)*X0 RETURN 9 CONTINUE C C RETURN C2 = 0 TO FLAG AN ERROR CONDITION. C X(3)=0 RETURN END C----------------------------------------------------------------------- SUBROUTINE POSDEF (Q,ILP) DIMENSION Q(3,3) R=Q(1,1) S=Q(1,1)*Q(2,2)-Q(2,1)*Q(1,2) T=Q(1,1)*Q(2,2)*Q(3,3)+Q(1,2)*Q(2,3)*Q(3,1)+Q(1,3)*Q(2,1)*Q(3,2) & -Q(3,1)*Q(2,2)*Q(1,3)-Q(3,2)*Q(2,3)*Q(1,1)-Q(3,3)*Q(2,1)*Q(1,2) IF (MIN(R,S,T).GT.0) RETURN WRITE (ILP,6000) ((Q(I,J),J=1,3),I=1,3) T=(Q(1,1)+Q(2,2)+Q(3,3))/3 DO 1 I=1,3 DO 1 J=1,3 IF (I.EQ.J) THEN Q(I,J)=T ELSE Q(I,J)=0 END IF 1 CONTINUE WRITE (ILP,6001) ((Q(I,J),J=1,3),I=1,3) 6000 FORMAT (/'0ATTENTION! FITTED PEAK-WIDTH ANISOTROPY TENSOR WAS NOT & POSITIVE-DEFINITE,'/' ----------'/'0Q11 Q12 Q13 ',3E10.3/' Q2', &'1 Q22 Q23 = ',3E10.3/' Q31 Q32 Q33 ',3E10.3/'0AND HAS BEEN REPL &ACED BY ITS ISOTROPIC EQUIVALENT') 6001 FORMAT ('0Q11 Q12 Q13 ',3E10.3/' Q21 Q22 Q23 = ',3E10.3/' Q31 Q3 &2 Q33 ',3E10.3) RETURN END C----------------------------------------------------------------------- SUBROUTINE READ1 (IFILE,IEND,II,IH,IK,IL,ANGLES,WIDTH,SPEED,XTIME, & NSTEP,YY) DIMENSION ANGLES(4),YY(1) INTEGER*2 JH,JK,JL,JY(96) READ (IFILE,END=9) II,JH,JK,JL,ANGLES,WIDTH,SPEED,(X,I=1,5),XTIME, & JY IH=JH IK=JK IL=JL NSTEP=96 DO 1 I=1,NSTEP YY(I)=JY(I) 1 CONTINUE IEND=0 RETURN 9 IEND=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE READR (IFILE,IEND,II,IH,IK,IL,ANGLES,NSTEP,XX,YY,VX,VY, & XTIME) C C SUBROUTINE READR MUST RETURN: C C INTEGER*4 VARIABLES: C C II - MEASUREMENT SERIAL NUMBER C IH,IK,IL - REFLECTION INDICES C NSTEP - NUMBER OF SCAN STEPS C C REAL*4 ARRAYS: C C ANGLES(4) - DIFFRACTOMETER SETTING ANGLES C XX(NSTEP) - ROCKING-ANGLE STEPS C YY(NSTEP) - COUNT-RATE STEPS C VX(NSTEP) - ROCKING-ANGLE VARIANCES C VY(NSTEP) - COUNT-RATE VARIANCES C C REAL*4 VARIABLE: C C XTIME - RADIATION EXPOSURE TIME C C THE DIFFRACTOMETER ANGLES MUST BE IN DEGREES, AND THEIR ORDER IS C IMPORTANT: C C DIFF. ANGLES(1) (2) (3) (4) C C INT.TAB. TWO-THETA OMEGA CHI PHI C BUS.LEV. TWO-THETA OMEGA CHI PHI C P3 TWO-THETA OMEGA PHI CHI C CAD4 THETA PHI OMEGA KAPPA C C THE ROCKING-ANGLE STEPS MUST INCREASE (FROM NEGATIVE TO POSITIVE) C IN ORDER OF INCREASING ABSOLUTE VALUE OF THE BRAGG ANGLE, THETA. C C RECOMMENDED UNITS: C C ROCKING ANGLE - DEGREES (THETA) C COUNT RATES - COUNTS PER SECOND C EXPOSURE TIME - HOURS C DIMENSION ANGLES(4),XX(1),YY(1),VX(1),VY(1) IEND=0 CALL READ1 (IFILE,IEND,II,IH,IK,IL,ANGLES,WIDTH,SPEED,XTIME,NSTEP, & YY) IF (IEND.NE.0) RETURN C C CONVERT FROM STEP-SCAN COUNT PROFILE TO COUNT-RATE VERSUS ROCKING- C ANGLE PROFILE. C CALL XYDATA (WIDTH,SPEED,NSTEP,XX,YY,VX,VY) RETURN END C----------------------------------------------------------------------- FUNCTION SINTHL(H,GINV) DIMENSION H(3),GINV(3,3) Q=0 DO 1 I=1,3 DO 1 J=1,3 Q=Q+H(J)*GINV(I,J)*H(I) 1 CONTINUE SINTHL=0.5*SQRT(Q) RETURN END C----------------------------------------------------------------------- SUBROUTINE SMOOTH (NSMOOTH,N,Y) C C SMOOTHING OF EQUALLY SPACED DATA BY FOURTH DIFFERENCES C C EQUIVALENT TO A FIVE-POINT LEAST-SQUARES FIT OF A PARABOLA, C C Y = A0 + A1*X + A2*X**2, C C TO THE (I - 2)ND, (I - 1)ST, ITH, (I + 1)ST, AND (I + 2)ND C POINTS, I.E., TO X = -2, -1, 0, 1, AND 2. C C C. LANCZOS (1956). APPLIED ANALYSIS, PP. 316-320. PRENTICE- C HALL, ENGLEWOOD CLIFFS, NEW JERSEY. C DIMENSION Y(N),T(-2:+2) IF (NSMOOTH.LE.0.OR.N.LT.6) RETURN DO 9 I=1,NSMOOTH DO 1 J=-2,+2 T(J)=Y(3+J) 1 CONTINUE D3=T(1)-3*T(0)+3*T(-1)-T(-2) D4=T(2)-4*T(1)+6*T(0)-4*T(-1)+T(-2) Y(1)=Y(1)+D3/5+3*D4/35 Y(2)=Y(2)-2*D3/5-D4/7 Y(3)=Y(3)-3*D4/35 DO 2 J=4,N-2 DO 3 K=-2,+1 T(K)=T(K+1) 3 CONTINUE T(2)=Y(J+2) D4=T(2)-4*T(1)+6*T(0)-4*T(-1)+T(-2) Y(J)=Y(J)-3*D4/35 2 CONTINUE D3=T(2)-3*T(1)+3*T(0)-T(-1) Y(N-1)=Y(N-1)+2*D3/5-D4/7 Y(N)=Y(N)-D3/5+3*D4/35 DO 9 J=1,N IF (Y(J).LT.0) Y(J)=0 9 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE TFIT C C LEAST-SQUARES FIT OF ISOTROPIC, SCALAR COEFFICIENTS OF TAN(THETA) C SPECTRAL DISPERSION TERMS C COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2 COMMON /BLOCK4/ IDIFF,JMODEL,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2, & Q1(3,3),Q2(3,3) CHARACTER MODEL*4,TYPE1*4,TYPE2*4 DIMENSION H(3),Y(3),Z(3),XX(3) DIMENSION AL(2,2),B1L(2),B2L(2),X1L(2),X2L(2) DIMENSION AG(2,2),B1G(2),B2G(2),X1G(2),X2G(2) DIMENSION AM(3,3),B1M(3),B2M(3),X1M(3),X2M(3) DOUBLE PRECISION AG,AL,AM C C ZERO NORMAL MATRICES AND VECTORS. C NPAR=2 DO 20 I=1,NPAR DO 21 J=I,NPAR AL(I,J)=0 AG(I,J)=0 21 CONTINUE B1L(I)=0 B2L(I)=0 B1G(I)=0 B2G(I)=0 20 CONTINUE NPAR=3 DO 30 I=1,NPAR DO 31 J=I,NPAR AM(I,J)=0 31 CONTINUE B1M(I)=0 B2M(I)=0 30 CONTINUE C C LOOP THROUGH DATA FILE OF LEAST-SQUARES OBSERVATIONS TO BUILD C NORMAL EQUATIONS. C NOBS=0 REWIND IO3 1 READ (IO3,END=9) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z, & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH C C ACCUMULATE NORMAL MATRICES AND VECTORS. C NOBS=NOBS+1 IF (TANTHM.EQ.0.OR.RHOM.EQ.90) THEN C C BETA-FILTER OR PERPENDICULAR MONOCHROMATOR GEOMETRY C NPAR=2 C C LORENTZIAN CONVOLUTION MODEL, W = SQRT(Q) + T*TANTH C XX(1)=1 XX(2)=TANTH Y1=W1 Y2=W2 CALL ACCUM2 (NPAR,XX,Y1,Y2,AL,B1L,B2L) C C GAUSSIAN CONVOLUTION MODEL, W**2 = Q + T**2*(TANTH)**2 C XX(1)=1 XX(2)=TANTH**2 Y1=W1**2 Y2=W2**2 CALL ACCUM2 (NPAR,XX,Y1,Y2,AG,B1G,B2G) ELSE C C PARALLEL MONOCHROMATOR GEOMETRY. GAUSSIAN CONVOLUTION MODEL, C W**2 = Q - 2*F*T*TANTH/TANTHM + T*(TANTH/TANTHM)**2 , 1 < F < 2 C NPAR=3 XX(1)=1 XX(2)=-2*TANTH/TANTHM XX(3)=(TANTH/TANTHM)**2 Y1=W1**2 Y2=W2**2 CALL ACCUM2 (NPAR,XX,Y1,Y2,AM,B1M,B2M) END IF C C LOOP BACK FOR NEXT REFLECTION. C GO TO 1 C C END LOOP THROUGH REFLECTION DATA. C 9 CONTINUE C C SOLVE LEAST-SQUARES NORMAL EQUATIONS, AND CALCULATE STATISTICS OF C FIT AND PARAMETER ERRORS. C IF (TANTHM.EQ.0.OR.RHOM.EQ.90) THEN C C BETA-FILTER OR PERPENDICULAR MONOCHROMATOR C F=0 NPAR=2 IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE TFIT' C C LORENTZIAN MODEL. Q = X(1)**2, T = X(2) C CALL MATINV (NPAR,NPAR,AL,DET) IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE TFIT' CALL MV (NPAR,AL,B1L,X1L) CALL MV (NPAR,AL,B2L,X2L) MODEL='LRNZ' Q1L=X1L(1)**2 Q2L=X2L(1)**2 T1L=X1L(2) T2L=X2L(2) CALL AGREE1 (MODEL,IO3,Q1L,Q2L,T1L,T2L,0.0,0.0,WL0, & NOBS,NPAR,CHSQ1L,CHSQ2L,RMSD1L,RMSD2L) SIGT1L=SQRT(CHSQ1L*AL(2,2)) SIGT2L=SQRT(CHSQ2L*AL(2,2)) C C GAUSSIAN MODEL. Q = X(1), T = SQRT(X(2)) C CALL MATINV (NPAR,NPAR,AG,DET) IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE TFIT' CALL MV (NPAR,AG,B1G,X1G) CALL MV (NPAR,AG,B2G,X2G) MODEL='GAUS' Q1G=X1G(1) Q2G=X2G(1) T1G=0 T2G=0 IF (X1G(2).GT.0) T1G=SQRT(X1G(2)) IF (X2G(2).GT.0) T2G=SQRT(X2G(2)) CALL AGREE1 (MODEL,IO3,Q1G,Q2G,T1G,T2G,0.0,0.0,WL0, & NOBS,NPAR,CHSQ1G,CHSQ2G,RMSD1G,RMSD2G) SIGT1G=0 SIGT2G=0 IF (T1G.NE.0) SIGT1G=SQRT(CHSQ1G*AG(2,2))/(2*T1G) IF (T2G.NE.0) SIGT2G=SQRT(CHSQ2G*AG(2,2))/(2*T2G) ELSE C C PARALLEL MONOCHROMATOR GEOMETRY. Q = X(1), F*T = X(2), T = X(3) C NPAR=3 IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE TFIT' CALL MATINV (NPAR,NPAR,AM,DET) IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE TFIT' CALL MV (NPAR,AM,B1M,X1M) CALL MV (NPAR,AM,B2M,X2M) MODEL='PLMC' Q1M=X1M(1) Q2M=X2M(1) F1=X1M(2)/X1M(3) F2=X2M(2)/X2M(3) F=(F1+F2)/2 IF (F.LT.1) F=1 IF (F.GT.2) F=2 T1=X1M(3) T2=X2M(3) CALL AGREE1 (MODEL,IO3,Q1M,Q2M,T1,T2,TANTHM,F,WL0, & NOBS,NPAR,CHISQ1,CHISQ2,RMSD1,RMSD2) SIGF1=F1*SQRT(CHISQ1*((AM(2,2)/(F1*T1))**2+(AM(3,3)/T1)**2)) SIGF2=F2*SQRT(CHISQ2*((AM(2,2)/(F2*T2))**2+(AM(3,3)/T2)**2)) SIGF=SQRT(SIGF1**2+SIGF2**2)/2 SIGT1=SQRT(CHISQ1*AM(3,3)) SIGT2=SQRT(CHISQ2*AM(3,3)) END IF C C CHOOSE MODELS. C IF (TANTHM.EQ.0.OR.RHOM.EQ.90) THEN IF (JMODEL.EQ.0) THEN C C FOR EACH HALF-WIDTH, CHOOSE THE MODEL WITH SMALLER ROOT-MEAN- C SQUARE ERROR-OF-FIT. C RMSD1=MIN(RMSD1L,RMSD1G) RMSD2=MIN(RMSD2L,RMSD2G) IF (RMSD1.EQ.RMSD1L) TYPE1='LRNZ' IF (RMSD1.EQ.RMSD1G) TYPE1='GAUS' IF (RMSD2.EQ.RMSD2L) TYPE2='LRNZ' IF (RMSD2.EQ.RMSD2G) TYPE2='GAUS' ELSE IF (JMODEL.EQ.1.) THEN TYPE1='LRNZ' TYPE2='LRNZ' ELSE IF (JMODEL.EQ.2.) THEN TYPE1='GAUS' TYPE2='GAUS' END IF ELSE TYPE1='PLMC' TYPE2='PLMC' END IF IF (TYPE1.EQ.'LRNZ') MODEL1=1 IF (TYPE1.EQ.'GAUS') MODEL1=2 IF (TYPE1.EQ.'PLMC') MODEL1=3 IF (TYPE2.EQ.'LRNZ') MODEL2=1 IF (TYPE2.EQ.'GAUS') MODEL2=2 IF (TYPE2.EQ.'PLMC') MODEL2=3 IF (TYPE1.EQ.'PLMC') RETURN C C EVALUATE AND STORE T AND SIGMA(T) VALUES. C C EXPECTED, LOWER LIMIT VALUES T1 AND T2 WERE CALCULATED FOR C TRUNCATED LORENTZIAN SHAPED SPECTRAL LINES IN SUBROUTINE INPUT. C T1E=T1 T2E=T2 IF (TYPE1.EQ.'LRNZ') THEN T1=T1L SIGT1=SIGT1L ELSE IF (TYPE1.EQ.'GAUS') THEN T1=T1G SIGT1=SIGT1G END IF IF (TYPE2.EQ.'LRNZ') THEN T2=T2L SIGT2=SIGT2L ELSE IF (TYPE2.EQ.'GAUS') THEN T2=T2G SIGT2=SIGT2G END IF IF (ABS(T1-T1E).GT.1.96*SIGT1.OR.ABS(T2-T2E).GT.1.96*SIGT2) & WRITE (ILP,6000) T1E,T1,SIGT1,T2E,T2,SIGT2 IF (T1.LT.0.OR.T2.LT.0) WRITE (ILP,6001) IF (T1.LT.0) T1=0 IF (T2.LT.0) T2=0 6000 FORMAT (/'0ATTENTION! FITTED SPECTRAL LINE-BROADENING COEFFICIENT &(S) OF T*TAN(THETA) DIFFER SIGNIFICANTLY FROM EXPECTED VALUE(S).'/ &' ----------'/'0EXPECTED: T1 = ',F8.4,' FITTED: T1 = ',E10.3, &' SIGMA(T1) = ',E10.3/' T2 = ',F8.4,' T2 = & ',E10.3,' SIGMA(T2) = ',E10.3) 6001 FORMAT (/'0NEGATIVE T*TAN(THETA) COEFFIFIENT(S) HAVE BEEN RESET TO & T = 0.') RETURN END C----------------------------------------------------------------------- SUBROUTINE UQFIT C C LEAST-SQUARES FIT OF ADJUSTED DIFFRACTOMETER ORIENTATION MATRIX, C U, AND ANISOTROPIC PEAK WIDTH TENSORS, Q. C CHARACTER ATIME*8,ADATE*9,TITLE*80 COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2,T1,T2,SIGT1,SIGT2 COMMON /BLOCK4/ IDIFF,JTYPE,TANTHM,RHOM,F,U(3,3),MODEL1,MODEL2, & Q1(3,3),Q2(3,3) DIMENSION H(3),Y(3),Z(3) DIMENSION AU1(3,3),BU1(3),AU2(3,3),BU2(3),AU3(3,3),BU3(3) DIMENSION U1(3),U2(3),U3(3) DIMENSION XX(6),AA(6,6),B1(6),B2(6),XX1(6),XX2(6) DIMENSION SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3) DIMENSION AU(3,3),QCRYST(3,3),UCRYST(3,3) DOUBLE PRECISION AU1,AU2,AU3,AA CHARACTER TYPE1*4,TYPE2*4 IF (MODEL1.EQ.1) TYPE1='LRNZ' IF (MODEL1.EQ.2) TYPE1='GAUS' IF (MODEL1.EQ.3) TYPE1='PLMC' IF (MODEL2.EQ.1) TYPE2='LRNZ' IF (MODEL2.EQ.2) TYPE2='GAUS' IF (MODEL2.EQ.3) TYPE2='PLMC' C C ZERO LEAST-SQUARES NORMAL MATRICES AND VECTORS. C NPAR=3 DO 10 I=1,NPAR DO 11 J=1,NPAR AU1(I,J)=0 AU2(I,J)=0 AU3(I,J)=0 11 CONTINUE BU1(I)=0 BU2(I)=0 BU3(I)=0 10 CONTINUE NPAR=6 DO 20 I=1,NPAR DO 21 J=I,NPAR AA(I,J)=0 21 CONTINUE B1(I)=0 B2(I)=0 20 CONTINUE C C LOOP THROUGH DATA FILE OF LEAST-SQUARES OBSERVATIONS TO BUILD C NORMAL EQUATIONS. C NOBS=0 REWIND IO3 1 READ (IO3,END=9) II,H,TWOTH,OMEGA,CHI,PHI,DSTAR,Y,Z, & I1,I2,XM,X0,X1,X2,L1,L2,W1,W2,TANTH C C ACCUMULATE LEAST-SQUARES NORMAL MATRICES AND VECTORS. C NOBS=NOBS+1 C C U-MATRIX C NPAR=3 CALL ACCUM1 (NPAR,H,DSTAR*Y(1),AU1,BU1) CALL ACCUM1 (NPAR,H,DSTAR*Y(2),AU2,BU2) CALL ACCUM1 (NPAR,H,DSTAR*Y(3),AU3,BU3) C C Q-TENSORS C NPAR=6 K=0 DO 5 I=1,3 DO 5 J=I,3 K=K+1 Q=Z(I)*Z(J) IF (I.NE.J) Q=2*Q XX(K)=Q 5 CONTINUE C C LORENTZIAN OR GAUSSIAN CONVOLUTION MODEL FOR BETA-FILTERED C CHARACTERISTIC X-RAYS OR PERPENDICULAR MONOCHROMATOR GEOMETRY, C GAUSSIAN CONVOLUTION MODEL FOR PARALLEL MONOCHROMATOR GEOMETRY C IF (TYPE1.EQ.'LRNZ') Y1=(W1-T1*TANTH)**2 IF (TYPE1.EQ.'GAUS') Y1=W1**2-(T1*TANTH)**2 IF (TYPE1.EQ.'PLMC') Y1=W1**2-(-2*F*T1*TANTH/TANTHM+T1*(TANTH/ & TANTHM)**2) IF (TYPE2.EQ.'LRNZ') Y2=(W2-T2*TANTH)**2 IF (TYPE2.EQ.'GAUS') Y2=W2**2-(T2*TANTH)**2 IF (TYPE2.EQ.'PLMC') Y2=W2**2-(-2*F*T2*TANTH/TANTHM+T2*(TANTH/ & TANTHM)**2) CALL ACCUM2 (NPAR,XX,Y1,Y2,AA,B1,B2) C C LOOP BACK FOR NEXT REFLECTION. C GO TO 1 C C END LOOP THROUGH REFLECTION DATA. C 9 CONTINUE C C SOLVE LEAST-SQUARES NORMAL EQUATIONS, CALCULATE STATISTICS OF FIT, C AND ESTIMATE PARAMETER ERRORS, ASSUMING COVARIANCES ARE C NEGLIGIBLE. C C DIFFRACTOMETER ORIENTATION MATRIX C NPAR=3 IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE UQFIT' CALL MATINV (NPAR,NPAR,AU1,DET1) CALL MATINV (NPAR,NPAR,AU2,DET2) CALL MATINV (NPAR,NPAR,AU3,DET3) IF (DET1.EQ.0.OR.DET2.EQ.0.OR.DET3.EQ.0) & STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE UQFIT' CALL MV (NPAR,AU1,BU1,U1) CALL MV (NPAR,AU2,BU2,U2) CALL MV (NPAR,AU3,BU3,U3) DO 40 I=1,3 U(I,1)=U1(I) U(I,2)=U2(I) U(I,3)=U3(I) 40 CONTINUE CALL AGREE0 (IO3,U,NOBS,NPAR,CHISQ1,CHISQ2,CHISQ3,RMSDX0,RMSDX1, & RMSDTH) DO 41 I=1,3 SIGU(I,1)=SQRT(AU1(I,I)*CHISQ1) SIGU(I,2)=SQRT(AU2(I,I)*CHISQ2) SIGU(I,3)=SQRT(AU3(I,I)*CHISQ3) 41 CONTINUE C C ANISOTROPIC PEAK WIDTH TENSORS C NPAR=6 IF (NOBS.LE.NPAR) STOP 'NOT ENOUGH DATA IN SUBROUTINE UQFIT' CALL MATINV (NPAR,NPAR,AA,DET) IF (DET.EQ.0) STOP 'ILL-CONDITIONED MATRIX IN SUBROUTINE UQFIT' CALL MV (NPAR,AA,B1,XX1) CALL MV (NPAR,AA,B2,XX2) K=0 DO 60 I=1,3 DO 60 J=I,3 K=K+1 Q1(I,J)=XX1(K) Q2(I,J)=XX2(K) Q1(J,I)=Q1(I,J) Q2(J,I)=Q2(I,J) 60 CONTINUE CALL POSDEF (Q1,ILP) CALL POSDEF (Q2,ILP) CALL AGREE1 (TYPE1,IO3,Q1,Q2,T1,T2,TANTHM,F,WL0, & NOBS,NPAR,CHISQ1,DUMMY1,RMSD1,DUMMY2) CALL AGREE1 (TYPE2,IO3,Q1,Q2,T1,T2,TANTHM,F,WL0, & NOBS,NPAR,DUMMY1,CHISQ2,DUMMY2,RMSD2) K=0 DO 61 I=1,3 DO 61 J=I,3 K=K+1 SIGQ1(I,J)=SQRT(AA(K,K)*CHISQ1) SIGQ2(I,J)=SQRT(AA(K,K)*CHISQ2) SIGQ1(J,I)=SIGQ1(I,J) SIGQ2(J,I)=SIGQ2(I,J) 61 CONTINUE C C WRITE RESULTS TO BGLP.DAT FILE. C WRITE (IO2,601) ((U(I,J),J=1,3),I=1,3) WRITE (IO2,601) ((SIGU(I,J),J=1,3),I=1,3) WRITE (IO2,602) TYPE1 WRITE (IO2,601) ((Q1(I,J),J=1,3),I=1,3) WRITE (IO2,601) ((SIGQ1(I,J),J=1,3),I=1,3) WRITE (IO2,601) T1 WRITE (IO2,601) SIGT1 WRITE (IO2,602) TYPE2 WRITE (IO2,601) ((Q2(I,J),J=1,3),I=1,3) WRITE (IO2,601) ((SIGQ2(I,J),J=1,3),I=1,3) WRITE (IO2,601) T2 WRITE (IO2,601) SIGT2 WRITE (IO2,601) F C C IF R.M.S. ERROR IN CALCULATED PEAK CENTROID POSITIONS EXCEEDS C HALF THE R.M.S. PEAK DISPLACEMENT, SET DEFAULT PEAK POSITION TO C SCAN MIDPOINT. C IF (RMSDX1.GT.0.5*RMSDX0) WRITE (IO2,603) 0.5,0.0,0.0 601 FORMAT (3(1X,3E15.8/)) 602 FORMAT (1X,A) 603 FORMAT (1X,3F10.2) CLOSE (UNIT=IO2,DISPOSE='SAVE') C C PRINT RESULTS. C WRITE (ILP,610) ATIME,ADATE,TITLE NPAR=9 WRITE (ILP,620) NOBS,NPAR,RMSDX0,RMSDX1,RMSDTH WRITE (ILP,610) ATIME,ADATE,TITLE NPAR=7 IF (TYPE1.EQ.'PLMC') NPAR=8 WRITE (ILP,634) NOBS,NPAR,TYPE1,RMSD1,TYPE2,RMSD2,RMSDTH WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,621) ((U(I,J),J=1,3),I=1,3),((SIGU(I,J),J=1,3),I=1,3) WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,637) IF (TYPE1.EQ.'LRNZ') WRITE (ILP,647) IF (TYPE1.EQ.'GAUS') WRITE (ILP,648) IF (TYPE1.EQ.'PLMC') WRITE (ILP,658) F WRITE (ILP,632) ((Q1(I,J),J=1,3),I=1,3) WRITE (ILP,633) ((SIGQ1(I,J),J=1,3),I=1,3) WRITE (ILP,643) T1 WRITE (ILP,644) SIGT1 WRITE (ILP,638) IF (TYPE2.EQ.'LRNZ') WRITE (ILP,647) IF (TYPE2.EQ.'GAUS') WRITE (ILP,648) IF (TYPE2.EQ.'PLMC') WRITE (ILP,658) F WRITE (ILP,632) ((Q2(I,J),J=1,3),I=1,3) WRITE (ILP,633) ((SIGQ2(I,J),J=1,3),I=1,3) WRITE (ILP,643) T2 WRITE (ILP,644) SIGT2 WRITE (ILP,639) C C TRANSFORM PEAK WIDTH TENSORS TO CRYSTALLOGRAPHIC RECIPROCAL C LATTICE AXES, AND FIND EIGENVALUES AND EIGENVECTORS. C WRITE (ILP,610) ATIME,ADATE,TITLE DO 19 I=1,3 AI=1/SQRT(GINV(I,I)) DO 19 J=1,3 AU(I,J)=AI*U(I,J) 19 CONTINUE CALL UQUT (AU,Q1,QCRYST) WRITE (ILP,641) ((QCRYST(I,J),J=1,3),I=1,3) CALL EIGEN (3,3,QCRYST,UCRYST) WRITE (ILP,642) (QCRYST(I,I),(UCRYST(J,I),J=1,3),I=1,3) CALL UQUT (AU,Q2,QCRYST) WRITE (ILP,653) ((QCRYST(I,J),J=1,3),I=1,3) CALL EIGEN (3,3,QCRYST,UCRYST) WRITE (ILP,642) (QCRYST(I,I),(UCRYST(J,I),J=1,3),I=1,3) WRITE (ILP,654) 610 FORMAT ('1'/'0PROGRAM REFPK. ',A,1X,A,'. ',A) 620 FORMAT ('0STATISTICS OF FIT FOR PEAK POSITIONS:'/ &'0RMSD0 = SQRT(SUM((X0 - XM)**2)/NOBS)'/ &' RMSD1 = SQRT(SUM((X0 - X0CALC)**2)/(NOBS - NPAR))'/ &'0XM IS THE SCAN MIDPOINT.'/ &' X0 IS THE INTENSITY-WEIGHTED CENTROID OF THE PEAK ABOVE BACKGROU &ND.'/ &' X0CALC IS THE PEAK POSITION CALCULATED FROM THE REFLECTION INDIC &ES AND THE RECALCULATED ORIENTATION MATRIX.'/ &'0NOBS = ',I6/ &' NPAR = ',I6/ &'0RMSD0 = ',F6.3,' DEGREES THETA'/ &' RMSD1 = ',F6.3,' DEGREES THETA'/ &'0ROOT-MEAN-SQUARE STEP WIDTH:'/ &'0RMSDX = ',F6.3,' DEGREES THETA') 634 FORMAT ('0STATISTICS OF FIT FOR PEAK WIDTHS:'/ &'0RMSD(I) = SQRT[SUM([WOBS(I) - WCALC(I)]**2)/(NOBS - NPAR)]'/ &'0I = 1 OR 2 FOR THE LOWER ANGLE OR HIGHER ANGLE HALF-PEAK, RESPEC &TIVELY.'/ &'0THE WOBS ARE FROM THE LEHMANN-LARSEN PEAK LIMITS.'/ &' THE WCALC ARE FROM THE FITTED PEAK WIDTH TENSORS.'/ &'0NOBS = ',I6/ &' NPAR = ',I6/ &'0LOWER ANGLE HALF-PEAK, W(1):'/ &'0MODEL(1) = ',1X,A/ &' RMSD(1) = ',F6.3,' DEGREES THETA'/ &'0HIGHER ANGLE HALF-PEAK, W(2):'/ &'0MODEL(2) = ',1X,A/ &' RMSD(2) = ',F6.3,' DEGREES THETA'/ &'0ROOT-MEAN-SQUARE STEP WIDTH:'/ &'0RMSDX = ',F6.3,' DEGREES THETA') 621 FORMAT ('0PEAK POSITION MATRIX'/ &'0DIFFRACTOMETER ORIENTATION MATRIX RECALCULATED BY LEAST-SQUARES' &/' FIT TO OMEGA VALUES FOR THE INTENSITY-WEIGHTED PEAK CENTROIDS'/ &' ALONG WITH THE DIFFRACTOMETER ANGLES CHI AND PHI.'/'0DIFFRACTOME &TER AXES AS DEFINED BY WALTER C. HAMILTON (1974).'/' INTERNATIONA &L TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV, PP. 273-284.'/ &'0U(1,1) U(1,2) U(1,3) ',3F12.8/ &' U(2,1) U(2,2) U(2,3) = ',3F12.8/ &' U(3,1) U(3,2) U(3,3) ',3F12.8/ &'0SIG(1,1) SIG(1,2) SIG(1,3) ',3F12.8/ &' SIG(2,1) SIG(2,2) SIG(2,3) = ',3F12.8/ &' SIG(3,1) SIG(3,2) SIG(3,3) ',3F12.8) 637 FORMAT ('0PEAK WIDTH TENSORS'/'0LOWER ANGLE HALF-PEAK, W1') 638 FORMAT ('0HIGHER ANGLE HALF-PEAK, W2') 632 FORMAT ( &'0Q(1,1) Q(1,2) Q(1,3) ',3F8.4/ &' Q(2,1) Q(2,2) Q(2,3) = ',3F8.4/ &' Q(3,1) Q(3,2) Q(3,3) ',3F8.4) 633 FORMAT ( &'0SIG(1,1) SIG(1,2) SIG(1,3) ',3F8.4/ &' SIG(2,1) SIG(2,2) SIG(2,3) = ',3F8.4/ &' SIG(3,1) SIG(3,2) SIG(3,3) ',3F8.4) 643 FORMAT ('0 T = ',F8.4) 644 FORMAT ('0 SIG = ',F8.4) 647 FORMAT ('0LORENTZIAN MODEL'/ &'0W = (Z(K)*Q(J,K)*Z(J))**(1/2) + T*TAN(THETA)') 648 FORMAT ('0GAUSSIAN MODEL'/ &'0W = (Z(K)*Q(J,K)*Z(J) + (T*TAN(THETA))**2)**(1/2)') 658 FORMAT ('0GAUSSIAN MODEL FOR PARALLEL MONOCHROMATOR GEOMETRY'/ &'0W = (Z(K)*Q(J,K)*Z(J) - 2*F*T*TAN(THETA)/TAN(THETA(M)) + T*(TAN( &THETA)/TAN(THETA(M)))**2)**(1/2), F = ',F4.1) 639 FORMAT ('0W1 AND W2 ARE THE BASE WIDTHS OF THE HALF-PEAKS BELOW TH &ETA(ALPHA-1)'/' AND ABOVE THETA(ALPHA-2), RESPECTIVELY, IN UNITS O &F DEGREES OF CRYSTAL'/' ROTATION (THETA).'/'0THE Z(3) ARE THE COMP &ONENTS OF A UNIT VECTOR NORMAL TO THE INCIDENT'/' AND DIFFRACTED B &EAM VECTORS.'/'0THE Q(3,3) ARE THE ELEMENTS OF SYMMETRIC TENSORS T &HAT ARE RELATED TO'/' THE ANISOTROPY OF THE SIZE AND MOSAICITY O', &'F THE SPECIMEN CRYSTAL.'/'0THE T ARE (SCALAR) COEFFICIENTS PROPOR &TIONAL TO THE SPECTRAL BAND'/' WIDTH OF THE CHARACTERISTIC ALPHA-1 & AND ALPHA-2 LINES.'/'0THE VECTORS Z AND TENSORS Q ARE REFERRED TO & CRYSTAL-FIXED CARTESIAN AXES'/' THAT ARE COINCIDENT WITH THE DIFF &RACTOMETER AXES AT OMEGA = PHI = CHI = 0.') 641 FORMAT ('0ANISOTROPIC PEAK WIDTH TENSORS TRANSFORMED TO NORMALIZED &'/' CRYSTALLOGRAPHIC RECIPROCAL LATTICE AXES'/'0 Q(CRYST) = (A*U) &*Q*(A*U)-TRANSPOSE'/'0 Z(CRYST) = Z*(A*U)-INVERSE,'/'0WHERE THE M &ATRIX'/'0 (1/ASTAR 0 0 )'/' A = (0 1/BSTAR & 0 )'/' (0 0 1/CSTAR)'/'0NORMALIZES TO BASI &S VECTORS OF UNIT LENGTH.'//'0LOWER ANGLE HALF-PEAK, Q1'/'0QCRYST( &I,J)'//3(3F8.4/)) 642 FORMAT (' EIGENVALUE(I) EIGENVECTOR(I,J), (I DOWN, J ACROSS)'// &3(' 'F8.4,6X,3F8.4/)) 653 FORMAT ('0HIGHER ANGLE HALF-PEAK, Q2'/'0QCRYST(I,J)'//3(3F8.4/)) 654 FORMAT ('0NARROWEST REFLECTION PEAK OCCURS WHEN THE EIGENVECTOR OF & THE SMALLEST'/' EIGENVALUE IS PERPENDICULAR TO THE EQUATORIAL PLA &NE, AND WIDEST PEAK'/' WHEN THE EIGENVECTOR OF THE LARGEST EIGENVA &LUE IS PERPENDICULAR.'/'0WITH EIGENVECTOR(I,J) PERPENDICULAR TO TH &E EQUATORIAL PLANE, THE WIDTH'/' OF THE REFLECTION HALF-PEAK IS SQ &RT(EIGENVALUE(I)) (DEGREES THETA).') 651 FORMAT ('0FAILURE. NOBS = ',I2,' .LE. NPAR = ',I2) 652 FORMAT ('0FAILURE. DET = 0. NORMAL MATRIX SINGULAR.') RETURN END C----------------------------------------------------------------------- SUBROUTINE UQUT (U,Q,T) C C MATRIX TRIPLE PRODUCT, T = U*Q*U-TRANSPOSE C DIMENSION U(3,3),Q(3,3),T(3,3) DO 1 I=1,3 DO 1 J=1,3 T(I,J)=0 DO 1 K=1,3 DO 1 L=1,3 T(I,J)=T(I,J)+U(I,K)*Q(K,L)*U(J,L) 1 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE XCHECK (X1,X0,X2,X,Y,I1,I2) C C FOR WEAK, HIGH-ANGLE REFLECTIONS, SUBROUTINE CENTROID SOMETIMES C FINDS X0 TOO CLOSE TO THE ALPHA-ONE PEAK. ENSURE THAT X1 IS AT C THE ALPHA-ONE PEAK. C DIMENSION X(1),Y(1) D=(X(I2)-X(I1))/(I2-I1+1) IF (X0-X1.LE.D) RETURN J1=NINDEX(X1) J0=NINDEX(X0) J1=MAX(J1,I1+1) J0=MIN(J0,I2-1) IF (J0-J1.LE.1) RETURN YMAX=0 DO 1 I=J1,J0 YI=Y(I-1)+Y(I)+Y(I+1) IF (YI.GT.YMAX) THEN YMAX=YI IMAX=I END IF 1 CONTINUE D=X(IMAX)-X1 IF (D.GT.0) THEN X1=X1+D X0=X0+D X2=X2+D END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE XYDATA (WIDTH,SPEED,NSTEP,XX,YY,VX,VY) C C CONVERT FROM STEP-SCAN COUNT PROFILE TO COUNT-RATE VERSUS ROCKING- C ANGLE PROFILE, AND APPLY COINCIDENCE CORRECTIONS FOR COUNTING C LOSSES. C C RETURN: C C ROCKING-ANGLE STEPS XX(NSTEP) IN DEGREES (THETA) INCREASING C (FROM NEGATIVE TO POSITIVE) IN ORDER OF INCREASING ABSOLUTE C VALUE OF THE BRAGG ANGLE, THETA; C C COUNT-RATE STEPS YY(NSTEP) IN COUNTS PER SECOND; C C VARIANCE STEPS VX(NSTEP) AND VY(NSTEP). C C THIS VERSION OF SUBROUTINE XYDATA TREATS OMEGA, OMEGA/THETA, OR C OMEGA/TWO-THETA STEP SCANS WITH: C C NSTEP EQUALLY SPACED STEPS, C ----- C TOTAL SCAN WIDTH IN DEGREES OF OMEGA ROTATION, C ----- C SCAN SPEED IN DEGREES (OMEGA) PER MINUTE. C ----- C C NEGATIVE SCAN WIDTH IS A FLAG THAT THE ATTENUATOR WAS USED. C DIMENSION XX(1),YY(1),VX(1),VY(1) COMMON /BLOCK7/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT C C ROCKING-ANGLE ABSCISSAE C X0=0.5*(1+NSTEP) DX=ABS(WIDTH)/NSTEP VARD=2*VARTHE/NSTEP**2 DO 1 I=1,NSTEP XX(I)=(I-X0)*DX VX(I)=(I-X0)**2*VARD 1 CONTINUE C C COUNT-RATE ORDINATES C SPEED=SPEED/60 T=(ABS(WIDTH)/NSTEP)/SPEED VART=2*(2*VARTHE/SPEED**2+VARTIM)/NSTEP**2 DO 2 I=1,NSTEP Y=YY(I) C C ALLOW FOR STEPS WITH VERY FEW OR NO COUNTS. C Y=Y+0.5 VARY=Y X=Y/T VARX=(VARY+X**2*VART)/T**2 C C DEAD TIME CORRECTION C Y=X/(1-TAU*X) VARY=Y**4*(VARX/X**4+VARTAU) C C ATTENUATOR CORRECTION C IF (WIDTH.LT.0) THEN X=Y VARX=VARY Y=ATT*X VARY=X**2*VARATT+ATT**2*VARX END IF YY(I)=Y VY(I)=VARY 2 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE YZ (DELTA,TWOTH,OMEGA,CHI,PHI,Y,Z) C C GET THE COMPONENTS ALONG CARTESIAN AXES FIXED IN THE CRYSTAL OF A C UNIT VECTOR, Y, PARALLEL TO THE RECIPROCAL LATTICE DIFFRACTION C VECTOR, DSTAR, AND OF A UNIT VECTOR, Z, NORMAL TO THE EQUATORIAL C PLANE OF THE INCIDENT AND DIFFRACTED BEAM VECTORS. C C AXES AS DEFINED BY WALTER C. HAMILTON (1974). INTERNATIONAL C TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV, PP. 273-284. C C THE CRYSTAL-FIXED CARTESIAN AXES ARE COINCIDENT WITH THE C DIFFRACTOMETER AXES AT OMEGA = PHI = CHI = 0. C C DELTA IS THE DISPLACEMENT OF THE PEAK CENTROID FROM THE SCAN C MIDPOINT. C C THE POSITIVE SENSE OF THE OMEGA ROTATION IS DEFINED SUCH THAT C DELTA(OMEGA) = -DELTA(X0). C DIMENSION Y(3),Z(3) DATA RAD /0.017453293/ IF (TWOTH.LT.0) DELTA=-DELTA COSOM=COS((OMEGA-DELTA)*RAD) SINOM=SIN((OMEGA-DELTA)*RAD) COSPH=COS(PHI*RAD) SINPH=SIN(PHI*RAD) COSCH=COS(CHI*RAD) SINCH=SIN(CHI*RAD) Y(1)= COSPH*SINOM+SINPH*COSCH*COSOM Y(2)=-SINPH*SINOM+COSPH*COSCH*COSOM Y(3)=-SINCH*COSOM IF (TWOTH.LT.0) THEN DO 1 I=1,3 Y(I)=-Y(I) 1 CONTINUE END IF Z(1)=SINPH*SINCH Z(2)=COSPH*SINCH Z(3)=COSCH RETURN END C-----------------------------------------------------------------------