PRO PLOTBPC,BPC_FILE,CHAN,OPTION,START  
;+
; NAME:
;   PLOTBPC
; PURPOSE:
;   Read the BPCxxxxxx.xxxx files written by BSDCAL and the BSD data file
;   specified within it and plot the spectral fits as given by parameters
;   in the BPC file over the raw data
;
; CALLING SEQUENCE:
;   PLOTBPC[,bpc_file,chan,option,start]
; 
;
; INPUTS:
;   BPC_FILE= Equal the BPC file name or I.D.
;        If BPC_FILE is undefined or the first character is
;        blank, PLOTBPC will prompt user for name.
;
;   CHAN    = required channel;
;
;   OPTION  = the fit option number required to display. If supplied
;             on the command line then only one option is allowed
;
;   START   = the time of first spectrum required ( default = all)
;
; COMMON BLOCKS:
;   None.
;
; SIDE EFFECTS:
;   None.
;
; MODIFICATION HISTORY:
; 0)  Written  May 1991   CDP   
; 1)  Revised to handle start time of spectra better   Sep 91  CDP
;
; 2)  Undo previous "bug-fix" and fix real one where rdmap variable
;     was incorrectly incremented between spectra      Nov 91  CDP
;
; 3)  Cut out bin-width scaling (not needed once BSDCAL sorted)  Nov 91 CDP
;
; 4)  Generalize and UNIXify handling of file names.   Mar 92 CDP
; 5)  Changed BGRND call to routine POLY.     June 92, RDB
; 6)  Use rd_bsd_head to access BSD file header     July  92   CDP
;
;--------------------------------------------------------------------------

   
; set up arrays required
;-----------------------

dummy    = lonarr(5)
rdummy   = fltarr(4)
time     = intarr(6)          ; spectrum start time (BPC file)
stime    = lonarr(7)          ; data start time (BSD file)
dw_etc   = fltarr(6)          ; dw, origwo, wo, wooff, sdgi, count rate
num      = intarr(2)          ; no. options, no. back
cs       = intarr(2)          ; channel & spectrum
bin      = intarr(2)          ; start & end bins
lparm    = fltarr(6)          ; line parameters
backparm = fltarr(6)          ; background parameters
optname  = ' '                ; character option name
data     = fltarr(26,200,10)  ; array for all information
                              ; dimensions are (parameter,spectrum,option)
sarray   = intarr(200)        ; array for spectrum numbers loaded
x        = fltarr(500)        ; array to hold bin values
y        = fltarr(500)        ; array to hold data values
temp     = fltarr(16)         ; fudge to get over read limitations
bcof     = fltarr(3)          ; background coefficients
vcof     = fltarr(9)          ; parameters for Voigt profile calculation
fwhml    = fltarr(4)          ; Lorentz width for each channel
wave     = fltarr(4)          ; wavelength of resonance line
ewid     = fltarr(4)          ; electronic resolution of each channel
tmp_val  = fltarr(6)          ; yet another dummy
xann     = fltarr(2)          ; x coordinates of annotation
yann     = fltarr(2)          ; y coordinates of annotation
maxopt   = intarr(4)          ; max number options in each channel
num_opt  = intarr(10)         ; array for input of options wanted
tchar    = '        '         ; character variable for input time
charchan = ' '                ; string for channel number input
ans      = ' '                ; one letter answer from terminal

; initialise integers used in FOR loops (version 2 much fussier than v 1!)
; ------------------------------------------------------------------------
nreq_opt = 0
lines    = 0
kspec    = 0
nrecs    = 0
kopt     = 0
nbinsp   = 0


; load known values 
;------------------
fwhml  = [ 2.08e-4, 2.06e-4, 3.89e-4, 13.71e-4 ]    ; Lorentzian width
ewid   = [ 0.0, 0.0, 0.0, 0.0 ]                     ; electronic width
wave   = [ 1.778, 1.8509, 3.1781, 5.0385 ]          ; wavelength resonance
maxopt = [5, 9, 3, 7]                               ; max number options 
                                                    ; in each channel
clear_screen = '\033[;H\033[2J'                     ; clear page
alpha        = '\x18'                               ; reset to alpha screen
bell         = '\007'                               ; warning
true         = 1                                    ; psuedo logical
false        = 0                                    ;  ditto
nparm = n_params()                                  ; set variable so can
                                                    ; be changed to loop

;
;  compile to get message away from plot
;
dumx = calfun()
;;dumx = bgrnd()

