write_xcat_Serkowski.pro
3.23 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
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