program RF_points

c  compile with (on the new laptop)
c   /usr/local/bin/g77 -ffixed-line-length-none -g -o RFpoints_30Nov05 RFpoints_30Nov05.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  -------------------

c  please acknowledge the NSF-supported Mesoscale Interface Mapping Project
c  based in the Materials Science & Engineering Dept 
c  at Carnegie Mellon University, Pittsburgh

c  new program May 05 for plotting points in RF space

c  based on MAPE.F

	implicit none

	integer m,iqcolor,kcount,iq0,nval
	parameter (m=100)
	real array(m,m,m)
c  ARRAY contains the data to be plotted

	integer mresolve,ndata,jndex,index,i,j,k,ii
	real tend,pi,pi2,degrad,deltae,recipr3
	integer ndiscrete,itmp(3),icount
	real vmax,vmin,vavg
	character fname*50,fname2*50,fname3*50
	real rtmp(3),rodr(3),rfsymm(3),rresult(3),weight,scale
	real quat(4),qresult(4)
	integer ir3val,ir2val,ir1val,type
	logical goodpoint,single,repeat
	character  crys_sym_type,sam_sym_type
	integer numsymm,numsamp
	real quatsymm,quatsamp

	common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp

	integer index1,nbox,idelta,isymm
	real garray(m,m),r3,xsize,ysize,rnorm,rnormmin
	integer idummy,jdummy
	real qtmp(4),qqsym(4),qint(4),delta,r1,r2
	real voxvol(m,m,m),volsum,avgvol,weight_total
	integer badpoints,ijk,kk
	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
	real xcentre, ycentre, XX1, YY1, WIDTH, HEIGHT, RAD, ANG
	real red,green,blue
	logical FILL,inside,plotted
	real txnew,tynew
	integer iq_input,iq_weight,num_pts
	real wgt_max,wgt_min
	integer iq_type,iq_skip,iq_header,ulimit,iq_cols_skip
	real dummy(20)
	logical plottable
	integer iq_scale
	real rescale
	real delta_r3, stmp   !  spacing of RF plots in R3, -.3,-.2,-.1,0,.1,.2,.3
	integer lim1,lim2,lim3  !  for the negative rotation "do"
	real angle
	integer num_angles  !  for reading one vs two orientations
	real q_input(8),rf_input(6),qtmp2(4),qtwo(4,2),rtemp
	integer exit_status  !  for call system

c  CODE::
c  following variables control plotting options

	print*,'This program plots discrete weighted misorientations or orientations in Rodrigues-Frank space'
	print*

 3	print*,'Plot misorientations (=0), or orientations (=1), or misorientations by combination (=2)?'
c  '
	read(*,*) iq_type
	if(iq_type.lt.0.or.iq_type.gt.2) goto 3

 5	print*,'Is input in quaternions (=0) or RF-vectors (=1)?'
	print*,'or, in Bunge Euler angles in radians (=2), or ditto in degrees (=3)'
	print*,'or, in angle + [u v w] (=4)?'
	read(*,*) iq_input
	if(iq_input.lt.0.or.iq_input.gt.4) goto 5

 10	print*,'How many columns to skip before orientations [2 in .mdf files, 1 in Maniatty files] ?'
c '
	read(*,*) iq_skip
	if(iq_skip.lt.0.or.iq_skip.gt.20) goto 10
	if(iq_type.eq.2) then
	   print*,'By the way:'
	   print*,' the two sets of orientations must be in contiguous columns (no data in between!)'
	endif

 15	print*,'Unweighted input (=0) / column with weights (=1)?'
	read(*,*) iq_weight
	if(iq_weight.lt.0.or.iq_weight.gt.1) goto 15

 20	print*,'How many header lines [1 in .mdf files; 2 in Maniatty files] ?'
c'
	read(*,*) iq_header
	if(iq_skip.lt.0) goto 20

 25	print*,'How many data columns to skip over [2 in Maniatty files] ?'
	read(*,*) iq_cols_skip
	if(iq_cols_skip.lt.0.or.iq_cols_skip.gt.20) goto 25

	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
	txold=0.
	tyold=0.
	fill = .true.
	width = 0.02 * xsize
	height = width
	rad = width/2.
	ang = 0.

	tend=sqrt(2.)-1.
