Changeset 31
- Timestamp:
- May 20, 2014 3:22:42 PM (10 years ago)
- File:
-
- 1 edited
-
trunk/doublefourier.py (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/doublefourier.py
r17 r31 2 2 3 3 import collections 4 import math 4 5 import numpy as np 5 6 import pp 6 7 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() 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 44 25 45 26 … … 48 29 """ 49 30 50 def __init__(self, parameters, previous_results ):31 def __init__(self, parameters, previous_results, job_server): 51 32 self.parameters = parameters 52 33 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() 34 self.job_server = job_server 59 35 60 36 self.result = collections.OrderedDict() … … 64 40 65 41 fts = self.previous_results['fts'] 66 f requency_axis= fts['fts_wn']42 fts_wn = fts['fts_wn'] 67 43 opd_max = fts['opd_max'] 68 44 fts_nsample = fts['ftsnsample'] … … 81 57 skymodel = skygenerator['sky model'] 82 58 spatial_axis = self.result['spatial axis'] = skygenerator['spatial axis'] 83 self.result['frequency axis'] = f requency_axis59 self.result['frequency axis'] = fts_wn 84 60 85 61 # assuming nx is even then transform has 0 freq at origin and [nx/2] is … … 100 76 # FTS path diff and possibly baseline itself vary with time 101 77 for tindex,t in enumerate(times): 102 # for tindex,t in enumerate(times[:1]): 103 78 79 # calculate the sky that the system is observing at this 80 # moment, incorporating various errors in turn 104 81 sky_now = skymodel 105 106 # what is the system looking at?107 # add 'errors' in order108 82 109 83 # 1. baseline should be perp to centre of field. … … 113 87 # for now assume 0 error but do full calculation for timing 114 88 # purposes 115 baseline_dx = 0.0 116 baseline_dy = 0.0 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 117 99 118 100 # calculate xpos, ypos in units of pixel - numpy arrays 119 101 # [row,col] 120 102 nx = len(spatial_axis) 121 colpos = float(nx-1) * float (baseline_dx - spatial_axis[0])/ \103 colpos = float(nx-1) * bx_error / \ 122 104 (spatial_axis[-1] - spatial_axis[0]) 123 rowpos = float(nx-1) * float (baseline_dy - spatial_axis[0])/ \105 rowpos = float(nx-1) * by_error / \ 124 106 (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 107 131 108 # calculate fourier phase shift to move point at [0,0] to 132 109 # [rowpos, colpos] 110 print 'shift', colpos, rowpos 133 111 shiftx = np.zeros([nx], np.complex) 134 112 shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex) … … 147 125 shift[:,i] *= shifty 148 126 127 jobs = {} 149 128 # 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) 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 165 139 166 140 if t == times[0]: 167 141 self.result['sky at time 0'] = sky_now 168 142 169 143 # 2. telescopes should be centred on centre of field 170 144 # Telescopes collect flux from the 'sky' and pass … … 176 150 177 151 # multiply sky by amplitude beam 1 * amplitude beam 2 178 for iwn,wn in enumerate(frequency_axis): 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): 179 156 # calculate shifted beams here 180 157 # for now assume no errors and just use beam 181 158 # 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 159 pass 160 161 # multiply sky by amplitude beams of 2 antennas 162 sky_now *= amp_beam_1 * amp_beam_2 188 163 189 164 if t == times[0]: … … 192 167 # 3. baseline error revisited 193 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 194 177 baseline_error = 0.0 195 178 baseline_now = baseline + baseline_error 196 # print 'baseline_now', baseline_now 197 198 # gamma_12[(baseline,t)] = 0.0 179 199 180 fft_now = np.zeros(np.shape(sky_now), np.complex) 200 spectrum = np.zeros(np.shape(f requency_axis), np.complex)201 for iwn,wn in enumerate(f requency_axis):181 spectrum = np.zeros(np.shape(fts_wn), np.complex) 182 for iwn,wn in enumerate(fts_wn): 202 183 203 184 # derive shift needed to place baseline at one of FFT coords 204 185 # this depends on physical baseline and frequency 205 186 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]) / \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]) / \ 211 192 (spatial_freq_axis[-1] - spatial_freq_axis[0]) 212 rowpos = float(nx-1) * float(baseline_now_lambdas[1] - spatial_freq_axis[0]) / \ 193 rowpos = float(nx-1) * \ 194 float(baseline_now_lambdas[1] - spatial_freq_axis[0]) / \ 213 195 (spatial_freq_axis[-1] - spatial_freq_axis[0]) 214 # print 'spatial colpos, rowpos', colpos, rowpos 196 colpos = 0.0 197 rowpos = 0.0 198 print 'spatial colpos, rowpos', colpos, rowpos 215 199 216 200 # calculate fourier phase shift to move point at [rowpos,colpos] to … … 219 203 shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex) 220 204 shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex) 221 shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / float(nx)) 205 shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / \ 206 float(nx)) 222 207 223 208 shifty = np.zeros([nx], np.complex) … … 238 223 # 2d fft 239 224 temp = np.fft.fft2(temp) 240 # shift 0 freq to centre of array241 fft_now[:,:,iwn] = np.fft.ifftshift(temp) 242 225 fft_now[:,:,iwn] = temp 226 227 # set amp/phase at this frequency 243 228 spectrum[iwn] = temp[0,0] 244 229 245 230 if t == times[0]: 246 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 247 238 248 239 # 3. FTS sampling should be accurate 249 240 # derive lag due to FTS path difference 250 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 251 260 mirror_error = 0.0 252 261 opd = 2.0 * (vdrive * t + mirror_error) … … 255 264 256 265 # make spectrum symmetric 257 nfreq = len(f requency_axis)266 nfreq = len(fts_wn) 258 267 symmetric_spectrum = np.zeros([2*nfreq]) 259 268 symmetric_spectrum[:nfreq] = spectrum … … 271 280 measurement[tindex] = spectrum_fft[0] 272 281 273 self.result['baseline interferograms'][tuple(baseline)] = measurement 282 self.result['baseline interferograms'][tuple(baseline)] = \ 283 measurement 274 284 275 285 return self.result
Note:
See TracChangeset
for help on using the changeset viewer.
