1 | from __future__ import absolute_import
|
---|
2 |
|
---|
3 | import collections
|
---|
4 | import math
|
---|
5 | import numpy as np
|
---|
6 | import pp
|
---|
7 |
|
---|
8 | import common.commonobjects as co
|
---|
9 |
|
---|
10 | def fftshift(data, shift):
|
---|
11 | """Method to shift a 2d complex data array by applying the
|
---|
12 | given phase shift to its Fourier transform.
|
---|
13 | """
|
---|
14 | # move centre of image to array origin
|
---|
15 | temp = numpy.fft.fftshift(data)
|
---|
16 | # 2d fft
|
---|
17 | temp = numpy.fft.fft2(temp)
|
---|
18 | # apply phase shift
|
---|
19 | temp *= shift
|
---|
20 | # transform and shift back
|
---|
21 | temp = numpy.fft.ifft2(temp)
|
---|
22 | temp = numpy.fft.fftshift(temp)
|
---|
23 |
|
---|
24 | return temp
|
---|
25 |
|
---|
26 |
|
---|
27 | class DoubleFourier(object):
|
---|
28 | """Class to compute interferograms.
|
---|
29 | """
|
---|
30 |
|
---|
31 | def __init__(self, parameters, previous_results, job_server):
|
---|
32 | self.parameters = parameters
|
---|
33 | self.previous_results = previous_results
|
---|
34 | self.job_server = job_server
|
---|
35 |
|
---|
36 | self.result = collections.OrderedDict()
|
---|
37 |
|
---|
38 | def run(self):
|
---|
39 | print 'DoubleFourier.run'
|
---|
40 |
|
---|
41 | fts = self.previous_results['fts']
|
---|
42 | fts_wn = fts['fts_wn']
|
---|
43 | opd_max = fts['opd_max']
|
---|
44 | fts_nsample = fts['ftsnsample']
|
---|
45 | vdrive = fts['vdrive']
|
---|
46 | delta_opd = fts['delta_opd']
|
---|
47 |
|
---|
48 | times = np.arange(int(fts_nsample), dtype=np.float)
|
---|
49 | times *= (opd_max / vdrive) / float(fts_nsample-1)
|
---|
50 |
|
---|
51 | beamsgenerator = self.previous_results['beamsgenerator']
|
---|
52 |
|
---|
53 | uvmapgenerator = self.previous_results['uvmapgenerator']
|
---|
54 | bxby = uvmapgenerator['bxby']
|
---|
55 |
|
---|
56 | skygenerator = self.previous_results['skygenerator']
|
---|
57 | skymodel = skygenerator['sky model']
|
---|
58 | spatial_axis = self.result['spatial axis'] = skygenerator['spatial axis']
|
---|
59 | self.result['frequency axis'] = fts_wn
|
---|
60 |
|
---|
61 | # assuming nx is even then transform has 0 freq at origin and [nx/2] is
|
---|
62 | # Nyquist frequency. Nyq freq = 0.5 * Nyquist sampling freq.
|
---|
63 | # Assume further that the fft is shifted so that 0 freq is at nx/2
|
---|
64 | nx = len(spatial_axis)
|
---|
65 | spatial_freq_axis = np.arange(-nx/2, nx/2, dtype=np.float)
|
---|
66 | sample_freq = (180.0 * 3600.0 / np.pi) / (spatial_axis[1] - spatial_axis[0])
|
---|
67 | spatial_freq_axis *= (sample_freq / nx)
|
---|
68 | self.result['spatial frequency axis'] = spatial_freq_axis
|
---|
69 |
|
---|
70 | self.result['baseline interferograms'] = collections.OrderedDict()
|
---|
71 | # for baseline in bxby:
|
---|
72 | for baseline in bxby[:3]:
|
---|
73 | print baseline
|
---|
74 | measurement = np.zeros(np.shape(times))
|
---|
75 |
|
---|
76 | # FTS path diff and possibly baseline itself vary with time
|
---|
77 | for tindex,t in enumerate(times):
|
---|
78 |
|
---|
79 | # calculate the sky that the system is observing at this
|
---|
80 | # moment, incorporating various errors in turn
|
---|
81 | sky_now = skymodel
|
---|
82 |
|
---|
83 | # 1. baseline should be perp to centre of field.
|
---|
84 | # If baseline is tilted then origin of sky map shifts.
|
---|
85 | # (I think effect could be corrected by changing
|
---|
86 | # FTS sample position to compensate.?)
|
---|
87 | # for now assume 0 error but do full calculation for timing
|
---|
88 | # purposes
|
---|
89 |
|
---|
90 | # perfect baseline is perpendicular to direction to 'centre'
|
---|
91 | # on sky. Add errors in x,y,z - z towards sky 'centre'
|
---|
92 | bz = 0.0
|
---|
93 | # convert bz to angular error
|
---|
94 | blength = math.sqrt(baseline[0]*baseline[0] +
|
---|
95 | baseline[1]*baseline[1])
|
---|
96 | bangle = (180.0 * 3600.0 / math.pi) * bz / blength
|
---|
97 | bx_error = bangle * baseline[0] / blength
|
---|
98 | by_error = bangle * baseline[1] / blength
|
---|
99 |
|
---|
100 | # calculate xpos, ypos in units of pixel - numpy arrays
|
---|
101 | # [row,col]
|
---|
102 | nx = len(spatial_axis)
|
---|
103 | colpos = float(nx-1) * bx_error / \
|
---|
104 | (spatial_axis[-1] - spatial_axis[0])
|
---|
105 | rowpos = float(nx-1) * by_error / \
|
---|
106 | (spatial_axis[-1] - spatial_axis[0])
|
---|
107 |
|
---|
108 | # calculate fourier phase shift to move point at [0,0] to
|
---|
109 | # [rowpos, colpos]
|
---|
110 | print 'shift', colpos, rowpos
|
---|
111 | shiftx = np.zeros([nx], np.complex)
|
---|
112 | shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
113 | shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
114 | shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / float(nx))
|
---|
115 |
|
---|
116 | shifty = np.zeros([nx], np.complex)
|
---|
117 | shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
118 | shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
119 | shifty = np.exp((-2.0j * np.pi * rowpos * shifty) / float(nx))
|
---|
120 |
|
---|
121 | shift = np.ones([nx,nx], np.complex)
|
---|
122 | for j in range(nx):
|
---|
123 | shift[j,:] *= shiftx
|
---|
124 | for i in range(nx):
|
---|
125 | shift[:,i] *= shifty
|
---|
126 |
|
---|
127 | jobs = {}
|
---|
128 | # go through freq planes and shift them
|
---|
129 | for iwn,wn in enumerate(fts_wn):
|
---|
130 | # submit jobs
|
---|
131 | indata = (sky_now[:,:,iwn], shift,)
|
---|
132 | jobs[wn] = self.job_server.submit(fftshift,
|
---|
133 | indata, (), ('numpy',))
|
---|
134 |
|
---|
135 | for iwn,wn in enumerate(fts_wn):
|
---|
136 | # collect and store results
|
---|
137 | temp = jobs[wn]()
|
---|
138 | sky_now[:,:,iwn] = temp
|
---|
139 |
|
---|
140 | if t == times[0]:
|
---|
141 | self.result['sky at time 0'] = sky_now
|
---|
142 |
|
---|
143 | # 2. telescopes should be centred on centre of field
|
---|
144 | # Telescopes collect flux from the 'sky' and pass
|
---|
145 | # it to the FTS beam combiner. In doing this each
|
---|
146 | # telescope multiplies the sky emission by its
|
---|
147 | # amplitude beam response - always real but with
|
---|
148 | # negative areas. Is this correct? Gives right
|
---|
149 | # answer for 'no error' case.
|
---|
150 |
|
---|
151 | # multiply sky by amplitude beam 1 * amplitude beam 2
|
---|
152 | amp_beam_1 = beamsgenerator['primary amplitude beam'].data
|
---|
153 | amp_beam_2 = beamsgenerator['primary amplitude beam'].data
|
---|
154 |
|
---|
155 | for iwn,wn in enumerate(fts_wn):
|
---|
156 | # calculate shifted beams here
|
---|
157 | # for now assume no errors and just use beam
|
---|
158 | # calculated earlier
|
---|
159 | pass
|
---|
160 |
|
---|
161 | # multiply sky by amplitude beams of 2 antennas
|
---|
162 | sky_now *= amp_beam_1 * amp_beam_2
|
---|
163 |
|
---|
164 | if t == times[0]:
|
---|
165 | self.result['sky*beams at time 0'] = sky_now
|
---|
166 |
|
---|
167 | # 3. baseline error revisited
|
---|
168 | # derive baseline at this time
|
---|
169 | #
|
---|
170 | # Perhaps baseline should be a continuous function in
|
---|
171 | # time, which would allow baselines that intentionally
|
---|
172 | # smoothly vary (as in rotating tethered assembly)
|
---|
173 | # and errors to be handled by one description.
|
---|
174 | #
|
---|
175 | # what follows assumes zero error
|
---|
176 |
|
---|
177 | baseline_error = 0.0
|
---|
178 | baseline_now = baseline + baseline_error
|
---|
179 |
|
---|
180 | fft_now = np.zeros(np.shape(sky_now), np.complex)
|
---|
181 | spectrum = np.zeros(np.shape(fts_wn), np.complex)
|
---|
182 | for iwn,wn in enumerate(fts_wn):
|
---|
183 |
|
---|
184 | # derive shift needed to place baseline at one of FFT coords
|
---|
185 | # this depends on physical baseline and frequency
|
---|
186 | baseline_now_lambdas = baseline_now * wn * 100.0
|
---|
187 |
|
---|
188 | # calculate baseline position in units of pixels of FFTed
|
---|
189 | # sky - numpy arrays [row,col]
|
---|
190 | colpos = float(nx-1) * \
|
---|
191 | float(baseline_now_lambdas[0] - spatial_freq_axis[0]) / \
|
---|
192 | (spatial_freq_axis[-1] - spatial_freq_axis[0])
|
---|
193 | rowpos = float(nx-1) * \
|
---|
194 | float(baseline_now_lambdas[1] - spatial_freq_axis[0]) / \
|
---|
195 | (spatial_freq_axis[-1] - spatial_freq_axis[0])
|
---|
196 | colpos = 0.0
|
---|
197 | rowpos = 0.0
|
---|
198 | print 'spatial colpos, rowpos', colpos, rowpos
|
---|
199 |
|
---|
200 | # calculate fourier phase shift to move point at [rowpos,colpos] to
|
---|
201 | # [0,0]
|
---|
202 | shiftx = np.zeros([nx], np.complex)
|
---|
203 | shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
204 | shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
205 | shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / \
|
---|
206 | float(nx))
|
---|
207 |
|
---|
208 | shifty = np.zeros([nx], np.complex)
|
---|
209 | shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
210 | shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
211 | shifty = np.exp((-2.0j * np.pi * rowpos * shifty) / float(nx))
|
---|
212 |
|
---|
213 | shift = np.ones([nx,nx], np.complex)
|
---|
214 | for j in range(nx):
|
---|
215 | shift[j,:] *= shiftx
|
---|
216 | for i in range(nx):
|
---|
217 | shift[:,i] *= shifty
|
---|
218 |
|
---|
219 | # move centre of sky image to origin
|
---|
220 | temp = np.fft.fftshift(sky_now[:,:,iwn])
|
---|
221 | # apply phase shift
|
---|
222 | temp *= shift
|
---|
223 | # 2d fft
|
---|
224 | temp = np.fft.fft2(temp)
|
---|
225 | fft_now[:,:,iwn] = temp
|
---|
226 |
|
---|
227 | # set amp/phase at this frequency
|
---|
228 | spectrum[iwn] = temp[0,0]
|
---|
229 |
|
---|
230 | if t == times[0]:
|
---|
231 | self.result['skyfft at time 0'] = fft_now
|
---|
232 |
|
---|
233 | axis = co.Axis(data=fts_wn, title='wavenumber',
|
---|
234 | units='cm-1')
|
---|
235 | temp = co.Spectrum(data=spectrum, axis=axis,
|
---|
236 | title='Detected spectrum', units='W sr-1 m-2 Hz-1')
|
---|
237 | self.result['skyfft spectrum at time 0'] = temp
|
---|
238 |
|
---|
239 | # 3. FTS sampling should be accurate
|
---|
240 | # derive lag due to FTS path difference
|
---|
241 | # 0 error for now
|
---|
242 |
|
---|
243 | if t == times[0]:
|
---|
244 | # test version
|
---|
245 | # inverse fft of emission spectrum at this point
|
---|
246 | temp = np.fft.ifft(spectrum)
|
---|
247 | pos = np.fft.rfftfreq(len(spectrum))
|
---|
248 |
|
---|
249 | # move 0 frequency to centre of array
|
---|
250 | temp = np.fft.fftshift(temp)
|
---|
251 | pos = np.fft.fftshift(pos)
|
---|
252 |
|
---|
253 | axis = co.Axis(data=pos, title='path difference',
|
---|
254 | units='cm')
|
---|
255 | temp = co.Spectrum(data=temp,
|
---|
256 | title='Detected interferogram', units='')
|
---|
257 |
|
---|
258 | self.result['test FTS at time 0'] = temp
|
---|
259 |
|
---|
260 | mirror_error = 0.0
|
---|
261 | opd = 2.0 * (vdrive * t + mirror_error)
|
---|
262 | print 'opd', opd, vdrive, t, mirror_error
|
---|
263 | opd_ipos = opd / delta_opd
|
---|
264 |
|
---|
265 | # make spectrum symmetric
|
---|
266 | nfreq = len(fts_wn)
|
---|
267 | symmetric_spectrum = np.zeros([2*nfreq])
|
---|
268 | symmetric_spectrum[:nfreq] = spectrum
|
---|
269 | symmetric_spectrum[nfreq:] = spectrum[::-1]
|
---|
270 |
|
---|
271 | # calculate shift needed to move point at opd to 0
|
---|
272 | shift = np.zeros([2*nfreq], dtype=np.complex)
|
---|
273 | shift[:nfreq] = np.arange(nfreq, dtype=np.complex)
|
---|
274 | shift[nfreq:] = np.arange(-nfreq, 0, dtype=np.complex)
|
---|
275 | shift = np.exp((-2.0j * np.pi * opd_ipos * shift) / float(nfreq))
|
---|
276 |
|
---|
277 | # apply phase shift and fft
|
---|
278 | symmetric_spectrum *= shift
|
---|
279 | spectrum_fft = np.fft.fft(symmetric_spectrum)
|
---|
280 | measurement[tindex] = spectrum_fft[0]
|
---|
281 |
|
---|
282 | self.result['baseline interferograms'][tuple(baseline)] = \
|
---|
283 | measurement
|
---|
284 |
|
---|
285 | return self.result
|
---|
286 |
|
---|
287 | def matlab_transform(self):
|
---|
288 | # readers should look at Izumi et al. 2006, Applied Optics, 45, 2576
|
---|
289 | # for theoretical background. Names of variables in the code correspond
|
---|
290 | # to that work.
|
---|
291 |
|
---|
292 | # For now, assume 2 light collectors giving one baseline at a time.
|
---|
293 | interferograms = {}
|
---|
294 |
|
---|
295 | for baseline in self.baselines:
|
---|
296 | interferogram = 0
|
---|
297 |
|
---|
298 | # baseline length (cm) and position angle
|
---|
299 | bu = baseline[0]
|
---|
300 | bv = baseline[1]
|
---|
301 | mod_b = np.sqrt(pow(bu,2) + pow(bv,2))
|
---|
302 | ang_b = np.arctan2(bv, bu)
|
---|
303 |
|
---|
304 | # loop over sky pixels covered by primary beam
|
---|
305 | nx = self.sky_s.shape[1]
|
---|
306 | ny = self.sky_s.shape[2]
|
---|
307 | for j in range(ny):
|
---|
308 | for i in range(nx):
|
---|
309 |
|
---|
310 | # inverse fft of emission spectrum at this point
|
---|
311 | temp = np.fft.ifft(self.sky_s[:,j,i])
|
---|
312 |
|
---|
313 | # move 0 frequency to centre of array
|
---|
314 | temp = np.fft.fftshift(temp)
|
---|
315 |
|
---|
316 | # length (radians) and position angle of theta vector
|
---|
317 | mod_theta = np.sqrt(pow(self.sky_x[i],2) + pow(self.sky_y[j],2))
|
---|
318 | ang_theta = np.arctan2(self.sky_y[j], self.sky_x[i])
|
---|
319 |
|
---|
320 | # calculate b.theta (the projection of b on theta)
|
---|
321 | # and the corresponding delay in units of wavelength at
|
---|
322 | # Nyquist frequency
|
---|
323 | delay = mod_theta * mod_b * np.cos(ang_b - ang_theta) * self.freqs[-1]
|
---|
324 |
|
---|
325 | # sampling is done at twice Nyquist freq so shift transformed
|
---|
326 | # spectrum by 2 * delay samples (approximated to nint)
|
---|
327 | # NOTE factor of 2 discrepency with matlab version! I think
|
---|
328 | # this is because there the variable 'Nyq' is the Nyquist
|
---|
329 | # sampling rate, not the Nyquist frequency.
|
---|
330 | temp = np.roll(temp, int(round(2.0 * delay)))
|
---|
331 |
|
---|
332 | # want only the real part of the result
|
---|
333 | interferogram += np.real(temp)
|
---|
334 |
|
---|
335 | def __repr__(self):
|
---|
336 | return 'DoubleFourier'
|
---|
337 |
|
---|