c....	spurion_routines.f

c===========================================================================
      subroutine init_vectors(a, aa, bmat, bmatuv, b, 
     &      uu, vv, uxv, uuv, vuv, u, v, uxvi, 
     &	    zetatitle, etatitle, c, angle, xstar, ystar)


      implicit none
      include 'io_files.inc'

c...	common blocks

c...	functions
      
c...	passed variables
      REAL*4 a(3),aa(3),bmat(3,3),bmatuv(3,3),b(3)
      REAL*4 uu(3),vv(3),uxv(3),uuv(3),vuv(3)
      INTEGER u(3),v(3),uxvi(3)
      character*10 zetatitle,etatitle
      real c
      real angle
C  The underlying plot and all calculations are done in AA^-1
C  the actual numbers on the graph may be changed to rlu so
C  that xworld= (x AA^-1)/xstar.  The plot will still be
C  square. 
      real xstar, ystar
      
c...	local variables
      INTEGER i
      CHARACTER*5 scale
      integer len
      character*80 aline

C********************************************************      
C  Generate and plot points in reciprocal space.

C  Convert u and v to floating point
      
 50   do i=1,3
        uu(i)=float(u(i))
        vv(i)=float(v(i))
      enddo

C  Calculate the reciprocal lattice vectors in
C  terms of Cartesian coordinates.

      call generateb(a,aa,b,c,bmat)
C      write(*,'(a31)') 'Reciprocal lattice vectors are:'
C      call displayvec(bmat)

C  Calculate the cross product (for main title)

C the answer uxv is in terms of the direct lattice vectors

      call cross(uu,vv,uxv)

      do i=1,3
         uxvi(i)=int(uxv(i))
      enddo

C put in an error trap

      if ((uxvi(1).eq.0).and.(uxvi(2).eq.0).and.(uxvi(3).eq.0)) then
         write(*,'(2a)') 'Orientation vectors are unsuitable!',
     &              ' (zero cross product)'
         call title(u,zetatitle,len)         
         write(*,'(4a)') 'Currently u=',zetatitle(1:len)
         call title(v,etatitle,len)
         write(*,'(4a)') 'Currently v=', etatitle(1:len)
c         write(*,'(4a)') 'Currently u=',zetatitle(1:len),' and v=',
c     &                    etatitle(1:len)
220      write(*,'(a,$)') 'Continue (y)/n? '
         read (*,'(a)',err=220) scale
         if (scale.eq.'n') stop
222      write(*,'(1x,a,$)') 'Enter new u (use spaces or commas): '
         read(*,*, err=222) (u(i),i=1,3)
224      write(*,'(1x,a,$)') 'Enter new v (use spaces or commas): '
         read(*,*, err=224) (v(i),i=1,3)
         go to 50 
      endif
C  Calculate the direct lattice vectors which really isn't
C  necessary.  Include amat in declarations if used.
C      call generatea(amat,bmat)
C      write(*,'(a27)') 'Direct lattice vectors are:'
C      call displayvec(amat)

C  Perform the rotational transformation to a coordinate
C  system where u // +x, u x v // +/- z, and v.y > 0.
      call transform(u,v,bmat,uuv,vuv,bmatuv)
C      write(*,'(a31)') 'Transformed reciprocal lattice vectors are:'
C      call displayvec(bmatuv)

C  What is the angle between v and eta which will appear in ytitle?

      angle=atan(vuv(1)/vuv(2))

      xstar=uuv(1)
      ystar=vuv(2)
 80   format(a80)
300   format (a,f5.3,a,f5.3,a)
      write(aline,300) 'Zeta* is ',uuv(1),' and eta* is ',
     &         vuv(2),' AA^-1.'
      write (*,80) aline
      write (unitlog,80) aline

310   format (a,a,f5.3,a,f5.3,a,  a, f5.3,a,f5.3,a)
      write (aline, 310)   'Orientation vector u is ', 
     &      '(', uuv(1),',', uuv(2), ') AA^-1  = ', 
     &      '(', uuv(1)/xstar,',', uuv(2)/ystar, ') rlu'  
      write (*,80) aline
      write (unitlog,80) aline

      write (aline, 310)   'Orientation vector v is ', 
     &      '(', vuv(1),',', vuv(2), ') AA^-1  = ', 
     &      '(', vuv(1)/xstar,',', vuv(2)/ystar, ') rlu'  
      write (*,80) aline
      write (unitlog,80) aline

      return
      end

c===========================================================================
      subroutine init_vectorlabels(rlu, u, v, zetatitle, etatitle, 
     &	xtitle, ytitle, x2title, y2title, angle, pi)

      implicit none

c...	common blocks

c...	functions

c...	passed variables
      logical rlu
      INTEGER u(3),v(3)
      character*10 zetatitle,etatitle
      CHARACTER*41 xtitle,ytitle
      CHARACTER*41 x2title,y2title
      real angle
      real pi

c...	local variables
      integer len

      
C  Produce axis labels with (hkl) and (h'k'l') of orientation vectors.
      
 60   if (rlu) then
         call title(u,zetatitle,len)
         write(xtitle,'(a22,a)') '\gz [rlu] || to ',zetatitle
         write(x2title,'(a22,a)') '\gz [\A\u-1\d] || to ',zetatitle
c         write(xtitle,'(a22,a)') '\gz (rlu) || to ',zetatitle
c         write(x2title,'(a22,a)') '\gz (\A\u-1\d) || to ',zetatitle

         call title(v,etatitle,len)
         write(ytitle,'(a15,i3,a10,a)') '\gy [rlu] ',
     &        int(180.0*angle/pi),' deg from ',etatitle
         write(y2title,'(a15,i3,a10,a)') '\gy [\A\u-1\d] ',
     &        int(180.0*angle/pi),' deg from ',etatitle
c         write(ytitle,'(a15,i3,a10,a)') '\gy (rlu) ',
c     &        int(180.0*angle/pi),' deg from ',etatitle
c         write(y2title,'(a15,i3,a10,a)') '\gy (\A\u-1\d) ',
c     &        int(180.0*angle/pi),' deg from ',etatitle
      else
         call title(u,zetatitle,len)
         write(xtitle,'(a22,a)') '\gz [\A\u-1\d] || to ',zetatitle
         write(x2title,'(a22,a)') '\gz [rlu] || to ',zetatitle
c         write(xtitle,'(a22,a)') '\gz (\A\u-1\d) || to ',zetatitle
c         write(x2title,'(a22,a)') '\gz (rlu) || to ',zetatitle
         call title(v,etatitle,len)
         write(ytitle,'(a15,i3,a10,a)') '\gy [\A\u-1\d] ',
     &        int(180.0*angle/pi),' deg from ',etatitle
         write(y2title,'(a15,i3,a10,a)') '\gy [rlu] ',
     &        int(180.0*angle/pi),' deg from ',etatitle
c          write(ytitle,'(a15,i3,a10,a)') '\gy (\A\u-1\d) ',
c     &        int(180.0*angle/pi),' deg from ',etatitle
c         write(y2title,'(a15,i3,a10,a)') '\gy (rlu) ',
c     &        int(180.0*angle/pi),' deg from ',etatitle
      endif

      return
      end

c===========================================================================
      subroutine init_scanvectors(constq, ec, es, npts, angle, 
     &      en, ei, ef, fe, h, hc, hs, nscan,
     &      mon, bmatuv, q, ki, kf, a4, a3, phimax, fail, pi, 
     &      qarr, kiarr, kfarr, qspur, 
     &      espur, qsparr, esparr, incohmon, incohana,
     &      incmonarr, incanaarr, prev)
cCommon     &      incmonarr, incanaarr, prev, xmin, xmax, ymin, ymax)
     
      implicit none
cDec01      IMPLICIT REAL (A-H,O-Z)
      include 'arraysize.inc'
      include 'spur_scan.inc'
      include 'io_files.inc'

c...	common blocks

c...	functions
      
c...	passed variables
      logical constq
      real ec, es
      integer npts
      real angle, en, ei, ef, fe
      REAL*4 h(3),hc(3),hs(3)
      integer nscan
      CHARACTER*5 mon
      REAL*4 bmatuv(3,3)
      real*4 q(3),ki(2),kf(2)
      real a4, a3
      real phimax
      logical fail
      real pi
      real*4 qarr(2,maxnpts),kiarr(2,maxnpts),kfarr(2,maxnpts)
      real qspur(2,n_spur_max,n_spur_max)
      real espur(n_spur_max,n_spur_max)
      real qsparr(2,n_spur_max,n_spur_max, maxnpts)
      real esparr(n_spur_max,n_spur_max,maxnpts)
      real*4 incohmon(2),incohana(2)
      real*4 incmonarr(2,maxnpts),incanaarr(2,maxnpts)
cCommon      real xmin, xmax, ymin, ymax
      logical prev
      
c...	local variables
      integer k, i, j, l
      real a4_1
      
C  Begin a loop for the points in the scan.

 63   do k=1,npts

C  Calculation of |ki| and |-kf| at (E,Q)
C     
C  value of en for point k         
         en=ec+(k-(float(npts)+1.0)/2.0)*es
C  value of q for point k         
         do i=1,3
            h(i)=hc(i)+(k-(float(npts)+1.0)/2.0)*hs(i)
            if (hs(i).ne.0) constq=.false.
         enddo

 65      if (mon.eq.'(MON)') then
            ei=fe
            ef=ei-en
         else
            ef=fe
            ei=ef+en
         endif

C  Calculation and text display of ki and Q and the scattering angles.

      if ((ei.le.0).or.(ef.le.0)) then
         write(*,'(2a,i2)') 'Non-positive energy at monochromator or ',
     &     'analyser.  Point# ',k
         write(*,'(2a)') 'This point and remaining points will not be',
     &     ' included.'
         write(unitlog,'(2a,i2)') 'Non-positive energy at ',
     &     'monochromator or analyser.  Point# ',k
         write(unitlog,'(2a)') 'This point and remaining points ',
     &     'will not be included.'
         nscan=k-1
         go to 510
      endif

c...	Note: calcki sets "a4_ok(pn) = .false." when scattering triangle
c...	 can not be closed. True otherwise.
      call calcki(bmatuv,h,ei,ef,q,ki,fail,a4,a3,k)

        angle4_store(k) = a4
cDebug	write (*,*) 'init_scanvectors: angle4=',k, 
cDebug     &		angle4_store(k) *180.0 /pi
	if (k .eq. 1) then
	  a4_1 = a4
	  angle4_varies = .false.
	else
	  if (a4_1 .ne. a4) angle4_varies = .true.
	endif

      if (a4.gt.phimax) then
        a4_ok(k) = .false.
 120    format (a,f6.2,a,i2,a,f6.2,a)
         write(*,120) 'A4=',180.0*a4/pi,' at point# '
     &           ,k,' exceeds limit=', (180.0*phimax/pi),'.'
         write(unitlog,120) 'A4=',180.0*a4/pi,' at point# '
     &           ,k,' exceeds limit=', (180.0*phimax/pi),'.'
c02Apr         write(*,'(a,f6.2,a,i2,a,f6.2,a)') 'A4=',180*a4/pi,' at point# '
c02Apr     &           ,k,' exceeds limit=',180*phimax/pi,'.'
c	A4_OK(k) = .false.
c      else
c        A4_OK(k) = .true.
      endif

C      write(*,'(f7.2,2xf7.2)') 180*a3/pi,180*a4/pi

C  Below is a routine that I added to calculate the Q value one would
C  obtain for a TOF type of scan (A3 and A4 fixed) as E changes 
C      phi=pi/180*88.951   ! the value of A4
C      ef=ef+6             ! the value of -E
C      call calckitof(bmat,ei,ef,h,q,phi,ki,u,v,uuv,vuv)
C  part of BT2 bkbo experiment, Nov. 1998

C  use vector arithmetic to calculate kf
      kf(1)=ki(1)-q(1)
      kf(2)=ki(2)-q(2)

C  store the scattering vectors in several arrays
      qarr(1,k)=q(1)
      qarr(2,k)=q(2)
      kiarr(1,k)=ki(1)
      kiarr(2,k)=ki(2)
      kfarr(1,k)=kf(1)
      kfarr(2,k)=kf(2)


C  Spurious scattering [first spurious peak qspur(*,1,1)=q(*)]
C  n is the number of multiple k's to calculate
      call calcspurious(ki,kf, qspur,espur)
cDec01      call calcspurious(ki,kf,n,qspur,espur)
      do i=1,n_spur_max
         do j=1,n_spur_max
            do l=1,2
               qsparr(l,i,j,k)=qspur(l,i,j)
            enddo
            esparr(i,j,k)=espur(i,j)
         enddo
      enddo

C  Spurious scattering from |ki|=|kf| i.e., scattering is
C  incoherent from the monochromator or analyser and
C  Bragg scattering from the other.  Effect will be large
C  when it matches a Bragg peak in the sample.

      call calcincoh(ki,kf,incohmon,incohana)

      do l=1,2
         incmonarr(l,k)=incohmon(l)
         incanaarr(l,k)=incohana(l)
      enddo

      if (.not. prev) then
c...	Already done at the main level. Don't override current values
C  Determine a reasonable scale for the plot using autoscale
C  routine.

C  The principal vectors which must be plotted are the Q and
C  ki and also u and v.

C      call autoscale(n*ki(1),n*ki(2),xmin,xmax,ymin,ymax)
C      call autoscale(uuv(1),uuv(2),xmin,xmax,ymin,ymax)
C      call autoscale(vuv(1),vuv(2),xmin,xmax,ymin,ymax)
      do i=1,n_spur_max
        do j=1,n_spur_max
          call autoscale(qspur(1,i,j),qspur(2,i,j))
cCommon          call autoscale(qspur(1,i,j),qspur(2,i,j),
cCommon     &                   xmin,xmax,ymin,ymax)
        enddo
      enddo

C  Prevents the plot from being excessively lopsided.
      call square()
cCommon      call square(xmin,xmax,ymin,ymax)

       endif
      enddo

      nscan=npts
 510  continue

      return
      end

c===========================================================================
      subroutine prep_elastic_spur(npts, ec, es,  fe, hc, hs, mon)
     
      implicit none
cDec01      IMPLICIT REAL (A-H,O-Z)
      include 'arraysize.inc'
      include 'spur_scan.inc'
      include 'io_files.inc'

c...	common blocks

c...	functions
      
c...	passed variables
      integer npts
      real ec, es
      real fe
      REAL*4 hc(3), hs(3)
      CHARACTER*5 mon
      
c...	local variables
      integer ii, jj, j
      real ei, ef
      real k_real, e_real
      real e_check
      real h(3)
      
