FUNCTION DUSTEM_GET_BAND_FLUX, dpath, inst, xs=xs, ys=ys, smi=smi, unit=unit, VERBOSE=verbose

  if n_params() EQ 0 then begin 
     print,'-------------------------------------------------------------------------------------------------------------------------'
     print,'FUNCTION GET_BAND_FLUX, dpath, inst, xs=xs, ys=ys, smi=smi, unit=unit, VERBOSE=verbose'
     print,'-------------------------------------------------------------------------------------------------------------------------'
     print,' returns new structure SMDAT with color corrected (cc) flux (MJy/sr) in instrument filter.'
     print,' or fills band flux in input structure SMI containing M_EMIS field'
     print,' Band flux integration always done in lambda, Deviations less than 1% for transmissions given vs. frequency.'
     print,''
     print,' SMDAT = { M_EMIS:  UNIT: for SED (erg/s/cm2/sr)' 
     print,'                    X:  xgrid for SED (microns) '
     print,'                    Y:  input SED '
     print,'                    YP: input polarized SED  '
     print,'           COM:  string for doc'
     print,'           NPAR: nr of parameters for overall CHI2'
     print,'           CHI2: chi-square on all data'
     print,'           INST: NAME: names of bands, '
     print,'                 X:    band center in microns, '                  
     print,'                 YD:   data points (erg/s/cm2/sr)'
     print,'                 ERR:  error on YD same unit'
     print,'                 YM:   SED of model in bands (erg/s/cm2/sr), '
     print,'                 FLX:  flux of model in bands (MJy/sr), '
     print,'                 CC:   color correction in bands, '
     print,'                 RR:   band sampling ratio flux/filt, '                   
     print,'                 ISEL: mask for band selection, 0:discard band, 1=take band, '
     print,'                 UNIT: info on units for X,YD,YM,FLX'
     print,'                 CHI2: chi-square on current bands masked with ISEL, '
     print,'                 TRANS: X: grid for band transmission, '
     print,'                        Y: band transmission } '
     print,''
     print,' DPATH  (I): path for FILTERS directory'
     print,' INST   (I): array of labels for band set, so far CAM, DIRBE, HFI, IRAC, IRAS, LFI, MIPS, PACS, SPIRE, WMAP'
     print,' XS     (I): SMDAT creation, input wave array in microns '
     print,' YS     (I): SMDAT creation, input SED in erg/s/cm2/sr, array(nwave,ntype) allows for several SEDs (grain type)'
     print,' SMI    (I): input SMDAT structure, must at least contain M_EMIS field. If set GET_BAND_FLUX will return filled SMI'
     print,''
     print,' Dependencies: '
     print,'  DUSTEM_BAND_CC, DUSTEM_ADD_MOD, HFI_READ_BANDPASS, GET_HFIBOLO_LIST, HFI_COLOR_CORRECTION '
     print,''
     print,' Examples: '
     print,' SMDAT creation: sm = GET_BAND_FLUX(dpath, [''iras'',''hfi''], xs=xs, ys=ys)'
     print,' SMDAT filling:  sm = GET_BAND_FLUX(dpath, [''spire''],smi=smi)               -> uses xs and ys from SMI (M_EMIS.X and Y)'
     print,' SMDAT filling:  sm = GET_BAND_FLUX(dpath, [''spire''],xs=xs,ys=ys,smi=smi)   -> force xs and ys in SMI'
     print,''
     print,' created by L Verstraete, IAS, Oct-dec 2010, Bandpass data from FILTERS'
     print,''
     print,'-------------------------------------------------------------------------------------------------------------------------'
     return,0
  endif

; init & check
  clight = 2.9979246D10  ; cm/s
  inst = STRUPCASE(STRTRIM(inst,2))
  n_inst = n_elements(inst)
  i_list = [ 'CAM', 'DIRBE', 'HFI', 'IRAC', 'IRAS', 'LFI', 'MIPS', 'PACS', 'SPIRE', 'WMAP' ]
  cx = 0
  for ik = 0, n_inst-1 do begin
     iix = WHERE( i_list EQ inst(ik), cxx )
     cx = cx + cxx
  endfor
  if cx EQ 0 then begin
     print,'(F) DUSTEM_GET_BAND_FLUX: no valid instrument name'
     return,-1
  endif
