[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
|
---|
| 16 | jnu = 2.0 * sc.h * pow(freq,3) / (pow(sc.c,2) *
|
---|
| 17 | (math.exp((sc.h * freq) / (sc.k * temperature)) - 1.0))
|
---|
| 18 |
|
---|
| 19 | return jnu
|
---|
| 20 |
|
---|
| 21 |
|
---|
| 22 | class SkyGenerator(object):
|
---|
| 23 | """Class to generate a model sky.
|
---|
| 24 | """
|
---|
| 25 |
|
---|
| 26 | def __init__(self, parameters, previous_results):
|
---|
| 27 | self.parameters = parameters
|
---|
| 28 | self.previous_results = previous_results
|
---|
| 29 | self.result = collections.OrderedDict()
|
---|
| 30 |
|
---|
| 31 | def run(self):
|
---|
| 32 | print 'SkyGenerator.run'
|
---|
| 33 |
|
---|
| 34 | fts = self.previous_results['fts']
|
---|
| 35 | frequency_axis = fts['fts_wn']
|
---|
| 36 | nspec = len(frequency_axis)
|
---|
| 37 |
|
---|
| 38 | beamsgenerator = self.previous_results['beamsgenerator']
|
---|
| 39 | npix = beamsgenerator['npix']
|
---|
| 40 | spatial_axis = beamsgenerator['spatial axis [arcsec]']
|
---|
| 41 |
|
---|
[33] | 42 | # skymodel is complex so that its fft can hold truncated version
|
---|
| 43 | # of infinitesimally sampled map - does that make sense?
|
---|
| 44 | skymodel = np.zeros([npix, npix, nspec], np.complex)
|
---|
[17] | 45 |
|
---|
| 46 | sky = self.parameters['substages']['Sky']
|
---|
| 47 | columns = sky['SourceNum'].keys()
|
---|
| 48 |
|
---|
| 49 | self.result['sources'] = collections.OrderedDict()
|
---|
| 50 |
|
---|
| 51 | for column in columns:
|
---|
| 52 | temp = sky['SourceNum'][column]
|
---|
| 53 | sourcenum = int(round(temp))
|
---|
| 54 | if sourcenum not in self.result['sources'].keys():
|
---|
| 55 | self.result['sources'][sourcenum] = {}
|
---|
| 56 |
|
---|
| 57 | type = sky['Type'][column]
|
---|
| 58 |
|
---|
| 59 | temp = sky['x pos [asec]'][column]
|
---|
| 60 | xpos = float(temp)
|
---|
| 61 |
|
---|
| 62 | temp = sky['y pos [asec]'][column]
|
---|
| 63 | ypos = float(temp)
|
---|
| 64 |
|
---|
| 65 | temp = sky['Temp'][column]
|
---|
| 66 | temperature = float(temp)
|
---|
| 67 |
|
---|
| 68 | temp = sky['cutoffmin'][column]
|
---|
| 69 | cutoffmin = float(temp)
|
---|
| 70 |
|
---|
| 71 | temp = sky['cutoffmax'][column]
|
---|
| 72 | cutoffmax = float(temp)
|
---|
| 73 |
|
---|
| 74 | temp = sky['emissivity'][column]
|
---|
| 75 | emissivity = float(temp)
|
---|
| 76 |
|
---|
| 77 | print 'generating source:%s type:%s xpos:%s ypos:%s temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
|
---|
| 78 | sourcenum, type, xpos, ypos, temperature, cutoffmin, cutoffmax,
|
---|
| 79 | emissivity)
|
---|
| 80 |
|
---|
| 81 | if type.upper().strip() == 'POINT':
|
---|
| 82 | source_spectrum = self._create_point_source(xpos, ypos,
|
---|
| 83 | temperature,
|
---|
| 84 | cutoffmin, cutoffmax, emissivity, skymodel,
|
---|
| 85 | spatial_axis, frequency_axis)
|
---|
| 86 | else:
|
---|
| 87 | source_spectrum = None
|
---|
| 88 | print "source type '%s' not yet implemented" % type
|
---|
| 89 |
|
---|
| 90 | self.result['sources'][sourcenum]['spectrum'] = source_spectrum
|
---|
| 91 |
|
---|
| 92 | self.result['sky model'] = skymodel
|
---|
| 93 | self.result['spatial axis'] = spatial_axis
|
---|
| 94 | self.result['frequency axis'] = frequency_axis
|
---|
| 95 |
|
---|
| 96 | return self.result
|
---|
| 97 |
|
---|
| 98 | def _create_point_source(self, xpos, ypos, temperature, cutoffmin,
|
---|
| 99 | cutoffmax, emissivity, skymodel, spatial_axis, frequency_axis):
|
---|
| 100 | """Create a point source.
|
---|
| 101 | """
|
---|
| 102 |
|
---|
| 103 | # calculate xpos, ypos in units of pixel - numpy arrays [row,col]
|
---|
| 104 | nx = len(spatial_axis)
|
---|
| 105 | colpos = float(nx-1) * float (xpos - spatial_axis[0]) / (spatial_axis[-1] - spatial_axis[0])
|
---|
| 106 | rowpos = float(nx-1) * float (ypos - spatial_axis[0]) / (spatial_axis[-1] - spatial_axis[0])
|
---|
| 107 |
|
---|
| 108 | if colpos < 0 or colpos > (nx-1) or rowpos < 0 or rowpos > (nx-1):
|
---|
| 109 | # point source is outside modelled area
|
---|
| 110 | return
|
---|
| 111 |
|
---|
[33] | 112 | # calculate fourier phase shift to move point at [0,0] to
|
---|
| 113 | # [rowpos, colpos]
|
---|
[17] | 114 | shiftx = np.zeros([nx], np.complex)
|
---|
| 115 | shiftx[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
| 116 | shiftx[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
| 117 | shiftx = np.exp((-2.0j * np.pi * colpos * shiftx) / float(nx))
|
---|
| 118 |
|
---|
| 119 | shifty = np.zeros([nx], np.complex)
|
---|
| 120 | shifty[:nx/2] = np.arange(nx/2, dtype=np.complex)
|
---|
| 121 | shifty[nx/2:] = np.arange(-nx/2, 0, dtype=np.complex)
|
---|
| 122 | shifty = np.exp((-2.0j * np.pi * rowpos * shifty) / float(nx))
|
---|
| 123 |
|
---|
| 124 | shift = np.ones([nx,nx], np.complex)
|
---|
| 125 | for j in range(nx):
|
---|
| 126 | shift[j,:] *= shiftx
|
---|
| 127 | for i in range(nx):
|
---|
| 128 | shift[:,i] *= shifty
|
---|
| 129 |
|
---|
| 130 | # calculate black-body spectrum, modified for cutoffs in
|
---|
| 131 | # sensitivity and source emissivity
|
---|
| 132 | bbmod = np.zeros(np.shape(frequency_axis))
|
---|
| 133 | # ignore floating point errors
|
---|
| 134 | old_settings = np.seterr(all='ignore')
|
---|
| 135 | for iwn,wn in enumerate(frequency_axis):
|
---|
| 136 | # simulate real-life 'rounded' cutoffs numerically
|
---|
| 137 | f1 = 1.0 / (1.0 + pow(cutoffmin/wn, 18) + pow(wn/cutoffmax, 24))
|
---|
| 138 | bbmod[iwn] = black_body(temperature, wn) * f1 * emissivity
|
---|
| 139 | # restore fp behaviour
|
---|
| 140 | ignore = np.seterr(**old_settings)
|
---|
| 141 |
|
---|
| 142 | # go through freq planes and add point source to each
|
---|
| 143 | for iwn,wn in enumerate(frequency_axis):
|
---|
| 144 | # create point in frequency space
|
---|
| 145 | temp = np.zeros([nx,nx])
|
---|
| 146 | temp[0,0] = bbmod[iwn]
|
---|
| 147 | # 2d fft
|
---|
| 148 | temp = np.fft.fft2(temp)
|
---|
| 149 | # apply phase shift to move point to required offset
|
---|
| 150 | temp *= shift
|
---|
| 151 | # transform back
|
---|
| 152 | temp = np.fft.ifft2(temp)
|
---|
| 153 | temp = np.real(temp)
|
---|
| 154 |
|
---|
| 155 | # add to sky model
|
---|
| 156 | skymodel[:,:,iwn] += temp
|
---|
| 157 |
|
---|
[33] | 158 | # return spectrum
|
---|
| 159 | axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
|
---|
| 160 | spectrum = co.Spectrum(data=bbmod, axis=axis, title='Source spectrum',
|
---|
| 161 | units='W sr-1 m-2 Hz-1')
|
---|
[17] | 162 |
|
---|
[33] | 163 | return spectrum
|
---|
| 164 |
|
---|
[17] | 165 | def __repr__(self):
|
---|
| 166 | return 'SkyGenerator'
|
---|
| 167 |
|
---|