PRO plotbft,BSDNAM,channel=channel,PLOT_PAGE=plot_page, $
    landscape=landscape,p2nd=p2nd,psym=psym,charsize=charsize
; Plot BCS cnts, Te,Td,EM,velocity
;+
; NAME:
;	plotbft
; PURPOSE:
;	Plot results of BSDFIT for BCS analysis.
;
; CALLING SEQUENCE:
;	plotbft[,BSDNAM,chan='1,2,3',/plot_page,/landscape,/p2nd]
;
; INPUTS:
;	BSDNAM		The BSD filename in the form 'bsdyymmdd.hhmm'
;
; OPTIONAL INPUT PARAMETERS:
;	CHANNEL=[2,3]	Scalar or vector containing required channels
;
;	/PLOT_PAGE	The presence of this parameter causes the light curve,
;			electron temperature, and emission measure to be plotted
;			on one page, and turbulent velocity (from Te-Td),
;			upflow  emission measure and upflow velocity
;			time histories to be plotted on another page.
;
;       /landscape	This switch should produce lanscape orientation plots
;			when the postscript output option is selected.
;
;	/p2nd		Plot the ratio of the first and second component
;			emission measures.
;
; OUTPUTS:
;	Only plots are produced.
;
; COMMON BLOCKS:
;	None.
;
; SIDE EFFECTS:
;	1) The BSD, BTH and BFT files are openned and read.
;	2) Calls procedures  rd_bsd_header,rd_bsd_rdmap,rd_bsd_data and READBFT
;
; RESTRICTIONS:
;	BSDNAM must be a string-type scalar.
;
;	MEGAPLOT reads the BSD and FIT files. The BSD
;	and FIT files must exist on the same directory.
;
; PROCEDURE:
;	The following plots are generated in turn for each channel:
;
;	1) Light curve  	-  All other time plots are plotted 
;			    		on the same scale.
;	2) Electron temperature	- In M K.
;	3) Doppler temperature  - In M K. - only when single plots on page
;	4) Emission measure	- * 1.e50 cm-3.
;	5) Ion (turbulent) velocity - In km s-1  (from Te - Td).
;	6) Upflow emission meas - * 1.e50 cm-3
;	7) Uplow velocity       - in km s-1
;
; MODIFICATION HISTORY:
;	Nov, 1985,	Written, J.R. Lemen, MSSL
;	Mar, 1986,	RDB+JRL - Improved titles; made compatible with EOF/MSSL
;	May, 1986,	AF - change to use new two-component FIT files
;       jul, 1991,      atp - cut out some rubbish + transfered to idl2.
;       Nov, 1991       AF - made it work again
;	May, 1992  RDB	Cleaned up routine a lot; sped it up; got rid of
;			redundant stuff to do with "transfer" from SMM
;       Jun, 1992       Atp - finished, converted to structures, much updates
;				much yet to do.
;	july 1992 vn2.1 atp - cleaned up + finished (!)
;	28-Feb-92  RDB	Changed VEL to BCS_IONVEL
;	01-Mar-93  RDB	Added CLEARPLOT, CLEAR_UTPLOT at end
;	02-Mar-92  RDB	Corrected PS_LONG and /LAND interaction
;			Channel select now by integer
;			Added prompt for pause when to screen
;	18-May-94  RDB  Fixed problem in creating tdx, and with plots
;-

;;;			String containing channels to be processed, 
;;;		    '*' means all, or use a range '1-3', or explicitly 
;;;		    '1,3,4', default is '*'


;	Extract the time and counts for the light curve:
;
;  modify to use set_utplot
;
;  first_time - internal - time of first spectrum in theory.
;  bsd_index(chan,0).time - internal - time of first spectrum in bsd file
;   bsd_index.time - internal- time of each spectrum in bsd file
;   data1(2,*) - time in minutes from first_time for each theory spectrum    
;
;	structure of data1 array:
;		0=SPC_NUM	  1=CHAN	    2=Relative time (min).
;		3=SDGI( sec )	  4=DW (A/bin)	    5=Max # iterations
;		6=# iterats.	  7=# degree free.  8=Reduced chi-square
;		9=cnvrg flag	  10=velocity flag  11=lithium flag
;		12=Te_6	    	  13=Td_6(w)        14=Td_6(k or j)
;		15=Background	  16=EM (1.e50)     17=EM2 (1.e50) - second comp
;		18=upflow velocit 19=N(Li)/N
;		20 to 27 are uncertainties of 12 to 19.
;
; do time axes
; which is first, first_time or bsd_index(chan,0).time
; then work out parameters for set_utplot
;

version = 'Plotbft Program  (idl2)  V3.1 18/05/94'
;
print,version

