pro plotbth,filename
;+
; NAME:                 plotbth v3.0 (IDL2) 
;
; PURPOSE:              Plot bsd file - Menu structure, quick look plotter 
;				overplotting of theory spectra.
;
; CALLING SEQUENCE:     plotbth[,bsdfilename]
;			(else prompted for [dir]filename[*]) 
;
; INPUTS:		Valid BSD file (see MKBSD.exe)
;			Matching SPC file (See BSDFIT.exe)
;
; OUTPUTS:              Plots to graphics devices.
;
; RESTRICTIONS:         This is a procedure, unlike plotbda, so data is
;                               not available at program termination. Use
;                               get_bsd to read in a BSD file to IDL.
;				use get_bth to read in the bth (SPC) file.
;                           You have to explicitly print postscript files after
;                               the program has finished. use PPRINT.
;
; PROCEDURE:             Read selected BSD file. Read selected BTH file.
;			 Setup plot parameters at main menu. Run series of 
;			 plots, with behaviour modified by plot submenu.
;
; MODIFICATION HISTORY:
;               author:       Andy Phillips , mssl 1991
;               vn 1.8:       atp, 15-jun-92
;		vn 2.1:	      atp, 07-jul-92 mods for spc structures
;		vn 3.0:	      atp, 23/11/92 - rationalise with ISAS
;		RDB  01-Mar-93	Changed "BTH" to "bth"
;		RDB  26-May-93	Modified oplot calls.
;
;-
;
; history: dont really want all this in doc_library
;
;  v0.0 atp dumped files contents
;  v1.0 atp M.S.S.L 15 - 5 - 91 rewrite, menu structure
;  v1.01 atp hardcopy
;  v1.1 atp idl1 demo version complete, still missing options, mssl specific
;  v1.2 atp idl1 non-mssl specific, missing options, added range.pro,bsdtime
;  v1.3 atp idl2 conversion functionally identical to v1.2 
;  v1.4 atp idl2 inserted error checking reorgainsed default plot variables
;	         into block at start, new bsdtime included.
;  v1.45 atp idl2 Now default plot unit at atart up is as defined externally.
;  v1.5  atp idl2 Fixed bugs/ relating to only one spectrum /channel 
;		 * scans if change channel etc..
;  v1.6  atp idl2 Fixed options 6 and 7. program is now complete aside from
;		  bugfixing.
;  v1.8  atp idl2 Copied plotbsd to plotbth. Overplotting of theory spectra.
;         af idl2 Two more routines - rdspc2_data, rdspc2_head 
;  v2.1  atp idl2 rewrote rdspc2_data etc... as rd_spc_head, rd_spc_rdmap,
;		  rd_spc_data, to give structures/consistency with bsd routines
;		  and to give unix compatibility. Also altered plotbth to give
;		  non - histogram primary component.
; 
;  19709::atp
;
;  if you intend to change anything please read comments.
;
; comments: 
;	Uncertainties should have option of error bars rather 
;	than +/- graphs.
;       DEFAULTS ARE BELOW
;       if you add another plotunit like an installable driver
;		eg dec sixel files, check defaults and code for reading them in
;		for soft plotunits eg tek search for string +++s+++
;		for hard copy search for string +++h+++
;	convention 1 = yes 0  = no.
;	spectra are selected by certain range criteria.
;	parsing of the selection string (spectstr) is accomplished in range.pro
;	it is not bullet-proof. The selection string contains a comma delinited
;	series of numbers, and/or two wild cards. 
;	* selects all spectra, the rest of the string is ignored.
;	- selects a range. eg 1-3 == 1,2,3.
;	bounds checking is performed, redundancy is not. it is perfectly
;	possible to replot the same spectrum several times eg.
;	1,2,3,1-2,3,3,3 will plot spectra 1,2,3,1,2,3,3,3.
;       Implicit assumption in range option that max(x) = x(max) , min(x) = x(0)
;       error handling uncertain.
;
; initialise local variables
; 1) PLOTTING DEFAULTS - no sanity checking. wrong values == crashed program.
;    plot unit now assumed to be initially defined externally.
;
 channel = 1;                     (int) Channel to plot  
 defpunit=!d.name;+++s+++h+++  (string) DEFAULT plot_unit (TEK,REGIS,PS,X)
 uncert = 0;                      (int) plot uncertainties flag
 waveyes = 0;                     (int) plot wavelength vs counts 
 theoryes = 1;			  (int) overplot theory spectrum
 chstring= '1';
 spectstr = '*';	       (string)	string for range of spectra. see above.
 orient = 'portrait';	       (string) orientation of postscript
