;+
; NAME:
;       bcs_emiss   (was bcs_counts)
; PURPOSE:
;       This routine plots the following graphs for the resonance 
;       lines of the Yohkoh BCS.
;         1) Temp vs G(Te) functions
;         2) Temp vs Fluxes
;         3) Temp vs Total counts
;         4) Temp vs Counts/S in peak bin. 
; CALLING SEQUENCE:
;       bcs_emiss
; INPUTS:
;       G(t) fundtions taken from $DIR_BCS_CAL/emmsbcs.xdr
; OPTIONAL INPUTS:
; KEYWORD INPUTS:
;       channel         Array of channel numbers  (Default is [1,2,3,4])
;                       Note: 1=FeXXVI, 2=FeXXV, 3=CaXIX, 4=SXV
;       abundance       String, or array.
;                       If array must be 4 elements long e.g.
;                          bcs_emiss,abun=[4.47e-5,4.47e-5,4.47e-6,8.13e-6]
;                       If string must be one of: 
;                          'GAD'  'MEWE'  'XRP'  'MEYER'  (def=MEYER)
;       emeasure        Specify value to be used  (def = 1.e50 cm-3)
;       summary         Creates a summary file called BCS_EMISS.OUT in
;                       current directory
; OUTPUTS:
;       Plots to current graphics device
; OPTIONAL OUTPUTS:
; RESTRICTIONS:
; PROCEDURE:
;       Read G(t) functions from file, calculate plots
; MODIFICATION HISTORY:
;       RDB  24-Sep-88  Written
;       ATP  16-Jul-93  Rewritten for IDL 2
;       RDB  18-Aug-94  Ported for Unix
;       RDB  26-Aug-94  Allowed EM other than 1.e50 (emeasure keyword)
;                       read instrument stuff from latest cal files
;                       Some mods to plots
;-

pro bcs_emiss,channel=channel,abundance=abundance,emeasure=emeasure, $
        summary=summary

;               BCS channels
elstr=['Fe XXVI','Fe XXV ','Ca XIX ','S XV   ']

;       Some atomic physics parameters
w_lam=[1.778,1.8508,3.1781,5.0382]              ;A

;       masses
z=[55.85,55.85,40.08,32.06]
;
;       some instrument parameters - The sensitivites come from
;       	$DIR_BCS_CAL:bcs*.00n
;
;sens=[0.145,0.096,0.198,0.037]                  ;cm2
;bin_wid=[0.253,0.355,0.146,0.497]*1.e-3         ;mA
sens = fltarr(4) & bin_wid=fltarr(4)
for j=1,4 do begin
  rd_bcscal,j,efarea=efarea,disp=disp
  sens(j-1)    = efarea				;cm2
  bin_wid(j-1) = abs(disp)			;mA                  
endfor

;
;       process options
;  
    if keyword_set(channel) then begin
        chan = fix(channel)
        nchans = n_elements(channel)
        endif else begin
        chan = [1,2,3,4]
        nchans = 4
        endelse

    if keyword_set(abundance) then begin
        fred=size(abundance)
        stype= fred(fred(0)+1)
        if stype eq 7 then begin
          case strupcase(abundance) of
          'GAD': begin
                abun_desc = 'GAD'
                abun=[4.47e-5,4.47e-5,4.47e-6,8.13e-6]
                end
          'MEWE': begin
                abun_desc = 'MEWE'
                abun=[3.98e-5,3.98e-5,2.00e-6,1.58e-5]
                end
          'XRP': begin
                abun_desc = 'JRL/XRP'
                abun=[3.98e-5,3.98e-5,3.98e-6,1.58e-5]
                end
          'MEYER': begin
                abun_desc = 'Meyer'
                abun=[3.92e-5,3.92e-5,2.94e-6,8.63e-6]
                end
          default: begin
                print,' '
                print,'Abundance name not recognised - using Meyer'
                abun_desc = 'Meyer'
                abun=[3.92e-5,3.92e-5,2.94e-6,8.63e-6]
                end
           endcase
        endif else begin
        abun_desc = 'Custom'
        abun = abundance
        endelse
      endif else begin
        print,'Using default abundance set  (Meyer)'
        abun_desc = 'Meyer'
        abun=[3.92e-5,3.92e-5,2.94e-6,8.63e-6]
      endelse

;	default emission measure if 10^50; allow other values
  emrq = 50
  if keyword_set(emeasure) then begin
    emrq = fix(emeasure)
  endif
  emstr = '(EM = 1.e' + string(emrq,format='(i2)') + ' cm-3)'
  print,emstr


datafilestring = concat_dir('$DIR_BCS_CAL','emmsbcs.xdr')
;;;  bcs_reademm2,datafilestring,data,text,elem

