diff --git a/src/idl/dustem_str_tools.pro b/src/idl/dustem_str_tools.pro new file mode 100755 index 0000000..603f4e3 --- /dev/null +++ b/src/idl/dustem_str_tools.pro @@ -0,0 +1,177 @@ +FUNCTION 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 + +FUNCTION 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 + +FUNCTION 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 + +FUNCTION FIL_CHI2, sf, ntype=ntype +; fills in the Chi-square fields of an input SMDAT structure +; (see format in 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) 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 = $ + 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 = CHI2( yy(it), mm(it), sf.npar, err=ee(it)) + + return, sf +END + +FUNCTION 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) ADD_MOD: array of sizes must be 2D' + endelse + if ctg EQ 0 then begin + if n_elements(ntrans) EQ 0 then st = CREATE_STRUCT( tag, STR_INST(nx,n2=ntype) ) $ + else st = CREATE_STRUCT( tag, 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) ADD_INST: INST already exists in structure' + endelse + return,sf +END -- libgit2 0.21.2