;
;
 ans = ' ';				Menu variable
 true =1;                              	if !true then quit
 first_time = 1;			if true then no hardcopy
 dum = ' '; 				dummy variable for get_kbrd
 ier = 0;				range error variable
 curspect = 0;				plot variable
 j = 0;					internal counter
 version = 'v3.0';			version of this program
 lskip = 0;				plot skipper
 rxlow = 0;				lower plotting x limit
 rxhi = 0;				higher plotting x limit
 temp = 0;				swapping variable
 numbsd = 0;				number of spectra in bsd file (theor)
 numthe = 0;				number of spectra in theory file
 firchn = 0;				spectrum no. of 1st theory spct /chan
 pointr = 0;				rdmap pointr for theor file
 spcnam = ' ';				name of spc file (theor)
 bpcnam = ' ';				name of bpc file used for bsdfit (theor)
 thevers = ' ';				text with version no. etc (theor)
 balcal = ' ';				ionization balance calcn used
 thetime =0;				time assoc with theory file
 thnpts =0;				number of points in theory spectra
 thespect =0;				theory spectrum
 velfit = 0;				theory variable
 primary = 0;				theory variable
 ier=0;					theory return status.
 compbounds = 0;			compute bounds for wavelength?
 got_data = 0;				have already got data via get_bsd
 nsumspec = 1;
 nsumstr = '1';
 kplotdirn = 1;
 khcflag = 0;

;
; print program header
;
print,'Plotbth  ',version
print,'Default plot unit - ',defpunit
print,' '
;
;
;  if !have filename read it in
;
get_file = 0;
if (n_elements(filename) eq 0) then get_file = 1 else if (filename eq '') then get_file=1
if (get_file eq 1) then begin
   bcs_getfile,filename,defstr="bsd*.*"
   if (filename eq '-1') then return
   endif
;
; strip out identifier 
;
break_file,filename,dsk_log,dir,bsdnam,ext,filver              
;
; extract identifier. 
;
ident = bsdnam+ext
;                                                                               
;
;  get theory spectra filename
;
ans = ' '
theorname = ' '
;
theorname = dsk_log+dir+"bth"+strmid(bsdnam,3,(strlen(bsdnam)-3))+ext
print,format = '("Looking for theory spectra in ",A)',theorname
print,format = '($,"Is this correct ? (y/n)[y] ")'
read,ans
if (ans eq '') then ans = 'y'
if (strlowcase(ans) eq 'n') then begin
	print,format='($,"Theory spectra file ?")'
	read,theorname
	endif

;
; read header info 
;
; rdbsd2head,filename,totspc,numchn,ptrdmp,floval,floset,flcorr,peak,
;  pknum,vernum,tdf,tdl
rd_bsd_header,filename,header
rd_bsd_rdmap,filename,header,rdmap
;
; read theory spectra header info
;
rd_spc_header,theorname,bth_header
rd_spc_rdmap,theorname,bth_header,bth_rdmap
;
;rdspc2_head,Theorname,numbsd,numthe,firchn,pointr,spcnam,bpcnam,balcal,thevers 
;
;
; Plotting defaults - do not alter below this line 
;
  plot_unit = defpunit; 
  plot_unit = strlowcase(plot_unit);  +++s+++
  if ((plot_unit eq 'tek') or (plot_unit eq 'regis')) then begin
	tekyes = 1; 
	no_stop = 0; 	                if true then wait for rtn b4 menu
  	endif else begin
	tekyes = 0;
	no_stop = 0;
	endelse
 
 if ((plot_unit eq 'ps')) then begin; +++h+++
	tekyes = 0;
	device,/portrait
	no_stop = 1;                 if true then is hardcopy+dont wait  
        endif

 if (plot_unit eq 'x') then device,retain=2

