PRO xss_sim, line_todo, NO_BND = no_bnd, NO_PLOT = no_plot, ITER=iter
; Generates an approximation to the XSS spectrum
; for a given CMDB file/line
; Requires that the CMDB file be loaded already
; (by ~dd/idl/cmdb/cmdb_load).
; Calling Sequence:
;
;     xss_sim[, line]
;
; Arguments:
;
;     line Integer line of the cmdb to be used, lines start
;          with line 1.
;          If no line_todo is specified then the internal
;          values are used to model a source.  Editing this
;          file can change those "manual" inputs.
;
; Keywords:
;            NO_BND    This suppresses the creation of the BND products
;                      and considerably speeds up execution! (Most time is
;                      spent blurring the BND...).
;
; Output products:
;
;     twr_id.spec      Spectrum file in MARX format
;     twr_id.spec.ps   Plot of spectrum
;     twr_id.BND.ps    Plot of spectrum as seen by a BND in 1000 sec.s
;     twr_id.BND.hist  The BND simulation in two column format (E, counts)
;     lx_es, lx_tubespec  These are in labx_common and can be used
;                         by other routines.
;
; Uses many of dd's idl routines in ~dd/idl/labxray and
; ~dd/idl/useful
; CMDB input/interpretation is provided by routines in ~dd/idl/cmdb
;
; 10/25/96 dd unbundled from toga_full_sim
; 1/8/97 dd input is line_todo
; 12/16/97 dd Modified for XRCF-HRMA source distance and NO_BND.
; 12/29/97 dd Modified to have a mono_effic parameter
; 1/12/98 dd Modified for ITER keyword: used to select a source
;            iteration.  It's value is ignored by all but the DCM and
;            HIREFS.  For these, the ITER selects among the scanned
;            energies/wavelengths.  ITER starts with value 0.
; 7/17/98 dd Add a x30 for the O-K contamination in Fe anode (for Fe-L).
; 4/9/99 dd Add Pb (2 microns per 37 microns) to Mn filters
; - - - - - - - - -
@labx_common  ; common for a variety of laboratory-related parameters
@cmdb_common  ; Use dd's CMDB reader


; There are a very finite set of parmeters that are used to define
; this source spectrum model
; If no row_todo is supplied the values below are used - useful
; for manual testing...

hrma_entrance = 52601.2  ; cm

; - - - - - Main parameters to change - - - - - 
if n_elements(line_todo) LE 0 then begin

  ; ** Manual entry - if no ria_todo array given - modify items below
  ;    in this file...

  ;All sources start with an EIPS except Penning
  ; spect_file used to select EIPS etc ('use_Dans_labxray_routines')
  ; OR Penning ('use_Penning')
  spect_file = 'use_Dans_labxray_routines' 


  ; EIPS, HIREFS, and DCM start with a target...
  lx_tname = 'Al' ; Select from 
                  ; Be,B,C,O,Mg,MgO,Al,Si,SiO,Ti,Cr,Fe,Ni,Cu,Zr,Nb,Mo,Ag,Sn,W
  lx_sV = 12.0 ; kV
  lx_sI = 1.0 ; mA

  ; Monochromators are implemented as a narrow band filter after the
  ; underlying EIPS source
  mono_energy = -1.0   ; Monochromator used if greater than 0.
  mono_res = 300.0     ; resolving power of monochromator
  mono_highorders = 1  ; 0 = no high orders, 1 = HIREFS approx
  mono_effic = 0.1     ; Monochromtor efficiency

  ; For ALL sources the filters are then applied:

  fmat1 = 'Al'
  fthk1 = 0.0
  fmat2 = 'NONE'
  fthk2 = 0.0

  run_name = 'xss_sim_out'

end else begin
  ; cmdb data uses idl row number (starts at 0)
  row_todo = line_todo - 1

  ; - - - This whole ELSE is for CMDB usage: it fills the above
  ;        variables based on CMDB line
  ; ** CMDB row index is specified, fill values 
  ; below using CMDB data

  ; Complain if there is no CMDB loaded
  if n_elements(cmdb_lines) LE 0 then begin
    print, ' xss_sim: CMDB not loaded! '
    RETURN
  end

  ; define the column integers to access the CMDB data
  ; CMDB can be either a .cal or a .trw file
  if cmdb_using_cal then begin
