OV-GSO-DC / dustem-wrapper_idl

Commit ad743c1b18f23b2cf248f66fda729538e4b607b3

Authored by Annie Hughes
1 parent da01e0aa
Exists in master

made function calls internally consistent

Showing 2 changed files with 34 additions and 186 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_create_rfield.pro
  Diff comments   View file @ad743c1
1   -FUNCTION HABING_FIELD, x, unit=unit
2   -; computes the Habing interstellar radiation field SED (ISRF)
3   -; in W/m2 (4*!pi*nu*I_nu)
4   -; (from Draine & Bertoldi 1996)
5   -; X (I): X-grid for the ISRF
6   -; UNIT (I): unit of the ISRF. Choices are
7   -; 'W': wave in um [Default]
8   -; 'F': frequency in cm-1
9   -; 'E': energy in eV
  1 +FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname, help=help
10 2  
11   - if n_elements( UNIT ) EQ 0 then unit = 'W' $
12   - else unit = strupcase( strcompress( unit,/rem) )
13   -
14   - x = double( x )
15   -
16   - CASE unit of
17   -
18   - 'W': x3 = 1.d1 * x
19   -
20   - 'F': x3 = 1.d5 / x
21   -
22   - 'E': x3 = 12.4 / x
23   -
24   - ENDCASE
25   -
26   - field = - x3^3*4.1667 + x3^2*12.5 - x3*4.3333
27   - field = 1.d-1 * 3.d8 * 1.d-14 * field
28   -
29   - i_neg = where( field LT 0.d0, c_neg )
30   - field( i_neg ) = 0.d0
31   -
32   - RETURN, field
33   -END ;FUNCTION HABING_FIELD
34   -
35   -
36   -FUNCTION MATHIS_FIELD, x, unit=unit
37   -; computes the Mathis interstellar radiation field SED (ISRF)
38   -; in W/m2 (4*!pi*nu*I_nu)
39   -; from Mathis et al. 1983, A&A 128, 212
40   -; X (I): X-grid for the ISRF
41   -; UNIT (I): unit of the ISRF. Choices are
42   -; 'W': wave in um [Default]
43   -; 'F': frequency in cm-1
44   -; 'E': energy in eV
45   -
46   - if n_elements( UNIT ) EQ 0 then unit = 'W' $
47   - else unit = strupcase( strcompress( unit,/rem) )
48   -
49   - ly_limit = 9.11267101235d-2
50   -
51   - x = double( x )
52   -
53   -;
54   -; visible part
55   -;
56   - wdil = 4.d0 * [ 4.d-13, 1.d-13, 1.d-14] ; Mathis definition
57   - ; 4*!pi*Wdil*B_nu
58   - field = !pi * x * BLACKBODY( x, [ 3.d3, 4.d3, 7.5d3 ], wdil=wdil, unit=unit )
59   -
60   -;
61   -; UV part
62   -;
63   -
64   -; first convert to (lambda / 1e3 AA)
65   - CASE unit of
66   -
67   - 'W': x3 = 1.d1 * x
68   -
69   - 'F': x3 = 1.d5 / x
70   -
71   - 'E': x3 = 12.4 / x
72   -
73   - ENDCASE
74   -
75   - il = where( x3 LE 2.46, cl )
76   - if cl GT 0 then field(il) = 0.D0
77   -
78   - il = where( x3 GE 1d1*ly_limit AND x3 LT 1.11, cl )
79   - if cl GT 0 then field(il) = 1.4817d-6 * x3(il)^(4.4172)
80   - il = where( x3 GE 1.11 AND x3 LT 1.34, cl )
81   - if cl GT 0 then field(il) = 2.0456d-6 * x3(il)
82   - il = where( x3 GE 1.34 AND x3 LT 2.46, cl )
83   - if cl GT 0 then field(il) = 3.3105d-6 * x3(il)^(-0.6678)
84   -
85   -
86   - RETURN, field
87   -END ;FUNCTION MATHIS_FIELD
88   -
89   -
90   -FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname
91   -
92   - if N_PARAMS() LT 1 then begin
  3 + if N_PARAMS() LT 1 or keyword_set(help) then begin
