PRO write_xcat_Mathis1990

; #################################################################################
;
; Write observed extinction curve from Mathis (19990) in IPAC format xcat
; 
; AUTHORS :  
;       Mathieu Compiegne & Vincent Guillet
;
; #################################################################################

; To read xcat file :
; -------------------
;st=read_xcat(file)
;help,st,/str

; #################################################################################
; Visible and NIR extinction curve according to Mathis (1999) 
; #################################################################################

Mathis90_Al_sur_Aj  = [0.0015, 0.0041, 0.0071, 0.013, 0.048, 0.075, 0.083, 0.053, 0.098, 0.192, 0.208, 0.157, 0.070, 0.095, 0.182, 0.382, 0.624, 1.00, 1.70, 2.66, 3.55, 4.70, 5.53, 5.87, 6.90, 7.63, 9.03, 11.29, 10.08, 8.93, 9.44, 11.09, 12.71, 17.2, 19.1, 9.15, 7.31, 3.39, 1.35]
lamb_Mathis = [250, 100, 60, 35, 25, 20, 18, 15, 12, 10, 9.7, 9, 7, 5, 3.4, 2.2, 1.65, 1.25, 0.9, 0.7, 0.55, 0.44, 0.365, 0.33, 0.28, 0.26, 0.24, 0.218, 0.2, 0.18, 0.15, 0.13, 0.12, 0.091, 0.073, 0.041, 0.023, 0.004, 0.002]
sigma_Mathis = Mathis90_Al_sur_Aj *  1.51d-22 / 1.086
sigma_Mathis = sigma_Mathis * 1d21

lambmin  = 1.5
good = WHERE(lamb_Mathis gt lambmin and lamb_Mathis lt 23.)
lamb_Mathis     = lamb_Mathis[good]
ext_mathis      = sigma_mathis[good]
extunc_mathis   =  ext_mathis * 0.15
lambunc_mathis  = lamb_mathis * 0.

next    = 40
xextmin = 1./lambmin
xextmax = 9.
xext=FINDGEN(next)*(xextmax-xextmin)/(next-1)+xextmin
lamb_ext = 1./xext
lambunc   = lamb_ext * 0.

; Complete with Fitzpatrick
R_V = 3.1
NHsurEBV = 5.8d21
fitz99_ext, R_V, NHsurEBV, xext, ext, extunc

extunc = extunc * 1.d21
ext    = ext    * 1.d21
lamb_ext = [lamb_ext, lamb_Mathis]
ext      = [ext, ext_mathis]
lambunc  = [lambunc, lambunc_mathis]
extunc   = [extunc, extunc_mathis]
ext      = ext[SORT(lamb_ext)]
lambunc  = lambunc[SORT(lamb_ext)]
extunc   = extunc[SORT(lamb_ext)]
lamb_ext = lamb_ext[SORT(lamb_ext)]
next     = N_ELEMENTS(ext)
extunc[WHERE(lamb_ext gt 0.5)] = ext[WHERE(lamb_ext gt 0.5)] * 0.15


one_ps = {instru:'',filter:'',wave:0.,spec:0.,error:0.,unit:''}
st = replicate(one_ps,next)

st.instru = 'EXTINCTION'
st.filter = 'SPECTRUM'
st.wave = lamb_ext
st.spec = ext
st.error = extunc
st.unit = 'cm2/H'

; Write out IPAC file
write_xcat,st,'Mathis1990_DISM.xcat'

; print and plot result
window,1
ploterror,1./lamb_ext,ext,0.*lamb_ext,extunc,psym=4
oplot,1./lamb_ext,ext

END