if max(header.numchn(*)) eq 0 then begin
 print,'no data for this file.'
 stop
 endif

range,spectstr,max(header.numchn(*)),spectra,ier,1; parse default specstr
ntodo = n_elements(spectra); number of spectra to plot in a run
;
; set up plot 
;
set_plot,plot_unit
!p.font = -1
!p.multi=0
; this is necessary to remove utplot probs
!x.ticks=0&!x.tickname=''&!x.tickv=0.0
!y.ticks=0&!y.tickname=''&!y.tickv=0.0
!x.range=0
!y.range=0
multiflag=0
!p.charsize = 0.9			; To fit title on page.
;
; output header info
;
print,'File - ',filename,' is a version',header.vernum,' file'
print,'Total number of spectra ',header.totspc
print,'First Spectrum starts at ',fmt_tim(header.tdf,/msec)
print,'Last spectrum starts at ',fmt_tim(header.tdl,/msec)
print,' '
print,'File - ',theorname,' version info:',string(bth_header.version)
print,'Number of BSD    spectra',bth_header.numbsd
print,'Number of Theory spectra',bth_header.numthe
print,'SPC filename ',string(bth_header.spcnam)
print,'BPC filename ',string(bth_header.bpcnam)
print,'Ion balance calculation ',string(bth_header.balcal)
print,' '
print,'press any key to continue.'
dum = get_kbrd(1)
;
; Commence loop
;
while (true eq 1) do begin
 print,' '
 print,' '
 print,' PlotTheory ',version,' Options '
 print,'  '
 print,' 1) Change channel (*=all)[',chstring,']'
 print,' 2) Change spectrum (*=all)[',spectstr,']'
 print,' 3) Plot Uncertainties if they exist (1/0)(toggle) [',uncert,']'
 if (!d.name eq 'PS') then begin
  print,' 4) Set plot device [',plot_unit,'-',orient,']'
 endif else begin
  print,' 4) Set plot device [',plot_unit,']'
 endelse
 print,' 5) Plot vs wavelength(=1) or bins(=0) (toggle) [',waveyes,']'
 print,' 6) Show Bsd_Header Structure '
 print,' 7) Set plot range (0,0)=all [',rxlow,',',rxhi,']'
 print,' 8) Overplot Theory Spectra if exist (1/0)(toggle) [',theoryes,']'
 print,' 9) Plot light curves for current channel/spectra'
 print,' Q) Quit '
 print,'  '
 print,' Any other key to plot current selection '
 print,'  Sequence Commands:'
 print,'   During a Sequence of Plots press Q or M to return to menu'
 print,'   type ? during plot mode for list of commands '
 print,' '
 ans = get_kbrd(1); 
 case ans of 
  '1': 	begin
	chstring= ' '; reinitialise
	print,'Enter new channel number (1...4)';
	print,'Or enter * for all channels on one page'
	print,format='($,"Channel")'

	read,chstring
	if (chstring eq '*') then begin
		channel = 1
		multiflag=1
		endif else begin
		channel = fix(chstring); cast to int
		multiflag=0
		endelse

	if ((channel gt 4) or (channel lt 1)) then begin; assume 4 channels
	    print,'Illegal choice - retry from menu'
	    channel = 1;
                    chstring = string (channel);
	    print,'press any key to continue. '
	    print,' '
	    dum = get_kbrd(1); halt to read error message
	    endif 
       	end
   '2': begin
        nsumstr='1';
	print,'Enter new spectrum number(1...',string(format='(i4)', $
                 header.numchn(channel-1)),')'
	print,' ("*" for all or a comma delimited series "1,3,5" or a range "1-3")'
	print,' '
	print,format='($,"Number(s)")'
	read,spectstr
	range,spectstr,header.numchn(channel-1),spectra,ier,1
	if (ier eq 2) then begin 
		print,'Error in choice - try again from menu'
		print,'press any key to continue. '
		print,' '
		spectra = 1
		spectstr = string(spectra)
		ntodo = 1
		dum = get_kbrd(1);
		endif else begin
		ntodo = n_elements(spectra)
		if (ntodo gt 1) then begin
			; more than one in sequence
			print,' During sequence press return for next plot '
			print,'  or M or Q for menu unless plotting to file '
			print,' '
			print,' press any key to continue.'
			dum = get_kbrd(1);
			endif else begin
			; only one in sequence
			endelse			
		endelse
       	end
   '3':	begin
	if (uncert eq 0) then begin
		uncert = 1 
		print,'Plotting uncertainties as +/- limits'
		print,'press any key to continue. '
		print,' '
		dum = get_kbrd(1)
		endif else uncert = 0
	end
   '4':	begin
	no_stop = 0; 	reset no_stop
	if (strlowcase(plot_unit) eq 'ps') then device,/close; +++h+++
	print,'Plot device      TEK - Tektronix type display e.g. kermit '
	print,'                 X - X workstations '
	print,'                 PS - PostScript file, for creating hardcopies'
	print,'                 REGIS - Regis terminals '
	print,' ';+++s+++h+++
	print,format='($,"Plot Device")'
	read,plot_unit
	plot_unit=strlowcase(plot_unit)
	case (plot_unit) of
	  'tek':begin; 
		tekyes = 1
		print,'Tektronix Type display - After plot press any key for menu'
		print,' '
		print,' press any key to continue.'
                                dum = get_kbrd(1)
		end
	    'x':begin
		tekyes = 0
                	print,'X Workstation display'
		print,' '
		print,' press any key to continue.'
		dum = get_kbrd(1)
		end
	   'ps':begin
		tekyes = 0
		no_stop = 1; dont stop in hardcopy output
		print,format='($,"Landscape or Portrait [P]")'
		read,dum
		if (dum eq '') then orient = 'portrait' 
		if (strlowcase(strmid(dum,0,1)) eq 'l') then orient='landscape' else orient = 'portrait'
		print,'PostScript file - plot(s) will be output without waiting'
		print,'  for keyboard input, selecting option 4 or quitting will '
		print,'  close the current file. Print with command pprint'
		print,' '
		print,' '
		print,' press any key to continue.'
		dum = get_kbrd(1)
		end
	'regis':begin
		tekyes = 1
		print,'Regis type display - After plot press key for menu'
		print,' ' 
		print,' press any key to continue.'
		dum = get_kbrd(1)
		end

