dustem_band_cc.pro
2.98 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
FUNCTION dustem_BAND_CC, w0, xs, ys, xt, yt, y0=y0, cc=cc, rr=rr, nint=nint
; returns flux in given instrument band (erg/cm2/sr)
; W0 (I): centroid of filter (any unit)
; XS (I): input x-grid for SED (same unit as W0)
; YS (I): input SED in erg/cm2/sr
; XT (I/O): x-grid of filter (same unit as W0)
; YT (I/O): filter transmission (arb. unit)
; RR (O): undersampling ratio: dx(SED)/dx(filt), usually 1.
; If RR<1 filter oversampled in flux integration
; Y0 (O): SED at band center
; CC (O): color correction
; NINT (O): nr of points in band integration
vlow = 1d-5
; remove bad and multiple values in filter
xtt=xt & ytt=yt
iq = UNIQ(xtt,sort(xtt))
xtt = xtt(iq) & ytt = ytt(iq)
in = where( (xtt GT 0d) AND (ytt GT 0d), cn )
if cn GT 0 then begin
xtt = xtt(in)
ytt = ytt(in)
endif else begin
; if no good point in filter, band flux and the rest is 0
sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt)
return,sed
endelse
nfilt = n_elements(xtt)
; sort arrays
ix = SORT(xtt)
xtt = xtt(ix) & ytt = ytt(ix)
ytt = ytt / MAX(ytt)
ix = SORT(xs)
xss = xs(ix) & yss = ys(ix)
xt=xtt & yt=ytt
; shrink SED array
ix = WHERE( xss GE MIN(xtt) AND xss LE MAX(xtt), cx)
if cx GT 0 then begin
; add points at edges for interpolation
i1 = 0
if ix(0) GT i1 then i1 = ix(0)-1
i2 = n_elements(xss)-1
if ix(cx-1) LT i2 then i2 = ix(cx-1)+1
xss = xss(i1:i2)
yss = yss(i1:i2)
endif else begin
; if no good point in SED, band flux and the rest is 0
sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt)
return,sed
endelse
nsed = n_elements(xss)
; get filter sampling f_res = dx(filt)
ix = WHERE( ytt/MAX(ytt) GT vlow, cx )
f_res = ABS( MEDIAN( xtt(ix(1:cx-1))-xtt(ix(0:cx-2)) ) )
; Interpolate to compute integration
IF NFILT GE NSED THEN BEGIN
; finer filter grid: integrate band flux on filter x-grid
xi = xtt
rr = 1d ; SED sampling is filter sampling
yi = INTERPOL(yss,xss,xi)
filt = ytt
ENDIF ELSE BEGIN
; finer SED grid: integrate band flux on SED x-grid
xi = xss & yi = yss
; get SED sampling: s_res = dx(SED) and sampling ratio
s_res = ABS( MEDIAN( xi(1:nsed-1)-xi(0:nsed-2) ) )
rr = s_res / f_res
filt = INTERPOL(ytt, xtt, xi)
io = WHERE( xi GT MAX(xtt) OR xi LT MIN(xtt), co)
if co GT 0 then filt(io) = 0d
in = where( filt LT 0d, cn)
if cn GT 0 then filt(in) = 0d
ENDELSE
xt=xi & yt=filt
; Integrate on filter bandpass (log-grid)
; trapezium safe and accurate enough -
np = n_elements(xi)
nint = np
dx = xi(1:np-1) - xi(0:np-2)
yint = yi*filt/xi
; a1 = INT_TABULATED(xi,yint, /double)
yint = (yint(0:np-2) + yint(1:np-1)) / 2d0
a1 = TOTAL( yint*dx )
yint = filt/xi
; a2 = INT_TABULATED(xi,yint, /double)
yint = (yint(0:np-2) + yint(1:np-1)) / 2d0
a2 = TOTAL( yint*dx )
; get color correction, color corrected band flux and sed
y0 = INTERPOL( yi, xi, w0)
cc = a1/a2/y0
sed = y0*cc
RETURN, sed
END