c...	indexsort.f
c..	From numerical recipes

c=====================================================
      subroutine indexx2d(m_rlv, arrin, indx1, indx2, 
     &	      usedescending, qmax, hh, kk, ll, x, y)
c...	2d extension to indexx from numberical recipes.
c...	Indexes a 2D array arrin(,) of length m_rlv by m_rlv 
c...	i.e. outputs the arrays indx1 and indx2 (each length m_rlv**2)
c...	such that arrin(indx1(j)), indx2(j)) is in ascending order
c...	for J = 1, ... m_rlv**2 .
c...	This uses a heap sort algorthm.

      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_layout.inc'
      
c...	common blocks
      integer ntiny, nabsent, npresent, nobserved
      common /nabsences/ ntiny, nabsent, npresent, nobserved

c...	passed variables      
      integer m_rlv
      real arrin(m_rlv_max,m_rlv_max)
      integer indx1(unique_rlv_max)
      integer indx2(unique_rlv_max)
      logical usedescending
      real 	qmax
      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 	x(m_rlv_max,m_rlv_max),y(m_rlv_max,m_rlv_max)

c...	local variables
      integer	j1, j2
      integer   j, i
      integer	l, ir
      integer	indx1temp, indx2temp
      real	atemp
      logical	debugthis
      parameter (debugthis = .false.)
            
c...	Initialize the index array, only using those 
c...	bragg peaks in the plotting range.
      j = 0
      do j1 = 1, m_rlv
         do j2 = 1, m_rlv
           if ((x(j1,j2).gt.xmin).and.(x(j1,j2).lt.xmax).and.
     &              (y(j1,j2).gt.ymin).and.(y(j1,j2).lt.ymax)) then
ccc              if (.not. absent(j1, j2)) then  
	       j = j + 1
	       indx1(j) = j1
	       indx2(j) = j2
ccc	     endif
           endif
	 enddo
      enddo
c02Jan      nobserved = j
      if (nobserved .eq. 0) then
        nobserved = j
      endif

c...	old way, using the whole array.
c      do j1 = 1, m_rlv
c         do j2 = 1, m_rlv
c	   j = (j1-1)*m_rlv + j2
c	   indx1(j) = j1
c	   indx2(j) = j2
ccDebug	   write(unitdebug, *) 'init', j, j1, j2, 
ccDebug     &	   	arrin(indx1(j), indx2(j))
c	 enddo
c      enddo

c	  WRITE (unitdebug,*) 'pre-indexx2d>>'
c	  do i = 1, m_rlv*m_rlv
c	    write (unitdebug,*) i, arrin(indx1(i), indx2(i))
cc            write (unitdebug,18) '  ', indx1(i), indx2(i), 
cc     &		hh(indx1(i), indx2(i)), kk(indx1(i), indx2(i)), 
cc     &		ll(indx1(i), indx2(i))
cc 18	    format(a, i3, i3, '(', i3, ',', i3, ',', i3,') ') 
c	  enddo
c	  write (unitdebug,*) 

cOld      if (m_rlv .eq. 1) return 
cOld         L = (m_rlv*m_rlv)/2 + 1
cOld      ir = m_rlv*m_rlv


      if (nobserved .eq. 1) return
      
      L = (nobserved)/2 + 1
      ir = nobserved
      
 10   continue
        if (L .gt. 1) then
	  L = L -1
	  indx1temp = indx1(L)
	  indx2temp = indx2(L)
	  atemp = arrin(indx1temp, indx2temp)
	else
	  indx1temp = indx1(ir)
	  indx2temp = indx2(ir)
	  atemp = arrin(indx1temp, indx2temp)
	  indx1(ir) = indx1(1)
	  indx2(ir) = indx2(1)
	  ir = ir -1
	  if (ir .eq. 1) then
	    indx1(1) = indx1temp 
	    indx2(1) = indx2temp 
	    goto 700
