program rexgbs ! compile with ! g77 -ffixed-line-length-none -O3 -o rexgbs rexgbs.f ! or ! g77 -ffortran-bounds-check -g -ffixed-line-length-none -o rexgbs rexgbs.f c calculate grain boundary type for various combinations of c texture components, analyze for CSL etc ! adr iv00, vi 03: last edit 29 Jul 07 ! to include HEX and ORTH ! last edit Jul-Aug 08 ! to include symmetry operation in between operations ! which allows twin chains to be calculated ! e.g. ! hex twins in Ti, Zr ! Compr twin: 0 35.1 0 ,or, 35.1 1 0 0 (Y-conv) like sig 11a ! TT1 : 0 84.8 30 ,or, 84.8 0.866 0.5 0 (Y-conv) like sig 11b ! TT2 : 0 64.6 0 ,or, 64.6 1 0 0 (Y-conv) like sig 7b ! note that, if the additional symmetry operator is selected ! then one must input two rotations ! and the output of Euler angles, quats is restricted to the forward set ! because the order in which the rotations are performed now matters c read orientations from TEXIN c symmetry elements from QUAT.SYMM.CUBIC c CSLs from QUAT.CSL.DATA ! allow for pairs of Rodrigues vectors, x 07 ! and Euler angles in radians, iii 08 ! entry of axis-angle, or pair of axis-angle specifications ! note that Q2ROD now correctly handles angles close to 180 ! Jan 08 - checking rotation of misorientation axis back to SAMPLE implicit none real pi integer i ,j ,k ,ijk ,jkl ,i2 integer index1 ,index2 ,index1a ,index2a real thetmin ,brandon ,crit ,theta ,al_min ,tmp real qr1(4),qr2(4),qtwo(4,2),qresult(4),axis(3) real angle(100,100),gbaxis(3,100,100),rnear(100,100) integer lsig(100,100) real p1(2) ,pp(2) ,p2(2) integer ac_min(3) real qresulta(100,100,4) integer ij,jk integer iq,iq2 real angle_input,misaxis(3),tmpvec(3),rfvector(3),tmpout(3) real rad , rad2 real rfvector2(3) character afile*100 real gbnorma(3),gbdirna(3),gbnormb(3),gbdirnb(3),Bunge(3) real matrixA(3,3),matrixB(3,3),matrix(3,3) real quatsymm(4,48) ,quatsamp(4,4) integer numsymm ,numsamp common/a1/ quatsymm,numsymm,quatsamp,numsamp real quatcsl(4,50) integer numcsl ,nsig(50) common/a2/quatcsl ,numcsl ,nsig c common/a3/aquat(4,100),ngrain,phi1(100),capphi(100),phi2(100) real aquat(4,100000) & ,phi1(100000),capphi(100000),phi2(100000),grwt(100000) integer ngrain common/a3/aquat ,phi1 ,capphi ,phi2 ,grwt ,ngrain integer iq_symm character inline*120 integer iq_extra , iex , ulimit ! for insertion of extra symm. oper. real qtmp(4) , qtmp1(4) , qtmp2(4) , d1 , d2 , d3 , rnorm real qdisq1(4) , qdisq2(4) real samp1(3) , samp2(3) ! for testing rotations of the misorient axis real qgrain1(4) , qgrain2(4) , qgrain3(4) , qgrain4(4) integer combine ! CODE: PI = 3.14159265 rad = 180. / pi rad2 = 360. / pi ! CALL GetArg(1, afile) write(*,*) 'Calculates a set of misorientations etc, from a list' write(*,*) ' of orientations in a file (TEXIN).' print*,'Be sure that you include ' write(*,*) ' all the sample-symmetry related variants for each' write(*,*) ' component or you will not get the full picture!' c write(*,*) 'Choice: read orientations from TEXIN [=0]' c write(*,*) 'or, enter orientations directly [=1]?' c read(*,*) iq iq = 0 CALL GetArg(1, afile) if(afile.eq.'') then print *,' No input file given as argument: input from terminal' print*,' Sample usage:' print*,'[/code/texture/]rexgbs data 0 ' iq = 1 endif iq_symm = 0 CALL GetArg(2, inline) if(inline.eq.'') then print *,' No symmetry as argument: input from terminal' print*,'Enter xtal symmetry: 0=cubic 1=Hex 2=Orth' read*,iq_symm else ! print*,'second argument = ',inline read(inline,"(i1)") iq_symm ! print*,'iq_symm = ',iq_symm if(iq_symm.ge.0 .and. iq_symm.lt.3) then print*,'Crystal symmetry: 0=cubic 1=Hex 2=Orth' print*,'Symmetry choice = ',iq_symm else stop 'invalid symmetry choice' endif endif call textur4(iq_symm) ! read in the symmetry operators print*,'CAUTION!!!' print*,'If you are combining 2 successive transformations' print*,'as in, e.g., twin chains, you must be v careful' print*,'about the order.' print*,'Remember that the second rotation is negated' print*,'whenever a pair of orientations is combined.' print*,'Accordingly, you should give the negative (inverse)' print*,'of the 2nd rotation to get the desired effect.' print*,'I.e. enter gA, then -gB.' print*,'To reverse the order, enter gB, then -gA' print*,'________________________________________' print* if(iq.lt.0.or.iq.gt.1) then stop 'invalid choice' elseif(iq.eq.0) then call textur2(afile) elseif(iq.eq.1) then call textur6 ngrain=2 ! limit to 2 grains c$$$ write(*,*) 'Choose also: enter Euler angles [=0]' c$$$ write(*,*) ' or an axis-angle misorientation [=1]?' c$$$ write(*,*) ' or as UVW-HKL-A + UVW-HKL-B [=2]?' c$$$ read(*,*) iq2 CALL GetArg(3, inline) if(inline.eq.'') then write(*,*) 'Choose a parameterization: ' write(*,*) 'A single entry defines the misorientation ' write(*,*) 'A paired entry defines two orientations:' write(*,*) 'Pair of Euler angles in radians [ = 0 ]' write(*,*) ' or a pair of Euler angles in degrees [ = 1 ]' print*,' or a single axis-angle misorientation [ = 2 ]?' write(*,*) ' or as UVW-HKL-A + UVW-HKL-B [ = 3 ]?' print*,' or as a single Rodrigues-Frank vector [ = 4 ]?' print*,' or as a PAIR of Rodrigues-Frank vectors [ = 5 ]?' print*,' or as a PAIR of axis-angle orientations [ = 6 ]?' print*,' or as a PAIR of quaternions [ = 7 ]?' print*,'BEWARE: remember that a pair with the same ' print*,'orientation is treated as zero misorientation!' print* read(*,*) iq2 else read(inline,"(i)") i2 if( iq2 .ge. 0 .and. iq2 .le. 7 ) then print* , 'input parameter selected = ' , iq2 else stop 'invalid input parameter choice' endif endif if(iq2.eq.0 .or. iq2.eq.1 .or. iq2.ge.5 .or. iq2.le.7 ) then print*,'Combine the two into one misor [ yes = 1 ]?' read*,combine end if if( iq2 .eq. 0 ) then write(*,*) 'Enter phi1A PhiA phi2A phi1B PhiB phi2B [rad]' read(*,*) phi1(1),capphi(1),phi2(1) & ,phi1(2),capphi(2),phi2(2) call quatB(phi1(1),capphi(1),phi2(1),qr1) call quatB(phi1(2),capphi(2),phi2(2),qr2) endif if( iq2 .eq. 1 ) then write(*,*) 'Enter phi1A PhiA phi2A phi1B PhiB phi2B [degr]' read(*,*) phi1(1),capphi(1),phi2(1) & ,phi1(2),capphi(2),phi2(2) call quatB(phi1(1)/rad,capphi(1)/rad,phi2(1)/rad,qr1) call quatB(phi1(2)/rad,capphi(2)/rad,phi2(2)/rad,qr2) endif if( iq2 .eq. 0 .or. iq2 .eq. 1 ) then do ijk=1,4 aquat(ijk,1)=qr1(ijk) aquat(ijk,2)=qr2(ijk) enddo ! endif elseif( iq2 .eq. 2 ) then write(*,*) 'Enter Angle(degr) axis1, axis2, axis3' read(*,*) angle_input,misaxis call donorm(misaxis) ! normalize phi1(1)=0. capphi(1)=0. phi2(1)=0. aquat(1,1)=0. aquat(2,1)=0. aquat(3,1)=0. aquat(4,1)=1. aquat(1,2)=misaxis(1)*sin(angle_input*pi/360.) aquat(2,2)=misaxis(2)*sin(angle_input*pi/360.) aquat(3,2)=misaxis(3)*sin(angle_input*pi/360.) aquat(4,2)=cos(angle_input*pi/360.) do ijk=1,4 qr2(ijk)=aquat(ijk,2) enddo write(*,"('quaternion = ',4f8.3,2i4)") qr2 call q2eulB(phi1(2),capphi(2),phi2(2),qr2) phi1(2)=phi1(2)*180./pi capphi(2)=capphi(2)*180./pi phi2(2)=phi2(2)*180./pi if(phi1(2).gt.360.) phi1(2)=phi1(2)-360. if(phi1(2).lt.0.) phi1(2)=phi2(2)+360. if(phi2(2).gt.360.) phi2(2)=phi1(2)-360. if(phi2(2).lt.0.) phi2(2)=phi2(2)+360. elseif( iq2 .eq. 3 ) then print*,'Enter H K L U V W for xtal A' read*,(gbnorma(i),i=1,3),(gbdirna(i),i=1,3) call uvwhkl2quat(gbnorma,gbdirna,qr1,Bunge,matrixA) do i = 1,4 aquat(i,1) = qr1(i) enddo phi1(1) = Bunge(1) capphi(1) = Bunge(2) phi2(1) = Bunge(3) print* print*,'Orientation matrix for grain A' do i = 1,3 print"('[ ',3(1x,f8.3),' ]')",(matrixA(i,j),j=1,3) enddo print*,'Enter H K L U V W for xtal B' read*,(gbnormb(i),i=1,3),(gbdirnb(i),i=1,3) c here we re-use code from rexgbs.f (aug 06) call uvwhkl2quat(gbnormb,gbdirnb,qr2,Bunge,matrixB) phi1(2) = Bunge(1) capphi(2) = Bunge(2) phi2(2) = Bunge(3) do i = 1,4 aquat(i,2) = qr2(i) enddo print*,'Orientation matrix for grain B' do i = 1,3 print"('[ ',3(1x,f8.3),' ]')",(matrixB(i,j),j=1,3) enddo print* print*,'Misorientation matrix from A to B' do i = 1,3 do j = 1,3 matrix(i,j) = 0. do k = 1,3 matrix(i,j) = 1 matrix(i,j) + matrixB(i,k)*matrixA(j,k) ! note inversion for matrix A enddo enddo print"('[ ',3(1x,f8.3),' ]')",(matrix(i,j),j=1,3) enddo tmp = (-1. + matrix(1,1)+matrix(2,2)+matrix(3,3)) * 0.5 if(tmp.lt.-1.) tmp = -1. if(tmp.gt.1.) tmp = 1. print*,'Angle = ',acos(tmp) * rad ! RF vector entry (single) elseif( iq2 .eq. 4 ) then write(*,*) 'Enter R1 R2 R3' read(*,*) rfvector call rod2q(rfvector,qr2) phi1(1) = 0. capphi(1) = 0. phi2(1) = 0. aquat(1,1) = 0. aquat(2,1) = 0. aquat(3,1) = 0. aquat(4,1) = 1. aquat(1,2) = qr2(1) aquat(2,2) = qr2(2) aquat(3,2) = qr2(3) aquat(4,2) = qr2(4) write(*,"('quaternion = ',4f8.3,2i4)") qr2 call q2eulB(phi1(2),capphi(2),phi2(2),qr2) phi1(2)=phi1(2)*180./pi capphi(2)=capphi(2)*180./pi phi2(2)=phi2(2)*180./pi if(phi1(2).gt.360.) phi1(2)=phi1(2)-360. if(phi1(2).lt.0.) phi1(2)=phi2(2)+360. if(phi2(2).gt.360.) phi2(2)=phi1(2)-360. if(phi2(2).lt.0.) phi2(2)=phi2(2)+360. ! pair of RF vectors elseif( iq2.eq. 5 ) then write(*,*) 'Enter R1 R2 R3 (for A), R1 R2 R3 (for B)' read(*,*) rfvector,rfvector2 call rod2q(rfvector,qr2) aquat(1,1) = qr2(1) aquat(2,1) = qr2(2) aquat(3,1) = qr2(3) aquat(4,1) = qr2(4) write(*,"('quaternion for A = ',4f8.3,2i4)") qr2 call q2eulB(phi1(1),capphi(1),phi2(1),qr2) phi1(1)=phi1(1)*180./pi capphi(1)=capphi(1)*180./pi phi2(1)=phi2(1)*180./pi if(phi1(1).gt.360.) phi1(1)=phi1(1)-360. if(phi1(1).lt.0.) phi1(1)=phi2(1)+360. if(phi2(1).gt.360.) phi2(1)=phi1(1)-360. if(phi2(1).lt.0.) phi2(1)=phi2(1)+360. call rod2q(rfvector2,qr2) aquat(1,2) = qr2(1) aquat(2,2) = qr2(2) aquat(3,2) = qr2(3) aquat(4,2) = qr2(4) write(*,"('quaternion for B = ',4f8.3,2i4)") qr2 call q2eulB(phi1(2),capphi(2),phi2(2),qr2) phi1(2) = phi1(2)*180./pi capphi(2) = capphi(2)*180./pi phi2(2) = phi2(2)*180./pi if(phi1(2).gt.360.) phi1(2) = phi1(2)-360. if(phi1(2).lt.0.) phi1(2) = phi2(2)+360. if(phi2(2).gt.360.) phi2(2) = phi1(2)-360. if(phi2(2).lt.0.) phi2(2) = phi2(2)+360. ! pair of axis-angle specs. elseif( iq2 .eq. 6 ) then write(*,*) 'Enter 1st Angle(degr) axis1, axis2, axis3' read(*,*) angle_input,misaxis call donorm(misaxis) ! normalize aquat(1,1) = misaxis(1)*sin(angle_input / rad2) aquat(2,1) = misaxis(2)*sin(angle_input / rad2) aquat(3,1) = misaxis(3)*sin(angle_input / rad2) aquat(4,1) = cos(angle_input / rad2 ) do ijk = 1 , 4 qr1(ijk) = aquat(ijk,1) enddo write(*,"('1st quaternion = ',4f8.3,2i4)") qr1 call q2eulB(phi1(1),capphi(1),phi2(1),qr1) phi1(1) = phi1(1) * rad capphi(1) = capphi(1) * rad phi2(1) = phi2(1) * rad if(phi1(1).gt.360.) phi1(2)=phi1(2)-360. if(phi1(1).lt.0.) phi1(2)=phi2(2)+360. if(phi2(1).gt.360.) phi2(2)=phi1(2)-360. if(phi2(1).lt.0.) phi2(2)=phi2(2)+360. write(*,*) 'Enter 2nd Angle(degr) axis1, axis2, axis3' read(*,*) angle_input,misaxis call donorm(misaxis) ! normalize aquat(1,2) = misaxis(1)*sin( angle_input / rad2 ) aquat(2,2) = misaxis(2)*sin( angle_input / rad2 ) aquat(3,2) = misaxis(3)*sin( angle_input / rad2 ) aquat(4,2) = cos( angle_input / rad2 ) do ijk=1,4 qr2(ijk) = aquat(ijk,2) enddo write(*,"('2nd quaternion = ',4f8.3,2i4)") qr2 call q2eulB(phi1(2),capphi(2),phi2(2),qr2) phi1(2)=phi1(2)* rad capphi(2)=capphi(2)* rad phi2(2)=phi2(2)* rad if(phi1(2).gt.360.) phi1(2)=phi1(2)-360. if(phi1(2).lt.0.) phi1(2)=phi2(2)+360. if(phi2(2).gt.360.) phi2(2)=phi1(2)-360. if(phi2(2).lt.0.) phi2(2)=phi2(2)+360. ! pair of quaternions elseif( iq2.eq. 7 ) then print*,'Enter 2 quaternions, with cos(t/2) as q4;' print*,'Enter q1 q2 q3 q4 (for A), q1 q2 q3 q4 (for B):' read(*,*) qr1 , qr2 aquat(1,1) = qr1(1) aquat(2,1) = qr1(2) aquat(3,1) = qr1(3) aquat(4,1) = qr1(4) write(*,"('check: quaternion for A = ',4f8.3,2i4)") qr1 call q2eulB(phi1(1),capphi(1),phi2(1),qr2) phi1(1)=phi1(1)*180./pi capphi(1)=capphi(1)*180./pi phi2(1)=phi2(1)*180./pi if(phi1(1).gt.360.) phi1(1)=phi1(1)-360. if(phi1(1).lt.0.) phi1(1)=phi2(1)+360. if(phi2(1).gt.360.) phi2(1)=phi1(1)-360. if(phi2(1).lt.0.) phi2(1)=phi2(1)+360. aquat(1,2) = qr2(1) aquat(2,2) = qr2(2) aquat(3,2) = qr2(3) aquat(4,2) = qr2(4) write(*,"('check: quaternion for B = ',4f8.3,2i4)") qr2 call q2eulB(phi1(2),capphi(2),phi2(2),qr2) phi1(2) = phi1(2)*180./pi capphi(2) = capphi(2)*180./pi phi2(2) = phi2(2)*180./pi if(phi1(2).gt.360.) phi1(2) = phi1(2)-360. if(phi1(2).lt.0.) phi1(2) = phi2(2)+360. if(phi2(2).gt.360.) phi2(2) = phi1(2)-360. if(phi2(2).lt.0.) phi2(2) = phi2(2)+360. elseif( iq2 .lt. 0 .or. iq2 .gt. 7 ) then stop 'invalid choice' endif ! if(iq2.eq.0 ... choice of input c everything should now be set up to do the calculation endif ! if(iq.lt.0 ... c get the orientations, symms etc. CALL GetArg(4, inline) if(inline.eq.'') then print*,'Insert additional symmetry operation [ yes=1 ]' read*,iq_extra else read(inline,"(i)") iq_extra if( iq_extra .eq. 1 ) then if ( iq2 .eq. 2 .or. iq2 .eq. 4 ) then print*,'Sorry, must enter TWO rotations' stop endif print*,'will insert extra symmetry operator ' endif endif !.......... COMBINE the two together ........... if(combine.eq.1) then do ijk = 1 , 3 qtwo( ijk , 1 ) = aquat( ijk , 1 ) qtwo( ijk , 2 ) = -1. * aquat( ijk , 2 ) ! reversal required to compensate for effect of comquat end do qtwo( 4 , 1 ) = aquat( 4 , 1 ) qtwo( 4 , 2 ) = aquat( 4 , 2 ) call comquat(qtwo,qr1) axis(1) = qr1(1) axis(2) = qr1(2) axis(3) = qr1(3) call vecnorm(axis) print*,'output from combination' print"('Quat from COMQUAT = ',4(2x,f8.3))" 1 ,(qr1( jkl ) , jkl = 1 , 4 ) print"('Misor. Axis= ',3(2x,f8.3))",(axis(ijk),ijk=1,3) do ijk = 1 , 4 aquat( ijk , 2 ) = qr1(ijk) aquat( ijk , 1 ) = 0. end do aquat( 4 , 1 ) = 1. iq2 = 2 phi1(1) = 0. capphi(1) = 0. phi2(1) = 0. call q2eulB (phi1(2) , capphi(2) , phi2(2) , qr1 ) phi1(2)=phi1(2)*180./pi capphi(2)=capphi(2)*180./pi phi2(2)=phi2(2)*180./pi if(phi1(2).gt.360.) phi1(2)=phi1(2)-360. if(phi1(2).lt.0.) phi1(2)=phi2(2)+360. if(phi2(2).gt.360.) phi2(2)=phi1(2)-360. if(phi2(2).lt.0.) phi2(2)=phi2(2)+360. end if ! ++++++ end of input +++++++++ ! ++++++++++++++++++++++++++++++ open(unit=1,file='rexgbs_output.txt',status='replace') print*,'output file rexgbs_output.txt contains list of misors.' write(1,*) ngrain,' = number of components' open(unit=2,file='rexgbs_misors.mdf',status='replace') print*, 'plot rexgbs_misors.mdf with RF_misor; e.g. ...' print*,'/code/RF_misor rexgbs_misors.mdf 1.5 0 0 999999999 0 0' print*,' for a cubic case, or this for HCP:' print*,'/code/RF_misor rexgbs_misors.mdf 1.5 0 1 999999999 0 0' write(2,*) ' i j q1 q2 q3 q4 E1 E2 E3 / for plotting with RFdist ' open(unit=3,file='rexgbs_Eulers.wts',status='replace') print*, 'process rexgbs_Eulers.wts with [/code/texture/]wts2pop' print*,'and then use Draw_stereograms to plot the discrete PFs' write(3,*) 'rexgbs E1 E2 E3 1.0 i j q1 q2 q3 q4 ' write(3,*) 'Evm F11 F12 F13 F21 F22 D23 F31 F32 F33' write(3,*) ' 0. 1. 0. 0. 0. 1. 0. 0. 0. 1. ' write(3,"('Bunge:Psi Theta phi ,gr.wt.')") open(unit=4,file='rexgbs_misors_forward.mdf',status='replace') print*, 'plot rexgbs_misors_forward.mdf with RF_misor' print*,'this contains misor-A followed by misor-B only' write(4,*) ' i j q1 q2 q3 q4 E1 E2 E3 / for plotting with RFdist ' ulimit = 1 if ( iq_extra .eq. 1 ) ulimit = numsymm ! used below in additional loop do 1000, i=1,(ngrain-1) p1(1)=phi1(i) pp(1)=capphi(i) p2(1)=phi2(i) write(*,*) write(1,*) write(1,*) '..................................................' write(1,*) 'Grain number ',i write(*,*) '1st Grain: Euler angles: ',p1(1),pp(1),p2(1) write(1,"('1st Grain: Euler angles: ',3(1x,f8.3))") $ p1(1),pp(1),p2(1) write(1,"(4f8.3,' = 1st gr quat')") (aquat(ijk,i),ijk=1,4) ! use this section when applying extra symmetry to 2nd g c$$$ do ijk=1,4 c$$$ qtwo( ijk , 1 ) = aquat( ijk , i ) c$$$ enddo do ijk = 1 , 4 qtmp1(ijk)= aquat( ijk , i ) enddo do 997 iex = 1 , ulimit if ( ulimit .gt. 1 ) then print*,'Inner symmetry element no. ',iex write(1,*) 'Inner symmetry element number ' , iex endif ! this is the extra, inserted symmetry element call presymm( qtmp1 , iex , qtmp2 ) ! call postsymm( qtmp1 , iex , qtmp2 ) ! write(1,*) ' 1st 2nd Angle Axis123 Quat. Sigma Near' do ijk = 1 , 4 qtwo( ijk , 1 ) = qtmp2(ijk) qdisq1(ijk) = aquat( ijk , i ) ! use w disquat2 enddo do 990, j = (i+1) , ngrain write(1,*) write(1,*) ' against grain number ' , j write(1,"(4f8.3,' = 2nd gr quat')") $ (aquat(ijk,j),ijk=1,4) c$$$ do ijk = 1 , 4 c$$$ qtmp1(ijk) = aquat( ijk , j ) c$$$ enddo c$$$ do 997 iex = 1 , ulimit c$$$ if ( ulimit .gt. 1 ) then c$$$ print*,'Inner symmetry element no. ',iex c$$$ endif c$$$! this is the extra, inserted symmetry element c$$$ call presymm( qtmp1 , iex , qtmp2 ) c$$$! call postsymm( qtmp1 , iex , qtmp2 ) c$$$ rnorm = sqrt(qtmp2(1)**2 + qtmp2(2)**2 c$$$ $ + qtmp2(3)**2 + qtmp2(4)**2) c$$$ if( rnorm .lt. 1.e-5 ) c$$$ $ stop 'Q2ROD: too small a quaternion!' c$$$ if( abs( rnorm - 1.0 ) .gt. 1.e-3 ) then c$$$ write(*,*) 'qtmp2 required normalization' c$$$ do ijk = 1 , 4 c$$$ qtmp2(ijk) = qtmp2(ijk) / rnorm c$$$ enddo c$$$ endif do ijk = 1 , 4 qtwo( ijk , 2 ) = aquat( ijk , j ) ! qtwo( ijk , 2 ) = qtmp2( ijk ) qdisq2(ijk) = aquat( ijk , j ) enddo call comquat(qtwo,qr1) axis(1) = qr1(1) axis(2) = qr1(2) axis(3) = qr1(3) call vecnorm(axis) print"('Quat from COMQUAT = ',4(2x,f8.3))" 1 ,(qr1( jkl ) , jkl = 1 , 4 ) print"('Misor. Axis= ',3(2x,f8.3))",(axis(ijk),ijk=1,3) do ijk = 1 , 4 qgrain1(ijk)= aquat( ijk , i ) qgrain2(ijk)= aquat( ijk , j ) enddo call qrvec ( qgrain1 , axis , samp1 ) print"('Qgrain1 = ',4(2x,f8.3))" 1 ,( qgrain1( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 1 = ',3(2x,f8.3))" $ , ( samp1(ijk) , ijk = 1 , 3 ) call qrvec ( qgrain2 , axis , samp2 ) print"('Qgrain2 = ',4(2x,f8.3))" 1 ,( qgrain2( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 2 = ',3(2x,f8.3))" $ , ( samp2(ijk) , ijk = 1 , 3 ) if(iq_symm.eq.0) then ! only for cubic c first use Dierk's subroutine MISORI p1(1)=phi1(i) pp(1)=capphi(i) p2(1)=phi2(i) p1(2)=phi1(j) pp(2)=capphi(j) p2(2)=phi2(j) write(*,*) '2nd Grain: Euler angles: ' $ , p1(2) , pp(2) , p2(2) write(1,*) '2nd Grain: Euler angles: ' $ , p1(2) , pp(2) , p2(2) call misori( p1 , pp , p2 , tmp , 1 , ac_min ) write(1,*) 'MISORI: angle= ',tmp $ ,' , axis= ',ac_min write(*,*) 'MISORI: angle= ',tmp $ ,' , axis= ',ac_min p1(1)=phi1(i) pp(1)=capphi(i) p2(1)=phi2(i) p1(2)=phi1(j) pp(2)=capphi(j) p2(2)=phi2(j) call misorinv(p1,pp,p2,tmp,1,ac_min) c write(1,*) 'MISORInv: angle= ',tmp,' axis= ',ac_min write(*,*) 'MISORInv: angle= ',tmp,' axis= ',ac_min call misquat(qtwo,tmp) write(1,*) 'MISQUAT: angle= ',tmp write(*,*) 'MISQUAT: angle= ',tmp ! & ,', IDs,angle,axis,quat,sig,nr' c call misquatinv(qtwo,tmp) c write(1,*) 'MISQUATINV: angle= ',tmp,' ,IDs,angle,axis,quat,sig,nr' c write(*,*) 'MISQUATINV: angle= ',tmp c endif ! cubic only if(iq_symm.eq.0) then call disquat( qr1 , index1 , index2 , qresult $ , index1a , index2a ) elseif(iq_symm.eq.1) then call disquat_hex & (qr1,index1,index2,qresult,index1a,index2a) elseif(iq_symm.eq.2) then call disquat_orth & (qr1,index1,index2,qresult,index1a,index2a) endif print* write(*,"('from DISQUAT, qresult = ',4(1x,f8.3)/4i4)") 1 qresult, index1 , index2 , index1a , index2a if(qresult(4).gt.1.0) qresult(4)=1.0 if(qresult(4).lt.-1.0) qresult(4)=-1.0 angle(i,j)=360./pi*acos(qresult(4)) write(*,"( 'from DISQUAT, Angle = ',f8.3)") angle(i,j) qresulta(i,j,1)=qresult(1) qresulta(i,j,2)=qresult(2) qresulta(i,j,3)=qresult(3) qresulta(i,j,4)=qresult(4) axis(1)=qresult(1) axis(2)=qresult(2) axis(3)=qresult(3) call vecnorm(axis) gbaxis(1,i,j)=axis(1) gbaxis(2,i,j)=axis(2) gbaxis(3,i,j)=axis(3) print"('Quat from DISQUAT = ',4(2x,f8.3))" 1 ,(qresulta(i,j,jkl),jkl=1,4) print"(' Axis= ',3(2x,f8.3))",(axis(ijk),ijk=1,3) do ijk = 1 , 4 qgrain1(ijk)= aquat( ijk , i ) enddo call qrvec ( qgrain1 , axis , samp1 ) print"('Qgrain1 = ',4(2x,f8.3))" 1 ,( qgrain1( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 1 = ',3(2x,f8.3))" $ , ( samp1(ijk) , ijk = 1 , 3 ) do ijk = 1 , 4 qgrain2(ijk)= aquat( ijk , j ) enddo call qrvec ( qgrain2 , axis , samp2 ) print"('Qgrain2 = ',4(2x,f8.3))" 1 ,( qgrain2( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 2 = ',3(2x,f8.3))" $ , ( samp2(ijk) , ijk = 1 , 3 ) call presymm( qgrain1 , index1 , qgrain3 ) call qrvec ( qgrain3 , axis , samp1 ) print"('Qgrain3 = ',4(2x,f8.3))" 1 ,( qgrain3( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 1, after symm. = ' $ ,3(2x,f8.3))" $ , ( samp1(ijk) , ijk = 1 , 3 ) call presymm( qgrain2 , index2 , qgrain4 ) call qrvec ( qgrain4 , axis , samp1 ) print"('Qgrain4 = ',4(2x,f8.3))" 1 ,( qgrain4( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 2, after symm. = ' $ ,3(2x,f8.3))" $ , ( samp2(ijk) , ijk = 1 , 3 ) do ijk=1,4 qr2(ijk) = qresulta(i,j,ijk) enddo call q2rod(rfvector,qr2) print"('RF_vec = ',3(2x,f8.3))" $ , ( rfvector(ijk) , ijk = 1 , 3 ) write(1,"('RF_vec = ',3(2x,f8.3))") $ ( rfvector(ijk) , ijk = 1 , 3 ) ! if you want to check the result (for cubics) uncomment this block if(iq_symm.eq.0) then call disquat2(qdisq1,qdisq2,index1 $ ,index2,qresult,index1a,index2a) if(qresult(4).gt.1.0) qresult(4)=1.0 if(qresult(4).lt.-1.0) qresult(4)=-1.0 angle(i,j)=360./pi*acos(qresult(4)) write(*,"( 'from DISQUAT2, Angle = ',f8.3)") $ angle(i,j) write(*,"('from DISQUAT2, qresult = ',4(1x,f8.3),4i4)") 1 qresult,index1,index2,index1a,index2a endif if(iq_symm.eq.0) then lsig(i,j)=0 rnear(i,j)=99. if(angle(i,j).lt.15.0) then lsig(i,j)=1 rnear(i,j)=angle(i,j) c call it sigma=1 if low angle bdy endif thetmin=66.0 if(lsig(i,j).eq.0) then do 150, ijk=1,numcsl brandon=15.0/sqrt(float(nsig(ijk))) call cslquat(qresult,ijk,theta,qr2) crit=theta/brandon ! write(*,"('sig,qr2,theta = ',i4,8f8.3)") nsig(ijk),qr2,theta,crit if(crit.lt.thetmin) then ! write(*,"('sig,qr2,theta = ',i4,8f8.3)") nsig(ijk),qr2,theta thetmin = crit lsig(i,j) = nsig(ijk) rnear(i,j) = crit endif 150 continue endif print* ! print"('indices: ',2i4,' angle= ',f8.3)",i,j,angle(i,j) print"('nearest Sigma, distance = 'i4,f8.3)" 1 ,lsig(i,j),rnear(i,j) ! write(1,"(' IDs, angle, axis, quat, sigma, nearness')") write(1,"('symm op indices: ',2(1x,i4))") i,j write(1,"('angle = ',f8.3)") angle(i,j) write(1,"('angle = ',3(1x,f8.3))") (axis(ijk),ijk=1,3) write(1,"('quat = ',4(1x,f8.3))") $ (qresulta(i,j,jkl),jkl=1,4) write(1,"('sig, near = ',i5,(1x,f8.3))") $ lsig(i,j),rnear(i,j) else ! non cubic, leave out Sigma write(1,*) " i j Angle(degr) Axis(3) Quaternion(4)" write(1,890) i,j,angle(i,j),(axis(ijk),ijk=1,3) & ,(qresulta(i,j,jkl),jkl=1,4) endif ! cubic? do jkl = 1 , 4 qtmp(jkl) = qresulta(i,j,jkl) enddo call q2eulB( d1 , d2 , d2 , qtmp ) write(2,"(2i6,8(1x,f10.4))") i,j $ ,(qtmp(jkl),jkl=1,4),d1*rad,d2*rad,d3*rad ! write out result for RFdist input goto 2100 ! jump around reverse order calculation if above not commented c ***************** now for the reverse order! call comquatinv(qtwo,qr1) call disquat(qr1,index1,index2,qresult,index1a,index2a) write(*,"('qresult= ',4f8.3,2i4)") qresult,index1a,index2a if(qresult(4).gt.1.0) qresult(4)=1.0 if(qresult(4).lt.-1.0) qresult(4)=-1.0 angle(i,j)=360./pi*acos(qresult(4)) write(*,"('Angle= ',f8.3)") angle(i,j) lsig(i,j)=0 rnear(i,j)=99. qresulta(i,j,1)=qresult(1) qresulta(i,j,2)=qresult(2) qresulta(i,j,3)=qresult(3) qresulta(i,j,4)=qresult(4) if(angle(i,j).lt.15.0) then lsig(i,j)=1 rnear(i,j)=angle(i,j) c call it sigma=1 if low angle bdy endif axis(1)=qresult(1) axis(2)=qresult(2) axis(3)=qresult(3) call vecnorm(axis) gbaxis(1,i,j)=axis(1) gbaxis(2,i,j)=axis(2) gbaxis(3,i,j)=axis(3) thetmin=66.0 if(lsig(i,j).eq.0) then do 250, ijk=1,numcsl brandon=15.0/sqrt(float(nsig(ijk))) call cslquat(qresult,ijk,theta,qr2) crit=theta/brandon if(crit.lt.thetmin) then c write(*,"('qr2,theta= ',8f8.3)") qr2,theta thetmin=crit lsig(i,j)=nsig(ijk) rnear(i,j)=crit endif 250 continue endif write(*,890) i,j,angle(i,j),(axis(ijk),ijk=1,3) & ,(qresulta(i,j,jkl),jkl=1,4),lsig(i,j),rnear(i,j) write(1,890) i,j,angle(i,j),(axis(ijk),ijk=1,3) & ,(qresulta(i,j,jkl),jkl=1,4),lsig(i,j),rnear(i,j) 890 format(2i4,8f8.3,i4,f8.3) 2100 continue ! jumps around reverse order calculation to here ! 997 enddo ! apply extra symmetry to 2nd quat 990 enddo 997 enddo ! apply extra symmetry to 1st quat 1000 enddo write(1,*) write(1,*) 'Next section intended as input for RFplot==rfp' write(1,*) 'which you can find in /dihedral/.' write(1,*) write(*,*) 'Last section of rexgbs.out intended as ' write(*,*) 'input for RFplot==rfp, ' print*,' which you can find in /dihedral/.' print*,'on the screen, R1,R2,R3; in the file R3,R2,R1 ' print* print*,'Output also available as rexgbs_misors.mdf, ' print*,' which you can plot with /code/MDF/RF_misor.' print* do 2000, ij = 1 , ngrain-1 do jk = (ij+1) , ngrain rfvector(1)=qresulta(ij,jk,1)/qresulta(ij,jk,4) rfvector(1)=qresulta(ij,jk,2)/qresulta(ij,jk,4) rfvector(1)=qresulta(ij,jk,3)/qresulta(ij,jk,4) ! print*,'qresulta: ',(qresulta(i,j,jkl),jkl=1,4) write(*,892) ij,jk,lsig(ij,jk) & ,angle(ij,jk) & ,qresulta(ij,jk,1)/qresulta(ij,jk,4) & ,qresulta(ij,jk,2)/qresulta(ij,jk,4) & ,qresulta(ij,jk,3)/qresulta(ij,jk,4) if(qresulta(ij,jk,4).ne.0.) then write(1,892) ij,jk,lsig(ij,jk) & ,angle(ij,jk) & ,qresulta(ij,jk,3)/qresulta(ij,jk,4) & ,qresulta(ij,jk,2)/qresulta(ij,jk,4) & ,qresulta(ij,jk,1)/qresulta(ij,jk,4) & ,1.0,0.,0.,1.0,0.,0. 892 format(3i4,4f8.3,6f5.1) endif c converts the quaternion to Rodrigues vector end do 2000 end do write(*,*) 'Rotate a vector from one lattice to the other [yes=1]?' read(*,*) iq if(iq.eq.1) then write(*,*) 'Enter the vector to be rotated (transformed)' read(*,*) tmpvec call donorm(tmpvec) call rfrvecR(rfvector,tmpvec,tmpout) write(*,"('input vector normalized',3(1x,f8.4))") tmpvec write(*,"('output vector',3(1x,f8.4))") tmpout endif close(1) close(2) call exit(0) end ! ! ............................................................. subroutine uvwhkl2quat( gbnorma , gbdirna , quatn , Bunge , a ) implicit none real scalar real gbnorma(3),gbdirna(3),quatn(4),Bunge(3) real pi,degrad,rtmp,t1(3),a(3,3),d1,d2,d3 integer i,iopt,kerr,ior,j character nomen*1 real tmp real axis(3) logical verbose common /general/ verbose c CODE:: print* pi = 3.14159 degrad = 180./pi call vecnorm(gbnorma) call vecnorm(gbdirna) rtmp = scalar(gbnorma,gbdirna) if(abs(rtmp).gt.1e-4) then stop 'uvwhkl2quat: sorry, uvw is not perp. to hkl' endif call vecpro2(gbnorma,gbdirna,t1) call vecnorm(t1) do 20, i = 1,3 a(3,i) = gbnorma(i) ! 3rd ROW a(2,i) = t1(i) ! 2nd ROW a(1,i) = gbdirna(i) ! 1st ROW 20 enddo ! gives the transpose of what you find in the books ! because this is what Euler needs as input iopt = 1 ! to get angles from matrix nomen = 'B' kerr = 0 ior = 1 ! grain number call euler(a,iopt,nomen,Bunge(1),Bunge(2),Bunge(3),ior,kerr) c note: should use Changsoo's more reliable code if(Bunge(1).lt.0.) Bunge(1) = Bunge(1) + 360. if(Bunge(2).lt.0.) Bunge(2) = Bunge(2) + 360. if(Bunge(3).lt.0.) Bunge(3) = Bunge(3) + 360. ! now we transpose A rtmp = a(1,3) a(1,3) = a(3,1) a(3,1) = rtmp rtmp = a(2,3) a(2,3) = a(3,2) a(3,2) = rtmp rtmp = a(1,2) a(1,2) = a(2,1) a(2,1) = rtmp rtmp = 1. + a(1,1) + a(2,2) + a(3,3) if ( rtmp .lt. 0. ) rtmp = 0. if ( abs(rtmp) .lt. 1.e-4 ) then ! these lines deal with the case that the angle = 180 degr axis(1) = sqrt( ( a(1,1) + 1.) * 0.5 ) axis(2) = sqrt( ( a(2,2) + 1.) * 0.5 ) axis(3) = sqrt( ( a(3,3) + 1.) * 0.5 ) quatn(4) = tmp * 0.25 quatn(3) = axis(3) quatn(2) = axis(2) quatn(1) = axis(1) else ! normal situation tmp = 2. * sqrt ( rtmp ) quatn(4) = tmp * 0.25 quatn(3) = ( a(1,2) - a(2,1) ) / tmp quatn(2) = ( a(3,1) - a(1,3) ) / tmp quatn(1) = ( a(2,3) - a(3,2) ) / tmp ! this implements Morawiec's formula to directly from Eulers to quat ! see the spreadsheet for 27-750, Hwk on orientations endif if (verbose) then print*, $ 'HKLUVW2QUAT: Bunge angles: ',Bunge(1),Bunge(2),Bunge(3) print*,'HKLUVW2QUAT: Quat : ',(quatn(i),i=1,4) endif if (verbose) then print*,'HKLUWV2QUAT: Orientation matrix:' do i = 1,3 print"('[ ',3(1x,f8.3),' ]')",(a(i,j),j=1,3) enddo endif return end 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 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 rod2q(d,quat) implicit none real d(3),quat(4),rtmp,rnorm,rnorm2 c converts Rodrigues vector to a quaternion c CODE:: rtmp=d(1)**2+d(2)**2+d(3)**2 c rnorm=sqrt(rtmp) rnorm2=sqrt(1.+rtmp) quat(4)=1./rnorm2 quat(3)=d(3)/rnorm2 quat(2)=d(2)/rnorm2 quat(1)=d(1)/rnorm2 c cute - this works for all rotation angles! return end c c ____________________________ c c subroutine q2rod(d,quat) implicit none integer i real d(3),quat(4),rtmp,rnorm,rnorm2,stheta,ttheta real large data large / 1.e33 / c converts a quaternion to a Rodrigues vector ! CODE:: rnorm=sqrt(quat(1)**2+quat(2)**2+quat(3)**2+quat(4)**2) if( rnorm .lt. 1.e-5 ) stop 'Q2ROD: too small a quaternion!' if( abs( rnorm - 1.0 ) .gt. 1.e-3 ) then write(*,*) 'Q2ROD: quat required normalization' do i=1,4 quat(i)=quat(i)/rnorm enddo endif if( abs(quat(4)) .gt. 1.e-2 ) then d(1) = quat(1)/quat(4) d(2) = quat(2)/quat(4) d(3) = quat(3)/quat(4) else c we calculate sine(theta/2) from quat1-3 c and scale by tangent/sine in order to make c the system deal with large numbers! c rtmp = sqrt(quat(1)**2+quat(2)**2+quat(3)**2) c this will be close to unity stheta = asin(rtmp) if ( abs(stheta) .gt. 1.e-5 ) then ttheta = tan(stheta) ! and this will be large d(1) = quat(1)*ttheta/stheta d(2) = quat(2)*ttheta/stheta d(3) = quat(3)*ttheta/stheta else d(1) = large d(2) = large d(3) = large endif endif return end c c _____________________________ c c subroutine comquat2(qq1,qq2,qresult) real qq1(4),qq2(4),qresult(4) c c PI=3.14159265 c c algorithm for forming resultant quaternion as Q1 x Q2^-1 c where Q2 follows Q1, i.e. is the second rotation; c note change of signs to get inverse of second orientation c qresult(1)=qq1(1)*qq2(4)-qq1(4)*qq2(1) & +qq1(2)*qq2(3)-qq1(3)*qq2(2) qresult(2)=qq1(2)*qq2(4)-qq1(4)*qq2(2) & +qq1(3)*qq2(1)-qq1(1)*qq2(3) qresult(3)=qq1(3)*qq2(4)-qq1(4)*qq2(3) & +qq1(1)*qq2(2)-qq1(2)*qq2(1) qresult(4)=qq1(4)*qq2(4)+qq1(1)*qq2(1) & +qq1(2)*qq2(2)+qq1(3)*qq2(3) c c write(*,*) 'qresult ',qresult c return end c c c _____________________________ c c subroutine comquat(qq,qresult) real qq(4,2),qresult(4) c c PI=3.14159265 c c algorithm for forming resultant quaternion as Q1 x Q2^-1 c where Q2 follows Q1, i.e. is the second rotation; c note change of signs to get inverse of second orientation c qresult(1)=qq(1,1)*qq(4,2)-qq(4,1)*qq(1,2) & +qq(2,1)*qq(3,2)-qq(3,1)*qq(2,2) qresult(2)=qq(2,1)*qq(4,2)-qq(4,1)*qq(2,2) & +qq(3,1)*qq(1,2)-qq(1,1)*qq(3,2) qresult(3)=qq(3,1)*qq(4,2)-qq(4,1)*qq(3,2) & +qq(1,1)*qq(2,2)-qq(2,1)*qq(1,2) qresult(4)=qq(4,1)*qq(4,2)+qq(1,1)*qq(1,2) & +qq(2,1)*qq(2,2)+qq(3,1)*qq(3,2) write(*,*) 'COMQUAT: qresult ',qresult return end c c _____________________________ c c subroutine comquatinv(qq,qresult) real qq(4,2),qresult(4) c c PI=3.14159265 c c algorithm for forming resultant quaternion as Q1^1 x Q2 c where Q2 follows Q1, i.e. is the second rotation; c note change of signs to get inverse of second orientation c qresult(1)=(-1.*qq(1,1)*qq(4,2))+qq(4,1)*qq(1,2) & +qq(2,1)*qq(3,2)-qq(3,1)*qq(2,2) qresult(2)=(-1.*qq(2,1)*qq(4,2))+qq(4,1)*qq(2,2) & +qq(3,1)*qq(1,2)-qq(1,1)*qq(3,2) qresult(3)=(-1.*qq(3,1)*qq(4,2))+qq(4,1)*qq(3,2) & +qq(1,1)*qq(2,2)-qq(2,1)*qq(1,2) qresult(4)=qq(4,1)*qq(4,2)+qq(1,1)*qq(1,2) & +qq(2,1)*qq(2,2)+qq(3,1)*qq(3,2) c c write(*,*) 'qresult ',qresult c return end c c _____________________________ c c subroutine presymm(qq,lindex,qresult) implicit none real qq(4),qresult(4) integer lindex ! include 'common.f' integer numsymm,numsamp real quatsymm(4,48),quatsamp(4,4) common/a1/ quatsymm,numsymm,quatsamp,numsamp ! CODE:: c$$$c algorithm for forming resultant quaternion c$$$c from applying a symmetry operator, QUATSYMM(n,lindex) c$$$c to the first quaternion/rotation c$$$ c$$$c If the symm oper == O, then c$$$c Qresult = (O x QQ) where QQ= Q1 x Q2^-1 c$$$c i.e. Qresult = (O x Q1) x Q2^-1 c$$$c note change of signs to get inverse of first orientation c$$$ c$$$ if(lindex.gt.numsymm) stop 'error in presymm, lindex>numsymm' c$$$ if(lindex.lt.1) stop 'error in presymm, lindex<1' c$$$ qresult(1) = quatsymm(4,lindex)*qq(1) + quatsymm(1,lindex)*qq(4) c$$$ & + quatsymm(3,lindex)*qq(2) - quatsymm(2,lindex)*qq(3) c$$$ qresult(2) = quatsymm(4,lindex)*qq(2) + quatsymm(2,lindex)*qq(4) c$$$ & + quatsymm(1,lindex)*qq(3) - quatsymm(3,lindex)*qq(1) c$$$ qresult(3) = quatsymm(4,lindex)*qq(3) + quatsymm(3,lindex)*qq(4) c$$$ & + quatsymm(2,lindex)*qq(1) - quatsymm(1,lindex)*qq(2) c$$$ qresult(4) = quatsymm(4,lindex)*qq(4) - quatsymm(1,lindex)*qq(1) c$$$ & - quatsymm(2,lindex)*qq(2) - quatsymm(3,lindex)*qq(3) ! old version 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.numsymm) stop 'error in PRESYMM, lindex>numsamp' if(lindex.lt.1) stop 'error in presymm, lindex<1' qresult(1)=qq(1)*quatsymm(4,lindex) + qq(4)*quatsymm(1,lindex) & - qq(2)*quatsymm(3,lindex) + qq(3)*quatsymm(2,lindex) qresult(2)=qq(2)*quatsymm(4,lindex) + qq(4)*quatsymm(2,lindex) & - qq(3)*quatsymm(1,lindex) + qq(1)*quatsymm(3,lindex) qresult(3)=qq(3)*quatsymm(4,lindex) + qq(4)*quatsymm(3,lindex) & - qq(1)*quatsymm(2,lindex) + qq(2)*quatsymm(1,lindex) qresult(4)=qq(4)*quatsymm(4,lindex) - qq(1)*quatsymm(1,lindex) & - qq(2)*quatsymm(2,lindex) - qq(3)*quatsymm(3,lindex) ! taken from MILL2EUL c write(*,*) 'qresult ',qresult return end c _____________________________ c c subroutine presamp(qq,lindex,qresult) real qq(4),qresult(4) integer lindex integer numsymm,numsamp real quatsymm(4,48),quatsamp(4,4) common/a1/ quatsymm,numsymm,quatsamp,numsamp ! CODE:: c$$$c algorithm for forming resultant quaternion c$$$c from applying a symmetry operator, QUATSYMM(n,lindex) c$$$c to the first quaternion/rotation c$$$ c$$$c If the symm oper == O, then c$$$c Qresult = (O x QQ) where QQ= Q1 x Q2^-1 c$$$c i.e. Qresult = (O x Q1) x Q2^-1 c$$$c note change of signs to get inverse of first orientation c$$$ c$$$ if( lindex .gt. numsamp) stop 'error in presymm, lindex>numsamp' c$$$ if( lindex .lt. 1 ) stop 'error in presamp, lindex<1' c$$$ qresult(1) = quatsamp(4,lindex)*qq(1) + quatsamp(1,lindex)*qq(4) c$$$ & + quatsamp(3,lindex)*qq(2) - quatsamp(2,lindex)*qq(3) c$$$ qresult(2) = quatsamp(4,lindex)*qq(2) + quatsamp(2,lindex)*qq(4) c$$$ & + quatsamp(1,lindex)*qq(3) - quatsamp(3,lindex)*qq(1) c$$$ qresult(3) = quatsamp(4,lindex)*qq(3) + quatsamp(3,lindex)*qq(4) c$$$ & + quatsamp(2,lindex)*qq(1) - quatsamp(1,lindex)*qq(2) c$$$ qresult(4) = quatsamp(4,lindex)*qq(4) - quatsamp(1,lindex)*qq(1) c$$$ & - quatsamp(2,lindex)*qq(2) - quatsamp(3,lindex)*qq(3) ! old version 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) ! taken from MILL2EUL c write(*,*) 'qresult ',qresult return end c c c _____________________________ c c subroutine postsymm(qq,lindex,qresult) real qq(4),qresult(4) integer lindex integer numsymm,numsamp real quatsymm(4,48),quatsamp(4,4) common/a1/ quatsymm,numsymm,quatsamp,numsamp c$$$c algorithm for forming resultant quaternion c$$$c from applying a symmetry operator, QUATSYMM(n,lindex) c$$$c to the second quaternion c$$$ c$$$c If the symm oper == O, then c$$$c Qresult = (QQ x O^-1) where QQ= Q1 x Q2^-1 c$$$c i.e. Qresult = Q1 x (O x Q2)^-1 c$$$c note change of signs to get inverse of first orientation c$$$ c$$$ if(lindex.gt.numsymm) stop 'error in postsymm, lindex>numsymm' c$$$ if(lindex.lt.1) stop 'error in postsymm, lindex<1' c$$$ qresult(1)=qq(1)*quatsymm(4,lindex) - qq(4)*quatsymm(1,lindex) c$$$ & + qq(2)*quatsymm(3,lindex) - qq(3)*quatsymm(2,lindex) c$$$ qresult(2)=qq(2)*quatsymm(4,lindex) - qq(4)*quatsymm(2,lindex) c$$$ & + qq(3)*quatsymm(1,lindex) - qq(1)*quatsymm(3,lindex) c$$$ qresult(3)=qq(3)*quatsymm(4,lindex) - qq(4)*quatsymm(3,lindex) c$$$ & + qq(1)*quatsymm(2,lindex) - qq(2)*quatsymm(1,lindex) c$$$ qresult(4)=qq(4)*quatsymm(4,lindex) + qq(1)*quatsymm(1,lindex) c$$$ & + qq(2)*quatsymm(2,lindex) + qq(3)*quatsymm(3,lindex) ! old version of POSTSYMM 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) ! taken from MILL2EUL ! write(*,*) 'postsymm: qresult ',qresult return end c c c _____________________________ c subroutine disquat2(qq1,qq2,index1,index2,qresult,index1a,index2a) implicit none real qq1(4),qq2(4),qqn(4),qresult(4),quint(4),quintn(4),qmax real q4neg integer index1,index2,index1a,index2a,index3,iq4neg,i,j,k integer ijk , kjl , jjm , jkl c algorithm for finding disorientation corresponding to c input quaternion, QQ (from COMQUAT): resultant quat. is in QRESULT c indices for symmetry operators are in INDEX1, INDEX2 logical keep integer ie , maxeul , eulcount parameter ( maxeul = 99000 ) real euler(maxeul,3) , epseul real d1 , d2 , d3 real eps , pi , rad real rodr(3) , qtemp(4) real qgrain1(4) , qgrain2(4) , axis(3) , samp1(3) , samp2(3) ! CODE:: eps = 1.0e-3 PI = 3.14159265 rad = pi / 180. index1 = 0 index2 = 0 index1a = 0 index2a = 0 index3 = 0 qmax = 0. epseul = 0.1 eulcount = 0 ! no. of new sets of Euler angles c$$$ print"('DISQUAT2 input 1: ',4(1x,f8.4))",qq1 c$$$ print"('DISQUAT2 input 2: ',4(1x,f8.4))",qq2 do 200, i=1,24 ! call presymm(qq1,i,quint) call postsymm(qq1,i,quint) ! now matches DISQUAT do 195, jjm=1,2 if(jjm.eq.2) then do 170, k=1,4 quint(k)=quint(k)*-1.0 170 continue c can take negative of Q because identity = +/-(0,0,0,1) endif do 190, j=1,24 c loop over the two sides of the grain boundary call postsymm(qq2,j,quintn) ! call presymm(qq2,j,quintn) call comquat2(quint,quintn,qqn) ! print"('DISQUAT2 qqn: ',4(1x,f8.4))",qqn axis(1) = qqn(1) axis(2) = qqn(2) axis(3) = qqn(3) call vecnorm(axis) print*,' i jjm j ',i , jjm , j print"('Quat from COMQUAT2 = ',4(2x,f8.3))" 1 ,(qqn( jkl ) , jkl = 1 , 4 ) print"('Misor. Axis= ',3(2x,f8.3))",(axis(ijk),ijk=1,3) do ijk = 1 , 4 qgrain1(ijk)= quint( ijk ) qgrain2(ijk)= quintn( ijk ) enddo call qrvec ( qgrain1 , axis , samp1 ) print"('Qgrain1 = ',4(2x,f8.3))" 1 ,( qgrain1( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 1 = ',3(2x,f8.3))" $ , ( samp1(ijk) , ijk = 1 , 3 ) call qrvec ( qgrain2 , axis , samp2 ) print"('Qgrain2 = ',4(2x,f8.3))" 1 ,( qgrain2( jkl ) , jkl = 1 , 4 ) print"('misor axis in sample 2 = ',3(2x,f8.3))" $ , ( samp2(ijk) , ijk = 1 , 3 ) do 185, kjl=1,2 if(kjl.eq.2) then do 180, k=1,4 qresult(k)=qqn(k)*-1.0 180 continue c can take negative of Q because identity = +/-(0,0,0,1) endif qresult(1)=qqn(1) qresult(2)=qqn(2) qresult(3)=qqn(3) qresult(4)=qqn(4) do 210, ijk=1,2 if(ijk.gt.1) qresult(4)=qresult(4)*-1. c take the inverse rotation, i.e. apply the switching symmetry ! print"('DISQUAT2 qresult: ',4(1x,f8.4))",qresult call q2rod(rodr,qresult) call q2eulB(d1,d2,d3,qresult) d1 = d1 / rad if( d1 .lt. 0. ) d1 = d1 + 360. d2 = d2 / rad if( d2 .lt. 0. ) d2 = d2 + 360. d3 = d3 / rad if( d3 .lt. 0. ) d3 = d3 + 360. c$$$ write(3,"(4(1x,g10.3),3i4,4(1x,g10.3))") c$$$ $ , d1 , d2 , d3 , 1. , i , j , kjl , qresult keep = .false. ! if( i .eq. 1 .and. j .eq. 1 .and. ijk .eq. 2 ) ! if( kjl .lt. 2 .and. ijk .lt. 2 ) if( ijk .eq. 2 ) $ keep = .true. ! debug ! keep = .true. if ( eulcount .gt. 0 ) then do ie = 1 , eulcount if ( abs(d1-euler(ie,1)) .lt. epseul $ .and. abs(d2-euler(ie,2)) .lt. epseul $ .and. abs(d3-euler(ie,3)).lt.epseul) then keep = .false. endif enddo endif if ( keep .and. (eulcount .lt. maxeul) ) then eulcount = eulcount + 1 euler(eulcount,1) = d1 euler(eulcount,2) = d2 euler(eulcount,3) = d3 c$$$ print"('DISQUAT2: I/J/KJL/RodrV/Euler: ' c$$$ $ ,3i4,6(1x,g10.3))" c$$$ $ ,i,j,kjl,rodr, d1 , d2 , d3 write(3,"(4(1x,g10.3),3i4,4(1x,g10.3))") $ d1 , d2 , d3 , 1. , i , j , kjl , qresult write(4,"(2(1x,i4),8(1x,g10.3))") $ i , j , qresult , 1. , d1 , d2 , d3 ! write(4,*) ' i j q1 q2 q3 q4 E1 E2 E3 / for plotting with RF_misor ' endif c$$$ write(*,"('i,j,kjl,q4neg,qresult ',3i3,5f8.3)") c$$$ $ i,j,kjl,q4neg,qresult if(qresult(1).ge.0.0) then c$$$ if((qresult(1)-qresult(2)).lt.eps) then c$$$ if((qresult(2)-qresult(3)).lt.eps) then c$$$ if((qresult(3)-qresult(4)).lt.eps) then if((qresult(1) .le. qresult(2))) then if((qresult(2) .le. qresult(3))) then if((qresult(3) .le. qresult(4))) then c these IFs ensure that 0<=q1R2>R3 return end c c c ----------------------------------------- c C subroutine ZUR BERECHNUNG VON ORIENTIERUNGSBEZIEHUNGEN ZWISCHEN C ZWEI VORGEGEBENEN ORIENTIERUNGEN c c subroutine provided by Dierk Raabe (to Paul Lee, 1996) c subroutine misori(p1,p,p2,al_min,nsymm,ac_min) implicit none c p1,p,p2 are (2 sets of) Bunge Euler angles c nsymm controls the sample symmetry c al_min is the disorientation angle c character meu,isy integer ac_min(3) integer i,j,k,i31max,nsymm,i31,i1,i24,i25,iz integer iomega,i_min,ii real D(3,3,2),SYM(3,3,24),P1(2),P(2),P2(2) real DM(3,3),DRR(3,3,24),A(3),IA(3) real DR(3,3) real tmp, tmp1, tmp2, tmp3 ,pi ,al_min ,fakt ,f real spur,sp,omega,alpha,x real c1,s1,pp1,pp,pp2,c,c2,s2,s PI=3.14159265 al_min=360. ac_min(1)=0 ac_min(2)=0 ac_min(3)=0 C ERSTELLEN DER SYMMETRIEMATRIZEN C ************************************************************** DO 2 I=1,3 DO 2 J=1,3 DO 2 K=1,24 2 SYM(I,J,K)=0. C 1 SYM(1,1,1)=1. SYM(2,2,1)=1. SYM(3,3,1)=1. C 5 SYM(1,1,2)=1. SYM(2,3,2)=-1. SYM(3,2,2)=1. C 2 SYM(1,1,3)=1. SYM(2,2,3)=-1. SYM(3,3,3)=-1. C 11 SYM(1,1,4)=1. SYM(2,3,4)=1. SYM(3,2,4)=-1. C 7 SYM(1,3,5)=-1. SYM(2,2,5)=1. SYM(3,1,5)=1. C 12 SYM(1,3,6)=1. SYM(2,2,6)=1. SYM(3,1,6)=-1. C 3 SYM(1,1,7)=-1. SYM(2,2,7)=1. SYM(3,3,7)=-1. C 4 SYM(1,1,8)=-1. SYM(2,2,8)=-1. SYM(3,3,8)=1. C 13 SYM(1,2,9)=1. SYM(2,1,9)=-1. SYM(3,3,9)=1. C 6 SYM(1,2,10)=-1. SYM(2,1,10)=1. SYM(3,3,10)=1. C 20 SYM(1,2,11)=-1. SYM(2,3,11)=1. SYM(3,1,11)=-1. C 23 SYM(1,3,12)=1. SYM(2,1,12)=-1. SYM(3,2,12)=-1. C 19 SYM(1,2,13)=-1. SYM(2,3,13)=-1. SYM(3,1,13)=1. C 21 SYM(1,3,14)=-1. SYM(2,1,14)=1. SYM(3,2,14)=-1. C 18 SYM(1,2,15)=1. SYM(2,3,15)=-1. SYM(3,1,15)=-1. C 22 SYM(1,3,16)=-1. SYM(2,1,16)=-1. SYM(3,2,16)=1. C 9 SYM(1,2,17)=1. SYM(2,3,17)=1. SYM(3,1,17)=1. C 8 SYM(1,3,18)=1. SYM(2,1,18)=1. SYM(3,2,18)=1. C 17 SYM(1,2,19)=1. SYM(2,1,19)=1. SYM(3,3,19)=-1. C 24 SYM(1,1,20)=-1. SYM(2,3,20)=1. SYM(3,2,20)=1. C 10 SYM(1,3,21)=1. SYM(2,2,21)=-1. SYM(3,1,21)=1. C 15 SYM(1,1,22)=-1. SYM(2,3,22)=-1. SYM(3,2,22)=-1. C 16 SYM(1,3,23)=-1. SYM(2,2,23)=-1. SYM(3,1,23)=-1. C 14 SYM(1,2,24)=-1. SYM(2,1,24)=-1. SYM(3,3,24)=-1. C ********************************** C ERSTELLEN DER BEIDEN DMATRIZEN c 205 FORMAT(A1) meu='E' isy='E' I31MAX=nsymm DO 31 I31=1,I31MAX DO 1 I1=1,2 FAKT=PI/180. PP1=P1(I1)*FAKT PP=P(I1)*FAKT PP2=P2(I1)*FAKT C1=COS(PP1) C=COS(PP) C2=COS(PP2) S1=SIN(PP1) S=SIN(PP) S2=SIN(PP2) D(1,1,I1)=C1*C2-S1*S2*C D(1,2,I1)=S1*C2+C1*S2*C D(1,3,I1)=S2*S D(2,1,I1)=-C1*S2-S1*C2*C D(2,2,I1)=-S1*S2+C1*C2*C D(2,3,I1)=C2*S D(3,1,I1)=S1*S D(3,2,I1)=-C1*S D(3,3,I1)=C ! GOTO 1 1 CONTINUE C ********************************** DO 24 I24=1,3 DO 25 I25=1,3 DM(I24,I25)=D(I24,I25,1) DR(I24,I25)=D(I24,I25,2) 25 CONTINUE 24 CONTINUE F=PI/180. print* print*,'1st matrix:' do i24 = 1,3 print"('[',3(2x,f8.3),' ]')",(DM(I24,I25), i25 = 1,3) enddo print* print*,'2nd matrix:' do i24 = 1,3 print"('[',3(2x,f8.3),' ]')",(dr(I24,I25), i25 = 1,3) enddo C **************************************************** C BESTIMMUNG DER INVERSEN MATRIX ZUR ORIENTIERUNG 1:DM C **************************************************** DO 3 I=1,3 DO 3 J=1,3 3 DM(I,J)=D(J,I,1) C ************************* C FELDER GLEICH NULL SETZEN C ************************* DO 4 I=1,3 DO 4 J=1,3 DO 4 K=1,24 DR(I,J)=0. 4 DRR(I,J,K)=0. C *********************************************************** C MATRIZENMULTIPLIKATION DER MATRIZEN D(I,J,2) UND DM=DR(I,J) C *********************************************************** DO 5 I=1,3 DO 5 J=1,3 DO 6 K=1,3 6 DR(I,J)=D(I,K,2)*DM(K,J)+DR(I,J) 5 CONTINUE c DM is associated with the first orientation c and has been transposed so that the misorientation is calculated properly c thus, we have product matrix as gB X gA^-1 C *********************************************** C MATRIZENMULITIPLIKATION DER MATRIZEN SYM UND DR C *********************************************** DO 7 IZ=1,24 DO 8 I=1,3 DO 8 J=1,3 DO 9 K=1,3 9 DRR(I,J,IZ)=SYM(I,K,IZ)*DR(K,J)+DRR(I,J,IZ) 8 CONTINUE c$$$ print* c$$$ print*,' Symmetry operator number ',iz c$$$ print*,'Product matrix for gA X gB^-1:' c$$$ do i24 = 1,3 c$$$ print"('[',3(2x,f8.3),' ]')",(drr(i24,i25,iz), i25 = 1,3) c$$$ enddo c$$$ tmp = (drr(1,1,iz)+drr(2,2,iz)+drr(3,3,iz)) c$$$ print*,'Trace = ',tmp c$$$ tmp2 = (tmp-1)/2. c$$$ if(tmp2.lt.-1.) tmp2 = -1. c$$$ if(tmp2.gt.1.) tmp2 = 1. c$$$ tmp3 = acos(tmp2)*180./pi c$$$ print*,' angle = ',tmp3 C ******************************* C BESTIMMUNG DES ROTATIONSWINKELS C ******************************* SPUR=0. DO 10 I=1,3 10 SPUR=SPUR+DRR(I,I,IZ) SP=(SPUR-1.)/2. IF(SP.GE.1.)SP=1.-1.E-5 IF(SP.LE.-1.)SP=-1.+1.E-5 OMEGA=PI/2.-ASIN(SP) IOMEGA=INT(OMEGA*180./PI+0.5) ALPHA=OMEGA*180./PI IF(IOMEGA.EQ.180.or.iomega.eq.0)GOTO 12 X=2.*SIN(OMEGA) C ******************************************************** C BESTIMMUNG DER ROTATIONSACHSE, INCL SONDERFALL OMEGA = 180 C ******************************************************** A(1)=100.*(DRR(2,3,IZ)-DRR(3,2,IZ))/X A(2)=100.*(DRR(3,1,IZ)-DRR(1,3,IZ))/X A(3)=100.*(DRR(1,2,IZ)-DRR(2,1,IZ))/X GOTO 13 12 A(1)=SQRT((DRR(1,1,IZ)+1.)/2.)*100. A(2)=SQRT((DRR(2,2,IZ)+1.)/2.)*100. A(3)=SQRT((DRR(3,3,IZ)+1.)/2.)*100. 13 DO 11 I=1,3 IF(A(I).EQ.0.)GOTO 19 IA(I)=INT(A(I)+0.5*A(I)/ABS(A(I))) GOTO 11 19 IA(I)=0 11 CONTINUE if(abs(alpha).lt.al_min) then al_min=abs(alpha) ac_min(1)=ia(1) ac_min(2)=ia(2) ac_min(3)=ia(3) end if 7 CONTINUE IF(I31.EQ.1)P1(2)=180.+P1(2) IF(I31.EQ.2)P1(2)=360.-P1(2) IF(I31.EQ.2)P2(2)=90.-P2(2) IF(I31.EQ.3)P1(2)=180.+P1(2) 31 CONTINUE i_min=100 if(abs(ac_min(1)).lt.i_min.and.abs(ac_min(1)).gt.0) $i_min=abs(ac_min(1)) if(abs(ac_min(2)).lt.i_min.and.abs(ac_min(2)).gt.0) $i_min=abs(ac_min(2)) if(abs(ac_min(3)).lt.i_min.and.abs(ac_min(3)).gt.0) $i_min=abs(ac_min(3)) do 5501 ii=1,3 5501 ac_min(ii)=nint(float(ac_min(ii))/float(i_min)) if(p1(1).gt.360.) p1(1)=p1(1)-360. if(p1(1).lt.0.) p1(1)=p2(1)+360. if(p2(1).gt.360.) p2(1)=p1(1)-360. if(p2(1).lt.0.) p2(1)=p2(1)+360. if(p1(2).gt.360.) p1(2)=p1(2)-360. if(p1(2).lt.0.) p1(2)=p2(2)+360. if(p2(2).gt.360.) p2(2)=p1(2)-360. if(p2(2).lt.0.) p2(2)=p2(2)+360. return end c c ------------------- c subroutine donorm(a) real a(3),b(3),rsum,rnorm rsum=0. do 100, i=1,3 rsum=rsum+a(i)**2 100 enddo if (rsum .lt. 1.e-5 ) stop 'zero length axis?!' rnorm=sqrt(rsum) do 200, i=1,3 b(i)=a(i)/rnorm a(i)=b(i) 200 enddo return end c c ______________________ c subroutine rfrvecR(rf,din,dout) c c uses RF-vector to rotate a vector from DIN to DOUT; c Eq. 1.36 in Sutton & Balluffi c REVERSE/NEGATIVE ROTATION; based on rfrvec; 7 vi 03 c real rf(3),din(3),dout(3) real tmp,tmp2,cross(3),rftmp(3) integer i c do i=1,3 rftmp(i)=-1.*rf(i) enddo tmp=2.*scalar(rftmp,din) tmp2=rftmp(1)**2+rftmp(2)**2+rftmp(3)**2 call vecpro2(din,rftmp,cross) dout(1)=(((tmp*rftmp(1))+(din(1)*(1.-tmp2))) &-(2.*cross(1)))/(1.+tmp2) dout(2)=(((tmp*rftmp(2))+(din(2)*(1.-tmp2))) &-(2.*cross(2)))/(1.+tmp2) dout(3)=(((tmp*rftmp(3))+(din(3)*(1.-tmp2))) &-(2.*cross(3)))/(1.+tmp2) return end c c ______________________ c subroutine q2eulB(p1,p,p2,qq) c converts quaternion to Bunge Euler angles c based on Altmann's solution for Euler->quat real qq(4),p1,p,p2 real sum,diff PI=3.14159265 if( ( abs(qq(2)) .lt. 1.e-35 ) $ .and. ( abs(qq(1)) .lt. 1.e-35 ) ) then ! diff=pi/4. diff = 0. ! looks better ... else diff = atan2( qq(2) , qq(1) ) endif c ATAN2 is unhappy is both arguments are exactly zero! if( ( abs(qq(3)) .lt. 1.e-35 ) $ .and. ( abs(qq(4)) .lt. 1.e-35 )) then ! sum=pi/4. sum = 0. else sum = atan2( qq(3) , qq(4) ) endif p1 = (diff+sum) p2 = (sum-diff) tmp = sqrt(qq(3)**2+qq(4)**2) if(tmp.gt.1.0) tmp=1.0 p = 2.*acos(tmp) 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 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 textur2(afile) c second version, for use in REXGBS, WTS2POP c not for use in REX2D implicit none integer i,j,k,ijk ,ng real rad ,d1 ,d2 ,d3 c real quatsymm,quatsamp,quatcsl,aquat,phi1,capphi,phi2,grwt real quatsymm(4,48) ,quatsamp(4,4) integer numsymm ,numsamp common/a1/ quatsymm,numsymm,quatsamp,numsamp real quatcsl(4,50) integer numcsl ,nsig(50) common/a2/quatcsl ,numcsl ,nsig c$$$ common/a1/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp c$$$ common/a2/quatcsl(4,50),numcsl,nsig(50) real aquat(4,100000) & ,phi1(100000),capphi(100000),phi2(100000),grwt(100000) integer ngrain common/a3/aquat ,phi1 ,capphi ,phi2 ,grwt ,ngrain c$$$ common/a3/aquat(4,100000),ngrain c$$$ & ,phi1(100000),capphi(100000),phi2(100000),grwt(100000) c common/a4/ textitl,eulan,iper c include 'common.f' c following code taken from LApp68.f for reading in TEXIN real evm,fk1(9),nstate character*80 textitl character eulan*25,nomen,iper*6 real ph,th,ps,dph,dth,dps,sph,cph,cth,sth,sps,cps,qr(4) real theta,r1,r2,r3 integer lmn character afile*100 rad=57.29578 ! symmetry now read in by textur4 write(*,*) 'reading CSL information from "quat.csl.data"' open(unit=3,name='quat.csl.data',status='old') read(3,*) ! skip over title line read(3,*) numcsl ! read number of CSL types if(numcsl.gt.50) then numcsl=50 write(*,*) 'too many CSL types, cutting back to 50!' else write(*,*) 'reading ',numcsl,' boundary types' endif read(3,*) ! skip over comment do 1900, i=1,numcsl read(3,*) nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4) c write(*,*) i,nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4) 1900 continue close(unit=3) write(*,*) 'reading grain orientations from ',afile c open(unit=3,name='texin',status='old') open(unit=3, file=afile, status='old') read(3,7) textitl print*,'Header Line: ',textitl 7 format(a) read(3,*) read(3,70) evm,fk1,nstate 70 format(10f7.3,i3) read(3,71) eulan,iper 71 format(a25,49x,a6) nomen=eulan(1:1) if(nomen.eq.'B') then write(*,*) 'detected Bunge angles' elseif(nomen.eq.'K') then write(*,*) 'detected Kocks angles' else print*,'Unknown Euler angle type, quitting' call exit(1) endif do 529 ng=1,100000 c do 520 i=1,18 c 520 state(i,ng)=0. c read(3,50,err=530,end=530) d1,d2,d3,grwt(ng),(state(i,ng),i=1,nstate) read(3,*,end=530,err=530) d1,d2,d3,grwt(ng) write(*,*) 'angles.. ',d1,d2,d3 c if((d1.eq.999.).and.(d2.eq.999.).and.(d3.eq.999.)) goto 530 50 format(6f8.2,24f6.2) c read(3,*,end=530) d1,d2,d3,grwt(ng) c .: for old unformatted TEXIN files c if(nomen.eq.'B') then phi1(ng)=d1 capphi(ng)=d2 phi2(ng)=d3 elseif(nomen.eq.'K') then phi1(ng)=d1+90. capphi(ng)=d2 phi2(ng)=90.-d3 endif if(phi1(ng).gt.360.) phi1(ng) = phi1(ng)-360. if(phi1(ng).lt.0.) phi1(ng) = phi1(ng)+360. if(phi2(ng).gt.360.) phi2(ng) = phi2(ng)-360. if(phi2(ng).lt.0.) phi2(ng)=phi2(ng)+360. 20 dth=d2 dps=d1 call quatB(phi1(ng)/rad,capphi(ng)/rad,phi2(ng)/rad,qr) c use the new conversion routine based on Altmann if Bunge c else, use longer routine do 527 j=1,4 527 aquat(j,ng)=qr(j) 529 continue 530 continue ngrain=ng-1 return end c c ________________________________________ c c subroutine textur4(iq_symm) c second version, for use in REXGBS, WTS2POP c not for use in REX2D integer numsymm,numsamp integer iq_symm,n_max real quatsymm(4,48),quatsamp(4,4) common/a1/ quatsymm,numsymm,quatsamp,numsamp ! CODE:: if(iq_symm.lt.0 .or. iq_symm.gt.2) stop 'bad symmetry class?' if(iq_symm.eq.0) then n_max = 24 elseif(iq_symm.eq.1) then n_max = 12 ! Hex elseif(iq_symm.eq.2) then n_max = 4 ! Orth endif ! confusing - wts2pop uses 1 to 3 c print*, 'reading symmetry operators from "quat.symm.cubic"' if(iq_symm.eq.0) then open(unit=3,name='quat.symm.cubic',status='old') elseif(iq_symm.eq.1) then open(unit=3,name='quat.symm.hex',status='old') elseif(iq_symm.eq.2) then open(unit=3,name='quat.symm.orth',status='old') endif read(3,*) ! skip over title line read(3,*) numsymm ! read number of symmetry operators if(numsymm.gt.n_max) then numsymm = n_max print*, 'too many symmetry operators, cutting back to ',n_max else c print*, 'reading ',numsymm,' symmetry operators' endif do 1800, i=1,numsymm read(3,*) (quatsymm(ijk,i),ijk=1,4) c print*, 'symm. no. ',i,(quatsymm(ijk,i),ijk=1,4) 1800 continue close(unit=3) c now we hard-wire the four operators into the code c so that choices by user are easier to deal with numsamp = 4 do i = 1,numsamp do ijk = 1,4 quatsamp(ijk,i) = 0. enddo enddo quatsamp(4,1) = 1. quatsamp(1,2) = 1. quatsamp(2,3) = 1. quatsamp(3,4) = 1. return end c c ________________________________________ c c subroutine textur6 implicit none c second version, for use in REXGBS, WTS2POP c not for use in REX2D ! edited because textur4 now reads in symmetry operators integer i,ijk,numsymm,numsamp real quatsymm(4,48),quatsamp(4,4) real quatcsl(4,50),theta,r1,r2,r3 integer numcsl,nsig(50),lmn common/a2/quatcsl,numcsl,nsig ! CODE:: write(*,*) 'reading CSL information from "quat.csl.data"' open(unit=3,name='quat.csl.data',status='old') read(3,*) ! skip over title line read(3,*) numcsl ! read number of CSL types if(numcsl.gt.50) then numcsl=50 write(*,*) 'too many CSL types, cutting back to 50!' else write(*,*) 'reading ',numcsl,' boundary types' endif read(3,*) ! skip over comment do 1900, i=1,numcsl read(3,*) nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4) c write(*,*) i,nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4) 1900 continue close(unit=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 character nomen save kor 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 elseif(nomen.eq.'B') then dps=d1-90. dph=90.-d3 else dth=d1 dph=d2-90. dps=90.-d3 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 ! note that this is the transpose of the standard definition ! because of the way in which Canova used the matrix return end c ----------------------------------------- c C subroutine ZUR BERECHNUNG VON ORIENTIERUNGSBEZIEHUNGEN ZWISCHEN C ZWEI VORGEGEBENEN ORIENTIERUNGEN c c subroutine provided by Dierk Raabe (to Paul Lee, 1996) c subroutine misorinv(p1,p,p2,al_min,nsymm,ac_min) c c p1,p,p2 are (2 sets of) Bunge Euler angles c nsymm controls the sample symmetry c al_min is the disorientation angle c character meu,isy integer ac_min(3) DIMENSION D(3,3,2),SYM(3,3,24),P1(2),P(2),P2(2) $,DM(3,3),DRR(3,3,24),A(3),IA(3) $,DR(3,3) real tmp, tmp1, tmp2, tmp3 real al_min PI=3.14159265 al_min=360. ac_min(1)=0 ac_min(2)=0 ac_min(3)=0 C ERSTELLEN DER SYMMETRIEMATRIZEN C ************************************************************** DO 2 I=1,3 DO 2 J=1,3 DO 2 K=1,24 2 SYM(I,J,K)=0. C 1 SYM(1,1,1)=1. SYM(2,2,1)=1. SYM(3,3,1)=1. C 5 SYM(1,1,2)=1. SYM(2,3,2)=-1. SYM(3,2,2)=1. C 2 SYM(1,1,3)=1. SYM(2,2,3)=-1. SYM(3,3,3)=-1. C 11 SYM(1,1,4)=1. SYM(2,3,4)=1. SYM(3,2,4)=-1. C 7 SYM(1,3,5)=-1. SYM(2,2,5)=1. SYM(3,1,5)=1. C 12 SYM(1,3,6)=1. SYM(2,2,6)=1. SYM(3,1,6)=-1. C 3 SYM(1,1,7)=-1. SYM(2,2,7)=1. SYM(3,3,7)=-1. C 4 SYM(1,1,8)=-1. SYM(2,2,8)=-1. SYM(3,3,8)=1. C 13 SYM(1,2,9)=1. SYM(2,1,9)=-1. SYM(3,3,9)=1. C 6 SYM(1,2,10)=-1. SYM(2,1,10)=1. SYM(3,3,10)=1. C 20 SYM(1,2,11)=-1. SYM(2,3,11)=1. SYM(3,1,11)=-1. C 23 SYM(1,3,12)=1. SYM(2,1,12)=-1. SYM(3,2,12)=-1. C 19 SYM(1,2,13)=-1. SYM(2,3,13)=-1. SYM(3,1,13)=1. C 21 SYM(1,3,14)=-1. SYM(2,1,14)=1. SYM(3,2,14)=-1. C 18 SYM(1,2,15)=1. SYM(2,3,15)=-1. SYM(3,1,15)=-1. C 22 SYM(1,3,16)=-1. SYM(2,1,16)=-1. SYM(3,2,16)=1. C 9 SYM(1,2,17)=1. SYM(2,3,17)=1. SYM(3,1,17)=1. C 8 SYM(1,3,18)=1. SYM(2,1,18)=1. SYM(3,2,18)=1. C 17 SYM(1,2,19)=1. SYM(2,1,19)=1. SYM(3,3,19)=-1. C 24 SYM(1,1,20)=-1. SYM(2,3,20)=1. SYM(3,2,20)=1. C 10 SYM(1,3,21)=1. SYM(2,2,21)=-1. SYM(3,1,21)=1. C 15 SYM(1,1,22)=-1. SYM(2,3,22)=-1. SYM(3,2,22)=-1. C 16 SYM(1,3,23)=-1. SYM(2,2,23)=-1. SYM(3,1,23)=-1. C 14 SYM(1,2,24)=-1. SYM(2,1,24)=-1. SYM(3,3,24)=-1. C ********************************** C ERSTELLEN DER BEIDEN DMATRIZEN c 205 FORMAT(A1) meu='E' isy='E' I31MAX=nsymm DO 31 I31=1,I31MAX DO 1 I1=1,2 FAKT=PI/180. PP1=P1(I1)*FAKT PP=P(I1)*FAKT PP2=P2(I1)*FAKT C1=COS(PP1) C=COS(PP) C2=COS(PP2) S1=SIN(PP1) S=SIN(PP) S2=SIN(PP2) D(1,1,I1)=C1*C2-S1*S2*C D(1,2,I1)=S1*C2+C1*S2*C D(1,3,I1)=S2*S D(2,1,I1)=-C1*S2-S1*C2*C D(2,2,I1)=-S1*S2+C1*C2*C D(2,3,I1)=C2*S D(3,1,I1)=S1*S D(3,2,I1)=-C1*S D(3,3,I1)=C GOTO 1 1 CONTINUE C ********************************** c DO 24 I24=1,3 c DO 25 I25=1,3 c DM(I24,I25)=D(I24,I25,1) c DR(I24,I25)=D(I24,I25,2) c 25 CONTINUE c24 CONTINUE F=PI/180. C **************************************************** C BESTIMMUNG DER INVERSEN MATRIX ZUR ORIENTIERUNG 1:DM C **************************************************** DO 3 I=1,3 DO 3 J=1,3 3 DM(I,J)=D(J,I,1) C ************************* C FELDER GLEICH NULL SETZEN C ************************* DO 4 I=1,3 DO 4 J=1,3 DO 4 K=1,24 DR(I,J)=0. 4 DRR(I,J,K)=0. C *********************************************************** C MATRIZENMULTIPLIKATION DER MATRIZEN D(I,J,2) UND DM=DR(I,J) C *********************************************************** c DO 5 I=1,3 c DO 5 J=1,3 c DO 6 K=1,3 c6 DR(I,J)=D(I,K,2)*DM(K,J)+DR(I,J) c5 CONTINUE C *********************************************** C MATRIZENMULITIPLIKATION DER MATRIZEN SYM UND DR C *********************************************** DO 7 IZ=1,24 DO I=1,3 DO J=1,3 dr(i,j) = 0. ! just to make sure DO K=1,3 c DR(I,J)=D(i,k,2)*SYM(k,j,IZ) + DR(I,J) DR(I,J) = SYM(i,k,IZ)*D(k,j,2) + DR(I,J) c apply symmetry before combining enddo enddo enddo c8 CONTINUE DO I=1,3 DO J=1,3 drr(i,j,iz) = 0. ! just to make sure DO K=1,3 drr(I,J,iz)=dr(k,i)*DM(j,k)+drr(I,J,iz) c and use transpose of DM, because it is already transposed c so as to have the product matrix as gB X gA^-1 (but both transposed!) enddo enddo enddo c5 CONTINUE c$$$ print* c$$$ print*,' Symmetry operator number ',iz c$$$ print*,'Product matrix for gB X gA^-1:' c$$$ do i24 = 1,3 c$$$ print"('[',3(2x,f8.3),' ]')",(drr(i24,i25,iz), i25 = 1,3) c$$$ enddo c$$$ tmp = (drr(1,1,iz)+drr(2,2,iz)+drr(3,3,iz)) c$$$ print*,'Trace = ',tmp c$$$ tmp2 = (tmp-1)/2. c$$$ if(tmp2.lt.-1.) tmp2 = -1. c$$$ if(tmp2.gt.1.) tmp2 = 1. c$$$ tmp3 = acos(tmp2)*180./pi c$$$ print*,' angle = ',tmp3 C ******************************* C BESTIMMUNG DES ROTATIONSWINKELS C ******************************* SPUR=0. DO 10 I=1,3 10 SPUR=SPUR+DRR(I,I,IZ) SP=(SPUR-1.)/2. IF(SP.GE.1.)SP=1.-1.E-5 IF(SP.LE.-1.)SP=-1.+1.E-5 OMEGA=PI/2.-ASIN(SP) IOMEGA=INT(OMEGA*180./PI+0.5) ALPHA=OMEGA*180./PI IF(IOMEGA.EQ.180.or.iomega.eq.0)GOTO 12 X=2.*SIN(OMEGA) C ******************************************************** C BESTIMMUNG DER ROTATIONSACHSE, INCL SONDERFALL OMEGA = 180 C ******************************************************** A(1)=100.*(DRR(2,3,IZ)-DRR(3,2,IZ))/X A(2)=100.*(DRR(3,1,IZ)-DRR(1,3,IZ))/X A(3)=100.*(DRR(1,2,IZ)-DRR(2,1,IZ))/X GOTO 13 12 A(1)=SQRT((DRR(1,1,IZ)+1.)/2.)*100. A(2)=SQRT((DRR(2,2,IZ)+1.)/2.)*100. A(3)=SQRT((DRR(3,3,IZ)+1.)/2.)*100. 13 DO 11 I=1,3 IF(A(I).EQ.0.)GOTO 19 IA(I)=INT(A(I)+0.5*A(I)/ABS(A(I))) GOTO 11 19 IA(I)=0 11 CONTINUE if(abs(alpha).lt.al_min) then al_min=abs(alpha) ac_min(1)=ia(1) ac_min(2)=ia(2) ac_min(3)=ia(3) end if 7 CONTINUE IF(I31.EQ.1)P1(2)=180.+P1(2) IF(I31.EQ.2)P1(2)=360.-P1(2) IF(I31.EQ.2)P2(2)=90.-P2(2) IF(I31.EQ.3)P1(2)=180.+P1(2) 31 CONTINUE i_min=100 if(abs(ac_min(1)).lt.i_min.and.abs(ac_min(1)).gt.0) $i_min=abs(ac_min(1)) if(abs(ac_min(2)).lt.i_min.and.abs(ac_min(2)).gt.0) $i_min=abs(ac_min(2)) if(abs(ac_min(3)).lt.i_min.and.abs(ac_min(3)).gt.0) $i_min=abs(ac_min(3)) do 5501 ii=1,3 5501 ac_min(ii)=nint(float(ac_min(ii))/float(i_min)) return end c ______________________ subroutine qrvec(qq,din,dout) c uses quaternion to rotate a vector from DIN to DOUT ! note the minus sign in the line with the permutation tensor ! we use this version to generate pole figures, which means that ! we are effectively performing the reverse transformation ! from crystal to sample axes (frames) ! cleaned up code vi 08, ADR c note that -q gives the same result as it should (represents same rotation) real qq(4),din(3),dout(3) real t1 , tmp c CODE:: ! print*,'qrvec: qq: ',qq t1 = qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2 do i = 1,3 dout(i) = t1*din(i) tmp = 0. do j = 1,3 tmp = tmp + (qq(j)*din(j)) enddo dout(i) = dout(i) + 2. * qq(i) * tmp do j = 1,3 if(i.ne.j) then do 80, ijk = 1,3 if(i.ne.ijk.and.j.ne.ijk) k=ijk 80 continue kl = j-i if( kl .eq. 2 ) kl = -1 if( kl .eq. -2 ) kl = 1 c poor man's permutation tensor! ! print*,'i,j,k,kl: ',i,j,k,kl dout(i) = dout(i) - (2.*qq(k)*qq(4)*float(kl)*din(j)) endif ! if(i.ne.j) then enddo enddo return end c c _____________________________ c }}}} subroutine invqrvec( qq , din , dout ) c uses quaternion to rotate a vector from DIN to DOUT; c but in the inverse sense so as to give sample to crystal axes c derived from qrvec, above routine c note that -q gives the same result as it should (represents same rotation) real qq(4),din(3),dout(3) real t1 c CODE:: t1 = qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2 do i = 1,3 dout(i) = t1*din(i) tmp = 0. do j = 1,3 tmp = tmp + (qq(j)*din(j)) enddo dout(i) = dout(i) + 2. * qq(i) * tmp do j = 1,3 if(i.ne.j) then do 80, ijk = 1,3 if(i.ne.ijk.and.j.ne.ijk) k=ijk 80 continue kl = j-i if( kl .eq. 2 ) kl = -1 if( kl .eq. -2 ) kl = 1 c poor man's permutation tensor! ! print*,'i,j,k,kl: ',i,j,k,kl dout(i) = dout(i) + (2.*qq(k)*qq(4)*float(kl)*din(j)) ! the ONLY difference between this and QRVEC is the + instead of - above endif ! if(i.ne.j) then enddo enddo return end c c c ________________________________________ c ! include 'misc.f'