if nparm lt 3 then option = 0                  ; if not on line, initialise
if nparm eq 0 then bpc_file = ' '              ; if not even the file 
                                               ; name given

;print, clear_screen                             ; clear screen

if (strmid(bpc_file,0,1) eq ' ') or $              ; if not valid BPC file
   (strlen(bpc_file) eq 0) then begin              ; name then manual input
   print,'$(a)', ' Enter the BPC file name (omit BPC prefix): ' 
   bpc_file = ' '  &  read,bpc_file
endif
if strlen(bpc_file) eq 0 then return               ; trap exit

if strpos(strupcase(bpc_file),'bpc') lt 0 then $    ; if BPC prefix not
   bpc_file = 'bpc'+bpc_file                        ; there, affix it

;------------------------------------------------------------------------
;  procedure loops back to here when request satisfied or user Quit     -
;  requested during plotting                                            -
;------------------------------------------------------------------------

continue:
data(*,*,*)  = 0                                    ; for safety
get_lun, lun                                        ; get free unit number 
                                                    ; to attach file
openr, lun, bpc_file                                ; open for readonly


bsd_name = ' '                                ; get the bsd file name from 
readf,lun,bsd_name                            ; the top line of the BPC   
bsd_name = strmid(bsd_name,$                  ; file
           strpos(bsd_name,'BSD file')+8,$
           strlen(bsd_name))
break_file,bsd_name,disk_log,dir_log,fil,ext,fver
bsd_name = strtrim(fil + ext,2)

buf = ' '                                     ; get the BPC file format   
readf,lun,buf                                 ; code
mode = fix(strtrim(strmid(buf,$
           strpos(buf,'mode:')+5,$
           strlen(buf)),2))

if mode eq 0 then begin                       ; if not 1, then error

   tex = ' Expected BPC format mode = 1.'$
         + ' Program aborted.'
   print, bell
   print,'$(/a)',tex
   return
endif

;-------------------------------------------------------------------
; read header of BSD file to get start time     
;-------------------------------------------------------------------
rd_bsd_header,bsd_name,bsd_h
;
; calculate day of the mission
;
dom,bsd_h.tdf(4),bsd_h.tdf(5),bsd_h.tdf(6),day_of_mission   

; time of first spectrum in mission seconds

bsd_stime = long(day_of_mission)*86400 + $            
            long(bsd_h.tdf(0))*3600 + bsd_h.tdf(1)*60 + bsd_h.tdf(2)    

;-------------------------------------------------------------------
;  get channel wanted
;-------------------------------------------------------------------

err_flag = 1                                  ; force at least one go
while err_flag eq 1 do begin
   if (nparm lt 2) then begin                 ; channel not supplied
      print,'$(a)', ' Enter the channel number: (Q to quit)'
      read,charchan                                  ; character input
      if strupcase(charchan) eq 'Q' then begin       ; null response quits
         close, lun
         return
      endif else begin
         chan = fix(charchan)                        ; otherwise translate
      endelse
   endif
   err_flag = 0                                ; reset flag
   if (chan lt 1 or chan gt 4) then begin      ; set sensible limits
      err_flag = 1
   endif
endwhile

;-----------------------------------------------------------------------
;  get options list wanted
;-----------------------------------------------------------------------

err_flag = 1                                    ; force at least one go
while err_flag eq 1 do begin
   if (nparm lt 3) or (option eq 0) then begin  ; option not on call line
      nreq_opt = fix(0)
      free_in, chan, num_opt, nreq_opt          ; free format input
   endif else begin
      nreq_opt = 1
      num_opt(0) = option
   endelse
   err_flag = 0                                         ; reset error flag
   for i = 1,nreq_opt do begin                          ; check each is OK
      if (num_opt(i-1) lt 1) or $
         (num_opt(i-1) gt maxopt(chan-1)) then begin    ; limit errors
         err_flag = 1              
         option   = 0                                   ; error, set flag
      endif
   endfor
endwhile

;-----------------------------------------------------------------------
;  get time of first spectrum wanted for plotting
;-----------------------------------------------------------------------

colpos1 = -1
colpos2 = -1

