Changeset 40

Timestamp:

Jun 10, 2014 10:41:57 AM (10 years ago)

Author:

JohnLightfoot

Message:

first working version

File:

• trunk/doublefourier.py (modified) (15 diffs)

trunk/doublefourier.py

@@ -39 +39 @@
         print 'DoubleFourier.run'
 
+        # convert lengths to m. Wavenumbers leave as cm-1.
         fts = self.previous_results['fts']
         fts_wn = fts['fts_wn']
-        opd_max = fts['opd_max']
+        fts_wn_truncated = fts['fts_wn_truncated']
+        opd_max = fts['opd_max'] / 100.0
         fts_nsample = fts['ftsnsample']
-        vdrive = fts['vdrive']
+        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 *= (opd_max / vdrive) / float(fts_nsample-1)
+        times *= delta_time
+        opd_start = -((fts_nsample + 1) / 2) * delta_opd
 
         beamsgenerator = self.previous_results['beamsgenerator']
@@ -57 +61 @@
         skymodel = skygenerator['sky model']
         spatial_axis = self.result['spatial axis'] = skygenerator['spatial axis']
-        self.result['frequency axis'] = fts_wn
+        self.result['frequency axis'] = fts_wn_truncated
 
         # assuming nx is even then transform has 0 freq at origin and [nx/2] is
@@ -70 +74 @@
         self.result['baseline interferograms'] = collections.OrderedDict()
 #        for baseline in bxby:
-        for baseline in bxby[:3]:
-            print baseline
+        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
@@ -78 +83 @@
 
                 # calculate the sky that the system is observing at this
-                # moment, incorporating various errors in turn
-                sky_now = skymodel
+                # 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.
@@ -108 +114 @@
                 # 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)
@@ -127 +132 @@
                 jobs = {}
                 # go through freq planes and shift them
-                for iwn,wn in enumerate(fts_wn):
+                for iwn,wn in enumerate(fts_wn_truncated):
                     # submit jobs
                     indata = (sky_now[:,:,iwn], shift,)
@@ -133 +138 @@
                       indata, (), ('numpy',))
                 
-                for iwn,wn in enumerate(fts_wn):
+                for iwn,wn in enumerate(fts_wn_truncated):
                     # collect and store results
                     temp = jobs[wn]()
@@ -139 +144 @@
 
                 if t == times[0]:
-                    self.result['sky at time 0'] = sky_now
+                    # 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
@@ -163 +169 @@
 
                 if t == times[0]:
-                    self.result['sky*beams at time 0'] = sky_now
+                    # 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
@@ -180 +187 @@
                 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):
+                for iwn,wn in enumerate(fts_wn_truncated):
 
                     # derive shift needed to place baseline at one of FFT coords
@@ -196 +203 @@
                     colpos = 0.0
                     rowpos = 0.0
-                    print 'spatial colpos, rowpos', colpos, rowpos
 
                     # calculate fourier phase shift to move point at [rowpos,colpos] to 
@@ -226 +232 @@
 
                     # set amp/phase at this frequency
-                    spectrum[iwn] = temp[0,0]
+                    spectrum[wn==fts_wn] = temp[0,0]
 
                 if t == times[0]:
-                    self.result['skyfft at time 0'] = fft_now
+                    self.result['skyfft at time 0'] = fft_now.copy()
 
                     axis = co.Axis(data=fts_wn, title='wavenumber',
@@ -242 +248 @@
 
                 if t == times[0]:
-                    # test version
+                    # test interferogram
                     # inverse fft of emission spectrum at this point
-                    temp = np.fft.ifft(spectrum)
-                    pos = np.fft.rfftfreq(len(spectrum))
+                    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
@@ -253 +264 @@
                     axis = co.Axis(data=pos, title='path difference',
                       units='cm')
-                    temp = co.Spectrum(data=temp,
+                    temp = co.Spectrum(data=temp, axis=axis,
                       title='Detected interferogram', units='')
 
@@ -259 +270 @@
 
                 mirror_error = 0.0
-                opd = 2.0 * (vdrive * t + mirror_error)
-                print 'opd', opd, vdrive, t, mirror_error
-                opd_ipos = opd / delta_opd
-
-                # make spectrum symmetric
+                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)
-                symmetric_spectrum = np.zeros([2*nfreq])
-                symmetric_spectrum[:nfreq] = spectrum
-                symmetric_spectrum[nfreq:] = spectrum[::-1]
-                
-                # calculate shift needed to move point at opd to 0
-                shift = np.zeros([2*nfreq], dtype=np.complex)
+                ndoublefreq = 2 * (nfreq - 1)
+                shift = np.zeros([ndoublefreq], dtype=np.complex)
                 shift[:nfreq] = np.arange(nfreq, dtype=np.complex)
-                shift[nfreq:] = np.arange(-nfreq, 0, dtype=np.complex)
-                shift = np.exp((-2.0j * np.pi * opd_ipos * shift) / float(nfreq))
+                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
-                symmetric_spectrum *= shift
-                spectrum_fft = np.fft.fft(symmetric_spectrum)
+                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)] = \
-              measurement
+              temp
 
         return self.result