Index: trunk/base.py
===================================================================
--- trunk/base.py	(revision 71)
+++ 	(revision )
@@ -1,23 +1,0 @@
-from __future__ import absolute_import
-
-from collections import defaultdict
-import os
-import pprint
-import xlrd
-
-from mako.template import Template
-
-class Base:
-    """
-    """
-
-    @staticmethod
-    def render(context):
-        print 'base.render'
-        link = os.path.join(context['dirname'], 't2-1-%s.html' % context['result'][0])
-        with open(link, 'w') as f:
-            t2_1 = Template(filename='t2-1empty.template')
-            f.write(t2_1.render(**context))
-
-        return os.path.basename(link)
-            
Index: trunk/beamsgenerator.py
===================================================================
--- trunk/beamsgenerator.py	(revision 71)
+++ 	(revision )
@@ -1,222 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import numpy as np
-import pp
-
-import common.commonobjects as co
-
-def calculate_dirty_beam(npix, rpix, mult, bmax, bxbylist, wnmax, wn):
-    lamb = 1.0 / (wn * 100.0)
-    lambda_min = 1.0 / (wnmax * 100.0)
-    maxbx = bmax / lambda_min
-
-    uv = numpy.zeros([mult*npix,mult*npix])
-
-    for bxby in bxbylist:
-        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 = numpy.fft.fftshift(uv)
-    temp = numpy.fft.fft2(temp)
-    temp = numpy.fft.fftshift(temp)
-    temp = numpy.power(temp, 2)
-    # convert to real, should be real anway as FT of symmetric function 
-    dirty_beam = numpy.abs(temp)
-    # and normalise
-    dirty_beam /= numpy.max(dirty_beam)
-
-    # truncate dirty beam pattern
-    dirty_beam = dirty_beam[
-      (mult-1)*rpix:(mult+1)*rpix, (mult-1)*rpix:(mult+1)*rpix]
-
-    # fit Gaussian to dirty beam
-    shape = numpy.shape(dirty_beam)
-    fitter = fitgaussian.FitGaussian()
-    dirty_centre = numpy.array(dirty_beam[shape[0]/2-5:shape[0]/2+5,
-      shape[1]/2-5:shape[1]/2+5])
-    p = fitter.fitgaussian(dirty_centre)
-
-    # clean beam
-    cp = (1.0, float(shape[0])/2.0, float(shape[1])/2.0, p[3], p[4], p[5])
-    rotgauss = fitter.gaussian(*cp)
-    bmax = numpy.max(dirty_beam)
-    clean_beam = numpy.fromfunction(rotgauss, numpy.shape(dirty_beam))
-
-    return dirty_beam, clean_beam
-
-def calculate_primary_beam(npix, rpix, mult, mx, my, bmax, m1_diameter,
-  wnmax, wn):
-    # assume uniform illumination of primary. uv goes out to
-    # +- bmax / lambda_min, fill uv plane appropriate to primary
-    # size
-    mirror = numpy.ones([mult*npix,mult*npix])
-    lambda_min = 1.0 / (wnmax * 100.0)
-    lamb = 1.0 / (wn * 100.0)
-    mirror[(numpy.sqrt(mx*mx + my*my) / (mult*rpix)) * 
-      (bmax / lambda_min) > (0.5 * m1_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 = numpy.fft.fftshift(mirror)
-    temp = numpy.fft.fft2(temp)
-    temp = numpy.fft.fftshift(temp)
-    # convert to real, should be real anyway as FT of symmetric function 
-    # square to give intensity response
-    primary_amplitude_beam = numpy.real(temp)
-    primary_intensity_beam = numpy.power(numpy.abs(temp), 2)
-    # and normalise
-    primary_amplitude_beam /= numpy.max(primary_amplitude_beam)
-    primary_intensity_beam /= numpy.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]
-
-    return primary_intensity_beam, primary_amplitude_beam
-
-
-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()
-
-        # start parallel python (pp), find the number of CPUS available
-        ppservers = ()
-        self.job_server = pp.Server(ppservers=ppservers)
-        print 'Sampler starting pp with %s workers' % \
-          self.job_server.get_ncpus()
-
-    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_wn_truncated = fts['fts_wn_truncated']
-      
-        # 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
-
-        lambda_min = 1.0 / (np.max(fts_wn) * 100.0)
-        max_pixsize = lambda_min / (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 that has flux (at wnmin), largest beam
-        
-        max_beam_radius = 1.22 * (1.0 / (np.min(fts_wn_truncated) * 100.0)) /\
-          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)
-        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
-        axis1 = co.Axis(data=-axis, title='RA offset', units='arcsec') 
-        axis2 = co.Axis(data=axis, title='Dec offset', units='arcsec') 
-        axis3 = co.Axis(data=fts_wn_truncated, title='Frequency', units='cm-1')
-
-        # calculate beams for each point on fts_wn_truncated
-        # 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()
-        self.result['clean beam'] = collections.OrderedDict()
-
-        wnmax = np.max(fts_wn_truncated)
-        m1_diameter = telescope['Primary mirror diameter']
-
-        # first calculate primary beam for each wn
-        jobs = {}
-        for wn in fts_wn_truncated:
-            # submit jobs
-            indata = (npix, rpix, mult, mx, my, bmax, m1_diameter, wnmax,
-              wn,)
-            jobs[wn] = self.job_server.submit(calculate_primary_beam,
-              indata, (), ('numpy',))
-
-        primary_amplitude_beam = np.zeros([npix,npix,len(fts_wn_truncated)],
-          np.complex)
-        for iwn,wn in enumerate(fts_wn_truncated):
-            # collect and store results
-            primary_intensity_beam, primary_amplitude_beam[:,:,iwn] = jobs[wn]()
-            image = co.Image(data=primary_intensity_beam, axes=[axis1, axis2],
-              title='Primary Beam %06.4g cm-1' % wn)
-            self.result['primary beam'][wn] = image
-
-        cube = co.Cube(data=primary_amplitude_beam, axes=[axis1, axis2, axis3],
-          title='Amplitude Primary Beam %06.4g cm-1' % wn)
-        self.result['primary amplitude beam'] = cube
-
-        # second calculate dirty beam, same use of pp
-        for wn in fts_wn_truncated:
-            #
-            # 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?
-
-            indata = (npix, rpix, mult, bmax, uvmapgen['bxby'], wnmax, wn,)
-            jobs[wn] = self.job_server.submit(calculate_dirty_beam,
-              indata, (), ('numpy','fitgaussian',))
-
-
-        # axes for dirty beam
-        axis = np.arange(-rpix, rpix, dtype=np.float)
-        axis *= (lambda_min / bmax)
-        axis *= 206264.806247
-        axis1 = co.Axis(data=-axis, title='RA offset', units='arcsec') 
-        axis2 = co.Axis(data=axis, title='Dec offset', units='arcsec') 
-
-        for wn in fts_wn_truncated:
-            dirty_beam,clean_beam = jobs[wn]()
-            image = co.Image(data=dirty_beam, axes=[axis1, axis2],
-              title='Dirty Beam %06.4g cm-1' % wn)
-            self.result['dirty beam'][wn] = image
-            image = co.Image(data=clean_beam, axes=[axis1, axis2],
-              title='Clean Beam %06.4g cm-1' % wn)
-            self.result['clean beam'][wn] = image
-
-        return self.result
-
-    def __repr__(self):
-        return 'PrimaryBeamGenerator'
-
Index: trunk/cleanimage.py
===================================================================
--- trunk/cleanimage.py	(revision 71)
+++ 	(revision )
@@ -1,67 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import math
-import numpy as np
-import pp
-
-import common.commonobjects as co
-import pythonclean
-
-
-class CleanImage(object):
-    """Class to compute clean image cube from uvspectra.
-    """
-
-    def __init__(self, previous_results, job_server):
-        self.previous_results = previous_results
-        self.job_server = job_server
-
-        self.result = collections.OrderedDict()      
-
-    def run(self):
-        print 'CleanImage.run'
-
-        # get dirty image and beam info
-        dirty = self.previous_results['dirtyimage']
-        dirtyimage = dirty['dirtyimage']
-        dirtybeam = dirty['dirtybeam']
-        spatial_axis = dirty['spatial axis [arcsec]']
-        wavenumber = dirty['wavenumber [cm-1]']
-
-        # clean image cube
-        cleanimage = np.zeros(np.shape(dirtyimage), np.float)
-        residualimage = np.zeros(np.shape(dirtyimage), np.float)
-
-        # calculate clean image for each wn
-        jobs = {}
-        # clean only central quarter
-        imsize = np.shape(dirtyimage)[:2]
-        window = [slice(imsize[0]/4, imsize[0]*3/4), slice(imsize[1]/4, imsize[1]*3/4)]
-        for iwn,wn in enumerate(wavenumber):
-            # submit jobs
-            immax = np.max(dirtyimage[:,:,iwn])
-            indata = (dirtyimage[:,:,iwn], dirtybeam[:,:,iwn], window, 0.1,
-              immax/100.0, 5000,)
-            jobs[wn] = self.job_server.submit(pythonclean.hogbom,
-              indata, (), ('numpy','pythonclean',))
-
-        for iwn,wn in enumerate(wavenumber):
-            # collect and store results
-            cleanimage[:,:,iwn], residualimage[:,:,iwn] = jobs[wn]()
-
-        self.result['dirtyimage'] = dirtyimage
-        self.result['cleanimage'] = cleanimage + residualimage
-        self.result['residualimage'] = residualimage
-        self.result['spatial axis [arcsec]'] = spatial_axis
-        self.result['spatial axis'] = spatial_axis
-        self.result['wavenumber [cm-1]'] = wavenumber
-
-        return self.result
-                
-    def __repr__(self):
-        return '''
-CleanImage:
-'''.format(
-          num_uvspectra=len(self.uvspectra))
-
Index: trunk/dirtyimage.py
===================================================================
--- trunk/dirtyimage.py	(revision 71)
+++ 	(revision )
@@ -1,117 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import math
-import numpy as np
-import pp
-
-import common.commonobjects as co
-
-def dirty(b_x_list, b_y_list, spectra, wn_spectra, iwn, wn, spatial_axis, npix):
-    image = numpy.zeros([npix, npix], numpy.float)
-    beam = numpy.zeros([npix, npix], numpy.float)
-
-    for ibx,b_x in enumerate(b_x_list):
-        b_x = b_x_list[ibx]
-        b_y = b_y_list[ibx]
-        spectrum = spectra[ibx]
-        wn_spectrum = wn_spectra[ibx]
-
-        # must be a better way of doing this
-        argx = numpy.radians(b_x * (wn * 100.0) * spatial_axis / 3600.0)
-        argy = numpy.radians(b_y * (wn * 100.0) * spatial_axis / 3600.0)
-        valx = numpy.exp(2.0j * math.pi * argx)
-        valy = numpy.exp(2.0j * math.pi * argy)
-
-        vis = spectrum[wn_spectrum==wn]
-
-        # numpy arrays are [row,col] or here [y,x]
-        for iy,y in enumerate(spatial_axis):
-            contribution = (vis * valx * valy[iy]) + \
-              numpy.conjugate(vis * valx * valy[iy])
-            image[iy,:] += contribution.real
-            contribution = (valx * valy[iy]) + \
-              numpy.conjugate(valx * valy[iy])
-            beam[iy,:] += contribution.real
-   
-    # normalise
-    # The sum of the dirty beam should equal the 0 baseline vis measurement: 0
-    # Normalising the volume would preserve flux under convolution but
-    # is problematic as the volume is 0, so normalise the peak to 1.
-    beam /= 2.0 * len(b_x_list)
-
-    # likewise for dirty map
-    image /= 2.0 * len(b_x_list)
-
-    return image, beam
-
-
-class DirtyImage(object):
-    """Class to compute dirty image cube from uvspectra.
-    """
-
-    def __init__(self, previous_results, job_server):
-        self.previous_results = previous_results
-        self.job_server = job_server
-
-        self.result = collections.OrderedDict()      
-
-    def run(self):
-        print 'DirtyImage.run'
-
-        # get relevant fts info
-        fts = self.previous_results['fts']
-        fts_wnmin = fts['wnmin']
-
-        # get primary beam info
-        beamsgenerator = self.previous_results['beamsgenerator']
-        npix = beamsgenerator['npix']
-        spatial_axis = beamsgenerator['spatial axis [arcsec]']
-
-        # get observation list
-        uvspectra = self.previous_results['uvspectra']['uvspectra']
-        wavenumber = uvspectra[0].wavenumber
-        wavenumber = wavenumber[wavenumber>fts_wnmin]
-
-        # dirty image cube
-        dirtyimage = np.zeros([npix, npix, len(wavenumber)], np.float)
-        dirtybeam = np.zeros([npix, npix, len(wavenumber)], np.float)
-
-        # calculate dirty image for each wn
-        # uvspectrum objects don't pickle so unpack them first
-        b_x_list = []
-        b_y_list = []
-        spectra = []
-        wn_spectra = []
-        for uvspectrum in uvspectra:
-            b_x_list.append(uvspectrum.baseline_x)
-            b_y_list.append(uvspectrum.baseline_y)
-            spectra.append(uvspectrum.spectrum)
-            wn_spectra.append(uvspectrum.wavenumber)
-
-        jobs = {}
-        for iwn,wn in enumerate(wavenumber):
-            # submit jobs
-            indata = (b_x_list, b_y_list, spectra, wn_spectra, iwn, wn, 
-              spatial_axis, npix,)
-            jobs[wn] = self.job_server.submit(dirty, indata, (),
-              ('numpy','math',))
-
-        for iwn,wn in enumerate(wavenumber):
-            # collect and store results
-            dirtyimage[:,:,iwn], dirtybeam[:,:,iwn] = jobs[wn]()
-
-        self.result['dirtyimage'] = dirtyimage
-        self.result['dirtybeam'] = dirtybeam
-        self.result['spatial axis [arcsec]'] = spatial_axis
-        self.result['spatial axis'] = spatial_axis
-        self.result['wavenumber [cm-1]'] = wavenumber
-
-        return self.result
-                
-    def __repr__(self):
-        return '''
-DirtyImage:
-'''.format(
-          num_uvspectra=len(self.uvspectra))
-
Index: trunk/doublefourier.py
===================================================================
--- trunk/doublefourier.py	(revision 71)
+++ 	(revision )
@@ -1,360 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import math
-import numpy as np
-import pp
-
-import common.commonobjects as co
-
-def fftshift(data, shift):
-    """Method to shift a 2d complex data array by applying the
-       given phase shift to its Fourier transform.
-    """
-    # move centre of image to array origin
-    temp = numpy.fft.fftshift(data)
-    # 2d fft
-    temp = numpy.fft.fft2(temp)
-    # apply phase shift
-    temp *= shift
-    # transform and shift back
-    temp = numpy.fft.ifft2(temp)
-    temp = numpy.fft.fftshift(temp)
-
-    return temp
-
-
-class DoubleFourier(object):
-    """Class to compute interferograms.
-    """
-
-    def __init__(self, parameters, previous_results, job_server):
-        self.parameters = parameters
-        self.previous_results = previous_results
-        self.job_server = job_server
-
-        self.result = collections.OrderedDict()      
-
-    def run(self):
-        print 'DoubleFourier.run'
-
-        # convert lengths to m. Wavenumbers leave as cm-1.
-        fts = self.previous_results['fts']
-        fts_wn = fts['fts_wn']
-        fts_wn_truncated = fts['fts_wn_truncated']
-        opd_max = fts['opd_max'] / 100.0
-        fts_nsample = fts['ftsnsample']
-        vdrive = fts['vdrive'] / 100.0
-        delta_opd = fts['delta_opd']
-        delta_time = delta_opd / vdrive
-
-        times = np.arange(int(fts_nsample), dtype=np.float)
-        times *= delta_time
-        opd_start = -((fts_nsample + 1) / 2) * delta_opd
-
-        beamsgenerator = self.previous_results['beamsgenerator']
-
-        uvmapgenerator = self.previous_results['uvmapgenerator']
-        bxby = uvmapgenerator['bxby']
-
-        skygenerator = self.previous_results['skygenerator']
-        skymodel = skygenerator['sky model']
-        spatial_axis = self.result['spatial axis'] = skygenerator['spatial axis']
-        self.result['frequency axis'] = fts_wn_truncated
-
-        # assuming nx is even then transform has 0 freq at origin and [nx/2] is
-        # Nyquist frequency. Nyq freq = 0.5 * Nyquist sampling freq. 
-        # Assume further that the fft is shifted so that 0 freq is at nx/2
-        nx = len(spatial_axis)
-        spatial_freq_axis = np.arange(-nx/2, nx/2, dtype=np.float)
-        sample_freq = (180.0 * 3600.0 / np.pi) / (spatial_axis[1] - spatial_axis[0])
-        spatial_freq_axis *= (sample_freq / nx)
-        self.result['spatial frequency axis'] = spatial_freq_axis
-        
-        self.result['baseline interferograms'] = collections.OrderedDict()
-#        for baseline in bxby:
-        for baseline in bxby[:1]:
-            print 'baseline', baseline
-            measurement = np.zeros(np.shape(times))
-            opd = np.zeros(np.shape(times))
-
-            # FTS path diff and possibly baseline itself vary with time
-            for tindex,t in enumerate(times):
-
-                # calculate the sky that the system is observing at this
-                # moment, incorporating various errors in turn. Be explicit
-                # about the copy otherwise a reference is made instead.
-                sky_now = skymodel.copy()
-
-                # 1. baseline should be perp to centre of field.
-                #    If baseline is tilted then origin of sky map shifts.
-                #    (I think effect could be corrected by changing
-                #    FTS sample position to compensate.?)
-                # for now assume 0 error but do full calculation for timing
-                # purposes
-
-                # perfect baseline is perpendicular to direction to 'centre'
-                # on sky. Add errors in x,y,z - z towards sky 'centre'
-                bz = 0.0
-                # convert bz to angular error
-                blength = math.sqrt(baseline[0]*baseline[0] +
-                  baseline[1]*baseline[1])
-                bangle = (180.0 * 3600.0 / math.pi) * bz / blength 
-                bx_error = bangle * baseline[0] / blength
-                by_error = bangle * baseline[1] / blength
-
-                # calculate xpos, ypos in units of pixel - numpy arrays 
-                # [row,col]
-                nx = len(spatial_axis)
-                colpos = float(nx-1) * bx_error / \
-                  (spatial_axis[-1] - spatial_axis[0])
-                rowpos = float(nx-1) * by_error / \
-                  (spatial_axis[-1] - spatial_axis[0])
-
-                # calculate fourier phase shift to move point at [0,0] to 
-                # [rowpos, colpos]
-                shiftx = np.zeros([nx], np.complex)
-                shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
-                shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-                shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / float(nx))
-
-                shifty = np.zeros([nx], np.complex)
-                shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
-                shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-                shifty = np.exp((-2.0j * np.pi * rowpos * shifty) / float(nx))
-
-                shift = np.ones([nx,nx], np.complex)
-                for j in range(nx):
-                    shift[j,:] *= shiftx
-                for i in range(nx):
-                    shift[:,i] *= shifty
-
-                jobs = {}
-                # go through freq planes and shift them
-                for iwn,wn in enumerate(fts_wn_truncated):
-                    # submit jobs
-                    indata = (sky_now[:,:,iwn], shift,)
-                    jobs[wn] = self.job_server.submit(fftshift,
-                      indata, (), ('numpy',))
-                
-                for iwn,wn in enumerate(fts_wn_truncated):
-                    # collect and store results
-                    temp = jobs[wn]()
-                    sky_now[:,:,iwn] = temp
-
-                if t == times[0]:
-                    # take copy of array to fix a snapshot of it now
-                    self.result['sky at time 0'] = sky_now.copy()
-
-                # 2. telescopes should be centred on centre of field
-                #    Telescopes collect flux from the 'sky' and pass
-                #    it to the FTS beam combiner. In doing this each
-                #    telescope multiplies the sky emission by its
-                #    amplitude beam response - always real but with
-                #    negative areas. Is this correct? Gives right
-                #    answer for 'no error' case.
-
-                # multiply sky by amplitude beam 1 * amplitude beam 2
-                amp_beam_1 = beamsgenerator['primary amplitude beam'].data
-                amp_beam_2 = beamsgenerator['primary amplitude beam'].data
-
-                for iwn,wn in enumerate(fts_wn):
-                    # calculate shifted beams here
-                    # for now assume no errors and just use beam
-                    # calculated earlier
-                    pass
-
-                # multiply sky by amplitude beams of 2 antennas
-                sky_now *= amp_beam_1 * amp_beam_2 
-
-                if t == times[0]:
-                    # take copy of array to fix a snapshot of it now
-                    self.result['sky*beams at time 0'] = sky_now.copy()
-                    
-                # 3. baseline error revisited
-                # derive baseline at this time
-                #
-                # Perhaps baseline should be a continuous function in
-                # time, which would allow baselines that intentionally
-                # smoothly vary (as in rotating tethered assembly)
-                # and errors to be handled by one description.
-                #
-                # what follows assumes zero error
-
-                baseline_error = 0.0
-                baseline_now = baseline + baseline_error
-
-                fft_now = np.zeros(np.shape(sky_now), np.complex)
-                spectrum = np.zeros(np.shape(fts_wn), np.complex)
-                for iwn,wn in enumerate(fts_wn_truncated):
-
-                    # derive shift needed to place baseline at one of FFT coords
-                    # this depends on physical baseline and frequency
-                    baseline_now_lambdas = baseline_now * wn * 100.0 
-
-                    # calculate baseline position in units of pixels of FFTed 
-                    # sky - numpy arrays [row,col]
-                    colpos = float(nx-1) * \
-                      float(baseline_now_lambdas[0] - spatial_freq_axis[0]) / \
-                      (spatial_freq_axis[-1] - spatial_freq_axis[0])
-                    rowpos = float(nx-1) * \
-                      float(baseline_now_lambdas[1] - spatial_freq_axis[0]) / \
-                      (spatial_freq_axis[-1] - spatial_freq_axis[0])
-                    colpos = 0.0
-                    rowpos = 0.0
-
-                    # calculate fourier phase shift to move point at [rowpos,colpos] to 
-                    # [0,0]
-                    shiftx = np.zeros([nx], np.complex)
-                    shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
-                    shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-                    shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / \
-                      float(nx))
-
-                    shifty = np.zeros([nx], np.complex)
-                    shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
-                    shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-                    shifty = np.exp((-2.0j * np.pi * rowpos * shifty) / float(nx))
-
-                    shift = np.ones([nx,nx], np.complex)
-                    for j in range(nx):
-                        shift[j,:] *= shiftx
-                    for i in range(nx):
-                        shift[:,i] *= shifty
-
-                    # move centre of sky image to origin
-                    temp = np.fft.fftshift(sky_now[:,:,iwn])
-                    # apply phase shift
-                    temp *= shift
-                    # 2d fft
-                    temp = np.fft.fft2(temp)
-                    fft_now[:,:,iwn] = temp
-
-                    # set amp/phase at this frequency
-                    spectrum[wn==fts_wn] = temp[0,0]
-
-                if t == times[0]:
-                    self.result['skyfft at time 0'] = fft_now.copy()
-
-                    axis = co.Axis(data=fts_wn, title='wavenumber',
-                      units='cm-1')
-                    temp = co.Spectrum(data=spectrum, axis=axis,
-                      title='Detected spectrum', units='W sr-1 m-2 Hz-1')
-                    self.result['skyfft spectrum at time 0'] = temp
-
-                # 3. FTS sampling should be accurate
-                #    derive lag due to FTS path difference
-                #    0 error for now
-
-                if t == times[0]:
-                    # test interferogram
-                    # inverse fft of emission spectrum at this point
-                    reflected_spectrum = np.zeros([2*(len(spectrum)-1)],
-                      np.complex)
-                    reflected_spectrum[:len(spectrum)].real = spectrum
-                    reflected_spectrum[len(spectrum):].real = spectrum[-2:0:-1]
-                    
-                    temp = np.fft.ifft(reflected_spectrum)
-                    pos = np.fft.fftfreq(len(reflected_spectrum), d=fts_wn[1]-fts_wn[0])
-
-                    # move 0 frequency to centre of array
-                    temp = np.fft.fftshift(temp)
-                    pos = np.fft.fftshift(pos)
-
-                    axis = co.Axis(data=pos, title='path difference',
-                      units='cm')
-                    temp = co.Spectrum(data=temp, axis=axis,
-                      title='Detected interferogram', units='')
-
-                    self.result['test FTS at time 0'] = temp
-
-                mirror_error = 0.0
-                opd[tindex] = opd_start + (vdrive * t + mirror_error)
-                opd_ipos = opd[tindex] / delta_opd
-
-                # calculate shift needed to move point at opd to 0
-                nfreq = len(fts_wn)
-                ndoublefreq = 2 * (nfreq - 1)
-                shift = np.zeros([ndoublefreq], dtype=np.complex)
-                shift[:nfreq] = np.arange(nfreq, dtype=np.complex)
-                shift[nfreq:] = np.arange(nfreq, dtype=np.complex)[-2:0:-1]
-                shift = np.exp((-2.0j * np.pi * opd_ipos * shift) / \
-                  float(ndoublefreq))
-
-                # reflect spectrum about 0 to give unaliased version
-                reflected_spectrum = np.zeros([2*(len(spectrum)-1)], np.complex)
-                reflected_spectrum[:len(spectrum)].real = spectrum
-                reflected_spectrum[len(spectrum):].real = spectrum[-2:0:-1]
-
-                # apply phase shift and fft
-                reflected_spectrum *= shift
-                spectrum_fft = np.fft.ifft(reflected_spectrum)
-                measurement[tindex] = spectrum_fft[0]
-
-#                if np.array_equal(baseline, bxby[0]):
-#                    if 'skyfft check' not in self.result.keys():
-#                        self.result['skyfft check'] = {}
-#                    temp = co.Spectrum(data=spectrum_fft, 
-#                      title='Test interferogram', units='')
-#                    self.result['skyfft check'][t] = temp
-
-            axis = co.Axis(data=np.array(opd), title='path difference',
-              units='cm')
-            temp = co.Spectrum(data=measurement, axis=axis,
-              title='Detected interferogram', units='')
-            self.result['baseline interferograms'][tuple(baseline)] = \
-              temp
-
-        return self.result
-
-    def matlab_transform(self):
-        # readers should look at Izumi et al. 2006, Applied Optics, 45, 2576
-        # for theoretical background. Names of variables in the code correspond
-        # to that work.
-
-        # For now, assume 2 light collectors giving one baseline at a time.
-        interferograms = {}
-
-        for baseline in self.baselines:
-            interferogram = 0
-
-            # baseline length (cm) and position angle
-            bu = baseline[0]
-            bv = baseline[1]
-            mod_b = np.sqrt(pow(bu,2) + pow(bv,2))
-            ang_b = np.arctan2(bv, bu)
-
-            # loop over sky pixels covered by primary beam
-            nx = self.sky_s.shape[1]
-            ny = self.sky_s.shape[2]
-            for j in range(ny):
-                for i in range(nx):
-
-                    # inverse fft of emission spectrum at this point
-                    temp = np.fft.ifft(self.sky_s[:,j,i])
-
-                    # move 0 frequency to centre of array
-                    temp = np.fft.fftshift(temp)
-
-                    # length (radians) and position angle of theta vector
-                    mod_theta = np.sqrt(pow(self.sky_x[i],2) + pow(self.sky_y[j],2))
-                    ang_theta = np.arctan2(self.sky_y[j], self.sky_x[i])
-
-                    # calculate b.theta (the projection of b on theta)
-                    # and the corresponding delay in units of wavelength at 
-                    # Nyquist frequency
-                    delay = mod_theta * mod_b * np.cos(ang_b - ang_theta) * self.freqs[-1]
-
-                    # sampling is done at twice Nyquist freq so shift transformed
-                    # spectrum by 2 * delay samples (approximated to nint)
-                    # NOTE factor of 2 discrepency with matlab version! I think
-                    # this is because there the variable 'Nyq' is the Nyquist
-                    # sampling rate, not the Nyquist frequency. 
-                    temp = np.roll(temp, int(round(2.0 * delay)))
-
-                    # want only the real part of the result
-                    interferogram += np.real(temp) 
-
-    def __repr__(self):
-        return 'DoubleFourier'
-
Index: trunk/fiins.py
===================================================================
--- trunk/fiins.py	(revision 71)
+++ 	(revision )
@@ -1,93 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import pp
-
-import beamsgenerator
-import cleanimage
-import dirtyimage
-import fts
-import loadparameters
-import observe
-import renderer
-import skygenerator
-import uvmapgenerator
-import uvspectra
-
-
-class Fiins(object):
-    """FISICA simulator.
-    """
-    def __init__(self, sky_spreadsheet='SkySparams.xlsx', sky_sheet='1point'):
-        self.sky_spreadsheet = sky_spreadsheet
-        self.sky_sheet = sky_sheet
-        self.result = collections.OrderedDict()
-
-        # start parallel python (pp), find the number of CPUS available
-        ppservers = ()
-        self.job_server = pp.Server(ppservers=ppservers)
-        print 'fiins starting pp with %s workers' % \
-          self.job_server.get_ncpus()
-
-    def simulate(self):
-        # read parameters
-        loadparams = loadparameters.LoadParameters(
-          sky_spreadsheet=self.sky_spreadsheet,
-          sky_sheet=self.sky_sheet)
-        obs_specification = loadparams.run()
-        self.result['loadparameters'] = obs_specification
-
-        # generate information on the FTS 
-        ftsd = fts.FTS(parameters=obs_specification)
-        self.result['fts'] = ftsd.run()
-        print ftsd
-
-        # generate UV map and observation framework
-        uvmapgen = uvmapgenerator.UVMapGenerator(parameters=obs_specification,
-          previous_results=self.result)
-        self.result['uvmapgenerator'] = uvmapgen.run()
-        print uvmapgen
-
-        # generate primary beams
-        beamsgen = beamsgenerator.BeamsGenerator(
-          previous_results=self.result)
-        self.result['beamsgenerator'] = beamsgen.run()
-
-        # construct sky
-        skygen = skygenerator.SkyGenerator(parameters=obs_specification,
-          previous_results=self.result)
-        self.result['skygenerator'] = skygen.run()
-
-        # calculate interferograms
-        obs = observe.Observe(parameters=obs_specification,
-          previous_results=self.result, job_server=self.job_server)
-        self.result['observe'] = obs.run()  
-
-        # recover spectra from interferograms
-        uvspec = uvspectra.UVspectra(previous_results=self.result,
-          job_server=self.job_server)
-        self.result['uvspectra'] = uvspec.run()  
-
-        # construct dirty image
-        dirty = dirtyimage.DirtyImage(previous_results=self.result,
-          job_server=self.job_server)
-        self.result['dirtyimage'] = dirty.run()  
-
-        # construct clean image
-        clean = cleanimage.CleanImage(previous_results=self.result,
-          job_server=self.job_server)
-        self.result['cleanimage'] = clean.run()  
-
-        # construct html description of result
-        self.render()
-        
-    def import_result(self):
-        print 'does nothing'
-
-    def render(self):
-        htmlrenderer = renderer.Renderer(result=self.result)
-        htmlrenderer.run()
-
-    def __repr__(self):
-        return 'FISICA simulator'
-
Index: trunk/fts.py
===================================================================
--- trunk/fts.py	(revision 71)
+++ 	(revision )
@@ -1,72 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import numpy as np
-
-
-class FTS(object):
-    """Class to describe the Fourier Transform Spectrometer.
-    """
-
-    def __init__(self, parameters):
-        self.parameters = parameters
-
-        # read params from FTSpectrograph sheet
-        ftsparm = parameters['substages']['FTSpectrograph']
-
-        self.result = collections.OrderedDict()
-
-        # slightly awkward getting values as they are keyed by 
-        # the row number in the parent spreadsheet
-        row = ftsparm['wn min [cm-1]'].keys()[0]
-        self.result['wnmin'] = wnmin = ftsparm['wn min [cm-1]'][row]
-        self.result['wnmax'] = wnmax = ftsparm['wn max [cm-1]'][row]
-        self.result['specres'] = specres = ftsparm['Spectral Res'][row]
-
-        # read params from FTSMechanical sheet
-        ftsparm = parameters['substages']['FTSMechanical']
-        row = ftsparm['V drive [cm/s]'].keys()[0]
-        self.result['vdrive'] = ftsparm['V drive [cm/s]'][row]
-
-        # derived parameters
-        self.result['delta_wn'] = delta_wn = wnmin / specres
-        self.result['delta_opd'] = 1.0 / (2.0 * wnmax * 100.0)
-        # number of unaliased spectral points
-        nspec = int(np.ceil(wnmax / delta_wn)) + 1
-        # ftsnsample symmetric about 0 opd
-        self.result['ftsnsample'] = ftsnsample = 2 * (nspec - 1)
-        self.result['opd_max'] = 1.0 / (2.0 * delta_wn)
-
-        # unaliased spectral points sampled by FTS (wn=wavenumber[cm-1]
-        i = np.arange(nspec)
-        fts_wn = wnmax * i / float(nspec-1)
-
-        # wn axis truncated at wnmin
-        fts_wn_truncated = fts_wn[fts_wn >= wnmin]
-
-        self.result['fts_wn'] = fts_wn
-        self.result['fts_wn_truncated'] = fts_wn_truncated
-
-    def run(self):
-        print 'FTS.run'        
-        return self.result
-
-    def __repr__(self):
-        return '''
-FTS:
-  wnmin    : {wnmin}
-  wnmax    : {wnmax}
-  specres  : {specres}
-  delta wn : {delta_wn}
-  delta opd: {delta_opd}
-  nsample  : {nsample}
-  max opd  : {opd_max}
-'''.format(
-          wnmin=self.result['wnmin'],
-          wnmax=self.result['wnmax'],
-          specres=self.result['specres'],
-          delta_wn=self.result['delta_wn'],
-          delta_opd=self.result['delta_opd'],
-          nsample=self.result['ftsnsample'],
-          opd_max=self.result['opd_max'])
-
Index: trunk/hogbom.py
===================================================================
--- trunk/hogbom.py	(revision 71)
+++ 	(revision )
@@ -1,130 +1,0 @@
-import numpy 
-
-def image(himage, psf, dirty, niter, thresh):
-    himage.open(psf)
-    psfdata = himage.getchunk()
-    print 'psf read, max val %s' % numpy.max(psfdata)
-    psfdata /= numpy.max(psfdata)
-    print 'psf normalised'
-
-#    if not dirty: 
-#        sky = numpy.zeros(numpy.shape(psf))
-#        dirty = numpy.zeros(numpy.shape(psf))
-    
-#        sky[70,70]=5000
-#        for i in range(numpy.shape(psf)[0]):
-#            for j in range(numpy.shape(psf)[1]):
-#                if pow((i-20),2) + pow((j-20),2) < 100:
-#                    sky[i,j] = 10000
-
-#        for i in range(numpy.shape(psf)[0]):
-#            for j in range(numpy.shape(psf)[1]):
-#                a1o, a2o=overlapIndices(sky, psf,
-#                                    i+1-sky.shape[0]/2,
-#                                    j+1-sky.shape[1]/2)
-#                dirty[a1o[0]:a1o[1],a1o[2]:a1o[3]]+=psf[a2o[0]:a2o[1],a2o[2]:a2o[3]]*sky[i,j]
-#    else:
-#        himage.open(dirty)
-#        dirtydata = himage.getchunk()
-#    himage.fromarray(outfile='dirty.im', pixels=dirty, overwrite=True) 
-#    print 'dirty read'
-
-    himage.open(dirty)
-    dirtydata = himage.getchunk()
-    print 'dirty read'
-
-    comps, res = hogbom(dirtydata, psfdata, True, 0.05, thresh, niter)
-    print 'clean completed'
-
-    himage.fromarray(outfile='residual.im', pixels=res, overwrite=True) 
-    himage.fromarray(outfile='comps.im', pixels=comps, overwrite=True) 
-    himage.fromarray(outfile='image.im', pixels=comps+res, overwrite=True) 
-    himage.close()
-
-
-# Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk> 
-# Initial version August 2010
-#
-# This file is part of pydeconv. This work is licensed under GNU GPL
-# V2 (http://www.gnu.org/licenses/gpl.html)
-"""
-Clean based deconvolution, using numpy
-"""
-
-def overlapIndices(a1, a2, 
-                   shiftx, shifty):
-    if shiftx >=0:
-        a1xbeg=shiftx
-        a2xbeg=0
-        a1xend=a1.shape[0]
-        a2xend=a1.shape[0]-shiftx
-    else:
-        a1xbeg=0
-        a2xbeg=-shiftx
-        a1xend=a1.shape[0]+shiftx
-        a2xend=a1.shape[0]
-
-    if shifty >=0:
-        a1ybeg=shifty
-        a2ybeg=0
-        a1yend=a1.shape[1]
-        a2yend=a1.shape[1]-shifty
-    else:
-        a1ybeg=0
-        a2ybeg=-shifty
-        a1yend=a1.shape[1]+shifty
-        a2yend=a1.shape[1]
-
-    return (a1xbeg, a1xend, a1ybeg, a1yend), (a2xbeg, a2xend, a2ybeg, a2yend)
-
-        
-
-def hogbom(dirty,
-           psf,
-           window,
-           gain,
-           thresh,
-           niter):
-    """
-    Hogbom clean
-
-    :param dirty: The dirty image, i.e., the image to be deconvolved
-
-    :param psf: The point spread-function
-
-    :param window: Regions where clean components are allowed. If
-    True, thank all of the dirty image is assumed to be allowed for
-    clean components
-
-    :param gain: The "loop gain", i.e., the fraction of the brightest
-    pixel that is removed in each iteration
-
-    :param thresh: Cleaning stops when the maximum of the absolute
-    deviation of the residual is less than this value
-
-    :param niter: Maximum number of components to make if the
-    threshold "thresh" is not hit
-    """
-    comps=numpy.zeros(dirty.shape)
-    res=numpy.array(dirty)
-    if window is True:
-        window=numpy.ones(dirty.shape,
-                          numpy.bool)
-    for i in range(niter):
-        mx, my=numpy.unravel_index(numpy.fabs(res[window]).argmax(), res.shape)
-        print 'mx, my', mx, my
-        mval=res[mx, my]*gain
-        print 'mval', mval
-        comps[mx, my]+=mval
-        a1o, a2o=overlapIndices(dirty, psf,
-                                mx+1-dirty.shape[0]/2,
-                                my+1-dirty.shape[1]/2)
-#                                mx-dirty.shape[0]/2,
-#                                my-dirty.shape[1]/2)
-        print 'a1o', a1o
-        print 'a2o', a2o
-        res[a1o[0]:a1o[1],a1o[2]:a1o[3]]-=psf[a2o[0]:a2o[1],a2o[2]:a2o[3]]*mval
-        print '59,56', res[59,56]
-        if numpy.fabs(res).max() < thresh:
-            break
-    return comps, res
Index: trunk/loadparameters.py
===================================================================
--- trunk/loadparameters.py	(revision 71)
+++ 	(revision )
@@ -1,235 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import os
-import pprint
-import xlrd
-
-
-class LoadParameters:
-    """Class to compute interferograms.
-    """
-
-    def __init__(self, sky_spreadsheet='Sky.xlsx', sky_sheet='1point',
-      instrument_spreadsheet='FIInS_Instrument_cor3.xlsx'):
-        self.result = collections.OrderedDict()
-        self.result['substages'] = collections.OrderedDict()
-
-        # access the relevant sheet of the spreadsheet
-        book = xlrd.open_workbook(sky_spreadsheet)
-        sheet = book.sheet_by_name(sky_sheet)
-
-        sheet_dict = collections.OrderedDict()
-
-        # Sheet has names in col 0, values in columns 1 to n
-        # Ignore row 0.
-
-        # read column names
-        for row in range(1,sheet.nrows):
-            colname = sheet.cell(row, 0)
-            sheet_dict[colname.value] = {}
-            for col in range(1,sheet.ncols):
-                sheet_dict[colname.value][col] = sheet.cell(row,col).value
-
-        self.result['substages']['Sky'] = sheet_dict
-
-        # the instrument spreadsheet
-        instrument_book = xlrd.open_workbook(instrument_spreadsheet)
-        for sheet in instrument_book.sheets():
-
-            if sheet.name == 'FTSpectrograph':
-                sheet_dict = collections.OrderedDict()
-
-                # ignore first row of sheet
-                # second row of sheet has column names
-                name_row = 1
-                # subsequent rows have entries, only one of which is
-                # 'selected'
-                select_col = None
-
-                # read column names
-                for col in range(sheet.ncols):
-                    sheet_dict[sheet.cell(name_row,col).value] = \
-                      collections.OrderedDict()
-                    if 'Selection' in sheet.cell(name_row,col).value:
-                        select_col = col
-
-                # sanity check 
-                if select_col is None:
-                    raise Exception, \
-                      'No Selection column in FTSpectrograph spreadsheet'
-        
-                # read entries from first 'selected' row
-                for row in range(name_row+1,sheet.nrows):
-                    if int(sheet.cell(row,select_col).value) > 0:
-                        for col in range(sheet.ncols):
-                            sheet_dict[sheet.cell(name_row,col).value][row] = \
-                              sheet.cell(row,col).value
-                        break
-
-                # resulting dict should contain the entries from the
-                # selected spreadsheet row
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'FTSMechanical':
-                sheet_dict = collections.OrderedDict()
-
-                # Sheet has names in col 0, values in col 1.
-                # Ignore row 0.
-                name_col = 0
-                entry_col = 1
-
-                # read column names
-                for row in range(1,sheet.nrows):
-                    sheet_dict[sheet.cell(row,name_col).value] = \
-                      collections.OrderedDict()
-                    sheet_dict[sheet.cell(row,name_col).value][entry_col] = \
-                      sheet.cell(row,entry_col).value
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'Telescope':
-                sheet_dict = collections.OrderedDict()
-
-                # Sheet has names in col 0, values in col 1, units in col 3.
-                # Ignore rows 0 and 1.
-                name_col = 0
-                entry_col = 1
-
-                # read column names
-                for row in range(2,sheet.nrows):
-                    sheet_dict[sheet.cell(row,name_col).value] = \
-                      sheet.cell(row,entry_col).value
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'Interferometer':
-                sheet_dict = collections.OrderedDict()
-
-                # ignore first 4 rows of sheet
-                # 5th row of sheet has column names
-                name_row = 4
-                # subsequent rows have entries, only one of which is
-                # 'selected'
-                select_col = None
-
-                # read column names
-                for col in range(sheet.ncols):
-                    sheet_dict[sheet.cell(name_row,col).value] = \
-                      collections.OrderedDict()
-                    if 'Select' in sheet.cell(name_row,col).value:
-                        select_col = col
-
-                # sanity check 
-                if select_col is None:
-                    raise Exception, \
-                      'No Selection column in Interferometer spreadsheet'
-        
-                # read entries from first 'selected' row
-                for row in range(name_row+1,sheet.nrows):
-                    select_cell = str(sheet.cell(row,select_col).value)
-                    select_val = 0
-                    try:
-                        select_val = int(float(select_cell))
-                    except:
-                        pass
-           
-                    if select_val == 1:
-                        for col in range(sheet.ncols):
-                            sheet_dict[sheet.cell(name_row,col).value][row] = \
-                              sheet.cell(row,col).value
-                        break
-
-                # resulting dict should contain the entries from the
-                # selected spreadsheet row
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'ColdOptics':
-                sheet_dict = collections.OrderedDict()
-
-                # Sheet has names in col 0, values in col 1.
-                # Ignore rows 0 and 1 and final row.
-                # Ignore cols 2 and higher - duplicate information.
-                name_col = 0
-                entry_col = 1
-
-                # read column names
-                for row in range(2,sheet.nrows-1):
-                    sheet_dict[sheet.cell(row,name_col).value] = \
-                      collections.OrderedDict()
-                    sheet_dict[sheet.cell(row,name_col).value][entry_col] = \
-                      sheet.cell(row,entry_col).value
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'Background':
-                sheet_dict = collections.OrderedDict()
-
-                # Sheet has names in col 0, values in col 1.
-                # Ignore rows 0 and 1.
-                name_col = 0
-                entry_col = 1
-
-                for row in range(2,sheet.nrows):
-                    sheet_dict[sheet.cell(row,name_col).value] = \
-                      collections.OrderedDict()
-                    sheet_dict[sheet.cell(row,name_col).value][entry_col] = \
-                      sheet.cell(row,entry_col).value
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'WarmOptics':
-                sheet_dict = collections.OrderedDict()
-
-                # Sheet has names in col 0, values in col 1.
-                # Ignore rows 0 and 1.
-                name_col = 0
-                entry_col = 1
-
-                for row in range(2,sheet.nrows):
-                    sheet_dict[sheet.cell(row,name_col).value] = \
-                      collections.OrderedDict()
-                    sheet_dict[sheet.cell(row,name_col).value][entry_col] = \
-                      sheet.cell(row,entry_col).value
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'Detectors':
-                sheet_dict = collections.OrderedDict()
-
-                # Sheet has names in col 0, values in col 1.
-                # Ignore rows 0 and 1.
-                name_col = 0
-                entry_col = 1
-
-                for row in range(2,sheet.nrows):
-                    sheet_dict[sheet.cell(row,name_col).value] = \
-                      collections.OrderedDict()
-                    sheet_dict[sheet.cell(row,name_col).value][entry_col] = \
-                      sheet.cell(row,entry_col).value
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            elif sheet.name == 'SimulatorControl':
-                sheet_dict = collections.OrderedDict()
-
-                # Sheet has names in col 0, values in col 1.
-                # Ignore row 0.
-                name_col = 0
-                entry_col = 1
-
-                for row in range(1,sheet.nrows):
-                    sheet_dict[sheet.cell(row,name_col).value] = \
-                      collections.OrderedDict()
-                    sheet_dict[sheet.cell(row,name_col).value][entry_col] = \
-                      sheet.cell(row,entry_col).value
-
-                self.result['substages'][sheet.name] = sheet_dict
-
-            else:
-                print 'loadparameters skipping', sheet.name
-
-    def run(self):
-        return self.result
Index: trunk/observe.py
===================================================================
--- trunk/observe.py	(revision 71)
+++ 	(revision )
@@ -1,297 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import math
-import numpy as np
-import pp
-
-import matplotlib.pyplot as plt
-
-import common.commonobjects as co
-
-def fftshift(data, shift):
-    """Method to shift a 2d complex data array by applying the
-       given phase shift to its Fourier transform.
-    """
-    # move centre of image to array origin
-    temp = numpy.fft.fftshift(data)
-    # 2d fft
-    temp = numpy.fft.fft2(temp)
-    # apply phase shift
-    temp *= shift
-    # transform and shift back
-    temp = numpy.fft.ifft2(temp)
-    temp = numpy.fft.fftshift(temp)
-
-    return temp
-
-def fields_match(previous_config, config, exclude_fields=[]):
-    """Function to compare 2 configs, ignoring any differences in
-    the components listed in 'exclude_fields'.
-    """
-    if previous_config is None:
-        return False
-
-    if not exclude_fields:
-        return previous_config==config
-    else:
-        # convert to dicts and eliminate exclude_fields
-        previous_dict = previous_config._asdict()
-        current_dict = config._asdict()
-        for field in exclude_fields:
-            ignore = previous_dict.pop(field, None)
-            ignore = current_dict.pop(field, None)
-
-        previous_set = set([k for k in previous_dict.items()])
-        current_set = set([k for k in current_dict.items()])
-        diff = previous_set.symmetric_difference(current_set)
-
-        debug = False
-        if(debug):
-            print diff
-
-        return not(bool(diff))
-
-
-class Observe(object):
-    """Class to compute interferograms.
-    """
-
-    def __init__(self, parameters, previous_results, job_server):
-        self.parameters = parameters
-        self.previous_results = previous_results
-        self.job_server = job_server
-
-        self.result = collections.OrderedDict()      
-
-    def run(self):
-        print 'Observe.run'
-
-        # access primary beam information
-        beamsgenerator = self.previous_results['beamsgenerator']
-
-        # access required FTS information
-        fts = self.previous_results['fts']
-        fts_wn = fts['fts_wn']
-        fts_wn_truncated = fts['fts_wn_truncated']
-        delta_opd = fts['delta_opd']
-
-        # access to model sky
-        skygenerator = self.previous_results['skygenerator']
-        sky_model = skygenerator['sky model']
-        spatial_axis = self.result['spatial axis'] = skygenerator['spatial axis']
-        nx = len(spatial_axis)
-
-        # assuming nx is even then transform of spatial axis has 0 freq at origin 
-        # and [nx/2] is Nyquist frequency [Nyq freq = 0.5 * Nyquist sampling freq].
-        # The fft is shifted so that 0 freq is at nx/2.
-        nx = len(spatial_axis)
-        spatial_freq_axis = np.arange(-nx/2, nx/2, dtype=np.float)
-        sample_freq = (180.0 * 3600.0 / np.pi) / (spatial_axis[1] - spatial_axis[0])
-        spatial_freq_axis *= (sample_freq / nx)
-        self.result['spatial frequency axis'] = spatial_freq_axis
-
-        # get list of instrument configuratons       
-        uvmapgenerator = self.previous_results['uvmapgenerator']
-        obs_framework = uvmapgenerator['obs_framework']
-
-        # and calculate the result for each configuration
-        previous_config = None
-        observed_obs_framework = []
-
-        debug_plotted = False
-        for config in obs_framework:
-#            if config.fts_start:
-#                print 'start FTS scan'
-#                print 'baseline', config.baseline_x, config.baseline_y,\
-#                  config.baseline_z
-  
-            if fields_match(previous_config, config,
-              exclude_fields=['scan_number',
-                              'fts_start',
-                              'fts_position',
-                              'fts_nominal_position',
-                              'baseline_x',
-                              'baseline_y',
-                              'baseline_z',
-                              'baseline_number',
-                              'data']):
-                pass
-            else:
-#                print 'previous', previous_config
-#                print 'current', config
-                # calculate new sky
-#                print 'calculating new sky'
-
-                # Be explicit about copies otherwise a reference is
-                # made which corrupts the original and leads to much
-                # confusion.
-                sky_now = sky_model.copy()
-
-                # calculate the sky that the system is observing at this
-                # moment, incorporating various errors in turn.
-
-                # 1. baseline should be perp to vector towards centre
-                #    of field. If baseline is tilted then origin of sky 
-                #    map shifts.
-                #    I think this error can be folded into that for the
-                #    FTS mirror, so leave handling it until then. Is
-                #    this correct?
-                
-                # 2. Flux-collecting telescopes should be pointing 
-                #    at centre of field. They each collect flux from 
-                #    the 'sky' and pass it to the FTS beam combiner.
-                #    In doing this each telescope multiplies the sky 
-                #    emission by its complex amplitude beam response.
-                #    So instead of measuring the correlation of the
-                #    'sky' with itself, we measure that x Beam1.Beam2*. 
-                #    Is this correct? Gives right answer for 'no error'
-                #    where Beam1.Beam2* == PSF.
-
-                # multiply sky by amp.beam 1 * conj(amp.beam 2)
-                amp_beam_1 = beamsgenerator['primary amplitude beam'].data
-                amp_beam_2 = beamsgenerator['primary amplitude beam'].data
-
-                # shift beams if pointing error
-                if config.point1_x_error or config.point1_y_error:
-                    # calculate shifted beam1
-#                    raise Exception, 'not implemented'
-                    current_amp_beam_1 = amp_beam_1
-                else:
-                    current_amp_beam_1 = amp_beam_1
-
-                if config.point2_x_error or config.point2_y_error:
-                    # calculate shifted beam2
-#                    raise Exception, 'not implemented'
-                    current_amp_beam_2 = amp_beam_2
-                else:
-                    current_amp_beam_2 = amp_beam_2
- 
-                sky_now *= current_amp_beam_1 * \
-                  np.conj(current_amp_beam_2)            
-
-            # 3. Get result for baseline at this time
-            if fields_match(previous_config, config,
-              exclude_fields=['fts_start',
-                              'fts_position',
-                              'fts_nominal_position',
-                              'data']):
-                pass
-            else:
-                # calculate new baseline spectrum
-#                print 'calculating new baseline spectrum'
- 
-                baseline = np.array([config.baseline_x, config.baseline_y])
-                fft_now = np.zeros(np.shape(sky_now), np.complex)
-                baseline_spectrum = np.zeros(np.shape(fts_wn), np.complex)
-
-                for iwn,wn in enumerate(fts_wn_truncated):
-
-                    # derive shift needed to place baseline at one of FFT
-                    # coords this depends on physical baseline and frequency
-                    baseline_lambda = baseline * wn * 100.0
-
-                    # calculate baseline position in units of pixels of 
-                    # FFTed sky - numpy arrays [row,col]
-                    colpos = float(nx-1) * baseline_lambda[0] / \
-                      (spatial_freq_axis[-1] - spatial_freq_axis[0])
-                    rowpos = float(nx-1) * baseline_lambda[1] / \
-                      (spatial_freq_axis[-1] - spatial_freq_axis[0])
-
-                    # calculate fourier phase shift to move point at 
-                    # [rowpos,colpos] to [0,0]
-                    shiftx = np.zeros([nx], np.complex)
-                    shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
-                    shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-                    shiftx = np.exp((2.0j * np.pi * colpos * shiftx) / \
-                      float(nx))
-
-                    shifty = np.zeros([nx], np.complex)
-                    shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
-                    shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-                    shifty = np.exp((2.0j * np.pi * rowpos * shifty) / \
-                      float(nx))
-
-                    shift = np.ones([nx,nx], np.complex)
-                    for j in range(nx):
-                        shift[j,:] *= shiftx
-                    for i in range(nx):
-                        shift[:,i] *= shifty
-
-                    # move centre of sky image to origin
-                    temp = np.fft.fftshift(sky_now[:,:,iwn])
-                    # apply phase shift
-                    temp *= shift
-                    # 2d fft
-                    temp = np.fft.fft2(temp)
-                    fft_now[:,:,iwn] = temp
-
-                    if wn==30.0 and not debug_plotted:
-
-                        freq = np.fft.fftfreq(n=np.shape(fft_now)[0],
-                          d=(spatial_axis[1]-spatial_axis[0]) / 206265.0)
-                        freq = np.arange(np.shape(fft_now)[0]) * (freq[1] - freq[0])
-
-                        ifreq = np.arange(np.shape(fft_now)[0])
-                        baseline_i = ifreq[abs(freq-abs(baseline_lambda[0])) < (freq[1]-freq[0])]
-                        baseline_j = ifreq[abs(freq-abs(baseline_lambda[1])) < (freq[1]-freq[0])]
-#                        print baseline, fft_now[baseline_i[0], baseline_j[0], iwn]
-
-#                        plt.figure()
-#                        plt.imshow(fft_now.real[:,:,iwn], interpolation='nearest', origin='lower',
-#                          aspect='equal', extent=[freq[0],freq[-1],freq[0],freq[-1]])
-#                        plt.colorbar(orientation='vertical')
-#                        plt.axis('image')
-#                        plt.savefig('debug_fft_real.png')
-#                        plt.close()
-
-#                        plt.figure()
-#                        plt.imshow(fft_now.imag[:,:,iwn], interpolation='nearest', origin='lower',
-#                          aspect='equal', extent=[freq[0],freq[-1],freq[0],freq[-1]])
-#                        plt.colorbar(orientation='vertical')
-#                        plt.axis('image')
-#                        plt.savefig('debug_fft_imag.png')
-#                        plt.close()
-
-#                        debug_plotted = True
-                  
-                    # set amp/phase at this frequency
-                    baseline_spectrum[wn==fts_wn] = temp[0,0]
-
-            # 4. Calculate interferogram value from this baseline on the 'sky'
-            #    at the current FTS mirror position.
-
-            opd = 2.0 * config.fts_position
-            opd_ipos = opd / delta_opd
-
-            # calculate shift needed to move point at opd to 0
-            nfreq = len(fts_wn)
-            ndoublefreq = 2 * (nfreq - 1)
-            shift = np.zeros([ndoublefreq], dtype=np.complex)
-            shift[:nfreq] = np.arange(nfreq, dtype=np.complex)
-            shift[nfreq:] = np.arange(nfreq, dtype=np.complex)[-2:0:-1]
-            shift = np.exp((-2.0j * np.pi * opd_ipos * shift) /\
-              float(ndoublefreq))
-
-            # reflect spectrum about 0 to give unaliased version
-            reflected_spectrum = np.zeros([2*(len(baseline_spectrum)-1)],
-              np.complex)
-            reflected_spectrum[:len(baseline_spectrum)] = baseline_spectrum
-            reflected_spectrum[len(baseline_spectrum):] = \
-              baseline_spectrum[-2:0:-1]
-
-            # apply phase shift and fft
-            reflected_spectrum *= shift
-            spectrum_fft = np.fft.ifft(reflected_spectrum)
-            config = config._replace(data = spectrum_fft[0].real)
-      
-            observed_obs_framework.append(config)
-
-            previous_config = config
-
-        self.result['observed_framework'] = observed_obs_framework
-        return self.result
-
-    def __repr__(self):
-        return 'Observe'
-
Index: trunk/pointingerrors.py
===================================================================
--- trunk/pointingerrors.py	(revision 71)
+++ 	(revision )
@@ -1,38 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import numpy as np
-
-
-class PointingErrors(object):
-    """Class to generate simulated pointing errors.
-    """
-
-    def __init__(self, parameters):
-        self.parameters = parameters
-        self.result = collections.OrderedDict()
-
-    def run(self):
-        print 'PointingErrors.run'
-        control = self.parameters['substages']['SimulatorControl']
-
-        # row number in the spreadsheet
-        row = control['PointingErrors'].keys()[0]
-        pointingerrors = control['PointingErrors'][row]
-        pointingerrors = int(round(pointingerrors))
-
-        if pointingerrors:
-            print 'ERROR: pointing errors not yet implemented'
-        else:
-            # dummy nbtot for now
-            nbtot = 10
-            pointing_err_rad = np.zeros([nbtot,2])
-            pointing_err_xy = np.zeros([nbtot,2])
-            pointing_err_t1 = np.zeros([nbtot,2])
-            pointing_err_t2 = np.zeros([nbtot,2])
-
-        return self.result
-
-    def __repr__(self):
-        return 'PointingErrors'
-
Index: trunk/pythonclean.py
===================================================================
--- trunk/pythonclean.py	(revision 71)
+++ 	(revision )
@@ -1,88 +1,0 @@
-#import numpy
-
-# Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk> 
-# Initial version August 2010
-#
-# This file is part of pydeconv. This work is licensed under GNU GPL
-# V2 (http://www.gnu.org/licenses/gpl.html)
-"""
-Clean based deconvolution, using numpy
-"""
-
-# Note added by jfl. I think a1 and a2 are assumed the same shape.
-
-def overlapIndices(a1, a2, 
-                   shiftx, shifty):
-    if shiftx >=0:
-        a1xbeg=shiftx
-        a2xbeg=0
-        a1xend=a1.shape[0]
-        a2xend=a1.shape[0]-shiftx
-    else:
-        a1xbeg=0
-        a2xbeg=-shiftx
-        a1xend=a1.shape[0]+shiftx
-        a2xend=a1.shape[0]
-
-    if shifty >=0:
-        a1ybeg=shifty
-        a2ybeg=0
-        a1yend=a1.shape[1]
-        a2yend=a1.shape[1]-shifty
-    else:
-        a1ybeg=0
-        a2ybeg=-shifty
-        a1yend=a1.shape[1]+shifty
-        a2yend=a1.shape[1]
-
-    return (a1xbeg, a1xend, a1ybeg, a1yend), (a2xbeg, a2xend, a2ybeg, a2yend)
-
-        
-
-def hogbom(dirty,
-           psf,
-           window,
-           gain,
-           thresh,
-           niter):
-    """
-    Hogbom clean
-
-    :param dirty: The dirty image, i.e., the image to be deconvolved
-
-    :param psf: The point spread-function
-
-    :param window: Regions where clean components are allowed. If
-    True, thank all of the dirty image is assumed to be allowed for
-    clean components
-
-    :param gain: The "loop gain", i.e., the fraction of the brightest
-    pixel that is removed in each iteration
-
-    :param thresh: Cleaning stops when the maximum of the absolute
-    deviation of the residual is less than this value
-
-    :param niter: Maximum number of components to make if the
-    threshold "thresh" is not hit
-    """
-    comps=numpy.zeros(dirty.shape)
-    res=numpy.array(dirty)
-    grid=numpy.indices(dirty.shape)
-    if window is True:
-        window = [slice(dirty.shape), slice(dirty.shape)]
-
-    for i in range(niter):
-        argmax=numpy.unravel_index(numpy.fabs(res[window]).argmax(), res[window].shape)
-        mx=grid[0][window][argmax]
-        my=grid[1][window][argmax]
-
-        mval=res[mx, my]*gain
-        comps[mx, my]+=mval
-        a1o, a2o=pythonclean.overlapIndices(dirty, psf,
-                                mx-dirty.shape[0]/2,
-                                my-dirty.shape[1]/2)
-        res[a1o[0]:a1o[1],a1o[2]:a1o[3]]-=psf[a2o[0]:a2o[1],a2o[2]:a2o[3]]*mval
-        if numpy.fabs(res).max() < thresh:
-            break
-    return comps, res
-
Index: trunk/renderer.py
===================================================================
--- trunk/renderer.py	(revision 71)
+++ 	(revision )
@@ -1,81 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import datetime
-import os
-import shutil
-
-from extern import mako
-from mako.template import Template
-from mako.lookup import TemplateLookup
-
-import templates
-import templates.resources as resources
-
-
-class Renderer(object):
-    """Class to render results into html.
-    """
-
-    def __init__(self, result):
-        self.result = result
-
-    def copy_resources(self, dirname):
-        outdir = os.path.join(dirname, 'resources')
-
-        # shutil.copytree complains if the output directory exists
-        if os.path.exists(outdir):
-            shutil.rmtree(outdir)
-
-        # copy all uncompressed non-python resources to output directory
-        src = os.path.dirname(resources.__file__)
-        dst = outdir
-        ignore_fn = shutil.ignore_patterns('*.zip','*.py','*.pyc','CVS*')
-        shutil.copytree(src,
-                        dst,
-                        symlinks=False,
-                        ignore=ignore_fn)
-
-#        # unzip fancybox to output directory
-#        infile = os.path.join(src, 'fancybox.zip')
-#        z = zipfile.ZipFile(infile, 'r')
-#        z.extractall(outdir)
-
-    def run(self):
-        # make a directory to contain this result
-        now = datetime.datetime.today()        
-        dirname = 'fisica-%s' % (now.strftime('%Y%m%dT%H%M%S'))
-        os.mkdir(dirname)
-
-        # make a resources directory and copy the Bootstrap stuff there
-        self.copy_resources(dirname)
-
-        # where to find the template files
-        templates_dir = os.path.dirname(templates.__file__)
-        mylookup = TemplateLookup(directories=[templates_dir])
-
-        # render the simulation steps first
-        for stage in self.result.items():
-            if 'substages' in stage[1].keys():
-                for substage in stage[1]['substages'].items():
-                    print 'rendering', stage[0], substage[0]
-                    context = {'dirname':dirname, 'renderstage':stage[0],
-                      'rendersubstage':substage[0], 'data':self.result}
-
-                    link = os.path.join(context['dirname'], '%s-%s.html' % (stage[0], substage[0]))
-                    with open(link, 'w') as f:
-                        template = mylookup.get_template('%s-%s.html' % (stage[0], substage[0]))
-                        f.write(template.render(**context))
-
-            else:
-                # no substages
-                print 'rendering', stage[0]
-                context = {'mylookup':mylookup, 'dirname':dirname,
-                  'renderstage':stage[0],
-                  'data':self.result}
-
-                link = os.path.join(context['dirname'], '%s.html' % (stage[0]))
-                with open(link, 'w') as f:
-                    template = mylookup.get_template('%s.html' % (stage[0]))
-                    f.write(template.render(**context))
-
Index: trunk/skygenerator.py
===================================================================
--- trunk/skygenerator.py	(revision 71)
+++ 	(revision )
@@ -1,262 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import math
-import numpy as np
-import scipy.constants as sc
-
-import common.commonobjects as co
-
-def black_body(temperature, wavenumber):
-    """Function to calculate Planck function. 
-       temperature - Kelvin
-       wavenumber - frequency in cm-1
-    """
-    freq = wavenumber * sc.c * 100.0
-    if freq > 0:
-        jnu = 2.0 * sc.h * pow(freq,3) / (pow(sc.c,2) *
-          (math.exp((sc.h * freq) / (sc.k * temperature)) - 1.0))
-    else:
-        jnu = 0.0
-
-    return jnu
-
-class BB_spectrum(object):
-    """Class to generate BB spectrum.
-    """
-    def __init__(self, temperature, frequency_axis, cutoffmin, cutoffmax,
-      emissivity):
-        self.temperature = temperature
-        self.frequency_axis = frequency_axis
-        self.cutoffmin = cutoffmin
-        self.cutoffmax = cutoffmax
-        self.emissivity = emissivity
-
-    def calculate(self):
-        spectrum = np.zeros(np.shape(self.frequency_axis))
-        # ignore floating point errors
-        old_settings = np.seterr(all='ignore')
-
-        for iwn,wn in enumerate(self.frequency_axis):
-            # simulate real-life 'rounded' cutoffs numerically  
-            f1 = 1.0 / (1.0 + pow(self.cutoffmin/wn, 18) + 
-              pow(wn/self.cutoffmax, 24))
-            spectrum[iwn] = black_body(self.temperature, wn) * f1 * \
-              self.emissivity
-
-        # restore fp behaviour
-        ignore = np.seterr(**old_settings)
-
-        return spectrum
-
-
-class linespectrum(object):
-    """Class to return spectrum with line at wn.
-    """
-    def __init__(self, temperature, linefreq, frequency_axis):
-        self.temperature = temperature
-        self.linefreq = linefreq
-        self.frequency_axis = frequency_axis
-
-    def calculate(self):
-        # spectrum is a sync function centred on the specified line freq
-        x = (self.frequency_axis - self.linefreq) / \
-          (self.frequency_axis[1] - self.frequency_axis[0])
-        spectrum = np.sinc(x) * self.temperature
-
-        return spectrum
-
-
-class SkyGenerator(object):
-    """Class to generate a model sky.
-    """
-
-    def __init__(self, parameters, previous_results):
-        self.parameters = parameters
-        self.previous_results = previous_results
-        self.result = collections.OrderedDict()
-
-    def run(self):
-        print 'SkyGenerator.run'
-
-        fts = self.previous_results['fts']
-        fts_wn_truncated = fts['fts_wn_truncated']
-        cutoffmin = fts['wnmin']
-        cutoffmax = fts['wnmax']
-
-        beamsgenerator = self.previous_results['beamsgenerator']
-        npix = beamsgenerator['npix']
-        spatial_axis = beamsgenerator['spatial axis [arcsec]']
-
-        # skymodel is complex so that its fft can hold truncated version
-        # of infinitesimally sampled map - does that make sense?
-        skymodel = np.zeros([npix, npix, len(fts_wn_truncated)], np.complex)
-
-        sky = self.parameters['substages']['Sky']
-        columns = sky['sourcenum'].keys()
-
-        self.result['sources'] = collections.OrderedDict()
-
-        for column in columns:
-            temp = sky['sourcenum'][column]
-            sourcenum = int(round(temp))
-            if sourcenum not in self.result['sources'].keys():
-                self.result['sources'][sourcenum] = {}
-
-            type = sky['type'][column]
-
-            temp = sky['x pos [asec]'][column]
-            xpos = float(temp)
-
-            temp = sky['y pos [asec]'][column]
-            ypos = float(temp)
-
-            temp = sky['xwidth'][column]
-            try:
-                xwidth = float(temp)
-            except:
-                xwidth = None
-
-            temp = sky['ywidth'][column]
-            try:
-                ywidth = float(temp)
-            except:
-                ywidth = None
-
-            spectrum = sky['spectrum'][column]
-
-            temp = sky['temperature'][column]
-            temperature = float(temp)
-
-            temp = sky['linefreq'][column]
-            try:
-                linefreq = float(temp)
-            except:
-                linefreq = None
-
-            temp = sky['emissivity'][column]
-            emissivity = float(temp)
-
-            print 'generating source:%s type:%s xpos:%s ypos:%s' % (sourcenum, 
-              type, xpos, ypos)
-            if type.upper().strip() == 'GAUSSIAN':
-                print '  fwhmx:%s fwhmy:%s' % (xwidth, ywidth)
-
-            if spectrum.upper().strip() == 'BLACKBODY':
-                print '  blackbody spectrum temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
-                  temperature, cutoffmin, cutoffmax, emissivity)
-                spectrum_func = BB_spectrum(temperature, fts_wn_truncated,
-                  cutoffmin, cutoffmax, emissivity)
-
-            elif spectrum.upper().strip() == 'LINE':
-                print '  line spectrum brightness temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
-                  temperature, cutoffmin, cutoffmax, emissivity)
-                spectrum_func = linespectrum(temperature, linefreq, fts_wn_truncated)
-
-            if type.upper().strip() == 'POINT':
-                source_spectrum = self._create_point_source(xpos, ypos,
-                  skymodel, spatial_axis, fts_wn_truncated, spectrum_func)
-            elif type.upper().strip() == 'GAUSSIAN':
-                source_spectrum = self._create_gaussian_source(xpos, ypos,
-                  xwidth, ywidth, skymodel, spatial_axis, fts_wn_truncated,
-                  spectrum_func)
-            else:
-                source_spectrum = None
-                print "source type '%s' not yet implemented" % type
-
-            self.result['sources'][sourcenum]['spectrum'] = source_spectrum
-
-        self.result['sky model'] = skymodel
-        self.result['spatial axis'] = spatial_axis
-        self.result['frequency axis'] = fts_wn_truncated 
-
-        return self.result
-
-    def _create_gaussian_source(self, xpos, ypos, xwidth, ywidth,
-      skymodel, spatial_axis, frequency_axis, spectrum_function):
-        """Create a point source.
-        """
-        # calculate spectrum
-        spectrum = spectrum_function.calculate()
-
-        # calculate Gaussian profile - a cunning use of array slicing found
-        # on the web
-        x = spatial_axis
-        y = x[:,np.newaxis]
-        profile = np.exp(-4*np.log(2) * ((x-xpos)**2/xwidth**2 + (y-ypos)**2/ywidth**2))
-
-#        # apodise
-#        profile[((x-xpos)**2 + (y-ypos)**2) < 4] = 1.0
-
-        # go through freq planes and add source to each
-        for iwn,wn in enumerate(frequency_axis):
-            # add to sky model
-            skymodel[:,:,iwn] += profile * spectrum[iwn]
-
-        # return spectrum
-        axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
-        spectrum = co.Spectrum(data=spectrum, axis=axis,
-          title='Source spectrum', units='W sr-1 m-2 Hz-1')
-
-        return spectrum
-
-    def _create_point_source(self, xpos, ypos, skymodel, spatial_axis,
-      frequency_axis, spectrum_function):
-        """Create a point source.
-        """
-
-        # calculate xpos, ypos in units of pixel - numpy arrays [row,col]
-        nx = len(spatial_axis)
-        colpos = float(nx-1) * float (xpos - spatial_axis[0]) / (spatial_axis[-1] - spatial_axis[0])
-        rowpos = float(nx-1) * float (ypos - spatial_axis[0]) / (spatial_axis[-1] - spatial_axis[0])
-
-        if colpos < 0 or colpos > (nx-1) or rowpos < 0 or rowpos > (nx-1):
-            # point source is outside modelled area 
-            return
-
-        # calculate fourier phase shift to move point at [0,0] to 
-        # [rowpos, colpos]
-        shiftx = np.zeros([nx], np.complex)
-        shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
-        shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-        shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / float(nx))
-
-        shifty = np.zeros([nx], np.complex)
-        shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
-        shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
-        shifty = np.exp((-2.0j * np.pi * rowpos * shifty) / float(nx))
-
-        shift = np.ones([nx,nx], np.complex)
-        for j in range(nx):
-            shift[j,:] *= shiftx
-        for i in range(nx):
-            shift[:,i] *= shifty
-
-        # calculate spectrum
-        spectrum = spectrum_function.calculate()
-
-        # go through freq planes and add point source to each
-        for iwn,wn in enumerate(frequency_axis):
-            # create point in frequency space
-            temp = np.zeros([nx,nx])
-            temp[0,0] = spectrum[iwn]
-            # 2d fft
-            temp = np.fft.fft2(temp)
-            # apply phase shift to move point to required offset
-            temp *= shift
-            # transform back
-            temp = np.fft.ifft2(temp)
-
-            # add to sky model
-            skymodel[:,:,iwn] += temp
-
-        # return spectrum
-        axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
-        spectrum = co.Spectrum(data=spectrum, axis=axis,
-          title='Source spectrum', units='W sr-1 m-2 Hz-1')
-
-        return spectrum
-
-    def __repr__(self):
-        return 'SkyGenerator'
-
Index: trunk/uvmapgenerator.py
===================================================================
--- trunk/uvmapgenerator.py	(revision 71)
+++ 	(revision )
@@ -1,138 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import numpy as np
-
-import pointingerrorgenerators
-
-
-class UVMapGenerator(object):
-    """Class to generate the UV map of the simulated observation.
-    """
-
-    def __init__(self, parameters, previous_results):
-        self.parameters = parameters
-        self.previous_results = previous_results
-        self.result = collections.OrderedDict()
-
-    def run(self):
-        print 'UVMapGenerator.run'
-
-        # The FTS mirror positions per scan. This is found from the 
-        # FTS results that should have already been calculated.
-        fts = self.previous_results['fts']
-        self.result['ftsnsample'] = ftsnsample = fts['ftsnsample']
-        delta_opd = fts['delta_opd']
-
-        interferometer = self.parameters['substages']['Interferometer']
-
-        # slightly awkward getting pattern value as it's keyed by its
-        # row number in the spreadsheet
-        row = interferometer['Pattern'].keys()[0]
-        pattern = interferometer['Pattern'][row]
-
-        if pattern.lower() == 'spiral':
-            bmax = interferometer['bmax [m]'][row]
-            bmin = interferometer['bmin [m]'][row]
-            n_baselines = int(interferometer['Num Baselines'][row])
-            bstep = interferometer['bstep [m]'][row]
-
-            # baseline increases by bstep for each circuit of spiral
-            n_laps = (bmax - bmin) / bstep
-
-#            n_baselines = 50
-            bxby = np.zeros([n_baselines,2])     
-#            bxby[0,0] = 10.0
-#            bxby[0,1] = 0.0
-#            bxby[1,0] = 20.0
-#            bxby[1,1] = 0.0
-#            bxby[2,0] = 30.0
-#            bxby[2,1] = 0.0
-#            bxby[3,0] = 40.0
-#            bxby[3,1] = 0.0
-#            bxby[4,0] = 0.0
-#            bxby[4,1] = 10.0
-#            bxby[5,0] = 0.0
-#            bxby[5,1] = 20.0
-#            bxby[6,0] = 0.0
-#            bxby[6,1] = 30.0
-#            bxby[7,0] = 0.0
-#            bxby[7,1] = 40.0
-#            bxby[8,0] = 40.0
-#            bxby[8,1] = 40.0
-#            bxby[9,0] = 20.0
-#            bxby[9,1] = 20.0
-#            bxby[10,0] = 30.0
-#            bxby[10,1] = 30.0
-
-
-            # n_baselines is total number of points along spiral
-            for ib in range(n_baselines):
-                # phi is angle around spiral, r its radius
-                phi = n_laps * 2 * np.pi * (n_baselines-ib) / n_baselines
-                r = bmin + bstep * phi / (2 * np.pi)
-
-                bxby[ib,:] = [r * np.cos(phi), r * np.sin(phi)]
-
-            self.result['pattern'] = pattern
-            self.result['bxby'] = bxby
-            self.result['bmin'] = bmin
-            self.result['bmax'] = bmax
-            self.result['n_baselines'] = n_baselines
-
-            # Construct the obs framework that will give one FTS
-            # scan at each baseline.
-     
-            # The framework will hold the baseline and FTS mirror position
-            # for each time that the detector is read. 
-
-            # The observation consists of a series of FTS scans. The 
-            # instrument baseline, pointing etc. can change continuously
-            # as the scan is performed.
- 
-            obs_framework = []
-            Config = collections.namedtuple('Config', [
-              'scan_number',
-              'baseline_x', 'baseline_y', 'baseline_z', 'baseline_number',
-              'fts_position', 'fts_nominal_position', 'fts_start',
-              'point1_x_error', 'point1_y_error',
-              'point2_x_error', 'point2_y_error',
-              'data'],
-              verbose=False)
-
-            for ib in range(n_baselines):
-                opd_start = -((ftsnsample + 1) / 2) * delta_opd
-
-                fts_start = True
-                for ifts in np.arange(ftsnsample):
-                    opd = opd_start + float(ifts) * delta_opd
-                    fts_mirror = opd / 2.0
-                    config = Config(ib, bxby[ib][0], bxby[ib][1], 0.0, ib,
-                      fts_mirror, fts_mirror, fts_start, 0.0, 0.0, 0.0,
-                      0.0, None)
-                    obs_framework.append(config) 
-                    fts_start = False
-
-            png = pointingerrorgenerators.GaussianPointingErrors(mu0=1.0,
-              sigma0=0.0, mu1=0.0, sigma1=0.0)
-            png.fill(obs_framework)
-
-            self.result['obs_framework'] = obs_framework
-
-        else:
-            raise Exception, 'unknown baseline pattern: %s' % pattern
-
-        return self.result
-
-    def __repr__(self):
-        return '''
-UVMapGenerator:
-  uv pattern      : {pattern}
-  FTS samples/scan: {ftsnsample}
-  obs.framework length: {len_framework}
-'''.format(
-          pattern=self.result['pattern'],
-          ftsnsample=self.result['ftsnsample'],
-          len_framework=len(self.result['obs_framework']))
-
-
Index: trunk/uvspectra.py
===================================================================
--- trunk/uvspectra.py	(revision 71)
+++ 	(revision )
@@ -1,102 +1,0 @@
-from __future__ import absolute_import
-
-import collections
-import math
-import numpy as np
-import pp
-
-import common.commonobjects as co
-
-
-class UVspectra(object):
-    """Class to compute spectra from each interferogram.
-    """
-
-    def __init__(self, previous_results, job_server):
-        self.previous_results = previous_results
-        self.job_server = job_server
-
-        self.uvspectra = []
-        self.result = collections.OrderedDict()      
-
-    def run(self):
-        print 'Uvspectra.run'
-
-        # get observation list
-        observe = self.previous_results['observe']
-        obs_framework = observe['observed_framework']
-
-        # and calculate the result for each configuration
-        UVspectrum = collections.namedtuple('uvspectrum', [
-          'scan_number',
-          'baseline_x', 'baseline_y', 'baseline_z', 'baseline_number',
-          'interferogram', 'mirror_position', 'spectrum', 'wavenumber'],
-          verbose=False)
-
-        fts_position = []
-        fts_data = []
-
-        for config in obs_framework:
-            if config.fts_start:
-                if fts_data:
-                    # reduce previous scan - assume symmetric -tau to +tau, so first
-                    # shift 0 freq from centre of scan
-                    temp = np.fft.fftshift(fts_data)
-                    # spectrum is complex
-                    spectrum = np.fft.fft(temp)
-                    wavenumber = np.fft.fftfreq(n=spectrum.size,
-                      d=2.0 * abs(fts_position[1] - fts_position[0])) / 100.0
-                    
-                    # save it
-                    uvspectrum = UVspectrum(config.scan_number,
-                      baseline[0], 
-                      baseline[1],
-                      0.0,
-                      config.baseline_number,
-                      np.array(fts_data),
-                      np.array(fts_position),
-                      spectrum,
-                      wavenumber)
-
-                    self.uvspectra.append(uvspectrum)
-
-                # start new scan
-                fts_position = [config.fts_nominal_position]
-                fts_data = [config.data]
-
-            else:
-
-                # add points to current scan
-                fts_position.append(config.fts_nominal_position)
-                fts_data.append(config.data)
-            baseline = (config.baseline_x, config.baseline_y)
-
-        # deal with last scan
-        if fts_data:
-            temp = np.fft.fftshift(fts_data)
-            spectrum = np.fft.fft(temp)
-            wavenumber = np.fft.fftfreq(n=spectrum.size,
-              d=2.0 * abs(fts_position[1] - fts_position[0])) / 100.0
-                    
-            uvspectrum = UVspectrum(config.scan_number,
-              baseline[0], 
-              baseline[1],
-              0.0,
-              config.baseline_number,
-              np.array(fts_data),
-              np.array(fts_position),
-              spectrum,
-              wavenumber)
-
-            self.uvspectra.append(uvspectrum)
-
-        self.result['uvspectra'] = self.uvspectra
-        return self.result
-                
-    def __repr__(self):
-        return '''
-UVspectra:
-  #uvspectra : {num_uvspectra}
-'''.format(
-          num_uvspectra=len(self.uvspectra))
-