93 4 print,'------------------------------------------------------------------------------------------------------'
94   - print,'FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname'
  5 + print,'FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname'
95 6 print,'------------------------------------------------------------------------------------------------------'
96 7 print,''
97 8 print,' generates the radiation field for DUSTEM (ISRF.DAT)'
... ... @@ -110,7 +21,7 @@ FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=c
110 21 print,' G0 (O): factor flux(6-13.6eV) wrt. Mathis field '
111 22 print,' CHI (O): scaling factor at 100nm wrt. Mathis field'
112 23 print,''
113   - print,'Example: tt = create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname=''ISRF.DAT'')'
  24 + print,'Example: tt = dustem_create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname=''ISRF.DAT'')'
114 25 print,''
115 26 print,' Created Aug. 2009, L. Verstraete, IAS'
116 27 print,' Force the WCUT point, May 2011, LV'
... ... @@ -137,40 +48,40 @@ FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=c
137 48 nx = n_elements(x)
138 49 if x(nx-1) NE xb(1) then x(nx-1)=xb(1) ; check rounding
139 50 endif else begin
140   - print,'(W) CREATE-RFIELD : using input wave x'
  51 + print,'(W) DUSTEM_CREATE-RFIELD : using input wave x'
141 52 x = double(x)
142 53 nx = n_elements(x)
143 54 ix = WHERE( ABS(x-wcut)/x LE 0.01, cx )
144 55 if cx EQ 0 then begin
145   - print,'(W) CREATE_RFIELD: your x-grid does not contain wcut.'
  56 + print,'(W) DUSTEM_CREATE_RFIELD: your x-grid does not contain wcut.'
146 57 print,' Should be included if radiation is 0 below wcut.'
147 58 endif
148 59 endelse
149 60  
150 61 rfield = dblarr(nx)
151 62 ; get Habing for normalization
152   - rhabing = HABING_FIELD(x,unit=unit)
  63 + rhabing = DUSTEM_HABING_FIELD(x,unit=unit)
153 64 rhabing = 1.d3*x*rhabing/3.d14 ; erg/cm2/s/Hz
154 65  
155 66 ; get isrf
156   - rmathis = MATHIS_FIELD(x,unit=unit)
  67 + rmathis = DUSTEM_MATHIS_FIELD(x,unit=unit)
157 68 rmathis = 1.d3*x*rmathis/3.d14 ; erg/cm2/s/Hz
158 69  
159 70 if n_elements(ISRF) EQ 0 then ISRF=1
160 71 if ISRF EQ 1 then begin
161   - print,'(W) CREATE_RFIELD : adding Mathis ISRF'
  72 + print,'(W) DUSTEM_CREATE_RFIELD : adding Mathis ISRF'
162 73 rfield = rmathis
163 74 endif else if ISRF EQ 2 then begin
164   - print,'(W) CREATE_RFIELD : adding Habing ISRF'
  75 + print,'(W) DUSTEM_CREATE_RFIELD : adding Habing ISRF'
165 76 rfield = rhabing
166 77 endif
167 78  
168 79 ; get blackbody
169 80 ntemp = n_elements(TEMP)
170 81 if ntemp GT 0 then begin
171   - bb = BLACKBODY( x, temp, unit=unit, wdil=wdil )
  82 + bb = DUSTEM_BLACKBODY( x, temp, unit=unit, wdil=wdil )
172 83 bb = bb * 1.d3 * !pi * x^2 / 3.d14 ; erg/cm2/s/Hz
173   - print,'(W) CREATE_RFIELD : adding BB with T= ',temp, format='(A38,10(1E10.4,1x))'
  84 + print,'(W) DUSTEM_CREATE_RFIELD : adding BB with T= ',temp, format='(A38,10(1E10.4,1x))'
