-
Notifications
You must be signed in to change notification settings - Fork 0
/
bisp_sim.py
138 lines (121 loc) · 12.6 KB
/
bisp_sim.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
#!/usr/bin/env python
#
# claw, 11jan13
#
# Simulating bispectrum for obs planning
# How to find f.o.v, number of beams, sensitivity
import numpy as n
import pylab as p
# antenna positions
#vla_cnb = n.array([[-1.60115e+06, -5.04173e+06, 3.55524e+06], [-1.60095e+06, -5.04213e+06, 3.55477e+06], [-1.59964e+06, -5.04295e+06, 3.5542e+06], [-1.60161e+06, -5.042e+06, 3.55465e+06], [-1.60078e+06, -5.03935e+06, 3.55876e+06], [-1.60107e+06, -5.04205e+06, 3.55482e+06], [-1.60111e+06, -5.04149e+06, 3.5556e+06], [-1.60115e+06, -5.042e+06, 3.55486e+06], [-1.601e+06, -5.0408e+06, 3.55661e+06], [-1.60259e+06, -5.04205e+06, 3.55414e+06], [-1.6008e+06, -5.04222e+06, 3.55471e+06], [-1.60171e+06, -5.04201e+06, 3.5546e+06], [-1.60204e+06, -5.04203e+06, 3.55443e+06], [-1.60291e+06, -5.04207e+06, 3.55398e+06], [-1.60117e+06, -5.0419e+06, 3.55499e+06], [-1.60042e+06, -5.04246e+06, 3.55454e+06], [-1.60018e+06, -5.04261e+06, 3.55443e+06], [-1.60123e+06, -5.04198e+06, 3.55486e+06], [-1.60069e+06, -5.03876e+06, 3.55963e+06], [-1.60062e+06, -5.04233e+06, 3.55463e+06], [-1.60093e+06, -5.04032e+06, 3.55733e+06], [-1.60086e+06, -5.03989e+06, 3.55797e+06], [-1.6023e+06, -5.04204e+06, 3.55429e+06], [-1.60132e+06, -5.04199e+06, 3.55481e+06], [-1.60145e+06, -5.04199e+06, 3.55474e+06], [-1.59993e+06, -5.04277e+06, 3.55432e+06], [-1.60106e+06, -5.04118e+06, 3.55606e+06]]) # 27 ants from sgra* obs
#vla_b = n.array([[ 4.74212121e+01, -1.23351515e+02, -6.83363636e+01], [ 1.54403030e+02, -4.05618182e+02, -2.25827273e+02], [ 3.09478788e+02, -8.13260606e+02, -4.52918182e+02], [ 5.05233333e+02, -1.33222727e+03, -7.43157576e+02], [ 7.41251515e+02, -1.95385455e+03, -1.08975455e+03], [ 1.01627576e+03, -2.67153939e+03, -1.48810303e+03], [ 1.33065455e+03, -3.48049697e+03, -1.93422121e+03], [ 1.65772727e+03, -4.37670606e+03, -2.44280303e+03], [ 2.02165758e+03, -5.35702121e+03, -2.99490606e+03], [ 3.59878788e+01, 1.35081818e+02, -5.16545455e+01], [ 1.15660606e+02, 4.43433333e+02, -1.71318182e+02], [ 2.31936364e+02, 8.88790909e+02, -3.43521212e+02], [ 3.81045455e+02, 1.45656061e+03, -5.62133333e+02], [ 5.66142424e+02, 2.13606667e+03, -8.19439394e+02], [ 7.73469697e+02, 2.92066667e+03, -1.12087273e+03], [ 1.00944545e+03, 3.80498485e+03, -1.45906061e+03], [ 1.26676970e+03, 4.78475152e+03, -1.83654545e+03], [ 1.55113636e+03, 5.85647576e+03, -2.24751515e+03], [ -7.56333333e+01, -1.18636364e+01, 1.12821212e+02],[ -2.46236364e+02, -3.84484848e+01, 3.93333333e+00], [ -4.94572727e+02, -7.71121212e+01, 7.29309091e+02], [ -8.10057576e+02, -1.26327273e+02, 1.19487879e+03], [ -1.18700909e+03, -1.85269697e+02, 1.75297576e+03], [ -1.67846364e+03, -2.62169697e+02, 2.48091515e+03], [ -2.11407273e+03, -3.30090909e+02, 3.12277576e+03], [ -2.65750303e+03, -4.15093939e+02, 3.92758182e+03], [ -3.25232121e+03, -5.08378788e+02, 4.80751212e+03]]) # from miriad/cat
wsrt = n.array([[ 0., 0. ], [ 120., 0. ], [ 240., 0. ], [ 360., 0. ], [ 480., 0. ], [ 600., 0. ], [ 720., 0. ], [ 840., 0. ], [ 960., 0. ], [ 1080., 0. ], [ 1200., 0. ], [ 1320., 0. ]] )
vla_d = 0.57*1e3*n.array([[ 0.00305045, 0.03486681], [ 0.00893224, 0.10209601], [ 0.01674565, 0.19140365], [ 0.02615514, 0.29895461], [ 0.03696303, 0.42248936], [ 0.04903413, 0.56046269], [ 0.06226816, 0.7117283 ], [ 0.07658673, 0.87539034], [ 0.09192633, 1.05072281], [ 0.02867032, -0.02007518], [ 0.08395162, -0.05878355], [ 0.1573876 , -0.11020398], [ 0.24582472, -0.17212832], [ 0.347405 , -0.2432556 ], [ 0.46085786, -0.32269615], [ 0.58524071, -0.40978995], [ 0.71981691, -0.50402122], [ 0.86398948, -0.60497195], [-0.03172077, -0.01479164], [-0.09288386, -0.04331245], [-0.17413325, -0.08119967], [-0.27197986, -0.12682629], [-0.38436803, -0.17923376], [-0.509892 , -0.23776654], [-0.64750886, -0.30193834], [-0.79640364, -0.37136912], [-0.95591581, -0.44575086]])
vla_c = 0.6*1e3*n.array([[ 0.00958713, 0.10958142], [ 0.02807276, 0.32087317], [ 0.05262918, 0.60155433], [ 0.08220187, 0.93957164], [ 0.11616952, 1.32782371], [ 0.15410727, 1.76145417], [ 0.19569992, 2.23686036], [ 0.24070115, 2.75122679], [ 0.28891134, 3.30227169], [ 0.09010672, -0.06309341], [ 0.26384793, -0.18474831], [ 0.49464674, -0.34635537], [ 0.77259197, -0.54097472], [ 1.0918443 , -0.76451761], [ 1.44841043, -1.0141879 ], [ 1.83932793, -1.28791128], [ 2.26228171, -1.58406671], [ 2.7153955 , -1.9013404 ], [-0.09969386, -0.04648801], [-0.2919207 , -0.13612486], [-0.54727592, -0.25519895], [-0.85479384, -0.39859691], [-1.20801382, -0.5633061 ], [-1.6025177 , -0.74726628], [-2.03502786, -0.94894907], [-2.50298286, -1.16716008], [-3.00430684, -1.40093129]])
vla_b = 0.58*1e3*n.array([[ 0.02309627, 0.26399159], [ 0.07504001, 0.85771122], [ 0.14950235, 1.70881966], [ 0.24380571, 2.78671198], [ 0.35627955, 4.07229384], [ 0.48573456, 5.55197147], [ 0.63126041, 7.21533946], [ 0.79212708, 9.05405395], [ 0.96772987, 11.06120308], [ 0.21707529, -0.15199776], [ 0.7052797 , -0.49384217], [ 1.40513007, -0.98388266], [ 2.29146052, -1.60449793], [ 3.34857014, -2.34469406], [ 4.56528105, -3.19664421], [ 5.93303706, -4.15435728], [ 7.44497719, -5.21302915], [ 9.09541793, -6.3686802 ], [ -0.24017156, -0.11199384], [ -0.78031971, -0.36386906], [ -1.55463241, -0.724937 ], [ -2.53526622, -1.18221406], [ -3.70484969, -1.72759978], [ -5.05101561, -2.35532726], [ -6.56429747, -3.06098218], [ -8.23710427, -3.8410248 ], [-10.0631478 , -4.69252289]])
vla_bP = 0.58*1e3*n.array([[ 0.02309627, 0.26399159], [ 0.07504001, 0.85771122], [ 0.14950235, 1.70881966], [ 0.24380571, 2.78671198], [ 0.35627955, 4.07229384], [ 0.48573456, 5.55197147], [ 0.63126041, 7.21533946], [ 0.79212708, 9.05405395], [ 0.96772987, 11.06120308], [ 0.21707529, -0.15199776], [ 3.34857014, -2.34469406], [ 4.56528105, -3.19664421], [ 5.93303706, -4.15435728], [ 7.44497719, -5.21302915], [ 9.09541793, -6.3686802 ], [ -0.24017156, -0.11199384], [ -1.55463241, -0.724937 ], [ -2.53526622, -1.18221406], [ -6.56429747, -3.06098218], [-10.0631478 , -4.69252289]])
vla_a = 0.57*1e3*n.array([[ 0.06100902, 0.69733629], [ 0.21244575, 2.428266 ], [ 0.44077028, 5.03802739], [ 0.73977906, 8.45571332], [ 1.10545242, 12.63537902], [ 1.53485127, 17.54343034], [ 2.02567851, 23.1536113 ], [ 2.5760603 , 29.44450398], [ 3.18442163, 36.3981058 ], [ 0.57340643, -0.40150351], [ 1.99671717, -1.39811641], [ 4.14267457, -2.90073196], [ 6.95297301, -4.86852412], [ 10.38983301, -7.27503939], [ 14.42563071, -10.10093536], [ 19.03877632, -13.3310947 ], [ 24.2116583 , -16.95318565], [ 29.92947346, -20.95684293], [ -0.63441545, -0.29583278], [ -2.20916291, -1.03014958], [ -4.58344485, -2.13729543], [ -7.69275207, -3.5871892 ], [-11.49528543, -5.36033963], [-15.96048198, -7.44249498], [-21.06445483, -9.8225166 ], [-26.7877186 , -12.49131833], [-33.11389509, -15.44126287]])
vla_cnb = 0.58*1e3*n.array([[ 0.02309627, 0.26399159], [ 0.07504001, 0.85771122], [ 0.14950235, 1.70881966], [ 0.24380571, 2.78671198], [ 0.35627955, 4.07229384], [ 0.48573456, 5.55197147], [ 0.63126041, 7.21533946], [ 0.79212708, 9.05405395], [ 0.96772987, 11.06120308], [ 0.09010672, -0.06309341], [ 0.26384793, -0.18474831], [ 0.49464674, -0.34635537], [ 0.77259197, -0.54097472], [ 1.0918443 , -0.76451761], [ 1.44841043, -1.0141879 ], [ 1.83932793, -1.28791128], [ 2.26228171, -1.58406671], [ 2.7153955 , -1.9013404 ], [-0.09969386, -0.04648801], [-0.2919207 , -0.13612486], [-0.54727592, -0.25519895], [-0.85479384, -0.39859691], [-1.20801382, -0.5633061 ], [-1.6025177 , -0.74726628], [-2.03502786, -0.94894907], [-2.50298286, -1.16716008], [-3.00430684, -1.40093129]])
vla_nb = 0.58*1e3*n.array([[ 0.02309627, 0.26399159], [ 0.07504001, 0.85771122], [ 0.14950235, 1.70881966], [ 0.24380571, 2.78671198], [ 0.35627955, 4.07229384], [ 0.48573456, 5.55197147], [ 0.63126041, 7.21533946], [ 0.79212708, 9.05405395], [ 0.96772987, 11.06120308]])
class receiver():
def __init__(self,freq,bw):
"""Receiver (band) specs.
"""
self.freq = freq
self.bw = bw
self.ll = 3e8/freq
class antenna():
def __init__(self,receiver,x,y,diameter):
"""Antenna contains a band, size and location.
"""
self.x = x
self.y = y
self.diameter = diameter
self.fov = n.degrees(receiver.ll/diameter)
self.ll = receiver.ll
self.freq = receiver.freq
self.bw = receiver.bw
class baseline():
def __init__(self, antenna1, antenna2):
"""Baseline contains two antennas
"""
# at zenith
self.u = (antenna1.x-antenna2.x)/antenna1.ll
self.v = (antenna1.y-antenna2.y)/antenna1.ll
self.delaybeamhwhm = (2 / n.sqrt(self.u**2+self.v**2)) * (antenna1.freq/antenna1.bw)
# self.delaybeamangle = n.angle(n.complex(self.u,self.v)) ??
# self.beam = gaussian2d(self.delaybeamhwhm, self.fov, self.delaybeamangle) ??
class triple():
def __init__(self, baseline1, baseline2, baseline3):
"""Triple made of three (could be two) baselines
"""
self.beam = baseline1.beam * baseline2.beam * baseline3.beam
def build_array(xyarr, freq=1.4e9, bw=128.e6, diameter=25., SEFD=400., tint=0.01, npol=2, thresh=5, uvlim=1e100):
"""create all baselines in array. produces:
-- max freq to avoid delay beam.
-- number of beams required to use all of bw
assumes circular delay beam for now.
xyarr is xy coords for each ant.
"""
rec = receiver(freq, bw)
na = len(xyarr)
blarr = []; blu = []; blv = []; uvdist = []
ntr = 0
for i in range(na):
ant1 = antenna(rec, xyarr[i,0], xyarr[i,1], diameter)
for j in range(i+1, na):
ant2 = antenna(rec, xyarr[j,0], xyarr[j,1], diameter)
bl = baseline(ant1, ant2)
if n.sqrt(bl.u**2 + bl.v**2) <= uvlim:
blarr.append(bl)
blu.append(bl.u)
blv.append(bl.v)
uvdist = n.sqrt(bl.u**2+bl.v**2)
for k in range(j+1,na):
ant3 = antenna(rec, xyarr[k,0], xyarr[k,1], diameter)
bl2 = baseline(ant2, ant3)
uvdist2 = n.sqrt( bl2.u**2 + bl2.v**2 )
uvdist3 = n.sqrt( (bl.u + bl2.u)**2 + (bl.v + bl2.v)**2 )
if ((uvdist <= uvlim) & (uvdist2 <= uvlim) & (uvdist3 <= uvlim)):
ntr+=1
blu = n.array(blu); blv = n.array(blv)
uvdist = n.sqrt(blu**2+blv**2)
uvmax = uvdist.max()
nbl = len(uvdist)
FWHM_delay = 4 * (n.degrees(1/uvmax)) * freq/bw
FWHM_primary = n.degrees(3e8/freq/diameter)
N_delaybeams = (FWHM_primary/FWHM_delay)**2
delay_bw = 4 * freq * diameter/(uvmax * 3e8/freq)
N_bw = bw/delay_bw
Q = SEFD / n.sqrt(2*tint*bw) # snr per bl per pol
sens_im = Q * thresh / n.sqrt(npol*nbl)
sens_bisp_allbw = Q * n.power(2 * thresh / n.sqrt(npol*ntr), 1/3.) # assuming all bw averaged. post-det comb of pols (like triples).
sens_bisp_allfov = Q * n.sqrt(bw/delay_bw) * n.power(2 * thresh / n.sqrt(N_bw*npol*ntr), 1/3.) # assuming delay beam must be primary beam. thus, does post-det comb of bw in segments of delay_bw (like triples)
print
print 'Observing configuration has %d antennas, %d baselines, and %d triples' % (na, nbl, ntr)
print 'Longest baseline (m): %d, Synth beam (asec, zenith): %.1f, Prim beam (deg) %.2f' % (uvmax*(3e8/freq), 3600*n.degrees(1/uvmax), FWHM_primary)
print 'Frequency (GHz): %.2f, Bandwidth (MHz): %d, SEFD (Jy): %d, tint (s): %.3f' % (freq/1e9, bw/1e6, SEFD, tint)
print
print 'Noise per bl-pol, Q (Jy, full BW): %.3f. Imaging sens (Jy, %d sigma): %.3f' % (Q, thresh, sens_im)
print
print 'Bispectrum parameters for whole BW (zenith):'
print 'FWHM_delay (deg), N_delaybeams, delay_bw (MHz), Bisp sens (Jy, %d sigma, ideal):' % (thresh)
print '%.3f, %.1f, %.3f' % (FWHM_delay, N_delaybeams, sens_bisp_allbw)
print
print 'Bispectrum parameters for whole FOV (zenith):'
if N_bw >= 1:
print 'delay_bw (MHz), N_bw, bisp sens (Jy, %d sigma, incoh)' % (thresh)
print '%.3f, %.1f, %.3f' % (delay_bw/1e6, N_bw, sens_bisp_allfov)
else:
print 'delay_bw (MHz), N_bw'
print '%.3f, %.1f' % (delay_bw/1e6, N_bw)
p.figure(1)
p.subplot(121)
p.plot(xyarr[:,0], xyarr[:,1],'.')
p.xlabel('Antenna x (m)')
p.ylabel('Antenna y (m)')
p.subplot(122)
p.plot(blu, blv, '.')
p.xlabel('Baseline u (lambda)')
p.ylabel('Baseline v (lambda)')
p.show()