if (nparm lt 4) then begin                    ; start time not on call line

   while colpos1 lt 0 do begin                 ; so get it here
      print,'$(a,i2,a,i2,a,i2,a,i2,a,i2,a,i4)',$
      ' First spectrum has start time (hh:mm:ss) of ',bsd_h.tdf(0),':',$
        bsd_h.tdf(1),':',bsd_h.tdf(2),' on date ',bsd_h.tdf(4),'/',$
        bsd_h.tdf(5),'/',bsd_h.tdf(6) 
      print,'$(a)',$
      ' Enter time of the first spectrum required (hh:mm[:ss])',$
      ' - default gives all -'
      read,tchar
      if strlen(tchar) eq 0 then begin
         tchar = string(bsd_h.tdf(0),':',bsd_h.tdf(1),':',bsd_h.tdf(2),$
         "(i2,a,i2,a,i2)")
      endif
      colpos1 = strpos(tchar,':')
   endwhile
   colpos2 = strpos(tchar,':',colpos1+1)      ; search for second colon for
                                              ; seconds field

endif

thr     = fix(strmid(tchar,0,colpos1))           ; extract hours
tmin    = fix(strmid(tchar,colpos1+1,2))         ; extract minutes
if colpos2 gt -1 then begin
   tsec = fix(strmid(tchar,colpos2+1,2))         ; extract seconds if there
endif else begin
   tsec = 0
endelse
start   = long(thr)*3600 + tmin*60 + tsec             ; start time in seconds
start   = start + long(day_of_mission)*86400          ; add start day

if start lt bsd_stime then begin                ; if requested start time
   start = start + 86400                        ; is less than first spectrum
endif                                           ; assume next day is wanted

print,'$(///2a)',' Spectral data from file: ',bsd_name
print,'$(2a)',' Reading results file ',strupcase(bpc_file)
print,'$(a20,i2,a10,10i2)',' Looking for channel ',chan,$
                           ' option #s',num_opt(0:nreq_opt-1)


erase                                            ; clear graphics screen
print, alpha                               ; return to alpha

readf,lun,buf                      ; read 1 more header record of ----'s  


kspec      = -1                          ; initialize spectrum total counter
sp_skip    = 0                          ; number of spectra time-skipped
inform_opt = 0                         ; has missing options been notified?

while not eof(lun) do begin         ; loop while not at end of file

   for i=0,1 do begin
      readf,lun,buf                 ; skip 2 lines below the -----'s
   endfor
   newdata = true

   readf,lun,cs,time,dw_etc,num     ; read header line

   dom,time(3),time(4),time(5),day_of_mission  ; get day of mission of spectrum

   tsec = long(time(0))*3600 + time(1)*60 + time(2)$
          + long(day_of_mission)*86400               ; time in seconds

   if cs(0) eq chan and (tsec ge start) then begin      ; is it the desired 
                                                        ; channel and is it
                                                        ; after specified
                                                        ; start time?

      for i=0,2 do begin
         readf,lun,buf              ; skip 3 lines of headers etc
      endfor

      kopt = -1                                    ; running total of options
                                                   ; found

      for iopt = 1,num(0) do begin              ; for all the options 
                                                ; in this channel
         extra = false                          ; not yet read extra lines
         bin = [0,0]
         lparm = [0.0,0.0,0.0,0.0,0.0,0.0]      ; zero out arrays
         backparm = [0.0,0.0,0.0,0.0,0.0,0.0]   

         tempback = fltarr(num(1)*2)            ; create temp array for backg

         readf,lun,opt,okflag,lines,bin,rchi,$
                                      lparm,tempback,optname ; read data line
         backparm(0:num(1)*2-1) = tempback                   ; store backgnd


         for ireq=0,nreq_opt-1 do begin         ; loop over requested options
            if opt eq num_opt(ireq) then begin  ; is this a correct option?   
                                                
               mode          = 1                ; set normal mode

               if newdata then begin
                  kspec      = kspec+1             ; increment spectrum total
                  newdata    = false               ; start of a new parameter
               endif

               kopt               = kopt + 1            ; options stored
               data(0,kspec,kopt) = cs(1)               ; save the spectrum number
   
               data(1,kspec,kopt) = bin(0)              ; save the start bin number
               data(2,kspec,kopt) = bin(1)              ; save the end bin number  
               data(3:8,kspec,kopt) = time(0:5)         ; save the time
   
               data(9, kspec,kopt)  = backparm(0)       ; save background constant
               data(10,kspec,kopt)  = backparm(2)       ; save background slope
               data(11,kspec,kopt)  = backparm(4)       ; save background curvature
                                          
               data(12,kspec,kopt)  = lparm(0)          ; save first line's height
               data(13,kspec,kopt)  = lparm(2)          ; save first line's width
               data(14,kspec,kopt)  = lparm(4)          ; save first line's position
             
               if lines gt 1 then begin             ; if there are more lines
                                                    ; in this fit
                  extra = true                      ; flag: have read them
                  for i=0,lines-2 do begin
                     if opt eq 2 then begin           ; is this opt 2 (vel)?
                        mode = 2
                        readf,lun,tmp_val              ; read additional lines 
                        data(15+i*3,kspec,kopt) = tmp_val(0) ; save next line's intensity
                        data(16+i*3,kspec,kopt) = tmp_val(4) ; save next line's vel
                        data(17+i*3,kspec,kopt) = 0.0      ; spare
                     endif else begin
                        readf,lun,lparm               ; read additional lines  
                        data(15+i*3,kspec,kopt) = lparm(0) ; save next line's height
                        data(16+i*3,kspec,kopt) = lparm(2) ; save next line's width
                        data(17+i*3,kspec,kopt) = lparm(4) ; save next line's posn
                     endelse
                  endfor
               endif
               data(21,kspec,kopt)  = mode                 ; save the mode for
                                                           ; this option
               data(22,kspec,kopt)  = lines                ; save # of lines
               ireq = fix(nreq_opt - 1)                    ; skip the rest
               data(23,kspec,0)     = kopt         ; save # options in spectrum
               data(24,kspec,kopt)  = okflag     ; was fit OK?
               data(25,kspec,kopt)  = dw_etc(0)  ; save dw for this one
            endif 

         endfor
         if not extra then begin
            for i=0,lines-2 do begin
               readf, lun,buf                 ; if necessary, skip other
            endfor                            ; lines
         endif
      endfor        
      
      if inform_opt eq 0 then begin
         if (kopt+1) lt nreq_opt then begin
            print, bell, ' NOTE: not all the options requested were found.'
            inform_opt = 1
         endif
      endif
         

      readf,lun,buf                             ; end of spectrum options so
                                                ; skip the ---'s line
   endif else begin
                                                ; not channel wanted or 
                                                ; spectrum too early

      if cs(0) eq chan and (tsec lt start) then begin ; count skipped spectra
         sp_skip = sp_skip + 1
      endif

      buf = ' ' 
      minibuf = ' '                             ; so skip till ------'s line

      while minibuf ne '----------' and not eof(lun) do begin
         readf,lun,buf                          ; search for next line of
         minibuf = strmid(buf,2,10)             ; --------'s
      endwhile
   endelse
