program RFfuncplot
c
c  program for reading in an energy distribution
c and plotting it out in triangular sections in RF space
c
c  adr  Feb 03; last edit mar 03; jul 06
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
c  based on MAPE.F
c
	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
c
	integer m,iqcolor,kcount,iq0,nval
	parameter (m=100)
	real array(m,m,m)
c  ARRAY contains the data to be plotted
c
	integer mresolve,ndata,jndex,index,i,j,k,ii
	real tend,pi,pi2,degrad,deltae
	integer ndiscrete,itmp(3),icount
	real vmax,vmin,vavg
	character fname*20
	real rtmp(3),rodr(3),rfsymm(3),rresult(3),weight,scale
	real quat(4),qresult(4)
	integer ir3val,ir2val,ir1val,type
	logical goodpoint
	character  crys_sym_type,sam_sym_type
	integer numsymm,numsamp
	real quatsymm,quatsamp
c
	common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c
	integer index1,nbox,idelta,isymm
	real garray(m,m),r1,r2,r3,xsize,ysize,rnorm,rnormmin,peak
        real exclude,cutoff

	real wdth,rmax,SMF(m,(-1*m):m)
	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
c
c  following variables control plotting options
c
	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.
	ysize=2.
	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
        print*,'RFfunc  adapted to show excluded zone'
        print*
        print*,'Enter size of zone in degrees:'
        read(*,*) exclude
        cutoff = tan(exclude/degrad/2.)

	call textur4              !reads in symmetry operators
	write(*,*) 'Found ',numsymm,' xtal symms',' and '
	write(*,*) numsamp,' sample symms'
c
	tx(1)=0.5
	ty(1)=1.
	tx(2)=0.5
	ty(2)=3.1
	tx(3)=0.5
	ty(3)=5.2
	tx(4)=0.5
	ty(4)=7.3
	tx(5)=4.0
	ty(5)=1.
	tx(6)=4.0
	ty(6)=3.1
	tx(7)=4.0
	ty(7)=5.2
	tx(8)=4.0
	ty(8)=7.3
c  7 triangles, ll corners
c
	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
c  10 triangles
c
	mresolve=10  !  one page for now!
	kcount=0
	vmax=0.
	vavg=0.
	vmin=999.
	deltae=2.*tend/float(m)
	write(*,*) 'Initializing array(s)'
	do 10, i=1,m
		if(mod(i,25).eq.0) write(*,*) 'Row number ',i
		do 10, j=1,m
			do 10, k=1,m
				array(i,j,k)=0.0
10	continue
c  this is inefficient for cubic systems, but does allow for less symmetric crystal symmetries
c
	crys_sym_type='C'
	sam_sym_type='O'
c
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.
c
	peak=8.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)
c                 print*,'r1/r2/r3/rnorm: ',r1,'/',r2,'/',r3,'/',rnorm
                 if(rnorm.lt.cutoff) then
                    array(i,j,k) = .3
c                    print*,'inside: cutoff: ',cutoff
c                    print*,'r1/r2/r3/rnorm: ',r1,'/',r2,'/',r3,'/',rnorm
                 else
                    array(i,j,k) = 6.
c                    array(i,j,k)=peak*(1./(1.+rnorm**2))**2 !  fixed 16Mar03
                 endif
c  for plotting an excluded zone in the center of the space

		 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)
	      enddo
	   enddo
	enddo

	vavg=vavg/float(kcount)
	write(*,*) 'min/avg/max ',vmin,vavg,vmax

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
	write(*,*) 'Enter 0 if this is OK, any other number to change'
	read(*,*) iq0
	if(iq0.ne.0) then
           do i=1,nval
              write(*,*) 'Enter contour value ',i
              read(*,*) cval(i)
           enddo
	endif

        txold=0.
        tyold=0.

	fname='RFfunc_GREY'//'.ps'
	call newdev (fname,15)
	call psinit(.true.)
	nbox=8
	idelta=m/nbox
	open(unit=3,file='RFfunc.output.txt',status='unknown')
	do i=1,nbox   ! this will be the loop over the sections
           index1 = idelta*(i-1) + 2
           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
           call keksym(-.2,.15,.15,trichar(i),0.,1,0)
           call keksym(2.1,1.8,.15,'R3 =',0.,4,0)
           call keknum(2.1,1.55,0.15,r3,0.,3,0)
           do j=1,m
              do k=1,m
                 garray(j,k) = array(j,k,index1) !  copy section 
              enddo
           enddo
