1 | from __future__ import absolute_import
|
---|
2 |
|
---|
3 | import collections
|
---|
4 | import math
|
---|
5 | import numpy as np
|
---|
6 | import pp
|
---|
7 |
|
---|
8 | import matplotlib.pyplot as plt
|
---|
9 |
|
---|
10 | import common.commonobjects as co
|
---|
11 |
|
---|
12 | def fftshift(data, shift):
|
---|
13 | """Method to shift a 2d complex data array by applying the
|
---|
14 | given phase shift to its Fourier transform.
|
---|
15 | """
|
---|
16 | # move centre of image to array origin
|
---|
17 | temp = numpy.fft.fftshift(data)
|
---|
18 | # 2d fft
|
---|
19 | temp = numpy.fft.fft2(temp)
|
---|
20 | # apply phase shift
|
---|
21 | temp *= shift
|
---|
22 | # transform and shift back
|
---|
23 | temp = numpy.fft.ifft2(temp)
|
---|
24 | temp = numpy.fft.fftshift(temp)
|
---|
25 |
|
---|
26 | return temp
|
---|
27 |
|
---|
28 | def fields_match(previous_config, config, exclude_fields=[]):
|
---|
29 | """Function to compare 2 configs, ignoring any differences in
|
---|
30 | the components listed in 'exclude_fields'.
|
---|
31 | """
|
---|
32 | if previous_config is None:
|
---|
33 | return False
|
---|
34 |
|
---|
35 | if not exclude_fields:
|
---|
36 | return previous_config==config
|
---|
37 | else:
|
---|
38 | # convert to dicts and eliminate exclude_fields
|
---|
39 | previous_dict = previous_config._asdict()
|
---|
40 | current_dict = config._asdict()
|
---|
41 | for field in exclude_fields:
|
---|
42 | ignore = previous_dict.pop(field, None)
|
---|
43 | ignore = current_dict.pop(field, None)
|
---|
44 |
|
---|
45 | previous_set = set([k for k in previous_dict.items()])
|
---|
46 | current_set = set([k for k in current_dict.items()])
|
---|
47 | diff = previous_set.symmetric_difference(current_set)
|
---|
48 |
|
---|
49 | debug = False
|
---|
50 | if(debug):
|
---|
51 | print diff
|
---|
52 |
|
---|
53 | return not(bool(diff))
|
---|
54 |
|
---|
55 |
|
---|
56 | class Observe(object):
|
---|
57 | """Class to compute interferograms.
|
---|
58 | """
|
---|
59 |
|
---|
60 | def __init__(self, parameters, previous_results, job_server):
|
---|
61 | self.parameters = parameters
|
---|
62 | self.previous_results = previous_results
|
---|
63 | self.job_server = job_server
|
---|
64 |
|
---|
65 | self.result = collections.OrderedDict()
|
---|
66 |
|
---|
67 | def run(self):
|
---|
68 | print 'Observe.run'
|
---|
69 |
|
---|
70 | # access primary beam information
|
---|
71 | beamsgenerator = self.previous_results['beamsgenerator']
|
---|
72 |
|
---|
73 | # access required FTS information
|
---|
74 | fts = self.previous_results['fts']
|
---|
75 | fts_wn = fts['fts_wn']
|
---|
76 | fts_wn_truncated = fts['fts_wn_truncated']
|
---|
77 | delta_opd = fts['delta_opd']
|
---|
78 |
|
---|
79 | # access to model sky
|
---|
80 | skygenerator = self.previous_results['skygenerator']
|
---|
81 | sky_model = skygenerator['sky model']
|
---|
82 | spatial_axis = self.result['spatial axis'] = skygenerator['spatial axis']
|
---|
83 | nx = len(spatial_axis)
|
---|
84 |
|
---|
85 | # assuming nx is even then transform of spatial axis has 0 freq at origin
|
---|
86 | # and [nx/2] is Nyquist frequency [Nyq freq = 0.5 * Nyquist sampling freq].
|
---|
87 | # The fft is shifted so that 0 freq is at nx/2.
|
---|
88 | nx = len(spatial_axis)
|
---|
89 | spatial_freq_axis = np.arange(-nx/2, nx/2, dtype=np.float)
|
---|
90 | sample_freq = (180.0 * 3600.0 / np.pi) / (spatial_axis[1] - spatial_axis[0])
|
---|
91 | spatial_freq_axis *= (sample_freq / nx)
|
---|
92 | self.result['spatial frequency axis'] = spatial_freq_axis
|
---|
93 |
|
---|
94 | # get list of instrument configuratons
|
---|
95 | uvmapgenerator = self.previous_results['uvmapgenerator']
|
---|
96 | obs_framework = uvmapgenerator['obs_framework']
|
---|
97 |
|
---|
98 | # and calculate the result for each configuration
|
---|
99 | previous_config = None
|
---|
100 | observed_obs_framework = []
|
---|
101 |
|
---|
102 | debug_plotted = False
|
---|
103 | for config in obs_framework:
|
---|
104 | # if config.fts_start:
|
---|
105 | # print 'start FTS scan'
|
---|
106 | # print 'baseline', config.baseline_x, config.baseline_y,\
|
---|
107 | # config.baseline_z
|
---|
108 |
|
---|
109 | if fields_match(previous_config, config,
|
---|
110 | exclude_fields=['scan_number',
|
---|
111 | 'fts_start',
|
---|
112 | 'fts_position',
|
---|
113 | 'fts_nominal_position',
|
---|
114 | 'baseline_x',
|
---|
115 | 'baseline_y',
|
---|
116 | 'baseline_z',
|
---|
117 | 'baseline_number',
|
---|
118 | 'data']):
|
---|
119 | pass
|
---|
120 | else:
|
---|
121 | # print 'previous', previous_config
|
---|
122 | # print 'current', config
|
---|
123 | # calculate new sky
|
---|
124 | # print 'calculating new sky'
|
---|
125 |
|
---|
126 | # Be explicit about copies otherwise a reference is
|
---|
127 | # made which corrupts the original and leads to much
|
---|
128 | # confusion.
|
---|
129 | sky_now = sky_model.copy()
|
---|
130 |
|
---|
131 | # calculate the sky that the system is observing at this
|
---|
132 | # moment, incorporating various errors in turn.
|
---|
133 |
|
---|
134 | # 1. baseline should be perp to vector towards centre
|
---|
135 | # of field. If baseline is tilted then origin of sky
|
---|
136 | # map shifts.
|
---|
137 | # I think this error can be folded into that for the
|
---|
138 | # FTS mirror, so leave handling it until then. Is
|
---|
139 | # this correct?
|
---|
140 |
|
---|
141 | # 2. Flux-collecting telescopes should be pointing
|
---|
142 | # at centre of field. They each collect flux from
|
---|
143 | # the 'sky' and pass it to the FTS beam combiner.
|
---|
144 | # In doing this each telescope multiplies the sky
|
---|
145 | # emission by its complex amplitude beam response.
|
---|
146 | # So instead of measuring the correlation of the
|
---|
147 | # 'sky' with itself, we measure that x Beam1.Beam2*.
|
---|
148 | # Is this correct? Gives right answer for 'no error'
|
---|
149 | # where Beam1.Beam2* == PSF.
|
---|
150 |
|
---|
151 | # multiply sky by amp.beam 1 * conj(amp.beam 2)
|
---|
152 | amp_beam_1 = beamsgenerator['primary amplitude beam'].data
|
---|
153 | amp_beam_2 = beamsgenerator['primary amplitude beam'].data
|
---|
154 |
|
---|
155 | # shift beams if pointing error
|
---|
156 | if config.point1_x_error or config.point1_y_error:
|
---|
157 | # calculate shifted beam1
|
---|
158 | # raise Exception, 'not implemented'
|
---|
159 | current_amp_beam_1 = amp_beam_1
|
---|
160 | else:
|
---|
161 | current_amp_beam_1 = amp_beam_1
|
---|
162 |
|
---|
163 | if config.point2_x_error or config.point2_y_error:
|
---|
164 | # calculate shifted beam2
|
---|
165 | # raise Exception, 'not implemented'
|
---|
166 | current_amp_beam_2 = amp_beam_2
|
---|
167 | else:
|
---|
168 | current_amp_beam_2 = amp_beam_2
|
---|
169 |
|
---|
170 | sky_now *= current_amp_beam_1 * \
|
---|
171 | np.conj(current_amp_beam_2)
|
---|
172 |
|
---|
173 | # 3. Get result for baseline at this time
|
---|
174 | if fields_match(previous_config, config,
|
---|
175 | exclude_fields=['fts_start',
|
---|
176 | 'fts_position',
|
---|
177 | 'fts_nominal_position',
|
---|
178 | 'data']):
|
---|
179 | pass
|
---|
180 | else:
|
---|
181 | # calculate new baseline spectrum
|
---|
182 | # print 'calculating new baseline spectrum'
|
---|
183 |
|
---|
184 | baseline = np.array([config.baseline_x, config.baseline_y])
|
---|
185 | fft_now = np.zeros(np.shape(sky_now), np.complex)
|
---|
186 | baseline_spectrum = np.zeros(np.shape(fts_wn), np.complex)
|
---|
187 |
|
---|
188 | for iwn,wn in enumerate(fts_wn_truncated):
|
---|
189 |
|
---|
190 | # derive shift needed to place baseline at one of FFT
|
---|
191 | # coords this depends on physical baseline and frequency
|
---|
192 | baseline_lambda = baseline * wn * 100.0
|
---|
193 |
|
---|
194 | # calculate baseline position in units of pixels of
|
---|
195 | # FFTed sky - numpy arrays [row,col]
|
---|
196 | colpos = float(nx-1) * baseline_lambda[0] / \
|
---|
197 | (spatial_freq_axis[-1] - spatial_freq_axis[0])
|
---|
198 | rowpos = float(nx-1) * baseline_lambda[1] / \
|
---|
199 | (spatial_freq_axis[-1] - spatial_freq_axis[0])
|
---|
200 |
|
---|
201 | # calculate fourier phase shift to move point at
|
---|
202 | # [rowpos,colpos] to [0,0]
|
---|
203 | shiftx = np.zeros([nx], np.complex)
|
---|
204 | shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
205 | shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
206 | shiftx = np.exp((2.0j * np.pi * colpos * shiftx) / \
|
---|
207 | float(nx))
|
---|
208 |
|
---|
209 | shifty = np.zeros([nx], np.complex)
|
---|
210 | shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
211 | shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
212 | shifty = np.exp((2.0j * np.pi * rowpos * shifty) / \
|
---|
213 | float(nx))
|
---|
214 |
|
---|
215 | shift = np.ones([nx,nx], np.complex)
|
---|
216 | for j in range(nx):
|
---|
217 | shift[j,:] *= shiftx
|
---|
218 | for i in range(nx):
|
---|
219 | shift[:,i] *= shifty
|
---|
220 |
|
---|
221 | # move centre of sky image to origin
|
---|
222 | temp = np.fft.fftshift(sky_now[:,:,iwn])
|
---|
223 | # apply phase shift
|
---|
224 | temp *= shift
|
---|
225 | # 2d fft
|
---|
226 | temp = np.fft.fft2(temp)
|
---|
227 | fft_now[:,:,iwn] = temp
|
---|
228 |
|
---|
229 | if wn==30.0 and not debug_plotted:
|
---|
230 |
|
---|
231 | freq = np.fft.fftfreq(n=np.shape(fft_now)[0],
|
---|
232 | d=(spatial_axis[1]-spatial_axis[0]) / 206265.0)
|
---|
233 | freq = np.arange(np.shape(fft_now)[0]) * (freq[1] - freq[0])
|
---|
234 |
|
---|
235 | ifreq = np.arange(np.shape(fft_now)[0])
|
---|
236 | baseline_i = ifreq[abs(freq-abs(baseline_lambda[0])) < (freq[1]-freq[0])]
|
---|
237 | baseline_j = ifreq[abs(freq-abs(baseline_lambda[1])) < (freq[1]-freq[0])]
|
---|
238 | # print baseline, fft_now[baseline_i[0], baseline_j[0], iwn]
|
---|
239 |
|
---|
240 | # plt.figure()
|
---|
241 | # plt.imshow(fft_now.real[:,:,iwn], interpolation='nearest', origin='lower',
|
---|
242 | # aspect='equal', extent=[freq[0],freq[-1],freq[0],freq[-1]])
|
---|
243 | # plt.colorbar(orientation='vertical')
|
---|
244 | # plt.axis('image')
|
---|
245 | # plt.savefig('debug_fft_real.png')
|
---|
246 | # plt.close()
|
---|
247 |
|
---|
248 | # plt.figure()
|
---|
249 | # plt.imshow(fft_now.imag[:,:,iwn], interpolation='nearest', origin='lower',
|
---|
250 | # aspect='equal', extent=[freq[0],freq[-1],freq[0],freq[-1]])
|
---|
251 | # plt.colorbar(orientation='vertical')
|
---|
252 | # plt.axis('image')
|
---|
253 | # plt.savefig('debug_fft_imag.png')
|
---|
254 | # plt.close()
|
---|
255 |
|
---|
256 | # debug_plotted = True
|
---|
257 |
|
---|
258 | # set amp/phase at this frequency
|
---|
259 | baseline_spectrum[wn==fts_wn] = temp[0,0]
|
---|
260 |
|
---|
261 | # 4. Calculate interferogram value from this baseline on the 'sky'
|
---|
262 | # at the current FTS mirror position.
|
---|
263 |
|
---|
264 | opd = 2.0 * config.fts_position
|
---|
265 | opd_ipos = opd / delta_opd
|
---|
266 |
|
---|
267 | # calculate shift needed to move point at opd to 0
|
---|
268 | nfreq = len(fts_wn)
|
---|
269 | ndoublefreq = 2 * (nfreq - 1)
|
---|
270 | shift = np.zeros([ndoublefreq], dtype=np.complex)
|
---|
271 | shift[:nfreq] = np.arange(nfreq, dtype=np.complex)
|
---|
272 | shift[nfreq:] = np.arange(nfreq, dtype=np.complex)[-2:0:-1]
|
---|
273 | shift = np.exp((-2.0j * np.pi * opd_ipos * shift) /\
|
---|
274 | float(ndoublefreq))
|
---|
275 |
|
---|
276 | # reflect spectrum about 0 to give unaliased version
|
---|
277 | reflected_spectrum = np.zeros([2*(len(baseline_spectrum)-1)],
|
---|
278 | np.complex)
|
---|
279 | reflected_spectrum[:len(baseline_spectrum)] = baseline_spectrum
|
---|
280 | reflected_spectrum[len(baseline_spectrum):] = \
|
---|
281 | baseline_spectrum[-2:0:-1]
|
---|
282 |
|
---|
283 | # apply phase shift and fft
|
---|
284 | reflected_spectrum *= shift
|
---|
285 | spectrum_fft = np.fft.ifft(reflected_spectrum)
|
---|
286 | config = config._replace(data = spectrum_fft[0].real)
|
---|
287 |
|
---|
288 | observed_obs_framework.append(config)
|
---|
289 |
|
---|
290 | previous_config = config
|
---|
291 |
|
---|
292 | self.result['observed_framework'] = observed_obs_framework
|
---|
293 | return self.result
|
---|
294 |
|
---|
295 | def __repr__(self):
|
---|
296 | return 'Observe'
|
---|
297 |
|
---|