from __future__ import absolute_import import collections import numpy as np class Axis(object): """Class to hold axis information. """ def __init__(self, data, title='', units=''): self.data = data self.title = title self.units = units class Image(object): """Class to hold 2d image, with axes, units, title, etc. """ def __init__(self, data, axes=[], title='', units=''): self.data = data self.axes = axes self.title = title self.units = units class BeamsGenerator(object): """Class to generate the primary beam(s) of the simulated observation. """ def __init__(self, previous_results): self.previous_results = previous_results self.result = collections.OrderedDict() def run(self): print 'BeamsGenerator.run' uvmapgen = self.previous_results['uvmapgenerator'] fts = self.previous_results['fts'] telescope = self.previous_results['loadparameters']['substages']\ ['Telescope'] # number of mirror positions sampled by FTS nsample = fts['ftsnsample'] fts_wn = fts['fts_wn'] fts_lambda = fts['fts_lambda'] # max pixel size (radians) that will fully sample the image. # Sampling freq is 2 * Nyquist freq. So sampling freq = 2 * b_max / lambda_min. bmax = uvmapgen['bmax'] self.result['max baseline [m]'] = bmax max_pixsize = np.min(fts_lambda) / (2.0 * bmax) self.result['max pixsize [rad]'] = max_pixsize # Number of pixels per beam # radius of first null of Airy disk is theta=1.22lambda/d. Number of # pixels is beam diameter/max_pixsize. Calculate this for the # longest wavelength, largest beam max_beam_radius = 1.22 * np.max(fts_lambda) /\ telescope['Primary mirror diameter'] self.result['beam diam [rad]'] = 2.0 * max_beam_radius # go out to radius of first null rpix = int(max_beam_radius / max_pixsize) # rpix = 32 npix = 2 * rpix self.result['npix'] = npix # spatial axes same for all wavelengths axis = np.arange(-rpix, rpix, dtype=np.float) axis *= max_pixsize axis *= 206264.806247 self.result['spatial axis [arcsec]'] = axis # calculate beams for each point on fts_wn # oversample uv grid by mult to minimise aliasing mult = 5 mx,my = np.meshgrid(np.arange(-mult*rpix, mult*rpix), np.arange(-mult*rpix, mult*rpix)) self.result['primary beam'] = collections.OrderedDict() self.result['primary amplitude beam'] = collections.OrderedDict() self.result['dirty beam'] = collections.OrderedDict() for wn in fts_wn: # first calculate primary beam # # assume uniform illumination of primary. uv goes out to # +- bmax / lambda_min, fill uv plane appropriate to primary # size mirror = np.ones([mult*npix,mult*npix]) lambda_min = 1.0 / (np.max(fts_wn) * 100.0) lamb = 1.0 / (wn * 100.0) mirror[(np.sqrt(mx*mx + my*my) / (mult*rpix)) * (bmax / lambda_min) > \ (0.5 * telescope['Primary mirror diameter'] / lamb)] = 0.0 # fftshift mirror to be centred at [0,0], # do a 2d ftt of it, # fftshift transform so that 0 freq is in centre of array temp = np.fft.fftshift(mirror) temp = np.fft.fft2(temp) temp = np.fft.fftshift(temp) # convert to real, should be real anyway as FT of symmetric function # square to give intensity response primary_amplitude_beam = np.real(temp) primary_intensity_beam = np.power(np.abs(temp), 2) # and normalise primary_amplitude_beam /= np.max(primary_amplitude_beam) primary_intensity_beam /= np.max(primary_intensity_beam) # truncate primary beam pattern primary_amplitude_beam = primary_amplitude_beam[ (mult-1)*rpix:(mult+1)*rpix, (mult-1)*rpix:(mult+1)*rpix] primary_intensity_beam = primary_intensity_beam[ (mult-1)*rpix:(mult+1)*rpix, (mult-1)*rpix:(mult+1)*rpix] # axes axis1 = Axis(data=-axis, title='RA offset', units='arcsec') axis2 = Axis(data=axis, title='Dec offset', units='arcsec') image = Image(data=primary_intensity_beam, axes=[axis1, axis2], title='Primary Beam %06.4g cm-1' % wn) self.result['primary beam'][wn] = image image = Image(data=primary_amplitude_beam, axes=[axis1, axis2], title='Amplitude Primary Beam %06.4g cm-1' % wn) self.result['primary amplitude beam'][wn] = image # second calculate dirty beam # # 1. uv plane covered by discrete points so should use a # slow dft rather than fft, which requires evenly spaced # samples. Perhaps a bit slow though as calcs done in Python # layer. # 2. Could use shift theorem to give even-spaced equivalent # of each discrete point, then use ffts. Equivalent to 1 # in speed. # 3. Oversample uv plane and accept small error in point # positions. Fastest but adequate? uv = np.zeros([mult*npix,mult*npix]) maxbx = bmax / np.min(fts_lambda) for bxby in uvmapgen['bxby']: ibx = int(round(((bxby[0] / lamb) / maxbx) * mult * rpix)) if ibx < 0: ibx = mult*npix + ibx iby = int(round(((bxby[1] / lamb) / maxbx) * mult * rpix)) if iby < 0: iby = mult*npix + iby uv[ibx,iby] = 1.0 temp = np.fft.fftshift(uv) temp = np.fft.fft2(temp) temp = np.fft.fftshift(temp) temp = np.power(temp, 2) # convert to real, should be real anway as FT of symmetric function dirty_beam = np.abs(temp) # and normalise dirty_beam /= np.max(dirty_beam) # truncate dirty beam pattern dirty_beam = dirty_beam[ (mult-1)*rpix:(mult+1)*rpix, (mult-1)*rpix:(mult+1)*rpix] # axes axis = np.arange(-rpix, rpix, dtype=np.float) axis *= (np.min(fts_lambda) / bmax) axis *= 206264.806247 axis1 = Axis(data=-axis, title='RA offset', units='arcsec') axis2 = Axis(data=axis, title='Dec offset', units='arcsec') image = Image(data=dirty_beam, axes=[axis1, axis2], title='Dirty Beam %06.4g cm-1' % wn) self.result['dirty beam'][wn] = image return self.result def __repr__(self): return 'PrimaryBeamGenerator'