first commit of fortran helper routines

Annie Hughes
1 parent 40772701
Showing 5 changed files with 143 additions and 0 deletions Show diff stats
src/idl/dustem_cut_off.pro
src/idl/dustem_habing_field.pro
src/idl/dustem_logn_vdist.pro
src/idl/dustem_mathis_field.pro
src/idl/dustem_plaw_vdist.pro
@@ -0,0 +1,16 @@
+FUNCTION DUSTEM_CUT_OFF, x, par
+; generates the large size cut-off 
+; from Weingartner & Draine 2001
+; x : grain size 
+; par(0) : threshold for cut-off (At)
+; par(1) : shape As (roughly the size where cut_off=0.5)
+
+ y = 0.d0*x + 1.d0
+ ix = WHERE( x GT par(0), cnt )
+ if CNT GT 0 then begin 
+     y(ix) = exp( -( (x(ix)-par(0)) / par(1))^3. )
+ endif else begin 
+     print,'(W)  DUSTEM_CUT_OFF: no size larger than At found'
+ endelse
+RETURN, y
+END
@@ -0,0 +1,33 @@
+FUNCTION DUSTEM_HABING_FIELD, x, unit=unit
+; computes the Habing interstellar radiation field SED (ISRF)
+; in W/m2 (4*!pi*nu*I_nu)
+; (from Draine & Bertoldi 1996)
+; X	(I): X-grid for the ISRF
+; UNIT	(I): unit of the ISRF. Choices are 
+;	     'W': wave in um [Default]
+;	     'F': frequency in cm-1
+;	     'E': energy in eV 
+
+ if n_elements( UNIT ) EQ 0 then unit = 'W' $
+ else unit = strupcase( strcompress( unit,/rem) )
+
+ x = double( x )
+
+ CASE unit of 
+
+   'W': x3 = 1.d1 * x 
+
+   'F': x3 = 1.d5 / x
+
+   'E': x3 = 12.4 / x
+
+ ENDCASE
+
+ field = - x3^3*4.1667 + x3^2*12.5 - x3*4.3333 
+ field = 1.d-1 * 3.d8 * 1.d-14 * field 
+
+ i_neg = where( field LT 0.d0, c_neg )
+ field( i_neg ) = 0.d0
+
+ RETURN, field
+END                             ;FUNCTION DUSTEM_HABING_FIELD
@@ -0,0 +1,20 @@
+FUNCTION DUSTEM_LOGN_VDIST, x, par
+; generates a volume normalized log-normal law 
+; returns distribution in nr of grains : dn/da
+; x : grain size 
+; par(0) : centroid of log-normal
+; par(1) : sigma of log-normal
+; par(2) : VOLUME normalization
+
+ np = n_elements(x)
+ x0 = par(0)
+ sigma = par(1)
+ y = exp(- 0.5 * ( alog(x/x0) / sigma )^2 ) / x
+ vy = x^4 * y
+ xm = par(0) * exp( -par(1)^2 ) ; x of max in dn/da
+ print,'(W) DUSTEM_LOGN_VDIST: dn/da max @',xm*1e7,' nm'
+ dx = alog(x(1:np-1)) - alog(x(0:np-2)) 
+ yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
+ y = par(2) * y / yi
+RETURN, y
+END 
@@ -0,0 +1,52 @@
+FUNCTION DUSTEM_MATHIS_FIELD, x, unit=unit
+; computes the Mathis interstellar radiation field SED (ISRF)
+; in W/m2 (4*!pi*nu*I_nu)
+; from Mathis et al. 1983, A&A 128, 212
+; X	(I): X-grid for the ISRF
+; UNIT	(I): unit of the ISRF. Choices are 
+;	     'W': wave in um [Default]
+;	     'F': frequency in cm-1
+;	     'E': energy in eV
+
+ if n_elements( UNIT ) EQ 0 then unit = 'W' $
+ else unit = strupcase( strcompress( unit,/rem) )
+
+ ly_limit = 9.11267101235d-2
+
+ x = double( x )
+
+;
+; visible part 
+;
+ wdil = 4.d0 * [ 4.d-13, 1.d-13, 1.d-14]	; Mathis definition
+						; 4*!pi*Wdil*B_nu
+ field = !pi * x * DUSTEM_BLACKBODY( x, [ 3.d3, 4.d3, 7.5d3 ], wdil=wdil, unit=unit )
+
+;
+; UV part
+;
+
+; first convert to (lambda / 1e3 AA)
+ CASE unit of 
+
+   'W': x3 = 1.d1 * x 
+
+   'F': x3 = 1.d5 / x
+
+   'E': x3 = 12.4 / x
+
+ ENDCASE
+ 
+ il = where( x3 LE 2.46, cl ) 
+ if cl GT 0 then field(il) = 0.D0
+
+ il = where( x3 GE 1d1*ly_limit AND x3 LT 1.11, cl )
+ if cl GT 0 then field(il) = 1.4817d-6 * x3(il)^(4.4172)
+ il = where( x3 GE 1.11 AND x3 LT 1.34, cl ) 
+ if cl GT 0 then field(il) = 2.0456d-6 * x3(il)
+ il = where( x3 GE 1.34 AND x3 LT 2.46, cl )
+ if cl GT 0 then field(il) = 3.3105d-6 * x3(il)^(-0.6678)
+  
+
+ RETURN, field
+END                             ;FUNCTION DUSTEM_MATHIS_FIELD
@@ -0,0 +1,22 @@
+FUNCTION DUSTEM_PLAW_VDIST, x, par 
+; generates a volume normalized power law x^par(0)
+; returns distribution in nr of grains : dn/da
+; x : grain size 
+; par(0) : power law index
+; par(1) : VOLUME normalization 
+; par(2) : curvature parameter beta 
+; par(3) : large size threshold At
+
+ np = n_elements(x)
+ y = x^par(0)
+; curvature term
+ if ((par(2) NE 0) AND (par(3) NE 0)) then begin 
+     psgn = par(2)/ABS(par(2))
+     y = y * ( 1.d0 + ABS(par(2))*x/par(3) )^psgn
+ endif 
+ vy = x^4 * y 
+ dx = ALOG(x(1:np-1)) - ALOG(x(0:np-2)) 
+ yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
+ y = par(1) * y / yi
+RETURN, y
+END