; find whether this is a creation
  ntg = N_TAGS(smi) 
  if N_TAGS(smi) EQ 0 then begin
     new_str = 1
  endif else begin
     tgs = TAG_NAMES(smi)
     ix = WHERE( tgs EQ 'M_EMIS', cx)
     if cx EQ 0 then new_str = 1 else new_str = 0
  endelse
; get start structure
  if new_str EQ 1 then begin        ; SMDAT creation
     tt = SIZE(ys) 
     if tt(0) GT 1 then ntype = tt(tt(0)) else ntype = 1
     if (n_elements(xs)+n_elements(ys)) EQ 0 OR (n_elements(xs) NE tt(1))  then begin
        print,'(F) DUSTEM_GET_BAND_FLUX: creating SMDAT must define xs(n1), ys(n1,n2)'
        return,0
     endif
     ix = SORT(xs)
     xw = xs(ix) & yw = ys(ix,0:ntype-1)
     smdat = { COM : '', NPAR : 0, CHI2 : 0d }
     smdat = DUSTEM_ADD_MOD(smdat,'emis',[n_elements(xw),ntype])
  endif else if new_str EQ 0 then begin  ; SMDAT filling 
     smdat = smi
     xw = smdat.m_emis.x
     yw = smdat.m_emis.y 
     if n_elements(xs) EQ n_elements(xw) then good = 1 
     if n_elements(ys) EQ n_elements(yw) then good = good + 1 
     if good EQ 2 then begin 
        ix = SORT(xs)
        tt = SIZE(ys)
        if tt(0) GT 1 then ntype = tt(tt(0)) else ntype = 1
        xw = xs(ix) & yw = ys(ix,0:ntype-1)
        if keyword_set(verbose) then print,'(W) DUSTEM_GET_BAND_FLUX: forcing xs, ys in SMI'
     endif 
     if n_elements(unit) EQ 0 then unit = smi.m_emis.unit
  endif
  if n_elements(unit) EQ 0 then unit = 'x(microns) SED(erg/s/cm2/sr)' 
  smdat.m_emis.unit = unit
  smdat.m_emis.x = xw
  smdat.m_emis.y = yw
  tag = TAG_NAMES(smdat)