174 85 print,' dilution factor wdil= ',wdil, format='(A38,10(1E10.4,1x))'
175 86 endif else bb=0.D0
176 87 rfield = rfield + bb
... ... @@ -189,7 +100,7 @@ FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=c
189 100 s2 = TOTAL( (rr(1:cg-1)+rr(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
190 101 g0 = s1/s2
191 102 endif else g0 = 0.d0
192   - print,'(W) CREATE_RFIELD : G0 =',g0
  103 + print,'(W) DUSTEM_CREATE_RFIELD : G0 =',g0
193 104  
194 105 ; get chi
195 106 ig = where( abs(2.*(x-x1)/(x+x1)) LE 1.d-1, cg)
... ... @@ -199,7 +110,7 @@ FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=c
199 110 ig = ig(0)
200 111 chi = rfield(ig) / rmathis(ig)
201 112 endif else chi = 0
202   - print,'(W) CREATE_RFIELD : chi =',chi
  113 + print,'(W) DUSTEM_CREATE_RFIELD : chi =',chi
203 114  
204 115 ; write file
205 116 if n_elements(FNAME) NE 0 then begin
... ... @@ -223,7 +134,7 @@ FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=c
223 134 printf, iu, x(i), rfield(i), format='(2(1E13.6,2x))'
224 135 endfor
225 136 FREE_LUN, iu
226   - print,'(W) CREATE_RFIELD: radiation field written in ', strtrim(fname,2)
  137 + print,'(W) DUSTEM_CREATE_RFIELD: radiation field written in ', strtrim(fname,2)
227 138 ENDIF
228 139  
229 140 RETURN, RFIELD
... ...
src/idl/dustem_create_vdist.pro
  Diff comments   View file @ad743c1
1   -FUNCTION PLAW, x, par
2   -; generates a volume normalized power law x^par(0)
3   -; returns distribution in nr of grains : dn/da
4   -; x : grain size
5   -; par(0) : power law index
6   -; par(1) : VOLUME normalization
7   -; par(2) : curvature parameter beta
8   -; par(3) : large size threshold At
9   -
10   - np = n_elements(x)
11   - y = x^par(0)
12   -; curvature term
13   - if ((par(2) NE 0) AND (par(3) NE 0)) then begin
14   - psgn = par(2)/ABS(par(2))
15   - y = y * ( 1.d0 + ABS(par(2))*x/par(3) )^psgn
16   - endif
17   - vy = x^4 * y
18   - dx = ALOG(x(1:np-1)) - ALOG(x(0:np-2))
19   - yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
20   - y = par(1) * y / yi
21   -RETURN, y
22   -END
23   -
24   -
25   -FUNCTION LOGN, x, par
26   -; generates a volume normalized log-normal law
27   -; returns distribution in nr of grains : dn/da
28   -; x : grain size
29   -; par(0) : centroid of log-normal
30   -; par(1) : sigma of log-normal
31   -; par(2) : VOLUME normalization
32   -
33   - np = n_elements(x)
34   - x0 = par(0)
35   - sigma = par(1)
36   - y = exp(- 0.5 * ( alog(x/x0) / sigma )^2 ) / x
37   - vy = x^4 * y
38   - xm = par(0) * exp( -par(1)^2 ) ; x of max in dn/da
39   - print,'(W) LOGN: dn/da max @',xm*1e7,' nm'
40   - dx = alog(x(1:np-1)) - alog(x(0:np-2))
41   - yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
42   - y = par(2) * y / yi
43   -RETURN, y
44   -END
45   -
46   -FUNCTION CUT_OFF, x, par
47   -; generates the large size cut-off
48   -; from Weingartner & Draine 2001
49   -; x : grain size
50   -; par(0) : threshold for cut-off (At)
51   -; par(1) : shape As (roughly the size where cut_off=0.5)
52   -
53   - y = 0.d0*x + 1.d0
54   - ix = WHERE( x GT par(0), cnt )
55   - if CNT GT 0 then begin
56   - y(ix) = exp( -( (x(ix)-par(0)) / par(1))^3. )
57   - endif else begin
58   - print,'(W) CUT_OFF: no size larger than At found'
59   - endelse
60   -RETURN, y
61   -END
62   -
63   -
64   -FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $
  1 +FUNCTION DUSTEM_CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $
65 2 norm=norm, fname=fname, C_WD=C_WD
66 3  
67 4 ;
... ... @@ -69,7 +6,7 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
69 6 ;
70 7 if N_PARAMS() LT 1 then begin
71 8 print,'----------------------------------------------------------------------------------------------------'
72   - print,'FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $'
  9 + print,'FUNCTION DUSTEM_CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $'
73 10 print,' norm=norm, fname=fname, C_WD=C_WD'
74 11 print,'----------------------------------------------------------------------------------------------------'
75 12 print,''
... ... @@ -82,8 +19,8 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
82 19 print,' MFRAC(I): mass fraction if composite'
83 20 print,' AG (O): size grid in cm (output)'
84 21 print,' NS (I): nr of sizes [Default = 10]'
85   - print,' SLAW (I): array of size dist. law of grains ''PLAW'' or ''LOGN'''
86   - print,' [Default = ''PLAW'' with index -3.5]'
  22 + print,' SLAW (I): array of size dist. law of grains ''DUSTEM_PLAW_VDIST'' or ''DUSTEM_LOGN_VDIST'''
  23 + print,' [Default = ''DUSTEM_PLAW_VDIST'' with index -3.5]'
87 24 print,' PAR (I): array(nlaw,npar) parameters for the size dist law: '
88 25 print,' [index,integral,beta,At] for PLAW and [center,width,integral] for LOGN'
89 26 print,' CUTOF(I): parameters for large size cut-off in microns. Default is [0.0107,0428].'
... ... @@ -94,11 +31,11 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
94 31 print,''
95 32 print,' Example 1: a power-law'
96 33 print,' par=dblarr(1,4) & par(0,0)=-3.21 & par(0,1)=1. & par(0,2)=0.3 & par(0,3)=0.164e-4'
97   - print,' sd = CREATE_VDIST( [3.e-8,2.5e-5], 3.3, slaw=''plaw'', par=par )'
  34 + print,' sd = DUSTEM_CREATE_VDIST( [3.e-8,2.5e-5], 3.3, slaw=''dustem_plaw_vdist'', par=par )'
98 35 print,''
99 36 print,' Example 2: a log-normal'
100 37 print,' par=dblarr(1,3) & par(0,0)=4.d-8 & par(0,1)=0.2 & par(0,2)=0.3 '
101   - print,' sd = CREATE_VDIST( [3.e-8,2.5e-5], 2.25, slaw=''logn'', par=par )'
  38 + print,' sd = DUSTEM_CREATE_VDIST( [3.e-8,2.5e-5], 2.25, slaw=''dustem_logn_vdist'', par=par )'
102 39 print,''
103 40 print,' Created March 2009, L. Verstraete, IAS'
104 41 print,''
... ... @@ -108,25 +45,25 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
108 45  
109 46 ; inits
110 47 if n_elements( AR ) EQ 0 then begin
111   - print,'(F) CREATE_VDIST: you must define a size range'
  48 + print,'(F) DUSTEM_CREATE_VDIST: you must define a size range'
112 49 RETURN,0
113 50 endif
114 51 if n_elements( RHO ) EQ 0 then begin
115   - print,'(F) CREATE_VDIST: you must define a grain mass density'
  52 + print,'(F) DUSTEM_CREATE_VDIST: you must define a grain mass density'
116 53 RETURN,0
117 54 endif
118 55 if n_elements( NS ) EQ 0 then ns=10
119 56 if n_elements(slaw) EQ 0 then begin
120   - slaw=['PLAW']
  57 + slaw=['DUSTEM_PLAW_VDIST']
121 58 par = dblarr(1,4)
122 59 par(0) = -3.5
123 60 par(1) = 1.d0
124 61 endif else begin
125 62 slaw=[ strupcase(strtrim(slaw,2)) ]
126 63 for i=0,n_elements( slaw )-1 do begin
127   - if (SLAW(i) NE 'PLAW') AND (SLAW(i) NE 'LOGN') then begin
128   - print,'(F) CREATE_VDIST: undefined law ',slaw(i)
129   - print,' only ''PLAW'' and ''LOGN'' '
  64 + if (SLAW(i) NE 'DUSTEM_PLAW_VDIST') AND (SLAW(i) NE 'DUSTEM_LOGN_VDIST') then begin
  65 + print,'(F) DUSTEM_CREATE_VDIST: undefined law ',slaw(i)
  66 + print,' only ''DUSTEM_PLAW_VDIST'' and ''DUSTEM_LOGN_VDIST'' '
130 67 return,0
131 68 endif
132 69 endfor
... ... @@ -142,7 +79,7 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
142 79  
143 80 ; WD01 style for carbon adapted to match cirrus emission
144 81 if keyword_set( C_WD ) then begin
145   - slaw = [ 'LOGN', 'LOGN', 'PLAW' ]
  82 + slaw = [ 'DUSTEM_LOGN_VDIST', 'DUSTEM_LOGN_VDIST', 'DUSTEM_PLAW_VDIST' ]
146 83 par = dblarr(3,4)
147 84 if TOTAL(CUTOF) NE 0 then cutof = [ 0.0107, 0.170 ] * 1.d-4
148 85  
... ... @@ -173,7 +110,7 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
173 110  
174 111 ; volume distribution
175 112 if n_elements( PAR ) EQ 0 then begin
176   - print,'(F) CREATE_VDIST: PAR array not defined'
  113 + print,'(F) DUSTEM_CREATE_VDIST: PAR array not defined'
177 114 RETURN, 0
178 115 endif
179 116 vdist = 0.d0
... ... @@ -182,20 +119,20 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
182 119 vdist = vdist + ag^(4.d0) * CALL_FUNCTION( slaw(i),ag,pp )
183 120 endfor
184 121 if TOTAL(CUTOF) NE 0 then begin
185   - print,'(W) CREATE_VDIST : cut-off applied, size and scale are ',cutof
186   - vdist = vdist * CUT_OFF(ag, cutof)
  122 + print,'(W) DUSTEM_CREATE_VDIST : cut-off applied, size and scale are ',cutof
  123 + vdist = vdist * DUSTEM_CUT_OFF(ag, cutof)
187 124 endif
188 125  
189 126 ; normalize
190 127 fac = TOTAL( vdist(1:ns-1)+vdist(0:ns-2) ) * 0.5 * da
191 128 vdist = norm * vdist / fac
192   - print,'(W) CREATE_VDIST: normalization factor ',fac / norm
  129 + print,'(W) DUSTEM_CREATE_VDIST: normalization factor ',fac / norm
193 130 yr = [1.d-4,1] * max(vdist)
194 131 plot_oo,1.e7*ag,vdist,xtit='a (nm)',yr=yr,ytit='Normalized a!u4!ndn/da', ps=-1
195 132  
196 133 ; effective density
197 134 rho_eff = TOTAL( fv*rho )
198   - print,'(W) CREATE_VDIST : effective mass density ',rho_eff,' g/cm3'
  135 + print,'(W) DUSTEM_CREATE_VDIST : effective mass density ',rho_eff,' g/cm3'
199 136  
200 137 ; write in fname
201 138 if n_elements(fname) NE 0 then begin
... ... @@ -217,7 +154,7 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par
217 154 printf, iu, ag(i), da, vdist(i,*), rho_eff, fv, format='(20(e11.4,1x))'
218 155 endfor
219 156 FREE_LUN, iu
220   - print,'(W) CREATE_VDIST: mass distribution written in ', strtrim(fname,2)
  157 + print,'(W) DUSTEM_CREATE_VDIST: mass distribution written in ', strtrim(fname,2)
221 158 endif
222 159  
223 160 RETURN, vdist
... ...