program trans c calculate orientation angles of reduced cell from old cell and phis, c and new cell and transformation matrix double precision car(3,3), ar(3,3), c(3,3) double precision at1(3,3),att1(3,3),at1inv(3,3),t(3,3),det character*80 fname pi=3.14159 cpi=pi/180.0 det=0.0 write(*,100) 100 format(1x,'name of input data file ?',$) read(*,110) fname 110 format(a) open(unit=8,file=fname,status='old') c c read in old cell and orientation read(8,*)as,bs,cs,alps,bets,gams write(6,*) write(6,*)'gn3 cell a*,b*,c*,alpha*,beta*,gamma*' write(6,200)as,bs,cs,alps,bets,gams alps=alps*cpi bets=bets*cpi gams=gams*cpi call setm(as,bs,cs,alps,bets,gams,ar) read(8,*)px,py,pz write(6,*)'gn3 orientation PX,PY,PZ' write(6,150)px,py,pz 150 format(3f10.3) call cwrite(px,py,pz,c) call matmult(c,ar,car) c c read in cell from tracer c read(8,*)as,bs,cs,alps,bets,gams write(6,*) write(6,*) write(6,*)'reduced (tracer) cell a*,b*,c*,alpha*,beta*,gamma*' write(6,200)as,bs,cs,alps,bets,gams 200 format(1x,6f10.5) alps=alps*cpi bets=bets*cpi gams=gams*cpi call setm(as,bs,cs,alps,bets,gams,at1) CD write(6,*)' matrix at1..standard orientation of tracer cell' CD call matprint(at1) c c read in transformation matrix and calc orientation angles c read(8,*)(t(1,i),i=1,3) read(8,*)(t(2,i),i=1,3) read(8,*)(t(3,i),i=1,3) write(6,*)'tracer transformation matrix (cell(1) -> cell(4))' call matprint(t) call matmult(at1,t,att1) call invert3(att1,at1inv,det) CD write(6,*)' matrix at1inv....inverse of [at1].[t]' CD call matprint(at1inv) call matmult(car,at1inv,c) CD write(6,*)' matrix c' CD call matprint(c) call solvps(c,px,py,pz) write(6,*) write(6,*)'missorientation of reduced cell ' write(6,900)px,py,pz 900 format(' PX, PY, PZ =',3f12.3) write(6,*)' ' write(6,*)'with respect to c* along spindle and a*' write(6,*)'nearest to x-ray beam in plane defined by' write(6,*)'the beam axis and rotation axis' c stop end C C ************************************************************* SUBROUTINE SOLVPS(AM,PX,PY,PZ) C GIVEN A ROTATION MATRIX AM, EXTRACT THE CORRESPONDING PX,PY,PZ DOUBLE PRECISION AM(3,3),det,aminv(3,3) PI=3.14159 CPI=PI/180 det=0.0 c c check determinant if .lt. 0 assume =-1 and multiply all c matrix elements by -1.0 to get correct orientation c c if not will get 180+px,-py,180+pz c call invert3(am,aminv,det) if (det.lt.0.0) then do 10 i=1,3 do 10 j=1,3 am(i,j)=-am(i,j) 10 continue endif c c SY=-AM(3,1) IF (SY.GE.1.000) SY=0.999999 IF (SY.LE.-1.00) SY=-0.999999 CY=SQRT(1-SY*SY) IF (CY.LT.0.000001) CY=0.000001 CZ=AM(1,1)/CY SZ=AM(2,1)/CY SX=AM(3,2)/CY CX=AM(3,3)/CY PX=ATAN2(SX,CX)/CPI PY=ATAN2(SY,CY)/CPI PZ=ATAN2(SZ,CZ)/CPI C RETURN END C C ******************************************************* SUBROUTINE SETM(AS,BS,CS,ALPS,BETS,GAMS,AM) C SETS UP ORTHOGONALISATION MATRIX FROM ASTAR TO GAMMASTAR C CSTAR PARALLEL TO Z, ASTAR NEAR x,angles in radians DOUBLE PRECISION AM(3,3) CA=COS(ALPS) SA=SIN(ALPS) CB=COS(BETS) SB=SIN(BETS) CG=COS(GAMS) SG=SIN(GAMS) COSALP=(CB*CG-CA)/(SB*SG) SINALP=SQRT(1-COSALP**2) COSBET=(CG*CA-CB)/(SG*SA) SINBET=SQRT(1-COSBET**2) COSGAM=(CA*CB-CG)/(SA*SB) SINGAM=SQRT(1-COSGAM**2) AM(1,1)=AS*SB AM(2,1)=0 AM(3,1)=AS*CB AM(1,2)=-BS*SA*COSGAM AM(2,2)=BS*SA*SINGAM AM(3,2)=BS*CA AM(1,3)=0 AM(2,3)=0 AM(3,3)=CS RETURN END C C *************************************************** SUBROUTINE MATPRINT(AM) DOUBLE PRECISION AM(3,3) DO 10 I=1,3 WRITE(6,900)(AM(I,J),J=1,3) 900 FORMAT(3F12.6) 10 CONTINUE RETURN END C C ********************************************************** SUBROUTINE MATMULT(A,B,C) C 3 X 3 MATRICES, AXB INTO C DOUBLE PRECISION A(3,3), B(3,3), C(3,3) DO 20 I=1,3 DO 10 J=1,3 C(I,J)=A(I,1)*B(1,J) + A(I,2)*B(2,J) + A(I,3)*B(3,J) 10 CONTINUE 20 CONTINUE RETURN END FUNCTION ANGLE( X,Y,Z,I,J,K) C C USES X,Y,Z ARRAYS, WORKS AND PRINTS ANGLE C ATOM I-ATOM J-ATOM K C (TAKEN FROM PROGRAM COCALC AND MODIFIED, JUNE 89) C DOUBLE PRECISION X(0:10),Y(0:10),Z(0:10) PI=3.14159 C WRITE(6,921) I, X(I),Y(I),Z(I) C WRITE(6,921) J, X(J),Y(J),Z(J) C WRITE(6,921) K, X(K),Y(K),Z(K) c 921 FORMAT(2X,I6,3F10.4) DISTIJ = (X(I)-X(J))**2+(Y(I)-Y(J))**2+(Z(I)-Z(J))**2 DISTIJ = SQRT(DISTIJ) DISTJK = (X(J)-X(K))**2+(Y(J)-Y(K))**2+(Z(J)-Z(K))**2 DISTJK = SQRT(DISTJK) DISTKI = (X(K)-X(I))**2+(Y(K)-Y(I))**2+(Z(K)-Z(I))**2 DISTKI = SQRT(DISTKI) C WRITE(6,930) DISTIJ,DISTJK,DISTKI C 930 FORMAT(3F10.3) COSANG = DISTIJ**2+DISTJK**2-DISTKI**2 COSANG = COSANG/(2*DISTIJ*DISTJK) IF (COSANG.LE.-1.0) COSANG=-0.999999 IF (COSANG.GE.1.0) COSANG=0.9999999 ANG = ACOS(COSANG) WRITE(6,900)I,J,K ,ANG*180/PI 900 FORMAT(10X,' ANGLE',2X,I2,'-',I2,'-',I2,'IS ',F6.1) WRITE(6,910) DISTIJ,DISTJK,DISTKI 910 FORMAT(' DISTANCES',3F10.3) C ANGLE=ANG RESULT=ANG RETURN END C ****************************************************** C subroutine invert3(a,ainv,det) c invert a 3x3 matrix ..return inverted matrix and determinant c double precision a(3,3),ainv(3,3),acof(3,3),det c c calculate cofactors: c acof(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2)) acof(1,2)=-1.0*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) acof(1,3)=(a(2,1)*a(3,2)-a(2,2)*a(3,1)) acof(2,1)=-1.0*(a(1,2)*a(3,3)-a(1,3)*a(3,2)) acof(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1)) acof(2,3)=-1.0*(a(1,1)*a(3,2)-a(1,2)*a(3,1)) acof(3,1)=(a(1,2)*a(2,3)-a(1,3)*a(2,2)) acof(3,2)=-1.0*(a(1,1)*a(2,3)-a(1,3)*a(2,1)) acof(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1)) c c calculate determinant c det=a(1,1)*acof(1,1)+a(1,2)*acof(1,2)+a(1,3)*acof(1,3) c c transpose and divide by determinant c do 10 i=1,3 do 10 j=1,3 ainv(i,j)=acof(j,i)/det 10 continue c return end C c*************************************************************** subroutine cwrite(px,py,pz,c) c writes the missetting matrix [c] from a knowledge of px,py,pz c real px,py,pz,cpx,spx,cpy,spy,cpz,spz double precision cx(3,3),cy(3,3),cz(3,3),c(3,3),czy(3,3) piby180=3.1415926536/180.0 c cpx=cos(px*piby180) spx=sin(px*piby180) cpy=cos(py*piby180) spy=sin(py*piby180) cpz=cos(pz*piby180) spz=sin(pz*piby180) c cx(1,1)=1.0 cx(1,2)=0.0 cx(1,3)=0.0 cx(2,1)=0.0 cx(2,2)=cpx cx(2,3)=-spx cx(3,1)=0.0 cx(3,2)=spx cx(3,3)=cpx c cy(1,1)=cpy cy(1,2)=0.0 cy(1,3)=spy cy(2,1)=0.0 cy(2,2)=1.0 cy(2,3)=0.0 cy(3,1)=-spy cy(3,2)=0.0 cy(3,3)=cpy c cz(1,1)=cpz cz(1,2)=-spz cz(1,3)=0.0 cz(2,1)=spz cz(2,2)=cpz cz(2,3)=0.0 cz(3,1)=0.0 cz(3,2)=0.0 cz(3,3)=1.0 c call matmult(cz,cy,czy) call matmult(czy,cx,c) c return end