      subroutine prepincoh(m_rlv,zeta,eta,incohmon,incohana,
     &  hh,kk,ll,circsize,ki,kf,pn,constq,h,flagincoh,absent,replot, 
     &  nincohmon, nincohana, flagincohmon,flagincohana)

C  Spurions from incoherent processes.

      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
c      include 'spur_color.inc'
      include 'spur_scan.inc'

c...	common blocks
      real tol, etol		!! spatial and energy tolerances       
      common /spur/ tol, etol
      
c...	passed variables
      integer m_rlv
      real zeta(m_rlv_max,m_rlv_max)
      real eta(m_rlv_max,m_rlv_max)
      real incohmon(2),incohana(2)
      integer hh(m_rlv_max,m_rlv_max)
      integer kk(m_rlv_max,m_rlv_max)
      integer ll(m_rlv_max,m_rlv_max)
      real circsize,ki(2),kf(2),h(3)
      logical constq,flagincoh,replot
      logical absent(m_rlv_max,m_rlv_max)
      integer pn, nincohmon, nincohana
      logical flagincohmon,flagincohana
      
c...	local variables
      integer i,j,k
c      logical flagincohmon,flagincohana
      integer sp(3),len
      character*10 spurtitle
      logical firstmatch
      integer iincoh, ii, jj 

c      real temp1, temp2
      
C  Check to see if there is a reciprocal lattice match.


c      write (*,*) 'Incoherent spurions matches at:'
c      write (unitspur,*) 'Incoherent spurions matches at:'
      
      flagincohmon=.false.
      flagincohana=.false.
cMoved      nincohmon = 0
cMoved      nincohana = 0
      do i=1,m_rlv
        do j=1,m_rlv

          if ((abs(zeta(i,j)-incohmon(1)).lt.tol).and.
     &      (abs(eta(i,j)-incohmon(2)).lt.tol).and.
     &       .not.absent(i,j)) then
            flagincohmon=.true.
	    nincohmon = nincohmon + 1
	    
            if (replot) goto 687

            sp(1)=hh(i,j)
            sp(2)=kk(i,j)
            sp(3)=ll(i,j)
            call title(sp,spurtitle,len)
            write(*,'(4a,i2,$)') 'Incoh. spurion from ',
     &        'MON at ',spurtitle(1:len),
     &        ', scan point#',pn
            write(unitspur,'(4a,i2,$)') 'Incoh. spurion from ',
     &        'MON at ',spurtitle(1:len),
     &        ', scan point#',pn
c            write(8,'(4a,i2,$)') 'Incoh. spurion from ',
c     &        'MON at ',spurtitle(1:len),
c     &        ', point#',pn
            if (constq) then
               write(*,'(a,f6.2)') ' E=',
     &              k2_to_E*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
               write(unitspur,'(a,f6.2)') ' E=',
     &              k2_to_E*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
c02Apr               write(*,'(a,f6.2)') ' E=',
c02Apr     &              2.072*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
c02Apr               write(unitspur,'(a,f6.2)') ' E=',
c02Apr     &              2.072*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
c               write(8,'(a,f6.2)') ' E=',
c     &              2.072*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
            else
               write(*,637) ' at Q=(',(h(k),k=1,3),')'
               write(unitspur,637) ' at Q=(',(h(k),k=1,3),')'
c               write(8,637) ' at Q=(',(h(k),k=1,3),')'
 637           format(a,f6.2,1x,f6.2,1x,f6.2,a)
            endif

	    lab_incmon(nincohmon) = spurtitle
	    ispot_incmon(nincohmon) = i
	    jspot_incmon(nincohmon) = j
	    scanpt_incmon(nincohmon) = pn
          endif

 687      if ((abs(zeta(i,j)-incohana(1)).lt.tol).and.
     &       (abs(eta(i,j)-incohana(2)).lt.tol).and.
     &        .not.absent(i,j)) then
            flagincohana=.true.
	    nincohana = nincohana + 1

            if (replot) goto 697

            sp(1)=hh(i,j)
            sp(2)=kk(i,j)
            sp(3)=ll(i,j)
            call title(sp,spurtitle,len)
            write(*,'(4a,i2,$)') 'Incoh. spurion from ',
     &        'ANA at ',spurtitle(1:len),
     &        ', scan point#',pn
            write(unitspur,'(4a,i2,$)') 'Incoh. spurion from ',
     &        'ANA at ',spurtitle(1:len),
     &        ', scan point#',pn
c            write(8,'(4a,i2,$)') 'Incoh. spurion from ',
c     &        'ANA at ',spurtitle(1:len),
c     &        ', point#',pn
            if (constq) then
               write(*,'(a,f6.2)') ' E=',
     &              k2_to_E*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
               write(unitspur,'(a,f6.2)') ' E=',
     &              k2_to_E*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
