diff --git a/src/idl/dustem_add_inst.pro b/src/idl/dustem_add_inst.pro new file mode 100644 index 0000000..85fe70c --- /dev/null +++ b/src/idl/dustem_add_inst.pro @@ -0,0 +1,30 @@ +FUNCTION DUSTEM_ADD_INST, sf, tag, nsz, ntrans=ntrans +; fills in data for existing instruments in SMDAT structure +; or add a new instrument to SMDAT structure +; +; SF (I): SMDAT input structure +; TAG (I): name of instrument to add in +; NSZ (I): array(2) of sizes [nr of data pts, nr of grain types] +; NTRANS (I): nr of points for transmission (bad flux data) + + stag = TAG_NAMES(sf) + ntag = N_TAGS(sf) + tag = STRUPCASE('I_'+STRTRIM(tag,2)) + itg = WHERE( stag EQ tag, ctg ) + if n_elements(nsz) EQ 2 then begin + nx = nsz(0) & ntype = nsz(1) + endif else begin + print,'(F) DUSTEM_ADD_INST: array of sizes must be 2D' + endelse + if ctg EQ 0 then begin + if n_elements(ntrans) EQ 0 then st = CREATE_STRUCT( tag, DUSTEM_STR_INST(nx,n2=ntype) ) $ + else st = CREATE_STRUCT( tag, DUSTEM_STR_INST(nx,n2=ntype,n3=ntrans) ) + sf = CREATE_STRUCT( sf, st ) + stag = TAG_NAMES(sf) + itg = WHERE( stag EQ tag ) + sf.(itg).isel = intarr(nx) + 1 + endif else begin + print,'(F) DUSTEM_ADD_INST: INST already exists in structure' + endelse + return,sf +END diff --git a/src/idl/dustem_add_mod.pro b/src/idl/dustem_add_mod.pro new file mode 100644 index 0000000..6791f0a --- /dev/null +++ b/src/idl/dustem_add_mod.pro @@ -0,0 +1,54 @@ +FUNCTION DUSTEM_ADD_MOD, sf, tag, nsz, unit=unit +; adds a model tag to a existing SMDAT structure +; SF (I): input SMDAT structure +; TAG (I): tag of model to be added +; NSZ (I): int array(2) of sizes [nr of data pts, nr of grain types] + + tag = STRLOWCASE(STRTRIM(tag,2)) + if n_elements(nsz) GT 1 then begin + n1 = nsz(0) & n2 = nsz(1) + if n_elements(nsz) EQ 3 then n3 = nsz(2) else n3 = 0 + endif else begin + print,'(F) ADD_MOD: array of sizes must > 1D' + endelse + if TAG EQ 'emis' then begin + if n_elements(unit) EQ 0 then unit='x(microns) SED(erg/s/cm2/sr)' + s1 = { UNIT : unit, $ + X : dblarr(n1), $ ; wave + Y : dblarr(n1,n2), $ ; SED(wave, grain type) (index ntype is total) + YP : dblarr(n1,n2) } ; polarized SED(wave, grain type) (index ntype is total) + ENDIF else if TAG EQ 'ext' then begin + if n_elements(unit) EQ 0 then unit='x(microns) sigma(cm2/H)' + s1 = {UNIT : unit, $ + X : dblarr(n1), $ ; wave + Y : dblarr(n1,n2), $ ; sigma_ext(wave, grain type) (index ntype is total) + ABS : dblarr(n1,n2), $ ; sigma_abs(wave, grain type) + SCA : dblarr(n1,n2), $ ; sigma_sca(wave, grain type) + ABS_P : dblarr(n1,n2), $ ; absorption sigma_pol(wave, grain type) + SCA_P : dblarr(n1,n2), $ ; scattering sigma_pol(wave, grain type) + ALB : dblarr(n1,n2), $ ; alb(wave, grain type) + XR : 0d, $ ; ref wave + YR_ABS : dblarr(n2), $ ; tau_abs/NH @ XR + YR_SCA : dblarr(n2), $ ; tau_sca/NH @ XR + RV : 0d } +; ENDIF else if TAG EQ 'pol' then begin +; if n_elements(unit) EQ 0 then unit='x(microns) SED(erg/s/cm2/sr)' +; s1 = {UNIT : unit, $ +; X: dblarr(n1), $ ; wave +; Y: dblarr(n1,n2), $ ; polarized SED(wave, grain type) (index ntype is total) +; POL: dblarr(n1,n2) } ; sigma_pol(wave, grain type) (ABS + SCAT) + ENDIF else if tag EQ 'sdist' then begin + if n3 EQ 0 then begin + print,'(F) ADD_MOD: 3d dimension missing ' + return,0 + endif + if n_elements(unit) EQ 0 then unit='x(cm) a^4*dn/da(cm3/H)' + s1 = {UNIT : unit, $ + XTOT : dblarr(n1), $ ; grain size + YTOT : dblarr(n1,n2), $ ; size distribution + XI : dblarr(n3,n2), $ ; grain size per type + YI : dblarr(n3,n2) } ; grain size dist per type + ENDIF + sm = CREATE_STRUCT( CREATE_STRUCT('M_'+tag,s1), sf ) + return, sm +END diff --git a/src/idl/dustem_band_cc.pro b/src/idl/dustem_band_cc.pro new file mode 100644 index 0000000..60c66e1 --- /dev/null +++ b/src/idl/dustem_band_cc.pro @@ -0,0 +1,101 @@ +FUNCTION dustem_BAND_CC, w0, xs, ys, xt, yt, y0=y0, cc=cc, rr=rr, nint=nint +; returns flux in given instrument band (erg/cm2/sr) +; W0 (I): centroid of filter (any unit) +; XS (I): input x-grid for SED (same unit as W0) +; YS (I): input SED in erg/cm2/sr +; XT (I/O): x-grid of filter (same unit as W0) +; YT (I/O): filter transmission (arb. unit) +; RR (O): undersampling ratio: dx(SED)/dx(filt), usually 1. +; If RR<1 filter oversampled in flux integration +; Y0 (O): SED at band center +; CC (O): color correction +; NINT (O): nr of points in band integration + + vlow = 1d-5 + +; remove bad and multiple values in filter + xtt=xt & ytt=yt + iq = UNIQ(xtt,sort(xtt)) + xtt = xtt(iq) & ytt = ytt(iq) + in = where( (xtt GT 0d) AND (ytt GT 0d), cn ) + if cn GT 0 then begin + xtt = xtt(in) + ytt = ytt(in) + endif else begin +; if no good point in filter, band flux and the rest is 0 + sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt) + return,sed + endelse + nfilt = n_elements(xtt) + +; sort arrays + ix = SORT(xtt) + xtt = xtt(ix) & ytt = ytt(ix) + ytt = ytt / MAX(ytt) + ix = SORT(xs) + xss = xs(ix) & yss = ys(ix) + xt=xtt & yt=ytt + +; shrink SED array + ix = WHERE( xss GE MIN(xtt) AND xss LE MAX(xtt), cx) + if cx GT 0 then begin +; add points at edges for interpolation + i1 = 0 + if ix(0) GT i1 then i1 = ix(0)-1 + i2 = n_elements(xss)-1 + if ix(cx-1) LT i2 then i2 = ix(cx-1)+1 + xss = xss(i1:i2) + yss = yss(i1:i2) + endif else begin +; if no good point in SED, band flux and the rest is 0 + sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt) + return,sed + endelse + nsed = n_elements(xss) + +; get filter sampling f_res = dx(filt) + ix = WHERE( ytt/MAX(ytt) GT vlow, cx ) + f_res = ABS( MEDIAN( xtt(ix(1:cx-1))-xtt(ix(0:cx-2)) ) ) + +; Interpolate to compute integration + IF NFILT GE NSED THEN BEGIN +; finer filter grid: integrate band flux on filter x-grid + xi = xtt + rr = 1d ; SED sampling is filter sampling + yi = INTERPOL(yss,xss,xi) + filt = ytt + ENDIF ELSE BEGIN +; finer SED grid: integrate band flux on SED x-grid + xi = xss & yi = yss +; get SED sampling: s_res = dx(SED) and sampling ratio + s_res = ABS( MEDIAN( xi(1:nsed-1)-xi(0:nsed-2) ) ) + rr = s_res / f_res + filt = INTERPOL(ytt, xtt, xi) + io = WHERE( xi GT MAX(xtt) OR xi LT MIN(xtt), co) + if co GT 0 then filt(io) = 0d + in = where( filt LT 0d, cn) + if cn GT 0 then filt(in) = 0d + ENDELSE + xt=xi & yt=filt + +; Integrate on filter bandpass (log-grid) +; trapezium safe and accurate enough - + np = n_elements(xi) + nint = np + dx = xi(1:np-1) - xi(0:np-2) + yint = yi*filt/xi +; a1 = INT_TABULATED(xi,yint, /double) + yint = (yint(0:np-2) + yint(1:np-1)) / 2d0 + a1 = TOTAL( yint*dx ) + yint = filt/xi +; a2 = INT_TABULATED(xi,yint, /double) + yint = (yint(0:np-2) + yint(1:np-1)) / 2d0 + a2 = TOTAL( yint*dx ) + +; get color correction, color corrected band flux and sed + y0 = INTERPOL( yi, xi, w0) + cc = a1/a2/y0 + sed = y0*cc + + RETURN, sed +END diff --git a/src/idl/dustem_blackbody.pro b/src/idl/dustem_blackbody.pro new file mode 100644 index 0000000..57857df --- /dev/null +++ b/src/idl/dustem_blackbody.pro @@ -0,0 +1,152 @@ +FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, check=check, g0=g0, chi=chi, $ + bb_to_hab=bb_to_hab, NORM=norm + if n_params() EQ 0 then begin + print,'------------------------------------------------------------------------------------------' + print,'FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, NORM=norm, check=check, g0=g0, chi=chi' + print,'------------------------------------------------------------------------------------------' + print,' computes a (or a sum of) blackbody brightness(es) B(x,T) in W/m2/(x unit)/sr ' + print,' X (I): blackbody abscissa in cm-1, GHz, microns or eV' + print,' TEMP (I): blackbody temperature (can be an array)' + print,' WDIL (I): dilution factor for each TEMP. [Default=1].' + print,' UNIT (I): blackbody x unit ''K'' is B_nu (W/m2/cm-1/sr), x in cm-1 ' + print,' ''F'' is B_nu (W/m2/Hz/sr), x in GHz [default]' + print,' ''W'' is B_lambda (W/m2/um/sr), x in microns.' + print,' ''E'' is B_E (W/m2/eV/sr), x in eV' + print,' NORM (I): overall normalizing factor, returns NORM*BB/MAX(BB)' + print,' CHECK (O): check energy conservation for the BB (power is sigma*T^4)' + print,' CHI (O): scaling factor @ 1000 A in ISRF unit (1.6e-3 erg/cm2/s)' + print,' G0 (O): scaling factor for integrated 5-13.6 eV flux in ISRF unit.' + print,'' + print,' Created 2001, L Verstraete IAS' + print,' More units and checks, Oct 2010 LV' + print,'---------------------------------------------------------------------------------------' + return,0 + endif + nx = n_elements( X ) + nt = n_elements( TEMP ) + nw = n_elements( WDIL ) + if nw EQ 0 then begin + wdil = dblarr(nt) + 1. + endif else begin + if nw LT nt then begin + print,'DUSTEM_BLACKBODY (F): WDIL and TEMP have different sizes' + return, 0 + endif + endelse + + if n_elements( UNIT ) EQ 0 then unit='F' $ + else begin + unit = strupcase(strcompress(unit,/rem)) + if unit NE 'K' AND unit NE 'F' AND unit NE 'W' AND unit NE 'E' then begin + print,'DUSTEM_BLACKBODY (W): not a valid BB unit --> make a B_nu' + unit = 'F' + endif + endelse + + if n_elements( CHECK ) GT 1 then check = 0 + + x = double(x) + temp = double(temp) + wdil = double(wdil) + + xpi = 3.1415926535897932385d + clight = 2.9979246d08 ; m/s + kb = 1.3806503d-23 + hp = 6.62606876d-34 + hck = hp*clight / kb + ze = 1.60217653d-19 ; electron charge + hev = hp/ze ; h in eV unit + + xly = 12.4 + isrf0 = 1.6d-6 ; 4*pi*nu*I_nu (W/m2) in solar neighborhood + +; +; loop on temperatures +; + bbt = dblarr(nx) + acheck = 0. + + FOR j = 0, nt-1 do begin + + xi = [6., 13.6] + + CASE UNIT OF + + 'K': begin + be = exp( -(1d2*x) * hck / temp(j) ) + bb = be * (1d2*x)^3 / (1.d0 - be) + bb = bb * 2d0 * hp * clight^2 + bb = bb * 1d2 ; m-1 to cm-1 + xi = xi / hev/clight/1d2 + x1 = xly / hev/clight/1d2 + end + + 'F': begin + be = exp( -(1d9*x) * hp/kb/ temp(j) ) + bb = be * (1d9*x)^3 / (1.d0 - be) + bb = bb * 2d0 * hp / clight^2 + xi = xi / hev / 1d9 + x1 = xly / hev / 1d9 + end + + 'W': begin + be = exp( -hp*clight/kb/ (1d-6*x) / temp(j) ) + bb = be / (1d-6*x)^5 / (1.d0 - be) + bb = bb * 2d0 * hp * clight^2 + bb = 1d-6 * bb ; m-1 to micron-1 + xi = hev*clight*1d6 / xi + x1 = hev*clight*1d6 / xly + end + + 'E': begin + be = exp( -(ze*x)/kb / temp(j) ) + bb = be * (ze*x)^3 / (1.d0 - be) + bb = bb * 2d0 / hp^3 / clight^2 + bb = bb * ze ; J-1 to eV-1 + xi = xi + x1 = xly + end + + ENDCASE + + ck = 0. + for jx = 1L,nx-1 do $ + ck = ck + ( bb(jx)+bb(jx-1) )*( x(jx)-x(jx-1) ) / 2.d0 + ck = xpi* ck / 5.6697d-8 / temp(j)^4 + + if n_elements( CHECK ) NE 0 then begin + if check EQ 1 then $ + print,'DUSTEM_BLACKBODY (W): integrated BB('+strtrim(temp(j),2)+$ + ' K) is '+ strtrim(ck,2)+ ' times Sigma*T^4' + endif + acheck = [ acheck, ck ] + + bbt = bbt + bb*wdil(j) + + ENDFOR + + check = acheck(1:*) + +; get G0 + ig = where( x GE MIN(xi) AND x LE MAX(xi), cg) + if cg GT 0 then begin + xx=x(ig) + yy=xpi*x(ig)*bbt(ig) + s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 ) + g0 = s1/isrf0 + endif else g0 = 0.d0 + +; get chi + ig = where( abs(2.*(x-x1)/(x+x1)) LE 1.d-1, cg) + if cg GT 0 then begin + tmp = min( abs(x(ig)-x1), i_min ) + ig = ig(i_min) + ig = ig(0) + chi = xpi * x(ig) * bbt(ig) / isrf0 + endif else chi = 0 + + if n_elements( NORM ) NE 0 then bbt = norm * bbt/max(bbt) + +the_end: + RETURN, BBT +END ;FUNCTION BLACKBODY diff --git a/src/idl/dustem_chi2.pro b/src/idl/dustem_chi2.pro new file mode 100644 index 0000000..4ed6c3c --- /dev/null +++ b/src/idl/dustem_chi2.pro @@ -0,0 +1,17 @@ +FUNCTION DUSTEM_CHI2, y, model, npar, err=err +; returns the chi-square value of fit MODEL to data Y +; +; NPAR (I): nr of parameters in the fit +; ERR (I): error of each data point. Default is ERR=0.1*Y + + ny = n_elements(y) + if n_elements(err) EQ 0 OR TOTAL(err) EQ 0 then begin + err = 0.1*y + print,'(W) CHI2: error missing, set to 10 %' + endif + ndof = ny - npar ; nr of degrees of freedom + chi = TOTAL( ((y-model)/err)^2 ) + if ndof GT 0 then chi = chi / ndof + + return, chi +END diff --git a/src/idl/dustem_create_rfield.pro b/src/idl/dustem_create_rfield.pro new file mode 100755 index 0000000..7aa4f8f --- /dev/null +++ b/src/idl/dustem_create_rfield.pro @@ -0,0 +1,230 @@ +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 + + 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 + print,'------------------------------------------------------------------------------------------------------' + print,'FUNCTION 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)' + print,' in erg/cm2/s/Hz (4*!pi*I_nu)' + print,' RFIELD = ISRF + WDIL*PI*BB(TEMP) ' + print,' ISRF from Mathis et al. 1983, A&A 128, 212' + print,' NB: if blackbody is isotropic set WDIL to 4 (to get 4*!pi factor)' + print,'' + print,' TEMP (I): blackbody temperature (can be an array). If 0 only ISRF.' + print,' X (I): X-grid for the ISRF in microns. Default is 200 pts over 0.01-10^5 microns.' + print,' (including WCUT point). You want to include WCUT for accurate edges.' + print,' ISRF (I): if set to 0 no ISRF is added, 1 is Mathis (default), 2 is Habing' + print,' WDIL (I): blackbody dilution factor (can be an array)' + print,' FNAME(I): filename for ISRF.DAT' + print,' WCUT (I): for wave < wcut radiation field is 0 ' + 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,'' + print,' Created Aug. 2009, L. Verstraete, IAS' + print,' Force the WCUT point, May 2011, LV' + print,'' + print,'------------------------------------------------------------------------------------------------------' + RETURN,0 + endif + +; inits + unit='W' + ly_limit = 9.11267101235d-2 + xi = [ ly_limit, 2.07d-1 ] ; for G0 6-13.6 eV in microns + x1 = 1.d-1 ; for CHI + + if n_elements(WCUT) EQ 0 then wcut=ly_limit + if n_elements(x) EQ 0 then begin + nx = 199 + xb = [ 1.d-2, 1.d5 ] ; wave boundaries in microns + dx = (alog(xb(1))-alog(xb(0))) / (nx-1) + x = alog(xb(0)) + dindgen(nx)*dx + x = exp(x) + x = [x, wcut] ; add wcut to avoid coarse edge + x = x(UNIQ(x,SORT(x))) + 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' + 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,' 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 = 1.d3*x*rhabing/3.d14 ; erg/cm2/s/Hz + +; get isrf + rmathis = 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' + rfield = rmathis + endif else if ISRF EQ 2 then begin + print,'(W) 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 = 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,' dilution factor wdil= ',wdil, format='(A38,10(1E10.4,1x))' + endif else bb=0.D0 + rfield = rfield + bb + +; apply cut + ix = WHERE( x LT WCUT, cx ) + if cx GT 0 then rfield(ix)=0.d0 + +; get G0 + ig = where( x GE xi(0) AND x LE xi(1), cg) + if cg GT 0 then begin + xx=x(ig) + yy=rfield(ig) + rr=rmathis(ig) + s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 ) + 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 + +; get chi + ig = where( abs(2.*(x-x1)/(x+x1)) LE 1.d-1, cg) + if cg GT 0 then begin + tmp = min( abs(x(ig)-x1), i_min ) + ig = ig(i_min) + ig = ig(0) + chi = rfield(ig) / rmathis(ig) + endif else chi = 0 + print,'(W) CREATE_RFIELD : chi =',chi + +; write file + if n_elements(FNAME) NE 0 then begin + OPENW, iu, fname, /get_lun + printf, iu, '# DUSTEM: exciting radiation field featuring ' + if ISRF EQ 1 then printf, iu, '# Mathis ISRF' + if NTEMP GT 0 then begin + a = '# Blackbody with T=' + a1 = '# dilution factor wdil=' + for i = 0, ntemp - 1 do begin + a = a + string( format='(2x,1E11.4)', temp(i) ) + a1 = a1 + string( format='(2x,1E11.4)', wdil(i) ) + endfor + printf, iu, a + printf, iu, a1 + endif + printf, iu, '# Nbr of points' + printf, iu, '# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)' + printf, iu, nx, format='(i4)' + for i=0,nx-1 do begin + 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) + ENDIF + + RETURN, RFIELD +END diff --git a/src/idl/dustem_create_vdist.pro b/src/idl/dustem_create_vdist.pro new file mode 100755 index 0000000..366a222 --- /dev/null +++ b/src/idl/dustem_create_vdist.pro @@ -0,0 +1,224 @@ +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, $ + norm=norm, fname=fname, C_WD=C_WD + +; +; Doc +; + 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,' norm=norm, fname=fname, C_WD=C_WD' + print,'----------------------------------------------------------------------------------------------------' + print,'' + print,'generates a dust size distribution in volume normalized to 1 or norm' + print,'uses power law or log-normal component' + print,'' + print,' AR (I): array(2) size range in cm ' + print,' RHO (I): density (g/cm3) of grain type (array if composite)' + print,' FV (I): volume fractions if composite' + 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,' 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].' + print,' NORM (I): normalization for the size distribution' + print,' FNAME(I): file name to write size dist. in' + print,' C_WD (I): keyword to generate a C dust size dist WD01 style' + print,' ( PAH, VSG as log-normal+power law a^(-2.54) ) ' + 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,'' + 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,'' + print,' Created March 2009, L. Verstraete, IAS' + print,'' + print,'----------------------------------------------------------------------------------------------------' + RETURN, 0 + endif + +; inits + if n_elements( AR ) EQ 0 then begin + print,'(F) 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' + RETURN,0 + endif + if n_elements( NS ) EQ 0 then ns=10 + if n_elements(slaw) EQ 0 then begin + slaw=['PLAW'] + 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'' ' + return,0 + endif + endfor + endelse + nmat = n_elements( RHO ) ; nr of bulk material in composite + if n_elements( FV ) EQ 0 then begin + if NMAT EQ 1 then fv = [1.] else $ + fv = (fltarr(nmat)+1.) / nmat + mfrac = fv + endif + if n_elements( CUTOF ) EQ 0 then cutof = [ 0.0107, 0.428] * 1.d-4 + if n_elements( NORM ) EQ 0 then norm = 1.d0 + +; WD01 style for carbon adapted to match cirrus emission + if keyword_set( C_WD ) then begin + slaw = [ 'LOGN', 'LOGN', 'PLAW' ] + par = dblarr(3,4) + if TOTAL(CUTOF) NE 0 then cutof = [ 0.0107, 0.170 ] * 1.d-4 + +; PAH log-normal + par(0,0) = 4.d-8 + par(0,1) = 0.2 + par(0,2) = 0.30 + +; VSG log-normal + par(1,0) = 7.d-8 + par(1,1) = 0.4 + par(1,2) = 0.08 + +; BG power law + par(2,0) = -2.54 + par(2,1) = 2.e-5 + par(2,2) = 0.0 ;-0.165 + par(2,3) = 0.0107e-4 + endif + nlaw = n_elements( SLAW ) + +; define size grid (log) + lar = alog(ar) + da = (lar(1)-lar(0)) / (ns-1) + a = lar(0) + dindgen(ns)*da + a(ns-1) = lar(1) ; roundup error + ag = exp(a) + +; volume distribution + if n_elements( PAR ) EQ 0 then begin + print,'(F) CREATE_VDIST: PAR array not defined' + RETURN, 0 + endif + vdist = 0.d0 + for i = 0, nlaw-1 do begin + pp = REFORM( par(i,*), n_elements(par(i,*)) ) + 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) + 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 + 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' + +; write in fname + if n_elements(fname) NE 0 then begin + OPENW, iu, fname, /get_lun + printf, iu, '# Size distribution of grain species' + printf, iu, '#' + printf, iu, '# Nbr of bulk materials' + printf, iu, '# Bulk densities in g/cm3 ' + printf, iu, '# Mass fractions for each bulk' + printf, iu, '# Nbr of size bins' + printf, iu, '# [ a (cm), dloga, a^4*dn/da, rho_eff, fv ] ' + printf, iu, '# fv: volume fraction of bulks, rho_eff: volume mean density' + printf, iu, '#' + printf, iu, nmat, format='(i2)' + printf, iu, rho, format='(10(e11.4,1x))' + printf, iu, mfrac, format='(10(e11.4,1x))' + printf, iu, ns, format='(i2,1x,e11.4)' + for i=0,ns-1 do begin + 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) + endif + +RETURN, vdist +END diff --git a/src/idl/dustem_fil_chi2.pro b/src/idl/dustem_fil_chi2.pro new file mode 100644 index 0000000..3439f16 --- /dev/null +++ b/src/idl/dustem_fil_chi2.pro @@ -0,0 +1,36 @@ +FUNCTION DUSTEM_FIL_CHI2, sf, ntype=ntype +; fills in the Chi-square fields of an input SMDAT structure +; (see format in DUSTEM_GET_BAND_FLUX) +; +; SF (I): input SMDAT structure +; NTYPE (I): index of SED to be used + + if ntype EQ 0 then begin + tt = SIZE(sf.sed.y) + ntype = tt(2) + endif + ntag = N_TAGS(sf) + stag = TAG_NAMES(sf) + itg = WHERE( STRPOS(stag,'I_') GE 0, ctg ) + if ctg EQ 0 then begin + print,'(F) DUSTEM_FIL_CHI2: no instrument to fill in' + return,0 + endif + yy=0d & mm=0d & ee=0d & ift=0 + for k = 0, ctg-1 do begin + i1 = WHERE( sf.(itg(k)).isel EQ 1 AND sf.(itg(k)).err GT 0, c1) + if c1 GT 0 then begin + sf.(itg(k)).chi2 = $ + DUSTEM_CHI2(sf.(itg(k)).yd(i1), sf.(itg(k)).ym(i1,ntype-1), sf.(itg(k)).npar, err=sf.(itg(k)).err(i1)) + endif + yy = [ yy, sf.(itg(k)).yd ] + ee = [ ee, sf.(itg(k)).err] + mm = [ mm, sf.(itg(k)).ym(*,ntype-1) ] + ift = [ ift, sf.(itg(k)).isel ] + endfor + yy=yy(1:*) & ee=ee(1:*) & mm=mm(1:*) & ift=ift(1:*) + it = WHERE( ift EQ 1 AND ee GT 0, ct) + if ct GT 0 then sf.chi2 = DUSTEM_CHI2( yy(it), mm(it), sf.npar, err=ee(it)) + + return, sf +END diff --git a/src/idl/dustem_serkowski.pro b/src/idl/dustem_serkowski.pro new file mode 100644 index 0000000..440022a --- /dev/null +++ b/src/idl/dustem_serkowski.pro @@ -0,0 +1,22 @@ +FUNCTION DUSTEM_SERKOWSKI, x, ka=ka, xmax=xmax, pmax=pmax + IF N_PARAMS() EQ 0 THEN BEGIN + print,'FUNCTION SERKOWSKI, x, ka=ka, xmax=xmax, pmax=pmax ' + print,' computes the Serkowski law (Draine & Fraisse 2009)' + print,'' + print,' X (I): array(n_qabs) inverse wavenumber in 1/microns' + print,' KA (I): K factor' + print,' XMAX (I): x-position of max' + print,' PMAX (I): max polar. fraction' + ENDIF + IF n_elements(KA) EQ 0 THEN ka = 0.92 ; Draine & Fraisse (2009) + IF n_elements(XMAX) EQ 0 THEN xmax = 1.82 + IF n_elements(pMAX) EQ 0 THEN pmax = 0.03 + + yy = pmax * EXP( -ka * ALOG(xmax/x)^2 ) + x0 = 1/1.39 + y0 = pmax * EXP( -ka * ALOG(xmax/x0)^2 ) + ix = WHERE( x LE x0, cx ) + IF cx GT 0 THEN yy(ix) = y0 * (x(ix)/x0)^1.7 + + RETURN, yy +END diff --git a/src/idl/dustem_show_fortran.pro b/src/idl/dustem_show_fortran.pro new file mode 100644 index 0000000..76e8435 --- /dev/null +++ b/src/idl/dustem_show_fortran.pro @@ -0,0 +1,1061 @@ +PRO dustem_show_fortran, model=model, data_dir=data_dir,dustem_dir=dustem_dir,nh=nh, fsed=fsed, sw=sw, inst=inst, smdat=smdat, com=com, wref=wref, bbpar=bbpar, wext=wext, $ + xr=xr, yr=yr, tit=tit, CMB=cmb, COSMO=COSMO, DATA=data, SHOW=show, wn=wn, HARD=hard + +;+ +; NAME: +; dustem_show_fortran +; PURPOSE: +; shows differents results of the DUSTEM fortran calculations, +; overlaid on fiducial diffuse ISM dust observables +; CATEGORY: +; DustEM, Distributed, User-Example +; CALLING SEQUENCE: +; dustem_show_results +; INPUTS: +; DATA_DIR : path to directory containing fiducial observational +; SED files. Defaults to Data/EXAMPLE_OBSDATA/ subdirectory. +; DUSTEM_DIR : path to directory containing DustEM data files. Default +; is the !dustem_dat directory (which is defined in +; a user's idl_startup file) +; model = specifies the interstellar dust mixture used by DustEM +; 'MC10' model from Compiegne et al 2010 (default) +; 'DBP90' model from Desert et al 1990 +; 'DL01' model from Draine & Li 2001 +; 'WD01_RV5p5B' model from Weingartner & Draine 2002 with Rv=5.5 +; 'DL07' model from Draine & Li 2007 +; 'J13' model from Jones et al 2013, as updated in +; Koehler et al 2014 +; 'G17_ModelA' model A from Guillet et al (2018). Includes +; polarisation. See Tables 2 and 3 of that paper for details. +; 'G17_ModelB' model B from Guillet et al (2018) +; 'G17_ModelC' model C from Guillet et al (2018) +; 'G17_ModelD' model A from Guillet et al (2018) +; NH : float, fiducial column density for emission [Default is 1e20 cm-2]' +; FSED : string, read SED from the file FSED +; SW : integer, width of window to SMOOTH the PAH emission +; INST : string array, instruments to be shown. Default is ['DIRBE','FIRAS','HFI','WMAP'] +; COM : string to describe SMDAT structure +; WEXT : normalize the UV-extinction to be 1 at wavelength WEXT (in microns). Also applied to IR-extinction +; WREF : get the dust cross-section per H at wavelength WREF (in microns) +; BBPAR : float tuple, T and beta for modified BB to get dust opacity from FIRAS +; XR,YR : plot ranges +; TIT : title for SMDAT and plots +; CMB : keyword to overlay the CMB +; COSMO : keyword to plot long wavelengths +; DATA : if DATA=0, only the model spectrum is shown +; if DATA=1, only the observational SED is shown +; if DATA=2, both the model and data are shown +; both are shown (default)' +; SHOW : string array, defines what to plot display. Possible +; options include +; ''emis'', ''extuv'', ''extir'', ''alb'', ''sdist'', ''polext'', ''polsed'', ''align'' +; Default is [''emis'']. SHOW=0 no plot +; WN : array of window numbers for plots +; HARD : array (as show) or keyword (/HARD) for hardcopy or filename of ps files, emission' +; extinction albedo and size dist. +; Default = [''show.ps''] +; +; OPTIONAL INPUT PARAMETERS: +; None +; OUTPUTS: +; OPTIONAL OUTPUT PARAMETERS: +; SMDAT : from DUSTEM_GET_BAND_FLUX (see format there), structure containing model projected on data points: I_INST field in SMDAT. +; SMDAT also contains model direct outputs in M_EMIS, M_EXT fields (see format in ADD_MOD.pro ) +; ACCEPTED KEY-WORDS: +; help = If set, print this help +; restart = If set, reinitialise DustEMWrap and use a new +; wavelength vector for integration. +; COMMON BLOCKS: +; None +; SIDE EFFECTS: +; None +; RESTRICTIONS: +; The DustEM fortran code must be installed +; The DustEMWrap IDL code must be installed +; PROCEDURES AND SUBROUTINES USED: +; *** COMMENT AH --> is this really NONE? **** +; EXAMPLES + ;; print,' Examples:' + ;; print,' path=''/Users/lverstra/DUSTEM_LOCAL/dustem4.0/''' + ;; print,' to show SED only: show_dustem,path' + ;; print,' SED+UV-extinction: show_dustem,path,show=''extuv'' ' + ;; print,' SED+IR extinction+size dist.: show_dustem,path,show=[''extir'',''sdist''] ' + ;; print,' ps files with default names: show_dustem,path,show=[''extir'',''alb''],/hard ' + ;; print,' ps files with given names: show_dustem,path,show=[''extir'',''sdist''],hard=[''f1.ps'',''f2.ps''] ' + ;; print,' no plot : show_dustem,path,show=0' + ;; print,' to get band fluxes in structure SM: show_dustem,path,smdat=sm' + ;; print,' see structure SM (overall): help,/str,sm' + ;; print,' see structure SM model field: help,/str,sm.m_emis' + ;; print,' see structure SM inst field: help,/str,sm.i_dirbe' + ;; print,' see structure SM inst flux field: help,/str,sm.i_dirbe.flx' + ;; print,' array(i,j+1) i:index of inst band, j:index of grain type (ntype+1 is total SED) ' + ;; print,' to select or add instrument field(s): show_dustem,path,smdat=sm,inst=[''spire''] (on existing SM to add instrumentx) ' +; MODIFICATION HISTORY: +; Written by L Verstraete and M Compiegne, Spring 2010 +; Modified by LV : change to cgs, add band fluxes in SMDAT, 2011 +; Modified by LV & V Guillet : add polarization part, 2017 +; Further 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_show_fortran' + goto,the_end +END + +; specify the grain model +IF keyword_set(model) THEN BEGIN + use_model=strupcase(model) +ENDIF ELSE BEGIN + use_model='MC10' ;Default is the Compiegne et al 2010 model +ENDELSE + +; run the fortran via the wrapper once +dustem_init,model=use_model +st_model=dustem_read_all(!dustem_soft_dir) +dustem_write_all,st_model,!dustem_dat +st=dustem_run() + + +data_path = !dustem_wrap_soft_dir+'/Data/' +fortran_path = !dustem_dat +if keyword_set(dustem_data_dir) then fortran_path = dustem_data_dir +if keyword_set(data_dir) then data_path = data_dir + + if n_elements(NH) EQ 0 then nh = 1.d20 + if n_elements(INST) EQ 0 then begin + inst = ['DIRBE','FIRAS','HFI','WMAP'] + endif else begin + inst = strupcase(inst) + inst = inst(UNIQ(inst, SORT(inst))) + endelse + n_inst = n_elements(inst) + if n_elements(WREF) EQ 0 then wref = 2.5d2 + if n_elements(BBPAR) EQ 0 then BBPAR = [17.75d, 1.8d0] + if n_elements(SW) EQ 0 then sw = 0 + if n_elements(TIT) EQ 0 then tit = 'DustEM' + if n_elements(COM) EQ 0 then com = tit + if n_elements(WN) EQ 0 then wn = [0, 2, 1, 3, 6, 4, 5, 7, 8, 9, 10, 11, 12, 13] ; absolute display + if n_elements( XR ) NE 2 then begin + if keyword_set( COSMO ) then xr = [1., 1.e5] else xr = [ 1, 1.e3 ] + endif else xr = [min(xr),max(xr)] + if n_elements( YR ) NE 2 then begin + if keyword_set( COSMO ) then yr = [1.e-12, 1.e-3] else yr = [1e-8, 1e-4] + endif else yr = [min(yr),max(yr)] + if n_elements( DATA ) EQ 0 then data = 0 + if n_elements( SHOW ) EQ 0 THEN BEGIN + show = ['emis'] + ENDIF ELSE BEGIN + show = [STRLOWCASE(STRTRIM(STRING(show),2))] + IF SHOW(0) NE '0' THEN BEGIN + show = ['emis',show] ; always show the SED + show = show(UNIQ(show)) + ENDIF + ENDELSE + pgrid = [ 'emis', 'extuv', 'extir', 'alb', 'sdist','polext', 'polsed', 'align' ] + nplots = n_elements( SHOW ) ; nr of plots + nhard = n_elements( HARD ) + if nhard EQ 0 OR show(0) EQ '0' THEN hard=['0'] + if nhard GT 0 then hard = [STRLOWCASE(STRTRIM(STRING(hard),2))] + if ((nhard GT 0) AND (nhard LT nplots)) OR (hard(0) EQ '1') then begin + ig = INTARR( nplots ) + for i = 1, nplots-1 do begin + ig(i) = WHERE( pgrid EQ show(i)) + endfor + hard = pgrid(ig(sort(ig))) + REPLICATE('.ps',nplots) + endif + hard = strtrim(hard,2) + nhard = n_elements( HARD ) + + +; constants + na = 6.02d23 + clight = 2.9979246d10 ; cm/s + +; +; plot inits +; +; set the rgb colors + red =[0,1,1,0,0,1] + green=[0,1,0,1,0,1] + blue =[0,1,0,0,1,0] + tvlct,255*red, 255*green, 255*blue + ls = [ [ 0,3,1,4,5 ], [ 0,3,1,4,5 ] ] + ihard = 0 + +; +; get the SED data +; +; New FarIR-mm SED from FBoulanger + RESTORE, data_path+'/EXAMPLE_OBSDATA/filters_ref_wave.xdr' + fact_em = 0.77 ; to account for 20% of H to be in ionized form + 3% H2 + to_sed_20 = 1d-17 * (1d20/1.82d18) ; MJy/sr->erg/s/cm2/sr and normalization to NH = 10^20 H/cm2 + +; FIRAS + READCOL, data_path+'/EXAMPLE_OBSDATA/diffuse_ISM_SED.dat', wfirasf, firasf, ufirasf, skipline=5, numline=156, /sil + nufirasf = clight/(wfirasf*1.d-4) + firasf = firasf * nufirasf * fact_em * to_sed_20 + ufirasf = ufirasf * nufirasf * to_sed_20 ; error + +; DIRBE 60 --> 240 microns + READCOL, data_path+'/EXAMPLE_OBSDATA/diffuse_ISM_SED.dat', wdirbef, dirbef, udirbef, skipline=162, numline=4, /sil + nudirbef = clight/ (wdirbef*1.d-4) + dirbef = dirbef * nudirbef * fact_em * to_sed_20 + udirbef = udirbef * nudirbef * to_sed_20 ; error + dwdirbef = dwdirbe[6:9]/2. + +; WMAP + READCOL,data_path+'/EXAMPLE_OBSDATA/diffuse_ISM_SED.dat', wwmapf, wmapf, uwmapf, skipline=167, numline=5,/sil + wwmapf = wwmapf[3:4] + wmapf = wmapf[3:4] + uwmapf = uwmapf[3:4] + nuwmapf = clight/(wwmapf*1.d-4) + wmapf = wmapf * nuwmapf * fact_em * to_sed_20 + uwmapf = uwmapf * nuwmapf * to_sed_20 ; error + +; DIRBE Arendt et al (1998) for |b|>25 +; numbers from Li & Draine, apj, 2001, 554, 778-802 + Wave_ARENDT = WDIRBE[2:*] + Dwave_ARENDT = DWDIRBE[2:*] / 2. + ARENDT = [0.97, 1.11, 7.16, 3.57, 5.30, 18.6, 22.5, 10.1] ; 10^-26 erg/s/H/sr + ARENDT = ARENDT * 1.d-26 * 1.d20 ; To erg s-1 cm-2 sr-1 for NH=10^20 H/cm2 + err_ARENDT = 0.2 * ARENDT ; from Dwek et al. 1997 DIRBE 1Sigma unc = 20% + nu_arendt = nudirbe[2:*] * 1.e9 + +; Normalization of the Arendt spectrum on the Boulanger 100 microns + wave_arendt_midIR = wave_arendt[0:3] + dwave_ARENDT_midir = Dwave_ARENDT[0:3] + nu_arendt_midIR = nu_arendt[0:3] + correl_coeff_midIR = [0.00183, 0.00291, 0.0462, 0.0480] ; for Inu + ucorrel_coeff_midIR = [0.00001, 0.00003, 0.0001, 0.0002] ; For Inu + arendt_midIR = correl_coeff_midIR * (dirbef[1]*100.*1e-4/clight*1e20) ; in (Inu) MJy sr-1 + err_arendt_midIR = FLTARR(N_ELEMENTS(arendt_midIR)) + for i=0,3 do begin + tmp= SQRT( ((udirbef[1]*100.*1e-4/clight*1e20)/(dirbef[1]*100.*1e-4/clight*1e20))^2. + (ucorrel_coeff_midIR[i]/correl_coeff_midIR[i])^2. ) + err_arendt_midIR[i] = arendt_midIR[i] * tmp + endfor + arendt_midIR = arendt_midIR / wave_arendt_midIR/1e-4*clight/1e20 ; to erg s-1 cm-2 sr-1 + err_arendt_midIR = err_arendt_midIR / wave_arendt_midIR/1e-4*clight/1e20 ; to erg s-1 cm-2 sr-1 + wdirbe_ful = [ wave_arendt_midir, wdirbef ] + dwdirbe_ful = [ dwave_arendt_midir, dwdirbef ] + dirbe_ful = [ arendt_midir, dirbef ] + udirbe_ful = [ err_arendt_midir, udirbef ] + +; isocam galactic spectrum + RESTORE, data_path+'/EXAMPLE_OBSDATA/spectre_gal.xdr' ; wgal in microns, spec_gal in MJy/sr + nuisocam = clight / wgal/1.d-4 + norm = 0.65*0.045/filters(9) ; normalization for Nh=1e20 cm-2 and I12/I100 =0.045 + iso = 1d-17 * spec_gal*nuisocam * norm ; galactic mid-IR SED in erg/s/cm2/sr + stars = 1d3 * wgal * DUSTEM_BLACKBODY(wgal, 3d3, unit='w') ; assumed spectrum for stars in cgs + stars = 1.3d-6 * stars/stars(0) ; normalization of stellar spectrum + err_cvf = 0.2 + isocam = iso-stars + smp = DUSTEM_GET_BAND_FLUX( data_path + 'FILTERS/', 'DIRBE', xs=wgal, ys=isocam ) + isocam = isocam / smp.i_dirbe.ym(4) * arendt_midIR[2] ; factor is now 1.013027 (1.0235712 before) + +; Arome 3.3 microns: Giard et al 1988 + xg = [3.3d] + yg = 1.1e-6 / 4./!pi * xg/0.05 ; SED erg/s/cm2/sr assumes a feature width of 0.05 mic + eg = 0.5e-6 / 4./!pi * xg/0.05 ; error + +; Modified BB overlaid to dust: default is 17.75 K and beta=1.8 (Miville-DeschĂȘnes et al 2013, Planck data) + np = 1000 + wbr = ALOG10( [20.,1.d5] ) + dwbr = (wbr(1)-wbr(0)) / (np-1) + lambdaf = 10.^(wbr(0) + findgen(np)*dwbr) + temp1 = bbpar(0) + beta = bbpar(1) + lamb_ref = 250. + qf = (lambdaf/ lamb_ref)^(-beta) + bbm = 1d3 * lambdaf * DUSTEM_BLACKBODY(lambdaf, temp1, unit='w') + SP1 = bbm*qf + firas_nuinu = firasf +; normalize to FIRAS @ wave WREF < wfiras_max and get dust cross-section per H + wfiras_max = 8d2 + if WREF LE wfiras_max then lamb_ref = wref + tt = MIN( ABS(lambdaf-lamb_ref), ibb ) + tt = MIN( ABS(wfirasf-lamb_ref), ifiras ) + nw = 1 + eps_firas = median(firasf(ifiras-nw:ifiras+nw)) / INTERPOL(sp1/qf, lambdaf, lamb_ref) / nh + SP1 = SP1 * median(firasf(ifiras-nw:ifiras+nw)) / INTERPOL(sp1, lambdaf, lamb_ref) + +; add HFI (Planck2011t, Dust in the diffuse interstellar medium and the Galactic halo) + nu_pep_dism = [5d3,3d3,857.,545.,353.] * 1d9 ; includes IRAS60 and 100 + w_pep_dism = 1d4*clight / nu_pep_dism + dw_pep_dism = dblarr(n_elements(nu_pep_dism)) + dw_pep_dism(0:1) = [31.,35.6] / 2. + dw_pep_dism(2:*) = w_pep_dism(2:*)/3. / 2. + y_pep_dism = [ 1.288, 6.522,5.624, 1.905, 0.465 ] * 1d-1 ; MJy/sr/1e20cm-2 + u_pep_dism = [ 0.436, 1.473, 1.015, 0.347, 0.100 ] * 1d-1 + y_pep_dism = fact_em * y_pep_dism * nu_pep_dism * 1d-17 ; MJy/sr/1e20cm-2 to erg/s/cm2/sr + u_pep_dism = fact_em * u_pep_dism * nu_pep_dism * 1d-17 + y_hfi_dism = y_pep_dism(0:2) & u_hfi_dism = u_pep_dism(0:2) + +; add cosmosomas + IF keyword_set( COSMO ) then begin +; Watson et al 2005 +; nu=1.e9*[ 0.408, 1.42, 10.9, 12.7, 14.7, 16.3, 22.8, 33.0, 40.9, 61.3, 93.8, 1250, 2143., 3.e3 ] +; i_jy=[6.3,7.3,17.,18.7,26.2,32.6,42.3,40.3,33.9,34.7,77.5,9.68e4,1.31e5,6.44e4 ] +; e_jy=[7.8,2.,0.1,0.4,0.6,1.5,0.2,0.4,0.7,1.8,4.3,7d2,1.25d3,100.] +; beam = 4.9e-4 ; steradians +; lref = 140d ; wave for normalization + +; nu in GHz and i_jy in Janskys +; Planck 2011p, "New light on anomalous microwave emission from spinning dust grains" + nu=[0.408,0.82,1.42,10.9,12.7,14.7,16.3,22.8,28.5,33.0,40.9,44.1,61.3,70.3,93.8,100.,143.,217.,353.,545.,857.,1250,2143.,3.e3] + iplck = where( (nu eq 33.) or (nu eq 40.9) or (nu eq 70.3) or (nu eq 100.) or (nu eq 143.) or (nu eq 217.) ) + nu=1d9*[0.408,0.82,1.42,10.9,12.7,14.7,16.3,22.8,28.5,33.0,40.9,44.1,61.3,70.3,93.8,100.,143.,217.,353.,545.,857.,1250,2143.,3.e3] + i_jy=[9.7,9.4,8.0,16.1,20.,28.4,35.8,39.8,38.1,36.7,30.8,28.5,27.4,32.2,63.,78.4,202.,1050.,3060.,1.53d4,4.87d4,9.3d4,1.17d5,5.36d4] + e_jy=[3.3,4.7,1.7,1.8,2.2,3.1,4.,4.2,4.,3.9,3.3,3.2,3.4,3.9,7.8,15.4,22.,128.,467.,2.09d3,6.11d3,1.29d4,1.45d4,6.67d3] + fwhm = 1.81d + beam = !pi*(!pi/1.8d2 * fwhm/2d0)^2 ; steradians + lref = 140d ; wave for normalization + + x_pep_ame = 1d4*clight/nu + nc = n_elements(x_pep_ame) + y_pep_ame = 1.e-23 * nu * i_jy / beam + e_pep_ame = 1.e-23 * nu * e_jy / beam + lref = 240d + iref = WHERE( wdirbe_ful EQ lref, cr) + if cr EQ 1 then begin + yic = INTERPOL(y_pep_ame,x_pep_ame,wdirbe_ful[iref]) + y_norm = dirbe_ful[iref]/yic + y_norm = y_norm[0] + endif else y_norm = 1d + is = WHERE( x_pep_ame/lref GE 0.9, cs) + x_pep_ame=x_pep_ame(is) & y_pep_ame=y_pep_ame(is)*y_norm & e_pep_ame=e_pep_ame(is)*y_norm + ENDIF + +; +; plot SED data +; + IF SHOW(0) NE '0' THEN BEGIN + if HARD(0) NE '0' then begin + !y.thick = 5 + !x.thick = 5 + !p.thick = 5 + !p.charsize = 1.3 + !p.charthick = 5 + !x.ticklen = 0.04 + set_plot,'ps' +; device,file=hard(0),xs=26, ys=20,/portrait,/color + device,file=hard(0),/portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(0), xs=900,ys=600,tit='SHOW_DUSTEM: SED '+strtrim(wn(0),2) + endelse + + xtit = 'Wavelength (!4l!3m)' + ytit = '!4m!3 I!d!4m!3!n (erg s!u-1!n cm!u-2!n sr!u-1!n)' + fine = 1 +; plot_oo, INDGEN(1),/NODAT,/xs,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit,tit=tit + plot_oo, INDGEN(1),/NODAT,xs=9,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit + axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)' + if DATA NE 0 then begin + if HARD(0) NE '0' then begin + oploterror,xg,yg,0.,eg,psym=5 + oploterror, wfirasf, firasf, wfirasf*0d, ufirasf, errthick=0.7, /nohat + oploterror, wdirbe_ful(0:3), dirbe_ful(0:3), dwdirbe_ful(0:3), udirbe_ful(0:3), psym=6 + oploterror, wdirbe_ful(6:*), dirbe_ful(6:*), dwdirbe_ful(6:*), udirbe_ful(6:*), psym=6 +; oploterror, wdirbe_ful, dirbe_ful, dwdirbe_ful, udirbe_ful, psym=6 + oploterror, wgal,isocam, wgal*0d, isocam*err_cvf, errthick=0.7, /nohat + oploterror, wwmapf, wmapf, wwmapf*0d, uwmapf, psym=4 + oploterror, w_pep_dism, y_pep_dism, dw_pep_dism, u_pep_dism, psym=5 + endif else begin + oploterror,xg,yg,0.,eg,psym=5 + oploterror, wfirasf, firasf, wfirasf*0d, ufirasf, errthick=0.4, /nohat + oploterror, wdirbe_ful(0:3), dirbe_ful(0:3), dwdirbe_ful(0:3), udirbe_ful(0:3), psym=6 + oploterror, wdirbe_ful(6:*), dirbe_ful(6:*), dwdirbe_ful(6:*), udirbe_ful(6:*), psym=6 +; oploterror, wdirbe_ful, dirbe_ful, dwdirbe_ful, udirbe_ful, psym=6 + oploterror, wgal,isocam, wgal*0d, isocam*err_cvf, errthick=0.4, /nohat + oploterror, wwmapf, wmapf, wwmapf*0d, uwmapf, psym=4 + oploterror, w_pep_dism, y_pep_dism, dw_pep_dism, u_pep_dism, psym=5 + endelse + endif + if DATA EQ 1 then begin + OPLOT, LAMBDAf, SP1, lines=1 ; FIR blackbody + sp1_firas = INTERPOL( sp1, lambdaf, wfirasf ) + firas_chi2 = DUSTEM_CHI2( firasf, sp1_firas, 3, err=ufirasf ) + print,'Grey body: Td = ',strtrim(bbpar(0),2),' and beta = ',strtrim(bbpar(1),2) + print,'Chi-squared of grey body to FIRAS: ', strtrim(firas_chi2,2) + endif + if keyword_set( CMB ) then begin + lambda_cmb = 1.d2^( 1.d0 + dindgen(500)*0.01 ) + sp_cmb = 1d3 * lambda_cmb * DUSTEM_BLACKBODY(LAMBDA_CMB, 2.728, unit='w') ; SED in erg/s/cm2/s/sr + OPLOT, LAMBDA_CMB, SP_CMB, lines=0 ; CMB + endif + if keyword_set( COSMO ) AND DATA NE 0 then begin + oploterror, x_pep_ame, y_pep_ame, y_pep_ame*0d, e_pep_ame, ps=4 +; oploterror, x_pep_ame(iplck), y_pep_ame(iplck), y_pep_ame(iplck)*0d, e_pep_ame(iplck), ps=5, syms=1.5, col=3 + endif + ENDIF + +; +; get the model SED +; + fname = fortran_path + 'data/GRAIN.DAT' + nlines = FILE_LINES( fname ) + OPENR, uu, fname, /get_lun + tmp = '#' + print,'(W) DUSTEM_SHOW_FORTRAN: GRAIN.DAT' + cnt = 0 + WHILE STRPOS(tmp,'#') EQ 0 do begin + READF, uu, tmp + cnt = cnt + 1 + ENDWHILE + r_opt = STRLOWCASE(STRTRIM(STRING(tmp),2)) + print,r_opt + READF, uu, g0 + ntype = nlines - (cnt+1) + nsize=intarr(ntype) + t_opt=strarr(ntype) + gtype = strarr(ntype) + propm = dblarr(ntype) + rho = dblarr(ntype) + for i=0,ntype-1 do begin + READF,uu,tmp + PRINT, tmp + tt = strsplit(tmp, /extract) + gtype(i) = strtrim(tt(0)) + nsize(i)=fix(tt(1)) + t_opt(i)=strlowcase(strtrim(tt(2))) + propm(i) = double(tt(3)) + rho(i) = double(tt(4)) + endfor + close,uu + free_lun,uu + nsz_max = MAX(nsize) + + IF n_elements(FSED) EQ 0 THEN fsed = fortran_path+'out/SED.RES' + OPENR, uu, fsed, /get_lun + tmp = '#' + WHILE (STRPOS(tmp,'#') EQ 0) do begin + READF, uu, tmp + ENDWHILE + tt = double( strsplit(tmp, ' ', /extract) ) + ntype_sed = fix(tt(0)) + if ntype_sed NE ntype then begin + print,'(F) DUSTEM_SHOW_FORTRAN: SED.RES & GRAIN.DAT have different NTYPE' + print,' data is not from present GRAIN.DAT' + return + endif + nlamb = fix(tt(1)) ; nr of wavelengths + x = dblarr(nlamb) + sedh = dblarr(nlamb,ntype+1) + for i=0,nlamb-1 do begin + READF, uu, tmp + tt = double( strsplit(tmp, ' ', /extract) ) + x(i) = tt(0) + sedh(i,*) = tt(1:*) + endfor + close,uu + free_lun,uu + xlamb= x + sed = sedh * nh / 4. / !pi ; in erg/s/cm2/sr + if keyword_set( SW ) then sed = SMOOTH( sed, sw) + + if (DATA NE 1) AND (SHOW(0) NE '0') then begin + dy = 10.^((ALOG10(yr(1))-ALOG10(yr(0))) / 25.) + dx = 10.^((ALOG10(xr(1))-ALOG10(xr(0))) / 50.) + yps= 10.^((ALOG10(yr(1))+ALOG10(yr(0))) / 2. + (ALOG10(yr(1))-ALOG10(yr(0))) / 2.3) + if keyword_set(COSMO) then begin + xpr=10.^((ALOG10(xr(1))+ALOG10(xr(0)))/3.) + yps= 10.^((ALOG10(yr(1))+ALOG10(yr(0))) / 2. - (ALOG10(yr(1))-ALOG10(yr(0))) / 4.) + endif else xpr=xr(0) + xpr = xpr*[dx,dx^4] +; overlay model SED in erg/s/cm2/sr + for i=0,ntype-1 do begin + oplot, x, sed(*,i), lin=ls(i+1) + oplot,xpr,yps*[1.,1.], lin=ls(i+1) + xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1 + yps = yps/dy + endfor + oplot, x, sed(*,ntype),lin=ls(0),col=2 + s1 = STRMID(STRTRIM(ROUND(1e1*g0)/1e1,2), 0, 4) + s2 = STRMID(STRTRIM(ROUND(1e2*alog10(nh))/1e2,2), 0, 5) + xyouts,/norm,0.6,0.85,'G!d0!n='+s1+' N!dH!n=10!u'+s2+'!ncm!u-2!n' +; overlay model polarized SED in erg/s/cm2/sr + ;IF (C_POLSED EQ 1) THEN oplot, x, sed_p(*,ntype),lin=ls(0),col=3 + endif + +; +; fill in SMDAT +; +; first get model fluxes in instrument bands + unit = 'x(microns) SED(erg/s/cm2/sr)' + if n_elements(SMDAT) GT 0 then begin + smdat = DUSTEM_GET_BAND_FLUX( data_path + 'FILTERS/', inst, xs=x, ys=sed, unit=unit, smi=smdat ) + endif else smdat = DUSTEM_GET_BAND_FLUX( data_path + 'FILTERS/', inst, xs=x, ys=sed, unit=unit ) + smdat.com = com + stag = TAG_NAMES(smdat) + ;if c_polsed EQ 1 then smdat.m_emis.yp = sed_p + +; then get diffuse data available here + ii = WHERE( inst EQ 'DIRBE', cic) + if cic GT 0 then begin + smdat.i_dirbe.unit = smdat.i_dirbe.unit + ' YD(erg/s/cm2/sr)' + smdat.i_dirbe.yd = [ 0d, 0d, dirbe_ful ] + smdat.i_dirbe.err = [ 0d, 0d, udirbe_ful] + nband = n_elements( smdat.i_dirbe.x) + smdat.i_dirbe.isel = [ 0,0,intarr(nband-2)+1 ] + smdat.i_dirbe.npar = ntype + endif + + ii = WHERE( inst EQ 'HFI', cic) + if cic GT 0 then begin + smdat.i_hfi.unit = smdat.i_hfi.unit + ' YD(erg/s/cm2/sr)' + ibx = indgen(3) ; only 3 first bands in PEP DISM + smdat.i_hfi.yd(ibx) = y_hfi_dism + smdat.i_hfi.err(ibx) = u_hfi_dism + smdat.i_hfi.isel = [1,1,1,0,0,0] + smdat.i_hfi.npar = 1 + endif + + ii = WHERE( inst EQ 'WMAP', cic) + if cic GT 0 then begin + smdat.i_wmap.unit = smdat.i_wmap.unit + ' YD(erg/s/cm2/sr)' + smdat.i_wmap.yd = [ 0d, 0d, 0d, wmapf ] ; only 61 and 94 GHz bands + smdat.i_wmap.err = [ 0d, 0d, 0d, uwmapf] + nband = n_elements( smdat.i_wmap.x) + smdat.i_wmap.isel = [ 0,0,0,intarr(nband-3)+1 ] + smdat.i_wmap.npar = 1 + endif + +; put FIRAS pts in SMDAT + ii = WHERE( inst EQ 'FIRAS', cic) + if cic GT 0 then begin + nband = n_elements(wfirasf) + smdat = DUSTEM_ADD_INST( smdat, 'FIRAS', [n_elements(firasf),ntype+1] ) + smdat.i_firas.unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr) YD(erg/s/cm2/sr)' + smdat.i_firas.yd = firasf + smdat.i_firas.err = ufirasf + sd_firas = dblarr(n_elements(wfirasf),ntype+1) + for itp=0,ntype do sd_firas(*,itp) = INTERPOL( sed(*,itp), x, wfirasf ) + smdat.i_firas.ym = sd_firas + smdat.i_firas.isel = intarr(nband)+1 + smdat.i_firas.npar = 2 + endif + +; overlay model band fluxes + if (DATA NE 1) AND (SHOW(0) NE '0') then begin + itg = WHERE( STRPOS(stag,'I_') GE 0 AND stag NE 'I_FIRAS', ctg ) +; then other instruments + if ctg GT 0 then begin + for k=0,ctg-1 do begin + oplot, smdat.(itg(k)).x, smdat.(itg(k)).ym(*,ntype), ps=6, syms=1.5, col=3 + endfor + endif + endif + +; +; get extinction +; + OPENR, uu, fortran_path+'out/EXT.RES', /get_lun + tmp = '#' + WHILE STRPOS(tmp,'#') EQ 0 do begin + READF, uu, tmp + ENDWHILE + tt = double( strsplit(tmp, ' ', /extract) ) + ntype_ext = fix(tt(0)) + if ntype_ext NE ntype then begin + print,'(F) DUSTEM_SHOW_FORTRAN: EXT.RES & GRAIN.DAT have different NTYPE' + print,' data is not from present GRAIN.DAT' + return + endif + nlamb = fix(tt(1)) ; nr of wavelengths + x = dblarr(nlamb) + stmp = dblarr(nlamb,2*ntype+1) + ssca = dblarr(nlamb,ntype+1) + sabs = dblarr(nlamb,ntype+1) + sext = dblarr(nlamb,ntype+1) + alb = dblarr(nlamb,ntype+1) + for i=0,nlamb-1 do begin + READF, uu, tmp + tt = double( strsplit(tmp, ' ', /extract) ) + x(i) = tt(0) + stmp(i,*) = tt(1:*) + endfor + CLOSE,uu + FREE_LUN,uu + wlm = x + x = 1. / wlm ; 1/micron + + mH = 1.67262158d-24 ; proton mass /gdust -> /gH + pm = [propm,propm] + for i=0,ntype-1 do begin + sabs(*,i) = stmp(*,i) + ssca(*,i) = stmp(*,ntype+i) + sext(*,i) = sabs(*,i) + ssca(*,i) + alb(*,i) = ssca(*,i) / sext(*,i) + endfor + for i=0,nlamb-1 do begin + sabs(i,ntype) = TOTAL(sabs(i,0:ntype-1)) + ssca(i,ntype) = TOTAL(ssca(i,0:ntype-1)) + sext(i,ntype) = TOTAL(sext(i,0:ntype-1)) + alb(i,ntype) = ssca(i,ntype)/sext(i,ntype) + endfor + +; get RV + av = INTERPOL( sext(*,ntype), wlm, 0.55 ) + ab = INTERPOL( sext(*,ntype), wlm, 0.44 ) + rv = av / (ab-av) + +; get 250 microns emissivity + eps_a = dblarr(ntype+1) & eps_e = eps_a & eps_s = eps_a + for i = 0,ntype-1 do begin + eps_a(i)=INTERPOL(sabs(*,i), wlm, wref) + eps_s(i)=INTERPOL(ssca(*,i), wlm, wref) + eps_e(i)=INTERPOL(sext(*,i), wlm, wref) + endfor + eps_a(ntype) = INTERPOL(sabs(*,ntype), wlm, wref) + eps_s(ntype) = INTERPOL(ssca(*,ntype), wlm, wref) + eps_e(ntype) = INTERPOL(sext(*,ntype), wlm, wref) + print,'(W) DUSTEM_SHOW_FORTRAN: model dust cross-section @ ',STRTRIM(wref,2),' per type and total (cm2/H)', $ + format='(A44,1E10.3,A27)' + print, ' Abs : ',eps_a, format='(A8,100(1E12.4))' + print, ' Sca : ',eps_s, format='(A8,100(1E12.4))' + print, ' Ext : ',eps_e, format='(A8,100(1E12.4))' + if WREF GT wfiras_max then begin + print,'(W) DUSTEM_SHOW_FORTRAN: WREF above longer FIRAS wave, using 250 microns' + wp = lamb_ref + endif else wp = wref + print,'(W) DUSTEM_SHOW_FORTRAN: dust cross-section @ ', STRTRIM(wp,2), ' microns from FIRAS',format='(A38,1E10.3,A20)' + print,' using T = ',bbpar(0), ' and beta = ',bbpar(1),' sigma(cm2/H) = ', eps_firas, $ + format='(A27,1E8.2,A13,1E8.2,A17,1E12.4)' + +; fill in model + it = WHERE( stag EQ 'M_EXT', ct) + if ct EQ 0 then begin + smdat = DUSTEM_ADD_MOD( smdat, 'EXT', [n_elements(wlm),ntype+1]) + smdat.m_ext.x = wlm + smdat.m_ext.y = sext + smdat.m_ext.abs = sabs + smdat.m_ext.sca = ssca + smdat.m_ext.alb = alb + smdat.m_ext.xr = wref + smdat.m_ext.yr_abs = eps_a + smdat.m_ext.yr_sca = eps_s + smdat.m_ext.rv = rv + endif + +; get standard Savage & Mathis 79 + Mathis 1990 extinction and fill in data + READCOL, data_path+'/EXAMPLE_OBSDATA/mean_ism_ext.dat', skip=3, x_sm, e_sm, /SIL + it = WHERE( stag EQ 'I_SM79', ct) + if ct EQ 0 then begin + smdat = DUSTEM_ADD_INST( smdat, 'SM79', [n_elements(x_sm),ntype+1] ) + smdat.i_sm79.x = x_sm + smdat.i_sm79.ym = INTERPOL( sext(*,ntype), wlm, x_sm) + smdat.i_sm79.yd = e_sm + smdat.i_sm79.err = smdat.i_sm79.yd*0.2 + smdat.i_sm79.unit = 'x(microns) sigma(cm2/H)' + smdat.i_sm79.npar = ntype + smdat.i_sm79.isel = intarr(n_elements(x_sm)) + 1 + endif + +; +; plot the vis-UV extinction +; + yscl = 1d0 + ip = WHERE( STRPOS(show,'extuv') GE 0, cp ) + IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(1), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Extinction UV '+strtrim(wn(1),2) + endelse + xtit = 'x (!4l!3m!u-1!n)' + xr = [0,11] + if n_elements(WEXT) then begin + xe = [x,1./wext] + xe = xe(SORT(XE)) + ei = INTERPOL(sext(*,ntype),x,xe) + tmp = MIN( ABS(xe - 1./wext),iw ) + yr = minmax( sext(*,ntype)/ei(iw) ) * [1.,1.] + yscl = 1.d0 / ei(iw) + ytit='Normalized !4r!3!dext!n ' + print,'(W) DUSTEM_SHOW_FORTRAN: extinction normalized to 1 at ',wext, ' microns yscl=',yscl + endif else begin + yr=[ 0, 2.5 ] + yscl = 1d0 + ytit='!4r!3!dext!n (10!u-21!n cm!u2!n per H)' + endelse + plot, [1d], [1d], lin=ls(0), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit + oplot, x, sext(*,ntype)*yscl, lin=ls(0), col=2 + oplot, 1d/x_sm, e_sm*1d21*yscl, ps=4 +; oploterror, 1d/x_sm, e_sm*1d21, smdat.i_sm79.err*1d21, psym = 4, /nohat + for i=0,ntype-1 do oplot, x, (sext(*,i))*yscl, lin=ls(i+1) + xyouts,/norm,0.2,0.8,'R!dV!n='+STRMID(STRTRIM(ROUND(1e1*rv)/1e1,2), 0, 4) + ENDIF + +; plot IR extinction + ip = WHERE( STRPOS(show,'extir') GE 0, cp ) + IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(2), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Extinction IR '+strtrim(wn(2),2) + endelse + xtit = textoidl('\lambda (\mum)') + xr = [1,400] + yr = [1e-6,1.] + if yscl NE 1d0 THEN ytit='Normalized !4r!3!dext!n' ELSE $ + ytit='!4r!3!dext!n (10!u-21!n cm!u2!n per H)' + plot_oo, [1d], [1d], xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit + oplot, wlm, sext(*,ntype)*yscl, lin=ls(0), col=2 + oplot, x_sm, e_sm*1d21*yscl, ps=4 +; oploterror, x_sm, e_sm*1d21, smdat.i_sm79.err*1d21, psym = 4, /nohat + for i=0,ntype-1 do oplot, wlm, sext(*,i)*yscl, lin=ls(i+1) + xyouts,/norm,0.8,0.85,'R!dV!n='+STRMID(STRTRIM(ROUND(1e1*rv)/1e1,2), 0, 4) + ENDIF + +; +; and the albedo +; + ip = WHERE( STRPOS(show,'alb') GE 0, cp ) + IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(3), xs=500,ys=350, tit='DUSTEM_SHOW_FORTRAN: Albedo '+strtrim(wn(3),2) + endelse + xtit = 'x (!4l!3m!u-1!n)' + xr = [0,11] + ytit='Albedo' + yr=[ 0, 1] + plot, x, alb(*,ntype), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit + oplot, x, alb(*,ntype), col=2 + for i=0,ntype-1 do oplot, x, alb(*,i), lin=ls(i+1) + ENDIF + + +; +; get the size distribution and plot +; + ip = WHERE( STRPOS(show,'sdist') GE 0, cp ) + ir = WHERE( STRPOS(r_opt,'sdist') GE 0, cr ) + IF CR EQ 0 THEN BEGIN + cp = 0 + print,'(W) DUSTEM_SHOW_FORTRAN: SDIST keyword not set in current run' + ENDIF ELSE IF ((SHOW(0) NE '0') AND (CP EQ 1)) THEN BEGIN + fn = fortran_path+'out/SDIST.RES' + OPENR, uu, fn, /get_lun + tt = '#' + WHILE STRPOS(tt,'#') EQ 0 do begin + READF, uu, tt + ENDWHILE + ax = dblarr(nsz_max,ntype) + ava = dblarr(nsz_max,ntype) + nsz_tot = TOTAL(nsize) + ava_tot = dblarr(nsz_tot,ntype+1) + ax_tot = 0d + FOR i=0,ntype-1 do begin + READF,uu,tt + for is=0,nsize(i)-1 do begin + READF, uu, tx, ty + ax(is,i) = tx + ava(is,i) = ty + endfor + ax_tot = [ax_tot, ax(0:nsize(i)-1,i)] + ENDFOR + close,uu + free_lun,uu + ax_tot = ax_tot(1:*) + ax_tot = ax_tot(SORT(ax_tot)) + FOR i=0,ntype-1 do begin + tt = INTERPOL( ava(0:nsize(i)-1,i), ax(0:nsize(i)-1,i), ax_tot ) + ix = WHERE( ax_tot LT MIN(ax(*,i)) OR ax_tot GT MAX(ax(*,i)), cx ) + if cx GT 0 then tt(ix) = 0d ; no extrapolation + ava_tot(*,i) = tt + endfor + for i=0, nsz_tot-1 do ava_tot(i,ntype) = TOTAL( ava_tot(i,0:ntype-1)) + +; fill in model + it = WHERE( stag EQ 'M_SDIST', ct) + if ct EQ 0 then begin + smdat = DUSTEM_ADD_MOD( smdat, 'SDIST', [nsz_tot,ntype+1,nsz_max]) + smdat.m_sdist.xtot = ax_tot + smdat.m_sdist.ytot = ava_tot + smdat.m_sdist.xi = ax + smdat.m_sdist.yi = ava + endif + + if HARD(0) NE '0' then begin + device, file=hard(nplots-1), /portrait,/color + endif else begin + window, wn(5), xs=500,ys=350, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: Size Distribution '+strtrim(wn(5),2) + endelse + yscl = 1.D29 + xr=[0.1,1e4] + yr=[-3,3] + xtit='a (nm)' + ytit='log( 10!u29!n n!dH!u-1!n a!u4!n dn/da (cm!u3!n/H)) ' + tit='DustEM' + plot_oi, ax(0:nsize(0)-1,0)*1.e7, alog10(ava(0:nsize(0)-1,0)*yscl), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit='DUSTEM' + FOR i=1,ntype-1 DO begin + oplot, ax(0:nsize(i)-1,i)*1.e7, alog10(ava(0:nsize(i)-1,i)*yscl), line=ls(i+1) + ENDFOR + + ENDIF + + +; Display polarized SED +; + ip = WHERE( STRPOS(show,'polsed') GE 0, c_polsed ) + IF (SHOW(0) NE '0') AND (C_POLsed EQ 1) THEN BEGIN + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(6), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Polarized SED '+strtrim(wn(6),2) + endelse + + ; get polarized SED + OPENR, uu, fortran_path+'out/SED_POL.RES', /get_lun + tmp = '#' + WHILE (STRPOS(tmp,'#') EQ 0) do begin + READF, uu, tmp + ENDWHILE + tt = double( strsplit(tmp, ' ', /extract) ) + sedh_p = dblarr(nlamb,ntype+1) + for i=0,nlamb-1 do begin + readf, uu, tmp + tt = double( strsplit(tmp, ' ', /extract) ) + sedh_p(i,*) = tt(1:*) + x(i) = tt(0) + endfor + close,uu + free_lun,uu + sed_p = sedh_p * nh / 4. / !pi ; in erg/s/cm2/sr + + yr = [1e-10,1e-5] + xr = [10,1e4] + dy = 10.^((ALOG10(yr(1))-ALOG10(yr(0))) / 25.) + dx = 10.^((ALOG10(xr(1))-ALOG10(xr(0))) / 50.) + yps= 10.^((ALOG10(yr(1))+ALOG10(yr(0))) / 2. + (ALOG10(yr(1))-ALOG10(yr(0))) / 2.3) + xpr=xr(0) + xpr = xpr*[dx,dx^4] + + xtit = 'Wavelength (!4l!3m)' + ytit = '!4m!3 P!d!4m!3!n (erg s!u-1!n cm!u-2!n sr!u-1!n)' + fine = 1 + plot_oo, INDGEN(1),/NODAT,xs=9,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit + axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)' + for i=0,ntype-1 do begin + oplot, x, sed_p(*,i), lin=ls(i+1) + oplot,xpr,yps*[1.,1.], lin=ls(i+1) + xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1 + yps = yps/dy + endfor + oplot,x,sed_p(*,ntype),lin=ls(0),col=3 + ENDIF + + +; +; get the polarisation extinction +; + ip = WHERE( STRPOS(show,'polext') GE 0, c_polext ) + IF (SHOW(0) NE '0') AND (C_POLEXT EQ 1) THEN BEGIN + + itp = WHERE( STRPOS(t_opt,'pol') GE 0, ctp) + IF ctp EQ 0 THEN BEGIN + c_pol = 0 + print,'(W) DUSTEM_SHOW_FORTRAN: no POL data in current run' + ENDIF ELSE BEGIN + OPENR, uu, fortran_path+'out/EXT_POL.RES', /get_lun + tmp = '#' + WHILE STRPOS(tmp,'#') EQ 0 do begin + READF, uu, tmp + ENDWHILE + tt = double( strsplit(tmp, ' ', /extract) ) + ntype_pol = fix(tt(0)) + if ntype_pol NE ntype then begin + print,'(F) DUSTEM_SHOW_FORTRAN: POL.RES & GRAIN.DAT have different NTYPE' + print,' data is not from present GRAIN.DAT' + return + endif + nlamb = fix(tt(1)) ; nr of wavelengths + x = dblarr(nlamb) + spabs = dblarr(nlamb,ntype+1) + spsca = dblarr(nlamb,ntype+1) + for i=0,nlamb-1 do begin + READF, uu, tmp + tt = double( strsplit(tmp, ' ', /extract) ) + x(i) = tt(0) + spabs(i,0:ntype-1) = tt(1:ntype) + spsca(i,0:ntype-1) = tt(ntype+1:2*ntype) + spabs(i,ntype) = TOTAL(spabs(i,0:ntype-1)) + spsca(i,ntype) = TOTAL(spsca(i,0:ntype-1)) + endfor + close,uu + free_lun,uu + +; fill in model + smdat.m_ext.abs_p = spabs + smdat.m_ext.sca_p = spsca + + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(7), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: Serkowski '+strtrim(wn(7),2) + endelse + +; plot Serkowski + yscl = 1.D23 + yscl2 = 1.D2 + xr=[0.2,10.] + yr=[0.1,5] + xtit = textoidl('1/\lambda (\mum^{-1})') + ytit = textoidl('\sigma_{pol} (10^{-23}cm^2/H)') + tit='DUSTEM' + plot_oo, 1./x, yscl2*(spabs(*,0)+spsca(*,0)), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit + FOR i=1,ntype-1 DO oplot, 1./x, yscl2*(spabs(*,i)+spsca(*,i)), line=ls(i+1) + oplot, 1./x, yscl2*(spabs(*,ntype)+spsca(*,ntype)), col=3 + ix = WHERE(x_sm LE 7d, csm) ; keep wave below 7 microns + if csm GT 0 then begin + xx = 1./x_sm(ix) + xe = e_sm(ix) + endif + ix = SORT(xx) + xx = xx(ix) & xe = xe(ix) + fpol = DUSTEM_SERKOWSKI(xx) + oplot,xx,yscl*fpol/5.8e21*3.1,ps=4 +; oplot,xx,25*fpol*3.1/5.8e21/xe, ps=5 + ENDELSE + + ENDIF + +; +; Polarization fraction in emission +; + IF (SHOW(0) NE '0') AND (C_POLSED EQ 1) THEN BEGIN +; plot fractional polarisation + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(8), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: P/I '+strtrim(wn(8),2) + endelse + xr=[10,1.e4] + yr=[0.,0.4] + xtit = textoidl('\lambda (\mum)') + ytit = textoidl('P/I') + tit='' + plot_oi, x, sed_p(*,0)/sed(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit + axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)' + FOR i=1,ntype-1 DO begin + oplot, x, sed_p(*,i)/sed(*,i), line=ls(i+1) + ENDFOR + oplot, x, sed_p(*,ntype)/sed(*,ntype),col=3 + ENDIF + +; +; Polarization fraction in extinction +; + IF (SHOW(0) NE '0') AND (C_POLEXT EQ 1) THEN BEGIN +; plot fractional polarisation + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(9), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: p/tau '+strtrim(wn(9),2) + endelse + yscl = 1. + xr=[1.e-1,1.e4] + yr=[0.,0.5] + xtit = textoidl('\lambda (\mum)') + ytit = textoidl('p/\tau') + tit='' + plot_oi, x, yscl*(spabs(*,0)+spsca(*,0))/sext(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit + FOR i=1,ntype-1 DO begin + oplot, x, yscl*(spabs(*,i)+spsca(*,i))/sext(*,i), line=ls(i+1) + ENDFOR + oplot, x, yscl*(spabs(*,ntype)+spsca(*,ntype))/sext(*,ntype), col=3 + ENDIF + +; +; Alignment function +; + ip = WHERE( STRPOS(show,'align') GE 0, c_align ) + IF (SHOW(0) NE '0') AND (C_ALIGN EQ 1) THEN BEGIN + fname = fortran_path + 'data/ALIGN.DAT' + nlines = FILE_LINES( fname ) + OPENR, uu, fname, /get_lun + tmp = '#' + print,'(W) DUSTEM_SHOW_FORTRAN: GRAIN.DAT' + cnt = 0 + WHILE STRPOS(tmp,'#') EQ 0 do begin + READF, uu, tmp + cnt = cnt + 1 + ENDWHILE + flags = tmp + readf, uu, anisG0 + + if HARD(0) NE '0' then begin + device, file=hard(ihard), /portrait,/color + ihard = ihard + 1 + endif else begin + window, wn(10), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: f_align '+strtrim(wn(10),2) + endelse + + ; Parametric + readf, uu, tmp + tt = strsplit(tmp, /extract) + align_type = tt(0) + if (align_type eq 'par') then begin + athresh = tt(1) + pstiff = tt(2) + plev = tt(3) + n = 100 + ; Grain radius in microns + radmin = 1d-3 + radmax = 10 + radius = radmin * exp(indgen(n)*alog(radmax/radmin)/n) + fpol = 0.5 * plev * (1 + TANH(ALOG(radius/athresh) / pstiff ) ) + xr=[radmin,radmax] + yr = [0,1.05] + xtit=textoidl('Grain radius (\mum)') + ytit='Alignment efficiency' + plot_oi,radius,fpol,xr=xr,yr=yr,/ys,xtit=xtit,ytit=ytit + endif + + ENDIF + + +; get all the chi2 + if n_elements(smdat) EQ 0 then npar = ntype + smdat = DUSTEM_FIL_CHI2( smdat,ntype=ntype+1 ) + stag = TAG_NAMES(smdat) + itx = WHERE( STRPOS(stag,'I_') GE 0, ctx ) + print,'Chi-square of model to :' + for ii=0,ctx-1 do print,stag(itx(ii))+': ', strtrim(smdat.(itx(ii)).chi2,2) + print,'Chi-square of model to all data : ', strtrim(smdat.chi2,2) + +; +; reset defaults +; + !x.tickname = '' + !y.tickname = '' + + if HARD(0) NE '0' then begin + device, /close + set_plot,'x' + print,'(W): hardcopies in ',hard + !y.thick = 0 + !x.thick = 0 + !p.thick = 0 + !p.charsize = 1 + !p.charthick = 0 + !x.ticklen = 0.02 + endif + + RETURN + the_end: + END + + + + diff --git a/src/idl/dustem_str_inst.pro b/src/idl/dustem_str_inst.pro new file mode 100644 index 0000000..3a315ad --- /dev/null +++ b/src/idl/dustem_str_inst.pro @@ -0,0 +1,36 @@ +FUNCTION DUSTEM_STR_INST, n1, n2=n2, n3=n3 +; returns structure containing flux and CC in band +; N1 (I): number of data points +; N2 (I): nr of grain types or different models +; N3 (I): nr of points in transmission (band flux data) + + if n_elements(n2) EQ 0 then n2 = 1 + + if n_elements(n3) NE 0 then begin + strct = { NAME : strarr(n1), $ + X : dblarr(n1), $ + YD : dblarr(n1), $ + ERR : dblarr(n1), $ + YM : dblarr(n1,n2), $ + FLX : dblarr(n1,n2), $ + CC : dblarr(n1,n2), $ + RR : dblarr(n1), $ + ISEL : intarr(n1)+1, $ + UNIT : '', $ + NPAR : 0, $ + CHI2 : 0d } + s1 = CREATE_STRUCT('TRANS', {X : dblarr(n1,n3), Y : dblarr(n1,n3) }) + strct = CREATE_STRUCT( strct, s1 ) + endif else if n_elements(BAND) EQ 0 then begin + strct = { X : dblarr(n1), $ + YD : dblarr(n1), $ + ERR : dblarr(n1) , $ + YM : dblarr(n1,n2), $ + ISEL : intarr(n1)+1, $ + UNIT : '', $ + NPAR : 0, $ + CHI2 : 0d } + endif + + RETURN, strct +END -- libgit2 0.21.2