Index: amsgenerator.py
===================================================================
--- /beamsgenerator.py	(revision 6)
+++ 	(revision )
@@ -1,183 +1,0 @@
-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'
-
