diff --git a/examples/README.md b/examples/README.md deleted file mode 100644 index 604040d..0000000 --- a/examples/README.md +++ /dev/null @@ -1,11 +0,0 @@ -# Running on NERSC - -Request a standard Haswell node (32 physical cores, 64 hyperthreads) -``` -salloc -N 1 -C haswell -q interactive -t 03:00:00 -``` - -Run Julia with 64 threads, -``` -julia -t 64 cib_planck2013_chunked.jl -``` diff --git a/examples/benchmark_fits_hdf5.jl b/examples/benchmark_fits_hdf5.jl deleted file mode 100644 index 3549551..0000000 --- a/examples/benchmark_fits_hdf5.jl +++ /dev/null @@ -1,36 +0,0 @@ -m = HealpixMap{Float64,RingOrder}(model.nside) - -using BenchmarkTools -using HDF5 -using JLD2 - -@btime Healpix.saveToFITS(m, "!test.fits", typechar="E") - -## - -function testHDF5() - h5open("test.h5", "w") do file - write(file, "A", m.pixels) # alternatively, say "@write file A" - end -end - -## -@btime testHDF5() - - -## -function testJLD2() - A = m.pixels - @save "test.jld2" A -end - - -## -@time testJLD2() -## -@time testHDF5() -## -@time Healpix.saveToFITS(m, "!test.fits", typechar="E") -## -@time Healpix.saveToFITS(m, "!test.fits", typechar="D") -## diff --git a/examples/cib_co_chunked.jl b/examples/cib_co_chunked.jl deleted file mode 100644 index 01a501e..0000000 --- a/examples/cib_co_chunked.jl +++ /dev/null @@ -1,208 +0,0 @@ -using XGPaint -using SparseArrays -using Healpix -using HDF5 -using JLD2, FileIO, CodecZlib -using DelimitedFiles -using Random -Random.seed!(3) - -halo_pos, halo_mass = read_halo_catalog_hdf5("/fs/lustre/cita/zack/projects/websky/websky_halos-light.hdf5") -cosmo = get_cosmology(h=0.677f0, OmegaM=0.310f0) -modelCIB = CIB_Scarfy{Float32}() -model_CO = CO_CROWNED{Float32}() -freqs = [ - "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4", - "63.9", "67.8", "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143", - "153", "164", "189", "210", "217", "232", "256", "275", "294", "306", "314", - "340", "353", "375", "409", "467", "525", "545", "584", "643", "729", "817", - "857", "906", "994", "1080" -] - -bandpass_PA6_090 = open(readdlm,"/home/dongwooc/lustrespace/tilec/data/PA6_avg_passband_f090_wErr_trunc_20200505.txt") -bandpass_PA6_150 = open(readdlm,"/home/dongwooc/lustrespace/tilec/data/PA6_avg_passband_f150_wErr_trunc_20200505.txt") -bandpass_PA4_220 = open(readdlm,"/home/dongwooc/lustrespace/tilec/data/PA4_avg_passband_f220_wErr_trunc_20200505.txt") - -function write_chunk( - output_dir, chunk_index, modelCIB, model_CO, cosmo, - pos, mass, freqs) - sourcesCIB = generate_sources(modelCIB, cosmo, pos, mass); - sources_CO = process_sources(model_CO, sourcesCIB, modelCIB); - fluxes_cen = Array{Float32,1}(undef,sourcesCIB.N_cen) - fluxes_sat = Array{Float32,1}(undef,sourcesCIB.N_sat) - println("writing info for ",sourcesCIB.N_cen," centrals") - h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index).h5"), "w") do file - write(file, "redshift", sourcesCIB.redshift_cen) - write(file, "theta", sourcesCIB.theta_cen) - write(file, "phi", sourcesCIB.phi_cen) - write(file, "LIR", sources_CO.LIR_cen) - write(file, "LCO", sources_CO.LcoJ_cen) - write(file, "LCI", sources_CO.LCI_cen) - end - println("writing info for ",sourcesCIB.N_sat," satellites") - h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index).h5"), "w") do file - write(file, "redshift", sourcesCIB.redshift_sat) - write(file, "theta", sourcesCIB.theta_sat) - write(file, "phi", sourcesCIB.phi_sat) - write(file, "LIR", sources_CO.LIR_sat) - write(file, "LCO", sources_CO.LcoJ_sat) - write(file, "LCI", sources_CO.LCI_sat) - end - m = HealpixMap{Float64,RingOrder}(modelCIB.nside) - mCO090 = HealpixMap{Float64,RingOrder}(model_CO.nside) - mCO150 = HealpixMap{Float64,RingOrder}(model_CO.nside) - mCO220 = HealpixMap{Float64,RingOrder}(model_CO.nside) - mCI090 = HealpixMap{Float64,RingOrder}(model_CO.nside) - mCI150 = HealpixMap{Float64,RingOrder}(model_CO.nside) - mCI220 = HealpixMap{Float64,RingOrder}(model_CO.nside) - println("*** 90 GHz ***") - fluxCO_090_cen, fluxCI_090_cen, fluxCO_090_sat, fluxCI_090_sat = XGPaint.paint!(mCO090, mCI090, bandpass_PA6_090, model_CO, sources_CO) - println("*** 150 GHz ***") - fluxCO_150_cen, fluxCI_150_cen, fluxCO_150_sat, fluxCI_150_sat = XGPaint.paint!(mCO150, mCI150, bandpass_PA6_150, model_CO, sources_CO) - println("*** 220 GHz ***") - fluxCO_220_cen, fluxCI_220_cen, fluxCO_220_sat, fluxCI_220_sat = XGPaint.paint!(mCO220, mCI220, bandpass_PA4_220, model_CO, sources_CO) - println("writing CO maps ...") - filename = joinpath(output_dir, "co_090.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCO090.pixels = mCO090.pixels + m0.pixels - end - Healpix.saveToFITS(mCO090, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_150.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCO150.pixels = mCO150.pixels + m0.pixels - end - Healpix.saveToFITS(mCO150, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_220.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCO220.pixels = mCO220.pixels + m0.pixels - end - Healpix.saveToFITS(mCO220, "!$(filename)", typechar="D") - println("writing CI maps ...") - filename = joinpath(output_dir, "ci_090.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCI090.pixels = mCI090.pixels + m0.pixels - end - Healpix.saveToFITS(mCI090, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "ci_150.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCI150.pixels = mCI150.pixels + m0.pixels - end - Healpix.saveToFITS(mCI150, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "ci_220.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCI220.pixels = mCI220.pixels + m0.pixels - end - Healpix.saveToFITS(mCI220, "!$(filename)", typechar="D") - println("writing CO fluxes ...") - for J in 1:7 - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCO$(J)_090.jld2"), "w") do file - write(file, "flux", sparse(fluxCO_090_cen[:,J])) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCO$(J)_150.jld2"), "w") do file - write(file, "flux", sparse(fluxCO_150_cen[:,J])) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCO$(J)_220.jld2"), "w") do file - write(file, "flux", sparse(fluxCO_220_cen[:,J])) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCO$(J)_090.jld2"), "w") do file - write(file, "flux", sparse(fluxCO_090_sat[:,J])) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCO$(J)_150.jld2"), "w") do file - write(file, "flux", sparse(fluxCO_150_sat[:,J])) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCO$(J)_220.jld2"), "w") do file - write(file, "flux", sparse(fluxCO_220_sat[:,J])) - end - end - println("writing CI fluxes ...") - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCI_090.jld2"), "w") do file - write(file, "flux", sparse(fluxCI_090_cen)) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCI_150.jld2"), "w") do file - write(file, "flux", sparse(fluxCI_150_cen)) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCI_220.jld2"), "w") do file - write(file, "flux", sparse(fluxCI_220_cen)) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCI_090.jld2"), "w") do file - write(file, "flux", sparse(fluxCI_090_sat)) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCI_150.jld2"), "w") do file - write(file, "flux", sparse(fluxCI_150_sat)) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCI_220.jld2"), "w") do file - write(file, "flux", sparse(fluxCI_220_sat)) - end - for freq in freqs - - XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, modelCIB, sourcesCIB, - fluxes_cen, fluxes_sat) - # save fluxes - h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - write(file, "flux", fluxes_cen) - end - h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - write(file, "flux", fluxes_sat) - end - # save maps - filename = joinpath(output_dir, "cib_$(freq).fits") - - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - m.pixels = m.pixels + m0.pixels - end - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - end -end - -function run_all_chunks(output_dir, halo_pos, halo_mass, freqs; N_chunks=4) - # provide views into halo positions and masses for chunks of the halos - N_halos = size(halo_mass, 1) - chunksize = trunc(Integer, N_halos / N_chunks + 1) - chunks = chunk(N_halos, chunksize) - for chunk_index in 1:length(chunks) - left_ind, right_ind = chunks[chunk_index] - println("Chunk ", chunk_index, "/", length(chunks), - " ", left_ind, " ", right_ind) - pos = @view halo_pos[:, left_ind:right_ind] - mass = @view halo_mass[left_ind:right_ind] - write_chunk(output_dir, chunk_index, modelCIB, model_CO, cosmo, - pos, mass, freqs) - end - # save maps - println("writing CO+CI maps ...") - m = HealpixMap{Float64,RingOrder}(model_CO.nside) - filename = joinpath(output_dir, "co_ci_090.fits") - m0 = Healpix.readMapFromFITS(joinpath(output_dir, "co_090.fits"), - 1, Float32) - m1 = Healpix.readMapFromFITS(joinpath(output_dir, "ci_090.fits"), - 1, Float32) - m.pixels = m0.pixels + m1.pixels - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_ci_150.fits") - m0 = Healpix.readMapFromFITS(joinpath(output_dir, "co_150.fits"), - 1, Float32) - m1 = Healpix.readMapFromFITS(joinpath(output_dir, "ci_150.fits"), - 1, Float32) - m.pixels = m0.pixels + m1.pixels - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_ci_220.fits") - m0 = Healpix.readMapFromFITS(joinpath(output_dir, "co_220.fits"), - 1, Float32) - m1 = Healpix.readMapFromFITS(joinpath(output_dir, "ci_220.fits"), - 1, Float32) - m.pixels = m0.pixels + m1.pixels - Healpix.saveToFITS(m, "!$(filename)", typechar="D") -end -## compute on all chunks, on all halos -scratch_dir = "/home/dongwooc/projectscratchspace/cib_co_sources_scarfy_refactor" -println("SCRATCH: ", scratch_dir) -mkpath(scratch_dir) -mkpath(joinpath(scratch_dir, "sources")) -run_all_chunks(scratch_dir, halo_pos, halo_mass, freqs) diff --git a/examples/cib_co_chunked_legacy.jl b/examples/cib_co_chunked_legacy.jl deleted file mode 100644 index 54862fc..0000000 --- a/examples/cib_co_chunked_legacy.jl +++ /dev/null @@ -1,378 +0,0 @@ -using XGPaint -using Healpix -using HDF5 -using JLD2, FileIO, CodecZlib -using DelimitedFiles -using Random -Random.seed!(3) - -halo_pos, halo_mass = read_halo_catalog_hdf5("/fs/lustre/cita/zack/projects/websky/websky_halos-light.hdf5") -cosmo = get_cosmology(h=0.677f0, OmegaM=0.310f0) -model = CIB_Planck2013{Float32}(z_evo="scarfy",m_evo="scarfy") -println(model.z_evo) -println(model.m_evo) -R_CO10_CI = 0.18*77.83 # r=0.18 times cubic frequency scaling -xRJ_GHz = 0.0176086 - -freqs = [ - "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4", - "63.9", "67.8", "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143", - "153", "164", "189", "210", "217", "232", "256", "275", "294", "306", "314", - "340", "353", "375", "409", "467", "525", "545", "584", "643", "729", "817", - "857", "906", "994", "1080" -] - -bandpass_PA6_090 = open(readdlm,"/home/dongwooc/lustrespace/tilec/data/PA6_avg_passband_f090_wErr_trunc_20200505.txt") -bandpass_PA6_150 = open(readdlm,"/home/dongwooc/lustrespace/tilec/data/PA6_avg_passband_f150_wErr_trunc_20200505.txt") -bandpass_PA4_220 = open(readdlm,"/home/dongwooc/lustrespace/tilec/data/PA4_avg_passband_f220_wErr_trunc_20200505.txt") - -bandpass_PA6_090_edges = [bandpass_PA6_090[1,1].-diff(bandpass_PA6_090[1:2,1])/2;bandpass_PA6_090[1:end-1,1]+diff(bandpass_PA6_090[:,1])/2;bandpass_PA6_090[end,1].+diff(bandpass_PA6_090[end-1:end,1])/2] -bandpass_PA6_150_edges = [bandpass_PA6_150[1,1].-diff(bandpass_PA6_150[1:2,1])/2;bandpass_PA6_150[1:end-1,1]+diff(bandpass_PA6_150[:,1])/2;bandpass_PA6_150[end,1].+diff(bandpass_PA6_150[end-1:end,1])/2] -bandpass_PA4_220_edges = [bandpass_PA4_220[1,1].-diff(bandpass_PA4_220[1:2,1])/2;bandpass_PA4_220[1:end-1,1]+diff(bandpass_PA4_220[:,1])/2;bandpass_PA4_220[end,1].+diff(bandpass_PA4_220[end-1:end,1])/2] - -bandpass_PA6_090_norm = bandpass_PA6_090[:,2]/sum(bandpass_PA6_090[:,2]) -bandpass_PA6_150_norm = bandpass_PA6_150[:,2]/sum(bandpass_PA6_150[:,2]) -bandpass_PA4_220_norm = bandpass_PA4_220[:,2]/sum(bandpass_PA4_220[:,2]) - -bandpass_PA6_090_diff = diff(bandpass_PA6_090_edges) -bandpass_PA6_150_diff = diff(bandpass_PA6_150_edges) -bandpass_PA4_220_diff = diff(bandpass_PA4_220_edges) - -bandpass_PA6_090_len = length(bandpass_PA6_090_norm) -bandpass_PA6_150_len = length(bandpass_PA6_150_norm) -bandpass_PA4_220_len = length(bandpass_PA4_220_norm) - -function write_chunk( - output_dir, chunk_index, model, cosmo, - pos, mass, freqs) - sources = generate_sources(model, cosmo, pos, mass); - fluxes_cen = Array{Float32,1}(undef,sources.N_cen) - fluxes_sat = Array{Float32,1}(undef,sources.N_sat) - LIR_cen = Array{Float32,1}(undef,sources.N_cen) - LIR_sat = Array{Float32,1}(undef,sources.N_sat) - LcoJ_cen = Array{Float32,2}(undef,sources.N_cen,7) - LcoJ_sat = Array{Float32,2}(undef,sources.N_sat,7) - Lnu_to_LIR(Td) = 9.521f8*(53.865 + 3.086*Td - 0.213*Td^2 + 0.00749*Td^3); - # calculate observing frequencies and quasi-temperatures for each halo - nuJ_cen = Array{Float32,2}(undef,(sources.N_cen,7)) - nuJ_sat = Array{Float32,2}(undef,(sources.N_sat,7)) - quasiTcoJ_cen = Array{Float32,2}(undef,(sources.N_cen,7)) - quasiTcoJ_sat = Array{Float32,2}(undef,(sources.N_sat,7)) - fluxCO_090_cen = Array{Float32,2}(undef,(sources.N_cen,7)) - fluxCO_090_sat = Array{Float32,2}(undef,(sources.N_sat,7)) - fill!(fluxCO_090_cen,zero(Float32)) - fill!(fluxCO_090_sat,zero(Float32)) - fluxCO_150_cen = Array{Float32,2}(undef,(sources.N_cen,7)) - fluxCO_150_sat = Array{Float32,2}(undef,(sources.N_sat,7)) - fill!(fluxCO_150_cen,zero(Float32)) - fill!(fluxCO_150_sat,zero(Float32)) - fluxCO_220_cen = Array{Float32,2}(undef,(sources.N_cen,7)) - fluxCO_220_sat = Array{Float32,2}(undef,(sources.N_sat,7)) - fill!(fluxCO_220_cen,zero(Float32)) - fill!(fluxCO_220_sat,zero(Float32)) - nuCI_cen = Array{Float32,1}(undef,sources.N_cen) - nuCI_sat = Array{Float32,1}(undef,sources.N_sat) - LCI_cen = Array{Float32,1}(undef,sources.N_cen) - LCI_sat = Array{Float32,1}(undef,sources.N_sat) - quasiTCI_cen = Array{Float32,1}(undef,sources.N_cen) - quasiTCI_sat = Array{Float32,1}(undef,sources.N_sat) - fluxCI_090_cen = Array{Float32,1}(undef,sources.N_cen) - fill!(fluxCI_090_cen,zero(Float32)) - fluxCI_150_cen = Array{Float32,1}(undef,sources.N_cen) - fill!(fluxCI_150_cen,zero(Float32)) - fluxCI_220_cen = Array{Float32,1}(undef,sources.N_cen) - fill!(fluxCI_220_cen,zero(Float32)) - fluxCI_090_sat = Array{Float32,1}(undef,sources.N_sat) - fill!(fluxCI_090_sat,zero(Float32)) - fluxCI_150_sat = Array{Float32,1}(undef,sources.N_sat) - fill!(fluxCI_150_sat,zero(Float32)) - fluxCI_220_sat = Array{Float32,1}(undef,sources.N_sat) - fill!(fluxCI_220_sat,zero(Float32)) - Threads.@threads :static for i in 1:sources.N_cen - # calculate dust temperatures and bolometric LIR for each source - Td = model.shang_Td * (1.f0 .+sources.redshift_cen[i])^model.shang_alpha; - LIR_cen[i] = Lnu_to_LIR(Td)*sources.lum_cen[i]*XGPaint.nu2theta(2.10833f12,sources.redshift_cen[i],model) - rnd = Float32(randn()) - r21 = 0.75f0+0.11f0*(rnd/2.0f0+373.f0/275.f0-sqrt(rnd^2.0f0/4.0f0-252.0f0/275.0f0*rnd+139129.0f0/75625.0f0)) - # generate CO SLED for each source with log-normal LCO(1-0) scatter - LcoJ_cen[i,1] = LIR_cen[i]^(1/1.37)*10^(1.74/1.37)*exp.((randn().-0.5*2.302585*0.3)*2.302585*0.3)*4.9e-5 - for J in 2:7 - LcoJ_cen[i,J] = LcoJ_cen[i,1]/(1.0f0+exp((log.(1.0f0/r21-1.0f0)-1.0f0)+Float32(J-1)))*J^3 - end - for J in 1:7 - nuJ_cen[i,J] = 115.27f0*J/(1.0f0+sources.redshift_cen[i]) - quasiTcoJ_cen[i,J] = LcoJ_cen[i,J]*(1.315415f0/4.0f0/pi/nuJ_cen[i,J]^2.0f0/sources.dist_cen[i]^2.0f0/(1.0f0+sources.redshift_cen[i])^2.0f0) # divide by dnu in GHz - xRJ = xRJ_GHz*nuJ_cen[i,J] - quasiTcoJ_cen[i,J]/= xRJ^2*exp(xRJ)/(exp(xRJ)-1)^2 - end - nuCI_cen[i] = 492.16f0/(1.0f0+sources.redshift_cen[i]) - LCI_cen[i] = LcoJ_cen[i,1]*R_CO10_CI*exp((randn()-0.5*2.302585*0.2)*2.302585*0.2) - quasiTCI_cen[i] = LCI_cen[i]*(1.315415f0/4.0f0/pi/nuCI_cen[i]^2.0f0/sources.dist_cen[i]^2.0f0/(1.0f0+sources.redshift_cen[i])^2.0f0) # divide by dnu in GHz - xRJ = xRJ_GHz*nuCI_cen[i] - quasiTCI_cen[i]/= xRJ^2*exp(xRJ)/(exp(xRJ)-1)^2 - end - Threads.@threads :static for i in 1:sources.N_sat - Td = model.shang_Td * (1.f0 .+sources.redshift_sat[i])^model.shang_alpha; - LIR_sat[i] = Lnu_to_LIR(Td)*sources.lum_sat[i]*XGPaint.nu2theta(2.10833f12,sources.redshift_sat[i],model) - rnd = Float32(randn()) - r21 = 0.75f0+0.11f0*(rnd/2.0f0+373.f0/275.f0-sqrt(rnd^2.0f0/4.0f0-252.0f0/275.0f0*rnd+139129.0f0/75625.0f0)) - LcoJ_sat[i,1] = LIR_sat[i]^(1/1.37)*10^(1.74/1.37)*exp((randn()-0.5*2.302585*0.3)*2.302585*0.3)*4.9e-5 - for J in 2:7 - LcoJ_sat[i,J] = LcoJ_sat[i,1]/(1.0f0+exp((log.(1.0f0/r21-1.0f0)-1.0f0)+Float32(J-1)))*J^3 - end - for J in 1:7 - nuJ_sat[i,J] = 115.27f0*J/(1.0f0+sources.redshift_sat[i]) - quasiTcoJ_sat[i,J] = LcoJ_sat[i,J]*(1.315415f0/4.0f0/pi/nuJ_sat[i,J]^2.0f0/sources.dist_sat[i]^2.0f0/(1.0f0+sources.redshift_sat[i])^2.0f0) # divide by dnu in GHz - xRJ = xRJ_GHz*nuJ_sat[i,J] - quasiTcoJ_sat[i,J]/= xRJ^2*exp(xRJ)/(exp(xRJ)-1)^2 - end - nuCI_sat[i] = 492.16f0/(1.0f0+sources.redshift_sat[i]) - LCI_sat[i] = LcoJ_sat[i,1]*R_CO10_CI*exp((randn()-0.5*2.302585*0.2)*2.302585*0.2) - quasiTCI_sat[i] = LCI_sat[i]*(1.315415f0/4.0f0/pi/nuCI_sat[i]^2.0f0/sources.dist_sat[i]^2.0f0/(1.0f0+sources.redshift_sat[i])^2.0f0) # divide by dnu in GHz - xRJ = xRJ_GHz*nuCI_sat[i] - quasiTCI_sat[i]/= xRJ^2*exp(xRJ)/(exp(xRJ)-1)^2 - end - println("writing info for ",sources.N_cen," centrals") - h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index).h5"), "w") do file - write(file, "redshift", sources.redshift_cen) - write(file, "theta", sources.theta_cen) - write(file, "phi", sources.phi_cen) - write(file, "LIR", LIR_cen) - write(file, "LCO", LcoJ_cen) - write(file, "LCI", LCI_cen) - end - println("writing info for ",sources.N_sat," satellites") - h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index).h5"), "w") do file - write(file, "redshift", sources.redshift_sat) - write(file, "theta", sources.theta_sat) - write(file, "phi", sources.phi_sat) - write(file, "LIR", LIR_sat) - write(file, "LCO", LcoJ_sat) - write(file, "LCI", LCI_sat) - end - m = HealpixMap{Float64,RingOrder}(model.nside) - mCO090 = HealpixMap{Float64,RingOrder}(model.nside) - mCO150 = HealpixMap{Float64,RingOrder}(model.nside) - mCO220 = HealpixMap{Float64,RingOrder}(model.nside) - mCO090_pixsr = nside2pixarea(mCO090.resolution.nside) - mCO150_pixsr = nside2pixarea(mCO150.resolution.nside) - mCO220_pixsr = nside2pixarea(mCO220.resolution.nside) - mCI090 = HealpixMap{Float64,RingOrder}(model.nside) - mCI150 = HealpixMap{Float64,RingOrder}(model.nside) - mCI220 = HealpixMap{Float64,RingOrder}(model.nside) - println("painting CO/CI from ",sources.N_cen," centrals") - Threads.@threads :static for i in 1:sources.N_cen - for J in 1:7 - j = searchsortedlast(bandpass_PA6_090_edges,nuJ_cen[i,J]) - if (j>0) && (j<=bandpass_PA6_090_len) - fluxCO_090_cen[i,J] = quasiTcoJ_cen[i,J]*bandpass_PA6_090_norm[j]/bandpass_PA6_090_diff[j]/mCO090_pixsr; - mCO090[sources.hp_ind_cen[i]]+= fluxCO_090_cen[i,J]; - end - j = searchsortedlast(bandpass_PA6_150_edges,nuJ_cen[i,J]) - if (j>0) && (j<=bandpass_PA6_150_len) - fluxCO_150_cen[i,J] = quasiTcoJ_cen[i,J]*bandpass_PA6_150_norm[j]/bandpass_PA6_150_diff[j]/mCO150_pixsr; - mCO150[sources.hp_ind_cen[i]]+= fluxCO_150_cen[i,J]; - end - j = searchsortedlast(bandpass_PA4_220_edges,nuJ_cen[i,J]) - if (j>0) && (j<=bandpass_PA4_220_len) - fluxCO_220_cen[i,J] = quasiTcoJ_cen[i,J]*bandpass_PA4_220_norm[j]/bandpass_PA4_220_diff[j]/mCO220_pixsr; - mCO220[sources.hp_ind_cen[i]]+= fluxCO_220_cen[i,J]; - end - end - j = searchsortedlast(bandpass_PA6_090_edges,nuCI_cen[i]) - if (j>0) && (j<=bandpass_PA6_090_len) - fluxCI_090_cen[i] = quasiTCI_cen[i]*bandpass_PA6_090_norm[j]/bandpass_PA6_090_diff[j]/mCO090_pixsr; - mCI090[sources.hp_ind_cen[i]]+= fluxCI_090_cen[i]; - end - j = searchsortedlast(bandpass_PA6_150_edges,nuCI_cen[i]) - if (j>0) && (j<=bandpass_PA6_150_len) - fluxCI_150_cen[i] = quasiTCI_cen[i]*bandpass_PA6_150_norm[j]/bandpass_PA6_150_diff[j]/mCO150_pixsr; - mCI150[sources.hp_ind_cen[i]]+= fluxCI_150_cen[i]; - end - j = searchsortedlast(bandpass_PA4_220_edges,nuCI_cen[i]) - if (j>0) && (j<=bandpass_PA4_220_len) - fluxCI_220_cen[i] = quasiTCI_cen[i]*bandpass_PA4_220_norm[j]/bandpass_PA4_220_diff[j]/mCO220_pixsr; - mCI220[sources.hp_ind_cen[i]]+= fluxCI_220_cen[i]; - end - end - println("painting CO/CI from ",sources.N_sat," satellites") - Threads.@threads :static for i in 1:sources.N_sat - for J in 1:7 - j = searchsortedlast(bandpass_PA6_090_edges,nuJ_sat[i,J]) - if (j>0) && (j<=bandpass_PA6_090_len) - fluxCO_090_sat[i,J] = quasiTcoJ_sat[i,J]*bandpass_PA6_090_norm[j]/bandpass_PA6_090_diff[j]/mCO090_pixsr; - mCO090[sources.hp_ind_sat[i]]+= fluxCO_090_sat[i,J]; - end - j = searchsortedlast(bandpass_PA6_150_edges,nuJ_sat[i,J]) - if (j>0) && (j<=bandpass_PA6_150_len) - fluxCO_150_sat[i,J] = quasiTcoJ_sat[i,J]*bandpass_PA6_150_norm[j]/bandpass_PA6_150_diff[j]/mCO150_pixsr; - mCO150[sources.hp_ind_sat[i]]+= fluxCO_150_sat[i,J]; - end - j = searchsortedlast(bandpass_PA4_220_edges,nuJ_sat[i,J]) - if (j>0) && (j<=bandpass_PA4_220_len) - fluxCO_220_sat[i,J] = quasiTcoJ_sat[i,J]*bandpass_PA4_220_norm[j]/bandpass_PA4_220_diff[j]/mCO220_pixsr; - mCO220[sources.hp_ind_sat[i]]+= fluxCO_220_sat[i,J]; - end - end - j = searchsortedlast(bandpass_PA6_090_edges,nuCI_sat[i]) - if (j>0) && (j<=bandpass_PA6_090_len) - fluxCI_090_sat[i] = quasiTCI_sat[i]*bandpass_PA6_090_norm[j]/bandpass_PA6_090_diff[j]/mCO090_pixsr; - mCI090[sources.hp_ind_sat[i]]+= fluxCI_090_sat[i]; - end - j = searchsortedlast(bandpass_PA6_150_edges,nuCI_sat[i]) - if (j>0) && (j<=bandpass_PA6_150_len) - fluxCI_150_sat[i] = quasiTCI_sat[i]*bandpass_PA6_150_norm[j]/bandpass_PA6_150_diff[j]/mCO150_pixsr; - mCI150[sources.hp_ind_sat[i]]+= fluxCI_150_sat[i]; - end - j = searchsortedlast(bandpass_PA4_220_edges,nuCI_sat[i]) - if (j>0) && (j<=bandpass_PA4_220_len) - fluxCI_220_sat[i] = quasiTCI_sat[i]*bandpass_PA4_220_norm[j]/bandpass_PA4_220_diff[j]/mCO220_pixsr; - mCI220[sources.hp_ind_sat[i]]+= fluxCI_220_sat[i]; - end - end - println("writing CO maps ...") - filename = joinpath(output_dir, "co_090.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCO090.pixels = mCO090.pixels + m0.pixels - end - Healpix.saveToFITS(mCO090, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_150.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCO150.pixels = mCO150.pixels + m0.pixels - end - Healpix.saveToFITS(mCO150, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_220.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCO220.pixels = mCO220.pixels + m0.pixels - end - Healpix.saveToFITS(mCO220, "!$(filename)", typechar="D") - println("writing CI maps ...") - filename = joinpath(output_dir, "ci_090.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCI090.pixels = mCI090.pixels + m0.pixels - end - Healpix.saveToFITS(mCI090, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "ci_150.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCI150.pixels = mCI150.pixels + m0.pixels - end - Healpix.saveToFITS(mCI150, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "ci_220.fits") - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - mCI220.pixels = mCI220.pixels + m0.pixels - end - Healpix.saveToFITS(mCI220, "!$(filename)", typechar="D") - println("writing CO fluxes ...") - for J in 1:7 - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCO$(J)_090.jld2"), "w") do file - write(file, "flux", fluxCO_090_cen[:,J]; compress=true) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCO$(J)_150.jld2"), "w") do file - write(file, "flux", fluxCO_150_cen[:,J]; compress=true) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCO$(J)_220.jld2"), "w") do file - write(file, "flux", fluxCO_220_cen[:,J]; compress=true) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCO$(J)_090.jld2"), "w") do file - write(file, "flux", fluxCO_090_sat[:,J]; compress=true) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCO$(J)_150.jld2"), "w") do file - write(file, "flux", fluxCO_150_sat[:,J]; compress=true) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCO$(J)_220.jld2"), "w") do file - write(file, "flux", fluxCO_220_sat[:,J]; compress=true) - end - end - println("writing CI fluxes ...") - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCI_090.jld2"), "w") do file - write(file, "flux", fluxCI_090_cen; compress=true) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCI_150.jld2"), "w") do file - write(file, "flux", fluxCI_150_cen; compress=true) - end - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_fluxCI_220.jld2"), "w") do file - write(file, "flux", fluxCI_220_cen; compress=true) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCI_090.jld2"), "w") do file - write(file, "flux", fluxCI_090_sat; compress=true) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCI_150.jld2"), "w") do file - write(file, "flux", fluxCI_150_sat; compress=true) - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_fluxCI_220.jld2"), "w") do file - write(file, "flux", fluxCI_220_sat; compress=true) - end - for freq in freqs - - XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources, - fluxes_cen, fluxes_sat) - # save fluxes - h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - write(file, "flux", fluxes_cen) - end - h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - write(file, "flux", fluxes_sat) - end - # save maps - filename = joinpath(output_dir, "cib_$(freq).fits") - - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - m.pixels = m.pixels + m0.pixels - end - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - end -end - -## the rest of this is more or less the same as cib_planck2013_chunked.jl -function run_all_chunks(output_dir, halo_pos, halo_mass, freqs; N_chunks=4) - # provide views into halo positions and masses for chunks of the halos - N_halos = size(halo_mass, 1) - chunksize = trunc(Integer, N_halos / N_chunks + 1) - chunks = chunk(N_halos, chunksize) - for chunk_index in 1:length(chunks) - left_ind, right_ind = chunks[chunk_index] - println("Chunk ", chunk_index, "/", length(chunks), - " ", left_ind, " ", right_ind) - pos = @view halo_pos[:, left_ind:right_ind] - mass = @view halo_mass[left_ind:right_ind] - write_chunk(output_dir, chunk_index, model, cosmo, - pos, mass, freqs) - end - # save maps - println("writing CO+CI maps ...") - m = HealpixMap{Float64,RingOrder}(model.nside) - filename = joinpath(output_dir, "co_ci_090.fits") - m0 = Healpix.readMapFromFITS(joinpath(output_dir, "co_090.fits"), - 1, Float32) - m1 = Healpix.readMapFromFITS(joinpath(output_dir, "ci_090.fits"), - 1, Float32) - m.pixels = m0.pixels + m1.pixels - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_ci_150.fits") - m0 = Healpix.readMapFromFITS(joinpath(output_dir, "co_150.fits"), - 1, Float32) - m1 = Healpix.readMapFromFITS(joinpath(output_dir, "ci_150.fits"), - 1, Float32) - m.pixels = m0.pixels + m1.pixels - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - filename = joinpath(output_dir, "co_ci_220.fits") - m0 = Healpix.readMapFromFITS(joinpath(output_dir, "co_220.fits"), - 1, Float32) - m1 = Healpix.readMapFromFITS(joinpath(output_dir, "ci_220.fits"), - 1, Float32) - m.pixels = m0.pixels + m1.pixels - Healpix.saveToFITS(m, "!$(filename)", typechar="D") -end -## compute on all chunks, on all halos -scratch_dir = "/home/dongwooc/scratchspace/cib_co_sources_scarfy" -println("SCRATCH: ", scratch_dir) -mkpath(scratch_dir) -mkpath(joinpath(scratch_dir, "sources")) -run_all_chunks(scratch_dir, halo_pos, halo_mass, freqs) diff --git a/examples/cib_nersc_planck.jl b/examples/cib_nersc_planck.jl deleted file mode 100644 index 36b80a2..0000000 --- a/examples/cib_nersc_planck.jl +++ /dev/null @@ -1,37 +0,0 @@ -using XGPaint -using Healpix -using Random -Random.seed!(3) - -## Load halos from HDF5 files, establish a CIB model and cosmology -halo_pos, halo_mass = read_halo_catalog_hdf5( - ENV["SCRATCH"] * "/websky_halos-light.hdf5") -# @load "/home/zequnl/testsamp.jld2" halo_mass halo_pos - -cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0) -model = CIB_Planck2013{Float32}() - -## Allocate some arrays and file them up for centrals and satellites -@time sources = generate_sources(model, cosmo, halo_pos, halo_mass); - -## Deposit the sources into maps -fluxes_cen = Array{Float32, 1}(undef, sources.N_cen) -fluxes_sat = Array{Float32, 1}(undef, sources.N_sat) -m = HealpixMap{Float64,RingOrder}(model.nside) - -for freq in [ - "18.7", "21.6", - "24.5", "27.3", "30.0", "35.9", - "41.7", "44.0", "47.4", "63.9", "67.8", - "70.0", "73.7", "79.6", "90.2", "100", "111", - "129", "143", "153", "164", "189", "210", "217", - "232", "256", "275", "294", "306", "314", "340", - "353", "375", "409", "467", "525", "545", "584", - "643", "729", "817", "857", "906", "994", "1080"] - - @time begin - XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources, - fluxes_cen, fluxes_sat) - Healpix.saveToFITS(m, "!/global/cscratch1/sd/xzackli/cib/cib$(freq).fits") - end -end diff --git a/examples/cib_planck2013.jl b/examples/cib_planck2013.jl deleted file mode 100644 index 73873fe..0000000 --- a/examples/cib_planck2013.jl +++ /dev/null @@ -1,29 +0,0 @@ -using XGPaint -using Healpix -using Random -Random.seed!(3) - -## Load halos from HDF5 files, establish a CIB model and cosmology -halo_pos, halo_mass = read_halo_catalog_hdf5( - "/tigress/zequnl/xgpaint/websky_halos-light.hdf5") -# @load "/home/zequnl/testsamp.jld2" halo_mass halo_pos - -cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0) -model = CIB_Planck2013{Float32}() - -## Allocate some arrays and file them up for centrals and satellites -@time sources = generate_sources(model, cosmo, halo_pos, halo_mass); - -## Deposit the sources into maps -fluxes_cen = Array{Float32, 1}(undef, sources.N_cen) -fluxes_sat = Array{Float32, 1}(undef, sources.N_sat) -m = HealpixMap{Float64,RingOrder}(model.nside) - -for freq in ["030", "090", "148", "219", "277", "350", - "143", "217", "353", "545", "857"] - @time begin - XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources, - fluxes_cen, fluxes_sat) - Healpix.saveToFITS(m, "!/tigress/zequnl/xgpaint/jl/cib/cib$(freq).fits") - end -end diff --git a/examples/cib_planck2013_chunked.jl b/examples/cib_planck2013_chunked.jl deleted file mode 100644 index 5db0ecf..0000000 --- a/examples/cib_planck2013_chunked.jl +++ /dev/null @@ -1,97 +0,0 @@ -using XGPaint -using Healpix -using HDF5 -using Random -Random.seed!(3) - -## Load halos from HDF5 files, establish a CIB model and cosmology -@time halo_pos, halo_mass = read_halo_catalog_hdf5( - "/global/cfs/cdirs/sobs/www/users/Radio_WebSky/websky_halos-light.hdf5"); -cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0) -model = CIB_Planck2013{Float32}(nside=4096) - -## Write one chunk to disk - -function write_chunk( - output_dir, chunk_index, model, cosmo, - pos, mass, freqs) - # Allocate some arrays and fill them up for centrals and satellites - @time sources = generate_sources(model, cosmo, pos, mass); - - # Deposit the sources into maps - fluxes_cen = Array{Float32, 1}(undef, sources.N_cen) - fluxes_sat = Array{Float32, 1}(undef, sources.N_sat) - m = HealpixMap{Float64,RingOrder}(model.nside) - - h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index).h5"), "w") do file - write(file, "redshift", sources.redshift_cen) - write(file, "theta", sources.theta_cen) - write(file, "phi", sources.phi_cen) - end - h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index).h5"), "w") do file - write(file, "redshift", sources.redshift_sat) - write(file, "theta", sources.theta_sat) - write(file, "phi", sources.phi_sat) - end - - # loop over all frequencies and paint sources to appropriate freq map - @time begin - for freq in freqs - - XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources, - fluxes_cen, fluxes_sat) - - # save sources with mass, redshift, angles - h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - write(file, "flux", fluxes_cen) - end - h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - write(file, "flux", fluxes_sat) - end - - filename = joinpath(output_dir, "cib_$(freq).fits") - - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - m.pixels = m.pixels + m0.pixels - end - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - end - end - -end - -## -function run_all_chunks(output_dir, halo_pos, halo_mass, freqs; N_chunks=2) - # provide views into halo positions and masses for chunks of the halos - N_halos = size(halo_mass, 1) - chunksize = trunc(Integer, N_halos / N_chunks + 1) - chunks = chunk(N_halos, chunksize) - for chunk_index in 1:length(chunks) - left_ind, right_ind = chunks[chunk_index] - println("Chunk ", chunk_index, "/", length(chunks), - " ", left_ind, " ", right_ind) - pos = @view halo_pos[:, left_ind:right_ind] - mass = @view halo_mass[left_ind:right_ind] - write_chunk(output_dir, chunk_index, model, cosmo, - pos, mass, freqs) - end -end -## compute on all chunks, on all halos - -freqs = [ - "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4", - "63.9", "67.8", "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143", - "153", "164", "189", "210", "217", "232", "256", "275", "294", "306", "314", - "340", "353", "375", "409", "467", "525", "545", "584", "643", "729", "817", - "857", "906", "994", "1080" -] -# freqs = ["143"] - -scratch_dir = joinpath(ENV["SCRATCH"], "cib_sources") -println("SCRATCH: ", scratch_dir) -mkpath(scratch_dir) -mkpath(joinpath(scratch_dir, "sources")) -run_all_chunks(scratch_dir, halo_pos, halo_mass, freqs) - -## diff --git a/examples/cib_planck2013_chunked_scarfy_zevo.jl b/examples/cib_planck2013_chunked_scarfy_zevo.jl deleted file mode 100644 index 109f4f1..0000000 --- a/examples/cib_planck2013_chunked_scarfy_zevo.jl +++ /dev/null @@ -1,98 +0,0 @@ -using XGPaint -using Healpix -using HDF5 -using Random -Random.seed!(3) - -## Load halos from HDF5 files, establish a CIB model and cosmology -@time halo_pos, halo_mass = read_halo_catalog_hdf5("/fs/lustre/cita/zack/projects/websky/websky_halos-light.hdf5") -cosmo = get_cosmology(h=0.677f0, OmegaM=0.310f0) -model = CIB_Planck2013{Float32}(nside=4096,z_evo="scarfy",quench=true) - -## Write one chunk to disk - -function write_chunk( - output_dir, chunk_index, model, cosmo, - pos, mass, freqs) - # Allocate some arrays and fill them up for centrals and satellites - @time sources = generate_sources(model, cosmo, pos, mass); - - # Deposit the sources into maps - fluxes_cen = Array{Float32, 1}(undef, sources.N_cen) - fluxes_sat = Array{Float32, 1}(undef, sources.N_sat) - m = HealpixMap{Float64,RingOrder}(model.nside) - - #h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index).h5"), "w") do file - # write(file, "redshift", sources.redshift_cen) - # write(file, "theta", sources.theta_cen) - # write(file, "phi", sources.phi_cen) - #end - #h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index).h5"), "w") do file - # write(file, "redshift", sources.redshift_sat) - # write(file, "theta", sources.theta_sat) - # write(file, "phi", sources.phi_sat) - #end - - # loop over all frequencies and paint sources to appropriate freq map - @time begin - for freq in freqs - - XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources, - fluxes_cen, fluxes_sat) - - # save sources with mass, redshift, angles - #h5open(joinpath(output_dir, "sources/cen_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - # write(file, "flux", fluxes_cen) - #end - #h5open(joinpath(output_dir, "sources/sat_chunk$(chunk_index)_flux_$(freq).h5"), "w") do file - # write(file, "flux", fluxes_sat) - #end - - filename = joinpath(output_dir, "cib_$(freq).fits") - - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - m.pixels = m.pixels + m0.pixels - end - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - end - end - -end - -## -function run_all_chunks(output_dir, halo_pos, halo_mass, freqs; N_chunks=2) - # provide views into halo positions and masses for chunks of the halos - N_halos = size(halo_mass, 1) - chunksize = trunc(Integer, N_halos / N_chunks + 1) - chunks = chunk(N_halos, chunksize) - for chunk_index in 1:length(chunks) - left_ind, right_ind = chunks[chunk_index] - println("Chunk ", chunk_index, "/", length(chunks), - " ", left_ind, " ", right_ind) - pos = @view halo_pos[:, left_ind:right_ind] - mass = @view halo_mass[left_ind:right_ind] - write_chunk(output_dir, chunk_index, model, cosmo, - pos, mass, freqs) - end -end -## compute on all chunks, on all halos - -""" -freqs = [ - "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4", - "63.9", "67.8", "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143", - "153", "164", "189", "210", "217", "232", "256", "275", "294", "306", "314", - "340", "353", "375", "409", "467", "525", "545", "584", "643", "729", "817", - "857", "906", "994", "1080" -] -""" -freqs = ["143", "857", "217", "353", "545"] - -scratch_dir = "/home/dongwooc/scratchspace/cib_scarfy_test" -println("SCRATCH: ", scratch_dir) -mkpath(scratch_dir) -mkpath(joinpath(scratch_dir, "sources")) -run_all_chunks(scratch_dir, halo_pos, halo_mass, freqs) - -## diff --git a/examples/data/3crr.dat b/examples/data/3crr.dat deleted file mode 100644 index 7fb1b64..0000000 --- a/examples/data/3crr.dat +++ /dev/null @@ -1,174 +0,0 @@ -3CRR_name IAU_name Redshift Flux178MHz_Jy Sp_index FR_class Flux5GHzcore_mJy Best_radio_map -4C12.03 0007+124 0.156 10.9 0.87 2 3.5 [601] -3C6.1 0013+790 0.8404 14.9 0.68 2 4.4 [530] -3C9 0017+154 2.012 19.4 1.12 2 4.9 [504] -3C13 0031+391 1.351 13.1 0.93 2 0.36 [503,537] -3C14 0033+183 1.469 11.3 0.81 2 10.6 [527] -3C16 0035+130 0.405 12.2 0.94 2 0.26 [533,601] -3C19 0038+328 0.482 13.2 0.63 2 0.47 [533] -3C20 0040+517 0.174 46.8 0.66 2 2.6 [605,505,522] -3C22 0048+509 0.938 13.2 0.78 2 7.83 [528] -3C28 0053+261 0.1952 17.8 1.06 2 <0.2 [506] -3C31 0104+321 0.0167 18.3 0.57 1 92 [606] -3C33 0106+130 0.0595 59.3 0.76 2 24 [601] -3C33.1 0106+729 0.181 14.2 0.62 2 12.3 [506] -3C34 0107+315 0.689 13.0 1.06 2 1.03 [529,505] -3C35 0109+492 0.0677 11.4 0.77 2 11.3 [506] -3C41 0123+329 0.795 11.6 0.51 2 0.64 [503,530] -3C42 0125+287 0.395 13.1 0.73 2 2.4 [601,505] -3C43 0127+233 1.47 12.6 0.75 0 <61 [531] -3C46 0132+376 0.4373 11.1 1.13 2 2.3 [506,533] -3C47 0133+207 0.425 28.8 0.98 2 73.6 [506,504] -3C48 0134+329 0.367 60.0 0.59 0 896 [607] -3C49 0138+136 0.6207 11.2 0.65 2 7.0 [531] -3C55 0154+286 0.735 23.4 1.04 2 5.32 [610,528] -3C61.1 0210+860 0.186 34.0 0.77 2 2.4 [601] -3C65 0220+397 1.176 16.6 0.75 2 0.52 [503] -3C66B 0220+427 0.0215 26.8 0.50 1 182 [76,608] -3C67 0221+276 0.3102 10.9 0.58 2 5.5 [510,609] -3C68.1 0229+341 1.238 14.0 0.80 2 1.1 [505] -3C68.2 0231+313 1.575 10.9 1.05 2 0.42 [503,610] -3C76.1 0300+162 0.0324 13.3 0.77 1 10 [601] -3C79 0307+169 0.2559 33.2 0.92 2 10 [140,522] -3C83.1B 0314+416 0.0255 29.0 0.62 1 40 [100] -3C84 0316+413 0.0177 66.8 0.78 1 59600 [512] -3C98 0356+102 0.0306 51.4 0.78 2 9.0 [613] -3C109 0410+110 0.3056 23.5 0.85 2 263 [614] -4C14.11 0411+141 0.206 12.1 0.84 2 29.69 [601,522] -3C123 0433+295 0.2177 206.0 0.70 2 100 [522] -3C132 0453+227 0.214 14.9 0.68 2 2.29 [506,522] -3C138 0518+165 0.759 24.2 0.46 0 94 [39] -3C147 0538+498 0.5450 65.9 0.46 0 2500 [531,615] -3C153 0605+480 0.2769 16.7 0.66 2 <1.7 [506,522] -3C171 0651+542 0.2384 21.3 0.87 2 2.2 [522] -3C172 0659+253 0.5191 16.5 0.86 2 0.43 [533] -3C173.1 0702+749 0.292 16.8 0.88 2 7.4 [601,522] -3C175 0710+118 0.768 19.2 0.98 2 23.5 [504] -3C175.1 0711+146 0.920 12.4 0.91 2 11 [530] -3C181 0725+147 1.382 15.8 1.0 2 6 [617] -3C184 0733+705 0.994 14.4 0.86 2 <0.2 [603] -3C184.1 0734+805 0.1187 14.2 0.68 2 6 [601,613] -3C186 0740+380 1.063 15.4 1.15 2 15 [531] -DA240 0745+560 0.0356 23.2 0.77 2 105 [506] -3C190 0758+143 1.197 16.4 0.93 2 73 [531] -3C191 0802+103 1.9523 14.2 0.98 2 42 [534] -3C192 0802+243 0.0598 23.0 0.79 2 8 [613] -3C196 0809+483 0.871 74.3 0.79 2 7 [618] -3C200 0824+294 0.458 12.3 0.84 2 35.1 [25] -4C14.27 0832+143 0.3920 11.2 1.15 2 <0.25 [601,533] -3C204 0833+654 1.112 11.4 1.08 2 26.9 [504] -3C205 0835+580 1.534 13.7 0.88 2 20 [535] -3C207 0838+133 0.684 14.8 0.90 2 510 [519] -3C208 0850+140 1.109 18.3 0.96 2 51.0 [504] -3C212 0855+143 1.049 16.5 0.92 2 150 [507] -3C215 0903+169 0.411 12.4 1.06 2 16.4 [504] -3C217 0905+380 0.8975 12.3 0.77 2 <0.6 [503] -3C216 0906+430 0.668 22.0 0.84 2 1050 [620] -3C219 0917+458 0.1744 44.9 0.81 2 51 [621] -3C220.1 0926+793 0.61 17.2 0.93 2 25 [25] -3C220.3 0931+836 0.685 17.1 0.75 2 <0.17 [57] -3C223 0936+361 0.1368 16.0 0.74 2 9 [601] -3C225B 0939+139 0.58 23.2 0.94 2 1.1 [533] -3C226 0941+100 0.82 16.4 0.88 2 7.5 [503] -4C73.08 0945+734 0.0581 15.6 0.85 2 11 [506] -3C228 0947+145 0.5524 23.8 1.0 2 13.3 [529] -3C231 0951+699 0.0013 15.9 0.28 1 -- [529] -3C234 0958+290 0.1848 34.2 0.86 2 90 [522] -3C236 1003+351 0.0989 15.7 0.51 2 84 [623] -3C239 1008+467 1.781 14.4 1.08 2 0.5 [503] -4C74.16 1009+748 0.81 12.8 0.87 2 1.5 [71] -3C241 1019+222 1.617 12.6 0.97 2 3 [38] -3C244.1 1030+585 0.428 22.1 0.82 2 1.98 [506,505] -3C245 1040+123 1.029 15.7 0.78 2 910 [626,527] -3C247 1056+432 0.7489 11.6 0.61 2 3.5 [503] -3C249.1 1100+772 0.311 11.7 0.81 2 71 [504] -3C252 1108+359 1.105 12.0 1.03 2 1.1 [505] -3C254 1111+408 0.734 21.7 0.96 2 19 [80] -3C263 1137+660 0.652 16.6 0.82 2 157 [504] -3C263.1 1140+223 0.824 19.8 0.87 2 3.2 [530] -3C264 1142+198 0.0208 28.3 0.75 1 200 [627,515] -3C265 1142+318 0.8108 21.3 0.96 2 2.89 [503,610,528] -3C266 1143+500 1.272 12.1 1.01 2 <0.17 [503] -3C267 1147+130 1.144 15.9 0.93 2 3.0 [503,610] -3C268.1 1157+732 0.9731 23.3 0.59 2 2.0 [610] -3C268.3 1203+645 0.371 11.7 0.50 2 1.0 [531] -3C268.4 1206+439 1.400 11.2 0.80 2 50 [80] -3C270.1 1218+339 1.519 14.8 0.75 2 190 [527] -3C272.1 1222+131 0.0029 21.1 0.60 1 180 [70] -A1552 1227+119 0.0837 12.5 0.94 1 27 [525] -3C274 1228+126 0.0041 1144.5 0.79 1 4e3 [630] -3C274.1 1232+216 0.422 18.0 0.87 2 3.0 [533] -3C275.1 1241+166 0.557 19.9 0.96 2 130 [533] -3C277.2 1251+159 0.766 13.1 1.02 2 0.48 [503] -3C280 1254+476 0.996 25.8 0.81 2 1.0 [80] -3C284 1308+277 0.2394 12.3 0.95 2 3.2 [629,522] -3C285 1319+428 0.0794 12.3 0.95 2 6 [632] -3C287 1328+254 1.055 17.8 0.42 0 -- [39] -3C286 1328+307 0.849 27.3 0.24 0 5554 [534,633] -3C288 1336+391 0.246 20.6 0.85 1 30 [634] -3C289 1343+500 0.9674 13.1 0.81 2 10.6 [80] -3C292 1349+647 0.71 11.0 0.80 2 1.0 [629] -3C293 1350+316 0.0452 13.8 0.45 1 100 [14] -3C294 1404+344 1.786 11.2 1.07 2 0.53 [80] -3C295 1409+524 0.4614 91.0 0.63 2 3 [517] -3C296 1414+110 0.0237 14.2 0.67 1 77 [636] -3C299 1419+419 0.367 12.9 0.65 2 1.0 [80] -3C300 1420+198 0.272 19.5 0.78 2 9 [522] -3C303 1441+522 0.141 12.2 0.76 2 150 [601] -3C305 1448+634 0.0417 17.1 0.85 1 <1.0 [506] -3C309.1 1458+718 0.904 24.7 0.53 2 2350 [534] -3C310 1502+262 0.0540 60.1 0.92 1 80 [637] -3C314.1 1510+709 0.1197 11.6 0.95 1 <1.0 [601] -3C315 1511+263 0.1083 19.4 0.72 1 <150 [506] -3C318 1517+204 1.574 13.4 0.78 2 <90 [638] -3C319 1522+546 0.192 16.7 0.90 2 <1.0 [522] -3C321 1529+242 0.096 14.7 0.60 2 30 [639] -3C322 1533+557 1.681 11.0 0.81 2 0.3 [537] -3C324 1547+215 1.2063 17.2 0.90 2 <0.14 [528,640] -3C326 1549+202 0.0895 22.2 0.88 2 13 [506,623] -3C325 1549+628 0.86 17.0 0.70 2 2.4 [505] -3C330 1609+660 0.5490 30.3 0.71 2 0.74 [505] -NGC6109 1615+351 0.0296 11.7 0.76 1 28 [511,506] -3C334 1618+177 0.555 11.9 0.86 2 111 [504] -3C336 1622+238 0.927 12.5 0.73 2 20.4 [504] -3C341 1626+278 0.448 11.8 0.85 2 1.0 [601] -3C338 1626+396 0.0303 51.1 1.19 1 105 [642] -3C340 1627+234 0.7754 11.0 0.73 2 1.1 [529] -3C337 1627+444 0.635 12.9 0.63 2 0.29 [503,106] -3C343 1634+628 0.988 13.5 0.37 0 <300 [38] -3C343.1 1637+626 0.750 12.5 0.32 2 <200 [38] -NGC6251 1637+826 0.024 11.6 0.72 1 350 [108] -3C346 1641+173 0.162 11.9 0.52 1 220 [506] -3C345 1641+399 0.594 11.8 0.27 0 8610 [644] -3C349 1658+471 0.205 14.5 0.74 2 25 [601,522] -3C351 1704+608 0.371 14.9 0.73 2 6.5 [504] -3C352 1709+460 0.806 12.3 0.88 2 3.18 [503] -3C356 1723+510 1.079 12.3 1.02 2 1.1 [528] -4C16.49 1732+160 1.296 11.4 1.0 2 16 [523] -4C13.66 1759+137 1.450 12.3 0.81 2 <2 [302] -3C368 1802+110 1.132 15.0 1.24 2 <0.21 [640] -3C380 1828+487 0.691 64.7 0.71 1 7447 [171] -3C381 1832+474 0.1605 18.1 0.81 2 5 [601,522] -3C382 1833+326 0.0578 21.7 0.59 2 188 [601,646] -3C386 1836+171 0.0177 26.1 0.59 1 120 [601] -3C388 1842+455 0.0908 26.8 0.70 2 62 [647] -3C390.3 1845+797 0.0569 51.8 0.75 2 330 [648] -3C401 1939+605 0.201 22.8 0.71 2 32 [506,522] -3C427.1 2104+763 0.572 29.0 0.97 2 0.8 [533] -3C432 2120+168 1.805 12.0 0.98 2 7.5 [504] -3C433 2121+248 0.1016 61.3 0.75 2 5 [646] -3C436 2141+279 0.2145 19.4 0.86 2 19 [506,522] -3C437 2145+151 1.480 15.9 0.79 2 <0.12 [503] -3C438 2153+377 0.290 48.7 0.88 2 7.1 [506,522] -3C441 2203+292 0.708 13.7 0.83 2 3.5 [505] -3C442A 2212+135 0.027 17.5 0.96 1 2 [649] -3C449 2229+390 0.0171 12.5 0.58 1 37 [506,650] -3C452 2243+394 0.0811 59.3 0.78 2 130 [646] -NGC7385 2247+113 0.0243 11.7 0.75 1 121 [506,539] -3C454 2249+185 1.757 12.6 0.90 0 <200 [531] -3C454.3 2251+158 0.859 14.2 0.04 0 12200 [651] -3C455 2252+129 0.5427 14.0 0.71 2 1.4 [519] -3C457 2309+184 0.428 14.3 1.01 2 2.1 [601,533] -3C465 2335+267 0.0293 41.2 0.75 1 270 [506,652] -3C469.1 2352+796 1.336 12.1 0.96 2 2.6 [653] -3C470 2356+438 1.653 11.0 0.77 2 2.0 [503] diff --git a/examples/data/sehgal_fig9_1.4GHz.txt b/examples/data/sehgal_fig9_1.4GHz.txt deleted file mode 100644 index f1d5497..0000000 --- a/examples/data/sehgal_fig9_1.4GHz.txt +++ /dev/null @@ -1,27 +0,0 @@ -0.00015486888004572868, 2.9190274700984284 -0.0002499478827893005, 3.861952354997236 -0.0003967951071321661, 5.052608741045338 -0.0006299167461912456, 6.684470970743401 -0.0009836294910426067, 9.042453561398359 -0.0015615230060004996, 12.232715901895336 -0.002478935534643514, 16.92172847690928 -0.003935338359604271, 24.20440436465614 -0.006247394411085373, 34.62135629035885 -0.00991780969288383, 49.521495894947364 -0.015744635704402734, 70.8342716300861 -0.02499478827893, 101.31951684191438 -0.03967951071321653, 140.15704829251132 -0.0640400427119727, 177.3426575209166 -0.0999999999999998, 216.99459808915134 -0.15875113751879782, 242.8628403388844 -0.2520192366351221, 259.95774268747783 -0.4000834049244464, 260.2488112085345 -0.635136956341498, 278.5674967997081 -1.0082871429943965, 282.00645800678825 -1.600667308959392, 291.92608035279227 -2.541077560864565, 264.3475917140384 -4.03398953310742, 273.6460606087317 -6.404004271197257, 292.90776683623056 -10.166429627278065, 224.3851567599586 -16.139322678351974, 313.8763345043145 -25.62135833971307, 179.93424204944841 diff --git a/examples/data/sehgal_fig9_145GHz.txt b/examples/data/sehgal_fig9_145GHz.txt deleted file mode 100644 index f262d79..0000000 --- a/examples/data/sehgal_fig9_145GHz.txt +++ /dev/null @@ -1,20 +0,0 @@ -0.0006350631327313502, 0.48496934285281956 -0.001, 0.6738627168030941 -0.0015746465956842633, 0.8576958985908937 -0.0025200458978585096, 1.0680004325145744 -0.003968181694030979, 1.245197084735032 -0.006350631327313502, 1.467799267622069 -0.01, 1.6559515234819164 -0.016003882425208067, 1.9952623149688788 -0.025200458978585097, 2.2758459260747887 -0.040330518255455486, 2.742175448122083 -0.06350631327313502, 2.742175448122083 -0.10163475708816792, 4.1595621630718425 -0.15746465956842634, 4.346055200262694 -0.25612425267988265, 5.3526818228471 -0.4033051825545549, 6.520572293428611 -0.655996328108289, 8.863781619257626 -0.9999999999999959, 4.903161412303259 -1.6265507027536286, 5.65432740962822 -2.5612425267988264, 11.15883992507748 -4.033051825545549, 5.531681197617221 diff --git a/examples/data/sehgal_fig9_33GHz.txt b/examples/data/sehgal_fig9_33GHz.txt deleted file mode 100644 index dc6e5e0..0000000 --- a/examples/data/sehgal_fig9_33GHz.txt +++ /dev/null @@ -1,21 +0,0 @@ -0.0006258315400152105, 0.9982701308273131 -0.000997538446877355, 1.4226286843916296 -0.0015900172640297115, 2.046119432823881 -0.0024835385036582587, 2.886991698970413 -0.003958613433958853, 3.893260033392083 -0.00624616886683321, 4.92182872528542 -0.009956023613324207, 5.828432831981491 -0.016030965671921877, 6.787884809066023 -0.025039678120311156, 7.267264700281397 -0.03991176542786095, 8.09698269862269 -0.0623404590827242, 8.809527433266034 -0.09936700334644462, 9.703087018653362 -0.15838512419278228, 8.782550063461155 -0.25245651695966786, 11.25928649094335 -0.39432620925119255, 9.526164064833942 -0.6285326468683717, 9.777771227262617 -1.0018438513877874, 7.662096252303133 -1.596879824086066, 4.136375730631529 -2.4942575269790654, 8.336370428761208 -6.337028857216717, 44.38519459717867 -9.998994096421844, 22.154492861369906 diff --git a/examples/data/sehgal_fig9_95GHz.txt b/examples/data/sehgal_fig9_95GHz.txt deleted file mode 100644 index 486fd31..0000000 --- a/examples/data/sehgal_fig9_95GHz.txt +++ /dev/null @@ -1,20 +0,0 @@ -0.0006195817341645596, 0.604104637652916 -0.0009771879850411914, 0.8400302889814686 -0.00155703867926312, 1.149190838317887 -0.0024818051496490992, 1.4538947716023338 -0.003916732764601186, 1.7449668785679036 -0.00618111290995378, 2.109431315454894 -0.00985795012311407, 2.3343870745268465 -0.015559773208435941, 2.7131749472584623 -0.02456250140963686, 3.0662866813146192 -0.03915813644219979, 3.7154348412066924 -0.061821310640063395, 4.095616752264089 -0.09858147929091977, 4.686379868342065 -0.15565503327296254, 5.025938135702978 -0.2480784921270924, 6.50335405551268 -0.38769411765212364, 7.187892563248658 -0.6181785008012204, 8.367829034270141 -0.9780606285591938, 5.6076218786688585 -1.5652892238489837, 2.781218966936857 -2.457448638834707, 11.143689679736001 -3.9141721226413044, 16.65056697504314 diff --git a/examples/data/sehgal_figure8_blue.txt b/examples/data/sehgal_figure8_blue.txt deleted file mode 100644 index c5beaeb..0000000 --- a/examples/data/sehgal_figure8_blue.txt +++ /dev/null @@ -1,6 +0,0 @@ -5.588713261387463e+25, 1.1155306488883394e-7 -1.7564009199951863e+26, 2.4448923533374074e-7 -5.519954321281483e+26, 7.737709943086199e-8 -1.7618451770855528e+27, 1.9325602208071095e-8 -5.537064338766726e+27, 7.750373260528731e-9 -1.7536851034159255e+28, 1.7266147069691356e-9 \ No newline at end of file diff --git a/examples/data/sehgal_figure8_green.txt b/examples/data/sehgal_figure8_green.txt deleted file mode 100644 index 728dd57..0000000 --- a/examples/data/sehgal_figure8_green.txt +++ /dev/null @@ -1,5 +0,0 @@ -1.891828311354444e+25, 0.00000858366364440464 -5.945570708544333e+25, 0.0000011155796205314259 -1.8976923481928418e+26, 8.054045337796261e-8 -5.963999994886519e+26, 2.8344321421904224e-8 -5.982486405870227e+27, 8.479173899075899e-10 \ No newline at end of file diff --git a/examples/data/sehgal_figure8_red.txt b/examples/data/sehgal_figure8_red.txt deleted file mode 100644 index 06e0f78..0000000 --- a/examples/data/sehgal_figure8_red.txt +++ /dev/null @@ -1,7 +0,0 @@ -5.848512932080513e+24, 0.000019905396851365134 -1.8374035259917507e+25, 0.000003612828459430432 -5.772495943042983e+25, 9.282548743463857e-7 -1.813521577654646e+26, 2.384993147124578e-7 -5.6974669970668993e+26, 1.1455265142927872e-7 -1.8494621790323921e+27, 3.753552880116237e-8 -6.003563256149004e+27, 8.995615307242464e-9 \ No newline at end of file diff --git a/examples/data/table.png b/examples/data/table.png deleted file mode 100644 index be89555..0000000 Binary files a/examples/data/table.png and /dev/null differ diff --git a/examples/data/wilman_151.txt b/examples/data/wilman_151.txt deleted file mode 100644 index 85664f8..0000000 --- a/examples/data/wilman_151.txt +++ /dev/null @@ -1,13 +0,0 @@ -0.00003608918289407965, 0.00024690129310615084 -0.0000828004671228646, 0.0011135993581813265 -0.00027709133868033484, 0.004006856031867952 -0.0008598569967363336, 0.019485862163414882 -0.0036089182894079646, 0.07011232879402685 -0.016334835079813696, 0.29328491084381414 -0.12541822663794908, 2.24113124144541 -1.302429121962331, 19.90973229455024 -8.59856996736331, 71.63746131066385 -38.919202673929355, 257.75966181331216 -404.1637673000104, 1000 -1572.9731573007896, 2123.7471239356023 -5263.946605735404, 2289.8818989855477 diff --git a/examples/demo_cib.jl b/examples/demo_cib.jl deleted file mode 100644 index 3ed8499..0000000 --- a/examples/demo_cib.jl +++ /dev/null @@ -1,79 +0,0 @@ -using XGPaint -using Healpix - -# example , ra and dec in radians, halo mass in M200c (Msun) -ra, dec, redshift, halo_mass = XGPaint.load_example_halos() - -# sort -ra, dec, redshift, halo_mass = sort_halo_catalog(ra, dec, redshift, halo_mass) - -print("Number of halos: ", length(halo_mass)) - -## -cosmo = get_cosmology(h=0.6774f0, OmegaM=0.3075f0) - -x, y, z = XGPaint.ra_dec_redshift_to_xyz(ra, dec, redshift, cosmo) -halo_pos = [x'; y'; z';] - -model = CIB_Planck2013{Float32}() - -## Allocate some arrays and file them up for centrals and satellites -@time sources = generate_sources(model, cosmo, halo_pos, halo_mass); - -## Deposit the sources into maps -fluxes_cen = Array{Float32, 1}(undef, sources.N_cen) -fluxes_sat = Array{Float32, 1}(undef, sources.N_sat) - -using Pixell -box = [4.5 -4.5; # RA - -3 3] * Pixell.degree # DEC -shape, wcs = geometry(CarClenshawCurtis{Float64}, box, 0.5 * Pixell.arcminute) -m = Enmap(zeros(Float32, shape), wcs) -XGPaint.paint!(m, 143.0f0 * 1.0f9, model, sources, fluxes_cen, fluxes_sat) - -## -using Plots -plot(log10.(m), c=:coolwarm) - - - - -## - - -## -function get_basic_halo_properties(halo_pos::Array{T,2}, - cosmo) where T - N_halos = size(halo_pos, 2) - redshift = Array{T}(undef, N_halos) - dist = Array{T}(undef, N_halos) - - r2z = XGPaint.build_r2z_interpolator( - 0.0, 4.5, cosmo) - Threads.@threads :static for i in 1:N_halos - dist[i] = sqrt(halo_pos[1,i]^2 + halo_pos[2,i]^2 + halo_pos[3,i]^2) - redshift[i] = r2z(dist[i]) - end - - return dist, redshift -end - -d̂, ẑ = get_basic_halo_properties([x'; y'; z'], cosmo) - -θ, ϕ = XGPaint.get_angles([x'; y'; z']) -ra_ = ϕ -dec_ = π/2 .- θ - -## -maximum((redshift .- ẑ)./redshift) - - -## -clip(x) = (x + 2π) % 2π - -maximum(abs.((clip.(ra_) .- clip.(ra)))), maximum(abs.(dec_ .- dec)) - -## - - -## diff --git a/examples/demo_tsz.jl b/examples/demo_tsz.jl deleted file mode 100644 index 27a29ed..0000000 --- a/examples/demo_tsz.jl +++ /dev/null @@ -1,40 +0,0 @@ -using XGPaint, Plots - -# example , ra and dec in radians, halo mass in M200c (Msun) -ra, dec, redshift, halo_mass = XGPaint.load_example_halos() -print("Number of halos: ", length(halo_mass)) - -## -b = 10.0 .^ (12:0.25:15) -histogram(halo_mass, bins=b, xaxis=(:log10, (1e12, 1e15)), yscale=:log10, label="", - xlabel="Halo mass (solar masses)", ylabel="number") - -## -print("Precomputing the model profile grid.\n") - -model = Battaglia16ThermalSZProfile(Omega_c=0.2589, Omega_b=0.0486, h=0.6774) -interp = build_interpolator(model, cache_file="cached_b16.jld2", overwrite=true) - - -ra, dec, redshift, halo_mass = sort_halo_catalog(ra, dec, redshift, halo_mass) - -## -using Pixell -box = [4.5 -4.5; # RA - -3 3] * Pixell.degree # DEC -shape, wcs = geometry(CarClenshawCurtis{Float64}, box, 0.5 * Pixell.arcminute) -m = Enmap(zeros(shape), wcs) - -workspace = profileworkspace(shape, wcs) - -@time paint!(m, model, workspace, interp, halo_mass, redshift, ra, dec) -# @time paint_map!(m, p, psa, sitp, halo_mass, redshift, ra, dec, 1:length(halo_mass)) - -## -plot(log10.(m), c = :thermal) - -## -plot((m), clim=(-1e-5, 1e-5)) - - - diff --git a/examples/lrg_chunked.jl b/examples/lrg_chunked.jl deleted file mode 100644 index 4b77681..0000000 --- a/examples/lrg_chunked.jl +++ /dev/null @@ -1,82 +0,0 @@ -using XGPaint -using Healpix -using JLD2, FileIO, CodecZlib -using HDF5 -using DelimitedFiles -using Random -Random.seed!(3) - -hdata = h5open("/home/dongwooc/projectscratchspace/websky_halos_rewrite/websky_halos-lesslight_20230612.h5","r") do file - (x = read(file,"x"), y = read(file,"y"), z = read(file,"z"), m = read(file,"M200m"), vrad = read(file,"vrad") ) -end -halo_pos = vcat(hdata.x',hdata.y',hdata.z') -halo_mass = hdata.m -halo_vrad = hdata.vrad -cosmo = get_cosmology(h=0.677f0, OmegaM=0.310f0) -model = LRG_Yuan23{Float32}(nside=4096) - -lrg_redshifts = [ - "0.35", "0.40", "0.45", "0.50", "0.55", - "0.60", "0.65", "0.70", "0.75", "0.80" -] - -function write_chunk( - output_dir, chunk_index, model, cosmo, - pos, mass, vrad, lrg_redshifts) - sources = generate_sources(model, cosmo, pos, mass, vrad); - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index).jld2"), "w") do file - file["redshift"]=sources.redshift_cen - file["theta"]=sources.theta_cen - file["phi"]=sources.phi_cen - write(file,"LRG",sources.lrg_cen) - file["m200c"]=sources.m200c_cen - file["hp_ind_cen"]=sources.hp_ind_cen - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index).jld2"), "w") do file - file["redshift"]=sources.redshift_sat - file["theta"]=sources.theta_sat - file["phi"]=sources.phi_sat - write(file,"LRG",sources.lrg_sat) - file["m200m"]=sources.m200m_sat - file["m200c"]=sources.m200c_sat - file["parent"]=sources.parent_sat - file["hp_ind_sat"]=sources.hp_ind_sat - end - m = HealpixMap{Float64,RingOrder}(model.nside) - for zlo in lrg_redshifts - XGPaint.paint!(m, model, sources, - parse(Float32, zlo), parse(Float32, zlo)+0.05f0) - # save maps - filename = joinpath(output_dir, "lrg_$(zlo).fits") - - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - m.pixels = m.pixels + m0.pixels - end - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - end -end - -## the rest of this is more or less the same as cib_planck2013_chunked.jl -function run_all_chunks(output_dir, halo_pos, halo_mass, halo_vrad, lrg_redshifts; N_chunks=4) - # provide views into halo positions and masses for chunks of the halos - N_halos = size(halo_mass, 1) - chunksize = trunc(Integer, N_halos / N_chunks + 1) - chunks = chunk(N_halos, chunksize) - for chunk_index in 1:length(chunks) - left_ind, right_ind = chunks[chunk_index] - println("Chunk ", chunk_index, "/", length(chunks), - " ", left_ind, " ", right_ind) - pos = @view halo_pos[:, left_ind:right_ind] - mass = @view halo_mass[left_ind:right_ind] - vrad = @view halo_vrad[left_ind:right_ind] - write_chunk(output_dir, chunk_index, model, cosmo, - pos, mass, vrad, lrg_redshifts) - end -end -## compute on all chunks, on all halos -scratch_dir = "/home/dongwooc/projectscratchspace/lrg_sources_vrad_m200c" -println("SCRATCH: ", scratch_dir) -mkpath(scratch_dir) -mkpath(joinpath(scratch_dir, "sources")) -run_all_chunks(scratch_dir, halo_pos, halo_mass, halo_vrad, lrg_redshifts) diff --git a/examples/lrg_chunked_cibmap.jl b/examples/lrg_chunked_cibmap.jl deleted file mode 100644 index c3a21bb..0000000 --- a/examples/lrg_chunked_cibmap.jl +++ /dev/null @@ -1,62 +0,0 @@ -# OBSOLETE! DO NOT USE WITH CURRENT XGPAINT -# hackish script for generating LRG/CIB maps after the fact -using XGPaint -using Healpix -using JLD2, FileIO, CodecZlib -using HDF5 -using DelimitedFiles -using Random -Random.seed!(3) - -cosmo = get_cosmology(h=0.677f0, OmegaM=0.310f0) -model = LRG_Yuan22{Float32}(nside=4096) - -freqs = ["353", "545", "857", - "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4", - "63.9", "67.8", "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143", - "153", "164", "189", "210", "217", "232", "256", "275", "294", "306", "314", - "340", "375", "409", "467", "525", "584", "643", "729", "817", - "906", "994", "1080" -] - -function write_chunk( - lrg_dir, cib_dir, chunk_index, model, cosmo, freqs) - m = HealpixMap{Float64,RingOrder}(model.nside) - jldopen(joinpath(lrg_dir, "sources/cen_chunk$(chunk_index).jld2"), "r") do file - println("loading centrals") - lrg_cen = file["LRG"] - println(sum(lrg_cen),"/",length(lrg_cen)) - hp_ind_lrgcen = file["hp_ind_cen"][lrg_cen] - # process centrals for this frequency - for freq in freqs - print(freq,"...") - fill!(m.pixels, zero(Float32)) - h5open(joinpath(cib_dir, "sources/cen_chunk$(chunk_index)_flux_$(freq).h5"), "r") do file2 - fluxes_lrgcen = read(file2["flux"])[lrg_cen] - Threads.@threads :static for i in 1:length(hp_ind_lrgcen) - m.pixels[hp_ind_lrgcen[i]] += fluxes_lrgcen[i] - end - end - # save maps - filename = joinpath(lrg_dir, "lrg_cen_cib_$(freq).fits") - - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - m.pixels = m.pixels + m0.pixels - end - Healpix.saveToFITS(m, "!$(filename)", typechar="D") - end - end - println("chunk done") -end - -function run_all_chunks(lrg_dir, cib_dir, freqs, N_chunks) - for chunk_index in 1:N_chunks - println("Chunk ", chunk_index, "/", N_chunks) - write_chunk(lrg_dir, cib_dir, chunk_index, model, cosmo, freqs) - end -end -## compute on all chunks, on all halos -scratch_dir = "/home/dongwooc/scratchspace/lrg_sources" -scratch2dir = "/home/dongwooc/scratchspace/cib_co_sources" -run_all_chunks(scratch_dir, scratch2dir, freqs, 4) diff --git a/examples/lrg_chunked_postmap.jl b/examples/lrg_chunked_postmap.jl deleted file mode 100644 index b9c9b57..0000000 --- a/examples/lrg_chunked_postmap.jl +++ /dev/null @@ -1,63 +0,0 @@ -# OBSOLETE! DO NOT USE WITH CURRENT XGPAINT -# hackish script for generating LRG maps after the fact -using XGPaint -using Healpix -using JLD2, FileIO, CodecZlib -using DelimitedFiles -using Random -Random.seed!(3) - -cosmo = get_cosmology(h=0.677f0, OmegaM=0.310f0) -model = LRG_Yuan22{Float32}(nside=4096) - -zlo = 0.64 -zhi = 0.86 - -function write_chunk( - output_dir, chunk_index, model, cosmo, zlo, zhi) - m = HealpixMap{Float64,RingOrder}(model.nside) - fill!(m.pixels, zero(Float32)) - jldopen(joinpath(output_dir, "sources/cen_chunk$(chunk_index).jld2"), "r") do file - redshift_cen = file["redshift"] - lrg_cen = file["LRG"] - hp_ind_cen = file["hp_ind_cen"] - # process centrals for this frequency - Threads.@threads :static for i in 1:length(redshift_cen) - if (lrg_cen[i]) && (redshift_cen[i] < zhi) && (redshift_cen[i] > zlo) - m.pixels[hp_ind_cen[i]] += 1 - end - end - end - jldopen(joinpath(output_dir, "sources/sat_chunk$(chunk_index).jld2"), "r") do file - redshift_sat = file["redshift"] - lrg_sat = file["LRG"] - hp_ind_sat = file["hp_ind_sat"] - # process satellites for this frequency - Threads.@threads :static for i in 1:length(redshift_sat) - if (lrg_sat[i]) && (redshift_sat[i] < zhi) && (redshift_sat[i] > zlo) - m.pixels[hp_ind_sat[i]] += 1 - end - end - end - # save maps - filename = joinpath(output_dir, "lrg_$(zlo)_$(zhi).fits") - - if chunk_index > 1 - m0 = Healpix.readMapFromFITS(filename, 1, Float32) - m.pixels = m.pixels + m0.pixels - end - Healpix.saveToFITS(m, "!$(filename)", typechar="D") -end - -function run_all_chunks(output_dir, zlo, zhi, N_chunks) - for chunk_index in 1:N_chunks - println("Chunk ", chunk_index, "/", N_chunks) - write_chunk(output_dir, chunk_index, model, cosmo, zlo, zhi) - end -end -## compute on all chunks, on all halos -scratch_dir = "/home/dongwooc/scratchspace/lrg_sources" -println("SCRATCH: ", scratch_dir) -mkpath(scratch_dir) -mkpath(joinpath(scratch_dir, "sources")) -run_all_chunks(scratch_dir, zlo, zhi, 4) diff --git a/examples/make_map_raw_car_martine.jl b/examples/make_map_raw_car_martine.jl deleted file mode 100644 index d97b176..0000000 --- a/examples/make_map_raw_car_martine.jl +++ /dev/null @@ -1,204 +0,0 @@ -# modified from Zack Li - -using Pixell, WCS, XGPaint -using Cosmology -using Interpolations -import XGPaint: AbstractProfile -using HDF5 -import JSON -using JLD2 -using ThreadsX -using FileIO -using Healpix - -print("Threads: ", Threads.nthreads(), "\n") -modeltype::String = "break" -healpix = true - - -if healpix - nside = 1024 - m = HealpixMap{Float64,RingOrder}(nside) -else - shape, wcs = fullsky_geometry(0.5 * Pixell.arcminute) - # for a box: - # box = [30 -30; # RA - # -25 25] * Pixell.degree # DEC - # shape, wcs = geometry(CarClenshawCurtis{Float64}, box, 0.5 * Pixell.arcminute) - m = Enmap(zeros(shape), wcs) - # precomputed sky angles - α_map, δ_map = posmap(shape, wcs) - psa = (sin_α=sin.(α_map), cos_α=cos.(α_map), sin_δ=sin.(δ_map), cos_δ=cos.(δ_map)) -end - -## -fid = h5open("/mnt/raid-cita/mlokken/buzzard/catalogs/halos/buzzard_halos.hdf5", "r") -# fid = h5open("/mnt/scratch-lustre/mlokken/pkpatch/halos_hdf5.h5", "r") -ra, dec = deg2rad.(fid["ra"]), deg2rad.(fid["dec"]) -redshift = collect(fid["z"]) -halo_mass = collect(fid["m200c"]) -# choose the cutoff for the pressure profiles here -cutoff = 4 # the standard -apply_beam = true -## -perm = sortperm(dec, alg=ThreadsX.MergeSort) -ra = ra[perm] -dec = dec[perm] -redshift = redshift[perm] -halo_mass = halo_mass[perm] - - - -## -print("Precomputing the model profile grid.\n") - -# set up a profile to paint -if modeltype=="break" - p = XGPaint.BreakModel(Omega_c=0.24, Omega_b=0.046, h=0.7, alpha_break=1.5) # Buzzard cosmology -elseif modeltype=="battaglia" - p = XGPaint.BattagliaProfile(Omega_c=0.24, Omega_b=0.046, h=0.7) # Buzzard cosmology -end - -# beam stuff -N_logθ = 512 -rft = RadialFourierTransform(n=N_logθ, pad=256) -beamstr = "" - -if apply_beam - ells = rft.l - fwhm = deg2rad(1.6/60.) # ACT beam - σ² = fwhm^2 / (8log(2)) - lbeam = @. exp(-ells * (ells+1) * σ² / 2) - beamstr = "_1p6arcmin" -end -model_file::String = "cached_$(modeltype)_$(beamstr).jld2" -if isfile(model_file) - print("Found cached Battaglia profile model. Loading from disk.\n") - model = load(model_file) - prof_logθs, prof_redshift, prof_logMs, prof_y = model["prof_logθs"], - model["prof_redshift"], model["prof_logMs"], model["prof_y"] -else - print("Didn't find a cached profile model. Computing and saving.\n") - logθ_min, logθ_max = log(minimum(rft.r)), log(maximum(rft.r)) - @time prof_logθs, prof_redshift, prof_logMs, prof_y = profile_grid(p; - N_logθ=length(rft.r), logθ_min=logθ_min, logθ_max=logθ_max) - - if apply_beam - # now apply the beam - @time XGPaint.transform_profile_grid!(prof_y, rft, lbeam) - prof_y = prof_y[begin+rft.pad:end-rft.pad, :, :] - prof_logθs = prof_logθs[begin+rft.pad:end-rft.pad] - XGPaint.cleanup_negatives!(prof_y) - end - jldsave(model_file; prof_logθs, prof_redshift, prof_logMs, prof_y) - -end - -if ~healpix # do this for CAR - itp = Interpolations.interpolate(log.(prof_y), BSpline(Cubic(Line(OnGrid())))) - sitp = scale(itp, prof_logθs, prof_redshift, prof_logMs); -end - - -## -function paint_map!(m::Enmap, p::XGPaint.AbstractProfile, psa, sitp, masses, - redshifts, αs, δs, irange; mult=4) - for i in irange - α₀ = αs[i] - δ₀ = δs[i] - mh = masses[i] - z = redshifts[i] - θmax = XGPaint.θmax(p, mh * XGPaint.M_sun, z, mult=mult) - profile_paint!(m, α₀, δ₀, psa, sitp, z, mh, θmax) - end -end - -function paint_map!(m::HealpixMap, p::XGPaint.AbstractProfile{T}, masses, - redshifts, αs, δs, irange; mult=4) where {T} - for i in irange - α₀ = αs[i] - δ₀ = δs[i] - mh = masses[i] - z = redshifts[i] - θmax = XGPaint.θmax(p, mh * XGPaint.M_sun, z, mult=mult) - θmin = 0.0 # don't know what this should be # error - beam_interp = XGPaint.realspacegaussbeam(T,fwhm)[0] #error - w = XGPaint.HealpixPaintingWorkspace(nside, θmin, θmax, beam_interp) - profile_paint!(m, α₀, δ₀, w, z, mh, θmax) - end -end - -function chunked_paint!(m::Enmap, p::XGPaint.AbstractProfile, psa, sitp, masses, - redshifts, αs, δs; mult=4) - m .= 0.0 - - N_sources = length(masses) - chunksize = ceil(Int, N_sources / (2Threads.nthreads())) - chunks = chunk(N_sources, chunksize); - - Threads.@threads for i in 1:Threads.nthreads() - chunk_i = 2i - i1, i2 = chunks[chunk_i] - paint_map!(m, p, psa, sitp, masses, redshifts, αs, δs, i1:i2, mult=mult) - end - - Threads.@threads for i in 1:Threads.nthreads() - chunk_i = 2i - 1 - i1, i2 = chunks[chunk_i] - paint_map!(m, p, psa, sitp, masses, redshifts, αs, δs, i1:i2, mult=mult) - end -end - -function chunked_paint!(m::HealpixMap, p::XGPaint.AbstractProfile, masses, - redshifts, αs, δs; mult=4) - m .= 0.0 - - N_sources = length(masses) - chunksize = ceil(Int, N_sources / (2Threads.nthreads())) - chunks = chunk(N_sources, chunksize); - - Threads.@threads for i in 1:Threads.nthreads() - chunk_i = 2i - i1, i2 = chunks[chunk_i] - paint_map!(m, p, masses, redshifts, αs, δs, i1:i2, mult=mult) - end - - Threads.@threads for i in 1:Threads.nthreads() - chunk_i = 2i - 1 - i1, i2 = chunks[chunk_i] - paint_map!(m, p, masses, redshifts, αs, δs, i1:i2, mult=mult) - end -end -## -# cut = eachindex(halo_mass)[begin:end-5] - -print(maximum(ra), minimum(ra), maximum(dec), minimum(dec)) -print("Painting map.\n") -if healpix - @time chunked_paint!(m, p, halo_mass, redshift, ra, dec, mult=cutoff) # for a Healpix -else - @time chunked_paint!(m, p, psa, sitp, halo_mass, redshift, ra, dec, mult=cutoff) # for an Enmap -end - -# -write_map( - "/mnt/raid-cita/mlokken/buzzard/ymaps/ymap_buzzard_$(modeltype)_bbps_car$(beamstr)_cutoff$(cutoff).fits", - m) - -using PyPlot -plt.clf() -plt.figure() -plt.imshow(log10.(m.data')) -plt.axis("off") -plt.savefig("test_Healpix.png", bbox_inches="tight",pad_inches = 0) -plt.gcf() - - - -# ## -# m .= 0 -# profile_paint!(m, 0., 0., psa, sitp, 0.1, 1e14, π/90) - -# ## -# using Plots -# Plots.plot(log10.(m)) diff --git a/examples/radio_ccat.jl b/examples/radio_ccat.jl deleted file mode 100644 index 9e82f5e..0000000 --- a/examples/radio_ccat.jl +++ /dev/null @@ -1,57 +0,0 @@ -using XGPaint -using Healpix -using HDF5 -import ThreadsX -using JLD2 - -# ENV["SCRATCH"] = "/home/zequnl/scratch/" -@time halo_pos, halo_mass = read_halo_catalog_hdf5( - joinpath(ENV["SCRATCH"],"websky_halos-light.hdf5")); - -## Load halos from HDF5 files, establish a CIB model and cosmology -cosmo = get_cosmology(Float32, h = 0.7, OmegaM = 0.25) -radio_model = Radio_Sehgal2009{Float32}(a_0 = -1) -@time begin - sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass); -end; - -@save joinpath(ENV["SCRATCH"], "CCAT/before_match/sources.h5") sources - -freqs = ["93", "145", "220", "280", "350", "405", "860"] - - -m = Map{Float64,RingOrder}(radio_model.nside) - -function generate_maps() - scratch_dir = ENV["SCRATCH"] - println("SCRATCH: ", scratch_dir) - mkpath(joinpath(scratch_dir, "CCAT/before_match/")) - mkpath(joinpath(scratch_dir, "CCAT/before_match/maps")) - - for freq in freqs - @time begin - - flux_I, redshift_I, θ_I, ϕ_I, flux_II, redshift_II, θ_II, ϕ_II = paint!( - m, parse(Float32, freq) * 1.0f9, radio_model, sources, - return_fluxes=true) - - flux = vcat(flux_I, flux_II) - redshift = vcat(redshift_I, redshift_II) - theta = vcat(θ_I, θ_II) - phi = vcat(ϕ_I, ϕ_II) - - h5open(joinpath(scratch_dir, "CCAT/before_match/catalog_$(freq).h5"), "w") do file - write(file, "flux", flux) - write(file, "redshift", redshift) - write(file, "theta", theta) - write(file, "phi", phi) - end - - map_path = joinpath(scratch_dir, "CCAT/before_match/maps/radio$(freq).fits") - Healpix.saveToFITS(m, "!$(map_path)") - end - end -end - -generate_maps() -## diff --git a/examples/radio_sehgal.jl b/examples/radio_sehgal.jl deleted file mode 100644 index 133fb54..0000000 --- a/examples/radio_sehgal.jl +++ /dev/null @@ -1,79 +0,0 @@ -using XGPaint -using Healpix -using HDF5 -import ThreadsX -using JLD2 - -# ENV["SCRATCH"] = "/home/zequnl/scratch/" -@time halo_pos, halo_mass = read_halo_catalog_hdf5( - joinpath(ENV["SCRATCH"],"websky_halos-light.hdf5")); - -## Load halos from HDF5 files, establish a CIB model and cosmology -cosmo = get_cosmology(Float32, h = 0.7, OmegaM = 0.25) -radio_model = Radio_Sehgal2009{Float32}(a_0 = -1) -@time begin - sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass); -end; - -@save joinpath(scratch_dir, "before_match/catalog_$(freq).h5") sources - - -freqs = [ - "18.7", "21.6", "24.5", "27.3", "30.0", "35.9", "41.7", "44.0", "47.4", "63.9", "67.8", - "70.0", "73.7", "79.6", "90.2", "100", "111", "129", "143", "153", "164", "189", "210", - "217", "232", "256", "275", "294", "306", "314", "340", "353", "375", "409", "467", - "525", "545", "584", "643", "729", "817", "857", "906", "994", "1080" -] - -# freqs = ["143.0"] - -m = HealpixMap{Float64,RingOrder}(radio_model.nside) - -function generate_maps() - scratch_dir = ENV["SCRATCH"] - println("SCRATCH: ", scratch_dir) - - for freq in freqs - @time begin -# flux_I, redshift_I, θ_I, ϕ_I, flux_II, redshift_II, θ_II, ϕ_II = paint!( -# m, parse(Float32, freq) * 1.0f9, radio_model, sources, -# return_fluxes=true) - -# flux = vcat(flux_I, flux_II) -# redshift = vcat(redshift_I, redshift_II) -# theta = vcat(θ_I, θ_II) -# phi = vcat(ϕ_I, ϕ_II) -# perm = sortperm(flux, rev=true, alg=ThreadsX.MergeSort) - -# h5open(joinpath(scratch_dir, "before_match/catalog_$(freq).h5"), "w") do file -# write(file, "flux", flux[perm]) -# write(file, "redshift", redshift[perm]) -# write(file, "theta", theta[perm]) -# write(file, "phi", phi[perm]) -# end - - flux_I, redshift_I, θ_I, ϕ_I, flux_II, redshift_II, θ_II, ϕ_II = paint!( - m, parse(Float32, freq) * 1.0f9, radio_model, sources, - return_fluxes=true) - - flux = vcat(flux_I, flux_II) - redshift = vcat(redshift_I, redshift_II) - theta = vcat(θ_I, θ_II) - phi = vcat(ϕ_I, ϕ_II) - # perm = sortperm(flux, rev=true, alg=ThreadsX.MergeSort) - - h5open(joinpath(scratch_dir, "before_match/catalog_$(freq).h5"), "w") do file - write(file, "flux", flux) - write(file, "redshift", redshift) - write(file, "theta", theta) - write(file, "phi", phi) - end - - map_path = joinpath(scratch_dir, "before_match/maps/radio$(freq).fits") - Healpix.saveToFITS(m, "!$(map_path)") - end - end -end - -generate_maps() -## diff --git a/examples/radio_sehgal_map.jl b/examples/radio_sehgal_map.jl deleted file mode 100644 index dcfe5e4..0000000 --- a/examples/radio_sehgal_map.jl +++ /dev/null @@ -1,25 +0,0 @@ -using XGPaint -using Healpix - -halo_pos, halo_mass = read_halo_catalog_hdf5( - "/tigress/zequnl/xgpaint/websky_halos-light.hdf5") - -## Load halos from HDF5 files, establish a CIB model and cosmology -cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0) -radio_model = Radio_Sehgal2009{Float32}(a_0=-1) - -## - -@time begin - sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass); -end; - -## -m = HealpixMap{Float64,RingOrder}(radio_model.nside) -for freq in ["030", "090", "148", "219", "277", "350"] - @time begin - paint!(m, parse(Float32, freq) * 1.0f9, radio_model, sources) - Healpix.saveToFITS(m, "!/tigress/zequnl/xgpaint/jl/radio$(freq).fits") - end -end -## diff --git a/examples/radio_sehgal_plots.jl b/examples/radio_sehgal_plots.jl deleted file mode 100644 index 0540fdb..0000000 --- a/examples/radio_sehgal_plots.jl +++ /dev/null @@ -1,178 +0,0 @@ -using XGPaint -using Healpix - -halo_pos, halo_mass = read_halo_catalog_hdf5( - "/home/zequnl/Data/websky_halos-light.hdf5") -## Load halos from HDF5 files, establish a CIB model and cosmology -cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0) -radio_model = Radio_Sehgal2009{Float32}() - -## - -@time begin - sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass); -end; - -## -@time begin -m = HealpixMap{Float64,RingOrder}(radio_model.nside) -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 151f6, radio_model, sources, return_fluxes=true) -end -## - - -using DelimitedFiles -using Unitful -using UnitfulAstro -using Cosmology -using PyPlot -using StatsBase - -figure(figsize=(7,7)) - -sehgal_red = readdlm("examples/data/sehgal_figure8_red.txt", ',', Float64, '\n') -sehgal_blue = readdlm("examples/data/sehgal_figure8_blue.txt", ',', Float64, '\n') -sehgal_green = readdlm("examples/data/sehgal_figure8_green.txt", ',', Float64, '\n') - -vol = ustrip(u"Mpc^3",comoving_volume(cosmo, 0.3)) -bins = (range(23, stop=29, step=0.5)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 -count_I = fit(Histogram, log10.(sources.L_I_151[redshift_I .< 0.3]), bins).weights - -count_II = fit(Histogram, log10.(sources.L_II_151[redshift_II .< 0.3]), bins).weights - -plot(mids, (count_I .+ count_II) ./ diff((bins)) ./ vol, label="websky total" ) -plot(mids, count_I ./ diff((bins)) ./ vol, label="websky FR I" ) -plot(mids, count_II ./ diff((bins)) ./ vol, label="websky FR II" ) - - -plot(log10.(sehgal_red[:, 1]), sehgal_red[:, 2], "rp", - label = "sehgal total", alpha = 0.6) -plot(log10.(sehgal_green[:, 1]), sehgal_green[:, 2], color = "green", - "^", label = "sehgal FR I", alpha = 0.6) -plot(log10.(sehgal_blue[:, 1]), sehgal_blue[:, 2], color = "blue", "s", - label = "sehgal FR II", alpha = 0.6) - -ylabel("counts / log bin / Mpc\$^3\$") -xlabel("log L [W/Hz]") -yscale("log") - -legend() -tight_layout() -savefig("sehgal_LF.pdf") -gcf() - -## - -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 1.4f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_1.4GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 1.4 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_1.4GHz.pdf") -gcf() - -## -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 33f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_33GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 33 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_33GHz.pdf") -gcf() -## -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 95f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_95GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 95 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_95GHz.pdf") -gcf() - -## -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 145f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_145GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 145 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_145GHz.pdf") -gcf() - -## diff --git a/examples/radio_spectra.jl b/examples/radio_spectra.jl deleted file mode 100644 index d4713e9..0000000 --- a/examples/radio_spectra.jl +++ /dev/null @@ -1,6 +0,0 @@ -using XGPaint -using Healpix - -nside = 4096 -m = HealpixMap{Float64,RingOrder}(nside) -XGPaint.catalog2map!(m, [1.0], [0.2], [0.3]) diff --git a/examples/radiocat_to_car.jl b/examples/radiocat_to_car.jl deleted file mode 100644 index 36e0965..0000000 --- a/examples/radiocat_to_car.jl +++ /dev/null @@ -1,110 +0,0 @@ -using HDF5 -using Pixell -using Healpix -using XGPaint - -cat_dir = "/global/project/projectdirs/sobs/www/users/Radio_WebSky/matched_catalogs_2" -# map_dir = "/global/project/projectdirs/sobs/www/users/Radio_WebSky/matched_maps/" -map_dir = "/global/cscratch1/sd/xzackli/radio_maps_will/" - - -shape, wcs = fullsky_geometry(deg2rad(0.5 / 60)) -m = Enmap(zeros(shape), wcs) - - -## -function enmap_py2jl(m_py) - m0 = pyconvert(Array{Float64, 2}, m_py) - permutedims(m0, (2,1)) -end - -using ThreadsX -using PythonCall -enmap = pyimport("pixell.enmap") - -shape_py, wcs_py = enmap.fullsky_geometry(res=deg2rad(0.5 / 60)) -pixsizes_py = enmap.pixsizemap(shape_py, wcs_py) -pixsizes = enmap_py2jl(pixsizes_py) - -## -function catalog2map!(m::Enmap{T}, flux, theta, phi, pixsizes) where T - N_halo = length(flux) - pixel_array = parent(m) - fill!(pixel_array, zero(T)) - α = phi - δ = π/2 .- theta - - ipix, jpix = sky2pix(m, α, δ) - - # try to prevent thread issues by sorting by theta - perm = sortperm(theta, rev=true, alg=ThreadsX.MergeSort) - Threads.@threads :static for i_perm in 1:N_halo - i_halo = perm[i_perm] - i = round(Int, ipix[i_halo]) - i = mod(i - 1, size(m, 1)) + 1 - j = round(Int, jpix[i_halo]) - pixel_array[i, j] += flux[i_halo] - end - - pixel_array ./= pixsizes - - top_pole = sum(pixel_array[:,end]) - bot_pole = sum(pixel_array[:,1]) - - pixel_array[:,end] .= top_pole - pixel_array[:,1] .= bot_pole - return m -end - -## - -# filename = "/global/project/projectdirs/sobs/www/users/Radio_WebSky/matched_catalogs_2/catalog_90.2.h5" -# sources = convert(Dict{String, Vector{Float64}}, read(h5open(filename, "r"))) -# catalog2map!(m, sources["flux"], sources["theta"], sources["phi"], pixsizes) -# write_map(map_path, m) - -## -write_map("test.fits", m) - -## -for (root, dirs, files) in walkdir(cat_dir) - # println("Directories in $root") - # for dir in dirs - # println(joinpath(root, dir)) # path to directories - # end - println("Files in $root") - for file in files - freq = replace(file, ".h5" => "", "catalog_" => "") - filename = joinpath(root, file) - println(filename) - sources = convert(Dict{String, Vector{Float64}}, read(h5open(filename, "r"))) - catalog2map!(m, sources["flux"], sources["theta"], sources["phi"], pixsizes) - map_path = joinpath(map_dir, "map_car_0.5arcmin_f$(freq).fits") - write_map(map_path, m) - end -end - -## - -wy, wx = enmap.calc_window(shape_py) - - - - -# function catalog2map!(m::HealpixMap{T,RingOrder}, flux, theta, phi) where T -# res = m.resolution -# pixel_array = m.pixels -# N_halo = length(flux) - -# # try to prevent thread issues by sorting by theta -# perm = sortperm(theta, rev=true, alg=ThreadsX.MergeSort) -# Threads.@threads :static for i_perm in 1:N_halo -# i_halo = perm[i_perm] -# hp_ind = Healpix.ang2pixRing(res, theta[i_halo], phi[i_halo]) -# pixel_array[hp_ind] += flux[i_halo] -# end - -# # divide by healpix pixel size -# per_pixel_steradian = 1 / nside2pixarea(res.nside) -# pixel_array .*= per_pixel_steradian -# end diff --git a/examples/rectpix.jl b/examples/rectpix.jl deleted file mode 100644 index c33e3b8..0000000 --- a/examples/rectpix.jl +++ /dev/null @@ -1,38 +0,0 @@ -using XGPaint -using PyCall -using PyPlot -using WCS -enmap = pyimport("pixell.enmap") - -wcs2 = WCSTransform(2; - cdelt = [1., 1.], #cdelt = [5/60., 5/60.], - ctype = ["RA---CAR", "DEC--CAR"], - crpix = [181.0, 91.0], - crval = [0., 0.]) -## -l = collect(1:1000) -cl = l .^ -2.5 -shape,wcs = enmap.geometry(shape=(180, 360), - res=deg2rad(1),pos=(0,0)) -modl = enmap.modlmap(shape, wcs); -cmb = enmap.rand_map(shape, wcs, - reshape(XGPaint.ellpad(cl),(1,1,length(cl)+1))); -## -imshow(cmb) -gcf() - -## - -coords = zeros(2,3) -coords = [ - 0.0 90.0 180.0; - 0.0 00.0 00.0 -] -(world_to_pix(wcs2, coords )) - -## -coords = [ - 1.0 5.0 180.0; - 1.0 1.0 1.0 -] -(pix_to_world(wcs2, coords )) diff --git a/examples/validate_radio.jl b/examples/validate_radio.jl deleted file mode 100644 index 0540fdb..0000000 --- a/examples/validate_radio.jl +++ /dev/null @@ -1,178 +0,0 @@ -using XGPaint -using Healpix - -halo_pos, halo_mass = read_halo_catalog_hdf5( - "/home/zequnl/Data/websky_halos-light.hdf5") -## Load halos from HDF5 files, establish a CIB model and cosmology -cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0) -radio_model = Radio_Sehgal2009{Float32}() - -## - -@time begin - sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass); -end; - -## -@time begin -m = HealpixMap{Float64,RingOrder}(radio_model.nside) -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 151f6, radio_model, sources, return_fluxes=true) -end -## - - -using DelimitedFiles -using Unitful -using UnitfulAstro -using Cosmology -using PyPlot -using StatsBase - -figure(figsize=(7,7)) - -sehgal_red = readdlm("examples/data/sehgal_figure8_red.txt", ',', Float64, '\n') -sehgal_blue = readdlm("examples/data/sehgal_figure8_blue.txt", ',', Float64, '\n') -sehgal_green = readdlm("examples/data/sehgal_figure8_green.txt", ',', Float64, '\n') - -vol = ustrip(u"Mpc^3",comoving_volume(cosmo, 0.3)) -bins = (range(23, stop=29, step=0.5)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 -count_I = fit(Histogram, log10.(sources.L_I_151[redshift_I .< 0.3]), bins).weights - -count_II = fit(Histogram, log10.(sources.L_II_151[redshift_II .< 0.3]), bins).weights - -plot(mids, (count_I .+ count_II) ./ diff((bins)) ./ vol, label="websky total" ) -plot(mids, count_I ./ diff((bins)) ./ vol, label="websky FR I" ) -plot(mids, count_II ./ diff((bins)) ./ vol, label="websky FR II" ) - - -plot(log10.(sehgal_red[:, 1]), sehgal_red[:, 2], "rp", - label = "sehgal total", alpha = 0.6) -plot(log10.(sehgal_green[:, 1]), sehgal_green[:, 2], color = "green", - "^", label = "sehgal FR I", alpha = 0.6) -plot(log10.(sehgal_blue[:, 1]), sehgal_blue[:, 2], color = "blue", "s", - label = "sehgal FR II", alpha = 0.6) - -ylabel("counts / log bin / Mpc\$^3\$") -xlabel("log L [W/Hz]") -yscale("log") - -legend() -tight_layout() -savefig("sehgal_LF.pdf") -gcf() - -## - -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 1.4f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_1.4GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 1.4 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_1.4GHz.pdf") -gcf() - -## -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 33f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_33GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 33 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_33GHz.pdf") -gcf() -## -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 95f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_95GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 95 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_95GHz.pdf") -gcf() - -## -plt.figure(figsize=(6,6)) - -bins = (range(-4.1, stop=2.5, step=0.2)); -mids = (bins[2:end] .+ bins[1:end-1]) ./ 2.0 - -flux_I, flux_II, redshift_I, redshift_II = paint!( - m, 145f9, radio_model, sources, return_fluxes=true) -count_I = fit(Histogram, log10.(flux_I), bins).weights -count_II = fit(Histogram, log10.(flux_II), bins).weights - -Sehgal = readdlm("examples/data/sehgal_fig9_145GHz.txt", ',', Float64, '\n'); -plot(Sehgal[:,1], Sehgal[:,2], "ro", label="sehgal 145 GHz") - -plot(10 .^ mids, (count_I .+ count_II) .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "total") -plot(10 .^ mids, count_I .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "I") -plot(10 .^ mids, count_II .* (10 .^ mids) .^ 2.5 ./ diff(10 .^ bins) / (4π), label = "II") - -yscale("log") -xscale("log") -xlabel("\$ S \$ [Jy]") -ylabel("\$ S^{2.5} d\\Sigma / dS \$ [Jy\$^{1.5}\$ sr\$^{-1}\$]") -xticks(10 .^ [-3.,-2.,-1.,0.,1.]) -legend() -tight_layout() -savefig("sehgal_145GHz.pdf") -gcf() - -##