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 |
|
---|
11 | class DirtyImage(object):
|
---|
12 | """Class to compute dirty image cube from uvspectra.
|
---|
13 | """
|
---|
14 |
|
---|
15 | def __init__(self, previous_results, job_server):
|
---|
16 | self.previous_results = previous_results
|
---|
17 | self.job_server = job_server
|
---|
18 |
|
---|
19 | self.result = collections.OrderedDict()
|
---|
20 |
|
---|
21 | def run(self):
|
---|
22 | print 'DirtyImage.run'
|
---|
23 |
|
---|
24 | # get relevant fts info
|
---|
25 | fts = self.previous_results['fts']
|
---|
26 | fts_wnmin = fts['wnmin']
|
---|
27 |
|
---|
28 | # get primary beam info
|
---|
29 | beamsgenerator = self.previous_results['beamsgenerator']
|
---|
30 | npix = beamsgenerator['npix']
|
---|
31 | spatial_axis = beamsgenerator['spatial axis [arcsec]']
|
---|
32 |
|
---|
33 | # get observation list
|
---|
34 | uvspectra = self.previous_results['uvspectra']['uvspectra']
|
---|
35 | wavenumber = uvspectra[0].wavenumber
|
---|
36 | wavenumber = wavenumber[wavenumber>fts_wnmin]
|
---|
37 |
|
---|
38 | # dirty image cube
|
---|
39 | dirtyimage = np.zeros([npix, npix, len(wavenumber)], np.float)
|
---|
40 | dirtybeam = np.zeros([npix, npix, len(wavenumber)], np.float)
|
---|
41 |
|
---|
42 | for uvspectrum in uvspectra:
|
---|
43 |
|
---|
44 | b_x = uvspectrum.baseline_x
|
---|
45 | b_y = uvspectrum.baseline_y
|
---|
46 | spectrum = uvspectrum.spectrum
|
---|
47 | wn_spectrum = uvspectrum.wavenumber
|
---|
48 |
|
---|
49 | for iwn,wn in enumerate(wavenumber):
|
---|
50 | # must be a better way of doing this
|
---|
51 | argx = b_x * (wn * 100.0) * spatial_axis / 206265.0
|
---|
52 | argy = b_y * (wn * 100.0) * spatial_axis / 206265.0
|
---|
53 | valx = np.exp(2j * math.pi * argx)
|
---|
54 | valy = np.exp(2j * math.pi * argy)
|
---|
55 |
|
---|
56 | vis = spectrum[wn_spectrum==wn]
|
---|
57 |
|
---|
58 | # numpy arrays are [row,col] or here [y,x]
|
---|
59 | for iy,y in enumerate(spatial_axis):
|
---|
60 | contribution = (vis * valx * valy[iy]) + \
|
---|
61 | np.conjugate(vis * valx * valy[iy])
|
---|
62 | dirtyimage[iy,:,iwn] += contribution.real
|
---|
63 | contribution = (valx * valy[iy]) + \
|
---|
64 | np.conjugate(valx * valy[iy])
|
---|
65 | dirtybeam[iy,:,iwn] += contribution.real
|
---|
66 |
|
---|
67 | # normalise
|
---|
68 | for iwn,wn in enumerate(wavenumber):
|
---|
69 | dirtyimage[:,:,iwn] /= np.sum(dirtybeam[:,:,iwn])
|
---|
70 | dirtybeam[:,:,iwn] /= np.sum(dirtybeam[:,:,iwn])
|
---|
71 | # dirtyimage /= 2.0 * float(len(uvspectra))
|
---|
72 | # dirtybeam /= 2.0 * float(len(uvspectra))
|
---|
73 | print 'sum image', np.sum(dirtyimage[:,:,0])
|
---|
74 | print 'sum beam', np.sum(dirtybeam[:,:,0])
|
---|
75 |
|
---|
76 | self.result['dirtyimage'] = dirtyimage
|
---|
77 | self.result['dirtybeam'] = dirtybeam
|
---|
78 | self.result['spatial axis [arcsec]'] = spatial_axis
|
---|
79 | self.result['spatial axis'] = spatial_axis
|
---|
80 | self.result['wavenumber [cm-1]'] = wavenumber
|
---|
81 |
|
---|
82 | return self.result
|
---|
83 |
|
---|
84 | def __repr__(self):
|
---|
85 | return '''
|
---|
86 | DirtyImage:
|
---|
87 | '''.format(
|
---|
88 | num_uvspectra=len(self.uvspectra))
|
---|
89 |
|
---|