;    +++s+++h+++ insert documentation and tekyes/nostop
;	flag conditions here

	ELSE:   begin
		print,'Error in choice - try again from menu'
		print,' resetting to default plot unit '
		print,' press any key to continue. '
		print,' '
		dum = get_kbrd(1)
		plot_unit = defpunit;
		end
	endcase
	set_plot,plot_unit
	if (plot_unit eq 'ps') then $
			if orient eq 'landscape' then device,/landscape $
			 else device,/portrait
	if (plot_unit eq 'x') then device,retain=2
        end

   '5': begin
	if (waveyes eq 0) then begin 
		print,'Plotting counts vs wavelength'
		waveyes = 1
		 if ((rxlow eq 0) and (rxhi eq 0)) then begin
		    ; no need to alter any thing
		     Compbounds = 0
		    endif else begin
		     ; need to recompute bounds
		     Compbounds =1
		     endelse
		endif else begin
		print,'Plotting counts vs bin number'
		waveyes = 0
		compbounds = 0
		endelse
		print,'  '
		print,' press any key to continue.'
		dum = get_kbrd(1)
	end
   '6': begin
	print,'dump of header info '
	print,'  '
	help,/struct,header
        print,'--Return-to-menu----'&dum = get_kbrd(1);
	end
   '7': begin
 	 print,' '
	 if (first_time eq 1) then begin
 		print,' Please plot spectrum before zooming in'
	                print,'press any key to continue. '
	                print,' '
	                dum = get_kbrd(1); halt to read error message
		endif else begin
	 	rxlow = ' '; reinitialise
	 	rxhi = ' '
		compbounds = 0
		if ((multiflag eq 1) and (waveyes eq 1)) then begin
			   ;choose a channel
			    print,'Multi - Plot and Wavelength scale selected,'
			    print,' Can only select range in terms of bins.'
			    print,' If you are interested in a part of a particular channel but want to'
			    print,' specify a wavelength range, and want all channels on the page,'
			    print,' try selecting one channel (option 1), setting the range (option 7),'
			    print,' then selecting all channels again. '
			    compbounds = 1
			    endif

	 	if ((waveyes eq 1) and (multiflag eq 0)) then begin
			; first ensure we have correct wo,dw
                        rd_bsd_data,filename,[curspect],header,rdmap, $
                           index,data,[1,2,3,4]
			wo=index(channel-1,0).wo
			dw=index(channel-1,0).dw
                     	wh=max(data(channel-1,0).wave)
			print,'Enter plot range between',wo,wh
			print,'(Enter 0 and 0 to reset to all)'
			print,format='($,"Low ")'
			read,rxlow
			print,format='($,"High")'
			read,rxhi
			rxlow = float(rxlow); cast to float
			rxhi = float(rxhi)
			if (rxlow gt rxhi) then begin ; lowest -> rxlow
				temp = rxlow;
				rxlow = rxhi
				rxhi = temp
				endif		
		
			if ((rxlow lt wo) or (rxhi gt wh)) then begin; bound check
				if ((rxlow eq 0) and (rxhi eq 0)) then begin
					print,'Resetting plot range'
					compbounds = 0
					endif else begin
		    			print,'Illegal choice - retry from menu'
	    	    			rxlow = 0&rxhi = 0
					compbounds = 0
	    	    			print,'press any key to continue. '
	    	    			print,' '
	    	    			dum = get_kbrd(1); halt to read error message
					endelse
	    	    		endif else begin; not illegal
				; we have to redo bounds every plot.
				compbounds = 1
				;convert to bins 
				rxlow = fix(((rxlow*1.0-wo)/dw))
				rxhi = fix(((rxhi*1.0-wo)/dw))
				endelse

			endif else begin ; not plotting wavelength or plotting many channels
			; implicit assumption in this is that max(x) = x(max)
			;				      min(x) = x(0)
			; This may not be true for some spectra?
			;
			print,'Enter plot range between',x(0),x(n_elements(x)-1) 
			print,'(Enter 0 and 0 to reset to all)'
			print,format='($,"Low ")'
			read,rxlow
			print,format='($,"High")'
			read,rxhi
			rxlow = float(rxlow); cast to float
			rxhi = float(rxhi)
			if (rxlow gt rxhi) then begin ; lowest -> rxlow
				temp = rxlow;
				rxlow = rxhi
				rxhi = temp
				endif		
		
			if ((rxlow lt x(0)) or (rxhi gt x(n_elements(x)-1))) then begin; bound check
				if ((rxlow eq 0) and (rxhi eq 0)) then begin
					print,'Resetting plot range'
					compbounds = 0
					endif else begin
			    		print,'Illegal choice - retry from menu'
		    	    		rxlow = 0&rxhi = 0
					compbounds = 0
	    		    		print,'press any key to continue. '
	    	    			print,' '
	    	    			dum = get_kbrd(1); halt to read error message
					endelse
	    	    		endif
      			endelse	;/* wavelength*/
		endelse; /*first_time*/
	end
   'Q': begin
	print,'ok - quitting'
        kktoggle = !d.name
        set_plot,'ps'
        device,/close
        set_plot,kktoggle
	true = 0
	end
   'q': begin
	print,'ok - quitting'
        kktoggle = !d.name
        set_plot,'ps'
        device,/close
        set_plot,kktoggle
	true = 0
	end

   '8':	begin
	if (theoryes eq 0) then begin
		theoryes = 1 
		print,'Overplotting theory spectra if they exist'
		print,'press any key to continue. '
		print,' '
		dum = get_kbrd(1)
		endif else theoryes = 0
	end
   '9': begin
	print,'Plotting Light Curve in ',orient,' orientation' 
	if (!d.name eq 'X') then begin
			winsav=!d.Window
			if (winsav eq -1) then winsav=0
			endif
	if (orient eq 'portrait') then begin
	  lcbsd,filename,/plot_page,/quiet,chan=chstring,spectra=spectstr,/stack
	 endif else begin
	  lcbsd,filename,/plot_page,/quiet,chan=chstring,spectra=spectstr,/lands,/stack
	 endelse
	if (!d.name eq 'X') then wset,winsav
	end 

   ELSE:begin

	first_time = 0; no longer true
        j = 0
