PROGRAM EVAL 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, EXTENSION 335 C ELECTRONIC MAIL: Blessing@HWI.Buffalo.Edu C C OCTOBER 1995,..., JANUARY 2001 C CHARACTER ATIME*8,ADATE*9,BTIME*8,BDATE*9,TITLE*80,FILE*80,FORM*11 CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5,XRAY*2,EL*2 DIMENSION CELL(6),ASTAR(6),GINV(3,3),B(3,3),C(3,3),COV(13,13), & U(3,3),V(3,3),AIAJ(6),VAR(13),H(3) DIMENSION EL(20),XI(20),ZI(20),A1(20),B1(20),A2(20),B2(20),A3(20), & B3(20),A4(20),B4(20),C0(20),FP(20),FPP(20) DATA B,C,COV /9*0,9*0,169*0/ DATA NLIST,CUTOFF,EMIN,EMAX /100,3,1.5,0.5/ DATA XMIN,XMAX,JMIN,JMAX,KMIN,KMAX,LMIN,LMAX /9,0,+99,-99,+99,-99, & +99,-99/ C C FORTRAN-90 DYNAMIC STORAGE ALLOCATION C REAL, ALLOCATABLE::DATA(:) INTEGER, ALLOCATABLE::INDEX(:) C C FORTRAN-90 TIME AND DATE TRANSLATIONS C CALL VTIME (ATIME) CALL VDATE (ADATE) c ATIME='hh:mm:ss' c ADATE='dd-mmm-yy' IO1=10 IO2=20 ILP=60 C C READ CONTROL DATA AND PARAMETERS FROM PROGRAM ROGERS OR PROGRAM LEVY. C OPEN (UNIT=IO1,FILE='eval.dat',STATUS='OLD') READ (IO1,999) BTIME READ (IO1,999) BDATE READ (IO1,999) TITLE READ (IO1,997) ITYPE READ (IO1,999) FILE READ (IO1,999) LATT READ (IO1,999) PTGP READ (IO1,998) CELL READ (IO1,999) XRAY READ (IO1,997) NEL DO 5 I=1,NEL READ (IO1,996) EL(I),XI(I),ZI(I),A1(I),B1(I),A2(I),B2(I) READ (IO1,995) A3(I),B3(I),A4(I),B4(I),C0(I),FP(I),FPP(I) 5 CONTINUE READ (IO1,998) SMIN,SMAX READ (IO1,993) ASOLV,BSOLV READ (IO1,994) SCALEK READ (IO1,994) B(1,1),B(1,2),B(1,3),B(2,2),B(2,3),B(3,3) READ (IO1,994,END=9) C(1,1),C(1,2),C(1,3),C(2,2),C(2,3),C(3,3) READ (IO1,993,END=9) ((COV(I,J),J=1,13),I=1,13) READ (IO1,992,END=9) NLIST,CUTOFF,EMAX,EMIN 999 FORMAT (1X,A) 998 FORMAT (1X,6F10.0) 997 FORMAT (1X,I6) 996 FORMAT (1X,A2,F10.2,F6.0,4(1X,F10.6)) 995 FORMAT (1X,5(1X,F10.6),2X,2(1X,F10.6)) 994 FORMAT (1X,6E12.5) 993 FORMAT (1X,7E10.3) 992 FORMAT (I10,3F10.0) 9 CLOSE (UNIT=IO1,STATUS='KEEP') CALL DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP) CALL METRIC (CELL,ASTAR,GINV) PI=3.14159 IF (B(2,2).EQ.0) THEN BISO=B(1,1) CISO=C(1,1) UISO=BISO/(8*PI**2) VISO=CISO/(8*PI**2) ELSE DO 6 I=1,3 DO 6 J=I,3 U(I,J)=B(I,J)/(2*PI**2*ASTAR(I)*ASTAR(J)) V(I,J)=C(I,J)/(2*PI**2*ASTAR(I)*ASTAR(J)) B(J,I)=B(I,J) C(J,I)=C(I,J) U(J,I)=U(I,J) V(J,I)=V(I,J) 6 CONTINUE END IF C C PRINT CONTROL DATA. C OPEN (UNIT=ILP,FILE='eval.lp',STATUS='NEW') WRITE (ILP,899) ATIME,ADATE,TITLE WRITE (ILP,898) BTIME,BDATE WRITE (ILP,897) FILE WRITE (ILP,896) SYST,LATT,LAUE,PTGP IF (PTGP.EQ.LAUE) THEN WRITE (ILP,8961) ELSE WRITE (ILP,8960) END IF WRITE (ILP,893) CELL WRITE (ILP,892) (EL(I),XI(I),I=1,NEL) WRITE (ILP,891) XRAY,(EL(I),FP(I),FPP(I),I=1,NEL) IF (SMIN.GT.0.OR.SMAX.LT.9) WRITE (ILP,895) SMIN,SMAX IF (ASOLV.NE.0) WRITE (ILP,613) ASOLV,BSOLV WRITE (ILP,899) ATIME,ADATE,TITLE IF (B(2,2).EQ.0) THEN WRITE (ILP,860) SCALEK,BISO,CISO DO 3 I=1,3 VAR(I)=COV(I,I) IF (VAR(I).EQ.0) VAR(I)=1 3 CONTINUE WRITE (ILP,861) (SQRT(COV(I,I)),I=1,3),(I,(COV(I,J)/SQRT(VAR(I)* & VAR(J)),J=1,3),I=1,3) ELSE WRITE (ILP,889) SCALEK, & B(1,1),B(1,2),B(1,3),B(2,2),B(2,3),B(3,3), & C(1,1),C(1,2),C(1,3),C(2,2),C(2,3),C(3,3), & U(1,1),U(1,2),U(1,3),U(2,2),U(2,3),U(3,3), & V(1,1),V(1,2),V(1,3),V(2,2),V(2,3),V(3,3) AIAJ(1)=ASTAR(1)**2 AIAJ(2)=ASTAR(1)*ASTAR(2) AIAJ(3)=ASTAR(1)*ASTAR(3) AIAJ(4)=ASTAR(2)**2 AIAJ(5)=ASTAR(2)*ASTAR(3) AIAJ(6)=ASTAR(3)**2 DO 13 I=1,13 VAR(I)=COV(I,I) IF (VAR(I).EQ.0) VAR(I)=1 13 CONTINUE WRITE (ILP,899) ATIME,ADATE,TITLE WRITE (ILP,888) (SQRT(COV(1+I,1+I))/(2*PI**2*AIAJ(I)),I=1,6), & (SQRT(COV(7+I,7+I))/(2*PI**2*AIAJ(I)),I=1,6), & (SQRT(COV(I,I)),I=1,13), & (I,(NINT(100*COV(I,J)/SQRT(VAR(I)*VAR(J))),J=1,13),I=1,13) END IF 899 FORMAT (/1H1,130('-')/1X,'PROGRAM EVAL. ',A,', ',A,'. ',A) 898 FORMAT (/1X, &'SCALE AND THERMAL PARAMETERS FROM PROGRAM ROGERS OR PROGRAM LEV', &'Y JOB:'/1X,A,', ',A) 897 FORMAT (/1X,'INPUT REFLECTION FILE: ',A) 896 FORMAT (/1X,A//1X,A,'-LATTICE'//1X, &'LAUE GROUP ',A//1X, &'POINT GROUP ',A) 8960 FORMAT (/1X,'NONCENTROSYMMETRIC STRUCTURE') 8961 FORMAT (/1X,'CENTROSYMMETRIC STRUCTURE') 893 FORMAT (/1X, &'LATTICE PARAMETERS'/1X, &' A B C ALPHA BETA GAMMA'/1X, & 3F10.4,3F10.3) 892 FORMAT (/1X,'UNIT CELL CONTENTS:'//(1X,A2,F10.2)) 891 FORMAT (/1X,'RADIATION: ',A2//1X,'ANOMALOUS DISPERSION CORRECTION & TERMS:'/1X,'--------------------------------------'/1X, &'EL FP FPP'/1X,'-- -- ---'/ &(1X,A2,2F10.3)) 895 FORMAT (/1X, &'E-VALUES ARE CALCULATED FOR ALL DATA, BUT'/1X, &'E-STATISTICS ARE COMPILED ONLY FOR THE DATA WITH SMIN .LE. SIN(', &'THETA)/LAMBDA .LE. SMAX,'/1X, &'SMIN = ',F6.3/1X,'SMAX = ',F6.3,' RECIPROCAL ANGSTROMS,'/1X, &'WHICH WERE USED TO ESTIMATE SCALEK AND B(I,J).') 613 FORMAT (/1X, &'BULK SOLVENT COMPENSATION PARAMETERS:'/1X, &'-------------------------------------'/1X, &' E(PROT) = E(MEAS)/[1 - ASOLV*EXP{-BSOLV*[SIN(THETA)/LAMBDA]**', &'2}'/1X, &' ASOLV = ',E10.3/1X, &' BSOLV = ',E10.3' A**2 (SQUARED ANGSTROMS)') 860 FORMAT (/1X, &'SCALEK, BISO, AND CISO ARE DEFINED ACCORDING TO:'//1X, &' F(ABSOLUTE) = SCALEK*F(RELATIVE)'/1X, &' f(T) = f(0)*, WHERE S = SIN(THETA)/LAMBDA'/1X, &' = EXP(-BISO*S**2 + CISO**2*S**4)'//1X, &'PARAMETERS:'//1X, &' SCALEK BISO CISO'/1X,3E12.4) 861 FORMAT (/1X, &'PARAMETER E.S.D.''S:'//1X, &' SIGMA(K) SIGMA(B) SIGMA(C)'/1X,3E12.2//1X, &'CORRELATION MATRIX:'//1X, &' K BISO CISO'//1X, &' 1 2 3'//(1X,I2,3F6.2)) 889 FORMAT (/1X, &'SCALEK, B(I,J), C(I,J), U(I,J), AND B(I,J) ARE DEFINED ACCORDIN', &'G TO:'//1X, &' F(ABSOLUTE) = SCALEK*F(RELATIVE)'/1X, &' f(T) = f(0)*EXP(-SUM(I)SUM(J) H(I)*H(J)*B(I,J))'/1X, &' f(T) = f(0)*EXP(-2*PI**2*SUM(I)SUM(J) H(I)*H(J)*ASTAR(I)*ASTA', &'R(J)*U(I,J))'/1X, &' = EXP(-2*(H(J)**H(I) - (H(', &'J)**H(I))**2))'//1X, &'SCALEK'/1X,E12.4//1X, &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENT PARAMETERS:'/1X, &' B(1,1) B(1,2) B(1,3) B(2,2) B(2,3) ', &' B(3,3)'/1X,6E12.4,' (DIMENSIONLESS)'//1X, &'B-DISTRIBUTION DISPERSION PARAMETERS:'/1X, &' C(1,1) C(1,2) C(1,3) C(2,2) C(2,3) ', &' C(3,3)'/' ',6E12.4,' (DIMENSIONLESS)'//1X, &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENTS:'/1X, &' U(1,1) U(1,2) U(1,3) U(2,2) U(2,3) ', &' U(3,3)'/1X,6E12.4,' SQUARED ANGSTROMS'//1X, &'U-DISTRIBUTION DISPERSION PARAMETERS:'/1X, &' V(1,1) V(1,2) V(1,3) V(2,2) V(2,3) ', &' V(3,3)'/1X,6E12.4,' SQUARED ANGSTROMS') 888 FORMAT (/1X, &'PARAMETER E.S.D.''S'//1X, &'NOTE THAT THE ORDER OF THE MEAN-SQUARE DISPLACEMENT PARAMETER E', &'RRORS IS'/1X, &'11, 12, 13, 22, 23, 33 NOT 11, 22, 33, 12, 13, 23.'/1X, &' ---'//1X, &' SIGMA(U11) SIGMA(U12) SIGMA(U13) SIGMA(U22) SIGMA(U23) S', &'IGMA(U33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X, &' SIGMA(V11) SIGMA(V12) SIGMA(V13) SIGMA(V22) SIGMA(V23) S', &'IGMA(V33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X, &' SIGMA(K)'/1X,E12.2//1X &' SIGMA(B11) SIGMA(B12) SIGMA(B13) SIGMA(B22) SIGMA(B13) S', &'IGMA(B33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X &' SIGMA(C11) SIGMA(C12) SIGMA(C13) SIGMA(C12) SIGMA(C23) S', &'IGMA(C33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X, &'CORRELATION MATRIX 100*CORR(I,J):'//1X, &' K B11 B12 B13 B22 B23 B33 C11 C12 C13 C22 C23 C33'//1X, &' 1 2 3 4 5 6 7 8 9 10 11 12 13'/ &/(1X,I2,2X,13I4)) C C CALCULATE INITIAL E-VALUES. C IF (ITYPE.LT.3) THEN FORM='FORMATTED' ELSE FORM='UNFORMATTED' END IF OPEN (UNIT=IO1,FILE=FILE,STATUS='OLD',FORM=FORM) OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='FORMATTED') SUMESQ=0 SUMM=0 IEND=0 10 CALL READ1 (ITYPE,IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF) IF (IEND.NE.0) GO TO 19 S=SINTHL(J,K,L,GINV) XMIN=MIN(XMIN,S) XMAX=MAX(XMAX,S) JMIN=MIN(JMIN,J) JMAX=MAX(JMAX,J) KMIN=MIN(KMIN,K) KMAX=MAX(KMAX,K) LMIN=MIN(LMIN,L) LMAX=MAX(LMAX,L) CALL FCALC (NEL,XI,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S,SUMFSQ) IF (B(2,2).EQ.0) THEN HTBH=BISO*S**2 HTCH=CISO*S**2 ELSE H(1)=J H(2)=K H(3)=L HTBH=VTMV(H,B) HTCH=VTMV(H,C) END IF Q=(1/SCALEK)*EXP(-HTBH+HTCH**2) VARQ=VARF(SCALEK,B,C,COV,S,H,HTBH,HTCH,Q) CALL ESYM (IPTGP,J,K,L,EPSILON,M) EPSILON=MLATT*EPSILON E=F/(Q*SQRT(EPSILON*SUMFSQ)) SIGE=E*SQRT((SIGF/F)**2+VARQ/Q**2) F=F*SCALEK SIGF=SIGF*SCALEK FSQ=FSQ*SCALEK**2 SIGFSQ=SIGFSQ*SCALEK**2 IF (ASOLV.NE.0) THEN Q=1/(1-ASOLV*EXP(-BSOLV*S**2)) E=E*Q SIGE=SIGE*Q END IF CALL WRITE1 (IO2,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) IF (SMIN.LE.S.AND.S.LE.SMAX) THEN CALL CSYM (IPTGP,J,K,L,IC) IF (IC.EQ.0) M=2*M SUMESQ=SUMESQ+M*E**2 SUMM=SUMM+M END IF GO TO 10 19 CONTINUE AVESQ=SUMESQ/SUMM RESCL=1/SQRT(AVESQ) ENDFILE IO2 CLOSE (UNIT=IO1,STATUS='KEEP') C C RE-SCALE E VALUES SO THAT = 1. C WRITE (ILP,899) ATIME,ADATE,TITLE WRITE (ILP,878) XMIN,XMAX,JMIN,JMAX,KMIN,KMAX,LMIN,LMAX,XMAX, &1/(2*XMAX) 878 FORMAT (/1X, &'SIN(THETA)/LAMBDA AND INDEX LIMITS:'//1X, &F7.4,' .LE. S .LE. ',F7.4,' ANGSTROM**(-1)'//1X, &I4,' .LE. H .LE. ',I4/1X, &I4,' .LE. K .LE. ',I4/1X, &I4,' .LE. L .LE. ',I4//1X, &'DATA RESOLUTION:'//1X, &'SMAX = ',F6.3,' ANGSTROM**(-1)'/1X,'DMIN = ',F5.2,' ANGSTROM') WRITE (ILP,880) AVESQ,RESCL WRITE (ILP,887) 880 FORMAT (/1X, &'MULTIPLICITY-WEIGHTED AVERAGE E**2 ',F12.4//1X, &'OVERALL RESCALING FACTOR FOR E-VALUES',F12.4) 887 FORMAT (/1X, &' H K L SIN(TH)/L D(HKL) FSQ SIGFSQ ', &' F SIGF E SIGE'/1X, &' - - - --------- ------ --- ------ ', &' - ---- - ----'/) 886 FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3) OPEN (UNIT=IO1,FILE='edata.hkl',STATUS='NEW',FORM='FORMATTED') REWIND IO2 IEND=0 NDATA=0 20 CALL READ2 (IO2,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) IF (IEND.NE.0) GO TO 29 NDATA=NDATA+1 E=E*RESCL SIGE=SIGE*RESCL CALL WRITE1 (IO1,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) CALL DPRINT (J,K,L,IPRINT) IF (IPRINT.NE.0) THEN S=SINTHL(J,K,L,GINV) D=1/(2*S) WRITE (ILP,886) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE END IF GO TO 20 29 CLOSE (UNIT=IO2,STATUS='DELETE') ENDFILE IO1 WRITE (ILP,881) NDATA 881 FORMAT (/1X,'NDATA = ',I6,' UNIQUE REFLECTION DATA') C C COMPILE E-STATISTICS. C WRITE (ILP,899) ATIME,ADATE,TITLE WRITE (ILP,800) SMIN,MIN(SMAX,XMAX) 800 FORMAT (/1X, &'E-STATISTICS FOR DATA WITH SMIN .LE. SIN(THETA)/LAMBDA .LE. SMA', &'X, WHERE'/1X, &'SMIN = ',F6.3,' RECIPROCAL ANGSTROMS'/1X, &'SMAX = ',F6.3) CALL ESTATS (IO1,ILP,SMIN,MIN(SMAX,XMAX),IPTGP,MLATT,GINV) C C ALLOCATE STORAGE FOR SORTING ARRAYS. C ALLOCATE (DATA(NDATA),INDEX(NDATA)) C C WRITE d(hkl)-SORTED OUTPUT FILE. C REWIND IO1 OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='UNFORMATTED', & ACCESS='DIRECT',RECL=36) DO I=1,NDATA CALL READ2 (IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) WRITE (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE DATA(I)=SINTHL(J,K,L,GINV) END DO CALL SORT (NDATA,DATA,INDEX) CLOSE (UNIT=IO1,STATUS='KEEP') OPEN (UNIT=IO1,FILE='edata.ddd',STATUS='NEW',FORM='FORMATTED') DO I=1,NDATA READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE CALL WRITE1 (IO1, J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) END DO C C LIST LOWEST RESOLUTION E-MAGNITUDES. C WRITE (ILP,899) ATIME,ADATE,TITLE WRITE (ILP,'(/1X,''LOWEST RESOLUTION E-MAGNITUDES:'')') WRITE (ILP,877) 877 FORMAT (/1X, &' H K L SIN(TH)/L D(HKL) FSQ SIGFSQ ', &' F SIGF E SIGE RANK'/1X, &' - - - --------- ------ --- ------ ', &' - ---- - ---- ----'/) 876 FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3,I10) NLIST=MIN(NLIST,NDATA) REWIND IO1 DO I=1,NLIST CALL READ2 (IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) S=SINTHL(J,K,L,GINV) D=1/(2*S) WRITE (ILP,876) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,I END DO C C WRITE E-SORTED OUTPUT FILE. C REWIND IO1 DO I=1,NDATA CALL READ2 (IO1,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) WRITE (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE DATA(I)=E END DO CALL SORT (NDATA,DATA,INDEX) CLOSE (UNIT=IO1,STATUS='KEEP') OPEN (UNIT=IO1,FILE='edata.eee',STATUS='NEW',FORM='FORMATTED') DO I=NDATA,1,-1 READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE CALL WRITE1 (IO1, J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) END DO CLOSE (UNIT=IO1,STATUS='KEEP') C C LIST LARGEST SIGNIFICANT E-MAGNITUDES. C WRITE (ILP,899) ATIME,ADATE,TITLE WRITE (ILP,'(/1X, &''LARGEST SIGNIFICANT E-MAGNITUDES:''//1X, &''E .GE. CUTOFF*SIGMA(E),''//1X, &''WHERE CUTOFF = '',F5.2)') CUTOFF WRITE (ILP,877) NLIST=MIN(NLIST,NDATA) N=0 DO 31 I=NDATA,1,-1 READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE IF (E.LT.CUTOFF*SIGE) GO TO 31 IF (E.LT.EMIN) GO TO 41 S=SINTHL(J,K,L,GINV) D=1/(2*S) WRITE (ILP,876) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,NDATA-I+1 N=N+1 IF (N.EQ.NLIST) GO TO 41 31 CONTINUE 41 CONTINUE C C LIST SMALLEST SIGNIFICANT E-MAGNITUDES. C WRITE (ILP,899) ATIME,ADATE,TITLE WRITE (ILP,'(/1X, &''SMALLEST SIGNIFICANT E-MAGNITUDES:''//1X, &''E .GE. CUTOFF*SIGMA(E),''//1X, &''WHERE CUTOFF = '',F5.2)') CUTOFF WRITE (ILP,877) N=0 DO 32 I=1,NDATA,+1 READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE IF (E.LT.CUTOFF*SIGE) GO TO 32 IF (E.GT.EMAX) GO TO 42 S=SINTHL(J,K,L,GINV) D=1/(2*S) WRITE (ILP,876) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,-I N=N+1 IF (N.EQ.NLIST) GO TO 42 32 CONTINUE 42 CONTINUE CLOSE (UNIT=IO2,STATUS='DELETE') WRITE (ILP,879) 879 FORMAT (/1X, &'OUTPUT E-FILES ARE FORMATTED, ASCII FILES (IH,IK,IL,FSQ,SIGFSQ,F, &SIGF,E,SIGE)'/1X, &'NAMED:'/1X, &' "edata.hkl" hkl-SORTED'/1X, &' "edata.ddd" d(hkl)-SORTED'/1X, &' "edata.eee" E(hkl)-SORTED') WRITE (6,879) STOP 'PROGRAM EVAL FINIS' END C----------------------------------------------------------------------- SUBROUTINE CSYM (PTGP,H,K,L,C) C C TABLE OF CENTROSYMMETRIC ROWS AND ZONES IN NON-CENTROSYMMETRIC C RECIPROCAL LATTICE POINT GROUPS C C TRICL. MONOCL. ORTHO. TETRAG. TRIG. RHOMB. CUB. C C (1) 1 (3) 2 (6) 222 (9) 4 (17) 3 (33) R3 (38) 23 C (2) -1 (4) M (7) MM2 (10) -4 (18) -3 (34) R-3 (39) M3 C (5) 2/M (8) MMM (11) 4/M C (19) 312 (35) R32 (40) 432 C (12) 422 (20) 321 (36) R3M (41) -43M C (13) 4MM (21) 31M (37) R-3M (42) M3M C (14) -42M (22) 3M1 C (15) -4M2 (23) -31M C (16) 4/MMM (24) -3M1 C C HEXAG. C C (25) 6 C (26) -6 C (27) 6/M C C (28) 622 C (29) 6MM C (30) -6M2 C (31) -62M C (32) 6/MMM C C A TOTAL OF 42 CASES IS TABULATED BECAUSE THERE ARE TEN CASES OF C ALTERNATIVE AXES FOR SEVEN OF THE 32 CRYSTALLOGRAPHIC POINT GROUPS: C C -42M -4M2 C 3 R3 C -3 R-3 C 312 321 R32 C 31M 3M1 R3M C -31M -3M1 R-3M C -62M -6M2 C C C = 0 FOR ACENTRIC DISTRIBUTIONS C C = 1 FOR CENTRIC DISTRIBUTIONS C INTEGER PTGP,H,C C=0 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,39,40,41,42) PTGP 1 RETURN 2 C=1 RETURN 3 IF (K.EQ.0) C=1 RETURN 4 IF (H.EQ.0.AND.L.EQ.0) C=1 RETURN 5 C=1 RETURN 6 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 RETURN 7 IF (L.EQ.0) C=1 RETURN 8 C=1 RETURN 9 IF (L.EQ.0) C=1 RETURN 10 IF (L.EQ.0) C=1 IF (H.EQ.0.AND.K.EQ.0) C=1 RETURN 11 C=1 RETURN 12 IF (L.EQ.0.OR.K.EQ.0.OR.H.EQ.0) C=1 IF (ABS(K).EQ.ABS(H)) C=1 RETURN 13 IF (L.EQ.0) C=1 RETURN 14 IF (L.EQ.0.OR.K.EQ.0.OR.H.EQ.0) C=1 RETURN 15 IF (L.EQ.0) C=1 IF (ABS(K).EQ.ABS(H)) C=1 RETURN 16 C=1 RETURN 17 RETURN 18 C=1 RETURN 19 IF (K.EQ.H.OR.K.EQ.-2*H) C=1 RETURN 20 IF (K.EQ.0.OR.H.EQ.0) C=1 RETURN 21 IF ((K.EQ.0.AND.L.EQ.0).OR.(H.EQ.0.AND.L.EQ.0)) C=1 RETURN 22 RETURN 23 C=1 RETURN 24 C=1 RETURN 25 IF (L.EQ.0) C=1 RETURN 26 IF (H.EQ.0.AND.K.EQ.0) C=1 RETURN 27 C=1 RETURN 28 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 IF (K.EQ.H.OR.K.EQ.-2*H) C=1 RETURN 29 IF (L.EQ.0) C=1 RETURN 30 IF (K.EQ.H.OR.K.EQ.-2*H) C=1 RETURN 31 IF (H.EQ.0.OR.K.EQ.0) C=1 RETURN 32 C=1 RETURN 33 RETURN 34 C=1 RETURN 35 IF (H.EQ.K.OR.K.EQ.L.OR.L.EQ.H) C=1 RETURN 36 IF (L.EQ.0.AND.K.EQ.-H) C=1 IF (K.EQ.0.AND.L.EQ.-H) C=1 IF (H.EQ.0.AND.L.EQ.-K) C=1 RETURN 37 C=1 RETURN 38 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 RETURN 39 C=1 RETURN 40 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 IF (ABS(H).EQ.ABS(K).OR.ABS(K).EQ.ABS(L).OR.ABS(L).EQ.ABS(H)) C=1 RETURN 41 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 RETURN 42 C=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP) C C CRYSTAL SYSTEM, LATTICE TYPE AND MULTIPLICITY, AND LAUE GROUP AND C POINT GROUP SYMMETRIES C C THERE ARE TEN CASES OF ALTERNATIVE AXES FOR SEVEN OF THE 32 C CRYSTALLOGRAPHIC POINT GROUPS (SEE SUBROUTINE ESYM), THUS A TOTAL OF C 42 CASES IS TABULATED. C C SIMILARLY, THERE ARE THREE CASES OF ALTERNATIVE AXES FOR TWO OF THE 11 C LAUE GROUPS, SO 14 CASES ARE TABULATED. C CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5,SYMBOL(42)*5 DATA SYMBOL &/'1 ','-1 ', & '2 ','M ','2/M ', & '222 ','MM2 ','MMM ', & '4 ','-4 ','4/M ', & '422 ','4MM ','-42M ','-4M2 ','4/MMM', & '3 ','-3 ', & '312 ','321 ','31M ','3M1 ','-31M ','-3M1 ', & '6 ','-6 ','6/M ', & '622 ','6MM ','-6M2 ','-62M ','6/MMM', & 'R3 ','R-3 ', & 'R32 ','R3M ','R-3M ', & '23 ','M3 ', & '432 ','-43M ','M3M '/ C C POINT GROUP INDEX NUMBERS: C C TRICLINIC 1, 2; C MONOCLINIC 3, 4, 5; C ORTHORHOMBIC 6, 7, 8; C TETRAGONAL 9, 10, 11, C 12, 13, 14, 15, 16; C TRIGONAL 17, 18, C 19, 20, 21, 22, 23, 24; C HEXAGONAL 25, 26, 27, C 28, 29, 30, 31, 32; C RHOMBOHEDRAL 33, 34, C 35, 36, 37; C CUBIC 38, 39, C 40, 41, 42. C DO 1 I=1,42 IF (PTGP.EQ.SYMBOL(I)) GO TO 2 1 CONTINUE STOP 'EVAL/DECODE: ERROR IN POINT GROUP SYMBOL' 2 IPTGP=I C C LAUE GROUP C IF (I.GE. 1) LAUE='-1 ' IF (I.GE. 3) LAUE='2/M ' IF (I.GE. 6) LAUE='MMM ' IF (I.GE. 9) LAUE='4/M ' IF (I.GE.12) LAUE='4/MMM' IF (I.GE.17) LAUE='-3 ' IF (I.GE.19) LAUE='-31M ' IF (I.EQ.20.OR.I.EQ.22.OR.I.EQ.24) LAUE='-3M1 ' IF (I.GE.25) LAUE='6/M ' IF (I.GE.28) LAUE='6/MMM' IF (I.GE.33) LAUE='R-3 ' IF (I.GE.35) LAUE='R-3M ' IF (I.GE.38) LAUE='M3 ' IF (I.GE.40) LAUE='M3M ' C C CRYSTAL SYSTEM C IF (I.GE. 1) SYST='TRICLINIC' IF (I.GE. 3) SYST='MONOCLINIC' IF (I.GE. 6) SYST='ORTHORHOMBIC' IF (I.GE. 9) SYST='TETRAGONAL' IF (I.GE.17) SYST='TRIGONAL' IF (I.GE.25) SYST='HEXAGONAL' IF (I.GE.33) SYST='RHOMBOHEDRAL' IF (I.GE.38) SYST='CUBIC' C C LATTICE CENTERING MULTIPLICITY C IF (LATT.EQ.'P') MLATT=1 IF (LATT.EQ.'A') MLATT=2 IF (LATT.EQ.'B') MLATT=2 IF (LATT.EQ.'C') MLATT=2 IF (LATT.EQ.'F') MLATT=4 IF (LATT.EQ.'I') MLATT=2 IF (LATT.EQ.'R') MLATT=1 IF (LATT.EQ.'R'.AND.SYST.NE.'RHOMBOHEDRAL') MLATT=3 RETURN END C----------------------------------------------------------------------- FUNCTION DELTA(I,J) C C KRONECKER DELTA C DELTA=0 IF (I.EQ.J) DELTA=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE DPRINT (IH,IK,IL,IPRINT) C C DECIDE TO PRINT OR NOT. C C IPRINT = 0 DO NOT C = 1 DO PRINT REFLECIOTN C C PRINT SPECIAL REFLECTIONS OF THE TYPE H00, 0K0, 00L, HH0, H0H,0KK, AND C HHH. C IPRINT=0 J=ABS(IH) K=ABS(IK) L=ABS(IL) IF (J.EQ.0.AND.K.EQ.0) GO TO 1 IF (J.EQ.0.AND.L.EQ.0) GO TO 1 IF (K.EQ.0.AND.L.EQ.0) GO TO 1 IF (J.EQ.K.AND.L.EQ.0) GO TO 1 IF (J.EQ.L.AND.K.EQ.0) GO TO 1 IF (K.EQ.L.AND.J.EQ.0) GO TO 1 IF (J.EQ.K.AND.K.EQ.L) GO TO 1 RETURN 1 IPRINT=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE ESTATS (IFILE,ILP,SMIN,SMAX,IPTGP,MLATT,GINV) C C E-STATISTICS THEORETICAL MOMENTS C ACENTRIC CENTRIC HYPERCENTRIC C C A1 = 0.866 0.798 0.718 C A2 = 1.000 1.000 1.000 C A3 = 1.329 1.596 1.916 C A4 = 2.000 3.000 4.500 C A5 = 3.323 6.383 12.260 C A6 = 6.000 15.000 37.500 C B1 = 0.736 0.968 1.145 C B2 = <(E**2 - 1)**2> 1.000 2.000 3.500 C B3 = <(E**2 - 1)**3> 2.000 8.000 26.000 C C3 = 2.415 8.691 26.903 C C F1 = FRACTION E .GE. 1 0.368 0.320 C F2 2 0.018 0.050 C F3 3 0.0001 0.003 C C F05 = FRACTION E .LT. 0.5 0.22 0.38 C F15 = .GT. 1.5 0.105 0.13 C DIMENSION F1(3),F2(3),F3(3),NI(7),A1(7),A2(7),A3(7),A4(7),A5(7), & A6(7),B1(7),B2(7),B3(7),C3(7),GINV(3,3),H(3),SI(10),NJ(10), & SJ(10),E1(10),E2(10),F05(10),F15(10) DATA F1,F2,F3,NI,A1,A2,A3,A4,A5,A6,B1,B2,B3,C3,NJ,SJ,E1,E2,F05,F15 & /146*0/ DIMENSION W(50),X(50),Y(50),Z(50) C C MOMENTS OF THE E-DISTRIBUTION C REWIND IFILE IEND=0 NREC=0 NA=0 NC=0 1 CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) IF (IEND.NE.0) GO TO 9 NREC=NREC+1 S=SINTHL(J,K,L,GINV) IF (SMIN.LE.S.AND.S.LE.SMAX) THEN CALL ESYM (IPTGP,J,K,L,EPSILON,M) CALL CSYM (IPTGP,J,K,L,IC) IF (IC.EQ.0) M=2*M DO 10 I=1,7 IF (I.EQ.2.AND.IC.NE.0) GO TO 10 IF (I.EQ.3.AND.IC.NE.1) GO TO 10 IF (I.EQ.4.AND.(J.EQ.0.OR.K.EQ.0.OR.L.EQ.0)) GO TO 10 IF (I.EQ.5.AND.J.NE.0) GO TO 10 IF (I.EQ.6.AND.K.NE.0) GO TO 10 IF (I.EQ.7.AND.L.NE.0) GO TO 10 NI(I)=NI(I)+M A1(I)=A1(I)+M*E A2(I)=A2(I)+M*E**2 A3(I)=A3(I)+M*E**3 A4(I)=A4(I)+M*E**4 A5(I)=A5(I)+M*E**5 A6(I)=A6(I)+M*E**6 B1(I)=B1(I)+M*ABS(E**2-1) B2(I)=B2(I)+M*(E**2-1)**2 B3(I)=B3(I)+M*(E**2-1)**3 C3(I)=C3(I)+M*ABS(E**2-1)**3 IF (I.LE.3) THEN IF (E.GE.1) F1(I)=F1(I)+M IF (E.GE.2) F2(I)=F2(I)+M IF (E.GE.3) F3(I)=F3(I)+M IF (E.LT.0.5) F05(I)=F05(I)+M IF (E.GT.1.5) F15(I)=F15(I)+M END IF 10 CONTINUE IF (IC.EQ.0) NA=NA+1 IF (IC.EQ.1) NC=NC+1 END IF GO TO 1 9 CONTINUE DO I=1,7 IF (NI(I).GT.0) THEN A1(I)=A1(I)/NI(I) A2(I)=A2(I)/NI(I) A3(I)=A3(I)/NI(I) A4(I)=A4(I)/NI(I) A5(I)=A5(I)/NI(I) A6(I)=A6(I)/NI(I) B1(I)=B1(I)/NI(I) B2(I)=B2(I)/NI(I) B3(I)=B3(I)/NI(I) C3(I)=C3(I)/NI(I) END IF END DO DO I=1,3 IF (NI(I).GT.0) THEN F1(I)=F1(I)/NI(I) F2(I)=F2(I)/NI(I) F3(I)=F3(I)/NI(I) F05(I)=F05(I)/NI(I) F15(I)=F15(I)/NI(I) END IF END DO WRITE (ILP,900) WRITE (ILP,901) A1, 0.866, 0.798, 0.718 WRITE (ILP,902) A2, 1.000, 1.000, 1.000 WRITE (ILP,903) A3, 1.329, 1.596, 1.916 WRITE (ILP,904) A4, 2.000, 3.000, 4.500 WRITE (ILP,905) A5, 3.323, 6.383, 12.260 WRITE (ILP,906) A6, 6.000, 15.000, 37.500 WRITE (ILP,907) B1, 0.736, 0.968, 1.145 WRITE (ILP,908) B2, 1.000, 2.000, 3.500 WRITE (ILP,909) B3, 2.000, 8.000, 26.000 WRITE (ILP,910) C3, 2.415, 8.691, 26.903 WRITE (ILP,911) NI WRITE (ILP,912) F1, 0.368, 0.320, & F2, 0.018, 0.050, & F3, 0.0001, 0.003, & (F05(I),I=1,3), 0.22, 0.38, & (F15(1),I=1,3), 0.105, 0.13 900 FORMAT (/1X, &' EXPERIMENTAL DATA AV', &'ERAGES THEORETICAL MOMENTS'/1X, &' -----------------------------------------', &'--------------------------- -------------------------------'/1X, &'AVERAGE ALL DATA ACENTRIC CENTRIC HKL ', &' 0KL H0L HK0 ACENTRIC CENTRIC HYPERCENTRIC'/1X, &' -------- -------- ------- --- ', &' --- --- --- -------- ------- ------------') 901 FORMAT (' ',10F10.3) 902 FORMAT (' ',10F10.3) 903 FORMAT (' ',10F10.3) 904 FORMAT (' ',10F10.3) 905 FORMAT (' ',10F10.3) 906 FORMAT (' ',10F10.3) 907 FORMAT (' ',10F10.3) 908 FORMAT (' <(E**2 - 1)**2> ',10F10.3) 909 FORMAT (' <(E**2 - 1)**3> ',10F10.3) 910 FORMAT (' ',10F10.3) 911 FORMAT (' --------------------------------', &'------------------------------------ -------------------------', &'------'/1X,'MULT.-WTD. NO. ',7I10) 912 FORMAT (/1X, &' EXPERIMENTAL ', &' THEORETICAL'/1X, &' ---------------------------- -', &'-----------------'/1X, &' ALL DATA ACENTRIC CENTRIC A', &'CENTRIC CENTRIC'/1X, &' -------- -------- ------- -', &'------- -------'/1X, &'FRACTION OF DATA WITH E .GE. 1',5F10.3/1X, &' 2',5F10.3/1X, &' 3',F11.4,F10.4,F9.3,F11.4,F9.3/1X, &' ---------------------------- ', &'------------------'/1X, &'FRACTION OF DATA WITH E .LT. 0.5',F7.2,1X,4(F9.2,1X)/1X, &' E .GT. 1.5',F7.2,1X,2(F9.2,1X),F10.3,F9.2/ &1X, &' ---------------------------- -', &'-----------------') C C E-STATISTICS IN EQUAL-VOLUME SHELLS OF SIN(THETA)/LAMBDA C N=10 DO I=1,N SI(I)=(FLOAT(I)/N)**(1/3.0)*SMAX END DO REWIND IFILE DO 2 IREC=1,NREC CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) S=SINTHL(J,K,L,GINV) IF (SMIN.LE.S.AND.S.LE.SMAX) THEN CALL ESYM (IPTGP,J,K,L,EPSILON,M) CALL CSYM (IPTGP,J,K,L,IC) IF (IC.EQ.0) M=2*M DO I=1,N IF (S.LE.SI(I)) THEN NJ(I)=NJ(I)+M SJ(I)=SJ(I)+M*S E1(I)=E1(I)+M*E E2(I)=E2(I)+M*E**2 IF (E.LT.0.5) F05(I)=F05(I)+M IF (E.GT.1.5) F15(I)=F15(I)+M GO TO 2 END IF END DO END IF 2 CONTINUE DO I=1,N IF (NJ(I).GT.0) THEN SJ(I)=SJ(I)/NJ(I) E1(I)=E1(I)/NJ(I) E2(I)=E2(I)/NJ(I) F05(I)=F05(I)/NJ(I) F15(I)=F15(I)/NJ(I) END IF END DO WRITE (ILP,920) (2*SJ(I),E2(I),E1(I),F05(I),F15(I),NJ(I),I=1,N) 920 FORMAT (/1X, &'E-STATISTICS IN EQUAL-VOLUME SHELLS OF DSTAR = 2*SIN(THETA)/LAM', &'BDA'//1X, &' F(E < 0.5) F(E > 1.5) MULT.-WTD.', &' NO.'/1X, &' ------- ------ --- ---------- ---------- ----------', &'----'/(1X,3F10.3,2F10.2,I10)) C C E-STATISTICS IN INDEX PARITY GROUPS C DO I=1,8 NJ(I)=0 E1(I)=0 E2(I)=0 F05(I)=0 F15(I)=0 END DO REWIND IFILE DO IREC=1,NREC CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) S=SINTHL(J,K,L,GINV) IF (SMIN.LE.S.AND.S.LE.SMAX) THEN CALL ESYM (IPTGP,J,K,L,EPSILON,M) CALL CSYM (IPTGP,J,K,L,IC) IF (IC.EQ.0) M=2*M CALL PARITY (I,J,K,L) NJ(I)=NJ(I)+M E1(I)=E1(I)+M*E E2(I)=E2(I)+M*E**2 IF (E.LT.0.5) F05(I)=F05(I)+M IF (E.GT.1.5) F15(I)=F15(I)+M END IF END DO DO I=1,8 IF (NJ(I).GT.0) THEN E1(I)=E1(I)/NJ(I) E2(I)=E2(I)/NJ(I) F05(I)=F05(I)/NJ(I) F15(I)=F15(I)/NJ(I) END IF END DO WRITE (ILP,931) (E2(I),E1(I),F05(I),F15(I),NJ(I),I=1,8) 931 FORMAT (/1X, &'E-STATISTICS FOR INDEX PARITY GROUPS'//1X, &' PARITY F(E < 0.5) F(E > 1.5) MULT.-WTD.', &' NO.'/1X, &' ------ ------ --- ---------- ---------- ----------', &'----'/ &' EEE',2F10.3,2F10.2,I10/ &' EEO',2F10.3,2F10.2,I10/ &' EOE',2F10.3,2F10.2,I10/ &' OEE',2F10.3,2F10.2,I10/ &' EOO',2F10.3,2F10.2,I10/ &' OEO',2F10.3,2F10.2,I10/ &' OOE',2F10.3,2F10.2,I10/ &' OOO',2F10.3,2F10.2,I10) C C PLOT THE E-DISTRIBUTIONS P(E) VS. E, WHERE 0 < P < 1 AND 0 < E < 5.0. C IF (NA.GT.0) THEN C C ACENTRIC DISTRIBUTION C N=MIN(50,NINT(0.1*NA)) DO I=1,N X(I)=0 Y(I)=0 END DO DX=5.0/N REWIND IFILE DO 20 IREC=1,NREC CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) S=SINTHL(J,K,L,GINV) IF (S.LT.SMIN.OR.S.GT.SMAX) GO TO 20 CALL CSYM (IPTGP,J,K,L,IC) IF (IC.EQ.1) GO TO 20 CALL ESYM (IPTGP,J,K,L,EPSILON,M) M=2*M DO I=1,N IF (E.LE.I*DX) THEN X(I)=X(I)+M*E Y(I)=Y(I)+M GO TO 20 END IF END DO 20 CONTINUE C C EVALUATE BIN AVERAGES, AND ELIMINATE ANY EMPTY BINS. C M=N N=0 DO I=1,M IF (Y(I).GT.0) THEN N=N+1 X(N)=X(I)/Y(I) Y(N)=Y(I) END IF END DO C C ENSURE INCREASING X VALUES. C M=N N=1 DO I=2,M IF (X(I).GT.X(I-1)) THEN N=N+1 X(N)=X(I) Y(N)=Y(I) END IF END DO C C NORMALIZE BIN AREAS. C A=Y(1)*(X(2)-X(1))+Y(N)*(X(N)-X(N-1)) DO I=2,N-1 A=A+0.5*Y(I)*(X(I+1)-X(I-1)) END DO DO I=1,N Y(I)=Y(I)/A END DO WRITE (ILP,600) CALL PLOTA (N,X,Y,ILP) END IF 600 FORMAT (/1H1,130('-')/1X,10X, &'PROGRAM EVAL. ACENTRIC E-DISTRIBUTION'/) IF (NC.GT.0) THEN C C CENTRIC DISTRIBUTION C N=MIN(50,NINT(0.1*NC)) DO I=1,N X(I)=0 Y(I)=0 END DO DX=5.0/N REWIND IFILE DO 21 IREC=1,NREC CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) S=SINTHL(J,K,L,GINV) IF (S.LT.SMIN.OR.S.GT.SMAX) GO TO 21 CALL ESYM (IPTGP,J,K,L,EPSILON,M) CALL CSYM (IPTGP,J,K,L,IC) IF (IC.EQ.0) GO TO 21 DO I=1,N IF (E.LE.I*DX) THEN X(I)=X(I)+M*E Y(I)=Y(I)+M GO TO 21 END IF END DO 21 CONTINUE C C EVALUATE BIN AVERAGES, AND ELIMINATE ANY EMPTY BINS. C M=N N=0 DO I=1,M IF (Y(I).GT.0) THEN N=N+1 X(N)=X(I)/Y(I) Y(N)=Y(I) END IF END DO C C ENSURE INCREASING X VALUES. C M=N N=1 DO I=2,M IF (X(I).GT.X(I-1)) THEN N=N+1 X(N)=X(I) Y(N)=Y(I) END IF END DO C C NORMALIZE BIN AREAS. C A=Y(1)*(X(2)-X(1))+Y(N)*(X(N)-X(N-1)) DO I=2,N-1 A=A+0.5*Y(I)*(X(I+1)-X(I-1)) END DO DO I=1,N Y(I)=Y(I)/A END DO WRITE (ILP,610) CALL PLOTC (N,X,Y,ILP) END IF 610 FORMAT (/1H1,130('-')/1X,10X, &'PROGRAM EVAL. CENTRIC E-DISTRIBUTION'/) C C PLOT VERSUS = 2*. C N=MIN(50,NINT(FLOAT(NREC)/50)) DO I=1,N W(I)=(FLOAT(I)/N)**(1/3.0)*SMAX END DO DO I=1,N X(I)=0 Y(I)=0 Z(I)=0 END DO REWIND IFILE DO 3 IREC=1,NREC CALL READ2 (IFILE,IEND,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) S=SINTHL(J,K,L,GINV) IF (S.LT.SMIN.OR.S.GT.SMAX) GO TO 3 CALL ESYM (IPTGP,J,K,L,EPSILON,M) CALL CSYM (IPTGP,J,K,L,IC) IF (IC.EQ.0) M=2*M DO I=1,N IF (S.LE.W(I)) THEN X(I)=X(I)+M*S Y(I)=Y(I)+M*E**2 Z(I)=Z(I)+M GO TO 3 END IF END DO 3 CONTINUE C C EVALUATE BIN AVERAGES, AND ELIMINATE ANY EMPTY BINS. C M=N N=0 DO I=1,M IF (Z(I).GT.0) THEN N=N+1 X(N)=X(I)/Z(I) Y(N)=Y(I)/Z(I) END IF END DO C C ENSURE INCREASING X VALUES. C M=N N=1 DO I=2,M IF (X(I).GT.X(I-1)) THEN N=N+1 X(N)=X(I) Y(N)=Y(I) END IF END DO WRITE (ILP,630) CALL PLOTS (N,X,Y,SMAX,ILP) 630 FORMAT (/1H1,130('-')/1X,10X, &'PROGRAM EVAL. VERSUS = 2* DI', &'STRIBUTION'/) RETURN END C----------------------------------------------------------------------- SUBROUTINE ESYM (PTGP,H,K,L,E,M) C C RECIPROCAL LATTICE POINT DEGENERACIES, E = EPSILON(HKL), AND C MULTIPLICITIES, M = M(HKL), FOR THE CRYSTALLOGRAPHIC POINT GROUPS C C AS TABULATED BY: C C HITOSHI IWASAKI AND TETSUZO ITO (1977). VALUES OF EPSILON FOR C OBTAINING NORMALIZED STRUCTURE FACTORS. ACTA CRYST. A33, 227-229. C C TRICL. MONOCL. ORTHO. TETRAG. TRIG. RHOMB. CUB. C C (1) 1 (3) 2 (6) 222 (9) 4 (17) 3 (33) R3 (38) 23 C (2) -1 (4) M (7) MM2 (10) -4 (18) -3 (34) R-3 (39) M3 C (5) 2/M (8) MMM (11) 4/M C (19) 312 (35) R32 (40) 432 C (12) 422 (20) 321 (36) R3M (41) -43M C (13) 4MM (21) 31M (37) R-3M (42) M3M C (14) -42M (22) 3M1 C (15) -4M2 (23) -31M C (16) 4/MMM (24) -3M1 C C HEXAG. C C (25) 6 C (26) -6 C (27) 6/M C C (28) 622 C (29) 6MM C (30) -6M2 C (31) -62M C (32) 6/MMM C C A TOTAL OF 42 CASES IS TABULATED BECAUSE THERE ARE TEN CASES OF C ALTERNATIVE AXES FOR SEVEN OF THE 32 CRYSTALLOGRAPHIC POINT GROUPS: C C -42M -4M2 C 3 R3 C -3 R-3 C 312 321 R32 C 31M 3M1 R3M C -31M -3M1 R-3M C -62M -6M2 C INTEGER PTGP,H E=1 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,39,40,41,42) PTGP 1 M=1 GO TO 99 2 M=2 GO TO 99 3 M=2 IF (H.EQ.0.AND.L.EQ.0) E=2 GO TO 99 4 M=2 IF (K.EQ.0) E=2 GO TO 99 5 M=4 IF (K.EQ.0) E=2 IF (H.EQ.0.AND.L.EQ.0) E=2 GO TO 99 6 M=4 IF (K.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=2 GO TO 99 7 M=4 IF (H.EQ.0) E=2 IF (K.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 8 M=8 IF (H.EQ.0) E=2 IF (K.EQ.0) E=2 IF (L.EQ.0) E=2 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 9 M=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 10 M=4 IF (H.EQ.0.AND.K.EQ.0) E=2 GO TO 99 11 M=8 IF (L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 12 M=8 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 13 M=8 IF (K.EQ.0.OR.H.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=8 GO TO 99 14 M=8 IF (K.EQ.H.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 15 M=8 IF (K.EQ.0.OR.H.EQ.0) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 16 M=16 IF (K.EQ.0.OR.H.EQ.0.OR.L.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=4 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=8 GO TO 99 17 M=3 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 18 M=6 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 19 M=6 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 20 M=6 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 21 M=6 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 22 M=6 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 23 M=12 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 24 M=12 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 25 M=6 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 26 M=6 IF (L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 27 M=12 IF (L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 28 M=12 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 29 M=12 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=12 GO TO 99 30 M=12 IF (L.EQ.0) E=2 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 31 M=12 IF (L.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 32 M=24 IF (L.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=4 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=12 GO TO 99 33 M=3 IF (K.EQ.H.AND.L.EQ.H) E=3 GO TO 99 34 M=6 IF (K.EQ.H.AND.L.EQ.H) E=3 GO TO 99 35 M=6 IF (L.EQ.0.AND.K.EQ.-H) E=2 IF (K.EQ.0.AND.L.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.-L) E=2 IF (K.EQ.H.AND.L.EQ.H) E=3 GO TO 99 36 M=6 IF (K.EQ.H.OR.L.EQ.H.OR.L.EQ.K) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=2 IF (K.EQ.H.AND.L.EQ.H) E=6 GO TO 99 37 M=12 IF (K.EQ.H.OR.L.EQ.H.OR.K.EQ.L) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=2 IF (K.EQ.-H.AND.L.EQ.0) E=2 IF (L.EQ.-H.AND.K.EQ.0) E=2 IF (L.EQ.-K.AND.H.EQ.0) E=2 IF (K.EQ.H.AND.L.EQ.H) E=6 GO TO 99 38 M=12 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3 IF (K.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=2 GO TO 99 39 M=24 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) E=2 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 40 M=24 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3 IF (L.EQ.0.AND.ABS(K).EQ.ABS(H)) E=2 IF (K.EQ.0.AND.ABS(L).EQ.ABS(H)) E=2 IF (H.EQ.0.AND.ABS(K).EQ.ABS(L)) E=2 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 41 M=24 IF (ABS(K).EQ.ABS(H).OR.ABS(L).EQ.ABS(H).OR.ABS(K).EQ.ABS(L)) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=6 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 42 M=48 IF (ABS(K).EQ.ABS(H).OR.ABS(L).EQ.ABS(H).OR.ABS(K).EQ.ABS(L)) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=6 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) E=2 IF (L.EQ.0.AND.ABS(K).EQ.ABS(H)) E=4 IF (K.EQ.0.AND.ABS(L).EQ.ABS(H)) E=4 IF (H.EQ.0.AND.ABS(K).EQ.ABS(L)) E=4 IF (K.EQ.0.AND.L.EQ.0) E=8 IF (H.EQ.0.AND.L.EQ.0) E=8 IF (H.EQ.0.AND.K.EQ.0) E=8 99 M=M/E RETURN END C----------------------------------------------------------------------- SUBROUTINE FCALC (N,X,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S, & SUMFSQ) DIMENSION X(1),A1(1),B1(1),A2(1),B2(1),A3(1),B3(1),A4(1),B4(1), & C0(1),FP(1),FPP(1) SUMFSQ=0 DO 1 I=1,N SUMFSQ=SUMFSQ+X(I)*((A1(I)*EXP(-B1(I)*S**2) & +A2(I)*EXP(-B2(I)*S**2) & +A3(I)*EXP(-B3(I)*S**2) & +A4(I)*EXP(-B4(I)*S**2)+C0(I)+FP(I))**2 & +FPP(I)**2) 1 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE FSQTOF (FSQ,SIGFSQ,F,SIGF) IF (FSQ.LE.0) THEN F=0 ELSE F=SQRT(FSQ) END IF IF (FSQ.LE.SIGFSQ) THEN SIGF=0.5*SQRT(SIGFSQ) ELSE SIGF=0.5*SIGFSQ/F END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE FTOFSQ (F,SIGF,FSQ,SIGFSQ) IF (F.LT.0) F=0 IF (SIGF.LT.0) SIGF=0 FSQ=F**2 SIGFSQ=MAX(2*F*SIGF,4*SIGF**2) RETURN END C----------------------------------------------------------------------- SUBROUTINE METRIC (A,B,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) RAD=PI/180 DEG=180/PI 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 PARITY (I,J,K,L) C C 1 EEE C 2 EEO C 3 EOE C 4 OEE C 5 EOO C 6 OEO C 7 OOE C 8 OOO C P=MOD(ABS(J),2) Q=MOD(ABS(K),2) R=MOD(ABS(L),2) IF (P.EQ.0.AND.Q.EQ.0.AND.R.EQ.0) I=1 IF (P.EQ.0.AND.Q.EQ.0.AND.R.EQ.1) I=2 IF (P.EQ.0.AND.Q.EQ.1.AND.R.EQ.0) I=3 IF (P.EQ.1.AND.Q.EQ.0.AND.R.EQ.0) I=4 IF (P.EQ.0.AND.Q.EQ.1.AND.R.EQ.1) I=5 IF (P.EQ.1.AND.Q.EQ.0.AND.R.EQ.1) I=6 IF (P.EQ.1.AND.Q.EQ.1.AND.R.EQ.0) I=7 IF (P.EQ.1.AND.Q.EQ.1.AND.R.EQ.1) I=8 RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOTA (N,X,Y,ILP) PARAMETER (NX=100,NY=50) CHARACTER AA*1 DIMENSION AA(0:NX,0:NY),X(N),Y(N) C C BLANK OUT THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C FILL IN THE GRID MARKS. C DO IX=0,NX AA(IX, 0)='.' AA(IX,NY)='.' END DO DO IX=0,NX,10 AA(IX, 0)='+' AA(IX,NY)='+' END DO DO IY=0,NY AA( 0,IY)='.' AA(NX,IY)='.' END DO DO IY=0,NY,10 AA( 0,IY)='+' AA(NX,IY)='+' END DO C C SET THE PLOT LIMITS, 0 < E < 5.0 AND 0 < P(E) < 1. C XMIN=0 XMAX=5 YMIN=0 YMAX=1 XDEL=XMAX-XMIN YDEL=YMAX-YMIN C C FILL IN THE DISTRIBUTION FUNCTION. C DX=XDEL/(2*NX) DO I=0,(2*NX) XI=I*DX YI=2*XI*EXP(-XI**2) IX=NINT(NX*(XI-XMIN)/XDEL) IY=NINT(NY*(YMAX-YI)/YDEL) IX=MAX(IX,0) IX=MIN(IX,NX) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='*' END DO C C FILL IN THE EXPERIMENTAL DATA POINTS. C DO I=1,N IX=NINT(NX*(X(I)-XMIN)/XDEL) IY=NINT(NY*(YMAX-Y(I))/YDEL) IX=MAX(IX,0) IX=MIN(IX,NX) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='O' END DO C C PRINT THE PLOT ARRAY. C DY=(YMAX-YMIN)/NY DO IY=0,NY IF (MOD(IY,10).EQ.0) THEN WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,101) (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,102) (XMIN+IX*DX,IX=0,10) 100 FORMAT (' ',F10.1,101A1) 101 FORMAT (' ',10X, 101A1) 102 FORMAT (' ',1X,11F10.1) WRITE (ILP,'(/1X,10X,''ACENTRIC DISTRIBUTION: P(E) = 2*E*EXP(-E** &2), "*" THEORETICAL, "O" EXPERIMENTAL'')') RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOTC (N,X,Y,ILP) PARAMETER (NX=100,NY=50) CHARACTER AA*1 DIMENSION AA(0:NX,0:NY),X(N),Y(N) C C BLANK OUT THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C FILL IN THE GRID MARKS. C DO IX=0,NX AA(IX, 0)='.' AA(IX,NY)='.' END DO DO IX=0,NX,10 AA(IX, 0)='+' AA(IX,NY)='+' END DO DO IY=0,NY AA( 0,IY)='.' AA(NX,IY)='.' END DO DO IY=0,NY,10 AA( 0,IY)='+' AA(NX,IY)='+' END DO C C SET THE PLOT LIMITS, 0 < E < 5.0 AND 0 < P(E) < 1. C XMIN=0 XMAX=5 YMIN=0 YMAX=1 XDEL=XMAX-XMIN YDEL=YMAX-YMIN C C FILL IN THE DISTRIBUTION FUNCTION. C DX=XDEL/(2*NX) C=SQRT(2/3.141593) DO I=0,(2*NX) XI=I*DX YI=C*EXP(-0.5*XI**2) IX=NINT(NX*(XI-XMIN)/XDEL) IY=NINT(NY*(YMAX-YI)/YDEL) IX=MAX(IX,0) IX=MIN(IX,NX) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='*' END DO C C FILL IN THE EXPERIMENTAL DATA POINTS. C DO I=1,N IX=NINT(NX*(X(I)-XMIN)/XDEL) IY=NINT(NY*(YMAX-Y(I))/YDEL) IX=MAX(IX,0) IX=MIN(IX,NX) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='O' END DO C C PRINT THE PLOT ARRAY. C DY=(YMAX-YMIN)/NY DO IY=0,NY IF (MOD(IY,10).EQ.0) THEN WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,101) (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,102) (XMIN+IX*DX,IX=0,10) 100 FORMAT (' ',F10.1,101A1) 101 FORMAT (' ',10X, 101A1) 102 FORMAT (' ',1X,11F10.1) WRITE (ILP,'(/1X,10X,''CENTRIC DISTRIBUTION: P(E) = SQRT(2/PI)*EX &P(-0.5*E**2), "*" THEORETICAL, "O" EXPERIMENTAL'')') RETURN END C----------------------------------------------------------------------- SUBROUTINE PLOTS (N,X,Y,SMAX,ILP) PARAMETER (NX=100,NY=50) CHARACTER AA*1 DIMENSION AA(0:NX,0:NY),X(N),Y(N),YPP(50) C C BLANK OUT THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C FILL IN THE GRID MARKS. C DO IX=0,NX AA(IX, 0)='.' AA(IX,NY)='.' END DO DO IX=0,NX,10 AA(IX, 0)='+' AA(IX,NY)='+' END DO DO IY=0,NY AA( 0,IY)='.' AA(NX,IY)='.' END DO DO IY=0,NY,10 AA( 0,IY)='+' AA(NX,IY)='+' END DO C C SET THE PLOT LIMITS, 0 < SIN(THETA)/LAMBDA < SMAX AND 0 < E(S) < 2.5. C XMIN=0 XMAX=SMAX YMIN=0 YMAX=2.5 XDEL=XMAX-XMIN YDEL=YMAX-YMIN C C MARK = 1. C DO IX=0,NX AA(IX,30)='.' END DO C C FILL IN THE SPLINE-INTERPOLATED DATA POINTS. C CALL SPLINE (N,X,Y,YPP) DX=(X(N)-X(1))/(2*NX) DO I=0,(2*NX) XI=X(1)+I*DX CALL SPLINT (N,X,Y,YPP,XI,YI) IX=NINT(NX*(XI-XMIN)/XDEL) IY=NINT(NY*(YMAX-YI)/YDEL) IX=MAX(IX,0) IX=MIN(IX,NX) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='*' END DO C C FILL IN THE EXPERIMENTAL DATA POINTS. C DO I=1,N IX=NINT(NX*(X(I)-XMIN)/XDEL) IY=NINT(NY*(YMAX-Y(I))/YDEL) IX=MAX(IX,0) IX=MIN(IX,NX) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='O' END DO C C PRINT THE PLOT ARRAY. C DY=(YMAX-YMIN)/NY DO IY=0,NY IF (MOD(IY,10).EQ.0) THEN WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,101) (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,102) (2*(XMIN+IX*DX),IX=0,10) 100 FORMAT (' ',F10.1,101A1) 101 FORMAT (' ',10X, 101A1) 102 FORMAT (' ',1X,11F10.4) RETURN END C----------------------------------------------------------------------- SUBROUTINE READ1 (ITYPE,IUNIT,IEND,IH,IK,IL,FSQ,SIGFSQ,F,SIGF) FSQ=0 SIGFSQ=0 F=0 SIGF=0 1 IF (ITYPE.EQ.0) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF IF (ITYPE.EQ.1) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ IF (ITYPE.EQ.2) READ (IUNIT,*,END=9) IH,IK,IL,F, SIGF IF (ITYPE.EQ.3) READ (IUNIT, END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF IF (ITYPE.EQ.4) READ (IUNIT, END=9) IH,IK,IL,FSQ,SIGFSQ IF (ITYPE.EQ.5) READ (IUNIT, END=9) IH,IK,IL,F, SIGF IF (SIGFSQ.LE.0.AND.SIGF.LE.0) GO TO 1 IF (FSQ.LE.0.AND.F.LE.0) GO TO 1 IF ((ITYPE.EQ.1.OR.ITYPE.EQ.4).AND.SIGFSQ.GT.0) & CALL FSQTOF (FSQ,SIGFSQ,F,SIGF) IF ((ITYPE.EQ.2.OR.ITYPE.EQ.5).AND.SIGF.GT.0) & CALL FTOFSQ (F,SIGF,FSQ,SIGFSQ) IF (F*SIGF.LE.0) GO TO 1 RETURN 9 IEND=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE READ2 (IFILE,IEND,IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE) READ (IFILE,*,END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE RETURN 9 IEND=1 RETURN END C----------------------------------------------------------------------- FUNCTION SINTHL(IH,IK,IL,GINV) DIMENSION H(3),GINV(3,3) H(1)=IH H(2)=IK H(3)=IL SINTHL=0.5*SQRT(VTMV(H,GINV)) RETURN END C----------------------------------------------------------------------- SUBROUTINE SORT (N,DATA,INDEX) C C INDEXES THE ARRAY DATA(N) AND RETURNS THE ARRAY INDEX(N) SORTED SUCH C THAT THE VALUES DATA(INDEX(I)) ARE IN ASCENDING ORDER FOR I = 1, 2, C ..., N. THE INPUT VARIABLES N AND DATA(N) ARE NOT CHANGED. C C EMPLOYS THE "HEAPSORT" ALGORITHM. C C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A. C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986). NUMERICAL RECIPIES: THE C --------- -------- --- C ART OF SCIENTIFIC COMPUTING, PP. 229-233. CAMBRIDGE, ENGLAND: C --- -- ---------- --------- C CAMBRIDGE UNIVERSITY PRESS. C DIMENSION DATA(N),INDEX(N) DO 10 I=1,N INDEX(I)=I 10 CONTINUE I=N/2+1 J=N 1 CONTINUE IF (I.GT.1) THEN I=I-1 INDEXT=INDEX(I) T=DATA(INDEXT) ELSE INDEXT=INDEX(J) T=DATA(INDEXT) INDEX(J)=INDEX(1) J=J-1 IF (J.EQ.1) THEN INDEX(1)=INDEXT RETURN END IF END IF K=I L=I+I 2 IF (L.LE.J) THEN IF (L.LT.J) THEN IF (DATA(INDEX(L)).LT.DATA(INDEX(L+1))) L=L+1 END IF IF (T.LT.DATA(INDEX(L))) THEN INDEX(K)=INDEX(L) K=L L=L+L ELSE L=J+1 END IF GO TO 2 END IF INDEX(K)=INDEXT GO TO 1 END C----------------------------------------------------------------------- SUBROUTINE SPLINE (N,X,Y,YPP) C C NATURAL CUBIC SPLINE FIT TO A TABULATED FUNCTION Y(I) = Y(X(I)), I = C 1, 2,..., N, WITH X(1) .LT. X(2) .LT. ... .LT. X(N). RETURNS THE C TABULATED SECOND DERIVATIVES y'' = d2y/dx2 = YPP(I), I = 1, 2,..., N, C WITH YPP(1) = YPP(N) = 0. C C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A. C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986). NUMERICAL RECIPIES: THE C ART OF SCIENTIFIC COMPUTING, PP. 86-89. CAMBRIDGE, ENGLAND: C CAMBRIDGE UNIVERSITY PRESS. C C MUST HAVE NMAX .GE. N. C PARAMETER (NMAX=1000) DIMENSION X(N),Y(N),YPP(N),U(NMAX) YPP(1)=0 YPP(N)=0 U(1)=0 U(N)=0 DO I=2,N-1 P=(X(I)-X(I-1))/(X(I+1)-X(I-1)) Q=P*YPP(I-1)+2 YPP(I)=(P-1)/Q U(I)=(6*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))/ & (X(I)-X(I-1)))/(X(I+1)-X(I-1))-P*U(I-1))/Q END DO DO I=N-1,1,-1 YPP(I)=YPP(I)*YPP(I+1)+U(I) END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE SPLINT (N,X,Y,YPP,XI,YI) C C CUBIC SPLINE INTERPOLATION OF A TABULATED FUNCTION Y(I) = Y(X(I)), I = C 1, 2,..., N, USING THE TABULATED SECOND DERIVATIVES YPP(I) FROM C SUBROUTINE SPLINE. C C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A. C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986). NUMERICAL RECIPIES: THE C ART OF SCIENTIFIC COMPUTING, PP. 86-89. CAMBRIDGE, ENGLAND: C CAMBRIDGE UNIVERSITY PRESS. C DIMENSION X(N),Y(N),YPP(N) IF (XI.LE.X(1)) THEN YI=Y(1) RETURN END IF IF (XI.GE.X(N)) THEN YI=Y(N) RETURN END IF I1=1 I2=N 1 IF (I2-I1.GT.1) THEN I=(I2+I1)/2 IF (X(I).GT.XI) THEN I2=I ELSE I1=I END IF GO TO 1 END IF H=X(I2)-X(I1) A=(X(I2)-XI)/H B=(XI-X(I1))/H YI=A*Y(I1)+B*Y(I2)+((A**3-A)*YPP(I1)+(B**3-B)*YPP(I2))*H**2/6 RETURN END C----------------------------------------------------------------------- FUNCTION VARF(SCALEK,B,C,COV,S,H,HTBH,HTCH,Q) DIMENSION B(3,3),C(3,3),COV(13,13),H(3),DQDP(13) DQDP(1)=-Q/SCALEK IF (B(2,2).EQ.0) THEN DQDP(2)=-Q*S**2 DQDP(3)=2*Q*C(1,1)*S**4 N=3 ELSE N=1 DO 1 I=1,3 DO 1 J=I,3 N=N+1 DQDP(N)=-Q*(2-DELTA(I,J))*H(I)*H(J) 1 CONTINUE DO 2 I=1,3 DO 2 J=I,3 N=N+1 DQDP(N)=2*Q*HTCH*(2-DELTA(I,J))*H(I)*H(J) 2 CONTINUE END IF VARF=0 DO 3 I=1,N DO 3 J=I,N VARF=VARF+(2-DELTA(I,J))*DQDP(I)*DQDP(J)*COV(I,J) 3 CONTINUE RETURN END C----------------------------------------------------------------------- FUNCTION VTMV(A,B) C C QUADRATIC FORM, VTMV, EQUALS ROW VECTOR, A-TRANSPOSE, TIMES SQUARE C MATRIX, B, TIMES COLUMN VECTOR, A. C DIMENSION A(3),B(3,3) VTMV=0 DO 1 I=1,3 DO 1 J=1,3 VTMV=VTMV+A(J)*B(I,J)*A(I) 1 CONTINUE RETURN END C----------------------------------------------------------------------- C===================================================================== C Machine-dependent system-service routines (such as TIME and DATE) C===================================================================== SUBROUTINE VFDATE(DAYTIME) C C This routine returns the date/time as a 20 charater string. C Truncation occurs if the length of the input string is less than 20. C The format of the returned string is "dd mmm yyyy hh:mm:ss". C Note: this routine is Y2K compliant. C C +++++++++++++++++++++++++++++++++++++++ C Machine dependent intel/Win32 version. C +++++++++++++++++++++++++++++++++++++++ C IMPLICIT NONE C input/output CHARACTER*(*) DAYTIME C local CHARACTER DATE*8, TIME*10, ZONE*5, MONTH(12)*3, TEMP*20 INTEGER VALUES(8) DATA MONTH/'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', & 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/ C CALL DATE_AND_TIME(DATE, TIME, ZONE, VALUES) TEMP = DATE(7:8) // ' ' // MONTH(VALUES(2)) // ' ' // & DATE(1:4) // ' ' // TIME(1:2) // ':' // TIME(3:4) // & ':' // TIME(5:6) DAYTIME = TEMP RETURN END C===================================================================== SUBROUTINE VDATE(ADATE) C C This routine returns the date as a 9 charater string. Truncation C occurs if the length of the input string is less than 9. The format C of the returned string is "dd-mmm-yy". C C Note: this routine is not Y2K compliant. C C +++++++++++++++++++++++++++++++++++++++ C Machine dependent intel/Win32 version. C +++++++++++++++++++++++++++++++++++++++ C IMPLICIT NONE C input/output CHARACTER*(*) ADATE C local CHARACTER TEMP*20 C CALL VFDATE(TEMP) ADATE = TEMP(1:2) // '-' // TEMP(4:6) // '-' // TEMP(10:11) RETURN END C===================================================================== SUBROUTINE VTIME(ATIME) C C This routine returns the time of day as an 8 character string. C Truncation occurs if the length of the input string is less than 8. C The format of the returned string is "hh:mm:ss". C C +++++++++++++++++++++++++++++++++++++++ C Machine dependent intel/Win32 version. C +++++++++++++++++++++++++++++++++++++++ C IMPLICIT NONE C input/output CHARACTER*(*) ATIME C local CHARACTER TEMP*20 C CALL VFDATE(TEMP) ATIME = TEMP(13:20) RETURN END C===================================================================== C----------------------------------------------------------------------- SUBROUTINE WRITE1 (IFILE,IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE) WRITE (IFILE,'(1X,3I5,4(E12.4,E10.2))') & IH,IK,IL,FSQ,SIGFSQ,F,SIGF,E,SIGE RETURN END C-----------------------------------------------------------------------