readcol, 'acis_s_xray_trans_1198.qdp', nrg_new, obf_new, /silent
nstep_new = n_elements(nrg_new)
set_plot, 'ps'
device, /land, filen='obf_decompose.ps'
plot_oo, nrg_new, -nrg_new^(2.7) * alog(obf_new)

; account for Al-L opacity at E > .28 keV
c_edge = 0.2838
al_l_edge = 0.075
n_edge = 0.40
o_edge = 0.529
al_k_edge = 1.55

; First, model the Al region without EXAFS and extrapolate to E > C-K
al_form = fltarr( nstep_new )
al_model = fltarr( nstep_new )

; define Al regions
al_region = where( nrg_new le c_edge )
above_al_l_edge = where( nrg_new gt al_l_edge )
above_al_region = where( nrg_new gt c_edge )

al_pow = 2.7
plot, nrg_new, -nrg_new^al_pow * alog(obf_new), xra=[.05, 0.5]
al_form[above_al_l_edge] = 0.042 * $
	(1.-exp(-(nrg_new[above_al_l_edge]-al_l_edge)/.13) )
oplot, nrg_new, al_form, lines=1
al_model[al_region] = -alog( obf_new[al_region] )	; use measurements below C-K
al_model[above_al_region] = nrg_new[above_al_region]^(-al_pow) * $
	al_form[above_al_region]	; use model above C-K
al_opacity_adj = obf_new * exp( al_model )

; Second, model the C-K region without EXAFS and extrapolate to E > N-K
c_region = where( nrg_new gt c_edge and nrg_new le n_edge )
above_c_edge = where( nrg_new gt c_edge )
above_c_region = where( nrg_new gt n_edge )

plot, nrg_new, -nrg_new^al_pow * alog(al_opacity_adj), xra=[.25, 0.8]

c_form = fltarr( nstep_new )
c_form[above_c_edge] = 0.042 * (1.-exp(-(nrg_new[above_c_edge]-c_edge)/.035) )
oplot, nrg_new, c_form, lines=1
c_model = fltarr( nstep_new )
c_model[c_region] = -alog(al_opacity_adj[c_region])
c_model[above_c_region] = nrg_new[above_c_region]^(-al_pow) * c_form[above_c_region]
c_opacity_adj = al_opacity_adj * exp( c_model )

; Next, model the N-K region without EXAFS and extrapolate to E > O-K
n_pow = 2.5
n_region = where( nrg_new gt n_edge and nrg_new le o_edge )
above_n_edge = where( nrg_new gt n_edge )
above_n_region = where( nrg_new gt o_edge )
plot, nrg_new, -nrg_new^n_pow * alog(c_opacity_adj), xra=[.4, 0.8], yra=[0,.05]

n_form = fltarr( nstep_new )
n_form[above_n_edge] = 0.0042 * (1.-exp(-(nrg_new[above_n_edge]-n_edge)/.02) )
oplot, nrg_new, n_form, lines=1
n_model = fltarr( nstep_new )
n_model[n_region] = -alog(c_opacity_adj[n_region])
n_model[above_n_region] = nrg_new[above_n_region]^(-n_pow) * n_form[above_n_region]
n_opacity_adj = c_opacity_adj * exp( n_model )

; Next, model the N-K region without EXAFS and extrapolate to E > O-K
o_pow = 2.7
o_region = where( nrg_new gt o_edge and nrg_new le al_k_edge )
above_o_edge = where( nrg_new gt o_edge )
above_o_region = where( nrg_new gt al_k_edge )
plot, nrg_new, -nrg_new^o_pow * alog(n_opacity_adj), xra=[.5, 1.6], yra=[0, .1]

o_form = fltarr( nstep_new )
o_form[above_o_edge] = 0.032 * (1.-exp(-(nrg_new[above_o_edge]-o_edge)/.02) )
oplot, nrg_new, o_form, lines=1
device, /close
o_model = fltarr( nstep_new )
o_model[o_region] = -alog(n_opacity_adj[o_region])
o_model[above_o_region] = nrg_new[above_o_region]^(-o_pow) * $
	o_form[above_o_region]
o_opacity_adj = n_opacity_adj * exp( o_model )

; Finally, all the rest is Al-K absorption
above_al_k_edge = where( nrg_new gt al_k_edge )
al_k_model = fltarr( nstep_new )
al_k_model[above_al_k_edge] = -alog(o_opacity_adj[above_al_k_edge])

; show all transmissions as contribute to the OBF
device, /land, filen='filter_components.ps'
plot_oo, nrg_new, -alog( obf_new ), xra=[.1,2.], /xsty, $
	xtit='Energy (keV)', ytit='-ln (Transmission)'
oplot, nrg_new, al_model, lines=1
oplot, nrg_new, c_model, lines=2
oplot, nrg_new, o_model, lines=3
oplot, nrg_new, al_k_model, lines=4
oplot, nrg_new, n_model, lines=5
legend, ['Full OBF', 'Al-L', 'C-K', 'N-K', 'O-K', 'Al-L'], $
	lines=[0, 1,2,5,3,4], /right
device, /close

forprint, nrg_new, al_model, c_model, n_model, o_model, al_k_model, $
	text='filter_components.txt', /silent