PRO plotbft,BSDNAM,channel=channel,PLOT_PAGE=plot_page, $
landscape=landscape,p2nd=p2nd,psym=psym,charsize=charsize
; Plot BCS cnts, Te,Td,EM,velocity
;;; 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