c		wdth=12.
c		CALL SMARRAY(SMF,WDTH,m)
c		CALL SMoothit(SMF,garray,rmax,m)
c  these were calls to 2D smoothing in each section plane
           write(3,*)
           write(3,*) 'Array for section ',i,'  R3= ',r3
           do k=1,m
              write(3,"(10(x,f8.4))") (garray(j,k), j=1,m)
           enddo

           write(*,*)
           write(*,*) 'Array for section ',i,'  R3= ',r3
           do k=1,m,10
              write(*,"(10(x,f8.4))") (garray(j,k), j=1,m,10)
           enddo

           call border(xsize,ysize,0000,0000,4,2,4,2)
c  no ticks, no borders

           call defclipbox(xclip,yclip,nclip,xsize)
           call GSAV
           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
           call GREST
           call DRWCRV (xclip, yclip, nclip, .03, .true.)

	enddo
	call plotnd
        call system('convert RFfunc_GREY.ps RFfunc_GREY.jpg')
        call system('rm RFfunc_GREY.ps ')

        txold=0.
        tyold=0.
	fname='RFfunc_CLR'//'.ps'
	call newdev (fname,15)
	call psinit(.true.)

	do i=1,nbox   ! this will be the loop over the sections
           index1 = idelta*(i-1) + 2
           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
           call keksym(-.2,.15,.15,trichar(i),0.,1,0)
           call keksym(2.1,1.8,.15,'R3 =',0.,4,0)
           call keknum(2.1,1.55,0.15,r3,0.,3,0)
           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 defclipbox(xclip,yclip,nclip,xsize)
           call GSAV
           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
           call GREST
           call DRWCRV (xclip, yclip, nclip, .03, .true.)

	enddo
	call plotnd
        call system('convert RFfunc_CLR.ps RFfunc_CLR.jpg')
        call system('rm RFfunc_CLR.ps ')
	close(3)

	call exit
c  ********** stop here for now!
c
c	vold=vavg/float(icount)
	write(*,*) 'Max., Min., Total = ',vmax,vmin,vavg
	if(vmin.le.0.0) vmin=1e-6
c  fix zero point for log scale!
	write(*,*) 'log10(Max.),log10( Min.) = ',alog10(vmax),alog10(vmin)
	write(*,*) 'Enter value for Average (change to rescale plot)'
	read(*,*) vavg
	vmax=vmax/vavg
	vmin=vmin/vavg
	write(*,*) 'What scale do you want to use (enter Max, e.g.)? '
	read(*,*) scale
	write(*,*) 'Color (=0) or black & white (=1)?'
	read(*,*) iqcolor
c
	type=3
	call rodplot(array,vmax,vmin,vavg,deltae,
     &  mresolve,type,scale,iqcolor,m)
c  plot the energy distribution
c  as a postscript file for viewing
1000	continue
c
	stop 'done!'
	end
c
c  ++++++++++++++++++++++++++++++
c
	subroutine defcol(colrval,cval,greyval)
c    making scale bar
c
	real colrval(3,0:10),cval(0:10)
	real greyval(0:10)
c
	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.1
	cval(2)=0.2
	cval(3)=0.5
	cval(4)=1.
	cval(5)=2.
	cval(6)=4.
	cval(7)=8.
c
	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
c
	return
	end
c
c  __________________________________
c
	subroutine rodplot(energy,vmax,vmin,vavg,deltae,
     &  mresolve,type,scale,iqcolor,m)
c
c  plots sets of Rodrigues triples in up to 7 triangles
c  output= postscript file
c
	integer type,ntriangs,iqcolor
	integer mmm
	parameter (mmm=400)
	real energy(mmm,mmm,mmm)
	real vmax,vmin,vavg,deltae