C  value of en for point k         
c         en=ec+(k-(float(npts)+1.0)/2.0)*es

      n_zeroelst = 0
c23456789_123456789_123456789_123456789_123456789_123456789_123456789_12
cDebug      write (*,*) 'prep_elastic_spur: Entered'
      write (*,*) '&&& Elastic spurions within scan range at:'
      write (unitspur,*) '&&& Elastic spurions within scan range at:'
      do ii=1,n_spur_max
        do jj=1,n_spur_max
	 if (.not. (ii .eq. jj)) then
c	Elastic spurion condition: ii*ei - jj*ef = 0	
c		From calcspurious()    
          if (mon.eq.'(MON)') then
            ei=fe
	    k_real =(float(npts)+1.0)/2.0 
     &	    	+(1.0-(float(ii)/float(jj))**2)*(ei/es) -ec/es
	    e_real = ec+(k_real-(float(npts)+1.0)/2.0)*es
	    ef = ei - e_real
          else
            ef=fe
	    k_real =(float(npts)+1.0)/2.0 
     &	     	-(1.0-(float(jj)/float(ii))**2)*(ef/es) -ec/es
	    e_real = ec+(k_real-(float(npts)+1.0)/2.0)*es
	    ei = ef + e_real
          endif
	  e_check = (float(ii)**2) * ei - (float(jj)**2)*ef

            do j=1,3
               h(j)=hc(j)+(k_real-(float(npts)+1.0)/2.0)*hs(j)
            enddo

cDebug	  write (*,*) 'prep_elastistic_spur: ', ii, '*ki -', jj, '*kf',
cDebug     &	    ' elast. at k=', k_real
cDebug          write (*,*) 'prep_elastistic_spur: ', 
cDebug     &	  	ii, '^2 *ei - ', jj, '^2 * ef = 0 = ', e_check
cDebug          write (*,*) 'prep_elastistic_spur: ', 
cDebug     &	  	'ei = ', ei,  ' ef = ', ef

	  if ((k_real .ge. 1.0) .and. (k_real .le. float(npts)) ) then
	    write (*,110) ii, jj,  k_real, 
     &		    h(1), h(2), h(3), e_real 
	    write (unitspur,110) ii, jj,  k_real, 
     &		    h(1), h(2), h(3), e_real 
110 	    format (i1,'*ki - ',  i1, '*kf at scan pt= ', f7.2,
     &		' Q=(',  f6.2, ',', f6.2, ',', f6.2, ') A^-1 dE=', f7.3) 
cDebug	    write (*,*) 'Elastic spurion ', ii, '*ki -', jj, '*kf',
cDebug     &	    ' at scan point ', k_real, 
cDebug     &	    ' dE=', e_real	!! ec+(k_real-(float(npts)+1.0)/2.0)*es
cDebug          write (*,*) 'prep_elastistic_spur: ', 
cDebug     &	  	ii, '^2 *ei - ', jj, '^2 * ef = 0 = ', e_check
cDebug          write (*,*) 'prep_elastistic_spur: ', 
cDebug     &	  	'ei = ', ei,  ' ef = ', ef
cDebug	  write (*,*)
	    n_zeroelst = n_zeroelst +1
	    iscan_zerospur(n_zeroelst) = k_real
	    ii_zerospur(n_zeroelst) = ii
	    jj_zerospur(n_zeroelst) = jj
	  endif
	 endif
        enddo
      enddo

      return
      end


c=====================================================================            
      subroutine init_qrange(qmin, qmax, qlowmin, qlowmax,
     &       nscan, kiarr, kfarr, phimax)
      
      implicit none
      include 'arraysize.inc'
c      include 'io_files.inc'

c...	common blocks

c...	functions
      
c...	passed variables
      real qmin,qmax
      real qlowmin,qlowmax
      integer nscan
      real*4 kiarr(2,maxnpts),kfarr(2,maxnpts)
      real phimax

c...	local variables
      real qlow1,qlow2
      real qex

C  Calculate the maximum obtainable Q (cosine law)
C  for the first and last point of the scan

      qex=sqrt(kiarr(1,1)**2.0+kiarr(2,1)**2.0+
     &         kfarr(1,1)**2.0+kfarr(2,1)**2.0-
     &  2*sqrt(kiarr(1,1)**2.0+kiarr(2,1)**2.0)*
     &    sqrt(kfarr(1,1)**2.0+kfarr(2,1)**2.0)*cos(phimax))

      qlow1=abs(sqrt(kiarr(1,1)**2.0+kiarr(2,1)**2.0)-
     &         sqrt(kfarr(1,1)**2.0+kfarr(2,1)**2.0))

      qmin=qex
      qmax=qex

      qex=sqrt(kiarr(1,nscan)**2.0+kiarr(2,nscan)**2.0+
     &         kfarr(1,nscan)**2.0+kfarr(2,nscan)**2.0-
     &  2*sqrt(kiarr(1,nscan)**2.0+kiarr(2,nscan)**2.0)*
     &    sqrt(kfarr(1,nscan)**2.0+kfarr(2,nscan)**2.0)*cos(phimax))

      qlow2=abs(sqrt(kiarr(1,nscan)**2.0+kiarr(2,nscan)**2.0)-
     &         sqrt(kfarr(1,nscan)**2.0+kfarr(2,nscan)**2.0))

      if (qex.gt.qmax) qmax=qex
      if (qex.lt.qmin) qmin=qex

      if (qlow1.lt.qlow2) then
         qlowmin=qlow1
         qlowmax=qlow2
      else
         qlowmin=qlow2
         qlowmax=qlow1
      endif

C      write(*,'(a,f5.2)')
C     & 'With this Ei and Ef the max attainable Q is ',qmax
     
      return
      end

