chianti.core.Spectrum

1 import copy 2 from datetime import datetime 3 import numpy as np 4 import chianti 5 import chianti.data as chdata 6 import chianti.constants as const 7 import chianti.filters as chfilters 8 import chianti.util as util 9 import chianti.io as chio 10 # 11 import chianti.Gui as chGui 12 # 13 from ._IonTrails import _ionTrails 14 from ._SpecTrails import _specTrails 15 16 defaults = chdata.Defaults 17

18 -class spectrum(_ionTrails, _specTrails):

19 ''' 20 Calculate the emission spectrum as a function of temperature and density. 21 22 one of the convenient things is that all of the instantiated ion classes, determined through such keywords as 'elementList', 23 'ionList', and 'minAbund' are kept in a dictionary self.IonInstances where self.IonInstances['mg_7'] is the class instance of 24 chianti.core.ion for 'mg_7'. All its methods and attributes are available. 25 26 includes elemental abundances and ionization equilibria 27 28 the set of abundances, a file in $XUVTOP/abundance, can be set with the keyword argument 'abundanceName' 29 30 temperature and density can be arrays but, unless the size of either is one (1), 31 the two must have the same size 32 33 the returned spectrum will be convolved with a filter of the specified width on the 34 specified wavelength array 35 36 the default filter is gaussianR with a resolving power of 1000. Other filters, 37 such as gaussian, box and lorentz, are available in chianti.filters. When using the box filter, 38 the width should equal the wavelength interval to keep the units of the continuum and line 39 spectrum the same. 40 41 Inherited methods include 'intensityList', 'intensityRatio' (between lines of different ions), 'intensityRatioSave' 42 and 'convolve' 43 44 A selection of elements can be make with elementList a list containing the names of elements 45 that are desired to be included, e.g., ['fe','ni'] 46 47 A selection of ions can be make with ionList containing the names of 48 the desired lines in Chianti notation, i.e. C VI = c_6 49 50 Both elementList and ionList can not be specified at the same time 51 52 a minimum abundance can be specified so that the calculation can be speeded up by excluding 53 elements with a low abundance. With solar photospheric abundances - 54 55 setting minAbund = 1.e-4 will include H, He, C, O, Ne 56 setting minAbund = 2.e-5 adds N, Mg, Si, S, Fe 57 setting minAbund = 1.e-6 adds Na, Al, Ar, Ca, Ni 58 59 Setting doContinuum =0 will skip the continuum calculation. 60 61 Setting em will multiply the spectrum at each temperature by the value of em. 62 63 em [for emission measure], can be a float or an array of the same length as the 64 temperature/density 65 '''

66 - def __init__(self, temperature, eDensity, wavelength, filter=(chfilters.gaussianR, 1000.), label=0, elementList = 0, ionList = 0, minAbund=0, doContinuum=1, em=0, keepIons=0, abundanceName=0, verbose=0, allLines=1):

