1 from datetime import datetime
2 import numpy as np
3 import pylab as pl
4 import chianti.filters as chfilters
5 import chianti.util as util
6 import chianti.io as chio
7 import chianti.data as chdata
8 import chianti.constants as const
9
10 defaults = chdata.Defaults
11
13 '''
14 a collection of methods for use in spectrum calculations
15 '''
16 - def __init__(self, temperature, density):
17 self.Temperature = temperature
18 self.EDensity = density
19 self.AbundanceName = defaults['abundfile']
20 self.AbundAll = chdata.Abundance[self.AbundanceName]['abundance']
21
22
23
24 - def convolve(self, wavelength=0, filter=(chfilters.gaussianR, 1000.), label=0, verbose=0):
25 '''
26 the first application of spectrum calculates the line intensities within the specified wavelength range and for set of ions specified
27
28 wavelength will not be used if applied to 'spectrum' objects
29
30 wavelength IS need for 'bunch' objects - in this case, the wavelength should not extend beyond the limits of the
31 wvlRange used for the 'bunch' calculation
32
33 '''
34 if not hasattr(self, 'IonInstances'):
35 print(' must set keepIons=1 in order to keep self.IonInstances')
36 return
37
38 if type(label)!= type(0):
39 if type(label) != str:
40 print(' label must either be zero or a string')
41 return
42
43 t1 = datetime.now()
44
45 if hasattr(self, 'Wavelength'):
46 nWvl = len(self.Wavelength)
47 wavelength = self.Wavelength
48 elif type(wavelength) == int:
49 print(' a wavelength array must be given')
50 return
51 else:
52 self.Wavelength = wavelength
53 nWvl = len(wavelength)
54 lineSpectrum = np.zeros((self.NTempDen, nWvl), 'float64').squeeze()
55 for akey in sorted(self.IonInstances.keys()):
56 if verbose:
57 print( ' trying ion = %s'%(akey))
58
59 if not 'errorMessage' in sorted(self.IonInstances[akey].Intensity.keys()):
60 if verbose:
61 print(' doing convolve on ion %s '%(akey))
62 self.IonInstances[akey].spectrum(wavelength, filter)
63
64 if 'errorMessage' in sorted(self.IonInstances[akey].Spectrum.keys()):
65 print(self.IonInstances[akey].Spectrum['errorMessage'])
66 else:
67 lineSpectrum += self.IonInstances[akey].Spectrum['intensity']
68
69
70
71
72
73 else:
74 if 'errorMessage' in sorted(self.IonInstances[akey].Intensity.keys()):
75 print(self.IonInstances[akey].Intensity['errorMessage'])
76
77 self.LineSpectrum = {'wavelength':wavelength, 'intensity':lineSpectrum.squeeze()}
78
79 total = self.LineSpectrum['intensity']
80
81
82
83 if hasattr(self, 'FreeFree'):
84 total += self.FreeFree['intensity']
85 if hasattr(self, 'FreeBound'):
86 total += self.FreeBound['intensity']
87 if hasattr(self, 'TwoPhoton'):
88 total += self.TwoPhoton['intensity']
89 self.Total = total
90
91
92 if self.NTempDen == 1:
93 integrated = total
94 else:
95 integrated = total.sum(axis=0)
96
97 t2 = datetime.now()
98 dt = t2 - t1
99 print(' elapsed seconds = %12.3e'%(dt.seconds))
100
101 if type(label) == type(''):
102 if hasattr(self, 'Spectrum'):
103 self.Spectrum[label] = {'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'em':self.Em, 'Abundance':self.AbundanceName}
104 else:
105 self.Spectrum = {label:{'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'em':self.Em, 'Abundance':self.AbundanceName}}
106 else:
107 self.Spectrum ={'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'Abundance':self.AbundanceName}
108 return
109
110
111
112 - def ionGate(self, elementList = 0, ionList = 0, minAbund=0, doContinuum=1, doWvlTest=1, verbose=0):
113 '''
114 creates a list of ions for free-free, free-bound, and line intensity calculations
115 if doing the radiative losses, accept all wavelength -> doWvlTest=0
116 '''
117
118 masterlist = chdata.MasterList
119 abundAll = self.AbundAll
120
121 nonzed = abundAll > 0.
122 minAbundAll = abundAll[nonzed].min()
123 if minAbund:
124 if minAbund < minAbundAll:
125 minAbund = minAbundAll
126 ionInfo = chio.masterListInfo()
127
128 if hasattr(self, 'Wavelength'):
129 wvlRange = [self.Wavelength.min(), self.Wavelength.max()]
130 elif hasattr(self, 'WvlRange'):
131 wvlRange = self.WvlRange
132 else:
133 print(' need a wavelength range in ionGate ')
134
135 temperature = self.Temperature
136
137
138 self.Todo = {}
139
140 if minAbund:
141 if doContinuum:
142 for iz in range(1, 31):
143 abundance = chdata.Abundance[self.AbundanceName]['abundance'][iz-1]
144 if abundance >= minAbund:
145 if verbose:
146 print(' %5i %5s abundance = %10.2e '%(iz, const.El[iz-1], abundance))
147
148 for ionstage in range(1, iz+2):
149 ionS = util.zion2name(iz, ionstage)
150 masterListTest = ionS in masterlist
151 masterListInfoTest = ionS in sorted(ionInfo.keys())
152 if masterListTest or masterListInfoTest:
153 if doWvlTest:
154 wvlTestMin = wvlRange[0] <= ionInfo[ionS]['wmax']
155 wvlTestMax = wvlRange[1] >= ionInfo[ionS]['wmin']
156 else:
157 wvlTestMin = 1
158 wvlTestMax = 1
159 ioneqTest = (temperature.max() >= ionInfo[ionS]['tmin']) and (temperature.min() <= ionInfo[ionS]['tmax'])
160
161 ionstageTest = ionstage > 1
162 if ionstageTest and ioneqTest and doContinuum:
163
164 if verbose:
165 print(' setting up continuum calculation for %s '%(ionS))
166 if ionS in sorted(self.Todo.keys()):
167 self.Todo[ionS] += '_ff_fb'
168 else:
169 self.Todo[ionS] = 'ff_fb'
170 if verbose:
171 print(' for ion %s do : %s'%(ionS, self.Todo[ionS]))
172
173 if elementList:
174 for i, one in enumerate(elementList):
175 elementList[i] = one.lower()
176 for one in masterlist:
177 stuff = util.convertName(one)
178 bare = stuff['Z'] == stuff['Ion']
179 if stuff['Element'] in elementList:
180 self.Todo[one] = 'line'
181 if doContinuum and not stuff['Dielectronic']:
182 self.Todo[one]+= '_ff'
183 if not bare:
184 self.Todo[one] += '_fb'
185 if ionList:
186 for one in ionList:
187 stuff = util.convertName(one)
188 bare = stuff['Z'] == stuff['Ion']
189 if masterlist.count(one):
190 self.Todo[one] = 'line'
191 if doContinuum and not stuff['Dielectronic']:
192 self.Todo[one]+= '_ff'
193 if not bare:
194 self.Todo[one] += '_fb'
195 else:
196 if verbose:
197 pstring = ' %s not in CHIANTI database'%(one)
198 print(pstring)
199
200
201
202 if minAbund:
203 for iz in range(1, 31):
204 abundance = chdata.Abundance[self.AbundanceName]['abundance'][iz-1]
205 if abundance >= minAbund:
206 if verbose:
207 print(' %5i %5s abundance = %10.2e '%(iz, const.El[iz-1], abundance))
208
209 for ionstage in range(1, iz+2):
210 ionS = util.zion2name(iz, ionstage)
211 masterListTest = ionS in masterlist
212 masterListInfoTest = ionS in sorted(ionInfo.keys())
213 if masterListTest or masterListInfoTest:
214 if masterListTest or masterListInfoTest:
215 if doWvlTest:
216 wvlTestMin = wvlRange[0] <= ionInfo[ionS]['wmax']
217 wvlTestMax = wvlRange[1] >= ionInfo[ionS]['wmin']
218 else:
219 wvlTestMin = 1
220 wvlTestMax = 1
221 ioneqTest = (temperature.max() >= ionInfo[ionS]['tmin']) and (temperature.min() <= ionInfo[ionS]['tmax'])
222
223 ionSd = util.zion2name(iz, ionstage, dielectronic=1)
224 masterListTestD = ionSd in masterlist
225 masterListInfoTestD = ionSd in sorted(ionInfo.keys())
226 if masterListTestD or masterListInfoTestD:
227 if doWvlTest:
228 wvlTestMinD = wvlRange[0] <= ionInfo[ionSd]['wmax']
229 wvlTestMaxD = wvlRange[1] >= ionInfo[ionSd]['wmin']
230 else:
231 wvlTestMinD = 1
232 wvlTestMaxD = 1
233 ioneqTestD = (temperature.max() >= ionInfo[ionSd]['tmin']) and (temperature.min() <=ionInfo[ionSd]['tmax'])
234
235 if masterListTest and wvlTestMin and wvlTestMax and ioneqTest:
236
237
238 if ionS in sorted(self.Todo.keys()):
239 self.Todo[ionS] += '_line'
240 else:
241 self.Todo[ionS] = 'line'
242
243 if verbose:
244 print(' for ion %s do : %s'%(ionS, self.Todo[ionS]))
245 if masterListTestD and wvlTestMinD and wvlTestMaxD and ioneqTestD:
246
247
248 if ionSd in sorted(self.Todo.keys()):
249 self.Todo[ionSd] += '_line'
250 else:
251 self.Todo[ionSd] = 'line'
252 if verbose:
253 print(' for ion %s do : %s'%(ionSd, self.Todo[ionSd]))
254 return
255
256
257
258 - def spectrumPlot(self, index=-1, integrated=0, saveFile=0, linLog = 'lin'):
259 '''
260 to plot the spectrum as a function of wavelength
261 '''
262 pl.figure()
263 mask = self.Em > 1.
264 if mask.sum() == 0:
265 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ ($\int\,$ N$_e\,$N$_H\,$d${\it l}$)$^{-1}$'
266 else:
267 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$'
268
269 xlabel = 'Wavelength ('+self.Defaults['wavelength'] +')'
270
271
272
273 pl.ion()
274
275 if integrated:
276 if 'wavelength' in sorted(self.Spectrum.keys()):
277 pl.plot(self.Spectrum['wavelength'], self.Spectrum['integrated'])
278 elif 'wvl' in sorted(self.Spectrum.keys()):
279 pl.plot(self.Spectrum['wvl'], self.Spectrum['integrated'])
280 pl.xlabel(xlabel)
281 pl.ylabel(ylabel)
282 pl.title('integrated spectrum')
283 else:
284 nTempDen = self.NTempDen
285 if nTempDen == 1:
286
287 if 'wavelength' in sorted(self.Spectrum.keys()):
288 pl.plot(self.Spectrum['wavelength'], self.Spectrum['intensity'])
289 elif 'wvl' in sorted(self.Spectrum.keys()):
290 pl.plot(self.Spectrum['wvl'], self.Spectrum['intensity'])
291 pl.title(' Temperature = %10.2e K'%(self.Temperature))
292 else:
293 if index < 0:
294 index = nTempDen/2
295 if 'wavelength' in sorted(self.Spectrum.keys()):
296 pl.plot(self.Spectrum['wavelength'], self.Spectrum['intensity'][index])
297 elif 'wvl' in sorted(self.Spectrum.keys()):
298 pl.plot(self.Spectrum['wvl'], self.Spectrum['intensity'][index])
299 pl.title(' Temperature = %10.2e K for index = %3i'%(self.Temperature[index], index))
300 pl.xlabel(xlabel)
301 pl.ylabel(ylabel)
302 if saveFile:
303 pl.savefig(saveFile)
304
305
306
307 - def lineSpectrumPlot(self, index = 0, integrated=0, saveFile=0, linLog = 'lin'):
308 '''
309 to plot the line spectrum as a function of wavelength
310 '''
311
312
313 pl.figure()
314 mask = self.Em > 1.
315 if mask.sum() == 0:
316 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ ($\int\,$ N$_e\,$N$_H\,$d${\it l}$)$^{-1}$'
317 else:
318 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$'
319
320
321 xlabel = 'Wavelength ('+self.Defaults['wavelength'].capitalize() +')'
322
323
324
325 pl.ion()
326
327 if integrated:
328 pl.plot(self.Spectrum['wavelength'], self.Spectrum['intensity'].sum(axis=0))
329 pl.xlabel(xlabel)
330 pl.ylabel(ylabel)
331 pl.title('integrated spectrum')
332 else:
333 nTempDen = self.NTempDen
334 if nTempDen == 1:
335
336 pl.plot(self.LineSpectrum['wavelength'], self.LineSpectrum['intensity'])
337 pl.title(' Temperature = %10.2e K '%(self.Temperature))
338 else:
339 pl.plot(self.LineSpectrum['wavelength'], self.LineSpectrum['intensity'][index])
340 pl.title(' Temperature = %10.2e K for index = %3i'%(self.Temperature[index], index))
341 pl.xlabel(xlabel)
342 pl.ylabel(ylabel)
343 if saveFile:
344 pl.savefig(saveFile)
345
346
347