      PROGRAM SPURION
C***********************************************************
C* SPURION CALCULATES THE LOCATION OF SPURIOUS NEUTRON 
C* SCATTERING PEAKS
C* WRITTEN BY CARL ADAMS, SEPT. 16, 1998
C***********************************************************

      implicit none
cDec01      IMPLICIT REAL (A-H,O-Z)
      include 'arraysize.inc'
      include 'io_files.inc'
      include 'spur_color.inc'
      include 'spur_menu.inc'		
      include 'spur_structmenu.inc'		
      include 'spur_layout.inc'
      include 'spur_scan.inc'
      include 'spur_buf.inc'
      include 'spur_warn.inc'

c...	data statements
cWrong      include 'spur_color_data.inc'

c...	common blocks
      logical noisy
      common /noisy/ noisy
 
      real branchdeg
      common /specparm/ branchdeg
      
      real tol, etol		!! spatial and energy tolerances       
      common /spur/ tol, etol

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

c...	functions
      logical bit_on
      logical read_qbuf, verify_qbuf, write_qbuf
      logical qbuf_consistent
      
c...	local variables
      character*80 version_stamp
      character*20 date_stamp

      real pi
      real ec, es
      integer npts
      real angle, en, ei, ef, fe
      integer nscan, nspurtot
      real a4, a3, a2, a6
      real qmin
      
c...	new variables
      integer m_rlv		! INteded as overrise of m_rlv_max
            
c...	DECLARED VARAIBLES

C  The number of aluminum peaks, the m x m size of the reciprocal 
C  lattice, and the number of n*useDescending-n*kf spurions to calculate.

c...	h(3) is the current h,k,l for point i of npts
c            h(i)=hc(i)+(k-(float(npts)+1.0)/2.0)*hs(i)
c	hc(3) is the center h,k,l of scan 
c	hs(3) is the step h,k,l of scan 

      REAL*4 a(3),aa(3),h(3),hc(3),hs(3),bmat(3,3),bmatuv(3,3),b(3)
      REAL*4 zeta(m_rlv_max,m_rlv_max), eta(m_rlv_max,m_rlv_max)
      real*4 q(3),ki(2),kf(2)
      real sfhkl(m_rlv_max,m_rlv_max)
      logical absent(m_rlv_max,m_rlv_max)
      real*4 qarr(2,maxnpts),kiarr(2,maxnpts),kfarr(2,maxnpts)
      real*4 incohmon(2),incohana(2)
      real*4 incmonarr(2,maxnpts),incanaarr(2,maxnpts)

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
      logical rlu,constq
      integer axis
      REAL*4 uu(3),vv(3),uxv(3),uuv(3),vuv(3)
      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 kistrt(2),kistop(2),kfstrt(2),kfstop(2)
cInc      real gal(nal)
      real minq,maxq,magq
      real qlowmin,qlowmax
cInc      integer hklal(3,nal)
cInc      logical alspur
      logical nearby
      INTEGER u(3),v(3),uxvi(3)
      integer hh(m_rlv_max,m_rlv_max)
      integer kk(m_rlv_max,m_rlv_max)
      integer ll(m_rlv_max,m_rlv_max)
      integer hklcen(3)
      integer sp(3)
      integer selct_sp(3)
cInc      integer ispur(n_spur_max*n_spur_max)
cInc      integer jspur(n_spur_max*n_spur_max)
cInc      integer ielst(n_spur_max*n_spur_max)
cInc      integer jelst(n_spur_max*n_spur_max)
cInc      integer hspur(n_spur_max*n_spur_max)
cInc      integer kspur(n_spur_max*n_spur_max)
cInc      integer lspur(n_spur_max*n_spur_max)
      integer nspur
      integer idev
      real phimax,qmax
      CHARACTER*5 mon,scale
      CHARACTER*41 xtitle,ytitle,maintitle
      CHARACTER*41 x2title,y2title
      character*10 zetatitle,etatitle,qztitle,spurtitle
      logical skip,right,gray,struct,fail,incohflag,prev
      logical rainbow
c02Mar      real circsize,c
      real circsize, relcircsize
      real ccircsize		!! Added TR, Nov 2001
      
      real dm,da
C this is an 8-bit binary number...
      integer spurhist,spurbit
      logical replot
      logical over
      logical overridesize
      real pointsize(m_rlv_max,m_rlv_max)

       integer indx1(unique_rlv_max), indx2(unique_rlv_max)
       integer indxunique1(unique_rlv_max)
       integer indxunique2(unique_rlv_max)
       integer multiplicity(unique_rlv_max)
       real	limlo(unique_rlv_max)
       real	limhi(unique_rlv_max)
       integer  nunique
       character*1 uniqueLabel(m_rlv_max,m_rlv_max)
       integer uniqueRank(m_rlv_max,m_rlv_max)
cChanged       real	tolerance
       real	sorttol, base
       logical	sortlogscale, plotunique, useDescending
       real	selct_sfhkl
       integer nincohmon, nincohana
cInc       integer nal_match
       integer  nelastictot
       integer	selected_bit
       logical hide_absences
       integer i, j
       character*1 whichqbuf, whichlpi
       logical showcalls
       parameter (showcalls = .false.) 
c       parameter (showcalls = .true.) 
       logical just_erase
       character*20 junkline

C Code for spurhist: (obsolete)
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
cccC     not used               2^7

      logical success, freshplot, reset_freshplot
cNope      real the_magq, temp_magq
cNope      integer iscan
      

c...	Initialization of variables
      data m_rlv /m_rlv_max/
      
      junkline = ' '
      
C  Display program introduction
      if (showcalls) write (*,*) 'main: call intro'
      call intro(version_stamp, date_stamp)

      if (showcalls) write (*,*) 'main: call load_menuitems'
      call load_menuitems()
      if (showcalls) write (*,*) 'main: call load_menu'
      call load_menu()
      if (showcalls) write (*,*) 'main: call load_stndannote' 
      call load_stndannote()

      if (showcalls) write (*,*) 'main: call load_modmenuitems'
      call load_modmenuitems()
      if (showcalls) write (*,*) 'main: call load_modmenu'
      call load_modmenu()

      if (showcalls) write (*,*) 'main: call load_structmenuitems'
      call load_structmenuitems()
      if (showcalls) write (*,*) 'main: call load_structmenu'
      call load_structmenu()

      if (showcalls) write (*,*) 'main: call load_constants'
      call load_constants()
      
