forked from dgary50/eovsa
-
Notifications
You must be signed in to change notification settings - Fork 0
/
delay_centers.py
123 lines (113 loc) · 5.12 KB
/
delay_centers.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
113
114
115
116
117
118
119
120
121
122
123
#
'''
Routine to capture packets while on satellite Galaxy 3C and analyze
them to check delay centers. Creates a plot of phase slopes, and
a table showing original delay centers (from acc's delay_centers.txt
file), and updated delay centers'''
#
# History:
# 2015-Jul-23 DG
# Started this history log. Fixed some errors in print-out of table,
# and improved formatting of the table.
# 2015-Aug-23 DG
# Changed to reflect new order of dimensions returned from
# pcapture.capture() routine.
#
import pcapture as p
import numpy as np
import matplotlib.pylab as plt
import urllib2
from util import Time
def delay_centers(overwrite=True):
''' Grabs raw packets while array is on the Galaxy 3C satellite,
and analyzes a specific range of channels in the 11-11.5 GHz
band. First it measures the phase slope and converts to delay
on each baseline, and then calculates the optimum delay
needed to correct for the phase slope.
TODO: It then optionally reads the delay_centers.txt file from
the ACC and updates it.
'''
# First see if we can read delay_centers.txt from the ACC (return if failure)
userpass = 'admin:observer@'
try:
# Read delay center file from ACC
dlafile = urllib2.urlopen('ftp://'+userpass+'acc.solar.pvt/parm/delay_centers.txt',timeout=1)
lines = dlafile.readlines()
delay_centers = []
dceny = []
for line in lines:
if line[0] != '#':
# Skip comment lines and take second number as delay center [ns]
delay_centers.append(float(line.strip().split()[1]))
dceny.append(float(line.strip().split()[2]))
except:
print Time().iso,'ACC connection for delay centers timed out.'
return [None]*5
# Capture 1 second of raw packets
out = p.capture('galaxy',1,'X',overwrite=overwrite)
# ch1,ch2 = 3600,3750
ch1,ch2 = 2050,3800
# Calculate the frequencies corresponding to selected channels ch1-ch2
f = np.linspace(ch1,ch2-1,ch2-ch1)*0.4/4096 + 12.15
# Determine which times in the capture contain valid (non-zero) data
good, = np.where(np.abs(out[0,0,100]) != 0)
outs = out[:,:,ch1:ch2,good]
good2, = np.where(np.abs(out[6,0,100]) != 0)
outs[6:] = out[6:,:,ch1:ch2,good2]
nbl,npol,nf,ngood = outs.shape
# NB: These settings hard-code the antenna assignments--probably needs
# to be aware of actual antenna assignments
bl = ['1-2','1-3','1-4','2-3','2-4','3-4','5-6','5-7','5-8','6-7','6-8','7-8']
pol = ['X','Y']
taus = []
c = np.zeros((nbl,npol))
for i in range(nbl):
# Determine the two antennas for this baseline (will break if ants > 9 -- i.e. not a single digit)
ant1 = int(bl[i][0])
ant2 = int(bl[i][2])
# For each baseline, determine phase slope by fitting a straight line
# to the mean phase over the relevant channels. Array pm are the
# two parameters of the fit [slope,offset] in radians
for j in range(npol):
phi = np.unwrap(np.angle(np.mean(outs[i,j],1))) # Mean over time axis
pm = np.polyfit(f,phi,1)
# plot the slopes and fits
plt.plot(f,phi,'.') # data points
y = np.polyval(pm,f)
plt.plot(f,y) # fit
chi_sq = np.sum((y - phi) ** 2)
# Convert slopes to delay in clock steps and in nsec
tau = int(pm[0]*0.8/(2*np.pi)+0.5) # delay in clock steps
tau_ns = pm[0]/(2*np.pi) # delay in nsec
plt.text(f[nf/2],y[nf/2],str(tau))
print bl[i]+pol[j],'=',tau,'(in ns:',tau_ns,'), chi_sq =',chi_sq
taus.append(tau_ns)
if j == 0:
c[i,j] = delay_centers[ant2-1] - delay_centers[ant1-1] - tau_ns
else:
c[i,j] = dceny[ant2-1] - dceny[ant1-1] - tau_ns
taus = np.array(taus)
# For now, analyze these in two sets of four
M = [[-1,1,0,0],[-1,0,1,0],[-1,0,0,1],[0,-1,1,0],[0,-1,0,1],[0,0,-1,1]]
x1,_,_,_ = np.linalg.lstsq(M,c[:6,0]) # Uses only X channel
x2,_,_,_ = np.linalg.lstsq(M,c[6:,0]) # Uses only X channel
y1,_,_,_ = np.linalg.lstsq(M,c[:6,1]) # Uses only Y channel
y2,_,_,_ = np.linalg.lstsq(M,c[6:,1]) # Uses only Y channel
# Calculate new delays for ants 1-4 based on Ant 1 as reference
newdlax1 = x1 - x1[0] + delay_centers[0]
newdlay1 = y1 - y1[0] + dceny[0]
# Calculate new delays for ants 5-8 based on Ant 5 as reference
newdlax2 = x2 - x2[0] + delay_centers[4]
newdlay2 = y2 - y2[0] + dceny[4]
for line in lines:
if line[0] == '#':
print line
else:
ant = int(line[:6])
if ant < 5:
print '{:4d} {:12.3f} {:12.3f} {:12.3f} {:12.3f}'.format(ant,delay_centers[ant-1],dceny[ant-1],newdlax1[ant-1],newdlay1[ant-1])
elif ant < 9:
print '{:4d} {:12.3f} {:12.3f} {:12.3f} {:12.3f}'.format(ant,delay_centers[ant-1],dceny[ant-1],newdlax2[ant-5],newdlay2[ant-5])
else:
print line[:-2]
return taus,c,M,x1,x2