cccc	    return
	  endif
	endif

        i = L
        j = L + L

 20     continue
        if (j .le. ir) then
          if (j .lt. ir) then
	    if (arrin(indx1(j), indx2(j)) .lt. 
     &		arrin(indx1(j+1), indx2(j+1)) ) then
     		j = j + 1
     	    endif
	  endif
	  if (atemp .lt. arrin(indx1(j), indx2(j)) ) then
	    indx1(i) = indx1(j)
	    indx2(i) = indx2(j)
	    i = j
	    j = j + j
	  else
	    j = ir + 1
	  endif
	  goto 20
        endif             
        indx1(i) =indx1temp
        indx2(i) =indx2temp
        goto 10
	
 700  continue
      if (usedescending) then	!! swap from ascending order to decending
        do j = 1, (nobserved)/2
cOld        do j = 1, (m_rlv*m_rlv)/2
	  indx1temp = indx1(j)
	  indx1(j) =  indx1(nobserved - j + 1)
cOld	  indx1(j) =  indx1(m_rlv*m_rlv - j + 1)
	  indx1(nobserved - j + 1) = indx1temp
cOld	  indx1(m_rlv*m_rlv - j + 1) = indx1temp
	  indx2temp = indx2(j)
	  indx2(j) =  indx2(nobserved - j + 1)
cOld	  indx2(j) =  indx2(m_rlv*m_rlv - j + 1)
	  indx2(nobserved - j + 1) = indx2temp
cOld	  indx2(m_rlv*m_rlv - j + 1) = indx2temp
	enddo
      endif

	  if (debugthis) then
	    WRITE (unitdebug,*) 'post-indexx2d>>'
	    WRITE (unitdebug,*) ' usedescending =', useDescending
	    do i = 1, (nobserved)
cOld	    do i = 1, (m_rlv*m_rlv)
              write (unitdebug,17) arrin(indx1(i), indx2(i)), i, 
     &	    	indx1(i), indx2(i), 
     &		hh(indx1(i), indx2(i)), 
     &		kk(indx1(i), indx2(i)), 
     &		ll(indx1(i), indx2(i))
 17	      format('sfhkl=', g12.5, 1x, 'i=', i4, 
     &		'. indx1=', i3, ' indx2=', i3, ' (', i3, ',', i3, ',', i3,') ') 
	    enddo
	    write (unitdebug,*) 
	  endif

      return
      end
      
c===============================================================================
      subroutine unique(m_rlv, arrin, absent,  indx1, indx2, 
     &	  	indxunique1, indxunique2, multiplicity, 
     &		limlo, limhi, nunique, uniqueLabel, uniqueRank,
     &		sorttol, base, uselogscale, usedecending, hh, kk, ll)
c02Mar     &		tolerance, uselogscale, usedecending, hh, kk, ll)
c...	Return indxunique1 and  indxunique2 as an index to
c...    the first occurance of each value in arrin   
c...	nunique ss the length of indxunique1 and  indxunique2  that is used.
c...	mulitplicity(j) contains the multiplicity of the value 
c...	arrin(indxunique1(j), indxunique2(j))
c...    If logscale is false, then all entries within tolerance will
c...	be lumped together.
c...	If logscale is true, then use int(tolerance) as the base for 
c...	the logscale. All entries within the range of base**n and base**(n+1)
c...	will be lumped together,

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

c...	common blocks
      integer natoms
      common /atoms/ natoms
      
      integer ntiny, nabsent, npresent, nobserved
      common /nabsences/ ntiny, nabsent, npresent, nobserved

c...	passed variables
      integer m_rlv
      real arrin(m_rlv_max,m_rlv_max)
      logical absent(m_rlv_max,m_rlv_max)
cCommon      integer nabsent
      integer indx1(unique_rlv_max), indx2(unique_rlv_max)
      integer indxunique1(unique_rlv_max)
      integer indxunique2(unique_rlv_max)
      integer multiplicity(unique_rlv_max)
      real	limlo(unique_rlv_max)
      real	limhi(unique_rlv_max)