c... Get filename base (no extension) for first time (get default values)
 22    continue
      if (showcalls) write (*,*) 'main: call cmd_filenames' 
      call cmd_filenames(version_stamp, struct,
     &		 whichqbuf, whichlpi, .true.)

c03Mar 22      call get_filenames(version_stamp, struct)
c	!! If can't opne LPI file struct=.false. and no intensities have been
c	!! calculated

c	whichqbuf = 'd' from cmd_filenames
      if (whichqbuf .eq. 'd' ) then
        if (showcalls) write (*,*) 'main: call qbuf_default'
        call qbuf_default()      ! default values in case no files read in.
      endif
      file_base = 'SPUR_DEFAULT'
c      file_bufin = file_base(1:ilen(file_base))//ext_bufin
c      file_lpi = file_base(1:ilen(file_base))//ext_lpi
cdebug      success = write_qbuf(0, 'All', junkline)	! to screen
      success = write_qbuf(unitdebug, 'All', junkline)	! debug
cDebug      write (*,*) 'main: Debug: Wrote out default qbuffer OK = ',success

      currentoptions = 0
      addoption = 0
      warn_options = 0
      warnstring = ' '
      current_structop = 0
      add_structop = 0
      noisy = .true.	     ! Extra debuggin info and IO help.
      freshplot = .true.
      just_erase = .false.
      hide_absences = .true.	!! Don't plot location of systematic absences
      scale=' '
      prev=.false.           ! Sfq^2 has not been calculated yet
      rainbow=.false.        ! indicate structure factor by gray or color
      aspect_ratio = ps_height/ ps_width ! set graphics window to 8.5 X 11
      write (*,*) 'main: aspect_ratio =', aspect_ratio
      
      struct = .true.		!! Plot structure factors
            
	nal_match = 0
        nincohmon = 0
	nincohana = 0
	nspurtot = 0
	nelastictot = 0
	
	draw_al_iscan  = 1
           
 24   continue

      do while (.not. bit_on(addoption,bit_quit)) 
      
      if (prev) then
        write (unitdebug,*) 'currentoptions =', currentoptions, 
     &		' addoption=', addoption, freshplot, prev
        if (showcalls) write (*,*) 'main: call load_itemvalue'
