PRO xss_sim, row_todo
; 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[, row]
;
; Arguments:
;
;     row  Integer row of the cmdb to be used, rows start
;          with row 0.
;          If no row_todo is specified then the internal
;          values are used to model a source.  Editing this
;          file can change those "manual" inputs.
;
; Keywords:
;            NONE (yet!)
;
; 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)
;
; 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
; 

@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...

; - - - - - Main parameters to change - - - - - 
if n_elements(row_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,C,O,Mg,MgO,Al,Si,SiO,Ti,Cr,Fe,Cu,Zr,Mo,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

  ; 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
  ; - - - 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','C','O','Mg','MgO','Al', 'Si', 'SiO', $
			'Ti','Cr','Fe','Cu','Mo','Zr','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 1.0 mA
      lx_sI = 1.0
      ; no monochromator
      mono_energy = -1.0
      mono_res = -1.0
      mono_highorders = 0
    end

    'HIREF-W': begin
      ct = 'W'
      lx_sI = 50.0 ; mA
      mono_energy = 0.0
      READS, cmdb_fields(row_todo,c_energy), mono_energy
      mono_res = 300.0
      mono_highorders = 1
    end

    'HIREF-C': begin
      ct = 'C'
      lx_sI = 50.0 ; mA
      mono_energy = 0.0
      READS, cmdb_fields(row_todo,c_energy), mono_energy
      mono_res = 300.0
      mono_highorders = 1
    end

    'DCM': begin
      ct = 'W'
      lx_sI = 50.0 ; mA
      mono_energy = 0.0
      READS, cmdb_fields(row_todo,c_energy), mono_energy
      mono_res = 300.0
      mono_highorders = 1
    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, ''
  print, '         xss_sim Parameters'
  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'
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 < 5.0)
    lx_estep = (5.0/2.5E4)
  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
  contam_spec = lx_tubespec
  ;
  lx_tname = 'O'
  contam_frac_I = 0.3
  lx_sI = contam_frac_I * save_curr
  labx_tubespec
  contam_spec = contam_spec + lx_tubespec
  ; put values back
  lx_tname = save_target
  lx_sI = save_curr

  ; real spectrum
  labx_tubespec

  ; 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 = 0.10 * 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
    cont_level = 1.E-6
    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 input
  lx_tubespec = ( lx_tubespec/(lx_es(1)-lx_es(0)) ) * (1./53175.4)^2

end ; CASE of 'use Dans...'

'use_Penning': begin
  ; 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

  ; 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./53175.4)^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., 130., 124.2, 117., $
	111.5, 108., 104., 99.5, 94., 90.8, 88.6, 85.6, $
	59., 56.]    
  p_flux = [2000.,8000., 50., 15., 12., $
	15., 15., 7., 1.5, 1.5, 1.5, 0.8, 0.8, $
	1.3, 4.]
  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

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
    '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
    'Mn': lx_tubespec = lx_tubespec * manganese(lx_es)^fthks(ifilter)
    'Mg': lx_tubespec = lx_tubespec * magnesium(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)
    'Cu': lx_tubespec = lx_tubespec * copper(lx_es)^fthks(ifilter)
    'CuL': lx_tubespec = lx_tubespec * copper(lx_es)^fthks(ifilter)
    'CuK': lx_tubespec = lx_tubespec * copper(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))
sim_flux_str = 'Total photons/cm^2 s = '+STRING(total_flux,FORMAT='(G8.3)')
print, '..........................................'
print, '  '+sim_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)
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, cmdb_flux_str

device, /close
set_plot, 'X'
SPAWN, 'cp idl.ps '+run_name+'.spec.ps'

; 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 * 48.0 * 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)
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'

!p.multi = 0

RETURN
END