      subroutine plotlattice(circsize,m_rlv,x,y,h,k,l,qmax,
     &	 gray, rainbow, struct,
     &   sfhkl,absent,replot, overridesize, pointsize,               
     &	 plotunique, uniqueLabel, hide_absences, nscan, qarr)
      
C  hmax is the maximum (h00) value corresponding to maximum Q for
C  the calculation.

C  reflmax is the maximum number of reflections (currently set for
C  a cubic material.
C  nref is the actual number
      
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_color.inc'
      include 'spur_layout.inc'
      include 'spur_scan.inc'
      
c...	common blocks
      logical warning
      common /reciplattice/ warning

      integer natoms
      common /atoms/ natoms

      integer ntiny, nabsent, npresent, nobserved
      common /nabsences/ ntiny, nabsent, npresent, nobserved

c...	functions
      integer select_color
      
c...	constants
c01Dec      integer nmax		! Now natoms_max in arraysize.inc
c01Dec     parameter (nmax=300)

c...	passed variables
      integer m_rlv
      logical struct, gray, rainbow, replot
      real x(m_rlv_max,m_rlv_max),y(m_rlv_max,m_rlv_max)
      integer h(m_rlv_max,m_rlv_max)
      integer k(m_rlv_max,m_rlv_max)
      integer l(m_rlv_max,m_rlv_max)
      real sfhkl(m_rlv_max,m_rlv_max)
      logical absent(m_rlv_max,m_rlv_max)
      real qmax
      logical overridesize
      real pointsize(m_rlv_max,m_rlv_max)

      logical 	plotUnique, hide_absences
      character*1 uniqueLabel(m_rlv_max,m_rlv_max)
      integer nscan
      real*4 qarr(2,maxnpts)
      
c...	local variables
      integer scale
c      real x1,x2,y1,y2
      integer i, j, nint
      integer hklunique(unique_rlv_max,3)
      real sfn(unique_rlv_max)
      real sfmax,sfmin
c      real zmin,zmax
c      integer mult,power,dec
      real circsize,size
      character*6 pgt
c      integer nsub

      integer itemp, this_color
c      integer color_temp
      integer pgt_len
      real	xtemp, ytemp

c      real	X1vp, X2vp, Y1vp, Y2vp
c      real	X1wrld, X2wrld, Y1wrld, Y2wrld
c      real 	xcenwrld, ycenwrld
      
      integer fill_solid, fill_hollow
      parameter (fill_solid = 1)
      parameter (fill_hollow = 2)

cDebug       write (*,*) 'plotlattice: Entry. replot=', replot
       
C      warning=.false.
C  Rescale size of circle to the q=0 scattering cross-section.

      call pgqch(size)
      if (.not.struct) go to 15

      call count_intensities(m_rlv, x, y, absent,  sfhkl, sfn, 
     &		h, k, l, hklunique, sfmin, sfmax, nint)

c	Only needed if will plot circles.
       call scale_circle(scale, sfmin, sfmax, replot, 
     &      m_rlv, pointsize, sfhkl)
       

C find the current limits of the plot window
c      call pgqwin(x1,x2,y1,y2)
c      write (*,*) 'plotlattice: x1,x2,y1,y2=', x1,x2,y1,y2

cDebug      write (*,*) 'plotlattice: xmin,xmax,ymin,ymax=', 
cDebug     &      xmin,xmax,ymin,ymax

      CALL PGSCH(1.0)

 15   continue

ccc      nobserved = 0
      do i=1,m_rlv
       do j=1,m_rlv

c	 !! only plot points inside of the plot's boundaries
         if ((x(i,j).gt.xmin).and.(x(i,j).lt.xmax).and.
     &              (y(i,j).gt.ymin).and.(y(i,j).lt.ymax)) then
ccc           nobserved = nobserved + 1
           if (.not. struct) then
C 		no struct calculated, use large white circles
cDebug              write (*,*) '*** if (.not. struct) then'
              overridesize = .true.
              call pgsci(color_ink)
	      call place_circle(x(i,j), y(i,j), pointsize(i,j), circsize,  
     &      		 fill_hollow, overridesize)

           else		!! if (.not. struct) then

	     if (absent(i,j)) then
cDebug                write (*,*) '*** if (absent(i,j)) then'
		if (.not. hide_absences) then
		  overridesize = .true.
                  call pgsci(color_yellowgreen)
	          call place_cross(x(i,j), y(i,j), 
     &			pointsize(i,j), circsize, overridesize)
		endif
             else
	      if (warning.and.
     &		     (sfhkl(i,j)/float(natoms).lt.tiny)) then
cDebug                write (*,*) '*** if (warning.and....'
                write(*,*) 'plotlattice: Nearly zero intensity marked',
     &	          ' blue at (', h(i,j), ',', k(i,j), ',', l(i,j), ')'
                write(unitlog,*) 'plotlattice: ', 
     &		  'Nearly zero intensity marked',
     &	          ' blue at (', h(i,j), ',', k(i,j), ',', l(i,j), ')'
		overridesize = .true.
                call pgsci(color_blue)
	        call place_circle(x(i,j), y(i,j), pointsize(i,j), circsize,  
     &      		 fill_hollow, overridesize)

	      else 	!! if (warning...) 
               if (gray) then		!! gray circles
cDebug                write (*,*) '*** if (gray) then	'
C		!! sets the pen colour to a grayscale
		overridesize = .false.
c02Jan		color_temp = 18+int(81.1*pointsize(i,j))
		this_color =  select_color(pointsize(i,j))
c		call select_color(pointsize(i,j))
		call place_circle(x(i,j), y(i,j), pointsize(i,j), circsize,  
     &      		 fill_solid, overridesize)
     	       else		!! grey numbers instead
cDebug                  write (*,*) '*** if (.not. gray) then	.. numbers'
     		  itemp = int(100.0*sfhkl(i,j))
		  pgt_len = 6
                  call PGNUMB (itemp, -2, 1, pgt, pgt_len)
		  this_color = select_color(pointsize(i,j))
c		  call select_color(pointsize(i,j))
c02Jan		  color_temp = 18+int(81.1*pointsize(i,j))
c02Jan     		  call pgsci(color_temp)
                  call pgptxt (x(i,j), y(i,j), 0.0, 0.5, pgt)   

	       endif

	      endif 	!! if (warning...) 

	     endif 	!! if (absent(i,j)) then

	     if (plotUnique. and. (gray)) then
	        if ((.not. absent(i,j)) .or. (.not. hide_absences)) then
c...		  Add labels for the intensities
		  xtemp = 1.0*x(i,j) + 2.5 * abs(circsize)
		  ytemp = 1.0*y(i,j) + 2.5 * abs(circsize)
     		  call pgsci(color_ink)
                  call pgptxt (xtemp, ytemp, 0.0, 0.5, 
     &			uniqueLabel(i,j) )
		endif
	     endif

           endif		!! if (.not. struct) then

         endif		!! in range of plot
       enddo
      enddo

      call pgupdt		!! flush the graphics buffer
      
      if (struct.and.warning.and.(.not.replot)) then
        write(*,*)
     &  'WARNING: # of Sf^2 < tiny (', (nabsent + ntiny), 
     &		') and systematic absences (',
     &   	nabsent, ') inconsistent.'
        write(*,*) '  Check smaller intensities marked blue.'
        write(unitlog,*)
     &  'WARNING: # of Sf^2 < tiny (', (nabsent + ntiny), 
     &		') and systematic absences (',
     &   	nabsent, ') inconsistent.'
        write(unitlog,*) '  Check smaller intensities marked blue.'
       endif
      call pgsch(size)
      return
      end



c======================================================
      subroutine count_intensities(m_rlv, x, y,absent,  sfhkl, sfn, 
     &		h, k, l, hklunique, sfmin, sfmax, nint)
c      subroutine count_intensities(xmin, xmax, ymin, ymax, m_rlv, x, y,
c     &       	absent,  sfhkl, sfn, 
cccc     &		nabsent, npresent, ntiny,
c     &		h, k, l, hklunique, sfmin, sfmax, nint)

c...	Some of this is redundent with routines in indexsort.f

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

c...	common blocks
      logical warning
      common /reciplattice/ warning

      integer natoms
      common /atoms/ natoms

      integer ntiny, nabsent, npresent, nobserved
      common /nabsences/ ntiny, nabsent, npresent, nobserved

c...	passed variables
      integer m_rlv
      real x(m_rlv_max,m_rlv_max),y(m_rlv_max,m_rlv_max)
      integer h(m_rlv_max,m_rlv_max)
      integer k(m_rlv_max,m_rlv_max)
      integer l(m_rlv_max,m_rlv_max)
      real sfhkl(m_rlv_max,m_rlv_max)
      logical absent(m_rlv_max,m_rlv_max)
cCommon      real 	xmin,xmax,ymin,ymax
      integer nint
ccc      integer nabsent, npresent, ntiny
      integer hklunique(unique_rlv_max,3)
      real sfn(unique_rlv_max)
      real sfmax,sfmin

c...	local variables
      integer i,j,n, len
      integer hkltmp(3)
      character*10 hkltitle

C  determine unique intensities and sort them
      nint=0
      nabsent = 0
      npresent = 0
      nobserved = 0
      ntiny = 0
      do i=1,m_rlv
        do j=1,m_rlv
	 if ( (xmin .le. x(i,j)) .and. (x(i,j) .le. xmax) .and.
     &	      	   (ymin .le. y(i,j)) .and. (y(i,j) .le. ymax) ) then
          nobserved = nobserved + 1
          if (absent(i,j)) then
C		!! absent reflections are not counted
	    nabsent = nabsent + 1
	    go to 12
	  else
	    npresent = npresent + 1
	  endif

C also if there was a warning then it is likely there are special conditions
C these points will still be plotted but they should not be counted in the
C scaling
          
          if (warning.and.(sfhkl(i,j)/float(natoms).lt. tiny)) then
c01Dec          if (warning.and.(sfhkl(i,j)/float(natoms).lt.1.0e-5)) then
C             write(*,*) 'warning',i,j,sfhkl(i,j)/float(natoms)
	     ntiny = ntiny + 1
             go to 12
          endif

C a WHILE loop
          n=1
  13      if (n.gt.nint) go to 14          ! add a unique reflection
          if (abs(sfhkl(i,j)/sfn(n)-1.0).le.0.001) then  ! not unique
             if (h(i,j)**2+k(i,j)**2+l(i,j)**2.gt.hklunique(n,1)**2+
     &       hklunique(n,2)**2+hklunique(n,3)**2) go to 12       
             hklunique(n,1)=h(i,j)
             hklunique(n,2)=k(i,j)
             hklunique(n,3)=l(i,j)
             go to 12
          endif

          n=n+1                            ! check next unique reflection
          go to 13
  14      nint=nint+1 
  	  if (nint .gt. unique_rlv_max) then
	    write (*,*) 'coutn_intensities: Overflow nint=', nint 
	  endif
  	  if (n .gt. unique_rlv_max) then
	    write (*,*) 'coutn_intensities: Overflow n=', n
	  endif

C          write(*,*) 'unique reflection',nint
          sfn(nint)=sfhkl(i,j)
          if (nint.eq.1) then 
             sfmin=sfn(1)
             sfmax=sfn(1)
          endif
          hklunique(nint,1)=h(i,j)
          hklunique(nint,2)=k(i,j)
          hklunique(nint,3)=l(i,j)
c          list(nint)=nint                  ! used for bubble sorting

	 endif		!! check that things are in the plot range

 12       enddo
       enddo


C      call bubblefloat(nint,sfn,list)

C      write(*,*) 'initial value',sfmin
C put a redundant check in here since this thing is not working very
C well, sfn(1) value is changing ?  CPA 8-02-00
      do n=1,nint
          if ((sfn(n).lt.sfmin).and.(sfn(n)/float(natoms).ge. tiny))
c01Dec          if ((sfn(n).lt.sfmin).and.(sfn(n)/float(natoms).ge.1.0e-5))
     &       then
             sfmin=sfn(n)
          endif
          if (sfn(n).gt.sfmax) sfmax=sfn(n)
      enddo

C  Output the unique intensities and reduced hkl to structinfo.out

      write (unitstruct, *) '>>  count_intensities()'
      write(unitstruct,*) 'There are ',nobserved,' intensities', 
     & 		' in the plot range.'
      write(unitstruct,*) 'There are ',nint,' unique intensities.'
      write(unitstruct,*) 'There are ',npresent,
     &			' inspected intensities.'
      write(unitstruct,*) 'There are ',nabsent,' absent intensities.'
      if ((npresent + nabsent) .ne. nobserved) then
        write(unitstruct,*) 'Error! (npresent + nabsent) .ne. ', 
     &		'nobserved)', npresent, nabsent, nobserved
      endif
      write(unitstruct,*) 'There are ',ntiny,' nearly zero intensities.'
      
      write (*,*) '>>  count_intensities()'
      write(*,*) 'There are ',nobserved,' intensities', 
     & 		' in the plot range.'
      write(*,*) 'There are ',nint,' unique intensities.'
      write(*,*) 'There are ',npresent,' inspected intensities.'
      write(*,*) 'There are ',nabsent,' absent intensities.'
      write(*,*) 'There are ',ntiny,' nearly zero intensities.'
      if ((npresent + nabsent) .ne. nobserved) then
        write(*,*) 'Error! (npresent + nabsent) .ne. ', 
     &		'nobserved)', npresent, nabsent, nobserved
      endif
      
C  Actually if there are any more than 10 unique intensities
C  outputting them really makes no sense.
      if (nint.le.10) then
         do n=1,nint
            do j=1,3
               hkltmp(j)=hklunique(n,j)
            enddo
            call title(hkltmp,hkltitle,len)
            write(unitstruct,'(f8.3,3a)') sfn(n),' at ',hkltitle(1:len),
     &           '-type peaks.'
c            write(7,'(f8.3,3a)') sfn(n),' at ',hkltitle(1:len),
c     &           '-type peaks.'
         enddo
         write(unitstruct,*)
c         write(7,*)
      endif

       return
       end

c======================================================
      subroutine scale_circle(scale, sfmin, sfmax, replot, 
     &      m_rlv, pointsize, sfhkl)

c...	Scale the size of the circle to represent the intensity either
c...	linear or logscale.

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

c...	common blocks
      logical warning
      common /reciplattice/ warning

      integer natoms
      common /atoms/ natoms

c...	passed variables
      real sfhkl(m_rlv_max,m_rlv_max)
      real pointsize(m_rlv_max,m_rlv_max)
      real sfmax,sfmin
      integer scale, m_rlv
      logical  replot

c...	local variables
      real zmin,zmax
      integer mult,power,dec, nsub
      integer  i, j

C  what sort of scaling?

      if (sfmax.eq.sfmin) then
        scale=0
C        write(*,*) sfmax,sfmin
        if (.not.replot) then 
	   write(*,*) 'All intensites identical: no scaling.'
	   write(unitlog,*) 'All intensites identical: no scaling.'
        endif
        do i=1,m_rlv
          do j=1,m_rlv
            pointsize(i,j)=sfhkl(i,j)/sfmax
          enddo
        enddo
        
      else if (sfmax/sfmin.le.100) then
        scale=1
        call linscale(sfmin,sfmax,zmin,zmax,nsub)
        if (.not.replot) then
	   write(*,'(a,f5.1,a,f5.1)') 
     &      'Truncated linear scaling from ',zmin,' to ',zmax
	   write(unitlog,'(a,f5.1,a,f5.1)') 
     &      'Truncated linear scaling from ',zmin,' to ',zmax
        endif
        do i=1,m_rlv
          do j=1,m_rlv
            pointsize(i,j)=sfhkl(i,j)/zmax
          enddo
        enddo

      else
        scale=2
        call logscale(sfmin,sfmax,mult,power,dec)
        zmax=float(mult)*10.0**float(power)
        zmin=float(mult)*10.0**float(power-dec)
        if (.not.replot) then
	  write(*,'(a,g7.2,a,g7.2)') 
     &    	'Truncated log scaling from ',zmin,' to ',zmax
	  write(unitlog,'(a,g7.2,a,g7.2)') 
     &    	'Truncated log scaling from ',zmin,' to ',zmax
        endif
        do i=1,m_rlv
          do j=1,m_rlv
            pointsize(i,j)=log10(sfhkl(i,j)/zmin)/log10(zmax/zmin)
C            if (sfhkl(i,j).gt.zmax) pointsize(i,j)=1.0
C            if (sfhkl(i,j).lt.zmin) pointsize(i,j)=0.0
          enddo
        enddo
      endif

      return
      end


C
C**********************************************************
C      
      subroutine linscale(x1,x2,xmin,xmax,nsub)

      implicit none
c...	functions
      integer myint
      
c...	passed variables
      real x1,x2,xmin,xmax
      integer nsub
      
c...	local variables
      real dif
      
C  how different are x1 and x2?  which power of 10 shows the
C  first difference?
      dif= float(myint(log10(x2-x1)))
C  first attempt at xmax is to round x2 up to the power 
C  of ten which shows up in dif
      xmax=myint(x2/10**dif+.99)*10**dif
C  and xmin is the same power rounded down
      xmin=myint(x1/10**dif)*10**dif
C  what is the best number of subdivisions?
      nsub=int((xmax-xmin)/10**dif+1e-4)

      return
      end

