PROGRAM BGLP 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 /BLOCKC/ NCOND,ICOND(38) COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2 COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2, & XX0,XW1,XW2 COMMON /BLOCK4/ C0KIN,C1KIN,C2KIN,C0DYN,C1DYN,C2DYN,FRACTD,VARFD, & INEUTRON COMMON /BLOCK5/ IDIFF,U(3,3),MODEL1,MODEL2,Q1(3,3),Q2(3,3),T1,T2, & F,SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3),SIGT1,SIGT2,TANTHM COMMON /BLOCK6/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS CALL INPUT CALL PROCESS CALL OUTPUT STOP 'PROGRAM BGLP FINIS!' END C----------------------------------------------------------------------- SUBROUTINE BGLP1 C C LINEAR BACKGROUND C PARAMETER (NMAX=100) DIMENSION X(NMAX),Y(NMAX) COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS DATA RAD /0.017453293/ C C FIT A STRAIGHT LINE TO THE PROFILE IN THE TWO BACKGROUND REGIONS, C REQUIRING THAT, FOR THE PURPOSE OF FITTING THE LINE, EACH C BACKGROUND SAMPLE INCLUDES AT LEAST TWO POINTS. C J1=MAX(L1-1,I1+1) J2=MIN(L2+1,I2-1) N=0 DO 1 I=I1,J1 N=N+1 X(N)=XX(I) Y(N)=YY(I) 1 CONTINUE DO 2 I=J2,I2 N=N+1 X(N)=XX(I) Y(N)=YY(I) 2 CONTINUE CALL FITLINE (N,X,Y,C0,C1,V0,V1,COV) C C SUBTRACT LINEAR BACKGROUND POINT-BY-POINT, DIVIDE OUT LORENTZ AND C POLARIZATION FACTORS POINT-BY-POINT, AND INTEGRATE THE PEAK C PROFILE BETWEEN PEAK LIMITS. C R=0 VR=0 DO 9 I=MAX(I1+1,L1),MIN(L2,I2-1) C C BACKGROUND SUBTRACTION C XI=XX(I) YI=YY(I)-C0-C1*XI VI=VY(I)+V0+V1*XI**2+2*XI*COV C C LORENTZ AND POLARIZATION CORRECTIONS C TH=(XX(I)+TH0)*RAD VT=VX(I)*RAD**2 CALL LZ (TH,VT,FL,VL) CALL PZ (TH,VT,FP,VP) CV=SQRT(VL*VP) YJ=YI IF (YI.EQ.0) YI=0.5*SQRT(VI) VI=(VI/YI**2+VL/FL**2+VP/FP**2+CV/(FL*FP))*(YI/(FL*FP))**2 YI=YJ/(FL*FP) C C INTEGRATED NET INTENSITY BY TRAPEZOIDAL INTEGRATION C DX=0.50*(XX(I+1)-XX(I-1)) VD=0.25*(VX(I+1)+VX(I-1)) R=R+YI*DX VR=VR+YI**2*VD+DX**2*VI 9 CONTINUE SIGR=SQRT(VR) RETURN END C----------------------------------------------------------------------- SUBROUTINE BGLP2 C C STRUCTURED BACKGROUND UNDER LOW ANGLE REFLECTIONS DUE TO THE C ABSORPTION EDGE OF THE BETA FILTER AND THE PEAK IN THE WHITE C RADIATION BACKGROUND C C R. J. NELMES (1975). ACTA CRYST. A31, 273-279. C PARAMETER (NMAX=100) DIMENSION X(NMAX),Y(NMAX) COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS DATA RAD /0.017453293/ C C FIT STRAIGHT LINE SEGMENTS TO THE STRUCTURED BACKGROUND. C A0=0 A1=0 VA0=0 VA1=0 COVA=0 E0=0 E1=0 VE0=0 VE1=0 COVE=0 B0=0 B1=0 VB0=0 VB1=0 COVB=0 J1=NINDEX(XE-EPS) J2=NINDEX(XE+EPS) IF (J2.LE.I1) GO TO 99 C C FIT ONE STRAIGHT LINE TO THE LOW ANGLE BACKGROUND BELOW THE TAIL C OF THE BETA FILTER ABSORPTION EDGE. C J1=MIN(J1,L1-1) J1=MAX(J1,I1) N=0 DO 1 I=I1,J1 N=N+1 X(N)=XX(I) Y(N)=YY(I) 1 CONTINUE IF (N.GE.3) THEN CALL FITLINE (N,X,Y,A0,A1,VA0,VA1,COVA) ELSE IF (N.EQ.2) THEN A0=0.5*(Y(1)+Y(2)) VA0=0.5*A0 ELSE A0=Y(1) VA0=A0 END IF C C FIT A SECOND STRAIGHT LINE TO THE ABSORPTION EDGE. C J2=MIN(J2,L1-1) J2=MAX(J2,J1) N=0 DO 2 I=J1,J2 N=N+1 X(N)=XX(I) Y(N)=YY(I) 2 CONTINUE IF (N.GE.3) THEN CALL FITLINE (N,X,Y,E0,E1,VE0,VE1,COVE) ELSE IF (N.EQ.2) THEN E0=0.5*(Y(1)+Y(2)) VE0=0.5*E0 ELSE E0=Y(1) VE0=E0 END IF C C THE ABSORPTION EDGE SHOULD HAVE POSITIVE INTERCEPT AND SLOPE. C IF (E0.LT.0.OR.E1.LT.0) GO TO 99 C C FIT A THIRD STRAIGHT LINE TO THE HIGH ANGLE, WHITE RADIATION C BACKGROUND ABOVE THE TAIL OF THE BRAGG PEAK. C J2=MAX(J2,L2+1) J2=MIN(J2,I2) N=0 DO 3 I=J2,I2 N=N+1 X(N)=XX(I) Y(N)=YY(I) 3 CONTINUE IF (N.GE.3) THEN CALL FITLINE (N,X,Y,B0,B1,VB0,VB1,COVB) ELSE IF (N.EQ.2) THEN B0=0.5*(Y(1)+Y(2)) VB0=0.5*B0 ELSE B0=Y(1) VB0=B0 END IF C C THE HIGH ANGLE BACKGROUND SHOULD HAVE POSITIVE INTERCEPT AND C NEGATIVE SLOPE. C IF (B0.LT.0.OR.B1.GT.0) GO TO 99 C C EVALUATE THE LOW ANGLE BACKGROUND AT XE - EPSILON AND THE HIGH C ANGLE BACKGROUND AT XE + EPSILON. C J1=NINDEX(XE-EPS) J2=NINDEX(XE+EPS) J1=MAX(J1,I1) J2=MIN(J2,I2) X1=XX(J1) X2=XX(J2) Y1=A0+A1*X1 Y2=B0+B1*X2 V1=VA0+VA1*X1**2+X1*COVA V2=VB0+VB1*X2**2+X2*COVB C C THE BACKGROUND AT XE - EPSILON SHOULD BE LOWER THAN THE C BACKGROUND AT XE + EPSILON. C IF (Y1.GE.Y2) GO TO 99 C C CALCULATE THE COEFFICIENTS OF THE LINE JOINING THE LOW ANGLE C BACKGROUND AT X1 = XE - EPSILON TO THE HIGH ANGLE BACKGROUND AT C X2 = XE + EPSILON. C IF (J1.LT.L1.AND.L1.LT.J2) THEN C C THE LOW ANGLE PEAK LIMIT FALLS ON THE ABSORPTION EDGE. C X1=XX(L1-1) Y1=E0+E1*X1 V1=VE0+VE1*X1**2+X1*COVE END IF E1=(Y2-Y1)/(X2-X1) E0=Y1-E1*X1 VE1=(V1+V2)/(X2-X1)**2 VE0=(V1*X2**2+V2*X1**2)/(X2-X1)**2 COVE=-(V1*X2+V2*X1)/(X2-X1)**2 C C SUBTRACT BACKGROUND POINT-BY-POINT, DIVIDE OUT LORENTZ AND C POLARIZATION FACTORS POINT-BY-POINT, AND INTEGRATE THE PEAK C PROFILE BETWEEN PEAK LIMITS. C R=0 VR=0 DO 9 I=MAX(I1+1,L1),MIN(L2,I2-1) C C BACKGROUND SUBTRACTION C XI=XX(I) IF (XI.LT.X1) THEN C0=A0 C1=A1 V0=VA0 V1=VA1 COV=COVA ELSE IF (X1.LE.XI.AND.XI.LE.X2) THEN C0=E0 C1=E1 V0=VE0 V1=VE1 COV=COVE ELSE C0=B0 C1=B1 V0=VB0 V1=VB1 COV=COVB END IF YI=YY(I)-C0-C1*XI VI=VY(I)+V0+V1*XI**2+2*XI*COV C C PREVENT TRUNCATION OF THE LOW ANGLE TAIL OF THE PEAK BY THE C ESTIMATED ABSORPTION EDGE. C IF (XI.LT.0.AND.YI.LT.0) YI=0 C C LORENTZ AND POLARIZATION CORRECTIONS C TH=(XX(I)+TH0)*RAD VT=VX(I)*RAD**2 CALL LZ (TH,VT,FL,VL) CALL PZ (TH,VT,FP,VP) CV=SQRT(VL*VP) YJ=YI IF (YI.EQ.0) YI=0.5*SQRT(VI) VI=(VI/YI**2+VL/FL**2+VP/FP**2+CV/(FL*FP))*(YI/(FL*FP))**2 YI=YJ/(FL*FP) C C INTEGRATED NET INTENSITY BY TRAPEZOIDAL INTEGRATION C DX=0.50*(XX(I+1)-XX(I-1)) VD=0.25*(VX(I+1)+VX(I-1)) R=R+YI*DX VR=VR+YI**2*VD+DX**2*VI 9 CONTINUE SIGR=SQRT(VR) IF (R.LT.CUTOFF*SIGR) GO TO 99 RETURN 99 CONTINUE C C IF STRUCTURED BACKGROUND MODEL IS NOT APPROPRIATE, RESORT TO C LINEAR BACKGROUND. C J2=MIN(NINDEX(XE+EPS),L1-1) I1=MAX(I1,J2) CALL BGLP1 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 DPRINT (JJ,JH,JK,JL,JPRINT) C C JPRINT = 1 FOR H00, 0K0, 00L, HH0, 0KK, H0H, AND HHH; C JPRINT = 0 OTHERWISE. C IF (JJ.LT.0) GO TO 10 J=ABS(JH) K=ABS(JK) L=ABS(JL) 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 JPRINT=0 RETURN 11 JPRINT=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE EXTREM (IXMIN,IXMAX,IX) PARAMETER (N1=4,N2=12,N3=50) DIMENSION IXMIN(N1,N2,N3),IXMAX(N1,N2,N3),ITYPE(N1),IX(N2) DATA ITYPE /5, 6, 8, 12/ 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(1),V(1) L=NINT(WINDOW*N) DO 1 I=L+1,N-L YB=0 YI=0 VB=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(1),Y(1) 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 C C 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 KAPPA IS MECHANICALLY LIMITED TO THE RANGE 0 TO 180 DEGREES. C WITH ALPHA = 50 DEGREES, THIS LIMITS CHI TO THE RANGE -100 TO C +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 I=0 ELSE H=A(I) END IF RETURN 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,L+H=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 IABSENT = 1 FOR A FORBIDDEN REFLECTION. C INTEGER H DIMENSION ICOND(38) IABSENT=0 IF (NCOND.EQ.0) RETURN 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.OR.MOD(L+H,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 RETURN 91 IABSENT=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE INPUT C C READ CONTROL DATA AND EVALUATE VARIABLES. C CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE1*80,FILE2*80,HCOND*18, & ZTIME*8,ZDATE*9,TYPE1*4,TYPE2*4 COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCKC/ NCOND,ICOND(38) COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2 COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2, & XX0,XW1,XW2 COMMON /BLOCK4/ C0KIN,C1KIN,C2KIN,C0DYN,C1DYN,C2DYN,FRACTD,VARFD, & INEUTRON COMMON /BLOCK5/ IDIFF,U(3,3),MODEL1,MODEL2,Q1(3,3),Q2(3,3),T1,T2, & F,SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3),SIGT1,SIGT2,TANTHM COMMON /BLOCK6/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT DIMENSION CELL(6) PI=ACOS(-1.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='bglp.dat',STATUS='OLD') READ (IO1,500) ZTIME READ (IO1,500) ZDATE READ (IO1,500) TITLE READ (IO1,500) FILE1 READ (IO1,500) FILE2 READ (IO1,501) CELL READ (IO1,502) IDIFF,INEUTRON READ (IO1,501) THM,RHOM,FRACTD,SIGFD READ (IO1,501) WLE,WL1,WL0,WL2 READ (IO1,501) SIGTHE,SIGTIM READ (IO1,501) TAU,SIGTAU,ATT,SIGATT READ (IO1,501) WINDOW READ (IO1,501) CUTOFF READ (IO1,502) NSMOOTH READ (IO1,505) SC1,SC2 READ (IO1,503) II1,II2,XT1,XT2 READ (IO1,504) ((U(I,J),J=1,3),I=1,3) READ (IO1,504) ((SIGU(I,J),J=1,3),I=1,3) READ (IO1,500) TYPE1 READ (IO1,504) ((Q1(I,J),J=1,3),I=1,3) READ (IO1,504) ((SIGQ1(I,J),J=1,3),I=1,3) READ (IO1,504) T1 READ (IO1,504) SIGT1 READ (IO1,500) TYPE2 READ (IO1,504) ((Q2(I,J),J=1,3),I=1,3) READ (IO1,504) ((SIGQ2(I,J),J=1,3),I=1,3) READ (IO1,504) T2 READ (IO1,504) SIGT2 F=0 XX0=0.5 XW1=0 XW1=0 READ (IO1,504,END=9) F READ (IO1,505,END=9) XX0,XW1,XW2 500 FORMAT (1X,A) 501 FORMAT (1X,6F10.5) 502 FORMAT (1X,2I10) 503 FORMAT (1X,2I10,2F10.3) 504 FORMAT (3(1X,3E15.8/)) 505 FORMAT (1X,4F10.2) 9 CLOSE (UNIT=IO1) 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 C C PRINT CONTROL DATA. C OPEN (UNIT=ILP,FILE='bglp.lp',STATUS='NEW',CARRIAGECONTROL= & 'FORTRAN') WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,600) ZTIME,ZDATE,TITLE WRITE (ILP,611) FILE1 WRITE (ILP,612) FILE2 WRITE (ILP,629) 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) IF (INEUTRON.EQ.0) WRITE (ILP,651) IF (INEUTRON.EQ.1) WRITE (ILP,652) IF (THM.NE.0) WRITE (ILP,628) THM,RHOM,FRACTD,SIGFD WRITE (ILP,631) WLE,WL1,WL0,WL2 WRITE (ILP,613) SIGTHE,SIGTIM WRITE (ILP,620) TAU,SIGTAU,ATT,SIGATT WRITE (ILP,622) II1,II2,XT1,XT2 IF (SC1.GT.0.OR.SC2.LT.1) WRITE (ILP,621) SC1,SC2 WRITE (ILP,634) WINDOW WRITE (ILP,626) CUTOFF IF (NSMOOTH.NE.2) WRITE (ILP,627) NSMOOTH WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,623) ((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,617) IF (TYPE1.EQ.'LRNZ') WRITE (ILP,609) IF (TYPE1.EQ.'GAUS') WRITE (ILP,608) 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,642) T1 WRITE (ILP,643) SIGT1 WRITE (ILP,618) IF (TYPE2.EQ.'LRNZ') WRITE (ILP,609) IF (TYPE2.EQ.'GAUS') WRITE (ILP,608) 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,642) T2 WRITE (ILP,643) SIGT2 C C CALCULATE RECIPROCAL SPACE METRIC TENSOR. C CALL METRIC (CELL,GINV) C C EVALUATE MONOCHROMATOR VARIABLES. C THM=THM*PI/180 RHOM=RHOM*PI/180 COSSQ=COS(RHOM)**2 SINSQ=SIN(RHOM)**2 C0KIN=COS(2*THM)**2 C1KIN=COSSQ+SINSQ*C0KIN C2KIN=SINSQ+COSSQ*C0KIN C0DYN=ABS(COS(2*THM)) C1DYN=COSSQ+SINSQ*C0DYN C2DYN=SINSQ+COSSQ*C0DYN IF (FRACTD.NE.0.AND.SIGFD.EQ.0) SIGFD=0.05*FRACTD TANTHM=TAN(THM) 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 VARFD=SIGFD**2 VARTAU=SIGTAU**2 VARATT=SIGATT**2 VARTHE=SIGTHE**2 VARTIM=SIGTIM**2 C C SET PEAK POSITION AND PEAK WIDTHS TO BE USED INSTEAD OF VALUES C CALCULATED FROM THE PEAK POSITION AND PEAK WIDTH PARAMETERS. C IF (XX0.NE.0.OR.XW1.NE.0.OR.XW2.NE.0) & WRITE (ILP,610) ATIME,ADATE,TITLE IF (XX0.NE.0) WRITE (ILP,647) XX0 IF (XW1.NE.0.OR.XW2.NE.0) WRITE (ILP,648) XW1,XW2 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,644) NCOND=0 7 READ (IO1,506,END=8) HCOND 506 FORMAT (A) WRITE (ILP,645) HCOND I=0 CALL HKLCOND (HCOND,I) NCOND=NCOND+1 ICOND(NCOND)=I GO TO 7 8 IF (NCOND.EQ.0) THEN WRITE (ILP,646) CLOSE (UNIT=IO1,DISPOSE='DELETE') ELSE CLOSE (UNIT=IO1) END IF C C OPEN INPUT AND OUTPUT REFLECTION DATA FILES. C OPEN (UNIT=IO1,FILE=FILE1,STATUS='OLD',FORM='UNFORMATTED', & READONLY) OPEN (UNIT=IO2,FILE=FILE2,STATUS='NEW',FORM='UNFORMATTED') OPEN (UNIT=IO3,FILE='OUTPUT.DAT',STATUS='SCRATCH',FORM= & 'UNFORMATTED') RETURN C C FORMATS FOR PRINTED OUTPUT C 600 FORMAT ('0INPUT DATA FROM PROGRAM REFPK JOB:'/ &'0 TIME DATE TITLE'/' ',A,1X,A,2X,A/) 610 FORMAT ('1'/'0PROGRAM BGLP. ',A,1X,A,'. 'A/) 611 FORMAT ('0INPUT FILE: ',A) 612 FORMAT ('0OUTPUT FILE: ',A/) 629 FORMAT ('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 ('0NICOLET (NEE SYNTEX) P3 DIFFRACTOMETER'/) 616 FORMAT ('0ENRAF-NONIUS CAD4 DIFFRACTOMETER'/) 651 FORMAT ('0X-RAY DATA'/) 652 FORMAT ('0NEUTRON DATA'/) 628 FORMAT ('0 THETA(M) RHO(M) FRAC(DYN) SIG(FRAC)'/' ',2F10.2, &2E10.3/) 631 FORMAT ('0FILT. EDGE ALPHA-ONE ALPHA-BAR ALPHA-TWO'/' ',4F10.5/) 613 FORMAT ('0SIG(THETA) SIG(TIME) (THETA IN DEGREES AND TIME IN MICRO &SECONDS)'/' ',F10.4,F10.2/) 620 FORMAT ('0 TAU SIG(TAU) ATT SIG(ATT) (TAU IN MICROSE &CONDS)'/' ',2F10.2,2F10.2/) 622 FORMAT ('0BEGINNING AND ENDING REFLECTION SERIAL NUMBERS AND X-RAY & EXPOSURE TIMES (HOURS):'/'0 IIMIN IIMAX XTMIN XTM &AX'/' ',2I10,2F10.2/) 621 FORMAT ('0SCAN LIMITS (DECIMAL FRACTION OF SCAN WIDTH):'/ &'0 I1 I2'/' ',2F10.2/) 634 FORMAT ('0WINDOW WIDTH FOR INITIAL SCAN FOR PEAK LIMITS (DECIMAL F &RACTION OF SCAN WIDTH)'/'0 WINDOW'/' ',F10.4) 626 FORMAT ('0E.S.D. MULTIPLIER FOR STATISTICAL SIGNIFICANCE TESTS'/'0 & CUTOFF'/' ',F10.2) 627 FORMAT ('0NUMBER OF PASSES OF DATA-SMOOTHING FOR PEAK LOCATION'/'0 & NSMOOTH'/' ',I10) 623 FORMAT ('0PEAK POSITION MATRIX'/ &'0ORIENTATION MATRIX RECALCULATED BY LEAST-SQUARES FIT TO OMEGA'/ &' VALUES FOR THE INTENSITY-WEIGHTED PEAK CENTROIDS ALONG WITH'/ &' THE DIFFRACTOMETER ANGLES CHI AND PHI.'/ &'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) 617 FORMAT ('0PEAK WIDTH TENSORS'/ &'0W1 AND W2 ARE THE BASE WIDTHS OF THE HALF-PEAKS BELOW'/ &' THETA(ALPHA-1) AND ABOVE THETA(ALPHA-2), RESPECTIVELY,'/ &' IN UNITS OF DEGREES OF CRYSTAL ROTATION (THETA).'/ &'0LOW ANGLE HALF-PEAK W1') 618 FORMAT ('0HIGH ANGLE HALF-PEAK W2') 609 FORMAT ('0LORENTZIAN MODEL'/ &'0W = (Z(K)*Q(J,K)*Z(J))**(1/2) + T*TAN(THETA)') 608 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 = ',F3.1) 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) 642 FORMAT ('0 T = ',F8.4) 643 FORMAT ('0 SIG = ',F8.4) 647 FORMAT ('0FOR WEAK REFLECTIONS WITH NO SIGNIFICANT PEAK ABOVE BACK &GROUND, THE'/' DEFAULT ESTIMATE OF PEAK POSITION WILL BE'/'0X0 =', &' ',F4.2,' (DECIMAL FRACTION OF SCAN WIDTH)'/'0INSTEAD OF THE VALU &E CALCULATED USING THE DIFFRACTOMETER ORIENTATION MATRIX.'/) 648 FORMAT ('0FOR EACH REFLECTION, THE PEAK WIDTH WILL BE FIXED AS W = & L2 - L1, WHERE'/'0L1 = ',F4.2,' (DECIMAL FRACTION OF SCAN WIDTH)' &/' L2 = ',F4.2,', AND THUS'/'0W = ',F4.2,' CENTERED ABOUT THE INTE &NSITY-WEIGHTED PEAK CENTROID,'/'0INSTEAD OF THE VALUES CALCULATE', &'D USING THE ANISOTROPIC PEAK WIDTH TENSORS.'/) 644 FORMAT ('0CONDITIONS LIMITING POSSIBLE REFLECTIONS'/) 645 FORMAT (' ',2X,A) 646 FORMAT (' NONE') END C----------------------------------------------------------------------- SUBROUTINE LZ (TH,VT,FL,VL) C C LORENTZ FACTOR. EQUATORIAL DIFFRACTION GEOMETRY C U=COS(2*TH) V=SIN(2*TH) FL=1/V VL=VT*(-2*U/V**2)**2 RETURN END C----------------------------------------------------------------------- SUBROUTINE METRIC (A,GB) C C CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND C RECIPROCAL LATTICE METRIC TENSORS FROM DIRECT LATTICE C 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----------------------------------------------------------------------- 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 /BLOCK7/ 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 CHARACTER ATIME*8,ADATE*9,TITLE*80,HEADING*120 COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2, & XX0,XW1,XW2 PARAMETER (M1=4,M2=12,M3=50,MI=26,MJ=13,MK=16) DIMENSION HEADING(M1),JX(M2),JXMIN(M1,M2,M3),JXMAX(M1,M2,M3) DIMENSION XI(MI),YI(MI),NI(MI),XJ(MJ),YJ(MJ),NJ(MJ), & XK(MK),YK(MK),NK(MK) C C DATA FOR CATALOGING REFLECTIONS C DATA XI /-1000, -300, -100, -50, -20, -10, -5, -2, -1, 0, 1, 2, 5, & 10, 20, 50, 100, 300, 1000, 3000, 10000, 30000, 100000, & 300000, 1000000, 0/ DATA XJ /-10, -5, -2, -1, 0, 1, 2, 5, 10, 20, 50, 100, 0/ DATA XK /0.5, 0.65, 0.8, 0.9, 1.0, 1.05, 1.1, 1.15, 1.175, 1.2, & 1.25, 1.3, 1.35, 1.4, 1.45, 1.5/ DATA NTOT,NREF,NCEN,NSIG /4*0/ C C INITIALIZE ARRAYS FOR CATALOGING VARIOUS CLASSES OF DATA. C DO 11 I=1,MI YI(I)=0 NI(I)=0 11 CONTINUE DO 12 J=1,MJ YJ(J)=0 NJ(J)=0 12 CONTINUE DO 13 K=1,MK YK(K)=0 NK(K)=0 13 CONTINUE DO 15 J1=1,M1 DO 15 J2=1,M2 DO 15 J3=1,M3 JXMIN(J1,J2,J3)=+1E5 JXMAX(J1,J2,J3)=-1E5 15 CONTINUE WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,635) C C LOOP THROUGH REFLECTION DATA FILE. C REWIND IO3 1 READ (IO3,END=9) II,IH,IK,IL,S,R,SIGR,I0,IE,L1,L2,W,WEAK NTOT=NTOT+1 IF (II.LT.0) GO TO 1 NREF=NREF+1 IF (WEAK.EQ.0) NCEN=NCEN+1 IF (R.GE.CUTOFF*SIGR) NSIG=NSIG+1 C C PRINT SELECTED REFLECTION RECORDS. C CALL DPRINT (II,IH,IK,IL,IPRINT) IF (IPRINT.EQ.1) WRITE (ILP,636) II,IH,IK,IL,I0,MAX(IE,0),L1,L2, & R,SIGR C C CATALOG DATA IN CLASSES OF R = I/LP, R/SIG(R), AND S = C SIN(THETA)/LAMBDA. C DO 31 I=1,MI-1 IF (R.LT.XI(I)) GO TO 32 31 CONTINUE I=MI 32 NI(I)=NI(I)+1 YI(I)=YI(I)+R/SIGR DO 33 J=1,MJ-1 IF (R/SIGR.LT.XJ(J)) GO TO 34 33 CONTINUE J=MJ 34 NJ(J)=NJ(J)+1 YJ(J)=YJ(J)+R/SIGR IF (R/SIGR.LT.CUTOFF) GO TO 1 DO 35 K=1,MK-1 IF (S.LT.XK(K)) GO TO 36 35 CONTINUE K=MK 36 NK(K)=NK(K)+1 YK(K)=YK(K)+R/SIGR C C CATALOG PROFILES WITH EXTREME FEATURES. C JX(1)=II JX(2)=IH JX(3)=IK JX(4)=IL JX(5)=1000*S JX(6)=1000*R JX(7)=1000*SIGR JX(8)=I0 JX(9)=MAX(IE,0) JX(10)=L1 JX(11)=L2 JX(12)=1000*W 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) DO 95 J1=1,M1 DO 95 J2=1,M2 DO 95 J3=1,M3 IF (JXMIN(J1,J2,J3).EQ.+1E5) JXMIN(J1,J2,J3)=0 IF (JXMAX(J1,J2,J3).EQ.-1E5) JXMAX(J1,J2,J3)=0 95 CONTINUE C C PRINT CATALOG OF CLASSES OF DATA. C WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,641) NTOT,NREF,NCEN,NSIG,CUTOFF ISUM=0 JSUM=0 KSUM=0 DO 37 I=1,MI IF (NI(I).EQ.0) GO TO 37 ISUM=ISUM+NI(I) YI(I)=YI(I)/NI(I) 37 CONTINUE DO 38 J=1,MJ IF (NJ(J).EQ.0) GO TO 38 JSUM=JSUM+NJ(J) YJ(J)=YJ(J)/NJ(J) 38 CONTINUE DO 39 K=1,MK IF (NK(K).EQ.0) GO TO 39 KSUM=KSUM+NK(K) YK(K)=YK(K)/NK(K) 39 CONTINUE WRITE (ILP,637) XI(1),NI(1),YI(1),(XI(I-1),XI(I),NI(I), &YI(I),I=2,(MI-1)),XI(MI-1),NI(MI),YI(MI) WRITE (ILP,640) ISUM WRITE (ILP,638) XJ(1),NJ(1),YJ(1),(XJ(J-1),XJ(J),NJ(J), &YJ(J),J=2,(MJ-1)),XJ(MJ-1),NJ(MJ),YJ(MJ) WRITE (ILP,640) JSUM WRITE (ILP,639) XK(1),NK(1),YK(1),(XK(K-1),XK(K),NK(K), &YK(K),K=2,(MK-1)),XK(MK-1),NK(MK),YK(MK) WRITE (ILP,640) KSUM C C PRINT CATALOGS OF EXTREME CASES. C WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,644) DO 49 J1=1,M1 WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,645) HEADING(J1) DO 59 J3=1,MIN(M3,NSIG) F5=JXMIN(J1,5,J3)*0.001 F6=JXMIN(J1,6,J3)*0.001 F7=JXMIN(J1,7,J3)*0.001 F12=JXMIN(J1,12,J3)*0.001 59 WRITE (ILP,646) (JXMIN(J1,J2,J3),J2=1,4),F5,F6,F7, &(JXMIN(J1,J2,J3),J2=8,11),F12 WRITE (ILP,610) ATIME,ADATE,TITLE WRITE (ILP,645) HEADING(J1) DO 69 J3=1,MIN(M3,NSIG) F5=JXMAX(J1,5,J3)*0.001 F6=JXMAX(J1,6,J3)*0.001 F7=JXMAX(J1,7,J3)*0.001 F12=JXMAX(J1,12,J3)*0.001 69 WRITE (ILP,646) (JXMAX(J1,J2,J3),J2=1,4),F5,F6,F7, &(JXMAX(J1,J2,J3),J2=8,11),F12 49 CONTINUE RETURN C C FORMAT STATEMENTS FOR PRINTED OUTPUT C 610 FORMAT ('1'/'0PROGRAM BGLP. ',A,1X,A,'. 'A/) 635 FORMAT ('0 I H K L X0 XE L1 L2 Y=I/LP', &' SIG(Y)'/) 636 FORMAT (' ',I6,2X,3I4,2X,2(I4,2X),2I4,2F10.2) 637 FORMAT ('0',12X'Y-INTERVAL NUMBER AVERAGE Y/SIG(Y), ' & 'WHERE Y = I/LP'/ & '0',9X, ' Y < ',F9.1,5X,I7,F9.1/ & 24(' ',F9.1,' < Y < ',F9.1,5X,I7,F9.1/), & ' ',F9.1,' < Y ',9X, 5X,I7,F9.1) 640 FORMAT (' ',32X,'TOTAL'/ ' ',30X,I7) 638 FORMAT ('0',12X,'R-INTERVAL NUMBER AVERAGE R, WHERE ' & 'R = Y/SIG(Y)'/ & '0',9X, ' R < ',F9.1,5X,I7,F9.1/ & 11(' ',F9.1,' < R < ',F9.1,5X,I7,F9.1/), & ' ',F9.1,' < R ',9X, 5X,I7,F9.1) 639 FORMAT ('1',12X,'S-INTERVAL NUMBER AVERAGE R, WHERE ' & 'S = SIN(THETA)/LAMBDA'/ & '0',9X, ' S < ',F9.3,5X,I7,F9.1/ & 14(' ',F9.3,' < S < ',F9.3,5X,I7,F9.1/), & ' ',F9.3,' < S ',9X, 5X,I7,F9.1) 641 FORMAT ('0NTOT = ',I6,' TOTAL REFLECTIONS (EXCLUDING SYMMETRY FOR &BIDDEN REFLECTIONS)'/ & ' NREF = ',I6,' NOT STANDARD REFERENCE REFLECTIONS'/ & ' NCEN = ',I6,' WITH A SIGNIFICANT INTENSITY-WEIGHTED PEA &K CENTROID'/ & ' NSIG = ',I6,' WITH Y > CUTOFF*SIGMA(Y), WHERE CUTOFF = & ',F5.2) 644 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 / &'S=SIN(THETA)/LAMBDA', &'Y=I/LP, LORENTZ AND POLARIZATION CORRECTED NET INTENSITY', &'X0, INTENSITY WEIGHTED PEAK CENTROID', &'W=(L2-L1+1)*DTH, BASE WIDTH OF REFLECTION (DEGREES THETA)'/ 645 FORMAT ('0',A/'0 I H K L S Y=I/LP SIG(Y)', &' X0 XE L1 L2 W'/) 646 FORMAT (' ',I6,2X,3I4,2XF5.3,2F10.2,2XI4,2XI4,2X,2I4,F10.2) END C----------------------------------------------------------------------- SUBROUTINE PROCESS PARAMETER (NMAX=100) CHARACTER ATIME*8,ADATE*9,TITLE*80 COMMON /BLOCK0/ IO1,IO2,IO3,ILP COMMON /BLOCKC/ NCOND,ICOND(38) COMMON /BLOCK1/ ATIME,ADATE,TITLE COMMON /BLOCK2/ GINV(3,3),WLE,WL1,WL0,WL2 COMMON /BLOCK3/ WINDOW,CUTOFF,NSMOOTH,II1,II2,XT1,XT2,SC1,SC2, & XX0,XW1,XW2 COMMON /BLOCK5/ IDIFF,U(3,3),MODEL1,MODEL2,Q1(3,3),Q2(3,3),T1,T2, & F,SIGU(3,3),SIGQ1(3,3),SIGQ2(3,3),SIGT1,SIGT2,TANTHM COMMON /BLOCK6/ VARTHE,VARTIM,TAU,VARTAU,ATT,VARATT COMMON /BLOCK7/ NSTEP,XX(NMAX),YY(NMAX),VX(NMAX),VY(NMAX) COMMON /BLOCK8/ TH0,X0,SIGX0,XE,X1,X2,W1,W2,I1,I2,L1,L2,R,SIGR,EPS DIMENSION H(3),ANGLES(4),YYTEMP(NMAX),VYTEMP(NMAX) PI=ACOS(-1.0) DEG=180/PI RAD=PI/180 C C LOOP THROUGH PROCESSING 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 SKIP SYMMETRY FORBIDDEN REFLECTIONS. C CALL HKLTEST (IH,IK,IL,NCOND,ICOND,IABSENT) IF (IABSENT.NE.0) GO TO 1 C C STORE INDICES IN 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 TRANSFORM TO SETTING ANGLES AS DEFINED BY WALTER HAMILTON, C INTERNATIONAL TABLES FOR X-RAY CRYSTALLOGRAPHY, VOL. IV, 1974, PP. C 273-284. C CALL GEOM (IDIFF,ANGLES,TWOTH,OMEGA,CHI,PHI) C C SET 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 SMOOTH PROFILE FOR PEAK CENTROID LOCATION. C DO 11 I=1,NSTEP YYTEMP(I)=YY(I) VYTEMP(I)=VY(I) 11 CONTINUE IF (NSMOOTH.NE.0) THEN CALL SMOOTH (NSMOOTH,NSTEP,YYTEMP) CALL SMOOTH (NSMOOTH,NSTEP,VYTEMP) END IF C C FIND INTENSITY-WEIGHTED CENTROID OF PEAK ABOVE BACKGROUND. C CALL CENTROID (WINDOW,CUTOFF,XX,YYTEMP,VX,VYTEMP,I1,I2,X0,SIGX0) C C TREAT SCANS THAT HAVE NO SIGNIFICANT PEAK ABOVE BACKGROUND. C IF (X0.EQ.90) THEN WEAK=1 CALL XCALC (XX0,XX,I1,I2,H,U,SIGU,GINV,TWOTH,OMEGA,CHI,PHI, & X0,SIGX0) ELSE WEAK=0 END IF C C CALCULATE POSITIONS OF SPECTRAL FEATURES. C XE=X0+(THE-TH0) X1=X0+(TH1-TH0) X2=X0+(TH2-TH0) IF (WEAK.EQ.0) CALL XCHECK (XE,X1,X0,X2,XX,YYTEMP,I1,I2) C C CALCULATE PEAK WIDTHS. C TANTH=TAN(TH0*RAD) CALL WCALC (XW1,XX,I1,I2,CHI,PHI,MODEL1,Q1,T1,F,TANTH,TANTHM, & SIGQ1,SIGT1,W1,SIGW1) CALL WCALC (XW2,XX,I1,I2,CHI,PHI,MODEL2,Q2,T2,F,TANTH,TANTHM, & SIGQ2,SIGT2,W2,SIGW2) C C LIMIT THE ESTIMATED UNCERTAINTIES IN PEAK POSITION AND PEAK WIDTHS C TO ONE STEP. C DX=(XX(I2)-XX(I1))/(I2-I1+1) SIGX0=MIN(SIGX0,DX) SIGW1=MIN(SIGW1,DX) SIGW2=MIN(SIGW2,DX) C C CALCULATE PEAK LIMITS. C L1=NINDEX(X1-SIGX0-W1-SIGW1) L2=NINDEX(X2+SIGX0+W2+SIGW2) C C PERFORM BACKGROUND SUBTRACTION AND LORENTZ AND POLARIZATION C CORRECTIONS. C IF (WLE.EQ.0.OR.WEAK.EQ.1) THEN CALL BGLP1 ELSE C C EPSILON IS AN ESTIMATE OF THE HALF-WIDTH OF THE BETA-FILTER C ABSORPTION EDGE BROADENED DUE TO FINITE INSTRUMENTAL RESOLUTION. C EPS=0.5*(W1+W2) J1=NINDEX(XE-EPS) J2=NINDEX(XE+EPS) IF (J2.LT.L1-1) THEN I1=MAX(I1,J2) CALL BGLP1 ELSE CALL BGLP2 END IF END IF C C WRITE REFLECTION RECORD TO OUTPUT FILES. C CALL WRITE1 (IO2,II,IH,IK,IL,ANGLES,R,SIGR,XTIME) I0=NINDEX(X0) IE=NINDEX(XE) IF (WLE.EQ.0) IE=0 W=XX(L2)-XX(L1) WRITE (IO3) II,IH,IK,IL,S,R,SIGR,I0,IE,L1,L2,W,WEAK C C LOOP BACK TO PROCESS NEXT REFLECTION. C GO TO 1 C C END LOOP THROUGH REFLECTION DATA. C 9 CLOSE (UNIT=IO2) ENDFILE IO3 RETURN END C----------------------------------------------------------------------- SUBROUTINE PZ (TH,VT,FP,VP) C C POLARIZATION FACTOR. WITH OR WITHOUT INCIDENT BEAM CRYSTAL C MONOCHROMATOR C COMMON /BLOCK4/ C0K,C1K,C2K,C0D,C1D,C2D,FD,VF,INEUTRON C C C0K = COS(2*THM)**2 KINEMATIC (MOSAIC) MONOCHROMATOR CRYSTAL C C0D = ABS(COS(2*THM)) DYNAMIC (PERFECT) MONOCHROMATOR CRYSTAL C C C1 = COS(RHO)**2 + SIN(RHO)**2*C0 C C2 = SIN(RHO)**2 + COS(RHO)**2*C0 C C RHO = 0 PARALLEL MONOCHROMATOR-DIFFRACTOMETER GEOMETRY C RHO = 90 DEG PERPENDICULAR GEOMETRY C C FD = 0 ALL KINEMATIC, NO DYNAMIC MONOCHROMATOR DIFFRACTION C FD = 1 ALL DYNAMIC, NO KINEMATIC MONOCHROMATOR DIFFRACTION C C VF = VAR(FD), THE VARIANCE OF THE DYNAMIC DIFFRACTION FRACTION C C FOR AN UNPOLARIZED INCIDENT BEAM (I.E., NO MONOCHROMATOR), C THM = 0, C0 = C1 = C2 = 1, AND FD = VF = 0. C C FOR NEUTRON DIFFRACTION DATA, INEUTRON = 1. C IF (INEUTRON.EQ.1) THEN FP=1 VP=0 ELSE U=COS(2*TH) V=SIN(2*TH) PK=(C1K+C2K*U**2)/(1+C0K) VK=VT*(-4*C2K*U*V/(1+C0K))**2 IF (FD.EQ.0) THEN FP=PK VP=VK ELSE PD=(C1D+C2D*U**2)/(1+C0D) VD=VT*(-4*C2D*U*V/(1+C0D))**2 FP=FD*PD+(1-FD)*PK VP=FD**2*VD+(1-FD)**2*VK+(PD**2+PK**2)*VF END IF END IF 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 WCALC (FW,X,I1,I2,CHI,PHI,MODEL,Q,T,F,TANTH,TANTHM, & SIGQ,SIGT,W,SIGW) C C CALCULATED PEAK WIDTHS C DIMENSION X(1),Q(3,3),SIGQ(3,3),Z(3) DATA RAD /0.017453293/ IF (FW.NE.0) THEN W=FW*(X(I2)-X(I1)) SIGW=0 RETURN END IF C C GET THE COMPONENTS ALONG CARTESIAN AXES FIXED IN THE CRYSTAL OF A C UNIT VECTOR, Z, NORMAL TO THE EQUATORIAL PLANE OF THE INCIDENT AND C DIFFRACTED BEAM VECTORS. C COSCH=COS(CHI*RAD) SINCH=SIN(CHI*RAD) COSPH=COS(PHI*RAD) SINPH=SIN(PHI*RAD) Z(1)=SINPH*SINCH Z(2)=COSPH*SINCH Z(3)=COSCH C C CALCULATE PEAK WIDTHS. C U=0 V=0 DO 10 I=1,3 DO 10 J=1,3 U=U+Z(I)*Z(J)*Q(I,J) V=V+(Z(I)*Z(J)*SIGQ(I,J))**2 10 CONTINUE IF (MODEL.EQ.1) THEN C C LORENTZIAN CONVOLUTION MODEL C W=SQRT(U)+T*TANTH SIGW=SQRT(V/(2*U)**2+(SIGT*TANTH)**2) ELSE IF (MODEL.EQ.2) THEN C C GAUSSIAN CONVOLUTION MODEL C W=SQRT(U+(T*TANTH)**2) SIGW=SQRT(V+(2*T*TANTH**2*SIGT)**2)/(2*W) ELSE IF (MODEL.EQ.3) THEN C C GAUSSIAN CONVOLUTION MODEL FOR PARALLEL MONOCHROMATOR GEOMETRY C W=SQRT(ABS(U-2*F*T*TANTH/TANTHM+T*(TANTH/TANTHM)**2)) SIGW=SQRT(V/(2*U)**2+(2*F*SIGT*TANTH/TANTHM)**2 & +(SIGT*(TANTH/TANTHM)**2)**2) END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE WRITE1 (IFILE,II,IH,IK,IL,ANGLES,R,SIGR,XTIME) DIMENSION ANGLES(4) WRITE (IFILE) II,IH,IK,IL,ANGLES,R,SIGR,XTIME RETURN END C----------------------------------------------------------------------- SUBROUTINE XCALC (FX,X,I1,I2,H,U,SIGU,GINV,TWOTH,OMEGA,CHI,PHI, & X0,SIGX0) C C EXPECTED PEAK POSITION C DIMENSION X(1),H(3),U(3,3),SIGU(3,3),GINV(3,3),Y(3),V(3) DATA RAD,DEG /0.017453293, 57.2957795/ IF (FX.NE.0) THEN X0=X(I1)+FX*(X(I2)-X(I1)) SIGX0=(X(I2)-X(I1))/(I2-I1+1) RETURN END IF C C GET THE COMPONENTS ALONG CARTESIAN AXES FIXED IN THE CRYSTAL OF A C UNIT VECTOR, Y, PARALLEL TO THE RECIPROCAL LATTICE VECTOR, DSTAR, C IN THE DIFFRACTING CONDITION. C DO 1 I=1,3 Y(I)=0 V(I)=0 DO 1 J=1,3 Y(I)=Y(I)+H(J)*U(J,I) V(I)=V(I)+(H(J)*SIGU(J,I))**2 1 CONTINUE DSTAR=2*SINTHL(H,GINV) DO 2 I=1,3 Y(I)=Y(I)/DSTAR V(I)=V(I)/DSTAR**2 2 CONTINUE IF (TWOTH.LT.0) THEN DO 3 I=1,3 Y(I)=-Y(I) 3 CONTINUE END IF C C CALCULATE EXPECTED PEAK DISPLACEMENT, DELTA(OMEGA). C 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 C C THE POSITIVE SENSE OF THE OMEGA ROTATION HAS BEEN DEFINED SUCH C THAT DELTA(X0) = -DELTA(OMEGA). C X0=0.5*(X(I1)+X(I2))-DELTA VARSIN=V(1)*COSPH**2+V(2)*SINPH**2 VARCOS=(V(1)*SINPH**2+V(2)*COSPH**2)*COSCH**2+V(3)*SINCH**2 OMEGA=ARCTAN(SINOM,COSOM) VAROM=(COS(OMEGA)**2/COSOM)**2*(TAN(OMEGA)**2*VARCOS+VARSIN) SIGX0=SQRT(VAROM)*DEG RETURN END C----------------------------------------------------------------------- SUBROUTINE XCHECK (XE,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 XE=XE+D 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 /BLOCK6/ 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-----------------------------------------------------------------------