cc=====================================================================            
c      subroutine init_qrange_alternative(qmin, qmax, qlowmin, qlowmax,
c     &       nscan, kiarr, kfarr, phimax)
c      
c      implicit none
c      include 'arraysize.inc'
cc      include 'io_files.inc'
c
cc...	common blocks
c
cc...	functions
c      
cc...	passed variables
c      real qmin,qmax
c      real qlowmin,qlowmax
c      integer nscan
c      real*4 kiarr(2,maxnpts),kfarr(2,maxnpts)
c      real phimax
c
cc...	local variables
c      real qlow1,qlow2
c      real qex
c      integer k
c      real q_k
c      
c      do k = 1, nscan
c        q_k = sqrt(kiarr(1,k)**2.0+kiarr(2,k)**2.0+
c     &         kfarr(1,k)**2.0+kfarr(2,k)**2.0-
c     &  2*sqrt(kiarr(1,k)**2.0+kiarr(2,k)**2.0)*
c     &    sqrt(kfarr(1,k)**2.0+kfarr(2,k)**2.0)*cos(angle4_store(k)))
c	if ((.not. A4_OK(k) .and. (angle4_store(k) .lt. acos(-1.0)) then
c	  qmin = max(qmin, q_k)
c	  qlowmin = min(qlowmin, q_k)
c	endif
c      enddo
c
c
cC  Calculate the maximum obtainable Q (cosine law)
cC  for the first and last point of the scan
c
c      qex=sqrt(kiarr(1,1)**2.0+kiarr(2,1)**2.0+
c     &         kfarr(1,1)**2.0+kfarr(2,1)**2.0-
c     &  2*sqrt(kiarr(1,1)**2.0+kiarr(2,1)**2.0)*
c     &    sqrt(kfarr(1,1)**2.0+kfarr(2,1)**2.0)*cos(phimax))
c
c      qlow1=abs(sqrt(kiarr(1,1)**2.0+kiarr(2,1)**2.0)-
c     &         sqrt(kfarr(1,1)**2.0+kfarr(2,1)**2.0))
c
c      qmin=qex
c      qmax=qex
c
c      qex=sqrt(kiarr(1,nscan)**2.0+kiarr(2,nscan)**2.0+
c     &         kfarr(1,nscan)**2.0+kfarr(2,nscan)**2.0-
c     &  2*sqrt(kiarr(1,nscan)**2.0+kiarr(2,nscan)**2.0)*
c     &    sqrt(kfarr(1,nscan)**2.0+kfarr(2,nscan)**2.0)*cos(phimax))
c
c      qlow2=abs(sqrt(kiarr(1,nscan)**2.0+kiarr(2,nscan)**2.0)-
c     &         sqrt(kfarr(1,nscan)**2.0+kfarr(2,nscan)**2.0))
c
c      if (qex.gt.qmax) qmax=qex
c      if (qex.lt.qmin) qmin=qex
c
c      if (qlow1.lt.qlow2) then
c         qlowmin=qlow1
c         qlowmax=qlow2
c      else
c         qlowmin=qlow2
c         qlowmax=qlow1
c      endif
c
cC      write(*,'(a,f5.2)')
cC     & 'With this Ei and Ef the max attainable Q is ',qmax
c     
c      return
c      end


c=====================================================================            
      subroutine plot_center(u, v, uuv, vuv, hklcen)
      implicit none
      include 'spur_layout.inc'

c...	common blocks

c...	functions
      
c...	passed variables
      INTEGER u(3),v(3)
      REAL*4 uuv(3),vuv(3)
      integer hklcen(3)

c...	local variables
      real xcen,ycen
      integer  nv, nu, i

C  When calculating the reciprocal lattice it should cover the range
C  of the spurions and the plot range.  This is easiest to do provided
C  that the reciprocal lattice and the plot have the same center.

      xcen=(xmax+xmin)/2.0
      ycen=(ymax+ymin)/2.0

C      xcen=0.0
C      ycen=0.0
      
C  the reciprocal lattice program will require an h-centre, k-centre,
C  and an l-centre to do this right.
C  v is the only thing that has a y-component
      nv=int(ycen/vuv(2))
C  which allows us to calculate the number of u vectors
      nu=int((xcen-vuv(1)*ycen/vuv(2))/uuv(1))

C  this is the hkl at roughly the centre of the plot

      do i=1,3
        hklcen(i)=nu*u(i)+nv*v(i)
      enddo
      
C      write(*,*) (hklcen(i),i=1,3)
      
      return
      end


cc==========================================================================
      subroutine plot_axes(uxvi, qztitle, rlu, right, axis,
cCommon     &      xmin, xmax, ymin, ymax,
     &      xtitle,ytitle, x2title,y2title, 
     &	    maintitle, xstar, ystar, 
     &      qmin, qmax, qlowmin, qlowmax, version_stamp)
C  setup the plot: axes, ticks, tick labels
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_color.inc'
      include 'spur_layout.inc'

c...	passed variables

      INTEGER uxvi(3)
      character*10 qztitle
      logical rlu, right
      integer axis
cCommon      real xmin, xmax, ymin, ymax
      CHARACTER*41 xtitle,ytitle,maintitle
      CHARACTER*41 x2title,y2title
      real xstar,ystar
      real qmin, qmax
      real qlowmin, qlowmax
      character*(80) version_stamp

c...	local variables
      integer len
      real temp1, temp2
      
C  setup the plot: axes, ticks, tick labels

 67   CALL PGSCH (1.5)
      call pgsls(1)  ! solid line style
      call pgsci(color_ink)  ! opposite to BG colour
c01Dec      call pgsci(1)  ! opposite to BG colour
      
      call title(uxvi,qztitle,len)  ! setup main title
      if (rlu) then
         axis=-1
      else
         axis=0
      endif


      if (right) then
        call init_pg2d(axis)
cCommon        call init_pg2d(xmin, xmax, ymin, ymax, axis)
c02Jan     &              1, axis)
c02Jan        CALL pgenv (xmin, xmax, ymin, ymax,
c02Jan     &              1, axis)
        write(maintitle,'(a31,a,a1)') 'RH Scattering Plane Normal to [',
     &    qztitle(2:len-1),']'
      else  
C  Flipping xmin and xmax in the viewport gives a left-handed spectrometer.
        call init_pg2d(axis)
cCommon        call init_pg2d(xmax, xmin, ymin, ymax, axis)
c02Jan        CALL pgenv (xmax, xmin, ymin, ymax, 
c02Jan     &              1, axis)
        write(maintitle,'(a31,a,a1)') 'LH Scattering Plane Normal to [',
     &    qztitle(2:len-1),']'
      endif

      if (rlu) then
	 temp1 = xmin/xstar
	 temp2 = xmax/xstar
c        call PGAXIS (OPT='N', X1, Y1, X2, Y2, V1, V2, 
c     :                   STEP=0.0, NSUB=1, DMAJL=0.0, DMAJR=0.5, 
c		FMIN=0.5, DISP=0.5, ORIENT=0)
         call pgaxis('N',xmin,ymin,xmax,ymin, temp1, temp2,
     &     0.0, 1, 0.5,0.5, 0.5, 0.5,0.0)

c	 ! now add the other units (A) on the top
         call pgsci(color_darkgrey)
c        call PGAXIS (OPT='N', X1, Y1, X2, Y2, V1, V2, 
c     :                   STEP=0.0, NSUB=1, DMAJL=0.0, DMAJR=0.5, 
c		FMIN=0.5, DISP=0.5, ORIENT=0)

c...	Let's put tick marks on both sides and the inside ones will be erased.
         call pgaxis('N',xmin,ymax,xmax,ymax, xmin, xmax,
     &     0.0, 1, 0.5, 0.5, 0.5, -0.5,0.0)
         call PGMTXT ('T', 1.5, 0.5, 0.5, x2title)
         call pgsci(color_ink)
	 
      	 temp1 = ymin/ystar
	 temp2 = ymax/ystar
c	 write (*,*) 'init_pg2d: ystar=', ystar, temp1, temp2, 
c     &	 xmin/xstar, xmax/xstar

         call pgaxis('N',xmin,ymin,xmin,ymax, temp1, temp2,
     &     0.0, 1, 0.5,0.5, 0.5, -0.75,0.0)

c	! now add the other units (A) on the right hand side
         call pgsci(color_darkgrey)
c...	Let's put tick marks on both sides and the inside ones will be erased.
         call pgaxis('N',xmax,ymin,xmax,ymax, ymin, ymax,
     &     0.0, 1, 0.5,0.5, 0.5, 0.5,0.0)
c         call PGMTXT (SIDE, DISP, COORD, FJUST, TEXT)
         call PGMTXT ('R', 2.25, 0.5, 0.5, y2title)
         call pgsci(color_ink)

c01Nov         call pgaxis('N',xmin,ymin,xmax,ymin,xmin/xstar,xmax/xstar,
c01Nov     &     0.0,1,0.0,0.5,0.5,0.5,0.0)
c01Nov         call pgaxis('N',xmin,ymin,xmin,ymax,ymin/ystar,ymax/ystar,
c01Nov     &     0.0,1,0.5,0.0,0.5,-0.75,0.0)
      else
c	 ! now add the other units (rlu) on the top
	 temp1 = xmin/xstar
	 temp2 = xmax/xstar
         call pgsci(color_darkgrey)
         call pgaxis('N',xmin,ymax,xmax,ymax, temp1, temp2,
     &     0.0, 1, 0.5, 0.5, 0.5, -0.5,0.0)
         call PGMTXT ('T', 1.5, 0.5, 0.5, x2title)
         call pgsci(color_ink)

c	! now add the other units (rlu) on the left hand side
      	 temp1 = ymin/ystar
	 temp2 = ymax/ystar
         call pgsci(color_darkgrey)
         call pgaxis('N',xmax,ymin,xmax,ymax, temp1, temp2,
     &     0.0, 1, 0.5,0.5, 0.5, 0.5,0.0)
         call PGMTXT ('R', 2.25, 0.5, 0.5, y2title)
         call pgsci(color_ink)

      endif

C  shade background in grey (note this overwrites the plot frame)
      call pgsfs(1)    ! a fill-style of 1 is solid
      call pgsci(color_lightgrey)   ! a grey colour 
c01Dec      call pgsci(60)   ! a grey colour 
      call PGRECT (xmin, xmax, ymin, ymax)

C  Blacken (or make it the background colour 0) a region of maximum
C  accessible Q (circle is still filled)
      
      call pgsci(color_bkgr)
c01Dec      call pgsci(0)
      call pgcirc(0.0,0.0,qmax)

      call pgsci(color_ink)     ! opposite BG color
c01Dec      call pgsci(1)     ! opposite BG color
      call pgsfs(2)     ! outline only
      call pgsls(2)     ! dashed line
      call pgcirc(0.0,0.0,qmin)
      call pgcirc(0.0,0.0,qlowmax)
      call pgsls(1)

C add something new, shade in grey the area of inaccessible low Q

      call pgsfs(1)    ! a fill-style of 1 is solid
      call pgsci(color_lightgrey)   ! a grey colour 
c01Dec      call pgsci(60)   ! a grey colour 
      call pgcirc(0.0,0.0,qlowmin)
      call pgsci(color_ink)
c01Dec      call pgsci(1)
      call pgsfs(2)

C  Have to redraw axes since they have overlaid by previous
C  commands.  The tick marks will be outside since PGPLOT_ENVOPT=I.

      call PGSCLP(0)    ! disable vprt clipping to ensure axes are visible 
      call pgrect (xmin,xmax,ymin,ymax)
      call pgsclp(1)    ! renable clipping
      call pgsfs(1)     ! solid fill style

C  Put on the previously determined labels.
      
c02Mar      CALL pglab (xtitle,ytitle,maintitle)
      CALL pglab (xtitle,ytitle,' ')
      call PGMTXT ('T', 3.0, 0.5, 0.5, maintitle)

c	Put on user name, time and date stamp
c02Jan      call pgiden()

ccc      call    sketch_layout()
      call    plot_stamp(version_stamp)    

C
      return
      end

cc==========================================================================
      subroutine plot_Q_ki_kf(qarr, kiarr, nscan)     
C  Plot Q, ki, kf  with a heavy line.

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

c...	passed variables
      real*4 qarr(2,maxnpts),kiarr(2,maxnpts)   !!,kfarr(2,maxnpts)
      integer nscan
      
c...	local variables
      integer itemp, i

C  Plot Q with a heavy line.
      
      call pgsci(color_ink)
c01Dec      call pgsci(1)
      call pgqlw(i)
      itemp = 3*i
      call pgslw(itemp)                   ! triple the thickness
c01Nov      call pgslw(3*i)                   ! triple the thickness
      CALL pgarro(0.0,0.0,qarr(1,1),qarr(2,1))
      CALL pgarro(0.0,0.0,qarr(1,nscan),qarr(2,nscan))
      call pgslw(i)


C  Plot ki and kf

      call pgsls(2)                     ! dashed line 
      CALL pgarro(0.0,0.0,kiarr(1,1),kiarr(2,1))
      CALL pgarro(kiarr(1,1),kiarr(2,1),qarr(1,1),qarr(2,1))
      call pgsls(1)

      call pgsls(2)                     ! dashed line 
      CALL pgarro(0.0,0.0,kiarr(1,nscan),kiarr(2,nscan))
      CALL pgarro(kiarr(1,nscan),kiarr(2,nscan),qarr(1,nscan),
     &            qarr(2,nscan))
      call pgsls(1)

C  End of default plot

      return
      end


c=============================================================================

      subroutine load_menuitems()
       
      implicit none
      include 'io_files.inc'
      include 'spur_menu.inc'

c...	passed variables 
      
      spurmenuitem(bit_spur) = char_spur
      spurmenuitem(bit_incoh) = char_incoh
      spurmenuitem(bit_al) = char_al
      spurmenuitem(bit_orientvec) = char_orientvec
      spurmenuitem(bit_labelvec) = char_labelvec
      spurmenuitem(bit_grid) = char_grid
      spurmenuitem(bit_mod) = char_mod

      spurmenuitem(bit_tol) = char_tol
      spurmenuitem(bit_bumpA4) = char_bumpA4
      spurmenuitem(bit_reflect) = char_reflect
cRemoved      spurmenuitem(bit_units) = char_units

      spurmenuitem(bit_zoomspur) = char_zoomspur
      spurmenuitem(bit_rescale) = char_rescale

      spurmenuitem(bit_structmenu) = char_structmenu	!! New menu
cMoved      spurmenuitem(bit_showintensities) = char_showintensities
cMoved      spurmenuitem(bit_circsize) = char_circsize
cMoved      spurmenuitem(bit_showsort) = char_showsort
cMoved      spurmenuitem(bit_sortlog) = char_sortlog
cMoved      spurmenuitem(bit_absent) = char_absent
      spurmenuitem(bit_annote) = char_annote

      spurmenuitem(bit_device) = char_device

      spurmenuitem(bit_default) = char_default
      spurmenuitem(bit_qbufout) = char_qbufout
      spurmenuitem(bit_ICP) = char_ICP
      spurmenuitem(bit_help) = char_help
      spurmenuitem(bit_quit) = char_quit
      
      return
      end

c=============================================================================

      subroutine load_menu()
       
      implicit none
      include 'io_files.inc'
      include 'spur_menu.inc'

c...	passed variables 
      
      spurmenu(bit_spur) = 
     & 'Toggle plot higher-order scattering spurions'
      spurmenu(bit_incoh) = 
     & 'Toggle plot incoherent scattering spurions'
      spurmenu(bit_al) = 
     & 'Toggle plot aluminum Ewald spheres'
      spurmenu(bit_orientvec) = 
     & 'Toggle plot u and v orientation vectors'
      spurmenu(bit_labelvec) = 
     & 'Toggle label Q, Ki and Kf wavevectors'
      spurmenu(bit_grid) = 
     & 'Toggle overlay grid on lattice plot'
      spurmenu(bit_mod) = 
     & 'Modify Q, E, and other Q-buffer values'
c     & 'Modify Q, E, and other buffer values'
      spurmenu(bit_tol) =
     & 'Change the rlv, energy spurion tolerances'
      spurmenu(bit_bumpA4) = 
     & 'BUMP! change maximum A4'
      spurmenu(bit_reflect) = 
     & 'Sf^2 for an (hkl) reflection'
cRemoved      spurmenu(bit_units) = 
cRemoved     & 'Toggle rlu - AA^-1 units'

c...	New menu item for the scan page; the spurion energy plot
      spurmenu(bit_zoomspur) = 
     & 'Zoom spurion energy plot (0,0,0,0 for auto)' 

      spurmenu(bit_rescale) = 
     & 'Rescale plot (reverse xmin and xmax for left)'

      spurmenu(bit_structmenu) = 
     & 'Presentation of structure factors menu'
c     & 'Presentation of lattice plot intensities menu'

cMoved      spurmenu(bit_showintensities) = 
cMoved     & 'Change display of Sf^2 (circles/#, grey/color)' 
ccc     & 'Change display of structure factor intensities (circles,#)' 

cMoved      spurmenu(bit_circsize) =
cMoved     & 'Change circ size (<0 vary, >0 uni, 0 no struct)'
ccc     & 'Change circ size '

cMoved      spurmenu(bit_showsort) = 
cMoved     & 'Toggle sorting of Sf^2 intensities display'
c     & 'Show sorting of intensities'

cMoved        spurmenu(bit_sortlog) = 
cMoved     & 'Select bunching S^2 = base^n or +_ tiny'
cc        spurmenu(bit_sortlog) = 
cc     & 'Bunch sorting by S^2 = base^n, rather than +_ tiny'
c      if (sortlogscale) then
c        write (spurmenu(bit_sortlog) , '(a, a, g12.5, a)') 
c     &   'Toggle logscale for S^2 sort multiplicity. ', 
c     &   '(', tolerance, ')^n.' 
c      else
c        write (spurmenu(bit_sortlog) , '(a, a, g12.5, a)') 
c     &   'Toggle logscale for S^2 sort multiplicity. ', 
c     &   '+- ', tolerance, '.' 
c      endif

cMoved      spurmenu(bit_absent) = 
cMoved     & 'Toggle hiding of the systematic absences'

      spurmenu(bit_annote) = 
     & 'Toggle display of notes next to plot'
c     & 'Toggle display of annotations next to plot'

      spurmenu(bit_device) = 
     & 'Printout/ change graphics device'

      spurmenu(bit_default) = 
     & 'Start over with defaults or new files'
c     & 'Start picture over with defaults'
      spurmenu(bit_qbufout) = 
     & 'Write out Q-buffer to a file'
c     & 'Write out to a buffer'
      spurmenu(bit_ICP) =
     & 'Generate spectrometer angles-ICP dry run'
      spurmenu(bit_help) = 
     & 'HELP!!'
      spurmenu(bit_quit) =
     & 'to quit'

      
      return
      end

c=============================================================================

      subroutine load_stndannote()
       
      implicit none
      include 'io_files.inc'
      include 'spur_menu.inc'
      include 'spur_annote.inc'

c...	passed variables 
c      character*80 stndannote(len_menu)
     
c...	local variables
      integer i
      
      nullstring = ''
      
      annoteoptions = 0
      do i = 1, bit_annote -1
        stndannote(i) = spurmenu(i)
      enddo

      do i = bit_annote, len_menu
        stndannote(i) = ' '
      enddo

      do i = len_menu + 1, max_annote
        stndannote(i) = ' '
      enddo

c...	Override with better descriptions

      stndannote(bit_spur) = 
     &'Spurions in yellow (red for elastic scattering, dE=0)'

      stndannote(bit_incoh) = 
     &'Incoh spurions from MON in cyan, ANA in blue'
c02Mar     &'Incoh spurions mon/ana(arrow) cyan(magneta)/blue(purple)'
ccc     &'Incoh spurions cyan(mon),blue(ana).Magenta/purple at peaks'

      stndannote(bit_al) = 
     &'Backgrd Al Ewald spheres in green. Red w/i tol of scan'

      stndannote(bit_orientvec) = 
     &'Xtl orientation vectors (u,v) are magenta'
     
      call set_bit(annoteoptions, bit_labelvec, 1)
      stndannote(bit_labelvec) = 
     &'Inked vectors are Q, ki, -kf at lims Ei(dash) to Ef(grey)'

c      call set_bit(annoteoptions, bit_showintensities, 1)
c      stndannote(bit_showintensities) = 
c     & 'Absences indicated by green circles. '

cccc    ' ...Show structure factor intensities (circles,#)' 
c      stndannote(bit_circsize) =
c     & '...Change circ size (<0 vary, >0 uni, 0 no struct)'

c...	debug
cDebug      do i = 1, max_annote
cDebug	write (*,*) 'load_stndannote:', i , stndannote(i)
cDebug      enddo
      
      return
      end

c=============================================================================
      subroutine load_itemvalue(tol, etol, 
cMoved     &		c, struct,
cMoved     &      gray, rainbow, plotunique, sortlogscale, tolerance, 
     &      phimax, rlu, right, selct_sp, selct_sfhkl, over,
     &	    m_rlv, nabsent, npresent, nspurtot, nelastic, nal_match,
     &	    nincohmon, nincohana, grid_style, grid_size1, grid_size2,
     &	    hide_absences)
       

      implicit none
      include 'io_files.inc'
      include 'spur_menu.inc'
      include 'spur_layout.inc'
c      include 'arraysize.inc'
c      include 'spur_scan.inc'

c...	functions
      logical bit_on
      integer ilen

c...	passed variables 
      real tol, etol, phimax
cMoved      real c, tolerance
cMoved      logical gray, rainbow
cMoved      logical plotunique, sortlogscale, struct
      logical rlu, right, over
      integer selct_sp(3)
      real	selct_sfhkl
      integer 	m_rlv, nabsent, npresent, nspurtot, nelastic
      integer  nal_match
      integer	nincohmon, nincohana
      character*1 grid_style		!!c: cartesian, o:orientation vectors
      real	grid_size1, grid_size2	!! step size in rlu, for x/u, y/v
      logical	hide_absences
      
c...    local variables
      integer i
      logical debug
      parameter (debug = .false.)
c      parameter (debug = .true.)
      
cDebug      if (debug) write (*,*) 'load_itemvalue: null strings'
      do i = 1, len_menu 
        spuritemvalue(i) = ''
      enddo
cDebug      if (debug) write (*,*) 'load_itemvalue: null strings loaded'

      write (spuritemvalue(bit_default), '(a)')  
     &      file_bufin(1:min(20,ilen(file_bufin)))
c02May      write (spuritemvalue(bit_default), '(a)')  
c02May     &      file_bufin(1:ilen(spuritemvalue))
c02Mar      write (spuritemvalue(bit_default), '(a)')  
c02Mar     &      file_bufin	
cDebug      if (debug) write (*,*)'load_itemvalue: ', 
cDebug     &      'report_menuitem(bit_default)'
      if (debug) call report_menuitem(bit_default)
      
c      if (debug) write (unitdebug, *) 'load_itemvalue: tol=', tol 
c      write (spuritemvalue(bit_tol), '(g12.5, a6)')  tol, ' AA^-1'
      write (spuritemvalue(bit_tol), '(f6.3, a5, f5.2, a4)')  
     &      tol, ' A^-1 ', etol, ' meV'	
c      write (spuritemvalue(bit_tol), '(g12.5)')  tol
cDebug      if (debug) write (*,*)'load_itemvalue: report_menuitem(bit_tol)'
      if (debug) call report_menuitem(bit_tol)
      
c      if (debug) write (unitdebug, *) 'load_itemvalue: gray=', gray 
cMoved      if (.not. struct) then
cMoved        write (spuritemvalue(bit_showintensities), *) 
cMoved     &		'No struct circles'
cMoved      else if ((gray) .and. (c .lt. 0)) then
cMoved	if (rainbow) then
cMoved          write (spuritemvalue(bit_showintensities), *)
cMoved     &		'varying colr circles'
cMoved	else
cMoved          write (spuritemvalue(bit_showintensities), *)
cMoved     &		'varying gray circles'
cMoved	endif
cMoved      else if ((gray) .and. (c .ge. 0)) then
cMoved	if (rainbow) then
cMoved          write (spuritemvalue(bit_showintensities), *) 
cMoved     &		'uniform colr circles'
cMoved        else
cMoved          write (spuritemvalue(bit_showintensities), *) 
cMoved     &		'uniform gray circles'
cMoved	endif
cMoved      else
cMoved	if (rainbow) then
cMoved          write (spuritemvalue(bit_showintensities), *) 'gray #s'
cMoved	else
cMoved          write (spuritemvalue(bit_showintensities), *) 'color #s'
cMoved	endif
cMoved      endif
cMoved      if (debug) call report_menuitem(bit_showintensities)
      
c      if (debug) write (unitdebug, *) 'bit_circsize=>', c, struct
cMoved      write (spuritemvalue(bit_circsize), 90) c, struct
cMoved 90   format( f10.3, ' struct=', L1)
c      write (spuritemvalue(bit_circsize), *) c
cMoved      if (debug) call report_menuitem(bit_circsize)

c      if (debug) write (unitdebug, *) 'load_itemvalue: plotunique', plotunique
cMoved      if (plotunique) then
cMoved        write (spuritemvalue(bit_showsort), *) 'ABC..cba labels'
cMoved      else
cMoved        write (spuritemvalue(bit_showsort), *) 'no sorting shown '
cMoved      endif
cMoved      if (debug) call report_menuitem(bit_showsort)

cMoved      if (hide_absences) then
cMoved        write (spuritemvalue(bit_absent), *) 'absences hidden'
cMoved      else
cMoved        write (spuritemvalue(bit_absent), *) 'absences crosses'
cMoved      endif
cMoved      if (debug) call report_menuitem(bit_absent)

cDebug      write (unitdebug, *) 'load_itemvalue: sort=', sortlogscale,
cDebug     &      	tolerance  
c        spurmenu(bit_sortlog) = 
c     & 'Bunch sorting by S^2 = base^n, rather than +_ tiny'
cMoved      if (sortlogscale) then
cMoved        write (spuritemvalue(bit_sortlog) , '(a, g12.5, a)') 
cMoved     &   '(', tolerance, ')^n' 
cMoved      else
cMoved        write (spuritemvalue(bit_sortlog) , '(a, g12.5)') 
cMoved     &   '+- ', tolerance 
cMoved      endif
cMoved      if (debug) call report_menuitem(bit_sortlog)

cDebug      write (unitdebug, *) 'load_itemvalue: phimax=', phimax
      write (spuritemvalue(bit_bumpa4), '(g12.5, a)') 
     &      phimax*180.0/acos(-1.0), ' degs' 
      if (debug) call report_menuitem(bit_bumpa4)

cDebug      write (unitdebug, *) 'load_itemvalue: rlu=', rlu 
cRemoved      if (rlu) then
cRemoved        write (spuritemvalue(bit_units),*) 'rlu'
cRemoved      else
cRemoved        write (spuritemvalue(bit_units),*) 'AA^-1'
cRemoved      endif
      
cDebug      write (unitdebug, *) 'load_itemvalue: right=', right
      if (auto_xmin_sc)  then 
        write (spuritemvalue(bit_zoomspur)(1:4), '(i3, 1x)') 0
      else
        write (spuritemvalue(bit_zoomspur)(1:4),'(i3,1x)') nint(xmin_sc)
cDebug	write (*,*) 'load_itemvalue: xmin_sc=', xmin_sc
      endif

      if (auto_xmax_sc)  then 
        write (spuritemvalue(bit_zoomspur)(5:8), '(i3, 1x)') 0
      else
        write (spuritemvalue(bit_zoomspur)(5:8),'(i3,1x)') nint(xmax_sc)
cDebug 	write (*,*) 'load_itemvalue: xmax_sc=', xmax_sc
      endif
      
      if (auto_ymin_sc_hs)  then 
        write (spuritemvalue(bit_zoomspur)(9:14),'(f5.1,1x)') 0.0
      else
	if (abs(ymin_sc_hs) .ge. 100.0) then
cDebug	  if (debug) write (*,*) 'load_itemvalue: ymin_sc_hs=', ymin_sc_hs
          write (spuritemvalue(bit_zoomspur)(9:14),'(i5,1x)') 
     &		  nint(ymin_sc_hs)
	else if (abs(ymax_sc_hs) .ge. 10.0) then
          write (spuritemvalue(bit_zoomspur)(9:14),'(f5.1,1x)')
     &		  ymin_sc_hs
	else
          write (spuritemvalue(bit_zoomspur)(9:14),'(f5.2,1x)')
     &	  	ymin_sc_hs
	endif
c02Mar        write (spuritemvalue(bit_zoomspur)(9:14),'(f5.1,1x)') ymin_sc_hs
      endif
      
      if (auto_ymax_sc_hs)  then 
        write (spuritemvalue(bit_zoomspur)(15:20),'(f5.1,1x)') 0.0
      else
	if (abs(ymax_sc_hs) .ge. 100.0) then
          write (spuritemvalue(bit_zoomspur)(15:20),'(i5,1x)') 
     &		  nint(ymax_sc_hs)
	else if (abs(ymax_sc_hs) .ge. 10.0) then
          write (spuritemvalue(bit_zoomspur)(15:20),'(f5.1,1x)')
     &		  ymax_sc_hs
	else
          write (spuritemvalue(bit_zoomspur)(15:20),'(f5.2,1x)')
     &	  	ymax_sc_hs
	endif
      endif
      if (debug) call report_menuitem(bit_zoomspur)

c      if      ((     auto_ymin_sc_hs) .and. (auto_ymin_sc_hs)) then
c        write (spuritemvalue(bit_zoomspur),*) 'auto auto'
c      else if ((.not.auto_ymin_sc_hs) .and. (auto_ymin_sc_hs)) then
c        write (spuritemvalue(bit_zoomspur),*) ymin_sc_hs, ' auto'
c      else if ((     auto_ymin_sc_hs) .and.(.not.auto_ymin_sc_hs)) then
c        write (spuritemvalue(bit_zoomspur),*) 'auto ', ymax_sc_hs
c      else if ((.not.auto_ymin_sc_hs).and.(.not.auto_ymin_sc_hs)) then
c        write (spuritemvalue(bit_zoomspur),*) ymin_sc_hs, ymax_sc_hs
c      endif


      if (right) then
c					      123456789_123456789_	
        write (spuritemvalue(bit_rescale),*) 'Right handed spectm'
      else
        write (spuritemvalue(bit_rescale),*) 'Left handed spectm'
      endif
      if (debug) call report_menuitem(bit_rescale)

cDebug      write (unitdebug, *) 'load_itemvalue: selct_sfhkl=', selct_sfhkl      
      write (spuritemvalue(bit_reflect), 101)
C     &      '( '(', I2, ',' I2, ',', I2, ')' )')
     &       (selct_sp(i), i = 1, 3), selct_sfhkl
