-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_diagnostics.py
94 lines (74 loc) · 2.26 KB
/
plot_diagnostics.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
91
92
93
94
import numpy as n
import matplotlib.pyplot as plt
import h5py
import sys
import glob
import scipy.signal as ss
fl=glob.glob("%s/lpi-*.h5"%(sys.argv[1]))
fl.sort()
n_t=len(fl)
h=h5py.File(fl[0],"r")
S=h["diagnostic_pwr_spec"][()]
R=h["retained_measurement_fraction"][()]
meas_delay=h["meas_delays_us"][()]
#plt.plot(meas_delay,R)
#plt.show()
dop=n.fft.fftshift(n.fft.fftfreq(1024,d=1/1e6))
dop_idx=n.where(n.abs(dop)<500e3)[0]
#plt.plot(dop,10.0*n.log10(S))
#plt.show()
n_fft=len(dop_idx)
h.close()
PS=n.zeros([n_t,n_fft],dtype=n.float32)
PS[:,:]=n.nan
RF=n.zeros([n_t,len(R)])
tv=n.zeros(n_t)
alpha=n.zeros(n_t)
ptx=n.zeros(n_t)
T_sys=n.zeros(n_t)
for fi,f in enumerate(fl):
h=h5py.File(f,"r")
if "retained_measurement_fraction" in h.keys():
RF[fi,:]=h["retained_measurement_fraction"][()]
if "T_sys" in h.keys():
T_sys[fi]=h["T_sys"][()]
if "P_tx" in h.keys():
ptx[fi]=h["P_tx"][()]
if "alpha" in h.keys():
alpha[fi]=h["alpha"][()]
if "diagnostic_pwr_spec" in h.keys():
# PS[fi,:]=n.convolve(h["diagnostic_pwr_spec"][()],n.repeat(1/5,5),mode="same")[dop_idx]
PS[fi,:]=h["diagnostic_pwr_spec"][()][dop_idx]#/alpha[fi]
PS[fi,:]=PS[fi,:]/n.nanmedian(PS[fi,:])
# PS[fi,:]=PS[fi,:]/n.nanmedian(PS[fi,:])
tv[fi]=h["i0"][()]
h.close()
#for fi in range(PS.shape[1]):
# PS[:,fi]=n.fft.ifft(n.fft.fft(PS[:,fi])*n.fft.fft(n.repeat(1/24,24),len(PS[:,fi]))).real
plt.subplot(311)
dB=10.0*n.log10(PS.T)
nf=n.nanmedian(dB)
print(nf)
plt.pcolormesh(tv,dop[dop_idx]/1e3,dB,vmin=nf-2,vmax=nf+2,cmap="plasma")
plt.title(sys.argv[1])
plt.xlabel("Time (unix)")
plt.ylabel("Frequency offset from 440.2 MHz (kHz)")
cb=plt.colorbar()
cb.set_label("Power at 1000 km (dB)")
plt.subplot(312)
plt.pcolormesh(tv,meas_delay*0.15,RF.T,vmin=0.7,vmax=1,cmap="jet")
plt.ylabel("Delay (km)")
plt.xlabel("Time (unix)")
cb=plt.colorbar()
cb.set_label("Retained measurement fraction")
plt.subplot(313)
plt.plot(tv,T_sys,label="T_sys")
plt.plot(tv,ptx/1000,label="P_tx/1000")
print(alpha)
plt.plot(tv,1000.0*(alpha-n.nanmedian(alpha))/n.nanmedian(alpha)+500,label="Receiver gain x 100")
plt.legend()
plt.ylim([0,2100])
plt.xlabel("Time (unix)")
#plt.ylabel("T$_{\\mathrm{sys}}$ (K)")
plt.tight_layout()
plt.show()