;David G. Whelan
;19 January, 2011

;a function that returns the passband flux densities in an evolutionary sequence
;structure, given the spectrum structure and the directory to all the filter 
;bandpasses

;INPUT PARAMETERS
;spec_struct: the structure of SEDs along the evolutionary sequence as returned by wrapper_read_spec.pro
;filter_dir: a string continaing the path to filters. Must contain a trailing `/'
;
;RETURNS:
;This function returns a structure with all of the band flux densities along every line of sight for the evolutionary sequence that spec_struct represents.

FUNCTION FILTERCON, SPEC_STRUCT, FILTER_DIR

;some constants
nlambda = 200
ntheta  = 10
nphi    = 20

;band info
bands=['2MASSJ','2MASSH','2MASSK','IRAC1','IRAC2','IRAC3','IRAC4',$
       'MIPS24','MIPS70','MIPS160','IRAS12','IRAS25','IRAS60','IRAS100',$
       'PACS70','PACS100','PACS160','SPIRE250','SPIRE350','SPIRE500']
nbands=n_elements(bands)
bandwave=[1.25,1.62,2.17,3.6,4.5,5.8,7.74,24,70,160,12,25,60,100,70,100,160,$
          250,350,500]
;all band flux densities will be in F_nu
flux_rin1=fltarr(ninc,nphi,nbands)
flux_rin2=fltarr(ninc,nphi,nbands)
flux_rin3=fltarr(ninc,nphi,nbands)
flux_rin4=fltarr(ninc,nphi,nbands)
flux_rin5=fltarr(ninc,nphi,nbands)
flux_rin6=fltarr(ninc,nphi,nbands)
flux_rin7=fltarr(ninc,nphi,nbands)
flux_rin8=fltarr(ninc,nphi,nbands)
flux_rin9=fltarr(ninc,nphi,nbands)
flux_rin10=fltarr(ninc,nphi,nbands)
flux_rin11=fltarr(ninc,nphi,nbands)
flux_rin12=fltarr(ninc,nphi,nbands)

;load the filters, and interpolate them to the wavelength array of the models
;adopted from sedmulti.pro:
nlamb=n_elements(spec_struct.lambda)
ftrans=fltarr(nbands,nlamb)
dlamb=[spec_struct.lambda[0:nlamb-2]-spec_struct.lambda[1:nlamb-1],$
       spec_struct.lambda[nw-2]-spec_struct.lambda[nw-1]]
for i=0,nbands-1,1 do begin
    readcol,filter_dir+bands[i]+'.filter',flamb0,ftrans0,/silent
    z=where(spec_struct.lambda gt min(flamb0) AND $
            spec_struct.lambda lt max(flamb0))
    z2=minmax(where(ftrans0/max(ftrans0) gt 0.5))
    ftrans[i,z]=interpol(ftrans0,flamb0,spec_struct.lambda[z])
    ftrans[i,z]=ftrans[i,z[/total(ftrans[i,z]*dw[z])
    bandwave[i]=total(spec_struct.lambda[z]*ftrans[i,z]*dw[z])
endfor

for nt=0,ntheta-1,1 do begin
    for np=0,nphi-1,1 do begin
        for i=0,nbands-1 do begin
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin1[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin1[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin2[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin2[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin3[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin3[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin4[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin4[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin5[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin5[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin6[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin6[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin7[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin7[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin8[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin8[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin9[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin9[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin10[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin10[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin11[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin11[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
            flux_lam=total(ftrans[i,*]*spec_struct.lamFlam_rin12[nt,np,*]/$
                           spec_struct.lambda*dw)
            flux_rin12[nt,np,i]=flux_lam*bandwave[i]^2/2.998e14
        endfor
    endfor
endfor

;the return structure:
ret_struct = { $
		SFE:spec_struct.sfe,$
		ROUT:spec_struct.rout,$
		SMTHFRACT:spec_struct.smthfract,$
		RIN:spec_struct.rin_vals,$
		NUMTHETA:ntheta,$
		NUMPHI:nphi,$
		BANDPASS_LAMBDAS:bandwave,$
		bandfluxes_rin1:flux_rin1,$
		bandfluxes_rin2:flux_rin2,$
		bandfluxes_rin3:flux_rin3,$
		bandfluxes_rin4:flux_rin4,$
		bandfluxes_rin5:flux_rin5,$
		bandfluxes_rin6:flux_rin6,$
		bandfluxes_rin7:flux_rin7,$
		bandfluxes_rin8:flux_rin8,$
		bandfluxes_rin9:flux_rin9,$
		bandfluxes_rin10:flux_rin10,$
		bandfluxes_rin11:flux_rin11,$
		bandfluxes_rin12:flux_rin12 $
	     }

return, ret_struct

END