@cmdb_cal_integers  ; include file of column integer def'ns.
  end else begin
@cmdb_trw_integers
  end

  spect_file = 'use_Dans_labxray_routines'
  ; Create spectrum file using Dan's X-GEF EIPS model...
  ; All the sources except Penning start with a EIPS source
  ; setup EIPS source
  CASE STRUPCASE(cmdb_fields(row_todo,c_source)) OF
    'EIPS': begin
      ; Set the EIPS target
      pieces = STR_SEP(cmdb_fields(row_todo,c_target),'-')
      if n_elements(pieces) LE 1 then begin
        print, ' * no valid target value'
        RETURN  ; can't simulate this
      end
      ct = pieces(0)  ; extract the target material
      ; Check if I have this target
      t_list = ['Be','B','C','O','Mg','MgO','Al', 'Si', 'SiO', $
			'Ti','Cr','Fe','Ni','Cu','Mo','Nb','Zr','Ag','Sn','W']
      ; if it is not in the list, return
      if MAX(ct EQ t_list) NE 1 then begin
        print, ' * target not in target list'
        RETURN
      end
      ; if it is 'O' or 'Si' use 'SiO'
      if ct EQ 'O' OR ct EQ 'Si' then ct = 'SiO'
      ; set EIPS current to value in CMDB
      temp_curr = 0.0
      READS, cmdb_fields(row_todo, c_src_current), temp_curr
      lx_sI = temp_curr
      ; no monochromator
      mono_energy = -1.0
      mono_res = -1.0
      mono_highorders = 0
      mono_effic = 0.1
    end

    ; For the monochromators the energy is selected (when it is a
    ; scan) by the value of ITER.
    ; Parameters determining monochromator energy selected:
    ;  * X_Ray_mono_initial     keV or A: num_value or N/A
    ;    X_Ray_mono_res_power
    ;  * X_Ray_mono_range       keV or A: num_value or N/A or 0
    ;  * X_Ray_mono_step_size   keV or A: num_value or N/A or 0
    ;    X_Ray_source_param    
    ; The function cmdb_mono_energies returns the list of energies.

    'HIREF-W': begin
      print, 'Source parameters: '+cmdb_fields(row_todo,c_X_ray_source_params)
      ct = 'W'
      lx_sI = 50.0 ; mA
      mono_energies = cmdb_mono_energies(row_todo)
      if n_elements(iter) EQ 0 then iter=0
      mono_energy = mono_energies((iter > 0) < (n_elements(mono_energies)-1))
      ; 12/29/97 dd Measured HSI data use to set the E/dE...
      mono_res = 60.0
      mono_highorders = 1
      ; default
      mono_effic = 0.001
      ; Set the efficiency based on the HIREFS crystal used:
      if STRPOS(cmdb_fields(row_todo,c_X_ray_source_params),',SX,') $
	GE 0 then begin
        mono_effic = 0.00015
        if mono_energy GT 1.3 then mono_effic = 0.0012
      end
      if STRPOS(cmdb_fields(row_todo,c_X_ray_source_params),',SA,') $
	GE 0 then mono_effic = 0.0002
      if STRPOS(cmdb_fields(row_todo,c_X_ray_source_params),',SB,') $
	GE 0 then begin
        mono_effic = 0.0007  ; set by E-LXH-EA-36.001
        if mono_energy LT 0.14 then begin
          ;mono_effic = 0.0006
          mono_effic = 0.00028
          mono_res = 80.
        end
      end
      if STRPOS(cmdb_fields(row_todo,c_X_ray_source_params),',SC,') $
	GE 0 then mono_effic = 0.0003
    end

    'HIREF-C': begin
      ct = 'C'
      lx_sI = 50.0 ; mA
      mono_energies = cmdb_mono_energies(row_todo)
      if n_elements(iter) EQ 0 then iter=0
      mono_energy = mono_energies((iter > 0) < (n_elements(mono_energies)-1))
      mono_res = 300.0
      mono_highorders = 1
      mono_effic = 0.1
    end

    'DCM': begin
      ct = 'W'
      lx_sI = 50.0 ; mA
      mono_energies = cmdb_mono_energies(row_todo)
      if n_elements(iter) EQ 0 then iter=0
      mono_energy = mono_energies((iter > 0) < (n_elements(mono_energies)-1))
      mono_res = 150.0
      mono_highorders = 1
      mono_effic = 0.01
    end

    'PENNING': begin
      spect_file = 'use_Penning'
    end

    ELSE: begin
      print, ' *** Unrecognized Source: '+cmdb_fields(row_todo,c_source)
      RETURN
    end

  ENDCASE

  if spect_file EQ 'use_Dans_labxray_routines' then begin
    ; EIPS getting source parameters
    lx_tname = ct ;
    ; get voltage
    voltage  = 0.0
    READS, cmdb_fields(row_todo, c_voltage), voltage
    lx_sV =  voltage; kV
  end

  fmat1 = cmdb_fields(row_todo,c_filter_mat_1)
  fthk1 = 0.0
  ; dont need thickness for no filter or OBF
  if NOT(fmat1 EQ 'NONE' OR fmat1 EQ 'N/A' OR fmat1 EQ '' OR $
					fmat1 EQ 'OBF') then begin  
    READS, cmdb_fields(row_todo, c_filter_thick_1), fthk1
  end
  fmat2 = cmdb_fields(row_todo,c_filter_mat_2)
  fthk2 = 0.0
  ; dont need thickness for no filter or OBF
  if NOT(fmat2 EQ 'NONE' OR fmat2 EQ 'N/A' OR fmat2 EQ '' OR $
					fmat2 EQ 'OBF') then begin  
    READS, cmdb_fields(row_todo, c_filter_thick_2), fthk2
  end

  run_name = cmdb_fields(row_todo, c_id)

