PROGRAM BAYES C C ROBERT H. BLESSING C HAUPTMAN-WOODWARD INSTITUTE C 73 HIGH STREET C BUFFALO, NEW YORK 14203, USA C TELEPHONE: (716) 856-9600, EXTENSION 335 C E-MAIL: blessing@hwi.buffalo.edu C C APRIL 1986,..., JANUARY 2001 C C GIVEN A SET OF MEASURED INTENSITY DATA, H, K, L, FSQ, AND SIGMA(FSQ), C THE PROGRAM USES A BAYESIAN PROBABILISTIC TREATMENT TO ESTIMATE C STATISTICAL EXPECTATION VALUES FOR FOBS**2, SIGMA(FOBS**2), FOBS, AND C SIGMA(FOBS). C C THE SUBROUTINES NEGIA AND NEGIC, WHICH ESTIMATE THE BAYESIAN POSTERIOR C MEANS AND VARIANCES, ARE FROM SUPPLEMENTARY PUBLICATION NO. SUP 33352 C FROM: C C SIMON FRENCH AND KEITH WILSON (1979). ON THE TREATMENT OF NEGATIVE C INTENSITY OBSERVATIONS. ACTA CRYST. A34, 517-525. C CHARACTER ATIME*8,ADATE*9,TITLE*80,FILEI*80,FILEO*80,FORM*11 CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5 DIMENSION CELL(6),GINV(3,3) 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 IO3=30 ILP=60 C C READ CONTROL DATA. C OPEN (UNIT=IO1,FILE='bayes.dat',STATUS='OLD') READ (IO1,'(A)') TITLE READ (IO1,*) ITYPE READ (IO1,'(A)') FILEI READ (IO1,'(A)') FILEO READ (IO1,'(A)') LATT READ (IO1,'(A)') PTGP READ (IO1,*) CELL NLOCAL=0 NOBAYES=0 READ (IO1,*,END=5) NLOCAL READ (IO1,*,END=5) NOBAYES 5 CLOSE (UNIT=IO1,STATUS='KEEP') CALL METRIC (CELL,GINV) CALL DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP,ILAUE) C C ECHO CONTROL DATA. C OPEN (UNIT=ILP,FILE='bayes.lp',STATUS='NEW') WRITE (ILP,6000) ATIME,ADATE,TITLE 6000 FORMAT (/1H1,130('-')/1X,'PROGRAM BAYES. ',A,', ',A,'. ',A) WRITE (ILP,'(/1X, &''INPUT REFLECTION FILE TYPE: ITYPE = '',I2)') ITYPE WRITE (ILP,'(/1X, &''INPUT REFLECTION FILE NAME: '',A//1X, &''OUTPUT REFLECTION FILE NAME: '',A)') & FILEI,FILEO WRITE (ILP,'(/1X,A//1X,A,''-LATTICE''//1X,''LAUE GROUP '',A//1X, &''POINT GROUP '',A)') SYST,LATT,LAUE,PTGP IF (PTGP.EQ.LAUE) THEN WRITE (ILP,'(/1X,''CENTROSYMMETRIC STRUCTURE'')') ELSE WRITE (ILP,'(/1X,''NONCENTROSYMMETRIC STRUCTURE'')') END IF WRITE (ILP,'(/1X, &''LATTICE PARAMETERS''/1X, &'' A B C ALPHA BETA GAMMA''/ &1X,3F10.4,3F10.3)') CELL IF (NOBAYES.NE.0) WRITE (ILP,'(//1X, &''BAYESIAN REPLACEMENT WILL NOT BE DONE. INSTEAD:''/1X, &'' FSQ = MAX(FSQ, 0.0)''/1X, &'' SIGMA_FSQ = MAX(SIGMA_FSQ, 0.0)''/1X, &'' F = SQRT(FSQ)''/1X, &'' IF (FSQ .GE. SIGMA_FSQ .AND. SIGMA_FSQ .GT. 0) THEN''/1X, &'' SIGMA_F = SIGMA_FSQ/(2*F)''/1X, &'' ELSE''/1X, &'' SIGMA_F = 0.5*SQRT(SIGMA_FSQ)''/1X, &'' END IF''/1X, &'' IF (SIGMA_FSQ .GT. 0 .AND. SIGMA_F .GT. 0) THEN''/1X, &'' IREJECT = 0''/1X, &'' ELSE''/1X, &'' IREJECT = 1''/1X, &'' END IF''/)') C C COPY DATA FROM THE INPUT REFLECTION DATA FILE TO A DIRECT ACCESS C WORKING DATA FILE THAT WILL ACCOMMODATE 15-WORD, 60-BYTE RECORDS: C C h, k, l, Fsq, sigma(Fsq), F, sigma(F), E, sigma(E), C icent, mult, epsilon, sin(theta)/lambda, , sigma() C IF (ITYPE.LT.2) THEN FORM='FORMATTED' ELSE FORM='UNFORMATTED' END IF OPEN (UNIT=IO1,FILE=FILEI,STATUS='OLD',FORM=FORM) OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='UNFORMATTED', & ACCESS='DIRECT',RECL=4*15) IF (IPTGP.LT.ILAUE) THEN JMIN=+999 JMAX=-999 KMIN=+999 KMAX=-999 LMIN=+999 LMAX=-999 END IF NDATA=0 IEND=0 1 CALL READ1 (ITYPE,IO1,IEND,J,K,L,FSQ,SIGFSQ) IF (IEND.NE.0) GO TO 9 NDATA=NDATA+1 WRITE (IO2,REC=NDATA) J,K,L,FSQ,SIGFSQ, & 0.0, 0.0, 0.0, 0.0, 0, 0, 0, 0.0, 0.0, 0.0 IF (IPTGP.LT.ILAUE) THEN C C STORE THE LAUE-UNIQUE INDEX LIMITS. C CALL EQUIV (ILAUE,J,K,L) JMIN=MIN(JMIN,J) JMAX=MAX(JMAX,J) KMIN=MIN(KMIN,K) KMAX=MAX(KMAX,K) LMIN=MIN(LMIN,L) LMAX=MAX(LMAX,L) END IF GO TO 1 9 CLOSE (UNIT=IO1,STATUS='KEEP') IF (IPTGP.LT.ILAUE) THEN C C STORE A TEMPORARY ARRAY OF THE BIJVOET- OR FRIEDEL-PAIR POPULATIONS C (0, 1, OR 2) INDEXED BY THE PACKED LAUE-UNIQUE INDICES. C NJ=JMAX-JMIN+1 NK=KMAX-KMIN+1 NL=LMAX-LMIN+1 NMAX=NJ*NK*NL+NK*NL+NL C C ALLOCATE STORAGE FOR SORTING ARRAYS. C ALLOCATE (DATA(NMAX),INDEX(NMAX)) DO I=1,NMAX DATA(I)=0 END DO DO I=1,NDATA READ (IO2,REC=I) J,K,L CALL EQUIV (ILAUE,J,K,L) JKL=(J-JMIN)*NK*NL+(K-KMIN)*NL+(L-LMIN) DATA(JKL)=DATA(JKL)+1 END DO END IF C C EVALUATE ADDITIONAL REFLECTION INFORMATION AND STORE IT IN THE WORKING C DATA FILE. C SMIN=9 SMAX=0 DO I=1,NDATA READ (IO2,REC=I) J,K,L C C GET ACENTRIC OR CENTRIC REFLECTION FLAG. C CALL CSYM (IPTGP,J,K,L,IC) C C GET POINT GROUP DEGENERACY AND MULTIPLICITY. C CALL ESYM (IPTGP,J,K,L,IE,M) C C ADJUST DEGENERACY FOR LATTICE MULTIPLICITY. C IE=MLATT*IE C C DOUBLE THE POINT GROUP MULTIPLICITY FOR ACENTRIC REFLECTIONS THAT ARE C REPRESENTED BY ONLY ONE MEMBER OF THE BIJVOET OR FRIEDEL PAIR. C IF (IC.EQ.0) THEN CALL EQUIV (ILAUE,J,K,L) JKL=(J-JMIN)*NK*NL+(K-KMIN)*NL+(L-LMIN) IF (DATA(JKL).EQ.1) M=2*M END IF C C GET SIN(THETA)/LAMBDA. C S=SINTHL(J,K,L,GINV) SMIN=MIN(SMIN,S) SMAX=MAX(SMAX,S) C C RECORD THE ADDITIONAL REFLECTION INFORMATION IN THE WORKING DATA FILE. C READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ WRITE (IO2,REC=I) J,K,L,FSQ,SIGFSQ, & 0.0, 0.0, 0.0, 0.0, IC,M,IE,S, 0.0, 0.0 END DO C C PRINT SOME DATA SET STATISTICS. C WRITE (ILP,'(/1X, &''N = '',I6,'' DATA FOR EVALUATING INTENSITY AVERAGES IN SHE'', &''LLS OF S = SIN(THETA)/LAMBDA''/1X, &''SMIN = '',F6.3,'' RECIPROCAL ANGSTROMS DMAX = '',F7.2, &'' ANGSTROMS''/1X, &''SMAX = '',F6.3,'' DMIN = '',F7.2)') & NDATA,SMIN,1/(2*SMIN),SMAX,1/(2*SMAX) C C ALLOCATE STORAGE FOR SORTING ARRAYS. C IF (IPTGP.LT.ILAUE) DEALLOCATE (DATA,INDEX) ALLOCATE (DATA(NDATA),INDEX(NDATA)) C C SORT ON INCREASING SIN(THETA)/LAMBDA. C DO I=1,NDATA READ (IO2,REC=I) (IDUMMY,J=1,3),(DUMMY,J=1,6),(IDUMMY,J=1,3),S DATA(I)=S END DO CALL SORT (NDATA,DATA,INDEX) C C EVALUATE LOCAL AVERAGE INTENSITIES IN SHELLS OF SIN(THETA)/LAMBDA C CONTAINING EQUAL NUMBERS OF REFLECTIONS. C C AT LEAST 10, AT MOST 500, REFLECTIONS PER SHELL. C IF (NLOCAL.EQ.0) THEN NLOCAL=MAX(NDATA/100,10) NLOCAL=MIN(NLOCAL,500,NDATA) END IF NHALF=NLOCAL/2 SUMM=0 SUM1=0 SUM2=0 YMAX=0 ZMAX=0 DO I=1,NHALF READ (IO2,REC=INDEX(I)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5), & IDUMMY,M,IE IF (FSQ.LT.0) FSQ=0 FSQ=FSQ/IE SUMM=SUMM+M SUM1=SUM1+M*FSQ SUM2=SUM2+M*FSQ**2 END DO DO I=1,NDATA I1=I-NHALF I2=I+NHALF IF (I1.LT.1) THEN I1=1 M1=0 Y1=0 ELSE READ (IO2,REC=INDEX(I1)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5), & IDUMMY,M,IE IF (FSQ.LT.0) FSQ=0 M1=M Y1=FSQ/IE END IF IF (I2.GT.NDATA) THEN I2=NDATA M2=0 Y2=0 ELSE READ (IO2,REC=INDEX(I2)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5), & IDUMMY,M,IE IF (FSQ.LT.0) FSQ=0 M2=M Y2=FSQ/IE END IF SUMM=SUMM-M1+M2 SUM1=SUM1-M1*Y1+M2*Y2 SUM2=SUM2-M1*Y1**2+M2*Y2**2 N=I2-I1+1 AVG=SUM1/SUMM AVGSQ=SUM2/SUMM RMSD=SQRT((SUMM/(SUMM-1))*(AVGSQ-AVG**2)) AVGFSQ=AVG SIGAVG=RMSD/SQRT(SUMM) READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ, & DUMMY,DUMMY,DUMMY,DUMMY,IC,M,IE,S WRITE (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ, & 0.0, 0.0, 0.0, 0.0, IC,M,IE,S,AVGFSQ,SIGAVG YMAX=MAX(YMAX,AVGFSQ) ZMAX=MAX(ZMAX,SIGAVG) END DO C C PRINT THE FIRST 100 REFLECTIONS SORTED ON SIN(THETA)/LAMBDA. C WRITE (ILP,6000) ATIME,ADATE,TITLE WRITE (ILP,'(/1X, &'' H K L FMEAS**2 SIGMA(FMEAS**2) SIN(TH)/L D(HKL) '', &'' SIGMA()''/1X, &'' - - - -------- --------------- --------- ------ '', &'' -------- ---------------'')') DO I=1,MIN(100,NDATA) READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,(DUMMY,N=1,4), & (IDUMMY,N=1,3),S,AVGFSQ,SIGAVG D=1/(2*S) WRITE (ILP,'(1X,3I4,2F12.2,4X,F10.4,F10.2,2F12.2)') & J,K,L,FSQ,SIGFSQ,S,D,AVGFSQ,SIGAVG END DO C C PLOT VERSUS DSTAR AND SIGMA() VERSUS DSTAR. C WRITE (ILP,6001) ATIME,ADATE,TITLE 6001 FORMAT (/1H1,130('-')/1X,10X,'PROGRAM BAYES. ',A,', ',A,'. ',A) CALL PLOT (ILP,IO2,NDATA,0.0,SMAX,0.0,1.1*YMAX,0) WRITE (ILP,'(/1X, &''INTENSITIES LOCALLY AVERAGED IN SHELLS OF S = SIN(THETA)/LAMB'', &''DA. Y = VERSUS X = DSTAR = 1/D = 2*S,''/1X, &''WHERE EPSILON = EPSILON(HKL) IS THE DEGENERACY OF THE RECIPRO'', &''CAL LATTICE POINT HKL.''/1X, &''NTOTAL = '',I6,'' TOTAL DATA SMAX = '', &F10.3,'' RECIPROCAL ANGSTROMS''/1X, &''NLOCAL = '',I6,'' LOCAL DATA PER SHELL DMIN = '', &F9.2,'' ANGSTROMS'')') & NDATA,SMAX,NLOCAL,1/(2*SMAX) WRITE (ILP,6001) ATIME,ADATE,TITLE CALL PLOT (ILP,IO2,NDATA,0.0,SMAX,0.0,1.1*ZMAX,1) WRITE (ILP,'(/1X, &''STANDARD UNCERTAINTIES OF THE LOCAL INTENSITY AVERAGES'',/1X, &''Y = SIGMA() VERSUS X = DSTAR = 2*SIN(THETA)'', &''/LAMBDA'')') C C TEST VALIDITY OF THE WILSON DISTRIBUTIONS. C CALL NZPLOT (NDATA,IO2,ILP,ATIME,ADATE,TITLE,PTGP,LAUE) C C OBTAIN BAYESIAN ESTIMATES OF POSTERIOR MEANS AND VARIANCES. C WRITE (ILP,6000) ATIME,ADATE,TITLE WRITE (ILP,6010) 6010 FORMAT (/1X, &' A PR', &'IORI A POSTERIORI'/1X, &' MEASUREMENT EXPECT', &'ATION BAYESIAN ESTIMATES'/1X, &' ----------------------- ------', &'----- -------------------------------------------'/1X, &' H K L F**2 SIGMA(F**2) SIN(TH)/L D(HKL) e* F**2 SIGMA(F**2) F SIGMA(F)'/1X, &' - - - ----------- ----------- --------- ------ ------', &'----- ----------- ----------- --------- ---------') 6011 FORMAT (1X,3I4,2F12.2,F10.4,F10.2,F12.2,2X,2F12.2,2F10.2) NN=0 NI=0 NJ=0 NF=0 N1REJ=0 N2REJ=0 N3REJ=0 OPEN (UNIT=IO3,STATUS='SCRATCH',FORM='UNFORMATTED', & ACCESS='DIRECT',RECL=4*11) IREC=0 DO 10 I=1,NDATA READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,DUMMY,DUMMY,DUMMY,DUMMY, & IC,M,IE,S,AVGFSQ,SIGAVG XI=FSQ SDI=SIGFSQ SIGMA=IE*AVGFSQ IF (NOBAYES.NE.0) THEN C C OPTIONALLY, DO NOT DO THE BAYESIAN REPLACEMENT. C IF (FSQ.LT.0) FSQ=0 IF (SIGFSQ.LT.0) SIGFSQ=0 F=SQRT(FSQ) IF (FSQ.GE.SIGFSQ.AND.SIGFSQ.GT.0) THEN SIGF=SIGFSQ/(2*F) ELSE SIGF=0.5*SQRT(SIGFSQ) END IF IF (SIGFSQ.GT.0.AND.SIGF.GT.0) THEN XJ=FSQ SDJ=SIGFSQ XF=F SDF=SIGF IREJ=0 ELSE IREJ=1 END IF ELSE C C CALL FRENCH AND WILSON SUBROUTINES TO GET BAYESIAN ESTIMATES. C IREJ=0 IF (IC.EQ.0) CALL NEGIA (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) IF (IC.EQ.1) CALL NEGIC (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) END IF C C ACCUMULATE STATISTICS FOR REJECTED REFLECTIONS. C IF (IREJ.NE.0) THEN IF (IREJ.EQ.+1) N1REJ=N1REJ+1 IF (IREJ.EQ.-1) N2REJ=N2REJ+1 GO TO 10 END IF C C REJECT REFLECTIONS WITH UNREASONABLY LARGE CHANGES IN LARGER-THAN- C AVERAGE FSQ. C IF (XJ.GT.SIGMA.AND.ABS(XJ-XI).GT.2*SDI) THEN N3REJ=N3REJ+1 WRITE (IO3,REC=N3REJ) J,K,L,XI,SDI,S,SIGMA,XJ,SDJ,XF,SDF GO TO 10 END IF C C RECORD THE ACCEPTED REFLECTIONS IN THE WORKING DATA FILE. C FSQ=XJ SIGFSQ=SDJ F=XF SIGF=SDF IREC=IREC+1 WRITE (IO2,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF, 0.0, 0.0, & IC,M,IE,S,AVGFSQ,SIGAVG C C PRINT SELECTED REFELCTIONS. C IF (IPRINT(J,K,L).EQ.1) WRITE (ILP,6011) J,K,L,XI,SDI,S,1/(2*S), & SIGMA,XJ,SDJ,XF,SDF C C COMPILE SOME DATA SET STATISTICS. C IF (XI.LT.0) NN=NN+1 IF (XI.LT.3*SDI) NI=NI+1 IF (XJ.LT.3*SDJ) NJ=NJ+1 IF (XF.LT.3*SDF) NF=NF+1 10 CONTINUE NDATA=NDATA-N1REJ-N2REJ-N3REJ C C PRINT DATA SET STATISTICS. C WRITE (ILP,6000) ATIME,ADATE,TITLE WRITE (ILP,'(/1X, &''N1 = '',I6 ,'' REFLECTIONS REJECTED BECAUSE MEASURED FSQ/SIG'', &''MA(FSQ) .LT. -4''/1X, &'' (UNREASONABLY NEGATIVE INTENSITY MEASUREMENTS)'')') & N1REJ WRITE (ILP,'(/1X, &''N2 = '',I6,'' REFLECTIONS REJECTED BECAUSE MEASURED FSQ/SIGM'', &''A(FSQ) - Q*SIGMA(FSQ)/ .LT. -4''/1X, &'' WHERE Q = 1 FOR ACENTRIC REFLECTIONS AND Q = 0.5'', &'' FOR CENTRIC REFLECTIONS''/1X, &'' (UNREASONABLY LARGE EXPERIMENTAL ERROR ESITMATES'', &'')'')') N2REJ WRITE (ILP,'(/1X, &''N3 = '',I6,'' REFLECTIONS REJECTED BECAUSE ABS(DELTA(FSQ)) .'', &''GT. 2*SIGMA(FSQ) .AND. FSQ .GT. ''/1X, &'' (UNREASONABLY LARGE CHANGES IN LARGER-THAN-LOCAL'', &''-AVERAGE INTENSITIES)'')') N3REJ WRITE (ILP,'(/1X, &''NP = '',I6,'' DATA PROCESSED''/1X, &''NN = '',I6,'' INPUT DATA WITH FSQ .LT. 0''/1X, &''NI = '',I6,'' INPUT DATA WITH FSQ .LT. 3*SIGMA(FSQ)''/1X, &''NJ = '',I6,'' OUTPUT DATA WITH FSQ .LT. 3*SIGMA(FSQ)''/1X, &''NF = '',I6,'' OUTPUT DATA WITH F .LT. 3*SIGMA(F)'')') & NDATA,NN,NI,NJ,NF IF (N3REJ.GT.0) THEN DO I=1,N3REJ READ (IO3,REC=I) J,K,L,XI,SDI,S,SIGMA,XJ,SDJ,XF,SDF DATA(I)=XJ/SIGMA END DO IF (N3REJ.GT.1) THEN CALL SORT (N3REJ,DATA,INDEX) ELSE INDEX(1)=1 END IF WRITE (ILP,'(//1X, &''REJECTED LARGER-THAN-LOCAL-AVERAGE INTENSITIES, LISTED IN ORD'', &''ER OF DECREASING E = SQRT[FSQ/(EPSILON*)]:'')') WRITE (ILP,6010) DO I=1,MIN(N3REJ,50) J=N3REJ-I+1 READ (IO3,REC=INDEX(J)) J,K,L,XI,SDI,S, SIGMA, & XJ,SDJ,XF,SDF WRITE (ILP,6011) J,K,L,XI,SDI,S,1/(2*S),SIGMA, & XJ,SDJ,XF,SDF END DO CLOSE (UNIT=IO3,STATUS='DELETE') END IF C C ESTIMATE LOCALLY NORMALIZED STRUCTURE FACTORS, C E(HKL) = F(HKL)/SQRT[EPSILON(HKL)*]. C CALL NORMAL (NDATA,IO2,IO1,ILP,ATIME,ADATE,TITLE,SMIN,SMAX,GINV, & DATA,INDEX,NLOCAL,NOBAYES) C C WRITE OUTPUT REFLECTION FILE. C OPEN (UNIT=IO1,FILE=FILEO,STATUS='NEW',FORM='FORMATTED') DO I=1,NDATA READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE IF (E.GT.0) CALL WRITE1 (IO1,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) END DO STOP 'PROGRAM BAYES 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,ILAUE) 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 I=1,42 IF (PTGP.EQ.SYMBOL(I)) GO TO 1 END DO STOP 'BAYES/DECODE: ERROR IN POINT GROUP SYMBOL' 1 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 ' DO I=1,42 IF (LAUE.EQ.SYMBOL(I)) GO TO 2 END DO 2 ILAUE=I 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----------------------------------------------------------------------- SUBROUTINE EQUIV (PTGP,H,K,L) C C EQUIVALENT, UNIQUE REFLECTION INDICES FOR THE CRYSTALLOGRAPHIC POINT C 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 SEE INTERNATIONAL TABLES FOR CRYSTALLOGRAPHY (1983). VOLUME A, PP. C 750-770. 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,X,Y,Z 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 C C 1 C 1 RETURN C C -1 C 2 IF (L.GT.0) RETURN H=-H K=-K L=-L IF (L.GT.0) RETURN IF (K.GT.0) RETURN H=-H K=-K IF (K.GT.0) RETURN IF (H.GE.0) RETURN H=-H RETURN C C 2 C 3 IF (L.GT.0) RETURN H=-H L=-L IF (L.GT.0) RETURN IF (H.GE.0) RETURN H=-H RETURN C C M C 4 K=ABS(K) RETURN C C 2/M C 5 K=ABS(K) GO TO 3 C C 222 C 6 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) THEN H=ABS(H) K=ABS(K) L=ABS(L) RETURN END IF IF (H.GE.0.AND.K.GE.0) RETURN X=H Y=K Z=L H=-X K=-Y L=Z IF (H.GE.0.AND.K.GE.0) RETURN H=-X K=Y L=-Z IF (H.GE.0.AND.K.GE.0) RETURN H=X K=-Y L=-Z RETURN C C MM2 C 7 H=ABS(H) K=ABS(K) RETURN C C MMM C 8 L=ABS(L) GO TO 7 C C 4 C 9 IF (H.EQ.0.AND.K.EQ.0) RETURN IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K H=-Y K=X IF (H.GT.0.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=-X RETURN C C -4 C 10 IF (H.EQ.0.AND.K.EQ.0) THEN L=ABS(L) RETURN END IF IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K H=-X K=-Y IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=-X L=-L IF (H.GT.0.AND.K.GE.0) RETURN H=-Y K=X RETURN C C 4/M C 11 L=ABS(L) GO TO 9 C C 422 C 12 IF (H.EQ.0.OR.K.EQ.0.OR.ABS(H).EQ.ABS(K)) L=ABS(L) 512 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=X IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=-X IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=Y L=-L GO TO 512 C C 4MM C 13 H=ABS(H) K=ABS(K) IF (H.GE.K) RETURN X=H H=K K=X RETURN C C -42M C 14 IF (H.EQ.0.OR.K.EQ.0) L=ABS(L) 514 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=-X L=-L IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=X IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=Y GO TO 514 C C -4M2 C 15 H=ABS(H) K=ABS(K) IF (H.EQ.K) L=ABS(L) IF (H.GE.K) RETURN X=H H=K K=X L=-L RETURN C C 4/MMM C 16 L=ABS(L) GO TO 13 C C 3 C 17 IF (H.EQ.0.AND.K.EQ.0) RETURN IF (H.LT.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.LT.0.AND.K.GE.0) RETURN H=Y K=I RETURN C C -3 C 18 IF (H.EQ.0.AND.K.EQ.0) THEN L=ABS(L) RETURN END IF 518 IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN X=H Y=K I=-H-K H=I K=X IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN H=Y K=I IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN H=-X K=-Y L=-L GO TO 518 C C 312 C 19 IF (H.EQ.0.OR.K.EQ.0.OR.-H.EQ.K) L=ABS(L) 519 IF (H.GE.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.0.AND.K.GE.0) RETURN H=-Y K=-X L=-L GO TO 519 C C 321 C 20 IF (H.EQ.K.OR.-2*H.EQ.K.OR.-H.EQ.2*K) L=ABS(L) 520 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=X L=-L GO TO 520 C C 31M C 21 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=X GO TO 21 C C 3M1 C 22 IF (H.GE.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.0.AND.K.GE.0) RETURN H=-Y K=-X GO TO 22 C C -31M C 23 IF (H.EQ.0.OR.K.EQ.0.OR.-H.EQ.K) L=ABS(L) 523 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=Y K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=Y K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=X K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=I K=Y IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=-X K=-Y L=-L GO TO 523 C C -3M1 C 24 IF (H.EQ.K.OR.-2*H.EQ.K.OR.-H.EQ.2*K) L=ABS(L) 524 IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=Y K=I IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=Y K=X L=-L IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=X K=I IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=I K=Y IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=-X K=-Y GO TO 524 C C 6 C 25 IF (H.EQ.0.AND.K.EQ.0) RETURN IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GT.0.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GT.0.AND.K.GE.0) RETURN H=-I K=-X IF (H.GT.0.AND.K.GE.0) RETURN H=-Y K=-I RETURN C C -6 C 26 L=ABS(L) GO TO 17 C C 6/M C 27 L=ABS(L) GO TO 25 C C 622 C 28 IF (H.EQ.0.OR.K.EQ.0.OR.ABS(H).EQ.ABS(K).OR.ABS(2*H).EQ.ABS(K).OR. & ABS(H).EQ.ABS(2*K)) L=ABS(L) 528 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=-I K=-X IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=-I IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=X L=-L GO TO 528 C C 6MM C 29 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=-I K=-X IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=-I IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=X GO TO 29 C C -6M2 C 30 L=ABS(L) IF (H.EQ.0.AND.K.EQ.0) RETURN 530 IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GT.0.AND.K.GE.0) RETURN H=-Y K=-X GO TO 530 C C -62M C 31 L=ABS(L) IF (H.EQ.0.AND.K.EQ.0) RETURN 531 IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=I IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=X GO TO 531 C C 6/MMM C 32 L=ABS(L) GO TO 29 C C R3 C 33 STOP 'BAYES/EQUIV: RE-INDEX ON HEXAGONAL AXES' C C R-3 C 34 GO TO 33 C C R32 C 35 GO TO 33 C C R3M C 36 GO TO 33 C C R-3M C 37 GO TO 33 C C 23 C 38 IF (H.EQ.K.AND.K.EQ.ABS(L)) RETURN 538 IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN X=H Y=K Z=L H=-X K=-Y L=Z IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN H=-X K=Y L=-Z IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN H=X K=-Y L=-Z IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN H=Z K=X L=Y GO TO 538 C C M3 C 39 H=ABS(H) K=ABS(K) L=ABS(L) IF (H.EQ.K.AND.K.EQ.L) RETURN 539 IF (H.GT.L.AND.K.GE.L) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GT.L.AND.K.GE.L) RETURN H=Y K=Z L=X GO TO 539 C C 432 C 40 IF (ABS(H).EQ.ABS(K).AND.ABS(K).EQ.ABS(L)) THEN H=ABS(H) K=ABS(K) L=ABS(L) RETURN END IF 540 IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Y K=Z L=X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-X K=-Y L=Z IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Z K=-X L=-Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Y K=Z L=-X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-X K=Y L=-Z IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Z K=-X L=Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Y K=-Z L=-X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=X K=-Y L=-Z IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Z K=X L=-Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Y K=-Z L=X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Y K=X L=-Z GO TO 540 C C -43M C 41 IF (H.GE.K.AND.K.GE.ABS(L)) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Y K=Z L=X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-X K=-Y L=Z IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Z K=-X L=-Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Y K=Z L=-X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-X K=Y L=-Z IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Z K=-X L=Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Y K=-Z L=-X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=X K=-Y L=-Z IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Z K=X L=-Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Y K=-Z L=X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Y K=X L=Z GO TO 41 C C M3M C 42 H=ABS(H) K=ABS(K) L=ABS(L) IF (H.GE.K.AND.K.GE.L) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GE.K.AND.K.GE.L) RETURN H=Y K=Z L=X IF (H.GE.K.AND.K.GE.L) RETURN H=Z K=Y L=X IF (H.GE.K.AND.K.GE.L) RETURN H=X K=Z L=Y IF (H.GE.K.AND.K.GE.L) RETURN H=Y K=X L=Z RETURN END C------------------------------------------------------------------------ FUNCTION ERRF(X) C C RETURNS THE ERROR FUNCTION C C ERF(X) = [2/SQRT(PI)]*INTEGRAL[0,ABS(X)] EXP(-T**2) DT C C RATIONAL APPROXIMATION FORMULA NO. 7.1.26 FROM: C C MILTON ABRAMOWITZ AND IRENE A. STEGUN (1965). HANDBOOK OF C MATHEMATICAL FUNCTIONS, P. 299. NEW YORK: DOVER PUBLICATIONS, INC. C REAL*8 Z,ONE,T,P,A1,A2,A3,A4,A5 DATA ONE,P,A1,A2,A3,A4,A5 /1.D0, 0.3275911D0, 0.254829592D0, & -0.284496736D0, 1.421413741D0, -1.43152027D0, 1.061405429D0/ Z=ABS(X) T=ONE/(ONE+P*Z) T=T*(A1+T*(A2+T*(A3+T*(A4+T*A5)))) ERRF=ONE-T*EXP(-Z**2) IF (X.LT.0) ERRF=-ERRF RETURN END C----------------------------------------------------------------------- SUBROUTINE ESTATS (NREC,IFILE,ILP,ATIME,ADATE,TITLE,SMIN,SMAX) 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) CHARACTER ATIME*8,ADATE*9,TITLE*80 WRITE (ILP,'(/1H1,130(''-'')/1X, &''SUBPROGRAM NORMAL/BAYES. '',A,'', '',A,''. '',A)') & ATIME,ADATE,TITLE WRITE (ILP,'(/1X, &''E-STATISTICS FOR NORMALIZED STRUCTURE FACTORS ESTIMATED AS E('', &''HKL) = F(HKL)/SQRT[EPSILON(HKL)*]'')') C C MOMENTS OF THE E-DISTRIBUTION C NA=0 NC=0 DO IREC=1,NREC READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (SMIN.LE.S.AND.S.LE.SMAX.AND.E.GT.0) THEN DO 1 I=1,7 IF (I.EQ.2.AND.IC.NE.0) GO TO 1 IF (I.EQ.3.AND.IC.NE.1) GO TO 1 IF (I.EQ.4.AND.(J.EQ.0.OR.K.EQ.0.OR.L.EQ.0)) GO TO 1 IF (I.EQ.5.AND.J.NE.0) GO TO 1 IF (I.EQ.6.AND.K.NE.0) GO TO 1 IF (I.EQ.7.AND.L.NE.0) GO TO 1 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 1 CONTINUE IF (IC.EQ.0) NA=NA+1 IF (IC.EQ.1) NC=NC+1 END IF END DO 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, &' -----------------------------------------', &'--------------------------- -------------------------------'/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 DSTAR = 2*SIN(THETA)/LAMBDA C N=10 DO I=1,N SI(I)=(FLOAT(I)/N)**(1/3.0)*SMAX END DO DO 2 IREC=1,NREC READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 2 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 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 DO IREC=1,NREC READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (SMIN.LE.S.AND.S.LE.SMAX.AND.E.GT.0) THEN 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, &' ------ ------ --- ---------- ---------- ----------', &'----'/1X, &' EEE',2F10.3,2F10.2,I10/1X, &' EEO',2F10.3,2F10.2,I10/1X, &' EOE',2F10.3,2F10.2,I10/1X, &' OEE',2F10.3,2F10.2,I10/1X, &' EOO',2F10.3,2F10.2,I10/1X, &' OEO',2F10.3,2F10.2,I10/1X, &' OOE',2F10.3,2F10.2,I10/1X, &' OOO',2F10.3,2F10.2,I10) C C PLOT 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 DO 10 IREC=1,NREC READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 10 IF (IC.EQ.1) GO TO 10 DO I=1,N IF (E.LE.I*DX) THEN X(I)=X(I)+M*E Y(I)=Y(I)+M GO TO 10 END IF END DO 10 CONTINUE DO I=1,N IF (Y(I).GT.0) X(I)=X(I)/Y(I) END DO 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,6000) ATIME,ADATE,TITLE CALL PLOTA (N,X,Y,ILP) END IF 6000 FORMAT (/1H1,130('-')/1X,10X, &'SUBPROGRAM NORMAL/BAYES. ',A,', ',A,'. ',A/) 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 DO 11 IREC=1,NREC READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 11 IF (IC.EQ.0) GO TO 11 DO I=1,N IF (E.LE.I*DX) THEN X(I)=X(I)+M*E Y(I)=Y(I)+M GO TO 11 END IF END DO 11 CONTINUE DO I=1,N IF (Y(I).GT.0) X(I)=X(I)/Y(I) END DO 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,6000) ATIME,ADATE,TITLE CALL PLOTC (N,X,Y,ILP) END IF 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 DO 3 IREC=1,NREC READ (IFILE,REC=IREC) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (S.LT.SMIN.OR.S.GT.SMAX.OR.E.LE.0) GO TO 3 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 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 WRITE (ILP,6000) ATIME,ADATE,TITLE CALL PLOTS (N,X,Y,SMAX,ILP) 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 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----------------------------------------------------------------------- FUNCTION IPRINT(IH,IK,IL) 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 METRIC (A,GB) C C CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND RECIPROCAL C LATTICE METRIC TENSORS FROM DIRECT LATTICE PARAMETERS. C DIMENSION A(6),B(6),GA(3,3),GB(3,3) PI=ACOS(-1.0) 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 NEGIA (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) C C FROM SUPPLEMENTARY PUBLICATION NO. SUP 33352 FROM: C C SIMON FRENCH AND KEITH WILSON (1979). ON THE TREATMENT OF C NEGATIVE INTENSITY OBSERVATIONS. ACTA CRYST. A34, 517-525. C C GIVEN AN OBSERVED INTENSITY AND ITS ASSOCIATED STANDARD DEVIATION, C THE SUBROUTINE CALCULATES AND RETURNS THE POSTERIOR MEAN AND C STANDARD DEVIATION OF BOTH THE TRUE INTENSITY AND THE TRUE C STRUCTURE FACTOR MODULUS. C C THE PRIOR DISTRIBUTION USED IS WILSON'S DISTRIBUTION FOR AN C ACENTRIC REFLECTION. FOR EXACT FORM SEE FRENCH AND WILSON (1977), C EQUATION (2.2). C C ACCURACY - PROBABLY BETTER THAN 5% C C ARGUMENTS: C C XI - OBSERVED INTENSITY C SDI - S.D. OF OBSERVATION C SIGMA - PARAMETER FOR WILSON'S DISTRIBUTION C XJ - POSTERIOR MEAN FOR THE TRUE INTENSITY C SDJ - POSTERIOR S.D. FOR THE TRUE INTENSITY C XF - POSTERIOR MEAN FOR THE TRUE STRUCTURE FACTOR MODULUS C SDF - POSTERIOR S.D. FOR THE TRUE STRUCTURE FACTOR MODULUS C IREJ - ON ENTRY: = -1 IF ONLY XF AND SDF REQUIRED C = 0 IF ALL XJ, SDJ, XF, AND SDJ REQUIRED C = +1 IF ONLY XJ AND SDJ REQUIRED C ON EXIT: = -1 IF XI/SDI - SDI/SIGMA < -4 C = 0 IF ALL O.K. C = +1 IF XI/SDI < -4, I.E., OBSERVATION IS C ESSENTIALLY IMPOSSIBLE C C ARRAYS ZJ, ZJSD, ZF, ZFSD CONTAIN THE TABULATED VALUES OF XJ, SDJ, C XF, AND SDF FOR XI/SDI - SDI/SIGMA = -4.0, (0.1), 3.0. LINEAR C INTERPOLATION IS USED IN THESE TABLES. C DIMENSION ZJ(71),ZJSD(71),ZF(71),ZFSD(71) DATA ZJ /0.226, 0.230, 0.235, 0.240, 0.246, 0.251, 0.257, & 0.263, 0.270, 0.276, 0.283, 0.290, 0.298, 0.306, & 0.314, 0.323, 0.332, 0.341, 0.351, 0.362, 0.373, & 0.385, 0.397, 0.410, 0.424, 0.439, 0.454, 0.470, & 0.487, 0.505, 0.525, 0.545, 0.567, 0.590, 0.615, & 0.641, 0.668, 0.698, 0.729, 0.762, 0.798, 0.835, & 0.875, 0.917, 0.962, 1.009, 1.059, 1.112, 1.167, & 1.226, 1.287, 1.352, 1.419, 1.490, 1.563, 1.639, & 1.717, 1.798, 1.882, 1.967, 2.055, 2.145, 2.236, & 2.329, 2.422, 2.518, 2.614, 2.710, 2.808, 2.906, & 3.004/ DATA ZJSD /0.217, 0.221, 0.226, 0.230, 0.235, 0.240, 0.245, & 0.250, 0.255, 0.261, 0.267, 0.273, 0.279, 0.286, & 0.292, 0.299, 0.307, 0.314, 0.322, 0.330, 0.339, & 0.348, 0.357, 0.367, 0.377, 0.387, 0.398, 0.409, & 0.421, 0.433, 0.446, 0.459, 0.473, 0.488, 0.503, & 0.518, 0.535, 0.551, 0.568, 0.586, 0.604, 0.622, & 0.641, 0.660, 0.679, 0.698, 0.718, 0.737, 0.757, & 0.776, 0.795, 0.813, 0.831, 0.848, 0.865, 0.881, & 0.895, 0.909, 0.921, 0.933, 0.943, 0.953, 0.961, & 0.968, 0.974, 0.980, 0.984, 0.988, 0.991, 0.994, & 0.996/ DATA ZF /0.423, 0.428, 0.432, 0.437, 0.442, 0.447, 0.453, & 0.458, 0.464, 0.469, 0.475, 0.482, 0.488, 0.495, & 0.502, 0.509, 0.516, 0.524, 0.532, 0.540, 0.549, & 0.557, 0.567, 0.576, 0.586, 0.597, 0.608, 0.619, & 0.631, 0.643, 0.656, 0.670, 0.684, 0.699, 0.714, & 0.730, 0.747, 0.765, 0.783, 0.802, 0.822, 0.843, & 0.865, 0.887, 0.911, 0.935, 0.960, 0.987, 1.014, & 1.042, 1.070, 1.100, 1.130, 1.161, 1.192, 1.224, & 1.257, 1.289, 1.322, 1.355, 1.388, 1.421, 1.454, & 1.487, 1.519, 1.551, 1.583, 1.615, 1.646, 1.676, & 1.706/ DATA ZFSD /0.216, 0.218, 0.220, 0.222, 0.224, 0.226, 0.229, & 0.231, 0.234, 0.236, 0.239, 0.241, 0.244, 0.247, & 0.250, 0.253, 0.256, 0.259, 0.262, 0.266, 0.269, & 0.272, 0.276, 0.279, 0.283, 0.287, 0.291, 0.295, & 0.298, 0.302, 0.307, 0.311, 0.315, 0.319, 0.324, & 0.328, 0.332, 0.337, 0.341, 0.345, 0.349, 0.353, & 0.357, 0.360, 0.364, 0.367, 0.369, 0.372, 0.374, & 0.375, 0.376, 0.377, 0.377, 0.377, 0.376, 0.374, & 0.372, 0.369, 0.366, 0.362, 0.358, 0.353, 0.348, & 0.343, 0.338, 0.332, 0.327, 0.321, 0.315, 0.310, & 0.304/ C IOP=IREJ IREJ=0 C C IS OBSERVATION POSSIBLE? C H=XI/SDI IF (H+4.0) 5,8,8 C 5 CONTINUE IREJ=1 GO TO 15 C C LOCATE RANGE OF DISTRIBUTION. C 8 CONTINUE H=H-SDI/SIGMA IF (H+4.0) 10,20,20 C C DISTRIBUTION IS ESSENTIALLY NEGATIVE. C 10 CONTINUE IREJ=-1 15 CONTINUE XJ=0.0 SDJ=1.0 XF=0.0 SDF=1.0 RETURN C 20 CONTINUE IF (H-3.0) 60,30,30 C C DISTRIBUTION IS ESSENTIALLY ALL POSTIVE. C 30 CONTINUE H=H*SDI IF (IOP) 40,35,35 35 CONTINUE XJ=H SDJ=SDI 40 CONTINUE IF (IOP) 45,45,50 45 CONTINUE XF=SQRT(H) SDF=0.5*SDI/XF 50 CONTINUE RETURN C C LOCATE XI IN THE TABLES. C 60 CONTINUE H=H+4.0 IH=H*10.0 A=(H-IH*0.1)*10.0 B=1.0-A IH=IH+1 IHP1=IH+1 C C INTERPOLATE IN THE TABLES. C IF (IOP) 70,65,65 65 CONTINUE XJ=(A*ZJ(IHP1)+B*ZJ(IH))*SDI SDJ=(A*ZJSD(IHP1)+B*ZJSD(IH))*SDI 70 CONTINUE IF (IOP) 75,75,80 75 CONTINUE SDF=SQRT(SDI) XF=(A*ZF(IHP1)+B*ZF(IH))*SDF SDF=(A*ZFSD(IHP1)+B*ZFSD(IH))*SDF 80 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE NEGIC (XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) C C FROM SUPPLEMENTARY PUBLICATION NO. SUP 33352 FROM: C C SIMON FRENCH AND KEITH WILSON (1979). ON THE TREATMENT OF C NEGATIVE INTENSITY OBSERVATIONS. ACTA CRYST. A34, 517-525. C C GIVEN AN OBSERVED INTENSITY AND ITS ASSOCIATED STANDARD DEVIATION, C THE SUBROUTINE CALCULATES AND RETURNS THE POSTERIOR MEAN AND C STANDARD DEVIATION OF BOTH THE TRUE INTENSITY AND THE TRUE C STRUCTURE FACTOR MODULUS. C C THE PRIOR DISTRIBUTION USED IS WILSON'S DISTRIBUTION FOR A C CENTRIC REFLECTION. FOR EXACT FORM SEE FRENCH AND WILSON (1977), C EQUATION (2.3). C C ACCURACY - PROBABLY BETTER THAN 5% C C ARGUMENTS: C C XI - OBSERVED INTENSITY C SDI - S.D. OF OBSERVATION C SIGMA - PARAMETER FOR WILSON'S DISTRIBUTION C XJ - POSTERIOR MEAN FOR THE TRUE INTENSITY C SDJ - POSTERIOR S.D. FOR THE TRUE INTENSITY C XF - POSTERIOR MEAN FOR THE TRUE STRUCTURE FACTOR MODULUS C SDF - POSTERIOR S.D. FOR THE TRUE STRUCTURE FACTOR MODULUS C IREJ - ON ENTRY: = -1 IF ONLY XF AND SDF REQUIRED C = 0 IF ALL XJ, SDJ, XF, AND SDJ REQUIRED C = +1 IF ONLY XJ AND SDJ REQUIRED C ON EXIT: = -1 IF XI/SDI - SDI/(2*SIGMA) < -4 C = 0 IF ALL O.K. C = +1 IF XI/SDI < -4, I.E., OBSERVATION IS C ESSENTIALLY IMPOSSIBLE C C ARRAYS ZJ, ZJSD, ZF, ZFSD CONTAIN THE TABULATED VALUES OF XJ, SDJ, C XF, AND SDF FOR XI/SDI - SDI/(2*SIGMA) = -4.0, (0.1), 4.0. LINEAR C INTERPOLATION IS USED IN THESE TABLES. C DIMENSION ZJ(81),ZJSD(81),ZF(81),ZFSD(81) DATA ZJ /0.114, 0.116, 0.119, 0.122, 0.124, 0.127, 0.130, 0.134, & 0.137, 0.141, 0.145, 0.148, 0.153, 0.157, 0.162, 0.166, & 0.172, 0.177, 0.183, 0.189, 0.195, 0.202, 0.209, 0.217, & 0.225, 0.234, 0.243, 0.253, 0.263, 0.275, 0.287, 0.300, & 0.314, 0.329, 0.345, 0.363, 0.382, 0.402, 0.425, 0.449, & 0.475, 0.503, 0.534, 0.567, 0.603, 0.642, 0.684, 0.730, & 0.779, 0.833, 0.890, 0.952, 1.018, 1.089, 1.164, 1.244, & 1.327, 1.416, 1.508, 1.603, 1.703, 1.805, 1.909, 2.015, & 2.123, 2.233, 2.343, 2.453, 2.564, 2.674, 2.784, 2.894, & 3.003, 3.112, 3.220, 3.328, 3.435, 3.541, 3.647, 3.753, & 3.962/ DATA ZJSD/ 0.158, 0.161, 0.165, 0.168, 0.172, 0.176, 0.179, 0.184, & 0.188, 0.192, 0.197, 0.202, 0.207, 0.212, 0.218, 0.224, & 0.230, 0.236, 0.243, 0.250, 0.257, 0.265, 0.273, 0.282, & 0.291, 0.300, 0.310, 0.321, 0.332, 0.343, 0.355, 0.368, & 0.382, 0.397, 0.412, 0.428, 0.445, 0.463, 0.481, 0.501, & 0.521, 0.543, 0.565, 0.589, 0.613, 0.638, 0.664, 0.691, & 0.718, 0.745, 0.773, 0.801, 0.828, 0.855, 0.881, 0.906, & 0.929, 0.951, 0.971, 0.989, 1.004, 1.018, 1.029, 1.038, & 1.044, 1.049, 1.052, 1.054, 1.054, 1.053, 1.051, 1.049, & 1.047, 1.044, 1.041, 1.039, 1.036, 1.034, 1.031, 1.029, & 1.028/ DATA ZF /0.269, 0.272, 0.276, 0.279, 0.282, 0.286, 0.289, 0.293, & 0.297, 0.301, 0.305, 0.309, 0.314, 0.318, 0.323, 0.328, & 0.333, 0.339, 0.344, 0.350, 0.356, 0.363, 0.370, 0.377, & 0.384, 0.392, 0.400, 0.409, 0.418, 0.427, 0.438, 0.448, & 0.460, 0.471, 0.484, 0.498, 0.512, 0.527, 0.543, 0.560, & 0.578, 0.597, 0.618, 0.639, 0.662, 0.687, 0.713, 0.740, & 0.769, 0.800, 0.832, 0.866, 0.901, 0.938, 0.976, 1.016, & 1.057, 1.098, 1.140, 1.183, 1.227, 1.270, 1.313, 1.356, & 1.398, 1.439, 1.480, 1.519, 1.558, 1.595, 1.632, 1.667, & 1.701, 1.735, 1.767, 1.799, 1.829, 1.859, 1.889, 1.917, & 1.945/ DATA ZFSD /0.203, 0.205, 0.207, 0.209, 0.211, 0.214, 0.216, 0.219, & 0.222, 0.224, 0.227, 0.230, 0.233, 0.236, 0.239, 0.243, & 0.246, 0.250, 0.253, 0.257, 0.261, 0.265, 0.269, 0.273, & 0.278, 0.283, 0.288, 0.293, 0.298, 0.303, 0.309, 0.314, & 0.320, 0.327, 0.333, 0.340, 0.346, 0.353, 0.361, 0.368, & 0.375, 0.383, 0.390, 0.398, 0.405, 0.413, 0.420, 0.427, & 0.433, 0.440, 0.445, 0.450, 0.454, 0.457, 0.459, 0.460, & 0.460, 0.458, 0.455, 0.451, 0.445, 0.438, 0.431, 0.422, & 0.412, 0.402, 0.392, 0.381, 0.370, 0.360, 0.349, 0.339, & 0.330, 0.321, 0.312, 0.304, 0.297, 0.290, 0.284, 0.278, & 0.272/ C IOP=IREJ IREJ=0 C C IS OBSERVATION POSSIBLE? C H=XI/SDI IF (H+4.0) 5,8,8 C 5 CONTINUE IREJ=1 GO TO 15 C 8 CONTINUE H=H-0.5*SDI/SIGMA IF (H+4.0) 10,20,20 C 10 CONTINUE IREJ=-1 15 CONTINUE XJ=0.0 SDJ=1.0 XF=0.0 SDF=1.0 RETURN C 20 CONTINUE IF (H-4.0) 60,30,30 C C DISTRIBUTION IS ESSENTIALLY ALL POSTIVE. C 30 CONTINUE XI2=1.0/(H*H) XI4=XI2*XI2 XXF=SQRT(H)*(1.0-0.375*XI2-0.6796875*XI4) XSDF=SQRT(H*(0.25*XI2+0.46875*XI4)) IF (IOP) 40,35,35 35 CONTINUE XJ=H*(1.0-0.5*XI2-0.75*XI4) SDJ=SDI*XSDF*2.0*XXF XJ=XJ*SDI 40 CONTINUE IF (IOP) 45,45,50 45 CONTINUE SDF=SQRT(SDI) XF=XXF*SDF SDF=XSDF*SDF 50 CONTINUE RETURN C C LOCATE XI IN THE TABLES. C 60 CONTINUE H=H+4.0 IH=H*10.0 A=(H-IH*0.1)*10.0 B=1.0-A IH=IH+1 IHP1=IH+1 C C INTERPOLATE IN THE TABLES. C IF (IOP) 70,65,65 65 CONTINUE XJ=(A*ZJ(IHP1)+B*ZJ(IH))*SDI SDJ=(A*ZJSD(IHP1)+B*ZJSD(IH))*SDI 70 CONTINUE IF (IOP) 75,75,80 75 CONTINUE SDF=SQRT(SDI) XF=(A*ZF(IHP1)+B*ZF(IH))*SDF SDF=(A*ZFSD(IHP1)+B*ZFSD(IH))*SDF 80 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE NORMAL (NDATA,IO2,IO1,ILP,ATIME,ADATE,TITLE,SMIN,SMAX, & GINV,DATA,INDEX,NLOCAL,NOBAYES) C C LOCALLY NORMALIZED STRUCTURE FACTOR MAGNITUDES C C E(HKL) = F(HKL)/SQRT[EPSILON(HKL)*] C CHARACTER ATIME*8,ADATE*9,TITLE*80,FILEO*80 DIMENSION GINV(3,3),DATA(1),INDEX(1) IF (NOBAYES.NE.0) GO TO 5 C C SORT REFLECTIONS ON INCREASING SIN(THETA)/LAMBDA. C DO I=1,NDATA READ (IO2,REC=I) (IDUMMY,J=1,3),(DUMMY,J=1,6),(IDUMMY,J=1,3),S DATA(I)=S END DO CALL SORT (NDATA,DATA,INDEX) C C EVALUATE LOCAL AVERAGE INTENSITIES IN SHELLS OF SIN(THETA)/LAMBDA C CONTAINING EQUAL NUMBERS OF REFLECTIONS. C IF (NLOCAL.EQ.0) THEN NLOCAL=MAX(NDATA/100,10) NLOCAL=MIN(NLOCAL,500,NDATA) END IF NHALF=NLOCAL/2 SUMM=0 SUM1=0 SUM2=0 DO I=1,NHALF READ (IO2,REC=INDEX(I)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5), & IDUMMY,M,IE IF (FSQ.LT.0) FSQ=0 FSQ=FSQ/IE SUMM=SUMM+M SUM1=SUM1+M*FSQ SUM2=SUM2+M*FSQ**2 END DO DO I=1,NDATA I1=I-NHALF I2=I+NHALF IF (I1.LT.1) THEN I1=1 M1=0 Y1=0 ELSE READ (IO2,REC=INDEX(I1)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5), & IDUMMY,M,IE IF (FSQ.LT.0) FSQ=0 M1=M Y1=FSQ/IE END IF IF (I2.GT.NDATA) THEN I2=NDATA M2=0 Y2=0 ELSE READ (IO2,REC=INDEX(I2)) (IDUMMY,J=1,3),FSQ,(DUMMY,J=1,5), & IDUMMY,M,IE IF (FSQ.LT.0) FSQ=0 M2=M Y2=FSQ/IE END IF SUMM=SUMM-M1+M2 SUM1=SUM1-M1*Y1+M2*Y2 SUM2=SUM2-M1*Y1**2+M2*Y2**2 N=I2-I1+1 AVG=SUM1/SUMM AVGSQ=SUM2/SUMM RMSD=SQRT((FLOAT(N)/(N-1))*(AVGSQ-AVG**2)) AVGFSQ=AVG SIGAVG=RMSD/SQRT(FLOAT(N)) READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,DUMMY,DUMMY, & IC,M,IE,S WRITE (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF, 0.0, 0.0, & IC,M,IE,S,AVGFSQ,SIGAVG END DO 5 CONTINUE C C EVALUATE LOCALLY NORMALIZED STRUCTURE FACTOR MAGNITUDES, AND FIND RE- C SCALING FACTOR TO FORCE = 1. C SUMM=0 SUM2=0 DO I=1,NDATA READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,DUMMY,DUMMY, & IC,M,IE,S,AVGFSQ IF (F.GT.0) THEN E=F/SQRT(IE*AVGFSQ) SUMM=SUMM+M SUM2=SUM2+M*E**2 END IF END DO AVESQ=SUM2/SUMM RESCL=1/SQRT(AVESQ) WRITE (ILP,6000) ATIME,ADATE,TITLE WRITE (ILP,6100) WRITE (ILP,6200) AVESQ,RESCL WRITE (ILP,6001) 6000 FORMAT (/1H1,130('-')/1X, &'SUBPROGRAM NORMAL/BAYES. ',A,', ',A,'. ',A) 6100 FORMAT (/1X, &'NORMALIZED AND RE-SCALED STRUCTURE FACTORS E(HKL) = F(HKL)/SQRT', &'[EPSILON(HKL)*)') 6200 FORMAT (/1X, &'MULTIPLICITY-WEIGHTED AVERAGE ESQ BEFORE RE-SCALING: = ', &' ',F6.3/1X, &'E-MAGNITUDE RE-SCALING FACTOR TO FORCE = 1: RESCl = ', &'1/SQRT() = ',F6.3) 6001 FORMAT (/1X, &' H K L SIN(TH)/L D(HKL) FSQ SIGFSQ ', &' F SIGF E SIGE'/1X, &' - - - --------- ------ --- ------ ', &' - ---- - ----') 6002 FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3) C C STORE RE-SCALED E-MAGNITUDES IN THE WORKING DATA FILE. C DO I=1,NDATA READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,DUMMY,DUMMY, & IC,M,IE,S,AVGFSQ,SIGAVG IF (F.GT.0) THEN E=F/SQRT(IE*AVGFSQ) SIGE=E*SQRT((SIGF**2/(IE*F**2))+(0.5*SIGAVG/AVGFSQ)**2) E=RESCL*E SIGE=RESCL*SIGE WRITE (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ,SIGAVG IF (IPRINT(J,K,L).EQ.1) THEN D=1/(2*S) WRITE (ILP,6002) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE END IF END IF END DO C C COMPILE AND PRINT E-STATISTICS. C CALL ESTATS (NDATA,IO2,ILP,ATIME,ADATE,TITLE,SMIN,SMAX) C C SORT ON d(hkl). C DO I=1,NDATA READ (IO2,REC=I) J,K,L DATA(I)=SINTHL(J,K,L,GINV) END DO CALL SORT (NDATA,DATA,INDEX) C C LIST LOWEST RESOLUTION E-MAGNITUDES. C WRITE (ILP,6000) ATIME,ADATE,TITLE WRITE (ILP,'(/1X,''LOWEST RESOLUTION E-MAGNITUDES:'')') WRITE (ILP,6005) DO I=1,MIN(100,NDATA) READ (IO2,REC=INDEX(I)) J,K,L, FSQ,SIGFSQ,F,SIGF,E,SIGE S=SINTHL(J,K,L,GINV) D=1/(2*S) WRITE (ILP,6006) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,I END DO C C WRITE E-SORTED OUTPUT FILE. C DO I=1,NDATA READ (IO2,REC=I) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE DATA(I)=E END DO CALL SORT (NDATA,DATA,INDEX) OPEN (UNIT=IO1,FILE='edata.bayes',STATUS='NEW',FORM='FORMATTED') DO I=NDATA,1,-1 READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE IF (E.GT.0) CALL WRITE1 (IO1,J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE) END DO CLOSE (UNIT=IO1,STATUS='KEEP') C C LIST LARGEST AND SMALLEST SIGNIFICANT E-MAGNITUDES. C WRITE (ILP,6000) ATIME,ADATE,TITLE WRITE (ILP,6051) WRITE (ILP,6005) I=NDATA N=0 DO WHILE (I.GE.1.AND.N.LT.100) READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (E.GT.3*SIGE) THEN D=1/(2*S) WRITE (ILP,6006) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,NDATA-I+1 N=N+1 END IF I=I-1 END DO WRITE (ILP,6000) ATIME,ADATE,TITLE WRITE (ILP,6052) WRITE (ILP,6005) I=1 N=0 DO WHILE (I.LE.NDATA.AND.N.LT.100) READ (IO2,REC=INDEX(I)) J,K,L,FSQ,SIGFSQ,F,SIGF,E,SIGE, & IC,M,IE,S,AVGFSQ IF (E.GT.3*SIGE) THEN D=1/(2*S) WRITE (ILP,6006) J,K,L,S,D,FSQ,SIGFSQ,F,SIGF,E,SIGE,-I N=N+1 END IF I=I+1 END DO WRITE (ILP,6009) 6051 FORMAT (/1X,'LARGEST E-MAGNITUDES WITH E .GT. 3*SIGMA(E)') 6052 FORMAT (/1X,'SMALLEST E-MAGNITUDES WITH E .GT. 3*SIGMA(E)') 6005 FORMAT (/1X, &' H K L SIN(TH)/L D(HKL) FSQ SIGFSQ ', &' F SIGF E SIGE RANK'/1X, &' - - - --------- ------ --- ------ ', &' - ---- - ---- ----') 6006 FORMAT (1X,3I4,F10.4,F10.2,2F12.2,2F10.2,2F10.3,I10) 6009 FORMAT (/1X, &'OUTPUT FILE OF E-VALUES ESTIMATED AS'/1X,' E(HKL) = F(HKL)/SQR', &'T[EPSILON(HKL)*]'/1X, &'IS AN E-SORTED, FORMATTED, ASCII FILE (IH,IK,IL,FSQ,SIGFSQ,F,SI', &'GF,E,SIGE)'/1X, &'NAMED: "edata.bayes"') RETURN END C----------------------------------------------------------------------- SUBROUTINE NZPLOT (NDATA,IFILE,ILP,ATIME,ADATE,TITLE,PTGP,LAUE) C C COMPILE STATISTICS ON AND PLOT THE CUMULATIVE WISLON DISTRIBUTIONS. C PARAMETER (NX=100,NY=50) CHARACTER AA*1 DIMENSION AA(0:NX,0:NY) CHARACTER ATIME*8,ADATE*9,TITLE*80,PTGP*5,LAUE*5 PARAMETER (NZ=10) DIMENSION ZI(NZ),FA(NZ),FC(NZ) WRITE (ILP,'(/1H1,130(''-'')/1X, &''PROGRAM BAYES. '',A,'', '',A,''. '',A)') ATIME,ADATE,TITLE WRITE (ILP,'(/1X, &''TESTS OF THE VALIDITY OF THE WILSON INTENSITY DISTRIBUTIONS''// &1X, &''POINT GROUP '',A/1X,''LAUE GROUP '',A)') PTGP,LAUE SUMMA=0 SUMFA=0 SUMZA=0 SUMASQ=0 SUMMC=0 SUMFC=0 SUMZC=0 SUMCSQ=0 DO I=1,NZ ZI(I)=FLOAT(I)/NZ FA(I)=0 FC(I)=0 END DO DO I=1,NDATA READ (IFILE,REC=I) J,K,L,FSQ,SIGFSQ,DUMMY,DUMMY,DUMMY,DUMMY, & IC,M,IE,S,AVGFSQ IF (FSQ.LT.-4*SIGFSQ) FSQ=-4*SIGFSQ Z=FSQ/(IE*AVGFSQ) F=SQRT(MAX(0.0,Z)) IF (IC.EQ.0) THEN C C ACENTRIC REFLECTION C SUMMA=SUMMA+M SUMFA=SUMFA+M*F SUMZA=SUMZA+M*Z SUMASQ=SUMASQ+M*Z**2 DO J=1,NZ IF (Z.LE.ZI(J)) FA(J)=FA(J)+M END DO ELSE C C CENTRIC REFLECTION C SUMMC=SUMMC+M SUMFC=SUMFC+M*F SUMZC=SUMZC+M*Z SUMCSQ=SUMCSQ+M*Z**2 DO J=1,NZ IF (Z.LE.ZI(J)) FC(J)=FC(J)+M END DO END IF END DO IF (SUMMA.EQ.0) THEN WRITE (ILP,'(/1X, &''NO ACENTRIC REFLECTIONS IN THE GIVEN POINT GROUP'')') ELSE AVGZA=SUMZA/SUMMA VARZA=SUMASQ/SUMMA-(SUMZA/SUMMA)**2 RHOA=(SUMFA/SUMMA)**2/AVGZA DO I=1,NZ FA(I)=FA(I)/SUMMA END DO WRITE (ILP,'(//1X, &''ACENTRIC WILSON DISTRIBUTION STATISTICS ON Z = FSQ/(EPSILON*<'', &''FSQ/EPSILON>''/1X, &''--------''/1X, &'' EXPERIMENTAL THEORETICAL''/1X, &'' ------------ -----------''/1X, &''MULT.-WTD. NO. OF REFLNS. '',I9/1X, &'' '',F9.3,5X,'' 1.000''/1X, &''<(Z - )**2> '',F9.3,5X,'' 1.000''/1X, &''RHO = **2/ '',F9.3,5X,'' 0.785 = PI/4'')') & NINT(SUMMA),AVGZA,VARZA,RHOA END IF IF (SUMMC.EQ.0) THEN WRITE (ILP,'(/1X, &''NO CENTRIC REFLECTIONS IN THE GIVEN POINT GROUP'')') ELSE AVGZC=SUMZC/SUMMC VARZC=SUMCSQ/SUMMC-(SUMZC/SUMMC)**2 RHOC=(SUMFC/SUMMC)**2/AVGZC DO I=1,NZ FC(I)=FC(I)/SUMMC END DO WRITE (ILP,'(//1X, &''CENTRIC WILSON DISTRIBUTION STATISTICS ON Z = FSQ/(EPSILON*''/1X, &''-------''/1X, &'' EXPERIMENTAL THEORETICAL''/1X, &'' ------------ -----------''/1X, &''MULT.-WTD. NO. OF REFLNS. '',I9/1X, &'' '',F9.3,5X,'' 1.000''/1X, &''<(Z - )**2> '',F9.3,5X,'' 2.000''/1X, &''RHO = **2/ '',F9.3,5X,'' 0.637 = 2/PI'')') & NINT(SUMMC),AVGZC,VARZC,RHOC END IF C C ERASE THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C PLOT 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 PLOT THE THEORETICAL ACENTRIC AND CENTRIC CUMULATIVE DISTRIBUTION C CURVES, N(Z) VS. Z FOR 0 .LE. Z .LE. 1. C XMIN=0 XMAX=1 YMIN=0 YMAX=1 XDEL=XMAX-XMIN YDEL=YMAX-YMIN DX=(XMAX-XMIN)/(2*NX) DO I=0,(2*NX) XI=-XMIN+I*DX IX=NINT(NX*(XI-XMIN)/XDEL) IX=MAX(IX,0) IX=MIN(IX,NX) C C ACENTRIC DISTRIBUTION N(Z) = 1 - EXP(-Z) C YI=1-EXP(-XI) IY=NINT(NY*(YMAX-YI)/YDEL) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='0' C C CENTRIC DISTRIBUTION N(Z) = ERF[SQRT(Z/2)] C YI=ERRF(SQRT(0.5*XI)) IY=NINT(NY*(YMAX-YI)/YDEL) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='1' END DO C C PLOT THE EXPERIMENTAL DATA POINTS. C DO I=1,NZ XI=ZI(I) IX=NINT(NX*(XI-XMIN)/XDEL) IX=MAX(IX,0) IX=MIN(IX,NX) C C ACENTRIC POINTS C IF (SUMMA.GT.0) THEN YI=FA(I) IY=NINT(NY*(YMAX-YI)/YDEL) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='A' END IF C C CENTRIC POINTS C IF (SUMMC.GT.0) THEN YI=FC(I) IY=NINT(NY*(YMAX-YI)/YDEL) IY=MAX(IY,0) IY=MIN(IY,NY) AA(IX,IY)='C' END IF END DO C C PRINT THE PLOT ARRAY. C WRITE (ILP,'(/1H1,130(''-'')/1X,10X, &''PROGRAM BAYES. '',A,'', '',A,''. '',A/)') ATIME,ADATE,TITLE DY=(YMAX-YMIN)/NY DO IY=0,NY IF (MOD(IY,10).EQ.0) THEN WRITE (ILP,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,'(1X,10X,101A1)') (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,'(2X,11F10.1)') (XMIN+IX*DX,IX=0,10) WRITE (ILP,'(/1X,10X, &''CUMULATIVE WILSON DISTRIBUTIONS N(Z) VS. Z = FSQ/(EPSILON*) FOR 0 .LE. Z .LE. 1''//1X,10X, &''ACENTRIC REFLECTIONS: THEORETICAL N(Z) = 1 - EXP(-Z) '', &'' "0", EXPERIMENTAL "A"''/1X,10X, &''CENTRIC " : " N(Z) = ERF[SQRT(Z/2)] '', &'' "1", " "C"'')') 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 PLOT (ILP,IFILE,NDATA,XMIN,XMAX,YMIN,YMAX,IFLAG) PARAMETER (NX=100,NY=50) CHARACTER AA*1 DIMENSION AA(0:NX,0:NY),D(0:10) C C ERASE THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C PLOT 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 PLOT DATA POINTS. C XDEL=XMAX-XMIN YDEL=YMAX-YMIN DO I=1,NDATA READ (IFILE,REC=I) (IDUMMY,J=1,3),(DUMMY,J=1,6),(IDUMMY,J=1,3), & S,YI,ZI XI=S IF (IFLAG.NE.0) YI=ZI 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 PRINT THE PLOT ARRAY. C DY=(YMAX-YMIN)/NY DO IY=0,NY IF (MOD(IY,10).EQ.0) THEN WRITE (ILP,'(1X,E10.3,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,'(1X,10X, 101A1)') (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,'(2X,11F10.3)') (2*(XMIN+IX*DX),IX=0,10) C C PRINT D-SPACINGS. C DO IX=0,10 XI=XMIN+IX*DX IF (XI.EQ.0) THEN D(IX)=9999.99 ELSE D(IX)=1/(2*ABS(XI)) END IF END DO WRITE (ILP,'(2X,11(F9.2,1X),'' d-SPACING (A)'')') (D(IX),IX=0,10) 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 ERASE THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C PLOT 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 PLOT 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 PLOT 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,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,'(1X,10X, 101A1)') (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,'(2X,11F10.1)') (XMIN+IX*DX,IX=0,10) WRITE (ILP,'(/1X,10X, &''ACENTRIC DISTRIBUTION: P(E) = 2*E*EXP(-E**2), "*" THEORETIC'', &''AL, "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 ERASE THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C PLOT 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 PLOT 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 PLOT 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,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,'(1X,10X, 101A1)') (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,'(2X,11F10.1)') (XMIN+IX*DX,IX=0,10) WRITE (ILP,'(/1X,10X, &''CENTRIC DISTRIBUTION: P(E) = SQRT(2/PI)*EXP(-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),D(0:10) C C ERASE THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C PLOT 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(TH)/L < 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 PLOT = 1. C DO IX=0,NX AA(IX,30)='.' END DO C C PLOT 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 PLOT 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,'(1X,F10.1,101A1)') YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,'(1X,10X, 101A1)') (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,'(2X,11F10.3)') (2*(XMIN+IX*DX),IX=0,10) C C PRINT D-SPACINGS. C DO IX=0,10 XI=XMIN+IX*DX IF (XI.EQ.0) THEN D(IX)=9999.99 ELSE D(IX)=1/(2*ABS(XI)) END IF END DO WRITE (ILP,'(2X,11(F9.2,1X),'' d-SPACING (A)'')') (D(IX),IX=0,10) WRITE (ILP,'(/1X,10X, &'' VERSUS = 2* DISTRIBUTION'')') RETURN END C----------------------------------------------------------------------- SUBROUTINE READ1 (ITYPE,IFILE,IEND,IH,IK,IL,FSQ,SIGFSQ) 1 IF (ITYPE.EQ.0) READ (IFILE,*,END=9) IH,IK,IL,FSQ,SIGFSQ IF (ITYPE.EQ.1) READ (IFILE,*,END=9) IH,IK,IL,F, SIGF IF (ITYPE.EQ.2) READ (IFILE, END=9) IH,IK,IL,FSQ,SIGFSQ IF (ITYPE.EQ.3) READ (IFILE, END=9) IH,IK,IL,F, SIGF IF (ITYPE.EQ.1.OR.ITYPE.EQ.3) THEN IF (F.LT.0) F=0 IF (SIGF.LT.0) SIGF=0 FSQ=F**2 SIGFSQ=MAX(2*F*SIGF,4*SIGF**2) END IF IF (SIGFSQ.LE.0) GO TO 1 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 Q=0 DO I=1,3 DO J=1,3 Q=Q+H(J)*GINV(I,J)*H(I) END DO END DO SINTHL=0.5*SQRT(Q) 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 ADAPTED 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. 229-233. CAMBRIDGE, ENGLAND: C CAMBRIDGE UNIVERSITY PRESS. C DIMENSION DATA(N),INDEX(N) DO I=1,N INDEX(I)=I END DO 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 DO WHILE (L.LE.J) 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 END DO 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 ADAPTED 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=100) 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 ADAPTED 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 DO WHILE (I2-I1.GT.1) I=(I2+I1)/2 IF (X(I).GT.XI) THEN I2=I ELSE I1=I END IF END DO 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----------------------------------------------------------------------- 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===================================================================== 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-----------------------------------------------------------------------