endwhile

if sp_skip gt 0 then begin                      ; inform of spectra skipped
   if sp_skip eq 1 then begin
      print, ' One spectrum skipped prior to start time.', bell 
   endif else begin
      print, sp_skip,' spectra skipped prior to start time.', bell 
   endelse
endif

if kspec eq -1 then begin                       ; inform if no spectra found
   print, bell, ' ** ERROR ** no spectra for that'
   print, ' channel and options-list located.', bell
   return
endif else begin
   if kspec eq 0 then begin
      print, ' Only one spectrum loaded.'       ; tidy in case only one
   endif else begin
      print, kspec+1,' spectra loaded.'         ; inform of success rate
   endelse                                      
endelse                                         


wait,3.0                                        ; wait for user digestion

close, lun                                      ; finished with BPC file
free_lun, lun                                   ; release unit
                                                
sarray = data(0,*,0)                            ; extract spectrum #s
ssarray  = intarr(kspec+1)                      ; sorted spectrum numbers
ssarray = sort(sarray(0:kspec))                 ; do the sort
;-----------------------------------------------------------------------
;   now access the BSD data file       
;-----------------------------------------------------------------------

get_lun, lun                 ; get a unit for the BSD data file
openr, lun, bsd_name         ; unformatted, readonly

numchan = lonarr(5)          ; total spectra and number per channel
ptrdmp  = lonarr(4)          ; roadmap pointers for each channel

readu, lun,numchan,ptrdmp    ; read first header record in file
numspec = numchan(chan)      ; total number of spectra in this channel
rdmap   = ptrdmp(chan-1)     ; the record where roadmap for this channel is
rdmap = (rdmap-1) * 64       ; byte location of start of roadmap

if kspec+1 gt numspec then begin      ; check number of spectra OK
   print, bell, ' ** ERRROR ** number of spectra referenced in BPC file'
   print, ' is greater than those present in BSD file', bell
   return
endif
   