c
	integer mresolve
c  MRESOLVE sets the resolution
c
	character*1 ten,one,ctype
c
	real pi,pi2,degrad,tri_size,rfincr
	real   tx(10),ty(10)
c
	tx(1)=1
	ty(1)=1
	tx(2)=1
	ty(2)=3.1
	tx(3)=1
	ty(3)=5.2
	tx(4)=1
	ty(4)=7.3
	tx(5)=3.5
	ty(5)=3.1
	tx(6)=3.5
	ty(6)=5.2
	tx(7)=3.5
	ty(7)=7.3
c  7 triangles, ll corners
c
	pi=4.*atan(1.)
	pi2=2.*atan(1.)
	tend=tan(pi/8.)
	degrad=180./pi
	
c
	ten=char(48+(mresolve/10))
	one=char(48+mod(mresolve,10))
	ctype=char(48+type)
	write(*,*) 'writing out rfdist'//'.'//ten//one//'.ps'
	open(unit=9,file='rfdist'//'.'//ten//one//'.ps',
     &	        status='unknown')
c
	write(9,4010) 
4010	format('%!PS-Adobe-3.0')
	write(9,4020) 
4020	format('%%BoundingBox: 36 36 576 792')
c  box for .5" margins on a 8.5x11 page
	write(9,4030) 
4030	format('%%Creator:RFdist.plot.f')
	write(9,4040) 
4040	format('%%EndComments')
c
	write(9,5010) 
5010	format('/rod2pix {347.646753 mul} def')
c  transforms RF units (0.41 = 45 degr about 100) to 2" or 144 pixels in PS space
	write(9,5020) 
5020	format('/inch {72 mul} def')
c  transforms inches to 72 pixels in PS space
	tri_size=2.0
c  sets size of triangle
	write(9,5030) tri_size,tri_size,-1*tri_size,-1*tri_size
5030	format('/triangle    %stack: x y == ll corner'/
     &	'{ newpath moveto'/
     &	' ',f8.3,' inch 0 rlineto'/
     &	' 0 ',f8.3,' inch rlineto'/
     &	' ',f8.3,' inch ',f8.3,' inch rlineto'/
     &	' closepath } def')
c  defines a triangle of size 2 inches, for plotting within
	box_size=0.2
	write(9,5040) box_size
5040	format('/bsz ',f8.3,' def')
	write(9,5050) 
5050	format('/box    %stack: x y == corner'/
     &	'{ newpath moveto'/
     &	' bsz 2 div inch 0 rmoveto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' bsz -1 mul inch 0 rlineto'/
     &	' 0 bsz -1 mul inch rlineto'/
     &	' bsz inch 0 rlineto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' closepath } def')
c  defines a box of size 0.2 inches, for use as a symbol
	write(9,5060) 
5060	format('/rbox    %draws a box wherever we are!'/
     &	'{ '/
     &	' bsz 2 div inch 0 rmoveto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' bsz -1 mul inch 0 rlineto'/
     &	' 0 bsz -1 mul inch rlineto'/
     &	' bsz inch 0 rlineto'/
     &	' 0 bsz 2 div inch rlineto'/
     &	' closepath } def')
c  defines a box of size 0.2 inches, for use as a symbol
	WRITE(9,4050) 
4050	format('/Times-Roman findfont 18 scalefont setfont')
	write(9,4090) 
4090	format('%%EndProlog')
	ntriangs=7
	    rfincr=1./3./7.
	if(mresolve.le.6) then
	    ntriangs=mresolve
    	    rfincr=tend/mresolve
    	  endif
c  number of sections to plot
c  increment of R3 between sections; MIN allows for fewer increments than seven
c
c	if(type.eq.0) call pltpts(disorient,ratio,ndata,rfincr,tx,ty,
c     &	   tri_size,ntriangs,scale)
c	if(type.eq.1) call pltpts2()
c  not yet implemented
	if(type.eq.2) call pltdist(energy,mmm,deltae,tx,ty,tri_size,
     &	   mresolve,ntriangs,rfincr,scale,iqcolor)
	if(type.eq.3) call pltdist(energy,mmm,deltae,tx,ty,tri_size,
     &	    mresolve,ntriangs,rfincr,scale,iqcolor)
	do 100, i=1,ntriangs
	rflabel=float(i)*rfincr-(rfincr/2.)
