first commit

Annie Hughes
1 parent 56f461b1
Showing 1 changed file with 177 additions and 0 deletions Show diff stats
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