101   format ('(', I2, ',', I2, ',', I2, ')=', g9.2 ) !! 5char + 6I +9=20
      if (debug) call report_menuitem(bit_reflect)

cDebug      write (unitdebug, *) 'load_itemvalue: over=', over
      if (over) then
        write (spuritemvalue(bit_labelvec),*) 'Overwrite labels'
      else
        write (spuritemvalue(bit_labelvec),*) ' '
      endif
      if (debug) call report_menuitem(bit_labelvec)

      if ( bit_on(currentoptions, bit_grid)) then
        write (spuritemvalue(bit_grid),150) grid_style, 
     &		grid_size1, grid_size2
 150    format (a1, ':', 1x,  f5.3, 1x, f5.3)
      else
        write (spuritemvalue(bit_grid),*) ' '
      endif

cDebug      write (unitdebug, *) 'load_itemvalue: qbufout='
      if ( bit_on(currentoptions, bit_qbufout)) then
        write (spuritemvalue(bit_qbufout),*) 
     &		file_bufout(1:min(20,ilen(file_bufout)))
c02May        write (spuritemvalue(bit_qbufout),*) file_bufout
      else
        write (spuritemvalue(bit_qbufout),*) 'no bufout yet'
      endif
      if (debug) call report_menuitem(bit_qbufout)
      
cDebug      write (unitdebug, *) 'load_itemvalue: icp'
      if ( bit_on(currentoptions, bit_icp)) then
        write (spuritemvalue(bit_icp),*)
     &		 file_icp(1:min(20,ilen(file_icp)))
c02May        write (spuritemvalue(bit_icp),*) file_icp
      else
        write (spuritemvalue(bit_icp),*) 'no icp yet'
      endif
      if (debug) call report_menuitem(bit_icp)

cDebug      write (unitdebug, *) 'load_itemvalue: nabsent=', nabsent
      write (spuritemvalue(bit_orientvec),205) nabsent, npresent
c02Jan      write (spuritemvalue(bit_orientvec),205) nabsent, m_rlv*m_rlv
 205   format ('nabsent=', i5, '/',   i5)     !! 8 + 1 + 5+5
c      write (spuritemvalue(bit_orientvec),*) 'nabsent=', nabsent, 
c     &         '/', m_rlv*m_rlv
cc     &		, 'ntiny=', ntiny
      if (debug) call report_menuitem(bit_orientvec)

cDebug      write (unitdebug, *) 'load_itemvalue: nspur=', nspurtot
      if ( bit_on(currentoptions, bit_spur)) then
        write (spuritemvalue(bit_spur),210) nspurtot, nelastic
 210    format ('spur=', i5, 1x,  'els=', i5)     !!5+5+1+4+5=20
c        write (spuritemvalue(bit_spur),*) 'nspur=', nspurtot, 
c     &		' nelast=', nelastic
      else
        write (spuritemvalue(bit_spur),*) ' '
      endif
      if (debug) call report_menuitem(bit_spur)

      if ( bit_on(currentoptions, bit_al)) then
        write (spuritemvalue(bit_al),*) 'nal_match=', nal_match
      else
        write (spuritemvalue(bit_al),*) ' '
      endif
cDebug      write  (unitdebug, *) 'load_itemvalue: ', spuritemvalue(bit_al)
      if (debug) call report_menuitem(bit_al)
      
      if ( bit_on(currentoptions, bit_incoh)) then
        write (spuritemvalue(bit_incoh),220) 'mon=', 
     &		nincohmon, 'ana=', nincohana
220     format (a4, i5, 1x, a4, i5)
c        write (spuritemvalue(bit_incoh),220) 'mon, ana=',
c     &		nincohmon, nincohana
c220     format (a9, i5, 1x, i5)
      else
        write (spuritemvalue(bit_incoh),*) ' '
      endif
cDebug      write  (unitdebug, *) 'load_itemvalue: ', spuritemvalue(bit_incoh)
      if (debug) call report_menuitem(bit_incoh)

      if ( bit_on(currentoptions, bit_mod)) then
        write (spuritemvalue(bit_mod),*) 'Modified' 
      else
        write (spuritemvalue(bit_mod),*) ' '
      endif

cDebug      write  (unitdebug, *) 'load_itemvalue: ', spuritemvalue(bit_mod)
      if (debug) call report_menuitem(bit_mod)


      if (debug) write (*,*) 'load_itemvalue: end of routine'

      return
      end

c=============================================================================
      subroutine report_menuitem(i)
      implicit none
      include 'io_files.inc'
      include 'spur_menu.inc'
      integer i
c      character*(*) astring
      
      write (unitdebug, *) 'report_menuitem: ', 
     &	      spurmenuitem(i),': ',  spuritemvalue(i)
      end
c=============================================================================

      subroutine display_menu(aversion, scale, spurhist) 
       
      implicit none
      include 'io_files.inc'
      include 'spur_menu.inc'
      include 'spur_warn.inc'

c...	functions
      logical bit_on
      integer ilen
      
c...	passed variables
      character*(*) aversion
      CHARACTER*2 scale
c      CHARACTER*5 scale
      integer spurhist

c...	local variables
      integer i, bit_value

c02Mar 220   write (*,'(A,$)') 'Type a command or hit <return> for options: '
c02Mar      read (*,'(a)',err=220) scale      ! read in a one character command
c      write(*,'(a,i3)') 'spurhist=',spurhist

c      call version(aversion)
c	Note: z/o in VMS fortran produces hexadecimal/octal output.      

c      write (*,'(a30, x, a24, x, a24, x, 0)') aversion,  
c     &		FILE_BUFIN, FILE_LPI,
c     &    	currentoptions
        
c02Mar      if ((scale(1:1).eq.' ') .or. (scale(1:1).eq.'?')) then

        write (*,*) aversion
        write (*,'(a25, x, a25, x, a25)') FILE_BUFIN, FILE_LPI,
     &		file_base(1:ilen(file_base))//'.<ext>'
        if (warn_options .ne. 0) write (*,'(a)') warnstring
c03Mar        write (*,'(a38, x, a38)') FILE_BUFIN, FILE_LPI  
ccc        write (*,'(a30, x, a24, x, a24)') aversion, FILE_BUFIN, FILE_LPI  
c      call display_spurhist(spurhist, 0)
c      call display_spurhist(spurhist, unitdebug)

c        write (*,*) 'currentoptions =', currentoptions, 
c     &		'addoption=', addoption
cc      write (unitdebug,*) 'currentoptions =', currentoptions 
      
        do i = 1, len_menu
	  write (*,*) spurmenuitem(i), ': ', spurmenu(i),
     &	   ' ', bit_on(currentoptions, i), bit_on(addoption, i), 
     &     ' ', spuritemvalue(i)
c	  write (*,*) spurmenuitem(i), ' :  ', spurmenu(i),
c     &	   ' ', bit_on(currentoptions, i),
c     &     ' ', spuritemvalue(i)
        enddo
222     write(*,'(a,$)') 'Command ("QQ" to quit)? '
c222     write(*,'(a,$)') 'Command ("q" to quit)? '
        read (*,'(a)',err=222) scale