c      integer indxunique1(0:unique_rlv_max)
c      integer indxunique2(0:unique_rlv_max)
c      integer multiplicity(0:unique_rlv_max)
c      real	limlo(0:unique_rlv_max)
c      real	limhi(0:unique_rlv_max)
      integer nunique
      character*1 uniqueLabel(m_rlv_max,m_rlv_max)
      integer	uniqueRank(m_rlv_max,m_rlv_max)
c02Mar      real 	tolerance
      real 	sorttol, base
      logical 	uselogscale
      logical	useDecending
      integer hh(m_rlv_max,m_rlv_max)
      integer kk(m_rlv_max,m_rlv_max)
      integer ll(m_rlv_max,m_rlv_max)
    
c...	local variables
      integer 	j, jq
      integer   i
c      integer j1, j2
      real 	value
c02Mar      integer 	base
cccc      real	limlo, limhi
      real 	pwrhi, pwrlo      
      integer	thisrank
      integer	absentRank
      parameter (absentRank = -1)
      logical 	debugthis, debugthis2
c      parameter (debugthis = .false.)
      parameter (debugthis = .true.)
      parameter (debugthis2 = .true.)
      logical inrange
      integer istart, nmult
           
c02Mar      if (uselogscale) then
c02Mar	base = nint(abs(tolerance))
c02Mar      endif

c      debugthis = .true.
      if (debugthis) then
	write (unitdebug,*) 'unique>>'
	write (unitdebug,*) 'uselogscale=',  uselogscale
	write (unitdebug,*) 'usedecending=', usedecending
	write (unitdebug,*) 'tolerance =', sorttol
	write (unitdebug,*) 'base=', base
	write (unitdebug,*) 'nobserved=', nobserved
      endif
	write (*,*) 'nobserved=', nobserved
      
      j = 1
      nabsent = 0

c...	Last item is always the absences 
      limlo(unique_rlv_max) = 0.0
      limhi(unique_rlv_max) = 0.0
      multiplicity(unique_rlv_max) = 0
      indxunique1(unique_rlv_max) = indx1(j)	!! Not correct if no absences
      indxunique2(unique_rlv_max) = indx2(j)
      uniqueRank(indx1(j), indx2(j))  = absentRank
      inrange  = .false.
      
      do while (absent(indx1(j), indx2(j)) .and. 
     &	      (j .le. min(nobserved, unique_rlv_max) ) )
cOld      do while (absent(indx1(j), indx2(j)) .and. 
cOld     &	      (j .le. min(m_rlv*m_rlv, unique_rlv_max) ) )
        nabsent = nabsent + 1
	multiplicity(unique_rlv_max) = multiplicity(unique_rlv_max) + 1
	if (nabsent .eq. 1) then
	  !! Stores the first absence found
          indxunique1(unique_rlv_max) = indx1(j) 
          indxunique2(unique_rlv_max) = indx2(j)
	endif
        uniqueRank(indx1(j), indx2(j))  = absentRank
        if (debugthis) then
	   write (unitdebug,*) 'absences at begining of list:'
	   write (unitdebug,*)
     &      j, indx1(j), indx2(j), '...', 
     &	    unique_rlv_max, 
     &		indxunique1(unique_rlv_max), indxunique2(unique_rlv_max), '...', 
     &	    ' mult=', multiplicity(unique_rlv_max),
     &      ' rank=',   uniqueRank(indx1(j), indx2(j)),
     &      ' nabsnt = ', nabsent, absent(indx1(j), indx2(j)) 
	   write (unitdebug,*) ' '
        endif
	j = j + 1

      enddo
      if (j .gt. min(nobserved, unique_rlv_max)) goto 666
cOld      if (j .gt. min(m_rlv*m_rlv, unique_rlv_max)) goto 666
      
      jq = j
      uniqueRank(indx1(j), indx2(j))  = jq
      indxunique1(jq) = indx1(j)
      indxunique2(jq) = indx2(j)
      multiplicity(jq) = 0
      value = arrin(indx1(j), indx2(j))

       if (uselogscale) then
	  pwrlo = int(log(value)/log(base) )