ind_Td = 13				; Doppler temperature from line Ca w
					; and from Fe w as well - 26-May-86

;--	Set up Month string:

month = ['   ','Jan','Feb','Mar','Apr','May','Jun',	$
		       'Jul','Aug','Sep','Oct','Nov','Dec']

te_tit = ['  ','Fe XXVI  ','Fe XXV  ','Ca XIX  ','S XV  '] 


N_args = N_params(0)			; Find out how many arguments

;	Default to all channels if sel_chan is not specified
;;kchan = 0
;;if (keyword_set(channel)) then begin
;;		range,channel,4,channels,ier,1
;;		if (ier ne 0) then begin
;;			print,'Error in channel selection'
;;			return
;;			endif
;;		chan = channels(0)
;;		kchan = n_elements(channels)-1
;;		endif else begin
;;		kchan = 3
;;		channels = [1,2,3,4]
;;		chan = 1
;;		endelse
kchan = 1
chans = [3]
if keyword_set(channel) then begin
   chans = [channel]
   kchan = n_elements(chans)
endif
 
;
;		Set up file names:
;
if  N_elements( bsdnam ) eq 0  then  ask='Y' else $
   if strpos(bsdnam,' ') eq 0  then  ask='Y' else ask = 'N'
if ask eq 'Y' then bcs_getfile,filename,defstr="bsd*.*" else filename=bsdnam 
if (filename eq '-1') then return

;		Pick off file id:
break_file,strlowcase(filename),dsk_log,dire,bsdnam,ext,filver
file_id = bsdnam
if strpos(bsdnam,'bsd') eq 0 then $
   file_id = strtrim(strmid(bsdnam,3,strlen(bsdnam)),2)

fitnam = dsk_log+dire+'bft' + file_id + ext		;Create BFT file name
bsdnam = dsk_log+dire+'bsd' + file_id + ext		;Create BSD file name


print, 'Opening files:  ',bsdnam
print, '         and    ',fitnam

;
;		Read the BSD file
rd_bsd_header,bsdnam,bsd_header
rd_bsd_rdmap,bsdnam,bsd_header,bsd_rdmap

chan_idx    = where(bsd_header.numchn gt 0) + 1		;/* active channels */
chan_num    = n_elements(chan_idx)			;/* number of active channels */

rd_bsd_data,bsdnam,indgen(max(bsd_header.numchn))+1,bsd_header,bsd_rdmap, $
 bsd_index,/nodata	


;;tdx         = intarr(7,3,bsd_header.totspc)
;;tdx(*,1,*)  = bsd_index.time
tdx = fix(bsd_index.time)

multi = 0
;
;
;
;--	Read the BFT file
readbft, fitnam, first_time, data1	; Read the BFTxxxx.xxx files.
;
; convert first_time to external rep
;
   f_time = intarr(7)
   f_time(0) = first_time(0:2)
   f_time(3) = 0
   f_time(4) = first_time(3:5)
  first_time = f_time
;
;	only do ps_long if not landsacpe
if (!d.name eq 'PS') then begin
   if (keyword_set(landscape)) then device,/landscape  $
   else  ps_long   			; stretch the page
endif

;
; Set up plotting conventions
;
!noeras = 0				; enable erasing of screen
set_xy
if (not keyword_set(psym)) then psym = 1 
!p.font = -1
chars = !p.charsize
if keyword_set(charsize) then !p.charsize = charsize else !p.charsize=1.2
;
;!x.style = (!x.style and 23)		; draw a box around plot ???
;!y.style = (!y.style and 23)		; draw a box around plot ???
;!p.psym = 0				; Plot as connected line
;!p.font=0
;
;--		Begin the loop on channels
;
for j=0,kchan-1 do begin

  if keyword_set(plot_page) then begin
     multi = 3
     !p.multi = [0,1,3,0,0]	; three plots on page - this is fixed **
  endif else !p.multi = 0

;;  if j gt 0 then chan = j+1		; next time through - next channel
  chan = chans(j)


;--		Plot the light curve:

  vv = where( bsd_index.chan eq chan )
  if (vv(0) eq -1 ) then begin
     print,'** No data for this channel **',chan
     avail = all_vals(bsd_index.chan) & avail = avail(where(avail gt 0))
     print,'** Data only available for',avail
     goto, nodata
  endif
   
  print,'$(3x,a,i4,a,i2,a)', 'There are',n_elements( vv ),' channel', $
			chan,' spectra'

  if (!d.name ne 'PS') and (!d.name ne 'SIXEL') then $
	print,'Next plot/page?: ',format='(a,$)'	;message for pause