end

; report parameters:'
  print, ''
  if n_elements(line_todo) GE 0 then begin
    print, '         xss_sim Parameters (cmdb line ='+$
	STRING(line_todo,FORMAT='(I5)')+')'
  end else begin
    print, '         xss_sim Parameters'
  end
  print, ''
  print, ' spect_file    : '+spect_file
if   spect_file EQ 'use_Dans_labxray_routines' then begin
  print, ' . . . . . EIPS and Monos . . . . . . . .'
  print, '  target       : '+lx_tname
  print, '  voltage      : ',lx_sV
  print, '  current      : ',lx_sI
  print, '  mono_energy  : ',mono_energy,' keV'
  print, '  mono_res     : ',mono_res
  print, '  mono_highord : ',mono_highorders,'  (0=DCM, 1=HIREFS)'
  print, '  mono_effic   : ',mono_effic
end
  print, ' . . . . . . . . . . . . . . . . . . . . .'
  print, '  filter 1     : '+fmat1, fthk1
  print, '  filter 2     : '+fmat2, fthk2
  print, ' . . . . . . . . . . . . . . . . . . . . .'
  print, '  run_name     : '+run_name

;.................................................................
; NOW do the real work...

CASE spect_file OF
'use_Dans_labxray_routines': begin
  ; Create spectrum
  ;lx_tname already setup
  ;lx_sV already setup
  ;lx_sI = 1.0 ; mA
  lx_sangle = 45.0

  ; Energy grid is linear and fine!
  lx_emin = 0.1
  lx_emax = (lx_sV < 9.9)
  lx_estep = 0.0005
  ; reduce the emin for C and Be
  if lx_tname EQ 'Be' OR lx_tname EQ 'C' then begin
    lx_emin = 0.050
    lx_emax = (lx_emax < 10.0)
    lx_estep = (10.0/2.5E4)
  end
  ; do the Mg-K at high resolution and cut out highest voltages
  if lx_tname EQ 'Mg' then begin
    lx_emin = 0.1
    lx_emax = (lx_sV < 5.0)
    lx_estep = 0.00025
  end

  ; fill lx_es and lx_tubespec
  ; lx_tubespec has units of photons/(sec bin steradian)

  ; add C and O contamination lines
  save_target = lx_tname
  save_curr = lx_sI
  ;
  lx_tname = 'C'
  contam_frac_I = 0.1
  lx_sI = contam_frac_I * save_curr
  labx_tubespec, /XRCF
  contam_spec = lx_tubespec
  ;
  ; don't add O to SiO !
  if save_target NE 'SiO' then begin
    if save_target EQ 'Fe' then begin
      ; The XRCF Fe anode has a ton of O-K...
      contam_frac_I = 0.3 * 28.0  ; for Fe-L !
    end else begin
      contam_frac_I = 0.3
    end
    lx_tname = 'O'
    lx_sI = contam_frac_I * save_curr
    labx_tubespec, /XRCF
    contam_spec = contam_spec + lx_tubespec
  end

  ; put values back
  lx_tname = save_target
  lx_sI = save_curr

  ; real spectrum
  labx_tubespec, /XRCF

  ; add in the contamination
  lx_tubespec = lx_tubespec + contam_spec

  ; Include monochromator
  ; Mono is a filter that has an overall efficiency factor followed by
  ; different order realative efficiencies: 100% at first order, 3% at second,
  ; order 0.5% third order, 0.7% 4th order, 0.3% 5th order
  ; and 1.E-5 at other energies.
  if mono_energy GT 0.0 then begin
    ; Overall effic (reflectivity of DCM, HIREFS effic)
    lx_tubespec = mono_effic * lx_tubespec
    ; Find the bandpass bins...
    mono_bins = where( approx_equal(lx_es, mono_energy, 0.5/mono_res), $
	nmonobins)
    if nmonobins LE 0 then RETURN
    print, ' nmonobins = ', nmonobins
    ; Adjusted for E-HXH-EA-35.002 (was 1.E-6) - need beeter measurement!
    cont_level = 1.E-4  
    lx_tubespec = cont_level * lx_tubespec  ; continuum level
    lx_tubespec(mono_bins) = (1./cont_level) * lx_tubespec(mono_bins)
    if mono_highorders EQ 1 then begin
      ; add higher orders
      order_fracs = [0.03, 0.005, 0.007, 0.003]  ; HIREFS ?
      mono_order = [2.,3.,4.,5.]
      for io=0,n_elements(order_fracs)-1 do begin
        mono_bins = where( approx_equal(lx_es, mono_order(io)*mono_energy, $
  		0.5/mono_res), nmonobins)
        if nmonobins GT 0 then begin
          lx_tubespec(mono_bins) = (1./cont_level) * $
		order_fracs(io)*lx_tubespec(mono_bins)
        end
      end
    end
  end
  ; convert to photons/cm^2/s/keV
  ; Real flux at XRCF at TMA/HRMA input
  lx_tubespec = ( lx_tubespec/(lx_es(1)-lx_es(0)) ) * (1./hrma_entrance)^2

