Index: /trunk/doublefourier.py
===================================================================
--- /trunk/doublefourier.py	(revision 30)
+++ /trunk/doublefourier.py	(revision 31)
@@ -2,44 +2,25 @@
 
 import collections
+import math
 import numpy as np
 import pp
 
-def test(nx=50, ny=50):
-    # NOT WORKING
-    # all wavenumbers are cm-1
-    wn_max = 100.0
-    wn_min = 50.0
-
-    resolution_min = 300.0
-
-    delta_wn = wn_min / resolution_min
-
-    # spectra are evenly sampled from 0 to wn_max with spacing delta_wn
-    nfreq = np.ceil(wn_max / delta_wn)
-    wn_max = delta_wn * nfreq
-
-    freqs = np.arange(nfreq)
-    freqs *= delta_wn
-
-    # construct a dummy 'sky'
-    # numpy default indexing is C-style so that the rightmost index
-    # 'changes the fastest'.
-    sky_s = np.zeros([nfreq, ny, nx], np.float)
-    sky_s[:,ny/2,nx/2] = 1.0
-
-    # x coords in radians
-    arcsec2rad = 1.0 / 206265
-    sky_x = (np.arange(nx) - nx/2) * 1.0 * arcsec2rad 
-    sky_y = (np.arange(ny) - ny/2) * 1.0 * arcsec2rad 
-
-    # dummy baselines (each baseline entry is [u,v] in centimetres)
-    baselines = []
-    for i in range(500):
-        baselines.append([20000.0, 20000.0])
-
-    df = DoubleFourier(sky_s=sky_s, freqs=freqs, sky_x=sky_x, sky_y=sky_y, 
-      baselines=baselines)
-
-    df.matlab_transform()
+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
 
 
@@ -48,13 +29,8 @@
     """
 
-    def __init__(self, parameters, previous_results):
+    def __init__(self, parameters, previous_results, job_server):
         self.parameters = parameters
         self.previous_results = previous_results
-
-        # paralle processing stuff
-        ppservers = ()
-        self.job_server = pp.Server(ppservers=ppservers)
-        print 'DoubleFourier starting pp with %s workers' % \
-          self.job_server.get_ncpus()
+        self.job_server = job_server
 
         self.result = collections.OrderedDict()      
@@ -64,5 +40,5 @@
 
         fts = self.previous_results['fts']
-        frequency_axis = fts['fts_wn']
+        fts_wn = fts['fts_wn']
         opd_max = fts['opd_max']
         fts_nsample = fts['ftsnsample']
@@ -81,5 +57,5 @@
         skymodel = skygenerator['sky model']
         spatial_axis = self.result['spatial axis'] = skygenerator['spatial axis']
-        self.result['frequency axis'] = frequency_axis
+        self.result['frequency axis'] = fts_wn
 
         # assuming nx is even then transform has 0 freq at origin and [nx/2] is
@@ -100,10 +76,8 @@
             # FTS path diff and possibly baseline itself vary with time
             for tindex,t in enumerate(times):
-#            for tindex,t in enumerate(times[:1]):
-
+
+                # calculate the sky that the system is observing at this
+                # moment, incorporating various errors in turn
                 sky_now = skymodel
-
-                # what is the system looking at?
-                # add 'errors' in order
 
                 # 1. baseline should be perp to centre of field.
@@ -113,22 +87,26 @@
                 # for now assume 0 error but do full calculation for timing
                 # purposes
-                baseline_dx = 0.0
-                baseline_dy = 0.0
+
+                # 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) * float (baseline_dx - spatial_axis[0]) / \
+                colpos = float(nx-1) * bx_error / \
                   (spatial_axis[-1] - spatial_axis[0])
-                rowpos = float(nx-1) * float (baseline_dy - spatial_axis[0]) / \
+                rowpos = float(nx-1) * by_error / \
                   (spatial_axis[-1] - spatial_axis[0])
-
-                if colpos < 0 or colpos > (nx-1) or rowpos < 0 or rowpos > \
-                  (nx-1):
-                    raise Exception, \
-                      'baseline centre outside limits of sky model'
 
                 # calculate fourier phase shift to move point at [0,0] to 
                 # [rowpos, colpos]
+                print 'shift', colpos, rowpos
                 shiftx = np.zeros([nx], np.complex)
                 shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
@@ -147,24 +125,20 @@
                     shift[:,i] *= shifty
 
+                jobs = {}
                 # go through freq planes and shift them
-                for iwn,wn in enumerate(frequency_axis):
-                    # move centre of image to array origin
-                    temp = np.fft.fftshift(sky_now[:,:,iwn])
-                    # 2d fft
-                    temp = np.fft.fft2(temp)
-                    # apply phase shift
-                    temp *= shift
-                    # transform and shift back
-                    temp = np.fft.ifft2(temp)
-                    temp = np.fft.fftshift(temp)
-#                    imag part should still be 0
-#                    print 'max real', np.max(np.real(temp))
-#                    print 'max imag', np.max(np.imag(temp))
-                    temp = np.abs(temp)
-                    sky_now[:,:,iwn] = np.fft.fftshift(temp)
+                for iwn,wn in enumerate(fts_wn):
+                    # submit jobs
+                    indata = (sky_now[:,:,iwn], shift,)
+                    jobs[wn] = self.job_server.submit(fftshift,
+                      indata, (), ('numpy',))
+                
+                for iwn,wn in enumerate(fts_wn):
+                    # collect and store results
+                    temp = jobs[wn]()
+                    sky_now[:,:,iwn] = temp
 
                 if t == times[0]:
                     self.result['sky at time 0'] = sky_now
-                
+
                 # 2. telescopes should be centred on centre of field
                 #    Telescopes collect flux from the 'sky' and pass
@@ -176,14 +150,15 @@
 
                 # multiply sky by amplitude beam 1 * amplitude beam 2
-                for iwn,wn in enumerate(frequency_axis):
+                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
-
-                    amplitude_beam = beamsgenerator['primary amplitude beam']\
-                      [wn]                
-                    amp_beam_1 = amplitude_beam.data
-                    amp_beam_2 = amplitude_beam.data
-                    sky_now[:,:,iwn] *= amp_beam_1 * amp_beam_2 
+                    pass
+
+                # multiply sky by amplitude beams of 2 antennas
+                sky_now *= amp_beam_1 * amp_beam_2 
 
                 if t == times[0]:
@@ -192,25 +167,34 @@
                 # 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
-#                print 'baseline_now', baseline_now
-
-#                gamma_12[(baseline,t)] = 0.0
+
                 fft_now = np.zeros(np.shape(sky_now), np.complex)
-                spectrum = np.zeros(np.shape(frequency_axis), np.complex)
-                for iwn,wn in enumerate(frequency_axis):
+                spectrum = np.zeros(np.shape(fts_wn), np.complex)
+                for iwn,wn in enumerate(fts_wn):
 
                     # 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 
-#                    print 'baseline lambda', wn, baseline_now_lambdas
-
-                    # 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]) / \
+
+                    # 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]) / \
+                    rowpos = float(nx-1) * \
+                      float(baseline_now_lambdas[1] - spatial_freq_axis[0]) / \
                       (spatial_freq_axis[-1] - spatial_freq_axis[0])
-#                    print 'spatial colpos, rowpos', colpos, rowpos
+                    colpos = 0.0
+                    rowpos = 0.0
+                    print 'spatial colpos, rowpos', colpos, rowpos
 
                     # calculate fourier phase shift to move point at [rowpos,colpos] to 
@@ -219,5 +203,6 @@
                     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))
+                    shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / \
+                      float(nx))
 
                     shifty = np.zeros([nx], np.complex)
@@ -238,15 +223,39 @@
                     # 2d fft
                     temp = np.fft.fft2(temp)
-                    # shift 0 freq to centre of array
-                    fft_now[:,:,iwn] = np.fft.ifftshift(temp)
-
+                    fft_now[:,:,iwn] = temp
+
+                    # set amp/phase at this frequency
                     spectrum[iwn] = temp[0,0]
 
                 if t == times[0]:
                     self.result['skyfft at time 0'] = fft_now
+
+                    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 version
+                    # inverse fft of emission spectrum at this point
+                    temp = np.fft.ifft(spectrum)
+                    pos = np.fft.rfftfreq(len(spectrum))
+
+                    # 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,
+                      title='Detected interferogram', units='')
+
+                    self.result['test FTS at time 0'] = temp
+
                 mirror_error = 0.0
                 opd = 2.0 * (vdrive * t + mirror_error)
@@ -255,5 +264,5 @@
 
                 # make spectrum symmetric
-                nfreq = len(frequency_axis)
+                nfreq = len(fts_wn)
                 symmetric_spectrum = np.zeros([2*nfreq])
                 symmetric_spectrum[:nfreq] = spectrum
@@ -271,5 +280,6 @@
                 measurement[tindex] = spectrum_fft[0]
 
-            self.result['baseline interferograms'][tuple(baseline)] = measurement
+            self.result['baseline interferograms'][tuple(baseline)] = \
+              measurement
 
         return self.result