;
; set plot range / clear plot range
;
if (compbounds eq 0) then kXrange = [rxlow,rxhi]
;
; start plotting loop
;
 	if (multiflag eq 1) then begin
		!p.multi = [0,2,2,0,0]
		chstart=1
		chend = 4
		!p.subtitle = ""
		endif else begin
		chstart = channel
		chend = channel
		!p.multi=0
		!p.subtitle = ""
		endelse

        nsummed = 1
        kplotdirn=1
	repeat begin 
          lskip =0;
	  tskip =0;
          if (j ge header.numchn(channel-1)) then lskip=1
	  if (lskip eq 0) then begin	
	  if (ntodo gt 1) then curspect = spectra(j) else curspect =fix(spectra,0)

;    rdbsd2data,filename,ptrdmp,numchn,curspect,channel,3,nbinsp,tsec79,td, $
;		sumdgi,crat,x,y,vy,wo,dw,ptoset,ptcorr,ochan,ier,uncspc
;rd_bsd_data,filename,header,rdmap,data

        rd_bsd_data,filename,[curspect],header,rdmap,index,data,[1,2,3,4]
        rd_spc_data,theorname,[curspect],bth_header,bth_rdmap,bth_index,bth_data,[1,2,3,4]

  
 	for chloop = chstart, chend do begin
	 tskip = 0
	 nsumspec = fix(nsumstr)	
         y = data(chloop-1,0).counts(0:index(chloop-1,0).nbinssp-2)

	  if (theoryes eq 1) then begin  
	   if (bth_rdmap(curspect-1).exist(chloop-1) eq 0) then tskip = 1 
	   endif else tskip = 1; no fit data
          


         if (multiflag eq 1) then begin
          if (chloop eq 1) then    $
  	   kptitle = 'Chan '+string(format='(i2)',chloop)+ $
		' Spect '+string(format='(i5)',curspect)+'   '+ $
                 'Sdgi '+string(format='(i4)',index(chloop-1,0).sdgi);

          if (chloop eq 2) then     $
           kptitle = 'Chan '+string(format='(i2)',chloop)+ $
                ' Time '+fmt_tim(index(chloop-1,0).time,/msec)

          if (chloop gt 2) then     $
           kptitle = 'Chan '+string(format='(i2)',chloop)

        endif else begin

	  kptitle = 'Chan '+string(format='(i2)',chloop)+ $
		' Spect '+string(format='(i5)',curspect)+'   '+ $
                 'Sdgi '+string(format='(i4)',index(chloop-1,0).sdgi) + $
                ' Time '+fmt_tim(index(chloop-1,0).time,/msec)
          endelse

          interval = index(chloop-1,0).sdgi