c  current R3 value; cells starting at zero (as opposed to centered on zero)
	rtmp=tri_size*rfincr
	write(9,7010) tx(i),ty(i)
7010	format(f8.3,' inch ',f8.3,' inch triangle '/
     &	' [] 0 setdash 2 setlinewidth stroke ')
	write(9,7020) tx(i)+0.2,ty(i)+1.6,rflabel
7020	format(f8.3,' inch ',f8.3,' inch moveto (R3 =',f6.3,') show')
c  should put the label above and to the right of the ll corner
c
	write(9,7030) tx(i),ty(i)
7030	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	write(9,7035) rflabel-(rfincr/2.),rflabel-(rfincr/2.)
7035	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point on the diagonal
	write(9,7040) tend-(rflabel-(rfincr/2.))
7040	format(f8.3,' rod2pix 0 rlineto',
     &	     ' [] 0 setdash 1 setlinewidth stroke')
c  draw solid line to the edge
	write(9,7050) tx(i),ty(i)
7050	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	write(9,7060) rflabel+(rfincr/2.),rflabel+(rfincr/2.)
7060	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point on the diagonal
	write(9,7070) tend-rflabel-(rfincr/2.)
7070	format(f8.3,' rod2pix 0 rlineto',
     &	     ' [2] 1 setdash 1 setlinewidth stroke [] 0 setdash')
c  draw dashed line to the edge
c
	if(rflabel.gt.0.1715) then
c  we need the cut at the top right corner
	write(9,7110) tx(i),ty(i)
7110	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	rod2=1.-(rflabel-(rfincr/2.))-tend
c   based on the cutoff plane being defined by R1+R2+R3=1
	write(9,7115) tend,rod2
7115	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point at edge of triangle
	rod2p=(1.-(rflabel-(rfincr/2.)))/2.
	write(9,7120) (rod2p-tend),(rod2p-rod2)
7120	format(f8.3,' rod2pix ',f8.3,' rod2pix rlineto',
     &	     ' [] 0 setdash 1 setlinewidth stroke')
c  draw solid line to the diagonal
	write(9,7130) tx(i),ty(i)
7130	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	rod2=1.-(rflabel+(rfincr/2.))-tend
	write(9,7135) tend,rod2
7135	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point at edge of triangle
	rod2p=(1.-(rflabel+(rfincr/2.)))/2.
	write(9,7140) rod2p-tend,rod2p-rod2
7140	format(f8.3,' rod2pix ',f8.3,' rod2pix rlineto',
     &	    ' [2] 1 setdash 1 setlinewidth stroke')
c  draw dashed line to the diagonal
	endif
c    draws lines in upper right corner of triangle
c	  
100	continue
c
	write(9,8010) tx(1)+2.7,ty(1)+.5,vmax
8010	format(f8.3,' inch ',f8.3,' inch moveto (max. E = ',f8.3,') show')
	write(9,8020) tx(1)+2.7,ty(1)+.9,vmin
8020	format(f8.3,' inch ',f8.3,' inch moveto (min. E = ',f8.3,') show')
	write(9,8030) tx(1)+2.7,ty(1)+1.3,vavg
8030	format(f8.3,' inch ',f8.3,' inch moveto (avg.=1 = ',f10.6,') show')
	write(9,9010) 
9010	format('showpage  '/'%%PageEnd'/'%%EOF')
	close(9)
	return
	end
c
c  _____________________________
c
	subroutine pltdist(array,m,delta,tx,ty,tri_size,mres
     &	    ,ntriangs,rfincr,scale,iqcolor)
	real array(m,m,m),delta,rod1,rod2,thue,tbright
	real tx(10),ty(10)
	integer m
