program mill2eul c program to convert Miller indices directly c to Bunge Euler angles c ADR i 02, i 05 ! further edits june 07 ! further edits oct 07, to allow for using larger sets as sample symmetry ! e.g. to calculate variants of a twin within a crystal ! minor change to output .WTS file, viii 08 implicit none real hkl(3),uvw(3),t1(3) real a(3,3),d1,d2,d3 integer i,j,k,iopt,ior,kerr,ijk character nomen logical result,result2 real qq1(4),quint(4),qend(4) integer iqbox real dd1,dd2,dd3,hkl_out(3),uvw_out(3),normhkl,normuvw real pi,degrad,twopi,rtmp real scalar integer iq0 real rodvec(3) real tend character header*100 character eulan*25,iper*6 real evm,fk1(9),nstate data fk1 / 1. , 0. , 0. , 0. , 1. , 0. , 0. , 0. , 1. / real quatsymm(4,48),quatsamp(4,48) integer numsymm,numsamp common/a1/ quatsymm,numsymm,quatsamp,numsamp ! CODE:: iper = ' 1 2 3' degrad = 57.29578 PI = 3.14159265 twopi = 2.0*PI nomen = 'B' tend = sqrt(2.) - 1. call textur4 !reads in symmetry operators write(*,*) 'Found ',numsymm,' xtal symms',' and ' write(*,*) numsamp,' sample symms' open(unit=2,file='Bunge.angles.txt',status='replace') write(2,"(' f1 F f2 h k l u v w flag')") open(unit=3,file='Rodrigues.vectors.txt',status='replace') write(3,"(' J K R1 R2 R3')") open(unit=4,file='mill2eul_2_popLA.wts',status='replace') ! changed file extension to .WTS so that wts2pop can read it write(4,"(' written by mill2eul ')") header='Evm F11 F12 F13 F21 ' &//'F22 F23 F31 F32 F33 nstate' write(4,"(a)") header 70 format(10f7.3,i3) evm = 0.0 nstate = 0.0 write(4,70) evm,fk1,nstate 71 format(a25,49x,a6) eulan(1:25)='Bunge phi1 Phi phi2:' ! fix angle type as Bunge (see below) write(4,71) eulan,iper print*,'Enter Miller indices [=0] or Bunge Euler angles [=1]? ' 1 ,'or as a Rodrigues vector [=2]?' read(*,*) iq0 if(iq0.eq.0) then write(*,*) 'To convert Miller indices to Bunge Euler angles' write(*,*) write(*,*) 'Enter (hkl)=sample.3 as 3 numbers' read(*,*) hkl(1),hkl(2),hkl(3) write(2,"( 'Plane HKL = ',3(f8.1))") hkl call vecnorm(hkl) 10 write(*,*) 'Enter [uvw]=sample.1 as 3 numbers' read(*,*) uvw(1),uvw(2),uvw(3) write(2,"( 'Plane UVW = ',3(f8.1))") uvw call vecnorm(uvw) rtmp=scalar(hkl,uvw) if(abs(rtmp).gt.1e-4) then write(*,*) 'Sorry, uvw is not perp. to hkl; try again' goto 10 endif call vecpro2(hkl,uvw,t1) call vecnorm(t1) do 20, i=1,3 a(1,i)=uvw(i) a(2,i)=t1(i) a(3,i)=hkl(i) 20 continue iopt=1 kerr=0 ior=1 ! grain number call euler(a,iopt,nomen,d1,d2,d3,ior,kerr) elseif(iq0.eq.1) then print*,'Enter 3 Euler angles in degrees (Bunge)' read(*,*) d1,d2,d3 elseif(iq0.eq.2) then print*,'Enter 3 Rodrigues parameters ' read(*,*) (rodvec(i),i=1,3) call rod2eul(d1,d2,d3,rodvec) d1 = d1 * degrad d2 = d2 * degrad d3 = d3 * degrad endif write(*,*) 'Euler angles = ',d1,d2,d3 call quatB(d1/degrad,d2/degrad,d3/degrad,qq1) print*,'How many crystal operators to use? [',numsymm,']' read(*,*) numsymm c allow for a smaller number of symms. print*,'How many sample operators to use? [',numsamp,']' read(*,*) numsamp c allow for a smaller number of symms. 23 print*,'Print all angles [0] or just in 90x90x90 [1]?' read(*,*) iqbox if(iqbox.lt.0.or.iqbox.gt.1) goto 23 do 900, j=1,numsymm ! crystal symmetry call postsymm(qq1,j,quint) do 880, k=1,numsamp ! sample symmetry call presamp(quint,k,qend) call q2eulB(d1,d2,d3,qend) if(d1.lt.0.0) d1=d1+twopi if(d2.lt.0.0) d2=d2+twopi if(d3.lt.0.0) d3=d3+twopi d1=d1*degrad d2=d2*degrad d3=d3*degrad c convert to Bunge Euler angles c write(*,*) d1, d2, d3 call testangs(d1,d2,d3,result) result2 = .false. if(d2.le.90..and.d3.le.90.) result2 = .true. dd1=d1/degrad dd2=d2/degrad dd3=d3/degrad hkl_out(1)=sin(dd2)*sin(dd3) hkl_out(2)=sin(dd2)*cos(dd3) hkl_out(3)=cos(dd2) uvw_out(1)=cos(dd1)*cos(dd3)-sin(dd1)*sin(dd3)*cos(dd2) uvw_out(2)=-1.*cos(dd1)*sin(dd3)-sin(dd1)*cos(dd3)*cos(dd2) uvw_out(3)=sin(dd1)*sin(dd2) c write(*,*) 'Plane HKL_OUT = ',hkl_out c write(*,*) 'Direction UVW = ',uvw_out normhkl = amax1(abs(hkl_out(1)),abs(hkl_out(2)) $ ,abs(hkl_out(3))) normuvw = amax1(abs(uvw_out(1)),abs(uvw_out(2)) $ ,abs(uvw_out(3))) do i=1,3 if(abs(hkl_out(i)).lt.normhkl.and.abs(hkl_out(i)).gt.0.01) & normhkl=abs(hkl_out(i)) if(abs(uvw_out(i)).lt.normuvw.and.abs(uvw_out(i)).gt.0.01) & normuvw=abs(uvw_out(i)) enddo do i=1,3 hkl_out(i)=hkl_out(i)/normhkl if(hkl_out(i) .gt. 999.) hkl_out(i) = 999. if(hkl_out(i) .lt. -999.) hkl_out(i) = -999. uvw_out(i)=uvw_out(i)/normuvw if(uvw_out(i) .gt. 999.) uvw_out(i) = 999. if(uvw_out(i) .lt. -999.) uvw_out(i) = -999. enddo call q2rod(qend,rodvec) c write(*,'("j,k, QEND = :",i6,i6,4(1x,f6.3))') c 1 j,k,(qend(i),i=1,4) c write(*,'("Rod vec = :",3(1x,g12.3))') (rodvec(i),i=1,3) ! write(3,'(i6,1x,i6,3(1x,g12.4))') j,k,(rodvec(i),i=1,3) c write(*,*) 'hkl norm, uvw norm = ',normhkl,normuvw c write(*,*) c write(*,*) 'hkl, uvw normalized' c write(*,'("(",3(1x,f6.1),") [",3(1x,f6.1),"]")') hkl_out,uvw_out if(result) then if(abs(rodvec(1)).le.tend .and. $ abs(rodvec(2)).le.tend .and. $ abs(rodvec(3)).le.tend $ ) then write(3,'(i6,1x,i6,3x,3(1x,g10.3),3x $ ,4(1x,g10.3)," all<=90 @ ")') $ j,k,(rodvec(i),i=1,3) ,(qend(ijk),ijk=1,4) else write(3,'(i6,1x,i6,3x,3(1x,g10.3) $ ,3x,4(1x,g10.3)," all<=90 ")') $ j,k,(rodvec(i),i=1,3) ,(qend(ijk),ijk=1,4) endif write(2, $ '(i6,1x,i6,3(1x,f6.2),3(1x,f6.1) $ ,3(1x,f6.1)," all<=90 ")') & j,k,d1, d2, d3, hkl_out,uvw_out write(*, $ '(i6,1x,i6,3(1x,f6.2),3(1x,f6.1) $ ,3(1x,f6.1)," all<=90 ")') & j,k,d1, d2, d3, hkl_out,uvw_out write(4, $ '(3(1x,f6.1) $ ,3(1x,f6.1))') & d1, d2, d3, 1.0 elseif(result2) then if(abs(rodvec(1)).le.tend .and. $ abs(rodvec(2)).le.tend .and. $ abs(rodvec(3)).le.tend $ ) then write(3,'(i6,1x,i6,3x,3(1x,g10.3) $ ,3x,4(1x,g10.3)," F,f2<=90 @ ")') $ j,k,(rodvec(i),i=1,3),(qend(ijk),ijk=1,4) else write(3,'(i6,1x,i6,3x,3(1x,g10.3) $ ,3x,4(1x,g10.3)," F,f2<=90 ")') $ j,k,(rodvec(i),i=1,3) ,(qend(ijk),ijk=1,4) endif write(2, $ '(i6,1x,i6,3(1x,f6.2),3(1x,f6.1) $ ,3(1x,f6.1)," F,f2<=90 ")') & j,k,d1, d2, d3, hkl_out,uvw_out write(*, $ '(i6,1x,i6,3(1x,f6.2),3(1x,f6.1) $ ,3(1x,f6.1)," F,f2<=90 ")') & j,k,d1, d2, d3, hkl_out,uvw_out write(4, $ '(3(1x,f6.1) $ ,3(1x,f6.1))') & d1, d2, d3, 1.0 else if(abs(rodvec(1)).le.tend .and. $ abs(rodvec(2)).le.tend .and. $ abs(rodvec(3)).le.tend) then write(3,'(i6,1x,i6,3x,3(1x,g10.3) $ ,3x,4(1x,g10.3)," ... @ ")') $ j,k,(rodvec(i),i=1,3) ,(qend(ijk),ijk=1,4) else write(3,'(i6,1x,i6,3x,3(1x,g10.3) $ ,3x,4(1x,g10.3)," ...")') $ j,k,(rodvec(i),i=1,3) ,(qend(ijk),ijk=1,4) endif write(2, $ '(i6,1x,i6,3(1x,f6.2), $ 3(1x,f6.1),3(1x,f6.1)," ...")') & j,k,d1, d2, d3, hkl_out,uvw_out write(4, $ '(3(1x,f6.1) $ ,3(1x,f6.1))') & d1, d2, d3, 1.0 if(iqbox.eq.0) then write(*, $ '(i6,1x,i6,3(1x,f6.2),3(1x,f6.1) $ ,3(1x,f6.1)," ...")') & j,k,d1, d2, d3, hkl_out,uvw_out endif endif c write(*,*) 'textangs result: ',result c if(result) then c this should yield 3 points per orientation c thereby enforcing the cubic crystal symmetry within the 90x90x90 space c c DEBUG LINES -> c write(*,*) c write(*,*) 'Angles for grain number ',i c write(*,*) phi1(i),capphi(i),phi2(i) c write(*,*) 'Angles returned from sort ' c endif 880 enddo 900 enddo close(unit=2) close(unit=3) close(unit=4) call exit end c c c ________________________________________ c subroutine vecnorm(vec) real vec(3),rnorm rnorm=0. do 10, i=1,3 rnorm=rnorm+vec(i)**2 10 continue if(rnorm.le.0.0) return rnorm=sqrt(rnorm) do 20, i=1,3 vec(i)=vec(i)/rnorm 20 continue return end c c -------------------- c function scalar(a,b) c return SCALAR PRODUCT of A and B real scalar,a(3),b(3) scalar=0. do 100, i=1,3 scalar=scalar+a(i)*b(i) 100 continue return end c c ____________________________ c subroutine vecpro2(v1,v2,vout) c vector product real v1(3),v2(3),vout(3) vout(3)=v1(1)*v2(2)-v1(2)*v2(1) vout(1)=v1(2)*v2(3)-v1(3)*v2(2) vout(2)=v1(3)*v2(1)-v1(1)*v2(3) return end c c ------------------- c subroutine euler(a,iopt,nomen,d1,d2,d3,ior,kerr) c Last revision 20nov90 UFK c common a(3,3),grvol(1152),epsga(5),ist1,ist2,sqr3,sqrh,ident(3,3) c SPECIAL VERSION WITHOUT COMMON BLOCK real a(3,3),d1,d2,d3,th,sth,cth,sph,cph,sps,cps,ps,ph,dth,dph,dps real pi,rad character nomen save kor PI=3.14159265 rad=57.29578 c c *** this subroutine calculates the euler angles associated with the c transformation matrix a(i,j) if iopt=1 and viceversa if iopt=2 c *** Note that a is sample (rows) in terms of crystal (columns); c -- opposite of standard g (e.g.Bunge) - this is Canova's c *** Note that in this version, the Euler angles are defined symmetrically: c so that interchanging phi and psi means transposing a. c ("Kocks" nomen: defined going from +X to +Y in both COD and SOD) c *** However, other angle conventions are translated, according to c nomen="K" - Kocks (as internally) -- also sometimes "N"... c "R" - Roe (Psi=psi, Phi=180-phi) c "B" - Bunge (phi1=90+psi, Phi=Theta, phi2=90-phi) c any other - Canova (Theta first, phiC=90+phi, omega=90-psi) c *** Note: only in symm. notation does a point with all Euler angles c between 0 and 90 deg appear in the +x,+y quadrant! c If you want to see an individual point in this quadrant and are: c using Roe, the third angle must be between 90 and 180; c Bunge, first ; c Canova, second . c *** Input and output Euler angles d1,d2,d3 in degrees c goto(5,20),iopt 5 if(abs(a(3,3)) .ge. 0.999) goto 10 th=acos(a(3,3)) sth=sin(th) c if((abs(a(2,3)/sth).lt.1e-35).and.(abs(a(1,3)/sth).lt.1e-35)) then ps=pi/4. else ps=atan2(a(2,3)/sth,a(1,3)/sth) endif c if((abs(a(3,2)/sth).lt.1e-35).and.(abs(a(3,1)/sth).lt.1e-35)) then ph=pi/4. else ph=atan2(a(3,2)/sth,a(3,1)/sth) endif ccccc if it bombs out here, both a-components in one arg. are zero c (this should not be possible, but has happened, probably fixed) c ADR: i 01 - above is the fix!! c go to 15 c 10 if((abs(a(1,2)/sth).lt.1e-35).and.(abs(a(1,1)/sth).lt.1e-35)) then ps=pi/4. else ps=0.5*atan2(a(1,2),-a(1,1)) endif c ph=-ps c The above still have the problem that they give too many c equivalents of the same grain in DIOROUT and density file. Therefore: if(kerr.eq.1.and.kor.ne.ior) then print*,'NOTE: grain',ior,' has Theta < 1 deg.: sometimes problems' print* kor=ior endif c 15 dth=th*rad dph=ph*rad dps=ps*rad d1=dps d2=dth if(nomen.eq.'K'.or.nomen.eq.'N') then d3=dph elseif(nomen.eq.'R') then d3=180.-dph elseif(nomen.eq.'B') then d1=dps+90. d3=90.-dph else d1=dth d2=dph+90. d3=90.-dps endif if(d1.ge.360.) d1=d1-360. if(d3.ge.360.) d3=d3-360. if(d1.lt.0.) d1=d1+360. if(d3.lt.0.) d3=d3+360. return c ************************************* 20 dth=d2 dps=d1 if(nomen.eq.'K') then dph=d3 elseif(nomen.eq.'R') then dph=180.-d3 if(dph.lt.0.) dph=dph+360. elseif(nomen.eq.'B') then dps=d1-90. dph=90.-d3 if(dps.lt.0.) dps=dps+360. if(dph.lt.0.) dph=dph+360. else dth=d1 dph=d2-90. if(dph.lt.0.) dph=dph+360. dps=90.-d3 if(dps.lt.0.) dps=dps+360. endif ph=dph/rad th=dth/rad ps=dps/rad sph=sin(ph) cph=cos(ph) sth=sin(th) cth=cos(th) sps=sin(ps) cps=cos(ps) a(1,1)=-sps*sph-cph*cps*cth a(2,1)=cps*sph-cph*sps*cth a(3,1)=cph*sth a(1,2)=sps*cph-sph*cps*cth a(2,2)=-cph*cps-sph*sps*cth a(3,2)=sth*sph a(1,3)=sth*cps a(2,3)=sps*sth a(3,3)=cth return end c c ________________________________________ c subroutine quatB(p1,p,p2,q) c convert Bunge Euler angles (radians) to quaternion; c direct conversion from angles, see Altmann's book c the rotation is a vector transformation (active rotation) c real p1,p,p2,q(4) double precision c1,c,c2,s1,s,s2,d(3,3),tmp(4),sin2,cos2,rtmp,rnorm c c PI=3.14159265 c form cosine, sine of Phi, and sum & diff of phi1, phi2 S=DSIN(0.5d0*P) C=DCOS(0.5d0*P) S1=DSIN(0.5d0*(P1-P2)) C1=DCOS(0.5d0*(P1-P2)) S2=DSIN(0.5d0*(P1+P2)) C2=Dcos(0.5d0*(P1+P2)) c write(*,*) 's1,c1,s,c,s2,c2' c write(*,*) s1,c1 c write(*,*) s,c c write(*,*) s2,c2 q(1)=s*c1 q(2)=s*s1 q(3)=c*s2 q(4)=c*c2 return end c c ------------------ c subroutine testangs(phi1, capphi, phi2,result) c tests to see if the angles fall within the specified range c which for now is just 0quat c real qq(4),p1,p,p2 double precision dp1,dp,dp2,dq1,dq2,dq3,dq4,sum,diff,pi,dtmp c PI=3.14159265d0 dq1=qq(1) dq2=qq(2) dq3=qq(3) dq4=qq(4) if((dabs(dq2).lt.1d-35).and.(dabs(dq1).lt.1d-35)) then diff=pi/4.d0 else diff=datan2(dq2,dq1) endif c ATAN2 is unhappy is both arguments are exactly zero! c if((dabs(dq3).lt.1d-35).and.(dabs(dq4).lt.1d-35)) then sum=pi/4.d0 else sum=datan2(dq3,dq4) endif dp1=(diff+sum) dp2=(sum-diff) dtmp=dsqrt(dq3**2+dq4**2) if(dtmp.gt.1.0d0) dtmp=1.0d0 dp=2.d0*dacos(dtmp) p1=dp1 p=dp p2=dp2 c write(*,*) ' quaternion input= ',qq c write(*,*) 'Bunge angles output= ',p1,p,p2 c write(*,*) ' angles [degrees]= ',180.*p1/pi,180.*p/pi,180.*p2/pi return end c c _____________________________ c c subroutine presamp(qq,lindex,qresult) real qq(4),qresult(4) integer lindex real quatsymm(4,48),quatsamp(4,48) integer numsymm,numsamp common/a1/ quatsymm,numsymm,quatsamp,numsamp c PI=3.14159265 c algorithm for forming resultant quaternion c from applying a symmetry operator, QUATSYMM(n,lindex) c to the first quaternion/rotation c If the symm oper == O, and the input quaternion, QQ, is an active c rotation, then Qresult = (O x QQ) c the "PRE" in the name refers to writing the operator before the c rotation/orientation in vector/tensor notation: Q' = O x Q c or in conventional quaternion notation, Qresult = QQ Ä Qsymm c Thus, for active rotations (standard definition of orientation in c mechanics) this is suitable for applying SAMPLE symmetry if(lindex.gt.numsamp) stop 'error in PRESAMP, lindex>numsamp' if(lindex.lt.1) stop 'error in presamp, lindex<1' qresult(1)=qq(1)*quatsamp(4,lindex) + qq(4)*quatsamp(1,lindex) & - qq(2)*quatsamp(3,lindex) + qq(3)*quatsamp(2,lindex) qresult(2)=qq(2)*quatsamp(4,lindex) + qq(4)*quatsamp(2,lindex) & - qq(3)*quatsamp(1,lindex) + qq(1)*quatsamp(3,lindex) qresult(3)=qq(3)*quatsamp(4,lindex) + qq(4)*quatsamp(3,lindex) & - qq(1)*quatsamp(2,lindex) + qq(2)*quatsamp(1,lindex) qresult(4)=qq(4)*quatsamp(4,lindex) - qq(1)*quatsamp(1,lindex) & - qq(2)*quatsamp(2,lindex) - qq(3)*quatsamp(3,lindex) ! write(*,*) 'qresult ',qresult return end c c _____________________________ c c subroutine postsymm(qq,lindex,qresult) real qq(4),qresult(4) integer lindex c real quatsymm(4,48),quatsamp(4,48) integer numsymm,numsamp common/a1/ quatsymm,numsymm,quatsamp,numsamp c PI=3.14159265 c algorithm for forming resultant quaternion/rotation c from applying a symmetry operator, QUATSYMM(n,lindex) c to the first quaternion/rotation, QQ c If the symm oper == O, and the input quaternion, QQ, is an active c rotation, then Qresult = (QQ x O) c the "POST" in the name refers to writing the operator after the c rotation/orientation in vector/tensor notation: Q' = Q x O c or in conventional quaternion notation, Qresult = Qsymm Ä QQ c Thus, for active rotations (standard definition of orientation in c mechanics) this is suitable for applying CRYSTAL symmetry if(lindex.gt.numsymm) stop 'error in presymm, lindex>numsymm' if(lindex.lt.1) stop 'error in presymm, lindex<1' qresult(1)=quatsymm(4,lindex)*qq(1) + quatsymm(1,lindex)*qq(4) & + quatsymm(3,lindex)*qq(2) - quatsymm(2,lindex)*qq(3) qresult(2)=quatsymm(4,lindex)*qq(2) + quatsymm(2,lindex)*qq(4) & + quatsymm(1,lindex)*qq(3) - quatsymm(3,lindex)*qq(1) qresult(3)=quatsymm(4,lindex)*qq(3) + quatsymm(3,lindex)*qq(4) & + quatsymm(2,lindex)*qq(1) - quatsymm(1,lindex)*qq(2) qresult(4)=quatsymm(4,lindex)*qq(4) - quatsymm(1,lindex)*qq(1) & - quatsymm(2,lindex)*qq(2) - quatsymm(3,lindex)*qq(3) ! write(*,*) 'qresult ',qresult return end !................................. subroutine q2rod(qq,rod) c converts quaternion to Rodrigues vector c ADR, 30 iii 02 implicit none real qq(4),rod(3) if(abs(qq(4)).gt.1e-6) then rod(1)=qq(1)/qq(4) rod(2)=qq(2)/qq(4) rod(3)=qq(3)/qq(4) else rod(1)=qq(1)*999999. rod(2)=qq(2)*999999. rod(3)=qq(3)*999999. endif return end !................................. subroutine rod2eul(d1,d2,d3,rodr) implicit none real d1,d2,d3,rodr(3),sum,diff c converts Rodrigues vector to Bunge Euler angles c d1== phi1, d2== Phi, d3== phi2 sum=atan(rodr(3)) diff=atan(rodr(2)/rodr(1)) d1=sum+diff d2=2.*atan(rodr(1)*cos(sum)/cos(diff)) d3=sum-diff return end