OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_plugin_mbbdy_isrf.pro 2.82 KB
Edit Raw Normal View History
759a527d   Ilyes Choubani   general update 
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
 
FUNCTION dustem_plugin_mbbdy_isrf, key=key,val=val,scope=scope,paramtag=paramtag,help=help

;+
; NAME:
;    dustem_plugin_mbbdy
; PURPOSE:
;    DUSTEM plugin to compute modified black-body emission
; CATEGORY:
;    DUSTEM Wrapper
; 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:
;    None
; PROCEDURE:
;    This is a dustem plugin
; EXAMPLES
;
; MODIFICATION HISTORY:
;    Written by JPB 
;-

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)
b0cabd3c   Ilyes Choubani   general update 
81
759a527d   Ilyes Choubani   general update 
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
 
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_QSED+ADD_USED'


the_end:
RETURN,output
  
END