| Home | Trees | Indices | Help |
|---|
|
|
1 from datetime import datetime
2 import numpy as np
3 import pylab as pl
4 np.seterr(over='ignore')
5 #import chianti.core as ch
6 import chianti
7 import chianti.data as chdata
8 import chianti.constants as const
9 import chianti.util as util
10 import chianti.Gui as chgui
11 from ._SpecTrails import _specTrails
12 #
13 defaults = chdata.Defaults
14 #chInteractive = chdata.chInteractive
15 #if chInteractive:
16 # import pylab as pl
17 #else:
18 # import matplotlib
19 # matplotlib.use('Agg')
20 # import matplotlib.pyplot as pl
21
23 '''
24 Calculate the emission spectrum as a function of temperature and density.
25
26 includes elemental abundances or ionization equilibria
27
28 temperature and density can be arrays but, unless the size of either is one (1),
29 the two must have the same size
30
31 the returned spectrum will be convolved with a filter of the specified width on the
32 specified wavelength array
33
34 the default filter is gaussianR with a resolving power of 1000. Other filters,
35 such as gaussian, box and lorentz, are available in chianti.filters. When using the box filter,
36 the width should equal the wavelength interval to keep the units of the continuum and line
37 spectrum the same.
38
39 A selection of ions can be make with ionList containing the names of
40 the desired lines in Chianti notation, i.e. C VI = c_6
41
42 a minimum abundance can be specified so that the calculation can be speeded up by excluding
43 elements with a low abundance. With solar photospheric abundances -
44
45 setting minAbund = 1.e-4 will include H, He, C, O, Ne
46 setting minAbund = 2.e-5 adds N, Mg, Si, S, Fe
47 setting minAbund = 1.e-6 adds Na, Al, Ar, Ca, Ni
48
49 Setting em will multiply the spectrum at each temperature by the value of em.
50
51 em [for emission measure], can be a float or an array of the same length as the
52 temperature/density.
53 '''
54 - def __init__(self, temperature, eDensity, elementList=0, ionList = 0, minAbund=0, doContinuum=1, abundanceName=0, verbose=0, allLines=1, keepIons=0):
55 t1 = datetime.now()
56 masterlist = chdata.MasterList
57 # use the ionList but make sure the ions are in the database
58 if ionList:
59 alist=[]
60 for one in ionList:
61 if masterlist.count(one):
62 alist.append(one)
63 else:
64 if verbose:
65 pstring = ' %s not in CHIANTI database'%(one)
66 print(pstring)
67 masterlist = alist
68 self.Defaults=defaults
69 self.Temperature = np.asarray(temperature, 'float64')
70 nTemp = self.Temperature.size
71 self.EDensity = np.asarray(eDensity, 'float64')
72 nDen = self.EDensity.size
73 nTempDen = max([nTemp, nDen])
74
75 if not abundanceName:
76 self.AbundanceName = self.Defaults['abundfile']
77 else:
78 if abundanceName in sorted(chdata.Abundance.keys()):
79 self.AbundanceName = abundanceName
80 else:
81 abundChoices = sorted(chdata.Abundance.keys())
82 # for one in wvl[topLines]:
83 # wvlChoices.append('%12.3f'%(one))
84 abundChoice = chgui.gui.selectorDialog(abundChoices,label='Select Abundance name')
85 abundChoice_idx = abundChoice.selectedIndex
86 self.AbundanceName = abundChoices[abundChoice_idx[0]]
87 abund = self.AbundanceName
88 print(' Abundance chosen: %s '%(self.AbundanceName))
89 #
90 abundAll = chdata.Abundance[self.AbundanceName]['abundance']
91 # needed by ionGate
92 self.AbundAll = abundAll
93 #
94 nonzed = abundAll > 0.
95 minAbundAll = abundAll[nonzed].min()
96 # if minAbund is even set
97 if minAbund:
98 if minAbund < minAbundAll:
99 minAbund = minAbundAll
100 # ionInfo = util.masterListInfo()
101 #
102 freeFreeLoss = np.zeros((nTempDen), 'float64').squeeze()
103 freeBoundLoss = np.zeros((nTempDen), 'float64').squeeze()
104 twoPhotonLoss = np.zeros((nTempDen), 'float64').squeeze()
105 boundBoundLoss = np.zeros((nTempDen), 'float64').squeeze()
106 twoPhotonLoss = np.zeros((nTempDen), 'float64').squeeze()
107 #
108 self.IonsCalculated = []
109 if keepIons:
110 self.IonInstances = {}
111 self.Finished = []
112 #
113 self.WvlRange = [0., 1.e+30]
114 #
115 self.ionGate(elementList = elementList, ionList = ionList, minAbund=minAbund, doContinuum=doContinuum, doWvlTest=0, verbose=verbose)
116 #
117 #
118 for akey in sorted(self.Todo.keys()):
119 zStuff = util.convertName(akey)
120 Z = zStuff['Z']
121 ionstage = zStuff['Ion']
122 dielectronic = zStuff['Dielectronic']
123 abundance = chdata.Abundance[self.AbundanceName]['abundance'][Z - 1]
124 if verbose:
125 print(' %5i %5s abundance = %10.2e '%(Z, const.El[Z-1], abundance))
126 if verbose:
127 print(' doing ion %s for the following processes %s'%(akey, self.Todo[akey]))
128 if ionstage != 1:
129 if verbose:
130 print(' calculating ff continuum for : %s'%(akey))
131 if 'ff' in self.Todo[akey]:
132 # need to skip the neutral
133 cont = chianti.core.continuum(akey, temperature, abundanceName=self.AbundanceName)
134 cont.freeFreeLoss()
135 freeFreeLoss += cont.FreeFreeLoss['rate']
136 # if nTempDen == 1:
137 # freeFree += cont.FreeFreeLoss['rate']
138 # else:
139 # for iTempDen in range(nTempDen):
140 # freeFree[iTempDen] += cont.FreeFree['rate'][iTempDen]*em[iTempDen]
141 # freeFree += cont.FreeFree['rate']
142 if 'fb' in self.Todo[akey]:
143 if verbose:
144 print(' calculating fb continuum for : %s'%(akey))
145 # try:
146 # does cont already exist - i.e did we create it for ff
147 if hasattr(cont, 'FreeFreeLoss'):
148 cont.freeBoundLoss()
149 else:
150 cont = chianti.core.continuum(akey, temperature, abundanceName=self.AbundanceName)
151 cont.freeBoundLoss()
152 if 'errorMessage' not in list(cont.FreeBoundLoss.keys()):
153 # an fblvl file exists for this ions
154 freeBoundLoss += cont.FreeBoundLoss['rate']
155 # except:
156 # cont = chianti.core.continuum(akey, temperature, abundanceName=self.AbundanceName)
157 # cont.freeBoundLoss()
158 # if 'errorMessage' not in list(cont.FreeBound.keys()):
159 # # an fblvl file exists for this ions
160 # freeBoundLoss += cont.FreeBoundLoss['rate']
161 # if nTempDen == 1:
162 # freeBoundLoss += cont.FreeBoundLoss['rate']
163 # else:
164 # for iTempDen in range(nTempDen):
165 # freeBound[iTempDen] += cont.FreeBound['rate'][iTempDen]*em[iTempDen]
166 # freeBound += cont.FreeBound['rate']
167 if 'line' in self.Todo[akey]:
168 if verbose:
169 print(' calculating spectrum for : %s'%(akey))
170 thisIon = chianti.core.ion(akey, temperature, eDensity, abundanceName=self.AbundanceName)
171 thisIon.intensity(allLines=allLines)
172 self.IonsCalculated.append(akey)
173 if 'errorMessage' not in list(thisIon.Intensity.keys()):
174 self.Finished.append(akey)
175 thisIon.boundBoundLoss()
176 # self.IonInstances[akey] = copy.deepcopy(thisIon)
177 # if setupIntensity:
178 # for bkey in self.Intensity:
179 # self.Intensity[bkey] = np.hstack((copy.copy(self.Intensity[bkey]), thisIon.Intensity[bkey]))
180 # else:
181 # setupIntensity = 1
182 # self.Intensity = thisIon.Intensity
183 boundBoundLoss += thisIon.BoundBoundLoss['rate']
184 # if nTempDen == 1:
185 # lineSpectrum += thisIon.Spectrum['intensity']
186 # else:
187 # for iTempDen in range(nTempDen):
188 # lineSpectrum[iTempDen] += thisIon.Spectrum['intensity'][iTempDen]
189 else:
190 if verbose:
191 print(thisIon.Intensity['errorMessage'])
192 # get 2 photon emission for H and He sequences
193 if (Z - ionstage) in [0, 1] and not dielectronic:
194 thisIon.twoPhotonLoss()
195 twoPhotonLoss += thisIon.TwoPhotonLoss['rate']
196
197 # for iz in range(31):
198 # abundance = chdata.Abundance['abundance'][iz-1]
199 # if abundance >= minAbund:
200 # print(' %5i %5s abundance = %10.2e '%(iz, const.El[iz-1], abundance))
201 # #
202 # for ionstage in range(1, iz+2):
203 # ionS = util.zion2name(iz, ionstage)
204 ## print ' ionS = ', ionS
205 # masterListTest = ionS in masterlist
206 # masterListInfoTest = ionS in sorted(ionInfo.keys())
207 # if masterListTest or masterListInfoTest:
208 # ioneqTest = (self.Temperature.max() >= ionInfo[ionS]['tmin']) and (self.Temperature.min() <= ionInfo[ionS]['tmax'])
209 # # construct similar test for the dielectronic files
210 # ionSd = util.zion2name(iz, ionstage, dielectronic=1)
211 # masterListTestD = ionSd in masterlist
212 # masterListInfoTestD = ionSd in sorted(ionInfo.keys())
213 # if masterListTestD or masterListInfoTestD:
214 # ioneqTestD = (self.Temperature.max() >= ionInfo[ionSd]['tmin']) and (self.Temperature.min() <=ionInfo[ionSd]['tmax'])
215 # ionstageTest = ionstage > 1
216 # if ionstageTest and ioneqTest and doContinuum:
217 # # ionS is the target ion, cannot be the neutral for the continuum
218 # print(' calculating continuum for %s'%(ionS))
219 # cont = chianti.core.continuum(ionS, temperature, abund=abund)
220 # cont.freeFreeLoss()
221 # # print dir(thisIon)
222 # # print ' wvl = ', thisIon.FreeFree['wvl']
223 ## if nTempDen ==1:
224 # freeFreeLoss += cont.FreeFreeLoss['rate']
225 ## else:
226 ## for iTempDen in range(nTempDen):
227 ## freeFreeLoss[iTempDen] += cont.FreeFreeLoss['rate'][iTempDen]
228 # #
229 # cont.freeBoundLoss()
230 # if 'errorMessage' not in sorted(cont.FreeBoundLoss.keys()):
231 # # an fblvl file exists for this ions
232 ## if nTempDen == 1:
233 # freeBoundLoss += cont.FreeBoundLoss['rate']
234 ## else:
235 ## freeBound[iTempDen] += cont.FreeBound['rate'][iTempDen]
236 # if masterListTest and ioneqTest:
237 # print(' calculating spectrum for %s '%(ionS))
238 # thisIon = chianti.core.ion(ionS, temperature, density, abund=abund)
239 ## print ' dir = ', dir(thisIon)
240 ## thisIon.emiss(wvlRange = wvlRange, allLines=allLines)
241 # thisIon.boundBoundLoss( allLines=allLines)
242 # # check that there are lines in this wavelength range
243 # if 'errorMessage' not in sorted(thisIon.BoundBoundLoss.keys()):
244 # thisIon.boundBoundLoss()
245 ## intensity = thisIon.Intensity['intensity']
246 ## if nTempDen == 1:
247 # boundBoundLoss += thisIon.BoundBoundLoss['rate']
248 ## else:
249 ## for iTempDen in range(nTempDen):
250 ## lineSpectrum[iTempDen] += thisIon.Spectrum['intensity'][iTempDen]
251 # # get 2 photon emission for H and He sequences
252 # if (iz - ionstage) in [0, 1]:
253 # thisIon.twoPhotonLoss()
254 # twoPhotonLoss += thisIon.TwoPhotonLoss['rate']
255 # # get dielectronic lines
256 # if masterListTestD and ioneqTestD:
257 # print(' calculating spectrum for %s '%(ionSd))
258 # thisIon = chianti.core.ion(ionSd, temperature, density, abund=abund)
259 ## print ' dir = ', dir(thisIon)
260 ## have to do all lines for the dielectronic satellites
261 ## thisIon.emiss(allLines=1)
262 # thisIon.intensity(allLines=allLines)
263 # # check that there are lines in this wavelength range - probably not redundant
264 # if 'errorMessage' not in sorted(thisIon.Intensity.keys()):
265 # thisIon.boundBoundLoss()
266 ## if nTempDen == 1:
267 # boundBoundLoss += thisIon.BoundBoundLoss['rate']
268 ## else:
269 ## for iTempDen in range(nTempDen):
270 ## lineSpectrum[iTempDen] += thisIon.Spectrum['intensity'][iTempDen]
271 self.FreeFreeLoss = freeFreeLoss
272 self.FreeBoundLoss = freeBoundLoss
273 self.BoundBoundLoss = boundBoundLoss
274 self.TwoPhotonLoss = twoPhotonLoss
275 #
276 total = freeFreeLoss + freeBoundLoss + boundBoundLoss + twoPhotonLoss
277 t2 = datetime.now()
278 dt=t2-t1
279 print(' elapsed seconds = %10.2e'%(dt.seconds))
280 self.RadLoss = {'rate':total, 'temperature':self.Temperature, 'density':self.EDensity, 'minAbund':minAbund, 'abundance':self.AbundanceName}
281 #
282 # -------------------------------------------------------------------
283 #
285 '''
286 to plot the radiative losses vs temperature
287 '''
288 fontsize = 16
289 temp = self.RadLoss['temperature']
290 rate = self.RadLoss['rate']
291 pl.loglog(temp, rate)
292 pl.xlabel('Temperature (K)',fontsize=fontsize)
293 pl.ylabel(r'erg s$^{-1}$ ($\int\,$ N$_e\,$N$_H\,$d${\it l}$)$^{-1}$',fontsize=fontsize)
294 # title = 'Radiative loss rate, minAbund = %6.1e'%(self.MinAbund)
295 # if self.Density.size == 1:
296 # title += ', density = %6.1e'%(self.Density)
297 # pl.title(title, fontsize=fontsize)
298
| Home | Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0.1 on Sun Feb 7 12:06:09 2016 | http://epydoc.sourceforge.net |