print,'Restoring data from: ',datafilestring
restore,datafilestring 
help,data,text,elem


;               Create storage arrays
;
tzaa=fltarr(150,4)
emfun=fltarr(150,4)
flux=fltarr(150,4)
counts=fltarr(150,4)
peaks=fltarr(150,4)
ntmp=n_elements(tz)-1


; Set up plot variables:

  ss=size(data) & num_cases=ss(2) & nmax=ss(1)

;	Get min and max temperatures:

;;  plot_xrange = fltarr(2) 
;; 
;;  plot_xrange(0) = min(10.^data(1,*)) 
;;  plot_xrange(1) = max(10.^(data(1,*)+(data(3,*)-1)*data(2,*)))
;;print,plot_xrange
plot_xrange = [1.e6,1.e8]

;;; Get min and max fluxes:
;;
;;  plot_yrange = fltarr(2)
;;
;;  npts = 8 + data(3,0)
;;  plot_yrange(0) = min(data(9:npts,0)) 
;;  plot_yrange(1) = max(data(9:npts,0))

;;  for ii=0,nchans-1 do begin
;;    k = chan(ii) - 1
;;    npts = 8 + data(3,k)
;;    plot_yrange(0) = plot_yrange(0) < min(data(9:npts,k)) 
;;    plot_yrange(1) = plot_yrange(1) > max(data(9:npts,k))
;;  endfor


  for ii=0,nchans-1 do begin
        
        k = chan(ii) - 1

        print,'$(/x,a,i4,4x,a)','Line',k,elstr(k)
;        print,data(0:8,k)

        tz1=(data(1,k)+data(2,k)+indgen((data(3,k)-1)*5.)*data(2,k)/5.)
;
;       Calculate temperature array
;
        Temp  = (data(1,k) + indgen(data(3,k))*data(2,k))
;
;       that corresponds to the emission function
;
        npts = 8 + data(3,k)
        emits = data(9:npts,k)
;
;       interpolate using spline
;
	emfact = 10.^(emrq-50)      
        print,'EM factor =',emfact
        emf=spline(temp,emits,tz1) * emfact	
        tz=10.^tz1                      ;expand interpolated temps

        flx=emf*abun(k)
        cps=flx*sens(k)

        lams=7.17e-7*w_lam(k)*sqrt(tz/z(k))
        peak=cps*bin_wid(k)/(sqrt(!pi)*lams)

        print,elstr(k),abun(k),w_lam(k),z(k),sens(k),bin_wid(k)

        insert,emfun,emf,0,k
        insert,flux,flx,0,k
        insert,counts,cps,0,k
        insert,tzaa,tz,0,k
        insert,peaks,peak,0,k
        endfor

print,'$(//x,14x,4(5x,a),/)',elstr
print,'$(x,a,4(e12.2))','Abundance     ',abun
print,'$(x,a,4f12.4)',  'Lamda w (A)   ',w_lam
print,'$(x,a,4f12.2)',  'Atomic weight ',z
print,'$(x,a,4f12.4)',  'Sensitivity   ',sens
print,'$(x,a,4f12.3)',  'Bin Width (mA)',bin_wid/1.e-3
print,' '
print,' Abundance Source:  ',abun_desc
print,' Instrument data from calibration files'
print,' Emission measure used:  1.e'+string(emrq,format='(i2)')+' cm-3'
print,' '
print,' '

if keyword_set(summary) then begin
  openw,lunw,'bcs_emiss.out',/get_lun
  printf,lunw,'$(//x,14x,4(5x,a),/)',elstr
  printf,lunw,'$(x,a,4(e12.2))','Abundance     ',abun
  printf,lunw,'$(x,a,4f12.4)','Lamda w (A)   ',w_lam
  printf,lunw,'$(x,a,4f12.2)','Atomic weight ',z
  printf,lunw,'$(x,a,4f12.4)','Sensitivity   ',sens
  printf,lunw,'$(x,a,4f12.3)','Bin Width (mA)',bin_wid/1.e-3
  printf,lunw,' '
  printf,lunw,' Abundance Source:  ',abun_desc
  printf,lunw,' Instrument data from calibration files'
  printf,lunw,' Emission measure used:  1.e'+string(emrq,format='(i2)')+' cm-3'
  printf,lunw,' '
  free_lun,lunw
endif

;
;-- Plot the lines
;
!p.linestyle=0
!noeras = 0
;


plot_xtitle='Electron Temperature'
plot_title='Yohkoh BCS G(Te) functions    ' + emstr    ;(EM = 1.e50 cm-3)'
plot_ytitle='G(Te)/Abundance'


;set_xy,1.e+6,9.e7,1.05e7,9.e11
    

