c f77 get_fluxes_SF_AGN.f -o get_fluxes_SF_AGN c Cannibalized get_fluxes_SF.f ccccccccccccccccc Version 2 cccccccccccccccccccccccccccccccccccccccccccc c Incorporated redshift ccccccccccccccccc Version 3 cccccccccccccccccccccccccccccccccccccccccccc c Computes total 5-1100um TIR energy integer nalpha,npts,nfilters,ntransmax parameter(nalpha=64,npts=1496,nfilters=20,ntransmax=700,nmix=21) parameter(n_all=3096) integer nf(nfilters) real lambda(npts),lognufnu(npts,nalpha),f(nfilters), : logf(nfilters) real lambdaIR(ntransmax,nfilters),transIR(ntransmax,nfilters), : dlambda(nfilters),delta_nu(nfilters) real lambda_all(n_all),trans_all(n_all),fsp1(6,nalpha),fsp(6), : ftot(nmix,nalpha),colorA(nmix,nalpha),colorB(nmix,nalpha) data nf/97,171,149,32,21,17,135,167,129,156,71,31,25,621, : 200,200,200,200,200,200/ data dlambda/0.05,0.05,0.05,0.5,3.,5.,0.01,0.01,0.02,0.03, : 0.2,2.0,4.0,0.5,2.,2.,2.,2.,2.,2./ data delta_nu/23.029e12,13.48e12,7.994e12,5.16e12,2.58e12,1.0e12, : 16.321e12,12.846e12,11.836e12,12.131e12,2.787e12,1.212e12, : 0.435e12,0.0824e12,1.296e12,1.016e12,8.150e11,3.702e11,2.592e11, : 2.290e11/ C These delta nu's correspond to the following bandwidths in um: C 3.500,6.47,6.000,10.76,30.98,33.36,0.69,0.87,1.28,2.57,5.30,21.32, c 35.77,198.70, c herschel dlambda/lambda=0.3086,0.3467,0.4515,0.3098,0.3026,0.3903 c And to the following central wavelengths: c iso 7; iras 12; iso 15; iras 25, 60, 100; sirtf 3.6, 4.5, 5.8, c 8.0, 24, 70, 160; scuba 850 c Then the output is for 20cm, 450u, and 200um c Then the final 6 fluxes are from Herschel: c HSO 71.3, 102.4, 166.1, 250.8, 350.0, 510.4 (from filter FWHMs). character filename(nfilters)*70 character fname(nmix)*32,sname(nmix)*27 character cAGN(nmix)*4,credshift*6 data cAGN/'0.00','0.05','0.10','0.15','0.20','0.25','0.30','0.35', : '0.40','0.45','0.50','0.55','0.60','0.65','0.70','0.75','0.80', : '0.85','0.90','0.95','1.00'/ open(10,file='/crow1/spectra/orig/flat.dat',status='old') do i=1,n_all read(10,*) lambda_all(i),trans_all(i) enddo close(10) open(40,file='Ldistrib/Ldistrib.SF_AGN.dat',status='unknown') open(41,file='figs/colors.fnu70fnu160.dat',status='unknown') open(42,file='figs/colors.fnu60fnu100.dat',status='unknown') redshift=0. credshift='0.0000' 99 print*,'0: run' print*,'11: redshift = ',redshift read*,nans if(nans.eq.11) then print*,'enter in new redshift (twice; 2nd in form V.WXYZ)' read*,redshift,credshift goto 99 endif do i=1,nmix sname(i)(1:16)='spectra/spectra.' sname(i)(17:20)=cAGN(i) sname(i)(21:27)='AGN.dat' fname(i)(1:13)='fluxes/model.' fname(i)(14:17)=cAGN(i) fname(i)(18:22)='AGN.z' fname(i)(23:28)=credshift fname(i)(29:32)='.dat' enddo ! done looping through AGN fractions ! read in the models; find row numbers for <20cm and <450um n20cm=1 n450um=1 n200um=1 nz20cm=1 nz60um=1 nz100um=1 filename(1)= '/crow1/scopes/iras/spectral_response_lw2.dat' filename(2)= '/crow1/scopes/iras/spectral_response_12mu.dat' filename(3)= '/crow1/scopes/iras/spectral_response_lw3.dat' filename(4)= '/crow1/scopes/iras/spectral_response_25mu.dat' filename(5)= '/crow1/scopes/iras/spectral_response_60mu.dat' filename(6)= '/crow1/scopes/iras/spectral_response_100mu.dat' filename(7)= '/crow1/scopes/sirtf/spectral_response_3.6mu.dat' filename(8)= '/crow1/scopes/sirtf/spectral_response_4.5mu.dat' filename(9)= '/crow1/scopes/sirtf/spectral_response_5.8mu.dat' filename(10)='/crow1/scopes/sirtf/spectral_response_8.0mu.dat' filename(11)='/crow1/scopes/sirtf/spectral_response_24mu.dat' filename(12)='/crow1/scopes/sirtf/spectral_response_70mu.dat' filename(13)='/crow1/scopes/sirtf/spectral_response_160mu.dat' filename(14)='/crow1/scopes/submm/spectral_response_850mu.dat' filename(15)='/crow1/scopes/hso/pacs070.txt' filename(16)='/crow1/scopes/hso/pacs100.txt' filename(17)='/crow1/scopes/hso/pacs160.txt' filename(18)='/crow1/scopes/hso/spire250.txt' filename(19)='/crow1/scopes/hso/spire350.txt' filename(20)='/crow1/scopes/hso/spire500.txt' do n=1,nfilters open(10,file=filename(n),status='old') do i=1,nf(n) read(10,*) lambdaIR(i,n),transIR(i,n) enddo close(10) enddo do nn=1,nmix ! loop thru different AGN fractions open(10,file=sname(nn),status='old') do k=1,npts read(10,803) lambda(k),(lognufnu(k,l),l=1,nalpha) 803 format(f10.3,64(1x,f8.4)) lambda(k)=lambda(k)*(1.+redshift) z20cm=200000.*(1.+redshift) z60um=60.*(1.+redshift) z100um=100.*(1.+redshift) if(lambda(k).lt.200000.) n20cm=k if(lambda(k).lt.450.) n450um=k if(lambda(k).lt.200.) n200um=k if(lambda(k).lt.z20cm) nz20cm=k if(lambda(k).lt.z60um) nz60um=k if(lambda(k).lt.z100um) nz100um=k enddo close(10) ! Compute fluxes under the bandpasses. open(10,file=fname(nn),status='unknown') do l=1,nalpha alpha=real(l)/(real(nalpha)/4.) do n=1,nfilters f(n)=0. do i=1,nf(n) do k=1,npts if(lambdaIR(i,n).lt.lambda(k)) then fb=10**lognufnu(k,l) fa=10**lognufnu(k-1,l) dy=fb-fa dx=lambda(k)-lambda(k-1) fmod=fa+(lambdaIR(i,n)-lambda(k-1))*dy/dx goto 200 endif enddo ! done looping thru lambda ! Here I divide by lambda to both convert nu*f_nu to f_lambda 200 f(n)=f(n)+fmod*transIR(i,n)*dlambda(n)/lambdaIR(i,n) ! W/m^2 enddo ! done looping thru filters lambda f(n)=f(n)/(1.+redshift) ! k correct enddo ! done looping thru filters do n=1,nfilters f(n)=f(n)*1e26/delta_nu(n) ! Jy logf(n)=log10(f(n)) enddo ! write results to file; append 20cm,450um,&200um flux densities write(10,300) alpha,(logf(n),n=1,14), : lognufnu(n20cm,l)+(lognufnu(n20cm+1,l)-lognufnu(n20cm,l))* : (200000.-lambda(n20cm))/(lambda(n20cm+1)-lambda(n20cm))- : log10(2.99792e14/200000.)+26., : lognufnu(n450um,l)+(lognufnu(n450um+1,l)-lognufnu(n450um,l))* : (450.-lambda(n450um))/(lambda(n450um+1)-lambda(n450um))- : log10(2.99792e14/450.)+26., : lognufnu(n200um,l)+(lognufnu(n200um+1,l)-lognufnu(n200um,l))* : (200.-lambda(n200um))/(lambda(n200um+1)-lambda(n200um))- : log10(2.99792e14/200.)+26.,(logf(n),n=15,20) 300 format(f7.5,1x,23(f8.4,1x)) colorA(nn,l)=logf(12)-logf(13) ! log(fnu70/fnu160) colorB(nn,l)=logf(5)-logf(6) ! log(fnu60/fnu100) enddo ! done looping thru alpha close(10) ! compute IR energy budget do l=1,nalpha ! do for all alpha values do ii=1,6 fsp(ii)=0. enddo C integrate flux from 5 to 13 microns dlambda_all=0.01 do i=462,1261 do j=1,npts if(lambda_all(i).lt.lambda(j)) then fb=10**lognufnu(j,l) fa=10**lognufnu(j-1,l) dy=fb-fa dx=lambda(j)-lambda(j-1) fj=fa+(lambda_all(i)-lambda(j-1))*dy/dx fsp(2)=fsp(2)+fj*trans_all(i)*dlambda_all/lambda_all(i)!W/m^2 goto 851 endif enddo 851 continue enddo C integrate flux from 13 to 20 microns do i=1262,1565 do j=1,npts if(lambda_all(i).lt.lambda(j)) then if(lambda_all(i).ge.15.7) then dlambda_all=0.1 else dlambda_all=0.01 endif fb=10**lognufnu(j,l) fa=10**lognufnu(j-1,l) dy=fb-fa dx=lambda(j)-lambda(j-1) fj=fa+(lambda_all(i)-lambda(j-1))*dy/dx fsp(3)=fsp(3)+fj*trans_all(i)*dlambda_all/lambda_all(i)!W/m^2 goto 852 endif enddo 852 continue enddo C integrate flux from 20 to 42 microns dlambda_all=0.1 do i=1566,1785 do j=1,npts if(lambda_all(i).lt.lambda(j)) then fb=10**lognufnu(j,l) fa=10**lognufnu(j-1,l) dy=fb-fa dx=lambda(j)-lambda(j-1) fj=fa+(lambda_all(i)-lambda(j-1))*dy/dx fsp(4)=fsp(4)+fj*trans_all(i)*dlambda_all/lambda_all(i)!W/m^2 goto 853 endif enddo 853 continue enddo C integrate flux from 42 to 122 microns do i=1786,2118 do j=1,npts if(lambda_all(i).lt.lambda(j)) then if(lambda_all(i).ge.71.) then dlambda_all=1. else dlambda_all=0.1 endif fb=10**lognufnu(j,l) fa=10**lognufnu(j-1,l) dy=fb-fa dx=lambda(j)-lambda(j-1) fj=fa+(lambda_all(i)-lambda(j-1))*dy/dx fsp(5)=fsp(5)+fj*trans_all(i)*dlambda_all/lambda_all(i)!W/m^2 goto 854 endif enddo 854 continue enddo C integrate flux from 122 to 1100 microns dlambda_all=1. do i=2119,3096 do j=1,npts if(lambda_all(i).lt.lambda(j)) then fb=10**lognufnu(j,l) fa=10**lognufnu(j-1,l) dy=fb-fa dx=lambda(j)-lambda(j-1) fj=fa+(lambda_all(i)-lambda(j-1))*dy/dx fsp(6)=fsp(6)+fj*trans_all(i)*dlambda_all/lambda_all(i) goto 855 endif enddo 855 continue enddo ftot(nn,l)=fsp(2)+fsp(3)+fsp(4)+fsp(5)+fsp(6) ! W/m^2 enddo ! done looping thru alpha enddo ! done looping through AGN fractions do l=1,nalpha alpha=real(l)/(real(nalpha)/4.) write(40,857) alpha,(log10(ftot(nn,l)),nn=1,nmix) 857 format(f6.4,21(1x,f8.4)) write(41,857) alpha,(colorA(nn,l),nn=1,nmix) ! log(fnu70/fnu160) write(42,857) alpha,(colorB(nn,l),nn=1,nmix) ! log(fnu60/fnu10) enddo close(40) close(41) close(42) end ! of main program