辅导案例-HWK3

欢迎使用51辅导,51作业君孵化低价透明的学长辅导平台,服务保持优质,平均费用压低50%以上! 51fudao.top
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 circulant, 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
51作业君

Email:51zuoyejun

@gmail.com

添加客服微信: abby12468