c02Mar	  pwrlo = int(log(value)/log(float(base)) )
	  if (value .le. 1.0) pwrlo = pwrlo - 1
	  pwrhi = pwrlo + 1
	  limlo(jq) = base**pwrlo
	  limhi(jq) = base**pwrhi
c02Mar	  limlo(jq) = float(base)**pwrlo
c02Mar	  limhi(jq) = float(base)**pwrhi
       else
           limlo(jq) = value - sorttol
           limhi(jq) = value + sorttol
c02Mar           limlo(jq) = value - tolerance
c02Mar           limhi(jq) = value + tolerance
       endif
       inrange = ( (limlo(jq) .le. arrin(indx1(j), indx2(j))) 
     &		.and. (arrin(indx1(j), indx2(j)) .le. limhi(jq)) )

       if (debugthis) then
	    write (unitdebug,*) '... rank=', jq, 
     &	    ' limits=', limlo(jq), limhi(jq)
       endif

      do while (j .le. min(nobserved, unique_rlv_max) )
cOld      do while (j .le. min(m_rlv*m_rlv, unique_rlv_max) )

  	  do while ((inrange) .and. 
     &		(j .le. min(nobserved, unique_rlv_max) ) )	  
cOld  	  do while ((inrange) .and. 
cOld     &		(j .le. min(m_rlv*m_rlv, unique_rlv_max) ) )	  
c  	  do while ( (limlo(jq) .le. arrin(indx1(j), indx2(j))) 
c     &		.and. (arrin(indx1(j), indx2(j)) .le. limhi(jq)) )
     	    multiplicity(jq) = multiplicity(jq) + 1
            if (debugthis) then
	      write (unitdebug,*)
     &      	j, ' val=', arrin(indx1(j), indx2(j)),
     &          indx1(j), indx2(j), '...', 
     &	    	jq, indxunique1(jq), indxunique2(jq), '...', 
     &	    	' mult=', multiplicity(jq),
     &	    	' rank=',  uniqueRank(indx1(j), indx2(j)),
     &          ' nabsnt = ', nabsent, absent(indx1(j), indx2(j))    
            endif

	    j = j + 1
	    if (j .le. min(nobserved, unique_rlv_max) ) then
cOld	    if (j .le. min(m_rlv*m_rlv, unique_rlv_max) ) then
              uniqueRank(indx1(j), indx2(j))  = jq
              inrange = ( (limlo(jq) .le. arrin(indx1(j), indx2(j))) 
     &		.and.   (arrin(indx1(j), indx2(j)) .le. limhi(jq)) )
	    endif
	  enddo
	
	if (j .le. min(nobserved, unique_rlv_max)) then
cOld	if (j .le. min(m_rlv*m_rlv, unique_rlv_max)) then
	
         if  (absent(indx1(j), indx2(j))) then 
	   if (nabsent .eq. 0) then
	     !! Stores the first absence found
             indxunique1(unique_rlv_max) = indx1(j)
             indxunique2(unique_rlv_max) = indx2(j)
	   endif
           nabsent = nabsent + 1
	   multiplicity(unique_rlv_max) = multiplicity(unique_rlv_max) + 1
	   uniqueRank(indx1(j), indx2(j)) = absentRank
           if (debugthis) then
	     if (nabsent .eq. 1) then
	       write (unitdebug,*) ' '
	       write (unitdebug,*) 'absences at end of list:'
	     endif
	     write (unitdebug,*)
     &      	j, ' val=', arrin(indx1(j), indx2(j)),
     &           indx1(j), indx2(j), '...', 
     &	    	unique_rlv_max, indxunique1(unique_rlv_max),
     &		 indxunique2(unique_rlv_max), '...', 
     &	    	' mult=', multiplicity(unique_rlv_max),
     &	    	' rank=',  uniqueRank(indx1(j), indx2(j)),
     &      	' nabsnt = ', nabsent, absent(indx1(j), indx2(j))
           endif

	   j = j + 1

         else

	  jq = jq + 1
	  uniqueRank(indx1(j), indx2(j))  = jq
          indxunique1(jq) = indx1(j)
          indxunique2(jq) = indx2(j)
          multiplicity(jq) = 0
          value = arrin(indx1(j), indx2(j))