end ; CASE of 'use Dans...'

'use_Penning': begin
  mono_energy = -1.
  ; use a Kramer continuum 'cause there is some!
  ; make sure variables are setup
  lx_tname = 'Al'
  labx_target
  ; Now change some values 
  lx_tZ = 0.1 ; Set to get continuum about right
  lx_tname = 'Penning'
  lx_sV = 0.9  ; kV
  lx_sI = 800.0 ; mA
  ; setup lx_ variables ...
  lx_emin = 0.050
  lx_emax = lx_sV
  lx_estep = 0.0001
  ; generate a grid
  lx_es = lx_emin+lx_estep*INDGEN(1+(lx_emax-lx_emin)/lx_estep)    
  ; and call the continuum model
  lx_tubespec = labx_tubeinuum(lx_es)
  if lx_sV GT 0.53 then begin
    ; add some O
    aboveO = where(lx_es GE 0.524)
    lx_tubespec(aboveO(0)) = lx_tubespec(aboveO(0))*1000.0
  end
  ; Reduce the continuum by factor of 4
  lx_tubespec = 0.25 * lx_tubespec

  ; convert to photons/cm^2/s/keV
  ; Real flux at XRCF at TMA input
  lx_tubespec = ( lx_tubespec/(lx_es(1)-lx_es(0)) ) * (1./hrma_entrance)^2

  ; Add in the Penning lines
  ; put in three place holders
  ; 160A , 10,000 ph/cm2sec
  ; 130A , 50 ph/cm2sec
  ; 56A , 4 ph/cm2sec
  ; p_lams = [160.,130.,56.]
  ; p_flux = [10000., 50., 4.]
  ; wavelengths from JeffK from Finley SPIE vol.18 no 5 (1979)p.649
  p_lams = [170., 160., 137.76, 137.6, 135.8, 133.0, 130.65, 130., $
	127.69, 127.16, 125.618, 124.75, 124.2, 117., $
	111.5, 108., 104., 99.5, 94., 90.8, 88.6, 85.6, $
	59.0, 59.27, 56.5, 56., 55.425]    
  p_flux = [2000.,8000., 5.0, 5.0, 10.0, 10.0, 18.0, 33., $
	 2., 2., 2., 6., 12., 12., $
	15., 15., 7., 1.5, 1.5, 1.5, 0.8, 0.8, $
	0.4, 0.4, 1.1, 3.5, 1.0]
  for il=0,n_elements(p_lams)-1 do begin
    e_lam = lx_akeV/p_lams(il)
    above_line = where(lx_es GE e_lam)
    lx_tubespec(above_line(0)) = lx_tubespec(above_line(0)) + $
		p_flux(il)/(lx_es(1)-lx_es(0))
  end
  ; ALWAYS have a 0.5 um parylene window
  lx_tubespec = lx_tubespec * (parylene(lx_es) > 1.E-10)^0.5

  ; Adjust the intensity to agree better with D-LXH-dF-7.001
  lx_tubespec = 0.05 * lx_tubespec