c  defines edge of the R-F space

	pi = 4.*atan(1.)
	pi2 = 2.*atan(1.)
	degrad = 180./pi
	recipr3 = 1./3.
	call textur4              !reads in symmetry operators
	write(*,*) 'Found ',numsymm,' xtal symms',' and '
	write(*,*) numsamp,' sample symms'

	tx(1) = 0.
	ty(1) = 0.
	tx(2) = 0.
	ty(2) = 2.1
	tx(3) = 0.
	ty(3) = 4.2
	tx(4) = 0.
	ty(4) = 6.3
	tx(5) = 3.5
	ty(5) = 0.
	tx(6) = 3.5
	ty(6) = 2.1
	tx(7) = 3.5
	ty(7) = 4.2
	tx(8) = 3.5
	ty(8) = 6.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.
	deltae = 2.*tend/float(m)
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

	icount=0
	volsum=0.0

c  now we read in the data
	write(*,*) 'Enter the filename for the input data, e.g. injy01.13:'
c$$$	if(iq_type.eq.0) then
c$$$	   write(*,*) 'An extension of .mdf is assumed'
c$$$	else
	   print*,'The extension (e.g. .txt) will be stripped off'
c$$$	endif
	read(*,'(a)') fname2
	index = 50
	do i = 50,1,-1
c		write(*,*) 'i,fname2(i:i) ',i,fname2(i:i)
	   if(fname2(i:i).eq.' ') index=i
	enddo
	index=index-1
	fname=fname2(1:index)
	fname2((index+1):(index+1))=char(0)
	write(*,*) 'fname2,index,fname'
	write(*,*) fname2,',   ',index,',   ',fname

	fname=fname2(1:index)
        single=.true.
