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