-
Notifications
You must be signed in to change notification settings - Fork 0
/
Plot_Residuals_Sky.py
90 lines (61 loc) · 3.09 KB
/
Plot_Residuals_Sky.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
import numpy
import argparse
from matplotlib import colors
from src.covariance import sky_covariance
from src.covariance import beam_covariance
from src.powerspectrum import compute_power
from src.powerspectrum import from_frequency_to_eta
from src.plottools import plot_2dpower_spectrum
def main(labelfontsize = 16, ticksize= 11):
k_perp_range = numpy.array([1e-4, 1.1e-1])
u_range = numpy.logspace(-1, numpy.log10(500), 100)
frequency_range = numpy.linspace(135, 165, 251) * 1e6
eta = from_frequency_to_eta(frequency_range)
position_error_power = calculate_sky_power_spectrum(u=u_range, nu=frequency_range)
beam_error_power = calculate_beam_power_spectrum(u=u_range, nu=frequency_range)
total_error_power = calculate_total_power_spectrum(u=u_range, nu=frequency_range)
figure, axes = pyplot.subplots(1, 3, figsize=(15, 5))
ps_norm = colors.LogNorm(vmin=1e3, vmax=1e15)
plot_2dpower_spectrum(u_range, eta, frequency_range, position_error_power, title="Sky Model Error", axes=axes[0],
axes_label_font=labelfontsize, tickfontsize=ticksize, colorbar_show=True,
xlabel_show=True, norm=ps_norm, ylabel_show=True)
plot_2dpower_spectrum(u_range, eta, frequency_range, beam_error_power, title="Beam Model Error", axes=axes[1],
axes_label_font=labelfontsize, tickfontsize=ticksize, colorbar_show=True,
xlabel_show=True, norm=ps_norm)
plot_2dpower_spectrum(u_range, eta, frequency_range, total_error_power, title="Total Modelling Error", axes=axes[2],
axes_label_font=labelfontsize, tickfontsize=ticksize, colorbar_show=True,
xlabel_show=True, norm=ps_norm, zlabel_show=True)
figure.tight_layout()
pyplot.show()
return
def calculate_sky_power_spectrum(u, nu):
variance = numpy.zeros((len(u), int(len(nu) / 2)))
print(f"Calculating covariances for all baselines")
for i in range(len(u)):
nu_cov = sky_covariance(u[i], v=0, nu=nu)
variance[i, :] = compute_power(nu, nu_cov)
return variance
def calculate_beam_power_spectrum(u, nu):
variance = numpy.zeros((len(u), int(len(nu) / 2)))
print(f"Calculating covariances for all baselines")
for i in range(len(u)):
nu_cov = beam_covariance(u[i], v=0, nu=nu, calibration_type='sky')
variance[i, :] = compute_power(nu, nu_cov)
return variance
def calculate_total_power_spectrum(u, nu):
variance = numpy.zeros((len(u), int(len(nu) / 2)))
print(f"Calculating covariances for all baselines")
for i in range(len(u)):
nu_cov = beam_covariance(u= u[i], v=0, nu=nu, calibration_type="sky") + \
sky_covariance(u=u[i], v=0, nu=nu)
variance[i, :] = compute_power(nu, nu_cov)
return variance
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--ssh", action="store_true", dest="ssh_key", default=False)
params = parser.parse_args()
import matplotlib
if params.ssh_key:
matplotlib.use("Agg")
from matplotlib import pyplot
main()