; instrument loop
  for ixs = 0, n_inst-1 do begin
     if (INST(ixs) EQ 'CAM') then begin
        Nband_sw = 11
        Nband_lw = 10
        nband = nband_sw + nband_lw  ; nr of bands
        ntrans = 160                 ; max nr of points in transmission grid
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif 
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        RESTORE, dpath+'ISO/cam_filter_transmission.save'
        smdat.(itx).x = [ sw_lambda_lc(0:nband_sw-1), lw_lambda_lc(0:nband_lw-1) ]
        smdat.(itx).name = [ strupcase(sw_name_arr(0:nband_sw-1)), strupcase(lw_name_arr(0:nband_lw-1)) ]
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        FOR i=0,Nband_sw-1 DO BEGIN
           ind = where(sw_lambda_arr(*,i) NE 0.)
           xt = sw_lambda_arr(ind,i)
           yt = sw_trans_arr(ind,i)
           nu0 = 1d4*clight/sw_lambda_lc[i]
           for itp = 0, ntype-1 do begin
              smdat.(itx).ym(i,itp) = DUSTEM_BAND_CC( sw_lambda_lc[i], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(i,itp) =  1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(i,itp) = tc
           endfor
           smdat.(itx).rr(i) = rb
           smdat.(itx).trans.x(i,0:nt-1) = xt
           smdat.(itx).trans.y(i,0:nt-1) = yt
        ENDFOR
        FOR i=0,Nband_lw-1 DO BEGIN
           ind=where(lw_lambda_arr(*,i) NE 0.)
           xt = lw_lambda_arr(ind,i)
           yt = lw_trans_arr(ind,i)/max(lw_trans_arr(ind,i))
           nu0 = 1d4*clight/lw_lambda_lc[i]
           for itp = 0, ntype-1 do begin
              smdat.(itx).ym(nband_sw+i,itp) = DUSTEM_BAND_CC( lw_lambda_lc[i], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(nband_sw+i,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(nband_sw+i,itp) = tc
           endfor
           smdat.(itx).rr(nband_sw+i) = rb
           smdat.(itx).trans.x(nband_sw+i,0:nt-1) = xt
           smdat.(itx).trans.y(nband_sw+i,0:nt-1) = yt
        ENDFOR ; CAM 
     ENDIF else if (INST(ixs) EQ 'DIRBE') then begin 
        nband = 10              ; nr of bands
        ntrans = 160            ; max nr of points in transmission grid
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband, n2=ntype, n3=ntrans) )
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        dir_band = dpath + 'DIRBE/'
        fname = dir_band + 'dirbe_response.asc'
        nwav = 800
        wav = dblarr(nwav)
        trans = dblarr(nwav,nband)
        OPENR, uu, fname, /get_lun
        tmp=''
        for i=0,11 do READF, uu, tmp
        for i=0,nwav-1 do begin 
           READF, uu, tmp
           tt = double( strsplit(tmp, ' ', /extract) ) ; first wave then 10 filters
           wav(i) = tt(0)
           trans(i,*) = tt(1:nband) 
        endfor
        CLOSE, uu
        free_lun, uu
        wbd = [1.25,2.2,3.5,4.9,12,25.,60.,100.,140.,240.]
        smdat.(itx).x = wbd
        smdat.(itx).name = 'DIRBE'+strtrim(indgen(10)+1,2)
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        for k=0,nband-1 do begin 
           nu0 = 1d4*clight/wbd[k]
           xt=wav & yt=trans(*,k)
           for itp = 0, ntype-1 do begin
              smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(k,itp) = tc
           endfor 
           smdat.(itx).rr(k) = rb
           smdat.(itx).trans.x(k,0:nt-1) = xt
           smdat.(itx).trans.y(k,0:nt-1) = yt
        endfor ; DIRBE
     ENDIF else if (INST(ixs) EQ 'HFI') then begin
        nband=6
        hfi_bp_version = 'v101'
        ntrans = 800
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        dir_band = dpath + 'HFI/OFFICIAL/'
        wbd = double([100.,143.,217.,353.,545.,857.]) ; band center in GHz
        bnm = STRTRIM(FIX(wbd),2)
        smdat.(itx).name = bnm
        smdat.(itx).x = 1d-5*clight/wbd
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
        bp = HFI_READ_BANDPASS(hfi_bp_version,path_rimo=dir_band,/rimo)
        NAME   = bp.name        ; name of the detector
        FREQ   = bp.freq/1d9    ; frequencies measured in GHz
        TRANS  = bp.trans       ; values of the transmission (normalized later)
        nu = 1d-5*clight/xw     ; microns to GHz
        ixw = SORT(nu) & nu = nu(ixw)
        FOR k=0,nband-1 do begin         
           BOLO_LIST = GET_HFIBOLO_LIST( channel = bnm(k) )
           N_BOLO_LIST = n_elements( bolo_list )  
           yt = 0d0
           For idet = 0 , n_bolo_list - 1 do begin ; ----- get transmission ------------
              Q_DET = where(bolo_list[ idet ] eq name)
              Q_DET = q_det(0)
              Q_TRANS = where( trans(*,q_det) gt 1d-6, nt )
              FREQ_D = freq( q_trans, q_det ) ; already in GHz
              TRANS_D = trans( q_trans, q_det )
              if idet EQ 0 then begin
                 freq_0 = freq_d
                 n0 = n_elements(freq_0)
                 yt = yt + trans_d
              endif else yt = yt + INTERPOL( trans_d, freq_d, freq_0)
           endfor
           smdat.(itx).trans.x(k,0:n0-1) = 1d-5*clight/freq_0
           smdat.(itx).trans.y(k,0:n0-1) = yt / MAX(yt)
           for itp=0,ntype-1 do begin
              inu = 1d17 * yw(ixw,itp)/nu/1d9 ; MJy/sr 
              y0 = INTERPOL(yw(ixw,itp),nu,wbd(k))
              cc = HFI_COLOR_CORRECTION( bp,bnm(k),bolo_cc, nu,inu )
              smdat.(itx).ym(k,itp)  = cc.cc * y0                  ; SED (erg/cm2/s/sr) @ band center
              smdat.(itx).flx(k,itp) = 1d17*cc.cc*y0/wbd(k)/1d9    ; in MJy/sr
              smdat.(itx).cc(k,itp)  = cc.cc
           endfor 
        ENDFOR       
     ENDIF ELSE if (INST(ixs) EQ 'IRAC') then begin
        nband = 4
        ntrans = 410
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        wbd = [ 3.550, 4.493, 5.731, 7.872 ] ; in microns 
        smdat.(itx).x = wbd
        bnm = ['IRAC1','IRAC2','IRAC3','IRAC4']
        smdat.(itx).name = bnm
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        dir_band = dpath + 'Spitzer/'
        for k = 0, nband-1 do begin 
           fn = dir_band + STRLOWCASE(bnm(k)) + '.dat'
           READCOL, fn, xt, yt, /sil ; !!! response in electrons per photon !!!
           yt = yt*xt
           nu0 = 1d4*clight/wbd[k]
           for itp=0,ntype-1 do begin
              smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(k,itp) = tc
           endfor 
           smdat.(itx).trans.x(k,0:nt-1) = xt
           smdat.(itx).trans.y(k,0:nt-1) = yt
           smdat.(itx).rr(k) = rb
        endfor ; IRAC
     ENDIF ELSE if (INST(ixs) EQ 'IRAS') then begin
        nband = 4
        ntrans = 140
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
           ttg = TAG_NAMES(smdat)
           itx = WHERE( ttg EQ 'I_'+inst(ixs))
           bnm = ['12','25','60','100']
           smdat.(itx).name = bnm
           wbd = DOUBLE(bnm)
           smdat.(itx).x = wbd
           smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
           dir_band = dpath + 'IRAS/IRAS'
           for k=0,nband-1 do begin 
              fn = dir_band + bnm(k)+'.txt'
              READCOL, fn, xt, yt,/SIL
              nu0 = 1d4*clight/wbd[k]
              for itp=0,ntype-1 do begin
                 smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
                 smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
                 smdat.(itx).cc(k,itp) = tc
              endfor
              smdat.(itx).trans.x(k,0:nt-1) = xt
              smdat.(itx).trans.y(k,0:nt-1) = yt
              smdat.(itx).rr(k) = rb
           endfor
        endif else begin ; IRAS no read
           ttg = TAG_NAMES(smdat)
           itx = WHERE( ttg EQ 'I_'+inst(ixs))
           wbd = smdat.(itx).x
           for k=0,nband-1 do begin 
              nu0 = 1d4*clight/wbd[k]
              xt = smdat.(itx).trans.x(k,*)
              yt = smdat.(itx).trans.y(k,*)
              for itp=0,ntype-1 do begin
                 smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
                 smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
                 smdat.(itx).cc(k,itp) = tc
                 smdat.(itx).rr(k) = rb
              endfor
           endfor    ; IRAS
        endelse
     ENDIF else if (INST(ixs) EQ 'LFI') then begin 
        wbd = [ 30.,44.,70. ]   ;band center in GHz
        nband = n_elements(wbd)
        ntrans = 210
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        smdat.(itx).x = 1d-5*clight/wbd 
        bnm = STRTRIM(FIX(wbd),2) + 'GHz'
        smdat.(itx).name = bnm
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        dir_band = dpath + 'LFI/LFI'
        FOR k=0L,Nband-1 DO BEGIN
           lfi_files = FINDFILE(dir_band+'*'+bnm(k)+'*.dat',count=Nfiles)
           FOR i=0L,Nfiles-1 DO BEGIN
              IF i EQ 0 THEN BEGIN
                 READCOL,lfi_files(i),nu,R0D0,R0D1,R1D0,R1D1, /SIL
              ENDIF ELSE BEGIN
                 READCOL,lfi_files(i),nu2,R0D0_2,R0D1_2,R1D0_2,R1D1_2, /SIL
                 R0D0 = R0D0 + INTERPOL(R0D0_2,nu2,nu)
                 R1D0 = R1D0 + INTERPOL(R1D0_2,nu2,nu)
                 R0D1 = R0D1 + INTERPOL(R0D1_2,nu2,nu)
                 R1D1 = R1D1 + INTERPOL(R1D1_2,nu2,nu)
              ENDELSE
           ENDFOR
           trans=sqrt(R0D0^2+R1D1^2+R1D0^2+R0D1^2) ; no idea if this is how it should be done
           xt = 1d-5*clight/nu
           yt = trans
           w0 = 1d-5*clight/wbd[k]
           nu0 = 1d9*wbd[k]
           for itp=0,ntype-1 do begin
              smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( w0, xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(k,itp) = tc
           endfor
           smdat.(itx).trans.x(k,0:nt-1) = xt
           smdat.(itx).trans.y(k,0:nt-1) = yt
           smdat.(itx).rr(k) = rb
        ENDFOR ; LFI
     ENDIF else if (INST(ixs) EQ 'MIPS') then begin
        bnm = [ '24', '70', '160' ] ; band name
        wbd = [ 23.68, 71.42, 155.90 ] ; band center in microns
        nband = n_elements(wbd)
        ntrans = 400
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        smdat.(itx).x = wbd
        smdat.(itx).name = bnm
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        dir_band = dpath + 'Spitzer/MIPSfiltsumm_'
        for k = 0, nband-1 do begin 
           fn = dir_band + bnm(k) + '.txt'
           READCOL, fn, xt, yt, /sil ; ; !!! response in electrons per photon energy !!!
           nu0 = 1d4*clight/wbd[k]
           for itp=0,ntype-1 do begin
              smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(k,itp) = tc
           endfor
           smdat.(itx).trans.x(k,0:nt-1) = xt
           smdat.(itx).trans.y(k,0:nt-1) = yt
           smdat.(itx).rr(k) = rb
        endfor ; MIPS
     ENDIF else if (INST(ixs) EQ 'PACS') then begin 
;        wbd=[75.,110.,170.]     ;microns
        wbd=[70.,100.,160.]     ;microns
        nband = n_elements(wbd)
        ntrans = 2030
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        smdat.(itx).x = wbd
        bnm = STRTRIM(FIX(wbd),2)
        smdat.(itx).name = bnm
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        dir_band = dpath + 'PACS/'
        file = dir_band+'PCalPhotometer_Absorption_FM_v2.fits'
        wave=mrdfits(file,1,h,/sil)
        Trans=mrdfits(file,2,h,/sil)
        file=dir_band+'PCalPhotometer_FilterTransmission_FM_v1.fits'
        st3=mrdfits(file,1,h,/sil) & lamb3=st3.wavelength & T3=st3.transmission
        st1=mrdfits(file,2,h,/sil) & lamb1=st1.wavelength & T1=st1.transmission
        st2=mrdfits(file,3,h,/sil) & lamb2=st2.wavelength & T2=st2.transmission
        mx = MAX([n_elements(lamb1),n_elements(lamb2),n_elements(lamb3)])
        yt  = dblarr(nband,mx)
        xt  = dblarr(nband,mx)
        TT1=interpol(Trans,wave,lamb1)
        TT2=interpol(Trans,wave,lamb2)
        TT3=interpol(Trans,wave,lamb3)
        yt(0,0:n_elements(lamb1)-1)=T1*TT1 & xt(0,0:n_elements(lamb1)-1) = lamb1
        yt(1,0:n_elements(lamb2)-1)=T2*TT2 & xt(1,0:n_elements(lamb2)-1) = lamb2 
        yt(2,0:n_elements(lamb3)-1)=T3*TT3 & xt(2,0:n_elements(lamb3)-1) = lamb3
        FOR k=0, nband-1 do begin 
           xtt = REFORM(xt(k,*),mx) & ytt = REFORM(yt(k,*),mx)
           nu0 = 1d4*clight/wbd[k]
           for itp=0,ntype-1 do begin
              smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xtt, ytt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(k,itp) = tc
           endfor
           smdat.(itx).trans.x(k,0:nt-1) = xtt
           smdat.(itx).trans.y(k,0:nt-1) = ytt
           smdat.(itx).rr(k) = rb
        endfor ; PACS
     ENDIF else if (INST(ixs) EQ 'SPIRE') then begin 
        wbd= [250., 350., 500.]
        nband = n_elements(wbd)
        ntrans = 510
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        smdat.(itx).x = wbd
        bnm = STRTRIM(FIX(wbd),2)
        smdat.(itx).name = bnm
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        tt  = dblarr(nband+1,512)   
        dir_band = dpath + 'SPIRE/'
        OPENR,unit,dir_band+'SPIRE_transmissions.txt',/get_lun
        sstr=''
        for i=0,2 do READF,unit,sstr
        cnt = 0 & ss=''
        WHILE not eof(unit) DO BEGIN
           READF,unit,ss
           tt(*,cnt) = double(strsplit(ss,/extract))
           cnt = cnt + 1
        ENDWHILE
        CLOSE, unit
        free_lun, unit
        FOR k=0, nband-1 do begin
           xt = 1d4/tt(0,*)     ; convert to microns
           yt = tt[k+1,*]
           w0 = wbd[k]
           nu0 = 1d4*clight/w0
           for itp=0,ntype-1 do begin
              smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( w0, xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
              smdat.(itx).cc(k,itp) = tc
           endfor
           smdat.(itx).trans.x(k,0:nt-1) = xt
           smdat.(itx).trans.y(k,0:nt-1) = yt
           smdat.(itx).rr(k) = rb
        endfor ; SPIRE
     ENDIF else if (INST(ixs) EQ 'WMAP') then begin 
        wbd = double([ 23., 33., 41., 61., 94. ])
        nband = n_elements(wbd)
        ntrans = 1000
        ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
        if new_str EQ 1 OR ctg EQ 0 then begin
           st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) ) 
           smdat = CREATE_STRUCT( smdat, st )
        endif
        ttg = TAG_NAMES(smdat)
        itx = WHERE( ttg EQ 'I_'+inst(ixs))
        smdat.(itx).x = 1d-5*clight/wbd
        bnm = STRTRIM(FIX(wbd),2)+'GHz'
        smdat.(itx).name = bnm
        smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)' 
        dir_band = dpath + 'WMAP/'
        wmap_files = FINDFILE(dir_band+'*v3.xcat', count=nfiles)
        FOR k=0, nband-1 do begin 
           READCOL, wmap_files(k), nu, trans, skip=4, /sil
           xt = 1d-5*clight/nu  ; !!!! integrate in microns
           yt = trans
           w0 = 1d-5*clight/wbd[k]
           nu0 = 1d9*wbd[k]
           for itp=0,ntype-1 do begin
              smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( w0, xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
              smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr 
              smdat.(itx).cc(k,itp) = tc
           endfor 
           smdat.(itx).trans.x(k,0:nt-1) = xt
           smdat.(itx).trans.y(k,0:nt-1) = yt
           smdat.(itx).rr(k) = rb
        endfor ; WMAP
     ENDIF

  ENDFOR

  RETURN, SMDAT
END