HWK3 November 18, 2019 1 Planck 1-d circular slices data project Note: I’m using weave to make a notebook from a script file. using Weave convert_doc(“notebook2.jl”, “notebook2.md”) # for markdown format convert_doc(“notebook2.jl”, “notebook2.ipynb”) # for jupyter notebook format Or directly from shell julia -e “using Weave; convert_doc(\”notebook2.jl\”, \”notebook2.ipynb\”)” julia -e “using Weave; convert_doc(\”notebook2.jl\”, \”notebook2.md\”)” 1.1 Start by defining a module to hold some constants and set the path to the data directory. [1]: module LocalMethods using GSL using LinearAlgebra using Statistics import JLD2 export cov_n1_dot_n2, cov_angle const data_dir = “/Users/ethananderes/Dropbox/Courses/STA250fall2019/Datasets/ ↪→Dataset-2-Planck_1d_transect/data/” const cls = JLD2.FileIO.load(data_dir*”cls.jld2″, “cls”) const w_on_Pl = cls[:cl_len_scalar][:,1] .* (2 .* cls[:ell] .+ 1) ./ (4) ./␣ ↪→cls[:factor_on_cl_cmb] const lmax = maximum(cls[:ell]) function cov_cos(cos,lmax::Int) @assert lmax >= 2 “lmax too small, must be greater than 1” dot(sf_legendre_Pl_array(lmax, cos)[3:end], w_on_Pl[3:lmax+1]) end cov_n1_dot_n2(n1_dot_n2) = cov_cos(n1_dot_n2,lmax) 1 cov_angle(angle) = cov_cos(cos(angle),lmax) function bin_mean(fk; bin=0) fkrtn = copy(fk) if bin > 1 Nm1 = length(fk) subNm1 = bin * (Nm1÷bin) fmk = reshape(fk[1:subNm1], bin, Nm1÷bin) fmk .= mean(fmk, dims=1) fkrtn[1:subNm1] .= vec(fmk) fkrtn[subNm1+1:end] .= mean(fk[subNm1+1:end]) end return fkrtn end end # module [1]: Main.LocalMethods 1.2 Load modules [2]: using PyPlot using PyCall using LinearAlgebra using Statistics using Random import JLD2 using FFTW using SparseArrays using ProgressMeter using .LocalMethods const prop_cycle = plt.rcParams[“axes.prop_cycle”] const colors = prop_cycle.by_key()[“color”] const mplot3d = pyimport(“mpl_toolkits.mplot3d”) [2]: PyObject ‘/Users/ethananderes/Software/anaconda3/lib/python3.7/site- packages/mpl_toolkits/mplot3d/__init__.py’> 1.3 Load the data The data loaded below contains 10 slices of the CMB at 10 different polar angles. 5 slices from the north pole, 5 slices from the south pole. The polar angles from which these slices are taken are given in variables _northp and _southp. The actual CMB measurements at each slice is stored 2 in variables cmb_northp and cmb_southp which are simply vectors the indices of which map to azmuthal coordinates stored in variable _ref. These slices are all taken at polar angles which have the same azmuthal coordinates, i.e. each CMB slice is the same length (denoted N below) the entries of which have asmuthal coordinates given in _ref. [3]: load_it = (jld2file, varname) -> JLD2.FileIO.load(LocalMethods. ↪→data_dir*jld2file, varname) _northp = load_it(“data_slices.jld2”, “_northp”) _southp = load_it(“data_slices.jld2”, “_southp”) cmb_northp = load_it(“data_slices.jld2”, “cmb_northp”) cmb_southp = load_it(“data_slices.jld2”, “cmb_southp”) _ref = load_it(“data_slices.jld2”, “_ref”) N = length(_ref) Nside = 2048 [3]: 2048 For the homework assignemnt, only analyze the first and fourth CMB slice from the south pole. Here is the code for extracting them. [4]: I_slice_1 = 1 I_slice_2 = 4 1 = _southp[I_slice_1] 2 = _southp[I_slice_2] cmb1 = cmb_southp[I_slice_1] cmb2 = cmb_southp[I_slice_2] @show 1, 2 @show summary(cmb1) @show summary(cmb2) @show summary(_ref); (1, 2) = (1.990924632438072, 2.2330667085254525) summary(cmb1) = “8192-element Array{Float64,1}” summary(cmb2) = “8192-element Array{Float64,1}” summary(_ref) = “8192-element Array{Float64,1}” [5]: ofset_cmb = 9000 .* vcat(_northp, _southp) ./ .- 0.5 |> reverse figure(figsize=(10,5)) for i = 1:length(_northp) plot(_ref, cmb_northp[i] .+ ofset_cmb[i], “k”, alpha = 0.1) end for i = 1:length(_southp) if i in [I_slice_1, I_slice_2] plot(_ref, cmb_southp[i] .+ ofset_cmb[i+length(_southp)]) else 3 plot(_ref, cmb_southp[i] .+ ofset_cmb[i+length(_southp)], “k”,␣ ↪→alpha = 0.1) end end # Note: The cmb is in units mK MicroKelvin (1 Kelvin = 1e+6 MicroKelvin) 1.4 Digression: Take a peak at the full sky data The full sky data file is also uploaded to Canvas (file data/data_full.jld2) but you don’t need it unless you want to look at it. In which case you will need healpy to make the following few plots. (install healpy at the command line with something like pip install –user healpy) [6]: const hp = pyimport(“healpy”) const ScS = pyimport(“scipy.spatial”) const ScI = pyimport(“scipy.interpolate”) cmb_full = load_it(“data_full.jld2”, “cmb_full”) full = load_it(“data_full.jld2”, “full”) full = load_it(“data_full.jld2”, “full”) hp.visufunc.mollview(cmb_full, xsize=2048, title=”full cmb map”) hp.visufunc.mollview(cmb_full, min=-250, max=500, xsize=Nside, title=”full cmb␣ ↪→map with restricted color scale”) 4 To get a visual where the CMB slices used in the homework come from, they are take at polar angles from within these strips. 5 [7]: custom_mask = (_northp[1] .<= full .<= _northp[end]) .| (_southp[1] .<=␣ ↪→full .<= _southp[end]) hp.visufunc.mollview(custom_mask .* cmb_full, min=-250, max=500, xsize=Nside) 1.5 1-st column of covariance covariance and cross-covariance matrix The isotropic nature of the CMB means we can compute the covariance between the pixel obser- vations as a function of the inner product of the three dimensional sphere coordinates. Here is a function that “eats polar coordaintes” and returns the tuple of three dim spherical coordinates. [8]: function n̂(,) nx = sin.() .* cos.() ny = sin.() .* sin.() nz = cos.() return (nx, ny, nz) end @show n̂(1, _ref[1]); @show n̂(2, _ref[1]); n̂(1, _ref[1]) = (0.9130366148294711, 0.0, -0.4078776041666667) n̂(2, _ref[1]) = (0.7885981873346127, 0.0, -0.6149088541666667) Now we can compute the inner product between a fixed pixel and all the rest on the same slice as follows 6 [9]: # we use `Ref` here since `n̂(1, _ref[1])` is a single 3-tuple that we want # broadcasted to the vector of 3-tuples `n̂.(1, _ref)` dot.(n̂.(1, _ref), Ref(n̂(1, _ref[1]))) [9]: 8192-element Array{Float64,1}: 1.0 0.9999997547967509 0.9999990191871476 0.9999977931716232 0.9999960767508986 0.9999938699259837 0.9999911726981767 0.9999879850690642 0.9999843070405215 0.9999801386147124 0.9999754797940889 0.9999703305813916 0.9999646909796499 0.9999646909796499 0.9999703305813916 0.9999754797940889 0.9999801386147124 0.9999843070405215 0.9999879850690642 0.9999911726981767 0.9999938699259837 0.9999960767508986 0.9999977931716232 0.9999990191871476 0.9999997547967509 Now we can use the methods given in LocalMethods (which “eats” inner products and returns cov) to compute the first column of the within group and accross group covariance matrix. [10]: # slice 1 within group covariance to first pixel Σ11 = dot.(n.̂(1, _ref), Ref(n̂(1, _ref[1]))) .|> cov_n1_dot_n2 # slice 2 within group covariance to first pixel Σ22 = dot.(n.̂(2, _ref), Ref(n̂(2, _ref[1]))) .|> cov_n1_dot_n2 # slice 2 / slice 1 across group covariance to first pixel Σ12 = dot.(n.̂(1, _ref), Ref(n̂(2, _ref[1]))) .|> cov_n1_dot_n2 # slice 1 / slice 2 across group covariance to first pixel Σ21 = dot.(n.̂(2, _ref), Ref(n̂(1, _ref[1]))) .|> cov_n1_dot_n2 [10]: 8192-element Array{Float64,1}: 748.7111817196152 7 748.708561033119 748.70069908468 748.6875962074554 748.66925295663 748.6456701093434 748.6168486645637 748.5827898430257 748.5434950871988 748.4989660613383 748.4492046516665 748.394212966757 748.3339933380622 748.3339933380622 748.394212966757 748.4492046516665 748.4989660613383 748.5434950871988 748.5827898430257 748.6168486645637 748.6456701093434 748.66925295663 748.6875962074554 748.70069908468 748.708561033119 [11]: figure(figsize=(10,4)) subplot(1,2,1) plot(_ref[1:200], Σ11[1:200], label=”slice 1″) plot(_ref[1:200], Σ22[1:200], “–“, label=”slice 2”) xlabel(“azmuth lag”) title(“within slice covariance”) legend() subplot(1,2,2) plot(_ref[1:200], Σ21[1:200], label=”slice 1″) plot(_ref[1:200], Σ12[1:200], “–“, label=”slice 1”) xlabel(“azmuth”) title(“across slice cross-covariance”) legend() 8 [11]: PyObject 1.6 Spectral analysis Using the fact that these 4 block matrices are ci
rculant, we take the (scaled) Discrete Fourier Transform to diagonalize them [12]: # Here are the pixel and grid parameters period = 2 Δ = period / N # pixel spacing Δ = 2 / period # Fourier spacing nyq = / Δ _ref = Δ .* (0:N-1) |> collect _nyq = @. ifelse(_ref <= nyq, _ref, _ref - 2nyq) [12]: 8192-element Array{Float64,1}: 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 9 -12.0 -11.0 -10.0 -9.0 -8.0 -7.0 -6.0 -5.0 -4.0 -3.0 -2.0 -1.0 [13]: # Here are the fft operators = plan_fft(similar(cmb1)) # unscaled fft = (1/sqrt(N)) * # unitary version = (Δ / √(2)) * # scaled fft which approximates the integral sc = (√(2)/Δ) * # scaling used for the eigenvalues from the first␣ ↪→col of the circulantmatrix [13]: 0.0007669903939428206 * FFTW forward plan for 8192-element array of Complex{Float64} [14]: # These model *cmb1 and *cmb2 scΣ11 = sc * Σ11 .|> real scΣ22 = sc * Σ22 .|> real scΣ21 = sc * Σ21 scΣ12 = sc * Σ12 [14]: 8192-element Array{Complex{Float64},1}: 154.62796786456275 + 2.786085204295846e-16im 599.569757606981 – 1.6779814513390042e-13im 779.6420483508 – 4.583568012836338e-13im 372.8710732243214 – 3.213635898096056e-13im 200.89061308534298 – 2.432443930829525e-13im 116.5589092472053 – 1.606333851540168e-13im 70.9989149940158 – 1.3088744033495058e-13im 44.79815710632035 – 7.668671487384125e-14im 28.98581845967752 – 6.384661237450354e-14im 19.12080169053594 – 4.618813712709086e-14im 12.844497256954483 – 4.475246054726254e-14im 8.730203340362827 – 3.0643558407250344e-14im 5.974482487946962 – 1.7940551232936005e-14im 5.974482487946963 + 1.9504201137134958e-14im 8.730203340362824 + 3.951261946811799e-14im 10 12.844497256954469 + 3.954397838664373e-14im 19.12080169053594 + 4.8859926922865755e-14im 28.98581845967752 + 6.384412657781292e-14im 44.79815710632037 + 7.424454965801742e-14im 70.9989149940158 + 1.2573664357393976e-13im 116.55890924720532 + 1.6118877847732418e-13im 200.890613085343 + 2.435562912903662e-13im 372.8710732243214 + 3.267591050513847e-13im 779.6420483508 + 4.569026892541495e-13im 599.5697576069812 + 1.681117049520985e-13im The theory of circulant matrices which predicts (√(2)/Δ) * Σ_ij is the expected power in *cmb_i * conj(*cmbj). [15]: cmb1 = * cmb1 cmb2 = * cmb2 [15]: 8192-element Array{Complex{Float64},1}: 133.90814399838476 + 2.1225283178226174e-16im 4.598592707638935 – 14.63128249815512im 11.703730351827819 + 19.90828970987386im -37.33998677029809 – 8.809719862382574im 37.647894786328564 – 24.40948634772287im -27.78505691144916 – 13.20704320742245im 3.7853845896430225 + 6.4362876647001475im 7.8393694627857045 – 11.794714499160008im 13.234914335385948 – 15.317386211599176im 1.6906889666819944 + 12.45987370140233im 8.969810355925363 – 7.910895730786762im 18.591384350538966 + 11.662364168115035im 17.818932365135325 – 23.169844939092187im 17.818932365135325 + 23.169844939092187im 18.59138435053897 – 11.662364168115039im 8.969810355925361 + 7.910895730786761im 1.6906889666819942 – 12.45987370140233im 13.234914335385948 + 15.31738621159918im 7.8393694627857045 + 11.794714499160008im 3.7853845896430203 – 6.436287664700148im -27.78505691144916 + 13.207043207422455im 37.647894786328564 + 24.40948634772287im -37.33998677029809 + 8.809719862382572im 11.703730351827819 – 19.90828970987386im 4.598592707638937 + 14.631282498155121im We can look at the “power” in these maps, i.e. |cmb1[k]|^2 and |cmb2[k]|^2, compared to what the signal covariance predicts 11 [16]: let plt_range = 10:N÷2+1, k = _ref[plt_range], km = k.^2, abs2cmb1 = abs2.(cmb1), abs2cmb2 = abs2.(cmb2), scΣ11 = scΣ11, scΣ22 = scΣ22 for myplot [semilogx, loglog] figure(figsize=(10,4)) sub = subplot(1,2,1) myplot(k, km .* abs2cmb1[plt_range], label=”abs2cmb1″,␣ ↪→color=colors[1], “.”, alpha=0.15) myplot(k, km .* scΣ11[plt_range], label=”scΣ11″,␣ ↪→color=colors[2]) xlabel(“wavenumber”) sub.set_ylim((km .* scΣ11[plt_range])[1], 2*maximum(km␣ ↪→.* scΣ11[plt_range])) sub.set_xlim(k[1], maximum(k)) legend() sub = subplot(1,2,2) myplot(k, km .* abs2cmb2[plt_range], label=”abs2cmb2″,␣ ↪→color=colors[1], “.”, alpha=0.15) myplot(k, km .* scΣ22[plt_range], label=”scΣ22″,␣ ↪→color=colors[2]) xlabel(“wavenumber”) sub.set_ylim((km .* scΣ22[plt_range])[1], 2*maximum(km␣ ↪→.* scΣ22[plt_range])) sub.set_xlim(k[1], maximum(k)) legend() end # for loop over myplot end 12 1.7 Spectral impact of the noise and beam From https://wiki.cosmos.esa.int/planckpla2015/index.php/Summary_of_HFI_data_characteristics Effective beam solid angle Ω [arcmin²] 60.44 Effecitive beam FWHM [arcmin] 7.3 Sensitivity per beam solid angle [K_CMB] 4.3 Temp sensitivity [Kdeg] 0.55 [17]: # Compute the standard deviation of white noise in each Healpix pixel. # Note: the Healpix pixel area 1.7² [arcmin²] with Nside = 2048. 13 # ############# # # # For the homework … assign the variable `sd_of_noise_per_1p7_arcmin_pixel`␣ ↪→the correct value. # Note: the Healpix pixel area 1.7² [arcmin²] with Nside = 2048. # # # ############## sd_of_noise_per_1p7_arcmin_pixel = ?????????? bFWHM_arcmin = 7.3 # arcmin [17]: 7.3 [18]: , = let n = sd_of_noise_per_1p7_arcmin_pixel, bfwhm = bFWHM_arcmin = fill(n^2, size(_ref)) = (Δ / Δ) .* ### note the scaling bFWHM_rad = deg2rad(bfwhm/60) = @. exp( – (bFWHM_rad^2) * abs2(_nyq) / (16*log(2)) ) , end [18]: ([0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252 … 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252], [1.0, 0.9999995934139979, 0.9999983736569836, 0.9999963407319328, 0.9999934946438046, 0.9999898353995421, 0.9999853630080723, 0.9999800774803053, 0.9999739788291351, 0.9999670670694394 … 0.999959342218079, 0.9999670670694394, 0.9999739788291351, 0.9999800774803053, 0.9999853630080723, 0.9999898353995421, 0.9999934946438046, 0.9999963407319328, 0.9999983736569836, 0.9999995934139979]) Now we can plot the theory beamed signal, the noise and binned bandpowers [19]: let bin_powers_width = 12, plt_range = 10:N÷2+1, k = _ref[plt_range], km = k.^2, abs2cmb1 = abs2.(cmb1), abs2cmb2 = abs2.(cmb2), scΣ11 = scΣ11, scΣ22 = scΣ22 14 beam_scΣ11 = scΣ11 .* abs2.() beam_scΣ22 = scΣ22 .* abs2.() tot_scΣ11 = scΣ11 .* abs2.() .+ tot_scΣ22 = scΣ22 .* abs2.() .+ for myplot [semilogx, loglog] figure(figsize=(10,4)) sub = subplot(1,2,1) myplot(k, km .* abs2cmb1[plt_range] |> x->LocalMethods. ↪→bin_mean(x; bin = bin_powers_width), label=”binned abs2cmb1″, color=colors[1], “.”,␣ ↪→alpha=0.15 ) myplot(k, km .* beam_scΣ11[plt_range],␣ ↪→label=”beam_scΣ11″, color=colors[2]) myplot(k, km .* [plt_range], label=”N”, “k”, alpha␣ ↪→= 0.8) myplot(k, km .* tot_scΣ11[plt_range], “–“,␣ ↪→label=”tot_scΣ11”, color=colors[7]) xlabel(“wavenumber”) sub.set_ylim((km .* beam_scΣ11[plt_range])[1],␣ ↪→2*maximum(km .* beam_scΣ11[plt_range])) sub.set_xlim(k[1], maximum(k)) legend() sub = subplot(1,2,2) myplot(k, km .* abs2cmb2[plt_range] |> x->LocalMethods. ↪→bin_mean(x; bin = bin_powers_width), label=”binned abs2cmb2″, color=colors[1], “.”,␣ ↪→alpha=0.15 ) myplot(k, km .* beam_scΣ22[plt_range],␣ ↪→label=”beam_scΣ22″, color=colors[2]) myplot(k, km .* [plt_range], label=”N”, “k”, alpha␣ ↪→= 0.8) myplot(k, km .* tot_scΣ22[plt_range], “–“,␣ ↪→label=”tot_scΣ22”, color=colors[7]) xlabel(“wavenumber”) sub.set_ylim((km .* beam_scΣ22[plt_range])[1],␣ ↪→2*maximum(km .* beam_scΣ22[plt_range])) sub.set_xlim(k[1], maximum(k)) legend() 15 end # for loop over myplot end Also the cross correlation. The circulant matrix theory predicts that E(cmb1 .* conj(cmb2)) = scΣ12 .* abs2.() and E(cmb2 .* conj(cmb1)) = scΣ21 .* abs2.() [20]: let bin_powers_width = 0, plt_range = 2:N÷2+1, k = _ref[plt_range], km = k.^0, crosscmb12 = cmb1 .* conj(cmb2), crosscmb21 = cmb2 .* conj(cmb1), 16 scΣ12 = scΣ12, scΣ21 = scΣ21 r_crosscmb = (crosscmb12
.+ crosscmb21) ./ 2 .|> real # this must be␣ ↪→real i_crosscmb = (crosscmb12 .- crosscmb21) ./ 2 .|> imag # this must be␣ ↪→purely imaginary r_beam_cross_scΣ = real.((scΣ12 .+ scΣ21) ./ 2) .* abs2.() i_beam_cross_scΣ = imag.((scΣ12 .- scΣ21) ./ 2) .* abs2.() myplot = semilogx #plot figure(figsize=(10,4)) sub = subplot(1,2,1) myplot(k, km .* r_crosscmb[plt_range] |> x->LocalMethods. ↪→bin_mean(x; bin = bin_powers_width), label=”binned r_crosscmb”, color=colors[1], “.”,␣ ↪→alpha=0.15 ) myplot(k, km .* r_beam_cross_scΣ[plt_range],␣ ↪→label=”r_beam_cross_scΣ”, color=colors[2]) xlabel(“wavenumber”) #sub.set_ylim((km .* r_beam_cross_scΣ[plt_range])[1],␣ ↪→2*maximum(km .* r_beam_cross_scΣ[plt_range])) sub.set_xlim(k[2], maximum(k)) legend() sub = subplot(1,2,2) myplot(k, km .* i_crosscmb[plt_range] |> x->LocalMethods. ↪→bin_mean(x; bin = bin_powers_width), label=”binned i_crosscmb”, color=colors[1], “.”,␣ ↪→alpha=0.15 ) myplot(k, km .* i_beam_cross_scΣ[plt_range],␣ ↪→label=”i_beam_cross_scΣ”, color=colors[2]) xlabel(“wavenumber”) #sub.set_ylim((km .* i_beam_cross_scΣ[plt_range])[1],␣ ↪→2*maximum(km .* i_beam_cross_scΣ[plt_range])) sub.set_xlim(k[2], maximum(k)) legend() end 17 [20]: PyObject 1.8 Frequency-by-frequency de-correlation of cmb1 and cmb2 [21]: tot_scΣ11 = scΣ11 .* abs2.() .+ tot_scΣ22 = scΣ22 .* abs2.() .+ tot_scΣ21 = scΣ21 .* abs2.() .|> real = tot_scΣ21 ./ sqrt.( tot_scΣ11 .* tot_scΣ22) nrm_cmb1 = cmb1 ./ sqrt.(tot_scΣ11) nrm_cmb2 = cmb2 ./ sqrt.(tot_scΣ22) z1 = nrm_cmb1 .+ exp.(im.*angle.()) .* nrm_cmb2 z2 = nrm_cmb1 .- exp.(im.*angle.()) .* nrm_cmb2 1 = 1 .- abs.() 2 = 1 .+ abs.() let bin_powers_width = 4, plt_range = 10:N÷2+1, k = _ref[plt_range], abs2_z1 = abs2.(z1), abs2_z2 = abs2.(z2), cross_z2_z1 = z2 .* conj.(z1), 1 = 1, 2 = 2 figure(figsize=(10,4)) sub = subplot(1,2,1) plot( 18 k, abs2_z1[plt_range] |> x->LocalMethods.bin_mean(x; bin =␣ ↪→bin_powers_width), label=”binned abs2(z1)”, color=colors[1], “.”, alpha=0. ↪→15 ) plot(k, 1[plt_range], label=”1″, color=colors[2]) sub = subplot(1,2,2) plot( k, abs2_z2[plt_range] |> x->LocalMethods.bin_mean(x; bin =␣ ↪→bin_powers_width), label=”binned abs2(z2)”, color=colors[1], “.”, alpha=0. ↪→15 ) plot(k, 2[plt_range], label=”2″, color=colors[2]) legend() figure(figsize=(10,4)) plot( k, real.(cross_z2_z1)[plt_range] |> x->LocalMethods. ↪→bin_mean(x; bin = bin_powers_width), label=”real(z1 * conj(z2))”, color=colors[1], “.”,␣ ↪→alpha=0.15 ) plot( k, imag.(cross_z2_z1)[plt_range] |> x->LocalMethods. ↪→bin_mean(x; bin = bin_powers_width), label=”imag(z1 * conj(z2))”, color=colors[2], “.”,␣ ↪→alpha=0.15 ) plot(k, fill(0, size(k)), “k”, alpha = 0.8) legend() end 19 [21]: PyObject 1.9 Wiener filtered south pole cap First we multiply the data (i.e. the two circular slices stacked) by the inverse cov matrix of the data. Now we make a function that computes the wiener filter on any unobserved circular slice. [22]: # ############# # # # For the homework … # Create a function `wf_cmb_fun()` that takes in a polar angle “ and 20 # returns the conditional expected value of `cmb` at the azmuth grid values␣ ↪→`_ref`. # # # ############## ???????????????? [22]: wf_cmb_fun (generic function with 1 method) Finally we create a grid of polar angles to predict on [23]: # if at north pole #s = range(0.001, max(2, 1) .+ 0.1, length=250) |> transpose # if at south pole s = range(min(2, 1) .- 0.1, , length=150) |> transpose grid = s .+ 0 .* _ref grid = 0 .* s .+ _ref function make_wf_cmb_mat(s, _ref) wf_cmb_mat_local = 0 .* s .+ 0 .* _ref @showprogress for icol = 1:length(s) #print(“$icol,”) wf_cmb_mat_local[:,icol] = wf_cmb_fun(s[icol]) end wf_cmb_mat_local end wf_cmb_mat = make_wf_cmb_mat(s, _ref) Progress: 100%|| Time: 0:06:09 [23]: 8192×150 Array{Float64,2}: 41.6078 43.7738 45.9311 48.05 50.0846 … -108.516 -108.587 -108.627 41.9043 44.0884 46.2649 48.4048 50.4616 -108.516 -108.587 -108.627 42.2009 44.4032 46.5991 48.7604 50.8397 -108.516 -108.587 -108.627 42.4976 44.7181 46.9337 49.1168 51.219 -108.515 -108.587 -108.627 42.7944 45.0332 47.2689 49.4739 51.5995 -108.515 -108.587 -108.627 43.0912 45.3486 47.6045 49.8318 51.9812 … -108.515 -108.586 -108.627 43.3882 45.664 47.9406 50.1903 52.3642 -108.515 -108.586 -108.627 43.6851 45.9797 48.2771 50.5495 52.7483 -108.515 -108.586 -108.627 43.9821 46.2956 48.6141 50.9095 53.1337 -108.515 -108.586 -108.627 44.279 46.6116 48.9516 51.2701 53.5204 -108.514 -108.586 -108.627 44.5759 46.9279 49.2894 51.6315 53.9081 … -108.514 -108.586 -108.627 21 44.8728 47.2443 49.6277 51.9935 54.297 -108.514 -108.586 -108.627 45.1696 47.5609 49.9663 52.3562 54.6871 -108.514 -108.586 -108.627 38.0591 40.0195 41.9633 43.8583 45.6599 … -108.518 -108.588 -108.627 38.3541 40.3307 42.2911 44.2028 46.0215 -108.518 -108.588 -108.627 38.6492 40.6422 42.6193 44.5481 46.3844 -108.518 -108.588 -108.627 38.9444 40.9541 42.9482 44.8944 46.7486 -108.517 -108.588 -108.627 39.2398 41.2662 43.2775 45.2415 47.1141 -108.517 -108.588 -108.627 39.5353 41.5787 43.6074 45.5895 47.481 … -108.517 -108.588 -108.627 39.831 41.8915 43.9378 45.9384 47.8491 -108.517 -108.587 -108.627 40.1268 42.2046 44.2688 46.2882 48.2185 -108.517 -108.587 -108.627 40.4228 42.5179 44.6003 46.6389 48.5892 -108.517 -108.587 -108.627 40.7189 42.8316 44.9323 46.9904 48.9612 -108.516 -108.587 -108.627 41.0151 43.1454 45.2647 47.3428 49.3344 … -108.516 -108.587 -108.627 41.3114 43.4595 45.5977 47.696 49.7089 -108.516 -108.587 -108.627 [24]: fig = figure() ax = mplot3d.Axes3D(fig) subplot(projection=”polar”) # if at north pole # pcolormesh(grid, grid, wf_cmb_mat, vmin=-125, vmax=225) # if at south pole pcolormesh(grid, .- grid, wf_cmb_mat, vmin=-135, vmax=235) # grid() 22 [24]: PyObject 23
