forked from rvweeren/lofar_facet_selfcal
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplot_tecandphase.py
112 lines (89 loc) · 3.36 KB
/
plot_tecandphase.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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#!/usr/bin/env python
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import os, sys
import numpy as np
import argparse
import tables
def phaseplot(phase):
phase = np.mod(phase+np.pi, 2.*np.pi) - np.pi
return phase
def TEC2phase(TEC,freq):
phase = -1.*(TEC*8.44797245e9/freq)
return phase
parser = argparse.ArgumentParser(description='Plot DPPP tecandphase solutions by combining the TEC and phase offset')
#imaging options
parser.add_argument('-h5','--H5file', help='H5 solution file, must contain tecandphase solutions', type=str, required=True)
parser.add_argument('-o','--outfile', help='Output figure name (png format)', type=str, required=True)
parser.add_argument('-p','--plotpoints', help='Plot points instead of lines', action='store_true')
#parser.add_argument('-f','--freq', help='Frequnecy at which the phase corrections are plotted', type=float, required=True)
args = vars(parser.parse_args())
H=tables.open_file(args['H5file'], mode='r')
try:
tec = H.root.sol000.tec000.val[:]
antennas = H.root.sol000.tec000.ant[:].tolist()
times= H.root.sol000.tec000.time[:]
freq = H.root.sol000.tec000.freq[:][0]
notec = False
except:
print('Your solutions contain do NOT contain TEC values')
#sys.exit()
notec = True
try:
phase = H.root.sol000.phase000.val[:]
antennas = H.root.sol000.phase000.ant[:].tolist()
times= H.root.sol000.phase000.time[:]
freq = H.root.sol000.phase000.freq[:][0]
containsphase = True
if notec:
freq = H.root.sol000.phase000.freq[:]
except:
print('Your solutions contain do NOT contain phase values')
containsphase = False
pass
print('Plotting at a frequency of:', freq/1e6, 'MHz')
times = (times-np.min(times))/3600. # in hrs since obs start
#print(tec.shape, phase.shape)
ysizefig = np.float(len(antennas))
refant = 0
antennas=[x.decode('utf-8') for x in antennas] # convert to proper string list
if 'ST001' in antennas:
refant = antennas.index('ST001')
print('ST001 detected, switching to refant:', refant)
fig, ax = plt.subplots(nrows=len(antennas), ncols=1, figsize=(9, 1.5*ysizefig), squeeze=True, sharex='col')
figcount = 0
if containsphase:
if notec:
freqidx = np.int(len(freq)/2)
refphase = phase[:,freqidx,refant,0,0]
else:
refphase = phase[:,refant,0,0] + TEC2phase(tec[:,refant,0,0], freq)
else:
refphase = TEC2phase(tec[:,refant,0,0], freq)
#print('Here')
for antenna_id, antenna in enumerate(antennas):
#if antenna_id != refant:
print(figcount)
if containsphase:
if notec:
phasep = phaseplot(phase[:,freqidx,antenna_id,0,0] - refphase)
else:
phasep = phaseplot(phase[:,antenna_id,0,0] + TEC2phase(tec[:,antenna_id,0,0], freq) - refphase)
else:
phasep = phaseplot(TEC2phase(tec[:,antenna_id,0,0], freq) - refphase)
if args['plotpoints']:
ax[figcount].plot(times, phasep, '.')
#ax[figcount].plot(times, phasep)
else:
ax[figcount].plot(times, phasep)
ax[figcount].set_ylabel('phase [rad]')
if figcount == len(antennas)-1:
ax[figcount].set_xlabel('time [hr]')
#print(type(antenna))
ax[figcount].set_title(antenna, position=(0.5, 0.75))
ax[figcount].set_ylim(-np.pi,np.pi)
figcount += 1
#plt.subplots_adjust(wspace=0, hspace=0)
plt.tight_layout(pad=1.0)
plt.savefig(str(args['outfile']))