c02Mar        call load_itemvalue(tol, etol, c, struct,
        call load_itemvalue(tol, etol, 
cMoved     &		relcircsize, struct,
cMoved     &	    gray, rainbow, plotunique, sortlogscale, tolerance, 
     &      phimax, rlu, right, selct_sp, selct_sfhkl, over,
     &	    m_rlv, nabsent, npresent, nspurtot, nelastictot, 
     &	    nal_match, nincohmon, nincohana, 
     &	    grid_style, grid_size1, grid_size2,hide_absences)
        write (unitdebug, *) 'calling display_menu. Prev scale=', scale
        if (showcalls) write (*,*) 'main: call display_menu'
        call display_menu(version_stamp, scale, spurhist)
        if (showcalls) write (*,*) 'main: call which_bit'
	call which_bit(addoption, selected_bit, 1)
c	Need to close up the graphics windows first before quiting.	
	if (bit_on(addoption,bit_quit)) go to 545

c...	Request to toggle off the option, so need a fresh plot
c...	Works for several options (see routine display_menu)
        if (selected_bit .ge. 1) then
         if (.not. bit_on(currentoptions, selected_bit)) then
           freshplot = .true.
c	  !! turn off this bit, so that the bits match for selected_bit	 
           if (showcalls) write (*,*) 'main: call set_bit'
	   call set_bit(addoption, selected_bit, 0)
           write (*,*) 'Toggle off ', spurmenuitem(selected_bit), ' : ', 
     &	   	spurmenu(selected_bit), ' ', 
     &		bit_on(currentoptions, selected_bit), 
     &      	bit_on(addoption, selected_bit),' ', 
     &		spuritemvalue(selected_bit)
        endif
       else
          go to 555		!! skip unrecognized commands
       endif
cDebug           write (*,*) 'Selected: ', spurmenuitem(selected_bit), ' : ', 
cDebug     &	   	spurmenu(selected_bit), ' ', 
cDebug     &		bit_on(currentoptions, selected_bit), 
cDebug     &      	bit_on(addoption, selected_bit),' ', 
cDebug     &		spuritemvalue(selected_bit)
cDebug           write (*,*) 'main: freshplot =', freshplot, ' prev=', prev

      endif		!!  if (prev) then

c----------------------------------------------------------------------
c 	commands which require that the plot be started over freshly.
c----------------------------------------------------------------------

C  zoom in/out of spurion energy plot on scan's page
      if (bit_on(addoption,bit_zoomspur) ) then 
	call cmd_zoomspur(nscan)
	freshplot = .true.
      endif

C  Rescale the plot
      if (bit_on(addoption,bit_rescale) ) then 
	call cmd_rescaleplot(xstar, ystar, rlu, right, circsize)
	freshplot = .true.
cDebug	write (*,*) 'after cmd_rescale ', xmin, xmax, ymin, ymax
      endif

C  Change the style of the display of the Bragg intensities
      if (bit_on(addoption,bit_structmenu)) then 
        call cmd_structmenu(relcircsize, sorttol, base, 
     &	      struct, gray, rainbow, plotunique, sortlogscale,
     &		hide_absences)
 	freshplot = .true.
      endif

C  Show numbers for the intensities instead of only circles.
cMoved      if (bit_on(addoption,bit_showintensities)) then 
cMoved 	call cmd_showintensities(gray, rainbow)
cMoved 	freshplot = .true.
cMoved       endif

C  Show numbers for the intensities instead of only circles.
cMoved       if (bit_on(addoption,bit_absent)) then 
cMoved 	hide_absences = .not. hide_absences 
cMoved 	freshplot = .true.
cMoved       endif

C  Change the circlesize for the structure factors. IF c=0, disables
C  the structure factor calculation completely in PLOTLATTICE      
cMoved      if (bit_on(addoption,bit_circsize)) then 
cMoved	call cmd_circsize(relcircsize, struct)
c02MarcMoved	call cmd_circsize(c, struct)
cMoved	freshplot = .true.
cMoved      endif
      
C  Sort the structure factors and then bunch by multiplicity 
c  if tolerance = 0.
c If tolerance = base, bunch by base^n (log scale).
cMoved      if (bit_on(addoption,bit_sortlog)) then 
cMoved	call cmd_sortlog(sortlogscale, tolerance)
cMoved	freshplot = .true.
cMoved      endif

      
      
C change the device number
      if (bit_on(addoption,bit_device)) then 
        call cmd_printout(idev)
c        write(*,'(a,i3)') 'current spurhist variable: ',spurhist
	just_erase = .false.
	freshplot = .true.
      endif

c...	Been removed from menu and spur_menu.inc
cC toggle the units between rlu and AA^-1
c      if (bit_on(addoption,bit_units)) then 
c         rlu=.not.rlu
c	 freshplot = .true.
c      endif

C modify Q, E, and npts
      if (bit_on(addoption,bit_mod)) then 
         warn_options = 0
	 warnstring = ' '
	 if (.not. qbuf_consistent()) then 
	   write (*,*) 'main: Something is inconsistent in the q-buffer!'
	 endif
         call cmd_modify()
c02mar         call cmd_modify(ec, es, hc, hs, fe, npts, mon)
	 if (.not. qbuf_consistent()) then 
	   write (*,*) 'main: Something is inconsistent in the q-buffer!'
	 endif
	 freshplot = .true.

c...	Eventually only use common block variables qb_*.
c...	For now copy the contents.
      do i = 1, 3
        a(i) =  qb_lattice_a(i)
        aa(i) =  qb_lattice_angle(i)
	hc(i) = qb_space_hc(i)
	hs(i) = qb_space_hs(i)
	u(i) = qb_orient_u(i)
	v(i) = qb_orient_v(i)
      enddo 
      ec = qb_energy_ec
      es = qb_energy_es
      fe = qb_energy_fe
      mon = qb_energy_fix
      dm = qb_energy_dM
      da = qb_energy_dA	
      npts = qb_time_scanpts	

         if (.not. verify_qbuf(a,aa,hc,hs,fe,ec,es,u,v,npts,mon,dm,da, 
     &		 .true.)) then
           write (*,*) 'main: Continuing anyway. #4'
	   stop
         endif

c...	Eventually only use common block variables qb_*.
c...	For now copy the contents. 

c      do i = 1, 3
cc        qb_lattice_a(i) = a(i)
cc        qb_lattice_angle(i) = aa(i)    
c	qb_space_hc(i) = hc(i)
c	qb_space_hs(i) = hs(i)
cc	qb_orient_u(i) = u(i)
cc	qb_orient_v(i) = v(i)
c      enddo 
c      qb_energy_ec = ec
c      qb_energy_es = es
c      qb_energy_fe = fe
c      qb_energy_fix = mon
cc      qb_energy_dM = dm
cc      qb_energy_dA = da	
c      qb_time_scanpts = npts   	
c         if (.not. verify_qbuf(a,aa,hc,hs,fe,ec,es,u,v,npts,mon,dm,da, 
c     &		 .true.)) then
c           write (*,*) 'main: Continuing anyway. #3'
c	   stop
c         endif

c	 !! energy resolution tol for spurion hitting the bragg peak
c	!! Using the step size of the scan is too small, practically
c	!! Not needed anyway, because the zero crossings of the 
c	!! elastic spurions are calculated exactly.
cNope      if (abs(qb_energy_es) .gt. tiny) then
cNope         etol= abs(qb_energy_es)
cNope       else
cNope        etol= 1.0 
cNope       endif

c...	The q tolerance is harder to figure out than the e-tol
c...	 from the scan step size.
c...	The step size is in terms of hkl, but the tolerance is used in
c...	the scattering plane, after bmat has bee used to transform it.
c...	Take out this attemp at cleverness.
cNope        the_magq = 0.0
cNope        do iscan = 1, nscan-1
cNope          temp_magq = sqrt((qarr(1,iscan) - qarr(1,iscan+1))**2 
cNope     &	  	+ (qarr(2,iscan) - qarr(2,iscan+1))**2)	!! in AA^-1
cNope	  if (temp_magq .gt. the_magq)  the_magq = temp_magq 
cNope        enddo
cNope        if (the_magq  .gt. tiny) then
cNope          tol= the_magq
cNope        else
cNope          tol= 0.05 		!! probably qb_space_hs(3) = 0 or nscan=1
cNope        endif

      endif

C set the AA^-1 tolerence for spurion matches
      if (bit_on(addoption,bit_tol)) then 
         call cmd_selecttolerance()
	 freshplot = .true.
      endif

C set the maximum A4 (BUMP, sorry, running out of letters)
      if (bit_on(addoption,bit_bumpA4)) then 
         call cmd_BumpA4( phimax, pi)
	 freshplot = .true.
      endif

C  Toggle Show the sorting of the intensities
cMoved      if (bit_on(addoption,bit_showsort)) then 
cMoved        plotunique = .not. plotunique      
cMovedc	call cmd_showsort(plotunique)
cMoved	freshplot = .true.
cMoved      endif

C  d: start picture over with default values or read in different files.
      if (bit_on(addoption,bit_default)) then 
c02Apr30        call pg_close(idev)  
cccc	just_erase = .true.
        call cmd_filenames(version_stamp, struct,
     &		 whichqbuf, whichlpi, .false.)
	 if (.not. qbuf_consistent()) then 
	   write (*,*) 'main: Something is inconsistent in the q-buffer!'
	 endif
	freshplot = .true.
	currentoptions = 0
	addoption = 0
	warn_options = 0
        warnstring = ' '
cc...	When open a new file, keep the old plotting options.
cNo!        current_structop = 0
        add_structop = 0
	prev = .false.
      endif


c------------------------------------------------------------
c 	prep for a fresh plot graphics
c------------------------------------------------------------

c     !! Only load default values at beginning or when requested by user. 
      if (currentoptions .eq. 0) then
        write (*,*) 'Main: (currentoptions .eq. 0). Start from scratch',
     &		' rainbow = false'
        skip=.false.           ! a leftover system variable
        gray=.true.            ! indicate structure factor by gray scale
        rainbow=.false.        ! indicate structure factor by gray or color
        struct=.true.          ! read in structural data from LPI     
cDebug        idev=2                 ! B&W PS device
        idev=0                 ! xserve device
        rlu=.true.             ! plot labels given in rlu 
        constq=.true.          ! constq or constE
        nearby=.true.          ! only plot Al Ewald spheres near Q
        prev=.false.           ! Sfq^2 has not been calculated yet
        replot=.false.         ! no plot has been made yet
        plotunique = .true.    ! Spots have been sorted by intensity
      			     ! and can be labeled in plot
        useDescending = .true.  
        auto_xmin_sc    = .true.
        auto_xmax_sc    = .true.
        auto_ymin_sc_hs = .true.
        auto_ymax_sc_hs = .true.

        spurhist=0
        xmin=-0.5              ! default plot limits
        xmax= 0.5
        ymin=-0.5
        ymax= 0.5
	
	do i = 1,n_spur_max
          do j = 1,n_spur_max
	     afloat_hkl_start(i,j)  = ''
	     afloat_hkl_end(i,j)  = ''
	  enddo
	enddo

        pi=4.0*atan(1.0)       ! useful for degrees/radians conversion
        phimax=115.0*pi/180    ! max value of a4

        branchdeg=-120.0
        tol= .05  !! rlv resolution tol for spurion hitting the bragg peak
        etol= 1.0 !! energy resolution tol for spurion hitting the bragg peak
	
	nal_match = 0
        nincohmon = 0
	nincohana = 0
	nspurtot = 0
	nelastictot = 0

        selct_sfhkl = 0.0
C  right=.false. for a left-handed spectrometer
        right=.true.
        freshplot = .true.

      if (whichqbuf .eq. 'd' ) then
        if (showcalls) write (*,*) 'main: call qbuf_default'
        call qbuf_default()      ! default values in case no files read in.
        success = .true.
      else if (whichqbuf .eq. 'c' ) then	!! keep current values
        success = .true.
      else 		!! read either from q-buffer file or data file
        success = read_qbuf(whichqbuf)
      endif
c02Mar        CALL QBUFIN(a,aa,hc,hs,fe,ec,es,u,v,npts,mon,dm,da, success)
        if (.not. success) then
	  write (*,*) 'main: Can not open file ',  FILE_BUFIN(1:50)
	  stop
	endif
c...     QBUFIN reports the symmetry, but it is ignored due to LPI file read later.
c        if (noisy) then
c    	    call reportsymmetry(ISYSTM)
c        endif

c	 !! energy resolution tol for spurion hitting the bragg peak
cNope      if (abs(qb_energy_es) .gt. tiny) then
cNope        etol= abs(qb_energy_es)
cNope      else
cNope        etol= 1.0 
cNope      endif

cNope        the_magq = 0.0
cNope        do iscan = 1, nscan-1
cNope          temp_magq = sqrt((qarr(1,iscan) - qarr(1,iscan+1))**2 
cNope     &	  	+ (qarr(2,iscan) - qarr(2,iscan+1))**2)	!! in AA^-1
cNope	  if (temp_magq .gt. the_magq)  the_magq = temp_magq 
cNope        enddo
cNope        if (the_magq  .gt. tiny) then
cNope          tol= the_magq
cNope        else
cNope          tol= 0.05 		!! probably qb_space_hs(3) = 0 or nscan=1
cNope        endif


c...	Eventually only use common block variables qb_*.
c...	For now copy the contents.
      do i = 1, 3
        a(i) =  qb_lattice_a(i)
        aa(i) =  qb_lattice_angle(i)
	hc(i) = qb_space_hc(i)
	hs(i) = qb_space_hs(i)
	u(i) = qb_orient_u(i)
	v(i) = qb_orient_v(i)
      enddo 
      ec = qb_energy_ec
      es = qb_energy_es
      fe = qb_energy_fe
      mon = qb_energy_fix
      dm = qb_energy_dM
      da = qb_energy_dA	
      npts = qb_time_scanpts	

         if (.not. verify_qbuf(a,aa,hc,hs,fe,ec,es,u,v,npts,mon,dm,da, 
     &		 .true.)) then
           write (*,*) 'main: Continuing anyway. #2'
	   stop
         endif

      WRITE(*,313) 'E-range: ',ec-es*(npts-1)/2.0,' to ',
     &              ec+es*(npts-1)/2.0,' with center (hkl): (',
     &              hc(1),hc(2),hc(3),').'
      WRITE(unitlog,313) 'E-range: ',ec-es*(npts-1)/2.0,' to ',
     &              ec+es*(npts-1)/2.0,' with center (hkl): (',
     &              hc(1),hc(2),hc(3),').'
 313  FORMAT(a9,f5.1,a4,f5.1,a21,f4.1,1x,f4.1,1x,f4.1,a2)

        write(unitspur,*)     'SPURION Report File'
        write(unitspur,'(a)') '---------------------'

cMoved        call version(aversion, date_stamp)
        write (unitspur,'(a39, x, a39)')   
     &      FILE_BUFIN, FILE_LPI  
        write (unitdebug,'( a39, x, a39)')
     &		FILE_BUFIN, FILE_LPI  
c        write (unitspur,'(a30, x, a24, x, a24)') version_stamp,  
c     &      FILE_BUFIN, FILE_LPI  
c        write (unitdebug,'(a30, x, a24, x, a24)') version_stamp, 
c     &		FILE_BUFIN, FILE_LPI  

      
        if (showcalls) write (*,*) 'main: call init_vectors'
        call init_vectors(a, aa, bmat, bmatuv, b, 
     &      uu, vv, uxv, uuv, vuv, u, v, uxvi, 
     &	    zetatitle, etatitle, relcircsize, angle, xstar, ystar)
c02Mar     &	    zetatitle, etatitle, c, angle, xstar, ystar)
      
