-
Notifications
You must be signed in to change notification settings - Fork 1
/
Simulator.py
160 lines (132 loc) · 4.64 KB
/
Simulator.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
"""
Created Wed 11 April
@author Landman Bester
Set of utilities to simulate corrupted visibility data
"""
import numpy as np
import itertools as it
import utils
def sim_uv(Na, Nt, umax, vmax, Autocor=False, rot_params=(1,1)):
"""
Simulates elliptical uv-tracks coverage
:param Na: number of antennae
:param Nt: number of times
:param umax: max u coordinate
:param vmax: max v coordinate
:param Autocor: whether to include auto-correlations or not
:param rot_params: tuple or list specifying major and minor axes of rotation
:returns:
upq: u coordinates for baseline pq
vpq: v coordinates for baseline pq
pqlist: list of antennae pairs ordered in same order as upq and vpq axis 0
N: number of baselines
"""
if Autocor:
pqlist = list(it.combinations_with_replacement(np.arange(Na), 2))
N = Na * (Na + 1) // 2
else:
pqlist = list(it.combinations(np.arange(Na), 2))
N = Na * (Na - 1) // 2
# choose random antennae locations
u = umax*np.random.random(Na)
v = vmax*np.random.random(Na)
# create baselines with time axis
upq = np.zeros([N, Nt])
vpq = np.zeros([N, Nt])
phi = np.linspace(0, 1.5*np.pi, Nt) # to simulate earth rotation
for i, pq in enumerate(iter(pqlist)):
p = pq[0]
q = pq[1]
upq[i, 0] = u[p] - u[q]
vpq[i, 0] = v[p] - v[q]
for j in xrange(1, Nt):
rottheta = np.array([[rot_params[0]*np.cos(phi[j]), -rot_params[1]*np.sin(phi[j])],
[rot_params[1]*np.sin(phi[j]), rot_params[0]*np.cos(phi[j])]])
upq[i, j], vpq[i, j] = np.dot(rottheta, np.array([upq[i, 0], vpq[i, 0]]))
return upq, vpq, pqlist, N
def sim_sky(Npix, Nsource, max_I, lmax, mmax, freqs, ref_freq):
"""
Simulates a sky randomly populated with sources
:param Npix:
:param Nsource:
:param max_I:
:param lmax:
:param mmax:
:return:
"""
Nnu = freqs.size
l = np.linspace(-lmax, lmax, Npix)
m = np.linspace(-mmax, mmax, Npix)
ll, mm = np.meshgrid(l, m)
lmsource = []
locs = []
alpha = []
IM = np.zeros([Nnu, Nsource], dtype=np.float64)
for i in xrange(Nsource):
locx = np.random.randint(2, Npix-2)
locy = np.random.randint(2, Npix-2)
locs.append((locx, locy))
alpha.append(-0.7 + 0.1*np.random.randn(1))
I0 = np.abs(max_I*np.random.randn())
IM[:, i] = I0*(freqs/ref_freq)**alpha[i]
lmsource.append((ll[locx, locy], mm[locx, locy]))
return IM.squeeze(), lmsource, locs
def sim_T_gains(Na, N, theta, bounds=None):
"""
Simulates DDE's
:param Na: number of antennae
:param Ns: [Nnu, Nt]
:param thetas: [theta_nu, theta_t]
:param bounds: [(nu_min, nu_max), (t_min, t_max)]
:return:
"""
if bounds is not None:
t = np.linspace(bounds[0][0], bounds[0][1], N[0])
else:
t = np.linspace(0.0, 1.0, N[0])
x = np.array([t])
meanf = lambda x: np.ones(x[0].size, dtype=np.complex128)
gains, Kmat = utils.draw_samples_ND_grid(x, theta, Na, meanf=meanf)
return gains, Kmat
def sim_DI_gains(Na, Ns, thetas, bounds=None):
"""
Simulates DDE's
:param Na: number of antennae
:param Ns: [Nnu, Nt]
:param thetas: [theta_nu, theta_t]
:param bounds: [(nu_min, nu_max), (t_min, t_max)]
:return:
"""
if bounds is not None:
nu = np.linspace(bounds[0][0], bounds[0][1], Ns[0])
t = np.linspace(bounds[1][0], bounds[1][1], Ns[1])
else:
nu = np.linspace(1.0, 2.0, Ns[0])
t = np.linspace(0.0, 1.0, Ns[1])
x = np.array([nu, t])
meanf = lambda x: np.ones([x[0].size, x[1].size], dtype=np.complex128)
gains, Kmat = utils.draw_samples_ND_grid(x, thetas, Na, meanf=meanf)
return gains, Kmat
def sim_DD_gains(Na, Ns, thetas, lm, bounds=None):
"""
Simulates DDE's
:param Na: number of antennae
:param Ns: [Nnu, Nt]
:param thetas: [theta_nu, theta_t, theta_l, theta_m]
:param lm: 2 x Nsource array of source coordinates (only doing point sources for now)
:param bounds: [(nu_min, nu_max), (t_min, t_max)]
:return:
"""
if bounds is not None:
nu = np.linspace(bounds[0][0], bounds[0][1], Ns[0])
t = np.linspace(bounds[1][0], bounds[1][1], Ns[1])
else:
nu = np.linspace(1.0, 2.0, Ns[0])
t = np.linspace(0,1, Ns[1])
l = lm[0]
m = lm[1]
Nsource = l.size
x = np.array([nu, t, l, m])
meanf = lambda x: np.ones([x[0].size, x[1].size, x[2].size, x[3].size], dtype=np.complex128)
gains = utils.draw_samples_ND_grid(x, thetas, Na, meanf=meanf)
return gains