for is=0,kspec do begin                 ; loop over spectra to be plotted
   sis = ssarray(is)                     ; order as sorted above
   numspec = fix( data(0,sis,0) ) - 1    ; spectrum to plot (zero start)
   nrecs = fix(numspec/16)        ; there are 16 spectral pointers per record
   point_lun, lun,rdmap+nrecs*64    ; go to where required spectrum pointer
                                    ; is located, ready for reading
   
   pointers  = lonarr(16)       ; temp storage for pointers
   readu, lun,pointers          ; read them
   
   spcrec = pointers(numspec-nrecs*16)   ; finally the record where 
                                        ; required spectrum is located
   spcbyt = (spcrec-1)*64        ; in bytes to beginning of record wanted
   point_lun, lun,spcbyt        ; go there ready to read
   
   readu, lun,pointers          ; read first housekeeping  record

   nbinsp = pointers(14)        ; to get the number bins per spectrum
   readu, lun,pointers          ; skip second housekeeping record
   
   nrecs = nbinsp/16            ; no. of complete records used by bins etc
   nleft = nbinsp - nrecs*16    ; any leftovers from complete records
   
   for i=1,nrecs do begin       ; read the bin values
      readu, lun,temp           ; one record at a time
      start = (i-1)*16          ; where to start loading temp into x
      x(start:start+15) = temp  ; load x
   endfor
   
   if nleft gt 0 then begin     ; have leftovers to deal with  
      readu, lun,temp           ; read 'em
      start = nrecs*16          ; where to load into x
      x(start:start+nleft-1) = temp(0:nleft-1)  ; load 'em
   endif
   
   for i=1,nrecs do begin       ; read the spectral data values
      readu, lun,temp           ; one record at a time
      start = (i-1)*16          ; where to start loading temp into y
      y(start:start+15) = temp  ; load y
   endfor
   
   if nleft gt 0 then begin     ; have leftovers to deal with  
      readu, lun,temp           ; read 'em
      start = nrecs*16          ; where to load into y
      y(start:start+nleft-1) = temp(0:nleft-1)  ; load 'em
   endif
   
;-----------------------------------------------------------------
;  attempt to plot data and overplot the fit
;-----------------------------------------------------------------

   kopt = fix(data(23,sis,0))                         ; max # of options 
                                                     ; in this spectrum
   minbin = 255                                      ; inverse limits to
   maxbin = 0                                        ; find real limits
   for i=0,kopt do begin
      if data(1,sis,i) lt minbin then begin           ; find min 
         minbin = data(1,sis,i)
      endif
      if data(2,sis,i) gt maxbin then begin           ; and max
         maxbin = data(2,sis,i)                       ; bins to be used
      endif
   endfor

   xs = minbin-1
   xe = maxbin-1                                     ; set data limits
                                                     ; -1 to account for
                                                     ; zero start arrays


;
;  not my idea! repeat this bit if want hardcopy after seeing terminal
;  version
;
   hardcopy:


;  set up plot title and information lines
;-----------------------------------------

   if nreq_opt gt 1 then begin                           ; more than 1 option?
      opt_str = string(' Options: ',num_opt(0),"(a,i1)") ; start string off
      for i=2,nreq_opt do begin
         temp_str = string(', ',num_opt(i-1),"(a,i1)")  ; next option
         opt_str = opt_str + temp_str                        ; add it
      endfor
   endif else begin
      opt_str = string(' Option: ',num_opt(0),"(a,i1)")     ; only one option
   endelse


   plt_title = string('BSDCAL results for Channel: ',chan,' Spectrum: ',$
                       numspec+1,opt_str,"(a,i1,a,i3,a)")
   time_tit  = string('Time: ',data(3,sis,0),'h ',data(4,sis,0),'m ',$
                      data(5,sis,0),'s on ',data(6,sis,0),'/',data(7,sis,0),$
                      '/',data(8,sis,0),"(a,i2,a,i2,a,i2,a,i2,a,i2,a,i4)")
; plot raw data
;--------------

   !p.linestyle = 0                                               ; solid line
   plot,x(xs:xe),y(xs:xe),title=plt_title,$
        xtitle='Bin value',ytitle='Counts/sec/unit bin position' ; plot data
   xwid = (!x.crange(1) - !x.crange(0))*0.05                     ; plot width
   xyouts, !x.crange(0)+xwid,0.95*!y.crange(1),bsd_name         ; bsd file
   xyouts, !x.crange(0)+xwid,0.9*!y.crange(1),time_tit         ; time info

   xann(0) = !x.crange(0)+xwid                                 ; legend    
   xann(1) = !x.crange(0)+xwid*2                               ; x posn
   yann(0) = 0.8*!y.crange(1)                                  ; and
   yann(1) = 0.8*!y.crange(1)                                  ; y posn
   oplot, xann,yann                                            ; output line
   xyouts, xann(1)+0.1*xwid,0.99*yann(1),'raw data'            ; annotate