cDebug	write (*,*) 'main: after init_vectors:'
cDebug	do i = 1, 3
cDebug	  write (*,*) 'main: bmat(', i, ', j)=',  (bmat(i,j), j=1,3)
cDebug	enddo
cDebug	do i = 1, 3
cDebug	  write (*,*) 'main: bmatuv(', i, ', j)=',  (bmatuv(i,j), j=1,3)
cDebug	enddo
	
      endif		!! if (currentoptions .eq. 0) 


      if (freshplot) then 
        if (showcalls) write (*,*) 'main: call init_vectorlabels'
 60     call init_vectorlabels(rlu, u, v, zetatitle, etatitle, 
     &	xtitle, ytitle, x2title, y2title,angle, pi)
        if (showcalls) write (*,*) 'main: call init_scanvectorlabels'
 63     call 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)
        call prep_elastic_spur(npts, ec, es, fe, hc, hs,  mon)
	
C 	suppose there were no suitable points
        if (nscan.eq.0) then
          write(*,'(2a)') 'PROBLEM: no points in scan!  Review and',
     &     ' modify these values.'
ccc          scale='m'
ccc          go to 781   ! shortcut to modify menu
          scale = char_mod
          addoption = 0
          call set_bit(addoption,bit_mod,1)
          call set_bit(currentoptions,bit_mod,1)
        endif
        if (showcalls) write (*,*) 'main: call init_qrange'
        call init_qrange(qmin, qmax, qlowmin, qlowmax, 
     &	      nscan, kiarr, kfarr, phimax)
        if (showcalls) write (*,*) 'main: call plot_center'
        call plot_center(u, v, uuv, vuv, hklcen)
