program RF_misor ! /usr/local/bin/g77 -ffixed-line-length-none -g -o RF_misor RF_misor-[DATE].f ! on the G4, change the "include" to read as follows: ! include '/Volumes/7DEC05/code/paul.lee.codes/psplot.txt' ! on most machines: ! include '/code/paul.lee.codes/psplot.txt' ! this warning is not a problem: !/Volumes/7DEC05/code/paul.lee.codes/psplot.txt: In subroutine `axis': !RF_misor-25Jun09.f:896: warning: ! call keksym( 4.0 , ycharpos ! 1 !/Volumes/7DEC05/code/paul.lee.codes/psplot.txt:241: (continued): ! call keksym (xa,ya,sizttl,iscr,theta,nac,1) ! 2 !Argument #4 of `keksym' is one type at (2) but is some other type at (1) [info -f g77 M GLOBALS] ! with gfortran, there are many warnings from the PSPLOT subroutines C program for reading in an energy distribution c and plotting it out in triangular sections in RF space ! adr Feb 03; last edit apr 03; for use with *.mdf output from REX3D ! sep 09, many changes obviously; also added axes in bottom RH corner c this program may be freely copied for use by others c please acknowledge the NSF-supported Mesoscle Interface Mapping Project c based in the Materials Science & Engineering Dept c at Carnegie Mellon University, Pittsburgh c based on MAPE.F implicit none c MRESOLVE sets the resolution, NDATA is the number of exptl TJs c DELTAR is used to index CORRECT, DELTAE is used to index ENERGY integer m,m3,iqcolor,kcount,iq0,nval parameter (m = 115) double precision array(m,m,m) c ARRAY contains the data to be plotted double precision dens_array(m,m,m) c DENS_ARRAY contains the normalization values double precision rr1,rr2,rr3 integer mresolve,ndata,jndex,i,j,k,ii double precision tend,pi,pi2,degrad,deltae integer ndiscrete,itmp(3),fz_count double precision vmax,vmin,vavg character fname*60,fname2*60,fname3*60 double precision rtmp(3),rodr(3),rfsymm(3),rresult(3),weight,scale double precision quat(4),qresult(4) integer ir3val,ir2val,ir1val,type logical goodpoint,single,repeat character crys_sym_type,sam_sym_type integer numsymm,numsamp double precision quatsymm,quatsamp common/a1symm/ quatsymm(4,48),quatsamp(4,4),numsamp,numsymm integer index1,nbox,idelta,isymm,iqnum,iq_read_type double precision r3,rnorm,rnormmin real xsize,ysize,garray(m,m) integer idummy,jdummy double precision qtmp(4),qqsym(4),qint(4),delta,r1,r2 double precision voxvol(m,m,m),volsum,avgvol,weight_total integer badpoints,ijk,kk c double precision rmax,SMF(m,(-1*m):m) double precision wdth real width real xclip(100),yclip(100) integer nclip character*1 trichar(10) real colrval(3,0:10),cval(0:10) real greyval(0:10) real tx(20),ty(20),txdiff,tydiff,txold,tyold real txc(20),tyc(20) COMMON/CONPAR/ ISPEC, IOFFP, SPVAL, ILEGG, ILABB, NHII, & NDECCN, NLBLL, LSCAL, LDASH, HGTLAB integer ISPEC, IOFFP, ILEGG, ILABB, NHII, & NDECCN, NLBLL, LSCAL, LDASH real spval, hgtlab real rtmp9,xx1,xx2,yy1,yy2,height,r3_real integer iq_symm character arg_line*40,arg_line_contour*40 double precision tan60,sin60,cos60,tan30,sin30,cos30 logical test_fz_hcp,test_fz_cubic,inside logical test_fz_orth real x_margin,y_margin,RFe integer dotpos,lastchar logical add_mdf integer data_points ! how many lines read in double precision angle ! for outputting MDF values integer jcount real Euler1(3),Euler2(3) double precision quat1(4),quat2(4),quint(4) real xshift , yshift logical header integer count character*160 title integer iq_weight ! to control use of weight real temp ! temporary result integer iq_print,iq_2,iq_debug ! control variables real rlimit,rl2 ! for size of Hex clip box integer iq_30_rotate ! control rotation of plot data integer vertical_section real ycharpos c following variables control plotting options c CODE:: PI = 3.14159265 print*,'Welcome to RF_misor' print*,'please acknowledge the NSF-supported ' print*,'Mesoscale Interface Mapping Project under the MRSEC' print*,'based in the Materials Science & Engineering Dept.' print*,'at Carnegie Mellon University, Pittsburgh' print* print*,' Usage = ' 1 //'./RFmisor input_data smooth_width[1.5] ' 2 //' quats[0]/RFvectors[1]/' 3 //'pair_Euler_angles[2]/1_set_Eulers[3]' 4 //' cubic[0]/hcp[1]/orth[2]/hcp-rotated-30[-1]' 5 //' max_pts[99999999]' 6 //' weighted[0=NO] accept_contours[0] ' 7 //' print_debug_files [1] ' print* ilabb = 0 ! contour labels are NOT plotted nhii = -1 ! suppresses "H" and "L" marks CALL GetArg(1,fname2) if(fname2.eq.'') then c detect an argument for the input file name; look for user input if not present write(*,*) 'Enter the filename for the input data' write(*,*) 'An extension of .mdf is assumed' read(*,'(a)') fname2 endif c now we read in the data add_mdf = .true. c lastchar = index(fname2,' ') lastchar = len_trim(fname2) dotpos = lastchar + 1 if(fname2(lastchar-3:lastchar).eq.'.mdf') then add_mdf = .false. dotpos = lastchar-3 endif c dotpos = index(fname2,'.') ! find the period c if(dotpos.eq.0) then c add_mdf = .true. c dotpos = lastchar c endif fname = fname2(1:(dotpos-1)) single=.true. if(add_mdf) then inquire(file=fname2,exist=repeat) else inquire(file=fname(1:(dotpos-1))//'.mdf',exist=repeat) endif if(.not.repeat) then write(*,*) 'Did NOT find ',fname2 write(*,*) 'Will look for series of .mdf files' single=.false. else write(*,*) 'Found ',fname2 endif if(add_mdf) then open(unit=2,file = fname2,status='old') else open(unit=2,file = fname(1:(dotpos-1))//'.mdf',status='old') endif read(2,"(a)") title i = len(title) print"(a)",'Header: '//title(1:i-1) ! assume there is at least one header line header = .true. count = 1 do while(header) read(2,"(a)") title if(title(1:1).eq.'#') then count = count + 1 i = len(title) print"(a)",'Header: '//title(1:i-1) else header = .false. endif enddo rewind(2) do i = 1,count read(2,*) enddo ! read past the header lines CALL GetArg(2,arg_line) c print*,'debug: arg_line = ',arg_line if(arg_line.eq.'') then write(*,*) 'Enter a smoothing width (1.5 suggested)' read(*,*) wdth c wdth=6. else read(arg_line,*) wdth endif CALL GetArg(3,arg_line) if(arg_line.eq.'') then print*,'read in quats from rex3d [=0] ' print*,'or from RF (eg genranRF) [=1]?' print*,'or from list of GB segments (2 sets Eulers RADS) [=2]?' print*,'or from list of GB pixels (1 set Eulers RADS) [=3]?' read(*,*) iq_read_type else read(arg_line,*) iq_read_type endif if(iq_read_type.lt.0.or.iq_read_type.gt.3) then stop 'cannot read that type' endif iq_30_rotate = 0 CALL GetArg(4,arg_line) if(arg_line.eq.'') then print*,'cubic [= 0] or hexagonal [= 1] or orthorhombic [=2]?' print*,' for hex rotated by 30 (to go from Y- ' $ //'to X-convention, e.g., use -1' read(*,*) iq_symm else read(arg_line,*) iq_symm endif if(iq_symm.eq.-1) then iq_30_rotate = 1 iq_symm = 1 endif if(iq_symm.lt.0.or.iq_symm.gt.2) then stop 'unknown symmetry type' endif CALL GetArg(5,arg_line) if(arg_line.eq.'') then print*,' number of lines (approx. no. points) = ' call system('wc -l '//fname2) print*,'Maximum points ?' read(*,*) iqnum else read(arg_line,*) iqnum endif CALL GetArg(6,arg_line) if(arg_line.eq.'') then write(*,*) 'Enter 0 to ignore weights ' print*,' Or any other number to use the weight information' read(*,*) iq_weight else read(arg_line,*) iq_weight endif if( iq_weight.eq.0 ) then print*,'UN-weighted MDF will be plotted' fname2 = fname2(1:(dotpos-1))//'_NoW' fname = fname(1:(dotpos-1))//'_NoW' dotpos = dotpos + 4 endif CALL GetArg(7,arg_line_contour) iq_debug = 0 CALL GetArg(8,arg_line) if(arg_line.eq.'1') then iq_debug = 1 endif if( iq_debug.eq.1 ) then print*,'will write debug info to files' else print*,'no debug info will be written' endif tan60 = tan (pi*60./180.) sin60 = sin (pi*60./180.) cos60 = cos (pi*60./180.) tan30 = tan (pi*30./180.) sin30 = sin (pi*30./180.) cos30 = cos (pi*30./180.) m3 = m*m*m ! number of cells ispec=1 ! must be non-zero for other values to take effect ioffp=0 ilegg=0 ilabb=1 nhii=-1 ! suppresses "H" and "L" marks lscal=1 nval=7 ! number of contour values c xsize=2. c ysize=2. tend = sqrt(2.)-1. c defines edge of the R-F space pi = 4.*atan(1.) pi2 = 2.*atan(1.) degrad = 180./pi call textur4(iq_symm) !reads in symmetry operators write(*,*) 'Found ',numsymm,' xtal symms',' and ' write(*,*) numsamp,' sample symms' if(iq_symm.eq.0) then xsize = 4. ! change this value if you plot a different area in each section x_margin = 0.5 y_margin = 1.3 xshift = 0.5 * xsize yshift = 2.1 elseif(iq_symm.eq.1) then xsize = 7. x_margin = 0.5 y_margin = 1.2 xshift = 0.45 * xsize yshift = 1.9 elseif(iq_symm.eq.2) then xsize = 4. x_margin = 0.5 y_margin = 1.2 xshift = 0.75 * xsize yshift = 2.1 endif ysize = xsize txold=0. tyold=0. tx(1) = x_margin ty(1) = y_margin tx(2) = x_margin ty(2) = ty(1) + yshift tx(3) = x_margin ty(3) = ty(2) + yshift tx(4) = x_margin ty(4) = ty(3) + yshift if( iq_symm.eq.1 ) then ! hex tx(5) = x_margin ty(5) = ty(4) + yshift tx(6) = 2.*x_margin + xshift ty(6) = y_margin elseif( iq_symm.eq.2 ) then ! ortho tx(5) = x_margin + xshift ty(5) = y_margin tx(6) = x_margin + xshift ty(6) = ty(5) + yshift else ! cubic tx(5) = 2.*x_margin + xshift ty(5) = y_margin tx(6) = tx(5) ty(6) = ty(5) + yshift end if tx(7) = tx(6) ty(7) = ty(6) + yshift tx(8) = tx(6) ty(8) = ty(7) + yshift tx(9) = tx(6) ty(9) = ty(8) + yshift tx(10) = tx(6) ty(10) = ty(9) + yshift ! this one for axes tx(11) = 3.*x_margin + 2.*xshift ty(11) = y_margin if( iq_symm.eq.2 ) then ! ortho tx(11) = -0.2*x_margin + 2.*xshift end if if(iq_debug.eq.1) then do i = 1,8 print*,'i,tx,ty ',i,tx(i),ty(i) enddo c 7 triangles, ll corners endif trichar(1)='a' trichar(2)='b' trichar(3)='c' trichar(4)='d' trichar(5)='e' trichar(6)='f' trichar(7)='g' trichar(8)='h' trichar(9)='i' trichar(10)='j' c 10 triangles mresolve = 10 ! one page for now! kcount = 0 vmax = 0. vavg = 0. vmin = 999999. ycharpos = 0.8 ! for positioning messages at bottom of plot ! the maximum value of any coordinate is DELTAE * M / 2 ! M = 115, M/2 = 57 if(iq_symm.eq.0) then c deltae = 2.0*tend/float(m) deltae = 2.0 * tend / 100. RFe = tend elseif(iq_symm.eq.1) then c deltae = 2.0/float(m) deltae = 2.0 / 100. RFe = 1. c need more space for HCP; deltae = .02; Rmax = 57*.02 = 1.14 ! about right for the outer edge of the FZ in this case elseif(iq_symm.eq.2) then deltae = 1.9 / 100. RFe = 1. c need slightly smaller space for orthorhombic ! deltae = .019; Rmax = 57*.019 ~ 1.1 endif c using a larger array allows for the smoothing to be done over a region larger c than the fundametal zone (FZ), so no loss in total intensity should occur write(*,*) 'Initializing array(s)' do i=1,m c if(mod(i,25).eq.0) write(*,*) 'Row number ',i do j=1,m do k=1,m array(i,j,k)=0.0 enddo enddo enddo c this is inefficient for cubic systems, but does allow for less symmetric crystal symmetries crys_sym_type='C' sam_sym_type='O' c now we estimate the volume of each cell according to the volume element size fz_count=0 volsum=0.0 do i=1,m r1 = dfloat(i-(m/2))*deltae c reverse of this is c i = int(r1/deltae) + m/2 do j = 1,m r2 = dfloat(j-(m/2))*deltae do k = 1,m r3 = dfloat(k-(m/2))*deltae rnorm = dsqrt(r1**2+r2**2+r3**2) ! voxvol(i,j,k) = 1. / (1.+rnorm**2)**2 voxvol(i,j,k) = $ 1.d0 * (1.d0/(1.d0 + rnorm**2))**2 ! fixed 16Mar03 c approx. volume of cell if(iq_symm.eq.0) then inside = test_fz_cubic(r1,r2,r3) elseif(iq_symm.eq.1) then inside = test_fz_hcp(r1,r2,r3) elseif(iq_symm.eq.2) then inside = test_fz_orth(r1,r2,r3) endif if(inside) then volsum = volsum+voxvol(i,j,k) fz_count = fz_count+1 endif enddo enddo enddo avgvol = volsum/float(fz_count) print*,'Total volume in FZ = ',volsum print*,'Average volume of cell in FZ = ',avgvol print*,'Number of cells in FZ (fz_count) = ',fz_count if(iq_debug.eq.1) then if(add_mdf) then open(unit=19,file = fname2(1:(dotpos-1))//'_RFvalues.txt' $ ,status='unknown') else open(unit=19,file = fname(1:(dotpos-1))//'_RFvalues.txt' $ ,status='unknown') endif write(19,*) ' R1 R2 R3 weight ' if(add_mdf) then open(unit=20 $ ,file = fname2(1:(dotpos-1))//'_RFvalues_cells.txt' $ ,status='unknown') else open(unit=20 & ,file = fname(1:(dotpos-1))//'_RFvalues_cells.txt' $ ,status='unknown') endif write(20,*) ' R1 R2 R3 weight x y z ' endif c icount=0 badpoints=0 weight_total=0.0 jcount = 0 do i = 1,iqnum ! label 580 if(mod(i,5000).eq.0) write(*,*) 'reading line no. ',i c goodpoint=.true. rnormmin=999999. isymm=0 c each input line should contain 2 indices and 5 floats: i,j,q1,q2,q3,q4,weight c$$$ if(iq_read_type.eq.0) then c$$$ read(2,*,err=580,end=30) idummy,jdummy c$$$ $ ,(qtmp(ii),ii=1,4),weight if(iq_read_type.eq.0) then read( 2 , * , end=30 ) idummy,jdummy $ ,(qtmp(ii),ii=1,4),weight elseif(iq_read_type.eq.1) then c$$$ if(iq_weight .eq. 0) then c$$$ read(2,*,err=580,end=30) (rodr(ii),ii=1,3) if(iq_weight .eq. 0) then read( 2 , * , end=30 ) (rodr(ii),ii=1,3) else c$$$ read(2,*,err=580,end=30) (rodr(ii),ii=1,3),weight read( 2 , * , end=30 ) (rodr(ii),ii=1,3),weight endif call rod2q(rodr,qtmp) elseif(iq_read_type.eq.2) then if(iq_weight .eq. 0) then c$$$ read(2,*,err=580,end=30) (Euler1(ii),ii=1,3) c$$$ $ ,(Euler2(ii),ii=1,3) read( 2 , * , end=30 ) (Euler1(ii),ii=1,3) $ ,(Euler2(ii),ii=1,3) else c$$$ read(2,*,err=580,end=30) (Euler1(ii),ii=1,3) c$$$ $ ,(Euler2(ii),ii=1,3),weight read( 2 , * , end=30 ) (Euler1(ii),ii=1,3) $ ,(Euler2(ii),ii=1,3),weight endif call quatB(Euler1(1),Euler1(2),Euler1(3),quat1) call quatB(Euler2(1),Euler2(2),Euler2(3),quat2) call comquat2(quat1,quat2,quint) do ii = 1,4 qtmp(ii) = quint(ii) enddo ! copy to single precision elseif(iq_read_type.eq.3) then if(iq_weight .eq. 0) then c$$$ read(2,*,err=580,end=30) (Euler1(ii),ii=1,3) read( 2 , * , end=30 ) (Euler1(ii),ii=1,3) else c$$$ read(2,*,err=580,end=30) (Euler1(ii),ii=1,3),weight read( 2 , * , end=30 ) (Euler1(ii),ii=1,3),weight endif call quatB(Euler1(1),Euler1(2),Euler1(3),qtmp) endif ! if(iq_read_type.eq.0) 575 format(4f12.4) c print*,i,'Q ',qtmp,' weight ',weight c$$$ write(*,*) c$$$ write(*,*) 'input line number ',i c$$$ print*,'input: ',(qtmp(ii),ii=1,4) if( weight.le.0.0 ) weight = 0. if( iq_weight.eq.0 ) weight = 1. ! the weight/no weight option rnorm=sqrt(qtmp(1)**2+qtmp(2)**2+qtmp(3)**2+qtmp(4)**2) if(rnorm.lt.1e-5) stop 'input: too small a quaternion!' if(abs(rnorm-1.0).gt.1.e-3) then c write(*,*) 'Q2ROD: quat required normalization' do ijk=1,4 qtmp(ijk)=qtmp(ijk)/rnorm enddo endif c now apply symmetry elements c note that the RF vector is an active rotation from ref to xtal frame do j=1,numsymm ! call presymm(qtmp,j,qint) call postsymm(qtmp,j,qint) do k = 1,numsymm ! call postsymm(qint,k,qresult) call presymm(qint,k,qresult) do kk = 1,-1,-2 ! generate negative rotation if(kk.lt.0) then do ijk = 1,3 qresult(ijk)=qresult(ijk)*float(kk) enddo endif c write(*,"('symm1/2, qresult= ',2i5,4(1x,f8.5))") j,k,qresult goodpoint=.true. c if(abs(qresult(4)).gt.0.85) then c don't bother to include unless angle <63 degrees call q2rod(rodr,qresult) ! convert the quaternion to a RF vector rresult(1) = rodr(1) rresult(2) = rodr(2) rresult(3) = rodr(3) ! write(*,"('rresult = ',i4,i4,3(1x,f8.4))") k,kk,rresult ! **** these next lines rotate the results 30 degrees if desired so as ! to convert Y-convention results to X-convention, for example if(iq_30_rotate .eq. 1) then rresult(1) = (0.8660254 * rodr(1)) + (rodr(2) *0.5) rresult(2) = (rodr(2)*0.8660254) - ( rodr(1) *0.5) endif if(rresult(1).lt.5. .and. & rresult(2).lt.5. .and. & rresult(3).lt.5.) then if(iq_debug.eq.1) $ write(19,"(4(1x,f8.4))") rresult, weight endif c i = int(r1/deltae) + m/2 c ir3val=1+int((rresult(3)+RFe-(0.5*deltae))/deltae) c reverse of this is: c r3 = float(i3)*deltae - RFe + (0.5*deltae) ! ir3val = int(rresult(3)/deltae) + m/2 temp = ((rresult(3)/deltae)+0.5) if( temp.lt.0.) then ir3val = int(temp) - 1 + m/2 else ir3val = int(temp) + m/2 endif if(ir3val.lt.1.or.ir3val.gt.m) then c write(*,*) 'bad point 3 at line ',i goodpoint=.false. endif c ir2val=1+int((rresult(2)+RFe-(0.5*deltae))/deltae) ! ir2val = int(rresult(2)/deltae) + m/2 temp = ((rresult(2)/deltae)+0.5) if( temp.lt.0.) then ir2val = int(temp) - 1 + m/2 else ir2val = int(temp) + m/2 endif if(ir2val.lt.1.or.ir2val.gt.m) then c write(*,*) 'bad point 2 at line ',i goodpoint=.false. endif c ir1val=1+int((rresult(1)+RFe-(0.5*deltae))/deltae) ! ir1val = int(rresult(1)/deltae) + m/2 temp = ((rresult(1)/deltae)+0.5) if( temp.lt.0.) then ir1val = int(temp) - 1 + m/2 else ir1val = int(temp) + m/2 endif if(ir1val.lt.1.or.ir1val.gt.m) then c write(*,*) 'bad point 1 at line ',i goodpoint=.false. endif c write(*,*) 'indices into ARRAY: ', ir1val,ir2val,ir3val if(goodpoint) then kcount=kcount+1 array(ir1val,ir2val,ir3val) = $ array(ir1val,ir2val,ir3val)+weight ! write(*,*) 'indices into ARRAY: ', ir1val,ir2val,ir3val c$$$ write(*,"('symm1/2, qresult= ',2i5,4(1x,f8.5))") j,k,qresult c$$$ write(*,"('rresult = ',i4,i4,i4,3(1x,f8.4))") j,k,kk,rresult c$$$ write(*,*) 'indices into ARRAY: ', ir1val,ir2val,ir3val if(iq_symm.eq.0) then inside = $ test_fz_cubic(rresult(1),rresult(2),rresult(3)) elseif(iq_symm.eq.1) then inside = $ test_fz_hcp(rresult(1),rresult(2),rresult(3)) elseif(iq_symm.eq.2) then inside = $ test_fz_orth(rresult(1),rresult(2),rresult(3)) endif if(inside) then weight_total = weight_total + weight jcount = jcount + 1 endif ! write(*,"('rresult = ',i4,i4,i4,3(1x,f8.4))") j,k,kk,rresult if(iq_debug.eq.1) write(20,"(4(1x,f8.4),3(i5),l6)") $ rresult, weight,ir1val,ir2val,ir3val,inside else badpoints=badpoints+1 c write(*,*) 'bad point 1 at line ',i c write(*,*) 'input line+symmetry:', (rresult(ii),ii=1,3),weight c write(*,*) 'symmetry element = ',isymm,' RNORM =',sqrt(rnorm) c write(*,*) 'indices into ARRAY: ', ir1val,ir2val,ir3val c rnorm=sqrt(rtmp(1)**2+rtmp(2)**2+rtmp(3)**2) c write(*,*) 'RTMP magnitude = ',rnorm endif c endif ! from qresult(4) test c don't use now that HCP is included enddo ! inversion enddo ! 2nd crystal symmetry loop enddo ! 1st xtal symm 580 enddo 30 continue data_points = i-1 if(iq_debug.eq.1) close(19) write(*,*) 'Generated ',kcount,' good points from ',i ! write(*,*) 'Number rejected (not v significant)= ',badpoints write(*,*) 'Total weight = ',weight_total print*,'No. points inside FZ = ',jcount print*,'You have approx. ',float(jcount)/float(fz_count) 1 ,' points per cell' print*,'This may help you decide how to choose ' 1 ,'the smoothing width' print*,'Now to normalize the values ' vavg = 0. do i=1,m r1 = float(i-(m/2))*deltae do j=1,m r2 = float(j-(m/2))*deltae do k=1,m r3 = float(k-(m/2))*deltae if(voxvol(i,j,k).gt.0.) then c array(i,j,k)=100.*array(i,j,k) ! array(i,j,k) = array(i,j,k) ! 1 *volsum/voxvol(i,j,k)/weight_total array(i,j,k) = array(i,j,k) 1 * volsum / voxvol(i,j,k) / weight_total ! to match RFdists_SQU else write(*,*) 'VOXVOL= ',voxvol(i,j,k),' at ',i,j,k endif if(iq_symm.eq.0) then inside = test_fz_cubic(r1,r2,r3) elseif(iq_symm.eq.1) then inside = test_fz_hcp(r1,r2,r3) elseif(iq_symm.eq.2) then inside = test_fz_orth(r1,r2,r3) endif if(inside) then if(array(i,j,k).gt.vmax) vmax=array(i,j,k) if(array(i,j,k).lt.vmin) vmin=array(i,j,k) vavg = vavg + (array(i,j,k)*voxvol(i,j,k)) endif enddo enddo enddo c vavg=vavg/float(fz_count)/volsum vavg = vavg / volsum write(*,*) 'before SMOOTH: min/avg/max ',vmin,vavg,vmax c write(*,*) 'Enter a smoothing width (1.5 suggested)' c read(*,*) wdth c wdth=6. CALL SMARRAY(SMF,WDTH,m) CALL SMoothit3D(SMF,array,rmax,m) write(*,*) 'SMoothit3D: rmax= ',rmax c$$$ print*,'print out a value?! yes=1' c$$$ read*,iq_print c$$$ if(iq_print.eq.1) then c$$$ 1011 print*,'enter i,j,k :' c$$$ read*,i,j,k c$$$ print*,'array val = ',array(i,j,k) c$$$ print*,'another value? yes=1' c$$$ read*,iq_2 c$$$ if(iq_2.eq.1) goto 1011 c$$$ endif vavg = 0. do i=1,m r1 = float(i-(m/2))*deltae do j=1,m r2 = float(j-(m/2))*deltae do k=1,m r3 = float(k-(m/2))*deltae c if(voxvol(i,j,k).gt.0.) then c array(i,j,k) = array(i,j,k) c 1 *volsum/voxvol(i,j,k)/weight_total c else c write(*,*) 'VOXVOL= ',voxvol(i,j,k),' at ',i,j,k c endif if(iq_symm.eq.0) then inside = test_fz_cubic(r1,r2,r3) elseif(iq_symm.eq.1) then inside = test_fz_hcp(r1,r2,r3) elseif(iq_symm.eq.2) then inside = test_fz_orth(r1,r2,r3) endif if(inside) then if(array(i,j,k).gt.vmax) vmax=array(i,j,k) if(array(i,j,k).lt.vmin) vmin=array(i,j,k) vavg=vavg+array(i,j,k) endif enddo enddo enddo c vavg=vavg/float(kcount) vavg=vavg/float(fz_count) write(*,*) 'after SMOOTH: min/avg/max ',vmin,vavg,vmax c write out various files for subsequent use if(add_mdf) then open(unit=17,file = fname2(1:(dotpos-1))//'_mdf.txt' $ ,status='unknown') else open(unit=17,file = fname(1:(dotpos-1))//'_mdf.txt' $ ,status='unknown') endif write(17,*) ' R1 R2 R3 MDF angle ' jcount = 0 do i = 1,m r1 = float(i-(m/2))*deltae do j = 1,m r2 = float(j-(m/2))*deltae do k = 1,m r3 = float(k-(m/2))*deltae rnorm = sqrt(r1**2+r2**2+r3**2) angle = degrad * 2. * atan(rnorm) write(17,"(5(2x,f15.4))") r1,r2,r3,array(i,j,k),angle if(iq_symm.eq.0) then inside = test_fz_cubic(r1,r2,r3) elseif(iq_symm.eq.1) then inside = test_fz_hcp(r1,r2,r3) elseif(iq_symm.eq.2) then inside = test_fz_orth(r1,r2,r3) endif if(inside) then vavg = vavg+array(i,j,k) jcount = jcount + 1 endif enddo enddo enddo close(17) vavg=vavg/float(jcount) print*, 'after WRITE: avg ',vavg open(unit=17,file = fname(1:(dotpos-1))//'_100axis.txt' $ ,status='unknown') write(17,*) ' R1 R2 R3 MDF_100 angle ' j = m/2 r2 = 0. k = m/2 r3 = 0. do i = 1,m r1 = float(i-(m/2))*deltae rnorm = sqrt(r1**2+r2**2+r3**2) angle = degrad * 2. * atan(rnorm) write(17,"(5(2x,f15.4))") r1,r2,r3,array(i,j,k),angle enddo open(unit=17,file = fname(1:(dotpos-1))//'_110axis.txt' $ ,status='unknown') write(17,*) ' R1 R2 R3 MDF_110 angle ' k = m/2 r3 = 0. do i = 1,m j = i r1 = float(i-(m/2))*deltae r2 = float(j-(m/2))*deltae rnorm = sqrt(r1**2+r2**2+r3**2) angle = degrad * 2. * atan(rnorm) write(17,"(5(2x,f15.4))") r1,r2,r3,array(i,j,k),angle enddo open(unit=17,file = fname(1:(dotpos-1))//'_111axis.txt' $ ,status='unknown') write(17,*) ' R1 R2 R3 MDF_111 angle ' do i = 1,m j = i r1 = float(i-(m/2))*deltae r2 = float(j-(m/2))*deltae k = i r3 = float(k-(m/2))*deltae rnorm = sqrt(r1**2+r2**2+r3**2) angle = degrad * 2. * atan(rnorm) write(17,"(5(2x,f15.4))") r1,r2,r3,array(i,j,k),angle enddo c now for contour plotting with psplot.txt call defcol(colrval,cval,greyval) c set the contour levels, colour values write(*,*) 'Contour levels set at:' do i=1,nval write(*,"('Contour no. ',i3,' value = ',f10.3)") i,cval(i) enddo ! CALL GetArg(7,arg_line) ! shifted this call to the beginning of the program to simplify the logic if(arg_line_contour.eq.'') then write(*,*) 'Enter 0 if this is OK: ' print*,' Or any other number to change the values' read(*,*) iq0 else read(arg_line_contour,*) iq0 endif if(iq0.ne.0) then do i = 1,nval write(*,*) 'Enter contour value ',i read(*,*) cval(i) enddo endif txold=0.+(xsize/2.) tyold=0.+(ysize/2.) ! txold = -2. ! tyold = -2. fname2 = fname(1:(dotpos-1))//'_mdf.ps' inquire(file=fname2,exist=repeat) if(repeat) then ! only delete the PS if conversion happened print*,'deleting PS file' call system('rm '//fname2) endif call newdev (fname2,15) call psinit(.true.) c nbox=8 nbox=7 if(iq_symm.eq.0) then ! delta = float(m)/2.345/float(nbox) delta = float(m)/2.59/float(nbox) elseif(iq_symm.eq.1) then delta = (0.28/deltae)/float(nbox) ! limit on R3 is set by tan(15degr) for Hexagonals = 0.268 elseif(iq_symm.eq.2) then delta = (1.08/deltae)/float(nbox) ! this will give us the upper half of the space for now endif if(iq_debug.eq.1) then open(unit=3,file='RFdist_output.txt',status='unknown',recl=200) endif txdiff = tx(1)-txold tydiff = ty(1)-tyold txold = tx(1) tyold = ty(1) c$$$ if(iq_symm.eq.2) y_margin = 3. * y_margin ! kludge for ortho case ! call plot(txdiff,tydiff,-3) if(iq_weight.eq.0 ) then call keksym( 4.0 , ycharpos $ ,.12,'Number Based/ No Weighting ',0.,27,0) endif if(iq_30_rotate.eq.1 ) then call keksym( 5.75 , ycharpos $ ,.12,'Data rotated 30 degrees',0.,23,0) end if call keksym( 0.75 , ycharpos $ , 0.15, 'No. of points = ',0.,16,0) call keknum( 2.85 , ycharpos $ , 0.15, float(data_points),0.,0,0) call keksym( 0.75 , ycharpos - 0.25 , .15 , 1 'Data file: '//fname(1:(dotpos-1))//'.mdf',0.,15+dotpos,0) call keksym( 0.75 , ycharpos - 0.5 $ , 0.15, 'Smoothing = ',0.,12,0) width = wdth call keknum( 2.75, ycharpos - 0.5 $ ,0.15,width,0.,2,0) if(iq_symm.eq.1) then call keksym( 0.75 , ycharpos - 0.7 $ , 0.12, 'Caution: R1 is Cartesian X,' $ //' = 2110 (the a-axis) for the X-convention' $ ,0.,69,0) call keksym( 0.85 , ycharpos - 0.9 $ , 0.12,' or 1010 for the Y-convention' ,0.,30,0) endif call plot(txdiff,tydiff,-3) do i = 1,nbox ! this will be the loop over the sections ! do i = 1,1 ! this will be the loop over the sections c index1=idelta*(i-1)+3 index1 = int(float(m/2) + (delta * float(i-1))) r3 = float(index1-(m/2))*deltae ! NOTE - this is exactly the same calculation of R3 as above for binning write(*,*) 'Section at index=, R3 = ',index1,r3 txdiff=tx(i)-txold tydiff=ty(i)-tyold txold=tx(i) tyold=ty(i) call plot(txdiff,tydiff,-3) write(*,*) write(*,*) 'Section no. ',i call keksym(-.2+(xsize/2.),.15+(ysize/2.) & ,.15,trichar(i),0.,1,0) r3_real = r3 if(iq_symm.eq.0) then call keksym(0.7+(xsize/2.),1.8+(ysize/2.) $ ,.15,'R3 =',0.,4,0) call keknum(0.7+(xsize/2.),1.55+(ysize/2.) $ ,0.15,r3_real,0.,3,0) elseif(iq_symm.eq.1) then call keksym(0.3+(xsize/2.),1.8+(ysize/2.7) $ ,.15,'R3 =',0.,4,0) call keknum(0.3+(xsize/2.),1.55+(ysize/2.7) $ ,0.15,r3_real,0.,3,0) elseif(iq_symm.eq.2) then call keksym(-0.9+(xsize/2.), 1.8+(ysize/2.7) $ ,.15,'R3 =',0.,4,0) call keknum(-0.9+(xsize/2.), 1.55+(ysize/2.7) $ ,0.15,r3_real,0.,3,0) endif do j=1,m do k=1,m garray(j,k)=array(j,k,index1) ! copy section if(garray(j,k).le.0.) garray(j,k) = 1.e-2 enddo enddo if(iq_debug.eq.1) then write(3,*) write(3,*) 'Array for section ',i $ ,',index= ',index1,' R3= ',r3 do k=1,m write(3,"(50(x,g8.2))") (garray(j,k), j=1,m,5) enddo endif call border(xsize,ysize,0000,0000,4,2,4,2) c no ticks, no borders call GSAV if(iq_symm.eq.0) then call defclipRFtri(xclip,yclip,nclip,xsize,r3) elseif(iq_symm.eq.1) then rlimit = 2. / deltae / float(m) rl2 = ((float((m/2) - 1 ) + 0.5) / float(m)) * xsize ! rl2 = (float((m/2) - 1 ) / float(m)) * xsize ! position of the point corresponding to x=0 print*,'scale fraction for Hex box = ',rlimit call defclipRFhcp(xclip,yclip,nclip,xsize,rlimit,rl2) elseif(iq_symm.eq.2) then ! call defclipsquare(xclip,yclip,nclip,xsize) call defclip_orth(xclip,yclip,nclip,xsize) endif call clipbox(xclip,yclip,nclip) call confill(garray,m,m,m,xsize,ysize,cval(1),greyval, & nval,0,0.) call conrec(garray,m,m,m,xsize,ysize,cval(1),nval) call clip ! clips off the FZ triangle call GREST if(iq_symm.eq.0) then call defclipbox(xclip,yclip,nclip,xsize) c$$$ elseif(iq_symm.eq.1) then c$$$ call defclipRFhcp(xclip,yclip,nclip,xsize,rlimit,r1) c$$$ elseif(iq_symm.eq.2) then c$$$ call defclipsquare(xclip,yclip,nclip,xsize) endif call DRWCRV (xclip, yclip, nclip, .01, .true.) c draws the full triangle around the section enddo ! loop over the boxes ! This section is special to HEX plots if(iq_symm.eq.1) then do vertical_section = 1 , 2 ! R1, R2 vertical sections print* if( vertical_section.eq.1 ) print*,'100-001 section' if( vertical_section.eq.2 ) print*,'010-001 section' do j=1,m do k=1,m if( vertical_section.eq.1 ) then garray(j,k) = array(j,(m/2),k) ! copy section else garray(j,k) = array((m/2),j,k) ! copy section end if if(garray(j,k).le.0.) garray(j,k) = 1.e-2 enddo enddo if(iq_debug.eq.1) then write(3,*) if( vertical_section.eq.1 ) then write(3,*) 'Array for 100-001 section ' else write(3,*) 'Array for 010-001 section ' end if do k=1,m write(3,"(50(x,g8.2))") (garray(j,k), j=1,m,5) enddo endif txdiff=tx(7+vertical_section)-txold tydiff=ty(7+vertical_section)-tyold txold=tx(7+vertical_section) tyold=ty(7+vertical_section) call plot(txdiff,tydiff,-3) call border(xsize,ysize,0000,0000,4,2,4,2) c no ticks, no borders if( vertical_section.eq.1 ) then call keksym(0.2+(xsize/2.), 2.0+(ysize/2.7) $ , .15, 'R2 = 0', 0., 6, 0) else call keksym(0.2+(xsize/2.), 2.0+(ysize/2.7) $ , .15, 'R1 = 0', 0., 6, 0) end if call keksym(-0.3+(xsize/2.), 0.15+(ysize/2.) & , 0.15, trichar(i), 0., 1, 0) call GSAV ! rlimit = 2. / deltae / float(m) ! print*,'scale fraction for Hex box = ',rlimit call defcliprect(xclip,yclip,nclip,xsize,rlimit,rl2) call clipbox(xclip,yclip,nclip) call confill(garray,m,m,m,xsize,ysize,cval(1),greyval, & nval,0,0.) call conrec(garray,m,m,m,xsize,ysize,cval(1),nval) call clip ! clips call GREST c call defcliprect(xclip,yclip,nclip,xsize,rlimit,rl2) call DRWCRV (xclip, yclip, nclip, .01, .true.) c draws the box around the section end do endif ! only for HEX plots c now to make a legend (still B&W) if(iq_symm.eq.1) then txdiff = tx(9)-txold + 0.99 tydiff = ty(9)-tyold + 0.9 ! make room for an additional plot for HEX else txdiff = tx(8)-txold tydiff = ty(8)-tyold endif txold = tx(8) tyold = ty(8) call plot(txdiff,tydiff,-3) call setlw(0.025) rtmp9 = 0.25 call keksym(-0.1+(xsize/2.),1.0+(ysize/2.),.15,'MRD',0.,3,0) do i=0,nval xx1 = 1.55 + (xsize/2.) xx2 = xx1 + 2.*rtmp9 yy1 = 0.2 + (ysize/2.) + float(i)*rtmp9 yy2 = yy1 height = rtmp9 c call keknum(-0.8+xx1, yy1-(height/3.), 0.15, cval(i)/100., 0., 2, 0) call keknum(-0.8+xx1, yy1-(height/3.), 0.15, cval(i), 0., 2, 0) CALL RECTFILg(XX1,YY1,XX2,YY2,HEIGHT,greyval(i)) CALL SLDLIN (xX1, yY1, xx1-0.10, yY1, 0.015) enddo ! *********** draw axes txdiff = tx(11)-txold tydiff = ty(11)-tyold txold = tx(11) tyold = ty(11) call plot(txdiff,tydiff,-3) call setlw(0.04) index1 = int( float(m/2) ) r3_real = float( index1 - (m/2) ) * deltae ! note that r3 is the value of the 3rd component, ! which for cubics controls the shape of the box if(iq_symm.eq.0) then ! call defclipRFtri(xclip,yclip,nclip,xsize,r3) xclip(1) = ( xsize / 2. ) + ( r3_real / tend * xsize / 2. ) yclip(1) = ( xsize / 2. ) + ( r3_real / tend * xsize / 2. ) xclip(2) = xsize * 0.8 yclip(2) = yclip(1) xclip(3) = xclip(1) yclip(3) = xsize * 0.8 elseif(iq_symm.eq.1) then rlimit = 2. / deltae / float(m) rl2 = ((float((m/2) - 1 ) + 0.5) / float(m)) * xsize ! call defclipRFhcp(xclip,yclip,nclip,xsize,rlimit,rl2) xclip(1) = rl2 yclip(1) = rl2 xclip(2) = xclip(1) + (( xsize / 2.) * rlimit * 0.8 ) yclip(2) = yclip(1) xclip(3) = xclip(1) yclip(3) = yclip(1) + (( xsize / 2.) * rlimit * 0.8 ) elseif(iq_symm.eq.2) then ! call defclip_orth(xclip,yclip,nclip,xsize) xclip(1) = xsize * 0.5 yclip(1) = xsize * 0.5 xclip(2) = xclip(1) + ( xsize * 0.3 ) yclip(2) = yclip(1) xclip(3) = xclip(1) yclip(3) = yclip(1) + ( xsize * 0.3 ) end if CALL ARROW ( xclip(1) , yclip(1) , xclip(2) , yclip(2) & , 0.07 , 15. , 0 ) CALL ARROW ( xclip(1) , yclip(1) , xclip(3) , yclip(3) & , 0.07 , 15. , 0 ) call keksym( xclip(1)+0.1 , yclip(3)-0.3 $ , 0.2 , 'R2' , 0. , 2 , 0 ) call keksym( xclip(2)-0.4 , yclip(1)+0.15 $ , 0.2 , 'R1' , 0. , 2 , 0 ) ! ****** end of axes call plotnd if(iq_debug.eq.1) close(3) fname3 = fname(1:(dotpos-1))//'_mdf_BW.gif' print*,'converting to gif' call system('convert '//fname2//' '//fname3) inquire(file=fname3,exist=repeat) c$$$ if(repeat) then ! only delete the PS if conversion happened c$$$ print*,'deleting PS file' c$$$ call system('rm '//fname2) c$$$ endif ! *********!!!!!!!!! end of greyscale ! repeat to get colour plot txold=0.+(xsize/2.) tyold=0.+(ysize/2.) c fname='RFdist.colour'//'.ps' fname2=fname(1:(dotpos-1))//'_mdf_clr.ps' inquire(file=fname2,exist=repeat) if(repeat) then ! only delete the PS if conversion happened print*,'deleting PS file' call system('rm '//fname2) endif call newdev (fname2,15) call psinit(.true.) txdiff = tx(1)-txold tydiff = ty(1)-tyold txold = tx(1) tyold = ty(1) ! call plot(txdiff,tydiff,-3) if(iq_weight.eq.0 ) then call keksym( 4.0 , ycharpos $ ,.12,'Number Based/ No Weighting ',0.,27,0) endif if(iq_30_rotate.eq.1 ) then call keksym( 5.75 , ycharpos $ ,.12,'Data rotated 30 degrees',0.,23,0) end if call keksym( 0.75 , ycharpos $ , 0.15, 'No. of points = ',0.,16,0) call keknum( 2.85 , ycharpos $ , 0.15, float(data_points),0.,0,0) call keksym( 0.75 , ycharpos - 0.25 , .15 , 1 'Data file: '//fname(1:(dotpos-1))//'.mdf',0.,15+dotpos,0) call keksym( 0.75 , ycharpos - 0.5 $ , 0.15, 'Smoothing = ',0.,12,0) width = wdth call keknum( 2.75, ycharpos - 0.5 $ ,0.15,width,0.,2,0) if(iq_symm.eq.1) then call keksym( 0.75 , ycharpos - 0.7 $ , 0.12, 'Caution: R1 is Cartesian X,' $ //' = 2110 (the a-axis) for the X-convention' $ ,0.,69,0) call keksym( 0.85 , ycharpos - 0.9 $ , 0.12,' or 1010 for the Y-convention' ,0.,30,0) endif call plot(txdiff,tydiff,-3) do i = 1,nbox ! this will be the loop over the sections c index1=idelta*(i-1)+3 index1 = int(float(m/2) + (delta * float(i-1))) r3 = float(index1-(m/2))*deltae write(*,*) 'Section (colour) at index=, R3 = ',index1,r3 txdiff=tx(i)-txold tydiff=ty(i)-tyold txold=tx(i) tyold=ty(i) call plot(txdiff,tydiff,-3) write(*,*) write(*,*) 'Section no. ',i call keksym(-.2+(xsize/2.),.15+(ysize/2.) & ,.15,trichar(i),0.,1,0) r3_real = r3 if(iq_symm.eq.0) then call keksym(0.7+(xsize/2.),1.8+(ysize/2.) $ ,.15,'R3 =',0.,4,0) call keknum(0.7+(xsize/2.),1.55+(ysize/2.) $ ,0.15,r3_real,0.,3,0) elseif(iq_symm.eq.1) then call keksym(0.3+(xsize/2.),1.8+(ysize/2.7) $ ,.15,'R3 =',0.,4,0) call keknum(0.3+(xsize/2.),1.55+(ysize/2.7) $ ,0.15,r3_real,0.,3,0) elseif(iq_symm.eq.2) then call keksym(-0.9+(xsize/2.), 1.8+(ysize/2.7) $ ,.15,'R3 =',0.,4,0) call keknum(-0.9+(xsize/2.), 1.55+(ysize/2.7) $ ,0.15,r3_real,0.,3,0) endif do j=1,m do k=1,m garray(j,k) = array(j,k,index1) ! copy section if(garray(j,k).le.0.) garray(j,k) = 1.e-2 enddo enddo call border(xsize,ysize,0000,0000,4,2,4,2) c no ticks, no borders call GSAV if(iq_symm.eq.0) then call defclipRFtri(xclip,yclip,nclip,xsize,r3) elseif(iq_symm.eq.1) then rlimit = 2. / deltae / float(m) rl2 = ((float((m/2) - 1 ) + 0.5) / float(m)) * xsize ! rl2 = (float((m/2) - 1 ) / float(m)) * xsize ! position of the point corresponding to x=0 call defclipRFhcp(xclip,yclip,nclip,xsize,rlimit,rl2) elseif(iq_symm.eq.2) then ! call defclipsquare(xclip,yclip,nclip,xsize) call defclip_orth(xclip,yclip,nclip,xsize) endif call clipbox(xclip,yclip,nclip) call concolr(garray,m,m,m,xsize,ysize,cval(1),colrval, & nval,0,0.) call conrec(garray,m,m,m,xsize,ysize,cval(1),nval) call clip ! clips off the FZ triangle call GREST if(iq_symm.eq.0) then call defclipbox(xclip,yclip,nclip,xsize) c$$$ elseif(iq_symm.eq.1) then c$$$ call defclipRFhcp(xclip,yclip,nclip,xsize,rlimit,rl2) c$$$ elseif(iq_symm.eq.2) then c$$$ call defclipsquare(xclip,yclip,nclip,xsize) endif call DRWCRV (xclip, yclip, nclip, .01, .true.) c draws the full triangle around the section enddo ! This section is special to HEX plots if(iq_symm.eq.1) then do vertical_section = 1 , 2 ! R1, R2 vertical sections print* if( vertical_section.eq.1 ) print*,'100-001 section' if( vertical_section.eq.2 ) print*,'010-001 section' do j=1,m do k=1,m if( vertical_section.eq.1 ) then garray(j,k) = array(j,(m/2),k) ! copy section else garray(j,k) = array((m/2),j,k) ! copy section end if if(garray(j,k).le.0.) garray(j,k) = 1.e-2 enddo enddo if(iq_debug.eq.1) then write(3,*) if( vertical_section.eq.1 ) then write(3,*) 'Array for 100-001 section ' else write(3,*) 'Array for 010-001 section ' end if do k=1,m write(3,"(50(x,g8.2))") (garray(j,k), j=1,m,5) enddo endif txdiff=tx(7+vertical_section)-txold tydiff=ty(7+vertical_section)-tyold txold=tx(7+vertical_section) tyold=ty(7+vertical_section) call plot(txdiff,tydiff,-3) call border(xsize,ysize,0000,0000,4,2,4,2) c no ticks, no borders if( vertical_section.eq.1 ) then call keksym(0.2+(xsize/2.), 2.0+(ysize/2.7) $ , .15, 'R2 = 0', 0., 6, 0) else call keksym(0.2+(xsize/2.), 2.0+(ysize/2.7) $ , .15, 'R1 = 0', 0., 6, 0) end if call keksym(-0.3+(xsize/2.), 0.15+(ysize/2.) & , 0.15, trichar(i), 0., 1, 0) call GSAV ! rlimit = 2. / deltae / float(m) ! print*,'scale fraction for Hex box = ',rlimit call defcliprect(xclip,yclip,nclip,xsize,rlimit,rl2) call clipbox(xclip,yclip,nclip) call concolr(garray,m,m,m,xsize,ysize,cval(1),colrval, & nval,0,0.) call conrec(garray,m,m,m,xsize,ysize,cval(1),nval) call clip ! clips call GREST c call defcliprect(xclip,yclip,nclip,xsize,rlimit,rl2) call DRWCRV (xclip, yclip, nclip, .01, .true.) c draws the box around the section end do endif ! only for HEX plots ! ***************** now to make a legend if(iq_symm.eq.1) then txdiff = tx(9)-txold + 0.99 tydiff = ty(9)-tyold + 0.9 ! make room for an additional plot for HEX txold = tx(9) + 0.99 tyold = ty(9) + 0.9 else txdiff = tx(8)-txold tydiff = ty(8)-tyold txold = tx(8) tyold = ty(8) endif call plot(txdiff,tydiff,-3) call setlw(0.025) rtmp9 = 0.25 call keksym(-0.1+(xsize/2.),1.0+(ysize/2.),.15,'MRD',0.,3,0) do i=0,nval xx1 = 1.55 + (xsize/2.) xx2 = xx1 + 2.*rtmp9 yy1 = 0.2 + (ysize/2.) + float(i)*rtmp9 yy2 = yy1 height = rtmp9 c call keknum(-0.8+xx1, yy1-(height/3.), 0.15, cval(i)/100., 0., 2, 0) call keknum(-0.8+xx1, yy1-(height/3.), $ 0.15, cval(i), 0., 2, 0) CALL RECTFILC(XX1,YY1,XX2,YY2,HEIGHT,colrval(1,i) $ ,colrval(2,i),colrval(3,i)) CALL SLDLIN (xX1, yY1, xx1-0.10, yY1, 0.015) enddo ! *********** draw axes txdiff = tx(11)-txold tydiff = ty(11)-tyold txold = tx(11) tyold = ty(11) call plot(txdiff,tydiff,-3) call setlw(0.04) index1 = int( float(m/2) ) r3_real = float( index1 - (m/2) ) * deltae ! note that r3 is the value of the 3rd component, ! which for cubics controls the shape of the box if(iq_symm.eq.0) then ! call defclipRFtri(xclip,yclip,nclip,xsize,r3) xclip(1) = ( xsize / 2. ) + ( r3_real / tend * xsize / 2. ) yclip(1) = ( xsize / 2. ) + ( r3_real / tend * xsize / 2. ) xclip(2) = xsize * 0.8 yclip(2) = yclip(1) xclip(3) = xclip(1) yclip(3) = xsize * 0.8 elseif(iq_symm.eq.1) then rlimit = 2. / deltae / float(m) rl2 = ((float((m/2) - 1 ) + 0.5) / float(m)) * xsize ! call defclipRFhcp(xclip,yclip,nclip,xsize,rlimit,rl2) xclip(1) = rl2 yclip(1) = rl2 xclip(2) = xclip(1) + (( xsize / 2.) * rlimit * 0.8 ) yclip(2) = yclip(1) xclip(3) = xclip(1) yclip(3) = yclip(1) + (( xsize / 2.) * rlimit * 0.8 ) elseif(iq_symm.eq.2) then xclip(1) = xsize * 0.5 yclip(1) = xsize * 0.5 xclip(2) = xclip(1) + ( xsize * 0.3 ) yclip(2) = yclip(1) xclip(3) = xclip(1) yclip(3) = yclip(1) + ( xsize * 0.3 ) end if CALL ARROW ( xclip(1) , yclip(1) , xclip(2) , yclip(2) & , 0.07 , 15. , 0 ) CALL ARROW ( xclip(1) , yclip(1) , xclip(3) , yclip(3) & , 0.07 , 15. , 0 ) call keksym( xclip(1)+0.1 , yclip(3)-0.3 $ , 0.2 , 'R2' , 0. , 2 , 0 ) call keksym( xclip(2)-0.4 , yclip(1)+0.15 $ , 0.2 , 'R1' , 0. , 2 , 0 ) ! ****** end of axes call plotnd fname3 = fname(1:(dotpos-1))//'_mdf_clr.gif' print*,'converting to gif (colour)' call system('convert '//fname2//' '//fname3) inquire(file=fname3,exist=repeat) c$$$ if(repeat) then ! only delete the PS if conversion happened c$$$ print*,'deleting PS file' c$$$ call system('rm '//fname2) c$$$ endif call exit(0) c ********** stop here for now! end c c ++++++++++++++++++++++++++++++ c function test_fz_cubic(r1,r2,r3) implicit none logical test_fz_cubic double precision tend,r1,r2,r3 c CODE:: tend = dsqrt(2.d0) - 1.d0 c defines edge of the R-F space test_fz_cubic = .false. if(r1.ge.0.d0 1 .and.r1.le.tend 2 .and.r2.ge.r3 3 .and.r2.le.r1 4 .and.r3.le.r1 5 .and.r3.ge.0.d0 6 .and.(r1+r2+r3).le.1.d0) then test_fz_cubic = .true. endif return end c c ++++++++++++++++++++++++++++++ c function test_fz_hcp(r1,r2,r3) implicit none logical test_fz_hcp double precision r1,r2,r3 double precision rtmp1,rtmp2,rtmp3 double precision pi,rtmp,eps double precision tan60,sin60,cos60,tan30,sin30,cos30,tan15 c CODE:: test_fz_hcp = .false. eps = 1.0d-2 PI = 3.14159265d0 ! print*,'testing for inside HCP FZ' ! print*,'input values: ', r1,r2,r3 tan60 = dtan (pi*60./180.) sin60 = dsin (pi*60./180.) cos60 = dcos (pi*60./180.) tan30 = dtan (pi*30./180.) sin30 = dsin (pi*30./180.) cos30 = dcos (pi*30./180.) tan15 = dtan (pi*15./180.) rtmp1 = r1*tan30 - r2 rtmp2 = r1*cos30 + r2*sin30 rtmp3 = r1*cos60 + r2*sin60 if(r2.ge.0. & .and.r1.le.1.d0 & .and.r3.ge.0.d0 & .and.r3.le.tan15 c these two IFs put the point inside the FZ for R3 & .and.rtmp1.ge.(0.d0 - eps) c line from origin at 30 degr & .and.rtmp2.le.1.d0 c 1st sector & .and.rtmp3.le.1.d0) then test_fz_hcp = .true. endif ! print*,'result = ',test_fz_hcp return end c c ++++++++++++++++++++++++++++++ c function test_fz_orth(r1,r2,r3) implicit none logical test_fz_orth double precision r1,r2,r3 ! FZ is a cube size 1, 0 to 1 all 3 directions c CODE:: test_fz_orth = .false. ! if(r1.ge.-1.d0 if(r1.ge.0.d0 1 .and.r1.le.1.d0 ! 2 .and.r2.ge.-1.d0 2 .and.r2.ge.0.d0 3 .and.r2.le.1.d0 4 .and.r3.le.1.d0 ! 5 .and.r3.ge.-1.d0) then 5 .and.r3.ge.0.d0) then test_fz_orth = .true. endif return end c c ++++++++++++++++++++++++++++++ c subroutine defcol(colrval,cval,greyval) c making scale bar real colrval(3,0:10),cval(0:10) real greyval(0:10) colrval(1,0)=0.0 colrval(2,0)=0.0 colrval(3,0)=0.0 ! black colrval(1,1)=0.0 colrval(2,1)=0.0 colrval(3,1)=1.0 ! blue colrval(1,2)=0.0 colrval(2,2)=1.0 colrval(3,2)=1.0 ! cyan colrval(1,3)=0.0 colrval(2,3)=1.0 colrval(3,3)=0.0 ! green colrval(1,4)=1.0 colrval(2,4)=1.0 colrval(3,4)=0.0 ! yellow colrval(1,5)=1.0 colrval(2,5)=0.0 colrval(3,5)=0.0 ! red colrval(1,6)=1.0 colrval(2,6)=0.0 colrval(3,6)=1.0 ! magenta colrval(1,7)=0.5 colrval(2,7)=0.5 colrval(3,7)=0.5 ! ? c add more values if you use more than 6 contour values! cval(0) = 0. cval(1) = 0.25 cval(2) = 0.50 cval(3) = 1. cval(4) = 1.50 cval(5) = 2.00 cval(6) = 4.00 cval(7) = 8.00 greyval(10)=1. greyval(9)=1. greyval(8)=1. greyval(7)=1. greyval(6)=0.95 greyval(5)=0.9 greyval(4)=0.8 greyval(3)=0.7 greyval(2)=0.6 greyval(1)=0.4 greyval(0)=0.1 return end c c ________________________________________ c subroutine defclipbox(xclip,yclip,nclip,side) c defines the clipping region c in this case a TRIANGLE c implicit none real side real xclip(100),yclip(100) integer nclip c nclip=4 c xclip(1)=0. c yclip(1)=0. c xclip(2)=side c yclip(2)=0. c xclip(3)=side c yclip(3)=side c xclip(4)=0. c yclip(4)=side xclip(1) = side/2. yclip(1) = side/2. xclip(2) = side yclip(2) = side/2. xclip(3) = side yclip(3) = side c xclip(4) =side/2. c yclip(4) =side xclip(4) = xclip(1) yclip(4) = yclip(1) c c do 6010, i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c 6010 continue c return end c c ________________________________________ c subroutine defclipRFhcp(xclip,yclip,nclip,side,rlimit,rl2) c defines the clipping region c in this case a polygon for HCP c corners of the FZ, in degrees c 0,0 c tan(45),0 c tan(45),tan(45)*tan(15) c tan(45)*cos(30),tan(45)*sin(30) implicit none double precision pi real xclip(100),yclip(100),side,rlimit integer nclip real tan45,sin45,cos45,tan30 $ ,sin30,cos30,tan15,sin15,cos15 real rside,r2,rl2 c CODE:: PI = 3.14159265 rside = (side / 2.) * rlimit r2 = side / 2. tan15 = tan (pi*15./180.) sin15 = sin (pi*15./180.) cos15 = cos (pi*15./180.) tan30 = tan (pi*30./180.) sin30 = sin (pi*30./180.) cos30 = cos (pi*30./180.) tan45 = tan (pi*45./180.) sin45 = sin (pi*45./180.) cos45 = cos (pi*45./180.) nclip = 5 xclip(1) = rl2 yclip(1) = rl2 xclip(2) = xclip(1) + (tan45 * rside) yclip(2) = yclip(1) xclip(3) = xclip(2) yclip(3) = rl2 + (tan45 * tan15 * rside) xclip(4) = rl2 + (tan45 * cos30 * rside) yclip(4) = rl2 + (tan45 * sin30 * rside) xclip(5) = xclip(1) yclip(5) = yclip(1) c do 6010, i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c 6010 continue return end c c ________________________________________ c subroutine defcliprect(xclip,yclip,nclip,side,rlimit,rl2) c defines the clipping region c in this case a RECTANGLE for a 100-001 plot in HEX implicit none real side,rlimit,rside,r2,tan15,tan45,pi,rl2 real xclip(100),yclip(100) integer nclip c CODE:: PI = 3.14159265 tan45 = tan (pi*45./180.) tan15 = tan (pi*15./180.) rside = (side / 2.) * rlimit r2 = side / 2. nclip=5 xclip(1) = rl2 yclip(1) = rl2 xclip(2) = rl2 + (tan45 * rside) yclip(2) = rl2 xclip(3) = xclip(2) yclip(3) = rl2 + (tan15 * rside) xclip(4) = xclip(1) yclip(4) = yclip(3) xclip(5) = xclip(1) yclip(5) = yclip(1) c do 6010, i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c 6010 continue return end c c ________________________________________ c subroutine defclipRFtri(xclip,yclip,nclip,side,r3) c defines the clipping region c in this case a TRIANGLE c implicit none double precision tend,r3,ax,ay,bx,by real xclip(100),yclip(100),side integer nclip logical corner,corner2 tend=sqrt(2.)-1. ax=(1.-r3)/2. ay=(1.-r3)/2. by=(1.-r3-tend) bx=1.-(2.*r3) corner=(ax.lt.tend) corner2=(bx.lt.tend) c special version for RF triangles c SIDE/2 assumed to be equivalent to sqrt(2)-1 xclip(1)=side/2.+(r3/tend*side/2.) yclip(1)=side/2.+(r3/tend*side/2.) xclip(2)=side yclip(2)=side/2.+(r3/tend*side/2.) if(corner) then if(corner2) then xclip(2)=side/2.+(bx/tend*side/2.) yclip(2)=side/2.+(r3/tend*side/2.) xclip(3)=side/2.+(ax/tend*side/2.) yclip(3)=side/2.+(ay/tend*side/2.) xclip(4)=xclip(1) yclip(4)=yclip(1) nclip=4 else xclip(3)=side yclip(3)=side/2.+(by/tend*side/2.) xclip(4)=side/2.+(ax/tend*side/2.) yclip(4)=side/2.+(ay/tend*side/2.) xclip(5)=xclip(1) yclip(5)=yclip(1) nclip=5 endif else xclip(3)=side yclip(3)=side xclip(4)=xclip(1) yclip(4)=yclip(1) nclip=4 endif c do 6010, i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c 6010 continue return end c ________________________________________ c subroutine defclip(xclip,yclip,nclip,radius) c defines the clipping region c in this case a circle c double precision angle,pi real xclip(100),yclip(100),radius integer nclip double precision tri_size,uparam,xu,yu,zu data pi/3.141592654/ c nclip=100 c THIS SECTION DEFINES A CIRCLE do 100, i=1,99 angle=float(i-1)/99.*2.*pi xclip(i)=radius+radius*cos(angle) yclip(i)=radius+radius*sin(angle) 100 continue xclip(100)=xclip(1) yclip(100)=yclip(1) c c do 6010, i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c 6010 continue c return end c c ________________________________________ c subroutine defclipsquare(xclip,yclip,nclip,side) c defines the clipping region c in this case a SQUARE of edge length SIDE implicit none real side real xclip(100),yclip(100) integer nclip c CODE:: nclip = 5 xclip(1) = 0. yclip(1) = 0. xclip(2) = side yclip(2) = 0. xclip(3) = side yclip(3) = side xclip(4) = 0. yclip(4) = side xclip(5) = xclip(1) yclip(5) = yclip(1) c do i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c enddo return end c c ________________________________________ c subroutine defclip_orth(xclip,yclip,nclip,side) c defines the clipping region c in this case a SQUARE of edge length SIDE/2 ! for orthorhombic misorientations implicit none real side real xclip(100),yclip(100) integer nclip c CODE:: nclip = 5 xclip(1) = side*.5 yclip(1) = side*.5 xclip(2) = side yclip(2) = side*.5 xclip(3) = side yclip(3) = side xclip(4) = side*.5 yclip(4) = side xclip(5) = xclip(1) yclip(5) = yclip(1) c do i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c enddo 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 iq_symm integer numsymm,numsamp double precision quatsymm,quatsamp common/a1symm/ quatsymm(4,48),quatsamp(4,4),numsamp,numsymm c common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp c CODE:: rad=57.29578 c write(*,*) 'reading symmetry operators from "quat.symm.cubic"' if(iq_symm.eq.0) then open(unit=3,file = 'quat.symm.cubic',status='old') elseif(iq_symm.eq.1) then open(unit=3,file = 'quat.symm.hex',status='old') elseif(iq_symm.eq.2) then open(unit=3,file = 'quat.symm.orth',status='old') endif 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 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 subroutine rodcom(r1,r2,rout) c c combines two rotations as Rodrigues vectors c where R1 is the first rotation, and R2 is the second rotation c ROUT= R1+R2- R1 cross R2 / 1 - (R1dotR2) double precision r1(3),r2(3),rout(3) double precision dot,cross(3) dot=1.-scalar(r1,r2) call vecpro2(r1,r2,cross) if(abs(dot).gt.1e-5) then rout(1)=(r1(1)+r2(1)-cross(1))/dot rout(2)=(r1(2)+r2(2)-cross(2))/dot rout(3)=(r1(3)+r2(3)-cross(3))/dot else rout(1)=1e38 rout(2)=1e38 rout(3)=1e38 c write(*,*) 'warning: RODCOM ->infinite result' endif return end c c ____________________________ c c subroutine rod2q(d,quat) double precision 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 ____________________________ c c subroutine q2rod(d,quat) implicit none double precision d(3),quat(4),rtmp,rnorm,rnorm2 integer i double precision stheta,ttheta c converts a quaternion to a Rodrigues vector 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-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! 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 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) real p1,p,p2 double precision q(4) double precision c1,c,c2,s1,s,s2,d(3,3),tmp(4),sin2,cos2,rtmp,rnorm c CODE:: 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 comquat2(qq1,qq2,qresult) double precision qq1(4),qq2(4),qresult(4) 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 ! CODE:: 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 write(*,*) 'qresult ',qresult return end c c _____________________________ c c subroutine presamp(qq,lindex,qresult) double precision qq(4),qresult(4) integer lindex integer numsymm,numsamp double precision quatsymm,quatsamp common/a1symm/ quatsymm(4,48),quatsamp(4,4),numsamp,numsymm c common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp 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 € 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 c subroutine presymm(qq,lindex,qresult) double precision qq(4),qresult(4) integer lindex c include 'common.f' integer numsymm,numsamp double precision quatsymm,quatsamp common/a1symm/ quatsymm(4,48),quatsamp(4,4),numsamp,numsymm c 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) double precision qq(4),qresult(4) integer lindex c include 'common.f' integer numsymm,numsamp double precision quatsymm,quatsamp common/a1symm/ quatsymm(4,48),quatsamp(4,4),numsamp,numsymm c 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 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 "dot" QQ c c Thus, for active rotations (standard definition of orientation in c mechanics) this is suitable for applying CRYSTAL symmetry c if(lindex.gt.numsymm) stop 'error in postsymm, lindex>numsymm' 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) ! very odd: this code below was the same as in presymm! ! picked up version from 3ddg ... 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 function scalar(a,b) c return SCALAR PRODUCT of A and B double precision 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 double precision 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 SMARRAY(SMF,WDTH,nnsize) c WDTH is the width of the smoothing filter c SMF is the array that contains the filter values c from JS Kallend's smoothing program for pole figures! c c IMPLICIT INTEGER*2 (I-N) implicit none integer i,k,nnsize c parameter (nnsize=100) DOUBLE PRECISION SMF((-1*nnsize):nnsize),dif,tot,wdth,rwdth,deg,sdeg c DO 100 I=2,19 c do 100 i=1,19 c DEG=(I-1)*0.087266 c SDEG=SIN(DEG) C 5*PI/180 c RWDTH=WDTH/SDEG c special version: adr: no adjustment of width for declination angle c dec 01 c c write(*,*) 'Enter smoothing width, e.g. 5' c read(*,*) wdth RWDTH=WDTH c c PRINT *,I,RWDTH,SDEG DO 200 K=0,nnsize c DIF=K*5./RWDTH DIF=float(K)/RWDTH ! take out 5 degree grid assumption IF(DIF.GT.8.)THEN SMF(K)=0. ELSE SMF(K)=EXP(-DIF*DIF) SMF(-K)=SMF(K) ENDIF c print *,'interm. result= ',k,smf(k) 200 continue c TOT=0. DO 210 K=1-nnsize,nnsize TOT=TOT+SMF(K) 210 CONTINUE c DO 220 K=-1*nnsize,nnsize SMF(K)=SMF(K)/TOT 220 CONTINUE 100 CONTINUE c now for special treatment of the first row c DO 230 K=-1*nnsize,nnsize c 230 SMF(1,K)=1./72 c c DO 300 I=1,nnsize SMF(-1*nnsize)=SMF(-1*nnsize)/2. SMF(nnsize)=SMF(-1*nnsize) c PRINT 1000,(SMF(K),K=(-1*nnsize),nnsize) 1000 FORMAT('smf=',1X, 10F8.3) c 300 CONTINUE c RETURN END c c ________________________________________ c C SUBROUTINE SMoothit(SMF,array,rmax,nnsize) c parameter (nnsize=100) c IMPLICIT INTEGER*2 (I-N) c INTEGER*2 ICOUNT(19,72) double precision array(nnsize,nnsize) DOUBLE PRECISION SMF((-1*nnsize):nnsize),smct(nnsize) c DOUBLE PRECISION SMF(19,-36:36),SMCT(72) double precision rmax,ccnt rmax=0. c write(*,*) 'ICOUNT: ' c write(*,*) 'SMoothit debug: 1st column, BEFORE 2nd index' c PRINT 1000,(array(I,1),i=1,nnsize) c DO 100 I=1,nnsize DO 150 K=1,nnsize CCNT=0. DO 160 IX=(-1*nnsize),nnsize if(smf(ix).gt.0.005)then INDX=K+IX IF(INDX.LT.1)INDX=INDX+nnsize IF(INDX.GT.nnsize)INDX=INDX-nnsize CCNT=CCNT+array(I,INDX)*SMF(IX) endif 160 CONTINUE SMCT(K)=CCNT 150 CONTINUE DO 110 K=1,nnsize c array(I,K)=SMCT(K)+.5 array(I,K)=SMCT(K) if(array(i,k).gt.rmax) rmax=array(i,k) 110 CONTINUE c 100 CONTINUE c c write(*,*) 'SMoothit debug: 1st column, 2nd index' c PRINT 1000,(array(I,1),i=1,nnsize) 1000 FORMAT(10(1x,f9.3)) c DO I=1,nnsize DO K=1,nnsize CCNT=0. DO IX=(-1*nnsize),nnsize if(smf(ix).gt.0.005)then INDX=K+IX IF(INDX.LT.1)INDX=INDX+nnsize IF(INDX.GT.nnsize)INDX=INDX-nnsize CCNT=CCNT+array(indx,I)*SMF(IX) endif enddo SMCT(K)=CCNT enddo DO K=1,nnsize c array(k,i)=SMCT(K)+.5 array(k,i)=SMCT(K) if(array(k,i).gt.rmax) rmax=array(k,i) enddo c PRINT 1000,(array(I,K),K=1,nnsize) c 1000 FORMAT(18I4) enddo c c write(*,*) 'SMoothit debug: 1st column, 1st index' c PRINT 1000,(array(I,1),i=1,nnsize) c RETURN END c c c ________________________________________ c C SUBROUTINE SMoothit3D(SMF,array,rmax,nnsize) c c adaptation of Kallend's smoothing program to 3D arrays c adr mar 03 c c parameter (nnsize=100) c IMPLICIT INTEGER*2 (I-N) implicit none c INTEGER*2 ICOUNT(19,72) integer i,j,k,mm,ix,indx,nnsize double precision array(nnsize,nnsize,nnsize) DOUBLE PRECISION SMF((-1*nnsize):nnsize),smct(nnsize) c DOUBLE PRECISION SMF(19,-36:36),SMCT(72) double precision rmax,ccnt rmax=0. c write(*,*) 'ICOUNT: ' c write(*,*) 'SMoothit debug: 1st column, BEFORE 2nd index' c PRINT 1000,(array(I,1),i=1,nnsize) c c smooth over the 1st index print*,'x-dimension smooth' DO mm=1,nnsize ! 110 c if(mod(mm,10).eq.0) write(*,*) '1st index; MM= ',mm DO I=1,nnsize ! 150 do K=1,nnsize ! 155 CCNT=0. DO IX=(-1*nnsize),nnsize ! 160 if(smf(ix).gt.0.005)then INDX=K+IX IF(INDX.LT.1)INDX=INDX+nnsize IF(INDX.GT.nnsize)INDX=INDX-nnsize CCNT=CCNT+array(mm,I,INDX)*SMF(IX) endif 160 enddo SMCT(K)=CCNT 155 enddo DO K=1,nnsize array(mm,I,K)=SMCT(K) if(array(mm,i,k).gt.rmax) rmax=array(mm,i,k) enddo 150 enddo 110 enddo 100 CONTINUE c c write(*,*) 'SMoothit debug: 1st column, 2nd index' c PRINT 1000,(array(I,1),i=1,nnsize) 1000 FORMAT(10(1x,f9.3)) c smooth over the 2nd index print*,'y-dimension smooth' DO mm=1,nnsize ! 110 c if(mod(mm,10).eq.0) write(*,*) '2nd index; MM= ',mm DO I=1,nnsize ! 150 do K=1,nnsize ! 155 CCNT=0. DO IX=(-1*nnsize),nnsize ! 160 if(smf(ix).gt.0.005)then INDX=K+IX IF(INDX.LT.1)INDX=INDX+nnsize IF(INDX.GT.nnsize)INDX=INDX-nnsize CCNT=CCNT+array(mm,INDX,I)*SMF(IX) endif enddo SMCT(K)=CCNT enddo DO K=1,nnsize array(mm,K,I)=SMCT(K) if(array(mm,k,i).gt.rmax) rmax=array(mm,k,i) enddo enddo enddo c smooth over the 3rd index print*,'z-dimension smooth' DO mm=1,nnsize ! 110 c if(mod(mm,10).eq.0) write(*,*) '3rd index; MM= ',mm DO I=1,nnsize ! 150 do K=1,nnsize ! 155 CCNT=0. DO IX=(-1*nnsize),nnsize ! 160 if(smf(ix).gt.0.005)then INDX=K+IX IF(INDX.LT.1)INDX=INDX+nnsize IF(INDX.GT.nnsize)INDX=INDX-nnsize CCNT=CCNT+array(INDX,mm,I)*SMF(IX) endif enddo SMCT(K)=CCNT enddo DO K=1,nnsize array(K,mm,I)=SMCT(K) if(array(k,mm,i).gt.rmax) rmax=array(k,mm,i) enddo enddo enddo c PRINT 1000,(array(I,K),K=1,nnsize) c 1000 FORMAT(18I4) c write(*,*) 'SMoothit debug: 1st column, 1st index' c PRINT 1000,(array(I,1),i=1,nnsize) RETURN END c ________________________________________ ! include '/code/paul.lee.codes/psplot.txt' ! include '/Volumes/7DEC05/code/paul.lee.codes/psplot.txt' include '/Volumes/7DEC05/code/paul.lee.codes/psplot-19Sep09.txt'