end  ; CASE of 'use_Penning'

ELSE: begin
  ; something else?!?
  print, ' spect file is unrecognized: '+spect_file
  RETURN
end

ENDCASE

; Include the filters for all sources
fmats = [fmat1,fmat2]
fthks = [fthk1,fthk2]
for ifilter = 0,1 do begin  
  CASE fmats(ifilter) OF
    'Al': begin
        ; Catch zero thickness values - put filter in anyway
        if fthks(ifilter) EQ 0 then begin
          if ifilter EQ 0 then fthks(ifilter)=18.38 else fthks(ifilter)=36.76
        end
        lx_tubespec = lx_tubespec * aluminum(lx_es)^fthks(ifilter)
      end
    'Fe': lx_tubespec = lx_tubespec * iron(lx_es)^fthks(ifilter)
    'Be': begin
        ; Catch zero thickness values - put filter in anyway
        if fthks(ifilter) EQ 0 then begin
          if ifilter EQ 0 then fthks(ifilter)=2.*0.903 else $
			fthks(ifilter)=4.*0.903
        end
        lx_tubespec = lx_tubespec * $
		(beryllium(lx_es,/MICRON) > 1.E-10)^fthks(ifilter)
      end
    'B': begin
        ; Catch zero thickness values - put filter in anyway
        if fthks(ifilter) EQ 0 then begin
          if ifilter EQ 0 then fthks(ifilter)=2.*1.25 else $
			fthks(ifilter)=4.*1.25
        end
        lx_tubespec = lx_tubespec * $
		(boron(lx_es,/MICRON) > 1.E-10)^fthks(ifilter)
      end
    'C': begin
        ; Catch zero thickness values - put filter in anyway
        if fthks(ifilter) EQ 0 then begin
          if ifilter EQ 0 then fthks(ifilter)=2.*5.16 else fthks(ifilter)=4.*5.16
        end
        lx_tubespec = lx_tubespec * (parylene(lx_es) > 1.E-10)^fthks(ifilter)
      end
    'C8H8': begin
        ; Catch zero thickness values - put filter in anyway
        if fthks(ifilter) EQ 0 then begin
          if ifilter EQ 0 then fthks(ifilter)=2.*5.16 else fthks(ifilter)=4.*5.16
        end
        lx_tubespec = lx_tubespec * (parylene(lx_es) > 1.E-10)^fthks(ifilter)
      end
    'Mn': begin
            ; The XRCF Mn filter includes Pb!  2 um Pb for each 37 um of Mn
            lx_tubespec = lx_tubespec * manganese(lx_es)^fthks(ifilter)
            lx_tubespec = lx_tubespec * lead(lx_es)^(fthks(ifilter)*2.0/37.0)
          end
    'Mg': lx_tubespec = lx_tubespec * magnesium(lx_es)^fthks(ifilter)
    'Mo': lx_tubespec = lx_tubespec * moly(lx_es)^fthks(ifilter)
    'Ag': lx_tubespec = lx_tubespec * silver(lx_es)^fthks(ifilter)
    'Sn': lx_tubespec = lx_tubespec * tin(lx_es)^fthks(ifilter)
    'Zr': begin
        if fthks(ifilter) EQ 0 then begin
          if ifilter EQ 0 then fthks(ifilter)=2.*1.98 else fthks(ifilter)=4.*1.98
        end
        lx_tubespec = lx_tubespec * zirconium(lx_es)^fthks(ifilter)
      end
    'Ti': lx_tubespec = lx_tubespec * titanium(lx_es)^fthks(ifilter)
    'TiK': lx_tubespec = lx_tubespec * titanium(lx_es)^fthks(ifilter)
    'TiL': lx_tubespec = lx_tubespec * titanium(lx_es)^fthks(ifilter)
    'Cr': lx_tubespec = lx_tubespec * chromium(lx_es)^fthks(ifilter)
    'CrK': lx_tubespec = lx_tubespec * chromium(lx_es)^fthks(ifilter)
    'CrL': lx_tubespec = lx_tubespec * chromium(lx_es)^fthks(ifilter)
    'V': lx_tubespec = lx_tubespec * vanadium(lx_es)^fthks(ifilter)
    ; Shift the CuL energy slightly to get better agreement between the
    ; Cu nofilter and Cu w/filter data
    'Cu': lx_tubespec = lx_tubespec * copper(0.998*lx_es)^fthks(ifilter)
    'CuL': lx_tubespec = lx_tubespec * copper(0.998*lx_es)^fthks(ifilter)
    'CuK': lx_tubespec = lx_tubespec * copper(0.998*lx_es)^fthks(ifilter)
    'Ni': lx_tubespec = lx_tubespec * nickel(lx_es)^fthks(ifilter)
    'NiK': lx_tubespec = lx_tubespec * nickel(lx_es)^fthks(ifilter)
    'NiL': lx_tubespec = lx_tubespec * nickel(lx_es)^fthks(ifilter)
    'OBF': lx_tubespec = lx_tubespec * acis2c_obf(lx_es)
    'NONE':
    'N/A':
    ELSE: print, '*** Unrecognized filter : '+fmats(ifilter)+'.'
  ENDCASE