67 # 68 t1 = datetime.now() 69 # creates Intensity dict from first ion calculated 70 setupIntensity = 0 71 # 72 masterlist = chdata.MasterList 73 # use the ionList but make sure the ions are in the database 74 #if elementList: 75 #for i, one in enumerate(elementList): 76 #elementList[i] = one.lower() 77 #alist = [] 78 #for one in masterlist: 79 #stuff = util.convertName(one) 80 #if stuff['Element'] in elementList: 81 #alist.append(one) 82 #masterlist = alist 83 #elif ionList: 84 #alist=[] 85 #for one in ionList: 86 #if masterlist.count(one): 87 #alist.append(one) 88 #else: 89 #if verbose: 90 #pstring = ' %s not in CHIANTI database'%(one) 91 #print(pstring) 92 #masterlist = alist 93 self.Defaults=defaults 94 self.Temperature = np.asarray(temperature, 'float64') 95 nTemp = self.Temperature.size 96 self.EDensity = np.asarray(eDensity, 'float64') 97 nDen = self.EDensity.size 98 self.NTempDen = max([nTemp, nDen]) 99 nTempDen = self.NTempDen 100 self.Wavelength = wavelength 101 # 102 if type(em) == int and em == 0: 103 em = np.ones(self.NTempDen, 'float64') 104 elif type(em) == float and em > 0.: 105 em = np.ones(self.NTempDen, 'float64')*em 106 elif type(em) == list or type(em) == tuple: 107 em = np.asarray(em, 'float64') 108 self.Em = em 109 # 110 if self.Em.any() > 0.: 111 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ $' 112 else: 113 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ ($\int\,$ N$_e\,$N$_H\,$d${\it l}$)$^{-1}$' 114 # 115 xlabel = 'Wavelength ('+self.Defaults['wavelength'] +')' 116 # 117 #self.AbundanceName = defaults['abundfile'] 118 #self.AbundanceAll = chdata.AbundanceAll 119 # 120 if abundanceName: 121 if abundanceName in list(chdata.Abundance.keys()): 122 self.AbundanceName = abundanceName 123 else: 124 abundChoices = list(chdata.Abundance.keys()) 125 # for one in wvl[topLines]: 126 # wvlChoices.append('%12.3f'%(one)) 127 abundChoice = chGui.gui.selectorDialog(abundChoices,label='Select Abundance name') 128 abundChoice_idx = abundChoice.selectedIndex 129 self.AbundanceName = abundChoices[abundChoice_idx[0]] 130 abundanceName = self.AbundanceName 131 print((' Abundance chosen: %s '%(self.AbundanceName))) 132 else: 133 self.AbundanceName = self.Defaults['abundfile'] 134 # 135 abundAll = chdata.Abundance[self.AbundanceName]['abundance'] 136 # needed by ionGate 137 self.AbundAll = abundAll 138 # 139 nonzed = abundAll > 0. 140 minAbundAll = abundAll[nonzed].min() 141 # if minAbund is even set 142 if minAbund: 143 if minAbund < minAbundAll: 144 minAbund = minAbundAll 145 else: 146 minAbund = minAbundAll 147 self.MinAbund = minAbund 148 # ionInfo = chio.masterListInfo() 149 wavelength = np.asarray(wavelength) 150 nWvl = wavelength.size 151 self.Wavelength = wavelength 152 wvlRange = [wavelength.min(), wavelength.max()] 153 # 154 freeFree = np.zeros((nTempDen, nWvl), 'float64').squeeze() 155 freeBound = np.zeros((nTempDen, nWvl), 'float64').squeeze() 156 twoPhoton = np.zeros((nTempDen, nWvl), 'float64').squeeze() 157 lineSpectrum = np.zeros((nTempDen, nWvl), 'float64').squeeze() 158 # 159 self.IonsCalculated = [] 160 if keepIons: 161 self.IonInstances = {} 162 self.Finished = [] 163 # 164 self.ionGate(elementList = elementList, ionList = ionList, minAbund=minAbund, doContinuum=doContinuum, verbose = verbose) 165 # 166 for akey in sorted(self.Todo.keys()): 167 zStuff = util.convertName(akey) 168 Z = zStuff['Z'] 169 ionstage = zStuff['Ion'] 170 dielectronic = zStuff['Dielectronic'] 171 abundance = chdata.Abundance[self.AbundanceName]['abundance'][Z - 1] 172 if verbose: 173 print(' %5i %5s abundance = %10.2e '%(Z, const.El[Z-1], abundance)) 174 if verbose: 175 print(' doing ion %s for the following processes %s'%(akey, self.Todo[akey])) 176 if 'ff' in self.Todo[akey]: 177 if verbose: 178 print(' calculating ff continuum for : %s'%(akey)) 179 FF = chianti.core.continuum(akey, temperature, abundanceName=self.AbundanceName, em=em) 180 FF.freeFree(wavelength) 181 # if nTempDen == 1: 182 # freeFree += cont.FreeFree['rate']*em[0] 183 # else: 184 # for iTempDen in range(nTempDen): 185 # freeFree[iTempDen] += cont.FreeFree['rate'][iTempDen]*em[iTempDen] 186 freeFree += FF.FreeFree['rate'] 187 if 'fb' in self.Todo[akey]: 188 if verbose: 189 print(' calculating fb continuum for : %s'%(akey)) 190 try: 191 FB = chianti.core.continuum(akey, temperature, abundanceName=self.AbundanceName, em=em) 192 FB.freeBound(wavelength) 193 if 'errorMessage' not in list(cont.FreeBound.keys()): 194 # an fblvl file exists for this ions 195 freeBound += FB.FreeBound['rate'] 196 except: 197 cont = chianti.core.continuum(akey, temperature, abundanceName=self.AbundanceName) 198 cont.freeBound(wavelength) 199 # if 'errorMessage' not in list(cont.FreeBound.keys()): 200 # # an fblvl file exists for this ions 201 # if nTempDen == 1: 202 # freeBound += cont.FreeBound['rate']*em[0] 203 # else: 204 # for iTempDen in range(nTempDen): 205 # freeBound[iTempDen] += cont.FreeBound['rate'][iTempDen]*em[iTempDen] 206 # freeBound += cont.FreeBound['rate'] 207 if 'line' in self.Todo[akey]: 208 if verbose: 209 print(' calculating spectrum for : %s'%(akey)) 210 thisIon = chianti.core.ion(akey, temperature, eDensity, abundanceName=self.AbundanceName) 211 thisIon.intensity(wvlRange=wvlRange, allLines=allLines, em=em) 212 self.IonsCalculated.append(akey) 213 if 'errorMessage' not in list(thisIon.Intensity.keys()): 214 self.Finished.append(akey) 215 thisIon.spectrum(wavelength, filter=filter, allLines=allLines, em=em) 216 if keepIons: 217 self.IonInstances[akey] = copy.deepcopy(thisIon) 218 if setupIntensity: 219 for bkey in self.Intensity: 220 self.Intensity[bkey] = np.hstack((copy.copy(self.Intensity[bkey]), thisIon.Intensity[bkey])) 221 else: 222 setupIntensity = 1 223 self.Intensity = thisIon.Intensity 224 lineSpectrum += thisIon.Spectrum['intensity'] 225 # if nTempDen == 1: 226 # lineSpectrum += thisIon.Spectrum['intensity'] 227 # else: 228 # for iTempDen in range(nTempDen): 229 # lineSpectrum[iTempDen] += thisIon.Spectrum['intensity'][iTempDen] 230 else: 231 if verbose: 232 print(thisIon.Intensity['errorMessage']) 233 # get 2 photon emission for H and He sequences 234 if (Z - ionstage) in [0, 1] and not dielectronic: 235 thisIon.twoPhoton(wavelength) 236 twoPhoton += thisIon.TwoPhoton['rate'] 237 238 self.FreeFree = {'wavelength':wavelength, 'intensity':freeFree.squeeze()} 239 self.FreeBound = {'wavelength':wavelength, 'intensity':freeBound.squeeze()} 240 self.LineSpectrum = {'wavelength':wavelength, 'intensity':lineSpectrum.squeeze()} 241 self.TwoPhoton = {'wavelength':wavelength, 'intensity':twoPhoton.squeeze()} 242 # 243 # 244 total = freeFree + freeBound + lineSpectrum + twoPhoton 245 self.Total = total 246 t2 = datetime.now() 247 dt=t2-t1 248 print(' elapsed seconds = %12.3f'%(dt.seconds)) 249 if nTempDen == 1: 250 integrated = total 251 else: 252 integrated = total.sum(axis=0) 253 # 254 if type(label) == type(''): 255 if hasattr(self, 'Spectrum'): 256 self.Spectrum[label] = {'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'em':em, 'ions':self.IonsCalculated, 257 'Abundance':self.AbundanceName, 'xlabel':xlabel, 'ylabel':ylabel} 258 else: 259 self.Spectrum = {label:{'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'em':em, 'ions':self.IonsCalculated, 260 'Abundance':self.AbundanceName, 'xlabel':xlabel, 'ylabel':ylabel}} 261 else: 262 self.Spectrum ={'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'integrated':integrated, 'ions':self.IonsCalculated, 263 'Abundance':self.AbundanceName, 'xlabel':xlabel, 'ylabel':ylabel}