cCommon        call plot_center(xmin, xmax, ymin, ymax, 
cCommon     &      u, v, uuv, vuv, hklcen)

cDebug	write (*,*) ' q-step from qb_space_hs= ',
cDebug     &	    sqrt( qb_space_hs(1)**2 + qb_space_hs(2)**2  + qb_space_hs(3)**2)
cNope        the_magq = 0.0
cNope        do iscan = 1, nscan-1
cNope          temp_magq = sqrt((qarr(1,iscan) - qarr(1,iscan+1))**2 
cNope     &	  	+ (qarr(2,iscan) - qarr(2,iscan+1))**2)	!! in AA^-1
cNope	   write (*,*) iscan, ' q-step = ', temp_magq
cNope	  if (temp_magq .gt. the_magq) then
cNope	    the_magq = temp_magq 
cNope	  endif
cNope	  	
cNope        enddo
cNope        if (the_magq  .gt. tiny) then
cNope          tol= the_magq
cNope        else
cNope          tol= 0.05 		!! probably qb_space_hs(3) = 0 or nscan=1
cNope        endif


      endif		!! if (freshplot)
       

      
C  Calculate the reciprocal lattice i.e., the (hkl) and the
C  zeta and eta coordinates assigned to each i,j which will be
C  placed in 5  m x m matrices.
      
C  only do these calculations once

      if ( .not.prev) then
c       !! AA^-1 resolution tol for spurion hitting the bragg peak

ccc      if ( (.not.prev) .or. freshplot) then
         if (showcalls) write (*,*) 'main: call calclattice'
         CALL calclattice(uuv,vuv,bmatuv, m_rlv, hh,kk,ll,
     &        zeta,eta, hklcen,u,v)
	  WRITE (unitdebug,*) 'post-calclattice>>'
cDec01         CALL calclattice(uuv,vuv,bmatuv,m,hh,kk,ll,zeta,eta,
cDec01     &        hklcen,u,v)
         if (showcalls) write (*,*) 'main: call fsqlattice'
         call fsqlattice(m_rlv,hh,kk,ll,struct,sfhkl,absent)
cDec01         call fsqlattice(m,hh,kk,ll,qmax,struct,sfhkl,absent)
	  WRITE (unitdebug,*) 'post-fsqlattice>>'
c	   !! Set default value of tolerance after fsqlattice()
c	  !!  because might need the value of natoms. 
c	  tolerance = 4.0	!! should be called base in this context.
	 			!!  Tol for matching intensities
c	  sortlogscale = .true.
	  sorttol = 100.0*tiny *float(natoms)
	  base = 2.0
c02Mar	  tolerance = 100.0*tiny *float(natoms)
	  sortlogscale = .false.
          prev=.true.
      endif

cChanged      if (freshplot .or. bit_on(addoption, bit_sortlog)
cChanged     &      .or. bit_on(addoption, bit_rescale)) then
      if (freshplot .or. bit_on(add_structop, sbit_logsort)
     &      .or. bit_on(addoption, bit_rescale)) then
c	 sorting the intensities in descending order

	  write (*,*) 'calling useDescending2'
          if (showcalls) write (*,*) 'main: call indexx2d'
	  call indexx2d(m_rlv, sfhkl, indx1, indx2, useDescending,
     &	   	qmax, hh, kk, ll, zeta,eta) 
cCommon	  call indexx2d(m_rlv, sfhkl, indx1, indx2, useDescending,
cCommon     &	   	qmax, hh, kk, ll, zeta,eta, xmin, xmax, ymin, ymax) 
c	  write (*,*) 'main: calling unique'
          if (showcalls) write (*,*) 'main: call unique'
	  call unique(m_rlv, sfhkl, absent, indx1, indx2,  
     &	  	indxunique1, indxunique2, multiplicity,  limlo, limhi,
     &		nunique, uniqueLabel,  uniqueRank,  sorttol, base,
c02Mar     &		nunique, uniqueLabel,  uniqueRank,  tolerance, 
     &		sortlogscale, useDescending, hh, kk, ll)
          if ((nabsent .gt. 0) .or. (ntiny .gt. 0)) then
	    spurbit = 1
            call bit(spurhist,7,spurbit,1)
	  else
	    spurbit = 0
            call bit(spurhist,7,spurbit,1)
	  endif
	endif		!! if (freshplot .or. bit_on(addoption, bit_sortlog)) 


cMoved      endif


      if (currentoptions .eq. 0) then
        relcircsize = 1.0