;	  kytitle = 'Counts/bin per '+string(format='(f4.1)',interval)+' secs'
	  if (waveyes eq 1) then kxtitle='Wavelength (angstroms)' else kxtitle='Bin no.'
	  kytitle = 'Counts per second'

  
          
	  if (compbounds eq 1) then begin
		; recompute bounds for wavelength
		 wlbound = (rxlow*index(chloop-1,0).dw)+index(chloop-1,0).wo
		 whbound = (rxhi*index(chloop-1,0).dw)+index(chloop-1,0).wo
		kXrange = [wlbound,whbound]
		endif

      	  if ((uncert eq 1) and (index(chloop-1,0).istat eq 0)) then begin
		print,'No uncertainties present -- hit a key to continue'
		dum = get_kbrd(1)
		uncert = 0;
		endif
	
	  if (waveyes eq 1) then begin
		px = data(chloop-1,0).wave(0:index(chloop-1,0).nbinssp-2)
		endif else begin
		px = data(chloop-1,0).bins(0:index(chloop-1,0).nbinssp-2)
		endelse


	  if (uncert eq 1) then begin
		top = y + vy
		bot = y - vy
                plot,px,top  ,Xrange=kxrange, font =-1, psym=10, $
                       xstyle=1,yrange=[0,0],linestyle=1,title=kptitle, $
                       xtitle=kxtitle,ytitle=kytitle
                oplot,px,y,psym=10,linestyle=0
                oplot,px,bot,psym=10,linestyle=1
		endif else begin
                plot,px,y  ,Xrange=kxrange, font =-1, psym=10, $
                       xstyle=1,yrange=[0,0],linestyle=0,title=kptitle, $
                       xtitle=kxtitle,ytitle=kytitle
	        endelse


	   if (tskip eq 0) then begin
		;stop
                ty = bth_data(chloop-1,0).counts(0:bth_index(chloop-1,0).nbinssp-2)
                tp = bth_data(chloop-1,0).primary(0:bth_index(chloop-1,0).nbinssp-2)
                oplot,px,ty,psym=0,linestyle=0
		if (bth_index(chloop-1,0).velfit ne 0) then begin
                  oplot,px,tp,psym=0,linestyle=1
                endif
           endif
	endfor

	  if (!d.name eq 'PS') then print,'Printed:',kptitle
	  if (khcflag eq 1) then begin
                    khcflag = 0
                    pprint
                    set_plot,kolddev
                    endif
		
	  if (ntodo gt 1)  then begin
            if (!d.name eq 'X') then print,kptitle
            if (no_stop eq 0) then dum = get_kbrd(1) else dum = get_kbrd(0);
                case strlowcase(dum) of 
                'c':  no_stop=1
                's':  no_stop=0
                'd':  kplotdirn=kplotdirn*(-1)
		'b':  j = j - 2*nsumspec*kplotdirn
		'!':  begin
			cmdstr = ''
			read,cmdstr
			dum = execute(cmdstr)
		        j = j - nsumspec*kplotdirn
			print,'+ [done]'
			;just in case ....
			set_plot,plot_unit
			end
		'm':  begin
			j = ntodo
                        kplotdirn=1
			end
                'q':  begin
			j = ntodo
			kplotdirn=1
			end
                'h': begin
                      khcflag = 1
                      kolddev = !d.name
                      set_plot,'ps'
		      j = j - nsumspec*kplotdirn
                     end
		'?': begin 
			print,'Plot Mode commands'
			print,' '
			print,'c for continous mode'
			print,'s for step by step'
			print,'q or m return to main menu '
			print,'d reverses time direction of plots'
			print,'b goes back one spectrum (depends on "d")'
			print,'h for hardcopy'
			print,'! for an IDL escape (can be used to call other programs)'
			print,'? for help'
			print,' '
		        j = j - nsumspec*kplotdirn
			end
                else:
                endcase
	     endif
	  endif; /*skip */
	  j = j + kplotdirn
          if (j lt 0) then j = ntodo
	  endrep until j ge ntodo

	if ((tekyes eq 1) and (ntodo eq 1)) then dum = get_kbrd(1);

       	end
   endcase
 endwhile
end