program RF_misor c /usr/local/bin/g77 -ffixed-line-length-none -g -o RF_misor RF_misor_13Oct06.f C program for reading in an energy distribution c and plotting it out in triangular sections in RF space c adr Feb 03; last edit apr 03; for use with *.mdf output from REX3D 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 c parameter (m=100) parameter (m = 150) 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*40,fname2*40,fname3*40 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(10),ty(10),txdiff,tydiff,txold,tyold real txc(10),tyc(10) 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 double precision tan60,sin60,cos60,tan30,sin30,cos30 logical test_fz_hcp,test_fz_cubic,inside double precision x_margin,RFe,y_margin integer dotpos,lastchar logical add_mdf integer data_points ! how many lines read in double precision angle ! for outputting MDF values integer jcount c following variables control plotting options c CODE:: PI=3.14159265 print*,'Welcome to MISOR. Usage = ', 1 './RFmisor input_data smooth_width[1.5] quats[0]/RFvectors[1]'// 2 ' cubic[0]/hcp[1] max_pts[99999999] accept_contours[0]' print* 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 c write(*,*) 'Enter the filename for the input data, e.g. injy01.13:' c write(*,*) 'An extension of .mdf is assumed' c read(*,'(a)') fname2 c$$$ index=40 c$$$ do 5, i=40,1,-1 c$$$c write(*,*) 'i,fname2(i:i) ',i,fname2(i:i) c$$$ if(fname2(i:i).eq.' ') index=i c$$$5 continue c$$$ index=index-1 c$$$ fname=fname2(1:index) c$$$ fname2((index+1):(index+1))=char(0) c$$$ write(*,*) 'fname2,index,fname' c$$$ write(*,*) fname2,', ',index,', ',fname c$$$ c$$$ fname=fname2(1:index) 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. inquire(file=fname(1:(dotpos-1))//'.mdf',exist=repeat) if(.not.repeat) then write(*,*) 'Did NOT find ',fname(1:(dotpos-1))//'.mdf' write(*,*) 'Will look for series of .mdf files' single=.false. else write(*,*) 'Found ',fname(1:(dotpos-1))//'.mdf' endif open(unit=2,name=fname(1:(dotpos-1))//'.mdf',status='old') read(2,*) c title line should have: c the crystal symmetry as a letter in the 1st column c so C for cubic, H for hexagonal, etc. c the sample symmetry as a letter in the 2nd column c so O for orthorhombic, X for monoclinic on 100, T for triclinic, etc. 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]?' read(*,*) iq_read_type else read(arg_line,*) iq_read_type endif if(iq_read_type.lt.0.or.iq_read_type.gt.1) then stop 'cannot read that type' endif CALL GetArg(4,arg_line) if(arg_line.eq.'') then print*,'cubic [= 0] or hexagonal [= 1]?' read(*,*) iq_symm else read(arg_line,*) iq_symm endif if(iq_symm.lt.0.or.iq_symm.gt.1) then stop 'unknown symmetry type' endif CALL GetArg(5,arg_line) if(arg_line.eq.'') then print*,' WC on input file > number of lines as 1st result' call system('wc '//fname(1:(dotpos-1))//'.mdf') print*,'Maximum points ?' read(*,*) iqnum else read(arg_line,*) iqnum 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. if(iq_symm.eq.0) then c xsize = 4. ! change this value if you plot a different area in each section xsize = 2. ! un-clipped version else xsize = 5. endif ysize = xsize txold=0. tyold=0. c tend=sqrt(2.)-1. c defines edge of the R-F space c 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' c x_margin = 0.5 x_margin = .25 + xsize y_margin = 1. tx(1) = x_margin ty(1) = y_margin + 1. tx(2) = x_margin ty(2) = y_margin + 3.1 tx(3) = x_margin ty(3) = y_margin + 5.2 tx(4) = x_margin ty(4) = y_margin + 7.3 c tx(5) = x_margin + (0.5 * xsize) + x_margin tx(5) = x_margin + (1.1 * xsize) + x_margin ty(5) = y_margin + 1. c tx(6) = x_margin + (0.5 * xsize) + x_margin tx(6) = x_margin + (1.1 * xsize) + x_margin ty(6) = y_margin + 3.1 c tx(7) = x_margin + (0.5 * xsize) + x_margin tx(7) = x_margin + (1.1 * xsize) + x_margin ty(7) = y_margin + 5.2 c tx(8) = x_margin + (0.5 * xsize) + x_margin tx(8) = x_margin + (1.1 * xsize) + x_margin ty(8) = y_margin + 7.3 c 7 triangles, ll corners 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=999. if(iq_symm.eq.0) then c deltae = 2.0*tend/float(m) deltae = 2.0 * tend/ 100. RFe = tend else c deltae = 2.0/float(m) deltae = 2.0 / 100. RFe = 1. c need more space for HCP endif print*,'Initializing array' 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=float(i-(m/2))*deltae c reverse of this is c i = int(r1/deltae) + m/2 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) voxvol(i,j,k) = 1. / (1.+rnorm**2)**2 c approx. volume of cell if(iq_symm.eq.0) then inside = test_fz_cubic(r1,r2,r3) else inside = test_fz_hcp(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 c icount=0 badpoints=0 weight_total=0.0 jcount = 0 do i = 1,iqnum ! label 580 if(mod(i,1000).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 if(iq_read_type.eq.0) then read(2,*,err=580,end=30) idummy,jdummy,(qtmp(ii),ii=1,4),weight else read(2,*,err=580,end=30) (rodr(ii),ii=1,3),weight call rod2q(rodr,qtmp) endif 575 format(4f12.4) c print*,i,'Q ',qtmp,' weight ',weight if(weight.le.0.0) weight=0. 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 c write(*,*) c write(*,*) 'input line number ',i c c now apply symmetry elements c note that the RF vector is an active rotation from ref to xtal frame c do j=1,numsymm c do k=1,4 c qqsym(k)=quatsymm(k,j) c enddo c call comquat2(qtmp,qqsym,qresult) call presymm(qtmp,j,qint) do k=1,numsymm call postsymm(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) c write(*,"('rresult = ',3(1x,f8.4))") rresult 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 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 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 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 c write(*,*) 'indices into ARRAY: ', ir1val,ir2val,ir3val if(iq_symm.eq.0) then inside = test_fz_cubic(rodr(1),rodr(2),rodr(3)) else inside = test_fz_hcp(rodr(1),rodr(2),rodr(3)) endif if(inside) then weight_total = weight_total+weight jcount = jcount + 1 endif 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 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*,'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 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) else inside = test_fz_hcp(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 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) else inside = test_fz_hcp(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 open(unit=17,name=fname(1:(dotpos-1))//'_mdf.txt',status='unknown') 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) else inside = test_fz_hcp(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,name=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,name=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,name=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(6,arg_line) if(arg_line.eq.'') then write(*,*) 'Enter 0 if this is OK: ' print*,' Or any other number to change the values' read(*,*) iq0 else read(arg_line,*) 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.) fname2 = fname(1:(dotpos-1))//'.mdf.ps' 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) else delta = (0.577/deltae)/float(nbox) endif open(unit=3,file='RFdist.output.txt',status='unknown') txdiff = tx(1)-txold tydiff = ty(1)-tyold txold = tx(1) tyold = ty(1) call plot(txdiff,tydiff,-3) c call keksym(0.2+(xsize/2.),-.3+(ysize/2.),.15,'No. of points = ',0.,16,0) c call keknum(2.2+(xsize/2.),-.3+(ysize/2.),0.15,float(data_points),0.,0,0) c call keksym(0.2+(xsize/2.),-.55+(ysize/2.),.15, call keksym(0.2-(xsize/1.3),-.3-(ysize/6.),.15 1 ,'No. of points = ',0.,16,0) call keknum(2.2-(xsize/1.3),-.3-(ysize/6.) 1 ,0.15,float(data_points),0.,0,0) call keksym(0.2-(xsize/1.3),-.55-(ysize/6.),.15, 1 'Data file: '//fname(1:(dotpos-1))//'.mdf',0.,15+dotpos,0) c call keksym(0.2+(xsize/2.),-.8+(ysize/2.),.15,'Smoothing = ',0.,12,0) c width = wdth c call keknum(2.2+(xsize/2.),-.8+(ysize/2.),0.15,width,0.,2,0) call keksym(0.2-(xsize/1.3),-.8-(ysize/6.),.15,'Smoothing = ',0.,12,0) width = wdth call keknum(2.2-(xsize/1.3),-.8-(ysize/6.),0.15,width,0.,2,0) 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 at R3 = ',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 c call keksym(-.2+(xsize/2.),.15+(ysize/2.) call keksym(-.2-(xsize/1.),.15+(ysize/2.) & ,.15,trichar(i),0.,1,0) r3_real = r3 if(iq_symm.eq.0) then c call keksym(0.7+(xsize/2.),1.8+(ysize/2.),.15,'R3 =',0.,4,0) c call keknum(0.7+(xsize/2.),1.55+(ysize/2.),0.15,r3_real,0.,3,0) call keksym(0.7-(xsize/1.),1.8+(ysize/4.),.15,'R3 =',0.,4,0) call keknum(0.7-(xsize/1.),1.55+(ysize/4.),0.15,r3_real,0.,3,0) else c call keksym(0.7+(xsize/2.),1.8+(ysize/2.7),.15,'R3 =',0.,4,0) c call keknum(0.7+(xsize/2.),1.55+(ysize/2.7),0.15,r3_real,0.,3,0) call keksym(0.7-(xsize/1.),1.8+(ysize/2.7),.15,'R3 =',0.,4,0) call keknum(0.7-(xsize/1.),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(index1,j,k) ! copy section enddo enddo write(3,*) write(3,*) 'Array for section ',i,' R3= ',r3 do k=1,m write(3,"(10(x,g8.2))") (garray(j,k), j=1,m) 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) else call defclipRFhcp(xclip,yclip,nclip,xsize) endif c 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) else call defclipRFhcp(xclip,yclip,nclip,xsize) endif call DRWCRV (xclip, yclip, nclip, .03, .true.) c draws the full triangle around the section enddo c now to make a legend txdiff = tx(8)-txold tydiff = ty(8)-tyold txold = tx(8) tyold = ty(8) call plot(txdiff,tydiff,-3) call setlw(0.025) rtmp9 = 0.25 c call keksym(-0.1+(xsize/2.),1.0+(ysize/2.),.15,'MRD',0.,3,0) call keksym(-0.1-(xsize/1.),1.0+(ysize/2.),.15,'MRD',0.,3,0) do i=0,nval xx1 = 1.55 + (xsize/1.) 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 call plotnd close(3) fname3 = fname(1:(dotpos-1))//'_mdf_noClip_BW.jpg' print*,'converting to jpg' call system('convert '//fname2//' '//fname3) print*,'deleting PS file' call system('rm '//fname2) c 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' 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) c call keksym(0.2+(xsize/2.),-.3+(ysize/2.),.15,'No. of points = ',0.,16,0) c call keknum(2.3+(xsize/2.),-.3+(ysize/2.),0.15,float(data_points),0.,0,0) c call keksym(0.2+(xsize/2.),-.55+(ysize/2.),.15, call keksym(0.2-(xsize/1.3),-.3-(ysize/6.),.15 1 ,'No. of points = ',0.,16,0) call keknum(2.3-(xsize/1.3) 1 ,-.3-(ysize/6.),0.15,float(data_points),0.,0,0) call keksym(0.2-(xsize/1.3),-.55-(ysize/6.),.15, 1 'Data file: '//fname(1:(dotpos-1))//'.mdf',0.,15+dotpos,0) c call keksym(0.2+(xsize/2.),-.8+(ysize/2.),.15,'Smoothing = ',0.,12,0) c width = wdth c call keknum(2.2+(xsize/2.),-.8+(ysize/2.),0.15,width,0.,2,0) call keksym(0.2-(xsize/1.3),-.8-(ysize/6.),.15,'Smoothing = ',0.,12,0) width = wdth call keknum(2.2-(xsize/1.3),-.8-(ysize/6.),0.15,width,0.,2,0) do i=1,nbox ! this will be the loop over the sections c index1=idelta*(i-1)+3 index1=(float(m)/2.)+(delta*float(i-1)) r3=float(index1-(m/2))*deltae write(*,*) 'Section (colour) at R3 = ',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 c call keksym(-.2+(xsize/2.),.15+(ysize/2.) call keksym(-.2-(xsize/1.),.15+(ysize/2.) & ,.15,trichar(i),0.,1,0) r3_real = r3 if(iq_symm.eq.0) then c call keksym(0.7+(xsize/2.),1.8+(ysize/2.),.15,'R3 =',0.,4,0) c call keknum(0.7+(xsize/2.),1.55+(ysize/2.),0.15,r3_real,0.,3,0) call keksym(0.7-(xsize/1.),1.8+(ysize/4.),.15,'R3 =',0.,4,0) call keknum(0.7-(xsize/1.),1.55+(ysize/4.),0.15,r3_real,0.,3,0) else c call keksym(0.7+(xsize/2.),1.8+(ysize/2.7),.15,'R3 =',0.,4,0) c call keknum(0.7+(xsize/2.),1.55+(ysize/2.7),0.15,r3_real,0.,3,0) call keksym(0.7-(xsize/1.),1.8+(ysize/2.7),.15,'R3 =',0.,4,0) call keknum(0.7-(xsize/1.),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(index1,j,k) ! copy section 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) else call defclipRFhcp(xclip,yclip,nclip,xsize) endif c 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) else call defclipRFhcp(xclip,yclip,nclip,xsize) endif call DRWCRV (xclip, yclip, nclip, .03, .true.) c draws the full triangle around the section enddo c now to make a legend txdiff = tx(8)-txold tydiff = ty(8)-tyold txold = tx(8) tyold = ty(8) call plot(txdiff,tydiff,-3) call setlw(0.025) rtmp9 = 0.25 c call keksym(-0.1+(xsize/2.),1.0+(ysize/2.),.15,'MRD',0.,3,0) call keksym(-0.1-(xsize/1.),1.0+(ysize/2.),.15,'MRD',0.,3,0) do i=0,nval c xx1 = 1.55 + (xsize/2.) xx1 = 1.55 - (xsize/1.) 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 call plotnd fname3 = fname(1:(dotpos-1))//'_mdf_noClip_clr.jpg' print*,'converting to jpg (colour)' call system('convert '//fname2//' '//fname3) print*,'deleting PS file' call system('rm '//fname2) 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=sqrt(2.)-1. c defines edge of the R-F space test_fz_cubic = .false. if(r1.ge.0. 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. 6 .and.(r1+r2+r3).le.1.) 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.0e-5 PI=3.14159265 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.) tan15 = tan (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.0 & .and.r3.ge.0. & .and.r3.le.tan15 c these two IFs put the point inside the FZ for R3 & .and.rtmp1.ge.0. c line from origin at 30 degr & .and.rtmp2.le.1. c 1st sector & .and.rtmp3.le.1.) then test_fz_hcp = .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) c6010 continue c return end c c ________________________________________ c subroutine defclipRFhcp(xclip,yclip,nclip,side) 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 integer nclip double precision tan45,sin45,cos45,tan30,sin30,cos30,tan15,sin15,cos15 c CODE:: PI = 3.14159265 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) = side/2. yclip(1) = side/2. xclip(2) = side/2. + (tan45 * side/2.) yclip(2) = side/2. xclip(3) = side/2. + (tan45 * side/2.) yclip(3) = side/2. + (tan45 * tan15 * side/2.) xclip(4) = side/2. + (tan45 * cos30 * side/2.) yclip(4) = side/2. + (tan45 * sin30 * side/2.) xclip(5) = xclip(1) yclip(5) = yclip(1) c do 6010, i=1,nclip c write(*,*) 'defclip: ',i,xclip(i),yclip(i) c6010 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 c special version for RF triangles c SIDE/2 assumed to be equivalent to Ã2-1 c 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) c6010 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) c6010 continue c 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,name='quat.symm.cubic',status='old') else open(unit=3,name='quat.symm.hex',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,name='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 c subroutine comquat2(qq1,qq2,qresult) double precision qq1(4),qq2(4),qresult(4) c c algorithm for forming resultant quaternion as Q1 x Q2^-1 c where Q2 follows Q1, i.e. is the second rotation; c note change of signs to get inverse of second orientation c qresult(1)=qq1(1)*qq2(4)-qq1(4)*qq2(1) & +qq1(2)*qq2(3)-qq1(3)*qq2(2) qresult(2)=qq1(2)*qq2(4)-qq1(4)*qq2(2) & +qq1(3)*qq2(1)-qq1(1)*qq2(3) qresult(3)=qq1(3)*qq2(4)-qq1(4)*qq2(3) & +qq1(1)*qq2(2)-qq1(2)*qq2(1) qresult(4)=qq1(4)*qq2(4)+qq1(1)*qq2(1) & +qq1(2)*qq2(2)+qq1(3)*qq2(3) c c write(*,*) 'qresult ',qresult c return end c c _____________________________ c c 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 € 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' 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 c230 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) c300 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)) 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 write(*,*) 'SMoothit debug: 1st column, 1st index' c PRINT 1000,(array(I,1),i=1,nnsize) 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 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) c1000 FORMAT(18I4) c c write(*,*) 'SMoothit debug: 1st column, 1st index' c PRINT 1000,(array(I,1),i=1,nnsize) RETURN END c c ________________________________________ c include '/code/paul.lee.codes/psplot.txt' c