end

; Report the resulting total flux, etc.
total_flux = TOTAL(lx_tubespec)*(lx_es(1)-lx_es(0))
if mono_energy GT 0. then begin
  line_energy = mono_energy
end else begin
  line_energy = 0.0
  READS, cmdb_fields(row_todo, c_energy), line_energy
end
line_region = where(lx_es GT 0.95*line_energy AND $
		lx_es LT 1.05*line_energy)
line_flux = TOTAL(lx_tubespec(line_region))*(lx_es(1)-lx_es(0))
sim_flux_str = 'Total photons/cm^2 s = '+STRING(total_flux,FORMAT='(G8.3)')
line_flux_str = STRING(line_energy,FORMAT='(F6.3)')+ ' keV' + $
	' phot/cm^2 s = '+STRING(line_flux,FORMAT='(G8.3)')
print, '..........................................'
print, '  '+sim_flux_str
print, '  '+line_flux_str
cmdb_flux_str = ''
if n_elements(row_todo) GE 1 then begin
  if cmdb_using_cal then begin
    cmdb_flux_str = 'CMDB  = '+cmdb_fields(row_todo, c_flux)+ ' [' + $
	cmdb_fields(row_todo, c_src_current)+' mA]'
  end else begin
    cmdb_flux_str = 'CMDB  = '+ $
	cmdb_fields(row_todo, c_X_Ray_Source_flux)+ ' [' + $
	cmdb_fields(row_todo, c_X_ray_source_current)+' mA]'
  end  
  print,   '           '+ cmdb_flux_str
end
print, ''

; write out the data in photons/cm^2/s/keV
spect_file = run_name + '.spec'
OPENW, spect_unit, spect_file, /GET
for ie=0,n_elements(lx_es)-1 do printf, spect_unit, $
	lx_es(ie), lx_tubespec(ie)
CLOSE, spect_unit
FREE_LUN, spect_unit

!p.multi = 0

; output a plot of the spectrum
; in photons/cm^2/s/keV
sourcexrange = [(MIN(lx_es) > 0.02),(MAX(lx_es) < 10.0)]
max_flux = MAX(lx_tubespec)
if NOT(KEYWORD_SET(NO_PLOT)) then begin
pre_print_sqr
plot_oo, lx_es, lx_tubespec, $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XTITLE = 'Energy (keV)', $
	XRANGE = sourcexrange, $
	YTITLE='Flux (photons/cm^2/s/keV, bin size = '+ $
		STRING(1000.0*(lx_es(1)-lx_es(0)),FORMAT='(F5.2)')+' eV)', $
	YRANGE = [1.E-8*max_flux,max_flux]