; calculate the fit for each option used   
;---------------------------------------
   for k=0,kopt do begin
      lines   = fix(data(22,sis,k))           ; lines in this option
      bcof    = data(9:11,sis,k)              ; load background coefficients
      okflag  = data(24,sis,k)                ; recall fit flag
      this_dw = data(25,sis,k)                ; recall dw for this fit

      for i=1,lines do begin                         ; the coefficient vector
         vcof((i-1)*3)   = data(12+(i-1)*3,sis,k)     ; height
         vcof((i-1)*3+1) = data(13+(i-1)*3,sis,k)     ; width
         vcof((i-1)*3+2) = data(14+(i-1)*3,sis,k)     ; position
      endfor
      pxmin = data(1,sis,k) - 1                      ; xmin for this option
      pxmax = data(2,sis,k) - 1                      ; xmax for this option
      if okflag gt 0 then begin         
         !p.linestyle = 2                           ; dashed line for OK fit
      endif else begin
         !p.linestyle = 1                           ; dotted for bad fit
      endelse

      if vcof(0) gt 0.0 then begin               ; check sensible               
         mode  = data(21,sis,k)                      ; set mode for this option
         yfit  = calfun(x,vcof,lines,fwhml(chan-1),$
                           this_dw,mode,wave(chan-1),ewid(chan-1))$
                           + poly(x,bcof)	;;bgrnd(x,bcof)
   
         oplot, x(pxmin:pxmax),yfit(pxmin:pxmax)       ; overlay calculated fit
      endif
      if(k eq 0) then begin
         yann(0) = 0.75*!y.crange(1)                      ; legend
         yann(1) = 0.75*!y.crange(1)                        ; y posn
         oplot, xann,yann                                     ; output line
         !p.linestyle = 0                                          ; solid line
         xyouts, xann(1)+0.1*xwid,0.99*yann(1),'calculated fit'   ; annotate
      endif
      if okflag eq 0 then begin
         yann(0) = 0.65*!y.crange(1)                      ; legend
         yann(1) = 0.65*!y.crange(1)                        ; y posn
         oplot, xann,yann                                     ; output line
         !p.linestyle = 0                                          ; solid line
         xyouts, xann(1)+0.1*xwid,0.99*yann(1),'fit not converged'
      endif
   ;
   ; plot the background used
   ;-------------------------
   
      yfit = poly(x,bcof)		;;bgrnd(x,bcof)
      !p.linestyle = 1
      oplot, x(pxmin:pxmax),yfit(pxmin:pxmax)
      if k eq 0 then begin
         yann(0) = 0.7*!y.crange(1)                               ; legend
         yann(1) = 0.7*!y.crange(1)                               ; y posn
         oplot, xann,yann                                         ; output line
         !p.linestyle = 0                                            ; solid line
         xyouts, xann(1)+0.1*xwid,0.99*yann(1),'background used'  ; annotate
      endif

   endfor          ;  end of options in this spectrum

 
   !p.linestyle = 0                                           ; solid line type
   if ans eq 'H' then begin                                ; back to terminal
      set_plot,save_dev
      ans = ' '
   endif else begin
      xyouts, 0.0,0.01,$
         ' Hit any key to proceed: (Q to quit: H for hard copy)',/normal 
      ans = get_kbrd(1)                                    ; wait for reply
   endelse

   ans = strupcase(ans)
   if ans eq 'Q' then begin
      is = kspec                                       ; kill loop by 
   endif                                               ; resetting counter

   if ans eq 'H' then begin                            ; like it! want 
      xyouts, 0.15,0.5 ,$                              ; hard copy
         ' HARD COPY being prepared in PostScript file (IDL.PS)',/normal
      save_dev = !d.name
      set_plot,'PS'
      goto, hardcopy
   endif

endfor

close, lun                   ; close BSD file
free_lun,lun                ; free unit
nparm = 0                  ; kid it that it needs channel number etc.
print, alpha              ; switch back to alpha screen
print, ' '
goto, continue           ; naughty but nice!

end