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 |
|
---|
8 | import common.commonobjects as co
|
---|
9 |
|
---|
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 | 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
|
---|
21 |
|
---|
22 | return jnu
|
---|
23 |
|
---|
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
|
---|
34 |
|
---|
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 |
|
---|
53 | class linespectrum(object):
|
---|
54 | """Class to return spectrum with line at wn.
|
---|
55 | """
|
---|
56 | def __init__(self, temperature, linefreq, frequency_axis):
|
---|
57 | self.temperature = temperature
|
---|
58 | self.linefreq = linefreq
|
---|
59 | self.frequency_axis = frequency_axis
|
---|
60 |
|
---|
61 | def calculate(self):
|
---|
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 |
|
---|
67 | return spectrum
|
---|
68 |
|
---|
69 |
|
---|
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']
|
---|
83 | fts_wn_truncated = fts['fts_wn_truncated']
|
---|
84 | cutoffmin = fts['wnmin']
|
---|
85 | cutoffmax = fts['wnmax']
|
---|
86 |
|
---|
87 | beamsgenerator = self.previous_results['beamsgenerator']
|
---|
88 | npix = beamsgenerator['npix']
|
---|
89 | spatial_axis = beamsgenerator['spatial axis [arcsec]']
|
---|
90 |
|
---|
91 | # skymodel is complex so that its fft can hold truncated version
|
---|
92 | # of infinitesimally sampled map - does that make sense?
|
---|
93 | skymodel = np.zeros([npix, npix, len(fts_wn_truncated)], np.complex)
|
---|
94 |
|
---|
95 | sky = self.parameters['substages']['Sky']
|
---|
96 | columns = sky['sourcenum'].keys()
|
---|
97 |
|
---|
98 | self.result['sources'] = collections.OrderedDict()
|
---|
99 |
|
---|
100 | for column in columns:
|
---|
101 | temp = sky['sourcenum'][column]
|
---|
102 | sourcenum = int(round(temp))
|
---|
103 | if sourcenum not in self.result['sources'].keys():
|
---|
104 | self.result['sources'][sourcenum] = {}
|
---|
105 |
|
---|
106 | type = sky['type'][column]
|
---|
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 |
|
---|
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]
|
---|
129 | temperature = float(temp)
|
---|
130 |
|
---|
131 | temp = sky['linefreq'][column]
|
---|
132 | try:
|
---|
133 | linefreq = float(temp)
|
---|
134 | except:
|
---|
135 | linefreq = None
|
---|
136 |
|
---|
137 | temp = sky['emissivity'][column]
|
---|
138 | emissivity = float(temp)
|
---|
139 |
|
---|
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)
|
---|
144 |
|
---|
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)
|
---|
150 |
|
---|
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 |
|
---|
156 | if type.upper().strip() == 'POINT':
|
---|
157 | source_spectrum = self._create_point_source(xpos, ypos,
|
---|
158 | skymodel, spatial_axis, fts_wn_truncated, spectrum_func)
|
---|
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)
|
---|
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
|
---|
171 | self.result['frequency axis'] = fts_wn_truncated
|
---|
172 |
|
---|
173 | return self.result
|
---|
174 |
|
---|
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 |
|
---|
203 | def _create_point_source(self, xpos, ypos, skymodel, spatial_axis,
|
---|
204 | frequency_axis, spectrum_function):
|
---|
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 |
|
---|
217 | # calculate fourier phase shift to move point at [0,0] to
|
---|
218 | # [rowpos, colpos]
|
---|
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 |
|
---|
235 | # calculate spectrum
|
---|
236 | spectrum = spectrum_function.calculate()
|
---|
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])
|
---|
242 | temp[0,0] = spectrum[iwn]
|
---|
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 |
|
---|
253 | # return spectrum
|
---|
254 | axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
|
---|
255 | spectrum = co.Spectrum(data=spectrum, axis=axis,
|
---|
256 | title='Source spectrum', units='W sr-1 m-2 Hz-1')
|
---|
257 |
|
---|
258 | return spectrum
|
---|
259 |
|
---|
260 | def __repr__(self):
|
---|
261 | return 'SkyGenerator'
|
---|
262 |
|
---|