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