From ad743c1b18f23b2cf248f66fda729538e4b607b3 Mon Sep 17 00:00:00 2001 From: Annie Hughes Date: Mon, 17 Oct 2022 19:04:46 +0100 Subject: [PATCH] made function calls internally consistent --- src/idl/dustem_create_rfield.pro | 119 +++++++++++++++-------------------------------------------------------------------------------------------------------- src/idl/dustem_create_vdist.pro | 101 +++++++++++++++++++---------------------------------------------------------------------------------- 2 files changed, 34 insertions(+), 186 deletions(-) diff --git a/src/idl/dustem_create_rfield.pro b/src/idl/dustem_create_rfield.pro index 7aa4f8f..9097180 100755 --- a/src/idl/dustem_create_rfield.pro +++ b/src/idl/dustem_create_rfield.pro @@ -1,97 +1,8 @@ -FUNCTION 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 +FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname, help=help - 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 HABING_FIELD - - -FUNCTION 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 * 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 MATHIS_FIELD - - -FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname - - if N_PARAMS() LT 1 then begin + if N_PARAMS() LT 1 or keyword_set(help) then begin print,'------------------------------------------------------------------------------------------------------' - print,'FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname' + print,'FUNCTION DUSTEM_CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname' print,'------------------------------------------------------------------------------------------------------' print,'' 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 print,' G0 (O): factor flux(6-13.6eV) wrt. Mathis field ' print,' CHI (O): scaling factor at 100nm wrt. Mathis field' print,'' - print,'Example: tt = create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname=''ISRF.DAT'')' + print,'Example: tt = dustem_create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname=''ISRF.DAT'')' print,'' print,' Created Aug. 2009, L. Verstraete, IAS' 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 nx = n_elements(x) if x(nx-1) NE xb(1) then x(nx-1)=xb(1) ; check rounding endif else begin - print,'(W) CREATE-RFIELD : using input wave x' + print,'(W) DUSTEM_CREATE-RFIELD : using input wave x' x = double(x) nx = n_elements(x) ix = WHERE( ABS(x-wcut)/x LE 0.01, cx ) if cx EQ 0 then begin - print,'(W) CREATE_RFIELD: your x-grid does not contain wcut.' + print,'(W) DUSTEM_CREATE_RFIELD: your x-grid does not contain wcut.' print,' Should be included if radiation is 0 below wcut.' endif endelse rfield = dblarr(nx) ; get Habing for normalization - rhabing = HABING_FIELD(x,unit=unit) + rhabing = DUSTEM_HABING_FIELD(x,unit=unit) rhabing = 1.d3*x*rhabing/3.d14 ; erg/cm2/s/Hz ; get isrf - rmathis = MATHIS_FIELD(x,unit=unit) + rmathis = DUSTEM_MATHIS_FIELD(x,unit=unit) rmathis = 1.d3*x*rmathis/3.d14 ; erg/cm2/s/Hz if n_elements(ISRF) EQ 0 then ISRF=1 if ISRF EQ 1 then begin - print,'(W) CREATE_RFIELD : adding Mathis ISRF' + print,'(W) DUSTEM_CREATE_RFIELD : adding Mathis ISRF' rfield = rmathis endif else if ISRF EQ 2 then begin - print,'(W) CREATE_RFIELD : adding Habing ISRF' + print,'(W) DUSTEM_CREATE_RFIELD : adding Habing ISRF' rfield = rhabing endif ; get blackbody ntemp = n_elements(TEMP) if ntemp GT 0 then begin - bb = BLACKBODY( x, temp, unit=unit, wdil=wdil ) + bb = DUSTEM_BLACKBODY( x, temp, unit=unit, wdil=wdil ) bb = bb * 1.d3 * !pi * x^2 / 3.d14 ; erg/cm2/s/Hz - print,'(W) CREATE_RFIELD : adding BB with T= ',temp, format='(A38,10(1E10.4,1x))' + print,'(W) DUSTEM_CREATE_RFIELD : adding BB with T= ',temp, format='(A38,10(1E10.4,1x))' print,' dilution factor wdil= ',wdil, format='(A38,10(1E10.4,1x))' endif else bb=0.D0 rfield = rfield + bb @@ -189,7 +100,7 @@ FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=c s2 = TOTAL( (rr(1:cg-1)+rr(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 ) g0 = s1/s2 endif else g0 = 0.d0 - print,'(W) CREATE_RFIELD : G0 =',g0 + print,'(W) DUSTEM_CREATE_RFIELD : G0 =',g0 ; get chi 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 ig = ig(0) chi = rfield(ig) / rmathis(ig) endif else chi = 0 - print,'(W) CREATE_RFIELD : chi =',chi + print,'(W) DUSTEM_CREATE_RFIELD : chi =',chi ; write file 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 printf, iu, x(i), rfield(i), format='(2(1E13.6,2x))' endfor FREE_LUN, iu - print,'(W) CREATE_RFIELD: radiation field written in ', strtrim(fname,2) + print,'(W) DUSTEM_CREATE_RFIELD: radiation field written in ', strtrim(fname,2) ENDIF RETURN, RFIELD diff --git a/src/idl/dustem_create_vdist.pro b/src/idl/dustem_create_vdist.pro index 366a222..d4861d6 100755 --- a/src/idl/dustem_create_vdist.pro +++ b/src/idl/dustem_create_vdist.pro @@ -1,67 +1,4 @@ -FUNCTION PLAW, 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 - - -FUNCTION LOGN, 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) LOGN: 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 - -FUNCTION 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) CUT_OFF: no size larger than At found' - endelse -RETURN, y -END - - -FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $ +FUNCTION DUSTEM_CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $ norm=norm, fname=fname, C_WD=C_WD ; @@ -69,7 +6,7 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par ; if N_PARAMS() LT 1 then begin print,'----------------------------------------------------------------------------------------------------' - print,'FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $' + print,'FUNCTION DUSTEM_CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $' print,' norm=norm, fname=fname, C_WD=C_WD' print,'----------------------------------------------------------------------------------------------------' print,'' @@ -82,8 +19,8 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par print,' MFRAC(I): mass fraction if composite' print,' AG (O): size grid in cm (output)' print,' NS (I): nr of sizes [Default = 10]' - print,' SLAW (I): array of size dist. law of grains ''PLAW'' or ''LOGN''' - print,' [Default = ''PLAW'' with index -3.5]' + print,' SLAW (I): array of size dist. law of grains ''DUSTEM_PLAW_VDIST'' or ''DUSTEM_LOGN_VDIST''' + print,' [Default = ''DUSTEM_PLAW_VDIST'' with index -3.5]' print,' PAR (I): array(nlaw,npar) parameters for the size dist law: ' print,' [index,integral,beta,At] for PLAW and [center,width,integral] for LOGN' 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 print,'' print,' Example 1: a power-law' 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' - print,' sd = CREATE_VDIST( [3.e-8,2.5e-5], 3.3, slaw=''plaw'', par=par )' + print,' sd = DUSTEM_CREATE_VDIST( [3.e-8,2.5e-5], 3.3, slaw=''dustem_plaw_vdist'', par=par )' print,'' print,' Example 2: a log-normal' print,' par=dblarr(1,3) & par(0,0)=4.d-8 & par(0,1)=0.2 & par(0,2)=0.3 ' - print,' sd = CREATE_VDIST( [3.e-8,2.5e-5], 2.25, slaw=''logn'', par=par )' + print,' sd = DUSTEM_CREATE_VDIST( [3.e-8,2.5e-5], 2.25, slaw=''dustem_logn_vdist'', par=par )' print,'' print,' Created March 2009, L. Verstraete, IAS' print,'' @@ -108,25 +45,25 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par ; inits if n_elements( AR ) EQ 0 then begin - print,'(F) CREATE_VDIST: you must define a size range' + print,'(F) DUSTEM_CREATE_VDIST: you must define a size range' RETURN,0 endif if n_elements( RHO ) EQ 0 then begin - print,'(F) CREATE_VDIST: you must define a grain mass density' + print,'(F) DUSTEM_CREATE_VDIST: you must define a grain mass density' RETURN,0 endif if n_elements( NS ) EQ 0 then ns=10 if n_elements(slaw) EQ 0 then begin - slaw=['PLAW'] + slaw=['DUSTEM_PLAW_VDIST'] par = dblarr(1,4) par(0) = -3.5 par(1) = 1.d0 endif else begin slaw=[ strupcase(strtrim(slaw,2)) ] for i=0,n_elements( slaw )-1 do begin - if (SLAW(i) NE 'PLAW') AND (SLAW(i) NE 'LOGN') then begin - print,'(F) CREATE_VDIST: undefined law ',slaw(i) - print,' only ''PLAW'' and ''LOGN'' ' + if (SLAW(i) NE 'DUSTEM_PLAW_VDIST') AND (SLAW(i) NE 'DUSTEM_LOGN_VDIST') then begin + print,'(F) DUSTEM_CREATE_VDIST: undefined law ',slaw(i) + print,' only ''DUSTEM_PLAW_VDIST'' and ''DUSTEM_LOGN_VDIST'' ' return,0 endif endfor @@ -142,7 +79,7 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par ; WD01 style for carbon adapted to match cirrus emission if keyword_set( C_WD ) then begin - slaw = [ 'LOGN', 'LOGN', 'PLAW' ] + slaw = [ 'DUSTEM_LOGN_VDIST', 'DUSTEM_LOGN_VDIST', 'DUSTEM_PLAW_VDIST' ] par = dblarr(3,4) if TOTAL(CUTOF) NE 0 then cutof = [ 0.0107, 0.170 ] * 1.d-4 @@ -173,7 +110,7 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par ; volume distribution if n_elements( PAR ) EQ 0 then begin - print,'(F) CREATE_VDIST: PAR array not defined' + print,'(F) DUSTEM_CREATE_VDIST: PAR array not defined' RETURN, 0 endif vdist = 0.d0 @@ -182,20 +119,20 @@ FUNCTION CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par vdist = vdist + ag^(4.d0) * CALL_FUNCTION( slaw(i),ag,pp ) endfor if TOTAL(CUTOF) NE 0 then begin - print,'(W) CREATE_VDIST : cut-off applied, size and scale are ',cutof - vdist = vdist * CUT_OFF(ag, cutof) + print,'(W) DUSTEM_CREATE_VDIST : cut-off applied, size and scale are ',cutof + vdist = vdist * DUSTEM_CUT_OFF(ag, cutof) endif ; normalize fac = TOTAL( vdist(1:ns-1)+vdist(0:ns-2) ) * 0.5 * da vdist = norm * vdist / fac - print,'(W) CREATE_VDIST: normalization factor ',fac / norm + print,'(W) DUSTEM_CREATE_VDIST: normalization factor ',fac / norm yr = [1.d-4,1] * max(vdist) plot_oo,1.e7*ag,vdist,xtit='a (nm)',yr=yr,ytit='Normalized a!u4!ndn/da', ps=-1 ; effective density rho_eff = TOTAL( fv*rho ) - print,'(W) CREATE_VDIST : effective mass density ',rho_eff,' g/cm3' + print,'(W) DUSTEM_CREATE_VDIST : effective mass density ',rho_eff,' g/cm3' ; write in fname 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 printf, iu, ag(i), da, vdist(i,*), rho_eff, fv, format='(20(e11.4,1x))' endfor FREE_LUN, iu - print,'(W) CREATE_VDIST: mass distribution written in ', strtrim(fname,2) + print,'(W) DUSTEM_CREATE_VDIST: mass distribution written in ', strtrim(fname,2) endif RETURN, vdist -- libgit2 0.21.2