c
	tend=sqrt(2.)-1.
	cell=tend/float(mres)
	deltax=tri_size/float(mres)
	do 150,  k=1,ntriangs
c  k indexes R3 by taking a section through RF space
	kindex=k
c  take value of k if mresolve is LE. 9.
	if(mresolve.gt.9) kindex=99
c  99 signals that we need to use interpolation to get the energy value
c
	rflabel=float(k)*rfincr-(rfincr/2.)
c  current R3 value; cells starting at zero (as opposed to centered on zero)
	write(9,*) ' gsave 0.5 setlinewidth '
c	write(9,5010) tx(k),ty(k)
c5010	format(f8.3,' inch ',f8.3,' inch triangle clip')
cc  clip the triangle to keep the plot clean since we are using filled boxes
	write(9,7030) tx(k),ty(k)
c	write(*,7030) tx(k),ty(k)
7030	format(' newpath ',f8.3,' inch ',f8.3,' inch moveto')
c  back to the corner of the triangle
	write(9,7035) rflabel-(rfincr/2.),rflabel-(rfincr/2.)
c	write(*,7035) rflabel-(rfincr/2.),rflabel-(rfincr/2.)
7035	format(f8.3,' rod2pix ',f8.3,' rod2pix rmoveto')
c  put point on the diagonal
	write(9,7040) tend-(rflabel-(rfincr/2.))
c	write(*,7040) tend-(rflabel-(rfincr/2.))
7040	format(f8.3,' rod2pix 0 rlineto')
c go to the RH edge
c
	rod2=amin1(tend,(1.-(rflabel-(rfincr/2.))-tend))
c   based on the cutoff plane being defined by R1+R2+R3=1
	write(9,7115) rod2-(rflabel-(rfincr/2.))
c	write(*,7115) rod2-(rflabel-(rfincr/2.))
7115	format(' 0 ',f8.3,' rod2pix rlineto')
c  put point at edge of triangle; have to subtract the R2 value that we got to in 7035
	rod2p=amin1(tend,(1.-(rflabel-(rfincr/2.)))/2.)
	write(9,7120) (rod2p-tend),(rod2p-rod2)
c	write(*,7120) (rod2p-tend),(rod2p-rod2)
7120	format(f8.3,' rod2pix ',f8.3,' rod2pix rlineto')
c  continue to the diagonal
	write(9,7130) tx(k),ty(k)
c	write(*,7130) tx(k),ty(k)
7130	format(f8.3,' inch ',f8.3,' inch lineto closepath clip')
c  back to the corner of the triangle
c
	write(9,*) '/bsz ',deltax,'  def'
	do 100, i=1,mres
c  i indexes R1
	rod1=float(i)*deltax-(deltax/2.)
c  calculate the x coordinate based on a triangle with 2" side
	do 100, j=1,i
c  j indexes R2
	rod2=float(j)*deltax-(deltax/2.)
	write(9,6020) tx(k)+rod1,ty(k)+rod2
6020	format(f8.3,' inch ',f8.3,' inch box')
c  boxes are centered on each point, but clipping takes care of spill-over
c
	if(kindex.ne.99) then
	    tmp=amax1(1e-6,array(i,j,kindex))
	    rztmp=alog10(tmp)
	else
	    r3incr=tend/float(mres)
	    kkk=int(rflabel/r3incr)
c	    kkk=int(0.5+rflabel/r3incr)
	    kk=max0(kkk,1)
	    kk1=max0(kk+1,mresolve)
	    tmp=(float(kk-1)+0.5)*r3incr
	    rtmp=(float(kk)+0.5)*r3incr
	    weight1=(rflabel-tmp)
	    weight2=(rtmp-rflabel)
	    write(*,*) 'i,j,kkk',i,j,kkk,'  rflabel ',rflabel
	    write(*,*) 'kk, kk1, weight1,2   ', kk,kk1,weight1,weight2
	    if(kk1.gt.i) then