264 # 265 # ----------------------------------------------------------------------- 266 # 267 # self.Spectrum ={'wavelength':wavelength, 'intensity':total.squeeze(), 'filter':filter[0].__name__, 'width':filter[1], 'ions':self.IonsCalculated, 'Abundance':self.AbundanceName} 268 # 269 # ---------------------------------------------------------------------------------------------- 270 #

271 -class bunch(_ionTrails, _specTrails):

272 ''' 273 Calculate the emission line spectrum as a function of temperature and density. 274 275 'bunch' is very similar to 'spectrum' except that continuum is not calculated and the spectrum 276 is not convolved over a filter. However, this can be done with the inherited convolve method 277 278 one of the convenient things is that all of the instantiated ion classes, determined through such keywords as 'elementList', 279 'ionList', and 'minAbund' are kept in a dictionary self.IonInstances where self.IonInstances['mg_7'] is the class instance of 280 chianti.core.ion for 'mg_7'. All its methods and attributes are available. 281 282 includes elemental abundances and ionization equilibria 283 284 the set of abundances, a file in $XUVTOP/abundance, can be set with the keyword argument 'abundanceName' 285 temperature and density can be arrays but, unless the size of either is one (1), 286 the two must have the same size 287 288 Inherited methods include 'intensityList', 'intensityRatio' (between lines of different ions), and 'intensityRatioSave' 289 and 'convolve'. 290 291 A selection of elements can be make with elementList a list containing the names of elements 292 that are desired to be included, e.g., ['fe','ni'] 293 294 A selection of ions can be make with ionList containing the names of 295 the desired lines in Chianti notation, i.e. C VI = c_6 296 297 Both elementList and ionList can not be specified at the same time 298 299 a minimum abundance can be specified so that the calculation can be speeded up by excluding 300 elements with a low abundance. With solar photospheric abundances - 301 302 setting minAbund = 1.e-4 will include H, He, C, O, Ne 303 setting minAbund = 2.e-5 adds N, Mg, Si, S, Fe 304 setting minAbund = 1.e-6 adds Na, Al, Ar, Ca, Ni 305 306 Setting em will multiply the spectrum at each temperature by the value of em. 307 308 em [for emission measure], can be a float or an array of the same length as the 309 temperature/density 310 ''' 311 # 312 # ------------------------------------------------------------------------------------ 313 #