cDebug	  write (*,*) j, indx1(j), indx2(j), jq, 
cDebug     &	     indxunique1(jq), indxunique2(jq),
cDebug     &	    multiplicity(jq), uniqueRank(indx1(j), indx2(j))  
          if (uselogscale) then
	    pwrlo = int(log(value)/log(base) )
c02Mar	    pwrlo = int(log(value)/log(float(base)) )
	    if (value .le. 1.0) pwrlo = pwrlo - 1
	    pwrhi = pwrlo + 1
	    limlo(jq) = base**pwrlo
	    limhi(jq) = base**pwrhi
c02Mar	    limlo(jq) = float(base)**pwrlo
c02Mar	    limhi(jq) = float(base)**pwrhi
           else
             limlo(jq) = value - sorttol
             limhi(jq) = value + sorttol
c02Mar             limlo(jq) = value - tolerance
c02Mar             limhi(jq) = value + tolerance
          endif
          if (debugthis) then
	    write (unitdebug,*) '... rank=', jq, 
     &	    ' limits=', limlo(jq), limhi(jq)
          endif
          inrange = ( (limlo(jq) .le. arrin(indx1(j), indx2(j))) 
     &		.and. (arrin(indx1(j), indx2(j)) .le. limhi(jq)) )
	  
	 endif
	endif
	if (jq .ge. min(nobserved, unique_rlv_max)) goto 666
cOld	if (jq .ge. min(m_rlv*m_rlv, unique_rlv_max)) goto 666

      enddo

      nunique = jq
      write (*,*) 'unique: nunique, nobserved =', nunique, nobserved

c      	!! Label for intensity range.
c	!! '*' for absent.
c	!! a..z for lowest intensities (a being the very lowest).
c	!! A...Z for the highest intensities (A being the highest intensity).
c	!! # for the middle intensities when there are more than 54. 
c	!! - for did something stupid

c...	Visual verification that I did not do something stupid.
      do j = 1 + min(nobserved, unique_rlv_max), unique_rlv_max
            uniqueLabel(indx1(j), indx2(j)) = '-'
      enddo


      do j = 1, min(nobserved, unique_rlv_max)
cOld      do j = 1, min(m_rlv*m_rlv, unique_rlv_max)
          thisrank = uniqueRank(indx1(j), indx2(j))
          if (thisrank .eq. absentRank) then
            uniqueLabel(indx1(j), indx2(j)) = '*'
ccc            uniqueLabel(indx1(j), indx2(j)) = '@'
          else if (thisrank .le. min(26, nunique/2)) then
	    if (useDecending) then
	      !! rank=1 is the largest number with useDecending = .true.
              uniqueLabel(indx1(j), indx2(j)) = 
     &	      	char(ichar('A') -1 + thisrank)
	    else
	      !! rank=1 is the smallest number with useDecending = .false.
              uniqueLabel(indx1(j), indx2(j)) = 
     &	      	char(ichar('a') -1 + thisrank)
	    endif
          else if ((nunique - thisrank + 1) .le.
     &		   min(26, (nunique-nunique/2))) then
	    if (useDecending) then
	      !! rank=nquniue is the smallest number with useDecending = .true.
              uniqueLabel(indx1(j), indx2(j)) = char(ichar('a') -1 + 
     &	      		(nunique - thisrank + 1))
	    else
	      !! rank=nunique is the largest number with useDecending = .false.
              uniqueLabel(indx1(j), indx2(j)) = char(ichar('A') -1 + 
     &	      		(nunique - thisrank + 1))
	    endif
          else
            uniqueLabel(indx1(j), indx2(j)) = '#'
          endif
c	  write (*,*) 'unique: rank ', j, indx1(j), indx2(j),
c     &		   arrin(indx1(j), indx2(j)), ' "', 
c     &		  uniqueLabel(indx1(j), indx2(j)), '" ', thisrank 
      enddo