c  if the second z index is larger than the x index, use only the
c lower point
c
	        weight1=0.
	        weight2=1.
	    endif
	    if(weight1.lt.0.) then
	        weight1=0.
	        weight2=1.
	    endif
	    if(weight2.lt.0.) then
	        weight2=0.
	        weight1=1.
	    endif
	    rtmp=weight1+weight2
	    weight1=1.-abs(weight1)/rtmp
	    weight2=1.-weight1
	    rtmp=(weight1*array(i,j,kk)+weight2*array(i,j,kk1))
	    write(*,*) 'weight1/weight2  rtmp ',weight1,weight2,rtmp
	    if(rtmp.le.0.0) rtmp=1e-6
	    rztmp=-1.*alog10(rtmp)
c	    rztmp=rtmp
	endif
c
	rtmp=rztmp/scale
	tmp=amin1(1.0,rtmp)
	rztmp=amax1(-1.0,tmp)
	thue=((1.-rztmp)/2.9)
c  colour (hue) ranges from about red to blue
c	tbright=(tmp+1.)/2.
	tsat=amax1(0.1,(2.*abs(tmp)/scale))
c	tbright=4.*abs(tmp)/scale
	if(tsat.gt.1.0) tsat=1.0
	if(iqcolor.eq.0) then
	    write(9,6030) thue,tsat,1.0
6030	format(' gsave ',3f8.3,' sethsbcolor fill grestore ')
	elseif(iqcolor.eq.1) then
	    write(9,6031) thue
6031	format(' gsave ',f8.3,' setgray fill grestore ')
	endif
c  outline a square symbol, set color and grey level, fill
	if(rztmp.gt.0.0) write(9,6040) tx(k)+rod1,ty(k)+rod2
6040	format(' newpath',f8.3,' inch ',f8.3,' inch ',
     &	'0.08 0 360 arc gsave 0 setgray fill grestore')
c   put a small dot in the box if the value is < 0
100	continue
	write(9,*) ' grestore'
150	continue
c
	write(9,*) '/bsz 0.5 def'
c  sets boxsize to 0.5" for legend
c
	do 200, i=-4,4
	rtmp=float(4-i)*0.5
	tmp=float(i)/4.
	if(mod(i,2).eq.0) write(9,8010) 
     &	     tx(4)+rtmp,ty(4)+2.15,exp(-1.*scale*tmp*2.3026)
8010	format(f8.3,' inch ',f8.3,' inch moveto (',f6.2,') show')
	write(9,6020) tx(4)+rtmp,ty(4)+2.6
	thue=((tmp+1.)/2.9)
c  colour (hue) ranges from about red to blue
	tsat=amax1(0.1,(2.*abs(tmp)/scale))
c	tbright=4.*abs(tmp)/scale
	if(tsat.gt.1.0) tsat=1.0
	if(iqcolor.eq.0) then
	    write(9,6030) thue,tsat,1.0
	elseif(iqcolor.eq.1) then
	    write(9,6031) thue
	endif
c  arrange for bright colors unless close to ratio=1
c  outline a square symbol, set color and grey level, fill

200	continue
	return
	end
c
c ________________________________________
c
	subroutine defclipbox(xclip,yclip,nclip,side)
c  defines the clipping region
c  in this case a simple square box
c
	implicit none
	real xclip(100),yclip(100),side
	integer nclip
c
	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
c	do 6010, i=1,nclip
c		write(*,*) 'defclip: ',i,xclip(i),yclip(i)
c6010	continue
c
	return
	end
c
c ________________________________________
c
	subroutine defclip(xclip,yclip,nclip,radius)
c  defines the clipping region
c  in this case a circle
c
	real xclip(100),yclip(100),radius,angle,pi
	integer nclip
	real 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
c
c  second version, for use in REXGBS, WTS2POP
c  not for use in REX2D
c
	integer numsymm,numsamp
c
	common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c
	rad=57.29578
c	write(*,*) 'reading symmetry operators from "quat.symm.cubic"'
	open(unit=3,name='quat.symm.cubic',status='old')
	read(3,*)	! skip over title line
	read(3,*) numsymm	! read number of symmetry operators
	if(numsymm.gt.24) then
		numsymm=24
		write(*,*) 'too many symmetry operators, cutting back to 24!'
		else
