PROGRAM GN3 C DERIVE RELATIVE UNIT CELL AND MISORIENTATION ANGLES, PX,PY,PZ C FROM SIX POINTS IN GNOMONIC PROJECTION DOUBLE PRECISION CAR(3,3), AR(3,3), ARINV(3,3),C(3,3),det DOUBLE PRECISION XG(0:6),YG(0:6),ZG(0:6) DOUBLE PRECISION XL(0:6),DC(0:6,0:6) DOUBLE PRECISION CELL(6),RCELL(6) character*80 fname character*1 spot(6) data spot /'A','F','B','D','C','E'/ PI=3.14159 CPI=PI/180.0 XG(0)=0 YG(0)=0 ZG(0)=0 det=0.0 write(*,100) 100 format('name of input data file ?',$) read(*,110) fname 110 format(a) open(unit=8,file=fname,status='old') READ(8,*)GPDIST write(6,*)' ' WRITE(6,6010)GPDIST 6010 FORMAT('Origin of Reciprocal Lattice to Gnomonic Plane = ', +f8.2,' mm') READ(8,*)(XG(I),YG(I),I=1,6) WRITE(6,6020) 6020 FORMAT('Coordinates of "Gnomonic Triangle" ...A,B,C are corners:') WRITE(6,6030) 6030 FORMAT('SPOT',2x,'XG',6x,'YG',6x,'ZG',6x,'XL') DO 10 I=1,6 ZG(I)=GPDIST XL(I)=SQRT(XG(I)**2 + YG(I)**2 + ZG(I)**2) WRITE(6,950)SPOT(I),XG(I),YG(I),ZG(I),XL(I) 950 FORMAT(A1,4F8.3) 10 CONTINUE DO 20 I=1,6 DC(I,1)=zG(I)/xl(I) DC(I,2)=xG(I)/XL(I) DC(I,3)=yG(I)/XL(I) 20 CONTINUE C THESE ARE DIRECTION COSINES OF LINES FROM ORIGIN TO EACH POINT C C USE 1ST, 3RD, 5TH TO GET ALPHASTAR, BETASTAR, GAMMASTAR GALPS=ACOS(DC(3,1)*DC(5,1)+DC(3,2)*DC(5,2)+DC(3,3)*DC(5,3)) GBETS=ACOS(DC(5,1)*DC(1,1)+DC(5,2)*DC(1,2)+DC(5,3)*DC(1,3)) GGAMS=ACOS(DC(1,1)*DC(3,1)+DC(1,2)*DC(3,2)+DC(1,3)*DC(3,3)) C CODE NEEDED HERE TO EXTRACT RATIOS AND RELAS, RELBS,RELGS ANG12=ANGLE(XG,YG,ZG,1,0,2) ANG23=ANGLE(XG,YG,ZG,2,0,3) ANG34=ANGLE(XG,YG,ZG,3,0,4) ANG45=ANGLE(XG,YG,ZG,4,0,5) ANG56=ANGLE(XG,YG,ZG,5,0,6) ANG61=ANGLE(XG,YG,ZG,6,0,1) CD WRITE(6,951)ANG12,ANG23,ANG34,ANG45,ANG56,ANG61 951 FORMAT(6F8.3) RABS=SIN(ANG23 )/SIN(ANG12 ) RBCS=SIN(ANG45 )/SIN(ANG34 ) RCAS=SIN(ANG61 )/SIN(ANG56 ) CD WRITE(6,941)RABS,RBCS,RCAS,RABS*RBCS*RCAS 941 FORMAT(' RATIOS AS/BS, BS/CS, CS/AS, PRODUCT(=1?)',4F9.4) C C C SET UP RELATIVE CELL (assuming a* = 0.2500) GAS=0.25 GBS=GAS/RABS GCS=(GBS/RBCS+GAS*RCAS)/2 WRITE(6,*)' ' WRITE(6,6900) WRITE(6,6910) WRITE(6,910)GAS,GBS,GCS,GALPS/CPI,GBETS/CPI,GGAMS/CPI 6900 FORMAT('Reciprocal Cell (assuming a* = 0.2500) ') 6910 FORMAT('a*,b*,c*,alpha*,beta*,gamma*') 910 FORMAT(3f8.4,3f10.3) C COPY INTO RCELL(1-6) AND CONVERT TO REAL CELL RCELL(1)=GAS RCELL(2)=GBS RCELL(3)=GCS RCELL(4)=GALPS RCELL(5)=GBETS RCELL(6)=GGAMS CALL RECIP(RCELL,CELL,1.0) write(6,6032) write(6,6034) write(6,6040)(cell(i),i=1,3),(cell(j)/CPI,j=4,6) 6032 format('Equivalent Real Cell') 6034 format('a,b,c,alpha,beta,gamma') 6040 format(3f8.4,3f10.3) rproduct=rabs*rbcs*rcas write(6,6090)rabs,rbcs,rcas,rproduct 6090 format('ratios a*/b*,b*/c*,c*/a*,(product=1 ?)',3f9.4, +' (',f6.4,')') if (rproduct.lt.0.9.or.rproduct.gt.1.1) then write(6,*) ' ' write(6,6100) 6100 format(' *** warning product not equal to 1 ...', + 'check spot coords ***') write(6,*) ' ' endif C C SET UP MATRIX CAR DO 30 I=1,3 CAR(I,1)=GAS*DC(1,I) CAR(I,2)=GBS*DC(3,I) CAR(I,3)=GCS*DC(5,I) 30 CONTINUE CD WRITE(6,*)' MATRIX C.AR' CD CALL MATPRINT(CAR) C IF CELL REDN NEEDED INSERT HERE C GAS -GGAMS WILL GIVE AS TO GAMS, THEN USE IN SETM C MR = MATRIX TO CONVERT OLD TO NEW INDICES C MATMULT CAR.MR INTO NEW CAR C C ************************************************** CALL SETM(GAS,GBS,GCS,GALPS,GBETS,GGAMS,AR) CD WRITE(6,*)' MATRIX A, STANDARD ORIENTATION FOR THIS CELL' CD CALL MATPRINT(AR) C NOW INVERT AR INTO ARINV CALL INVERT3(AR,ARINV,det) CD WRITE(6,*)' MATRIX ARINV' CD CALL MATPRINT(ARINV) CALL MATMULT(CAR,ARINV,C) CD WRITE(6,*)' MATRIX C' CD CALL MATPRINT(C) CALL SOLVPS(C,PX,PY,PZ) WRITE(6,*) WRITE(6,900) WRITE(6,901) WRITE(6,902)PX,PY,PZ 900 FORMAT('Missetting angles PX,PY,PZ with respect', +' to spindle along c*,') 901 FORMAT('a* nearest axis to x-ray beam in plane defined', +' by the beam axis and rotation axis.') 902 FORMAT('PX,PY,PZ = ',3f12.3) c write(6,*)' ' write(6,6200) write(6,6210) write(6,6220) write(6,6230) write(6,6240) write(6,6250) write(6,6260) write(6,6270) 6200 format(' **********************************************') 6210 format(' ** NOW REFINE THIS CELL IN GENLAUE **') 6220 FORMAT(' ** ONCE YOU HAVE A GOOD RMS THEN REDUCE TO **') 6230 FORMAT(' ** THE CONVENTIONAL CELL USING TRACER **') 6240 FORMAT(' ** ORIENTATION OF TRACER REDUCED CELL **') 6250 FORMAT(' ** OBTAINED FROM PROGRAM "TRANS" **') 6260 FORMAT(' ** HAPPY REFINING ! **') 6270 FORMAT(' **********************************************') 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) CD write (6,18) 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 CD write (6,19) CD write (6,20) endif c 18 format(1x,'determinant = ',f6.3) 19 format(1x,' ') 20 format(1x,' determinant less than 0 elements mulipied by -1') 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) CD WRITE(6,900)I,J,K ,ANG*180/PI 900 FORMAT(10X,' ANGLE',2X,I2,'-',I2,'-',I2,'IS ',F6.1) CD 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 RECIP(CELL,RCELL,CONST) DOUBLE PRECISION CELL(6),RCELL(6) C A TO GAMMA, A AND RADIANS, INVERT TO RECIP CELL PI=3.14159 SA=SIN(CELL(4)) SB=SIN(CELL(5)) SG=SIN(CELL(6)) CA=COS(CELL(4)) CB=COS(CELL(5)) CG=COS(CELL(6)) CD WRITE(6,900)(CELL(I),I=1,3),(CELL(I)*180/PI,I=4,6) 900 FORMAT(' INPUT CELL: ',3F10.5,3F10.3) CASTAR=(CB*CG-CA)/(SB*SG) CBSTAR=(CG*CA-CB)/(SG*SA) CGSTAR=(CA*CB-CG)/(SA*SB) X=0.9999999 IF (CASTAR.GE.1.0000) CASTAR=X IF (CBSTAR.GE.1.0000) CBSTAR=X IF (CGSTAR.GE.1.0000) CGSTAR=X IF (CAS TAR.LE.-1.0000) CASTAR=-X IF (CBSTAR.LE.-1.0000) CBSTAR=-X IF (CGSTAR.LE.-1.0000) CGSTAR=-X SASTAR=SQRT(1-CASTAR*CASTAR) SBSTAR=SQRT(1-CBSTAR*CBSTAR) SGSTAR=SQRT(1-CGSTAR*CGSTAR) RCELL(4)=ACOS(CASTAR) RCELL(5)=ACOS(CBSTAR) RCELL(6)=ACOS(CGSTAR) RCELL(1)=CONST/(CELL(1)*SB*SGSTAR) RCELL(2)=CONST/(CELL(2)*SG*SASTAR) RCELL(3)=CONST/(CELL(3)*SA*SBSTAR) CD WRITE(6,*)' CONSTANT ',CONST CD WRITE(6,910)(RCELL(I),I=1,3),(RCELL(I)*180/PI,I=4,6) 910 FORMAT(' RECIPROCAL:', 3F12.6,3F 8.3) RETURN END 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