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