OV-GSO-DC / dustem-wrapper_idl

Commit aed288deff4171c3cbb35f4cf18bf34878a923a8

Authored by Annie Hughes
1 parent 4da6e627
Exists in master

first commit of IAS routines

Showing 1 changed file with 517 additions and 0 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_get_band_flux.pro 0 → 100644
  Diff comments   View file @aed288d
... ... @@ -0,0 +1,517 @@
  1 +FUNCTION DUSTEM_GET_BAND_FLUX, dpath, inst, xs=xs, ys=ys, smi=smi, unit=unit, VERBOSE=verbose
  2 +
  3 + if n_params() EQ 0 then begin
  4 + print,'-------------------------------------------------------------------------------------------------------------------------'
  5 + print,'FUNCTION GET_BAND_FLUX, dpath, inst, xs=xs, ys=ys, smi=smi, unit=unit, VERBOSE=verbose'
  6 + print,'-------------------------------------------------------------------------------------------------------------------------'
  7 + print,' returns new structure SMDAT with color corrected (cc) flux (MJy/sr) in instrument filter.'
  8 + print,' or fills band flux in input structure SMI containing M_EMIS field'
  9 + print,' Band flux integration always done in lambda, Deviations less than 1% for transmissions given vs. frequency.'
  10 + print,''
  11 + print,' SMDAT = { M_EMIS: UNIT: for SED (erg/s/cm2/sr)'
  12 + print,' X: xgrid for SED (microns) '
  13 + print,' Y: input SED '
  14 + print,' YP: input polarized SED '
  15 + print,' COM: string for doc'
  16 + print,' NPAR: nr of parameters for overall CHI2'
  17 + print,' CHI2: chi-square on all data'
  18 + print,' INST: NAME: names of bands, '
  19 + print,' X: band center in microns, '
  20 + print,' YD: data points (erg/s/cm2/sr)'
  21 + print,' ERR: error on YD same unit'
  22 + print,' YM: SED of model in bands (erg/s/cm2/sr), '
  23 + print,' FLX: flux of model in bands (MJy/sr), '
  24 + print,' CC: color correction in bands, '
  25 + print,' RR: band sampling ratio flux/filt, '
  26 + print,' ISEL: mask for band selection, 0:discard band, 1=take band, '
  27 + print,' UNIT: info on units for X,YD,YM,FLX'
  28 + print,' CHI2: chi-square on current bands masked with ISEL, '
  29 + print,' TRANS: X: grid for band transmission, '
  30 + print,' Y: band transmission } '
  31 + print,''
  32 + print,' DPATH (I): path for FILTERS directory'
  33 + print,' INST (I): array of labels for band set, so far CAM, DIRBE, HFI, IRAC, IRAS, LFI, MIPS, PACS, SPIRE, WMAP'
  34 + print,' XS (I): SMDAT creation, input wave array in microns '
  35 + print,' YS (I): SMDAT creation, input SED in erg/s/cm2/sr, array(nwave,ntype) allows for several SEDs (grain type)'
  36 + print,' SMI (I): input SMDAT structure, must at least contain M_EMIS field. If set GET_BAND_FLUX will return filled SMI'
  37 + print,''
  38 + print,' Dependencies: '
  39 + print,' DUSTEM_BAND_CC, DUSTEM_ADD_MOD, HFI_READ_BANDPASS, GET_HFIBOLO_LIST, HFI_COLOR_CORRECTION '
  40 + print,''
  41 + print,' Examples: '
  42 + print,' SMDAT creation: sm = GET_BAND_FLUX(dpath, [''iras'',''hfi''], xs=xs, ys=ys)'
  43 + print,' SMDAT filling: sm = GET_BAND_FLUX(dpath, [''spire''],smi=smi) -> uses xs and ys from SMI (M_EMIS.X and Y)'
  44 + print,' SMDAT filling: sm = GET_BAND_FLUX(dpath, [''spire''],xs=xs,ys=ys,smi=smi) -> force xs and ys in SMI'
  45 + print,''
  46 + print,' created by L Verstraete, IAS, Oct-dec 2010, Bandpass data from FILTERS'
  47 + print,''
  48 + print,'-------------------------------------------------------------------------------------------------------------------------'
  49 + return,0
  50 + endif
  51 +
  52 +; init & check
  53 + clight = 2.9979246D10 ; cm/s
  54 + inst = STRUPCASE(STRTRIM(inst,2))
  55 + n_inst = n_elements(inst)
  56 + i_list = [ 'CAM', 'DIRBE', 'HFI', 'IRAC', 'IRAS', 'LFI', 'MIPS', 'PACS', 'SPIRE', 'WMAP' ]
  57 + cx = 0
  58 + for ik = 0, n_inst-1 do begin
  59 + iix = WHERE( i_list EQ inst(ik), cxx )
  60 + cx = cx + cxx
  61 + endfor
  62 + if cx EQ 0 then begin
  63 + print,'(F) DUSTEM_GET_BAND_FLUX: no valid instrument name'
  64 + return,-1
  65 + endif
  66 +; find whether this is a creation
  67 + ntg = N_TAGS(smi)
  68 + if N_TAGS(smi) EQ 0 then begin
  69 + new_str = 1
  70 + endif else begin
  71 + tgs = TAG_NAMES(smi)
  72 + ix = WHERE( tgs EQ 'M_EMIS', cx)
  73 + if cx EQ 0 then new_str = 1 else new_str = 0
  74 + endelse
  75 +; get start structure
  76 + if new_str EQ 1 then begin ; SMDAT creation
  77 + tt = SIZE(ys)
  78 + if tt(0) GT 1 then ntype = tt(tt(0)) else ntype = 1
  79 + if (n_elements(xs)+n_elements(ys)) EQ 0 OR (n_elements(xs) NE tt(1)) then begin
  80 + print,'(F) DUSTEM_GET_BAND_FLUX: creating SMDAT must define xs(n1), ys(n1,n2)'
  81 + return,0
  82 + endif
  83 + ix = SORT(xs)
  84 + xw = xs(ix) & yw = ys(ix,0:ntype-1)
  85 + smdat = { COM : '', NPAR : 0, CHI2 : 0d }
  86 + smdat = DUSTEM_ADD_MOD(smdat,'emis',[n_elements(xw),ntype])
  87 + endif else if new_str EQ 0 then begin ; SMDAT filling
  88 + smdat = smi
  89 + xw = smdat.m_emis.x
  90 + yw = smdat.m_emis.y
  91 + if n_elements(xs) EQ n_elements(xw) then good = 1
  92 + if n_elements(ys) EQ n_elements(yw) then good = good + 1
  93 + if good EQ 2 then begin
  94 + ix = SORT(xs)
  95 + tt = SIZE(ys)
  96 + if tt(0) GT 1 then ntype = tt(tt(0)) else ntype = 1
  97 + xw = xs(ix) & yw = ys(ix,0:ntype-1)
  98 + if keyword_set(verbose) then print,'(W) DUSTEM_GET_BAND_FLUX: forcing xs, ys in SMI'
  99 + endif
  100 + if n_elements(unit) EQ 0 then unit = smi.m_emis.unit
  101 + endif
  102 + if n_elements(unit) EQ 0 then unit = 'x(microns) SED(erg/s/cm2/sr)'
  103 + smdat.m_emis.unit = unit
  104 + smdat.m_emis.x = xw
  105 + smdat.m_emis.y = yw
  106 + tag = TAG_NAMES(smdat)
  107 +
  108 +
  109 +; instrument loop
  110 + for ixs = 0, n_inst-1 do begin
  111 + if (INST(ixs) EQ 'CAM') then begin
  112 + Nband_sw = 11
  113 + Nband_lw = 10
  114 + nband = nband_sw + nband_lw ; nr of bands
  115 + ntrans = 160 ; max nr of points in transmission grid
  116 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  117 + if new_str EQ 1 OR ctg EQ 0 then begin
  118 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  119 + smdat = CREATE_STRUCT( smdat, st )
  120 + endif
  121 + ttg = TAG_NAMES(smdat)
  122 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  123 + RESTORE, dpath+'ISO/cam_filter_transmission.save'
  124 + smdat.(itx).x = [ sw_lambda_lc(0:nband_sw-1), lw_lambda_lc(0:nband_lw-1) ]
  125 + smdat.(itx).name = [ strupcase(sw_name_arr(0:nband_sw-1)), strupcase(lw_name_arr(0:nband_lw-1)) ]
  126 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  127 + FOR i=0,Nband_sw-1 DO BEGIN
  128 + ind = where(sw_lambda_arr(*,i) NE 0.)
  129 + xt = sw_lambda_arr(ind,i)
  130 + yt = sw_trans_arr(ind,i)
  131 + nu0 = 1d4*clight/sw_lambda_lc[i]
  132 + for itp = 0, ntype-1 do begin
  133 + 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 )
  134 + smdat.(itx).flx(i,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  135 + smdat.(itx).cc(i,itp) = tc
  136 + endfor
  137 + smdat.(itx).rr(i) = rb
  138 + smdat.(itx).trans.x(i,0:nt-1) = xt
  139 + smdat.(itx).trans.y(i,0:nt-1) = yt
  140 + ENDFOR
  141 + FOR i=0,Nband_lw-1 DO BEGIN
  142 + ind=where(lw_lambda_arr(*,i) NE 0.)
  143 + xt = lw_lambda_arr(ind,i)
  144 + yt = lw_trans_arr(ind,i)/max(lw_trans_arr(ind,i))
  145 + nu0 = 1d4*clight/lw_lambda_lc[i]
  146 + for itp = 0, ntype-1 do begin
  147 + 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 )
  148 + smdat.(itx).flx(nband_sw+i,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  149 + smdat.(itx).cc(nband_sw+i,itp) = tc
  150 + endfor
  151 + smdat.(itx).rr(nband_sw+i) = rb
  152 + smdat.(itx).trans.x(nband_sw+i,0:nt-1) = xt
  153 + smdat.(itx).trans.y(nband_sw+i,0:nt-1) = yt
  154 + ENDFOR ; CAM
  155 + ENDIF else if (INST(ixs) EQ 'DIRBE') then begin
  156 + nband = 10 ; nr of bands
  157 + ntrans = 160 ; max nr of points in transmission grid
  158 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  159 + if new_str EQ 1 OR ctg EQ 0 then begin
  160 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband, n2=ntype, n3=ntrans) )
  161 + smdat = CREATE_STRUCT( smdat, st )
  162 + endif
  163 + ttg = TAG_NAMES(smdat)
  164 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  165 + dir_band = dpath + 'DIRBE/'
  166 + fname = dir_band + 'dirbe_response.asc'
  167 + nwav = 800
  168 + wav = dblarr(nwav)
  169 + trans = dblarr(nwav,nband)
  170 + OPENR, uu, fname, /get_lun
  171 + tmp=''
  172 + for i=0,11 do READF, uu, tmp
  173 + for i=0,nwav-1 do begin
  174 + READF, uu, tmp
  175 + tt = double( strsplit(tmp, ' ', /extract) ) ; first wave then 10 filters
  176 + wav(i) = tt(0)
  177 + trans(i,*) = tt(1:nband)
  178 + endfor
  179 + CLOSE, uu
  180 + free_lun, uu
  181 + wbd = [1.25,2.2,3.5,4.9,12,25.,60.,100.,140.,240.]
  182 + smdat.(itx).x = wbd
  183 + smdat.(itx).name = 'DIRBE'+strtrim(indgen(10)+1,2)
  184 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  185 + for k=0,nband-1 do begin
  186 + nu0 = 1d4*clight/wbd[k]
  187 + xt=wav & yt=trans(*,k)
  188 + for itp = 0, ntype-1 do begin
  189 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  190 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  191 + smdat.(itx).cc(k,itp) = tc
  192 + endfor
  193 + smdat.(itx).rr(k) = rb
  194 + smdat.(itx).trans.x(k,0:nt-1) = xt
  195 + smdat.(itx).trans.y(k,0:nt-1) = yt
  196 + endfor ; DIRBE
  197 + ENDIF else if (INST(ixs) EQ 'HFI') then begin
  198 + nband=6
  199 + hfi_bp_version = 'v101'
  200 + ntrans = 800
  201 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  202 + if new_str EQ 1 OR ctg EQ 0 then begin
  203 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  204 + smdat = CREATE_STRUCT( smdat, st )
  205 + endif
  206 + ttg = TAG_NAMES(smdat)
  207 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  208 + dir_band = dpath + 'HFI/OFFICIAL/'
  209 + wbd = double([100.,143.,217.,353.,545.,857.]) ; band center in GHz
  210 + bnm = STRTRIM(FIX(wbd),2)
  211 + smdat.(itx).name = bnm
  212 + smdat.(itx).x = 1d-5*clight/wbd
  213 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  214 + bp = HFI_READ_BANDPASS(hfi_bp_version,path_rimo=dir_band,/rimo)
  215 + NAME = bp.name ; name of the detector
  216 + FREQ = bp.freq/1d9 ; frequencies measured in GHz
  217 + TRANS = bp.trans ; values of the transmission (normalized later)
  218 + nu = 1d-5*clight/xw ; microns to GHz
  219 + ixw = SORT(nu) & nu = nu(ixw)
  220 + FOR k=0,nband-1 do begin
  221 + BOLO_LIST = GET_HFIBOLO_LIST( channel = bnm(k) )
  222 + N_BOLO_LIST = n_elements( bolo_list )
  223 + yt = 0d0
  224 + For idet = 0 , n_bolo_list - 1 do begin ; ----- get transmission ------------
  225 + Q_DET = where(bolo_list[ idet ] eq name)
  226 + Q_DET = q_det(0)
  227 + Q_TRANS = where( trans(*,q_det) gt 1d-6, nt )
  228 + FREQ_D = freq( q_trans, q_det ) ; already in GHz
  229 + TRANS_D = trans( q_trans, q_det )
  230 + if idet EQ 0 then begin
  231 + freq_0 = freq_d
  232 + n0 = n_elements(freq_0)
  233 + yt = yt + trans_d
  234 + endif else yt = yt + INTERPOL( trans_d, freq_d, freq_0)
  235 + endfor
  236 + smdat.(itx).trans.x(k,0:n0-1) = 1d-5*clight/freq_0
  237 + smdat.(itx).trans.y(k,0:n0-1) = yt / MAX(yt)
  238 + for itp=0,ntype-1 do begin
  239 + inu = 1d17 * yw(ixw,itp)/nu/1d9 ; MJy/sr
  240 + y0 = INTERPOL(yw(ixw,itp),nu,wbd(k))
  241 + cc = HFI_COLOR_CORRECTION( bp,bnm(k),bolo_cc, nu,inu )
  242 + smdat.(itx).ym(k,itp) = cc.cc * y0 ; SED (erg/cm2/s/sr) @ band center
  243 + smdat.(itx).flx(k,itp) = 1d17*cc.cc*y0/wbd(k)/1d9 ; in MJy/sr
  244 + smdat.(itx).cc(k,itp) = cc.cc
  245 + endfor
  246 + ENDFOR
  247 + ENDIF ELSE if (INST(ixs) EQ 'IRAC') then begin
  248 + nband = 4
  249 + ntrans = 410
  250 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  251 + if new_str EQ 1 OR ctg EQ 0 then begin
  252 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  253 + smdat = CREATE_STRUCT( smdat, st )
  254 + endif
  255 + ttg = TAG_NAMES(smdat)
  256 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  257 + wbd = [ 3.550, 4.493, 5.731, 7.872 ] ; in microns
  258 + smdat.(itx).x = wbd
  259 + bnm = ['IRAC1','IRAC2','IRAC3','IRAC4']
  260 + smdat.(itx).name = bnm
  261 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  262 + dir_band = dpath + 'Spitzer/'
  263 + for k = 0, nband-1 do begin
  264 + fn = dir_band + STRLOWCASE(bnm(k)) + '.dat'
  265 + READCOL, fn, xt, yt, /sil ; !!! response in electrons per photon !!!
  266 + yt = yt*xt
  267 + nu0 = 1d4*clight/wbd[k]
  268 + for itp=0,ntype-1 do begin
  269 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  270 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  271 + smdat.(itx).cc(k,itp) = tc
  272 + endfor
  273 + smdat.(itx).trans.x(k,0:nt-1) = xt
  274 + smdat.(itx).trans.y(k,0:nt-1) = yt
  275 + smdat.(itx).rr(k) = rb
  276 + endfor ; IRAC
  277 + ENDIF ELSE if (INST(ixs) EQ 'IRAS') then begin
  278 + nband = 4
  279 + ntrans = 140
  280 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  281 + if new_str EQ 1 OR ctg EQ 0 then begin
  282 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  283 + smdat = CREATE_STRUCT( smdat, st )
  284 + ttg = TAG_NAMES(smdat)
  285 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  286 + bnm = ['12','25','60','100']
  287 + smdat.(itx).name = bnm
  288 + wbd = DOUBLE(bnm)
  289 + smdat.(itx).x = wbd
  290 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  291 + dir_band = dpath + 'IRAS/IRAS'
  292 + for k=0,nband-1 do begin
  293 + fn = dir_band + bnm(k)+'.txt'
  294 + READCOL, fn, xt, yt,/SIL
  295 + nu0 = 1d4*clight/wbd[k]
  296 + for itp=0,ntype-1 do begin
  297 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  298 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  299 + smdat.(itx).cc(k,itp) = tc
  300 + endfor
  301 + smdat.(itx).trans.x(k,0:nt-1) = xt
  302 + smdat.(itx).trans.y(k,0:nt-1) = yt
  303 + smdat.(itx).rr(k) = rb
  304 + endfor
  305 + endif else begin ; IRAS no read
  306 + ttg = TAG_NAMES(smdat)
  307 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  308 + wbd = smdat.(itx).x
  309 + for k=0,nband-1 do begin
  310 + nu0 = 1d4*clight/wbd[k]
  311 + xt = smdat.(itx).trans.x(k,*)
  312 + yt = smdat.(itx).trans.y(k,*)
  313 + for itp=0,ntype-1 do begin
  314 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  315 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  316 + smdat.(itx).cc(k,itp) = tc
  317 + smdat.(itx).rr(k) = rb
  318 + endfor
  319 + endfor ; IRAS
  320 + endelse
  321 + ENDIF else if (INST(ixs) EQ 'LFI') then begin
  322 + wbd = [ 30.,44.,70. ] ;band center in GHz
  323 + nband = n_elements(wbd)
  324 + ntrans = 210
  325 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  326 + if new_str EQ 1 OR ctg EQ 0 then begin
  327 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  328 + smdat = CREATE_STRUCT( smdat, st )
  329 + endif
  330 + ttg = TAG_NAMES(smdat)
  331 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  332 + smdat.(itx).x = 1d-5*clight/wbd
  333 + bnm = STRTRIM(FIX(wbd),2) + 'GHz'
  334 + smdat.(itx).name = bnm
  335 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  336 + dir_band = dpath + 'LFI/LFI'
  337 + FOR k=0L,Nband-1 DO BEGIN
  338 + lfi_files = FINDFILE(dir_band+'*'+bnm(k)+'*.dat',count=Nfiles)
  339 + FOR i=0L,Nfiles-1 DO BEGIN
  340 + IF i EQ 0 THEN BEGIN
  341 + READCOL,lfi_files(i),nu,R0D0,R0D1,R1D0,R1D1, /SIL
  342 + ENDIF ELSE BEGIN
  343 + READCOL,lfi_files(i),nu2,R0D0_2,R0D1_2,R1D0_2,R1D1_2, /SIL
  344 + R0D0 = R0D0 + INTERPOL(R0D0_2,nu2,nu)
  345 + R1D0 = R1D0 + INTERPOL(R1D0_2,nu2,nu)
  346 + R0D1 = R0D1 + INTERPOL(R0D1_2,nu2,nu)
  347 + R1D1 = R1D1 + INTERPOL(R1D1_2,nu2,nu)
  348 + ENDELSE
  349 + ENDFOR
  350 + trans=sqrt(R0D0^2+R1D1^2+R1D0^2+R0D1^2) ; no idea if this is how it should be done
  351 + xt = 1d-5*clight/nu
  352 + yt = trans
  353 + w0 = 1d-5*clight/wbd[k]
  354 + nu0 = 1d9*wbd[k]
  355 + for itp=0,ntype-1 do begin
  356 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( w0, xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  357 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  358 + smdat.(itx).cc(k,itp) = tc
  359 + endfor
  360 + smdat.(itx).trans.x(k,0:nt-1) = xt
  361 + smdat.(itx).trans.y(k,0:nt-1) = yt
  362 + smdat.(itx).rr(k) = rb
  363 + ENDFOR ; LFI
  364 + ENDIF else if (INST(ixs) EQ 'MIPS') then begin
  365 + bnm = [ '24', '70', '160' ] ; band name
  366 + wbd = [ 23.68, 71.42, 155.90 ] ; band center in microns
  367 + nband = n_elements(wbd)
  368 + ntrans = 400
  369 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  370 + if new_str EQ 1 OR ctg EQ 0 then begin
  371 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  372 + smdat = CREATE_STRUCT( smdat, st )
  373 + endif
  374 + ttg = TAG_NAMES(smdat)
  375 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  376 + smdat.(itx).x = wbd
  377 + smdat.(itx).name = bnm
  378 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  379 + dir_band = dpath + 'Spitzer/MIPSfiltsumm_'
  380 + for k = 0, nband-1 do begin
  381 + fn = dir_band + bnm(k) + '.txt'
  382 + READCOL, fn, xt, yt, /sil ; ; !!! response in electrons per photon energy !!!
  383 + nu0 = 1d4*clight/wbd[k]
  384 + for itp=0,ntype-1 do begin
  385 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  386 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  387 + smdat.(itx).cc(k,itp) = tc
  388 + endfor
  389 + smdat.(itx).trans.x(k,0:nt-1) = xt
  390 + smdat.(itx).trans.y(k,0:nt-1) = yt
  391 + smdat.(itx).rr(k) = rb
  392 + endfor ; MIPS
  393 + ENDIF else if (INST(ixs) EQ 'PACS') then begin
  394 +; wbd=[75.,110.,170.] ;microns
  395 + wbd=[70.,100.,160.] ;microns
  396 + nband = n_elements(wbd)
  397 + ntrans = 2030
  398 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  399 + if new_str EQ 1 OR ctg EQ 0 then begin
  400 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  401 + smdat = CREATE_STRUCT( smdat, st )
  402 + endif
  403 + ttg = TAG_NAMES(smdat)
  404 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  405 + smdat.(itx).x = wbd
  406 + bnm = STRTRIM(FIX(wbd),2)
  407 + smdat.(itx).name = bnm
  408 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  409 + dir_band = dpath + 'PACS/'
  410 + file = dir_band+'PCalPhotometer_Absorption_FM_v2.fits'
  411 + wave=mrdfits(file,1,h,/sil)
  412 + Trans=mrdfits(file,2,h,/sil)
  413 + file=dir_band+'PCalPhotometer_FilterTransmission_FM_v1.fits'
  414 + st3=mrdfits(file,1,h,/sil) & lamb3=st3.wavelength & T3=st3.transmission
  415 + st1=mrdfits(file,2,h,/sil) & lamb1=st1.wavelength & T1=st1.transmission
  416 + st2=mrdfits(file,3,h,/sil) & lamb2=st2.wavelength & T2=st2.transmission
  417 + mx = MAX([n_elements(lamb1),n_elements(lamb2),n_elements(lamb3)])
  418 + yt = dblarr(nband,mx)
  419 + xt = dblarr(nband,mx)
  420 + TT1=interpol(Trans,wave,lamb1)
  421 + TT2=interpol(Trans,wave,lamb2)
  422 + TT3=interpol(Trans,wave,lamb3)
  423 + yt(0,0:n_elements(lamb1)-1)=T1*TT1 & xt(0,0:n_elements(lamb1)-1) = lamb1
  424 + yt(1,0:n_elements(lamb2)-1)=T2*TT2 & xt(1,0:n_elements(lamb2)-1) = lamb2
  425 + yt(2,0:n_elements(lamb3)-1)=T3*TT3 & xt(2,0:n_elements(lamb3)-1) = lamb3
  426 + FOR k=0, nband-1 do begin
  427 + xtt = REFORM(xt(k,*),mx) & ytt = REFORM(yt(k,*),mx)
  428 + nu0 = 1d4*clight/wbd[k]
  429 + for itp=0,ntype-1 do begin
  430 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( wbd[k], xw, yw(*,itp), xtt, ytt, y0=y0, cc=tc, rr=rb, ni=nt )
  431 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  432 + smdat.(itx).cc(k,itp) = tc
  433 + endfor
  434 + smdat.(itx).trans.x(k,0:nt-1) = xtt
  435 + smdat.(itx).trans.y(k,0:nt-1) = ytt
  436 + smdat.(itx).rr(k) = rb
  437 + endfor ; PACS
  438 + ENDIF else if (INST(ixs) EQ 'SPIRE') then begin
  439 + wbd= [250., 350., 500.]
  440 + nband = n_elements(wbd)
  441 + ntrans = 510
  442 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  443 + if new_str EQ 1 OR ctg EQ 0 then begin
  444 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  445 + smdat = CREATE_STRUCT( smdat, st )
  446 + endif
  447 + ttg = TAG_NAMES(smdat)
  448 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  449 + smdat.(itx).x = wbd
  450 + bnm = STRTRIM(FIX(wbd),2)
  451 + smdat.(itx).name = bnm
  452 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  453 + tt = dblarr(nband+1,512)
  454 + dir_band = dpath + 'SPIRE/'
  455 + OPENR,unit,dir_band+'SPIRE_transmissions.txt',/get_lun
  456 + sstr=''
  457 + for i=0,2 do READF,unit,sstr
  458 + cnt = 0 & ss=''
  459 + WHILE not eof(unit) DO BEGIN
  460 + READF,unit,ss
  461 + tt(*,cnt) = double(strsplit(ss,/extract))
  462 + cnt = cnt + 1
  463 + ENDWHILE
  464 + CLOSE, unit
  465 + free_lun, unit
  466 + FOR k=0, nband-1 do begin
  467 + xt = 1d4/tt(0,*) ; convert to microns
  468 + yt = tt[k+1,*]
  469 + w0 = wbd[k]
  470 + nu0 = 1d4*clight/w0
  471 + for itp=0,ntype-1 do begin
  472 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( w0, xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  473 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  474 + smdat.(itx).cc(k,itp) = tc
  475 + endfor
  476 + smdat.(itx).trans.x(k,0:nt-1) = xt
  477 + smdat.(itx).trans.y(k,0:nt-1) = yt
  478 + smdat.(itx).rr(k) = rb
  479 + endfor ; SPIRE
  480 + ENDIF else if (INST(ixs) EQ 'WMAP') then begin
  481 + wbd = double([ 23., 33., 41., 61., 94. ])
  482 + nband = n_elements(wbd)
  483 + ntrans = 1000
  484 + ix = WHERE(tag EQ 'I_'+inst(ixs), ctg)
  485 + if new_str EQ 1 OR ctg EQ 0 then begin
  486 + st = CREATE_STRUCT( 'I_'+inst(ixs), DUSTEM_STR_INST(nband,n2=ntype,n3=ntrans) )
  487 + smdat = CREATE_STRUCT( smdat, st )
  488 + endif
  489 + ttg = TAG_NAMES(smdat)
  490 + itx = WHERE( ttg EQ 'I_'+inst(ixs))
  491 + smdat.(itx).x = 1d-5*clight/wbd
  492 + bnm = STRTRIM(FIX(wbd),2)+'GHz'
  493 + smdat.(itx).name = bnm
  494 + smdat.(itx).unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr)'
  495 + dir_band = dpath + 'WMAP/'
  496 + wmap_files = FINDFILE(dir_band+'*v3.xcat', count=nfiles)
  497 + FOR k=0, nband-1 do begin
  498 + READCOL, wmap_files(k), nu, trans, skip=4, /sil
  499 + xt = 1d-5*clight/nu ; !!!! integrate in microns
  500 + yt = trans
  501 + w0 = 1d-5*clight/wbd[k]
  502 + nu0 = 1d9*wbd[k]
  503 + for itp=0,ntype-1 do begin
  504 + smdat.(itx).ym(k,itp) = DUSTEM_BAND_CC( w0, xw, yw(*,itp), xt, yt, y0=y0, cc=tc, rr=rb, ni=nt )
  505 + smdat.(itx).flx(k,itp) = 1.d17*tc*y0/nu0 ; in MJy/sr
  506 + smdat.(itx).cc(k,itp) = tc
  507 + endfor
  508 + smdat.(itx).trans.x(k,0:nt-1) = xt
  509 + smdat.(itx).trans.y(k,0:nt-1) = yt
  510 + smdat.(itx).rr(k) = rb
  511 + endfor ; WMAP
  512 + ENDIF
  513 +
  514 + ENDFOR
  515 +
  516 + RETURN, SMDAT
  517 +END
... ...