;
;	form the time for the banner
;       this needs fixing
;
  a1 = tdx(4,1,vv(0)) 
  a2 = tdx(5,1,vv(0))
  a3 = tdx(6,1,vv(0))
  a2 = month(a2)
  if a3 gt 1900 then a3 = a3-1900		;get year in 91 format
  tdate = string(a1,a2,a3,format='(i2,1h-,a,1h-,i2)')
;
;  work out offset between bsd data set and fit data set, label for utplot
;
  ktmp = int2secarr([[bsd_index(vv(0)).time],[first_time]])
;;  print,ktmp
  t_offset = ktmp(1)
  if (t_offset ge 0) then first_ptime=bsd_index(vv(0)).time $
	else first_ptime=first_time
;
; create xarray for utplot
;
  xtime = int2secarr(bsd_index(vv).time(*))
;
; make a secarr for the fit parameters
; 
  fv = where( data1(1,*) eq chan )
  if (fv(0) eq -1) then begin
    !noeras=0
    nofit=1
  endif else nofit = 0

 if (nofit eq 0) then begin
  fitime = data1(2,fv) ; time in minutes
  nfitpar = n_elements(fitime)
  fxtime = lonarr(nfitpar) + fitime*60.0 		;+ t_offset 
  endif
;
; light curve
;
  kytitle = 'Counts per sec'
  kptitle = string('BCS Chan #',chan,te_tit(chan),tdate,file_id,$
	"(a,i1,2x,a,3x,a,2x,'(',a,')')")

  if multi ne 0 then begin
    kytitle = 'CNTS/SEC'
    !p.linestyle = 0
    endif

  cnts = bsd_index(vv).crate
 
  utplot,xtime,cnts,first_ptime,ytitle=kytitle,title=kptitle

  xtit = 'Start Time ('+fmt_tim(first_ptime)+')'

 if (multi eq 0) then begin
  ans = '' 
  if (!d.name ne 'PS') and (!d.name ne 'SIXEL') then ans=get_kbrd(1)	; allow a pause
  endif 

