1 | from __future__ import absolute_import
|
---|
2 |
|
---|
3 | import collections
|
---|
4 | import numpy as np
|
---|
5 | import pp
|
---|
6 |
|
---|
7 | import common.commonobjects as co
|
---|
8 |
|
---|
9 | def calculate_dirty_beam(npix, rpix, mult, bmax, bxbylist, wn_max, wn):
|
---|
10 | lamb = 1.0 / (wn * 100.0)
|
---|
11 | lambda_min = 1.0 / (wn_max * 100.0)
|
---|
12 | maxbx = bmax / lambda_min
|
---|
13 |
|
---|
14 | uv = numpy.zeros([mult*npix,mult*npix])
|
---|
15 |
|
---|
16 | for bxby in bxbylist:
|
---|
17 | ibx = int(round(((bxby[0] / lamb) / maxbx) * mult * rpix))
|
---|
18 | if ibx < 0:
|
---|
19 | ibx = mult*npix + ibx
|
---|
20 | iby = int(round(((bxby[1] / lamb) / maxbx) * mult * rpix))
|
---|
21 | if iby < 0:
|
---|
22 | iby = mult*npix + iby
|
---|
23 | uv[ibx,iby] = 1.0
|
---|
24 |
|
---|
25 | temp = numpy.fft.fftshift(uv)
|
---|
26 | temp = numpy.fft.fft2(temp)
|
---|
27 | temp = numpy.fft.fftshift(temp)
|
---|
28 | temp = numpy.power(temp, 2)
|
---|
29 | # convert to real, should be real anway as FT of symmetric function
|
---|
30 | dirty_beam = numpy.abs(temp)
|
---|
31 | # and normalise
|
---|
32 | dirty_beam /= numpy.max(dirty_beam)
|
---|
33 |
|
---|
34 | # truncate dirty beam pattern
|
---|
35 | dirty_beam = dirty_beam[
|
---|
36 | (mult-1)*rpix:(mult+1)*rpix, (mult-1)*rpix:(mult+1)*rpix]
|
---|
37 |
|
---|
38 | return dirty_beam
|
---|
39 |
|
---|
40 | def calculate_primary_beam(npix, rpix, mult, mx, my, bmax, m1_diameter,
|
---|
41 | wn_max, wn):
|
---|
42 | # assume uniform illumination of primary. uv goes out to
|
---|
43 | # +- bmax / lambda_min, fill uv plane appropriate to primary
|
---|
44 | # size
|
---|
45 | mirror = numpy.ones([mult*npix,mult*npix])
|
---|
46 | lambda_min = 1.0 / (wn_max * 100.0)
|
---|
47 | lamb = 1.0 / (wn * 100.0)
|
---|
48 | mirror[(numpy.sqrt(mx*mx + my*my) / (mult*rpix)) *
|
---|
49 | (bmax / lambda_min) > (0.5 * m1_diameter / lamb)] = 0.0
|
---|
50 |
|
---|
51 | # fftshift mirror to be centred at [0,0],
|
---|
52 | # do a 2d ftt of it,
|
---|
53 | # fftshift transform so that 0 freq is in centre of array
|
---|
54 | temp = numpy.fft.fftshift(mirror)
|
---|
55 | temp = numpy.fft.fft2(temp)
|
---|
56 | temp = numpy.fft.fftshift(temp)
|
---|
57 | # convert to real, should be real anyway as FT of symmetric function
|
---|
58 | # square to give intensity response
|
---|
59 | primary_amplitude_beam = numpy.real(temp)
|
---|
60 | primary_intensity_beam = numpy.power(numpy.abs(temp), 2)
|
---|
61 | # and normalise
|
---|
62 | primary_amplitude_beam /= numpy.max(primary_amplitude_beam)
|
---|
63 | primary_intensity_beam /= numpy.max(primary_intensity_beam)
|
---|
64 |
|
---|
65 | # truncate primary beam pattern
|
---|
66 | primary_amplitude_beam = primary_amplitude_beam[
|
---|
67 | (mult-1)*rpix:(mult+1)*rpix,
|
---|
68 | (mult-1)*rpix:(mult+1)*rpix]
|
---|
69 | primary_intensity_beam = primary_intensity_beam[
|
---|
70 | (mult-1)*rpix:(mult+1)*rpix,
|
---|
71 | (mult-1)*rpix:(mult+1)*rpix]
|
---|
72 |
|
---|
73 | return primary_intensity_beam, primary_amplitude_beam
|
---|
74 |
|
---|
75 |
|
---|
76 | class BeamsGenerator(object):
|
---|
77 | """Class to generate the primary beam(s) of the simulated observation.
|
---|
78 | """
|
---|
79 |
|
---|
80 | def __init__(self, previous_results):
|
---|
81 | self.previous_results = previous_results
|
---|
82 | self.result = collections.OrderedDict()
|
---|
83 |
|
---|
84 | # start parallel python (pp), find the number of CPUS available
|
---|
85 | ppservers = ()
|
---|
86 | self.job_server = pp.Server(ppservers=ppservers)
|
---|
87 | print 'Sampler starting pp with %s workers' % \
|
---|
88 | self.job_server.get_ncpus()
|
---|
89 |
|
---|
90 | def run(self):
|
---|
91 | print 'BeamsGenerator.run'
|
---|
92 |
|
---|
93 | uvmapgen = self.previous_results['uvmapgenerator']
|
---|
94 | fts = self.previous_results['fts']
|
---|
95 | telescope = self.previous_results['loadparameters']['substages']\
|
---|
96 | ['Telescope']
|
---|
97 |
|
---|
98 | # number of mirror positions sampled by FTS
|
---|
99 | nsample = fts['ftsnsample']
|
---|
100 | fts_wn = fts['fts_wn']
|
---|
101 | fts_lambda = fts['fts_lambda']
|
---|
102 |
|
---|
103 | # max pixel size (radians) that will fully sample the image.
|
---|
104 | # Sampling freq is 2 * Nyquist freq. So sampling freq = 2 * b_max / lambda_min.
|
---|
105 | bmax = uvmapgen['bmax']
|
---|
106 | self.result['max baseline [m]'] = bmax
|
---|
107 | max_pixsize = np.min(fts_lambda) / (2.0 * bmax)
|
---|
108 | self.result['max pixsize [rad]'] = max_pixsize
|
---|
109 |
|
---|
110 | # Number of pixels per beam
|
---|
111 | # radius of first null of Airy disk is theta=1.22lambda/d. Number of
|
---|
112 | # pixels is beam diameter/max_pixsize. Calculate this for the
|
---|
113 | # longest wavelength, largest beam
|
---|
114 | max_beam_radius = 1.22 * np.max(fts_lambda) /\
|
---|
115 | telescope['Primary mirror diameter']
|
---|
116 | self.result['beam diam [rad]'] = 2.0 * max_beam_radius
|
---|
117 |
|
---|
118 | # go out to radius of first null
|
---|
119 | rpix = int(max_beam_radius / max_pixsize)
|
---|
120 | npix = 2 * rpix
|
---|
121 | self.result['npix'] = npix
|
---|
122 |
|
---|
123 | # spatial axes same for all wavelengths
|
---|
124 | axis = np.arange(-rpix, rpix, dtype=np.float)
|
---|
125 | axis *= max_pixsize
|
---|
126 | axis *= 206264.806247
|
---|
127 | self.result['spatial axis [arcsec]'] = axis
|
---|
128 | axis1 = co.Axis(data=-axis, title='RA offset', units='arcsec')
|
---|
129 | axis2 = co.Axis(data=axis, title='Dec offset', units='arcsec')
|
---|
130 | axis3 = co.Axis(data=fts_wn, title='Frequency', units='cm-1')
|
---|
131 |
|
---|
132 | # calculate beams for each point on fts_wn
|
---|
133 | # oversample uv grid by mult to minimise aliasing
|
---|
134 | mult = 5
|
---|
135 | mx,my = np.meshgrid(np.arange(-mult*rpix, mult*rpix),
|
---|
136 | np.arange(-mult*rpix, mult*rpix))
|
---|
137 | self.result['primary beam'] = collections.OrderedDict()
|
---|
138 | self.result['primary amplitude beam'] = collections.OrderedDict()
|
---|
139 | self.result['dirty beam'] = collections.OrderedDict()
|
---|
140 |
|
---|
141 | wn_max = np.max(fts_wn)
|
---|
142 | m1_diameter = telescope['Primary mirror diameter']
|
---|
143 |
|
---|
144 | # first calculate primary beam for each wn
|
---|
145 | jobs = {}
|
---|
146 | for wn in fts_wn:
|
---|
147 | # submit jobs
|
---|
148 | indata = (npix, rpix, mult, mx, my, bmax, m1_diameter, wn_max,
|
---|
149 | wn,)
|
---|
150 | jobs[wn] = self.job_server.submit(calculate_primary_beam,
|
---|
151 | indata, (), ('numpy',))
|
---|
152 |
|
---|
153 | primary_amplitude_beam = np.zeros([npix,npix,len(fts_wn)], np.complex)
|
---|
154 | for iwn,wn in enumerate(fts_wn):
|
---|
155 | # collect and store results
|
---|
156 | primary_intensity_beam, primary_amplitude_beam[:,:,iwn] = jobs[wn]()
|
---|
157 | image = co.Image(data=primary_intensity_beam, axes=[axis1, axis2],
|
---|
158 | title='Primary Beam %06.4g cm-1' % wn)
|
---|
159 | self.result['primary beam'][wn] = image
|
---|
160 |
|
---|
161 | cube = co.Cube(data=primary_amplitude_beam, axes=[axis1, axis2, axis3],
|
---|
162 | title='Amplitude Primary Beam %06.4g cm-1' % wn)
|
---|
163 | self.result['primary amplitude beam'] = cube
|
---|
164 |
|
---|
165 | # second calculate dirty beam, same use of pp
|
---|
166 | for wn in fts_wn:
|
---|
167 | #
|
---|
168 | # 1. uv plane covered by discrete points so should use a
|
---|
169 | # slow dft rather than fft, which requires evenly spaced
|
---|
170 | # samples. Perhaps a bit slow though as calcs done in Python
|
---|
171 | # layer.
|
---|
172 | # 2. Could use shift theorem to give even-spaced equivalent
|
---|
173 | # of each discrete point, then use ffts. Equivalent to 1
|
---|
174 | # in speed.
|
---|
175 | # 3. Oversample uv plane and accept small error in point
|
---|
176 | # positions. Fastest but adequate?
|
---|
177 |
|
---|
178 | indata = (npix, rpix, mult, bmax, uvmapgen['bxby'], wn_max, wn,)
|
---|
179 | jobs[wn] = self.job_server.submit(calculate_dirty_beam,
|
---|
180 | indata, (), ('numpy',))
|
---|
181 |
|
---|
182 |
|
---|
183 | # axes for dirty beam
|
---|
184 | axis = np.arange(-rpix, rpix, dtype=np.float)
|
---|
185 | axis *= (np.min(fts_lambda) / bmax)
|
---|
186 | axis *= 206264.806247
|
---|
187 | axis1 = co.Axis(data=-axis, title='RA offset', units='arcsec')
|
---|
188 | axis2 = co.Axis(data=axis, title='Dec offset', units='arcsec')
|
---|
189 |
|
---|
190 | for wn in fts_wn:
|
---|
191 | dirty_beam = jobs[wn]()
|
---|
192 | image = co.Image(data=dirty_beam, axes=[axis1, axis2],
|
---|
193 | title='Dirty Beam %06.4g cm-1' % wn)
|
---|
194 | self.result['dirty beam'][wn] = image
|
---|
195 |
|
---|
196 | return self.result
|
---|
197 |
|
---|
198 | def __repr__(self):
|
---|
199 | return 'PrimaryBeamGenerator'
|
---|
200 |
|
---|