c02Mar        c=1.0
C  	Calculate a reasonable circle size for the reciprocal
C  	lattice points.
        circsize= (xmax-xmin)/120.0
      endif      !!  if (currentoptions .eq. 0) then

C --------------BEGINNING OF PLOTTING LOOP---------------

      if (freshplot)  then
C	  Begin pgplot routine, set character size, and axes.
C 	 Call pg_init with idev=
c		0 for XSERVE
c		1 for XWIN
c		2 for grey PS
c		3 for colour PS
        if (showcalls) write (*,*) 'main: call pg_close'
        call pg_close(idev, just_erase)             ! close PGPLOT
        if (showcalls) write (*,*) 'main: call pg_init'
 68     CALL pg_init(idev, rainbow, just_erase)
cccc 68     call plot_begin(idev)
        just_erase = .true.
      endif 		!! if (freshplot) 
      
C  setup the plot: axes, ticks, tick labels

      if (freshplot) then
	if (showcalls) write (*,*) 'main: calling plot_axes'
cDebug	write (*,*) 'plot_axes range: ', xmin, xmax, ymin, ymax
        if (showcalls) write (*,*) 'main: call plot_axes'
 67     CALL 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)
cDebug	write (*,*) 'plot_axes range: ', xmin, xmax, ymin, ymax
      endif
	
      call pgupdt		!! flush the graphics buffer

      if (freshplot) then

c02Mar        ccircsize = c*circsize
        ccircsize = relcircsize*circsize

c        if (makescanplot) then
c          if (showcalls) write (*,*) 'main: call plot_scan'
c          call plot_scan(constq, ec, es, hc, hs, mon, nscan, 
c     &		  esparr, pointsize, relcircsize, ispur, jspur,
c     &		nspurtot, sfhkl, whichqbuf, version_stamp)
cc02Mar          call plot_scan(constq, ec, es, hc, hs, mon, nscan, whichqbuf)
c	endif

C 	 Plot the reciprocal lattice.
        if (showcalls) write (*,*) 'main: call plotlattice'
        CALL plotlattice(ccircsize,m_rlv,zeta,eta,hh,kk,ll,qmax,
     &  	gray, rainbow, struct,sfhkl,absent, 
     &		replot,  overridesize, pointsize, plotunique, uniqueLabel, 
     &  	 hide_absences, nscan, qarr)
cCommon     &  	 hide_absences, xmin, xmax, ymin, ymax)
c01Nov      CALL plotlattice(c*circsize,m,zeta,eta,hh,kk,ll,qmax,
c01Nov     &  gray,struct,sfhkl,absent)

        if (makescanplot) then
          if (showcalls) write (*,*) 'main: call plot_scan'
c	  ! first load in spurion info without plotting on LATTICE page.
	  call 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, .true.,  bmatuv)
cDebug         write (*,*) 'main: returned from cmd_spurion'
         call cmd_aluminum(nscan, npts, minq, maxq, magq, 
     &   	nearby, qarr, kiarr, kfarr, 
     &		h, hc, hs, incmonarr, incanaarr,
     &	 	constq, spurhist, mon, pi, fe, phimax, .false.)
	 call 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, .false.)
          call plot_scan(constq, ec, es, hc, hs, mon, nscan, 
     &		  esparr, pointsize, relcircsize, 
     &		nspurtot, nelastictot, nincohmon, nincohana, 
     &		sfhkl, phimax, whichqbuf, version_stamp)
c02Mar          call plot_scan(constq, ec, es, hc, hs, mon, nscan, whichqbuf)

cDebug         write (*,*) 'main: post plot_scan: nspurtot=', nspurtot
cDebug	 do i = 1, nspurtot
cDebug          write (*,*) 'main: post plot_scan: ', i, scanpt_spot(i),
cDebug     &		  ispot(i), jspot(i), ispur(i), jspur(i)
cDebug         enddo   

	endif

          if ((nabsent .gt. 0) .or. (ntiny .gt. 0)) then
	    spurbit = 1
            call bit(spurhist,7,spurbit,1)
	  else
	    spurbit = 0
            call bit(spurhist,7,spurbit,1)
	  endif

cDebug       write (*,*) 'returned from plotlattice'

C	  Plot Q, ki, kf  with a heavy line.
        if (showcalls) write (*,*) 'main: call plot_Q_ki_kf'
        call plot_Q_ki_kf(qarr, kiarr, nscan)
	
	reset_freshplot = .true.

C	  End of default minimalist plot
      else
	reset_freshplot = .false.
      endif 		!! if (freshplot) 
      call pgupdt		!! flush the graphics buffer

C  Beginning of command sequence

c      call annotate_spurhist(spurhist, annotation, unitdebug) 
cDebug      write (*,*) 'finished calling annotate_spurhist ',unitdebug
c      call annotate_spurhist(spurhist, annotation, 0) 
cDebug      write (*,*) 'finished calling annotate_spurhist ',0

ccc 781   continue


c------------------------------------------------------------
c 	Commands which super-impose on the current graphics
c	The program remembers that these commands have been requested.
c------------------------------------------------------------

C  Plot aluminum Ewald spheres that are within qmax in green      
      if (bit_on(addoption,bit_al) .or. 
     &	      (freshplot .and. bit_on(currentoptions,bit_al))) then 
	 spurbit = 1
         call bit(spurhist,2,spurbit,1)
	 if (bit_on(addoption,bit_al)) then
c           if (.not. constq) then
c             write (*,*) 'For clarity, in case of a q-scan, ', 
c     &		'this routine only reports some of the Al Aleward spheres.'
             write (*,*) 'For clarity, ', 
     &		'this routine only reports some of the Al Aleward spheres.'
120          WRITE (*,'(a, i3, a, $)') 'Enter which scan point: [1-', 
     &		nscan, ', or 0 for all]: '
	     read (*,*, err=120) draw_al_iscan
	     if ((draw_al_iscan .lt. 0) .or. 
     &		 (draw_al_iscan .gt. nscan)) goto 120