c      do j1 = 1, m_rlv
c       do j2 = 1, m_rlv
c          thisrank = uniqueRank(j1, j2)
c          if (thisrank .le. min(26, nunique/2)) then
c	    if (useDecending) then
c	      !! rank=1 is the largest number with useDecending = .true.
c              uniqueLabel(j1, j2) = char(ichar('A') -1 + thisrank)
c	    else
c	      !! rank=1 is the smallest number with useDecending = .false.
c              uniqueLabel(j1, j2) = char(ichar('a') -1 + thisrank)
c	    endif
c          else if ((nunique - thisrank + 1) .le.
c     &		   min(26, (nunique-nunique/2))) then
c	    if (useDecending) then
c	      !! rank=nquniue is the smallest number with useDecending = .true.
c              uniqueLabel(j1, j2) = char(ichar('a') -1 + 
c     &	      		(nunique - thisrank + 1))
c	    else
c	      !! rank=nunique is the largest number with useDecending = .false.
c              uniqueLabel(j1, j2) = char(ichar('A') -1 + 
c     &	      		(nunique - thisrank + 1))
c	    endif
c          else
c            uniqueLabel(j1, j2) = '#'
c          endif
c	  write (*,*) 'unique: rank ', j1, j2, arrin(j1,j2), '"', 
c     &		  uniqueLabel(j1, j2), '"', thisrank 
c        enddo
c      enddo

	  if (debugthis2) then
	    WRITE (unitdebug,*) 'post-unique>>'
	    WRITE (unitdebug,*) ' tolerance =', sorttol
	    WRITE (unitdebug,*) ' base =', base
	    WRITE (unitdebug,*) ' uselogscale = ', uselogscale
     	    WRITE (unitdebug,*) ' nunique =', nunique
	    nmult = 0
	    do i = 1, nunique
	      nmult = nmult + multiplicity(i)
	      write (unitdebug,*) multiplicity(i), ' in range ', 
     &	      		limlo(i), limhi(i)
	    enddo
	    WRITE (unitdebug,*) 'total  # of non-absent intensisties=', 
     &		    nmult
     	    WRITE (unitdebug,*) ' nabsent =', nabsent
	    WRITE (unitdebug,*) 'total  # of non-absent and ', 
     &		    'absent intensisties=', nmult+nabsent
            WRITE (unitdebug,*)
	    
	    if (nabsent .gt. 0) then 
	      istart = 0
	    else
	      istart = 1
	    endif

	    do i = istart, nunique
cNothere	     if (debugthis) then
cDebug	    write (*,*) 'unique: unformatted write to unitdebug 1'
	      write (unitdebug,*) 
     &		uniqueLabel(indxunique1(i), indxunique2(i)), i, 
     &		arrin(indxunique1(i), indxunique2(i)),
     &		' limlo =', limlo(i), ' limhi=', limhi(i), 
     &		' mult=', multiplicity(i)
cNothere	     endif
c		maybe next line causes SGI to crash.      
c	inwrite: off end of record
c	apparent state: unit 20 named SPUR_DEFAULT.debug
c	last format: (a)
c	Unit 20 is a sequential formatted internal file
cDebug	    write (*,*) 'unique: formatted write to unitdebug 2'
              write (unitdebug,18) '              ', 
     &	    	indxunique1(i), indxunique2(i), 
     &		hh(indxunique1(i), indxunique2(i)), 
     &		kk(indxunique1(i), indxunique2(i)), 
     &		ll(indxunique1(i), indxunique2(i))
 18	      format(a, i3, i3, ' (', i3, ',', i3, ',', i3,') ') 
	    enddo
	    write (unitdebug,*) 
cDebug	    write (*,*) 'unique: finished write to unitdebug 3'
	  endif
	  

      goto 999
      
 666  write (*,*) 'unique: Bug in algorthm. j, jq = ', j, jq
      stop
           
 999  return
      end
