1 | from __future__ import absolute_import
|
---|
2 |
|
---|
3 | import collections
|
---|
4 | import math
|
---|
5 | import numpy as np
|
---|
6 | import pp
|
---|
7 |
|
---|
8 | import common.commonobjects as co
|
---|
9 |
|
---|
10 | def dirty(b_x_list, b_y_list, spectra, wn_spectra, iwn, wn, spatial_axis, npix):
|
---|
11 | image = numpy.zeros([npix, npix], numpy.float)
|
---|
12 | beam = numpy.zeros([npix, npix], numpy.float)
|
---|
13 |
|
---|
14 | for ibx,b_x in enumerate(b_x_list):
|
---|
15 | b_x = b_x_list[ibx]
|
---|
16 | b_y = b_y_list[ibx]
|
---|
17 | spectrum = spectra[ibx]
|
---|
18 | wn_spectrum = wn_spectra[ibx]
|
---|
19 |
|
---|
20 | # must be a better way of doing this
|
---|
21 | argx = numpy.radians(b_x * (wn * 100.0) * spatial_axis / 3600.0)
|
---|
22 | argy = numpy.radians(b_y * (wn * 100.0) * spatial_axis / 3600.0)
|
---|
23 | valx = numpy.exp(2.0j * math.pi * argx)
|
---|
24 | valy = numpy.exp(2.0j * math.pi * argy)
|
---|
25 |
|
---|
26 | vis = spectrum[wn_spectrum==wn]
|
---|
27 |
|
---|
28 | # numpy arrays are [row,col] or here [y,x]
|
---|
29 | for iy,y in enumerate(spatial_axis):
|
---|
30 | contribution = (vis * valx * valy[iy]) + \
|
---|
31 | numpy.conjugate(vis * valx * valy[iy])
|
---|
32 | image[iy,:] += contribution.real
|
---|
33 | contribution = (valx * valy[iy]) + \
|
---|
34 | numpy.conjugate(valx * valy[iy])
|
---|
35 | beam[iy,:] += contribution.real
|
---|
36 |
|
---|
37 | # normalise
|
---|
38 | # The sum of the dirty beam should equal the 0 baseline vis measurement: 0
|
---|
39 | # Normalising the volume would preserve flux under convolution but
|
---|
40 | # is problematic as the volume is 0, so normalise the peak to 1.
|
---|
41 | beam /= 2.0 * len(b_x_list)
|
---|
42 |
|
---|
43 | # likewise for dirty map
|
---|
44 | image /= 2.0 * len(b_x_list)
|
---|
45 |
|
---|
46 | return image, beam
|
---|
47 |
|
---|
48 |
|
---|
49 | class DirtyImage(object):
|
---|
50 | """Class to compute dirty image cube from uvspectra.
|
---|
51 | """
|
---|
52 |
|
---|
53 | def __init__(self, previous_results, job_server):
|
---|
54 | self.previous_results = previous_results
|
---|
55 | self.job_server = job_server
|
---|
56 |
|
---|
57 | self.result = collections.OrderedDict()
|
---|
58 |
|
---|
59 | def run(self):
|
---|
60 | print 'DirtyImage.run'
|
---|
61 |
|
---|
62 | # get relevant fts info
|
---|
63 | fts = self.previous_results['fts']
|
---|
64 | fts_wnmin = fts['wnmin']
|
---|
65 |
|
---|
66 | # get primary beam info
|
---|
67 | beamsgenerator = self.previous_results['beamsgenerator']
|
---|
68 | npix = beamsgenerator['npix']
|
---|
69 | spatial_axis = beamsgenerator['spatial axis [arcsec]']
|
---|
70 |
|
---|
71 | # get observation list
|
---|
72 | uvspectra = self.previous_results['uvspectra']['uvspectra']
|
---|
73 | wavenumber = uvspectra[0].wavenumber
|
---|
74 | wavenumber = wavenumber[wavenumber>fts_wnmin]
|
---|
75 |
|
---|
76 | # dirty image cube
|
---|
77 | dirtyimage = np.zeros([npix, npix, len(wavenumber)], np.float)
|
---|
78 | dirtybeam = np.zeros([npix, npix, len(wavenumber)], np.float)
|
---|
79 |
|
---|
80 | # calculate dirty image for each wn
|
---|
81 | # uvspectrum objects don't pickle so unpack them first
|
---|
82 | b_x_list = []
|
---|
83 | b_y_list = []
|
---|
84 | spectra = []
|
---|
85 | wn_spectra = []
|
---|
86 | for uvspectrum in uvspectra:
|
---|
87 | b_x_list.append(uvspectrum.baseline_x)
|
---|
88 | b_y_list.append(uvspectrum.baseline_y)
|
---|
89 | spectra.append(uvspectrum.spectrum)
|
---|
90 | wn_spectra.append(uvspectrum.wavenumber)
|
---|
91 |
|
---|
92 | jobs = {}
|
---|
93 | for iwn,wn in enumerate(wavenumber):
|
---|
94 | # submit jobs
|
---|
95 | indata = (b_x_list, b_y_list, spectra, wn_spectra, iwn, wn,
|
---|
96 | spatial_axis, npix,)
|
---|
97 | jobs[wn] = self.job_server.submit(dirty, indata, (),
|
---|
98 | ('numpy','math',))
|
---|
99 |
|
---|
100 | for iwn,wn in enumerate(wavenumber):
|
---|
101 | # collect and store results
|
---|
102 | dirtyimage[:,:,iwn], dirtybeam[:,:,iwn] = jobs[wn]()
|
---|
103 |
|
---|
104 | self.result['dirtyimage'] = dirtyimage
|
---|
105 | self.result['dirtybeam'] = dirtybeam
|
---|
106 | self.result['spatial axis [arcsec]'] = spatial_axis
|
---|
107 | self.result['spatial axis'] = spatial_axis
|
---|
108 | self.result['wavenumber [cm-1]'] = wavenumber
|
---|
109 |
|
---|
110 | return self.result
|
---|
111 |
|
---|
112 | def __repr__(self):
|
---|
113 | return '''
|
---|
114 | DirtyImage:
|
---|
115 | '''.format(
|
---|
116 | num_uvspectra=len(self.uvspectra))
|
---|
117 |
|
---|