C A MINI INTERACTIVE INTRAB DIMENSION RLC(3),X(3),Y(3),Z(3),XA(3,3),D(3,3) CHARACTER*1 ANS DPR=57.29578 OPEN (7,FILE='MINITRAB.DAT',STATUS='OLD') READ(7,*)RLC,CA,CB,CG WRITE(*,*)'THESE ARE THE CELL PARAMETERS ON MINITRAB.DAT' WRITE(*,21)RLC,CA,CB,CG WRITE(*, '(2X,''ARE THESE CORRECT [Y/N]?'')') READ(*,10)ANS 10 FORMAT(A1) IF(ANS.EQ.'Y')GO TO 101 50 WRITE(*, '(2X,''PROVIDE THEM IN FREE FORMAT'')') READ(5,*,ERR=50)RLC,CA,CB,CG BACKSPACE 7 WRITE(7,21)RLC,CA,CB,CG 21 FORMAT(3F9.3,3F9.2) 101 CA=COS(CA/DPR) CB=COS(CB/DPR) CG=COS(CG/DPR) SG=SQRT(1.-CG*CG) CR=(CA-CB*CG)/SG CP=SQRT(1.-CB*CB-CR*CR) 100 NUM=2 WRITE(*, '(2X,''FOR ANGLE AT B GIVE ATOMS IN ORDER A-B-C'')') 150 WRITE(*, '(2X,''GIVE FRACT COORDS OF FIRST ATOM'')') READ(5,*,ERR=150)X(1),Y(1),Z(1) 160 WRITE(*, '(2X,''GIVE FRACT COORDS OF SECOND ATOM'')') READ(5,*,ERR=160)X(2),Y(2),Z(2) WRITE(*, '(2X,''DO YOU WANT THE BOND ANGLE [Y/N]?'')') READ(5,10)ANS IF(ANS.EQ.'N')GO TO 200 NUM=3 170 WRITE(*, '(2X,''GIVE COORDS OF THIRD ATOM'')') READ(5,*,ERR=170)X(3),Y(3),Z(3) 200 DO 3 M=1,NUM XA(M,1)=X(M)*RLC(1)+Y(M)*RLC(2)*CG+Z(M)*RLC(3)*CB XA(M,2)=Y(M)*RLC(2)*SG+Z(M)*RLC(3)*CR XA(M,3)=Z(M)*RLC(3)*CP 3 CONTINUE NM=NUM-1 DO 6 I=1,NM MM=I+1 DO 6 J=MM,NUM DX=XA(I,1)-XA(J,1) DY=XA(I,2)-XA(J,2) DZ=XA(I,3)-XA(J,3) D(I,J)=SQRT(DX*DX+DY*DY+DZ*DZ) WRITE(*,20)I,J,D(I,J) 20 FORMAT(2X,'BOND LENGTH ',I3,'-'I1,' IS ',F8.3,' A') 6 CONTINUE IF(NUM.EQ.3)GO TO 250 275 WRITE(*, '(2X,''DO YOU WISH TO QUIT [Y/N]?'')') READ(5,10)ANS IF(ANS.EQ.'Y')GO TO 300 GO TO 100 250 A=DPR*ACOS((D(1,2)**2+D(2,3)**2-D(1,3)**2)/(2.*D(1,2)*D(2,3))) WRITE(*,25)A 25 FORMAT(2X,'BOND ANGLE IS ',F6.2,' DEG.') GO TO 275 300 STOP END C2345678 DSPACING FROM KEYED-IN HKL CHARACTER*1 ANS DPR=57.296 OPEN(7,FILE='CELL.OUT') WRITE(*,*)' THE CURRENT PARAMETERS & WAVELENGTH ARE:' READ(7,40)A,B,C,AL,BE,GA,W 40 FORMAT(7F9.4) WRITE(*,*)A,B,C,AL,BE,GA,W WRITE(*,*)' CORRECT?' READ(*,20)ANS IF(ANS.EQ.'Y')GO TO 50 WRITE(*,*)' GIVE A,B, ETC., AND WAVELENGTH' READ(*,*)A,B,C,AL,BE,GA,W REWIND 7 WRITE(7,40)A,B,C,AL,BE,GA,W 50 CA=COS(AL/DPR) SA=SIN(AL/DPR) CB=COS(BE/DPR) SB=SIN(BE/DPR) CG=COS(GA/DPR) SG=SIN(GA/DPR) TOP=1.-CA*CA-CB*CB-CG*CG+2.*CA*CB*CG 30 WRITE(*,*)' GIVE H,K,L' READ(*,*)IH,IK,IL FH=FLOAT(IH) FK=FLOAT(IK) FL=FLOAT(IL) BOT=(FH*SA/A)**2+(FK*SB/B)**2+(FL*SG/C)**2 1+2.*FK*FL*((CB*CG-CA)/(B*C)) 2+2.*FL*FH*((CG*CA-CB)/(C*A)) 3+2.*FH*FK*((CA*CB-CG)/(A*B)) D=SQRT(TOP/BOT) WRITE(*,*)D TT=2.*DPR*ASIN(W/(2.*D)) WRITE(*,10)IH,IK,IL,D,TT 10 FORMAT(3I4,2F8.2) WRITE(*,*)' ANOTHER Y/N?' READ(*,20)ANS 20 FORMAT(A1) IF(ANS.EQ.'Y')GO TO 30 STOP END C234567PGM TO GIVE CELL AXIAL LENGTH FROM OSC/ROT PHOTO CHARACTER*1 ANS WRITE(*,*)' PGM TO CALCULATE AXIAL LENGTH FROM POSITIONS OF' WRITE(*,*)' LAYER LINES' WRITE(*,*)' GIVE WAVELENGTH:' READ(*,*)W WRITE(*,*)' GIVE RADIUS OF CAMERA. WEISSENBERG IS 28.65 mm' READ(*,*)R 30 WRITE(*,*)' GIVE SERIAL NUMBER OF LAYER LINE. EQUATOR =0' READ(*,*)N WRITE(*,*)' GIVE DISTANCE OF THIS LAYER LINE FROM EQUATOR (mm)' READ(*,*)S A=(W*N)/SIN(ATAN(S/R)) WRITE(*,10)A 10 FORMAT(' AXIAL LENGTH = ',F9.2,' A') WRITE(*,*)' ANOTHER LAYER LINE Y/N?' READ(*,20)ANS 20 FORMAT(A1) IF(ANS.EQ.'Y')GO TO 30 STOP END C ORTH.FOR ADAPTED FROM MINITRAB.FOR DIMENSION RLC(3),X(3),Y(3),Z(3),XA(3,3),D(3,3) CHARACTER*1 ANS DPR=57.29578 OPEN (7,FILE='MINITRAB.DAT',STATUS='OLD') READ(7,*)RLC,CA,CB,CG WRITE(*,*)'THESE ARE THE CELL PARAMETERS ON MINITRAB.DAT' WRITE(*,21)RLC,CA,CB,CG WRITE(*, '(2X,''ARE THESE CORRECT [Y/N]?'')') READ(*,10)ANS 10 FORMAT(A1) IF(ANS.EQ.'Y')GO TO 101 50 WRITE(*, '(2X,''PROVIDE THEM IN FREE FORMAT'')') READ(5,*,ERR=50)RLC,CA,CB,CG BACKSPACE 7 WRITE(7,21)RLC,CA,CB,CG 21 FORMAT(3F9.3,3F9.2) 101 CA=COS(CA/DPR) CB=COS(CB/DPR) CG=COS(CG/DPR) SG=SQRT(1.-CG*CG) CR=(CA-CB*CG)/SG CP=SQRT(1.-CB*CB-CR*CR) 150 WRITE(*, '(2X,''GIVE FRACT COORDS OF ATOM'')') READ(5,*,ERR=150)X(1),Y(1),Z(1) M=1 XA(M,1)=X(M)*RLC(1)+Y(M)*RLC(2)*CG+Z(M)*RLC(3)*CB XA(M,2)=Y(M)*RLC(2)*SG+Z(M)*RLC(3)*CR XA(M,3)=Z(M)*RLC(3)*CP WRITE(*,3) XA(M,1),XA(M,2),XA(M,3) 3 FORMAT(10X,3F9.4) 275 WRITE(*, '(2X,''DO YOU WISH TO QUIT [Y/N]?'')') READ(5,10)ANS IF(ANS.EQ.'Y')GO TO 300 GO TO 150 300 STOP END C2345678 REAL TO RECIP CELL DPR=57.296 WRITE(*,*)' GIVE A,B, ETC.' READ(*,*)A,B,C,AL,BE,GA CA=COS(AL/DPR) SA=SIN(AL/DPR) CB=COS(BE/DPR) SB=SIN(BE/DPR) CG=COS(GA/DPR) SG=SIN(GA/DPR) V=SQRT(1.-CA*CA-CB*CB-CG*CG+2.*CA*CB*CG) V=A*B*C*V AS=B*C*SA/V BS=A*C*SB/V CS=A*B*SG/V ALS=DPR*ACOS((CB*CG-CA)/(SB*SG)) BES=DPR*ACOS((CA*CG-CB)/(SA*SG)) GAS=DPR*ACOS((CA*CB-CG)/(SA*SB)) WRITE(*,*)AS,BS,CS,ALS,BES,GAS STOP END C234567 GIVEN THE ENDS OF TWOO CHORDS FIND CENTRE OF CIRCLE DIMENSION X(4),Y(4) CHARACTER*1 ANS 20 WRITE(*,*)' ENTER FOUR LINES. X,Y VALUES FOR' WRITE(*,*)' THE PAIRS OF TERMINI OF TWO CHORDS' DO 100 I=1,4 READ(*,*)X(I),Y(I) 100 CONTINUE C NORMAL TO FIRST CHORD SL1=-(X(2)-X(1))/(Y(2)-Y(1)) C MIDPOINT XM1=(X(1)+X(2))/2. YM1=(Y(1)+Y(2))/2. C SECOND CHORD SL2=-(X(4)-X(3))/(Y(4)-Y(3)) XM2=(X(3)+X(4))/2. YM2=(Y(3)+Y(4))/2. C1=YM1-SL1*XM1 C2=YM2-SL2*XM2 XC=(C2-C1)/(SL1-SL2) YC=SL1*XC+C1 WRITE(*,10)XC,YC 10 FORMAT(2F7.1) WRITE(*,*)' ANOTHER PAIR OF CHORDS, Y/N?' READ(*,11)ANS 11 FORMAT(A1) IF(ANS.EQ.'Y')GO TO 20 STOP END C2345678 REFORMAT A SHELX .HKL FILE OPEN(7,FILE=' ') OPEN(8,FILE=' ') 20 READ(7,10,END=999)IH,IK,IL,FS,SFS 10 FORMAT(3I4,2F8.2) WRITE(8,11)IH,IK,IL,FS,SFS 11 FORMAT(3I4,2F9.2) GO TO 20 999 STOP END C234567 BIN\SP3.FOR FIND TWO METHYLENE HYDROGEN ATOMS DIMENSION RLC(3),RAN(3),FX(3),FY(3),FZ(3),X(3),Y(3),Z(3) DATA RLC/4.992,7.503,67.448/ DATA RAN/90.,90.,90./ DPR=57.2958 CA=COS(RAN(1)/DPR) CB=COS(RAN(2)/DPR) CG=COS(RAN(3)/DPR) SG=SIN(RAN(3)/DPR) CR=(CA-CB*CG)/SG CP=SQRT(1.-CB*CB-CR*CR) WRITE(*,*)' GIVE FRACT COORDS OF THREE CONSEC ATOMS' READ(*,*)FX(1),FY(1),FZ(1) READ(*,*)FX(2),FY(2),FZ(2) READ(*,*)FX(3),FY(3),FZ(3) DO 100 N=1,3 X(N)=FX(N)*RLC(1)+FY(N)*RLC(2)*CG+FZ(N)*RLC(3)*CB Y(N)=FY(N)*RLC(2)*SG+FZ(N)*RLC(3)*CR 100 Z(N)=FZ(N)*RLC(3)*CP C SP3 ROUTINE AB=SQRT((X(1)-X(2))**2+(Y(1)-Y(2))**2+(Z(1)-Z(2))**2) BC=SQRT((X(2)-X(3))**2+(Y(2)-Y(3))**2+(Z(2)-Z(3))**2) R=BC/AB XG=X(2)-((X(2)-X(3))/R) YG=Y(2)-((Y(2)-Y(3))/R) ZG=Z(2)-((Z(2)-Z(3))/R) XP=(X(1)+XG)/2. YP=(Y(1)+YG)/2. ZP=(Z(1)+ZG)/2. PB=SQRT(( C *** B M F I T *** BEST MOLECULAR FIT *** S. C. NYBURG ** C ************************* UPDATE OCT 1981 *********** IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 TYPE CHARACTER*8 PTYPE INTEGER*2 ORD COMMON TYPE(2,99),N,J,DS,NA,NB COMMON D(99), U(3,3),W(99),X(2,99),Y(2,99),Z(2,99),PTYPE(99) DIMENSION CQ(3,3),A(9),R(9),ORD(99),XP(99),YP(99),ZP(99),EPX(99) DIMENSION CA(2),CB(2),CG(2),SG(2),CR(2),CP(2),EX(2,99) DIMENSION EY(2,99),EZ(2,99),FRX(2,99),FRY(2,99),FRZ(2,99),RLC(2,3) DIMENSION XS(2,99),YS(2,99),ZS(2,99),RAN(2,3),EPY(99),EPZ(99) DIMENSION TITLE(18) OPEN(1, FILE=' ') OPEN(2, FILE=' ') OPEN(7, FILE=' ') OPEN(8, FILE=' ') 1000 READ(7,79,END=999)TITLE 79 FORMAT(18A4) DO 100 NY=1,2 DO 100 NZ=1,50 X(NY,NZ)=0. Y(NY,NZ)=0. 100 Z(NY,NZ)=0. DS=0. J=0 C ITAP IS TH INPUT FILE HNIT FOR MOLECULE A IF=1,FOR B IF = 2 READ(7,*)N,NA,NB,NFIX,NFIT,NW,ITEMP,NDISCA,NDISCB,NCUTA,NCUTB,NNN C** 1 FORMAT(12I5) IF(NA-NB)101,101,103 101 NTOT=NB GO TO 10 103 NTOT=NA 10 CONTINUE DPR=57.2957795 DO 33 K=1,2 C ****** RAN IN EITHER COSINES OR DEGREES READ(7,* )RLC(K,1),RLC(K,2),RLC(K,3),RAN(K,1),RAN(K,2),RAN(K,3) C**12 FORMAT(6F10.6) IF (RAN(K,1).LT.1.0) GO TO 102 RAN(K,1)=COS(RAN(K,1)/DPR) RAN(K,2)=COS(RAN(K,2)/DPR) RAN(K,3)=COS(RAN(K,3)/DPR) 102 CA(K) = RAN(K,1) CB(K) = RAN(K,2) CG(K) = RAN(K,3) SG(K) = SQRT(1.0 - RAN(K,3)**2) CR(K)=(CA(K)-CB(K)*CG(K))/SG(K) CP(K)=SQRT(1.-CB(K)*CB(K)-CR(K)*CR(K)) IF(K-1)22,35,22 35 NOW=NA NOWW=NDISCA ITAP=1 IF(NCUTA.EQ.0)GO TO 2333 II=0 2033 II=II+1 READ(ITAP,*) C2133 FORMAT(A4) IF(II.EQ.NCUTA)GO TO 2333 GO TO 2033 2333 CONTINUE GO TO 25 22 NOW=NB NOWW=NDISCB ITAP=2 IF(NCUTB.EQ.0)GO TO 3333 II=0 3033 II=II+1 READ(ITAP,*) C3133 FORMAT(A4) IF(II.EQ.NCUTB)GO TO 3333 GO TO 3033 3333 CONTINUE 25 DO 33 KJ=1,NOW IF(KJ.GT.NOWW)GO TO 4033 READ(ITAP,*)TYPE(K,KJ),FRX(K,KJ),FRY(K,KJ),FRZ(K,KJ) 3 FORMAT (A6,21X,3F9.6 ) IF(ITEMP.NE.0) READ(ITAP,*) GO TO 4032 4033 READ(7,*)TYPE(K,KJ),FRX(K,KJ),FRY(K,KJ),FRZ(K,KJ) IF(ITEMP.NE.0) READ(7,3) 4032 CONTINUE X(K,KJ)=FRX(K,KJ)*RLC(K,1)+FRY(K,KJ)*RLC(K,2)*CG(K)+FRZ(K,KJ)* 1RLC(K,3)*CB(K) Y(K,KJ)=FRY(K,KJ)*RLC(K,2)*SG(K)+FRZ(K,KJ)*RLC(K,3)*CR(K) Z(K,KJ)=FRZ(K,KJ)*RLC(K,3)*CP(K) IF(NW)44,44,45 44 W(KJ)=1. GO TO 33 45 READ(7,* )EX(K,KJ),EY(K,KJ),EZ(K,KJ) EX(K,KJ)=EX(K,KJ)*RLC(K,1)+EY(K,KJ)*RLC(K,2)*CG(K)+EZ(K,KJ)* 1RLC(K,3)*CB(K) EY(K,KJ)=EY(K,KJ)*RLC(K,2)*SG(K)+EZ(K,KJ)*RLC(K,3)*CR(K) EZ(K,KJ)=EZ(K,KJ)*RLC(K,3)*CP(K) C**46 FORMAT( 27X,3F9.6) EX(K,KJ)=SQRT((EX(K,KJ)**2+EY(K,KJ)**2+EZ(K,KJ)**2)/3.) 33 CONTINUE IF(NNN.EQ.0)GO TO 7 READ(7,* )ORD C*105 FORMAT(24I3) DO 8 NORD=1,NB IF(ORD(NORD).EQ.0)ORD(NORD)=NORD PTYPE(NORD)=TYPE(2,ORD(NORD)) XP(NORD)=X(2,ORD(NORD)) YP(NORD)=Y(2,ORD(NORD)) ZP(NORD)=Z(2,ORD(NORD)) EPX(NORD)=EX(2,ORD(NORD)) EPY(NORD)=EY(2,ORD(NORD)) EPZ(NORD)=EZ(2,ORD(NORD)) 8 CONTINUE DO 9 M=1,NB TYPE(2,M)=PTYPE(M) X(2,M)=XP(M) Y(2,M)=YP(M) Z(2,M)=ZP(M) EX(2,M)=EPX(M) EY(2,M)=EPY(M) EZ(2,M)=EPZ(M) 9 CONTINUE C FIND WEIGHT FOR EACH MATCHED ATOM PAIR 7 IF(NW.EQ.0) GO TO 49 SES=0. DO 47 KJ=1,N EX(1,KJ)=SQRT((EX(1,KJ))**2+(EX(2,KJ))**2) 47 SES=SES+ (1./(EX(1,KJ)**2)) AN=N AN=AN/SES DO 48 KJ=1,N 48 W(KJ)=AN/(EX(1,KJ)**2) 49 WRITE(8,179)TITLE 179 FORMAT(1H1,20A4,//) IF(NW.EQ.0) GO TO 75 WRITE(8,76) 76 FORMAT(31X,'THIS IS A WEIGHTED FIT'//) GO TO 80 75 WRITE(8,77) 77 FORMAT(30X,'THIS IS AN UNWEIGHTED FIT'//) 80 WRITE(8,23) 23 FORMAT(72H ORTHOGONAL COORDINATES OF ATOMS 1 BEFORE MATCHING//) WRITE(8,24) 24 FORMAT(21X,'FIRST MOLECULE',31X,'SECOND MOLECULE'//) WRITE(8,20)RLC(1,1),RLC(1,2),RLC(1,3),RLC(2,1),RLC(2,2),RLC(2,3) 20 FORMAT(3X,'A=',F9.4,' B=',F9.4,' C=',F9.4,10X,'A=',F9.4,' B= 1',F9.4,' C=',F9.4) WRITE(8,43)RAN(1,1),RAN(1,2),RAN(1,3),RAN(2,1),RAN(2,2),RAN(2,3) 43 FORMAT(3X,'COSA=',F9.6,',COSB=',F9.6,',COSG=',F9.6,5X,'COSA=',F9.6 1,',COSB=',F9.6,',COSG=',F9.6//) CALL NATOMS C TYPE OF FORMAT 60 IF((NA-N).EQ.0)GO TO 50 IF((NB-N).EQ.0)GO TO 54 IF(NA-NB)51,52,53 C OUTPUT ROUTINE F5 51 NP=N+1 NAP=NA+1 IF(J.NE.0 )GO TO 84 WRITE(8,26)(MI,TYPE(1,MI),X(1,MI),Y(1,MI),Z(1,MI),TYPE(2,MI),X(2,M 1I),Y(2,MI),Z(2,MI),MI=NP,NA) 26 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6) WRITE(8,74)(MP,TYPE(2,MP),X(2,MP),Y(2,MP),Z(2,MP),MP=NAP,NB) 74 FORMAT(I4,39X,A6,3F10.6) GO TO 34 84 WRITE(8,82)(MM,TYPE(1,MM),X(1,MM),Y(1,MM),Z(1,MM),TYPE(2,MM),X(2,M 1M),Y(2,MM),Z(2,MM),D(MM),MM=NP,NA) 82 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) 85 WRITE(8,42)(MR,TYPE(2,MR),X(2,MR),Y(2,MR),Z(2,MR),MR=NAP,NB) 42 FORMAT(I4,39X,A6,3F10.6) 34 CONTINUE GO TO 55 C OUTPUT ROUTINE F6 52 NP=N+1 IF(J.NE.0 )GO TO 86 WRITE(8,27)(MS,TYPE(1,MS),X(1,MS),Y(1,MS),Z(1,MS),TYPE(2,MS),X(2,M 1S),Y(2,MS),Z(2,MS),MS=NP,NA) 27 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6) GO TO 87 86 WRITE(8,36)(MT,TYPE(1,MT),X(1,MT),Y(1,MT),Z(1,MT),TYPE(2,MT),X(2,M 1T),Y(2,MT),Z(2,MT),D(MT),MT=NP,NA) 36 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) 87 CONTINUE GO TO 55 C OUTPUT ROUTINE F4 53 NP=N+1 NBP=NB+1 IF(J.NE.0 )GO TO 81 WRITE(8,37)(MG,TYPE(1,MG),X(1,MG),Y(1,MG),Z(1,MG),TYPE(2,MG),X(2,M 1G),Y(2,MG),Z(2,MG),MG=NP,NB) 37 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6) WRITE(8,58)(MH,TYPE(1,MH),X(1,MH),Y(1,MH),Z(1,MH),MH=NBP,NA) GO TO 135 81 WRITE(8,39)(MN,TYPE(1,MN),X(1,MN),Y(1,MN),Z(1,MN),TYPE(2,MN),X(2,M 1N),Y(2,MN),Z(2,MN),D(MN),MN=NP,NB) 39 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) 83 WRITE(8,58)(ML,TYPE(1,ML),X(1,ML),Y(1,ML),Z(1,ML),ML=NBP,NA) 58 FORMAT(I4,A6,3F10.6) 135 CONTINUE GO TO 55 50 IF(NB-N)66,65,66 65 WRITE(8,29) 29 FORMAT(33H NONE//) GO TO 55 C OUTPUT ROUTINE F3 66 NBP=N+1 WRITE(8,40)(MF,TYPE(2,MF),X(2,MF),Y(2,MF),Z(2,MF),MF=NBP,NB) 40 FORMAT(I4,39X,A6,3F10.6) GO TO 55 C OUTPUT ROUTINE F2 54 NAP=N+1 WRITE(8,41)(ME,TYPE(1,ME),X(1,ME),Y(1,ME),Z(1,ME),ME=NAP,NA) 41 FORMAT(I4,A6,3F10.6) 55 IF(J.NE.0)GO TO 61 WRITE(8,90) 90 FORMAT(1H0) WRITE(8,93) 93 FORMAT(69H CENTROID 1 CENTROID/) C FIND PAIR-WEIGHTED CENTROIDS OF MOLECULES IF(NFIX.NE.0)GO TO 62 XA=0. YA=0. ZA=0. XB=0. YB=0. ZB=0. DO 4 NF=1,N XA=XA+X(1,NF)*W(NF) YA=YA+Y(1,NF)*W(NF) ZA=ZA+Z(1,NF)*W(NF) XB=XB+X(2,NF)*W(NF) YB=YB+Y(2,NF)*W(NF) 4 ZB=ZB+Z(2,NF)*W(NF) AN=N XA=XA/AN YA=YA/AN ZA=ZA/AN XB=XB/AN YB=YB/AN ZB=ZB/AN GO TO 63 62 XA=X(1,NFIX) YA=Y(1,NFIX) ZA=Z(1,NFIX) XB=X(2,NFIX) YB=Y(2,NFIX) ZB=Z(2,NFIX) 63 WRITE(8,94)XA,YA,ZA,XB,YB,ZB 94 FORMAT(10X,3F10.6,9X,3F10.6,//) WRITE(8,78) 78 FORMAT(1H1) IF(NW.EQ.0) GO TO 175 WRITE(8,76) GO TO 180 175 WRITE(8,77) C MOVE CENTROIDS TO ORIGIN 180 DO 5 NG=1,NTOT X(1,NG)=X(1,NG)-XA Y(1,NG)=Y(1,NG)-YA Z(1,NG)=Z(1,NG)-ZA X(2,NG)=X(2,NG)-XB Y(2,NG)=Y(2,NG)-YB Z(2,NG)=Z(2,NG)-ZB XS(2,NG)=X(2,NG) YS(2,NG)=Y(2,NG) 5 ZS(2,NG)=Z(2,NG) DO 997 I=1,3 DO 997 JJJ=1,3 U(I,JJJ)=0. 997 CQ(I,JJJ)=0. DO 17 I=1,N CQ(1,1)=CQ(1,1)+X(1,I )*X(2,I )*W(I) CQ(1,2)=CQ(1,2)+Y(1,I )*X(2,I )*W(I) CQ(1,3)=CQ(1,3)+Z(1,I )*X(2,I )*W(I) CQ(2,1)=CQ(2,1)+X(1,I )*Y(2,I )*W(I) CQ(2,2)=CQ(2,2)+Y(1,I )*Y(2,I )*W(I) CQ(2,3)=CQ(2,3)+Z(1,I )*Y(2,I )*W(I) CQ(3,1)=CQ(3,1)+X(1,I )*Z(2,I )*W(I) CQ(3,2)=CQ(3,2)+Y(1,I )*Z(2,I )*W(I) CQ(3,3)=CQ(3,3)+Z(1,I )*Z(2,I )*W(I) 17 CONTINUE A(1) =CQ(1,1)*CQ(1,1)+CQ(1,2)*CQ(1,2)+CQ(1,3)*CQ(1,3) A(2) =CQ(1,1)*CQ(2,1)+CQ(1,2)*CQ(2,2)+CQ(1,3)*CQ(2,3) A(4) =CQ(1,1)*CQ(3,1)+CQ(1,2)*CQ(3,2)+CQ(1,3)*CQ(3,3) A(3) =CQ(2,1)*CQ(2,1)+CQ(2,2)*CQ(2,2)+CQ(2,3)*CQ(2,3) A(5) =CQ(2,1)*CQ(3,1)+CQ(2,2)*CQ(3,2)+CQ(2,3)*CQ(3,3) A(6) =CQ(3,1)*CQ(3,1)+CQ(3,2)*CQ(3,2)+CQ(3,3)*CQ(3,3) J=1 CALL EIGEN(A,R) NWR=0 NMIR=0 CALL VECTOR(A,R,CQ,NWR,NMIR,DETA) 213 IF(NWR.NE.0) GO TO 204 IF(NMIR.EQ.1) GO TO 209 IF(DETA.LT.0.) GO TO 203 WRITE(8,211) 211 FORMAT(20X,'THE BEST FIT IS GIVEN WITH THE SECOND MOLECULE UN-MIRR 1ORED',///) GO TO 207 203 WRITE(8,212) 212 FORMAT(20X,'THE BEST FIT IS GIVEN WITH THE SECOND MOLECULE MIRRORE 1D',///) GO TO 207 204 IF(DDET.LT.0.) GO TO 205 WRITE(8,206) 206 FORMAT(1H1,35X,'THIS IS THE MIRRORED FIT',///) GO TO 2 205 WRITE(8,208) 208 FORMAT(1H1,20X,'THIS IS THE FIT WITH THE SECOND MOLECULE UN-MIRROR 1ED',///) GO TO 2 209 WRITE(8,210) 210 FORMAT( 1X,'EITHER ONE OF THE MOLECULES IS FLAT OR THE FIT IS THE 1SAME WHETHER OR NOT THE 2ND MOLECULE IS MIRRORED',///) 207 DDET=DETA 2 DO 994 I=1,NB XN=A(1)*XS(2,I)+A(2)*YS(2,I)+A(3)*ZS(2,I) YN=A(4)*XS(2,I)+A(5)*YS(2,I)+A(6)*ZS(2,I) ZN=A(7)*XS(2,I)+A(8)*YS(2,I)+A(9)*ZS(2,I) X(2,I)=XN Y(2,I)=YN Z(2,I)=ZN 994 CONTINUE CALL DSUM WRITE(8,11) 11 FORMAT(87H 1 DELTA/) CALL NATOMS GO TO 60 61 CONTINUE C THIS FINAL SECTION PERFORMS THE TRANSFORMATION OF 2ND MOLECULE INT C THE FRACTIONAL COORDINATES OF THE FIRST. C IF THIS IS REQUIRED, PUT NFIT=1. OTHERWISE NFIT=0. IF(NFIT.EQ.0) GO TO 161 C BRING 2ND MOLECULE BACK TO SAME CRYSTAL FRAME AS 1ST MOLECULE C MOVE 2ND MOLECULE TO CENTROID OF 1ST MOLECULE IN ITS ORIGINAL POSI DO 16 NG=1,NB X(2,NG)=X(2,NG)+XA Y(2,NG)=Y(2,NG)+YA Z(2,NG)=Z(2,NG)+ZA 16 CONTINUE C INVERSE MATRIX TO BRING TO FRACTIONAL COORDS DO 30 K=1,NB X(2,K)=(X(2,K)-(Y(2,K)*CG(1)/SG(1)) 1 +(Z(2,K)/CP(1))*((CR(1)*CG(1)/SG(1)) -CB(1)))/RLC(1,1) Y(2,K)= (Y(2,K)/SG(1)-((Z(2,K)*CR(1))/(CP(1)*SG(1))))/RLC(1,2) 30 Z(2,K)=Z(2,K)/(CP(1)*RLC(1,3)) WRITE(8,250) 250 FORMAT ( 1H0,///,' NO', 4X,' ATOM ',6X,'FRACT. COORDINATES', 1 11X,' ATOM FRACT. COORDS. OF 2ND MOLECULE',/23X,'FIRST MOLECULE 2',26X,'IN CELL OF FIRST') IF(NA-NB)67,67,68 67 NUM=NA GO TO 70 68 NUM=NB 70 WRITE(8,59)(M,TYPE(1,M),FRX(1,M),FRY(1,M),FRZ(1,M),TYPE(2,M),X(2,M 1),Y(2,M),Z(2,M),M=1,NUM) 59 FORMAT(I4,5X,A6,3F10.6,5X,A6,3F10.6) IF(NA-NB)71,161,72 71 NAP=NA+1 WRITE(8,99)(M,TYPE(2,M),X(2,M),Y(2,M),Z(2,M),M=NAP,NB) 99 FORMAT(I4,46X,A6,3F10.6) GO TO 161 72 NBP=NB+1 WRITE(8,73)(M,TYPE(1,M),FRX(1,M),FRY(1,M),FRZ(1,M),M=NBP,NA) 73 FORMAT(I4,5X,A6,3F10.6) 161 IF(NMIR.EQ.1) GO TO 1000 NWR=NWR+1 IF(NWR.EQ.2)GO TO 1000 CALL VECTOR(A,R,CQ,NWR,NMIR,DETA) GO TO 204 999 STOP END SUBROUTINE EIGEN(A,R) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(6),R(9) RANGE=1.0E-6 IQ=-3 DO 20 J=1,3 IQ=IQ+3 DO 20 I=1,3 IJ=IQ+I R(IJ)=0.0 IF(I-J) 20,15,20 15 R(IJ)=1.0 20 CONTINUE 25 ANORM=0.0 DO 35 I=1,3 DO 35 J=I,3 IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM)165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*RANGE/3. IND=0 THR=ANORM 45 THR=THR/3. 50 L=1 55 M=L+1 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(ABS(A(LM))-THR) 130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X =0.5*(A(LL)-A(MM)) 68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y)))) SINX2=SINX*SINX 78 COSX=SQRT(1.0-SINX2) COSX2=COSX*COSX SINCS=SINX*COSX ILQ=3*(L-1) IMQ=3*(M-1) DO 125 I=1,3 IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GO TO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GO TO 110 105 IL=L+IQ 110 X =A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 ILR=ILQ+I IMR=IMQ+I X =R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X =2.0*A(LM)*SINCS Y =A(LL)*COSX2+A(MM)*SINX2-X X =A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X 130 IF(M-3) 135,140,135 135 M=M+1 GO TO 60 140 IF(L-2) 145,150,145 145 L=L+1 GO TO 55 150 IF(IND-1) 160,155,160 155 IND=0 GO TO 50 160 IF(THR-ANRMX) 165,165,45 165 IQ=-3 DO 185 I=1,3 IQ=IQ+3 LL=I+(I*I-I)/2 JQ=3*(I-2) DO 185 J=I,3 JQ=JQ+3 MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X =A(LL) A(LL)=A(MM) A(MM)=X DO 180 K=1,3 ILR=IQ+K IMR=JQ+K X =R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE RETURN END SUBROUTINE VECTOR(A,R,CQ,NWR,NMIR,DETA) IMPLICIT REAL*8 (A-H, O-Z) DIMENSION CQ(3,3) DIMENSION A(9),R(9) DIMENSION U(3,3) DIMENSION WV(3),ZZ(3,3),RC(3,3),B(3,3) DO 28 I=1,3 DO 28 J=1,3 28 U(I,J)=0. IF(NWR.NE.0) GO TO 52 WV(1)=A(1) WV(2)=A(3) WV(3)=A(6) ZZ(1,1)=R(1) ZZ(1,2)=R(2) ZZ(1,3)=R(3) ZZ(2,1)=R(4) ZZ(2,2)=R(5) ZZ(2,3)=R(6) ZZ(3,1)=R(7) ZZ(3,2)=R(8) ZZ(3,3)=R(9) DO 21 I=1,3 DO 21 J=1,3 B(I,J)=0. 21 RC(I,J)=CQ(J,I) C MAKE SURE DET OF RC IS NOT EXACTLY ZERO 999 FORMAT(/' VALUE OF DETERMINANT',E14.4) DETR= RC(1,1)*RC(2,2)*RC(3,3) $ -RC(1,1)*RC(2,3)*RC(3,2) $ +RC(1,2)*RC(2,3)*RC(3,1) $ -RC(1,2)*RC(2,1)*RC(3,3) $ +RC(1,3)*RC(2,1)*RC(3,2) $ -RC(1,3)*RC(2,2)*RC(3,1) WRITE(8,999)DETR IF(ABS(DETR).LT.0.0001 )GO TO 17 DO 19 K=1,3 DO 19 I=1,3 DO 19 J=1,3 B(K,I)=B(K,I)+RC(I,J)*ZZ(K,J) 19 CONTINUE GO TO 38 17 NMIR=1 IF(WV(2).LT.1.0E-6) GO TO 39 DO 40 K=1,2 DO 40 I=1,3 DO 40 J=1,3 B(K,I)=B(K,I)+RC(I,J)*ZZ(K,J) 40 CONTINUE B (3,1)=B (1,2)*B (2,3)-B (1,3)*B (2,2) B (3,2)=B (1,3)*B (2,1)-B (2,3)*B (1,1) B (3,3)=B (1,1)*B (2,2)-B (2,1)*B (1,2) 38 BS= B(1,1)**2+B(1,2)**2+B(1,3)**2 IF (BS.EQ.0.0) GO TO 100 FACS=1./(B(1,1)**2+B(1,2)**2+B(1,3)**2) DO 13 J=1,3 13 B(1,J)=SQRT(FACS)*B(1,J) 100 BS= B(2,1)**2+B(2,2)**2+B(2,3)**2 IF(BS.EQ.0.0) GO TO 101 FACS=1./(B(2,1)**2+B(2,2)**2+B(2,3)**2) DO 14 J=1,3 14 B(2,J)=SQRT(FACS)*B(2,J) 101 BS= B(3,1)**2+B(3,2)**2+B(3,3)**2 IF(BS.EQ.0.0) GO TO 102 FACS=1./(B(3,1)**2+B(3,2)**2+B(3,3)**2) DO 15 J=1,3 15 B(3,J)=SQRT(FACS)*B(3,J) GO TO 102 39 DO 41 I=1,3 DO 41 J=1,3 B(1,I)=B(1,I)+RC(I,J)*ZZ(1,J) 41 CONTINUE FACS=1./(B(1,1)**2+B(1,2)**2+B(1,3)**2) DO 42 J=1,3 B(1,J)=SQRT(FACS)*B(1,J) 42 CONTINUE IF(B(1,1).EQ.0.) GO TO 43 IF(B(1,2).EQ.0.) GO TO 37 B(2,2)=B(1,1)/SQRT(B(1,1)**2+B(1,2)**2) B(2,1)=-(B(1,2)*B(2,2))/B(1,1) B(2,3)=0. GO TO 44 43 B(2,1)=1. B(2,2)=0. B(2,3)=0. GO TO 44 37 B(2,1)=0. B(2,2)=1. B(2,3)=0. 44 B(3,1)=B(1,2)*B(2,3)-B(1,3)*B(2,2) B(3,2)=B(1,3)*B(2,1)-B(1,1)*B(2,3) B(3,3)=B(1,1)*B(2,2)-B(1,2)*B(2,1) GO TO 102 52 B(3,1)=-B(3,1) B(3,2)=-B(3,2) B(3,3)=-B(3,3) 102 DO 20 I=1,3 DO 20 J=1,3 DO 20 K=1,3 20 U(I,J)=U(I,J)+B(K,I)*ZZ(K,J) A(1)=U(1,1) A(2)=U(1,2) A(3)=U(1,3) A(4)=U(2,1) A(5)=U(2,2) A(6)=U(2,3) A(7)=U(3,1) A(8)=U(3,2) A(9)=U(3,3) C IS A MATRIX ORTHOGONAL? A11=A(1)*A(1)+A(2)*A(2)+A(3)*A(3) A12=A(1)*A(4)+A(2)*A(5)+A(3)*A(6) A13=A(1)*A(7)+A(2)*A(8)+A(3)*A(9) A21=A(4)*A(1)+A(5)*A(2)+A(6)*A(3) A22=A(4)*A(4)+A(5)*A(5)+A(6)*A(6) A23=A(4)*A(7)+A(5)*A(8)+A(6)*A(9) A31=A(7)*A(1)+A(8)*A(2)+A(9)*A(3) A32=A(7)*A(4)+A(8)*A(5)+A(9)*A(6) A33=A(7)*A(7)+A(8)*A(8)+A(9)*A(9) WRITE(8,998)A11,A12,A13 WRITE(8,998)A21,A22,A23 WRITE(8,998)A31,A32,A33 998 FORMAT(/20X,3E14.4,/) DETA=A(1)*A(5)*A(9)-A(1)*A(6)*A(8)+A(2)*A(6)*A(7)-A(2)*A(4)*A(9)+A 1(3)*A(4)*A(8)-A(3)*A(5)*A(7) RETURN END SUBROUTINE DSUM IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 TYPE COMMON TYPE(2,99),N,J,DS,NA,NB COMMON D(99), U(3,3),W(99),X(2,99),Y(2,99),Z(2,99) DIMENSION A(9),R(9) DS=0. IF(NA-NB)96,96,98 96 NSUM=NA GO TO 99 98 NSUM=NB 99 DO 60 MB=1,NSUM 60 D(MB)= (X(1,MB)-X(2,MB))**2+(Y(1,MB)-Y(2,MB))**2+(Z(1,MB 1)-Z(2,MB))**2 DO 160 MU=1,N 160 DS=DS+D(MU) DO 161 MV=1,NSUM 161 D(MV)=SQRT(D(MV)) RETURN END SUBROUTINE NATOMS IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 TYPE COMMON TYPE(2,99),N,J,DS,NA,NB COMMON D(99), U(3,3),W(99),X(2,99),Y(2,99),Z(2,99) IF(J.NE.0 )GO TO 56 WRITE(8,26)(MA,TYPE(1,MA),X(1,MA),Y(1,MA),Z(1,MA),TYPE(2,MA),X(2,M 1A),Y(2,MA),Z(2,MA),MA=1,N) 26 FORMAT(I4,A6,3F10.6,12X,A6,3F10.6) GO TO 57 56 WRITE(8,27)(MD,TYPE(1,MD),X(1,MD),Y(1,MD),Z(1,MD),TYPE(2,MD),X(2,M 1D),Y(2,MD),Z(2,MD),D(MD),MD=1,N) 27 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) WRITE(8,164) DS 164 FORMAT(/,' SUM OF DIFFERENCES SQUAR 1ED ',F10.6,//) 57 WRITE (8,29) 29 FORMAT(2H0 ,///) WRITE (8,28) 28 FORMAT(56H UNMATCHED ATOMS 1//) RETURN END C23456789 MATINV.FOR INVERSE OF 3 X 3 WRITE(*,*)' READ IN THE THREE ROWS OF ELEMENTS' READ(*,*)A11,A12,A13 READ(*,*)A21,A22,A23 READ(*,*)A31,A32,A33 DET=A11*A22*A33-A11*A23*A32+A12*A23*A31 1-A12*A21*A33+A13*A21*A32-A13*A22*A31 IF(DET.EQ.0.0)GO TO 90 B11=( A22*A33-A23*A32)/DET B12=(-A12*A33+A13*A32)/DET B13=( A12*A23-A13*A22)/DET B21=(-A21*A33+A23*A31)/DET B22=( A11*A33-A13*A31)/DET B23=(-A11*A23+A13*A21)/DET B31=( A21*A32-A22*A31)/DET B32=(-A11*A32+A12*A31)/DET B33=( A11*A22-A12*A21)/DET WRITE(*,*)' INVERSE IS' WRITE(*,10)B11,B12,B13 WRITE(*,10)B21,B22,B23 WRITE(*,10)B31,B32,B33 10 FORMAT(3F9.2,/) GO TO 100 90 WRITE(*,*)' DETERMINANT IS ZERO' 100 STOP END C ROTATE.FOR C TO TAKE A SET OF POINTS IN 3D (SAY RECIP LATTICE) C AND ROTATE ABOUT EITHER X,Y OR Z THRU THE ORIGIN C TO GENERATE NEW COORDINATES. ORGINAL IN UNIT 7 C IS LEFT UNCHANGED. OUTPUT GOES TO UNIT 8 AND THREE C234567 EXCEL PLOT FILES PLX.DAT,PLY.DAT,PLZ.DAT IN 9,10,11 COMMON X(100),Y(100),Z(100),XN(100),YN(100),ZN(100) CHARACTER*1 ANS WRITE(*,*)' GIVE NAME OF FILE CONTAINING COORDS TO BE ROTATED' OPEN(7,FILE=' ') WRITE(*,*)' GIVE NAME OF ROTATED DATA' OPEN(8,FILE=' ') OPEN(9,FILE='PLX.DAT') OPEN(10,FILE='PLY.DAT') OPEN(11,FILE='PLZ.DAT') NUM=0 21 FORMAT(3F8.4) DO 10 I=1,100 READ(7,21,END=999)X(I),Y(I),Z(I) NUM=NUM+1 10 CONTINUE 999 DO 12 I=1,NUM WRITE(8,21)X(I),Y(I),Z(I) 12 CONTINUE GO TO 99 30 REWIND 8 DO 11 I=1,NUM READ(8,21)X(I),Y(I),Z(I) 11 CONTINUE 99 WRITE(*,*)' GIVE AXIS X,Y, OR Z ABOUT WHICH ROTN MADE' READ(*,20)ANS 20 FORMAT(A1) WRITE(*,*)' AXES ARE RIGHT-HANDED; LOOKING ALONG Z' WRITE(*,*)' FROM -VE TO +VE, CLOCKWISE BRINGS X TO Y' WRITE(*,*)' GIVE SIZE OF ANGLE, +VE FOR CLOCK, -VE FOR ANTI' READ(*,*)ANG ANG=ANG/57.296 IF(ANS.EQ.'X')CALL ROT1(NUM,ANG) IF(ANS.EQ.'Y')CALL ROT2(NUM,ANG) IF(ANS.EQ.'Z')CALL ROT3(NUM,ANG) WRITE(*,*)' ANOTHER ROTATION?' READ(*,20)ANS IF(ANS.EQ.'Y')GO TO 30 C MAKE PLX.DAT, ETC., COVER THE RANGE +- 0.6,+-0.6 (SQUARE PLOT) X1=.6 Y1=0. X2=-.6 Y2=0. X3=0. Y3=.6 X4=0. Y4=-.6 WRITE(9,21)X1,Y1 WRITE(10,21)X1,Y1 WRITE(11,21)X1,Y1 WRITE(9,21)X2,Y2 WRITE(10,21)X2,Y2 WRITE(11,21)X2,Y2 WRITE(9,21)X3,Y3 WRITE(10,21)X3,Y3 WRITE(11,21)X3,Y3 WRITE(9,21)X4,Y4 WRITE(10,21)X4,Y4 WRITE(11,21)X4,Y4 DO 40 I=1,NUM WRITE(9,22)ZN(I),YN(I) WRITE(10,22)XN(I),ZN(I) WRITE(11,22)XN(I),YN(I) 40 CONTINUE 22 FORMAT(2F8.4) STOP END SUBROUTINE ROT1(NUM,ANG) COMMON X(100),Y(100),Z(100),XN(100),YN(100),ZN(100) REWIND 8 DO 10 I=1,NUM READ(8,20)X(I),Y(I),Z(I) XN(I)=X(I) YN(I)=Y(I)*COS(ANG)-Z(I)*SIN(ANG) ZN(I)=Y(I)*SIN(ANG)+Z(I)*COS(ANG) 10 CONTINUE REWIND 8 DO 11 I=1,NUM WRITE(8,20)XN(I),YN(I),ZN(I) 11 CONTINUE 20 FORMAT(3F8.4) RETURN END SUBROUTINE ROT2(NUM,ANG) COMMON X(100),Y(100),Z(100),XN(100),YN(100),ZN(100) REWIND 8 DO 10 I=1,NUM READ(8,20)X(I),Y(I),Z(I) XN(I)=X(I)*COS(ANG)-Z(I)*SIN(ANG) YN(I)=Y(I) ZN(I)=X(I)*SIN(ANG)+Z(I)*COS(ANG) 10 CONTINUE REWIND 8 DO 11 I=1,NUM WRITE(8,20)XN(I),YN(I),ZN(I) 11 CONTINUE 20 FORMAT(3F8.4) RETURN END SUBROUTINE ROT3(NUM,ANG) COMMON X(100),Y(100),Z(100),XN(100),YN(100),ZN(100) REWIND 8 DO 10 I=1,NUM READ(8,20)X(I),Y(I),Z(I) XN(I)=X(I)*COS(ANG)-Y(I)*SIN(ANG) YN(I)=X(I)*SIN(ANG)+Y(I)*COS(ANG) ZN(I)=Z(I) 10 CONTINUE REWIND 8 DO 11 I=1,NUM WRITE(8,20)XN(I),YN(I),ZN(I) 11 CONTINUE 20 FORMAT(3F8.4) RETURN END