c02Mar      endif

      addoption = 0
      do i =1, len_menu
c        if (spurmenuitem(i) .eq. scale(1:1)) then
        if (spurmenuitem(i) .eq. scale) then
	   call set_bit(addoption,i,1)
	   if (i .eq. bit_spur) then
	     call toggle_bit(currentoptions,i, bit_value)
	   else if (i .eq. bit_al) then
	     call toggle_bit(currentoptions,i, bit_value)
	   else if (i .eq. bit_incoh) then
	     call toggle_bit(currentoptions,i, bit_value)
	   else if (i .eq. bit_orientvec) then
	     call toggle_bit(currentoptions,i, bit_value)
	   else if (i .eq. bit_labelvec) then
	     call toggle_bit(currentoptions,i, bit_value)
	   else if (i .eq. bit_grid) then
	     call toggle_bit(currentoptions,i, bit_value)
ccc	   else if (i .eq. bit_showsort) then
ccc	     call toggle_bit(currentoptions,i, bit_value)
	   else if (i .eq. bit_annote) then
	     call toggle_bit(currentoptions,i, bit_value)
	   else
	     call set_bit(currentoptions,i,1)
	   endif
	endif
      enddo

      return
      end
     
c=============================================================================

      subroutine display_spurhist(spurhist, unit)
c...	Use UPPERCASE letters to indicate that the condition is true (bit=1)
c....   Use lowercase letters to indicate that the condition is false (bit=0)      

C     spurion on or off      2^0
C     incoherent on or off   2^1
C     aluminum on or off     2^2
C     prox (0) or all (1)    2^3
C     orient vec on or off   2^4
C     label wavevectors      2^5
C     overwrite (yes/no)     2^6
C     Systematic absences    2^7
c
c123456789_123456789_123456789_123456789_123456789_123456789_123456789_123456789
cspurion   incohernt aluminum  nearby    orient_vc label_wav overwrite absences 


      implicit none
c...	passed variables
      integer spurhist, unit
      
c...	local variables
      integer ipwr, spurbit
      character*10 ahist(0:7)
      
      data ahist /'spurion   ', 'incohernt ', 'aluminum  ', 
     &      	  'nearby    ', 'orient_vc ', 'label_wav ', 
     & 		  'overwrite ', 'absences  '/
     
      do ipwr = 0, 7
        call bit(spurhist,ipwr,spurbit,0)
	 if (spurbit.eq.1) then
	   call supper(ahist(ipwr)) 
	 else
	   call lower(ahist(ipwr)) 
	 endif
      enddo
      if (unit .eq. 0) then
        write (*,'(8(a))') (ahist(ipwr), ipwr = 0, 7)
      else
        write (unit,'(8(a))') (ahist(ipwr), ipwr = 0, 7)
      endif

      return
      end
      
      
c=========================================================================

      subroutine cmd_spurion( m_rlv, npts, nscan, nspurtot, nspur, 
     &      qspur, espur, qsparr, esparr, kiarr, kfarr,
     &      circsize, ki, kf, zeta, eta, hh, kk, ll, 
     &      absent, replot, sp, spurtitle, 
     &      ec, es, h, hc, hs,  constq, nelastictot, noplot, bmatuv)
C  Plot, count, and record spurious Q's

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

c...	common blocks
      real tol, etol		!! spatial and energy tolerances		!! implicit       
      common /spur/ tol, etol

c...	functions
      logical	invert
cUnused      integer 	ilen
      
c...	passed variables
      integer  m_rlv, npts, nscan, nspurtot, nspur
      real qspur(2,n_spur_max,n_spur_max)
      real espur(n_spur_max,n_spur_max)
      real qsparr(2,n_spur_max,n_spur_max, maxnpts)
      real esparr(n_spur_max,n_spur_max,maxnpts)
      real*4 kiarr(2,maxnpts),kfarr(2,maxnpts)
      real circsize
      real ki(2),kf(2)
      real zeta(m_rlv_max,m_rlv_max),eta(m_rlv_max,m_rlv_max)
      integer hh(m_rlv_max,m_rlv_max)
      integer kk(m_rlv_max,m_rlv_max)
      integer ll(m_rlv_max,m_rlv_max)
      logical absent(m_rlv_max,m_rlv_max)
      logical replot
      integer sp(3)
      character*10 spurtitle
      real ec, es
      REAL*4 h(3),hc(3),hs(3)
      logical constq
      integer nelastictot
c      real bmat(3,3)
      real bmatuv(3,3)
      logical noplot
            
c...	local variables
      integer iscan, i, j, l
      integer ii, jj
      real bmatuv_inv(3,3)
      integer nelastic

cUnused      integer len
cUnused      character*13 a_hkl		! only integers for q
cReport      character*19 afloat_hkl
cReport      character*(80)	line_scan, line_spur, line_info
cReport      logical firstbragg
      
      write(*,*) '&&& Elastic spurions and ', 
     &	      'higher order spurions matching Bragg peaks at:'
      write(unitspur,*) '&&& All higher order spurions:'

      if (noplot) then
        do i =1, 3
	  do j = 1, 3
	    bmatuv_inv(i,j) = dble(bmatuv(i,j))
          enddo
        enddo
cOpps      if (.not. invert(3, bmatuv)) then
        if (.not. invert(3, bmatuv_inv)) then
	  write (*,*) 'cmd_spur: Error inverting bmatuv'
        endif

         nspurtot=0
         nelastictot=0
      endif
      
         do iscan=1,nscan

cReport           do j=1,3
cReport             h(j)=hc(j)+(k-(float(npts)+1.0)/2.0)*hs(j)
cReport           enddo
cReport	   write (line_scan, 115) k, ec+(k-(float(npts)+1.0)/2.0)*es, 
cReport     &		   (h(j),j=1,3) 	   
cReport115	   format ('Scan point #=', i3, ' at dE=', f7.3, 
cReport     &		' at Q=(',f6.2,1x,f6.2,1x,f6.2,')')
cReport           write (unitspur, '(a)') line_scan
cReport           write (*, '(a)') line_scan
           

C setup variables again
            do ii=1,n_spur_max
             do jj=1,n_spur_max
                  do l=1,2
                     qspur(l,ii,jj)=qsparr(l,ii,jj,iscan)
                  enddo
                  espur(ii,jj)=esparr(ii,jj,iscan)
	      enddo
	    enddo


            do l=1,2
               ki(l)=kiarr(l,iscan)
               kf(l)=kfarr(l,iscan)
            enddo

	    if (noplot) then
              call prepspur(n_spur_max, nspurtot,m_rlv,qspur,espur,
     &		 nspur, zeta, eta, hh,kk,ll,
     &		 ki,kf,absent, replot, nelastic, nelastictot, iscan)

c...		New location. Not used by plotspur or reportspur.
              nelastictot=nelastictot + nelastic
              nspurtot=nspurtot+nspur
c	      write (*,*) 'cmd_spurion: increasing nelastictot ',
c     &		       nelastictot , nelastic
c	      write (*,*) 'cmd_spurion: increasing nspurtot ',
c     &		       nspurtot , nspur
	    endif

	    if (.not. noplot) then
             call plotspur(n_spur_max, m_rlv,qspur,espur,
     &		 circsize,zeta,eta, 
     &		 ki,kf,absent, replot,  iscan)
            endif
	    
	    if (noplot) then
	      call reportspur(n_spur_max, nspurtot, nelastictot, 
     &		qspur,espur, bmatuv_inv, ec, es, h, hc, hs, nscan, iscan)
     
     	    endif

cReportC  		Print out the spurions
cReport
cReportc            if ((nspur.ne.0).and.(nspurtot.eq.0)) 
cReportc     &           write(*,'(a19)') 'Spurion matches at:'
cReport
cReport	    firstbragg = .true.
cReport            do ii=1,n_spur_max
cReport             do jj=1,n_spur_max
cReportc	        if (( ii .ne. 1) .and. (jj .ne. 1)) then
cReport	          call hkl_float(qspur(1,ii,jj), qspur(2,ii,jj), 
cReport     &		    bmatuv_inv, afloat_hkl)
cReport     		  write (line_spur, 124) ii, jj, afloat_hkl, espur(ii,jj)
cReport124		  format (1x,i2, '*Ki - ', i2, '*Kf =',1x,
cReport     &			a,' at E_spur=', f8.2)

cReport		  if ( k .eq. nint(xmin_sc)) then
cReport		    afloat_hkl_start(ii,jj) = afloat_hkl
cReport		  endif
cReport		  if ( k .eq. nint(xmax_sc)) then
cReport		    afloat_hkl_end(ii,jj) = afloat_hkl
cReport		  endif
ccReport		endif

cReport                  if (abs(espur(ii,jj)).lt. etol) then
cReport		    write (line_info,*) ' Elastic red. '
cReportc...		    Override  line_info, if there was match to a Bragg peak.		    
cReport		    if (nspur .ne. 0) then
cReport                     do i=nspurtot+1,nspurtot + nspur
cReport		      if ((ispur(i+nspurtot).eq.ii) 
cReport     &			 .and.(jspur(i+nspurtot).eq.jj)) then
cReport     			if (firstbragg) then
cReport     			  firstbragg = .false.
cReport			  write (*,*)
cReport			  write (*, '(a80)') line_scan
cReport			endif
cReport		        write (line_info,*) ' Elastic red.   Bragg'
cReport		      endif
cReport		     enddo
cReport		    endif
ccReport		report all those higher order spurions which are eleastic
cReport	            if (( ii .ne. 1) .and. (jj .ne. 1)) then
cReport                      write (*, *) line_spur(1:ilen(line_spur)), 
cReport     &		  	line_info(1:ilen(line_info))
cReport     		    endif
cReport
cReport		  else
cReport		    write (line_info,*) ' Inelst yellow.'
cReportc...		    Override  line_info, if there was match to a Bragg peak.		    
cReport		    if (nspur .ne. 0) then
cReport                     do i=nspurtot+1,nspurtot + nspur
cReport		      if ((ispur(i+nspurtot).eq.ii) 
cReport     &			 .and.(jspur(i+nspurtot).eq.jj)) then
cReport     			if (firstbragg) then
cReport     			  firstbragg = .false.
cReport			  write (*,*)
cReport			  write (*, '(a)') line_scan
cReport			endif
cReport		        write (line_info,*) ' Inelst yellow. Bragg'
ccReport		        write (line_info,*) ' Inelastic (yellow) Bragg'
cReportc		Only report those higher order spurions which hit Bragg peaks
cReport	                if (( ii .ne. 1) .and. (jj .ne. 1)) then
cReport                          write (*, *) line_spur(1:ilen(line_spur)), 
cReport     &			  	line_info(1:ilen(line_info))
cReport     			endif
cReport		      endif
cReport		     enddo
cReport		    endif
cReport		  endif
cReport
cReportc		  Write all higher order spurions to the file *.spur		  
cReport                  write (unitspur, *) line_spur(1:ilen(line_spur)), 
cReport     &		  	line_info(1:ilen(line_info))
cReport
cReport
cReport 26                 format(a, f7.3, $)
cReport 637                 format(a,f6.2,1x,f6.2,1x,f6.2,a, $)
cReport
cReport
cReport
cReport             enddo
cReport            enddo
cReport	    
cReport 24         format(1x,i2, '*Ki - ', i2, '*Kf =', $)
cReport 25         format(1x,a,a,f8.2,a,i2,$)
cMoved            nspurtot=nspurtot+nspur


cReport         write(unitspur,*) 		!! blank line between scan points
cReport         write(*,*) 		!! blank line between scan points

         enddo		!! iscan = 1, nscan
	 
cMoved         write(*,*) 		!! blank line at end of scan points

C  No spurions?
cReport	write(*,*) nspurtot, ' matches of higher order spurions ', 
cReport     &		'to Bragg peaks, out of ', n_spur_max*n_spur_max*nscan 
cReport	write(unitspur,*) nspurtot, ' matches of higher order spurions ', 
cReport     &		'to Bragg peaks, out of ', n_spur_max*n_spur_max*nscan 

cReport	write(*,*) nelastictot, ' elastic higher order spurions, ', 
cReport     &		'out of ', n_spur_max*n_spur_max*nscan 
cReport	write(unitspur,*) nelastictot, ' matches of higher order spurions, ', 
cReport     &		'out of ', n_spur_max*n_spur_max*nscan 


cReport      write(*,*) ''
cReport      write(unitspur,*) ''

c      write (*,*)  'cmd_spurion: goal points: ',
c     &	      nint(xmin_sc), nint(xmax_sc)
c            do ii=1,n_spur_max
c             do jj=1,n_spur_max
c	       write (*,*) 'cmd_spurion: labels =', ii, jj, 
c     &		 '"', afloat_hkl_start(ii,jj),
c     &		 '" "', afloat_hkl_end(ii,jj), '"'
c	     enddo
c	    enddo


      call rank_higher(nspurtot, nelastictot)
      
      return
      end


c==========================================================================
      subroutine cmd_incoherent(nscan, m_rlv, zeta, eta,
     &	kiarr, kfarr, incohmon,incohana, incmonarr, incanaarr, 
     &  hh,kk,ll,circsize,ki,kf,constq,h,incohflag,absent,replot, 
     &  npts, hc, hs, nincohmon, nincohana, pleaseplot)
C  Plot and report spurious Q's due to incoherent scattering.
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	passed variables
      integer nscan, m_rlv
      real zeta(m_rlv_max,m_rlv_max)
      real eta(m_rlv_max,m_rlv_max)
      real incohmon(2),incohana(2)
      real*4 incmonarr(2,maxnpts),incanaarr(2,maxnpts)
      real*4 kiarr(2,maxnpts),kfarr(2,maxnpts)
      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,incohflag,replot
      logical absent(m_rlv_max,m_rlv_max)
      integer npts
      REAL*4 hc(3),hs(3)
      integer nincohmon, nincohana
      logical pleaseplot
      
c...	local variables
      integer k, l, j
      logical  flagincohmon,flagincohana
            
      write (*,*) '&&& Incoherent spurions matches at:'
      write (unitspur,*) 'Incoherent spurions matches at:'
      
         incohflag=.false.
         nincohmon = 0
         nincohana = 0
         do k=1,nscan
            do l=1,2
               incohmon(l)=incmonarr(l,k)
               incohana(l)=incanaarr(l,k)
               ki(l)=kiarr(l,k)
               kf(l)=kfarr(l,k)
            enddo
            do j=1,3
               h(j)=hc(j)+(k-(float(npts)+1.0)/2.0)*hs(j)
            enddo
            call prepincoh(m_rlv,zeta,eta,incohmon,incohana, hh,kk,ll,
     &  	circsize,ki,kf,k,constq,h,incohflag,absent, replot,
     &		nincohmon, nincohana,  flagincohmon,flagincohana)
     	    if (pleaseplot) then
              call plotincoh(m_rlv,zeta,eta,incohmon,incohana, hh,kk,ll,
     &  	circsize,ki,kf,k,constq,h,incohflag,absent, replot,
     &		nincohmon, nincohana,  flagincohmon,flagincohana)
     	    endif
        enddo
        if (.not.incohflag) then
           write(*,*) 'No incoherent spurions found.'
           write(unitspur,*) 'No incoherent spurions found.'
