[17] | 1 | from __future__ import absolute_import
|
---|
| 2 |
|
---|
| 3 | import collections
|
---|
| 4 | import math
|
---|
| 5 | import numpy as np
|
---|
| 6 | import scipy.constants as sc
|
---|
| 7 |
|
---|
[33] | 8 | import common.commonobjects as co
|
---|
| 9 |
|
---|
[17] | 10 | def black_body(temperature, wavenumber):
|
---|
| 11 | """Function to calculate Planck function.
|
---|
| 12 | temperature - Kelvin
|
---|
| 13 | wavenumber - frequency in cm-1
|
---|
| 14 | """
|
---|
| 15 | freq = wavenumber * sc.c * 100.0
|
---|
[38] | 16 | if freq > 0:
|
---|
| 17 | jnu = 2.0 * sc.h * pow(freq,3) / (pow(sc.c,2) *
|
---|
| 18 | (math.exp((sc.h * freq) / (sc.k * temperature)) - 1.0))
|
---|
| 19 | else:
|
---|
| 20 | jnu = 0.0
|
---|
[17] | 21 |
|
---|
| 22 | return jnu
|
---|
| 23 |
|
---|
[38] | 24 | class BB_spectrum(object):
|
---|
| 25 | """Class to generate BB spectrum.
|
---|
| 26 | """
|
---|
| 27 | def __init__(self, temperature, frequency_axis, cutoffmin, cutoffmax,
|
---|
| 28 | emissivity):
|
---|
| 29 | self.temperature = temperature
|
---|
| 30 | self.frequency_axis = frequency_axis
|
---|
| 31 | self.cutoffmin = cutoffmin
|
---|
| 32 | self.cutoffmax = cutoffmax
|
---|
| 33 | self.emissivity = emissivity
|
---|
[17] | 34 |
|
---|
[38] | 35 | def calculate(self):
|
---|
| 36 | spectrum = np.zeros(np.shape(self.frequency_axis))
|
---|
| 37 | # ignore floating point errors
|
---|
| 38 | old_settings = np.seterr(all='ignore')
|
---|
| 39 |
|
---|
| 40 | for iwn,wn in enumerate(self.frequency_axis):
|
---|
| 41 | # simulate real-life 'rounded' cutoffs numerically
|
---|
| 42 | f1 = 1.0 / (1.0 + pow(self.cutoffmin/wn, 18) +
|
---|
| 43 | pow(wn/self.cutoffmax, 24))
|
---|
| 44 | spectrum[iwn] = black_body(self.temperature, wn) * f1 * \
|
---|
| 45 | self.emissivity
|
---|
| 46 |
|
---|
| 47 | # restore fp behaviour
|
---|
| 48 | ignore = np.seterr(**old_settings)
|
---|
| 49 |
|
---|
| 50 | return spectrum
|
---|
| 51 |
|
---|
| 52 |
|
---|
[45] | 53 | class linespectrum(object):
|
---|
[38] | 54 | """Class to return spectrum with line at wn.
|
---|
| 55 | """
|
---|
[45] | 56 | def __init__(self, temperature, linefreq, frequency_axis):
|
---|
| 57 | self.temperature = temperature
|
---|
| 58 | self.linefreq = linefreq
|
---|
[38] | 59 | self.frequency_axis = frequency_axis
|
---|
| 60 |
|
---|
| 61 | def calculate(self):
|
---|
[45] | 62 | # spectrum is a sync function centred on the specified line freq
|
---|
| 63 | x = (self.frequency_axis - self.linefreq) / \
|
---|
| 64 | (self.frequency_axis[1] - self.frequency_axis[0])
|
---|
| 65 | spectrum = np.sinc(x) * self.temperature
|
---|
| 66 |
|
---|
[38] | 67 | return spectrum
|
---|
| 68 |
|
---|
| 69 |
|
---|
[17] | 70 | class SkyGenerator(object):
|
---|
| 71 | """Class to generate a model sky.
|
---|
| 72 | """
|
---|
| 73 |
|
---|
| 74 | def __init__(self, parameters, previous_results):
|
---|
| 75 | self.parameters = parameters
|
---|
| 76 | self.previous_results = previous_results
|
---|
| 77 | self.result = collections.OrderedDict()
|
---|
| 78 |
|
---|
| 79 | def run(self):
|
---|
| 80 | print 'SkyGenerator.run'
|
---|
| 81 |
|
---|
| 82 | fts = self.previous_results['fts']
|
---|
[38] | 83 | fts_wn_truncated = fts['fts_wn_truncated']
|
---|
[45] | 84 | cutoffmin = fts['wnmin']
|
---|
| 85 | cutoffmax = fts['wnmax']
|
---|
[17] | 86 |
|
---|
| 87 | beamsgenerator = self.previous_results['beamsgenerator']
|
---|
| 88 | npix = beamsgenerator['npix']
|
---|
| 89 | spatial_axis = beamsgenerator['spatial axis [arcsec]']
|
---|
| 90 |
|
---|
[33] | 91 | # skymodel is complex so that its fft can hold truncated version
|
---|
| 92 | # of infinitesimally sampled map - does that make sense?
|
---|
[38] | 93 | skymodel = np.zeros([npix, npix, len(fts_wn_truncated)], np.complex)
|
---|
[17] | 94 |
|
---|
| 95 | sky = self.parameters['substages']['Sky']
|
---|
[45] | 96 | columns = sky['sourcenum'].keys()
|
---|
[17] | 97 |
|
---|
| 98 | self.result['sources'] = collections.OrderedDict()
|
---|
| 99 |
|
---|
| 100 | for column in columns:
|
---|
[45] | 101 | temp = sky['sourcenum'][column]
|
---|
[17] | 102 | sourcenum = int(round(temp))
|
---|
| 103 | if sourcenum not in self.result['sources'].keys():
|
---|
| 104 | self.result['sources'][sourcenum] = {}
|
---|
| 105 |
|
---|
[45] | 106 | type = sky['type'][column]
|
---|
[17] | 107 |
|
---|
| 108 | temp = sky['x pos [asec]'][column]
|
---|
| 109 | xpos = float(temp)
|
---|
| 110 |
|
---|
| 111 | temp = sky['y pos [asec]'][column]
|
---|
| 112 | ypos = float(temp)
|
---|
| 113 |
|
---|
[45] | 114 | temp = sky['xwidth'][column]
|
---|
| 115 | try:
|
---|
| 116 | xwidth = float(temp)
|
---|
| 117 | except:
|
---|
| 118 | xwidth = None
|
---|
| 119 |
|
---|
| 120 | temp = sky['ywidth'][column]
|
---|
| 121 | try:
|
---|
| 122 | ywidth = float(temp)
|
---|
| 123 | except:
|
---|
| 124 | ywidth = None
|
---|
| 125 |
|
---|
| 126 | spectrum = sky['spectrum'][column]
|
---|
| 127 |
|
---|
| 128 | temp = sky['temperature'][column]
|
---|
[17] | 129 | temperature = float(temp)
|
---|
| 130 |
|
---|
[45] | 131 | temp = sky['linefreq'][column]
|
---|
| 132 | try:
|
---|
| 133 | linefreq = float(temp)
|
---|
| 134 | except:
|
---|
| 135 | linefreq = None
|
---|
[17] | 136 |
|
---|
| 137 | temp = sky['emissivity'][column]
|
---|
| 138 | emissivity = float(temp)
|
---|
| 139 |
|
---|
[45] | 140 | print 'generating source:%s type:%s xpos:%s ypos:%s' % (sourcenum,
|
---|
| 141 | type, xpos, ypos)
|
---|
| 142 | if type.upper().strip() == 'GAUSSIAN':
|
---|
| 143 | print ' fwhmx:%s fwhmy:%s' % (xwidth, ywidth)
|
---|
[17] | 144 |
|
---|
[45] | 145 | if spectrum.upper().strip() == 'BLACKBODY':
|
---|
| 146 | print ' blackbody spectrum temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
|
---|
| 147 | temperature, cutoffmin, cutoffmax, emissivity)
|
---|
| 148 | spectrum_func = BB_spectrum(temperature, fts_wn_truncated,
|
---|
| 149 | cutoffmin, cutoffmax, emissivity)
|
---|
[38] | 150 |
|
---|
[45] | 151 | elif spectrum.upper().strip() == 'LINE':
|
---|
| 152 | print ' line spectrum brightness temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
|
---|
| 153 | temperature, cutoffmin, cutoffmax, emissivity)
|
---|
| 154 | spectrum_func = linespectrum(temperature, linefreq, fts_wn_truncated)
|
---|
| 155 |
|
---|
[17] | 156 | if type.upper().strip() == 'POINT':
|
---|
| 157 | source_spectrum = self._create_point_source(xpos, ypos,
|
---|
[38] | 158 | skymodel, spatial_axis, fts_wn_truncated, spectrum_func)
|
---|
[45] | 159 | elif type.upper().strip() == 'GAUSSIAN':
|
---|
| 160 | source_spectrum = self._create_gaussian_source(xpos, ypos,
|
---|
| 161 | xwidth, ywidth, skymodel, spatial_axis, fts_wn_truncated,
|
---|
| 162 | spectrum_func)
|
---|
[17] | 163 | else:
|
---|
| 164 | source_spectrum = None
|
---|
| 165 | print "source type '%s' not yet implemented" % type
|
---|
| 166 |
|
---|
| 167 | self.result['sources'][sourcenum]['spectrum'] = source_spectrum
|
---|
| 168 |
|
---|
| 169 | self.result['sky model'] = skymodel
|
---|
| 170 | self.result['spatial axis'] = spatial_axis
|
---|
[38] | 171 | self.result['frequency axis'] = fts_wn_truncated
|
---|
[17] | 172 |
|
---|
| 173 | return self.result
|
---|
| 174 |
|
---|
[45] | 175 | def _create_gaussian_source(self, xpos, ypos, xwidth, ywidth,
|
---|
| 176 | skymodel, spatial_axis, frequency_axis, spectrum_function):
|
---|
| 177 | """Create a point source.
|
---|
| 178 | """
|
---|
| 179 | # calculate spectrum
|
---|
| 180 | spectrum = spectrum_function.calculate()
|
---|
| 181 |
|
---|
| 182 | # calculate Gaussian profile - a cunning use of array slicing found
|
---|
| 183 | # on the web
|
---|
| 184 | x = spatial_axis
|
---|
| 185 | y = x[:,np.newaxis]
|
---|
| 186 | profile = np.exp(-4*np.log(2) * ((x-xpos)**2/xwidth**2 + (y-ypos)**2/ywidth**2))
|
---|
| 187 |
|
---|
| 188 | # # apodise
|
---|
| 189 | # profile[((x-xpos)**2 + (y-ypos)**2) < 4] = 1.0
|
---|
| 190 |
|
---|
| 191 | # go through freq planes and add source to each
|
---|
| 192 | for iwn,wn in enumerate(frequency_axis):
|
---|
| 193 | # add to sky model
|
---|
| 194 | skymodel[:,:,iwn] += profile * spectrum[iwn]
|
---|
| 195 |
|
---|
| 196 | # return spectrum
|
---|
| 197 | axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
|
---|
| 198 | spectrum = co.Spectrum(data=spectrum, axis=axis,
|
---|
| 199 | title='Source spectrum', units='W sr-1 m-2 Hz-1')
|
---|
| 200 |
|
---|
| 201 | return spectrum
|
---|
| 202 |
|
---|
[38] | 203 | def _create_point_source(self, xpos, ypos, skymodel, spatial_axis,
|
---|
| 204 | frequency_axis, spectrum_function):
|
---|
[17] | 205 | """Create a point source.
|
---|
| 206 | """
|
---|
| 207 |
|
---|
| 208 | # calculate xpos, ypos in units of pixel - numpy arrays [row,col]
|
---|
| 209 | nx = len(spatial_axis)
|
---|
| 210 | colpos = float(nx-1) * float (xpos - spatial_axis[0]) / (spatial_axis[-1] - spatial_axis[0])
|
---|
| 211 | rowpos = float(nx-1) * float (ypos - spatial_axis[0]) / (spatial_axis[-1] - spatial_axis[0])
|
---|
| 212 |
|
---|
| 213 | if colpos < 0 or colpos > (nx-1) or rowpos < 0 or rowpos > (nx-1):
|
---|
| 214 | # point source is outside modelled area
|
---|
| 215 | return
|
---|
| 216 |
|
---|
[33] | 217 | # calculate fourier phase shift to move point at [0,0] to
|
---|
| 218 | # [rowpos, colpos]
|
---|
[17] | 219 | shiftx = np.zeros([nx], np.complex)
|
---|
| 220 | shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
| 221 | shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
| 222 | shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / float(nx))
|
---|
| 223 |
|
---|
| 224 | shifty = np.zeros([nx], np.complex)
|
---|
| 225 | shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
| 226 | shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
| 227 | shifty = np.exp((-2.0j * np.pi * rowpos * shifty) / float(nx))
|
---|
| 228 |
|
---|
| 229 | shift = np.ones([nx,nx], np.complex)
|
---|
| 230 | for j in range(nx):
|
---|
| 231 | shift[j,:] *= shiftx
|
---|
| 232 | for i in range(nx):
|
---|
| 233 | shift[:,i] *= shifty
|
---|
| 234 |
|
---|
[38] | 235 | # calculate spectrum
|
---|
| 236 | spectrum = spectrum_function.calculate()
|
---|
[17] | 237 |
|
---|
| 238 | # go through freq planes and add point source to each
|
---|
| 239 | for iwn,wn in enumerate(frequency_axis):
|
---|
| 240 | # create point in frequency space
|
---|
| 241 | temp = np.zeros([nx,nx])
|
---|
[38] | 242 | temp[0,0] = spectrum[iwn]
|
---|
[17] | 243 | # 2d fft
|
---|
| 244 | temp = np.fft.fft2(temp)
|
---|
| 245 | # apply phase shift to move point to required offset
|
---|
| 246 | temp *= shift
|
---|
| 247 | # transform back
|
---|
| 248 | temp = np.fft.ifft2(temp)
|
---|
| 249 |
|
---|
| 250 | # add to sky model
|
---|
| 251 | skymodel[:,:,iwn] += temp
|
---|
| 252 |
|
---|
[33] | 253 | # return spectrum
|
---|
| 254 | axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
|
---|
[38] | 255 | spectrum = co.Spectrum(data=spectrum, axis=axis,
|
---|
| 256 | title='Source spectrum', units='W sr-1 m-2 Hz-1')
|
---|
[17] | 257 |
|
---|
[33] | 258 | return spectrum
|
---|
| 259 |
|
---|
[17] | 260 | def __repr__(self):
|
---|
| 261 | return 'SkyGenerator'
|
---|
| 262 |
|
---|