PRO write_xcat_Serkowski,file,lambda,K=K,lambda_max=lambda_max,psAv_max=psAv_max

;##################################################################
; NAME:
;	write_xcat_Serkowski
; PURPOSE:
;    Computes observed sed and color correction for a given spectrum in a given filter
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    write_xcat_Serkowski,lambda,K=K,lambda_max=lambda_max,psAv_max=psAv_max
; INPUTS:
;    file: output IPAC file
;    lambda: array of wavenumbers 
; OPTIONAL INPUT PARAMETERS:
;    K: Serkowski K factor
;    lambda_max: wavelength of maximal polarization
;    psAv_max: maximal fraction of polarization p/Av
; OUTPUTS:
;    sed: SED as observed in filters provided in filter_names
; OPTIONAL OUTPUT PARAMETERS:
;    cc        = color correction coefficients (spec*cc=sed)
; ACCEPTED KEY-WORDS:
;    fluxconv  = if set, these are taken for flux conventions
;              possible values: 'nuInu=cste', 'IRAC', 'MIPS', 'CMB'
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;
; MODIFICATION HISTORY:
;    Written by Vincent Guillet (Jan. 2012)
;
;##################################################################

; Optional parameters

; Extrapolate to NIR with a power-law (Martin etal .)
; Draine & Fraisse (2009) : Power law coef -1.7 from 1.39 to 5 microns
;IF n_elements(K) EQ 0 THEN K = 0.92
;lmin=1.39
;beta=-1.7
; Martin & Whittet (1990) : Power law coef -1.8 from 0.9 to 5 microns
IF n_elements(K) EQ 0 THEN K = 1.15
lmin=0.9
beta=-1.8

IF n_elements(lambda_max) EQ 0 THEN lambda_max = 0.55

; Max % polarization by extinction observed in the diffuse medium is 
;       p = 0.09 * E(B-V)
; E(B-V) = NH / (5.8d21 H cm-2 mag-1) (Bohlin, Savage & Drake 1978)
; A(V) = NH / (1.87d21 H cm-2 mag-1) 
; Av = E(B-V) * Rv
; Av = tau(V) * 1.086
; Therefore p/Av = 0.09 / Rv = 0.029 with Rv=3.1
;           p/tau = 0.09/Rv*1.086 = 0.31 (and not 0.0267 like in Martin 2007)
;----------------------------------------------------------------------------
; dem:	p 	= 0.09 E(B-V)
; 	tau 	= E(B-V) * Rv / 1.086
;  donc	p/tau 	= 0.09 / Rv * 1.086 = 0.031, not 0.0267 = 0.09 / Rv / 1.086
;----------------------------------------------------------------------------

IF n_elements(psAv_max) EQ 0 THEN psAv_max = 0.029
; The corresponding maximal dust polarization cross-section is:
; sigma_pol_max = (p/Av)_max / 1.87d21 H/cm2 = psAv_max/1.87 (x10^-21 cm2/H)
sigma_pol_max_per_H = psAV_max / 1.87 

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

nlambda = n_elements(lambda)

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

; Serkowski's law
extpol = sigma_pol_max_per_H * EXP( -K * ALOG(lambda/lambda_max)^2 )

im = WHERE(lambda gt lmin, count)
if count GT 0 then extpol(im) = extpol(im(0)) * (lambda(im)/lambda(im(0)))^beta
uncextpol = extpol * 0.15

st.instru = 'POLAR_EXT'
st.filter = 'SPECTRUM'
st.wave = lambda
st.spec = extpol
st.error = uncextpol
st.unit = 'cm2/H'

; Write out IPAC file
write_xcat,st,file

; Plot result
window,0
ploterror,lambda,extpol,0*lambda,uncextpol,psym=4,/xlog,/ylog,xrange=[min(lambda),max(lambda)]
oplot,lambda,extpol

stop

END