314 - def __init__(self, temperature, eDensity, wvlRange, elementList=0, ionList=0, minAbund=0, keepIons=0, em=0, abundanceName=0, verbose=0, allLines=1):

315 # 316 t1 = datetime.now() 317 # creates Intensity dict from first ion calculated 318 setupIntensity = 0 319 # 320 self.Defaults=defaults 321 temperature = np.asarray(temperature, 'float64') 322 self.Temperature = temperature 323 eDensity = np.asarray(eDensity, 'float64') 324 self.EDensity = eDensity 325 326 # 327 self.EDensity = np.asarray(eDensity,'float64') 328 self.NEDens = self.EDensity.size 329 ndens = self.EDensity.size 330 ntemp = self.Temperature.size 331 tst1 = ndens == ntemp 332 tst1a = ndens != ntemp 333 tst2 = ntemp > 1 334 tst3 = ndens > 1 335 tst4 = ndens > 1 and ntemp > 1 336 if tst1 and ntemp == 1: 337 self.NTempDen = 1 338 elif tst1a and (tst2 or tst3) and not tst4: 339 self.NTempDen = ntemp*ndens 340 if ntemp == self.NTempDen and ndens != self.NTempDen: 341 self.EDensity = np.ones_like(self.Temperature)*self.EDensity 342 elif ndens == self.NTempDen and ntemp != self.NTempDen: 343 self.Temperature = np.ones_like(self.EDensity)*self.Temperature 344 elif tst1 and tst4: 345 self.NTempDen = ntemp 346 # 347 if type(em) == int and em == 0: 348 em = np.ones(self.NTempDen, 'float64') 349 elif type(em) == float and em > 0.: 350 em = np.ones(self.NTempDen, 'float64')*em 351 elif type(em) == list or type(em) == tuple: 352 em = np.asarray(em, 'float64') 353 self.Em = em 354 # 355 #if em != 0: 356 #if type(em) == type(float): 357 #if nTempDen > 1: 358 #em = np.ones_like(self.Temperature)*em 359 #nEm = nTempDen 360 #else: 361 #nEm = 1 362 #else: 363 #em = np.asarray(em, 'float64') 364 #nEm = em.size 365 #if nEm != nTempDen: 366 #print(' the emission measure array must be the same size as the temperature/density array') 367 #return 368 #else: 369 #nEm = 0 370 #self.NEm = nEm 371 #self.AbundanceName = defaults['abundfile'] 372 #self.AbundanceAll = chdata.AbundanceAll 373 # 374 if abundanceName: 375 if abundanceName in list(chdata.Abundance.keys()): 376 self.AbundanceName = abundanceName 377 else: 378 abundChoices = list(chdata.Abundance.keys()) 379 # for one in wvl[topLines]: 380 # wvlChoices.append('%12.3f'%(one)) 381 abundChoice = chGui.gui.selectorDialog(abundChoices,label='Select Abundance name') 382 abundChoice_idx = abundChoice.selectedIndex 383 self.AbundanceName = abundChoices[abundChoice_idx[0]] 384 abundanceName = self.AbundanceName 385 print((' Abundance chosen: %s '%(self.AbundanceName))) 386 else: 387 self.AbundanceName = self.Defaults['abundfile'] 388 # 389 abundAll = chdata.Abundance[self.AbundanceName]['abundance'] 390 # needed by ionGate 391 self.AbundAll = abundAll 392 # 393 # nonzed = abundAll > 0. 394 # minAbundAll = abundAll[nonzed].min() 395 # # if minAbund is even set 396 # if minAbund: 397 # if minAbund < minAbundAll: 398 # minAbund = minAbundAll 399 # self.minAbund = minAbund 400 ionInfo = chio.masterListInfo() 401 # # 402 # self.Intensity = {'ionS':[], 'lvl1':[], 'lvl2':[], 'wvl':np.ndarray, 'pretty1':np.ndarray, 'pretty2':np.ndarray, 'intensity':np.zeros((nTempDen, 0),'float64'), 'obs':np.ndarray } 403 self.IonsCalculated = [] 404 if keepIons: 405 self.IonInstances = {} 406 self.Finished = [] 407 # 408 # also needed by ionGate 409 self.WvlRange = np.asarray(wvlRange, 'float64') 410 # 411 self.ionGate(elementList = elementList, ionList = ionList, minAbund=minAbund, doContinuum=0, verbose = verbose) 412 # 413 for ionS in sorted(self.Todo.keys()): 414 nameStuff = util.convertName(ionS) 415 Z = nameStuff['Z'] 416 417 if verbose: 418 print(' calculating %s'%(ionS)) 419 thisIon = chianti.core.ion(ionS, temperature, eDensity, abundance=abundAll[Z-1]) 420 thisIon.intensity(wvlRange=wvlRange, allLines = allLines, em=em) 421 self.IonsCalculated.append(ionS) 422 # 423 if 'errorMessage' not in list(thisIon.Intensity.keys()): 424 self.Finished.append(ionS) 425 # thisIon.spectrum(wavelength, filter=filter) 426 if keepIons: 427 self.IonInstances[ionS] = copy.deepcopy(thisIon) 428 if setupIntensity: 429 for akey in self.Intensity: 430 self.Intensity[akey] = np.hstack((copy.copy(self.Intensity[akey]), thisIon.Intensity[akey])) 431 else: 432 setupIntensity = 1 433 # print(' creating Intensity dict from ion %s'%(ionS)) 434 self.Intensity = thisIon.Intensity 435 else: 436 if verbose: 437 print(thisIon.Intensity['errorMessage']) 438 # 439 # 440 t2 = datetime.now() 441 dt=t2-t1 442 print(' elapsed seconds = %12.3f'%(dt.seconds)) 443 return

444 # 445 # ------------------------------------------------------------------------- 446 # 447

Source Code for Module chianti.core.Spectrum