dustem_plugin_mbbdy_isrf.pro
3.05 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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
FUNCTION dustem_plugin_mbbdy_isrf, key=key,val=val,scope=scope,paramtag=paramtag,help=help
;+
; NAME:
; dustem_plugin_mbbdy_isrf
;
; PURPOSE:
; DUSTEM plugin to compute modified black-body emission
;
; CATEGORY:
; DUSTEM Wrapper, Plugins, Lowlevel
;
; CALLING SEQUENCE:
; output=dustem_plugin_mbbdy([,key=][,val=])
;
; INPUTS:
; None
;
; OPTIONAL INPUT PARAMETERS:
; key = input parameter number
; val = input parameter value
;
; OUTPUTS:
; output = MBB spectrum (on dustem wavelengths) (output[*,0] is Stokes I, output[*,1] Stokes Q and output[*,2] Stokes U)
;
; OPTIONAL OUTPUT PARAMETERS:
; None
;
; ACCEPTED KEY-WORDS:
; help = if set, print this help
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; None
;
; RESTRICTIONS:
; The DustEM fortran code must be installed
; The DustEMWrap idl code must be installed
;
; PROCEDURE:
; This is a dustem plugin
;
; EXAMPLES:
;
; MODIFICATION HISTORY:
; Written by JPB
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_plugin_mbbdy_isrf'
output=0.
goto,the_end
ENDIF
;default values of input parameters
amp = 0.5 ; Amplitude of the modified black-body emission
temp = 40 ; Temperature of the dust species
beta = 1.8 ; Emissivity index
lambref = 250 ;microns
;default polarization values
smallp=0.
psi=0.
paramtag=['Amp','beta (Plaw_index)','p','Psi (deg)']
IF keyword_set(key) THEN BEGIN
ind1=where(key EQ 1,count1) ; reading amp
ind2=where(key EQ 2,count2) ; reading beta
;In the case this plugin is used in polarization
ind3=where(key EQ 3,count3) ; reading smallp
ind4=where(key EQ 4,count4) ; reading psi
IF count1 NE 0 then amp=val[ind1[0]] ;else default value
IF count2 NE 0 then beta=val[ind2[0]] ; // //
IF count3 NE 0 then smallp=val[ind3[0]] ; // //
IF count4 NE 0 then psi=val[ind4[0]] ; // //
ENDIF
;********YOU NEEED TO MAKE SURE THAT LAMBIR AND MA_ISRF.LAMBISRF ARE THE SAME********
;************************************************************************************
;integrating the mathis isrf
ma_isrf_dat=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
ma_isrf=dustem_read_isrf(ma_isrf_dat)
un=uniq(ma_isrf.isrf)
int_mathis=INT_TABULATED((ma_isrf.lambisrf)[un],(ma_isrf.isrf)[un])
;integrating the composite isrf
int_isrf=INT_TABULATED((ma_isrf.lambisrf)[un],((*(*!dustem_plugin).stellar_population))[un])
;lambir=dustem_get_wavelengths()
lambir=dustem_get_wavelengths()
Nwavs=n_elements(lambir)
tau_o_nh = 1.00E-25*(*!dustem_HCD)*(lambir/lambref)^(-1.*beta)
;fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/ma_isrf.lambisrf)*1.e20/1.e7
tempref=19.00 ;K
temp=tempref*(int_isrf/int_mathis)^(1./(4+beta))
;stop
spec=amp*tau_o_nh*dustem_planck_function(temp,lambir)
;stop
output=fltarr(Nwavs,3)
output[*,0]=spec
;polarization this actually need to
polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
output[*,1]=Q
output[*,2]=U
scope='ADD_SED+ADD_POLSED'
the_end:
RETURN,output
END