c           write(8,*) 'No incoherent spurions found.'
        endif

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


      return
      end


c==========================================================================
      subroutine cmd_aluminum(nscan, npts, minq, maxq, magq,
     &   nearby, qarr, kiarr, kfarr, 
     &	 h, hc, hs, incmonarr, incanaarr,
     &	 constq, spurhist, mon, pi, fe, phimax, plotlattice)
C  Plot aluminum Ewald spheres that are within qmax in green
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_buf.inc'
      include 'spur_color.inc'
      include 'spur_scan.inc'

c...	common blocks
      real tol, etol		!! spatial and energy tolerances		!! implicit       
      common /spur/ tol, etol

c...	functions
      character*10	name_color
      
c...	passed variables
      integer nscan, npts
c02Apr      real gal(nal)
      real minq,maxq,magq
c02Apr      integer hklal(3,nal)
      logical nearby
      real*4 qarr(2,maxnpts),kiarr(2,maxnpts),kfarr(2,maxnpts)
      REAL*4 h(3),hc(3),hs(3)
      real*4 incmonarr(2,maxnpts),incanaarr(2,maxnpts)
      logical constq
      integer spurhist
      CHARACTER*5 mon
      real pi
      real fe
      real phimax
      logical plotlattice

cInc      logical alspur
cInc      integer nal_match
      
c...	local variables
      integer i, j, jj, iscan
      integer spurbit
      real angle4, asin_angle4		!! A4 angle for Al (not sample!)
      character*10 thiscolor
      character*80 aline
c      logical angle4_OK
      logical thismatches
c      character*1 beforeOrAfter


C  Calculate the aluminum Ewald spheres.  Use text output
C  to give the one in closest proximity to current Q.

         call alewald()
c         call alewald(gal,hklal)

c...	Check min and max of q's in the scan range. (So can test Al q's later.)
 	iscan = 1
        magq=sqrt(qarr(1,iscan)**2+qarr(2,iscan)**2)
	minq=magq
	maxq=magq
        do iscan=2,nscan
          magq=sqrt(qarr(1,iscan)**2+qarr(2,iscan)**2)
          if (magq.lt.minq) minq=magq
          if (magq.gt.maxq) maxq=magq
	enddo


        write (*,*) '&&& Al Eward spheres at:'
        write (unitspur,*)'&&& Al Eward spheres at:'
C  report locations of Al Alewald spheres first, before looking for matches
        do iscan=1,nscan
          magq=sqrt(qarr(1,iscan)**2+qarr(2,iscan)**2)
	  do i = 1, nal

	    angle4 = 190.0		!! Nonsense value
c...	     AL Alewald spheres fixed in .rl.u.
c		But their A4 angle depends on Ei.	    
            if (mon.eq.'(MON)') then 
	      asin_angle4 =		!!  Laue scattering law
     &	      	0.5*gal(i)/sqrt(kiarr(1,iscan)**2+kiarr(2,iscan)**2)
            else
	      asin_angle4 =		!!  Laue scattering law
     &	      	0.5*gal(i)/sqrt(kfarr(1,iscan)**2+kfarr(2,iscan)**2)
            endif

c 	    if (i .eq. 1) then		!! debug
c              magq=sqrt(qarr(1,iscan)**2+qarr(2,iscan)**2)
c	      write (*,*) 'cmd_aluminum: i, ki= ', iscan,
c     &	     	kiarr(1,iscan), kiarr(2,iscan), '|q|=',magq 
c	      write (*,*) 'cmd_aluminum: i, kf= ', iscan,
c     &	     	kfarr(1,iscan), kfarr(2,iscan), '|q|=',magq 
c	    endif			!! debug
	    
	    if (abs(asin_angle4  ) .le. 1.0) then
	      angle4 = 360/pi* asin(asin_angle4 )
	    endif


c02Apr 541        format(a,3i1,2a,f5.2,a,$)
c02Apr            if ((gal(i).ge.(minq-.5)).and.(gal(i).le.(maxq+.5)))

	     if ( (draw_al_iscan .eq. 0) .or. 
     &	         (draw_al_iscan .eq. iscan) )  then
c	     if ( ((draw_al_iscan .eq. 0) .and. (iscan .eq. 1)) .or. 
c     &	         ((draw_al_iscan .eq. 0) .and. (iscan .eq. nscan)) .or. 
c     &	         (draw_al_iscan .eq. iscan) )  then
 
              if (angle4 .le. phimax*180.0/pi)  then
cOld                if ((gal(i) .ge. minq) .and. (gal(i) .le. maxq)) then
cOld 540             format ('Al(', 3i1, ')', ' at |Q|=', f5.2, 
cOld     &			' AA^-1 or A4= ', f6.2, 
cOld     &			' at fixed E', f4.1, 1x, a5,  ' pt#=', i3, a)
c 540             format ('Al(', 3i1, ')', ' at |Q|=', f5.2, 
c     &			' AA^-1 or A4= ', f6.2, 
c     &			' at fixed E', f4.1, 1x, a,  ' pt#=', i3, a)

c		 beforeOrAfter = '~'
cOld		 write (aline, 540) (hklal(j,i),j=1,3), gal(i), 
cOld     &		 	angle4,fe, mon, iscan,' InsideScan' 
c		 write (aline, 540) (hklal(j,i),j=1,3), gal(i), 
c     &		 	angle4,fe, mon, iscan,' Inside scan' 
		
cOld                else if (gal(i) .lt. minq)  then
c		 beforeOrAfter = '<'
cOld		 write (aline, 540) (hklal(j,i),j=1,3), gal(i), 
cOld     &		 	angle4,fe, mon,iscan, ' BeforeScan' 
c     &		 	angle4,fe, mon,iscan, ' Before scan' 
cOld                else if (gal(i) .gt. maxq) then
c		 beforeOrAfter = '>'
cOld		 write (aline, 540) (hklal(j,i),j=1,3), gal(i), 
cOld     &		 	angle4,fe, mon,iscan, ' After scan' 
cOld     		endif

cOld      	        write (*,'(a80)') aline
cOld      	        write (unitspur,'(a80)') aline

		 write (aline, 560) (hklal(j,i),j=1,3), gal(i), 
     &		 	angle4,fe, mon, 	!!beforeOrAfter,
     &			iscan, magq, angle4_store(iscan)*180.0/acos(-1.0)
 560             format ('Al(', 3i1, ')', ' |Q|=', f5.2, 
     &			'A^-1 A4=', f6.2,  ' fix E', f4.1, 
     &			 a5,  '; pt#=', i3, ' |Q|=', f6.2, 'A^-1 A4=', f6.2)

      	        write (*,'(a80)') aline
      	        write (unitspur,'(a80)') aline
              endif	      

	      if (i .eq. nal) then
	        write (*,*) ' '
	        write (unitspur,*) ' '
	      endif
	    endif
	    
	  enddo		!! i = 1, nal
	enddo		!! iscan=1,nscan
 


      write (*,*) '&&& Al spurions at:'
      write (unitspur,*)'&&& Al spurions at:'

      if ( plotlattice) then
         call pgsfs(2)          ! outline only
         call pgsls(2)          ! dashed line
      endif 		!! if ( plotlattice) then

C  check for Q matches with while loop
         nal_match = 0
         alspur=.false.
	do i = 1, nal
	  thismatches = .false.
c	  Too hard to see diff between yellowgreen and green. Try cyangreen.
c          call pgsci(color_greencyan)
c	  thiscolor = name_color(color_greencyan)
          call pgsci(color_green)
	  thiscolor = name_color(color_green)
          do iscan=1,nscan
            magq=sqrt(qarr(1,iscan)**2+qarr(2,iscan)**2)
            if (abs(gal(i)-magq).lt.tol) then
               do j=1,3
                  h(j)=hc(j)+(iscan-(float(npts)+1.0)/2.0)*hs(j)
               enddo
	       nal_match = nal_match + 1
               alspur=.true.
	       thismatches = .true.
	       iqal_ial(nal_match) = i
	       iqal_iscan(nal_match) = iscan

c                call pgsci(color_green)
c		thiscolor = name_color(color_green)
                call pgsci(color_red)
		thiscolor = name_color(color_red)


	         write (aline, 539) (hklal(j,i),j=1,3), 
     &		 iscan, (h(jj),jj=1,3), magq, thiscolor
     		 write (*,'(a80)') aline
     		 write (unitspur,'(a80)') aline
 539             format ('Al(', 3i1,
     &		   ') spurion at pt#', 
     &		   i2, ' Q=(', f5.2, ',', f5.2, ',', f5.2, ') [rlu] ', 
     &		   '|Q|=', f5.2, ' AA^-1 ', a10)
c 539             format ('Al(', 3i1,
c     &		   ') reflection matches scan pt#', 
c     &		   i2, ' Q=(', f5.2, ',', f5.2, ',', f5.2, ') [rlu] ', a)
		 
	    endif


	    angle4 = 190.0		!! Nonsense value
            if (mon.eq.'(MON)') then 
	      asin_angle4 =
     &	      	0.5*gal(i)/sqrt(kiarr(1,iscan)**2+kiarr(2,iscan)**2)
c	           angle4 = 360/pi* 		!!  Laue scattering law
c     &        asin(0.5*gal(i)/sqrt(kiarr(1,iscan)**2+kiarr(2,iscan)**2))
            else
	      asin_angle4 =
     &	      	0.5*gal(i)/sqrt(kfarr(1,iscan)**2+kfarr(2,iscan)**2)
            endif

	    if (abs(asin_angle4  ) .le. 1.0) then
	      angle4 = 360/pi* asin(asin_angle4 )
	    endif


 541        format(a,3i1,2a,f5.2,a,$)

              if (.not. (angle4 .le. phimax*180.0/pi)) then
                call pgsci(color_darkgrey)		!! Override color
		thiscolor = name_color(color_darkgrey)
 	      endif
c	      write (*,*) 'cmd_aluminum: angle4=', angle4, ' ', 
c     &		      thiscolor, phimax*180.0/pi

	  enddo		!! i = 1, nal
          if ( plotlattice) then
            call pgcirc(0.0,0.0,gal(i))
	  endif		!! if ( plotlattice) then
	enddo		!! iscan=1,nscan

            if (alspur) then 
	       spurbit = 1
	    else
	       spurbit = 0
	    endif
            call bit(spurhist,3,spurbit,1)

         if (.not.alspur) then
            write(*,*) 'Q does not match any Al scattering.'
            write(unitspur,*) 'Q does not match any Al scattering.'
         endif

ccc         nearby=.not.nearby
         nearby= alspur

      write (*,*) ' '
      write (unitspur,*) ' '
	 
      return
      end


c==========================================================================
      subroutine cmd_orientvecs(uuv, vuv)
C  Draw the u and v orientation vectors with solid, magenta lines.
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_color.inc'

c...	passed variables
      REAL*4 uuv(3),vuv(3)

c...	local variables
      real xch, ych, size
      real temp1


        call pgsci(color_magenta)
c01dec        call pgsci(6)
        call pgsls(1)
C  It is possible that these vectors are much smaller than the rest
C  of the plot so check that the arrow head size is suitable.
        call PGQCS(4, XCH, YCH)         ! character size in world coor
        call pgqch(size)                ! character size in 1/40 of viewport
        if (4.0*xch.gt.sqrt(uuv(1)**2+uuv(2)**2)) then
	  temp1 = size*sqrt(uuv(1)**2+uuv(2)**2)/(2.0*xch)
          call pgsch(temp1)
	endif
c01Nov        if (4.0*xch.gt.sqrt(uuv(1)**2+uuv(2)**2))
c01Nov     &    call pgsch(size*sqrt(uuv(1)**2+uuv(2)**2)/(2.0*xch))
        CALL pgarro(0.0,0.0,uuv(1),uuv(2))
        CALL pgarro(0.0,0.0,vuv(1),vuv(2))
        call pgsci(color_ink)
c01Dec        call pgsci(1)
        call pgsch(size)                ! back to original size

      return
      end


c==========================================================================
      subroutine cmd_zoomspur(nscan) 
C  Zoom in/out of the spurion energy plot on the scan page
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_layout.inc'

c...	passed variables
      integer nscan


c220     write(*,'(a,a,$)') 'Enter spurion energy range for plot ', 
c     &		'(0,0 for auto):'
220         write(*,'(a,a,$)') 'Enter scan pnts and ', 
     &	 'spurion energy ranges (0,0,0,0 for auto):'
        read(*,*,err=220) xmin_sc, xmax_sc, ymin_sc_hs, ymax_sc_hs

      if (xmin_sc .gt. xmax_sc) goto 220
      if ((xmin_sc .lt. 1.0) .and. (xmin_sc .ne. 0.0) ) then
	write (*,*) 'cmd_zoomspur: Try again! ', 
     &		'Only positve scan ranges allowed.'
        goto 220
      endif
      if (xmax_sc .gt. nscan) then
	write (*,*) 'cmd_zoomspur: Try again! ', 
     &		'Largest allowed scan range is ', nscan
        goto 220
      endif

      if (ymin_sc_hs .gt. ymax_sc_hs) goto 220
      
      if (xmin_sc .eq. 0.0) then 
        auto_xmin_sc = .true.
      else
         auto_xmin_sc = .false.
      endif
      
      if (xmax_sc .eq. 0.0) then 
        auto_xmax_sc = .true.
      else
         auto_xmax_sc = .false.
      endif
      
      if (ymin_sc_hs .eq. 0.0) then 
        auto_ymin_sc_hs = .true.
      else
         auto_ymin_sc_hs = .false.
      endif
      
      if (ymax_sc_hs .eq. 0.0) then
        auto_ymax_sc_hs = .true.
      else
        auto_ymax_sc_hs = .false.
      endif
        
      return
      end

c==========================================================================
      subroutine cmd_rescaleplot(xstar, ystar, rlu, right, circsize)
C  Rescale the plot
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_layout.inc'