c           else
c      	     write (*,*) ' Constant q scan, so only 1 set of ', 
c     &		'Al Alewalds spheres to draw'
c             draw_al_iscan = 1
c           endif      
	 endif
	 
         call cmd_aluminum(nscan, npts, minq, maxq, magq, 
     &   	nearby, qarr, kiarr, kfarr, 
     &		h, hc, hs, incmonarr, incanaarr,
     &	 	constq, spurhist, mon, pi, fe, phimax, .true.)
c         call cmd_aluminum(nscan, npts, gal, minq, maxq, magq, hklal, 
c     &   	alspur,nearby, qarr, kiarr, kfarr, 
c     &		h, hc, hs, incmonarr, incanaarr,
c     &	 	constq, spurhist, mon, pi, fe, nal_match)

      else if (freshplot) then
c...	Refresh the al spurion info for the scan plot but don't plot on Lattice.
         call cmd_aluminum(nscan, npts, minq, maxq, magq, 
     &   	nearby, qarr, kiarr, kfarr, 
     &		h, hc, hs, incmonarr, incanaarr,
     &	 	constq, spurhist, mon, pi, fe, phimax, .false.)

      endif

C  Plot, count, and record spurious Q's
      if ((bit_on(addoption,bit_spur) .and. 
     &		 bit_on(currentoptions,bit_spur)).or. 
     &	      (freshplot .and. bit_on(currentoptions,bit_spur))) then 
	 spurbit = 1
         call bit(spurhist,0,spurbit,1)

cDebug	write (*,*) 'main: before cmd_spur:'
cDebug	do i = 1, 3
cDebug	  write (*,*) 'main: bmatuv(', i, ', j)=',  (bmatuv(i,j), j=1,3)
cDebug	enddo

	 call cmd_spurion(  m_rlv, npts, nscan, nspurtot, nspur, 
     &      qspur, espur, qsparr, esparr, kiarr, kfarr,
     &      circsize, ki, kf, zeta, eta,
     &      hh, kk, ll, 
cInc     &      ispur, jspur, hspur, kspur, lspur, ielst, jelst,
     &      absent, replot, sp, spurtitle, 
     &      ec, es, h, hc, hs, constq, nelastictot, .false.,  bmatuv)
cDebug         write (*,*) 'main: returned from cmd_spurion'
      endif

cDebug      write (*,*) 'main: post cmd_spurion: nspurtot=', nspurtot
cDebug      do i = 1, nspurtot
cDebug        write (*,*) 'main: post cmd_spurion : ', i, scanpt_spot(i),
cDebug     &		  ispot(i), jspot(i), ispur(i), jspur(i)
cDebug      enddo   
cDebug
cDebug	write (*,*) 'main: after cmd_spur:'
cDebug	do i = 1, 3
cDebug	  write (*,*) 'main: bmatuv(', i, ', j)=',  (bmatuv(i,j), j=1,3)
cDebug	enddo
      
C  Plot and report spurious Q's due to incoherent scattering.
      if (bit_on(addoption,bit_incoh) .or. 
     &	      (freshplot .and. bit_on(currentoptions,bit_incoh))) then 
	 spurbit = 1
         call bit(spurhist,1,spurbit,1)
	 call 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, .true.)
      else if (freshplot) then
	 call 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, .false.)
      endif
      
C  Draw the u and v orientation vectors with solid, magenta lines.
      if (bit_on(addoption, bit_orientvec) .or. 
     &	 (freshplot .and. bit_on(currentoptions,bit_orientvec))) then 
	 spurbit = 1
        call bit(spurhist,4,spurbit,1)
	call cmd_orientvecs(uuv, vuv)
      endif

C  label Q, ki, and kf
      if (bit_on(addoption,bit_labelvec) .or. 
     &   (freshplot .and. bit_on(currentoptions,bit_labelvec))) then 
         spurbit = 1
         call bit(spurhist,5,spurbit,1)
         call cmd_labelvec(kiarr, kfarr, kistrt, kistop, kfstrt, kfstop,
     &      nscan, over, spurhist)
cCommon         call cmd_labelvec(kiarr, kfarr, kistrt, kistop, kfstrt, kfstop,
cCommon     &      nscan, over, spurhist, xmin, xmax, ymin, ymax)
      endif

C  place a grid on Lattice plot
      if (bit_on(addoption,bit_grid) ) then 
c	  Ask for style and step size      
         call cmd_grid(.true., .true., uuv, vuv, rlu, xstar, ystar)
      endif
      if (freshplot .and. bit_on(currentoptions,bit_grid)) then 
c	  Keep previous values for style and step size      
         call cmd_grid(.false., .true., uuv, vuv, rlu, xstar, ystar)
      endif
      

C  Add annotation text to graph
      if (bit_on(addoption, bit_annote) .or. 
     &	 (freshplot .and. bit_on(currentoptions,bit_annote))) then 
