1 from datetime import datetime
2 import copy
3
4
5 import numpy as np
6 import chianti
7 import chianti.data as chdata
8 import chianti.constants as const
9 import chianti.filters as chfilters
10 import chianti.util as util
11 import chianti.io as chio
12 import chianti.Gui as chgui
13
14 from ._IonTrails import _ionTrails
15 from ._SpecTrails import _specTrails
16
17
18
19
20
21
22
23
24
25
26
27
28
29 from ipyparallel import Client
30
31
32
33
34 defaults = chdata.Defaults
35
36
37
38
39
41 '''
42 this is the multiprocessing version of spectrum for using inside an IPython Qtconsole or notebook.
43
44 be for creating an instance, it is necessary to type something like the following into a console
45 > ipcluster start --profile=notebook --n=3
46 this is the way to invoke things under the IPython 2.1 notation, a bit different in 2.0
47
48 Calculate the emission spectrum as a function of temperature and density.
49
50 temperature and density can be arrays but, unless the size of either is one (1),
51 the two must have the same size
52
53 the returned spectrum will be convolved with a filter of the specified width on the
54 specified wavelength array
55
56 the default filter is gaussianR with a resolving power of 100. Other filters,
57 such as gaussian, box and lorentz, are available in chianti.filters. When using the box filter,
58 the width should equal the wavelength interval to keep the units of the continuum and line
59 spectrum the same.
60
61 A selection of elements can be make with elementList a list containing the names of elements
62 that are desired to be included, e.g., ['fe','ni']
63
64 A selection of ions can be make with ionList containing the names of
65 the desired lines in Chianti notation, i.e. C VI = c_6
66
67 Both elementList and ionList can not be specified at the same time
68
69 a minimum abundance can be specified so that the calculation can be speeded up by excluding
70 elements with a low abundance. With solar photospheric abundances -
71
72 setting minAbund = 1.e-4 will include H, He, C, O, Ne
73 setting minAbund = 2.e-5 adds N, Mg, Si, S, Fe
74 setting minAbund = 1.e-6 adds Na, Al, Ar, Ca, Ni
75
76 Setting em will multiply the spectrum at each temperature by the value of em.
77
78 em [for emission measure], can be a float or an array of the same length as the
79 temperature/density.
80 allLines = 1 will include lines with either theoretical or observed wavelengths. allLines=0 will
81 include only those lines with observed wavelengths
82
83 proc = the number of processors to use
84 timeout - a small but non-zero value seems to be necessary
85 '''
86 - def __init__(self, temperature, eDensity, wavelength, filter=(chfilters.gaussianR, 1000.), label=0, elementList = 0, ionList = 0, minAbund=0, keepIons=0, doContinuum=0, allLines = 1, em=0, abundanceName=0, verbose=0, timeout=0.1):
87
88 t1 = datetime.now()
89
90 rcAll = Client()
91
92 lbvAll = rcAll.load_balanced_view()
93
94
95
96
97 setupIntensity = 0
98
99 masterlist = chdata.MasterList
100 self.Defaults = defaults
101
102 self.Temperature = np.asarray(temperature,'float64')
103 self.EDensity = np.asarray(eDensity,'float64')
104 self.NEDens = self.EDensity.size
105 ndens = self.EDensity.size
106 ntemp = self.Temperature.size
107 tst1 = ndens == ntemp
108 tst1a = ndens != ntemp
109 tst2 = ntemp > 1
110 tst3 = ndens > 1
111 tst4 = ndens > 1 and ntemp > 1
112 if tst1 and ntemp == 1:
113 self.NTempDen = 1
114 elif tst1a and (tst2 or tst3) and not tst4:
115 self.NTempDen = ntemp*ndens
116 if ntemp == self.NTempDen and ndens != self.NTempDen:
117 self.EDensity = np.ones_like(self.Temperature)*self.EDensity
118 elif ndens == self.NTempDen and ntemp != self.NTempDen:
119 self.Temperature = np.ones_like(self.EDensity)*self.Temperature
120 elif tst1 and tst4:
121 self.NTempDen = ntemp
122
123
124 if type(em) == int and em == 0:
125 em = np.ones(self.NTempDen, 'float64')
126 elif type(em) == float and em > 0.:
127 em = np.ones(self.NTempDen, 'float64')*em
128 elif type(em) == list or type(em) == tuple:
129 em = np.asarray(em, 'float64')
130 self.Em = em
131
132
133 if self.Em.any() > 0.:
134 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ $'
135 else:
136 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ ($\int\,$ N$_e\,$N$_H\,$d${\it l}$)$^{-1}$'
137
138 xlabel = 'Wavelength ('+self.Defaults['wavelength'] +')'
139
140 self.AllLines = allLines
141
142 if not abundanceName:
143 self.AbundanceName = self.Defaults['abundfile']
144 else:
145 if abundanceName in chdata.Abundance:
146 self.AbundanceName = abundanceName
147 else:
148 abundChoices = list(chdata.Abundance.keys())
149
150
151 abundChoice = chgui.gui.selectorDialog(abundChoices,label='Select Abundance name')
152 abundChoice_idx = abundChoice.selectedIndex
153 self.AbundanceName = abundChoices[abundChoice_idx[0]]
154 abundanceName = self.AbundanceName
155 print(' Abundance chosen: %s '%(self.AbundanceName))
156
157
158 abundAll = chdata.Abundance[self.AbundanceName]['abundance']
159 self.AbundAll = abundAll
160 self.MinAbund = minAbund
161
162
163
164
165
166
167
168
169 wavelength = np.asarray(wavelength)
170 nWvl = wavelength.size
171 self.Wavelength = wavelength
172
173
174
175 freeFree = np.zeros((self.NTempDen, nWvl), 'float64').squeeze()
176 freeBound = np.zeros((self.NTempDen, nWvl), 'float64').squeeze()
177 twoPhoton = np.zeros((self.NTempDen, nWvl), 'float64').squeeze()
178 lineSpectrum = np.zeros((self.NTempDen, nWvl), 'float64').squeeze()
179
180
181 allInpt = []
182
183 if keepIons:
184 self.IonInstances = {}
185 self.FbInstances = {}
186 self.FfInstances = {}
187
188
189
190 self.ionGate(elementList = elementList, ionList = ionList, minAbund=minAbund, doContinuum=doContinuum, verbose = verbose)
191
192 for akey in sorted(self.Todo.keys()):
193 zStuff = util.convertName(akey)
194 Z = zStuff['Z']
195 abundance = chdata.Abundance[self.AbundanceName]['abundance'][Z - 1]
196 if verbose:
197 print(' %5i %5s abundance = %10.2e '%(Z, const.El[Z-1], abundance))
198 if verbose:
199 print(' doing ion %s for the following processes %s'%(akey, self.Todo[akey]))
200 if 'ff' in self.Todo[akey]:
201
202
203 allInpt.append([akey, 'ff', temperature, wavelength, abundance, em])
204 if 'fb' in self.Todo[akey]:
205
206
207 allInpt.append([akey, 'fb', temperature, wavelength, abundance, em])
208 if 'line' in self.Todo[akey]:
209
210
211 allInpt.append([akey, 'line', temperature, eDensity, wavelength, filter, allLines, abundance, em, doContinuum])
212
213 for anInpt in allInpt:
214 lbvAll.apply(doAll, anInpt)
215 lbvAll.wait()
216 lbvAll.get_result()
217 if verbose:
218 print(' got all ff, fb, line results')
219 ionsCalculated = []
220
221 for ijk in range(len(list(lbvAll.results.values()))):
222 out = list(lbvAll.results.values())[ijk]
223 if type(out) != list:
224 print(' a problem has occured - this can be caused by')
225 print('running Python3 and not using ipcluster3')
226 return
227 ionS = out[0]
228 if verbose:
229 print(' collecting calculation for %s'%(ionS))
230 ionsCalculated.append(ionS)
231 calcType = out[1]
232 if verbose:
233 print(' processing %s results'%(calcType))
234
235 if calcType == 'ff':
236 thisFf = out[2]
237 if 'errorMessage' not in sorted(thisFf.keys()):
238 if keepIons:
239 self.FfInstances[ionS] = thisFf
240
241
242
243
244
245 freeFree += thisFf['rate']
246 elif type(thisFf) == str:
247 print(' error in FfCont %s'%(thisFf))
248 else:
249 print(thisFf['errorMessage'])
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268 elif calcType == 'fb':
269 thisFb = out[2]
270 if verbose:
271 print(' fb ion = %s'%(ionS))
272 if hasattr(thisFb, 'FreeBound'):
273 if 'errorMessage' not in sorted(thisFb.FreeBound.keys()):
274 if keepIons:
275 self.FbInstances[ionS] = thisFb
276
277
278
279
280
281 freeBound += thisFb.FreeBound['rate']
282 else:
283 print(thisFb.FreeBound['errorMessage'])
284
285 elif calcType == 'line':
286 thisIon = out[2]
287 if not 'errorMessage' in sorted(thisIon.Intensity.keys()):
288 if keepIons:
289 self.IonInstances[ionS] = thisIon
290 thisIntensity = thisIon.Intensity
291
292 if setupIntensity:
293 for akey in sorted(self.Intensity.keys()):
294 self.Intensity[akey] = np.hstack((self.Intensity[akey], thisIntensity[akey]))
295 else:
296 setupIntensity = 1
297 self.Intensity = thisIntensity
298
299 lineSpectrum += thisIon.Spectrum['intensity']
300
301
302
303
304
305
306 if len(out) == 4:
307 tp = out[3]
308 if self.NTempDen == 1:
309 twoPhoton += tp['rate']
310 else:
311 for iTempDen in range(self.NTempDen):
312 twoPhoton[iTempDen] += tp['rate'][iTempDen]
313 else:
314 if 'errorMessage' in sorted(thisIon.Intensity.keys()):
315 print(thisIon.Intensity['errorMessage'])
316
317
318 self.IonsCalculated = ionsCalculated
319
320
321 self.FreeFree = {'wavelength':wavelength, 'intensity':freeFree.squeeze()}
322 self.FreeBound = {'wavelength':wavelength, 'intensity':freeBound.squeeze()}
323 self.LineSpectrum = {'wavelength':wavelength, 'intensity':lineSpectrum.squeeze()}
324 self.TwoPhoton = {'wavelength':wavelength, 'intensity':twoPhoton.squeeze()}
325
326 total = freeFree + freeBound + lineSpectrum + twoPhoton
327
328 t2 = datetime.now()
329 dt=t2-t1
330 print(' elapsed seconds = %12.3e'%(dt.seconds))
331 rcAll.purge_results('all')
332
333 if self.NTempDen == 1:
334 integrated = total
335 else:
336 integrated = total.sum(axis=0)
337
338 if type(label) == type(''):
339 if hasattr(self, 'Spectrum'):
340 print(' hasattr = true')
341 self.Spectrum[label] = {'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'em':em, 'Abundance':self.AbundanceName,
342 'xlabel':xlabel, 'ylabel':ylabel}
343 else:
344 self.Spectrum = {label:{'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'em':em, 'Abundance':self.AbundanceName,
345 'xlabel':xlabel, 'ylabel':ylabel}}
346 else:
347 self.Spectrum = {'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'em':em, 'Abundance':self.AbundanceName,
348 'xlabel':xlabel, 'ylabel':ylabel}
349
350
351
353 '''
354 to process ff, fb and line inputs
355 '''
356 ionS = inpt[0]
357 calcType = inpt[1]
358 if calcType == 'ff':
359 temperature = inpt[2]
360 wavelength = inpt[3]
361 abund = inpt[4]
362 em = inpt[5]
363 FF = chianti.core.continuum(ionS, temperature, abundance=abund, em=em)
364 FF.freeFree(wavelength)
365
366
367 return [ionS, calcType, copy.copy(FF.FreeFree)]
368 elif calcType == 'fb':
369 temperature = inpt[2]
370 wavelength = inpt[3]
371 abund = inpt[4]
372 em = inpt[5]
373 cont = chianti.core.continuum(ionS, temperature, abundance=abund, em=em)
374 cont.freeBound(wavelength)
375 return [ionS, calcType, copy.deepcopy(cont)]
376 elif calcType == 'line':
377 temperature = inpt[2]
378 density = inpt[3]
379 wavelength = inpt[4]
380 wvlRange = [wavelength.min(), wavelength.max()]
381 filter = inpt[5]
382 allLines = inpt[6]
383 abund = inpt[7]
384 em = inpt[8]
385 doContinuum = inpt[9]
386 thisIon = chianti.core.ion(ionS, temperature, density, abundance=abund)
387 thisIon.intensity(wvlRange = wvlRange, allLines = allLines, em=em)
388 if 'errorMessage' not in thisIon.Intensity.keys():
389 thisIon.spectrum(wavelength, filter=filter, allLines=allLines)
390
391
392
393
394 outList = [ionS, calcType, copy.deepcopy(thisIon)]
395 if not thisIon.Dielectronic and doContinuum:
396 if (thisIon.Z - thisIon.Ion) in [0, 1]:
397 thisIon.twoPhoton(wavelength, em=em)
398 outList.append(thisIon.TwoPhoton)
399 return outList
400
401
402
403