| 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 |
|
|---|
