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