program voxel2avs ! compile: g77 -O3 -o voxel2avs voxel2avs-[DATE].f ! or, g77 -O0 -g -ffortran-bounds-check -o voxel2avs voxel2avs-[DATE].f c this program converts a voxel file (as in 3d Monte Carlo) c to an AVS format, suitable for reading in to LaGrit c ADR June 04 c added VTK output, Jan 07, Nov 07 c both basic point output and unstructured grid ! added inverse pole figure color scheme, similar to TSL ! leave spin <= 0 empty ! Jan 08, add FFT output (list of Euler angles, coords, material) ! Feb 08, add GB analysis in the form of a .vtk file ! for use with FFT2dx if you like ! May08 fixed problems in textur3 ! Jul 08, added DX output, with padding, for Sukbin's meshing code ! Oct 08, experimented with COLOR_SCALARS; found out that if you ! write a set of triples (one per point) with range 0..1, paraview ! interprets them as unsigned chars in the range 0..255 (huh?!) ! and you can use them to color points as you wish ! tested this with simple.spn and with combination of ! Cleaned_Micro_Random_11Dec06.dx and ! randomized_grain_orientation.txt as input ! Mar 09 - can now read Groeber files (Structure + list or orientations) ! Apr 09 - flipping orientations read from Groeber files according to TSL ! convention for reference frames ! May 09 - fixed problem with logic of colour scheme ! June 09 - minor adjustments to the FFT output implicit none integer i,j,k,ii,jj integer size , nbors , xsize , ysize , zsize parameter (size = 390, nbors = 26 , xsize = 1024 $ , ysize = 512 , zsize = 32 ) integer nnbors(nbors,3) character fname*40 integer isite,jsite,ksite real dum1,dum2,dum3 integer xx,yy,zz integer*2 spins( xsize , ysize , zsize ) common / structure / xx , yy , zz , nnbors , spins integer npoints,nelemavs,nvalues,nvalues_e character e_type*3 real xcoord,ycoord,zcoord ! coordinates integer i1,i2,j1,j2,k1,k2 ! coords integer n1,n2,n3,n4,n5,n6,n7,n8 ! node numbers integer icount integer spin_min,spin_max ! min, max values of spin integer irange,irange2 real range,range2 real red,green,blue integer n_type ! node type; 0=interior 2=interface 10=boundary 12=interface+boundary integer n_mtl ! node material (taken from voxel above node) integer ix,iy,iz ! voxel coords integer mi,mj,mk,ni,nj,nk ! coords for voxels above, below a node integer numsymm,numsamp,isymm real pi, pi2,degrad,tend,phi1,capphi,phi2 integer norients,numcsl,nsig,ngrain real quatsymm,quatsamp,quatcsl,aquat,grwt,tayf parameter(norients=100000) character*80 textitl character eulan*25,nomen,iper*6 integer numsigmas parameter ( numsigmas = 300 ) common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp ! common/a2/quatcsl(4,50),numcsl,nsig(50) common/a2/quatcsl(4,numsigmas),numcsl,nsig(numsigmas) common/a3/aquat(4,0:norients),ngrain,phi1(0:norients) 1 ,capphi(0:norients),phi2(0:norients) common/a4/ grwt(0:norients),tayf(0:norients),textitl,eulan,iper real quat(4),qresult(4),rodr(3),rresult(3),qsave(4) real qtwo(4,2),qtmp(4) integer ijk,iq,dotpos real tmp,cutoff real rnorm, xtmp,ytmp,ztmp integer date_time(8) character stamp*99 character tim*8,dat*9 character cell*3 real clr_red(-3:norients),clr_green(-3:norients) & ,clr_blue(-3:norients) integer clrrange , r13rd , r23rd , clr13 , clr23 character input_line*130,dumm*6 logical flag integer number,itmp,imx,imy,imz,dump(40) integer iq_colour ! choice of colour scheme ! for IVP colour method real hkl(3),dout(3) ! use for RGB colour real d001(3,3),d110(6,3),d111(4,3),vec(3),qr(4) real sqrt2,sqrt3,rsqrt2,rsqrt3,correct,drmax real scalar character inline*80,inp_name*100 integer output_type integer count ! no of args integer input_type logical found_grid integer max_size ! grid size integer threshold ! choice of spin number, below which no display integer kk,idum,matl,ioffx,ioffy,ioffz ! for FFT grids real p1,pp,p2 ! Euler angles integer jtmp,lagb,hagb,sigma3,sigma5,sigma7,sigma9 logical color_on integer sig_out , sig_val ! to return a sigma value real nearness ! and how near it is integer inbr,jnbr,knbr,nnspin real angle,theta integer csl_count(numsigmas) , csl_total ! for counting CSL fractions real fraction1 , fraction2 integer uxlimit,uylimit,uzlimit,fft_size ! for FFT structure size integer iqheader , iqphasenum , igasnumber ! CODE:: call newdate(dat,tim) stamp = 'Written by voxel2avs Date: '//DAT//' '//TIM count = iargc() if(count.le.0) then print*, 1 'Usage = ./Voxel2avs data output control texture-list scheme' print* print*,' where DATA contains the structure to be drawn; ' 2 ,' OUTPUT is a choice as follows: ' print*,' 0 = write AVS file ' print*,' 1 = write VTK file' print*,' 2 = write unstructured grid VTK file ' & ,'which is required for volume rendering' print*,' 3 = output FFT structure' print*,' 4 = output DX file (for Sukbin meshing)' print* 2 ,'CONTROL contains at least one character, as follows' print*,'0|1|2|3|4 means use texture as follows:' print*,' 0 => just use spin numbers' print*,' 1 => use info from texin' print*,' 2 => use info from texin1 and texin2' print*,' 3 => use info from Sukbin type list' print*,' 4 => use .dx as microstructure input and Sukbin list' print*,' 5 => read input from Groeber files' print* print*,'texture-list = file of orients' print*,' - needed for Sukbin/Groeber' print* print*,'SCHEME means colour according to:' print*,' Rodrigues-Frank vectors =1, ' print*,'or IVP with sample 3 (ND) =2,' print*,'or IVP with sample 2 (TD) =3,' print*,'or IVP with sample 1 (RD) =4' print*,'NB: you may have to UNcheck "Map Scalars" in paraview' print*,' in order to see the color-by-orientation' print* print*,'NB: spin numbers below zero will not be plotted' print*,'Adjust THRESHOLD to alter this behavior' print*,'But, you may have to use Threshold in Paraview' print* endif threshold = 1 ! smallest spin value plotted CALL GetArg(1, fname) if(fname.eq.'') then print *,' No file type given as argument:' print*,'Enter filename containing data' read(*,*) fname endif CALL GetArg(2, inline) if(inline.eq.'') then print* print*,'output is a choice as follows ' print*,' 0 = write AVS file ' print*,' 1 = write VTK file' print*,' 2 = write unstructured grid VTK file ' print*,' 3 = write FFT .txt file ' print*,' 4 = write .dx file ' read*,output_type else read(inline,"(i2)") output_type endif print*,'output type = ',output_type CALL GetArg(3, inline ) if(inline .eq. '') then print* print*,'Color by number [=0] or color by orientation; ' print*,'Orientations in texin [1] or texin1+texin2 [2] ?' print*,' or in a Sukbin-type texture list [3] ?' print*,' or .dx as microstructure input and Sukbin list [4]' print*,' or Groeber lists [5]' read * , iq else ! extract control parameters from input iq = -1 jtmp = index(inline,'0') if(jtmp.gt.0) iq = 0 jtmp = index(inline,'1') if(jtmp.gt.0) iq = 1 jtmp = index(inline,'2') if(jtmp.gt.0) iq = 2 jtmp = index(inline,'3') if(jtmp.gt.0) iq = 3 jtmp = index(inline,'4') if(jtmp.gt.0) iq = 4 jtmp = index(inline,'5') if(jtmp.gt.0) iq = 5 if(iq.lt.0) stop 'bad value of IQ' endif if( iq.lt.0 .or. iq.gt.5 ) then print*,' not an option for the input option' call exit(1) endif if( iq.eq.0 ) print*,'color type by spin number ' if(iq.eq.1) print*,'color type from texin' if(iq.eq.2) print*,'color type from texin1 and texin2' if(iq.eq.3) print*,'color type from Sukbin orient. list' if(iq.eq.4) print*,'.DX structure and Sukbin orient. list' if(iq.eq.5) print*,'Groeber structure and orient. list' pi = 4.*atan(1.) pi2 = 2.*atan(1.) degrad = 180./pi tend = sqrt(2.) - 1. cutoff = 1.0 ! criterion for new orientation sqrt2=sqrt(2.0) sqrt3=sqrt(3.0) rsqrt2=1./sqrt2 rsqrt3=1./sqrt(3.0) d001(1,1)=0. d001(1,2)=0. d001(1,3)=1. d001(2,1)=1. d001(2,2)=0. d001(2,3)=0. d001(3,1)=0. d001(3,2)=1. d001(3,3)=0. d110(1,1)=0. d110(1,2)=rsqrt2 d110(1,3)=rsqrt2 d110(2,1)=rsqrt2 d110(2,2)=0 d110(2,3)=rsqrt2 d110(3,1)=rsqrt2 d110(3,2)=rsqrt2 d110(3,3)=0. d110(4,1)=0. d110(4,2)=-1.*rsqrt2 d110(4,3)=rsqrt2 d110(5,1)=-1.*rsqrt2 d110(5,2)=0 d110(5,3)=rsqrt2 d110(6,1)=rsqrt2 d110(6,2)=-1.*rsqrt2 d110(6,3)=0. d111(1,1)=rsqrt3 d111(1,2)=rsqrt3 d111(1,3)=rsqrt3 d111(2,1)=-1.*rsqrt3 d111(2,2)=rsqrt3 d111(2,3)=rsqrt3 d111(3,1)=rsqrt3 d111(3,2)=-1.*rsqrt3 d111(3,3)=rsqrt3 d111(4,1)=-1.*rsqrt3 d111(4,2)=-1.*rsqrt3 d111(4,3)=rsqrt3 call prep(nnbors) ! set up NN array e_type = 'hex' ! print*,'Enter filename containing data' ! read(*,*) fname open(17,file=fname,status='old') if(iq.eq.4) then ! this for .DX file found_grid = .true. do while(found_grid) read(17,"(a)") input_line ! print*,'checking input:',input_line jtmp = 0 jtmp = index(input_line,'counts') if(jtmp.gt.0) then read(input_line(jtmp+6:jtmp+60),*) n1,n2,n3 c$$$ xx = n1-1 c$$$ yy = n2-1 c$$$ zz = n3-1 zz = n1-1 yy = n2-1 xx = n3-1 ! changed the order 5 Feb 08 to get consistency ! between the various descriptions found_grid = .false. endif enddo ! read a few lines to get the grid size else ! it's a PH file read(17,*) xx,yy,zz ! also works for Groeber file endif print*,'Dimensions of grid = ',xx,yy,zz max_size = max0(xx,yy,zz) ! if (( xx .gt. size ) .or. ( yy .gt. size ) ! & .or. ( zz .gt. size )) then if (( xx .gt. xsize ) .or. ( yy .gt. ysize ) & .or. ( zz .gt. zsize )) then stop 'too big!' endif imx = xx imy = yy imz = zz if( iq.eq.4 ) then ! this for .DX file found_grid = .true. do while(found_grid) read(17,"(a)") input_line ! print*,'checking input:',input_line jtmp = 0 jtmp = index(input_line,'3 class') if(jtmp.gt.0) found_grid = .false. enddo elseif( iq.eq.0 .or. iq.eq.1 .or. iq.eq.2 .or. iq.eq.3 ) then ! reading PH file ! print*,'reading line 2 ' read(17,"(a)") input_line ! print*,'checking input:',input_line ! dummy read; had keyword, time, , temperature read(17,*) dum1,dum2,dum3 ! dummy read; had radius,linevol,nucvol ! print*, '3rd line: ', dum1,dum2,dum3 ! dummy read; had radius,linevol,nucvol endif ! no extra reads needed for Groeber file spin_min = 9999999 spin_max = -9999 flag = .true. isite = 0 jsite = 1 ksite = 1 icount = 0 spin_max = 0 input_type = 0 ! 0 means basic list directed; 1 means reading one line at a time if(input_type.eq.0) then ! if(iq.eq.4) then ! read(17,*) (((spins(isite,jsite,ksite),ksite=1,zz) ! & ,jsite=1,yy),isite=1,xx) ! else read(17,*) (((spins(isite,jsite,ksite),isite=1,xx) & ,jsite=1,yy),ksite=1,zz) ! endif ! 5 Feb 08 - the DX files have Z varying fastest do isite = 1,xx do jsite = 1,yy do ksite = 1,zz if ( iq.eq.5 ) then ! Groeber, add 1 spins(isite,jsite,ksite) $ = spins(isite,jsite,ksite) + 1 endif if(spins(isite,jsite,ksite).lt.spin_min) then spin_min = spins(isite,jsite,ksite) endif if(spins(isite,jsite,ksite).gt.spin_max) then spin_max = spins(isite,jsite,ksite) endif enddo enddo enddo c 992 format(10i6/) 992 format(10i6) elseif(input_type.eq.1) then do while(flag) if(mod(icount,1000).eq.0) then print*,'reading line ',icount print*,'x,y,z: ',isite,jsite,ksite print*,'max spin = ', spin_max endif read(17,"(a)") input_line ! print*,input_line icount = icount + 1 itmp = len_trim(input_line) ! print* ! print*,'length of input = ',itmp number = itmp/6 ! print*,number,' numbers found in line',icount read(input_line,*) (dump(i), i = 1,number) ! print*,' tmp array: ',(dump(i), i = 1,number) c tells us number of i6 numbers available in the line do i = 1,number isite = isite + 1 if(isite.gt.imx) then isite = 1 jsite = jsite + 1 if(jsite.gt.imy) then jsite = 1 ksite = ksite + 1 endif endif spins(isite,jsite,ksite) = dump(i) if(dump(i).gt.spin_max) spin_max = dump(i) if(dump(i).lt.spin_min) spin_min = dump(i) c$$$ dumm = input_line(((i-1)*6+1):(i*6)) c$$$ read(dumm,"(i6)") spins(isite,jsite,ksite) ! print*, 'i,j,k,spin',isite,jsite,ksite,spins(isite,jsite,ksite) if(isite.ge.imx .and. jsite.ge.imy .and. ksite.ge.imz) & flag = .false. enddo ! do i = 1,number c pause enddo ! do while(flag endif ! if(input_type.eq.0) close(17) print*,'max. spin = ',spin_max print*,'min. spin = ',spin_min if( spin_min.lt.0 ) then print*,' negative spins will be ignored' endif if( spin_min .eq. 0 ) then print*,' min_spin = 0; adding 1 to all spins ' do ksite = 1 , zz do jsite = 1 , yy do isite = 1 , xx spins(isite,jsite,ksite) = spins(isite,jsite,ksite) + 1 end do end do end do endif if ( spin_min .eq. 0 ) then print*,'spin_min=0; adjusting spin_min, spin_max' spin_min = 1 spin_max = spin_max + 1 print*,'spin_min= ',spin_min,', spin_max= ',spin_max end if ! was below, after reading in orientations c$$$ print* c$$$ print*,'Color by number [=0] or color by orientation; ' c$$$ print*,'Orientations in texin [1] or texin1+texin2 [2] ?' c$$$ read*,iq if( iq.lt.0 .or. iq.gt.5 ) then print*,iq,' is not a valid value for control parameter' call exit(1) endif print*,'MAIN: reading in symmetry operators' if( iq.ge.1 .and. iq.le.5 ) then if( iq.eq.3 .or. iq.eq.4 .or. iq.eq.5 ) then CALL GetArg(4, inp_name) if(inp_name.eq.'') then print* print*,'Enter the filename for the texture list' read"(a)",inp_name ! else ! read(inline,"(a)") inp_name end if print*,'reading orientations from: ',inp_name end if ! if(iq.eq.3 call textur3( iq , spin_min , spin_max , inp_name ) ! call textur3(iq,spin_max,inp_name) ! call textur3(iq,spin_max) print*,'MAIN: Found ',numsymm,' xtal symms',' and ' print*,numsamp,' sample symms' CALL GetArg(5, inline) ! print*,'inline: ',inline if(inline.eq.'') then print* print*,'Colour scheme? RF-color=1,' print*,'or IVP with sample 3 (ND) =2,' print*,'or IVP with sample 2 (TD) =3,' print*,'or IVP with sample 1 (RD) =4' read*,iq_colour else read(inline,"(i1)") iq_colour end if if(iq_colour.lt.1 .or. iq_colour.gt.4) then print*,iq_colour,' = invalid colour scheme choice' call exit(1) end if if(iq_colour.gt.1) then jtmp = 5 - iq_colour hkl(1)=0. hkl(2)=0. hkl(3)=0. hkl(jtmp) = 1. ! set the desired sample axis for making the IVP map print*,'check HKL values: ',hkl end if end if if( iq.ge.1 .and. iq.le.5 ) then ! now we figure out colors according to values of the Rodrigues parameters range = 2.*tend + 1.e-4 do i = 1 , spin_max c$$$ print*,'no. of orients in av = ',orients_count(i) do j=1,4 c qtmp(j) = orients(i,j) qtmp(j) = aquat(j,i) enddo !DEBUG: print*,'#, orient. = ' , i , ( qtmp(j) , j = 1 , 4 ) rnorm=sqrt(qtmp(1)**2+qtmp(2)**2+qtmp(3)**2+qtmp(4)**2) if(abs(rnorm-1.0).gt.1.e-2) then write(*,"('MAIN: i, rnorm,qtmp ',i5,6(1x,f10.3))") 1 i,rnorm,qtmp endif ! ^^^^^^^^^^^^^^^^^^^ ! colour by RF vector values if( iq_colour.eq.1 ) then do j = 1,numsymm c do it in quaternions call postsymm(qtmp,j,qresult) call q2rod(rodr,qresult) if(abs(rodr(1)).le.tend.and. & abs(rodr(2)).le.tend.and. & abs(rodr(3)).le.tend.and. & (abs(rodr(1))+abs(rodr(2))+abs(rodr(3))).le.1.) then rresult(1)=rodr(1) rresult(2)=rodr(2) rresult(3)=rodr(3) isymm=j clr_red(i) = (rresult(1)+tend) / range clr_green(i) = (rresult(2)+tend) / range clr_blue(i) = (rresult(3)+tend) / range endif enddo ! j = 1, numsymm c call q2rod(rresult,qtmp) c$$$ print*,'symm_# R G B = ',isymm c$$$ 1 ,clr_red(i),clr_green(i),clr_blue(i) endif ! ^^^^^^^^^^^^^^^^^^^ ! end of colour by RF vector values ! &&&&&&&&&&&&&&&&&&& ! colour by inv. PF map if(iq_colour.ge.2 .and. iq_colour.le.4) then do j = 1,3 qr(j) = qtmp(j) ! reverse the rotation enddo qr(4) = -1.*qtmp(4) ! reverse the rotation call qrvec(qr,hkl,dout) ! write(*,"('%%INVPFMAP: DOUT= ',3(1x,f6.3))") ! $ (dout(ijk),ijk=1,3) c red=0. do j=1,3 do k=1,3 vec(k) = d001(j,k) enddo drmax = scalar(dout,vec) ! write(*,"('%%INVPFMAP: RED/drmax= ',3(1x,f6.3))") drmax if(abs(drmax).gt.red) then red = abs(drmax) if(red.gt.1.0) red = 1.0 drmax = acos(red) red = amax1(0.0,1.0-(drmax/1.1)) endif enddo c green = 0. do j = 1,6 do k = 1,3 vec(k) = d110(j,k) enddo drmax = scalar(dout,vec) ! write(*,"('%%INVPFMAP: GREEN/drmax= ',3(1x,f6.3))") drmax if(abs(drmax).gt.green) then green = abs(drmax) if(green.gt.1.0) green = 1.0 drmax = acos(green) green = amax1(0.0,1.0-(drmax/0.9)) endif enddo blue = 0. do j = 1,4 do k = 1,3 vec(k) = d111(j,k) enddo drmax = scalar(dout,vec) ! write(*,"('%%INVPFMAP: BLUE/drmax= ',3(1x,f6.3))") drmax if(abs(drmax).gt.blue) then blue = abs(drmax) if(blue.gt.1.0) blue = 1.0 drmax = acos(blue) blue = amax1(0.0,1.0-(drmax/0.9)) endif enddo c do j = 1,4 c drmax=red+green+blue drmax = amax1(red,green,blue) correct = 1./drmax ! write(*,"('%%INVPFMAP: j,correct= ' ! $,i6,3(1x,f6.3))") j,correct red = red*correct green = green*correct blue = blue*correct c enddo ! very crude way to get L1 norm=1 enddo ! very crude way to get largest value=1 ! print*,'i,R,G,B= ',i,red,green,blue c$$$ clr_red(i) = 255. * red c$$$ clr_green(i) = 255. * green c$$$ clr_blue(i) = 255. * blue clr_red(i) = red clr_green(i) = green clr_blue(i) = blue ! remember to write out opacity also endif ! &&&&&&&&&&&&&&&&&&& ! end of colour by inv. PF map ! print*,i,clr_red(i),clr_green(i),clr_blue(i) enddo ! do i=1,spin_max endif ! if(iq.ge.0 .or. iq.lt.4) if ( iq .eq. 0 ) then clrrange = spin_max - spin_min clr13 = clrrange / 3 clr23 = 2 * ( clrrange / 3 ) r13rd = spin_min + clrrange / 3 r23rd = spin_min + clr23 print*,'min, r13rd, r23rd max' print*, spin_min, r13rd, r23rd , spin_max do i = spin_min , spin_max if ( i .ge. spin_min .and. i .le. r13rd ) then clr_red(i) = 1.-(abs(float(i-spin_min)) & / float( r13rd - spin_min )) clr_green(i) = 1.-(abs(float(i-r13rd)) & / float( r13rd - spin_min )) clr_blue(i) = 0. endif if ( i .gt. r13rd .and. i .le. r23rd ) then clr_red(i) = 0. clr_green(i) = 1.-(abs(float(i-r13rd)) & / float( r23rd - r13rd )) clr_blue(i) = 1.-(abs(float(i-r23rd)) & / float( r23rd - r13rd )) endif if ( i .gt. r23rd .and. i .le. spin_max ) then clr_red(i) = 1.-(abs(float(i-spin_max)) & / float( spin_max - r23rd)) clr_green(i) = 0. clr_blue(i) = 1.-(abs(float(i-r23rd)) & / float( spin_max - r23rd)) endif print*,'i,red,green,blue' print*,i,clr_red(i), clr_green(i), clr_blue(i) c$$$ clr_red(i) = 255. * clr_red(i) c$$$ clr_green(i) = 255. * clr_green(i) c$$$ clr_blue(i) = 255. * clr_blue(i) ! color scalars enddo endif ! pause ! +++++++++++++++++++++++++++++++++++++++++++++++++ ! AVS output if(output_type.eq.0) then dotpos = 0 dotpos = index(fname,'.') if(dotpos.eq.0) then dotpos = index(fname,' ') endif ! open(16,file='micro.avs',status='unknown') open(16,file=fname(1:dotpos-1)//'.avs',status='unknown') npoints = (xx+1)*(yy+1)*(zz+1) nelemavs = xx*yy*zz nvalues = 0 ! no atttributes needed for nodes nvalues_e = 3 ! three colour values, or orientation values for elements write(16,*) npoints,nelemavs,nvalues,nvalues_e icount = 0 do i = 0,xx xcoord = float(i)/float(max_size) ix = min0(i+1,xx) do j = 0,yy ycoord = float(j)/float(max_size) iy = min0(j+1,yy) do k = 0,zz icount = icount+1 zcoord = float(k)/float(max_size) iz = min0(k+1,zz) n_mtl = spins(ix,iy,iz) ! print*,'ix,iy,iz,spins= ',ix,iy,iz,n_mtl if(i.eq.0.or.i.eq.xx.or.j.eq.0.or. & j.eq.yy.or.k.eq.0.or.k.eq.zz) then n_type = 10 else n_type = 0 do mi = -1,1,2 ni = mi * (-1) do mj = -1,1,2 nj = mj * (-1) c gives the opposite corner for the x and y coords c do mk = -1,1,1 c print*,'lower voxel:',i+max0(0,mi),j+max0(0,mj),k-1 c print*,'spin:',spins(i+max0(0,mi),j+max0(0,mj),k-1) c print*,'upper voxel:',i+max0(0,ni),j+max0(0,nj),k c print*,'spin:',spins(i+max0(0,ni),j+max0(0,nj),k) i1 = i+max0(0,mi) j1 = j+max0(0,mj) c k1 = k-1 k1 = k i2 = i+max0(0,ni) j2 = j+max0(0,nj) ! print*,'i1,j1,k1,i2,j2,k',i1,j1,k1,i2,j2,k ! if(spins(i1,j1,k1).ne.spins(i2,j2,k)) then ! print*,'i1,j1,k1,i2,j2,k',i1,j1,k1,i2,j2,k+1 if(spins(i1,j1,k1).ne.spins(i2,j2,k+1)) then n_type = 2 ! to get the code to run, I had to change the indexing in the z direction ! this needs to be checked ..... endif ! if(spins(i+ etc enddo enddo c if any pair of spins separated by the node are dissimilar, c then set the node type to be "interface" endif ! if(i.eq.0 etc write(16,"(3(1x,i8),3(1x,f10.4))") & icount,n_type,n_mtl,xcoord,ycoord,zcoord ! write(*,"(3(1x,i8),3(1x,f10.4))") ! & icount,n_type,n_mtl,xcoord,ycoord,zcoord enddo enddo enddo c writes out the node coords. icount = 0 do i = 1,xx i1 = i-1 i2 = i do j = 1,yy j1 = j-1 j2 = j do k = 1,zz icount = icount+1 k1 = k k2 = k+1 n1 = (yy+1)*(zz+1)*i1 + (zz+1)*j1 + k1 n2 = (yy+1)*(zz+1)*i2 + (zz+1)*j1 + k1 n3 = (yy+1)*(zz+1)*i2 + (zz+1)*j2 + k1 n4 = (yy+1)*(zz+1)*i1 + (zz+1)*j2 + k1 n5 = (yy+1)*(zz+1)*i1 + (zz+1)*j1 + k2 n6 = (yy+1)*(zz+1)*i2 + (zz+1)*j1 + k2 n7 = (yy+1)*(zz+1)*i2 + (zz+1)*j2 + k2 n8 = (yy+1)*(zz+1)*i1 + (zz+1)*j2 + k2 write(16,"(i8.8,1x,i6,1x,a3,8(1x,i6))") & icount,spins(i,j,k),e_type,n1,n2,n3,n4,n5,n6,n7,n8 enddo enddo enddo c write out the element type c no need for a section for node values write(16,*) '# no node values required' c section for element values write(16,*) ' 3 1 1 1 1' ! again, 3 attributes for each element write(16,*) ' 1, real' write(16,*) ' 1, real' write(16,*) ' 1, real' c not quite sure what needs to be written out here icount = 0 irange = spin_max - spin_min irange2 = irange/2 range2 = 0.5*float(spin_max - spin_min) do i=1,xx do j = 1,yy do k = 1,zz icount = icount + 1 red = float(max0(0,spins(i,j,k)-irange2))/range2 blue = float(max0(irange2,int(spins(i,j,k))))/range2 green = float(iabs(spins(i,j,k)-irange2))/range2 c$$$ red = red * 255. c$$$ green = green * 255. c$$$ blue = blue * 255. write(16,"(i8.8,1x,i6,3(1x,f7.3))") & icount , spins(i,j,k) , red $ , green , blue , 1.0 c$$$ write(16,"(i8.8,1x,i6,3(1x,i4))") c$$$ & icount,spins(i,j,k) c$$$ & , int(red) , int(green) , int(blue) enddo enddo enddo endif ! if(output_type.eq.0 c +++++++++++++++++++++++++++++++++++++++++++++++ ! now for VTK output if(output_type.eq.1) then dotpos = 0 dotpos = index(fname,'.') if(dotpos.eq.0) then dotpos = index(fname,' ') endif open(3,file=fname(1:(dotpos-1))//'.vtk',status='unknown') c for output to VTK, paraview write(3,"('# vtk DataFile Version 2.0')") write(3,"(' data set from ',a,a)") fname(1:dotpos+4),stamp write(3,"('ASCII')") write(3,"('DATASET STRUCTURED_POINTS')") c write(3,"('DIMENSIONS ',3(2x,i8))") mx,my,num_z_values write(3,"('DIMENSIONS ',3(2x,i8))") xx,yy,zz write(3,"('ORIGIN 0.0 0.0 0.0')") write(3,"('SPACING 1 1 1')") write(3,"('POINT_DATA ',i12)")xx*yy*zz write(3,"('SCALARS myScalars int 1') ") write(3,"('LOOKUP_TABLE default')") c do k = 1,num_z_values do k = 1,zz do j = 1,yy do i = 1,xx ! VTK data types states that X varies fastest ! (C++ order, as opposed to Fortran ?) c$$$ if(point(i,j,k).eq.0) then c$$$ point(i,j,k) = 99999 c$$$ endif write(3,"(1x,i6,$)") spins(i,j,k) c write(3,*) point(i,j,k) c write(3,*) correspond(point(i,j,k)) enddo enddo enddo write(3,*) if( iq_colour.gt.0 ) then ! no point in writing this block if all we have is the spins! write(3,*) ! write(3,"('VECTORS Rodrigues float 3') ") ! write(3,"('VECTORS Rodrigues int 3') ") ! write(3,"('COLOR_SCALARS Rodrigues 4') ") if( iq_colour.eq.1 ) then write(3,"('COLOR_SCALARS Rodrigues 3') ") endif if( iq_colour.eq.4 ) then write(3,"('COLOR_SCALARS IVP-X 3') ") endif if( iq_colour.eq.3 ) then write(3,"('COLOR_SCALARS IVP-Y 3') ") endif if( iq_colour.eq.2 ) then write(3,"('COLOR_SCALARS IVP-Z 3') ") endif do k = 1,zz do j = 1,yy do i = 1,xx ! VTK data types states that X varies fastest ! (C++ order, as opposed to Fortran ?) ii = spins(i,j,k) write(3,"(3(1x,f8.1))") $ clr_red(ii) , clr_green(ii) $ , clr_blue(ii) end do end do end do write(3,*) end if close(3) endif ! if(output_type.eq.1 ! +++++++++++++++++++++++++++++++ if(output_type.eq.2) then print*,'writing VTK unstructured grid to' print*,fname(1:(dotpos-1))//'_UNSTRUC_GRID.vtk' dotpos = 0 dotpos = index(fname,'.') if(dotpos.eq.0) then dotpos = index(fname,' ') endif open(3,file=fname(1:(dotpos-1))//'_UNSTRUC_GRID.vtk' 1 ,status='unknown') c for output to VTK, paraview write(3,"('# vtk DataFile Version 2.0')") write(3,"(' data set from ',a,' ',a)") fname(1:dotpos+4),stamp write(3,"('ASCII')") write(3,"(/'DATASET UNSTRUCTURED_GRID')") npoints = (xx+1)*(yy+1)*(zz+1) write(3,"(/'POINTS ',i10,' float')") npoints c write(3,"('DIMENSIONS ',3(2x,i8))") mx,my,num_z_values icount = 0 do i = 0,xx xcoord = float(i)/float(max_size) ix = min0(i+1,xx) do j = 0,yy ycoord = float(j)/float(max_size) iy = min0(j+1,yy) do k = 0,zz icount = icount+1 zcoord = float(k)/float(max_size) iz = min0(k+1,zz) write(3,"(3(1x,f10.4))") & xcoord,ycoord,zcoord ! write(*,"(3(1x,i8),3(1x,f10.4))") ! & icount,n_type,n_mtl,xcoord,ycoord,zcoord enddo enddo enddo print*,'wrote ',icount,' points' c writes out the node coords. nelemavs = xx*yy*zz ! this is the default number of elements, obviously do i = 1,xx do j = 1,yy do k = 1,zz if(spins(i,j,k).lt.threshold) then nelemavs = nelemavs - 1 ! remove an edge element from the count endif enddo enddo enddo write(3,"(/'CELLS ',i10,2x,i10)") nelemavs,nelemavs*9 icount = 0 do i = 1,xx i1 = i-1 i2 = i do j = 1,yy j1 = j-1 j2 = j do k = 1,zz if(spins(i,j,k).ge.threshold) then ! converse condition from above icount = icount+1 k1 = k-1 k2 = k n1 = (yy+1)*(zz+1)*i1 + (zz+1)*j1 + k1 n2 = (yy+1)*(zz+1)*i2 + (zz+1)*j1 + k1 n3 = (yy+1)*(zz+1)*i2 + (zz+1)*j2 + k1 n4 = (yy+1)*(zz+1)*i1 + (zz+1)*j2 + k1 n5 = (yy+1)*(zz+1)*i1 + (zz+1)*j1 + k2 n6 = (yy+1)*(zz+1)*i2 + (zz+1)*j1 + k2 n7 = (yy+1)*(zz+1)*i2 + (zz+1)*j2 + k2 n8 = (yy+1)*(zz+1)*i1 + (zz+1)*j2 + k2 write(3,"(' 8 ',8(1x,i7))") & n1,n2,n3,n4,n5,n6,n7,n8 ! only write out an element if it's in desired range endif ! (spins(i,j,k).ge.0) then enddo enddo enddo ! writes out the element nodes in the correct order for the hexahedron type (=12) cell = ' 12' write(3,"(/'CELL_TYPES ',i10)") nelemavs write(3,"(10(a))") (cell,i=1,nelemavs) c$$$ print*,'no of lines to print = ',(nelemavs/10)+1 c$$$ do i = 1,(nelemavs/10)+1 c$$$ write(3,"('12 12 12 12 12 12 12 12 12 12')") c$$$ enddo write(3,*) ! writes out the element type, which is 12 for a hexahedron ! a few extra entries is not a problem if(iq.ge.1 .and. iq.le.4) then ! use texture info write(3,"(/'POINT_DATA ',i10)") npoints write(3,"('COLOR_SCALARS colors 3')") do i = 0,xx i1 = min0(i+1,xx) do j = 0,yy j1 = min0(j+1,yy) do k = 0,zz k1 = min0(k+1,zz) ! coordinates of corresponding cell c$$$ red = float(max0(0,spins(i,j,k)-irange2))/range2 c$$$ blue = float(max0(irange2,spins(i,j,k)))/range2 c$$$ green = float(iabs(spins(i,j,k)-irange2))/range2 ! red = float(spins(i1,j1,k1)) / float(spin_max) ! green = 1.0 - (float(spins(i1,j1,k1)) / float(spin_max)) ! blue = 0. red = clr_red(spins(i1,j1,k1)) green = clr_green(spins(i1,j1,k1)) blue = clr_blue(spins(i1,j1,k1)) ! write(3,"(4(1x,f10.4))") ! & red,green,blue,0.5 write(3,"(4(1x,f10.4))") & red,green,blue ! we have to color the nodes which is weird, ! coloring the elements would make more sense enddo enddo enddo elseif(iq.eq.0) then ! or just spin numbers print* print*,'Note that you will have to re-set the color scale max' print*,' after reading the image into Paraview ' print*,' in order to see the range of colors that you want' print*,'Reminder, maximum spin value = ',spin_max write(3,"(/'POINT_DATA ',i10)") npoints write(3,"('SCALARS scalars int'/'LOOKUP_TABLE default')") ! write(3,"(i6)") (i,i = 1,npoints) do i = 0,xx i1 = min0(i+1,xx) do j = 0,yy j1 = min0(j+1,yy) do k = 0,zz k1 = min0(k+1,zz) ! coordinates of corresponding cell write(3,"(i6)") spins(i1,j1,k1) enddo enddo enddo endif ! if(iq.ge.1 .or. iq.le.4) endif ! if(output_type.eq.2 ! end of unstructured grid for VTK ! +++++++++++++++++++++++++++++++ if( output_type .eq. 3 ) then csl_total = 0 if( numsigmas .lt. 99 ) stop 'make number of CSLs >99!' do i = 1,numsigmas csl_count(i) = 0 enddo ! see below for CSL counts dotpos = 0 dotpos = index(fname,'.') if(dotpos.eq.0) then dotpos = index(fname,' ') endif print* print*,'++++++++++++++++' print*,'output type 3 / writing FFT structure to:' print*,fname(1:(dotpos-1))//'_FFT.txt' print*,'Also writing a .vtk file with CSL grain boundary info' print*,fname(1:(dotpos-1))//'_FFT_GBvoxels.vtk' print"(' and providing CSL analysis in ',$)" print*,fname(1:(dotpos-1))//'_CSL_analysis.txt' open(3,file=fname(1:(dotpos-1))//'_FFT.txt' 1 ,status='unknown') open(4,file=fname(1:(dotpos-1))//'_FFT_GBvoxels.vtk' 1 ,status='unknown') write(4,"('# vtk DataFile Version 2.0')") write(4,"(' data set from ',a,a)") fname(1:dotpos+4),stamp write(4,"('ASCII')") write(4,"('DATASET STRUCTURED_POINTS')") c write(3,"('DIMENSIONS ',3(2x,i8))") mx,my,num_z_values write(4,"('DIMENSIONS 128 128 128 ')") write(4,"('ORIGIN 0.0 0.0 0.0')") write(4,"('SPACING 1 1 1')") write(4,"('POINT_DATA ',i12)") 128*128*128 write(4,"('SCALARS Sigma_Values int 1') ") write(4,"('LOOKUP_TABLE default')") print* print*,'Some questions about the FFT file ...' print*,'Note that including a buffer layer is automatic' print*,'do you want the 2 dummy lines at the top? [1=yes]' read* , iqheader if ( iqheader.eq.1 ) then write(3,"(' 1')") write(3,"(' 0.00 0.00 0.00 0.00000 0')") ! dummy lines end if print*,'Make the Phase = Grain number ? [1=yes]' read* , iqphasenum print*,'What value for the gas/buffer ? [3 for Metz steels]' read* , igasnumber ! did this for the Metz steel specimens print*,'Grid size as x, y, z =',xx,yy,zz $ , ' , total points = ' , xx*yy*zz print*,'Enter dimensions of FFT grid (power of 2 !) :' read*,uxlimit,uylimit,uzlimit fft_size = uxlimit * uylimit * uzlimit print*,'You can have ',fft_size,', so what offset in x y z ?' read*,ioffx,ioffy,ioffz ! if( xx-ioffx .lt. 128 ) stop 'too much in x!' ! if( yy-ioffy .lt. 128 ) stop 'too much in y!' ! if( zz-ioffz .lt. 128 ) stop 'too much in z!' print*, 'x limits: ',ioffx+1 , min0( xx , ioffx+uxlimit ) print*, 'y limits: ',ioffy+1 , min0( yy , ioffy+uylimit ) print*, 'z limits: ',ioffz+1 , min0( zz , ioffz+uzlimit ) do k = ioffz+1 , ioffz + uzlimit if( mod(k,10) .eq. 0 ) print*,'z= ',k do j = ioffy+1 , ioffy + uylimit do i = ioffx+1 , ioffx + uxlimit ! if( mod(i,10) .eq. 0 ) print*,'x= ',i ! decide if real point or not ! as of 5 feb 08, X varies fastest for this type of VTK file idum = spins(i,j,k) if( i .gt. xx .or. j .gt.yy .or. k .gt. zz $ .or. idum .le. 0 ) then p1 = 0. pp = 0. p2 = 0. idum = spin_max + 2 matl = igasnumber else do kk = 1 , 4 quat(kk) = aquat(kk,spins(i,j,k)) enddo call q2eulB(p1,pp,p2,quat) p1 = p1 * degrad if ( p1 .lt. 0. ) p1 = p1 + 360. pp = pp * degrad ! if ( pp .lt. 0. ) pp = pp + 360. p2 = p2 * degrad if ( p2 .lt. 0. ) p2 = p2 + 360. ! idum = spins(i,j,k) matl = 1 if ( iqphasenum.eq.1 ) matl = idum ! did this for the Metz steel specimens c$$$ if( idum .le. 0 ) then c$$$ idum = spin_max+3 c$$$ matl = 2 c$$$ endif endif write(3,"(3(1x,f8.3),6(1x,i6))") & p1,pp,p2,(i-ioffx),(j-ioffy),(k-ioffz) & ,idum,matl ! last number is material number ! now we make the complementary grid for GB analysis if( i .gt. xx .or. j .gt.yy .or. k .gt. zz ) then sig_out = 0 else if( spins(i,j,k) .gt. 0 ) then ! guard against negative/zero spins qtwo(1,1)=aquat(1,spins(i,j,k)) qtwo(2,1)=aquat(2,spins(i,j,k)) qtwo(3,1)=aquat(3,spins(i,j,k)) qtwo(4,1)=aquat(4,spins(i,j,k)) else qtwo(1,1)= 0. qtwo(2,1)= 0. qtwo(3,1)= 0. qtwo(4,1)= 1. endif sig_out = 0 theta = 0. nearness = 2. ! set the values to defaults in case no GB found at this voxel ! do kk = 1,6 ! 1st NN only do kk = 1,nbors ! 26 NN call neighs(i,j,k,kk,inbr,jnbr,knbr) nnspin = spins(inbr,jnbr,knbr) if (spins(i,j,k) .ne. nnspin $ .and. nnspin .gt. 1 ) then qtwo(1,2) = aquat(1,nnspin) qtwo(2,2) = aquat(2,nnspin) qtwo(3,2) = aquat(3,nnspin) qtwo(4,2) = aquat(4,nnspin) call calc_gb(qtwo,theta,sig_val,nearness & ,color_on,red,green,blue & ,lagb,hagb,sigma3,sigma5,sigma7,sigma9) if( nearness .lt. 1.0 ) sig_out = sig_val endif enddo ! kk=1,6 if( theta .lt. 15. ) then sig_out = 1 ! flag low angle boundaries endif if( theta .lt. 0.9 ) then sig_out = 0 ! eliminate very low angle boundaries (probably not in most cleaned data sets endif if( theta .gt. 15. .and. nearness .gt. 1.0 ) then sig_out = 99 ! flag the general HAGB endif endif ! in range for the actual grid if( sig_out .gt. 0) then csl_total = csl_total + 1 csl_count(sig_out) = csl_count(sig_out) + 1 endif write (4,"(i3)") sig_out enddo enddo enddo close(3) close(4) open(5,file=fname(1:(dotpos-1))//'_CSL_analysis.txt' 1 ,status='unknown') write(5,"(' Sigma Fraction Fraction_no_LAGB ')") do i = 1 , 99 , 2 fraction1 = float(csl_count(i)) / float(csl_total) fraction2 = $ float(csl_count(i)) / float(csl_total - csl_count(1)) if( fraction1 .gt. 0. ) then write(5,*) i , fraction1 endif enddo close(5) print*,'completed FFT output: max. spin value= ',spin_max+2 endif ! if(output_type.eq.3 ! end of unstructured grid for VTK ! +++++++++++++++++++++++++++++++ if( output_type .eq. 4 ) then dotpos = 0 dotpos = index(fname,'.') if(dotpos.eq.0) then dotpos = index(fname,' ') endif print* print*,'++++++++++++++++' print*,'output type 4 / writing DX structure to:' print*,fname(1:(dotpos-1))//'_padded.dx' open(3,file=fname(1:(dotpos-1))//'_padded.dx' 1 ,status='unknown') write(3 , "('object 1 class gridpositions counts ',3i9)") & xx + 3 , yy + 3 , zz+3 write(3, "('origin 0 0 0')") write(3, "('delta 0 0 1')") write(3, "('delta 0 1 0')") write(3, "('delta 1 0 0')") write(3,*) write(3 , "('object 2 class gridconnections counts',3i9)") & xx + 3 , yy + 3 , zz+3 write(3,*) write(3,*) & "object 3 class array type float rank 0 items " & , ( (xx+2) * (yy+2) * (zz+2) ) , " data follows" do i = 0 , xx+1 do j = 0 , yy+1 do k = 0 , zz+1 ! pad the lattice with "-3" on all sides if( i .eq. 0 .or. i .eq. xx+1 & .or. j .eq. 0 .or. j .eq. yy+1 & .or. k .eq. 0 .or. k .eq. zz+1 ) then write(3,"(i3)") -3 else write(3,"(i6)") spins( i , j , k ) endif enddo enddo enddo write( 3 , * ) write( 3 , * ) 'attribute ','"dep"',' string ','"connections"' write( 3, * ) write( 3 , * ) 'object ','"3d voxel2avs-output"',' class field' write( 3 , * ) 'component ','"positions"',' value 1' write( 3 , * ) 'component ','"connections"',' value 2' write( 3 , * ) 'component ','"Spins"',' value 3' write(3 , * ) write(3 , * ) 'end' endif ! if( output_type .eq. 4 ) ! end of unstructured grid for padded DX output ! +++++++++++++++++++++++++++++++ call exit(0) 999 print*, 'could not find ',fname call exit(1) end ! ++++++++++++++++++++++++++++++++ C*****NEWDATE / borrowed from psplot.txt subroutine newdate(dat,tim) integer date_time(8) character*10 dummyarg character*(*)dat,tim character*36 monnam data monnam/'JANFEBMARAPRMAYJUNJULAUGSEPOCTNOVDEC'/ call date_and_time(dummyarg,dummyarg,dummyarg,date_time) c Convert new date/time to the old way write(dat(1:2),'(i2.2)')date_time(3) mind = (date_time(2)-1)*3+1 dat(4:6) = monnam(mind:mind+2) iyronly = date_time(1)-(date_time(1)/100)*100 write(dat(8:9),'(i2.2)')iyronly dat(3:3) = '-' dat(7:7) = '-' write(tim(1:2),'(i2.2)')date_time(5) write(tim(4:5),'(i2.2)')date_time(6) write(tim(7:8),'(i2.2)')date_time(7) tim(3:3) = ':' tim(6:6) = ':' return end c +++++++++++++++++++++++ subroutine prep(nnbors) integer nbors parameter (nbors = 26) integer nnbors(nbors,3) do i=1,nbors do j=1,3 nnbors(i,j)=0 enddo enddo nnbors(1,1)=1 nnbors(2,2)=1 nnbors(3,3)=1 nnbors(4,1)=-1 nnbors(5,2)=-1 nnbors(6,3)=-1 nnbors(7,1)=1 nnbors(7,2)=1 nnbors(8,1)=1 nnbors(8,3)=1 nnbors(9,2)=1 nnbors(9,3)=1 nnbors(10,1)=-1 nnbors(10,2)=1 nnbors(11,1)=-1 nnbors(11,3)=1 nnbors(12,1)=-1 nnbors(12,2)=-1 nnbors(13,1)=-1 nnbors(13,3)=-1 nnbors(14,2)=-1 nnbors(14,1)=1 nnbors(15,2)=-1 nnbors(15,3)=1 nnbors(16,2)=-1 nnbors(16,3)=-1 nnbors(17,3)=-1 nnbors(17,1)=1 nnbors(18,3)=-1 nnbors(18,2)=1 nnbors(19,1)=1 nnbors(19,2)=1 nnbors(19,3)=1 nnbors(20,1)=-1 nnbors(20,2)=1 nnbors(20,3)=1 nnbors(21,1)=1 nnbors(21,2)=-1 nnbors(21,3)=1 nnbors(22,1)=1 nnbors(22,2)=1 nnbors(22,3)=-1 nnbors(23,1)=-1 nnbors(23,2)=-1 nnbors(23,3)=1 nnbors(24,1)=-1 nnbors(24,2)=-1 nnbors(24,3)=-1 nnbors(25,1)=1 nnbors(25,2)=-1 nnbors(25,3)=-1 nnbors(26,1)=-1 nnbors(26,2)=1 nnbors(26,3)=-1 return end ! ++++++++++++++++++++++++++++++ c c NEIGHS : Finds site indices of the nearest neighbors of a site. c --------------------------------------------------------------- c This routine supplies periodic boundary conditions. It relies c on the matrix of increments for nn sites 'nnbors' (see LINEPREP). subroutine neighs(isite,jsite,ksite,nbor,inbr,jnbr,knbr) integer isite,jsite,ksite,nbor,inbr,jnbr,knbr integer size , nbors , xsize , ysize , zsize parameter (size = 390, nbors = 26 , xsize = 1024 $ , ysize = 512 , zsize = 32 ) integer nnbors(nbors,3) integer xx,yy,zz integer*2 spins( xsize , ysize , zsize ) common / structure / xx , yy , zz , nnbors , spins ! CODE:: inbr = isite+nnbors(nbor,1) jnbr = jsite+nnbors(nbor,2) knbr = ksite+nnbors(nbor,3) inbr = mod((inbr+xx-1),xx)+1 jnbr = mod((jnbr+yy-1),yy)+1 knbr = mod((knbr+zz-1),zz)+1 return end ! +++++++++++++++++++++++++++++++++++++++++++++ c LIKES : Finds the number of like neighbors of a c site with a given spin c -------------------------------------------------------------- c function likes(isite,jsite,ksite,spin) integer i,isite,jsite,ksite,spin,inbr,jnbr,knbr integer size , nbors , xsize , ysize , zsize parameter (size = 390, nbors = 26 , xsize = 1024 $ , ysize = 512 , zsize = 32 ) integer nnbors(nbors,3) integer xx,yy,zz integer*2 spins( xsize , ysize , zsize ) common / structure / xx , yy , zz , nnbors , spins likes = 0 do 10 i = 1,nbors call neighs(isite,jsite,ksite,i,inbr,jnbr,knbr) c$$$ print*,'inbr,jnbr,knbr spin ' c$$$ & ,inbr,jnbr,knbr,spins(inbr,jnbr,knbr) if(spins(inbr,jnbr,knbr).eq.spin) likes = likes+1 10 enddo return end !******************************************** subroutine calc_gb(qq,angle,sig_val,rnear & ,color_on,red,green,blue & ,lagb,hagb,sigma3,sigma5,sigma7,sigma9) implicit none integer size , nbors , xsize , ysize , zsize parameter (size = 390, nbors = 26 , xsize = 1024 $ , ysize = 512 , zsize = 32 ) integer nnbors(nbors,3) integer xx,yy,zz integer*2 spins( xsize , ysize , zsize ) common / structure / xx , yy , zz , nnbors , spins real qq(4,2),rnear,qresult(4),qr1(4),angle,qr2(4) real axis(3),crit,thetmin,pi,rad,brandon,theta integer i,ijk,jkl,sig_val,lsig integer index1,index2,index1a,index2a integer numsymm,numsamp,isymm integer norients,numcsl,nsig,ngrain real quatsymm,quatsamp,quatcsl,aquat,grwt,tayf integer numsigmas parameter ( numsigmas = 300 ) common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp common/a2/quatcsl(4,numsigmas),numcsl,nsig(numsigmas) logical color_on real red,green,blue integer lagb,hagb,sigma3,sigma5,sigma7,sigma9 ! CODE:: pi = 3.14159265 rad = 180./pi c$$$ print"('Quaternion 1 = ',4(2x,f8.3))" c$$$ 1 ,(qq(jkl,1),jkl=1,4) c$$$ print"('Quaternion 2 = ',4(2x,f8.3))" c$$$ 1 ,(qq(jkl,2),jkl=1,4) call comquat(qq,qr1) call disquat(qr1,index1,index2,qresult,index1a,index2a) ! print* ! write(*,"('from DISQUAT, qresult = ',4(1x,f8.3),2i4)") ! 1 qresult,index1a,index2a if( qresult(4) .gt. 1.0 ) qresult(4) = 1.0 if( qresult(4) .lt. -1.0 ) qresult(4) = -1.0 angle = 360./pi*acos(qresult(4)) ! write(*,"( 'from DISQUAT, Angle = ',f8.3,' degrees')") angle lsig = 0 rnear = 99. if(angle.lt.15.0) then lsig = 1 rnear = angle 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) thetmin = 66.0 if(lsig.eq.0) then do 150, 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 = nsig(ijk) rnear = crit endif 150 continue endif ! print"('indices: ',2i4,' angle= ',f8.3)",i,j,angle ! print"(' Axis= ',3(2x,f8.3))",(axis(ijk),ijk=1,3) ! print"('Quaternion= ',4(2x,f8.3))" ! 1 ,(qresult(jkl),jkl=1,4) ! print"('nearest Sigma:',i4,', distance = ',f8.3)" ! 1 ,lsig,rnear sig_val = lsig ! in v3, figure out colour and whether to plot in here c$$$ if(sig_val.gt.1 .and. rnear.lt.1.) then c$$$ print*,'sigma ',sig_val,', rnear = ',rnear c$$$ endif ! print*,'return disorientation angle in DEGREES' if(angle.lt.15.0 .and. lagb.gt.0) then !ASSIGN color to RED for LAGB color_on = .true. red = 1.0 green = 0 blue = 0 endif !ASSIGN color to BLACK for HAGB if(angle.ge.15.0 .and. hagb.gt.0) then color_on = .true. red = 0.1 green = 0.1 blue = 0.1 endif if(sigma3.gt.0 .and. sig_val.eq.3 & .and. rnear.lt.1.) then ! print*,'found sig-3!' color_on = .true. red = 1.0 green = 1.0 blue = 0. endif if(sigma5.gt.0 .and. sig_val.eq.5 & .and. rnear.lt.1.) then color_on = .true. red = 0.0 green = 1.0 blue = 0. endif if(sigma7.gt.0 .and. sig_val.eq.7 & .and. rnear.lt.1.) then color_on = .true. red = 0.0 green = 1.0 blue = 1. endif if(sigma9.gt.0 .and. sig_val.eq.9 & .and. rnear.lt.1.) then color_on = .true. red = 0.0 green = 0.0 blue = 1. endif return end 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) c c write(*,*) 'qresult ',qresult c return end c c _____________________________ c c subroutine disquat(qq,index1,index2,qresult,index1a,index2a) implicit none real qq(4),qqn(4),qresult(4),quint(4),quintn(4),qmax integer index1,index2,index1a,index2a,index3,iq4neg,i,j,k,ip,jp integer ijk,kjl real switch integer iqindex,iq1index real q1max,disor 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 index1=0 index2=0 index1a=0 index2a=0 call qnorm(qq) ! write(*,"('DISQUAT: qq normed = ',4f8.3)") qq if(abs(qq(4)).gt.0.999) then ! very low angle boundary, return zero misorientation ! print*,'almost zero or 180' qresult(1) = 0. qresult(2) = 0. qresult(3) = 0. qresult(4) = 1. index1 = 1 index2 = 1 index1a = 0 index2a = 0 return endif ! reproduce misquat code in order to estimate the angle qmax = 0. iqindex = 0 do i = 1,4 qq(i) = abs(qq(i)) if(qq(i).gt.qmax) then qmax = qq(i) iqindex = i endif enddo q1max = 0. iq1index = 0 c find the next highest q component do i=1,4 if(i.ne.iqindex) then if(qq(i).gt.q1max) then q1max = qq(i) iq1index = i endif endif enddo disor = amax1(qmax,(qmax+q1max)/sqrt(2.), & (qq(1)+qq(2)+qq(3)+qq(4))/2.) if(disor.gt.1.0) disor = 1.0 disor = disor - .05 if(disor.lt.-1.0) disor = -1.0 ! a small offset to avoid round off problems ! from here, the standard calculation qmax=0. ip=1 jp=1 do 200, i=1,24 call presymm(qq,i,quint) ! write(*,"('after presymm, quint = ',4f8.3)") quint do 190, j=1,24 c loop over the two sides of the grain call postsymm(quint,j,quintn) ! write(*,"('after postsymm, quint = ',4f8.3)") quintn if(abs(quintn(4)).gt.disor) then ! only go into the loops if the angle is known to be smaller ! than the previously determined minimum do 185, kjl=1,2 switch=1.0 if(kjl.eq.2) then switch=-1.0 endif do k=1,4 qqn(k)=quintn(k)*switch enddo c can take negative of Q because identity = +/-(0,0,0,1) do 210, ijk=1,2 qresult(1)=qqn(1) qresult(2)=qqn(2) qresult(3)=qqn(3) qresult(4)=qqn(4) iq4neg=ijk if(ijk.eq.2) qresult(4)=-1.*qqn(4) c take the inverse rotation, i.e. apply the switching symmetry ! print*,'i,j,kjl,ijk ',i,j,kjl,ijk ! write(*,"('DISQUAT: qresult = ',4f8.3)") qresult if(qresult(1).ge.0.0) then ! print*,'q1 > 0' c if((qresult(1)-qresult(2)).lt.1e-5) then c if((qresult(2)-qresult(3)).lt.1e-5) then c if((qresult(3)-qresult(4)).lt.1e-5) then if(qresult(1).le.qresult(2)) then ! print*,'q1 <= q2' if(qresult(2).le.qresult(3)) then ! print*,'q2 <= q3' if(qresult(3).le.qresult(4)) then ! print*,'q3 <= q4' c these IFs ensure that 0<=q1 qmax' c so now we test to see if we have a smaller angle qmax=qresult(4) index1=i index2=j c index3=0 index1a=0 if(ijk.eq.2) index1a=1 index2a=0 if(kjl.eq.2) index2a=1 c record the result if we have found a quat in the fund. zone c , and a smaller angle endif endif endif endif endif 210 enddo 185 enddo endif ! test for the angle 190 enddo ! j 200 enddo ! i if(index1.eq.0.or.index2.eq.0) then print*,'input quat.= ',qq stop 'DISQUAT failed!' endif c call presymm(qq,index1,quint) c call postsymm(quint,index2,qqn) ! write(*,"('DISQUAT: qqn = ',4f8.3)") qqn do 280, k=1,4 if(index2a.eq.1) then c if we had to take the negative of the Q, then do so c to the answer c qresult(k) = qqn(k) * (-1.0) c can take negative of Q because identity = +/-(0,0,0,1) else qresult(k) = qqn(k) endif 280 continue c if(index1a.eq.1) then qresult(4) = (-1.) * qresult(4) c restore the 4th component if disorient found for ijk=1 endif c c write(*,*) 'DISQUAT: qmax= ',qmax return end c c ____________________________ c subroutine cslquat(qq,numsig,theta,qresult) real qq(4),theta,qresult(4),tmp(2) real pi integer numsig ! index of CSL boundary - must be supplied integer numcsl,nsig real quatcsl integer numsigmas parameter ( numsigmas = 300 ) common/a2/quatcsl(4,numsigmas),numcsl,nsig(numsigmas) PI=3.14159265 c c algorithm for forming resultant quaternion c and determining angle based on fourth component c c note change of signs to get inverse of second orientation c qresult(1)=qq(1)*quatcsl(4,numsig)-qq(4)*quatcsl(1,numsig) & +qq(2)*quatcsl(3,numsig)-qq(3)*quatcsl(2,numsig) qresult(2)=qq(2)*quatcsl(4,numsig)-qq(4)*quatcsl(2,numsig) & +qq(3)*quatcsl(1,numsig)-qq(1)*quatcsl(3,numsig) qresult(3)=qq(3)*quatcsl(4,numsig)-qq(4)*quatcsl(3,numsig) & +qq(1)*quatcsl(2,numsig)-qq(2)*quatcsl(1,numsig) qresult(4)=qq(4)*quatcsl(4,numsig)+qq(1)*quatcsl(1,numsig) & +qq(2)*quatcsl(2,numsig)+qq(3)*quatcsl(3,numsig) c c write(*,*) 'qresult ',qresult if(qresult(4).gt.1.0) qresult(4)=1.0 if(qresult(4).lt.-1.0) qresult(4)=-1.0 theta=acos(qresult(4))*360./pi c write(*,*) 'theta= ',theta c return end c c _____________________________ c c subroutine presymm(qq,lindex,qresult) real qq(4),qresult(4) integer lindex c c include 'common.f' common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp c PI=3.14159265 c c algorithm for forming resultant quaternion c from applying a symmetry operator, QUATSYMM(n,lindex) c to the first quaternion/rotation c c If the symm oper == O, then c Qresult = (O x QQ) where QQ= Q1 x Q2^-1 c i.e. Qresult = (O x Q1) x Q2^-1 c note change of signs to get inverse of first orientation c 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) c c write(*,*) 'qresult ',qresult c return end c c _____________________________ c c subroutine postsymm(qq,lindex,qresult) real qq(4),qresult(4) integer lindex c include 'common.f' common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp 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 c$$$ print*,'POSTSYMM: lindex, symmop' c$$$ print"(i3,4(1x,f5.2))",lindex,(quatsymm(i,lindex),i=1,4) if(lindex.gt.numsymm .or. lindex.lt.1) then print*,'error in postsymm, lindex out of range' print*,'lindex, q ',lindex,qq call exit(1) endif ! this chunk of code came from one of the programs but is incorrect! 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) ! this is the correct code 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) rnorm=sqrt(qresult(1)**2+qresult(2)**2+qresult(3)**2 & +qresult(4)**2) if(abs(rnorm-1.0).gt.1.e-2) then write(*,*) 'qresult ',qresult,' index= ',lindex write(*,"('POSTSYMM: rnorm,quat ',6(1x,f10.3))")rnorm,qresult endif return end c ________________________________________ subroutine rod2q(d,quat) real d(3),quat(4),rtmp,rnorm,rnorm2 c converts Rodrigues vector to a quaternion c 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 subroutine qnorm(q) implicit none real q(4),norm integer i c CODE:: norm = q(1)**2+q(2)**2+q(3)**2+q(4)**2 do i = 1,4 q(i) = q(i) / norm enddo return end c c ______________________ c subroutine q2eulB(p1,p,p2,qq) c c converts quaternion to Bunge Euler angles c based on Altmann's solution for Euler->quat c real qq(4),p1,p,p2 real sum,diff c PI=3.14159265 c if((abs(qq(2)).lt.1e-35).and.(abs(qq(1)).lt.1e-35)) then diff=pi/4. else diff=atan2(qq(2),qq(1)) endif c ATAN2 is unhappy is both arguments are exactly zero! c if((abs(qq(3)).lt.1e-35).and.(abs(qq(4)).lt.1e-35)) then sum=pi/4. else sum=atan2(qq(3),qq(4)) endif c 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 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) real p1,p,p2,q(4) double precision c1,c,c2,s1,s,s2,d(3,3),tmp(4),sin2,cos2,rtmp,rnorm c CODE:: 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 ! ________________________________________ subroutine textur3( iq , spin_min , qvalue , inp_name ) implicit none ! second version, for use in WIRETEX.COLOUR ! not for use in REX3D ! updated Mar 09, to read in Groeber files (iq=5) ! updated Apr 09, to flip orientations from Groeber input integer spin_min integer iq,norients real quatsymm,quatsamp,quatcsl,aquat,grwt,tayf parameter(norients=100000) integer numsigmas parameter ( numsigmas = 300 ) common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp common/a2/quatcsl(4,numsigmas),numcsl,nsig(numsigmas) common/a3/aquat(4,0:norients),ngrain,phi1(0:norients) 1 ,capphi(0:norients),phi2(0:norients) common/a4/ grwt(0:norients),tayf(0:norients),textitl,eulan,iper integer numsymm,numcsl,nsig,ngrain,numsamp real phi1,capphi,phi2 ! include 'common.f' ! 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 integer qvalue,q2 ! max orient number, must be even integer i,j,k,ijk,ng real rad,d1,d2,d3,dummy character inp_name*100 real rtmp1,rtmp2,rtmp3 integer tmp1,tmp2,tmp3 real rot180a110(4) ! 180 degr about 110 data rot180a110 / 0.707107 , 0.707107 , 0. , 0. / real quats(4,2),qresult(4) ! for use with comquat ! CODE:: rad = 57.29578 q2 = qvalue/2 quats(1,2) = rot180a110(1) quats(2,2) = rot180a110(2) quats(3,2) = rot180a110(3) quats(4,2) = rot180a110(4) write(*,*) 'TEXTUR3: symmetry operators from "quat.symm.cubic"' open(unit=3,file ='quat.symm.cubic',status='old') read(3,*) ! skip over title line read(3,*) numsymm ! read number of symmetry operators if(numsymm.gt.24) then numsymm=24 write(*,*) 'too many symmetry operators, cutting back to 24!' else ! write(*,*) 'reading ',numsymm,' symmetry operators' endif do 1800, i=1,numsymm read(3,*) (quatsymm(ijk,i),ijk=1,4) ! write(*,*) 'symm. no. ',i,(quatsymm(ijk,i),ijk=1,4) 1800 continue close(unit=3) write(*,*) 1 'TEXTUR3: sample symmetry operators from "quat.symm.ort"' open(unit=3,file ='quat.symm.ort',status='old') read(3,*) ! skip over title line read(3,*) numsamp ! read number of sample symmetry operators if(numsamp.gt.4) then numsamp=4 ! write(*,*) 'too many symmetry operators, cutting back to 4!' else ! write(*,*) 'reading ',numsamp,' symmetry operators' endif do 1810, i=1,numsamp read(3,*) (quatsamp(ijk,i),ijk=1,4) ! write(*,*) 'samp. symm. no. ',i,(quatsamp(ijk,i),ijk=1,4) 1810 continue close(unit=3) write(*,*) 'TEXTUR3: CSL information from "quat.csl.data"' open(unit=3,file ='quat.csl.data',status='old') read(3,*) ! skip over title line read(3,*) numcsl ! read number of CSL types if(numcsl.gt.numsigmas) then numcsl=numsigmas ! 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) ! write(*,*) i,nsig(i),theta,lmn,r1,r2,r3,(quatcsl(ijk,i),ijk=1,4) 1900 continue close(unit=3) ! print*,'iq = ',iq if(iq.eq.2) then write(*,*) 'reading grain orientations from TEXIN1 and TEXIN2' elseif(iq.eq.1) then write(*,*) 'reading grain orientations from TEXIN' elseif(iq.eq.3) then write(*,*) ' orientations from Sukbin-type texture list' elseif(iq.eq.4) then write(*,*) ' structure from .DX and orientations & from Sukbin-type texture list' elseif( iq .eq. 5 ) then write(*,*) ' structure+orientations input from Groeber files' elseif(iq.le.0) then print*, 'iq LE 0' call exit(1) elseif( iq .gt. 5 ) then print*,'iq GT 5' call exit(1) endif do i = 1,qvalue do j = 1,4 aquat(j,i) = 0. enddo enddo do 529 ng = 1 , qvalue if(ng.eq.1.or.(iq.eq.2 .and. ng.eq.(q2+1))) then c adr, 22 vii 00; actually need independent orientations for all Q if(iq.eq.2) then if(ng.eq.1) open(unit=3,file ='texin1',status='old') if(ng.eq.(q2+1)) open(unit=3,file ='texin2',status='old') elseif(iq.eq.1) then ! texin source if(ng.eq.1) open(unit=3,file ='texin',status='old') elseif( iq .eq. 3 $ .or. iq .eq. 4 ! Sukbin-type source $ .or. iq .eq. 5 ) then ! Groeber-type source if(ng.eq.1) open(unit=3,file =inp_name,status='old') endif if(iq.eq.1 .or. iq.eq.2) then read(3,7) textitl 7 format(a) read(3,*) ! read(3,70) evm,fk1,nstate read( 3 , * ) 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(*,"('%%TEXTUR3: detected Bunge angles')") elseif(nomen.eq.'K') then write(*,"('%%TEXTUR3: detected Kocks angles')") endif c print 76,textitl 76 format(/ & ' TEXINx file= '/1x,a/) elseif( iq.eq.3 .or. iq.eq.4) then ! dealing with Sukbin input read(3,*) ! 1 header line read(3,*) ! discard 1st zero quat ! no header lines in Groeber files elseif( iq.eq.5 ) then ! dealing with Groeber input print*,'TEXTUR3:: reading Groeber texture list: ' print*,' if this fails, try not reading any header lines' read(3,*) ! 1 header line ! 1 header lines in Groeber files (but some of them have none!) endif endif if( iq .eq. 1 .or. iq .eq. 2 ) then read(3,*,end=530,err=530) d1,d2,d3 c$$$ read(3,*,end=530) d1,d2,d3,dummy,tayf(ng) c$$$ if(tayf(ng).eq.0.) tayf(ng)=1. c now let's read in the taylor factor for each orientation (skip wt.) ! print*,'no, Euler angles.. ',ng,d1,d2,d3 50 format(6f8.2,24f6.2) if(nomen.eq.'B') then phi1(ng)=d1 capphi(ng)=d2 phi2(ng)=d3 elseif(nomen.eq.'K') then phi1(ng)=d1+90. if(phi1(ng).gt.360.) phi1(ng)=phi1(ng)-360. capphi(ng)=d2 phi2(ng)=90.-d3 if(phi2(ng).lt.0.) phi2(ng)=phi2(ng)+360. endif 20 dth=d2 dps=d1 call quatB(phi1(ng)/rad,capphi(ng)/rad,phi2(ng)/rad,qr) ! use the new conversion routine based on Altmann if Bunge ! else, use longer routine do j = 1,4 aquat(j,ng)=qr(j) ! note use of ng as index end do ! print*,'no, quat.. ' , ng , (aquat(j,ng),j=1,4) elseif ( iq .eq. 3 .or. iq .eq. 4 ) then ! **** Sukbin lists read(3,*,end=530) tmp1,tmp2,tmp3,(qr(j),j=1,4) c$$$ print* c$$$ print"(4(i7),' q: ',4(1x,f8.3))" c$$$ & ,ng,tmp1,tmp2,tmp3,(qr(j),j=1,4) If( tmp2.lt.1 .or. tmp2.gt.qvalue ) then print*,'TEXTUR3: spin number out of range !' call exit(1) endif call qnorm(qr) ! 3 integers (sequence_no, spin_number, volume) then quaternion do j = 1,4 aquat(j,tmp2)=qr(j) ! note use of tmp2 as index end do elseif( iq .eq. 5 ) then ! ***** Groeber lists read( 3 , * , end=530 ) tmp1 , d1 , d2 , d3 if( spin_min.eq.0 ) tmp1 = tmp1 + 1 ! add 1 to match SPINS input ! grain (spin) number, 3 Eulers in radians if( tmp1 .le. 0 .or. tmp1 .gt. qvalue ) then ! note that Groeber structures start with spin = 0 print*,'TEXTUR3: spin number out of range !' call exit(1) endif call quatB( d1 , d2, d3 , qr ) quats(1,1) = qr(1) quats(2,1) = qr(2) quats(3,1) = qr(3) quats(4,1) = qr(4) call comquat(quats,qresult) ! flips orientation do j = 1,4 aquat( j , tmp1 ) = qresult( j ) ! note use of tmp1 as index end do ! print*,'no, quat.. ',tmp1,(aquat(j,tmp1),j=1,4) endif ! (iq.eq.3 .or. i ! print*,'no, quat.. ',tmp2,(aquat(j,tmp2),j=1,4) 529 continue 530 continue ngrain=ng-1 do i = 1,qvalue if(aquat(1,i).eq.0. .and. & aquat(2,i).eq.0. .and. & aquat(3,i).eq.0. .and. & aquat(4,i).eq.0.) then print"(i7,' q: ',4(1x,f8.3))",i,(aquat(j,i),j=1,4) print*,'zero quat found at i = ',i,' replacing' aquat(4,i) = 1.0 endif enddo return end c ________________________________________ c c subroutine textur4 c c second version, for use in REXGBS, WTS2POP c not for use in REX2D c integer numsymm,numsamp c common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp rad=57.29578 c write(*,*) 'reading symmetry operators from "quat.symm.cubic"' open(unit=3,file='quat.symm.cubic',status='old') read(3,*) ! skip over title line read(3,*) numsymm ! read number of symmetry operators if(numsymm.gt.24) then numsymm=24 write(*,*) 'too many symmetry operators, cutting back to 24!' else c write(*,*) 'reading ',numsymm,' symmetry operators' endif do 1800, i=1,numsymm read(3,*) (quatsymm(ijk,i),ijk=1,4) c write(*,*) 'symm. no. ',i,(quatsymm(ijk,i),ijk=1,4) 1800 continue close(unit=3) c write(*,*) 'reading sample symmetry operators from "quat.symm.ort"' open(unit=3,file='quat.symm.ort',status='old') read(3,*) ! skip over title line read(3,*) numsamp ! read number of sample symmetry operators if(numsamp.gt.4) then numsamp=4 write(*,*) 'too many symmetry operators, cutting back to 4!' else c write(*,*) 'reading ',numsamp,' symmetry operators' endif do 1810, i=1,numsamp read(3,*) (quatsamp(ijk,i),ijk=1,4) c write(*,*) 'samp. symm. no. ',i,(quatsamp(ijk,i),ijk=1,4) 1810 continue close(unit=3) return end c c _____________________________ c c subroutine misquat(qq,thetamin) real qq(4,2),thetamin,qresult(4),tmp(2),rquat(4) real disor,pi real qmax,q1max,q2max PI=3.14159265 c algorithm for forming resultant quaternion c and determining minimum angle taken from Sutton & Baluffi c note that the resultant quaternion is not returned c because it is not in the fundamental zone c note change of signs to get inverse of second orientation 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) c write(*,*) 'qresult ',qresult qmax=0. iqindex=0 do 10, i=1,4 qresult(i)=abs(qresult(i)) if(qresult(i).gt.qmax) then qmax=qresult(i) iqindex=i endif 10 continue q1max=0. iq1index=0 c find the next highest q component do 20, i=1,4 if(i.eq.iqindex) goto 20 if(qresult(i).gt.q1max) then q1max=qresult(i) iq1index=i endif 20 continue disor=amax1(qmax,(qmax+q1max)/sqrt(2.), & (qresult(1)+qresult(2)+qresult(3)+qresult(4))/2.) if(disor.gt.1.0) disor=1.0 if(disor.lt.-1.0) disor=-1.0 thetamin=acos(disor)*360./pi c write(*,*) 'thetamin ',thetamin q2max=0. iq2index=0 c find the next highest q component do 30, i=1,4 if(i.eq.iqindex.or.i.eq.iq1index) goto 30 if(qresult(i).gt.q2max) then q2max=qresult(i) iq2index=i endif 30 continue rquat(4)=qmax rquat(3)=q1max rquat(2)=q2max do 40, i=1,4 if(i.ne.iqindex.and. & i.ne.iq1index.and. & i.ne.iq2index) rquat(1)=qresult(i) 40 continue c supply the sorted quaternion c CAUTION: note that q1R2>R3 return end c c ____________________________ c c subroutine q2rod(d,quat) real d(3),quat(4),rtmp,rnorm,rnorm2 c converts a quaternion to a Rodrigues vector c rnorm=sqrt(quat(1)**2+quat(2)**2+quat(3)**2+quat(4)**2) if(rnorm.lt.1e-5) stop 'Q2ROD: too small a quaternion!' if(abs(rnorm-1.0).gt.1.e-2) then write(*,"('Q2ROD: quat to be normalized: ',6(1x,f10.3))") $ rnorm,quat 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! rtmp=sqrt(quat(1)**2+quat(2)**2+quat(3)**2) c this will be close to unity stheta=asin(rtmp) 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 endif return end c c _____________________________ c c subroutine presamp(qq,lindex,qresult) real qq(4),qresult(4) integer lindex c common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp c c algorithm for forming resultant quaternion c from applying a symmetry operator, QUATSYMM(n,lindex) c to the first quaternion/rotation c 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 with Qsymm c c Thus, for active rotations (standard definition of orientation in c mechanics) this is suitable for applying SAMPLE symmetry c 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) c c write(*,*) 'qresult ',qresult c 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 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 subroutine qrvec(qq,din,dout) c c uses quaternion to rotate a vector from DIN to DOUT; c note similarity to above routine c note that -q gives the same result as it should (represents same rotation) c real qq(4),din(3),dout(3) real t1 c t1=qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2 do 100, i=1,3 dout(i)=t1*din(i) c do 90, j=1,3 dout(i)=dout(i)+(2.*qq(i)*qq(j)*din(j)) c 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! dout(i)=dout(i)-(2.*qq(k)*qq(4)*float(kl)*din(j)) endif 90 continue 100 continue return end