c02Mar        call load_itemvalue(tol, etol, c, struct,
        call load_itemvalue(tol, etol, 
cMoved     &		relcircsize, struct,
cMoved     &	    gray, rainbow, plotunique, sortlogscale, tolerance, 
     &      phimax, rlu, right, selct_sp, selct_sfhkl, over,
     &	    m_rlv, nabsent, npresent, nspurtot, nelastictot, nal_match,
     &	    nincohmon, nincohana, grid_style, grid_size1, grid_size2,
     &	    hide_absences)
cc	call cmd_annotation(unitdebug, c, version_stamp) !! to the file unitdebug
c	call cmd_annotation(0, c, version_stamp)	 !! to the screen
cc	call cmd_annotation(-1, c, version_stamp)	!! no ASCII printout
cDebug      write (*,*) 'main:1 i, j', indxunique1(1), indxunique2(1)
c02Mar	call cmd_annotation(0, c, version_stamp,
	call cmd_annotation(0, version_stamp,
     &		sfhkl, nunique, uniqueLabel, 
     &      indxunique1, indxunique2, multiplicity, 
     &	    pointsize, circsize, overridesize, rainbow)	 !! to the screen
      endif


c------------------------------------------------------------
c 	Now the plot already has stuff on it, so no longer fresh.
c------------------------------------------------------------
      if (bit_on(currentoptions,bit_al)) then 
	freshplot = .false.
      endif
      if (bit_on(currentoptions,bit_spur)) then 
	freshplot = .false.
      endif
      if (bit_on(currentoptions,bit_incoh)) then 
	freshplot = .false.
      endif
      if (bit_on(currentoptions,bit_orientvec)) then 
	freshplot = .false.
      endif
      if (bit_on(currentoptions,bit_labelvec)) then 
	freshplot = .false.
      endif

      if (reset_freshplot) then		!! basic plot made.
	freshplot = .false.
      endif
c------------------------------------------------------------
c 	commands which only output to TTY or text files
c------------------------------------------------------------

C  Sf^2 printout
      if (bit_on(addoption,bit_reflect)) then 
         call  cmd_selectReflection( selct_sp, spurtitle, m_rlv, 
     &      hh, kk, ll, zeta, eta, absent, sfhkl, selct_sfhkl,
     &		qlowmin, qmax)
cWrong	freshplot = .false.
      endif

C generate angles like an ICP dry run
      if (bit_on(addoption,bit_ICP)) then 
	 call cmd_genICPrun(mon, npts, bmatuv, h, q, hc, hs, ki, 
     &      ec, es, en, ei, ef, fe, da, dm, a2, a3, a4, a6, pi)
cWrong	freshplot = .false.
      endif

C output values to a buffer
      if (bit_on(addoption,bit_qbufout)) then 
c	 Double check that we have not forgotten to update something.
         if (.not. verify_qbuf(a,aa,hc,hs,fe,ec,es,u,v,npts,mon,dm,da, 
     &		 .true.)) then
           write (*,*) 'main: Continuing anyway.'
	   stop
         endif
         success = write_qbuf(unitbufout, 'All', junkline)
         success = write_qbuf(0, 'All', junkline)	! also to screen
	 write (*,*) ' '
c02Mar         CALL QBUFOUT(a,aa,hc,hs,fe,ec,es,u,v,npts,mon,dm,da)
cWrong	freshplot = .false.
      endif


C  HELP!!
      if (bit_on(addoption,bit_help)) then 
c	Now using help-paul-unx.c, and ASCII file spurion.hlp
c  With empty path, use logical MSRHELP to determine where the help file is.
c	On VMS, "assign <helpdir> MSRHELP:"
c	On Unix,"setenv MSRHELP <helpdir>"
c	On Unix, "export MSRHELP='<helpdir>'"
c	 CALL HELP(HLP_PATH, HLP_FIL, KEY_WD)
c		In *.hlp file, spurion is currently only 1 level available.
c		So use key word to skip to it.
c		If have libraries (like ASK, DOPLOT, MINUIT, etc.)
c		then they should all be level 1 and key_wd = ''
	 call help('', 'spurion', 'spurion')
cWrong	freshplot = .false.
      endif

 545  if (bit_on(addoption,bit_quit)) then 
C 	---- closing sequence--------
        if (showcalls) write (*,*) 'main: call pg_close 545'
        call pg_close(idev, .false.)             ! close PGPLOT

        close(unit=unitspur)
        close(unit=unitstruct)
        close(unit=unitlog)
        close(unit=unitdebug)
        WRITE(*,'(a)') 'Ending SPURION.'
      endif	!!  if (bit_on(addoption,bit_quit)) then 

c...	Ignore unrecognised commands.
c...	This piece of code is now complicated by toggling of 
c...	some commands via settting addoption = 0.
c      if ((addoption .eq. 0) .and. (scale .ne. ' ') 
c     &      .and. (.not. bit_on(addoption,bit_default))) then 
c	write (*,*) 'Unrecognized command: "', scale, '"'
c      endif	!! if (addoption .eq. 0) then 

      
C ---- loop back to menu ----

555     continue
cDebug        write (*,*) 'main: end of loop. freshplot =',
cDebug     &		 freshplot, ' prev=', prev
      enddo 		!! do while (.not. bit_on(addoption,bit_quit)) 

      END		!! main

C*******************************************************
C
C      Subroutines
C
C*******************************************************
      
      subroutine intro(version_stamp, date_stamp)
      implicit none

c...	passed variables
      character*(*) version_stamp
      character*(*) date_stamp
      
c...	local variables
      integer i, j 

C     Program introduction

	call version(version_stamp, date_stamp)
	

      WRITE(*,*)
      WRITE(*,105) 
      write(*,110) version_stamp(1:72)
      WRITE(*,110) 'Written by Jeff Lynn, the dark ages, ca. 1970'
      WRITE(*,110) 'Graphical version written by Carl Adams, Feb. 2001'
      WRITE(*,110) 'Modified by Tanya Riseman, Nov, 2001 - May, 2002'
      WRITE(*,105) 
      WRITE(*,*)
      WRITE(*,*)

c      !! To see if screen is set to 80 characters wide
      write (*, *) 'Check that screen is 80 characters wide', 
     &      ' and 24 lines long:'
      WRITE(*,'(80(a1))') ((char(mod(i,10) +ichar('0')), i = 1, 10), 
     &	      j = 1, 8)
      do i = 2, 23
        write (*,*) i, '= line number' 
      enddo
      WRITE(*,'(80(a1))') ((char(mod(i,10) +ichar('0')), i = 1, 10), 
     &	      j = 1, 8)
      

105   format (2x, 76('*'), 2x) 
110   format (2x, '* ', a72, ' *', 2x) 

      return
      end
      
C*******************************************************
      subroutine version(astring, date_stamp)
      implicit none
     
      character*(*) astring
      character*(*) date_stamp

c...			      123456789_123456789_123456789_
      call datestamp(date_stamp)
CCC      write (astring, '(a, a, a, a)') 'spur_0201_$(OS=g77,sgi).'	!! 25
c     &      		     ' V0.4.4 Jan, 2002.',			!! 19
c     &      		     ' Used on ', date_stamp			!! 10+20
      write (astring, 110) date_stamp	
110   format ( 'spur_030115_$(g77,sgi). V0.4.16 Jan, 2003.',
     &		 ' Used on ', a20)
c110   format ( 'spur_020514_$(g77,sgi). V0.4.15 May, 2002.',
c     &		 ' Used on ', a20)
         

      return
      end

C*******************************************************
     