c...	passed variables
cCommon      real xmin, xmax, ymin, ymax
      real xstar,ystar
      logical rlu, right
      real circsize

c...	local variables
      real tmp

220   if (rlu) then
        write(*,'(a,$)') 'Enter xmin,xmax,ymin,ymax (in rlu): '
      else 
        write(*,'(a,$)') 'Enter xmin,xmax,ymin,ymax (in A^-1): '
      endif

        read(*,*,err=220) xmin,xmax,ymin,ymax
        if (rlu) then		!! convert to A^-1
           xmin=xmin*xstar
           xmax=xmax*xstar
           ymin=ymin*ystar
           ymax=ymax*ystar
        endif
        if (xmin.gt.xmax) then    ! if xmin > xmax then change orientation
          right=.not.right        ! of the spectrometer
          tmp=xmin
          xmin=xmax
          xmax=tmp
        endif
	
ccc      call square(xmin, xmax, ymin, ymax)
      circsize = (xmax- xmin)/120.0

cDebug      write (*,*) 'cmd_rescale: ', xmin, xmax, ymin, ymax
      
      return
      end

c==========================================================================
      subroutine cmd_showintensities_old(gray, rainbow)
C  Show numbers for the intensities instead of only circles.
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	passed variables
      logical gray, rainbow
      
c...	local variables
       integer temp

 20     write(*,'(a,a,$)') 
     &	  'SF^2 display: 0 for gray/colored circles, 1 for numbers: '
ccc     &	  'SF^2 display: Enter 0 for circles, 1 for numbers: '
        read (*,*, err=20) temp
        if (temp .le. 0) then
          gray= .true.
        else 
	  gray = .false.
        endif

120     write(*,'(a,$)') 
     &	  'SF^2 display: Enter 0 for gray, 1 for color: ' 
        read (*,*, err=120) temp
        if (temp.le.0) then
          rainbow= .false.
        else 
          rainbow= .true.
        endif

      return
      end

c==========================================================================
      subroutine cmd_circsize_old(c, struct)
C  Change the circlesize.  Also toggles 'grayscale' or can disable
C  the structure factor calculation completely in PLOTLATTICE      
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	common blocks
      integer natoms
      common /atoms/ natoms
      
c...	passed variables
      real c
      logical struct
      
c...	local variables

c 20     write(*,'(a,$)') 
c     &	 'Enter circle size scale factor (0 no structure calc): '
 20     continue
        write (*,*) 'Scheme for the circle size: '
        write (*,*) '  <0 for varying circles of varying gray... '
        write (*,*) '   0 for no structure factor ', 
     &		'(identical circles everywhere)... ' 
        write (*,*) '  >0 for uniform circles of varying gray...'
        write (*,*) 

220     write (*,'(a,$)') 'Enter circle size scale factor: '
        read (*,*,err=220) c
	
c        if (c.lt.0) then
c	   go to 20
c        endif

        if (abs(c) .le. 1.0e-5) then
          struct = .false.
cTR          struct=.not.struct
          c=1.0
	  write (*,*) 'Not using structure calculation'
cTR          write(*,'(2a)') 'Toggling Sf^2, sytematic absences ',
cTR     &      'calculation.'
        else
          struct = .true.
        endif

      return
      end

c==========================================================================
      subroutine cmd_sortlog_old(sortlogscale, tolerance)
C  Change the circlesize.  Also toggles 'grayscale' or can disable
C  the structure factor calculation completely in PLOTLATTICE      
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	common blocks
      integer natoms
      common /atoms/ natoms
      
c...	passed variables
      logical  sortlogscale
      real tolerance
      
c...	local variables

	  write (*,*) 'Intensities between base^n and base^(n+1) will',
     &	    ' have the same label.' 
          write (*,*) 'If a base of zero is entered, intensities ',
     &			' within tolerance will have the same label'
220       write(*,'(a,$)') 'Enter base (0 for none): '
          read(*,*,err=220) tolerance
	  if (tolerance .lt. 1.0) then
	     	  sortlogscale = .false.
		  tolerance = 100.0*tiny *float(natoms)
          else
	      sortlogscale = .true.
	  endif

      return
      end


c==========================================================================
      subroutine cmd_printout(idev )
C change the device number
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	passed variables
      integer idev

c...	local variables

 771 	 continue 
       write(*,*) 
       write(*,*) 'Device #   =0 for XSERVE (screen,',
     &        ' no printout)'
       write(*,*) '           =1 for XWIN (screen,',
     &        ' no printout)'

       write (*,*) '           =2 for B&W postscript file:',
     &  ' <base>_bw.ps   = lattice plot'
       write (*,*) '           =3 for color postscript file:', 
     &  ' <base>_cl.ps   = lattice plot'
       write (*,*) '           =4 for B&W postscript file:' ,
     &  ' <base>_scbw.ps = scan plot'
       write (*,*) '           =5 for color postscript file:' ,
     &  ' <base>_sccl.ps = scan plot'
       write (*,*) ' Apologies: PGPLOT can only process', 
     &       ' one PS "window" at a time.'
       write (*,*) ' To properly close a PS file, a new device ',
     &       'must be selected or quit program.'
c       write (*,*) ' Apology: You must switch to another', 
c     &       ' device or end the program to properly close a PS file.'
 220    write(*,'(a, $)') 'Enter device number (0-5)' 
       read(*,*,err=220) idev
        if ((idev.lt.0) .or. (idev.gt.5)) then
          write(*,*) 'Enter 0, 1, 2, 3, 4 or 5.'
          goto 771
        endif
   
   
   
c 771 	 continue
c         write(*,'(2a)') 'Enter device number: 0 for XSERVE (screen,',
c     &        ' no printout)'
c         write(*,'(2a)') '                     1 for XWIN (screen,',
c     &        ' no printout)'
c         write(*,'(2a)') '                     2 for BWPS (black',
c     &        ' and white postcript) of Lattice plot'
c         write(*,'(2a)') '                     3 for CPS (color',
c     &        ' postscript) of Lattice plot'
c         write(*,'(2a)') '                     4 for BWPS (black',
c     &        ' and white postcript) of Scan plot'
c         write(*,'(2a)') '                     5 for CPS (color',
c     &        ' postscript) of Scan plot'
cc         write(*,'(2a)') '                     2 for BWPS (black',
cc     &        ' and white postcript)'
cc         write(*,'(2a)') '                     3 for CPS (color',
cc     &        ' postscript)'
c220      write(*,'(a,$)') 'Number of selection: '
c         read(*,*,err=220) idev
c        if ((idev.lt.0) .or. (idev.gt.5)) then
c          write(*,*) 'Enter 0, 1, 2, 3, 4, or 5.'
c          goto 771
c        endif
cc02Jan        idev=idev-1

c        if ((idev.eq.2).or.(idev.eq.3)) then
c          write(*,*) 'Output to spurion.ps.  Use: lpr -Pn spurion.ps'
c          write(*,*) '(n is printer number) to printout after quitting'
c          write(*,*) 'spurion.'
c          write(*,*)
c        endif

      return
      end

c==========================================================================
      subroutine cmd_labelvec(kiarr, kfarr, kistrt, kistop,
     &       kfstrt, kfstop, nscan, over, spurhist)
C  label Q, ki, and kf
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_menu.inc'
      include 'spur_layout.inc'
      
c...    functions
      logical	bit_on

c...	passed variables
      real*4 kiarr(2,maxnpts),kfarr(2,maxnpts)
      real kistrt(2),kistop(2),kfstrt(2),kfstop(2)
      logical over
      integer nscan, spurhist
cCommon      real xmin, xmax, ymin, ymax

c...	local variables
      character*2 key
      integer spurbit, l
      real temp1, temp2

         do l=1,2
            kistrt(l)=kiarr(l,1)
            kistop(l)=kiarr(l,nscan)
            kfstrt(l)=kfarr(l,1)
            kfstop(l)=kfarr(l,nscan)
         enddo
	 
         if (bit_on(addoption,bit_labelvec)) then
           write(*,'(a)') 'Labelling Q, ki, and -kf.'
220        write (*,'(A,$)') 'Okay to overwrite on plot ((y)/n)? '
           read (*,'(A)',err=220) key
           if ((key.eq.'n') .or. (key.eq.'N') )then
             over=.false.      
	     spurbit = 0
             call bit(spurhist,6,spurbit,1)
           else if ((key.eq.'y') .or. (key.eq.'Y')) then
             over=.true.
	     spurbit = 1
             call bit(spurhist,6,spurbit,1)
           else
	     go to 220
           endif
         endif
	 
	 temp1 = xmax-xmin
	 temp2 = ymax-ymin
         call label(kistrt,kistop,kfstrt,kfstop, temp1, temp2,
     &        over)
c01Nov         call label(kistrt,kistop,kfstrt,kfstop,xmax-xmin,ymax-ymin,
c01Nov     &        over)


      return
      end

c==========================================================================
      subroutine cmd_annotation(unitout, version_stamp,
     &		arrin, nunique, uniqueLabel, 
     &      indxunique1, indxunique2, multiplicity,
     &	    pointsize, circsize, overridesize, rainbow)
C  Put/add sorting labels text on right hand side of graph
    
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_menu.inc'
      include 'spur_annote.inc'
      
c...    functions
      logical bit_on
      logical write_qbuf

c...	passed variables
       integer 	unitout
c02Mar       real 	c	!! now relcircsize
c02Mar       real 	relcircsize
       character*(80) version_stamp
c      character*80 annotation(*)
c      integer 	nannote
      real arrin(m_rlv_max,m_rlv_max)
      integer nunique
      character*1 uniqueLabel(m_rlv_max,m_rlv_max)
      integer indxunique1(unique_rlv_max)
      integer indxunique2(unique_rlv_max)
      integer multiplicity(unique_rlv_max)
c      integer indxunique1(0:unique_rlv_max)
c      integer indxunique2(0:unique_rlv_max)
c      integer multiplicity(0:unique_rlv_max)
      logical overridesize, rainbow
      real pointsize(m_rlv_max,m_rlv_max)
      real circsize



c...	local variables
      integer i, iwhere
      logical success

      len_annote = 0

c...	Annotate the q-buffer file
        iwhere = 1
        success =  write_qbuf( -1, 'Qbuffer', annotation(iwhere))
        iwhere = iwhere + 1
        success =  write_qbuf( -1, 'Lattice', annotation(iwhere))
        iwhere = iwhere + 1
        success =  write_qbuf( -1, 'Energy', annotation(iwhere))
        iwhere = iwhere + 1
        success =  write_qbuf( -1, 'Space', annotation(iwhere))
        iwhere = iwhere + 1
        success =  write_qbuf( -1, 'Orient', annotation(iwhere))
        iwhere = iwhere + 1
        success =  write_qbuf( -1, 'Temp/H', annotation(iwhere))
        iwhere = iwhere + 1
        success =  write_qbuf( -1, 'Time', annotation(iwhere))

        write (*,*) 'annotate_spurhist: Debug: loaded up qbuf. ',
     &      ' iwhere = ', iwhere
c        do i = 1, 7
c        	write (*,'(a79)') annotation(i)
c        enddo

      len_annote = len_annote + iwhere
      
      do i = 1, min(max_annote, len_menu)
c02Mar      do i = 1, min(max_annote, len_menu)
	if (len_annote .eq. max_annote) then
	  write (*,*) 'cmd_annote: Programmer! ', 
     &	  	'Need to lengthen array stdannote via max_annote'
	  goto 444
	endif

      
        if (bit_on(currentoptions, i) .and. 
     &		.not.(i.eq. bit_annote) ) then
	  call set_bit(annoteoptions, i, 1)
cDebug          write (*,*) 'cmd_annote: setting bit on:', i
	endif
	
 120	format (a58,2x, a20)
	if (bit_on(annoteoptions, i)) then
          len_annote = len_annote + 1

cMoved?	  if (i .eq. bit_circsize) then
cccc    ' ...Show structure factor intensities (circles,#)' 
c      stndannote(bit_circsize) =
c     & '...Change circ size (<0 vary, >0 uni, 0 no struct)'
cMoved?	    if (abs(c) .lt. 1.0e-5) then 
cMoved?              write (annotation(len_annote), 120)
cMoved?     &	     'No structure factors shown', 
cMoved?     &			  spuritemvalue(i)
cMoved?	    else if (c .lt. 0.0) then 
cMoved?              write (annotation(len_annote),120)
cMoved?     &	     'Intensity proportional to circle size and grey tone', 
cMoved?     &			  spuritemvalue(i)
cMoved?	    else	!! if (c .gt. 0.0) then 
cMoved?              write (annotation(len_annote),'(a56,4x, a20)')
cMoved?     &	     'Intensity proportional to grey tone of uniform circles', 
cMoved?     &			  spuritemvalue(i)
cMoved?            endif
cMoved?          else
            write (annotation(len_annote),'(a56,4x, a20)') 
     &			  stndannote(i), spuritemvalue(i)
cMoved?          endif

cDebug           write (*,*) 'cmd_annote: bit on:', i, 
cDebug      &		  'writing: len_annote=', len_annote
cDebug 	  write (*,*) '   ', annotation(len_annote) 
	endif 

c	if ( (stndannote(i) .ne. nullstring) .and. 
c     &		 	 bit_on(currentoptions, i) ) then
c          write (annotation(i),'(a56,4x, a20)') stndannote(i),
c     &			  spuritemvalue(i)
c        endif

      enddo

 444  continue
 
 

      if (unitout .gt. 0) then
       write (unitout,*) 'cmd_annote>>' 
        do i = 1, len_annote
          write (unitout,'(a80)') annotation(i)
        enddo
        write (unitout,*) 

      else if (unitout .eq. 0) then
cDebug       write (*,*) 'cmd_annote: len_annote=',len_annote
cDebug     &		, annoteoptions, currentoptions
        do i = 1, len_annote
          write (*,'(a80)') annotation(i)
c          write (*,'(a, i2, a60)') 'cmd_annote:', i, annotation(i)
        enddo
        write (*,*) 

      else if (unitout .lt. 0) then

      endif

cDebug      write (*,*) 'cmd_annotation:1 i, j', indxunique1(1),indxunique2(1)
      call plot_txtinfo(annotation, len_annote)
cMoved      call plot_txtinfo(annotation, len_annote, version_stamp, .true.)
cDebug      write (*,*) 'cmd_annotation:1 i, j', indxunique1(1),indxunique2(1)
      call plot_labelinfo(arrin, nunique, uniqueLabel, 
     &      indxunique1, indxunique2, multiplicity,
     &	    pointsize, circsize, overridesize, rainbow)

      return
      end


c==========================================================================
      subroutine cmd_modify()