c		write(*,*) 'reading ',numsymm,' symmetry operators'
		endif
c
	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
c
	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)
	real r1(3),r2(3),rout(3)
	real 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)
	real d(3),quat(4),rtmp,rnorm,rnorm2
c  converts Rodrigues vector to a quaternion 
c
	rtmp=d(1)**2+d(2)**2+d(3)**2
c	rnorm=sqrt(rtmp)
	rnorm2=sqrt(1.+rtmp)
	quat(4)=1./rnorm2
	quat(3)=d(3)/rnorm2
	quat(2)=d(2)/rnorm2
	quat(1)=d(1)/rnorm2
c  cute - this works for all rotation angles!
	return
	end
c
c
c  ____________________________
c
c
	subroutine q2rod(d,quat)
	real d(3),quat(4),rtmp,rnorm,rnorm2
c  converts a quaternion to a Rodrigues vector
c
	rnorm=sqrt(quat(1)**2+quat(2)**2+quat(3)**2+quat(4)**2)
	if(rnorm.lt.1e-5) stop 'Q2ROD: too small a quaternion!'
	if(abs(rnorm-1.0).gt.1.e-3) then
		write(*,*) 'Q2ROD: quat required normalization'
		do i=1,4
			quat(i)=quat(i)/rnorm
		enddo
	endif
c
	if(quat(4).gt.1.e-2) then
		d(1)=quat(1)/quat(4)
		d(2)=quat(2)/quat(4)
		d(3)=quat(3)/quat(4)
	else
c  we calculate sine(theta/2) from quat1-3
c  and scale by tangent/sine in order to make
c  the system deal with large numbers!
c
		rtmp=sqrt(quat(1)**2+quat(2)**2+quat(3)**2)
c  this will be close to unity
		stheta=asin(rtmp)
		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
c
	return
	end
c
c
c _____________________________
c
c
	subroutine presamp(qq,lindex,qresult)
	real qq(4),qresult(4)
	integer lindex
c
	common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c
c  algorithm for forming resultant quaternion
c  from applying a symmetry operator, QUATSYMM(n,lindex)
c  to the first quaternion/rotation
c
c  If the symm oper == O, and the input quaternion, QQ, is an active
c  rotation, then Qresult = (O x QQ) 
c  the "PRE" in the name refers to writing the operator before the
c  rotation/orientation in vector/tensor notation:  Q' = O x Q
c  or in conventional quaternion notation, Qresult = QQ € 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 postsymm(qq,lindex,qresult)
        real qq(4),qresult(4)
        integer lindex
c
c        include 'common.f'
        common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c
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 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
	function scalar(a,b)
c  return SCALAR PRODUCT of A and B
	real scalar,a(3),b(3)
	scalar=0.
	do 100, i=1,3
	scalar=scalar+a(i)*b(i)
100	continue
	return
	end
c
c ____________________________
c
      subroutine vecpro2(v1,v2,vout)
c  vector product
	real v1(3),v2(3),vout(3)
	vout(3)=v1(1)*v2(2)-v1(2)*v2(1)
	vout(1)=v1(2)*v2(3)-v1(3)*v2(2)
	vout(2)=v1(3)*v2(1)-v1(1)*v2(3)
	return
	end
c
c ________________________________________
c
        SUBROUTINE 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)
	REAL 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)
	real array(nnsize,nnsize)
	REAL SMF((-1*nnsize):nnsize),smct(nnsize)
c        REAL SMF(19,-36:36),SMCT(72)
	real 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)
c1000   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
	real array(nnsize,nnsize,nnsize)
	REAL SMF((-1*nnsize):nnsize),smct(nnsize)
c        REAL SMF(19,-36:36),SMCT(72)
	real rmax,ccnt
c
	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
	DO mm=1,nnsize     !  110
		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
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  smooth over the 2nd index
	DO mm=1,nnsize     !  110
		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
c  smooth over the 3rd index
	DO mm=1,nnsize     !  110
		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
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)
c
        RETURN
        END
c
c ________________________________________
c
	include 'psplot.txt'
c