C..... C .. PROGRAM ABSORB2 C..... C .. APPLY THE ABSORPTION CORRECTION TO THE DATA C .. USING THE CONDITIONS SET UP IN THE ABSORB.DAT C .. FILE. PRINT OUT SUMMARIES OF THE CORRECTION. C..... PROGRAM ABSORB2 PARAMETER (M1=9,M2=15,M3=50) CHARACTER TITLE*80, ATIME*8, ADATE*9 CHARACTER FIL01*80, FIL02*80 DIMENSION CELL(6) CHARACTER*4 WAVE CHARACTER*2 ATOMS(10) DIMENSION STOIC(10) DIMENSION FACES(50,4) DIMENSION NGP(3), D1(3),D2(3),D3(3) DIMENSION ORMAT(3,3) DIMENSION NTRP(50) COMMON /INPT01/ TITLE, ATIME, ADATE COMMON /INPT02/ FIL01, FIL02, IPATH COMMON /INPT03/ CELL COMMON /INPT04/ WAVE COMMON /LAMBDA/ FLAM COMMON /INPT05/ ATOMS COMMON /INPT06/ ZCELL, STOIC COMMON /INPT07/ FACES COMMON /INPT08/ NGP, D1,D2,D3 COMMON /INPT09/ ORMAT, NDIFF COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP COMMON /INPT11/ NBPR, A0, A1, A2, A3 COMMON /INPT12/ NSOL, EMUS, NTRP COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C- DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3) DIMENSION FS(50,5) DIMENSION VR(150,3), MM(150,3) PARAMETER (NGMAX=32**3) DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3) DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM COMMON /CRYSTL/ TMOLW, RHO, EMUX COMMON /MATRIX/ UBF, UBI, UBGF, UBGI COMMON /PLANES/ NF, FS COMMON /VERTEX/ NV, VR, MM COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB C- CHARACTER HEADING*120 DIMENSION HEADING(M1),JX(M2),JXMIN(M1,M2,M3),JXMAX(M1,M2,M3) DIMENSION ANGLS(4) DIMENSION SIP(3), SDP(3), SIT(3), SDT(3) DIMENSION SIR(3), SDR(3), SIS(3), SDS(3) DIMENSION PMT(3), PMP(3), PMS(3) DATA NFL1/30/, NFL2/40/ C- DOUBLE PRECISION ABSCO1,ABSCO2,ABSCO3,ABSCO4,ABSCOT,TBAR DIMENSION S0(3),S1(3),COSIR(3),COSDR(3),COSIS(3),COSDS(3) C- NIN=IO5 NOU=ILP CALL PAGE1 CALL CELL2 CALL MATX2 CALL EMUX2 CALL FACE2 CALL GRID2 CALL KAPP2 C- OPEN(UNIT=NFL1,FILE=FIL01,STATUS='OLD',FORM='UNFORMATTED',ERR=997) OPEN(UNIT=NFL2,FILE=FIL02,STATUS='NEW',FORM='UNFORMATTED') C- DO 10 J1=1,M1 DO 10 J2=1,M2 DO 10 J3=1,M3 JXMIN(J1,J2,J3)=+1E7 JXMAX(J1,J2,J3)=-1E7 10 CONTINUE NTOTL = 0 1000 READ(NFL1,END=1500,ERR=997) * NREF, NH,NK,NL, ANGLS, FSQ,SIG, XTM IF (NDIFF.EQ.1) CALL P3ROT (ANGLS, SIP,SDP,PMP) IF (NDIFF.EQ.2) CALL C4ROT (ANGLS, SIP,SDP,PMP) IF (NDIFF.EQ.3) CALL BLROT (ANGLS, SIP,SDP,PMP) C- GET REAL SPACE COMPONENTS OF SI AND SD: SIR AND SDR CALL MTV1 (UBGI,SIP, SIR) CALL MTV1 (UBGI,SDP, SDR) C- GET RECIPROCAL SPACE COMPONENTS OF SI AND SD: SIS AND SDS CALL MTV1 (UBI,SIP, SIS) CALL MTV1 (UBI,SDP, SDS) C- GET RECIPROCAL SPACE COMPONENTS OF PM: PMS CALL MTV1 (UBI,PMP, PMS) ABSCO1 = 0. ABSCO2 = 0. ABSCO3 = 0. ABSCO4 = 0. ABSCOT = 0. TBAR = 0. DO 100 N=1,NPTS W = WABC(N) C- CRYSTAL CONTRIBUTION CALL PATH1 (SIR,N, PXI,NPI) CALL PATH1 (SDR,N, PXD,NPD) PXTL = PXI + PXD AXTL = EXP( -EMUX * PXTL ) TXTL = PXTL * AXTL C- CAPILLARY CONTRIBUTION AKAP = 1. IF(NKAP.EQ.0) GO TO 50 CALL MTV1 (UPTSI,SIP, SIT) CALL MTV1 (UPTSI,SDP, SDT) CALL PATH2 (SIT,N, W1I,W2I,PKI) CALL PATH2 (SDT,N, W1D,W2D,PKD) PKAP = PKI + PKD AKAP = EXP( -EMUK * PKAP ) 50 CONTINUE C- BEAM NONHOMOGENIETY CONTRIBUTION APRO = 1. IF(NBPR.EQ.0) GO TO 60 CALL BPRO3 (PMS,N, DEL,APRO) 60 CONTINUE C- TRAPPED SOLVENT CONTRIBUTION ASOL = 1. IF(NSOL.EQ.0) GO TO 70 CALL PATH3 (W1I, PSI) CALL PATH3 (W1D, PSD) PSOL = PSI + PSD ASOL = EXP( -EMUS * PSOL ) 70 CONTINUE C- TOTAL UP THE VARIOUS CONTRIBUTIONS ABSCO1 = ABSCO1 + W * AXTL ABSCO2 = ABSCO2 + W * APRO ABSCO3 = ABSCO3 + W * AKAP ABSCO4 = ABSCO4 + W * ASOL ABSCOT = ABSCOT + W * AXTL * APRO * AKAP * ASOL TBAR = TBAR + W * TXTL 100 CONTINUE C- AFSQ = FSQ / ABSCOT ASIG = SIG / ABSCOT TBAR = TBAR / ABSCO1 C C WRITE OUT THE ABSORPTION CORRECTED DATA. C C FOR OUTPUT OF THE REVERSE-INCIDENT AND DIFFRACTED BEAM DIRECTIONS IN C THE CRYSTAL, CONVERT FROM UNIT DIRECTION VECTOR COMPONENTS TO C DIRECTION COSINES REFERRED TO THE CRYSTALLOGRAPHIC DIRECT SPACE AXES. C CALL DIRCOS (SIR,GR,COSIR) CALL DIRCOS (SDR,GR,COSDR) CALL DIRCOS (SIS,GS,COSIS) CALL DIRCOS (SDS,GS,COSDS) DO I=1,3 S0(I)=COSIR(I) S1(I)=COSDR(I) END DO C IF(IPATH.EQ.0) WRITE(NFL2) NREF, NH,NK,NL, ANGLS, AFSQ,ASIG, XTM IF(IPATH.EQ.1) WRITE(NFL2) NREF, NH,NK,NL, ANGLS, AFSQ,ASIG, XTM, & SNGL(TBAR),(S0(I),I=1,3),(S1(I),I=1,3) NTOTL = NTOTL + 1 C C CATALOG PROFILES WITH EXTREME FEATURES. C IF(NREF.LT.0) GO TO 1000 IF(FSQ.LE.2.0*SIG) GO TO 1000 TWOTH=ANGLS(1) IF(NDIFF.EQ.2) TWOTH=2.0*TWOTH JX(1)=NREF JX(2)=NH JX(3)=NK JX(4)=NL JX(5)= 1000*FSQ JX(6)= 1000*SIG JX(7)= 1000*AFSQ JX(8)= 1000*ASIG JX(9)= 1000*ABSCO1 JX(10)=1000*ABSCO2 JX(11)=1000*ABSCO3 JX(12)=1000*ABSCO4 JX(13)=1000*ABSCOT JX(14)=1000*TBAR JX(15)=1000*TWOTH CALL EXTREM (JXMIN,JXMAX,JX) GO TO 1000 C- 1500 CLOSE(UNIT=NFL1,STATUS='KEEP') CLOSE(UNIT=NFL2,STATUS='KEEP') C- C- START OUTPUT WITH SUMMARY OF INPUT DATA C- CALL OUTP1 C C PRINT CATALOGS OF EXTREME CASES. C OPEN(UNIT=NLP,FILE='absorb.lp',STATUS='UNKNOWN') DO I=1,1E9 READ (NLP,'(A1)',END=1) DUMMY END DO 1 WRITE (NLP,'(/1X)') DO 49 J1=1,M1 IF(J1.EQ.4 .AND. NBPR.EQ.0) GO TO 49 IF(J1.EQ.5 .AND. NKAP.EQ.0) GO TO 49 IF(J1.EQ.6 .AND. NSOL.EQ.0) GO TO 49 WRITE (NLP,610) ATIME,ADATE,TITLE WRITE (NLP,645) HEADING(J1) DO 59 J3=1,MIN(M3,NTOTL) I1 =JXMIN(J1,1,J3) I2 =JXMIN(J1,2,J3) I3 =JXMIN(J1,3,J3) I4 =JXMIN(J1,4,J3) F5 =JXMIN(J1,5,J3)*0.001 F6 =JXMIN(J1,6,J3)*0.001 F7 =JXMIN(J1,7,J3)*0.001 F8 =JXMIN(J1,8,J3)*0.001 F9 =JXMIN(J1,9,J3)*0.001 F10=JXMIN(J1,10,J3)*0.001 F11=JXMIN(J1,11,J3)*0.001 F12=JXMIN(J1,12,J3)*0.001 F13=JXMIN(J1,13,J3)*0.001 F14=JXMIN(J1,14,J3)*0.001 F15=JXMIN(J1,15,J3)*0.001 WRITE (NLP,646) I1,I2,I3,I4, * F5,F6,F7,F8,F9,F10,F11,F12,F13,F14,F15 59 CONTINUE WRITE (NLP,610) ATIME,ADATE,TITLE WRITE (NLP,645) HEADING(J1) DO 69 J3=1,MIN(M3,NTOTL) I1 =JXMAX(J1,1,J3) I2 =JXMAX(J1,2,J3) I3 =JXMAX(J1,3,J3) I4 =JXMAX(J1,4,J3) F5 =JXMAX(J1,5,J3)*0.001 F6 =JXMAX(J1,6,J3)*0.001 F7 =JXMAX(J1,7,J3)*0.001 F8 =JXMAX(J1,8,J3)*0.001 F9 =JXMAX(J1,9,J3)*0.001 F10=JXMAX(J1,10,J3)*0.001 F11=JXMAX(J1,11,J3)*0.001 F12=JXMAX(J1,12,J3)*0.001 F13=JXMAX(J1,13,J3)*0.001 F14=JXMAX(J1,14,J3)*0.001 F15=JXMAX(J1,15,J3)*0.001 WRITE (NLP,646) I1,I2,I3,I4, * F5,F6,F7,F8,F9,F10,F11,F12,F13,F14,F15 69 CONTINUE 49 CONTINUE C C FORMAT STATEMENTS FOR PRINTED OUTPUT C DATA HEADING / 1' FSQ: F SQUARE MAGNITUDE UNCORRECTED', 2' FSQ/A: F SQUARE MAGNITUDE CORRECTED', 3' AXTL: TRANSMISSION FACTOR FOR CRYSTAL ONLY', 4' BPRO: BEAM NON-HOMOGENIETY PROFILE CORRECTION', 5' AKAP: TRANSMISSION FACTOR FOR CAPILLARY ONLY', 6' ASOL: TRANSMISSION FACTOR FOR TRAPPED SOLVENT ONLY', 7' ATOT: TRANSMISSION FACTOR FOR ALL ABSORPTION EFFECTS', 8' TBAR: TRANSMISSION PATH LENGTH THROUGH CRYSTAL', 9' 2THETA: TWO THETA MAGNITUDE' / 610 FORMAT (1H1/1X,'PROGRAM ABSORB. 'A,', ',A'. 'A) 645 FORMAT (/1X,A//1X, &' NREF',' H K L', &' FSQ',' SIG ', &' FSQ/A',' SIG/A ', &' AXTL',' BPRO',' AKAP',' ASOL',' ATOT', &' TBAR',' 2THETA',/) 646 FORMAT (1X,I7,3I4,2(F9.2,F6.2,2X),5F7.3,F6.3,F8.3) CLOSE(UNIT=NLP,STATUS='KEEP') STOP 'ABSORB2 IS DONE' C- 997 WRITE(IO6,9003) 9003 FORMAT(/1X,'ERROR READING OR OPENING INPUT FILE') STOP END C----------------------------------------------------------------------- SUBROUTINE DIRCOS (U,G,COSU) C C DIRECTION COSINES OF VECTOR U REFERRED TO AXES WITH METRIC MATRIX G. C C U(I)*V(J)*G(I,J) C COS PHI = ------------------------------------------- C SQRT{[U(I)*U(J)*G(I,J)]*[V(I)*V(J)*G(I,J)]} C DIMENSION U(3),G(3,3),COSU(3) DIMENSION X(3),Y(3),Z(3) DATA X /1, 0, 0/ DATA Y /0, 1, 0/ DATA Z /0, 0, 1/ DO I=1,3 COSU(I)=0 END DO DO I=1,3 DO J=1,3 COSU(1)=COSU(1)+U(I)*X(J)*G(I,J) COSU(2)=COSU(2)+U(I)*Y(J)*G(I,J) COSU(3)=COSU(3)+U(I)*Z(J)*G(I,J) END DO END DO USQ=0 DO I=1,3 DO J=1,3 USQ=USQ+U(I)*U(J)*G(I,J) END DO END DO DO I=1,3 COSU(I)=COSU(I)/SQRT(USQ*G(I,I)) END DO RETURN END C----------------------------------------------------------------------- C..... C .. SUBROUTINE EXTREME C..... SUBROUTINE EXTREM (IXMIN,IXMAX,IX) PARAMETER (N1=9,N2=15,N3=50) DIMENSION IXMIN(N1,N2,N3),IXMAX(N1,N2,N3),ITYPE(N1),IX(N2) DATA ITYPE /5, 7, 9, 10, 11, 12, 13, 14, 15/ C C TABULATE IN RANKED ORDER THE N3 SMALLEST AND N3 LARGEST C VALUES OF THE VARIABLE OF TYPE I1. C C I1 TYPE OF VARIABLE C I2 VALUE OF VARIABLE C I3 RANK IN LIST C DO 10 I1=1,N1 I2=ITYPE(I1) DO 11 I3=1,N3 IF (IX(I2).GE.IXMIN(I1,I2,I3)) GO TO 11 DO 12 J3=N3,(I3+1),-1 DO 12 I2=1,N2 12 IXMIN(I1,I2,J3)=IXMIN(I1,I2,J3-1) DO 13 I2=1,N2 13 IXMIN(I1,I2,I3)=IX(I2) GO TO 20 11 CONTINUE 20 CONTINUE I2=ITYPE(I1) DO 21 I3=1,N3 IF (IX(I2).LE.IXMAX(I1,I2,I3)) GO TO 21 DO 22 J3=N3,(I3+1),-1 DO 22 I2=1,N2 22 IXMAX(I1,I2,J3)=IXMAX(I1,I2,J3-1) DO 23 I2=1,N2 23 IXMAX(I1,I2,I3)=IX(I2) GO TO 10 21 CONTINUE 10 CONTINUE RETURN END C..... C .. BLOCKDATA ABSORB3 C..... BLOCKDATA ABSORB3 C- CHARACTER TITLE*80, ATIME*8, ADATE*9 CHARACTER FIL01*80, FIL02*80 DIMENSION CELL(6) CHARACTER*4 WAVE CHARACTER*2 ATOMS(10) DIMENSION STOIC(10) DIMENSION FACES(50,4) DIMENSION NGP(3), D1(3),D2(3),D3(3) DIMENSION ORMAT(3,3) DIMENSION NTRP(50) COMMON /INPT01/ TITLE, ATIME, ADATE COMMON /INPT02/ FIL01, FIL02, IPATH COMMON /INPT03/ CELL COMMON /INPT04/ WAVE COMMON /LAMBDA/ FLAM COMMON /INPT05/ ATOMS COMMON /INPT06/ ZCELL, STOIC COMMON /INPT07/ FACES COMMON /INPT08/ NGP, D1,D2,D3 COMMON /INPT09/ ORMAT, NDIFF COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP COMMON /INPT11/ NBPR, A0, A1, A2, A3 COMMON /INPT12/ NSOL, EMUS, NTRP COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C- DATA TITLE/'ABSORB1...VERSION OF 10/II/82'/ DATA ATIME/'hh:mm:ss'/, ADATE/'dd-mmm-yy'/ DATA FIL01/'FIL01.DATA'/, FIL02/'FIL02.DATA'/, IPATH/0/ DATA CELL/1.,1.,1.,90.,90.,90./ DATA WAVE/'MOKA'/ DATA FLAM/0.7107/ DATA ATOMS/10*' '/ DATA ZCELL/1./, STOIC/10*0./ DATA FACES/200*0./ DATA NGP/3*2/ DATA D1/1.,0.,0./, D2/0.,1.,0./, D3/0.,0.,1./ DATA ORMAT/9*0./ DATA NDIFF/1/ DATA NKAP/0/, DKAP/1./, TKAP/0.01/, TAU/0./, SGM/0./ DATA XDSP/0./, YDSP/0./ DATA NBPR/0/, A0/1./, A1/0./, A2/0./, A3/0./ DATA NSOL/0/, EMUS/0./, NTRP/50*0/ DATA IO5/5/, IO6/6/, NLP/20/, NAB/30/, NIN/1/, NOU/1/ END C..... C .. SUBROUTINE AMI3 C..... C .. INVERT AN ASYMMETRIC 3 X 3 MATRIX C..... SUBROUTINE AMI3(A,B,DET) DIMENSION A(3,3), B(3,3) B(1,1) = A(2,2)*A(3,3) - A(3,2)*A(2,3) B(2,1) = A(2,1)*A(3,3) - A(3,1)*A(2,3) B(3,1) = A(2,1)*A(3,2) - A(3,1)*A(2,2) DET = A(1,1)*B(1,1) - A(1,2)*B(2,1) + A(1,3)*B(3,1) IF (DET .EQ. 0.0) RETURN B(1,2) = A(1,2)*A(3,3) - A(3,2)*A(1,3) B(2,2) = A(1,1)*A(3,3) - A(3,1)*A(1,3) B(3,2) = A(1,1)*A(3,2) - A(3,1)*A(1,2) B(1,3) = A(1,2)*A(2,3) - A(2,2)*A(1,3) B(2,3) = A(1,1)*A(2,3) - A(2,1)*A(1,3) B(3,3) = A(1,1)*A(2,2) - A(2,1)*A(1,2) DO 10 I=1,3 DO 15 J=1,3 B(I,J) = (-1.0)**(I+J)*B(I,J)/DET 15 CONTINUE 10 CONTINUE RETURN END C..... C .. SUBROUTINE BLRMT C..... C .. GET THE BUSING-LEVY ROTATION MATRIX FOR THE C .. BUSING-LEVY DIFFRACTOMETER ACTA CRYST 22,457-64(1967) C .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI C..... SUBROUTINE BLRMT(ANGLS, RMTF,ST,CT) DIMENSION ANGLS(4), RMTF(3,3) DATA CONS/57.29577951/ C- TTH = ANGLS(1)/CONS OME = ANGLS(2)/CONS PHI = ANGLS(3)/CONS CHI = ANGLS(4)/CONS THT = 0.5*TTH ST = SIN(THT) CT = COS(THT) SO = SIN(OME) CO = COS(OME) SP = SIN(PHI) CP = COS(PHI) SC = SIN(CHI) CC = COS(CHI) C- EXPLICIT GENERATION OF ROTATION MATRIX RMTF(1,1) = +CO*CC*CP-SO*SP RMTF(1,2) = +CO*CC*SP+SO*CP RMTF(1,3) = +CO*SC RMTF(2,1) = -SO*CC*CP-CO*SP RMTF(2,2) = -SO*CC*SP+CO*CP RMTF(2,3) = -SO*SC RMTF(3,1) = -SC*CP RMTF(3,2) = -SC*SP RMTF(3,3) = +CC C- RETURN END C..... C .. SUBROUTINE BLROT C..... C .. GET THE [PHI] FRAME COMPONENTS OF THE UNIT C .. VECTORS FOR THE BUSING-LEVY DIFFRACTOMETER C .. -SI- REVERSE INCIDENT BEAM C .. -SD- DIFFRACTED BEAM C .. -PM- MONOCHROMATOR VERTICAL PROFILE C .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI C..... SUBROUTINE BLROT(ANGLS, SIP,SDP,PMP) DIMENSION ANGLS(4), SIP(3), SDP(3), PMP(3) DIMENSION RMTF(3,3), RMTI(3,3), SDT(3), SIT(3), PMT(3) C- CALL BLRMT (ANGLS, RMTF,ST,CT) C- C- TRANSPOSE OF ROTATION MATRIX EQUALS ITS INVERSE DO 100 I=1,3 DO 50 J=1,3 RMTI(I,J) = RMTF(J,I) 50 CONTINUE 100 CONTINUE C- SIT(1) = +ST SIT(2) = -CT SIT(3) = 0.0 SDT(1) = +ST SDT(2) = +CT SDT(3) = 0.0 PMT(1) = 0.0 PMT(2) = 0.0 PMT(3) = 1.0 C- CALL MTV1(RMTI,SIT, SIP) CALL MTV1(RMTI,SDT, SDP) CALL MTV1(RMTI,PMT, PMP) RETURN END C..... C .. SUBROUTINE BPRO2 C..... C .. WRITE THE ANALYTICAL FORM OF THE BEAM PROFILE C .. THE EQUATION IS OF THE FORM C .. F = A0 + A1*DEL +A2*DEL*DEL +A3*DEL*DEL*DEL C..... SUBROUTINE BPRO2 COMMON /INPT11/ NBPR, A0, A1, A2, A3 COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU IF(NBPR.EQ.0) GO TO 100 IF(NBPR.EQ.1) GO TO 200 IF(NBPR.EQ.2) GO TO 300 IF(NBPR.EQ.3) GO TO 400 IF(NBPR.EQ.4) GO TO 500 100 WRITE(NOU,9001) 9001 FORMAT(/1X,'NO MONOCROMATOR / NO BEAM PROFILE') RETURN 200 WRITE(NOU,9002) A0 9002 FORMAT(/1X,'BEAM: SCALE CORRECTION / A0: ',F8.4) RETURN 300 WRITE(NOU,9003) A0, A1 9003 FORMAT(/1X,'BEAM: LINEAR CORRECTION / A0,A1: ',2F8.4) RETURN 400 WRITE(NOU,9004) A0, A1, A2 9004 FORMAT(/1X,'BEAM: QUADRATIC CORRECTION / A0,A1,A2: ',3F8.4) RETURN 500 WRITE(NOU,9005) A0, A1, A2, A3 9005 FORMAT(/1X,'BEAM: CUBIC CORRECTION /A0,A1,A2,A3: ',4F8.4) RETURN END C..... C .. SUBROUTINE BPRO3 C..... C .. CALCULATE THE BEAM PROFILE CONTRIBUTION ASSUMING A C .. PERFECTLY PARALLEL BEAM AND A NONHOMOGENIETY IN ONE C .. DIRECTION ONLY C .. DIFFRACTOMETER NONHOMOGENIETY DIRECTION C .. NICOLET P3 HORIZONTAL C .. NONIUS CAD4 VERTICAL C .. BUSING-LEVY VERTICAL C..... C .. BEAM INTENSITY GIVEN IN ANALYTICAL FORM FOR A LEAST-SQUARES FIT C .. OF COUNT RATE VERSUS DISPLACEMENT FROM THE DIFFRACTOMETER CENTER C .. COUNT RATE IS IN FRACTIONAL FORM C .. DISPLACEMENT IS IN MILLIMETERS C .. FIT IS TO AN EQUATION OF THE FORM C .. BPF = A0 + A1*DEL + A2*DEL**2 + A3*DEL**3 C .. BPF IS THE FRACTIONAL COUNT RATE C .. DEL IS THE PROJECTION OF POINT RABC ON NONHOMOGENIETY AXIS C .. A0, A1, A2, A3 ARE FIT TO EXPERIMENTAL POINTS C..... SUBROUTINE BPRO3 (PM,N, DEL,BPF) COMMON /INPT11/ NBPR, A0, A1, A2, A3 PARAMETER (NGMAX=32**3) DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3) COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC DIMENSION PM(3) C- DEL = PM(1)*RABC(N,1) + PM(2)*RABC(N,2) + PM(3)*RABC(N,3) BPF = 1. IF(NBPR.EQ.0) RETURN BPF = A0 IF(NBPR.EQ.1) RETURN BPF = BPF + A1*DEL IF(NBPR.EQ.2) RETURN BPF = BPF + A2*DEL*DEL IF(NBPR.EQ.3) RETURN BPF = BPF + A3*DEL*DEL*DEL RETURN END C..... C .. SUBROUTINE C4RMT C..... C .. GET THE BUSING-LEVY ROTATION MATRIX FOR C .. THE ENRAF-NONIUS CAD-4 DIFFRACTOMETER C .. ANGLES ARE: THETA, PHIK, OMGK, KAPPA C..... SUBROUTINE C4RMT(ANGLS, RMTF,ST,CT) DIMENSION ANGLS(4), RMTF(3,3) DIMENSION A(3,3), B(3,3), C(3,3), D(3,3), E(3,3) DIMENSION DE(3,3), CDE(3,3), BCDE(3,3), ABCDE(3,3) CON = 57.2957795 ALPHA = 49.96681/CON C- THETA = ANGLS(1)/CON PHIK = ANGLS(2)/CON OMGK = ANGLS(3)/CON FKAP = ANGLS(4)/CON C C CORRECTION OF CODE 10 MARCH 1987 (GDT) C OMEGA(KAPPA) IS REALLY OMEGA(BUSING-LEVY) PLUS THETA C OMGK = OMGK - THETA C- SA = SIN(ALPHA) CA = COS(ALPHA) ST = SIN(THETA) CT = COS(THETA) SP = SIN(PHIK) CP = COS(PHIK) SO = SIN(OMGK) CO = COS(OMGK) SK = SIN(FKAP) CK = COS(FKAP) C- A(1,1) = +CO A(1,2) = +SO A(1,3) = +0.0 A(2,1) = -SO A(2,2) = +CO A(2,3) = +0.0 A(3,1) = +0.0 A(3,2) = +0.0 A(3,3) = +1.0 C- B(1,1) = +CA B(1,2) = +0.0 B(1,3) = +SA B(2,1) = +0.0 B(2,2) = +1.0 B(2,3) = +0.0 B(3,1) = -SA B(3,2) = +0.0 B(3,3) = +CA C- C(1,1) = +CK C(1,2) = +SK C(1,3) = +0.0 C(2,1) = -SK C(2,2) = +CK C(2,3) = +0.0 C(3,1) = +0.0 C(3,2) = +0.0 C(3,3) = +1.0 C- D(1,1) = +CA D(1,2) = +0.0 D(1,3) = -SA D(2,1) = +0.0 D(2,2) = +1.0 D(2,3) = +0.0 D(3,1) = +SA D(3,2) = +0.0 D(3,3) = +CA C- E(1,1) = +CP E(1,2) = +SP E(1,3) = +0.0 E(2,1) = -SP E(2,2) = +CP E(2,3) = +0.0 E(3,1) = +0.0 E(3,2) = +0.0 E(3,3) = +1.0 C- CALL MTM1 (D,E, DE) CALL MTM1 (C,DE, CDE) CALL MTM1 (B,CDE, BCDE) CALL MTM1 (A,BCDE, ABCDE) DO 100 I=1,3 DO 50 J=1,3 RMTF(I,J) = ABCDE(I,J) 50 CONTINUE 100 CONTINUE C- RETURN END C..... C .. SUBROUTINE C4ROT C..... C .. GET THE [PHI] FRAME COMPONENTS OF THE UNIT C .. VECTORS FOR THE ENRAF-NONIUS CAD-4 DIFFRACTOMETER C .. -SI- REVERSE INCIDENT BEAM C .. -SD- DIFFRACTED BEAM C .. -PM- MONOCHROMATOR VERTICAL PROFILE C .. ANGLES ARE: THETA, PHIK, OMGK, KAPPA C..... SUBROUTINE C4ROT(ANGLS, SIP,SDP,PMP) DIMENSION ANGLS(4), SIP(3), SDP(3), PMP(3) DIMENSION RMTF(3,3), RMTI(3,3), SIT(3), SDT(3), PMT(3) C- CALL C4RMT (ANGLS, RMTF,ST,CT) C- INVERSE OF ROTATION MATRIX EQUALS ITS TRANSPOSE DO 100 I=1,3 DO 50 J=1,3 RMTI(I,J) = RMTF(J,I) 50 CONTINUE 100 CONTINUE C SIT(1) = +CT SIT(2) = +ST SIT(3) = +0.0 SDT(1) = -CT SDT(2) = +ST SDT(3) = +0.0 C PMT(1) = +0.0 PMT(2) = +0.0 PMT(3) = +1.0 C- CALL MTV1 (RMTI,SIT, SIP) CALL MTV1 (RMTI,SDT, SDP) CALL MTV1 (RMTI,PMT, PMP) C- RETURN END C..... C .. SUBROUTINE CELL2 C..... C .. GIVEN THE UNIT CELL CONSTANTS C .. GENERATE THE REAL AND RECIPROCAL CELL CONSTANTS C .. THE REAL AND RECIPROCAL METRIC TENSORS C .. AND THE VOLUMES. C..... SUBROUTINE CELL2 DIMENSION CELL(6) COMMON /INPT03/ CELL DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DATA CON/57.29577951/ C- WRITE OUT TERMS EXPLICITLY A = CELL(1) B = CELL(2) C = CELL(3) ALP = CELL(4) BET = CELL(5) GAM = CELL(6) C- CONVERT ANGLES TO RADIANS CALP = COS (ALP/CON) SALP = SIN (ALP/CON) CBET = COS (BET/CON) SBET = SIN (BET/CON) CGAM = COS (GAM/CON) SGAM = SIN (GAM/CON) IF( ALP.EQ.90.) CALP = 0. IF( ALP.EQ.90.) SALP = 1. IF( BET.EQ.90.) CBET = 0. IF( BET.EQ.90.) SBET = 1. IF( GAM.EQ.90.) CGAM = 0. IF( GAM.EQ.90.) SGAM = 1. C- BUILD THE GR MATRIX EXPLICITLY GR(1,1) = A*A GR(1,2) = A*B*CGAM GR(1,3) = A*C*CBET GR(2,1) = GR(1,2) GR(2,2) = B*B GR(2,3) = B*C*CALP GR(3,1) = GR(1,3) GR(3,2) = GR(2,3) GR(3,3) = C*C C- INVERT GR TO GET GS AND VOL**2 CALL AMI3 (GR, GS, V2) CVR = SQRT(V2) CVS = 1.0/CVR C- GET RECIPROCAL CELL EXPLICITLY AST = B*C*SALP*CVS BST = C*A*SBET*CVS CST = A*B*SGAM*CVS CALS = (CBET*CGAM - CALP) / (SBET*SGAM) SALS = SQRT(1. - CALS*CALS) CBES = (CALP*CGAM - CBET) / (SALP*SGAM) SBES = SQRT(1. - CBES*CBES) CGAS = (CALP*CBET - CGAM) / (SALP*SBET) SGAS = SQRT(1. - CGAS*CGAS) IF( ABS(CALS).LE.1E-8 ) CALS = 0. IF( ABS(CBES).LE.1E-8 ) CBES = 0. IF( ABS(CGAS).LE.1E-8 ) CGAS = 0. IF( CALS.EQ.0. ) SALS = 1. IF( CBES.EQ.0. ) SBES = 1. IF( CGAS.EQ.0. ) SGAS = 1. C- STUFF RR AND RS RR(1) = A RR(2) = B RR(3) = C RR(4) = CALP RR(5) = CBET RR(6) = CGAM RS(1) = AST RS(2) = BST RS(3) = CST RS(4) = CALS RS(5) = CBES RS(6) = CGAS C- BUILD THE B MATRIX EXPLICITLY BM(1,1) = AST BM(1,2) = BST*CGAS BM(1,3) = CST*CBES BM(2,1) = 0.0 BM(2,2) = BST*SGAS BM(2,3) = -CST*SBES*CALP BM(3,1) = 0.0 BM(3,2) = 0.0 BM(3,3) = CST*SBES*SALP RETURN END C..... C .. SUBROUTINE CELL3 C..... C .. PRINT OR TYPE OUT THE CELL INFORMATION. C..... SUBROUTINE CELL3 DIMENSION CELL(6) COMMON /INPT03/ CELL DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C- WRITE(NOU,102) CELL 102 FORMAT(/1X,'A= ',3X,F7.3,4X,'B=',3X,F7.3,3X,'C=',4X,F7.3,/, * 1X,'ALPHA=',F7.3,3X,' BETA=',F7.3,3X,'GAMMA=',F7.3) WRITE(NOU,103) 103 FORMAT(/1X,'REAL METRIC TENSOR') WRITE(NOU,104) ( (GR(I,J),J=1,3), I=1,3 ) 104 FORMAT(1X,3E15.7) WRITE(NOU,105) 105 FORMAT(1X,'RECIPROCAL METRIC TENSOR') WRITE(NOU,104) ( (GS(I,J),J=1,3), I=1,3 ) WRITE(NOU,106) 106 FORMAT(1X,'BUSING-LEVY B MATRIX') WRITE(NOU,104) ( (BM(I,J),J=1,3), I=1,3) WRITE(NOU,107) CVR 107 FORMAT(/1X,'REAL CELL VOL(A**3): ',E15.6) WRITE(NOU,108) CVS 108 FORMAT(1X,'RECL. CELL VOL(A**-3): ',E12.6) RETURN END C..... C .. SUBROUTINE LREAL C..... C .. GET THE LENGTH OF A VECTOR GIVEN C .. ITS REAL SPACE COMPONENTS C..... SUBROUTINE LREAL(FX,FY,FZ, DREL) DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM C- S1 = GR(1,1)*FX + GR(1,2)*FY + GR(1,3)*FZ S2 = GR(2,1)*FX + GR(2,2)*FY + GR(2,3)*FZ S3 = GR(3,1)*FX + GR(3,2)*FY + GR(3,3)*FZ D2 = FX*S1 + FY*S2 + FZ*S3 IF (D2.LE.0.) D2 = 0. DREL = SQRT(D2) RETURN END C..... C .. SUBROUTINE LSTAR C..... C .. CALCULATE THE LENGTH OF A RECIPROCAL VECTOR C..... SUBROUTINE LSTAR(FH,FK,FL, DSTR) DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM R1 = GS(1,1)*FH + GS(1,2)*FK + GS(1,3)*FL R2 = GS(2,1)*FH + GS(2,2)*FK + GS(2,3)*FL R3 = GS(3,1)*FH + GS(3,2)*FK + GS(3,3)*FL D2 = FH*R1 + FK*R2 + FL*R3 IF (D2.LE.0.) D2 = 0. DSTR = SQRT(D2) RETURN END C..... C .. SUBROUTINE EDGE1 C..... C .. GENERATE THE EDGES OF A CRYSTAL C..... SUBROUTINE EDGE1 DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DIMENSION VR(150,3), MM(150,3) COMMON /VERTEX/ NV, VR, MM COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU CHARACTER*1 STAR, BLNK, LABL, LABLS(250) DIMENSION JS(250), ELENS(250) DATA STAR/'*'/, BLNK/' '/ C- WRITE(NOU,9001) 9001 FORMAT(/1X,'DIAGONALS AND EDGES OF CRYSTAL IN MILLIMETERS') DO 200 I=1,NV-1 MI1 = MM(I,1) MI2 = MM(I,2) MI3 = MM(I,3) NS = 0 DO 150 J=I+1,NV NS = NS+1 MJ1 = MM(J,1) MJ2 = MM(J,2) MJ3 = MM(J,3) CALL TWOCO (MI1,MI2,MI3, MJ1,MJ2,MJ3, NN) IF (NN.EQ.1) LABL = STAR IF (NN.EQ.0) LABL = BLNK DELA = VR(I,1) - VR(J,1) DELB = VR(I,2) - VR(J,2) DELC = VR(I,3) - VR(J,3) CALL LREAL (DELA,DELB,DELC, ELEN) ELENS(NS) = ELEN LABLS(NS) = LABL JS(NS) = J 150 CONTINUE WRITE(NOU,9003) I,(JS(N),ELENS(N),LABLS(N),N=1,NS) 9003 FORMAT(1X,'VERTEX # ',I2,' TO',/, * 5(' #',I2,F7.3,A)) 200 CONTINUE C- RETURN END C..... C .. SUBROUTINE EMUX2 C..... C .. GIVEN THE CHEMICAL COMPOSITION OF THE CRYSTAL CHEMICAL UNIT C .. AND THE UNIT CELL VOLUME, CALCULATE THE MOLECULAR WEIGHT OF C .. THE CRYSTAL CHEMICAL UNIT, THE DENSITY OF THE CRYSTAL, AND C .. THE LINEAR ABSORPTION COEFFICIENT OF THE CRYSTAL. THE VALUES C .. OF THE MASS ABSORPTION COEFFICIENTS FOR ELEMENTS 1 TO 96 ARE C .. TAKEN FROM THE INTERNATIONAL TABLES VOLUME IV AND ARE TABULATED C .. FOR COPPER, MOLYBDENUM AND SILVER RADIATION (BOTH K-ALPHA AND C .. K-BETA). THE VALUES FOR THE ATOMIC WEIGHTS ARE FROM THE SARGENT C .. AND WELCH SCIENTIFIC COMPANY "PERIODIC TABLE OF THE ELEMENTS" C .. COPYRIGHT 1968 (WITH THE ATOMIC WEIGHT OF LITHIUM CORRECTED). C .. VALUES OF THE UNIT CELL VOLUME COME FROM SUBROUTINE CELL2. C..... SUBROUTINE EMUX2 CHARACTER*4 WAVE COMMON /INPT04/ WAVE COMMON /LAMBDA/ FLAM CHARACTER*2 ATOMS(10) COMMON /INPT05/ ATOMS DIMENSION STOIC(10) COMMON /INPT06/ ZCELL, STOIC DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM C- COMMON /CRYSTL/ TMOLW, RHO, EMUX CHARACTER*2 SYMBL,BLANK,ATLAB(96) CHARACTER*4 CUKA,CUKB, MOKA,MOKB, AGKA,AGKB DIMENSION ECUKA(96), ECUKB(96) DIMENSION EMOKA(96), EMOKB(96) DIMENSION EAGKA(96), EAGKB(96) DIMENSION ATWHT(96), EWORK(96) C- DATA CUKA/'CUKA'/, CUKB/'CUKB'/ DATA MOKA/'MOKA'/, MOKB/'MOKB'/ DATA AGKA/'AGKA'/, AGKB/'AGKB'/ DATA ATLAB/'H ' ,'HE' ,'LI' ,'BE' ,'B ' ,'C ' , * 'N ' ,'O ' ,'F ' ,'NE' ,'NA' ,'MG' , * 'AL' ,'SI' ,'P ' ,'S ' ,'CL' ,'AR' , * 'K ' ,'CA' ,'SC' ,'TI' ,'V ' ,'CR' , * 'MN' ,'FE' ,'CO' ,'NI' ,'CU' ,'ZN' , * 'GA' ,'GE' ,'AS' ,'SE' ,'BR' ,'KR' , * 'RB' ,'SR' ,'Y ' ,'ZR' ,'NB' ,'MO' , * 'TC' ,'RU' ,'RH' ,'PD' ,'AG' ,'CD' , * 'IN' ,'SN' ,'SB' ,'TE' ,'I ' ,'XE' , * 'CS' ,'BA' ,'LA' ,'CE' ,'PR' ,'ND' , * 'PM' ,'SM' ,'EU' ,'GD' ,'TB' ,'DY' , * 'HO' ,'ER' ,'TM' ,'YB' ,'LU' ,'HF' , * 'TA' ,'W ' ,'RE' ,'OS' ,'IR' ,'PT' , * 'AU' ,'HG' ,'TL' ,'PB' ,'BI' ,' ' , * ' ' ,'RA' ,' ' ,' ' ,' ' ,'TH' , * ' ' ,'U ' ,' ' ,'PU' ,' ' ,' ' / DATA AVOG /1.660565/ DATA ECUKA/ 0.3912, 0.2835, 0.4770, 1.0070, 2.1420, 4.2190, * 7.1420, 11.0300, 15.9500, 22.1300, 30.3000, 40.8800, * 50.2300, 65.3200, 77.2800, 92.5300,109.2000,119.5000, * 148.4000,171.4000,186.0000,202.4000,222.6000,252.3000, * 272.5000,304.4000,338.6000, 48.8300, 51.5400, 59.5100, * 62.1300, 67.9200, 75.6500, 82.8900, 90.2900, 97.0200, * 106.3000,115.3000,127.1000,136.8000,148.8000,158.3000, * 167.7000,180.8000,194.1000,205.0000,218.1000,229.3000, * 242.1000,253.3000,266.5000,273.4000,291.7000,309.8000, * 325.4000,336.1000,353.5000,378.8000,402.2000,417.9000, * 441.1000,453.5000,417.9000,426.7000,321.9000,336.6000, * 128.4000,134.3000,140.2000,144.7000,152.0000,157.7000, * 161.5000,170.5000,178.3000,183.8000,192.2000,198.2000, * 207.8000,216.2000,222.2000,232.1000,242.9000, 0.0000, * 0.0000,263.7000, 0.0000, 0.0000, 0.0000,306.8000, * 0.0000,305.7000, 0.0000,352.9000, 0.0000, 0.0000/ DATA ECUKB/ 0.3882, 0.2623, 0.3939, 0.7742, 1.5900, 3.0930, * 5.2150, 8.0620, 11.6600, 16.2400, 22.2300, 30.0800, * 37.1400, 48.3700, 57.4400, 68.9000, 81.7900, 89.3400, * 111.7000,129.0000,140.8000,153.2000,168.0000,191.1000, * 206.7000,233.6000,258.7000,282.8000, 38.7400, 45.3000, * 46.6500, 51.4400, 57.0100, 62.3200, 68.0700, 73.2200, * 80.1600, 86.7700, 96.1900,103.3000,112.3000,119.7000, * 126.9000,137.0000,147.1000,155.6000,165.8000,174.0000, * 183.3000,193.1000,203.6000,208.4000,221.9000,235.9000, * 247.5000,256.0000,270.8000,289.6000,306.8000,319.8000, * 336.4000,346.6000,369.0000,377.2000,399.9000,360.2000, * 272.4000,291.7000,288.5000,110.9000,116.5000,121.0000, * 123.9000,131.5000,137.3000,141.3000,147.6000,151.2000, * 160.6000,166.6000,171.1000,178.6000,187.5000, 0.0000, * 0.0000,203.0000, 0.0000, 0.0000, 0.0000,236.6000, * 0.0000,236.2000, 0.0000,271.2000, 0.0000, 0.0000/ DATA EMOKA/ 0.3727, 0.2019, 0.1968, 0.2451, 0.3451, 0.5348, * 0.7898, 1.1470, 1.5840, 2.2090, 2.9390, 3.9790, * 5.0430, 6.5330, 7.8700, 9.6250, 11.6400, 12.6200, * 16.2000, 19.0000, 21.0400, 23.2500, 25.2400, 29.2500, * 31.8600, 37.7400, 41.0200, 47.2400, 49.3400, 55.4600, * 56.9000, 60.4700, 65.9700, 68.8200, 74.6800, 79.1000, * 83.0000, 88.0400, 97.5600, 16.1000, 16.9600, 18.4400, * 19.7800, 21.3300, 23.0500, 24.4200, 26.3800, 27.7300, * 29.1300, 31.1800, 33.0100, 33.9200, 36.3300, 38.3100, * 40.4400, 42.3700, 45.3400, 48.5600, 50.7800, 53.2800, * 55.5200, 57.9600, 61.1800, 62.7900, 66.7700, 68.8900, * 72.1400, 75.6100, 78.9800, 80.2300, 84.1800, 86.3300, * 89.5100, 95.7600, 98.7400,100.2000,103.4000,108.6000, * 111.3000,114.7000,119.4000,122.8000,125.9000, 0.0000, * 0.0000,117.2000, 0.0000, 0.0000, 0.0000, 99.4600, * 0.0000, 96.6700, 0.0000, 48.8400, 0.0000, 0.0000/ DATA EMOKB/ 0.3699, 0.1972, 0.1866, 0.2216, 0.2928, 0.4285, * 0.6054, 0.8545, 1.1540, 1.5970, 2.0980, 2.8250, * 3.5850, 4.6240, 5.5690, 6.8350, 8.2610, 8.9490, * 11.5100, 13.5600, 15.0000, 16.6500, 18.0700, 20.9900, * 22.8900, 27.2100, 29.5100, 34.1800, 35.7700, 40.2600, * 41.6900, 44.2600, 48.5700, 51.2000, 55.5600, 58.6400, * 62.0700, 65.5900, 72.5700, 75.2000, 81.2200, 13.2900, * 14.3000, 15.4000, 16.6500, 17.6300, 19.1000, 20.1300, * 21.1800, 22.6200, 23.9100, 24.6700, 26.5300, 27.8600, * 29.5100, 31.0000, 33.1000, 35.5400, 37.0900, 38.8800, * 40.5200, 42.4000, 44.7400, 45.9500, 48.8800, 50.3800, * 52.7600, 55.0700, 57.9400, 59.2200, 62.0400, 64.1500, * 66.0700, 70.5700, 72.4700, 74.1300, 77.2000, 80.2300, * 82.3300, 85.3000, 88.2500, 90.5500, 93.5000, 0.0000, * 0.0000,100.7000, 0.0000, 0.0000, 0.0000, 73.3400, * 0.0000, 72.6300, 0.0000, 78.9900, 0.0000, 0.0000/ DATA EAGKA/ 0.3663, 0.1931, 0.1788, 0.2047, 0.2558, 0.3537, * 0.4756, 0.6479, 0.8507, 1.1630, 1.5010, 2.0040, * 2.5400, 3.2540, 3.9110, 4.8160, 5.8030, 6.2800, * 8.0800, 9.5730, 10.5400, 11.7600, 12.7800, 14.8800, * 16.2300, 19.3100, 20.9200, 24.3200, 25.5200, 28.7200, * 29.9300, 31.7800, 35.0000, 37.2300, 40.4100, 42.5300, * 45.3300, 47.8900, 52.9600, 55.2400, 59.6300, 62.5300, * 64.9900, 11.0000, 11.9100, 12.5800, 13.6700, 14.4500, * 15.2500, 16.1900, 17.0500, 17.7200, 19.1800, 20.0100, * 21.3100, 22.4500, 23.8500, 25.7200, 26.8000, 28.0200, * 29.2500, 30.6900, 32.4000, 33.3000, 35.4200, 36.5200, * 38.3100, 39.8300, 42.1200, 43.2100, 45.2700, 46.9300, * 48.2500, 51.4000, 52.7800, 54.2400, 56.7500, 58.5800, * 60.2100, 62.6300, 64.5600, 66.1400, 68.7300, 0.0000, * 0.0000, 74.1400, 0.0000, 0.0000, 0.0000, 86.2900, * 0.0000, 86.9200, 0.0000, 58.7600, 0.0000, 0.0000/ DATA EAGKB/ 0.3626, 0.1894, 0.1726, 0.1920, 0.2294, 0.3011, * 0.3851, 0.5041, 0.6401, 0.8605, 1.0880, 1.4360, * 1.8120, 2.3010, 2.7550, 3.4010, 4.0790, 4.4090, * 5.6670, 6.7520, 7.4020, 8.2980, 9.0240, 10.5200, * 11.4900, 13.6700, 14.7900, 17.2500, 18.1500, 20.4200, * 21.4100, 22.7300, 25.1100, 26.9400, 29.2400, 30.7100, * 32.9400, 34.8100, 38.4900, 40.4000, 43.6200, 45.7300, * 47.7300, 50.3400, 52.0300, 57.0300, 9.7800, 10.3700, * 10.9700, 11.5700, 12.1400, 12.7000, 13.8500, 14.3600, * 15.3700, 16.2500, 17.1500, 18.5900, 19.3400, 20.1600, * 21.0800, 22.1800, 23.4300, 24.1000, 25.6300, 26.4400, * 27.7900, 28.7700, 30.5800, 31.4600, 32.9800, 34.2400, * 35.1700, 37.3500, 38.3800, 39.6000, 74.5900, 42.6800, * 43.9300, 45.8700, 47.1500, 48.2200, 50.4100, 0.0000, * 0.0000, 54.4400, 0.0000, 0.0000, 0.0000, 62.8800, * 0.0000, 64.2400, 0.0000, 69.7500, 0.0000, 0.0000/ DATA ATWHT/ 1.0080, 4.0026, 6.9410, 9.0122, 10.8110, 12.0112, * 14.0067, 15.9994, 18.9984, 20.1830, 22.9898, 24.3120, * 26.9815, 28.0860, 30.9738, 32.0640, 35.4530, 39.9480, * 39.1020, 40.0800, 44.9560, 47.9000, 50.9420, 51.9960, * 54.9380, 55.8430, 58.9330, 58.7100, 63.5400, 65.3700, * 69.7200, 72.5900, 74.9220, 78.9600, 79.9090, 83.8000, * 85.4700, 87.6200, 88.9050, 91.2200, 92.9060, 95.9400, * 98.0000,101.0700,102.9050,106.4000,107.8700,112.4000, * 114.8200,118.6900,121.7500,127.6000,126.9040,131.3000, * 132.9050,137.3400,138.9100,140.1200,140.9070,144.2400, * 147.0000,150.3500,151.9600,157.2500,158.9240,162.5000, * 164.9300,167.2600,168.9340,173.0400,174.9700,178.4900, * 180.9480,183.8500,186.2000,190.2000,192.2000,195.0900, * 196.9670,200.5900,204.3700,207.1900,208.9800,000.0000, * 000.0000,222.0000,000.0000,000.0000,000.0000,232.0380, * 000.0000,238.0300,000.0000,242.0000,000.0000,000.0000/ DATA WCUKA/1.5418/, WCUKB/1.3922/ DATA WMOKA/0.7107/, WMOKB/0.6323/ DATA WAGKA/0.5608/, WAGKB/0.4970/ DATA EWORK/96*0.0/ C..... C .. NO CALCULATIONS IF NO ATOMS SPECIFIED C..... IF (STOIC(1).LE.0.0) RETURN C..... C- GET THE MASS-ABSORPTION COEFFICIENTS AND WAVELENGTHS APPROPRIATE C- TO THE RADIATION USED IF (WAVE.EQ.CUKA) GO TO 100 IF (WAVE.EQ.CUKB) GO TO 110 IF (WAVE.EQ.MOKA) GO TO 120 IF (WAVE.EQ.MOKB) GO TO 130 IF (WAVE.EQ.AGKA) GO TO 140 IF (WAVE.EQ.AGKB) GO TO 150 GO TO 160 100 FLAM = WCUKA DO 101 J=1,96 EWORK(J) = ECUKA(J) 101 CONTINUE GO TO 200 110 FLAM = WCUKB DO 111 J=1,96 EWORK(J) = ECUKB(J) 111 CONTINUE GO TO 200 120 FLAM = WMOKA DO 121 J=1,96 EWORK(J) = EMOKA(J) 121 CONTINUE GO TO 200 130 FLAM = WMOKB DO 131 J=1,96 EWORK(J) = EMOKB(J) 131 CONTINUE GO TO 200 140 FLAM = WAGKA DO 141 J=1,96 EWORK(J) = EAGKA(J) 141 CONTINUE GO TO 200 150 FLAM = WAGKB DO 151 J=1,96 EWORK(J) = EAGKB(J) 151 CONTINUE GO TO 200 160 DO 161 J=1,96 EWORK(J) = 0.0 161 CONTINUE 200 CONTINUE C- ACCUMULATE THE SUMS NECESSARY SUM1 = 0. SUM2 = 0. DO 500 N=1,10 SYMBL = ATOMS(N) FRACT = STOIC(N) DO 450 J=1,96 IF( SYMBL.NE.ATLAB(J) ) GO TO 450 ATW = ATWHT(J) FMU = EWORK(J) 450 CONTINUE FAT = FRACT*ATW SUM1 = SUM1 + FAT SUM2 = SUM2 + FAT*FMU 500 CONTINUE C- GET THE VARIOUS DERIVED DATA TMOLW = SUM1 RHO = AVOG*ZCELL*SUM1/CVR IF( SUM2.EQ.0. ) RETURN EMUX = RHO*SUM2/SUM1 C- CONVERT EMUX TO 1/MM EMUX = EMUX/10. RETURN END C..... C .. SUBROUTINE EMUX3 C..... C .. WRITE OUT THE LINEAR ABSORPTION COEFFICIENT AND C .. RELATED PERTINENT DATA C..... SUBROUTINE EMUX3 CHARACTER*4 WAVE COMMON /INPT04/ WAVE COMMON /LAMBDA/ FLAM CHARACTER*2 ATOMS(10) COMMON /INPT05/ ATOMS DIMENSION STOIC(10) COMMON /INPT06/ ZCELL,STOIC COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU COMMON /CRYSTL/ TMOLW, RHO, EMUX C- NTYPE = 0 DO 10 J=1,10 IF(STOIC(J).NE.0.0) NTYPE=NTYPE+1 10 CONTINUE WRITE(NOU,200) ZCELL 200 FORMAT(/1X,'NUMBER OF C.C.U./UNIT CELL: ',F5.1) WRITE(NOU,201) 201 FORMAT(1X,'CONTENTS OF C.C.U.') DO 700 J=1,NTYPE WRITE(NOU,202) ATOMS(J), STOIC(J) 202 FORMAT(1X,A,3X,F5.1) 700 CONTINUE WRITE(NOU,203) TMOLW 203 FORMAT(1X,'MOLECULAR WEIGHT: ',F7.2) WRITE(NOU,204) RHO 204 FORMAT(1X,'RHO (MG/MM**3): ',F7.4) WRITE(NOU,205) WAVE 205 FORMAT(1X,'WAVE SYMBOL: ',2X,A) WRITE(NOU,206) FLAM 206 FORMAT(1X,'LAMBDA (A): ',F6.4) WRITE(NOU,207) EMUX 207 FORMAT(1X,'CRYSTAL MU (1/MM): ',F8.5) RETURN END C..... C .. SUBROUTINE FACE2 C..... C .. GIVEN THE MILLERIAN SYMBOLS OF THE CRYSTAL PLANES C .. GENERATE THE UNIT VECTORS IN THE [STAR] FRAME C .. WHICH ARE PERPENDICULAR TO THE PLANES. C..... C .. FS(N,1)...A* COMPONENT OF UNIT VECTOR C .. FS(N,2)...B* " " " " C .. FS(N,3)...C* " " " " C .. FS(N,4)...PERPENDICULAR DISTANCE IN MM. C .. FS(N,5)...DIST TO FACE FOR SOLVENT CALC C..... SUBROUTINE FACE2 DIMENSION FACES(50,4) COMMON /INPT07/ FACES DIMENSION FS(50,5) COMMON /PLANES/ NF, FS C- NF = 0 J = 1 10 FH = FACES(J,1) FK = FACES(J,2) FL = FACES(J,3) PD = FACES(J,4) CALL LSTAR(FH,FK,FL, DSTR) IF( DSTR.LE.0. ) RETURN C- FS(J,1) = FH/DSTR FS(J,2) = FK/DSTR FS(J,3) = FL/DSTR FS(J,4) = PD FS(J,5) = 0.0 NF = NF + 1 J = J + 1 IF( J .LE. 50 ) GO TO 10 RETURN END C..... C .. SUBROUTINE FACE3 C..... C .. PRINT OR TYPE OUT THE CRYSTAL FACE INFORMATION C..... SUBROUTINE FACE3 DIMENSION FACES(50,4) COMMON /INPT07/ FACES DIMENSION FS(50,5) COMMON /PLANES/ NF, FS COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C- WRITE(NOU,9001) 9001 FORMAT(/1X,'CRYSTAL FACE DATA') WRITE(NOU,9002) 9002 FORMAT(1X,' #',7X,'MILLER INDICES',5X, *' LENGTH',5X,'[STAR] COMPONENTS OF UNIT VECTOR') DO 10 I=1,NF WRITE(NOU,9003) I,(FACES(I,J),J=1,4),(FS(I,J),J=1,3) 9003 FORMAT(1X,I3,3F8.3,2X,F8.5,1X,3E13.4) 10 CONTINUE RETURN END C..... C .. SUBROUTINE FILS2 C..... C .. OUTPUT THE NAMES OF THE INPUT AND OUTPUT FILES C..... SUBROUTINE FILS2 CHARACTER FIL01*80, FIL02*80 COMMON /INPT02/ FIL01, FIL02, IPATH COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C.... WRITE(NOU,9001) FIL01 WRITE(NOU,9002) FIL02 WRITE(NOU,9003) IPATH 9001 FORMAT(/1X,'INPUT FILE NAME: ',A) 9002 FORMAT( 1X,'OUTPUT FILE NAME: ',A) 9003 FORMAT(/1X, &'IPATH = ',I1//1X, &'IPATH = 0 DO NOT'/1X, &' = 1 DO'/1X, &' WRITE TBAR AND S0 AND S1 COMPONENTS TO OUTPUT FILE') RETURN END C..... C .. SUBROUTINE FLIPS C..... C .. GENERATE THE DUMMY FACE VECTOR TO BE C .. USED IN THE SECTIONING OF THE CRYSTAL. C .. IF THE DISTANCE IS NEGATIVE THEN FLIP C .. THE SIGNS OF THE INDICES OF THE PLANE. C .. THAT IS, AN HKL PLANE WITH A NEGATIVE C .. DISTANCE IS ACTUALLY AN -H-K-L PLANE C .. WITH A POSITIVE DISTANCE. C..... SUBROUTINE FLIPS (DVEC,EPSL, FVEC) DIMENSION DVEC(3), FVEC(4) C- FVEC(1) = DVEC(1) FVEC(2) = DVEC(2) FVEC(3) = DVEC(3) FVEC(4) = EPSL IF (EPSL.GE.0.0) RETURN DO 10 N=1,4 FVEC(N) = -FVEC(N) 10 CONTINUE RETURN END C..... C .. SUBROUTINE GAUTAB C..... C .. GIVEN NGP, THE GAUSSIAN GRID NUMBER, TRANSFER THE C .. GAUSSIAN FRACTIONS AND GAUSSIAN WEIGHTS TO WORKING C .. ARRAYS PASSED BACK TO CALLING ROUTINE. THE VALUES C .. IN THESE TABLES ARE CALCULATED FROM THE TABLES OF C .. ABRAMOWITZ AND STEGUN. THE TABLES WERE PUNCHED BY C .. THAT STAR OF THE STAGE, OUR OWN CATHERINE DEVINE, C .. AND CONVERTED TO THE PRESENT FORM BY GDT. C..... SUBROUTINE GAUTAB (NGP, FGAUSS,WGAUSS) DIMENSION FGAUSS(96), WGAUSS(96) DIMENSION NALLOW(19) DIMENSION F02( 2), W02( 2) DIMENSION F03( 3), W03( 3) DIMENSION F04( 4), W04( 4) DIMENSION F05( 5), W05( 5) DIMENSION F06( 6), W06( 6) DIMENSION F07( 7), W07( 7) DIMENSION F08( 8), W08( 8) DIMENSION F09( 9), W09( 9) DIMENSION F10(10), W10(10) DIMENSION F12(12), W12(12) DIMENSION F16(16), W16(16) DIMENSION F20(20), W20(20) DIMENSION F24(24), W24(24) DIMENSION F32(32), W32(32) DIMENSION F40(40), W40(40) DIMENSION F48(48), W48(48) DIMENSION F64(64), W64(64) DIMENSION F80(80), W80(80) DIMENSION F96(96), W96(96) C NGAUSS = 2 DATA F02/.2113248,.7886752/ DATA W02/.5000000,.5000000/ C NGAUSS = 3 DATA F03/.1127017,.5000000,.8872983/ DATA W03/.2777778,.4444444,.2777778/ C NGAUSS = 4 DATA F04/.0694318,.3300095,.6699905,.9305682/ DATA W04/.1739274,.3260726,.3260726,.1739274/ C NGAUSS = 5 DATA F05/.0469101,.2307653,.5000000,.7692347,.9530899/ DATA W05/.1184634,.2393143,.2844445,.2393143,.1184634/ C NGAUSS = 6 DATA F06/.0337653,.1693953,.3806904,.6193096,.8306047,.9662347/ DATA W06/.0856623,.1803808,.2339569,.2339569,.1803808,.0856623/ C NGAUSS = 7 DATA F07/.0254461,.1292344,.2970774,.5000000,.7029226,.8707656, * .9745539/ DATA W07/.0647425,.1398527,.1909150,.2089796,.1909150,.1398527, * .0647425/ C NGAUSS = 8 DATA F08/.0198551,.1016667,.2372338,.4082827,.5917173,.7627662, * .8983333,.9801450/ DATA W08/.0506142,.1111905,.1568533,.1813419,.1813419,.1568533, * .1111905,.0506142/ C NGAUSS = 9 DATA F09/.0159199,.0819845,.1933143,.3378733,.5000000,.6621267, * .8066857,.9180155,.9840801/ DATA W09/.0406372,.0903241,.1303053,.1561735,.1651197,.1561735, * .1303053,.0903241,.0406372/ C NGAUSS = 10 DATA F10/.0130467,.0674683,.1602952,.2833023,.4255629,.5744371, * .7166977,.8397048,.9325317,.9869533/ DATA W10/.0333356,.0747257,.1095432,.1346333,.1477621,.1477621, * .1346333,.1095432,.0747257,.0333356/ C NGAUSS = 12 DATA F12/.0092197,.0479413,.1150486,.2063410,.3160843,.4373833, * .5626167,.6839157,.7936590,.8849514,.9520587,.9907803/ DATA W12/.0235876,.0534697,.0800392,.1015837,.1167462,.1245735, * .1245735,.1167462,.1015837,.0800392,.0534697,.0235876/ C NGAUSS = 16 DATA F16/.0052995,.0277125,.0671844,.1222978,.1910619,.2709916, * .3591982,.4524938,.5475063,.6408018,.7290084,.8089381, * .8777022,.9328156,.9722875,.9947005/ DATA W16/.0135763,.0311268,.0475793,.0623145,.0747980,.0845783, * .0913017,.0947253,.0947253,.0913017,.0845783,.0747980, * .0623145,.0475793,.0311268,.0135763/ C NGAUSS = 20 DATA F20/.0034357,.0180140,.0438828,.0804415,.1268341,.1819731, * .2445665,.3131469,.3861071,.4617367,.5382633,.6138930, * .6868531,.7554335,.8180269,.8731660,.9195585,.9561172, * .9819860,.9965643/ DATA W20/.0088070,.0203007,.0313360,.0416383,.0509650,.0590972, * .0658443,.0710481,.0745865,.0763767,.0763767,.0745865, * .0710481,.0658443,.0590972,.0509650,.0416383,.0313360, * .0203007,.0088070/ C NGAUSS = 24 DATA F24/.0024064,.0126357,.0308627,.0567923,.0899990,.1299379, * .1759531,.2272892,.2831033,.3424787,.4044406,.4679716, * .5320284,.5955595,.6575214,.7168968,.7727108,.8240469, * .8700621,.9100010,.9432077,.9691373,.9873644,.9975936/ DATA W24/.0061706,.0142657,.0221387,.0296493,.0366732,.0430951, * .0488093,.0537221,.0577528,.0608352,.0629188,.0639691, * .0639691,.0629188,.0608352,.0577528,.0537221,.0488093, * .0430951,.0366732,.0296493,.0221387,.0142657,.0061706/ C NGAUSS = 32 DATA F32/.0013680,.0071943,.0176188,.0325469,.0518394,.0753162, * .1027581,.1339090,.1684778,.2061421,.2465501,.2893244, * .3340657,.3803563,.4277640,.4758461,.5241538,.5722360, * .6196437,.6659343,.7106757,.7534500,.7938579,.8315222, * .8660911,.8972420,.9246838,.9481606,.9674531,.9823812, * .9928058,.9986320/ DATA W32/.0035093,.0081372,.0126961,.0171370,.0214180,.0254990, * .0293420,.0329111,.0361729,.0390969,.0416560,.0438260, * .0455869,.0469222,.0478193,.0482700,.0482700,.0478193, * .0469222,.0455869,.0438260,.0416560,.0390969,.0361729, * .0329111,.0293420,.0254990,.0214180,.0171370,.0126961, * .0081372,.0035093/ C NGAUSS = 40 DATA F40/.0008811,.0046369,.0113701,.0210416,.0335936,.0489506, * .0670202,.0876939,.1108471,.1363408,.1640216,.1937231, * .2252665,.2584621,.2931104,.3290030,.3659239,.4036512, * .4419580,.4806138,.5193862,.5580420,.5963488,.6340761, * .6709970,.7068896,.7415379,.7747335,.8062770,.8359784, * .8636592,.8891529,.9123061,.9329798,.9510494,.9664064, * .9789584,.9886299,.9953631,.9991189/ DATA W40/.0022607,.0052491,.0082106,.0111229,.0139685,.0167301, * .0193911,.0219354,.0243479,.0266139,.0287199,.0306531, * .0324020,.0339560,.0353058,.0364433,.0373616,.0380552, * .0385199,.0387530,.0387530,.0385199,.0380552,.0373616, * .0364433,.0353058,.0339560,.0324020,.0306531,.0287199, * .0266139,.0243479,.0219354,.0193911,.0167301,.0139685, * .0111229,.0082106,.0052491,.0022607/ C NGAUSS = 48 DATA F48/.0006145,.0032349,.0079377,.0147042,.0235062,.0343066, * .0470605,.0617140,.0782059,.0964669,.1164205,.1379830, * .1610638,.1855663,.2113877,.2384195,.2665485,.2956568, * .3256221,.3563188,.3876181,.4193888,.4514976,.4838099, * .5161901,.5485023,.5806112,.6123819,.6436812,.6743780, * .7043433,.7334515,.7615805,.7886124,.8144338,.8389362, * .8620170,.8835795,.9035331,.9217942,.9382860,.9529396, * .9656934,.9764938,.9852958,.9920623,.9967651,.9993855/ DATA W48/.0015766,.0036638,.0057386,.0077897,.0098081,.0117854, * .0137132,.0155836,.0173886,.0191207,.0207726,.0223373, * .0238084,.0251795,.0264451,.0275998,.0286386,.0295574, * .0303522,.0310197,.0315571,.0319621,.0322331,.0323688, * .0323688,.0322331,.0319621,.0315571,.0310197,.0303522, * .0295574,.0286386,.0275998,.0264451,.0251795,.0238084, * .0223373,.0207726,.0191207,.0173886,.0155836,.0137132, * .0117854,.0098081,.0077897,.0057386,.0036638,.0015766/ C NGAUSS = 64 DATA F64/.0003475,.0018300,.0044933,.0083318,.0133366,.0194956, * .0267943,.0352154,.0447389,.0553423,.0670003,.0796853, * .0933673,.1080138,.1235901,.1400591,.1573819,.1755172, * .1944223,.2140522,.2343602,.2552985,.2768170,.2988649, * .3213899,.3443386,.3676564,.3912882,.4151778,.4392686, * .4635035,.4878249,.5121751,.5364966,.5607314,.5848222, * .6087118,.6323436,.6556615,.6786101,.7011351,.7231830, * .7447016,.7656398,.7859478,.8055777,.8244828,.8426182, * .8599409,.8764099,.8919863,.9066327,.9203147,.9329997, * .9446577,.9552611,.9647846,.9732057,.9805045,.9866634, * .9916682,.9955067,.9981701,.9996525/ DATA W64/.0008916,.0020735,.0032523,.0044234,.0055840,.0067315, * .0078630,.0089759,.0100687,.0111351,.0121763,.0131888, * .0141699,.0151174,.0160289,.0169026,.0177361,.0185276, * .0192751,.0199769,.0206313,.0212367,.0217918,.0222953, * .0227458,.0231424,.0234841,.0237701,.0239997,.0241724, * .0242878,.0243455,.0243455,.0242878,.0241724,.0239997, * .0237701,.0234841,.0231424,.0227458,.0222953,.0217918, * .0212367,.0206313,.0199769,.0192751,.0185276,.0177361, * .0169026,.0160289,.0151174,.0141699,.0131888,.0121763, * .0111351,.0100687,.0089759,.0078630,.0067315,.0055840, * .0044234,.0032523,.0020735,.0008916/ C NGAUSS = 80 DATA F80/.0002231,.0011750,.0028862,.0053543,.0085757,.0125454, * .0172575,.0227046,.0288786,.0357700,.0433685,.0516622, * .0606387,.0702843,.0805842,.0915230,.1030836,.1152488, * .1279998,.1413174,.1551812,.1695700,.1844621,.1998347, * .2156643,.2319270,.2485979,.2656517,.2830623,.3008033, * .3188476,.3371678,.3557360,.3745238,.3935027,.4126438, * .4319180,.4512958,.4707478,.4902443,.5097557,.5292522, * .5487042,.5680820,.5873562,.6064972,.6254762,.6442640, * .6628322,.6811524,.6991967,.7169377,.7343483,.7514021, * .7680730,.7843357,.8001653,.8155379,.8304300,.8448188, * .8586826,.8720002,.8847512,.8969164,.9084771,.9194158, * .9297157,.9393613,.9483379,.9566315,.9642300,.9711214, * .9772954,.9827425,.9874546,.9914243,.9946457,.9971138, * .9988250,.9997769/ DATA W80/.0005725,.0013317,.0020902,.0028455,.0035964,.0043420, * .0050809,.0058120,.0065344,.0072467,.0079481,.0086374, * .0093134,.0099753,.0106220,.0112526,.0118660,.0124612, * .0130376,.0135791,.0141299,.0146442,.0151361,.0156051, * .0160503,.0164710,.0168666,.0172365,.0175803,.0178972, * .0181869,.0184489,.0186829,.0188882,.0190648,.0192125, * .0193309,.0194199,.0194792,.0195089,.0195089,.0194792, * .0194199,.0193309,.0192125,.0190648,.0188882,.0186829, * .0184489,.0181869,.0178972,.0175803,.0172365,.0168666, * .0164710,.0160503,.0156051,.0151361,.0146442,.0141299, * .0135791,.0130376,.0124612,.0118660,.0112526,.0106220, * .0099753,.0093134,.0086374,.0079481,.0072467,.0065344, * .0058120,.0050809,.0043420,.0035964,.0028455,.0020902, * .0013317,.0005725/ C NGAUSS = 96 DATA F96/.0001552,.0008178,.0020091,.0037281,.0059730,.0087413, * .0120304,.0158366,.0201558,.0249836,.0303149,.0361437, * .0424643,.0492697,.0580527,.0643058,.0725205,.0811883, * .0902998,.0998456,.1098155,.1201988,.1309847,.1421616, * .1537178,.1656409,.1779183,.1905371,.2034838,.2167448, * .2303059,.2441529,.2582710,.2726453,.2872605,.3021012, * .3171516,.3323958,.3478176,.3634006,.3791284,.3949844, * .4109516,.4270132,.4431520,.4593512,.4755935,.4918616, * .5081384,.5244065,.5406488,.5568479,.5729868,.5890484, * .6050156,.6208716,.6365994,.6521825,.6676043,.6828485, * .6978988,.7127395,.7273547,.7417290,.7558471,.7696941, * .7832552,.7965162,.8094629,.8220817,.8343592,.8462822, * .8578384,.8690153,.8798012,.8901845,.9001544,.9097002, * .9188117,.9274795,.9356943,.9419473,.9507303,.9575357, * .9638563,.9696851,.9750164,.9798442,.9841634,.9879696, * .9912587,.9940271,.9962720,.9979909,.9991822,.9998448/ DATA W96/.0003984,.0009270,.0014553,.0019823,.0025071,.0030293, * .0035483,.0040634,.0045743,.0050804,.0055811,.0060758, * .0065641,.0070455,.0075193,.0079853,.0084428,.0088913, * .0093303,.0097595,.0101784,.0105864,.0109833,.0113686, * .0117417,.0121024,.0124503,.0127850,.0131061,.0134135, * .0137065,.0139850,.0142487,.0144973,.0147306,.0149482, * .0151500,.0153357,.0155051,.0156582,.0157946,.0159144, * .0160173,.0161031,.0161719,.0162236,.0162580,.0162753, * .0162753,.0162580,.0162236,.0161719,.0161031,.0160173, * .0159144,.0157946,.0156582,.0155051,.0153357,.0151500, * .0149482,.0147306,.0144973,.0142487,.0139850,.0137065, * .0134135,.0131061,.0127850,.0124503,.0121024,.0117417, * .0113686,.0109833,.0105864,.0101784,.0097595,.0093303, * .0088913,.0084428,.0079853,.0075193,.0070455,.0065641, * .0060758,.0055811,.0050804,.0045743,.0040634,.0035483, * .0030293,.0025071,.0019823,.0014553,.0009270,.0003984/ DATA NALLOW/2,3,4,5,6,7,8,9,10,12,16,20,24,32,40,48,64,80,96/ C- ERROR CONDITIONS REQUIRING RETURN DO 200 J=1,96 FGAUSS(J) = 0. WGAUSS(J) = 0. 200 CONTINUE IF (NGP.LT.2 .OR. NGP.GT.96) RETURN NOK = 0 DO 210 J=1,19 IF( NGP .EQ. NALLOW(J) ) NOK = 1 210 CONTINUE IF( NOK .NE. 1 ) RETURN C- FIND THE PROPER GAUSS POINT TABLE DO 220 J=1,19 IF( NGP .EQ. NALLOW(J) ) NGPNT = J 220 CONTINUE GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)NGPNT C- 1 DO 1000 J=1,2 FGAUSS(J) = F02(J) WGAUSS(J) = W02(J) 1000 CONTINUE RETURN C- 2 DO 2000 J=1,3 FGAUSS(J) = F03(J) WGAUSS(J) = W03(J) 2000 CONTINUE RETURN C- 3 DO 3000 J=1,4 FGAUSS(J) = F04(J) WGAUSS(J) = W04(J) 3000 CONTINUE RETURN C- 4 DO 4000 J=1,5 FGAUSS(J) = F05(J) WGAUSS(J) = W05(J) 4000 CONTINUE RETURN C- 5 DO 5000 J=1,6 FGAUSS(J) = F06(J) WGAUSS(J) = W06(J) 5000 CONTINUE RETURN C- 6 DO 6000 J=1,7 FGAUSS(J) = F07(J) WGAUSS(J) = W07(J) 6000 CONTINUE RETURN C- 7 DO 7000 J=1,8 FGAUSS(J) = F08(J) WGAUSS(J) = W08(J) 7000 CONTINUE RETURN C- 8 DO 8000 J=1,9 FGAUSS(J) = F09(J) WGAUSS(J) = W09(J) 8000 CONTINUE RETURN C- 9 DO 9000 J=1,10 FGAUSS(J) = F10(J) WGAUSS(J) = W10(J) 9000 CONTINUE RETURN C- 10 DO 10000 J=1,12 FGAUSS(J) = F12(J) WGAUSS(J) = W12(J) 10000 CONTINUE RETURN C- 11 DO 11000 J=1,16 FGAUSS(J) = F16(J) WGAUSS(J) = W16(J) 11000 CONTINUE RETURN C- 12 DO 12000 J=1,20 FGAUSS(J) = F20(J) WGAUSS(J) = W20(J) 12000 CONTINUE RETURN C- 13 DO 13000 J=1,24 FGAUSS(J) = F24(J) WGAUSS(J) = W24(J) 13000 CONTINUE RETURN C- 14 DO 14000 J=1,32 FGAUSS(J) = F32(J) WGAUSS(J) = W32(J) 14000 CONTINUE RETURN C- 15 DO 15000 J=1,40 FGAUSS(J) = F40(J) WGAUSS(J) = W40(J) 15000 CONTINUE RETURN C- 16 DO 16000 J=1,48 FGAUSS(J) = F48(J) WGAUSS(J) = W48(J) 16000 CONTINUE RETURN C- 17 DO 17000 J=1,64 FGAUSS(J) = F64(J) WGAUSS(J) = W64(J) 17000 CONTINUE RETURN C- 18 DO 18000 J=1,80 FGAUSS(J) = F80(J) WGAUSS(J) = W80(J) 18000 CONTINUE RETURN C- 19 DO 19000 J=1,96 FGAUSS(J) = F96(J) WGAUSS(J) = W96(J) 19000 CONTINUE RETURN C- END C..... C .. SUBROUTINE GRID2 C..... C .. GENERATE THE GAUSSIAN GRID POINTS AND WEIGHTS. C .. THE POINTS ARE SPECIFIED BY THE COMPONENTS OF C .. VECTORS IN THE [REAL] FRAME AND THE WEIGHTS C .. ARE PROPORTIONAL TO THE VOLUME ELEMENTS SUR- C .. ROUNDING THE POINTS. C .. CALCULATE COORDINATES OF POINTS W/R/T CENTER OF MASS C .. PABC(I,3) W/R/T ABSORPTION ORIGIN OA C .. RABC(I,3) W/R/T CENTER OF MASS C..... SUBROUTINE GRID2 DIMENSION NGP(3), D1(3),D2(3),D3(3) COMMON /INPT08/ NGP, D1,D2,D3 DIMENSION FS(50,5) COMMON /PLANES/ NF, FS PARAMETER (NGMAX=32**3) DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3) COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC DIMENSION F1(4), F2(4), F3(4), POR(3) DIMENSION U1(3), U2(3), U3(3) DIMENSION GF1(96),GW1(96), GF2(96),GW2(96), GF3(96),GW3(96) C- NPTS = 0 VXTL = 0.0 DO 10 I=1,3 VCOM(I) = 0. 10 CONTINUE DO 30 I=1,NGMAX WABC(I) = 0. DO 20 J=1,3 PABC(I,J) = 0. RABC(I,J) = 0. 20 CONTINUE 30 CONTINUE C- IF (NF .LE. 4) RETURN NP1 = NGP(1) NP2 = NGP(2) NP3 = NGP(3) CALL GAUTAB (NP1, GF1,GW1) CALL GAUTAB (NP2, GF2,GW2) CALL GAUTAB (NP3, GF3,GW3) CALL VSTAR (VSTR) CALL UNIVS (D1, U1) CALL UNIVS (D2, U2) CALL UNIVS (D3, U3) C- CALL VERT1 CALL MINMAX (U1, D1MIN,D1MAX) D1DEL = D1MAX - D1MIN C- DO 1000 N1 = 1,NP1 EPS1 = D1MIN + GF1(N1)*D1DEL PIE1 = GW1(N1)*D1DEL CALL FLIPS (U1,EPS1, F1) CALL VERT2 (F1) CALL MINMAX (U2, D2MIN,D2MAX) D2DEL = D2MAX - D2MIN C- DO 900 N2 = 1,NP2 EPS2 = D2MIN + GF2(N2)*D2DEL PIE2 = GW2(N2)*D2DEL CALL FLIPS (U2,EPS2, F2) CALL VERT3 (F1,F2) CALL MINMAX (U3, D3MIN,D3MAX) D3DEL = D3MAX - D3MIN C- DO 800 N3 = 1,NP3 EPS3 = D3MIN + GF3(N3)*D3DEL PIE3 = GW3(N3)*D3DEL CALL FLIPS (U3,EPS3, F3) CALL VERT4 (F1,F2,F3, POR) VOL = PIE1*PIE2*PIE3*VSTR NPTS = NPTS + 1 PABC(NPTS,1) = POR(1) PABC(NPTS,2) = POR(2) PABC(NPTS,3) = POR(3) WABC(NPTS) = VOL VXTL = VXTL + VOL DO 700 I = 1,3 VCOM(I) = VCOM(I) + VOL*POR(I) 700 CONTINUE 800 CONTINUE 900 CONTINUE 1000 CONTINUE C- DO 1100 I =1,3 VCOM(I) = VCOM(I)/VXTL 1100 CONTINUE DO 1200 I =1,NPTS WABC(I) = WABC(I)/VXTL 1200 CONTINUE DO 1300 I =1,NPTS DO 1300 J =1,3 RABC(I,J) = PABC(I,J) - VCOM(J) 1300 CONTINUE C- CALL VERT1 RETURN END C..... C .. SUBROUTINE GRID3 C..... C .. OUTPUT INFORMATION RELEVANT TO GRID SECTIONING C..... SUBROUTINE GRID3 DIMENSION NGP(3), D1(3),D2(3),D3(3) COMMON /INPT08/ NGP, D1,D2,D3 PARAMETER (NGMAX=32**3) DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3) COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC COMMON /CRYSTL/ TMOLW, RHO, EMUX COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU DIMENSION U1(3),U2(3),U3(3) DATA CONS/57.2957795/ C- CALL UNIVS (D1, U1) CALL UNIVS (D2, U2) CALL UNIVS (D3, U3) CALL VERT1 CALL MINMAX (U1, DMIN1,DMAX1) CALL MINMAX (U2, DMIN2,DMAX2) CALL MINMAX (U3, DMIN3,DMAX3) CALL SDS1 (U1,U2, DOT12) CALL SDS1 (U1,U3, DOT13) CALL SDS1 (U2,U3, DOT23) DOT12 = CONS*ACOS(DOT12) DOT13 = CONS*ACOS(DOT13) DOT23 = CONS*ACOS(DOT23) CMASS = VXTL*RHO C- WRITE(NOU,8001) 8001 FORMAT(/1X,'GRID SECTIONING INDICES AND EXTREMES') WRITE(NOU,8002) 8002 FORMAT(1X,5X,'INDICES',8X,'DMIN',5X,'DMAX',3X,'#PTS') WRITE(NOU,8003) D1,DMIN1,DMAX1,NGP(1) WRITE(NOU,8003) D2,DMIN2,DMAX2,NGP(2) WRITE(NOU,8003) D3,DMIN3,DMAX3,NGP(3) 8003 FORMAT(1X,3F5.0,4X,F6.3,3X,F6.3,I5) WRITE(NOU,8004) DOT12,DOT13,DOT23 8004 FORMAT(/1X,'ANGLE D1^D2: ',F6.2,/, * 1X,'ANGLE D1^D3: ',F6.2,/, * 1X,'ANGLE D2^D3: ',F6.2 ) WRITE(NOU,8005) VXTL 8005 FORMAT(/1X,'CRYSTAL VOLUME(MM**3): ',E12.5) WRITE(NOU,8006) CMASS 8006 FORMAT(1X,'CRYSTAL MASS(MG): ',E12.5) WRITE(NOU,8007) VCOM 8007 FORMAT(1X,'C.O.M. [REAL] VECTOR: ',3E12.4) RETURN END C..... C .. SUBROUTINE INNIE C..... C .. ANSWERS THE AGE-OLD QUESTION: IS THE C .. GENERATED POINT WITHIN THE CRYSTAL? C .. NIN = +1 YES IT IS C .. NIN = -1 NO IT IS NOT C..... SUBROUTINE INNIE (POR, NIN) DIMENSION FS(50,5) COMMON /PLANES/ NF, FS DIMENSION POR(3) DATA EPSIL/0.0001/ C- NIN = +1 DO 100 M=1,NF DOT = POR(1)*FS(M,1) * +POR(2)*FS(M,2) * +POR(3)*FS(M,3) * -EPSIL IF( DOT.GE.FS(M,4) ) NIN = -1 100 CONTINUE RETURN END C..... C .. SUBROUTINE KAPP2 C..... C .. CALCULATE THE U(PHI,TAU,SGM) MATRIX AND CONVERT C .. THE POINTS PABC INTO THEIR [TUBE] COMPONENTS C..... SUBROUTINE KAPP2 CHARACTER*4 WAVE COMMON /INPT04/ WAVE COMMON /LAMBDA/ FLAM DIMENSION ORMAT(3,3) COMMON /INPT09/ ORMAT, NDIFF COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP PARAMETER (NGMAX=32**3) DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3) COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3) COMMON /MATRIX/ UBF, UBI, UBGF, UBGI C- DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3) COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB C- DIMENSION UPHI1(3,3), UPHI2(3,3), UTS(3,3) DIMENSION VKPH(3), VKTB(3), VNTB(3), P(3), D(3) DIMENSION ETAB(3,5) DATA CON/57.2957795/ CHARACTER*4 CUKA,CUKB, MOKA,MOKB, AGKA,AGKB DATA CUKA/'CUKA'/, CUKB/'CUKB'/ DATA MOKA/'MOKA'/, MOKB/'MOKB'/ DATA AGKA/'AGKA'/, AGKB/'AGKB'/ C- C- RETURN IF NO CAPILLARY IF(NKAP.EQ.0) RETURN C- C- DIMENSIONS OF THE CAPPILARY RAD2 = (0.5*DKAP)**2. RPT2 = (0.5*DKAP + TKAP)**2. C- C- LOOK UP THE CAPILLARY LINEAR ABSORPTION COEFFICIENT IF(WAVE.EQ.CUKA .OR. WAVE.EQ.CUKB) NRAD = 1 IF(WAVE.EQ.MOKA .OR. WAVE.EQ.MOKB) NRAD = 2 IF(WAVE.EQ.AGKA .OR. WAVE.EQ.AGKB) NRAD = 3 C- TABLE OF EMUK CALCULATED BY HAND C- ETAB(NRADIATION,NKAPPILARY) IN RECIPROCAL MM C NRAD = 1 CU KA NKAP = 1 "GLASKAPILLAREN" C 2 MO KA 2 "PYREX 7740" C 3 AG KA 3 "AR GLASS" C 4 FUSED SILICA C 5 LINDEMANN GLASS ETAB(1,1) = 10.870 ETAB(1,2) = 7.600 ETAB(1,3) = 13.298 ETAB(1,4) = 8.021 ETAB(1,5) = 1.484 ETAB(2,1) = 1.058 ETAB(2,2) = 0.733 ETAB(2,3) = 1.360 ETAB(2,4) = 0.807 ETAB(2,5) = 0.161 ETAB(3,1) = 0.538 ETAB(3,2) = 0.375 ETAB(3,3) = 0.700 ETAB(3,4) = 0.411 ETAB(3,5) = 0.094 EMUK = ETAB(NRAD,NKAP) C- NICOLET P3/F UPHI MATRIX UPHI1(1,1) = +1. UPHI1(1,2) = +0. UPHI1(1,3) = +0. UPHI1(2,1) = +0. UPHI1(2,2) = +1. UPHI1(2,3) = +0. UPHI1(3,1) = +0. UPHI1(3,2) = +0. UPHI1(3,3) = +1. C- NONIUS CAD-4 AND BUSING-LEVY UPHI MATRIX UPHI2(1,1) = +0. UPHI2(1,2) = -1. UPHI2(1,3) = +0. UPHI2(2,1) = +1. UPHI2(2,2) = +0. UPHI2(2,3) = +0. UPHI2(3,1) = +0. UPHI2(3,2) = +0. UPHI2(3,3) = +1. C- BUILD THE UTS MATRIX EXPLICITLY CTAU = COS (TAU/CON) STAU = SIN (TAU/CON) CSGM = COS (SGM/CON) SSGM = SIN (SGM/CON) UTS(1,1) = +CTAU UTS(1,2) = +STAU*SSGM UTS(1,3) = -STAU*CSGM UTS(2,1) = +0. UTS(2,2) = +CSGM UTS(2,3) = +SSGM UTS(3,1) = +STAU UTS(3,2) = -CTAU*SSGM UTS(3,3) = +CTAU*CSGM C- GET UPTSF AND ITS INVERSE UPTSI IF(NDIFF.EQ.1) CALL MTM1 (UPHI1,UTS, UPTSF) IF(NDIFF.EQ.2) CALL MTM1 (UPHI2,UTS, UPTSF) IF(NDIFF.EQ.3) CALL MTM1 (UPHI2,UTS, UPTSF) CALL AMI3 (UPTSF,UPTSI,UDET) IF(UDET.EQ.0.0) RETURN C- GET UBIG CALL MTM1 (UPTSI,UBGF, UBIG) C- CONVERT KAPILLARY DISPLACEMENT FROM [PHI] TO [TUBE] VKPH(1) = XDSP VKPH(2) = YDSP VKPH(3) = 0.0 CALL MTV1 (UPTSI,VKPH, VKTB) C- CONVERT CENTER-OF-MASS VECTOR FROM [REAL] TO [TUBE] CALL MTV1 (UBIG,VCOM, VNTB) C- FORM THE VECTOR { -N -K } VNK1 = -VNTB(1) -VKTB(1) VNK2 = -VNTB(2) -VKTB(2) VNK3 = -VNTB(3) -VKTB(3) C- CONVERT THE PABC VECTORS FROM [REAL] TO [TUBE] C- AND FORM THE DTUB VECTORS P - N - K DO 1000 N = 1,NPTS P(1) = PABC(N,1) P(2) = PABC(N,2) P(3) = PABC(N,3) CALL MTV1 (UBIG,P, D) DTUB(N,1) = D(1) + VNK1 DTUB(N,2) = D(2) + VNK2 DTUB(N,3) = D(3) + VNK3 1000 CONTINUE C- RETURN END C..... C .. SUBROUTINE KAPP3 C..... C .. CALCULATE AND OUTPUT THE [TUBE] COORDINATES OF THE C .. CRYSTAL VERTICES. CHECK TO SEE THAT ALL THE VERTICES C .. ARE WITHIN THE TUBE RADIUS. C..... SUBROUTINE KAPP3 COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DIMENSION VR(150,3), MM(150,3) COMMON /VERTEX/ NV, VR, MM PARAMETER (NGMAX=32**3) DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3) COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3) COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB C- DIMENSION VKPH(3), VKTB(3), VNTB(3), P(3), D(3) DIMENSION VTUB(150,3) C- RETURN IF NO KAPILLARY IF(NKAP.EQ.0) RETURN C- CONVERT KAPILLARY DISPLACEMENT FROM [PHI] TO [TUBE] VKPH(1) = XDSP VKPH(2) = YDSP VKPH(3) = 0.0 CALL MTV1 (UPTSI,VKPH, VKTB) C- CONVERT CENTER-OF-MASS VECTOR FROM [REAL] TO [TUBE] CALL MTV1 (UBIG,VCOM, VNTB) C- FORM THE VECTOR { -N -K } VNK1 = -VNTB(1) -VKTB(1) VNK2 = -VNTB(2) -VKTB(2) VNK3 = -VNTB(3) -VKTB(3) C- CONVERT THE VERTEX VECTORS FROM [REAL] TO [TUBE] DO 1000 N = 1,NV P(1) = VR(N,1) P(2) = VR(N,2) P(3) = VR(N,3) CALL MTV1 (UBIG,P, D) VTUB(N,1) = D(1) + VNK1 VTUB(N,2) = D(2) + VNK2 VTUB(N,3) = D(3) + VNK3 1000 CONTINUE C- IF(NKAP.EQ.1) WRITE(NOU,9101) IF(NKAP.EQ.2) WRITE(NOU,9102) IF(NKAP.EQ.3) WRITE(NOU,9103) IF(NKAP.EQ.4) WRITE(NOU,9104) IF(NKAP.EQ.5) WRITE(NOU,9105) 9101 FORMAT(/1X,'CRYSTAL IN A "GLASKAPILLAREN" CAPILLARY') 9102 FORMAT(/1X,'CRYSTAL IN A "PYREX 7740" CAPILLARY') 9103 FORMAT(/1X,'CRYSTAL IN A "AR GLASS" CAPILLARY') 9104 FORMAT(/1X,'CRYSTAL IN A FUSED SILICA CAPILLARY') 9105 FORMAT(/1X,'CRYSTAL IN A LINDEMANN GLASS CAPILLARY') WRITE(NOU,9004) EMUK 9004 FORMAT(1X,'CAPILLARY LINEAR ABSORPTION COEFF(MM): ',F6.3) WRITE(NOU,9005) DKAP, TKAP 9005 FORMAT(1X,'CAPILLARY DIAMETER: ',F6.3,/, * 1X,'CAPILLARY THICKNESS:',F6.3 ) WRITE(NOU,9006) SGM, TAU 9006 FORMAT(1X,'TOP ARC(SGM): ',F7.3,/, * 1X,'LOW ARC(TAU): ',F7.3 ) WRITE(NOU,9007) XDSP, YDSP 9007 FORMAT(1X,'X[PHI] DISPLACEMENT(MM): ',F8.5,/, * 1X,'Y[PHI] DISPLACEMENT(MM): ',F8.5 ) C- WRITE OUT THE DISTANCES FROM THE [TUBE] ORIGIN C- TO THE VERTICES IN THE CAPILLARY X,Y PLANE WRITE(NOU,9008) 9008 FORMAT(/1X,'DISTANCES OF VERTICES FROM CAPILLARY WALL IN MM') WRITE(NOU,9009) 9009 FORMAT(1X,'NVERT',4X,'DIST',20X,'[TUBE] COMPONENTS') NWARN = 0 RAD1 = SQRT(RAD2) DO 2000 N=1,NV DIST = SQRT(VTUB(N,1)**2.+VTUB(N,2)**2.) IF(DIST.GE.RAD1) NWARN = 1 DIST = RAD1 - DIST WRITE(NOU,9010) N, DIST, (VTUB(N,I),I=1,3) 9010 FORMAT(1X,I4,2X,F9.6,5X,3E15.6) 2000 CONTINUE C- IF(NWARN.NE.0) WRITE(NOU,9011) 9011 FORMAT(/1X,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!',/, * 1X,'ONE OR MORE VERTICES CALCULATED NOT IN TUBE',/, * 1X,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' ) RETURN END C..... C .. SUBROUTINE LIMS1 C..... C .. CALCULATE AND OUTPUT THE LIMITS OF THE C .. CRYSTAL IN VARIOUS PRINCIPAL DIRECTIONS C..... SUBROUTINE LIMS1 COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU DIMENSION G1(3),U1(3) DIMENSION GG(13,3) C- GG( 1,1) = +1. GG( 1,2) = +0. GG( 1,3) = +0. C- GG( 2,1) = +0. GG( 2,2) = +1. GG( 2,3) = +0. C- GG( 3,1) = +0. GG( 3,2) = +0. GG( 3,3) = +1. C- GG( 4,1) = +1. GG( 4,2) = +1. GG( 4,3) = +0. C- GG( 5,1) = +1. GG( 5,2) = +0. GG( 5,3) = +1. C- GG( 6,1) = +0. GG( 6,2) = +1. GG( 6,3) = +1. C- GG( 7,1) = +1. GG( 7,2) = -1. GG( 7,3) = +0. C- GG( 8,1) = +1. GG( 8,2) = +0. GG( 8,3) = -1. C- GG( 9,1) = +0. GG( 9,2) = +1. GG( 9,3) = -1. C- GG(10,1) = +1. GG(10,2) = +1. GG(10,3) = +1. C- GG(11,1) = +1. GG(11,2) = +1. GG(11,3) = -1. C- GG(12,1) = +1. GG(12,2) = -1. GG(12,3) = +1. C- GG(13,1) = +1. GG(13,2) = -1. GG(13,3) = -1. C- WRITE(NOU,9001) 9001 FORMAT(/1X,'LIMITS OF CRYSTAL IN MILLIMETERS') WRITE(NOU,9002) 9002 FORMAT(1X,2X,'DIRECTION',12X,'DMIN',9X,'DMAX',8X,'DELT') DO 100 I=1,13 G1(1) = GG(I,1) G1(2) = GG(I,2) G1(3) = GG(I,3) CALL UNIVS (G1, U1) CALL MINMAX (U1, DMIN, DMAX) DDEL = ABS(DMAX-DMIN) WRITE(NOU,9003) G1,DMIN,DMAX,DDEL 9003 FORMAT(1X,3F4.0,5X,3F12.6) 100 CONTINUE RETURN END C..... C .. SUBROUTINE MATX2 C..... C .. CALCULATE THE UBF UBI UBGF AND UBGI MATRICES C..... SUBROUTINE MATX2 DIMENSION ORMAT(3,3) COMMON /INPT09/ ORMAT, NDIFF DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3) COMMON /MATRIX/ UBF, UBI, UBGF, UBGI C- DO 1200 I=1,3 DO 1100 J=1,3 UBF(I,J) = ORMAT(I,J) 1100 CONTINUE 1200 CONTINUE CALL AMI3 (UBF, UBI,UDET) CALL MTM1 (UBF,GR, UBGF) CALL AMI3 (UBGF, UBGI,UDET) RETURN END C..... C .. SUBROUTINE MATX3 C .. OUTPUT THE DATA ON THE ORIENTATION MATRIX AND DIFFRACTOMETER TYPE C..... SUBROUTINE MATX3 DIMENSION ORMAT(3,3) COMMON /INPT09/ ORMAT, NDIFF DIMENSION UBF(3,3), UBI(3,3), UBGF(3,3), UBGI(3,3) COMMON /MATRIX/ UBF, UBI, UBGF, UBGI COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C- IF(NDIFF.EQ.1) GO TO 100 IF(NDIFF.EQ.2) GO TO 200 IF(NDIFF.EQ.3) GO TO 300 100 WRITE(NOU,8001) GO TO 500 200 WRITE(NOU,8002) GO TO 500 300 WRITE(NOU,8003) GO TO 500 8001 FORMAT(/1X,'NICOLET P3/F DIFFRACTOMETER') 8002 FORMAT(/1X,'NONIUS CAD-4 DIFFRACTOMETER') 8003 FORMAT(/1X,'BUSING-LEVY DIFFRACTOMETER') 500 CONTINUE C- WRITE(NOU,9001) 9001 FORMAT(1X,'ORIENTATION MATRIX UBF(I,J)') DO 1000 I=1,3 WRITE(NOU,9002) (UBF(I,J),J=1,3) 1000 CONTINUE 9002 FORMAT(1X,3E15.6) RETURN END C..... C .. SUBROUTINE MINMAX C..... C .. GIVEN A LIST OF COMPONENTS OF VERTEX VECTORS C .. IN THE [REAL] FRAME SEARCH FOR THE MINIMUM C .. AND MAXIMUM OF THE DOT-PRODUCT U*-DOT-V C .. GIVEN UNIT VECTOR U* IN THE [STAR] FRAME C..... SUBROUTINE MINMAX(UVECT, DMIN,DMAX) DIMENSION VR(150,3), MM(150,3) COMMON /VERTEX/ NV, VR, MM DIMENSION UVECT(3) C- U1 = UVECT(1) U2 = UVECT(2) U3 = UVECT(3) DMIN = +1E10 DMAX = -1E10 DO 10 J=1,NV V1 = VR(J,1) V2 = VR(J,2) V3 = VR(J,3) DOT = U1*V1 + U2*V2 + U3*V3 IF (DOT.LE.DMIN) DMIN = DOT IF (DOT.GE.DMAX) DMAX = DOT 10 CONTINUE RETURN END C..... C .. SUBROUTINE MTM1 C..... C .. MULTIPLY MATRIX B BY MATRIX A C .. TO GET MATRIX C C = A B C..... SUBROUTINE MTM1 (A,B, C) DIMENSION A(3,3), B(3,3), C(3,3) DO 10 I=1,3 DO 5 J=1,3 C(I,J) = 0. 5 CONTINUE 10 CONTINUE C- DO 30 I=1,3 DO 25 J=1,3 DO 20 N=1,3 C(I,J) = C(I,J) + A(I,N)*B(N,J) 20 CONTINUE 25 CONTINUE 30 CONTINUE C- RETURN END C..... C .. SUBROUTINE MTV1 C..... C .. MULTIPLY VECTOR -R- BY MATRIX -G- TO GET C .. VECTOR -S- C..... SUBROUTINE MTV1(G, R, S) DIMENSION G(3,3), R(3), S(3) S(1) = G(1,1)*R(1) + G(1,2)*R(2) + G(1,3)*R(3) S(2) = G(2,1)*R(1) + G(2,2)*R(2) + G(2,3)*R(3) S(3) = G(3,1)*R(1) + G(3,2)*R(2) + G(3,3)*R(3) RETURN END C..... C .. SUBROUTINE OUTP1 C..... C .. PRINT OUT A SUMMARY OF INPUT AND PRIMARY OUTPUT DATA C..... SUBROUTINE OUTP1 COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C- OPEN(UNIT=NLP,FILE='absorb.lp',STATUS='UNKNOWN') DO I=1,1E9 READ (NLP,'(A1)',END=1) DUMMY END DO 1 WRITE (NLP,'(/1X)') I=NOU NOU=NLP CALL TITL2 CALL FILS2 CALL BPRO2 CALL CELL2 CALL CELL3 CALL FACE2 CALL FACE3 CALL VERT1 CALL VERT5 CALL LIMS1 CALL EDGE1 CALL MATX2 CALL MATX3 CALL EMUX2 CALL EMUX3 CALL GRID2 CALL GRID3 CALL KAPP2 CALL KAPP3 CALL SOLN4 CLOSE(UNIT=NLP,STATUS='KEEP') NOU=I RETURN END C..... C .. SUBROUTINE P3RMT C..... C .. GET THE BUSING-LEVY ROTATION MATRIX FOR THE C .. NICOLET P3/F DIFFRACTOMETER C .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI C..... SUBROUTINE P3RMT(ANGLS, RMTF,ST,CT) DIMENSION ANGLS(4), RMTF(3,3) DATA CONS/57.29577951/ C- TTH = ANGLS(1)/CONS OME = ANGLS(2)/CONS PHI = ANGLS(3)/CONS CHI = ANGLS(4)/CONS THT = 0.5*TTH OMT = OME - THT ST = SIN(THT) CT = COS(THT) SO = SIN(OMT) CO = COS(OMT) SP = SIN(PHI) CP = COS(PHI) SC = SIN(CHI) CC = COS(CHI) C- EXPLICIT GENERATION OF ROTATION MATRIX RMTF(1,1) = +CO*CP + SO*CC*SP RMTF(1,2) = +CO*SP - SO*CC*CP RMTF(1,3) = -SO*SC RMTF(2,1) = +SO*CP - CO*CC*SP RMTF(2,2) = +SO*SP + CO*CC*CP RMTF(2,3) = +CO*SC RMTF(3,1) = +SC*SP RMTF(3,2) = -SC*CP RMTF(3,3) = +CC C- RETURN END C..... C .. SUBROUTINE P3ROT C..... C .. GET THE [PHI] FRAME COMPONENTS OF THE UNIT C .. VECTORS FOR THE NICOLET P3/F DIFFRACTOMETER C .. -SI- REVERSE INCIDENT BEAM C .. -SD- DIFFRACTED BEAM C .. -PM- MONOCHROMATOR HORIZONTAL PROFILE C .. ANGLES ARE: TWO THETA, OMEGA, PHI, CHI C..... SUBROUTINE P3ROT(ANGLS, SIP,SDP,PMP) DIMENSION ANGLS(4), SIP(3), SDP(3), PMP(3) DIMENSION RMTF(3,3), RMTI(3,3), SDT(3), SIT(3), PMT(3) C- CALL P3RMT (ANGLS, RMTF,ST,CT) C- C- TRANSPOSE OF ROTATION MATRIX EQUALS ITS INVERSE DO 100 I=1,3 DO 50 J=1,3 RMTI(I,J) = RMTF(J,I) 50 CONTINUE 100 CONTINUE C- SIT(1) = -CT SIT(2) = +ST SIT(3) = 0.0 SDT(1) = +CT SDT(2) = +ST SDT(3) = 0.0 PMT(1) = +CT PMT(2) = +ST PMT(3) = +0.0 C- CALL MTV1(RMTI,SIT, SIP) CALL MTV1(RMTI,SDT, SDP) CALL MTV1(RMTI,PMT, PMP) RETURN END C..... C .. SUBROUTINE PAGE1 C..... C .. READ IN FROM FILE NAB THE ARCHIVAL QUALITY DATA C..... SUBROUTINE PAGE1 CHARACTER TITLE*80, ATIME*8, ADATE*9 CHARACTER FIL01*80, FIL02*80 DIMENSION CELL(6) CHARACTER*4 WAVE CHARACTER*2 ATOMS(10) DIMENSION STOIC(10) DIMENSION FACES(50,4) DIMENSION NGP(3), D1(3),D2(3),D3(3) DIMENSION ORMAT(3,3) DIMENSION NTRP(50) COMMON /INPT01/ TITLE, ATIME, ADATE COMMON /INPT02/ FIL01, FIL02, IPATH COMMON /INPT03/ CELL COMMON /INPT04/ WAVE COMMON /LAMBDA/ FLAM COMMON /INPT05/ ATOMS COMMON /INPT06/ ZCELL, STOIC COMMON /INPT07/ FACES COMMON /INPT08/ NGP, D1,D2,D3 COMMON /INPT09/ ORMAT, NDIFF COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP COMMON /INPT11/ NBPR, A0, A1, A2, A3 COMMON /INPT12/ NSOL, EMUS, NTRP COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU COMMON /CRYSTL/ TMOLW, RHO, EMUX C- OPEN(UNIT=NAB,FILE='absorb.dat',STATUS='OLD',ERR=999) C- READ(NAB,7001) TITLE READ(NAB,7001) ATIME READ(NAB,7001) ADATE READ(NAB,7001) FIL01 READ(NAB,7001) FIL02 READ(NAB,7010) IPATH READ(NAB,7002) CELL READ(NAB,7001) WAVE READ(NAB,7002) FLAM READ(NAB,7003) ATOMS READ(NAB,7004) ZCELL READ(NAB,7005) STOIC READ(NAB,7002) EMUX READ(NAB,7006) ( (FACES(I,J),J=1,4), I=1,50) READ(NAB,7007) NGP READ(NAB,7008) D1,D2,D3 READ(NAB,7009) ( (ORMAT(I,J),J=1,3), I=1,3) READ(NAB,7010) NDIFF READ(NAB,7011) NKAP, DKAP, TKAP READ(NAB,7012) TAU, SGM, XDSP, YDSP READ(NAB,7013) NBPR, A0, A1, A2, A3 READ(NAB,7014) NSOL, EMUS READ(NAB,7015) NTRP CLOSE(UNIT=NAB,STATUS='KEEP') RETURN C- 999 WRITE(IO6,8001) RETURN C- 7001 FORMAT(A) 7002 FORMAT(6F10.5) 7003 FORMAT(10A) 7004 FORMAT(F5.2) 7005 FORMAT(10F7.2) 7006 FORMAT(3F10.3,F10.5,/) 7007 FORMAT(3I5) 7008 FORMAT(3F6.2,/) 7009 FORMAT(3E15.7) 7010 FORMAT(I5) 7011 FORMAT(I5,F7.4,F7.4) 7012 FORMAT(4F10.5) 7013 FORMAT(I5,4F10.6) 7014 FORMAT(I5,F10.6) 7015 FORMAT(50I1) C- 8001 FORMAT(1X,'ERROR OPENING THE ABSORB.DAT FILE') END C..... C .. SUBROUTINE PAGE2 C..... C .. WRITE OUT ON FILE NAB THE ARCHIVAL QUALITY DATA C..... SUBROUTINE PAGE2 CHARACTER TITLE*80, ATIME*8, ADATE*9 CHARACTER FIL01*80, FIL02*80 DIMENSION CELL(6) CHARACTER*4 WAVE CHARACTER*2 ATOMS(10) DIMENSION STOIC(10) DIMENSION FACES(50,4) DIMENSION NGP(3), D1(3),D2(3),D3(3) DIMENSION ORMAT(3,3) DIMENSION NTRP(50) COMMON /INPT01/ TITLE, ATIME, ADATE COMMON /INPT02/ FIL01, FIL02, IPATH COMMON /INPT03/ CELL COMMON /INPT04/ WAVE COMMON /LAMBDA/ FLAM COMMON /INPT05/ ATOMS COMMON /INPT06/ ZCELL, STOIC COMMON /INPT07/ FACES COMMON /INPT08/ NGP, D1,D2,D3 COMMON /INPT09/ ORMAT, NDIFF COMMON /INPT10/ NKAP, DKAP, TKAP, TAU, SGM, XDSP, YDSP COMMON /INPT11/ NBPR, A0, A1, A2, A3 COMMON /INPT12/ NSOL, EMUS, NTRP COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU COMMON /CRYSTL/ TMOLW, RHO, EMUX OPEN(UNIT=NAB,FILE='absorb.dat.new',STATUS='NEW') WRITE(NAB,7001) TITLE WRITE(NAB,7001) ATIME WRITE(NAB,7001) ADATE WRITE(NAB,7001) FIL01 WRITE(NAB,7001) FIL02 WRITE(NAB,7010) IPATH WRITE(NAB,7002) CELL WRITE(NAB,7001) WAVE WRITE(NAB,7002) FLAM WRITE(NAB,7003) ATOMS WRITE(NAB,7004) ZCELL WRITE(NAB,7005) STOIC WRITE(NAB,7002) EMUX WRITE(NAB,7006) ( (FACES(I,J),J=1,4), I=1,50) WRITE(NAB,7007) NGP WRITE(NAB,7008) D1,D2,D3 WRITE(NAB,7009) ( (ORMAT(I,J),J=1,3), I=1,3) WRITE(NAB,7010) NDIFF WRITE(NAB,7011) NKAP, DKAP, TKAP WRITE(NAB,7012) TAU, SGM, XDSP, YDSP WRITE(NAB,7013) NBPR, A0, A1, A2, A3 WRITE(NAB,7014) NSOL, EMUS WRITE(NAB,7015) NTRP CLOSE(UNIT=NAB,STATUS='KEEP') 7001 FORMAT(A) 7002 FORMAT(6F10.5) 7003 FORMAT(10A) 7004 FORMAT(F5.2) 7005 FORMAT(10F7.2) 7006 FORMAT(3F10.3,F10.5,/) 7007 FORMAT(3I5) 7008 FORMAT(3F6.2,/) 7009 FORMAT(3E15.7) 7010 FORMAT(I5) 7011 FORMAT(I5,F7.4,F7.4) 7012 FORMAT(4F10.5) 7013 FORMAT(I5,4F10.6) 7014 FORMAT(I5,F10.6) 7015 FORMAT(50I1) RETURN END C..... C .. SUBROUTINE PATH1 C..... C .. FIND THE ENTRY/EXIT FACE OF THE BEAM C .. WITH [REAL] SPACE COMPONENTS -SR- FOR C .. THE POINT WITH [REAL] COMPONENTS -P- C .. AND CALCULATE THE PATH LENGTH -PL- C .. IN MILLIMETERS C..... SUBROUTINE PATH1 (SR,JPT, PL,NX) DIMENSION FS(50,5) COMMON /PLANES/ NF, FS PARAMETER (NGMAX=32**3) DIMENSION VCOM(3), PABC(NGMAX,3), WABC(NGMAX), RABC(NGMAX,3) COMMON /FGAUSS/ VXTL, VCOM, NPTS, PABC, WABC, RABC DIMENSION SR(3) C- PL = 1E12 NX = 0 S1 = SR(1) S2 = SR(2) S3 = SR(3) P1 = PABC(JPT,1) P2 = PABC(JPT,2) P3 = PABC(JPT,3) C- DO 1000 N=1,NF PN = 0. U1 = FS(N,1) U2 = FS(N,2) U3 = FS(N,3) DV = FS(N,4) SU = U1*S1 + U2*S2 + U3*S3 IF( SU.LE.0.0 ) GO TO 999 PU = U1*P1 + U2*P2 + U3*P3 PN = ( DV - PU ) / SU IF( PN.GE.PL ) GO TO 999 PL = PN NX = N 999 FS(N,5) = PN 1000 CONTINUE C- RETURN END C..... C .. SUBROUTINE PATH2 C..... C .. CALCULATE THE PATH LENGTH THROUGH THE CAPILLARY C .. WINN IS DISTANCE IN MM TO INNER WALL C .. WOUT IS DISTANCE IN MM TO OUTER WALL C .. PCAP IS DISTANCE IN MM THROUGH CAPILLARY C..... SUBROUTINE PATH2 (S,JPT, WINN,WOUT,PCAP) PARAMETER (NGMAX=32**3) DIMENSION UPTSF(3,3), UPTSI(3,3), UBIG(3,3), DTUB(NGMAX,3) COMMON /KAPCOR/ EMUK, RAD2, RPT2, UPTSF, UPTSI, UBIG, DTUB DIMENSION S(3) S1 = S(1) S2 = S(2) D1 = DTUB(JPT,1) D2 = DTUB(JPT,2) PT1 = D1*S1 + D2*S2 PT2 = S1*S1 + S2*S2 PT3 = D1*D1 + D2*D2 WINN = (-PT1 + SQRT(PT1*PT1 - PT2*(PT3-RAD2))) / PT2 WOUT = (-PT1 + SQRT(PT1*PT1 - PT2*(PT3-RPT2))) / PT2 PCAP = WOUT - WINN RETURN END C..... C .. SUBROUTINE PATH3 C..... C .. CALCULATE PATH LENGTH THROUGH SOLVENT C..... SUBROUTINE PATH3 (WINN, PSOL) DIMENSION NTRP(50) COMMON /INPT12/ NSOL, EMUS, NTRP DIMENSION FS(50,5) COMMON /PLANES/ NF, FS C- PSOL = 0.0 DO 1000 N=1,NF C- DOES PLANE N TRAP SOLVENT? IF(NTRP(N).EQ.0) GO TO 1000 C- IS RAY ACUTE TO FACE NORMAL? F = FS(N,5) IF(F.LE.0.) GO TO 1000 C- DOES RAY INTERSECT FACE WITHIN CAPILLARY? IF(F.GE.WINN) GO TO 1000 C- RAY GOES THROUGH SOLVENT PF = WINN - F C- MAXIMUM DISTANCE IS CORRECT IF( PF .GE. PSOL ) PSOL = PF 1000 CONTINUE RETURN END C..... C .. SUBROUTINE RTOS1 C..... C .. CONVERT VECTOR COMPONENTS IN THE [REAL] FRAME C .. TO COMPONENTS IN THE [STAR] FRAME. C..... SUBROUTINE RTOS1(R,S) DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DIMENSION R(3), S(3) S(1) = GR(1,1)*R(1) + GR(1,2)*R(2) + GR(1,3)*R(3) S(2) = GR(2,1)*R(1) + GR(2,2)*R(2) + GR(2,3)*R(3) S(3) = GR(3,1)*R(1) + GR(3,2)*R(2) + GR(3,3)*R(3) RETURN END C..... C .. SUBROUTINE SDS1 C..... C .. GIVEN TWO VECTORS -S1- AND -S2- GET THEIR DOT PRODUCT C..... SUBROUTINE SDS1 (S1, S2, DOT) DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DIMENSION S1(3), S2(3) C- X1 = GS(1,1)*S1(1) + GS(1,2)*S1(2) + GS(1,3)*S1(3) X2 = GS(2,1)*S1(1) + GS(2,2)*S1(2) + GS(2,3)*S1(3) X3 = GS(3,1)*S1(1) + GS(3,2)*S1(2) + GS(3,3)*S1(3) C- DOT = S2(1)*X1 + S2(2)*X2 + S2(3)*X3 C- RETURN END C..... C .. SUBROUTINE SOLN4 C..... C .. OUTPUT THE SOLVENT INFORMATION C..... SUBROUTINE SOLN4 DIMENSION NTRP(50) COMMON /INPT12/ NSOL, EMUS, NTRP COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU DIMENSION FS(50,5) COMMON /PLANES/ NF, FS C- IF (NSOL .NE. 0) GO TO 1001 WRITE(NOU,9001) 9001 FORMAT(/1X,'NO TRAPPED SOLVENT') RETURN 1001 WRITE(NOU,9002) 9002 FORMAT(/1X,'SOLVENT TRAPPED BY ONE OR MORE FACES') WRITE(NOU,9003) EMUS 9003 FORMAT(1X,'MU OF SOLVENT(1/MM): ',F7.3) WRITE(NOU,9004) 9004 FORMAT(1X,'STATUS: NOT TRAPPING(0); TRAPPING(1)') WRITE(NOU,9005) 9005 FORMAT(1X,'FACE#',3X,'HKL',12X,'STATUS') DO 1002 J=1,NF WRITE(NOU,9006) J, (FS(J,I),I=1,3), NTRP(J) 9006 FORMAT(1X,I5,3F5.0,5X,I1) 1002 CONTINUE RETURN END C..... C .. SUBROUTINE TITL2 C..... C .. OUTPUT THE TITLE DATE AND PROGRAM HEADER C..... SUBROUTINE TITL2 CHARACTER TITLE*80, ATIME*8, ADATE*9 COMMON /INPT01/ TITLE, ATIME, ADATE COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU CALL TIME(ATIME) CALL DATE(ADATE) WRITE(NOU,9001) ATIME, ADATE 9001 FORMAT(/1X,'PROGRAM ABSORB1...VERSION 1',5X,A,3X,A) WRITE(NOU,9002) TITLE 9002 FORMAT(1X,A) RETURN END C..... C .. SUBROUTINE TWOCO C..... C .. "OH. I GET IT. NOBODY CARES!" C .. THE FOONEWAL DEWECTOR, 1982 C .. FIND OUT IF TWO SETS OF THREE NUMBERS C .. HAVE ANY TWO IN COMMON. A ROUTINE TO C .. DISTINGUISH TRUE EDGES FROM GENERAL C .. DIAGONALS. C..... C .. NTWO = 0 NO MATCH C .. NTWO = 1 A MATCH C..... SUBROUTINE TWOCO (I1,I2,I3, J1,J2,J3, NTWO) NTWO = 0 C- IF( I1.EQ.J1 .AND. I2.EQ.J2 ) NTWO = 1 IF( I1.EQ.J1 .AND. I2.EQ.J3 ) NTWO = 1 IF( I1.EQ.J1 .AND. I3.EQ.J2 ) NTWO = 1 IF( I1.EQ.J1 .AND. I3.EQ.J3 ) NTWO = 1 C- IF( I1.EQ.J2 .AND. I2.EQ.J1 ) NTWO = 1 IF( I1.EQ.J2 .AND. I2.EQ.J3 ) NTWO = 1 IF( I1.EQ.J2 .AND. I3.EQ.J1 ) NTWO = 1 IF( I1.EQ.J2 .AND. I3.EQ.J3 ) NTWO = 1 C- IF( I1.EQ.J3 .AND. I2.EQ.J1 ) NTWO = 1 IF( I1.EQ.J3 .AND. I2.EQ.J2 ) NTWO = 1 IF( I1.EQ.J3 .AND. I3.EQ.J1 ) NTWO = 1 IF( I1.EQ.J3 .AND. I3.EQ.J2 ) NTWO = 1 C- IF( I2.EQ.J1 .AND. I1.EQ.J2 ) NTWO = 1 IF( I2.EQ.J1 .AND. I1.EQ.J3 ) NTWO = 1 IF( I2.EQ.J1 .AND. I3.EQ.J2 ) NTWO = 1 IF( I2.EQ.J1 .AND. I3.EQ.J3 ) NTWO = 1 C- IF( I2.EQ.J2 .AND. I1.EQ.J1 ) NTWO = 1 IF( I2.EQ.J2 .AND. I1.EQ.J3 ) NTWO = 1 IF( I2.EQ.J2 .AND. I3.EQ.J1 ) NTWO = 1 IF( I2.EQ.J2 .AND. I3.EQ.J3 ) NTWO = 1 C- IF( I2.EQ.J3 .AND. I1.EQ.J1 ) NTWO = 1 IF( I2.EQ.J3 .AND. I1.EQ.J2 ) NTWO = 1 IF( I2.EQ.J3 .AND. I3.EQ.J1 ) NTWO = 1 IF( I2.EQ.J3 .AND. I3.EQ.J2 ) NTWO = 1 C- IF( I3.EQ.J1 .AND. I1.EQ.J2 ) NTWO = 1 IF( I3.EQ.J1 .AND. I1.EQ.J3 ) NTWO = 1 IF( I3.EQ.J1 .AND. I2.EQ.J2 ) NTWO = 1 IF( I3.EQ.J1 .AND. I2.EQ.J3 ) NTWO = 1 C- IF( I3.EQ.J2 .AND. I1.EQ.J1 ) NTWO = 1 IF( I3.EQ.J2 .AND. I1.EQ.J3 ) NTWO = 1 IF( I3.EQ.J2 .AND. I2.EQ.J1 ) NTWO = 1 IF( I3.EQ.J2 .AND. I2.EQ.J3 ) NTWO = 1 C- IF( I3.EQ.J3 .AND. I1.EQ.J1 ) NTWO = 1 IF( I3.EQ.J3 .AND. I1.EQ.J2 ) NTWO = 1 IF( I3.EQ.J3 .AND. I2.EQ.J1 ) NTWO = 1 IF( I3.EQ.J3 .AND. I2.EQ.J2 ) NTWO = 1 C- RETURN END C..... C .. SUBROUTINE UNIVR C..... C .. GIVEN THE COMPONENTS OF A REAL VECTOR DVEC C .. GET THE COMPONENTS OF THE UNIT VECTOR C .. UVEC PARALLEL TO DVEC C..... SUBROUTINE UNIVR (DVEC, UVEC) DIMENSION DVEC(3), UVEC(3) C- UVEC(1) = 0. UVEC(2) = 0. UVEC(3) = 0. D1 = DVEC(1) D2 = DVEC(2) D3 = DVEC(3) CALL LREAL (D1,D2,D3, DREL) IF (DREL.LE.0.0) RETURN UVEC(1) = D1/DREL UVEC(2) = D2/DREL UVEC(3) = D3/DREL RETURN END C..... C .. SUBROUTINE UNIVS C..... C .. GIVEN THE COMPONENTS OF A RECIPROCAL VECTOR C .. DVEC, GET THE COMPONENTS OF THE UNIT VECTOR C .. UVEC PARALLEL TO DVEC C..... SUBROUTINE UNIVS (DVEC, UVEC) DIMENSION DVEC(3), UVEC(3) C- UVEC(1) = 0. UVEC(2) = 0. UVEC(3) = 0. D1 = DVEC(1) D2 = DVEC(2) D3 = DVEC(3) CALL LSTAR (D1,D2,D3, DSTR) IF (DSTR.LE.0.0) RETURN UVEC(1) = D1/DSTR UVEC(2) = D2/DSTR UVEC(3) = D3/DSTR RETURN END C..... C .. SUBROUTINE VERT1 C..... C .. FROM A LIST OF UNIT VECTORS PERPENDICULAR TO THE C .. CRYSTAL FACES AND THEIR PERPENDICULAR DISTANCES C .. GENERATE THE [REAL] COMPONENTS OF THE VERTEX C .. VECTORS. C..... SUBROUTINE VERT1 DIMENSION FS(50,5) COMMON /PLANES/ NF, FS DIMENSION VR(150,3), MM(150,3) COMMON /VERTEX/ NV, VR, MM DIMENSION FF(3,3),FI(3,3), DT(3),POR(3) C- NV = 0 IF (NF .LE. 3) RETURN DO 30 I=1,NF-2 FF(1,1) = FS(I,1) FF(1,2) = FS(I,2) FF(1,3) = FS(I,3) DT(1) = FS(I,4) C- DO 20 J=I+1,NF-1 FF(2,1) = FS(J,1) FF(2,2) = FS(J,2) FF(2,3) = FS(J,3) DT(2) = FS(J,4) C- DO 10 K=J+1,NF FF(3,1) = FS(K,1) FF(3,2) = FS(K,2) FF(3,3) = FS(K,3) DT(3) = FS(K,4) C- CALL AMI3 (FF, FI, DET) IF (DET.EQ.0.) GO TO 10 CALL MTV1 (FI, DT, POR) CALL INNIE (POR, NIN) IF (NIN.EQ.-1) GO TO 10 C- NV = NV + 1 VR(NV,1) = POR(1) VR(NV,2) = POR(2) VR(NV,3) = POR(3) C- MM(NV,1) = I MM(NV,2) = J MM(NV,3) = K C- 10 CONTINUE 20 CONTINUE 30 CONTINUE C- RETURN END C..... C .. SUBROUTINE VERT2 C..... C .. GIVEN A PARTICULAR CRYSTAL PLANE -F1- GENERATE C .. THE COMPONENTS IN THE [REAL] FRAMES C .. OF ALL THE VERTICES FORMED BY -F1-. C..... SUBROUTINE VERT2(F1) DIMENSION FS(50,5) COMMON /PLANES/ NF, FS DIMENSION VR(150,3), MM(150,3) COMMON /VERTEX/ NV, VR, MM DIMENSION FF(3,3),FI(3,3), DT(3),POR(3) DIMENSION F1(4) C- NV = 0 IF (NF .LE. 2) RETURN FF(1,1) = F1(1) FF(1,2) = F1(2) FF(1,3) = F1(3) DT(1) = F1(4) C- DO 20 J=1,NF-1 FF(2,1) = FS(J,1) FF(2,2) = FS(J,2) FF(2,3) = FS(J,3) DT(2) = FS(J,4) C- DO 10 K=J+1,NF FF(3,1) = FS(K,1) FF(3,2) = FS(K,2) FF(3,3) = FS(K,3) DT(3) = FS(K,4) C- CALL AMI3 (FF, FI, DET) IF (DET.EQ.0.) GO TO 10 CALL MTV1 (FI, DT, POR) CALL INNIE (POR, NIN) IF (NIN.EQ.-1) GO TO 10 C- NV = NV + 1 VR(NV,1) = POR(1) VR(NV,2) = POR(2) VR(NV,3) = POR(3) C- 10 CONTINUE 20 CONTINUE C- RETURN END C..... C .. SUBROUTINE VERT3 C..... C .. GIVEN TWO PARTICULAR CRYSTAL PLANES -F1- AND -F2- C .. GENERATE THE [REAL] COMPONENTS OF ALL C .. THE VERTICES FORMED WITH -F1- AND -F2-. C..... SUBROUTINE VERT3(F1,F2) DIMENSION FS(50,5) COMMON /PLANES/ NF, FS DIMENSION VR(150,3), MM(150,3) COMMON /VERTEX/ NV, VR, MM DIMENSION FF(3,3),FI(3,3), DT(3),POR(3) DIMENSION F1(4), F2(4) C- NV = 0 IF (NF .LE. 1) RETURN FF(1,1) = F1(1) FF(1,2) = F1(2) FF(1,3) = F1(3) DT(1) = F1(4) C- FF(2,1) = F2(1) FF(2,2) = F2(2) FF(2,3) = F2(3) DT(2) = F2(4) C- DO 10 K=1,NF FF(3,1) = FS(K,1) FF(3,2) = FS(K,2) FF(3,3) = FS(K,3) DT(3) = FS(K,4) C- CALL AMI3 (FF, FI, DET) IF (DET.EQ.0.) GO TO 10 CALL MTV1 (FI, DT, POR) CALL INNIE (POR, NIN) IF (NIN.EQ.-1) GO TO 10 C- NV = NV + 1 VR(NV,1) = POR(1) VR(NV,2) = POR(2) VR(NV,3) = POR(3) C- 10 CONTINUE C- RETURN END C..... C .. SUBROUTINE VERT4 C..... C .. GIVEN THREE PARTICULAR FACE VECTORS -F1- -F2- -F3- C .. WITH COMPONENTS IN THE [STAR] FRAME GET THEIR POINT C .. OF INTERSECTION (VERTEX) IN THE [REAL] FRAMES. C .. DO NOT CHECK IF POINT IS IN CRYSTAL SINCE C .. IT MUST BE BY THE METHOD OF GENERATION. C..... SUBROUTINE VERT4(F1,F2,F3, POR) DIMENSION FF(3,3),FI(3,3), DT(3),POR(3) DIMENSION F1(4), F2(4), F3(4) C- FF(1,1) = F1(1) FF(1,2) = F1(2) FF(1,3) = F1(3) DT(1) = F1(4) C- FF(2,1) = F2(1) FF(2,2) = F2(2) FF(2,3) = F2(3) DT(2) = F2(4) C- FF(3,1) = F3(1) FF(3,2) = F3(2) FF(3,3) = F3(3) DT(3) = F3(4) C- CALL AMI3 (FF, FI, DET) IF (DET.EQ.0.) RETURN CALL MTV1 (FI, DT, POR) C- RETURN END C..... C .. SUBROUTINE VERT5 C..... C .. GIVEN A SET OF VERTEX VECTORS AND THEIR CONNECTION C .. TABLE, PRINT OR TYPE OUT THE LIST OF VECTORS. C..... SUBROUTINE VERT5 DIMENSION VR(150,3), MM(150,3) COMMON /VERTEX/ NV, VR, MM COMMON /LOCAL1/ IO5, IO6, NLP, NAB, NIN, NOU C- WRITE(NOU,9001) 9001 FORMAT(/1X,'VERTEX VECTORS') WRITE(NOU,9002) 9002 FORMAT(1X,' #',' FACES ',3X, * ' [REAL] COMPONENTS OF VERTICES') DO 10 I=1,NV WRITE(NOU,9003) I,(MM(I,J),J=1,3),(VR(I,J),J=1,3) 9003 FORMAT(1X,I3,3I4,3X,3E13.4) 10 CONTINUE RETURN END C..... C .. SUBROUTINE VSTAR C..... C .. GET THE SCALE FACTOR WHICH YIELDS THE VOLUME C .. OF THE GRID POINT WHEN MULTIPLIED BY THE THREE C .. GAUSSIAN FRACTIONAL LENGTHS EPS1, EPS2, EPS3 C .. SCALE = D*1 X D*2 X D*3 / DET X V* C..... SUBROUTINE VSTAR (VSTR) DIMENSION NGP(3), D1(3),D2(3),D3(3) COMMON /INPT08/ NGP, D1,D2,D3 DIMENSION RR(6), RS(6), GR(3,3), GS(3,3), BM(3,3) COMMON /UNITCL/ RR, RS, GR, GS, CVR, CVS, BM DIMENSION FHKL(3,3), GHKL(3,3) C- FH = D1(1) FK = D1(2) FL = D1(3) CALL LSTAR (FH,FK,FL, D1STR) FH = D2(1) FK = D2(2) FL = D2(3) CALL LSTAR (FH,FK,FL, D2STR) FH = D3(1) FK = D3(2) FL = D3(3) CALL LSTAR (FH,FK,FL, D3STR) FHKL(1,1) = D1(1) FHKL(1,2) = D1(2) FHKL(1,3) = D1(3) FHKL(2,1) = D2(1) FHKL(2,2) = D2(2) FHKL(2,3) = D2(3) FHKL(3,1) = D3(1) FHKL(3,2) = D3(2) FHKL(3,3) = D3(3) CALL AMI3 (FHKL, GHKL,VDET) VSTR = (D1STR*D2STR*D3STR) / (VDET*CVS) RETURN END