C modify Q, E, and npts
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_modmenu.inc'
      include 'spur_buf.inc'

c...	functions
      logical bit_on

c...	passed variables
c02Mar      real ec, es
c02Mar      REAL*4 hc(3), hs(3)
c02Mar      real fe
c02Mar      integer npts
c02Mar      CHARACTER*5 mon

c...	local variables
      character*5 key
      integer i
      character*72 aline
      integer newbit

      aline = 'Modify buffer parameters.'
c     &'Modify parameters. Remember: these are not written to *.bufout!'

      current_modop = 0
      add_modop = 0

      do while (.not. bit_on(add_modop,bit_modquit)) 
        call load_mod_itemvalue()
c02Mar        call load_mod_itemvalue(ec, es, hc, hs, fe, npts, mon)
        call display_submenu(aline, key, current_modop, add_modop,
     &	      len_modmenu, modmenu, moditemvalue, modmenuitem, .true.)
c        call display_modmenu(aline, key, current_modop, add_modop,
c     &	      len_modmenu, modmenu, moditemvalue, modmenuitem)

	 call which_bit(add_modop,newbit,1)
 220	 if ((newbit .eq. 0) .or. (newbit .gt. len_modmenu)) then
	   write (*,*) 'cmd_modify: newbit is out of range ', newbit
	   goto 555
	 else
	   write (*,'(a,a)') modmenu(newbit), ' Enter value(s):'
c	   write (*,'(a,a,$)') modmenu(newbit), ' Enter value(s):'
	 endif

c...	Lattice:
         if (bit_on(add_modop,bit_lp) ) then
	   read(*,*,err=220) (qb_lattice_a(i),i=1,3)

         else if (bit_on(add_modop,bit_la) ) then
	   read(*,*,err=220) (qb_lattice_angle(i),i=1,3)

c...	Energy
         else if (bit_on(add_modop,bit_ec) ) then
	   read(*,*,err=220) qb_energy_ec

         else if (bit_on(add_modop,bit_es) ) then
	   read(*,*,err=220) qb_energy_es

         else if (bit_on(add_modop,bit_fe) ) then
	   read(*,*,err=220) qb_energy_fe

         else if (bit_on(add_modop,bit_mo) ) then
ccc	   read(*,*,err=220) qb_energy_fix;	
c	!! Toggle instead
	  if (qb_energy_fix .eq. '(MON)') then
	     qb_energy_fix  = '(ANA)'
	  else
	     qb_energy_fix  = '(MON)'
	  endif

         else if (bit_on(add_modop,bit_dm) ) then
	   read(*,*,err=220) qb_energy_dm

         else if (bit_on(add_modop,bit_da) ) then
	   read(*,*,err=220) qb_energy_da

c...	Space:
         else if (bit_on(add_modop,bit_qc) ) then
	   read(*,*,err=220) (qb_space_hc(i),i=1,3)

         else if (bit_on(add_modop,bit_qs) ) then
	   read(*,*,err=220) (qb_space_hs(i),i=1,3)

c...	Orient:
         else if (bit_on(add_modop,bit_uu) ) then
	   read(*,*,err=220) (qb_orient_u(i),i=1,3)

         else if (bit_on(add_modop,bit_vv) ) then
	   read(*,*,err=220) (qb_orient_v(i),i=1,3)

c...	Temp:
c         else if (bit_on(add_modop,bit_tc) ) then
c	   read(*,*,err=220) qb_temp_tc

         else if (bit_on(add_modop,bit_t0) ) then
	   read(*,*,err=220) qb_temp_t0
	   if (qb_temp_t0 .lt. 0) goto 220

         else if (bit_on(add_modop,bit_ts) ) then
	   read(*,*,err=220) qb_temp_incT

c...	Time:
         else if (bit_on(add_modop,bit_np) ) then
	   read(*,*,err=220) qb_time_scanpts
	   if ((qb_time_scanpts .lt. 0) .or. 
     &	       (qb_time_scanpts .gt. maxnpts)) goto 220

         else if (bit_on(add_modop,bit_modhelp) ) then 
	   call help('', 'spurion', 
     &	   'spurion MainMenu mm_ModifyMenu')
         else if (bit_on(add_modop,bit_modquit) ) then 
	 else
	   write (*,*) 'cmd_modify: Unrecognized command: "', key, '"'
         endif

c...	old menu style
c         if (bit_on(add_modop,bit_ec) ) then 
c220         write(*,'(a,$)') 'New E-center: '
c            read(*,*,err=220) ec
cc         endif

c         else if (bit_on(add_modop,bit_es) ) then 
c222         write(*,'(a,$)') 'New E-step: '
c            read(*,*, err=222) es
cc         endif

c         else if (bit_on(add_modop,bit_qc) ) then 
c224         write(*,'(a,$)') 'New Q-center: '
c            read(*,*, err=224) (hc(i),i=1,3)
cc         endif

c         else if (bit_on(add_modop,bit_qs) ) then 
c226         write(*,'(a,$)') 'New Q-step: '
c            read(*,*,err=226) (hs(i),i=1,3)
cc         endif

c         else if (bit_on(add_modop,bit_mo) ) then 
c            if (mon.eq.'(MON)') then
c               mon='(ANA)'
c            elsec
c               mon='(MON)'
c            endif
cc         endif

c         else if (bit_on(add_modop,bit_fe) ) then 
c228        write(*,'(a,$)') 'Fixed energy?  '
c           read(*,*,err=228) fe
cc         endif

c         else if (bit_on(add_modop,bit_np) ) then 
c230         write(*,'(a,$)') 'How many points? '
c            read(*,*,err=230) npts
cc         endif

c         else if (bit_on(add_modop,bit_modhelp) ) then 
cc	   call help('', 'spurion', 'spurion')
cc         endif

c         else if (bit_on(add_modop,bit_modquit) ) then 
cc	   do i = 1, len_modmenu
cc             write(unitspur,*) modmenuitem(i), ' : ', modmenu(i),
cc     &     '   ', moditemvalue(i)
cc	   enddo

c	 else
c	   write (*,*) 'cmd_modify: Unrecognized command: "', key, '"'
c         endif

 555    continue
      enddo
      
      return
      end


c==========================================================================
      subroutine cmd_selecttolerance()
C set the AA^-1 tolerence for spurion matches
    
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	common blocks
      real tol, etol		!! spatial and energy tolerances       
      common /spur/ tol, etol

c...	passed variables
c...	local variables

 480     write(*,'(a,f5.3,a,$)') 'Current rlv tolerance ',tol,
     &        ' AA^-1.  Enter new value. '
         read(*,*, err=480) tol
         if (.not.(tol .gt. 0.0)) then
            write(*,*) 'Must be a positive value.'
            go to 480
         endif
         write(unitspur,*) 'New spatial tolerance is ',tol,' AA^-1.'

 580     write(*,'(a,f5.3,a,$)') 'Current energy  tolerance ', etol,
     &        ' meV?.  Enter new value. '
         read(*,*, err=580) etol
         if (.not.(etol .gt. 0.0)) then
            write(*,*) 'Must be a positive value.'
            go to 580
         endif
         write(unitspur,*) 'New energy tolerance is ', etol,' meV?.'

      return
      end

c==========================================================================
      subroutine cmd_BumpA4( phimax, pi)
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	passed variables
      real phimax, pi

c...	local variables

 580     write(*,'(a,f5.1,a,$)') 'Current max A4 ',phimax*180.0/pi,
     &        ' degrees.  Enter new value. '
         read(*,*,err=580) phimax
         if (.not.(phimax.gt.0).or.(phimax.gt.180)) then
            write(*,*) 'Must be between 0 and 180.'
            go to 580
         endif
         write(unitspur,*) 'New A4 max is ',phimax,' degrees.'
c         write(8,*) 'New A4 max is ',phimax,' degrees.'
         phimax=pi*phimax/180.0

      return
      end


c==========================================================================
      subroutine cmd_genICPrun(mon, npts, bmatuv, h, q, hc, hs, ki, 
     &      ec, es, en, ei, ef, fe, da, dm, a2, a3, a4, a6, pi)
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	passed variables
      CHARACTER*5 mon
      integer npts
      real bmatuv(3,3),h(3),q(3)
      REAL*4 hc(3), hs(3)
      real ki(2)
      real ec, es
      real en, ei, ef, fe, da, dm
      real a2, a3, a4, a6
      real pi

c...	local variables
      integer k, j
      logical fail

       OPEN(UNIT=uniticp,FILE=file_icp,STATUS='UNKNOWN', err=660)

        write(*,*) 'Pt#  E(meV)   Qh    Qk    Ql     A2       ',
     &        'A3      A4       A6'
         write(*,*) '---  ------ ----- ----- ----- --------',
     &        ' -------- -------- --------'
        write(uniticp,*,err=666) 'Pt#  E(meV)   Qh    Qk    Ql',
     &        '     A2       A3      A4       A6'
         write(uniticp,*,err=666) '---  ------ ----- ----- ----- ',
     &        '-------- -------- -------- --------'
 484     format(i3,2x,f6.2,1x,f5.2,1x,f5.2,1x,f5.2,1x,f8.3,1x,f8.3,
     &        1x,f8.3,1x,f8.3)
c 484     format(1x,i2,2x,f6.2,1x,f5.2,1x,f5.2,1x,f5.2,1x,f8.3,1x,f8.3,
c     &        1x,f8.3,1x,f8.3)
         do k=1,npts
            en=ec+(k-(float(npts)+1.0)/2.0)*es
            do j=1,3
               h(j)=hc(j)+(k-(float(npts)+1.0)/2.0)*hs(j)
            enddo

            if (mon.eq.'(MON)') then
               ei=fe
               ef=ei-en
            else
               ef=fe
               ei=ef+en
            endif

            a2=180.0/pi*2.0*asin(pi*sqrt(k2_to_E/ei)/dm)
            a6=180.0/pi*2.0*asin(pi*sqrt(k2_to_E/ef)/da)
c02Apr            a2=180.0/pi*2.0*asin(pi*sqrt(2.072/ei)/dm)
c02Apr            a6=180.0/pi*2.0*asin(pi*sqrt(2.072/ef)/da)
            call calcki(bmatuv,h,ei,ef,q,ki,fail,a4,a3,k)
            a3=180/pi*a3
            a4=180/pi*a4
            write(*,484) k,en,(h(j),j=1,3),a2,a3,a4,a6
            write(uniticp,484,err=666) k,en,(h(j),j=1,3),a2,a3,a4,a6
         enddo

      go to 999
 660  write (*,*) 'cmd_genICPrun: Error opening file:', file_icp

      go to 999
 666  write (*,*) 'cmd_genICPrun: Error writing to file:', file_icp
 
 999  CLOSE(uniticp)

      return
      end

c==========================================================================
      subroutine cmd_selectReflection( sp, spurtitle, m_rlv, 
     &      hh, kk, ll, zeta, eta, absent, sfhkl, selct_sfhkl, 
     &	    qlow, qhigh)
C  Sf^2 printout
     
      implicit none
      include 'arraysize.inc'
      include 'io_files.inc'

c...	functions
c02May      logical inqrange
      
c...	passed variables
      integer sp(3)
      character*10 spurtitle
      integer m_rlv
      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*4 zeta(m_rlv_max,m_rlv_max), eta(m_rlv_max,m_rlv_max)
      logical absent(m_rlv_max,m_rlv_max)
      real sfhkl(m_rlv_max,m_rlv_max), selct_sfhkl
      real qlow, qhigh

c...	local variables
      integer len, i, j
      logical firsttime

      firsttime = .true.
      
100      write(*,'(a,$)') 'h k l of desired reflection? '
C just use sp variable for double duty
         read(*,*,err=100) (sp(i),i=1,3)
         call title(sp,spurtitle,len)
C  try to match with one of the calculated points
         do i=1,m_rlv
            do j=1,m_rlv
               if ((hh(i,j).eq.sp(1)).and.(kk(i,j).eq.sp(2)).and.
     &              (ll(i,j).eq.sp(3))) goto 421
            enddo
         enddo

	 if (firsttime) then
           write (*,*) 'Possible reflections:'
           do i=1,m_rlv
             do j=1,m_rlv
c02May		if (inqrange(qlow, qhigh, zeta(i,j), eta(i,j))) then
		  if (absent(i,j)) then
		    write (*, 110) hh(i,j), kk(i,j), ll(i,j),
     &			     sfhkl(i,j), ' in range'
		  else
		    write (*, 112) hh(i,j), kk(i,j), ll(i,j), ' absent'
		  endif
cDebug		else
cDebug		  write (*, 110) hh(i,j), kk(i,j), ll(i,j), ' out range'
c02May		endif
              enddo
           enddo
110       format ('Sf^2(', I3, ',', I3, ',', I3, ')=', g12.5, a )
112       format ('Sf^2(', I3, ',', I3, ',', I3, ')=', 12x, a )
           write (*,*) ''
	   firsttime = .false.
	   goto 100
	 else
	   selct_sfhkl = -666.6
           write(*,'(3a)') 'Sorry, ',spurtitle(1:len),
     &              ' reflection does not match those generated.'
           goto 999
	 endif
	 
 421     if (absent(i,j)) then
            write(*,'(2a)') spurtitle(1:len),
     &              ' reflection is systematically absent.'
            write(unitstruct,'(2a)') spurtitle(1:len),
     &              ' reflection is systematically absent.'
c            write(7,'(2a)') spurtitle(1:len),
c     &              ' reflection is systematically absent.'
	     selct_sfhkl = 0.0
         else
            write(*,'(3a,f8.3)') 'Sf^2',spurtitle(1:len),' = ',
     &                sfhkl(i,j)
            write(unitstruct,'(3a,f8.3)') 'Sf^2',spurtitle(1:len),' = ',
     &                sfhkl(i,j)
c            write(7,'(3a,f8.3)') 'Sf^2',spurtitle(1:len),' = ',
c     &                sfhkl(i,j)
	     selct_sfhkl = sfhkl(i,j)
         endif
	 
999   return
      end



c==========================================================================
c      logical function inqrange(qlow, qhigh, zzeta, eeta)
cc...	Removed May 2002
c      implicit none
cc      include 'arraysize.inc'
c      
cc..	passed variables      
c      real qlow, qhigh
c      real zzeta
c      real eeta
c     
cc...	local variables
c      real qq
c      
c      qq = zzeta*zzeta + eeta*eeta 
c      if ((qq .ge. qlow) .and. (qq. le. qhigh)) then
c        inqrange = .true.
c      else
c        inqrange = .false.
c      endif
c      
c      return
c      end