c$$$	if(iq_type.eq.0) then
c$$$	   inquire(file=fname(1:index)//'.mdf',exist=repeat)
c$$$	   if(.not.repeat) then
c$$$	      write(*,*) 'Did NOT find ',fname(1:index)//'.mdf'
c$$$	      write(*,*) 'Will look for series of .mdf files'
c$$$	      single=.false.
c$$$	   else
c$$$	      write(*,*) 'Found ',fname(1:index)//'.mdf'
c$$$	   endif
c$$$	   
c$$$	   open(unit=2,name=fname(1:index)//'.mdf',status='old')
c$$$	else
	   open(unit=2,name=fname2,status='old')
c$$$	endif
	if(iq_header.gt.0) then
	   do i = 1,iq_header
	      read(2,*) 
	   enddo
	endif

	icount=0
	badpoints=0
	weight_total=0.0

	if(iq_type.eq.0.or.iq_type.eq.2) then
	   fname2=fname(1:index)//'_mdf_discrete.ps'
	   call newdev (fname2,index+17)
	   xsize=4.		!  change this value if you plot a different area in each section
	   ysize=4.
	   nbox=7
	   lim1 = 1
	   lim2 = -1
	   lim3 = -2
	else
	   fname2=fname(1:index)//'_odf_RF.ps'
	   call newdev (fname2,index+10)
	   xsize=2.
	   ysize=2.
	   nbox = 7
	   idelta = m/nbox
	   delta_r3 = 0.125
	   lim1 = 1
	   lim2 = 1
	   lim3 = 1
	endif
	call psinit(.true.)

	delta=float(m)/2.345/float(nbox)
	open(unit=3,file='RFpoints_output.txt',status='unknown')
	do i=1,nbox   ! this will be the loop over the sections
c		index1=idelta*(i-1)+3
c	   index1=(float(m)/2.)+(delta*float(i-1))
c	   r3=float(index1-(m/2))*deltae

	   if(iq_type.eq.0.or.iq_type.eq.2) then
	      r3 = 0.05 * float(i-1) * 7. / float(nbox)	!  works best for nbox = 7
	   else
	      index1 = idelta*(i-1)+3
c	      r3 = float(index1-(m/2))*deltae
	      r3 = float(i-4)*delta_r3
	   endif

	   write(*,*) 'Section at R3 = ',r3
	   txdiff=tx(i)-txold
	   tydiff=ty(i)-tyold
	   txold=tx(i)
	   tyold=ty(i)
	   print*,'before CALL PLOT, txold,tyold= ',txold,tyold
	   call plot(txdiff,tydiff,-3)
	   write(*,*)
	   write(*,*) 'Section no. ',i
	   if(iq_type.eq.0.or.iq_type.eq.2) then
	      call keksym(-.2+(xsize/2.),0.15+(ysize/2.)
     &						,.15,trichar(i),0.,1,0)
	      call keksym(0.3 +(xsize/2.),1.3+(ysize/2.),.15,'R3 =',0.,4,0)
	      call keknum(0.3 +(xsize/2.),1.05+(ysize/2.),0.15,r3,0.,3,0)
	      call border(xsize,ysize,0000,0000,4,2,4,2)
	   else
	      call keksym(-1.2+(xsize/2.),-0.15+(ysize/2.)
     &						,.15,trichar(i),0.,1,0)
	      call keksym(1.1 +(xsize/2.),0.75+(ysize/2.),.15,'R3 =',0.,4,0)
	      call keknum(1.1 +(xsize/2.),0.5+(ysize/2.),0.15,r3,0.,3,0)
	      call border(xsize,ysize,1111,0000,2,1,2,1)
	   endif

	   call GSAV
c	   call defclipRFtri(xclip,yclip,nclip,xsize,r3)
c	   call clipbox(xclip,yclip,nclip)
c	   call confill(garray,m,m,m,xsize,ysize,cval(1),greyval,
c     &        nval,0,0.)
c	   call conrec(garray,m,m,m,xsize,ysize,cval(1),nval)
c	   call clip		!  clips off the FZ triangle
c	   call GREST
	   if(iq_type.eq.0.or.iq_type.eq.2) then
	      call defclipbox(xclip,yclip,nclip,xsize)
	   else
	      call defclip_squ(xclip,yclip,nclip,xsize)
	   endif
	   call DRWCRV (xclip, yclip, nclip, .03, .true.)
c  draws the full triangle around the section
	enddo  !  do i=1,nbox
	txdiff = txold*-1.
	tydiff = tyold*-1.
	call plot(txdiff,tydiff,-3)  !  get us back to the start
	if(iq_type.eq.1) then
	   call keksym(0.5,9.3,.15,'input:',0.,6,0)
	   call keksym(0.3,9.0,.15,fname(1:index),0.,index,0)
	else
	   call keksym(0.5,0.9,.15,'input:',0.,6,0)
	   call keksym(0.3,0.6,.15,fname(1:index),0.,index,0)
	endif


	txold=0.
	tyold=0.
c  begin the loop to read in points
	wgt_max = 0.
	wgt_min = 1.e30
	if(iq_input.eq.0.or.iq_input.eq.4) then
	   num_angles = 4
	else
	   num_angles = 3
	endif
	if(iq_type.eq.2) then
	   num_angles = num_angles*2
	endif

 	do i=1,10000000		!  label 580
	   if(mod(i,10).eq.0) write(*,*) 'reading line no. ',i
c     goodpoint=.true.
	   isymm=0
	   weight = 0.
	   
	   if(iq_input.eq.0.or.iq_input.eq.4) then
	      if(iq_weight.eq.0) then
		 if(iq_skip.gt.0) then
		    read(2,*,err=580,end=30) (dummy(ii),ii=1,iq_skip)
     $                   ,(q_input(ii),ii=1,num_angles)
		 else
		    read(2,*,err=580,end=30) (q_input(ii),ii=1,num_angles)
		 endif
		 weight = 1.
	      else
		 if(iq_skip.gt.0.and.iq_cols_skip.gt.0) then
		    read(2,*,err=580,end=30) (dummy(ii),ii=1,iq_skip)
     $                   ,(q_input(ii),ii=1,num_angles),(dummy(ii)
     $                   ,ii=1,iq_cols_skip),weight
		 elseif(iq_skip.gt.0.and.iq_cols_skip.le.0) then
		    read(2,*,err=580,end=30) (dummy(ii),ii=1,iq_skip)
     $                   ,(q_input(ii),ii=1,num_angles),weight
		 elseif(iq_skip.le.0.and.iq_cols_skip.le.0) then
		    read(2,*,err=580,end=30) (q_input(ii),ii=1,num_angles),weight
		 endif
	      endif
	   endif

	   if(iq_input.ge.1.and.iq_input.le.3) then
	      if(iq_weight.eq.0) then
		 if(iq_skip.gt.0) then
		    read(2,*,err=580,end=30) (dummy(ii),ii=1,iq_skip)
     $                   ,(rf_input(ii),ii=1,num_angles)
		 else
		    read(2,*,err=580,end=30) (rf_input(ii),ii=1,num_angles)
		 endif
		 weight = 1.
	      else
		 if(iq_skip.gt.0.and.iq_cols_skip.gt.0) then
		    read(2,*,err=580,end=30) (dummy(ii),ii=1,iq_skip)
     $                   ,(rf_input(ii),ii=1,num_angles)
     $                   ,(dummy(ii),ii=1,iq_cols_skip),weight
		 elseif(iq_skip.gt.0.and.iq_cols_skip.le.0) then
		    read(2,*,err=580,end=30) (dummy(ii),ii=1,iq_skip)
     $                   ,(rf_input(ii),ii=1,num_angles),weight
		 elseif(iq_skip.le.0.and.iq_cols_skip.le.0) then
		    read(2,*,err=580,end=30) (rf_input(ii),ii=1,num_angles),weight
		 endif
	      endif
c$$$  if(iq_input.eq.1) then
c$$$  call rod2q(rodr,qtmp)
c$$$  elseif(iq_input.eq.2) then
c$$$  call quatB(rf_input(1),rf_input(2),rf_input(3),qtmp)
c$$$  elseif(iq_input.eq.3) then
c$$$  call quatB(rf_input(1)/degrad,rf_input(2)/degrad,rf_input(3)/degrad,qtmp)
c$$$  endif
c     orientations do not have to be calculated here - see below
	   endif

	   if(weight.gt.wgt_max) wgt_max = weight
	   if(weight.lt.wgt_min) wgt_min = weight
 580    enddo                   !  i =
 30     num_pts = i-1
        print*,'no. of points = ',num_pts
        print*,'min, max weight = ',wgt_min,wgt_max
c$$$  print*,'Apply a scaling factor [0 = YES]'
c$$$  read(*,*) iq_scale
c$$$  if(iq_scale.eq.0) then
c$$$	      print*,'Enter scaling factor'
c$$$	      read(*,*) rescale
c$$$	   endif

	   if(iq_type.eq.1) then
	   call keksym(5.,9.3,.15,'max. weight: ',0.,13,0)
	   call keknum(6.7,9.3,0.15,wgt_max,0.,2,0)
	   call keksym(5.,9.,.15,'min. weight: ',0.,13,0)
	   call keknum(6.7,9.,0.15,wgt_min,0.,2,0)
	   else
	   call keksym(5.,0.9,.15,'max. weight: ',0.,13,0)
	   call keknum(6.7,0.9,0.15,wgt_max,0.,2,0)
	   call keksym(5.,0.6,.15,'min. weight: ',0.,13,0)
	   call keknum(6.7,0.6,0.15,wgt_min,0.,2,0)
	   endif

	   rewind(2)
	   if(iq_header.gt.0) then
	      do i = 1,iq_header
		 read(2,*) 
	      enddo
	   endif

	   do i = 1,num_pts
	      if(mod(i,1).eq.0) then
c		 print*
		 write(*,*) '\015 reading line no. ',i
	      endif
c     goodpoint=.true.
	      isymm=0
	      
	      weight = 0.
	      if(iq_input.eq.0.or.iq_input.eq.4) then
		 if(iq_weight.eq.0) then
		    if(iq_skip.gt.0) then
		       read(2,*) (dummy(ii),ii=1,iq_skip),(q_input(ii),ii=1,num_angles)
		    else
		       read(2,*) (q_input(ii),ii=1,num_angles)
		    endif
		    weight = 1.
		 else
		    if(iq_skip.gt.0.and.iq_cols_skip.gt.0) then
		       read(2,*) (dummy(ii),ii=1,iq_skip),(q_input(ii),ii=1,num_angles),(dummy(ii),ii=1,iq_cols_skip),weight
		    elseif(iq_skip.gt.0.and.iq_cols_skip.le.0) then
		       read(2,*) (dummy(ii),ii=1,iq_skip),(q_input(ii),ii=1,num_angles),weight
		    elseif(iq_skip.le.0.and.iq_cols_skip.le.0) then
		       read(2,*) (q_input(ii),ii=1,num_angles),weight
		    endif
		 endif  !  if(iq_weight
		 do ii = 1,4
		    qtmp(ii) = q_input(ii)
		    if(iq_type.eq.2) then
		       qtmp2(ii) = q_input(ii+4)
		    endif
		 enddo  !  xfer the numbers
	      endif
	      
	      if(iq_input.ge.1.and.iq_input.le.3) then
		 if(iq_weight.eq.0) then
		    if(iq_skip.gt.0) then
		       read(2,*) (dummy(ii),ii=1,iq_skip),(rf_input(ii),ii=1,num_angles)
		    else
		       read(2,*) (rf_input(ii),ii=1,num_angles)
		    endif
		    weight = 1.
		 else
		    if(iq_skip.gt.0.and.iq_cols_skip.gt.0) then
		       read(2,*) (dummy(ii),ii=1,iq_skip),(rf_input(ii),ii=1,num_angles),(dummy(ii),ii=1,iq_cols_skip),weight
		    elseif(iq_skip.gt.0.and.iq_cols_skip.le.0) then
		       read(2,*) (dummy(ii),ii=1,iq_skip),(rf_input(ii),ii=1,num_angles),weight
		    elseif(iq_skip.le.0.and.iq_cols_skip.le.0) then
		       read(2,*) (rf_input(ii),ii=1,num_angles),weight
		    endif
		 endif
		 if(iq_input.eq.1) then
		    do ii = 1,3
		       rodr(ii) = rf_input(ii)
		    enddo	!  xfer the numbers
		    call rod2q(rodr,qtmp)
		    if(iq_type.eq.2) then
		       do ii = 1,3
			  rodr(ii) = rf_input(ii+3)
		       enddo	!  xfer the numbers
		       call rod2q(rodr,qtmp2)
		    endif
		 elseif(iq_input.eq.2) then
		    call quatB(rf_input(1),rf_input(2),rf_input(3),qtmp)
		    if(iq_type.eq.2) then
		       call quatB(rf_input(4),rf_input(5),rf_input(6),qtmp2)
		    endif
		 elseif(iq_input.eq.3) then
		    call quatB(rf_input(1)/degrad,rf_input(2)/degrad,rf_input(3)/degrad,qtmp)
		    if(iq_type.eq.2) then
		       call quatB(rf_input(4)/degrad,rf_input(5)/degrad,rf_input(6)/degrad,qtmp2)
		    endif
		 elseif(iq_input.eq.4) then
c  angle + axis input
                    do ijk = 1,3
                       rodr(ijk) = q_input(ijk+1)
                    enddo
                    call donorm(rodr)
                    rtemp = sin(0.5*q_input(1)/degrad)
                    do ijk = 1,3
                       qtmp(ijk) =  rtemp * rodr(ijk)
                    enddo
                    qtmp(4) = cos(0.5*q_input(1)/degrad)
		    if(iq_type.eq.2) then
                       do ijk = 1,3
                          rodr(ijk) = q_input(ijk+5)
                       enddo
                       call donorm(rodr)
                       rtemp = sin(0.5*q_input(5)/degrad)
                       do ijk = 1,3
                          qtmp(ijk) =  rtemp * rodr(ijk)
                       enddo
                       qtmp(4) = cos(0.5*q_input(5)/degrad)
		    endif
		 endif
	      endif

575	      format(4f12.4)

              print*,'Quat input = ',(qtmp(ijk),ijk=1,4)
	      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

	      if(iq_type.eq.2) then
		 if(weight.le.0.0) weight=0.
		 rnorm=sqrt(qtmp2(1)**2+qtmp2(2)**2+qtmp2(3)**2+qtmp2(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: quat2 required normalization'
		    do ijk=1,4
		       qtmp2(ijk)=qtmp2(ijk)/rnorm
		    enddo
		 endif
		 do ijk=1,4
		    qtwo(ijk,1)= qtmp(ijk)
		    qtwo(ijk,2)= qtmp2(ijk)
 		 enddo
		 call comquat(qtwo,qtmp)
c  makes new quat with misorientation from 2 orientations		 
	      endif  !  if (iq_type
	      
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

	      plotted = .false.
	      ulimit = 1
	      if(iq_type.eq.0.or.iq_type.eq.2) ulimit = numsymm
c	      if(iq_type.eq.1) ulimit = 1 !  force triclinic sample symmetry for now
	      if(iq_type.eq.1) ulimit = numsamp
	      do j = 1,numsymm
		 if(iq_type.eq.0.or.iq_type.eq.2) then
		    call presymm(qtmp,j,qint)
		 else
		    call postsymm(qtmp,j,qint)
		 endif
		 do k = 1,ulimit
		    if(iq_type.eq.0.or.iq_type.eq.2) then
		       call postsymm(qint,k,qresult)
		    else
		       call presamp(qint,k,qresult)
		    endif
		    do kk = lim1,lim2,lim3     ! generate negative rotation
		       if(kk.lt.0) then
			  do ijk=1,3
			     qresult(ijk)=qresult(ijk)*float(kk)
			  enddo
		       endif
c      write(*,"('symm 1 2, qresult = ',2i5,4(1x,f8.5))") j,k,qresult
		       goodpoint=.true.
		       angle = degrad*2.*acos(qresult(4))
c		       if(abs(qresult(4)).gt.0.45) 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, rotation angle = ',4(1x,g11.4))") rresult,angle
		       
c  and now to plot it
c		       ir3val = 1 + int((float(nbox) * rresult(3) / recipr3) - 0.01)
			  plottable = .false.
			  if(iq_type.eq.0.or.iq_type.eq.2) then  !  misorientations
			     ir3val = 1 + int((rresult(3) / 0.05) * (float(nbox)/7.))
			     if(ir3val.ge.0.and.ir3val.le.nbox
     1		             .and.rresult(2).le.rresult(1)
     2		             .and.rresult(2).ge.rresult(3)
     3		             .and.rresult(1).ge.0.
     4			     .and.rresult(2).ge.0.
     5			     .and.rresult(1).le.tend
     6			     .and.rresult(3).ge.0.
     7		             .and.(.not.plotted)) then
				plottable = .true.
			     endif
			  elseif(iq_type.eq.1) then    !  orientations
c  the "0.10" is the spacing of the boxes in R3
c			     ir3val = 3 + int((rresult(3) / 0.10) * (float(nbox)/8.))
			     stmp = rresult(3) / delta_r3
			     if(stmp.gt.0.) then
				ir3val = 4 + int(stmp + 0.5)
			     else
				ir3val = 4 + int(stmp - 0.5)
			     endif
c			     print*,' ir3val = ',ir3val
			     if(ir3val.lt.1) ir3val = 1
			     if(ir3val.gt.nbox) ir3val = nbox
			     if(
     1		             rresult(3).ge.(-1.*tend)
     2	                     .and.rresult(3).le.tend
     3		             .and.rresult(1).ge.(-1.*tend)
     4		             .and.rresult(2).ge.(-1.*tend)
     5	                     .and.rresult(1).le.tend
     6	                     .and.rresult(2).le.tend
     7			     .and.(((.not.plotted).and.(iq_type.eq.0.or.iq_type.eq.2)).or.iq_type.eq.1)) then
c  if you want sample symmetries other than triclinic, then must change the FZ here
				plottable = .true.
			     endif  !  test within a box
			  endif  !  misorient vs orientation

			  if(plottable) then
c			  print*
c			  print*,'Plottable: rresult: ',rresult
c			  print*,'ir3val,tx,ty ',ir3val,tx(ir3val),ty(ir3val)
c			  txnew = tx(ir3val)+(xsize/2.)
c			  tynew = ty(ir3val)+(xsize/2.)
			  txnew = tx(ir3val) !  back to standard approach in order
			  tynew = ty(ir3val) !  to use defclipRFtri
			  txdiff=txnew-txold
			  tydiff=tynew-tyold
			  txold=txnew
			  tyold=tynew
c			  print*,'before CALL PLOT, txold,tyold= ',txold,tyold
			  call plot(txdiff,tydiff,-3) !  get into the correct triangle
			  red = 0.99
			  green = 0.1
			  blue = 0.1
			  CALL SETCOLR (RED, GREEN, BLUE)
			  xcentre = xsize/2. + (xsize/2. * rresult(1) / tend)
			  ycentre = xsize/2. + (xsize/2. * rresult(2) / tend)
			  xx1 = xcentre - (width/2.)
			  yy1 = ycentre - (width/2.)
			  r3 = 0.05 * float(ir3val-1) * 7. / float(nbox) !  works best for nbox = 7
c       call testRFtri(xcentre,ycentre,xsize,r3,inside)
			  if(iq_weight.eq.1) then
			     width = 0.02 * xsize * amax1(0.8,(3.5 * (weight-wgt_min) / (wgt_max-wgt_min)))
			     height = width
			     rad = width/3.
			     xx1 = xcentre - (width/2.)
			     yy1 = ycentre - (width/2.)
			     red = 2. * amax1(0.,(((weight-wgt_min) / (wgt_max-wgt_min))-0.5))
			     green = 1. - (2. * abs(0.5 - ((weight-wgt_min) / (wgt_max-wgt_min))))
c       blue = 0.99 * (1. - (weight / wgt_max))
			     blue = 1.0 - (2. * amin1(0.5,((weight-wgt_min) / (wgt_max-wgt_min))))
			     CALL SETCOLR (RED, GREEN, BLUE)
                          else
                             width = 0.02 * xsize
			     height = width
			     rad = width/3.
			     xx1 = xcentre - (width/2.)
			     yy1 = ycentre - (width/2.)
                             red = 0.2
                             green = 0.2
                             blue = 0.2
			     CALL SETCOLR (RED, GREEN, BLUE)
			  endif
			  inside = .false.
			  inside = ((rresult(1)+rresult(2)+rresult(3)).le.1.)
			  if(inside) then
			     plotted = .true.
			     CALL RRECT (XX1, YY1, WIDTH, HEIGHT, RAD, ANG, FILL)
c			     print*,'XX1, YY1, WIDTH, HEIGHT, RAD, ANG, FILL'
c			     print*,XX1, YY1, WIDTH, HEIGHT, RAD, ANG, FILL
			  endif
		       endif	!  if(ir3val
		       
c		    endif	!  from qresult(4) test
		 enddo		!  inversion 
	      enddo		!  2nd crystal symmetry loop
	   enddo		!  1st xtal symm
	enddo                   !  each point
     
	call plotnd

	if(iq_type.eq.0.or.iq_type.eq.2) then
	   fname3 = fname(1:index)//'_mdf_discrete.jpg'
	else
	   fname3 = fname(1:index)//'_odf_RF.jpg'
	endif

	call system('convert '//fname2//' '//fname3,exit_status)
        print*,'exit_status from CONVERT = ',exit_status
        if(exit_status.eq.0) then
           call system('rm '//fname2)
        endif

	Close(3)  !  end of plotting
	close(2)  !  end of plotting

	call exit(0)   !  successful ending

	end

c
c ________________________________________
c
	subroutine testRFtri(x,y,side,r3,inside)
c  tests if a point is within the clipping region
c  in this case a TRIANGLE

	implicit none

	integer i
	real xclip(10),yclip(10),side,tend,r3,ax,ay,bx,by
	integer nclip
	logical corner,corner2,inside
	real x,y,slope

c  CODE::

	inside = .true.
	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 Ã2-1

	xclip(1) = side/2.+(r3/tend*side/2.)
	yclip(1) = side/2.+(r3/tend*side/2.)
	xclip(2) = side
	yclip(2) = yclip(1)
	if(y.gt.x) inside = .false.
	if(x.lt.xclip(1)) inside = .false.
	if(x.gt.xclip(2)) inside = .false.
	if(y.lt.yclip(1)) inside = .false.
	if(corner) then
	   if(corner2) then
c  this is a simple triangle
	      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
c  this has a notch cut off the top
	      xclip(3)=side
	      yclip(3)=side/2.+(by/tend*side/2.)
c  3rd point is up the RHS of the triangle
	      xclip(4)=side/2.+(ax/tend*side/2.)
	      yclip(4)=side/2.+(ay/tend*side/2.)
	      slope = (yclip(3)-yclip(4))/(xclip(4)-xclip(3))
	      if((yclip(4)-y)/(x-xclip(4)).le.slope) inside = .false.
	      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 i=1,nclip
c	   write(*,"('defclip: ',i4,2(2x,f6.2))") i,xclip(i),yclip(i)
c	enddo

	return
	end

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)=25.
	cval(2)=50.
	cval(3)=100.
	cval(4)=150.
	cval(5)=200.
	cval(6)=400.
	cval(7)=800.

	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
	subroutine defclip_squ(xclip,yclip,nclip,side)
c  defines the clipping region
c  in this case a simple square box

	implicit none
	real xclip(100),yclip(100),side
	integer nclip

	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
	subroutine defclipbox(xclip,yclip,nclip,side)
c  defines the clipping region
c  in this case a TRIANGLE
c
	implicit none
	real xclip(100),yclip(100),side
	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 defclipRFtri(xclip,yclip,nclip,side,r3)
c  defines the clipping region
c  in this case a TRIANGLE

	implicit none

	integer i
	real xclip(100),yclip(100),side,tend,r3,ax,ay,bx,by
	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 Ã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 i=1,nclip
c	   write(*,"('defclip: ',i4,2(2x,f6.2))") i,xclip(i),yclip(i)
c	enddo

	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
	subroutine quatB(p1,p,p2,q)
c  convert Bunge Euler angles (radians) to quaternion; 
c  direct conversion from angles, see Altmann's book
c  the rotation is a vector transformation (active rotation)
c
	real p1,p,p2,q(4)
	double precision c1,c,c2,s1,s,s2,d(3,3),tmp(4),sin2,cos2,rtmp,rnorm
c
c
	PI=3.14159265
c  form cosine, sine of Phi, and sum & diff of phi1, phi2
	S=DSIN(0.5d0*P)
	C=DCOS(0.5d0*P)
	S1=DSIN(0.5d0*(P1-P2))
	C1=DCOS(0.5d0*(P1-P2))
	S2=DSIN(0.5d0*(P1+P2))
	C2=Dcos(0.5d0*(P1+P2))
c      write(*,*) 's1,c1,s,c,s2,c2'
c      write(*,*) s1,c1
c      write(*,*) s,c
c      write(*,*) s2,c2
	q(1)=s*c1
	q(2)=s*s1
	q(3)=c*s2
	q(4)=c*c2
	return
	end

c
c _____________________________
c
c
	subroutine comquat(qq,qresult)
	real qq(4,2),qresult(4)
c
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion as Q1 x Q2^-1
c  where Q2 follows Q1, i.e. is the second rotation;
c  note change of signs to get inverse of second orientation
c
	qresult(1)=qq(1,1)*qq(4,2)-qq(4,1)*qq(1,2)
     &	+qq(2,1)*qq(3,2)-qq(3,1)*qq(2,2)
	qresult(2)=qq(2,1)*qq(4,2)-qq(4,1)*qq(2,2)
     &	+qq(3,1)*qq(1,2)-qq(1,1)*qq(3,2)
	qresult(3)=qq(3,1)*qq(4,2)-qq(4,1)*qq(3,2)
     &	+qq(1,1)*qq(2,2)-qq(2,1)*qq(1,2)
	qresult(4)=qq(4,1)*qq(4,2)+qq(1,1)*qq(1,2)
     &	+qq(2,1)*qq(2,2)+qq(3,1)*qq(3,2)
c
c	write(*,*) 'qresult ',qresult
c
	return
	end
c
c ------------------------------------------------
c
	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

	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)

	implicit none

	real d(3),quat(4),rtmp,rnorm,rnorm2
	integer i
	real 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)
	real 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)
	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 presymm(qq,lindex,qresult)
	real qq(4),qresult(4)
	integer lindex
c
c	include 'common.f'
	common/a1symm/ quatsymm(4,48),numsymm,quatsamp(4,4),numsamp
c       PI=3.14159265
c
c  algorithm for forming resultant quaternion
c  from applying a symmetry operator, QUATSYMM(n,lindex)
c  to the first quaternion/rotation
c
c  If the symm oper == O, then
c  Qresult = (O x QQ) where QQ= Q1 x Q2^-1
c  i.e. Qresult = (O x Q1) x Q2^-1
c  note change of signs to get inverse of first orientation
c
	if(lindex.gt.numsymm) stop 'error in presymm, lindex>numsymm'
	if(lindex.lt.1) stop 'error in presymm, lindex<1'
	qresult(1)=quatsymm(4,lindex)*qq(1)+quatsymm(1,lindex)*qq(4)
     &	-quatsymm(3,lindex)*qq(2)+quatsymm(2,lindex)*qq(3)
	qresult(2)=quatsymm(4,lindex)*qq(2)+quatsymm(2,lindex)*qq(4)
     &	-quatsymm(1,lindex)*qq(3)+quatsymm(3,lindex)*qq(1)
	qresult(3)=quatsymm(4,lindex)*qq(3)+quatsymm(3,lindex)*qq(4)
     &	-quatsymm(2,lindex)*qq(1)+quatsymm(1,lindex)*qq(2)
	qresult(4)=quatsymm(4,lindex)*qq(4)-quatsymm(1,lindex)*qq(1)
     &	-quatsymm(2,lindex)*qq(2)-quatsymm(3,lindex)*qq(3)
c
c	write(*,*) 'qresult ',qresult
c
	return
	end
c
c _____________________________
c
c
        subroutine postsymm(qq,lindex,qresult)
        real qq(4),qresult(4)
        integer lindex
c
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 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
	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 ________________________________________

      subroutine donorm(a)

      implicit none

      real a(3),b(3),rsum,rnorm
      integer i

c     CODE::
      rsum=0.
      do 100, i=1,3
         rsum=rsum+a(i)**2
 100  continue
      rnorm=sqrt(rsum)
      if(rnorm.gt.1.e-6) then
         do 200, i=1,3
            b(i)=a(i)/rnorm
            a(i)=b(i)
 200     continue
      endif
      return
      end

c     +++++++++++++++++
      include '/code/paul.lee.codes/psplot.txt'