| 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 |
|
|---|