c02Apr               write(*,'(a,f6.2)') ' E=',
c02Apr     &              2.072*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
c02Apr               write(unitspur,'(a,f6.2)') ' E=',
c02Apr     &              2.072*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
c               write(8,'(a,f6.2)') ' E=',
c     &              2.072*(ki(1)**2+ki(2)**2-kf(1)**2-kf(2)**2)
            else
               write(*,637) ' at Q=(',(h(k),k=1,3),')'
               write(unitspur,637) ' at Q=(',(h(k),k=1,3),')'
c               write(8,637) ' at Q=(',(h(k),k=1,3),')'
            endif
	    lab_incana(nincohana) = spurtitle
	    ispot_incana(nincohana) = i
	    jspot_incana(nincohana) = j
	    scanpt_incana(nincohana) = pn
          endif

 697    enddo
      enddo
      if (flagincohana.or.flagincohmon) flagincoh=.true.


c...	Load up "rankings"
	iincoh = 0
	do ii = 1, m_rlv
	  do jj = 1, m_rlv
	    firstmatch = .true.
	    do i = 1, nincohmon
	      if ((ii .eq. ispot_incmon(i))
     &	      .and. (jj .eq. jspot_incmon(i)) ) then
		if (firstmatch) then
	          iincoh = iincoh + 1
		  firstmatch = .false.
		endif
	        rankincmon(i) = iincoh
	      endif
	    enddo
	  enddo
	enddo
	n_rankincmon = iincoh


	iincoh = 0
	do ii = 1, m_rlv
	  do jj = 1, m_rlv
	    firstmatch = .true.
	    do i = 1, nincohana
	      if ((ii .eq. ispot_incana(i))
     &	      .and. (jj .eq. jspot_incana(i)) ) then
		if (firstmatch) then
	          iincoh = iincoh + 1
		  firstmatch = .false.
		endif
	        rankincana(i) = iincoh
	      endif
	    enddo
	  enddo
	enddo
	n_rankincana = iincoh


      return
      end

c============================================================================
      subroutine plotincoh(m_rlv,zeta,eta,incohmon,incohana,
     &  hh,kk,ll,circsize,ki,kf,pn,constq,h,flagincoh,absent,replot, 
     &  nincohmon, nincohana,  flagincohmon,flagincohana)

C  Spurions from incoherent processes.

      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_color.inc'

c...	common blocks
      real tol, etol		!! spatial and energy tolerances       
      common /spur/ tol, etol
      
c...	passed variables
      integer m_rlv
      real zeta(m_rlv_max,m_rlv_max)
      real eta(m_rlv_max,m_rlv_max)
      real incohmon(2),incohana(2)
      integer hh(m_rlv_max,m_rlv_max)
      integer kk(m_rlv_max,m_rlv_max)
      integer ll(m_rlv_max,m_rlv_max)
      real circsize,ki(2),kf(2),h(3)
      logical constq,flagincoh,replot
      logical absent(m_rlv_max,m_rlv_max)
      integer pn, nincohmon, nincohana
      logical flagincohmon,flagincohana
      
c...	local variables
c      integer i,j,k
c      integer sp(3),len
c      character*10 spurtitle

      real temp1, temp2
      

c...	Now plot
      call pgsfs(1)
      call pgsls(1)

      if (flagincohmon) then
        call pgsci(color_cyan)   ! cyan
c02Mar        call pgsci(color_magenta)   ! pink/purple
c01Dec        call pgsci(6)   ! pink/purple
	temp1 = 1.1*circsize
        call pgcirc(incohmon(1),incohmon(2), temp1)
c01Nov        call pgcirc(incohmon(1),incohmon(2),1.1*circsize)
	temp1 = incohmon(1)+kf(1)
	temp2 = incohmon(2)+kf(2)
        call pgarro(temp1, temp2, incohmon(1),incohmon(2))
c01Nov        call pgarro(incohmon(1)+kf(1),incohmon(2)+kf(2),
c01Nov     &                  incohmon(1),incohmon(2))
      else
        call pgsci(color_cyan)   ! cyan
c01Dec        call pgsci(5)   ! cyan
	temp1 = 1.1*circsize
        call pgcirc(incohmon(1),incohmon(2), temp1)
c01Nov        call pgcirc(incohmon(1),incohmon(2),1.1*circsize)
      endif

      if (flagincohana) then
        call pgsci(color_bluecyan)
c02Mar        call pgsci(color_bluemagenta)
c01Dec        call pgsci(6)
	temp1 = 1.1*circsize
        call pgcirc(incohana(1),incohana(2), temp1)
c01Nov        call pgcirc(incohana(1),incohana(2),1.1*circsize)
        call pgarro(ki(1),ki(2),incohana(1),incohana(2))
      else
        call pgsci(color_bluecyan)
c01Dec        call pgsci(5)
	temp1 = 1.1*circsize
        call pgcirc(incohana(1),incohana(2),temp1)
c01Nov        call pgcirc(incohana(1),incohana(2),1.1*circsize)
      endif

c      write (*,*) 
c      write (unitspur,*)

      return
      end