; Add the flux strings
XYOUTS, sourcexrange(0)*1.5, 0.3*max_flux, sim_flux_str
XYOUTS, sourcexrange(0)*1.5, 0.1*max_flux, line_flux_str
XYOUTS, sourcexrange(0)*1.5, 0.03*max_flux, cmdb_flux_str

device, /close
set_plot, 'X'
SPAWN, 'cp idl.ps '+run_name+'.spec.ps'
end ; NO_PLOT
;- - - - - - - BND output - - - -
if NOT(KEYWORD_SET(NO_BND)) then begin

; output a plot of the spectrum as seen by BND
; lx_tubespec is in photons/cm^2/s/keV
bnd_time = 1000.0
bnd_spec = lx_tubespec * 35.75 * bnd_time * (lx_es(1)-lx_es(0))  ; photons per bin
; 
; For BND detector we need to add efficiency and energy blur
; include the FPC efficiency, hxda_fpc(E).
bnd_temp = hxda_fpc(lx_es) * bnd_spec
bnd_es = lx_es

; Add the blur
;  ... first ...
; Increase the bin size by factor of 20
bnd_estep = 20.0 * lx_estep
ipoint = LONG(0)
iout = 0
while ipoint+19 LE n_elements(bnd_temp)-1 do begin
  bnd_es(iout) = TOTAL(lx_es(ipoint:ipoint+19))/20.0
  bnd_spec(iout) = TOTAL(bnd_temp(ipoint:ipoint+19))
  ipoint = ipoint + 20
  iout = iout + 1
end
bnd_es = bnd_es(0:iout-1)
bnd_spec = bnd_spec(0:iout-1)

; add Ar escape too!
save_tube = lx_tubespec
save_es = lx_es
lx_tubespec = bnd_spec
lx_es = bnd_es
lx_dname = 'SMD'
labx_escape
bnd_spec = lx_tubespec
lx_tubespec = save_tube
lx_es = save_es

; detector resolution vs energy, FWHM in keV
; use resolution of the X-GEF SMD...
res = labx_smd_res(bnd_es)
res = res/(2.35*bnd_estep) ; one-sigma in bins

npts = n_elements(bnd_spec)
nptsm = npts-1
out_spec = 0. * bnd_spec
dout = 0. * bnd_spec

; go through the input spectrum
for in = 0, nptsm do begin
  res_e = res(in)  ; one_sigma in bins at this energy
  amp = bnd_spec(in)
  dout = 0. * dout
  for io = 0, nptsm do begin
    sigmas = (in-io)/res_e
    if(sigmas GT -6.0 AND sigmas LT 6.0) then $
	      dout(io) = EXP(-(sigmas^2)/2.)
  end
  out_spec = out_spec + (amp/res_e)*dout
end

; Add Poisson counting stats
bnd_spec = g_poiss(out_spec/SQRT(2.*!PI))

; write out the BND histogram data
hist_file = run_name + '.BND.hist'
OPENW, hist_unit, hist_file, /GET
for ie=0,n_elements(bnd_es)-1 do printf, hist_unit, $
	bnd_es(ie), bnd_spec(ie)
CLOSE, hist_unit
FREE_LUN, hist_unit

max_flux = MAX(bnd_spec)
if NOT(KEYWORD_SET(NO_PLOT)) then begin
pre_print_sqr
plot_oo, bnd_es, bnd_spec, PSYM=10, $
	TITLE = run_name + $
	' ( BND Counts,'+STRING(bnd_time,FORMAT='(F7.0)')+' sec, '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XTITLE = 'Energy (keV)', XRANGE = sourcexrange, $
	YTITLE='Flux (BND events/bin, bin size ='+$
		STRING(bnd_estep*1000.0,FORMAT='(F6.1)')+' eV)', $
	YRANGE = [0.1,max_flux]
device, /close
set_plot, 'X'
SPAWN, 'cp idl.ps '+run_name+'.BND.ps'
SPAWN, 'rm idl.ps'
end ; NO_PLOT

end  ; of NO_BND
; - - - -

!p.multi = 0

RETURN
END