if (nofit eq 0) then begin
 ;
 ;  Plot electron temperature:
 ;  
 Te_6 = data1( 12, * )
 Te_6 = Te_6( fv )			; Electron temperature in M K
 d_Te = data1( 21, * )
 d_Te = d_Te( fv )			; Uncertainties in electron temperature
 ;
 ;  Set up plot titles
 ;
 kytitle = 'Electron temperature!C(x10!E6!N K)'
 kptitle = string('BCS ',te_tit(chan),' electron temp ',tdate, $
 	file_id,"(3a,3x,a,2x,'(',a,')')")
 
 if multi ne 0 then begin
    kptitle = ' '
    kytitle='Te (x10!E6!N K)'
 endif
 ;
 ; do plot
 ;
 set_utplot,xrange=xtime
 plot,fxtime,Te_6,title=kptitle,ytitle=kytitle,xtitle=xtit
 errplot,fxtime,Te_6-d_Te,Te_6+d_Te
 ;
 if multi eq 0 then begin
  ans = '' 
  if (!d.name ne 'PS') and (!d.name ne 'SIXEL') then ans=get_kbrd(1)	; allow a pause
  endif 
 ;
 ; end electron temperature begin emission measure 
 ; Plot Emission (*1.e50) cm-3:
 ;
 EM_50 = data1( 16, * )		; Emission meas. (* 1.e50 cm-3)
 EM_50 = EM_50( fv )
 d_EM  = data1( 25, * )
 d_EM  = d_EM( fv )			; Uncertainties
 ;
 kytitle = 'Emission measure (x10!E50!N cm!E-3!N)'
 kptitle = string('BCS ',te_tit(chan),' EM * 1.e50 cm-3',tdate,	$
	file_id,"(3a,3x,a,2x,'(',a,')')")

 if multi ne 0 then begin
  kptitle = ' '
  kytitle='EM (x10!E50!N cm!E-3!N)'
  endif
 ;
 ; do plot
 ;
 set_utplot,xrange=xtime
 plot,fxtime,EM_50,ytitle=kytitle,title=kptitle,xtitle=xtit
 errplot,fxtime,EM_50-d_EM,EM_50+d_EM
 ;
 if (multi eq 0 or multi eq 3) then begin
  ans = '' 
  if (!d.name ne 'PS') and (!d.name ne 'SIXEL') then ans=get_kbrd(1)	; allow a pause
  endif 
 ;
 ; End emission measure begin turbulent velocity
 ; Plot Velocity (km s-1) from Te-Td
 ;
 ; get doppler temperatures
 ;
 Td_6 = data1( ind_Td, * )
 Td_6 = Td_6( fv )			; Doppler temperature in M K
 d_Td = data1( ind_Td + 9, * )
 d_Td = d_Td( fv )			; Uncertainties in doppler temperature
 ;
 ; work out turbulent velocities
 ;
 del_T = Td_6-Te_6
 ij = where( del_T lt 0. )
 if !err ne 0 then del_T(ij) = 0.	; If Td-Te is negative, set to 0.
 ;
 velo = bcs_ionvel(del_T,strmid(te_tit(chan),0,2))	; Calculate velocities
 d_te = d_te < 3e7
 d_td = d_td < 3e7
 d_vel = velo * 0.5 * (sqrt( double(d_Te)^2 + double(d_Td)^2 ) / (Td_6-Te_6))	; uncertainties
 ;
 kytitle = 'Turbulent velocity km per sec'
 kptitle = string('BCS ',te_tit(chan),' velocity',tdate,	$
 	file_id,"(3a,3x,a,2x,'(',a,')')")
 ;
 if multi ne 0 then begin
  ;the first plot on a new page:
  kptitle = string('BCS Chan #',chan,te_tit(chan),tdate,file_id,$
	"(a,i1,2x,a,3x,a,2x,'(',a,')')")
  kytitle='Turbulent velocity (km s!E-1!N)'
  endif
 ;
 ; do plot
 ; 
 set_utplot,xrange=xtime
 plot,fxtime,velo,ytitle=kytitle,title=kptitle,xtitle=xtit
 errplot,fxtime,velo-d_vel,velo+d_vel
 ;
 if multi eq 0 then begin
  ans = '' 
  if (!d.name ne 'PS') and (!d.name ne 'SIXEL') then ans=get_kbrd(1); allow a pause
  endif 
 ;
 ; end turbulent velocity 
 ; check second component:
 ;
 vel_flag = data1( 10,*)
 vel_flag = vel_flag( fv )		; velocity flag for this channel
 vnz = where(vel_flag gt 0)
 if (vnz(0) eq -1) then novel = 1 else novel = 0
 ;
 if (novel eq 0) then begin
  ;
  ; plot second component emission measure
  ;
  EM2_50 = data1(17,*)
  EM2_50 = EM2_50( fv )		; second component emission measure
  d_em2 = data1( 26 , * )
  d_em2 = d_em2( fv )		; uncertainty of second comp EM
  ;
  kytitle = 'Upflow emission measure (x10!E50!N cm!E-3!N)'
  kptitle = string('BCS ',te_tit(chan),' EM2 * 1.e50 cm-3',tdate,	$
 	file_id,"(3a,3x,a,2x,'(',a,')')")
  ;
  if multi ne 0 then begin
     kptitle = ' '
     kytitle='EM2 (x10!E50!N cm!E-3!N)'
     if keyword_set(p2nd) then !y.title='EM2/EM1'
     endif
  ;
  ;
  set_utplot,xrange=xtime
  if not keyword_set(p2nd) then begin
   plot,fxtime,EM2_50,ytitle=kytitle,title=kptitle,xtitle=xtit
   errplot,fxtime,EM2_50-d_EM2,EM2_50+d_EM2
  endif else begin
   plot,fxtime,em2_50/em_50,psym=psym,xtitle=xtit
   errplot,fxtime,(em2_50-d_em2)/em_50,(em2_50+d_em2)/em_50
  endelse

  ;
  if multi eq 0 then begin
   ans = '' 
   if (!d.name ne 'PS') and (!d.name ne 'SIXEL') then ans=get_kbrd(1)	; allow a pause
   endif 
  ;
  ; end upflow emission measure
  ; plot upflow velocity
  ;
  vup = data1( 18, * )
  vup = - vup( fv )		; upflow velocity ( > 0. )
  d_vup = data1( 27, *)
  d_vup = d_vup( fv )		; uncertainty of upflow velocity

  kytitle = 'Upflow velocity km per sec'
  kptitle = string('BCS ',te_tit(chan),' upflow velocity',tdate,	$
         file_id,"(3a,3x,a,2x,'(',a,')')")

   if multi ne 0 then begin
     kptitle = ' '
     kytitle='Upflow velocity (km s!E-1!N)'
     endif
  ;
  set_utplot,xrange=xtime
  plot,fxtime,vup,ytitle=kytitle,title=kptitle,xtitle=xtit
  errplot,fxtime,vup-d_vup,vup+d_vup
  ;
  ; end vel plot
  ; end
  ;
  endif ;/* no vel */
 endif ;/* no fit */
ans = '' 
if (!d.name ne 'PS') and (!d.name ne 'SIXEL') then ans=get_kbrd(1)	; allow a pause
;
nodata:
;
endfor	
;
!x.range = 0 

ps_reset   ; reset the page to the default length
!p.charsize = chars
clearplot
clear_utplot

return
end