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