;    plot_yrange = [(min(emfun(*,chan-1))>1),max(emfun(*,chan-1))]
;pmm,emfun(*,chan-1)
cymax = fix(alog10(max(emfun(*,chan-1)))+1.5)
plot_yrange = [10.^(cymax-5),10.^cymax]
 
plot_over = 0 

for ii=0,nchans-1 do begin

   k = chan(ii) - 1

   qqq=where(tzaa(*,k) gt 0)
   qqqm=max(qqq)
   !c = 0 

   if not plot_over then begin
     plot_oo,tzaa(0:qqqm,k),emfun(*,k),xrange=plot_xrange,yrange=plot_yrange, $
        xtitle=plot_xtitle,ytitle=plot_ytitle,title=plot_title,charsiz=1.2
   endif else begin
     oplot,tzaa(0:qqqm,k),emfun(*,k)
   endelse
   
   plot_over=1
 
   endfor

ans=''
if !d.name  ne 'PS' then read,'Next plot ',ans

;;plot_title='Yohkoh BCS Fluxes    (EM = 1.e50 cm-3) '+abun_desc+' abundances'
plot_title='Yohkoh BCS Fluxes    '+emstr+'  '+abun_desc+' abundances'
plot_ytitle='Photons/cm2/sec in line'

;set_xy,1.e+6,9.e7,100.,9.e+6

;   plot_yrange = [(min(flux(*,chan-1))>1),max(flux(*,chan-1))]
;pmm,flux(*,chan-1)
cymax = fix(alog10(max(flux(*,chan-1)))+1.5)
plot_yrange = [10.^(cymax-5),10.^cymax]


plot_over = 0

 for ii=0,nchans-1 do begin

   k = chan(ii) - 1

   qqq=where(tzaa(*,k) gt 0)
   qqqm=max(qqq)
   !c = 0 

   if not plot_over then begin
     plot_oo,tzaa(1:qqqm,k),flux(1:qqqm,k),xrange=plot_xrange,yrange=plot_yrange, $
        xtitle=plot_xtitle,ytitle=plot_ytitle,title=plot_title,charsiz=1.2
   endif else begin
     oplot,tzaa(1:qqqm,k),flux(1:qqqm,k)
   endelse
   
   plot_over=1
 
   endfor

if !d.name ne 'PS' then read,'Next plot ',ans


  plot_title='Yohkoh BCS Total Counts    '+emstr	;(EM = 1.e50 cm-3) 
  plot_ytitle='Counts/sec in line'

;   plot_yrange = [(min(counts(*,chan-1))>1),max(counts(*,chan-1))]
;pmm,counts(*,chan-1)
cymax = fix(alog10(max(counts(*,chan-1)))+1.5)
plot_yrange = [10.^(cymax-5),10.^cymax]

   plot_over = 0

 for ii=0,nchans-1 do begin

   k = chan(ii) - 1

   qqq=where(tzaa(*,k) gt 0)
   qqqm=max(qqq)
   !c = 0 


   if not plot_over then begin
     plot_oo,tzaa(1:qqqm,k),counts(1:qqqm,k),xrange=plot_xrange,yrange=plot_yrange, $
        xtitle=plot_xtitle,ytitle=plot_ytitle,title=plot_title,charsiz=1.2
   endif else begin
     oplot,tzaa(1:qqqm,k),counts(1:qqqm,k)
   endelse
   
   plot_over=1
 
   endfor

if !d.name ne 'PS' then read,'Next plot ',ans

plot_title='Yohkoh BCS Peak Counts    '+emstr	;(EM = 1.e50 cm-3)'
plot_ytitle='Counts/sec in peak bin'
;   plot_yrange = [(min(peaks(*,chan-1))>1),max(peaks(*,chan-1))]
;pmm,peaks(*,chan-1)
cymax = fix(alog10(max(peaks(*,chan-1)))+1.5)
plot_yrange = [10.^(cymax-5),10.^cymax]

plot_over = 0

 for ii=0,nchans-1 do begin

   k = chan(ii) - 1

   qqq=where(tzaa(*,k) gt 0)
   qqqm=max(qqq)
   !c = 0 


   if not plot_over then begin
     plot_oo,tzaa(1:qqqm,k),peaks(1:qqqm,k),xrange=plot_xrange,yrange=plot_yrange, $
        xtitle=plot_xtitle,ytitle=plot_ytitle,title=plot_title,charsiz=1.2
   endif else begin
     oplot,tzaa(1:qqqm,k),peaks(1:qqqm,k)
   endelse
   
   plot_over=1
 
   endfor

if !d.name ne 'PS' then read,'Next plot ',ans

!c = 0 
  
end