chianti.core.Ion

22 ''' 23 The top level class for performing spectral calculations for an ion in the CHIANTI database. 24 25 ionStr is a string corresponding such as 'c_5' that corresponds to the C V ion. 26 temperature in Kelvin 27 eDensity in cm^-3 28 radTemperature, the radiation black-body temperature in Kelvin 29 rPlot, the distance from the center of the star in stellar radii 30 31 Inherited methods include 'intensityList', 'intensityRatio' (between lines of different ions), 'intensityRatioSave' 32 and 'convolve' 33 '''

34 - def __init__(self, ionStr, temperature=None, eDensity=None, pDensity='default', radTemperature=0, rStar=0, abundanceName=0, abundance=0, verbose=0, setup=True, em=0):

35 ''' 36 The top level class for performing spectral calculations for an ion in the CHIANTI database. 37 38 ionStr is a string corresponding such as 'c_5' that corresponds to the C V ion. 39 temperature in Kelvin 40 eDensity in cm^-3 41 radTemperature, the radiation black-body temperature in Kelvin 42 rStar, the distance from the center of the star in stellar radii 43 44 the elemental abundance values can be set to one of the names in the XUVTOP/abundance directory 45 without the '.abund' suffix, such as abundanceName = 'sun_photospheric_1998_grevesse' 46 47 ''' 48 # 49 # 50 # self.__version__ = chianti.__version__ 51 self.IonStr=ionStr 52 self.Z=util.convertName(ionStr)['Z'] 53 self.Ion=util.convertName(ionStr)['Ion'] 54 self.Dielectronic=util.convertName(ionStr)['Dielectronic'] 55 self.Spectroscopic=util.zion2spectroscopic(self.Z,self.Ion) 56 self.FileName=util.zion2filename(self.Z, self.Ion,dielectronic=self.Dielectronic ) 57 # 58 self.Defaults=chdata.Defaults 59 # 60 if abundance: 61 self.Abundance=abundance 62 elif abundanceName: 63 if abundanceName in chdata.AbundanceList: 64 self.AbundanceName = abundanceName 65 self.Abundance = chdata.Abundance[self.AbundanceName]['abundance'][self.Z-1] 66 else: 67 abundChoices = chdata.AbundanceList 68 # for one in wvl[topLines]: 69 # wvlChoices.append('%12.3f'%(one)) 70 abundChoice = chGui.gui.selectorDialog(abundChoices,label='Select Abundance name') 71 abundChoice_idx = abundChoice.selectedIndex 72 self.AbundanceName = abundChoices[abundChoice_idx[0]] 73 self.Abundance = chdata.Abundance[self.AbundanceName]['abundance'][self.Z-1] 74 print(' Abundance chosen: %s '%(self.AbundanceName)) 75 else: 76 self.AbundanceName = self.Defaults['abundfile'] 77 self.Abundance = chdata.Abundance[self.AbundanceName]['abundance'][self.Z-1] 78 # 79 # 80 self.IoneqName = self.Defaults['ioneqfile'] 81 #MasterList = chdata.MasterList 82 # 83 # 84 self.RadTemperature = radTemperature 85 self.RStar = rStar 86 # 87 # ip in eV, but don't read for bare ions 88 if self.Ion <= self.Z: 89 self.Ip=chdata.Ip[self.Z-1, self.Ion-1-self.Dielectronic] 90 self.FIP = chdata.Ip[self.Z-1, 0] 91 if self.Dielectronic: 92 self.UpperIp=chdata.Ip[self.Z-1, self.Ion-1] 93 # 94 if type(temperature) != type(None): 95 self.Temperature = np.array(temperature,'float64') 96 self.NTemp = self.Temperature.size 97 # 98 # 99 # 100 # self.IoneqAll = io.ioneqRead(ioneqname = self.Defaults['ioneqfile']) 101 # if chInteractive: 102 self.IoneqAll = chdata.IoneqAll 103 self.ioneqOne() 104 # 105 # this needs to go after setting temperature and reading ionization equilibria 106 if pDensity == 'default': 107 self.p2eRatio() 108 # 109 if type(eDensity) != type(None): 110 self.EDensity = np.asarray(eDensity,'float64') 111 self.NEDens = self.EDensity.size 112 ndens = self.EDensity.size 113 ntemp = self.Temperature.size 114 tst1 = ndens == ntemp 115 tst1a = ndens != ntemp 116 tst2 = ntemp > 1 117 tst3 = ndens > 1 118 tst4 = ndens > 1 and ntemp > 1 119 if tst1 and ntemp == 1: 120 self.NTempDen = 1 121 elif tst1a and (tst2 or tst3) and not tst4: 122 self.NTempDen = ntemp*ndens 123 if ntemp == self.NTempDen and ndens != self.NTempDen: 124 self.EDensity = np.ones_like(self.Temperature)*self.EDensity 125 elif ndens == self.NTempDen and ntemp != self.NTempDen: 126 self.Temperature = np.ones_like(self.EDensity)*self.Temperature 127 elif tst1 and tst4: 128 self.NTempDen = ntemp 129 # 130 if pDensity == 'default' and type(eDensity) != type(None): 131 if tst1 and tst2 and tst3: 132 self.PDensity = np.zeros((ntemp), 'float64') 133 for itemp in range(ntemp): 134 self.PDensity[itemp] = self.ProtonDensityRatio[itemp]*self.EDensity[itemp] 135 elif tst2 and tst3 and not tst1: 136 print(' if both temperature and eDensity are arrays, they must be of the same size') 137 return 138 else: 139 self.PDensity = self.ProtonDensityRatio*self.EDensity 140 else: 141 self.PDensity = pDensity 142 if setup: 143 if self.IonStr in chdata.MasterList: 144 self.setup() 145 else: 146 self.setupIonrec() 147 # 148 if type(em) == int and em == 0: 149 # if hasattr(self, 'Em'): 150 # em = self.Em 151 # else: 152 # em = np.ones(self.NTempDen, 'float64') 153 # self.Em = em 154 pass 155 elif type(em) == float and em > 0.: 156 em = np.ones(self.NTempDen, 'float64')*em 157 self.Em = em 158 elif type(em) == list or type(em) == tuple or type(em) == np.ndarray: 159 em = np.asarray(em, 'float64') 160 self.Em = em

161 # 162 # ------------------------------------------------------------------------------ 163 #

164 - def diCross(self, energy=None, verbose=False):

165 ''' 166 Calculate the direct ionization cross section. 167 168 Given as a function of the incident electron energy in eV, puts values into DiCross 169 ''' 170 iso = self.Z - self.Ion + 1 171 if type(energy) == type(None): 172 btenergy=0.1*np.arange(10) 173 btenergy[0]=0.01 174 dum=np.ones(len(btenergy)) 175 [energy, dum]=util.descale_bti(btenergy, dum, 2., self.Ip) 176 energy=np.asarray(energy, 'float64') 177 # 178 if iso == 1 and self.Z >= 6: 179 # hydrogenic sequence 180 ryd=27.2113845/2. 181 u=energy/self.Ip 182 ev1ryd=self.Ip/ryd 183 a0=0.5291772108e-8 184 a_bohr=const.pi*a0**2 # area of bohr orbit 185 if self.Z >= 20: 186 ff = (140.+(self.Z/20.)**3.2)/141. 187 else: 188 ff = 1. 189 qr = util.qrp(self.Z,u)*ff 190 bb = 1. # hydrogenic 191 qh = bb*a_bohr*qr/ev1ryd**2 192 self.DiCross = {'energy':energy, 'cross':qh} 193 elif iso == 2 and self.Z >= 10: 194 # use 195 ryd=27.2113845/2. 196 u=energy/self.Ip 197 ev1ryd=self.Ip/ryd 198 a0=0.5291772108e-8 199 a_bohr=const.pi*a0**2 # area of bohr orbit 200 if self.Z >= 20: 201 ff=(140.+(self.Z/20.)**3.2)/141. 202 else: 203 ff=1. 204 qr=util.qrp(self.Z,u)*ff 205 bb=2. # helium-like 206 qh=bb*a_bohr*qr/ev1ryd**2 207 self.DiCross={'energy':energy, 'cross':qh} 208 else: 209 # try: 210 # diparams = self.DiParams 211 # except: 212 # self.DiParams = io.diRead(self.IonStr) 213 # diparams = self.DiParams 214 if not hasattr(self, 'DiParams'): 215 self.DiParams = io.diRead(self.IonStr) 216 217 cross=np.zeros(len(energy), 'Float64') 218 219 for ifac in range(self.DiParams['info']['nfac']): 220 # prob. better to do this with masked arrays 221 goode=energy > self.DiParams['ev1'][ifac] 222 if goode.sum() > 0: 223 dum=np.ones(len(energy)) 224 btenergy, btdum=util.scale_bti(energy[goode],dum[goode], self.DiParams['btf'][ifac], self.DiParams['ev1'][ifac] ) 225 # these interpolations were made with the scipy routine used here 226 y2=interpolate.splrep(self.DiParams['xsplom'][ifac], self.DiParams['ysplom'][ifac], s=0) 227 btcross=interpolate.splev(btenergy, y2, der=0) 228 energy1, cross1=util.descale_bti(btenergy, btcross,self.DiParams['btf'][ifac], self.DiParams['ev1'][ifac] ) 229 offset=len(energy)-goode.sum() 230 if verbose: 231 pl.plot(self.DiParams['xsplom'][ifac], self.DiParams['ysplom'][ifac]) 232 pl.plot(btenergy, btcross) 233 if offset > 0: 234 seq=[np.zeros(offset, 'Float64'), cross1] 235 cross1=np.hstack(seq) 236 cross += cross1*1.e-14 237 self.DiCross={'energy':energy, 'cross':cross}

238 # 239 #----------------------------------------------------------- 240 #

241 - def diRate(self):

242 ''' 243 Calculate the direct ionization rate coefficient as a function of temperature (K) 244 ''' 245 # if hasattr(self, 'DiParams'): 246 # DiParams = self.DiParams 247 # else: 248 # DiParams = io.diRead(self.IonStr) 249 # 250 if hasattr(self, 'Temperature'): 251 temperature=self.Temperature 252 else: 253 print(' temperature is not defined') 254 return 255 # 256 # gauss laguerre n=12 257 # 258 ngl=12 259 xgl=np.asarray([0.115722117358021,0.611757484515131,1.512610269776419,2.833751337743509 260 ,4.599227639418353,6.844525453115181,9.621316842456871,13.006054993306350 261 ,17.116855187462260,22.151090379396983,28.487967250983992,37.099121044466926], 'float64') 262 263 264 wgl=np.asarray([2.647313710554435e-01,3.777592758731382e-01,2.440820113198774e-01,9.044922221168074e-02 265 ,2.010238115463406e-02,2.663973541865321e-03,2.032315926629993e-04,8.365055856819753e-06 266 ,1.668493876540914e-07,1.342391030515027e-09,3.061601635035012e-12,8.148077467426124e-16], 'float64') 267 # 268 alpha=5.287e+13 269 tev=const.boltzmannEv*temperature 270 # 271 ntemp=temperature.size 272 # 273 # 274 if ntemp == 1: 275 x0=self.Ip/tev # Ip in eV 276 beta=np.sqrt(const.boltzmann*temperature) 277 egl=self.Ip+xgl*tev 278 self.diCross(energy=egl) 279 crossgl=self.DiCross['cross'] 280 term1=wgl*xgl*crossgl 281 term2=wgl*crossgl 282 newcross=alpha*beta*np.exp(-x0)*(term1.sum()+x0*term2.sum()) 283 rate=newcross 284 else: 285 rate=np.zeros(ntemp, 'float64') 286 for itemp in range(ntemp): 287 x0=self.Ip/tev[itemp] # Ip in eV 288 beta=np.sqrt(const.boltzmann*temperature[itemp]) 289 egl=self.Ip+xgl*tev[itemp] 290 self.diCross(energy=egl) 291 crossgl=self.DiCross['cross'] 292 term1=wgl*xgl*crossgl 293 term2=wgl*crossgl 294 newcross=alpha*beta*np.exp(-x0)*(term1.sum()+x0*term2.sum()) 295 rate[itemp]=newcross 296 self.DiRate={'temperature':temperature, 'rate':rate}

297 # 298 #----------------------------------------------------------- 299 #

300 - def eaDescale(self):

301 """ 302 Calculates the effective collision strengths (upsilon) 303 for excitation-autoionization as a function of temperature. 304 """ 305 # 306 # xt=kt/de 307 # 308 # need to make sure elvl is >0, except for ground level 309 # 310 if hasattr(self, 'EaParams'): 311 eaparams=self.EaParams 312 else: 313 self.EaParams = io.eaRead(self.IonStr) 314 eaparams=self.EaParams 315 # 316 if hasattr(self, 'Temperature'): 317 temperature=self.Temperature 318 else: 319 print(' temperature is not defined') 320 return 321 ntemp=temperature.size 322 nsplups=len(eaparams['de']) 323 if ntemp > 1: 324 ups=np.zeros((nsplups,ntemp),"Float64") 325 else: 326 ups=np.zeros(nsplups,"Float64") 327 # 328 for isplups in range(0,nsplups): 329 l1=self.EaParams["lvl1"][isplups]-1 330 l2=self.EaParams["lvl2"][isplups]-1 331 ttype=self.EaParams["ttype"][isplups] 332 cups=self.EaParams["cups"][isplups] 333 nspl=self.EaParams["nspl"][isplups] 334 de=self.EaParams["de"][isplups] 335 dx=1./(float(nspl)-1.) 336 ## print self.EaParams["EaParams"][l1,l2] 337 splups=self.EaParams["splups"][isplups,0:nspl] 338 kte=const.boltzmannEv*temperature/(const.ryd2Ev*de) 339 # 340 if ttype ==1: 341 st=1.-np.log(cups)/np.log(kte+cups) 342 xs=dx*np.arange(nspl) 343 y2=interpolate.splrep(xs,splups,s=0) #allow smoothing,s=0) 344 sups=interpolate.splev(st,y2,der=0) 345 ups[isplups]=sups*np.log(kte+np.exp(1.)) 346 # 347 if ttype == 2: 348 st=kte/(kte+cups) 349 xs=dx*np.arange(nspl) 350 y2=interpolate.splrep(xs,splups,s=0) 351 sups=interpolate.splev(st,y2,der=0) 352 ups[isplups]=sups 353 # 354 if ttype == 3: 355 st=kte/(kte+cups) 356 xs=dx*np.arange(nspl) 357 y2=interpolate.splrep(xs,splups,s=0) 358 sups=interpolate.splev(st,y2,der=0) 359 ups[isplups]=sups/(kte+1.) 360 # 361 if ttype == 4: 362 st=1.-np.log(cups)/np.log(kte+cups) 363 xs=dx*np.arange(nspl) 364 y2=interpolate.splrep(xs,splups,s=0) 365 sups=interpolate.splev(st,y2,der=0) 366 ups[isplups]=sups*np.log(kte+cups) 367 # 368 if ttype == 5: 369 st=kte/(kte+cups) 370 xs=dx*np.arange(nspl) 371 y2=interpolate.splrep(xs,splups,s=0) #allow smoothing,s=0) 372 sups=interpolate.splev(st,y2,der=0) 373 ups[isplups]=sups/(kte+0.) 374 # 375 # 376 elif ttype > 5: print(' t_type ne 1,2,3,4,5 = %5i %5i %5i'%(ttype,l1,l2)) 377 # 378 # 379 ups=np.where(ups > 0.,ups,0.) 380 # 381 self.EaParams['ups']=ups 382 return ups

383 # 384 # ------------------------------------------------------------------------------------- 385 #

386 - def eaCross(self, energy=0, verbose=False):

387 ''' 388 Provide the excitation-autoionization cross section. 389 390 Energy is given in eV. 391 ''' 392 # get neaev from diparams file 393 # 394 # try: 395 # diparams = self.DiParams 396 # except: 397 # self.DiParams = io.diRead(self.IonStr) 398 # diparams = self.DiParams 399 # 400 if self.DiParams['info']['neaev'] == 0: 401 # print ' no EA rates' 402 return 403 else: 404 if not energy.all(): 405 energy=self.Ip*10.**(0.05*np.arange(31)) 406 try: 407 easplom=self.Easplom 408 except: 409 # self.splomRead() 410 self.Easplom = util.splomRead(self.IonStr, ea=1) 411 easplom =self.Easplom 412 # 413 # multiplicity of ground level already included 414 # 415 # splomDescale takes care of when energy < threshold 416 # 417 omega = util.splomDescale(easplom, energy) 418 # 419 # need to replicate neaev 420 ntrans=len(easplom['deryd']) 421 nsplom=easplom['splom'].shape[1] 422 x=0.25*np.arange(nsplom) 423 eaev=self.DiParams['eaev'] 424 if len(eaev) ==1: 425 for itrans in range(ntrans): 426 eaev.append(eaev[0]) 427 428 totalCross = np.zeros_like(energy) 429 ntrans = omega.shape[0] 430 for itrans in range(ntrans): 431 lvl1 = self.Easplom['lvl1'][itrans] 432 # the collision strengths have already by divided by the 433 # statistical weight of the ground level 2j+1 434 cross = eaev[itrans]*const.bohrCross*omega[itrans]/(energy/const.ryd2Ev) 435 totalCross += cross 436 self.EaCross = {'energy':energy, 'cross':totalCross} 437 return

438 # 439 # ------------------------------------------------------------------------------------- 440 #

441 - def eaRate(self):

442 '''Calculate the excitation-autoionization rate coefficient.''' 443 # get neaev from diparams file 444 # 445 if not hasattr(self, 'DiParams'): 446 self.DiParams = io.diRead(self.IonStr) 447 # 448 if self.DiParams['info']['neaev'] == 0: 449 # print ' no EA rates' 450 return 451 else: 452 if hasattr(self, 'Temperature'): 453 temperature=self.Temperature 454 else: 455 bte=0.1*np.arange(10) 456 bte[0]=0.01 457 dum=np.ones(10, 'Float64') 458 [temperature, dum]=util.descale_bt(bte, dum, self.EaParams['cups'][0], self.DiParams['de'][0]) 459 self.Temperature=temperature 460 if hasattr(self, 'EaParams'): 461 eaparams=self.EaParams 462 else: 463 self.eaParams = io.eaRead(self.IonStr) 464 self.eaDescale() 465 eaparams=self.EaParams 466 # 467 # need to replicate neaev 468 nups=len(eaparams['de']) 469 tev=const.boltzmannEv*temperature 470 earate=np.zeros(temperature.size, 'Float64') 471 eaev=self.DiParams['eaev'] 472 if len(eaev) ==1: 473 for iups in range(nups): 474 eaev.append(eaev[0]) 475 476 for iups in range(nups): 477 x0=const.ryd2Ev*eaparams['de'][iups]/tev 478 # upsilon has already been divided by the statistical weight of the 479 # ground level 2j+1 480 earate+=eaev[iups]*8.63e-6*eaparams['ups'][iups]*np.exp(-x0)/(np.sqrt(temperature)) 481 self.EaRate={'rate':earate, 'temperature':temperature} 482 return

483 # 484 # ------------------------------------------------------------------------------------- 485 #

486 - def ionizCross(self, energy=0):

487 ''' 488 Provides the total ionization cross section. 489 490 uses diCross and eaCross. 491 ''' 492 if self.Z < self.Ion: 493 # print ' this is a bare nucleus and has no ionization rate' 494 self.IonizRate = {'rate':np.zeros_like(self.Temperature), 'temperature':self.Temperature} 495 return 496 if not energy.all(): 497 energy = self.Ip*10.**(0.025*np.arange(401)) 498 # 499 self.diCross(energy) 500 self.eaCross(energy) 501 if self.DiParams['info']['neaev'] == 0: 502 ionizCross=self.DiCross['cross'] 503 else: 504 ionizCross=self.DiCross['cross']+self.EaCross['cross'] 505 self.IonizCross = {'cross':ionizCross, 'energy':energy} 506 return

507 # 508 # ------------------------------------------------------------------------------------- 509 #

510 - def ionizRate(self):

511 ''' 512 Provides the total ionization rate. 513 514 Calls diRate and eaRate. 515 ''' 516 if self.Z < self.Ion: 517 # print ' this is a bare nucleus and has no ionization rate' 518 self.IonizRate = {'rate':np.zeros_like(self.Temperature), 'temperature':self.Temperature} 519 return 520 self.diRate() 521 self.eaRate() 522 if self.DiParams['info']['neaev'] == 0: 523 ionizrate=self.DiRate['rate'] 524 else: 525 ionizrate=self.DiRate['rate']+self.EaRate['rate'] 526 self.IonizRate = {'rate':ionizrate, 'temperature':self.DiRate['temperature']} 527 return

528 # 529 # ------------------------------------------------------------------------------------- 530 # 531 # 532 # ------------------------------------------------------------------------------------- 533 #

534 - def rrRate(self):

535 '''Provide the radiative recombination rate coefficient as a function of temperature (K).''' 536 if hasattr(self, 'Temperature'): 537 temperature=self.Temperature 538 else: 539 print(' temperature is not defined') 540 return 541 rrparamsfile = util.ion2filename(self.IonStr) + '.rrparams' 542 if hasattr(self, 'RrParams'): 543 rrparams=self.RrParams 544 elif os.path.isfile(rrparamsfile): 545 self.RrParams = io.rrRead(self.IonStr) 546 rrparams=self.RrParams 547 else: 548 self.RrRate={'temperature':temperature, 'rate':np.zeros_like(temperature)} 549 return 550 # 551 # print ' rr params type = ', rrparams['rrtype'] 552 # 553 if rrparams['rrtype'] == 1: 554 a=rrparams['params'][3] 555 b=rrparams['params'][4] 556 t0=rrparams['params'][5] 557 t1=rrparams['params'][6] 558 rate=a/(np.sqrt(temperature/t0)*(1.+np.sqrt(temperature/t0))**(1.-b)*(1.+np.sqrt(temperature/t1))**(1.+b)) 559 self.RrRate={'temperature':temperature, 'rate':rate} 560 elif rrparams['rrtype'] == 2: 561 a=rrparams['params'][3] 562 b=rrparams['params'][4] 563 t0=rrparams['params'][5] 564 t1=rrparams['params'][6] 565 c=rrparams['params'][7] 566 t2=rrparams['params'][8] 567 b += c*np.exp(-t2/temperature) 568 rate=a/(np.sqrt(temperature/t0)*(1.+np.sqrt(temperature/t0))**(1.-b)*(1.+np.sqrt(temperature/t1))**(1.+b)) 569 self.RrRate={'temperature':temperature, 'rate':rate} 570 elif rrparams['rrtype'] == 3: 571 a=rrparams['params'][2] 572 b=rrparams['params'][3] 573 rate=a/(temperature/1.e+4)**b 574 self.RrRate={'temperature':temperature, 'rate':rate} 575 else: 576 self.RrRate={'temperature':temperature, 'rate':np.zeros_like(temperature)}

577 578 # 579 # ------------------------------------------------------------------------------------- 580 #

581 - def drRate(self):

582 '''Provide the dielectronic recombination rate coefficient as a function of temperature (K). 583 ''' 584 # 585 if hasattr(self, 'Temperature'): 586 temperature=self.Temperature 587 else: 588 print(' temperature is not defined') 589 return {'errorMessage':' temperature is not defined'} 590 drparamsfile = util.ion2filename(self.IonStr) + '.drparams' 591 if hasattr(self, 'DrParams'): 592 drparams=self.DrParams 593 elif os.path.isfile(drparamsfile): 594 self.DrParams = io.drRead(self.IonStr) 595 drparams=self.DrParams 596 else: 597 self.DrRate = {'rate':np.zeros_like(temperature), 'temperature':temperature} 598 return 599 # 600 if drparams['drtype'] == 1: 601 # badnell type 602 drenergy=drparams['eparams'] 603 drcoef=drparams['cparams'] 604 gcoef = drenergy > 0. 605 ncoef=gcoef.sum() 606 # print ' ncoef = ', gcoef.sum() 607 rate=np.zeros(temperature.size, 'float64') 608 for icoef in range(ncoef): 609 rate += drcoef[icoef]*np.exp(-drenergy[icoef]/temperature) 610 rate=rate/temperature**1.5 611 self.DrRate={'temperature':temperature, 'rate':rate} 612 elif drparams['drtype'] == 2: 613 # shull type 614 params = drparams['params'] 615 adi = params[0] 616 bdi = params[1] 617 t0 = params[2] 618 t1 = params[3] 619 rate=adi*np.exp(-t0/temperature)*(1.+bdi*np.exp(-t1/temperature))/temperature**1.5 620 self.DrRate={'temperature':temperature, 'rate':rate}

621 # 622 # ------------------------------------------------------------------------------------- 623 #

624 - def cireclvlDescale(self, lvlType):

625 '''Interpolate and extrapolate cilvl and reclvl rates. 626 lvltype must be either 'reclvl', 'cilvl' or 'rrlvl' 627 Used in level population calculations. 628 ''' 629 if hasattr(self, 'Temperature'): 630 temperature=self.Temperature 631 else: 632 print(' temperature is not defined') 633 return {'errorMessage':' temperature is not defined'} 634 lvlfile = util.ion2filename(self.IonStr)+'.' + lvlType 635 if lvlType == 'reclvl': 636 if hasattr(self, 'Reclvl'): 637 lvl = self.Reclvl 638 elif os.path.isfile(lvlfile): 639 # print ' reading reclvl file' 640 self.Reclvl = io.cireclvlRead(self.IonStr, '.'+lvlType) 641 lvl = self.Reclvl 642 else: 643 self.ReclvlRate = {'rate':np.zeros_like(temperature)} 644 return 645 elif lvlType == 'cilvl': 646 if hasattr(self, 'Cilvl'): 647 lvl = self.Cilvl 648 elif os.path.isfile(lvlfile): 649 # print ' reading reclvl file' 650 self.Cilvl = io.cireclvlRead(self.IonStr, '.'+lvlType) 651 lvl = self.Cilvl 652 else: 653 # self.CilvlRate = {'rate':np.zeros_like(temperature)} 654 return 655 elif lvlType == 'rrlvl': 656 if hasattr(self, 'Rrlvl'): 657 lvl = self.Rrlvl 658 elif os.path.isfile(lvlfile): 659 # print ' reading reclvl file' 660 self.Rrlvl = io.cireclvlRead(self.IonStr, '.'+lvlType) 661 lvl = self.Rrlvl 662 else: 663 self.RrlvlRate = {'rate':np.zeros_like(temperature)} 664 return 665 # 666 # the rates and temperatures in reclvl are not all the same 667 # 668 ntemp = temperature.size 669 nlvl = len(lvl['lvl1']) 670 if ntemp == 1: 671 rate = np.zeros(( nlvl), 'float64') 672 # previous takes care of temperatures below reclvl['temperature'].min() 673 if temperature > lvl['temperature'].max(): 674 # extrapolate as 1/temperature 675 for itrans in range(nlvl): 676 # lvl2 = self.Reclvl['lvl2'][itrans] 677 # lvlTemp = lvl['ntemp'][itrans] 678 rate[itrans] = lvl['rate'][itrans,-1]*(lvl['temperature'][itrans, -1]/temperature) 679 elif temperature < lvl['temperature'].min(): 680 # rate is already set to zero 681 pass 682 else: 683 for itrans in range(nlvl): 684 lvl2 = lvl['lvl2'][itrans] 685 nrecTemp = lvl['ntemp'][itrans] 686 y2 = interpolate.splrep(np.log(lvl['temperature'][itrans][:nrecTemp]), np.log(lvl['rate'][itrans][:nrecTemp])) 687 cirec = np.exp(interpolate.splev(np.log(temperature),y2)) 688 rate[itrans] = cirec.squeeze() 689 else: 690 # ntemp > 1 691 rate = np.zeros(( nlvl, temperature.size), 'float64') 692 # 693 for itrans in range(nlvl): 694 lvl2 = lvl['lvl2'][itrans] 695 nTemp = lvl['ntemp'][itrans] 696 y2 = interpolate.splrep(np.log(lvl['temperature'][itrans, :nTemp]), np.log(lvl['rate'][itrans, :nTemp])) 697 goodLow = temperature < lvl['temperature'][itrans].min() 698 if goodLow.sum() >0: 699 # print ' number of low temperatures = ', goodLow.sum() 700 lowT = temperature[goodLow] 701 good1 = temperature >= lvl['temperature'][itrans].min() 702 good2 = temperature <= lvl['temperature'][itrans].max() 703 realgood = np.logical_and(good1,good2) 704 if realgood.sum() > 0: 705 # print ' number of mid temperatures = ', realgood.sum() 706 midT = temperature[realgood] 707 goodHigh = temperature > lvl['temperature'][itrans].max() 708 if goodHigh.sum() > 0: 709 # print ' number of high temperatures = ', goodHigh.sum() 710 highT = temperature[goodHigh] 711 lvl2 = lvl['lvl2'][itrans] 712 nTemp = lvl['ntemp'][itrans] 713 newRate = np.zeros(ntemp, 'float64') 714 index = 0 715 if goodLow.sum() == 1: 716 # lowRec = np.exp(interpolate.splev(np.log(lowT),y2)) 717 # newRec[index] = lowRec 718 newRate[index] = 0. 719 index += 1 720 elif goodLow.sum() > 1: 721 # lowRec = np.exp(interpolate.splev(np.log(lowT),y2)) 722 for idx in range(goodLow.sum()): 723 # newRec[index] = lowRec[idx] 724 newRate[index] = 0. 725 index += 1 726 if realgood.sum() == 1: 727 midRec = np.exp(interpolate.splev(np.log(midT),y2)) 728 newRate[index] = midRec 729 index += 1 730 elif realgood.sum() > 1: 731 midRec = np.exp(interpolate.splev(np.log(midT),y2)) 732 for idx in range(realgood.sum()): 733 newRate[index] = midRec[idx] 734 index += 1 735 if goodHigh.sum() == 1: 736 # highRec = np.exp(interpolate.splev(np.log(highT),y2)) 737 # newRec[index] = highRec 738 newRate[index] = 0. 739 index += 1 740 elif goodHigh.sum() > 1: 741 # highRec = np.exp(interpolate.splev(np.log(highT),y2)) 742 for idx in range(goodHigh.sum()): 743 # print ' index, idx = ', index, idx 744 # newRec[index] = highRec[idx] 745 newRate[index] = 0. 746 index += 1 747 rate[itrans] = newRate 748 if lvlType == 'reclvl': 749 self.ReclvlRate = {'rate':rate, 'lvl1':lvl['lvl1'], 'lvl2':lvl['lvl2'], 'temperature':temperature} 750 elif lvlType == 'cilvl': 751 self.CilvlRate = {'rate':rate, 'lvl1':lvl['lvl1'], 'lvl2':lvl['lvl2'], 'temperature':temperature} 752 elif lvlType == 'rrlvl': 753 self.RrlvlRate = {'rate':rate, 'lvl1':lvl['lvl1'], 'lvl2':lvl['lvl2'], 'temperature':temperature}

754 755 # 756 # ------------------------------------------------------------------------------------- 757 # 758 # def reclvlDescale(self): 759 # ''' 760 # Interpolate and extrapolate reclvl rates. 761 # 762 # Used in level population calculations. 763 # I think this is made redundant by cireclvlDescale 764 # ''' 765 # if hasattr(self, 'Temperature'): 766 # temperature=self.Temperature 767 # else: 768 # print ' temperature is not defined' 769 # self.ReclvlRate = None 770 # return 771 # reclvlfile = util.ion2filename(self.IonStr)+'.reclvl' 772 # if hasattr(self, 'Reclvl'): 773 # reclvl = self.Reclvl 774 # elif os.path.isfile(reclvlfile): 775 ## print ' reading reclvl file' 776 # reclvl = io.cireclvlRead(self.IonStr, 'reclvl') 777 # else: 778 # self.ReclvlRate = {'rate':zeros_like(temperature)} 779 # return 780 # # 781 # # the rates and temperatures in reclvl are not all the same 782 # # 783 # ntemp = temperature.size 784 # if ntemp == 1: 785 # recRate = np.zeros(( len(reclvl['lvl1'])), 'float64') 786 # # previous takes care of temperatures below reclvl['temperature'].min() 787 # if temperature > reclvl['temperature'].max(): 788 # # extrapolate as 1/temperature 789 # for itrans in range(len(reclvl['lvl1'])): 790 ## lvl2 = self.Reclvl['lvl2'][itrans] 791 # nrecTemp = reclvl['ntemp'][itrans] 792 # recRate[itrans] = self.Reclvl['rate'][itrans,nrecTemp-1]*(reclvl['temperature'][itrans, nrecTemp-1]/temperature) 793 # else: 794 # for itrans in range(len(self.Reclvl['lvl1'])): 795 # lvl2 = self.Reclvl['lvl2'][itrans] 796 # nrecTemp = self.Reclvl['ntemp'][itrans] 797 # y2 = interpolate.splrep(np.log(self.Reclvl['temperature'][itrans, :nrecTemp]), np.log(self.Reclvl['rate'][itrans, :nrecTemp])) 798 # rec = np.exp(interpolate.splev(np.log(temperature),y2)) 799 # recRate[itrans] = rec.squeeze() 800 # else: 801 # # ntemp > 1 802 # recRate = np.zeros(( len(reclvl['lvl1']), temperature.size), 'float64') 803 # # 804 # for itrans in range(len(reclvl['lvl1'])): 805 # lvl2 = self.Reclvl['lvl2'][itrans] 806 # nrecTemp = self.Reclvl['ntemp'][itrans] 807 # y2 = interpolate.splrep(np.log(self.Reclvl['temperature'][itrans, :nrecTemp]), np.log(self.Reclvl['rate'][itrans, :nrecTemp])) 808 # goodLow = temperature < self.Reclvl['temperature'][itrans].min() 809 # if goodLow.sum() >0: 810 ## print ' number of low temperatures = ', goodLow.sum() 811 # lowT = temperature[goodLow] 812 # good1 = temperature >= self.Reclvl['temperature'][itrans].min() 813 # good2 = temperature <= self.Reclvl['temperature'][itrans].max() 814 # realgood = np.logical_and(good1,good2) 815 # if realgood.sum() > 0: 816 ## print ' number of mid temperatures = ', realgood.sum() 817 # midT = temperature[realgood] 818 # goodHigh = temperature > self.Reclvl['temperature'][itrans].max() 819 # if goodHigh.sum() > 0: 820 ## print ' number of high temperatures = ', goodHigh.sum() 821 # highT = temperature[goodHigh] 822 # lvl2 = self.Reclvl['lvl2'][itrans] 823 # nrecTemp = self.Reclvl['ntemp'][itrans] 824 # newRec = np.zeros(ntemp, 'float64') 825 # index = 0 826 # if goodLow.sum() == 1: 827 ## lowRec = np.exp(interpolate.splev(np.log(lowT),y2)) 828 ## newRec[index] = lowRec 829 # newRec[index] = 0. 830 # index += 1 831 # elif goodLow.sum() > 1: 832 ## lowRec = np.exp(interpolate.splev(np.log(lowT),y2)) 833 # for idx in range(goodLow.sum()): 834 ## newRec[index] = lowRec[idx] 835 # newRec[index] = 0. 836 # index += 1 837 # if realgood.sum() == 1: 838 # midRec = np.exp(interpolate.splev(np.log(midT),y2)) 839 # newRec[index] = midRec 840 # index += 1 841 # elif realgood.sum() > 1: 842 # midRec = np.exp(interpolate.splev(np.log(midT),y2)) 843 # for idx in range(realgood.sum()): 844 # newRec[index] = midRec[idx] 845 # index += 1 846 # if goodHigh.sum() == 1: 847 ## highRec = np.exp(interpolate.splev(np.log(highT),y2)) 848 ## newRec[index] = highRec 849 # newRec[index] = 0. 850 # index += 1 851 # elif goodHigh.sum() > 1: 852 ## highRec = np.exp(interpolate.splev(np.log(highT),y2)) 853 # for idx in range(goodHigh.sum()): 854 ## print ' index, idx = ', index, idx 855 ## newRec[index] = highRec[idx] 856 # newRec[index] = 0. 857 # index += 1 858 # recRate[itrans] = newRec 859 # self.ReclvlRate = {'rate':recRate, 'lvl2':reclvl['lvl2'], 'temperature':temperature} 860 # 861 # ------------------------------------------------------------------------------------- 862 #

863 - def drRateLvl(self, verbose=0):

864 ''' 865 to calculate the level resolved dielectronic rate from the higher ionization stage to 866 the ion of interest 867 rates are determined from autoionizing A-values 868 the dictionary self.DrRateLvl contains 869 rate = the dielectronic rate into an autoionizing level 870 effRate = the dielectronic rate into an autoionizing level mutilplied by the branching 871 ratio for a stabilizing transition 872 totalRate = the sum of all the effRates 873 ''' 874 if not hasattr(self, 'Higher'): 875 nameStuff = util.convertName(self.IonStr) 876 z = nameStuff['Z'] 877 stage = nameStuff['Ion'] 878 higherStr = util.zion2name(z, stage+1) 879 self.Higher = ion(higherStr, self.Temperature, self.EDensity) 880 # 881 # (4 pi a0^2)^(3/2) = 6.6011e-24 (Badnell et al, 2003, A&A 406, 1151 882 # coef1 = 6.6011e-24*(const.hartree/(2.*const.boltzmann*self.Temperature))**1.5 883 coef2 = (const.planck)**3/(2.*const.pi*const.emass*const.boltzmann*self.Temperature)**1.5 884 # next from Aped 885 # coef3 = (4.*const.pi/(const.boltzmannEv*self.Temperature/const.ryd2Ev))**(1.5)*(const.bohr)**3 886 coef = coef2 887 # print ' coefs = ', coef1, coef2, coef3 888 nt = self.Temperature.size 889 allRate = [] 890 effRate = [] 891 de = [] 892 erg = [] 893 ipErg = self.Ip*const.ev2Erg 894 lvl = [] 895 branch = [] 896 dekt = [] 897 totalRate = np.zeros(nt, 'float64') 898 # lvlformat = '%7i%7i%10.2e%10.2e' 899 for i, avalue in enumerate(self.Auto['avalue']): 900 elvl1idx = self.Elvlc['lvl'].index(self.Auto['lvl1'][i]) 901 elvl2idx = self.Elvlc['lvl'].index(self.Auto['lvl2'][i]) 902 gUpper = float(self.Elvlc['mult'][elvl2idx]) 903 gLower = float(self.Higher.Elvlc['mult'][elvl1idx]) 904 # print i, autoa['lvl2'][i], gLower, gUpper, avalue 905 ecm2 = self.Elvlc['ecm'][elvl2idx] 906 if ecm2 < 0.: 907 ecm2 = self.Elvlc['ecmth'][elvl2idx] 908 de1 = ecm2*const.invCm2Erg - self.Ip*const.ev2Erg 909 erg.append(ecm2*const.invCm2Erg) 910 #de1 = ecm2*const.invCm2Ev - self.Ip 911 de.append(de1) 912 dekt1 = de1/(const.boltzmann*self.Temperature) 913 dekt.append(dekt1) 914 #dekt = de1/(const.boltzmannEv*self.Temperature) 915 #expkt = np.exp(-de1/(const.boltzmann*self.Temperature)) 916 expkt = np.exp(-dekt1) 917 rate = coef*gUpper*expkt*avalue/(2.*gLower) 918 branch1 = self.Wgfa['avalueLvl'][elvl2idx]/(avalue + self.Wgfa['avalueLvl'][elvl2idx]) 919 branch.append(branch1) 920 lvl.append(self.Auto['lvl2'][i]) 921 # print i, self.Auto['lvl2'][i], cnt 922 # print i, j, self.Wgfa['lvl1'][idx[-1 ]], self.Wgfa['lvl2'][idx[-1]], self.Auto['lvl2'][i] 923 # print ' lvl2, rate = ', autoa['lvl2'][i], rate 924 # lvlstr = lvlformat%(self.Auto['lvl1'][i], self.Auto['lvl2'][i], avalue, branch) 925 # if verbose: 926 # tstr = lvlstr 927 # rstr = lvlstr 928 # if nt == 1: 929 # tstr += '%10.2e'%(self.Temperature) 930 # rstr += '%10.2e%10.2e'%(rate, rate*branch) 931 # else: 932 # for it in range(nt): 933 # tstr += '%10.2e'%(self.Temperature[it]) 934 # rstr += '%10.2e'%(rate[it]) 935 # print tstr 936 # print rstr 937 allRate.append(rate) 938 effRate.append(rate*branch1) 939 totalRate += rate*branch1 940 # outpt.write(tstr +'\n') 941 # outpt.write(rstr + '\n') 942 self.DrRateLvl = {'rate':allRate, 'effRate':effRate, 'totalRate':totalRate, 'de':de, 'avalue':self.Auto['avalue'], 'lvl':lvl, 'branch':branch, 'dekt':dekt, 'erg':erg, 'ipErg':ipErg} #, 'lvl1':lvl1, 'lvl2':lvl2} - in self.Auto

943 # 944 # ------------------------------------------------------------------------------------- 945 #

946 - def recombRate(self):

947 '''Provides the total recombination rate coefficient. 948 949 Calls drRate and rrRate''' 950 # 951 if hasattr(self, 'Temperature'): 952 temperature=self.Temperature 953 else: 954 print(' temperature is not defined') 955 self.RecombRate = {'errorMessage':' temperature is not defined'} 956 if self.Ion == 1: 957 # print ' this is a neutral and has no recombination rate' 958 self.RecombRate = {'rate':np.zeros_like(temperature), 'temperature':temperature} 959 return 960 self.rrRate() 961 self.drRate() 962 if not hasattr(self, 'DrRate'): 963 rate=self.RrRate['rate'] 964 else: 965 rate=self.RrRate['rate']+self.DrRate['rate'] 966 self.RecombRate={'rate':rate, 'temperature':temperature} 967 return

968 # 969 # ------------------------------------------------------------------------------------- 970 #

971 - def p2eRatio(self):

972 '''Calculates the proton density to electron density ratio. 973 974 Uses the abundance and ionization equilibrium.''' 975 # version 7 now get the abundancd from chianti.data/chdata 976 #if hasattr(self, 'Abundance'): 977 #ab=self.Abundance 978 #else: 979 #abundName = self.Defaults['abundfile'] 980 #io.abundanceRead(abundancename = abundName) 981 # 982 #if not hasattr(self, 'IoneqAll'): 983 #self.IoneqAll = io.ioneqRead(ioneqname = ioneqname) 984 # 985 if hasattr(self, 'Temperature'): 986 temperature=self.Temperature 987 else: 988 temperature = self.IoneqAll['ioneqTemperature'] 989 if not hasattr(self, 'AbundanceName'): 990 self.AbundanceName = self.Defaults['abundfile'] 991 # 992 # temperature=self.IoneqTemperature 993 # else: temperature=np.asarray(temperature,'float32') 994 # 995 nTemp=temperature.size 996 #nEl = chdata.Abundance[self.AbundanceName]['abundance'].size 997 # 998 eDensity=np.zeros(nTemp,'Float64') 999 pDensity=np.zeros(nTemp,'Float64') 1000 ionDensity=0. 1001 # zDensity=np.zeros(nTemp, 'float32') 1002 # 1003 # 1004 # only hydrogen contributes to the proton density 1005 anEl = 0 1006 ion = 1 1007 good = self.IoneqAll['ioneqAll'][anEl,ion] > 0. 1008 y2 = interpolate.splrep(np.log(self.IoneqAll['ioneqTemperature'][good]),np.log(self.IoneqAll['ioneqAll'][anEl,ion,good]),s=0) 1009 bad1 = np.log(temperature) < np.log(self.IoneqAll['ioneqTemperature'][good].min()) 1010 bad2 = np.log(temperature) > np.log(self.IoneqAll['ioneqTemperature'][good].max()) 1011 bad=np.logical_or(bad1,bad2) 1012 goodt=np.logical_not(bad) 1013 thisIoneq=np.where(goodt,10.**interpolate.splev(np.log(temperature),y2,der=0),0.) 1014 pDensity+=chdata.Abundance[self.AbundanceName]['abundance'][anEl]*thisIoneq 1015 # ionDensity=self.Abundance[anEl] 1016 # 1017 # all the rest do contribute to the electron and ion densities 1018 El=[iEl for iEl in range(50) if chdata.Abundance[self.AbundanceName]['abundance'][iEl] > 0.] 1019 for anEl in El: 1020 ionDensity+=chdata.Abundance[self.AbundanceName]['abundance'][anEl] 1021 for ion in range(1,anEl+2): 1022 good = self.IoneqAll['ioneqAll'][anEl,ion] > 0. 1023 y2 = interpolate.splrep(np.log(self.IoneqAll['ioneqTemperature'][good]),np.log(self.IoneqAll['ioneqAll'][anEl,ion,good]),s=0) 1024 bad1 = np.log(temperature) < np.log(self.IoneqAll['ioneqTemperature'][good].min()) 1025 bad2 = np.log(temperature) > np.log(self.IoneqAll['ioneqTemperature'][good].max()) 1026 bad = np.logical_or(bad1,bad2) 1027 goodt = np.logical_not(bad) 1028 # good1 = temperature >= self.IoneqAll['ioneqTemperature'][good].min() 1029 # good2 = temperature <= self.IoneqAll['ioneqTemperature'][good].min() 1030 # goodt = np.logical_and(good1, good2) 1031 # thisIoneq=np.where(goodt,10.**interpolate.splev(np.log(temperature),y2,der=0),0.) 1032 # for temperatures outside the range of the ioneq, the default is set to 1, since 1033 # this is the denominator 1034 thisIoneq=np.where(goodt,10.**interpolate.splev(np.log(temperature),y2,der=0),1.) 1035 eDensity+=float(ion)*chdata.Abundance[self.AbundanceName]['abundance'][anEl]*thisIoneq 1036 self.ProtonDensityRatio=pDensity/eDensity 1037 self.EDensity=eDensity 1038 self.IonDensity=ionDensity 1039 self.IonDensityRatio=ionDensity/eDensity

1040 # return # {'PDensityRatio':self.PDensityRatio,'PDensity':pDensity,'EDensity':eDensity} 1041 # 1042 # ------------------------------------------------------------------------------------- 1043 #

1044 - def upsilonDescale(self, prot=0, diel=0):

1045 """ 1046 Provides the temperatures and effective collision strengths (upsilons) 1047 set prot for proton rates 1048 otherwise, ce will be set for electron collision rates 1049 uses the new format "scups" files 1050 """ 1051 # 1052 # xt=kt/de 1053 # 1054 # 1055 if prot: 1056 ce = 0 1057 try: 1058 nscups=len(self.Psplups["lvl1"]) 1059 except: 1060 self.Psplups=io.splupsRead(self.IonStr,prot=1) 1061 if type(self.Psplups) == type(None): 1062 self.PUpsilon = None 1063 return 1064 else: 1065 nscups = len(self.Cilvl["lvl1"]) 1066 elif diel: 1067 ce = 0 1068 try: 1069 nsplups = len(self.DielSplups["lvl1"]) 1070 except: 1071 self.DielSplups = io.scupsRead(self.IonStr) 1072 if type(self.DielSplups) == type(None): 1073 self.DielUpsilon = None 1074 return 1075 else: 1076 nscups = len(self.Scups["lvl1"]) 1077 else: 1078 ce=1 1079 try: 1080 nscups=len(self.Scups["lvl1"]) 1081 except: 1082 self.Scups = io.scupsRead(self.IonStr) 1083 if not self.Scups['status']: 1084 self.Upsilon = None 1085 return 1086 else: 1087 nscups = len(self.Scups["lvl1"]) 1088 # 1089 # 1090 if hasattr(self, 'Temperature'): 1091 temperature=self.Temperature 1092 else: 1093 print(' Temperature undefined') 1094 return {'errorMessage':' Temperature undefined'} 1095 # 1096 # if hasattr(self, 'Elvlc'): 1097 # nlvls=len(self.Elvlc["lvl"]) 1098 # else: 1099 # self.elvlcRead() 1100 if not hasattr(self, 'Elvlc'): 1101 self.elvlcRead() 1102 1103 # nlvls=len(self.Elvlc["lvl"]) 1104 # 1105 # need to make sure elvl is >0, except for ground level 1106 eryd=np.asarray(self.Elvlc["eryd"]) 1107 erydth=np.asarray(self.Elvlc["erydth"]) 1108 elvlc=np.where(eryd >= 0.,eryd,erydth) 1109 ## de=self.Elvlc["de"] 1110 temp=np.asarray(temperature) 1111 ntemp=temp.size 1112 if ntemp > 1: 1113 ups = np.zeros((nscups,ntemp),"Float64") 1114 exRate = np.zeros((nscups,ntemp),"Float64") 1115 dexRate = np.zeros((nscups,ntemp),"Float64") 1116 else: 1117 ups = np.zeros(nscups,"Float64") 1118 exRate = np.zeros((nscups,ntemp),"Float64") 1119 dexRate = np.zeros((nscups,ntemp),"Float64") 1120 deAll = [] 1121 # 1122 for iscups in range(nscups): 1123 if prot: 1124 # for proton rates 1125 l1=self.Psplups["lvl1"][iscups]-1 1126 l1idx = self.Elvlc['lvl'].index(self.Psplups['lvl1'][iscups]) 1127 l2=self.Psplups["lvl2"][iscups]-1 1128 l2idx = self.Elvlc['lvl'].index(self.Psplups['lvl2'][iscups]) 1129 ttype=self.Psplups["ttype"][iscups] 1130 cups=self.Psplups["cups"][iscups] 1131 nspl=self.Psplups["nspl"][iscups] 1132 dx=1./(float(nspl)-1.) 1133 xs=dx*np.arange(nspl) 1134 # splups=self.Psplups["splups"][isplups,0:nspl] 1135 splups=self.Psplups["splups"][iscups] 1136 de=elvlc[l2idx]-elvlc[l1idx] 1137 # de=self.Psplups['de'][isplups] # these are generally 0. 1138 kte = const.boltzmann*temp/(de*const.ryd2erg) 1139 # elif diel: 1140 # # 1141 # l1 = self.Scups["lvl1"][isplups]-1 1142 # l1idx = self.Elvlc['lvl'].index(self.Scups['lvl1'][isplups]) 1143 # l2 = self.Scups["lvl2"][isplups]-1 1144 # l2idx = self.Elvlc['lvl'].index(self.Scups['lvl2'][isplups]) 1145 # ttype = self.Scups["ttype"][isplups] 1146 # cups = self.Scups["cups"][isplups] 1147 # nspl = self.Scups["nspl"][isplups] 1148 # ttype = self.Scups["ttype"][isplups] 1149 # dx = 1./(float(nspl)-1.) 1150 # xs=dx*np.arange(nspl) 1151 ## splups = self.Scups["splups"][isplups,0:nspl] 1152 # splups = self.Scups["splups"][isplups] 1153 # de=self.Scups['de'][isplups] 1154 # kte = const.boltzmann*temp/(de*const.ryd2erg) 1155 else: 1156 # electron collisional excitation or dielectronic excitation 1157 l1=self.Scups["lvl1"][iscups]-1 1158 l1idx = self.Elvlc['lvl'].index(self.Scups['lvl1'][iscups]) 1159 l2=self.Scups["lvl2"][iscups]-1 1160 l2idx = self.Elvlc['lvl'].index(self.Scups['lvl2'][iscups]) 1161 ttype=self.Scups["ttype"][iscups] 1162 cups=self.Scups["cups"][iscups] 1163 nspl=self.Scups["ntemp"][iscups] 1164 xs = self.Scups['btemp'][iscups] 1165 # splups=self.Scups["splups"][isplups,0:nspl] 1166 scups=self.Scups["bscups"][iscups] 1167 de=self.Scups['de'][iscups] 1168 kte = const.boltzmann*temp/(de*const.ryd2erg) 1169 # 1170 der=0 1171 if ttype == 1: 1172 st=1.-np.log(cups)/np.log(kte+cups) 1173 y2=interpolate.splrep(xs,scups,s=0) 1174 sups=interpolate.splev(st,y2,der=der) 1175 # sups=interpolate.spline(xs, splups, st) 1176 ups[iscups]=sups*np.log(kte+np.exp(1.)) 1177 # 1178 if ttype == 2: 1179 st=kte/(kte+cups) 1180 y2=interpolate.splrep(xs,scups,s=0) 1181 sups=interpolate.splev(st,y2,der=der) 1182 ups[iscups]=sups 1183 # 1184 if ttype == 3: 1185 st=kte/(kte+cups) 1186 y2=interpolate.splrep(xs,scups,s=0) 1187 sups=interpolate.splev(st,y2,der=der) 1188 ups[iscups]=sups/(kte+1.) 1189 # 1190 if ttype == 4: 1191 st=1.-np.log(cups)/np.log(kte+cups) 1192 y2=interpolate.splrep(xs,scups,s=0) 1193 sups=interpolate.splev(st,y2,der=der) 1194 ups[iscups]=sups*np.log(kte+cups) 1195 # 1196 if ttype == 5: 1197 # dielectronic rates 1198 st=kte/(kte+cups) 1199 # xs=dx*np.arange(nspl) 1200 y2=interpolate.splrep(xs,scups,s=0) 1201 sups=interpolate.splev(st,y2,der=der) 1202 ups[iscups]=sups/(kte+0.) 1203 # 1204 # descale proton values 1205 if ttype == 6: 1206 st=kte/(kte+cups) 1207 # xs=dx*np.arange(nspl) 1208 y2=interpolate.splrep(xs,scups,s=0) 1209 sups=interpolate.splev(st,y2,der=der) 1210 # ups[iscups] = sups 1211 ups[iscups]=10.**sups 1212 # 1213 elif ttype > 6: print(' t_type ne 1,2,3,4,5 = %5i %5i %5i '%(ttype,l1,l2)) 1214 # 1215 if ce: 1216 if self.Dielectronic: 1217 # the dielectronic ions will eventually be discontinued 1218 de = np.abs((elvlc[l2idx] - self.UpperIp/const.ryd2Ev) - elvlc[l1idx]) 1219 else: 1220 de = np.abs(elvlc[l2idx] - elvlc[l1idx]) 1221 deAll.append(de) 1222 # print ' ce lvl1 %5i lvl2 %5i de %10.2e'%(l1, l2, de) 1223 ekt = (de*const.ryd2erg)/(const.boltzmann*temp) 1224 fmult1 = float(self.Elvlc["mult"][l1idx]) 1225 fmult2 = float(self.Elvlc["mult"][l2idx]) 1226 dexRate[iscups] = const.collision*ups[iscups]/(fmult2*np.sqrt(temp)) 1227 exRate[iscups] = const.collision*ups[iscups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 1228 elif diel: 1229 # print ' diel lvl1 %5i lvl2 %5i de %10.2e'%(l1, l2, de) 1230 de = np.abs((elvlc[l2idx] - self.Ip/const.ryd2Ev) - elvlc[l1idx]) 1231 ekt = (de*const.ryd2erg)/(const.boltzmann*temp) 1232 fmult1 = float(self.Elvlc["mult"][l1idx]) 1233 fmult2 = float(self.Elvlc["mult"][l2idx]) 1234 exRate[iscups] = const.collision*ups[iscups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 1235 elif prot: 1236 de = np.abs(elvlc[l2idx]- elvlc[l1idx]) 1237 ekt = (de*1.57888e+5)/temp 1238 fmult1 = float(self.Elvlc["mult"][l1idx]) 1239 fmult2 = float(self.Elvlc["mult"][l2idx]) 1240 dexRate[iscups] = const.collision*ups[iscups]/(fmult2*np.sqrt(temp)) 1241 exRate[iscups] = const.collision*ups[iscups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 1242 # 1243 ups=np.where(ups > 0.,ups,0.) 1244 # 1245 if prot == 1: 1246 self.PUpsilon = {'upsilon':ups, 'temperature':temperature, 'exRate':exRate, 'dexRate':dexRate} 1247 # elif diel == 1: 1248 # self.DielUpsilon = {'upsilon':ups, 'temperature':temperature, 'exRate':exRate} 1249 else: 1250 self.Upsilon = {'upsilon':ups, 'temperature':temperature, 'exRate':exRate, 'dexRate':dexRate, 'de':deAll}

1251 # 1252 # ------------------------------------------------------------------------------------- 1253 #

1254 - def upsilonDescaleSplups(self, prot=0, diel=0):

1255 """ 1256 Provides the temperatures and effective collision strengths (upsilons) 1257 set prot for proton rates 1258 otherwise, ce will be set for electron collision rates 1259 """ 1260 # 1261 # xt=kt/de 1262 # 1263 # 1264 if prot: 1265 ce = 0 1266 try: 1267 nsplups=len(self.Psplups["lvl1"]) 1268 except: 1269 self.Psplups=io.splupsRead(self.IonStr,prot=1) 1270 if type(self.Psplups) == type(None): 1271 self.PUpsilon = None 1272 return 1273 else: 1274 nsplups = len(self.Cilvl["lvl1"]) 1275 elif diel: 1276 ce = 0 1277 try: 1278 nsplups = len(self.DielSplups["lvl1"]) 1279 except: 1280 self.DielSplups = io.splupsRead(self.IonStr,diel=1) 1281 if type(self.DielSplups) == type(None): 1282 self.DielUpsilon = None 1283 return 1284 else: 1285 nsplups = len(self.DielSplups["lvl1"]) 1286 else: 1287 ce=1 1288 try: 1289 nsplups=len(self.Splups["lvl1"]) 1290 except: 1291 self.Splups = io.splupsRead(self.IonStr) 1292 if type(self.Splups) == type(None): 1293 self.Upsilon = None 1294 return 1295 else: 1296 nsplups = len(self.Splups["lvl1"]) 1297 # 1298 # 1299 if hasattr(self, 'Temperature'): 1300 temperature=self.Temperature 1301 else: 1302 print(' Temperature undefined') 1303 return {'errorMessage':' Temperature undefined'} 1304 # 1305 # if hasattr(self, 'Elvlc'): 1306 # nlvls=len(self.Elvlc["lvl"]) 1307 # else: 1308 # self.elvlcRead() 1309 # nlvls=len(self.Elvlc["lvl"]) 1310 if not hasattr(self, 'Elvlc'): 1311 self.elvlcRead() 1312 # 1313 # need to make sure elvl is >0, except for ground level 1314 eryd=np.asarray(self.Elvlc["eryd"]) 1315 erydth=np.asarray(self.Elvlc["erydth"]) 1316 elvlc=np.where(eryd > 0.,eryd,erydth) 1317 ## de=self.Elvlc["de"] 1318 temp=np.asarray(temperature) 1319 ntemp=temp.size 1320 if ntemp > 1: 1321 ups = np.zeros((nsplups,ntemp),"Float64") 1322 exRate = np.zeros((nsplups,ntemp),"Float64") 1323 dexRate = np.zeros((nsplups,ntemp),"Float64") 1324 else: 1325 ups = np.zeros(nsplups,"Float64") 1326 exRate = np.zeros((nsplups,ntemp),"Float64") 1327 dexRate = np.zeros((nsplups,ntemp),"Float64") 1328 deAll = [] 1329 # 1330 for isplups in range(nsplups): 1331 if prot: 1332 # for proton rates 1333 l1=self.Psplups["lvl1"][isplups]-1 1334 l2=self.Psplups["lvl2"][isplups]-1 1335 fmult1 = self.Elvlc['mult'][l1] 1336 fmult2 = self.Elvlc['mult'][l2] 1337 ttype=self.Psplups["ttype"][isplups] 1338 cups=self.Psplups["cups"][isplups] 1339 nspl=self.Psplups["nspl"][isplups] 1340 dx=1./(float(nspl)-1.) 1341 # splups=self.Psplups["splups"][isplups,0:nspl] 1342 splups=self.Psplups["splups"][isplups] 1343 de=elvlc[l2]-elvlc[l1] 1344 # de=self.Psplups['de'][isplups] # these are generally 0. 1345 kte = const.boltzmann*temp/(de*const.ryd2erg) 1346 elif diel: 1347 # 1348 l1 = self.DielSplups["lvl1"][isplups]-1 1349 l2 = self.DielSplups["lvl2"][isplups]-1 1350 ttype = self.DielSplups["ttype"][isplups] 1351 cups = self.DielSplups["cups"][isplups] 1352 nspl = self.DielSplups["nspl"][isplups] 1353 ttype = self.DielSplups["ttype"][isplups] 1354 dx = 1./(float(nspl)-1.) 1355 # splups = self.DielSplups["splups"][isplups,0:nspl] 1356 splups = self.DielSplups["splups"][isplups] 1357 de=self.DielSplups['de'][isplups] 1358 kte = const.boltzmann*temp/(de*const.ryd2erg) 1359 else: 1360 # electron collisional excitation 1361 l1=self.Splups["lvl1"][isplups]-1 1362 l2=self.Splups["lvl2"][isplups]-1 1363 ttype=self.Splups["ttype"][isplups] 1364 cups=self.Splups["cups"][isplups] 1365 nspl=self.Splups["nspl"][isplups] 1366 dx=1./(float(nspl)-1.) 1367 ## print self.Splups["splups"][l1,l2] 1368 # splups=self.Splups["splups"][isplups,0:nspl] 1369 splups=self.Splups["splups"][isplups] 1370 # de=elvlc[l2]-elvlc[l1] 1371 de=self.Splups['de'][isplups] 1372 kte = const.boltzmann*temp/(de*const.ryd2erg) 1373 # 1374 der=0 1375 if ttype == 1: 1376 st=1.-np.log(cups)/np.log(kte+cups) 1377 xs=dx*np.arange(nspl) 1378 y2=interpolate.splrep(xs,splups,s=0) 1379 sups=interpolate.splev(st,y2,der=der) 1380 # sups=interpolate.spline(xs, splups, st) 1381 ups[isplups]=sups*np.log(kte+np.exp(1.)) 1382 # 1383 if ttype == 2: 1384 st=kte/(kte+cups) 1385 xs=dx*np.arange(nspl) 1386 y2=interpolate.splrep(xs,splups,s=0) 1387 sups=interpolate.splev(st,y2,der=der) 1388 ups[isplups]=sups 1389 # 1390 if ttype == 3: 1391 st=kte/(kte+cups) 1392 xs=dx*np.arange(nspl) 1393 y2=interpolate.splrep(xs,splups,s=0) 1394 sups=interpolate.splev(st,y2,der=der) 1395 ups[isplups]=sups/(kte+1.) 1396 # 1397 if ttype == 4: 1398 st=1.-np.log(cups)/np.log(kte+cups) 1399 xs=dx*np.arange(nspl) 1400 y2=interpolate.splrep(xs,splups,s=0) 1401 sups=interpolate.splev(st,y2,der=der) 1402 ups[isplups]=sups*np.log(kte+cups) 1403 # 1404 if ttype == 5: 1405 # dielectronic rates 1406 st=kte/(kte+cups) 1407 xs=dx*np.arange(nspl) 1408 y2=interpolate.splrep(xs,splups,s=0) 1409 sups=interpolate.splev(st,y2,der=der) 1410 ups[isplups]=sups/(kte+0.) 1411 # 1412 # descale proton values 1413 if ttype == 6: 1414 st=kte/(kte+cups) 1415 xs=dx*np.arange(nspl) 1416 y2=interpolate.splrep(xs,splups,s=0) 1417 sups=interpolate.splev(st,y2,der=der) 1418 # ups[isplups] = sups 1419 ups[isplups]=10.**sups 1420 # 1421 elif ttype > 6: print(' t_type ne 1,2,3,4,5 = %5i %5i %5i'%(ttype,l1,l2)) 1422 # 1423 if ce: 1424 if self.Dielectronic: 1425 # the dielectronic ions will eventually be discontinued 1426 de = np.abs((elvlc[l2] - self.UpperIp/const.ryd2Ev) - elvlc[l1]) 1427 else: 1428 de = np.abs(elvlc[l2] - elvlc[l1]) 1429 deAll.append(de) 1430 # print ' ce lvl1 %5i lvl2 %5i de %10.2e'%(l1, l2, de) 1431 ekt = (de*const.ryd2erg)/(const.boltzmann*temp) 1432 fmult1 = float(self.Elvlc["mult"][l1]) 1433 fmult2 = float(self.Elvlc["mult"][l2]) 1434 dexRate[isplups] = const.collision*ups[isplups]/(fmult2*np.sqrt(temp)) 1435 exRate[isplups] = const.collision*ups[isplups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 1436 elif diel: 1437 # print ' diel lvl1 %5i lvl2 %5i de %10.2e'%(l1, l2, de) 1438 de = np.abs((elvlc[l2] - self.Ip/const.ryd2Ev) - elvlc[l1]) 1439 ekt = (de*const.ryd2erg)/(const.boltzmann*temp) 1440 fmult1 = float(self.Elvlc["mult"][l1]) 1441 fmult2 = float(self.Elvlc["mult"][l2]) 1442 exRate[isplups] = const.collision*ups[isplups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 1443 elif prot: 1444 de = np.abs(elvlc[l2]- elvlc[l1]) 1445 ekt = (de*1.57888e+5)/temp 1446 fmult1 = float(self.Elvlc["mult"][l1]) 1447 fmult2 = float(self.Elvlc["mult"][l2]) 1448 dexRate[isplups] = ups[isplups] 1449 exRate[isplups] = (fmult1/fmult2)*ups[isplups]*np.exp(-ekt) 1450 # 1451 ups=np.where(ups > 0.,ups,0.) 1452 # 1453 if prot: 1454 self.PUpsilon = {'temperature':temperature, 'exRate':exRate, 'dexRate':dexRate} 1455 elif diel: 1456 self.DielUpsilon = {'upsilon':ups, 'temperature':temperature, 'exRate':exRate} 1457 else: 1458 self.Upsilon = {'upsilon':ups, 'temperature':temperature, 'exRate':exRate, 'dexRate':dexRate, 'de':deAll}

1459 # 1460 # ------------------------------------------------------------------------- 1461 #

1462 - def setup(self, dir=0, verbose=0):

1463 ''' 1464 if ion is initiated with setup=0, this allows the setup to be done at a later point 1465 perhaps, more importantly, by setting dir to a directory cotaining the necessary files 1466 for a ChiantiPy ion, it allows one to setup an ion with files not in the current 1467 Chianti directory 1468 ''' 1469 # 1470 # read in all data if in masterlist 1471 # if not, there should still be ionization and recombination rates 1472 # 1473 MasterList = chdata.MasterList 1474 # 1475 if self.IonStr in MasterList: 1476 if dir: 1477 fileName = os.path.join(dir, self.IonStr) 1478 self.Elvlc = io.elvlcRead('',filename=fileName+'.elvlc') 1479 self.Wgfa = io.wgfaRead('',filename=fileName+'.wgfa', elvlcname=fileName+'.elvlc') 1480 self.Nwgfa=len(self.Wgfa['lvl1']) 1481 nlvlWgfa = max(self.Wgfa['lvl2']) 1482 nlvlList =[nlvlWgfa] 1483 # print 'fileName = ', fileName 1484 scupsfile = fileName + '.scups' 1485 # read the splups file 1486 if os.path.isfile(scupsfile): 1487 # happens the case of fe_3 and prob. a few others 1488 self.Scups = io.scupsRead('', filename=scupsfile) 1489 self.Nscups=len(self.Scups['lvl1']) 1490 nlvlScups = max(self.Scups['lvl2']) 1491 nlvlList.append(nlvlScups) 1492 else: 1493 self.Nscups = 0 1494 nlvlScups = 0 1495 else: 1496 fileName = util.ion2filename(self.IonStr) 1497 self.Elvlc = io.elvlcRead(self.IonStr) 1498 self.Wgfa = io.wgfaRead(self.IonStr) 1499 self.Nwgfa=len(self.Wgfa['lvl1']) 1500 nlvlWgfa = max(self.Wgfa['lvl2']) 1501 nlvlList =[nlvlWgfa] 1502 # print 'fileName = ', fileName 1503 scupsfile = fileName + '.scups' 1504 if os.path.isfile(scupsfile): 1505 # happens the case of fe_3 and prob. a few others 1506 self.Scups = io.scupsRead(self.IonStr) 1507 self.Nscups=len(self.Scups['lvl1']) 1508 nlvlScups = max(self.Scups['lvl2']) 1509 nlvlList.append(nlvlScups) 1510 else: 1511 self.Nscups = 0 1512 nlvlScups = 0 1513 ## self.Nlvls = nlvlElvlc 1514 # 1515 file = fileName +'.cilvl' 1516 if os.path.isfile(file): 1517 self.Cilvl = io.cireclvlRead('',filename = fileName, cilvl=1) 1518 self.Ncilvl=len(self.Cilvl['lvl1']) 1519 nlvlCilvl = max(self.Cilvl['lvl2']) 1520 nlvlList.append(nlvlCilvl) 1521 else: 1522 self.Ncilvl = 0 1523 # .reclvl file may not exist 1524 reclvlfile = fileName +'.reclvl' 1525 if os.path.isfile(reclvlfile): 1526 self.Reclvl = io.cireclvlRead('',filename=fileName, reclvl=1) 1527 self.Nreclvl = len(self.Reclvl['lvl1']) 1528 nlvlReclvl = max(self.Reclvl['lvl2']) 1529 nlvlList.append(nlvlReclvl) 1530 else: 1531 self.Nreclvl = 0 1532 # .dielsplups files not longer exist 1533 # dielsplupsfile = fileName +'.splups' 1534 # if self.Dielectronic and os.path.isfile(dielsplupsfile): 1535 # self.DielSplups = io.splupsRead('', filename=dielsplupsfile, diel=1) 1536 # self.Ndielsplups=len(self.DielSplups["lvl1"]) 1537 # nlvlDielSplups = max(self.DielSplups['lvl2']) 1538 # nlvlList.append(nlvlDielSplups) 1539 # else: 1540 # self.Ndielsplups = 0 1541 # 1542 # psplups file may not exist 1543 psplupsfile = fileName +'.psplups' 1544 if os.path.isfile(psplupsfile): 1545 self.Psplups = io.splupsRead('', filename=psplupsfile, prot=True) 1546 self.Npsplups=len(self.Psplups["lvl1"]) 1547 else: 1548 self.Npsplups = 0 1549 # 1550 drparamsFile = fileName +'.drparams' 1551 if os.path.isfile(drparamsFile): 1552 self.DrParams = io.drRead(self.IonStr) 1553 # 1554 rrparamsFile = fileName +'.rrparams' 1555 if os.path.isfile(rrparamsFile): 1556 self.RrParams = io.rrRead(self.IonStr) 1557 1558 # not needed for ion, only phion 1559 # photoxfile = util.ion2filename(self.IonStr)+'.photox' 1560 # if os.path.isfile(photoxfile): 1561 # self.Photox = util.photoxRead(self.IonStr) 1562 # 1563 # need to determine the number of levels that can be populated 1564 nlvlElvlc = len(self.Elvlc['lvl']) 1565 # print ' nlvlElvlc = ', nlvlElvlc 1566 # print ' other nlvls = ', nlvlList 1567 # nlvlWgfa = max(self.Wgfa['lvl2']) 1568 # elvlc file can have more levels than the rate level files 1569 self.Nlvls = min([nlvlElvlc, max(nlvlList)]) 1570 else: 1571 try: 1572 self.Elvlc = io.elvlcRead(self.IonStr, verbose=verbose) 1573 except: 1574 print(' the ion %s is not in the CHIANTI masterlist '%(self.IonStr)) 1575 print(' elvlc file NOT available for %s'%(self.IonStr)) 1576 return

1577 # 1578 # ------------------------------------------------------------------------- 1579 #

1580 - def setupIonrec(self, dir=0, verbose=0):

1581 ''' 1582 this allows a bare-bones ion object to be setup up with just the ionization and recombination rates 1583 mainly for ions without a complete set of files - one that is not in the MasterList 1584 ''' 1585 # 1586 # 1587 if dir: 1588 fileName = os.path.join(dir, self.IonStr) 1589 else: 1590 fileName = util.ion2filename(self.IonStr) 1591 # 1592 elvlcname=fileName+'.elvlc' 1593 if os.path.isfile(elvlcname): 1594 self.Elvlc = io.elvlcRead('',elvlcname) 1595 else: 1596 zstuff = util.convertName(self.IonStr) 1597 if zstuff['Ion'] - zstuff['Z'] != 1: 1598 # don't expect one for the bare ion 1599 print(' Elvlc file missing for '+self.IonStr) 1600 return 1601 # 1602 file = fileName +'.cilvl' 1603 if os.path.isfile(file): 1604 self.Cilvl = io.cireclvlRead('',filename = fileName, cilvl=1) 1605 self.Ncilvl=len(self.Cilvl['lvl1']) 1606 # nlvlCilvl = max(self.Cilvl['lvl2']) 1607 # nlvlList.append(nlvlCilvl) 1608 else: 1609 self.Ncilvl = 0 1610 # .reclvl file may not exist 1611 reclvlfile = fileName +'.reclvl' 1612 if os.path.isfile(reclvlfile): 1613 self.Reclvl = io.cireclvlRead('',filename=fileName, reclvl=1) 1614 self.Nreclvl = len(self.Reclvl['lvl1']) 1615 # nlvlReclvl = max(self.Reclvl['lvl2']) 1616 # nlvlList.append(nlvlReclvl) 1617 else: 1618 self.Nreclvl = 0 1619 # # 1620 # # 1621 drparamsFile = fileName +'.drparams' 1622 if os.path.isfile(drparamsFile): 1623 self.DrParams = io.drRead(self.IonStr) 1624 # 1625 rrparamsFile = fileName +'.rrparams' 1626 if os.path.isfile(rrparamsFile): 1627 self.RrParams = io.rrRead(self.IonStr)

1628 # 1629 # ------------------------------------------------------------------------- 1630 #

1631 - def spectrum(self, wavelength, filter=(chfilters.gaussianR,1000.), label=0, allLines=1, em=0):

1632 ''' 1633 Calculates the line emission spectrum for the specified ion. 1634 1635 Convolves the results of intensity to make them look like an observed spectrum 1636 the default filter is the gaussianR filter with a resolving power of 1000. Other choices 1637 include chianti.filters.box and chianti.filters.gaussian. When using the box filter, 1638 the width should equal the wavelength interval to keep the units of the continuum and line 1639 spectrum the same. 1640 1641 includes ionization equilibrium and elemental abundances 1642 1643 can be called multiple times to use different filters and widths 1644 uses label to keep the separate applications of spectrum sorted by the label 1645 for example, do .spectrum( .... labe='test1') 1646 and do .spectrum( .... label = 'test2') 1647 then will get self.Spectrum.keys() = test1, test2 and 1648 self.Spectrum['test1'] = {'intensity':aspectrum, 'wvl':wavelength, 'filter':useFilter.__name__, 'filterWidth':useFactor} 1649 1650 Note: scipy.ndimage.filters also includes a range of filters. 1651 ''' 1652 # aspectrum = np.zeros_like(wavelength) 1653 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ ($\int\,$ N$_e\,$N$_H\,$d${\it l}$)$^{-1}$' 1654 # 1655 xlabel = 'Wavelength ('+self.Defaults['wavelength'] +')' 1656 # 1657 nTemp = self.Temperature.size 1658 nDens = self.EDensity.size 1659 useFilter = filter[0] 1660 useFactor= filter[1] 1661 # 1662 wvlRange = [wavelength.min(), wavelength.max()] 1663 if hasattr(self, 'Intensity'): 1664 intensity = self.Intensity 1665 else: 1666 self.intensity(wvlRange=wvlRange, allLines=allLines, em=em) 1667 intensity = self.Intensity 1668 # 1669 # if intensity had been called with em, then the intensities are already multiply by em 1670 if hasattr(self, 'Em'): 1671 em = self.Em 1672 useEm = 0 1673 elif type(em) == int and em == 0: 1674 em = np.ones(self.NTempDen, 'float64') 1675 self.Em = em 1676 useEm = 1 1677 elif type(em) == float and em > 0.: 1678 em = np.ones(self.NTempDen, 'float64')*em 1679 self.Em = em 1680 useEm = 1 1681 elif type(em) == list or type(em) == tuple or type(em) == np.ndarray: 1682 em = np.asarray(em, 'float64') 1683 self.Em = em 1684 useEm = 1 1685 # 1686 if self.Em.any() > 0.: 1687 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ ' 1688 else: 1689 ylabel = r'erg cm$^{-2}$ s$^{-1}$ sr$^{-1} \AA^{-1}$ ($\int\,$ N$_e\,$N$_H\,$d${\it l}$)$^{-1}$' 1690 # 1691 xlabel = 'Wavelength ('+self.Defaults['wavelength'] +')' 1692 # 1693 if self.NTempDen == 1: 1694 aspectrum = np.zeros_like(wavelength) 1695 if not 'errorMessage' in self.Intensity.keys(): 1696 idx = util.between(self.Intensity['wvl'], wvlRange) 1697 if len(idx) == 0: 1698 print(' no lines in wavelength range %12.2f - %12.2f'%(wavelength.min(), wavelength.max())) 1699 self.Spectrum = {'errorMessage':' no lines in wavelength range %12.2f - %12.2f'%(wavelength.min(), wavelength.max())} 1700 return 1701 for iwvl in idx: 1702 wvlCalc = self.Intensity['wvl'][iwvl] 1703 aspectrum += useFilter(wavelength, wvlCalc, factor=useFactor)*intensity['intensity'][iwvl] 1704 if useEm: 1705 aspectrum *= em 1706 elif self.NTempDen > 1: 1707 aspectrum = np.zeros((self.NTempDen, wavelength.size), 'float64') 1708 if not 'errorMessage' in self.Intensity.keys(): 1709 idx = util.between(self.Intensity['wvl'], wvlRange) 1710 if len(idx) == 0: 1711 print(' no lines in wavelength range %12.2f - %12.2f'%(wavelength.min(), wavelength.max())) 1712 self.Spectrum = {'errorMessage':' no lines in wavelength range %12.2f - %12.2f'%(wavelength.min(), wavelength.max())} 1713 return 1714 for itemp in range(self.NTempDen): 1715 for iwvl in idx: 1716 wvlCalc = self.Intensity['wvl'][iwvl] 1717 aspectrum[itemp] += useFilter(wavelength, wvlCalc, factor=useFactor)*self.Intensity['intensity'][itemp, iwvl] 1718 if useEm: 1719 aspectrum[itemp] *= em[itemp] 1720 # aspectrum[itemp] += useFilter(wavelength, wvlCalc, factor=useFactor)*self.Intensity['intensity'][iwvl, itemp] 1721 # for iwvl, wvlCalc in enumerate(self.Intensity['wvl']): 1722 # aspectrum[itemp] += useFilter(wavelength, wvlCalc, factor=useFactor)*self.Intensity['intensity'][itemp, iwvl] 1723 if type(label) == type(''): 1724 if hasattr(self, 'Spectrum'): 1725 self.Spectrum[label] = {'intensity':aspectrum, 'wvl':wavelength, 'filter':useFilter.__name__, 'filterWidth':useFactor, 'allLines':allLines, 'em':em, 'xlabel':xlabel, 'ylabel':ylabel} 1726 else: 1727 self.Spectrum = {label:{'intensity':aspectrum, 'wvl':wavelength, 'filter':useFilter.__name__, 'filterWidth':useFactor, 'allLines':allLines, 'em':em, 'xlabel':xlabel, 'ylabel':ylabel}} 1728 1729 else: 1730 self.Spectrum = {'intensity':aspectrum, 'wvl':wavelength, 'filter':useFilter.__name__, 'filterWidth':useFactor, 'allLines':allLines, 'em':em, 'xlabel':xlabel, 'ylabel':ylabel}

1731 # 1732 # ------------------------------------------------------------------------------------- 1733 #

1734 - def populate(self, popCorrect=1, verbose=0, **kwargs):

1735 """ 1736 Calculate level populations for specified ion. 1737 possible keyword arguments include temperature, eDensity, pDensity, radTemperature and rStar 1738 """ 1739 # 1740 # 1741 for one in kwargs.keys(): 1742 if one not in chdata.keywordArgs: 1743 print(' following keyword is not understood - %20s '%(one)) 1744 # 1745 nlvls=self.Nlvls 1746 nwgfa=self.Nwgfa 1747 nscups=self.Nscups 1748 npsplups=self.Npsplups 1749 # 1750 if 'temperature' in kwargs.keys(): 1751 self.Temperature = np.asarray(kwargs['temperature']) 1752 temperature = self.Temperature 1753 elif hasattr(self, 'Temperature'): 1754 temperature=self.Temperature 1755 else: 1756 print(' no temperature values have been set') 1757 return 1758 # 1759 if 'eDensity' in kwargs.keys(): 1760 self.EDensity = np.asarray(kwargs['eDensity']) 1761 eDensity = self.EDensity 1762 elif hasattr(self, 'EDensity'): 1763 eDensity = self.EDensity 1764 else: 1765 print(' no eDensity values have been set') 1766 return 1767 # 1768 if 'pDensity' in kwargs.keys(): 1769 if kwargs['pDensity'] == 'default': 1770 self.p2eRatio() 1771 protonDensity = self.ProtonDensityRatio*self.EDensity 1772 else: 1773 try: 1774 self.PDensity = np.asarray(kwargs['pDensity']) 1775 except: 1776 print(' could not interpret value for keyword pDensity') 1777 print(' should be either "default" or a number or array') 1778 return 1779 else: 1780 if hasattr(self, 'PDensity'): 1781 protonDensity = self.PDensity 1782 else: 1783 self.p2eRatio() 1784 self.PDensity = self.ProtonDensityRatio*self.EDensity 1785 protonDensity = self.PDensity 1786 print(' proton density not specified, set to \"default\" ') 1787 # 1788 if 'radTemperature' in kwargs.keys() and 'rStar' in kwargs.keys(): 1789 self.RadTemperature = np.asarray(kwargs['radTemperature']) 1790 radTemperature = np.array(self.RadTemperature) 1791 self.RStar = np.asarray(kwargs['rStar']) 1792 rStar = np.asarray(self.RStar) 1793 elif hasattr(self, 'RadTemperature') and hasattr(self, 'RStar'): 1794 radTemperature = self.RadTemperature 1795 rStar = self.RStar 1796 # 1797 # 1798 # the Dielectronic test should eventually go away 1799 rec = 0 1800 ci = 0 1801 if popCorrect and (not self.Dielectronic): 1802 if self.Ncilvl: 1803 ci = 1 1804 cilvl = self.Cilvl 1805 if hasattr(self, 'CilvlRate'): 1806 cilvlRate = self.CilvlRate 1807 else: 1808 self.cireclvlDescale('cilvl') 1809 cilvlRate = self.CilvlRate 1810 self.recombRate() 1811 # 1812 lowers = util.zion2name(self.Z, self.Ion-1) 1813 # get the lower ionization stage 1814 lower = ion(lowers, temperature=self.Temperature, eDensity = self.EDensity) 1815 lower.ionizRate() 1816 # need to get multiplicity of lower ionization stage 1817 lowMult = lower.Elvlc['mult'] 1818 # else: 1819 # ci = 0 1820 # try: 1821 rec = 0 1822 if self.Nreclvl: 1823 rec = 1 1824 reclvl = self.Reclvl 1825 if hasattr(self, 'ReclvlRate'): 1826 reclvlRate = self.ReclvlRate 1827 else: 1828 self.cireclvlDescale('reclvl') 1829 reclvlRate = self.ReclvlRate 1830 # if self.Ndielsplups: 1831 # rec = 1 1832 # if hasattr(self, 'DielUpsilon'): 1833 # dielexRate = self.DielUpsilon['exRate'] 1834 # else: 1835 # print(' doing upsilonDescale') 1836 # self.upsilonDescale(diel=1) 1837 # dielexRate = self.DielUpsilon['exRate'] 1838 1839 # except: 1840 # self.Reclvl = io.cireclvlRead(self.IonStr,'reclvl' ) 1841 # reclvl = self.Reclvl 1842 # self.reclvlDescale() 1843 # if type(self.Reclvl) != type(None): 1844 # rec = 1 1845 # self.ionizRate() 1846 # # get the higher ionization stage 1847 # highers = util.zion2name(self.Z, self.Ion+1) 1848 ## print ' highers = ', highers 1849 # higher = ion(highers, temperature=self.Temperature, eDensity=self.EDensity) 1850 # higher.recombRate() 1851 # else: 1852 # rec = 0 1853 # 1854 # elif self.Dielectronic: 1855 ## self.Ndielsplups and self.Dielectronic: 1856 # rec = 0 1857 # if hasattr(self, 'DielUpsilon'): 1858 # dielexRate = self.DielUpsilon['exRate'] 1859 # else: 1860 ## print(' doing upsilonDescale') 1861 # self.upsilonDescaleSplups(diel=1) 1862 # dielexRate = self.DielUpsilon['exRate'] 1863 # 1864 if rec: 1865 # if self.Ndielsplups: 1866 # self.upsilonDescale(diel=1) 1867 # dielexRate = self.DielUpsilon['exRate'] 1868 # get ionization rate of this iion 1869 self.ionizRate() 1870 # get the higher ionization stage 1871 highers = util.zion2name(self.Z, self.Ion+1) 1872 higher = ion(highers, temperature=self.Temperature, eDensity=self.EDensity) 1873 higher.recombRate() 1874 # print ' nlvls, ci, rec = ', nlvls, ci, rec 1875 # 1876 rad=np.zeros((nlvls+ci+rec,nlvls+ci+rec),"float64") # the populating matrix for radiative transitions 1877 # 1878 # 1879 for iwgfa in range(nwgfa): 1880 l1 = self.Wgfa["lvl1"][iwgfa]-1 1881 l2 = self.Wgfa["lvl2"][iwgfa]-1 1882 rad[l1+ci,l2+ci] += self.Wgfa["avalue"][iwgfa] 1883 rad[l2+ci,l2+ci] -= self.Wgfa["avalue"][iwgfa] 1884 # photo-excitation and stimulated emission 1885 if self.RadTemperature: 1886 if not self.RStar: 1887 dilute = 0.5 1888 else: 1889 dilute = util.dilute(self.RStar) 1890 # next - don't include autoionization lines 1891 if abs(self.Wgfa['wvl'][iwgfa]) > 0.: 1892 de = const.invCm2Erg*(self.Elvlc['ecm'][l2] - self.Elvlc['ecm'][l1]) 1893 dekt = de/(const.boltzmann*self.RadTemperature) 1894 # photoexcitation 1895 phexFactor = dilute*(float(self.Elvlc['mult'][l2])/float(self.Elvlc['mult'][l1]))/(np.exp(dekt) -1.) 1896 rad[l2+ci,l1+ci] += self.Wgfa["avalue"][iwgfa]*phexFactor 1897 rad[l1+ci,l1+ci] -= self.Wgfa["avalue"][iwgfa]*phexFactor 1898 # stimulated emission 1899 stemFactor = dilute/(np.exp(-dekt) -1.) 1900 rad[l1+ci,l2+ci] += self.Wgfa["avalue"][iwgfa]*stemFactor 1901 rad[l2+ci,l2+ci] -= self.Wgfa["avalue"][iwgfa]*stemFactor 1902 1903 # 1904 # 1905 if self.Nscups: 1906 self.upsilonDescale(diel=self.Dielectronic) 1907 #ups = self.Upsilon['upsilon'] 1908 exRate = self.Upsilon['exRate'] 1909 dexRate = self.Upsilon['dexRate'] 1910 # 1911 if npsplups: 1912 self.upsilonDescaleSplups(prot=1) 1913 # pups = self.PUpsilon['upsilon'] 1914 pexRate = self.PUpsilon['exRate'] 1915 pdexRate = self.PUpsilon['dexRate'] 1916 1917 # 1918 temp = temperature 1919 ntemp = temp.size 1920 dens = self.EDensity 1921 ndens = dens.size 1922 cc=const.collision*self.EDensity 1923 # 1924 if npsplups: 1925 cp=const.collision*protonDensity 1926 if ntemp > 1 and ndens >1 and ntemp != ndens: 1927 print(' unless temperature or eDensity are single values') 1928 print(' the number of temperatures values must match the ') 1929 print(' the number of eDensity values') 1930 return 1931 # 1932 # get corrections for recombination and excitation 1933 # 1934 nscups = self.Nscups 1935 # 1936 # 1937 # the way temperature and density are now (9/2015) handled as arrays of the same size 1938 # one the ndens == ntemp =1 case and the ndens >1 and ntemp>1 case are really needed 1939 # 1940 # first, for ntemp=ndens=1 1941 if ndens == 1 and ntemp == 1: 1942 popmat=np.copy(rad) 1943 if not self.Dielectronic: 1944 for iscups in range(0,nscups): 1945 l1=self.Scups["lvl1"][iscups]-1 1946 l2=self.Scups["lvl2"][iscups]-1 1947 # 1948 popmat[l1+ci,l2+ci] += self.EDensity*dexRate[iscups] 1949 popmat[l2+ci,l1+ci] += self.EDensity*exRate[iscups] 1950 popmat[l1+ci,l1+ci] -= self.EDensity*exRate[iscups] 1951 popmat[l2+ci,l2+ci] -= self.EDensity*dexRate[iscups] 1952 # 1953 for isplups in range(0,npsplups): 1954 l1=self.Psplups["lvl1"][isplups]-1 1955 l2=self.Psplups["lvl2"][isplups]-1 1956 # 1957 popmat[l1+ci,l2+ci] += self.PDensity*pdexRate[isplups] 1958 popmat[l2+ci,l1+ci] += self.PDensity*pexRate[isplups] 1959 popmat[l1+ci,l1+ci] -= self.PDensity*pexRate[isplups] 1960 popmat[l2+ci,l2+ci] -= self.PDensity*pdexRate[isplups] 1961 # now include ionization rate from 1962 if ci: 1963 # print ' ci = ', ci 1964 # 1965 # the ciRate can be computed for all temperatures 1966 # 1967 ciTot = 0. 1968 for itrans in range(len(cilvl['lvl1'])): 1969 lvl1 = cilvl['lvl1'][itrans]-1 1970 lvl2 = cilvl['lvl2'][itrans]-1 1971 # de = cilvl['de'][itrans] 1972 # ekt = (de*1.57888e+5)/temperature 1973 # mult = lowMult[lvl1-1] 1974 # this is kind of double booking the ionization rate components 1975 popmat[lvl2+ci, lvl1] += self.EDensity*self.CilvlRate['rate'][itrans] 1976 popmat[lvl1, lvl1] -= self.EDensity*self.CilvlRate['rate'][itrans] 1977 ciTot += self.EDensity*self.CilvlRate['rate'][itrans] 1978 # 1979 popmat[1, 0] += (self.EDensity*lower.IonizRate['rate'] - ciTot) 1980 popmat[0, 0] -= (self.EDensity*lower.IonizRate['rate'] - ciTot) 1981 popmat[0, 1] += self.EDensity*self.RecombRate['rate'] 1982 popmat[1, 1] -= self.EDensity*self.RecombRate['rate'] 1983 if rec: 1984 # 1985 # print ' rec, dielTot = ', rec, dielTot 1986 # 1987 for itrans in range(self.Nreclvl): 1988 # lvl1 = reclvl['lvl1'][itrans]-1 1989 lvl2 = reclvl['lvl2'][itrans]-1 1990 popmat[lvl2+ci, -1] += self.EDensity*reclvlRate['rate'][itrans] 1991 popmat[-1, -1] -= self.EDensity*reclvlRate['rate'][itrans] 1992 if self.Nreclvl: 1993 recTot = reclvlRate['rate'].sum(axis=0) 1994 else: 1995 recTot = 0. 1996 1997 # 1998 popmat[-1, ci] += self.EDensity*self.IonizRate['rate'] 1999 popmat[ci, ci] -= self.EDensity*self.IonizRate['rate'] 2000 # next 2 line take care of overbooking 2001 popmat[ci, -1] += self.EDensity*(higher.RecombRate['rate']- recTot) 2002 popmat[-1, -1] -= self.EDensity*(higher.RecombRate['rate']- recTot) 2003 # 2004 if self.Dielectronic: 2005 # dielTot = 0. 2006 for iscups in range(0,nscups): 2007 l1=self.Scups["lvl1"][iscups]-1 2008 l2=self.Scups["lvl2"][iscups]-1 2009 # 2010 popmat[l2+ci,-1] += self.EDensity*exRate[iscups] 2011 popmat[-1, -1] -= self.EDensity*exRate[iscups] 2012 # for dielectronic ions, l1 = ground level of the ion itself 2013 # branch = np.zeros(self.Ndielsplups, 'float64') 2014 # for isplups in range(0,self.Ndielsplups): 2015 # l1 = self.DielSplups["lvl1"][isplups]-1 + nlvls 2016 # l2 = self.DielSplups["lvl2"][isplups]-1 2017 # auto = rad[l1+ci, l2+ci] 2018 # avalue = rad[:, l2+ci] 2019 # good = avalue > 0. 2020 # avalueTot = avalue[good].sum() 2021 # branch[isplups] = (avalueTot-auto)/avalueTot 2022 ## print ' l1 %4i l2 %4i auto %10.2e avalue %10.2e tot %10.3f'%( l1, l2, auto, avalueTot, branch[isplups]) 2023 # self.DielUpsilon['branch'] = branch 2024 # # 2025 # dielTot = 0. 2026 # print ' Ndielsplups > 0 ' 2027 # for isplups in range(0,self.Ndielsplups): 2028 # l1 = self.DielSplups["lvl1"][isplups]-1 2029 # l2 = self.DielSplups["lvl2"][isplups]-1 2030 # # 2031 ## print ' l1, l2, dielexRate = ', l1, l2, dielexRate[isplups] 2032 # popmat[l2+ci,l1+ci] += self.EDensity*dielexRate[isplups] 2033 # popmat[l1+ci,l1+ci] -= self.EDensity*dielexRate[isplups] 2034 # 2035 # dielTot += self.EDensity*exRate[isplups]*branch[isplups] 2036 # else: 2037 # dielTot = 0. 2038 2039 # print ' higher, rec , dieltot = ', self.EDensity*higher.RecombRate['rate'], self.EDensity*reclvlRate['rate'].sum(axis=0), dielTot 2040 # normalize to unity 2041 # print(' rec = %5i ci = %5i'%(rec, ci)) 2042 norm=np.ones(nlvls+ci+rec,'float64') 2043 if ci: 2044 norm[0] = 0. 2045 if rec: 2046 norm[nlvls+ci+rec-1] = 0. 2047 if self.Dielectronic: 2048 norm[nlvls-1] = 0. 2049 popmat[nlvls+ci+rec-1]=norm 2050 # popmata = np.copy(popmat) 2051 # popmata[nlvls+ci+rec-1]=norm 2052 #popmata[nlvls+ci+rec-1]=norm 2053 b=np.zeros(nlvls+ci+rec,'float64') 2054 b[nlvls+ci+rec-1]=1. 2055 # print ' norm = ', norm 2056 # print 'popmat, last line', popmat[-1] 2057 # print ' b = ', b 2058 # popmat[nlvls/2]=norm 2059 # b=np.zeros(nlvls+ci+rec,'float64') 2060 # b[nlvls/2]=1. 2061 # if rec: 2062 # fullpop = np.linalg.solve(popmat,b) 2063 # pop = fullpop[ci:ci+nlvls+rec-1] 2064 # else: 2065 # fullpop = np.linalg.solve(popmat,b) 2066 # pop = fullpop[ci:] 2067 # fullpop = np.linalg.solve(popmat,b) 2068 try: 2069 fullpop=np.linalg.solve(popmat,b) 2070 pop = fullpop[ci:ci+nlvls] 2071 except np.linalg.LinAlgError: 2072 pop = np.zeros(nlvls, 'float64') 2073 # print ' error in matrix inversion, setting populations to zero at T = ', ('%8.2e')%(temperature) 2074 # 2075 # ---------------------------------------------------------------------------------- 2076 # next, in case of a single eDensity value 2077 # pop = np.linalg.solve(popmat,b) 2078 elif ndens == 1: 2079 # print(' doing ndens = 1 size of EDensity = %5i'%(self.EDensity.size)) 2080 pop=np.zeros((ntemp, nlvls),"float64") 2081 # pop=np.zeros((ntemp,ci + nlvls + rec),"float64") 2082 for itemp in range(0,ntemp): 2083 popmat=np.copy(rad) 2084 for iscups in range(0,nscups): 2085 l1=self.Scups["lvl1"][iscups]-1 2086 l2=self.Scups["lvl2"][iscups]-1 2087 popmat[l1+ci,l2+ci] += self.EDensity*dexRate[iscups, itemp] 2088 popmat[l2+ci,l1+ci] += self.EDensity*exRate[iscups, itemp] 2089 popmat[l1+ci,l1+ci] -= self.EDensity*exRate[iscups, itemp] 2090 popmat[l2+ci,l2+ci] -= self.EDensity*dexRate[iscups, itemp] 2091 for isplups in range(0,npsplups): 2092 l1=self.Psplups["lvl1"][isplups]-1 2093 l2=self.Psplups["lvl2"][isplups]-1 2094 # for proton excitation, the levels are all below the ionization potential 2095 # 2096 popmat[l1+ci,l2+ci] += self.PDensity[itemp]*pdexRate[isplups, itemp] 2097 popmat[l2+ci,l1+ci] += self.PDensity[itemp]*pexRate[isplups, itemp] 2098 popmat[l1+ci,l1+ci] -= self.PDensity[itemp]*pexRate[isplups, itemp] 2099 popmat[l2+ci,l2+ci] -= self.PDensity[itemp]*pdexRate[isplups, itemp] 2100 # now include ionization rate from 2101 if ci: 2102 # print ' ci = ', ci 2103 # 2104 # the ciRate can be computed for all temperatures 2105 # 2106 ciTot = 0. 2107 for itrans in range(len(cilvl['lvl1'])): 2108 lvl1 = cilvl['lvl1'][itrans]-1 2109 lvl2 = cilvl['lvl2'][itrans]-1 2110 # de = cilvl['de'][itrans] 2111 # ekt = (de*1.57888e+5)/temperature 2112 #mult = lowMult[lvl1-1] 2113 # this is kind of double booking the ionization rate components 2114 popmat[lvl2+ci, lvl1] += self.EDensity*self.CilvlRate['rate'][itrans, itemp] 2115 popmat[lvl1, lvl1] -= self.EDensity*self.CilvlRate['rate'][itrans, itemp] 2116 ciTot += self.EDensity*self.CilvlRate['rate'][itrans, itemp] 2117 # popmat[lvl2, lvl1-1] += self.EDensity*cirate[itemp] 2118 # popmat[lvl1-1, lvl1-1] -= self.EDensity*cirate[itemp] 2119 popmat[1, 0] += (self.EDensity*lower.IonizRate['rate'][itemp] - ciTot) 2120 popmat[0, 0] -= (self.EDensity*lower.IonizRate['rate'][itemp] - ciTot) 2121 popmat[0, 1] += self.EDensity*self.RecombRate['rate'][itemp] 2122 popmat[1, 1] -= self.EDensity*self.RecombRate['rate'][itemp] 2123 if rec: 2124 # 2125 if self.Nreclvl: 2126 recTot = self.ReclvlRate['rate'][:, itemp].sum() 2127 else: 2128 recTot = 0. 2129 # 2130 popmat[-1, ci] += self.EDensity*self.IonizRate['rate'][itemp] 2131 popmat[ci, ci] -= self.EDensity*self.IonizRate['rate'][itemp] 2132 popmat[ci, -1] += self.EDensity*(higher.RecombRate['rate'][itemp]- recTot) 2133 popmat[-1, -1] -= self.EDensity*(higher.RecombRate['rate'][itemp]- recTot) 2134 # popmat[ci, -1] += self.EDensity*(higher.RecombRate['rate'][itemp]- self.ReclvlRate['rate'][:, itemp].sum()) - dielTot 2135 # popmat[-1, -1] -= self.EDensity*(higher.RecombRate['rate'][itemp]- self.ReclvlRate['rate'][:, itemp].sum()) + dielTot 2136 # popmat[ci, -1] += self.EDensity*higher.RecombRate['rate'][itemp] 2137 # popmat[-1, -1] -= self.EDensity*higher.RecombRate['rate'][itemp] 2138 # 2139 # for itrans in range(len(reclvl['lvl1'])): 2140 for itrans in range(self.Nreclvl): 2141 lvl1 = reclvl['lvl1'][itrans]-1 2142 lvl2 = reclvl['lvl2'][itrans]-1 2143 popmat[lvl2+ci, -1] += self.EDensity*self.ReclvlRate['rate'][itrans, itemp] 2144 popmat[-1, -1] -= self.EDensity*self.ReclvlRate['rate'][itrans, itemp] 2145 # 2146 # if self.Dielectronic: 2147 # branch = np.zeros(self.Ndielsplups, 'float64') 2148 # for isplups in range(0,self.Ndielsplups): 2149 # l1 = self.DielSplups["lvl1"][isplups]-1 + nlvls 2150 # l2 = self.DielSplups["lvl2"][isplups]-1 2151 # auto = rad[l1+ci, l2+ci] 2152 # avalue = rad[:, l2+ci] 2153 # good = avalue > 0. 2154 # avalueTot = avalue[good].sum() 2155 # branch[isplups] = (avalueTot-auto)/avalueTot 2156 ## print ' l1 %4i l2 %4i auto %10.2e avalue %10.2e tot %10.3f'%( l1, l2, auto, avalueTot, branch[isplups]) 2157 # self.DielUpsilon['branch'] = branch 2158 # # 2159 # dielTot = 0. 2160 # print ' Ndielsplups > 0 ' 2161 # for isplups in range(0,self.Ndielsplups): 2162 # l1 = self.DielSplups["lvl1"][isplups]-1 2163 # l2 = self.DielSplups["lvl2"][isplups]-1 2164 # # 2165 # popmat[l2+ci,l1+ci] += self.EDensity*dielexRate[isplups, itemp] 2166 # popmat[l1+ci,l1+ci] -= self.EDensity*dielexRate[isplups, itemp] 2167 # 2168 # dielTot += self.EDensity*dielexRate[isplups, itemp]*branch[isplups] 2169 # else: 2170 # dielTot = 0. 2171 # normalize to unity 2172 norm=np.ones(nlvls+ci+rec,'float64') 2173 if ci: 2174 norm[0] = 0. 2175 if rec: 2176 norm[-1] = 0. 2177 if self.Dielectronic: 2178 norm[-1] = 0. 2179 popmat[nlvls+ci+rec-1]=norm 2180 b=np.zeros(nlvls+ci+rec,'float64') 2181 b[nlvls+ci+rec-1]=1. 2182 try: 2183 thispop=np.linalg.solve(popmat,b) 2184 pop[itemp] = thispop[ci:ci+nlvls] 2185 except np.linalg.LinAlgError: 2186 pop[itemp] = np.zeros(nlvls, 'float64') 2187 # print ' error in matrix inversion, setting populations to zero at T = ', ('%8.2e')%(temperature[itemp]) 2188 # 2189 elif ntemp == 1: 2190 # pop=np.zeros((ndens,nlvls),"float64") 2191 # print(' doing ntemp = 1 size of Temperature = %5i'%(self.Temperature.size)) 2192 pop=np.zeros((ndens,nlvls),"float64") 2193 for idens in range(0,ndens): 2194 popmat=np.copy(rad) 2195 for isplups in range(0,nscups): 2196 l1=self.Scups["lvl1"][isplups]-1 2197 l2=self.Scups["lvl2"][isplups]-1 2198 # if self.Dielectronic: 2199 # de=np.abs((self.Elvlc["eryd"][l2]-self.Ip/const.ryd2Ev)-self.Elvlc["eryd"][l1]) 2200 # else: 2201 # de=np.abs(self.Elvlc["eryd"][l2]-self.Elvlc["eryd"][l1]) 2202 # ekt=(de*1.57888e+5)/temp 2203 # fmult1=float(self.Elvlc["mult"][l1]) 2204 # fmult2=float(self.Elvlc["mult"][l2]) 2205 # popmat[l1+ci,l2+ci]+=cc[idens]*ups[isplups]/(fmult2*np.sqrt(temp)) 2206 # popmat[l2+ci,l1+ci]+=cc[idens]*ups[isplups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2207 # popmat[l1+ci,l1+ci]-=cc[idens]*ups[isplups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2208 # popmat[l2+ci,l2+ci]-=cc[idens]*ups[isplups]/(fmult2*np.sqrt(temp)) 2209 # 2210 popmat[l1+ci,l2+ci] += self.EDensity[idens]*dexRate[isplups] 2211 popmat[l2+ci,l1+ci] += self.EDensity[idens]*exRate[isplups] 2212 popmat[l1+ci,l1+ci] -= self.EDensity[idens]*exRate[isplups] 2213 popmat[l2+ci,l2+ci] -= self.EDensity[idens]*dexRate[isplups] 2214 # 2215 for isplups in range(0,npsplups): 2216 l1=self.Psplups["lvl1"][isplups]-1 2217 l2=self.Psplups["lvl2"][isplups]-1 2218 # # for proton excitation, the levels are all below the ionization potential 2219 # de=np.abs(self.Elvlc["eryd"][l2]-self.Elvlc["eryd"][l1]) 2220 # ekt=(de*1.57888e+5)/temp 2221 # fmult1=float(self.Elvlc["mult"][l1]) 2222 # fmult2=float(self.Elvlc["mult"][l2]) 2223 # popmat[l1+ci,l2+ci]+=cp[idens]*pups[isplups]/(fmult2*np.sqrt(temp)) 2224 # popmat[l2+ci,l1+ci]+=cp[idens]*pups[isplups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2225 # popmat[l1+ci,l1+ci]-=cp[idens]*pups[isplups]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2226 # popmat[l2+ci,l2+ci]-=cp[idens]*pups[isplups]/(fmult2*np.sqrt(temp)) 2227 # 2228 popmat[l1+ci,l2+ci] += self.PDensity[idens]*pdexRate[isplups] 2229 popmat[l2+ci,l1+ci] += self.PDensity[idens]*pexRate[isplups] 2230 popmat[l1+ci,l1+ci] -= self.PDensity[idens]*pexRate[isplups] 2231 popmat[l2+ci,l2+ci] -= self.PDensity[idens]*pdexRate[isplups] 2232 # now include ionization rate from 2233 if ci: 2234 # print ' ci = ', ci 2235 # 2236 # 2237 ciTot = 0. 2238 for itrans in range(len(cilvl['lvl1'])): 2239 lvl1 = cilvl['lvl1'][itrans] -1 2240 lvl2 = cilvl['lvl2'][itrans] -1 2241 # de = cilvl['de'][itrans] 2242 # ekt = (de*1.57888e+5)/temperature 2243 # this is kind of double booking the ionization rate components 2244 # popmat[lvl2, lvl1-1] += self.EDensity[idens]*cirate 2245 # popmat[lvl1-1, lvl1-1] -= self.EDensity[idens]*cirate 2246 popmat[lvl2+ci, lvl1] += self.EDensity[idens]*self.CilvlRate['rate'][itrans] 2247 popmat[lvl1, lvl1] -= self.EDensity[idens]*self.CilvlRate['rate'][itrans] 2248 ciTot += self.EDensity[idens]*self.CilvlRate['rate'][itrans] 2249 popmat[1, 0] += (self.EDensity[idens]*lower.IonizRate['rate'] -ciTot) 2250 popmat[0, 0] -= (self.EDensity[idens]*lower.IonizRate['rate'] -ciTot) 2251 popmat[0, 1] += self.EDensity[idens]*self.RecombRate['rate'] 2252 popmat[1, 1] -= self.EDensity[idens]*self.RecombRate['rate'] 2253 if rec: 2254 # print ' rec = ', rec 2255 if self.Nreclvl: 2256 # print ' ReclvlRate.shape = ', self.ReclvlRate['rate'].shape 2257 recTot = self.ReclvlRate['rate'].sum() 2258 else: 2259 recTot = 0. 2260 # 2261 popmat[-1, ci] += self.EDensity[idens]*self.IonizRate['rate'] 2262 popmat[ci, ci] -= self.EDensity[idens]*self.IonizRate['rate'] 2263 popmat[ci, -1] += self.EDensity[idens]*(higher.RecombRate['rate'] - recTot) 2264 popmat[-1, -1] -= self.EDensity[idens]*(higher.RecombRate['rate'] - recTot) 2265 # popmat[ci, -1] += self.EDensity[idens]*higher.RecombRate['rate'] 2266 # popmat[-1, -1] -= self.EDensity[idens]*higher.RecombRate['rate'] 2267 # 2268 # for itrans in range(len(reclvl['lvl1'])): 2269 for itrans in range(self.Nreclvl): 2270 lvl1 = reclvl['lvl1'][itrans]-1 2271 lvl2 = reclvl['lvl2'][itrans]-1 2272 popmat[lvl2+ci, -1] += self.EDensity[idens]*self.ReclvlRate['rate'][itrans] 2273 popmat[-1, -1] -= self.EDensity[idens]*self.ReclvlRate['rate'][itrans] 2274 # 2275 # if self.Dielectronic: 2276 # branch = np.zeros(self.Ndielsplups, 'float64') 2277 # for isplups in range(0,self.Ndielsplups): 2278 # l1 = self.DielSplups["lvl1"][isplups]-1 + nlvls 2279 # l2 = self.DielSplups["lvl2"][isplups]-1 2280 # auto = rad[l1+ci, l2+ci] 2281 # avalue = rad[:, l2+ci] 2282 # good = avalue > 0. 2283 # avalueTot = avalue[good].sum() 2284 # branch[isplups] = (avalueTot-auto)/avalueTot 2285 ## print ' l1 %4i l2 %4i auto %10.2e avalue %10.2e tot %10.3f'%( l1, l2, auto, avalueTot, branch[isplups]) 2286 # self.DielUpsilon['branch'] = branch 2287 # # 2288 # dielTot = 0. 2289 ## print ' Ndielsplups > 0 ' 2290 # for isplups in range(0,self.Ndielsplups): 2291 # l1 = self.DielSplups["lvl1"][isplups]-1 2292 # l2 = self.DielSplups["lvl2"][isplups]-1 2293 # # 2294 # popmat[l2+ci,l1+ci] += self.EDensity[idens]*dielexRate[isplups] 2295 # popmat[l1+ci,l1+ci] -= self.EDensity[idens]*dielexRate[isplups] 2296 # 2297 # dielTot += self.EDensity[idens]*dielexRate[isplups]*branch[isplups] 2298 # else: 2299 # dielTot = 0. 2300 2301 # normalize to unity 2302 norm=np.ones(nlvls+ci+rec,'float64') 2303 if ci: 2304 norm[0] = 0. 2305 if rec: 2306 norm[-1] = 0. 2307 if self.Dielectronic: 2308 norm[-1] = 0. 2309 popmat[nlvls+ci+rec-1]=norm 2310 b=np.zeros(nlvls+ci+rec,'float64') 2311 b[nlvls+ci+rec-1]=1. 2312 try: 2313 thispop=np.linalg.solve(popmat,b) 2314 pop[idens] = thispop[ci:ci+nlvls] 2315 except np.linalg.LinAlgError: 2316 pop[idens] = np.zeros(nlvls, 'float64') 2317 # print ' error in matrix inversion, setting populations to zero at eDensity = ', ('%8.2e')%(eDensity[idens]) 2318 # thispop=np.linalg.solve(popmat,b) 2319 # if rec: 2320 # pop[idens] = thispop[ci:ci+nlvls+rec-1] 2321 # else: 2322 # pop[idens] = thispop[ci:] 2323 # pop[idens] = thispop[ci:ci+nlvls] 2324 # 2325 elif ntemp>1 and ntemp==ndens: 2326 # print(' doing ntemp > 1 and ndens > 1') 2327 pop=np.zeros((ntemp,nlvls),"float64") 2328 # pop=np.zeros((ntemp,ci+nlvls+rec),"float64") 2329 for itemp in range(0,ntemp): 2330 temp=self.Temperature[itemp] 2331 popmat=np.copy(rad) 2332 for isplups in range(0,nscups): 2333 l1=self.Scups["lvl1"][isplups]-1 2334 l2=self.Scups["lvl2"][isplups]-1 2335 # if self.Dielectronic: 2336 # de=np.abs((self.Elvlc["eryd"][l2]-self.Ip/const.ryd2Ev)-self.Elvlc["eryd"][l1]) 2337 # else: 2338 # de=np.abs(self.Elvlc["eryd"][l2]-self.Elvlc["eryd"][l1]) 2339 # ekt=(de*1.57888e+5)/temp 2340 # fmult1=float(self.Elvlc["mult"][l1]) 2341 # fmult2=float(self.Elvlc["mult"][l2]) 2342 # popmat[l1+ci,l2+ci]+=cc[itemp]*ups[isplups,itemp]/(fmult2*np.sqrt(temp)) 2343 # popmat[l2+ci,l1+ci]+=cc[itemp]*ups[isplups,itemp]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2344 # popmat[l1+ci,l1+ci]-=cc[itemp]*ups[isplups,itemp]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2345 # popmat[l2+ci,l2+ci]-=cc[itemp]*ups[isplups,itemp]/(fmult2*np.sqrt(temp)) 2346 # 2347 popmat[l1+ci,l2+ci] += self.EDensity[itemp]*dexRate[isplups, itemp] 2348 popmat[l2+ci,l1+ci] += self.EDensity[itemp]*exRate[isplups, itemp] 2349 popmat[l1+ci,l1+ci] -= self.EDensity[itemp]*exRate[isplups, itemp] 2350 popmat[l2+ci,l2+ci] -= self.EDensity[itemp]*dexRate[isplups, itemp] 2351 # proton rates 2352 for isplups in range(0,npsplups): 2353 l1=self.Psplups["lvl1"][isplups]-1 2354 l2=self.Psplups["lvl2"][isplups]-1 2355 # for proton excitation, the levels are all below the ionization potential 2356 # de=np.abs(self.Elvlc["eryd"][l2]-self.Elvlc["eryd"][l1]) 2357 # ekt=(de*1.57888e+5)/temp 2358 # fmult1=float(self.Elvlc["mult"][l1]) 2359 # fmult2=float(self.Elvlc["mult"][l2]) 2360 # popmat[l1+ci,l2+ci]+=cp[itemp]*pups[isplups,itemp]/(fmult2*np.sqrt(temp)) 2361 # popmat[l2+ci,l1+ci]+=cp[itemp]*pups[isplups,itemp]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2362 # popmat[l1+ci,l1+ci]-=cp[itemp]*pups[isplups,itemp]*np.exp(-ekt)/(fmult1*np.sqrt(temp)) 2363 # popmat[l2+ci,l2+ci]-=cp[itemp]*pups[isplups,itemp]/(fmult2*np.sqrt(temp)) 2364 # 2365 popmat[l1+ci,l2+ci] += self.PDensity[itemp]*pdexRate[isplups, itemp] 2366 popmat[l2+ci,l1+ci] += self.PDensity[itemp]*pexRate[isplups, itemp] 2367 popmat[l1+ci,l1+ci] -= self.PDensity[itemp]*pexRate[isplups, itemp] 2368 popmat[l2+ci,l2+ci] -= self.PDensity[itemp]*pdexRate[isplups, itemp] 2369 # now include ionization rate from 2370 if ci: 2371 # print ' ci = ', ci 2372 # 2373 # the ciRate can be computed for all temperatures 2374 # 2375 ciTot = 0. 2376 for itrans in range(len(cilvl['lvl1'])): 2377 lvl1 = cilvl['lvl1'][itrans] -1 2378 lvl2 = cilvl['lvl2'][itrans] -1 2379 # this is kind of double booking the ionization rate components 2380 # popmat[lvl2, lvl1-1] += self.EDensity[itemp]*cirate[itemp] 2381 # popmat[lvl1-1, lvl1-1] -= self.EDensity[itemp]*cirate[itemp] 2382 popmat[lvl2+ci, lvl1] += self.EDensity[itemp]*self.CilvlRate['rate'][itrans, itemp] 2383 popmat[lvl1, lvl1] -= self.EDensity[itemp]*self.CilvlRate['rate'][itrans, itemp] 2384 ciTot += self.EDensity[itemp]*self.CilvlRate['rate'][itrans, itemp] 2385 popmat[1, 0] += (self.EDensity[itemp]*lower.IonizRate['rate'][itemp] - ciTot) 2386 popmat[0, 0] -= (self.EDensity[itemp]*lower.IonizRate['rate'][itemp] - ciTot) 2387 popmat[0, 1] += self.EDensity[itemp]*self.RecombRate['rate'][itemp] 2388 popmat[1, 1] -= self.EDensity[itemp]*self.RecombRate['rate'][itemp] 2389 if rec: 2390 # 2391 if self.Nreclvl: 2392 recTot = self.ReclvlRate['rate'][:, itemp].sum() 2393 else: 2394 recTot = 0. 2395 # 2396 # print ' rec = ', rec 2397 popmat[-1, ci] += self.EDensity[itemp]*self.IonizRate['rate'][itemp] 2398 popmat[ci, ci] -= self.EDensity[itemp]*self.IonizRate['rate'][itemp] 2399 popmat[ci, -1] += self.EDensity[itemp]*(higher.RecombRate['rate'][itemp] - recTot) 2400 popmat[-1, -1] -= self.EDensity[itemp]*(higher.RecombRate['rate'][itemp] - recTot) 2401 # popmat[ci, -1] += self.EDensity[itemp]*higher.RecombRate['rate'][itemp] 2402 # popmat[-1, -1] -= self.EDensity[itemp]*higher.RecombRate['rate'][itemp] 2403 # 2404 for itrans in range(self.Nreclvl): 2405 lvl1 = reclvl['lvl1'][itrans]-1 2406 lvl2 = reclvl['lvl2'][itrans]-1 2407 popmat[lvl2+ci, -1] += self.EDensity[itemp]*self.ReclvlRate['rate'][itrans, itemp] 2408 popmat[-1, -1] -= self.EDensity[itemp]*self.ReclvlRate['rate'][itrans, itemp] 2409 # normalize to unity 2410 # if self.Dielectronic: 2411 # branch = np.zeros(self.Ndielsplups, 'float64') 2412 # for isplups in range(0,self.Ndielsplups): 2413 # l1 = self.DielSplups["lvl1"][isplups]-1 + nlvls 2414 # l2 = self.DielSplups["lvl2"][isplups]-1 2415 # auto = rad[l1+ci, l2+ci] 2416 # avalue = rad[:, l2+ci] 2417 # good = avalue > 0. 2418 # avalueTot = avalue[good].sum() 2419 # branch[isplups] = (avalueTot-auto)/avalueTot 2420 ## print ' l1 %4i l2 %4i auto %10.2e avalue %10.2e tot %10.3f'%( l1, l2, auto, avalueTot, branch[isplups]) 2421 # self.DielUpsilon['branch'] = branch 2422 # # 2423 # dielTot = 0. 2424 # for isplups in range(0,self.Ndielsplups): 2425 # l1 = self.DielSplups["lvl1"][isplups]-1 2426 # l2 = self.DielSplups["lvl2"][isplups]-1 2427 # # 2428 # popmat[l2+ci,l1+ci] += self.EDensity[itemp]*dielexRate[isplups, itemp] 2429 # popmat[l1+ci,l1+ci] -= self.EDensity[itemp]*dielexRate[isplups, itemp] 2430 # 2431 # dielTot += self.EDensity[itemp]*dielexRate[isplups, itemp]*branch[isplups] 2432 # else: 2433 # dielTot = 0. 2434 2435 norm=np.ones(nlvls+ci+rec,'float64') 2436 if ci: 2437 norm[0] = 0. 2438 if rec: 2439 norm[-1] = 0. 2440 if self.Dielectronic: 2441 norm[-1] = 0. 2442 popmat[nlvls+ci+rec-1]=norm 2443 b=np.zeros(nlvls+ci+rec,'float64') 2444 b[nlvls+ci+rec-1]=1. 2445 try: 2446 thispop=np.linalg.solve(popmat,b) 2447 pop[itemp] = thispop[ci:ci+nlvls] 2448 except np.linalg.LinAlgError: 2449 pop[itemp] = np.zeros(nlvls, 'float64') 2450 # print ' error in matrix inversion, setting populations to zero at T = ', ('%8.2e')%(temperature[itemp]) 2451 # thispop=np.linalg.solve(popmat,b) 2452 # if rec: 2453 # pop[itemp] = thispop[ci:ci+nlvls+rec-1] 2454 # else: 2455 # pop[itemp] = thispop[ci:] 2456 # pop[itemp] = thispop[ci:ci+nlvls] 2457 # 2458 pop=np.where(pop >0., pop,0.) 2459 self.Population={"temperature":temperature,"eDensity":eDensity,"population":pop, "protonDensity":protonDensity, "ci":ci, "rec":rec} 2460 # 2461 return

2462 # 2463 # ------------------------------------------------------------------------------------- 2464 #

2465 - def popPlot(self,top=10, plotFile=0, outFile=0, pub=0):

2466 """ 2467 Plots populations vs temperature or eDensity. 2468 2469 top specifies the number of the most highly populated levels to plot 2470 if pub is set, the want publication plots (bw, lw=2). 2471 """ 2472 #self.Population={"temperature":temperature,"eDensity":eDensity,"population":pop} 2473 if pub: 2474 fontsize=16 2475 else: 2476 fontsize=14 2477 # 2478 if hasattr(self, 'Population'): 2479 temperature=self.Population["temperature"] 2480 eDensity=self.Population["eDensity"] 2481 pop=self.Population["population"] 2482 else: 2483 self.populate() 2484 temperature=self.Population["temperature"] 2485 eDensity=self.Population["eDensity"] 2486 pop=self.Population["population"] 2487 # 2488 # for case of only a single density and temperature 2489 if len(pop.shape) == 1: 2490 spop = np.sort(pop) 2491 idx = np.argsort(pop) 2492 minPop = spop[-top:].min()/2. 2493 if top > pop.size: 2494 top = pop.size 2495 for itop in range(1, top+1): 2496 x = [idx[-itop], idx[-itop], idx[-itop]+1, idx[-itop]+1] 2497 y = [minPop, spop[-itop], spop[-itop], minPop] 2498 pl.semilogy(x, y, 'k') 2499 pl.axis([0, max(idx[-top:])+1, minPop, 1.]) 2500 pl.xlabel('Level', fontsize=fontsize) 2501 pl.ylabel('Population', fontsize=fontsize) 2502 return 2503 # 2504 # find the top most populated levels 2505 # 2506 lvl=self.Elvlc["lvl"] 2507 # nlvls=len(lvl) 2508 nlvls = self.Nlvls 2509 if top > nlvls: 2510 top = nlvls 2511 maxpop=np.zeros(nlvls,'Float64') 2512 for ilvl in range(nlvls): 2513 maxpop[ilvl]=pop[:,ilvl].max() 2514 # 2515 lvlsort=np.take(lvl,np.argsort(maxpop)) 2516 toplvl=lvlsort[-top:] 2517 # 2518 # temp=np.asarray(temperature,'Float32') 2519 ntemp = temperature.size 2520 if ntemp > 0: 2521 # if temperature.all() == temperature[0]: 2522 if temperature[0] == temperature[-1]: 2523 ntemp = 1 2524 # 2525 ndens = eDensity.size 2526 if ndens > 0: 2527 # if eDensity.all() == eDensity[0]: 2528 if eDensity[0] == eDensity[-1]: 2529 ndens = 1 2530 # 2531 # 2532 ylabel='Population' 2533 title=self.Spectroscopic 2534 # 2535 pl.figure() 2536 # 2537 pl.ion() 2538 # 2539 # 2540 if ndens == 1: 2541 toppops = np.zeros((top, ntemp), 'float64') 2542 for ilvl in range(top): 2543 toppops[ilvl] = pop[:, toplvl[ilvl]-1] 2544 nonzero = toppops > 0. 2545 ymin = min(toppops[nonzero]) 2546 for lvl in toplvl: 2547 # for some low temperature, populations can not be calculated 2548 good = pop[:, lvl-1] > 0 2549 if pub: 2550 pl.loglog(temperature[good],pop[good,lvl-1], 'k',lw=2) 2551 else: 2552 pl.loglog(temperature[good],pop[good,lvl-1]) 2553 skip=3 2554 if good.sum() == ntemp: 2555 start=divmod(lvl,ntemp)[1] 2556 for itemp in range(start,ntemp,ntemp//skip): 2557 pl.text(temperature[itemp],pop[itemp,lvl-1],str(lvl)) 2558 else: 2559 newtemp=[] 2560 for i, one in enumerate(temperature): 2561 if good[i]: 2562 newtemp.append(one) 2563 start = divmod(lvl, len(newtemp))[1] + ntemp - good.sum() 2564 for itemp in range(start,ntemp,ntemp//skip): 2565 pl.text(temperature[itemp],pop[itemp,lvl-1],str(lvl)) 2566 xlabel='Temperature (K)' 2567 pl.xlabel(xlabel,fontsize=fontsize) 2568 pl.ylabel(ylabel,fontsize=fontsize) 2569 dstr=' - Density = %10.2e (cm$^{-3}$)' % eDensity[0] 2570 pl.title(title+dstr,fontsize=fontsize) 2571 pl.xlim(temperature.min(),temperature.max()) 2572 # nonzero = pop 2573 # yl=pl.ylim() 2574 pl.ylim(ymin,1.2) 2575 elif ntemp == 1: 2576 xlabel=r'Electron Density (cm$^{-3}$)' 2577 # for lvl in toplvl: 2578 # pl.loglog(eDensity,pop[:,lvl-1]) 2579 # skip=min(3, ndens) 2580 # start=divmod(lvl,ndens)[1] 2581 # for idens in range(start,ndens,ndens//skip): 2582 # pl.text(eDensity[idens],pop[idens,lvl-1],str(lvl)) 2583 toppops = np.zeros((top, ndens), 'float64') 2584 for ilvl in range(top): 2585 toppops[ilvl] = pop[:, toplvl[ilvl]-1] 2586 nonzero = toppops > 0. 2587 ymin = min(toppops[nonzero]) 2588 for lvl in toplvl: 2589 # for some low temperature, populations can not be calculated 2590 good = pop[:, lvl-1] > 0 2591 if pub: 2592 pl.loglog(eDensity[good],pop[good,lvl-1], 'k', lw=2) 2593 else: 2594 pl.loglog(eDensity[good],pop[good,lvl-1]) 2595 skip=3 2596 if good.sum() == ndens: 2597 start=divmod(lvl,ndens)[1] 2598 for idens in range(start,ndens,ndens//skip): 2599 pl.text(eDensity[idens],pop[idens,lvl-1],str(lvl)) 2600 else: 2601 newdens=[] 2602 for i, one in enumerate(eDensity): 2603 if good[i]: 2604 newdens.append(one) 2605 start = divmod(lvl, len(newdens))[1] + ndens - good.sum() 2606 for idens in range(start,ndens,ndens//skip): 2607 pl.text(eDensity[idens],pop[idens, lvl-1],str(lvl)) 2608 pl.xlabel(xlabel,fontsize=fontsize) 2609 pl.ylabel(ylabel,fontsize=fontsize) 2610 tstr=' - T = %10.2e (K)' % temperature[0] 2611 pl.title(title+tstr,fontsize=fontsize) 2612 pl.xlim(eDensity[eDensity.nonzero()].min(),eDensity.max()) 2613 yl=pl.ylim() 2614 pl.ylim(yl[0],1.2) 2615 else: 2616 # pl.figure() 2617 ax = pl.subplot(111) 2618 # for lvl in toplvl: 2619 # pl.loglog(temperature,pop[:,lvl-1]) 2620 # skip = min(3, ntemp) 2621 # start=divmod(lvl,ntemp)[1] 2622 # for itemp in range(start,ntemp,ntemp//skip): 2623 # pl.text(temperature[itemp],pop[itemp,lvl-1],str(lvl)) 2624 toppops = np.zeros((top, ntemp), 'float64') 2625 for ilvl in range(top): 2626 toppops[ilvl] = pop[:, toplvl[ilvl]-1] 2627 nonzero = toppops > 0. 2628 ymin = min(toppops[nonzero]) 2629 for lvl in toplvl: 2630 # for some low temperature, populations can not be calculated 2631 good = pop[:, lvl-1] > 0 2632 if pub: 2633 pl.loglog(temperature[good],pop[good,lvl-1], 'k', lw=2) 2634 else: 2635 pl.loglog(temperature[good],pop[good,lvl-1]) 2636 skip=3 2637 if good.sum() == ntemp: 2638 start=divmod(lvl,ntemp)[1] 2639 for itemp in range(start,ntemp,ntemp//skip): 2640 pl.text(temperature[itemp],pop[itemp,lvl-1],str(lvl)) 2641 else: 2642 newtemp=[] 2643 for i, one in enumerate(temperature): 2644 if good[i]: 2645 newtemp.append(one) 2646 start = divmod(lvl, len(newtemp))[1] + ntemp - good.sum() 2647 for itemp in range(start,ntemp,ntemp//skip): 2648 pl.text(temperature[itemp],pop[itemp,lvl-1],str(lvl)) 2649 xlabel='Temperature (K)' 2650 pl.xlabel(xlabel,fontsize=fontsize) 2651 pl.ylabel(ylabel,fontsize=fontsize) 2652 # pl.title(title,fontsize=fontsize) 2653 # pl.xlim(temperature.min(),temperature.max()) 2654 # yl=pl.ylim() 2655 # pl.ylim(ymin,1.2) 2656 pl.axis([temperature.min(),temperature.max(), ymin, 1.2]) 2657 pl.text(0.1, 0.5,title, horizontalalignment='center', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 2658 # 2659 ax2 = pl.twiny() 2660 xlabel=r'Electron Density (cm$^{-3}$)' 2661 pl.xlabel(xlabel, fontsize=fontsize) 2662 pl.loglog(eDensity,pop[:,toplvl[0]], visible=False) 2663 ax2.xaxis.tick_top() 2664 # pl.figure() 2665 # for lvl in toplvl: 2666 # pl.loglog(eDensity,pop[:,lvl-1]) 2667 # skip = min(3, ntemp) 2668 # start=divmod(lvl,ndens)[1] 2669 # for idens in range(start,ndens,ndens//skip): 2670 # pl.text(eDensity[idens],pop[idens,lvl-1],str(lvl)) 2671 # xlabel=r'Electron Density (cm$^{-3}$)' 2672 # pl.xlabel(xlabel,fontsize=fontsize) 2673 # pl.ylabel(ylabel,fontsize=fontsize) 2674 # pl.title(title,fontsize=fontsize) 2675 # pl.xlim(eDensity.min(),eDensity.max()) 2676 # yl=pl.ylim() 2677 # pl.ylim(yl[0],1.2) 2678 if outFile: 2679 pl.savefig(outFile) 2680 self.Population['toplvl'] = toplvl 2681 return

2682 # 2683 # ------------------------------------------------------------------------------------- 2684 #

2685 - def emiss(self, wvlRange = 0, allLines=1):

2686 """ 2687 Calculate the emissivities for lines of the specified ion. 2688 2689 wvlRange can be set to limit the calculation to a particular wavelength range 2690 2691 units: ergs cm^-3 s^-1 str^-1 2692 2693 Does not include elemental abundance or ionization fraction 2694 2695 Wavelengths are sorted 2696 set allLines = 1 to include unidentified lines 2697 """ 2698 # 2699 # 2700 if hasattr(self, 'Population'): 2701 # doPopulate=False 2702 pop=self.Population['population'] 2703 else: 2704 self.populate() 2705 pop=self.Population["population"] 2706 # 2707 # 2708 # nlvls=len(self.Elvlc['lvl']) 2709 ## good=self.Wgfa['avalue'] > 0. 2710 # using [:] to make a copy things don't change elsewhere 2711 wvl = np.asarray(self.Wgfa["wvl"], 'float64') 2712 obs = np.where(wvl > 0., 'Y', 'N') 2713 if allLines: 2714 wvl=np.abs(wvl) 2715 l1 = np.asarray(self.Wgfa['lvl1'], 'int64') 2716 l2 = np.asarray(self.Wgfa["lvl2"], 'int64') 2717 avalue = np.asarray(self.Wgfa["avalue"], 'float64') 2718 if 'pretty1' in self.Wgfa.keys(): 2719 pretty1 = np.asarray(self.Wgfa['pretty1']) 2720 if 'pretty2' in self.Wgfa.keys(): 2721 pretty2 = np.asarray(self.Wgfa['pretty2']) 2722 # 2723 # make sure there are lines in the wavelength range, if specified 2724 2725 if wvlRange: 2726 realgood = util.between(wvl, wvlRange) 2727 l1 = l1[realgood] 2728 l2 = l2[realgood] 2729 wvl = wvl[realgood] 2730 avalue = avalue[realgood] 2731 obs = obs[realgood] 2732 if 'pretty1' in self.Wgfa.keys(): 2733 pretty1 = pretty1[realgood] 2734 if 'pretty2' in self.Wgfa.keys(): 2735 pretty2 = pretty2[realgood] 2736 # 2737 # two-photon decays have wvl=0 and nonzero avalues 2738 # zed = wvl.count(0.) 2739 nonzed = wvl != 0. 2740 wvl = wvl[nonzed] 2741 l1 = l1[nonzed] 2742 l2 = l2[nonzed] 2743 avalue = avalue[nonzed] 2744 pretty1 = pretty1[nonzed] 2745 pretty2 = pretty2[nonzed] 2746 obs = obs[nonzed] 2747 nwvl=len(wvl) 2748 # 2749 if nwvl == 0: 2750 self.Emiss = {'errorMessage':self.Spectroscopic+' no lines in this wavelength range'} 2751 return 2752 # 2753 # 2754 try: 2755 ntempden,nlvls = pop.shape 2756 em=np.zeros((nwvl, ntempden),'Float64') 2757 # if self.Temperature.size < ntempden: 2758 # temperature = np.repeat(self.Temperature, ntempden) 2759 # else: 2760 # temperature = self.Temperature 2761 if self.EDensity.size < ntempden: 2762 eDensity = np.repeat(self.EDensity, ntempden) 2763 else: 2764 eDensity = self.EDensity 2765 except: 2766 # nlvls=len(pop) 2767 ntempden=1 2768 em=np.zeros(nwvl,'Float64') 2769 eDensity = self.EDensity 2770 # temperature = self.Temperature 2771 # 2772 plotLabels={} 2773 # 2774 if self.Defaults['wavelength'] == 'angstrom': 2775 plotLabels["xLabel"]="Angstroms" 2776 elif self.Defaults['wavelength'] == 'nm': 2777 plotLabels["xLabel"]="nanometers" 2778 elif self.Defaults['wavelength'] == 'kev': 2779 plotLabels["xLabel"] = "kev" 2780 # 2781 if self.Defaults['flux'] == 'energy': 2782 factor=const.planck*const.light/(4.*const.pi*1.e-8*wvl) 2783 plotLabels["yLabel"]="ergs cm^-3 s^-1" 2784 elif self.Defaults['flux'] == 'photon': 2785 factor=np.ones((nwvl),'Float64')/(4.*const.pi) 2786 plotLabels["yLabel"]="photons cm^-3 s^-1" 2787 # 2788 if ntempden > 1: 2789 for itempden in range(ntempden): 2790 for iwvl in range(nwvl): 2791 p = pop[itempden,l2[iwvl]-1] 2792 # p = pop[l2[iwvl]-1, itempden] 2793 em[iwvl, itempden] = factor[iwvl]*p*avalue[iwvl] 2794 if self.Defaults['wavelength'] == 'kev': 2795 wvl = const.ev2Ang/np.asarray(wvl) 2796 elif self.Defaults['wavelength'] == 'nm': 2797 wvl = wvl/10. 2798 # em = em.take(wvl.argsort(),axis=0) 2799 # wvl.sort() 2800 # idx = np.argsort(wvl) 2801 # l1 = l1[idx] 2802 # l2 = l2[idx] 2803 else: 2804 for iwvl in range(0,nwvl): 2805 p=pop[l2[iwvl]-1] 2806 em[iwvl]=factor[iwvl]*p*avalue[iwvl] 2807 if self.Defaults['wavelength'] == 'kev': 2808 wvlE=const.ev2Ang/np.asarray(wvl) 2809 elif self.Defaults['wavelength'] == 'nm': 2810 wvl=wvl/10. 2811 # idx = np.argsort(wvl) 2812 # wvl = wvl[idx] 2813 # em = em[idx] 2814 # l1 = l1[idx] 2815 # l2 = l2[idx] 2816 # lvl1 = l1.tolist() 2817 # lvl2 = l2.tolist() 2818 nlvl = len(l1) 2819 ionS = np.asarray([self.IonStr]*nlvl) 2820 Emiss = {'ionS':ionS,"wvl":wvl, "emiss":em, "plotLabels":plotLabels, 'lvl1':l1, 'lvl2':l2, 'avalue':avalue, 'obs':obs, 'pretty1':pretty1, 'pretty2':pretty2} 2821 # if self.Wgfa.has_key('pretty1'): 2822 # Emiss['pretty1'] = np.asarray(p1) 2823 # if self.Wgfa.has_key('pretty2'): 2824 # Emiss['pretty2'] = np.asarray(p2) 2825 # everything in Emiss should be a numpy array 2826 self.Emiss = Emiss 2827 return

2828 # 2829 # --------------------------------------------------------------------------- 2830 #

2831 - def emissList(self, index=-1, wvlRange=None, wvlRanges=None, top=10, relative=0, outFile=0 ):

2832 ''' 2833 List the emissivities. 2834 2835 wvlRange, a 2 element tuple, list or array determines the wavelength range 2836 2837 Top specifies to plot only the top strongest lines, default = 10 2838 2839 normalize = 1 specifies whether to normalize to strongest line, default = 0 2840 ''' 2841 # 2842 if outFile: 2843 output = open(outFile, 'w') 2844 # 2845 if not hasattr(self, 'Emiss'): 2846 try: 2847 self.emiss() 2848 except: 2849 print(' emissivities not calculated or emiss() is unable to calculate them') 2850 print(' perhaps the temperature and/or eDensity are not set') 2851 return 2852 # 2853 emissivity = copy.copy(self.Emiss) 2854 emiss = emissivity['emiss'] 2855 ionS = emissivity['ionS'] 2856 wvl = emissivity['wvl'] 2857 lvl1 = emissivity['lvl1'] 2858 lvl2 = emissivity['lvl2'] 2859 avalue = emissivity['avalue'] 2860 obs = emissivity['obs'] 2861 pretty1 = emissivity['pretty1'] 2862 pretty2 = emissivity['pretty2'] 2863 # plotLabels = emissivity['plotLabels'] 2864 # 2865 temperature = self.Temperature 2866 eDensity = self.EDensity 2867 # 2868 ndens = eDensity.size 2869 ntemp = temperature.size 2870 # 2871 if ndens == 1 and ntemp == 1: 2872 dstr = ' - Density = %10.2e (cm$^{-3}$)' %(eDensity) 2873 tstr = ' - T = %10.2e (K)' %(temperature) 2874 elif ndens == 1 and ntemp > 1: 2875 if index < 0: 2876 index = ntemp/2 2877 print('using index = %5i specifying temperature = %10.2e'%(index, temperature[index])) 2878 self.Message = 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 2879 2880 emiss=emiss[:, index] 2881 elif ndens > 1 and ntemp == 1: 2882 if index < 0: 2883 index = ntemp/2 2884 print('using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index])) 2885 self.Message = 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 2886 emiss=emiss[:, index] 2887 elif ndens > 1 and ntemp > 1: 2888 if index < 0: 2889 index = ntemp/2 2890 print('using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index])) 2891 self.Message = 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index]) 2892 emiss=emiss[:, index] 2893 # 2894 if wvlRange: 2895 wvlIndex=util.between(wvl,wvlRange) 2896 elif wvlRanges: 2897 wvlIndex = [] 2898 for awvlRange in wvlRanges: 2899 wvlIndex.extend(util.between(wvl,awvlRange)) 2900 else: 2901 wvlIndex = range(wvl.size) 2902 # 2903 # 2904 # get those in the right wavelength range 2905 # 2906 emiss = emiss[wvlIndex] 2907 ionS = ionS[wvlIndex] 2908 wvl = wvl[wvlIndex] 2909 lvl1 = lvl1[wvlIndex] 2910 lvl2 = lvl2[wvlIndex] 2911 avalue = avalue[wvlIndex] 2912 obs = obs[wvlIndex] 2913 pretty1 = pretty1[wvlIndex] 2914 pretty2 = pretty2[wvlIndex] 2915 # 2916 self.Error = 0 2917 if wvl.size == 0: 2918 print('No lines in this wavelength interval') 2919 self.Error = 1 2920 self.Message = 'No lines in this wavelength interval' 2921 return 2922 # 2923 elif top == 0: 2924 top = wvl.size 2925 elif top > wvl.size: 2926 top = wvl.size 2927 # 2928 # now sort by intensity/emissivity 2929 isrt = np.argsort(emiss) 2930 ionS = ionS[isrt[-top:]] 2931 wvl = wvl[isrt[-top:]] 2932 lvl1 = lvl1[isrt[-top:]] 2933 lvl2 = lvl2[isrt[-top:]] 2934 obs = obs[isrt[-top:]] 2935 emiss = emiss[isrt[-top:]] 2936 avalue = avalue[isrt[-top:]] 2937 pretty1 = pretty1[isrt[-top:]] 2938 pretty2 = pretty2[isrt[-top:]] 2939 # 2940 # must follow setting top 2941 # 2942 if relative: 2943 emiss = emiss/emiss[:top].max() 2944 # 2945 # 2946 idx = np.argsort(wvl) 2947 # 2948 fmt = '%5s %5i %5i %25s - %25s %12.3f %12.3e %12.2e %1s' 2949 print( ' ') 2950 print( '------------------------------------------') 2951 print(' ') 2952 print(' Ion lvl1 lvl2 lower upper Wvl(A) Emissivity A value Obs') 2953 for kdx in idx: 2954 print(fmt%(ionS[kdx], lvl1[kdx], lvl2[kdx], pretty1[kdx], pretty2[kdx], wvl[kdx], emiss[kdx], avalue[kdx], obs[kdx])) 2955 print(' ') 2956 print(' ------------------------------------------') 2957 print(' ') 2958 # 2959 self.Emiss['wvlTop'] = wvl[idx] 2960 self.Emiss['emissTop'] = emiss[idx] 2961 if outFile: 2962 fmt = '%5s %5i %5i %25s - %25s %12.3f %12.3e %12.2e %1s' 2963 output.write(' \n') 2964 output.write(' ------------------------------------------ \n') 2965 output.write(' ') 2966 output.write(' Ion lvl1 lvl2 lower upper Wvl(A) Emissivity A value Obs \n') 2967 for kdx in idx: 2968 output.write(fmt%(ionS[kdx], lvl1[kdx], lvl2[kdx], pretty1[kdx], pretty2[kdx], wvl[kdx], emiss[kdx], avalue[kdx], obs[kdx])) 2969 output.write(' \n') 2970 output.write(' ------------------------------------------ \n') 2971 output.write(' \n ') 2972 output.close()

2973 # 2974 # --------------------------------------------------------------------------- 2975 #

2976 - def emissPlot(self, index=-1, wvlRange=None, top=10, linLog='lin', relative=0, verbose=0, plotFile = 0, saveFile=0 ):

2977 '''Plot the emissivities. 2978 2979 wvlRange, a 2 element tuple, list or array determines the wavelength range 2980 2981 Top specifies to plot only the top strongest lines, default = 10 2982 2983 linLog specifies a linear or log plot, want either lin or log, default = lin 2984 2985 normalize = 1 specifies whether to normalize to strongest line, default = 0''' 2986 # 2987 title=self.Spectroscopic 2988 # 2989 # doEmiss=False 2990 if hasattr(self, 'Emiss'): 2991 em = self.Emiss 2992 else: 2993 try: 2994 self.emiss() 2995 em = self.Emiss 2996 except: 2997 print(' emissivities not calculated and emiss() is unable to calculate them') 2998 print(' perhaps the temperature and/or eDensity are not set') 2999 return 3000 emiss = em['emiss'] 3001 wvl = em['wvl'] 3002 temperature = self.Temperature 3003 eDensity = self.EDensity 3004 # 3005 ndens = eDensity.size 3006 ntemp = temperature.size 3007 # 3008 if ndens == 1 and ntemp == 1: 3009 dstr = ' - Density = %10.2e (cm$^{-3}$)' %(eDensity) 3010 tstr = ' - T = %10.2e (K)' %(temperature) 3011 elif ndens == 1 and ntemp > 1: 3012 if index < 0: 3013 index = ntemp/2 3014 print('using index = %5i specifying temperature = %10.2e'%(index, temperature[index])) 3015 self.Message = 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 3016 # if chInteractive: 3017 # print 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 3018 # else: 3019 # self.Message = 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 3020 emiss=emiss[:, index] 3021 dstr=' - Density = %10.2e (cm$^{-3}$)' % eDensity 3022 tstr=' - T = %10.2e (K)' % temperature[index] 3023 elif ndens > 1 and ntemp == 1: 3024 if index < 0: 3025 index = ntemp/2 3026 print('using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index])) 3027 self.Message = 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 3028 # if chInteractive: 3029 # print 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 3030 # else: 3031 # self.Message = 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 3032 emiss=emiss[:, index] 3033 dstr=' - Density = %10.2e (cm$^{-3}$)' % eDensity[index] 3034 tstr=' - T = %10.2e (K)' % temperature 3035 elif ndens > 1 and ntemp > 1: 3036 if index < 0: 3037 index = ntemp/2 3038 print('using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index])) 3039 self.Message = 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index]) 3040 # if chInteractive: 3041 # print 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index]) 3042 # else: 3043 # self.Message = 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index]) 3044 emiss=emiss[:, index] 3045 dstr=' - Density = %10.2e (cm$^{-3}$)' % eDensity[index] 3046 tstr=' - T = %10.2e (K)' % temperature[index] 3047 if type(wvlRange) != type(None): 3048 wvlIndex = util.between(wvl, wvlRange) 3049 else: 3050 wvlIndex = range(wvl.size) 3051 emiss = emiss[wvlIndex] 3052 wvl = wvl[wvlIndex] 3053 # 3054 self.Error = 0 3055 if wvl.size == 0: 3056 print('No lines in this wavelength interval') 3057 self.Error = 1 3058 self.Message = 'No lines in this wavelength interval' 3059 # if chInteractive: 3060 # print 'No lines in this wavelength interval' 3061 # else: 3062 # self.Error = 1 3063 # self.Message = 'No lines in this wavelength interval' 3064 return 3065 elif top == 0: 3066 top = wvl.size 3067 elif wvl.size > top: 3068 isrt = np.argsort(emiss) 3069 wvl = wvl[isrt[-top:]] 3070 emiss = emiss[isrt[-top:]] 3071 else: 3072 top = wvl.size 3073 # must follow setting top 3074 # 3075 pl.figure() 3076 ylabel = 'Emissivity' 3077 if relative: 3078 emiss = emiss/emiss[:top].max() 3079 ylabel += ' (Relative)' 3080 # 3081 xlabel = 'Wavelength ('+self.Defaults['wavelength'] +')' 3082 # 3083 ymin = 10.**(np.log10(emiss.min()).round(0)-0.5 ) 3084 # 3085 pl.ion() 3086 # if chInteractive: 3087 # pl.ion() 3088 # else: 3089 # pl.ioff() 3090 # 3091 for idx in range(top): 3092 xx=[wvl[idx], wvl[idx]] 3093 if linLog == 'lin': 3094 yy=[0., emiss[idx]] 3095 pl.plot(xx, yy) 3096 else: 3097 yy=[ymin/10., emiss[idx]] 3098 pl.semilogy(xx, yy) 3099 pl.xlabel(xlabel) 3100 pl.ylabel(ylabel) 3101 pl.title(title+tstr+dstr) 3102 if wvlRange: 3103 pl.axis([wvlRange[0], wvlRange[1], ymin, emiss.max()]) 3104 if plotFile: 3105 pl.savefig(plotFile) 3106 # 3107 idx = np.argsort(wvl) 3108 self.Emiss['wvlTop'] = wvl[idx] 3109 self.Emiss['emissTop'] = emiss[idx]

3110 # 3111 # ------------------------------------------------------------------------------------- 3112 #

3113 - def emissRatio(self,wvlRange=None, wvlRanges=None,top=10):

3114 """ 3115 Plot the ratio of 2 lines or sums of lines. 3116 Shown as a function of density and/or temperature. 3117 For a single wavelength range, set wvlRange = [wMin, wMax] 3118 For multiple wavelength ranges, set wvlRanges = [[wMin1,wMax1],[wMin2,wMax2], ...] 3119 A plot of relative emissivities is shown and then a dialog appears for the user to 3120 choose a set of lines. 3121 """ 3122 # 3123 # self.Emiss={"temperature":temperature,"density":density,"wvl":wvl,"emiss":em, 3124 # "plotLabels":plotLabels} 3125 # 3126 if hasattr(self, 'Emiss'): 3127 doEmiss = False 3128 em = self.Emiss 3129 else: 3130 doEmiss = True 3131 # 3132 # 3133 if doEmiss: 3134 # new values of temperature or eDensity 3135 self.emiss() 3136 em = self.Emiss 3137 # 3138 # 3139 fontsize=14 3140 # 3141 # temperature = self.Temperature 3142 eDensity = self.EDensity 3143 emiss = em['emiss'] 3144 ionS = em['ionS'] 3145 wvl = em["wvl"] 3146 lineLabel = [] 3147 for iline, ions in enumerate(ionS): 3148 lineLabel.append(ions+' '+str(wvl[iline])) 3149 # 3150 plotLabels = em["plotLabels"] 3151 xLabel = plotLabels["xLabel"] 3152 yLabel = plotLabels["yLabel"] 3153 # 3154 # find which lines are in the wavelength range if it is set 3155 # 3156 # 3157 if wvlRange: 3158 igvl=util.between(wvl,wvlRange) 3159 elif wvlRanges: 3160 igvl = [] 3161 for awvlRange in wvlRanges: 3162 igvl.extend(util.between(wvl,awvlRange)) 3163 else: 3164 igvl=range(len(wvl)) 3165 # 3166 nlines=len(igvl) 3167 # 3168 # print ' nlines = ',nlines 3169 # print ' iglv = ',igvl 3170 igvl=np.take(igvl,wvl[igvl].argsort()) 3171 # find the top most intense lines 3172 # 3173 if top > nlines: 3174 top=nlines 3175 # 3176 maxEmiss = np.zeros(nlines,'Float64') 3177 print(' maxEmiss.shape = %s'%(str(maxEmiss.shape))) 3178 for iline in range(nlines): 3179 maxEmiss[iline] = emiss[igvl[iline]].max() 3180 for iline in range(nlines): 3181 if maxEmiss[iline] == maxEmiss.max(): 3182 maxAll = emiss[igvl[iline]] 3183 # line=range(nlines) 3184 igvlsort=np.take(igvl,np.argsort(maxEmiss)) 3185 # print 'igvlsort = ', igvlsort 3186 topLines=igvlsort[-top:] 3187 # print ' topLines = ', topLines 3188 maxWvl='%5.3f' % wvl[topLines[-1]] 3189 # maxline=topLines[-1] 3190 # 3191 topLines=topLines[wvl[topLines].argsort()] 3192 # 3193 # 3194 # need to make sure there are no negative values before plotting 3195 good = np.where(emiss > 0.) 3196 emissMin=emiss[good].min() 3197 bad=np.where(emiss <= 0.) 3198 emiss[bad]=emissMin 3199 # 3200 # 3201 ntemp=self.Temperature.size 3202 # 3203 ndens=self.EDensity.size 3204 # 3205 ylabel='Emissivity relative to '+maxWvl 3206 title=self.Spectroscopic 3207 # 3208 # 3209 if ndens==1 and ntemp==1: 3210 print(' only a single temperature and eDensity') 3211 return 3212 elif ndens == 1: 3213 xlabel='Temperature (K)' 3214 xvalues=self.Temperature 3215 outTemperature=self.Temperature 3216 outDensity=np.zeros(ntemp,'Float64') 3217 outDensity.fill(self.EDensity) 3218 desc_str=' at Density = %10.2e (cm)$^{-3}$' % self.EDensity 3219 elif ntemp == 1: 3220 xvalues=self.EDensity 3221 outTemperature=np.zeros(ndens,'Float64') 3222 outTemperature.fill(self.Temperature) 3223 outDensity=self.EDensity 3224 xlabel=r'$\rm{Electron Density (cm)^{-3}}$' 3225 desc_str=' at Temp = %10.2e (K)' % self.Temperature 3226 else: 3227 outTemperature=self.Temperature 3228 outDensity=self.EDensity 3229 xlabel='Temperature (K)' 3230 xvalues=self.Temperature 3231 desc_str=' for variable Density' 3232 # 3233 # put all actual plotting here 3234 # 3235 pl.ion() 3236 # if chInteractive: 3237 # pl.ion() 3238 # else: 3239 # pl.ioff() 3240 # 3241 # maxAll is an array 3242 ymax = np.max(emiss[topLines[0]]/maxAll) 3243 ymin = ymax 3244 pl.figure() 3245 ax = pl.subplot(111) 3246 nxvalues=len(xvalues) 3247 for iline in range(top): 3248 tline=topLines[iline] 3249 pl.loglog(xvalues,emiss[tline]/maxAll) 3250 if np.min(emiss[tline]/maxAll) < ymin: 3251 ymin = np.min(emiss[tline]/maxAll) 3252 if np.max(emiss[tline]/maxAll) > ymax: 3253 ymax = np.max(emiss[tline]/maxAll) 3254 skip=2 3255 start=divmod(iline,nxvalues)[1] 3256 for ixvalue in range(start,nxvalues,nxvalues//skip): 3257 pl.text(xvalues[ixvalue],emiss[tline,ixvalue]/maxAll[ixvalue],str(wvl[tline])) 3258 pl.xlim(xvalues.min(),xvalues.max()) 3259 # pl.ylim(ymin, ymax) 3260 pl.xlabel(xlabel,fontsize=fontsize) 3261 pl.ylabel(ylabel,fontsize=fontsize) 3262 if ndens == ntemp and ntemp > 1: 3263 pl.text(0.07, 0.5,title, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 3264 # 3265 ax2 = pl.twiny() 3266 xlabelDen=r'Electron Density (cm$^{-3}$)' 3267 pl.xlabel(xlabelDen, fontsize=fontsize) 3268 pl.loglog(eDensity,emiss[topLines[top-1]]/maxAll, visible=False) 3269 ax2.xaxis.tick_top() 3270 pl.ylim(ymin/1.2, 1.2*ymax) 3271 else: 3272 pl.ylim(ymin/1.2, 1.2*ymax) 3273 pl.title(title+desc_str,fontsize=fontsize) 3274 pl.draw() 3275 # need time to let matplotlib finish plotting 3276 time.sleep(0.5) 3277 # 3278 # get line selection ************************************************************ 3279 # 3280 selectTags = [] 3281 for itop in topLines: 3282 selectTags.append(ionS[itop]+ ' '+ str(wvl[itop])) 3283 # 3284 numden = chGui.gui.choice2Dialog(selectTags) 3285 # numden = gui.choice2Dialog(wvl[topLines]) 3286 # 3287 # num_idx and den_idx are tuples 3288 # 3289 num_idx=numden.numIndex 3290 if len(num_idx) == 0: 3291 print(' no numerator lines were selected') 3292 return 3293 # 3294 den_idx=numden.denIndex 3295 if len(den_idx) == 0: 3296 print(' no denominator lines were selected') 3297 return 3298 # 3299 numEmiss=np.zeros(len(xvalues),'Float64') 3300 for aline in num_idx: 3301 numEmiss+=emiss[topLines[aline]] 3302 # 3303 denEmiss=np.zeros(len(xvalues),'Float64') 3304 for aline in den_idx: 3305 denEmiss+=emiss[topLines[aline]] 3306 # 3307 # plot the desired ratio 3308 # maxAll is an array 3309 pl.figure() 3310 ax = pl.subplot(111) 3311 pl.loglog(xvalues,numEmiss/denEmiss) 3312 pl.xlim(xvalues.min(),xvalues.max()) 3313 pl.xlabel(xlabel,fontsize=fontsize) 3314 pl.ylabel('Emissivity Ratio ('+self.Defaults['flux']+')',fontsize=fontsize) 3315 desc = '' 3316 for aline in num_idx: 3317 desc += ' ' + selectTags[aline] 3318 # desc += ' ' + str(wvl[topLines[aline]]) 3319 desc +=' / ' 3320 for aline in den_idx: 3321 desc += ' ' + selectTags[aline] 3322 # desc += ' ' + str(wvl[topLines[aline]]) 3323 if ndens == ntemp and ntemp > 1: 3324 pl.text(0.07, 0.5,desc, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 3325 # 3326 ax2 = pl.twiny() 3327 xlabelDen=r'Electron Density (cm$^{-3}$)' 3328 pl.xlabel(xlabelDen, fontsize=fontsize) 3329 pl.loglog(eDensity,numEmiss/denEmiss, visible=False) 3330 ax2.xaxis.tick_top() 3331 else: 3332 # pl.ylim(ymin, ymax) 3333 pl.title(desc,fontsize=fontsize) 3334 # desc=title+' '+str(wvl[num_line])+' / '+str(wvl[den_line])+' '+desc_str 3335 # pl.title(desc, fontsize=fontsize) 3336 # pl.title(title+' '+str(wvl[num_line])+' / '+str(wvl[den_line])+' '+desc_str,fontsize=fontsize) 3337 # pl.draw() 3338 # pl.ioff() 3339 # pl.show() 3340 # 3341 intensityRatioFileName=self.IonStr 3342 for aline in num_idx: 3343 intensityRatioFileName+= '_%3i'%(wvl[topLines[aline]]) 3344 intensityRatioFileName+='_2' 3345 for aline in den_idx: 3346 intensityRatioFileName+= '_%3i'%(wvl[topLines[aline]]) 3347 intensityRatioFileName+='.rat' 3348 self.IntensityRatio={'ratio':numEmiss/denEmiss,'desc':desc, 3349 'temperature':outTemperature,'eDensity':outDensity,'filename':intensityRatioFileName, 'numIdx':num_idx, 'denIdx':den_idx}

3350 # 3351 # --------------------------------------------------------------------------- 3352 #

3353 - def intensity(self, wvlRange = None, allLines=1, em=0):

3354 """ 3355 Calculate the intensities for lines of the specified ion. 3356 3357 wvlRange, a 2 element tuple, list or array determines the wavelength range 3358 3359 units: ergs cm^-3 s^-1 str^-1 3360 3361 includes elemental abundance and ionization fraction. 3362 3363 the emission measure 'em' is included if specified 3364 """ 3365 # emiss ={"wvl":wvl, "emiss":em, "plotLabels":plotLabels} 3366 # 3367 if type(em) == int and em == 0: 3368 if hasattr(self, 'Em'): 3369 em = self.Em 3370 else: 3371 em = np.ones(self.NTempDen, 'float64') 3372 self.Em = em 3373 elif type(em) == float and em > 0.: 3374 em = np.ones(self.NTempDen, 'float64')*em 3375 self.Em = em 3376 elif type(em) == list or type(em) == tuple or type(em) == np.ndarray: 3377 em = np.asarray(em, 'float64') 3378 self.Em = em 3379 # so we know that it has been applied 3380 # 3381 if not hasattr(self, 'Emiss'): 3382 self.emiss(wvlRange = wvlRange, allLines=allLines) 3383 emiss = copy.copy(self.Emiss) 3384 else: 3385 emiss = copy.copy(self.Emiss) 3386 if 'errorMessage' in emiss.keys(): 3387 self.Intensity = {'errorMessage': self.Spectroscopic+' no lines in this wavelength region'} 3388 return 3389 # 3390 # everything in emiss should be a numpy array 3391 emissivity = emiss['emiss'] 3392 ionS = emiss['ionS'] 3393 wvl = emiss['wvl'] 3394 lvl1 = emiss['lvl1'] 3395 lvl2 = emiss['lvl2'] 3396 obs = emiss['obs'] 3397 pretty1 = emiss['pretty1'] 3398 pretty2 = emiss['pretty2'] 3399 avalue = emiss['avalue'] 3400 # 3401 # print(' shape of emissivity %5i'%(len(emissivity.shape))) 3402 # if len(emissivity.shape) > 1: 3403 # print( ' emissivity.shape %5i %5i '%(emissivity.shape)) 3404 # else: 3405 # print(' len of emissivity'%(len(emissivity))) 3406 # 3407 if hasattr(self, 'Abundance'): 3408 ab=self.Abundance 3409 else: 3410 self.Abundance = io.abundanceRead() 3411 ab=self.Abundance 3412 if hasattr(self, 'IoneqOne'): 3413 thisIoneq=self.IoneqOne 3414 else: 3415 self.ioneqOne() 3416 thisIoneq=self.IoneqOne 3417 # try: 3418 # nwvl, ntempden = emissivity.shape 3419 # intensity = np.zeros((ntempden, nwvl),'Float64') 3420 # if thisIoneq.size == 1: 3421 # thisIoneq = np.ones(ntempden, 'float64')*thisIoneq 3422 # for it in range(ntempden): 3423 # intensity[it] = ab*thisIoneq[it]*emissivity[:, it] 3424 # except: 3425 # nwvl=len(emissivity) 3426 # ntempden=1 3427 # intensity = ab*thisIoneq*emissivity 3428 3429 if len(emissivity.shape) > 1: 3430 nwvl, ntempden= emissivity.shape 3431 intensity = np.zeros((ntempden, nwvl),'Float64') 3432 if thisIoneq.size == 1: 3433 thisIoneq = np.ones(ntempden, 'float64')*thisIoneq 3434 for it in range(ntempden): 3435 intensity[it] = ab*thisIoneq[it]*emissivity[:, it]*em[it]/self.EDensity[it] 3436 else: 3437 nwvl=len(emissivity) 3438 ntempden=1 3439 intensity = ab*thisIoneq*emissivity*em/self.EDensity 3440 if ntempden == 1: 3441 integrated = intensity 3442 else: 3443 integrated = intensity.sum(axis=0) 3444 Intensity = {'intensity':intensity, 'integrated':integrated,'ionS':ionS, 'wvl':wvl, 'lvl1':lvl1, 'lvl2':lvl2, 'pretty1':pretty1, 'pretty2':pretty2, 'obs':obs, 'avalue':avalue, 'em':em} 3445 self.Intensity = Intensity

3446 # 3447 # --------------------------------------------------------------------------- 3448 #

3449 - def boundBoundLoss(self, wvlRange = None, allLines=1):

3450 """ 3451 Calculate the summed radiative loss rate for all lines of the specified ion. 3452 3453 wvlRange, a 2 element tuple, list or array determines the wavelength range 3454 3455 units: ergs cm^-3 s^-1 3456 3457 includes elemental abundance and ionization fraction. 3458 """ 3459 # emiss ={"wvl":wvl, "emiss":em, "plotLabels":plotLabels} 3460 # 3461 self.emiss(wvlRange = wvlRange, allLines=allLines) 3462 emiss = self.Emiss 3463 if 'errorMessage' in emiss.keys(): 3464 self.Intensity = {'errorMessage': self.Spectroscopic+' no lines in this wavelength region'} 3465 return 3466 em = emiss['emiss'] 3467 wvl = emiss['wvl'] 3468 eDensity = self.EDensity 3469 if hasattr(self, 'Abundance'): 3470 ab=self.Abundance 3471 else: 3472 self.Abundance = io.abundanceRead() 3473 ab=self.Abundance 3474 if hasattr(self, 'IoneqOne'): 3475 thisIoneq=self.IoneqOne 3476 else: 3477 self.ioneqOne() 3478 thisIoneq=self.IoneqOne 3479 try: 3480 nwvl, ntempden = em.shape 3481 intensity = np.zeros((ntempden, nwvl),'Float64') 3482 if thisIoneq.size == 1: 3483 thisIoneq = np.ones(ntempden, 'float64')*thisIoneq 3484 for it in range(ntempden): 3485 if self.Defaults['flux'] != 'energy': 3486 intensity[it] = 4.*const.pi*(const.planck*const.light*1.e+8/wvl)*ab*thisIoneq[it]*em[:, it] 3487 else: 3488 intensity[it] = 4.*const.pi*ab*thisIoneq[it]*em[:, it]/eDensity[it] 3489 loss = intensity.sum(axis=1) 3490 except: 3491 nwvl=len(em) 3492 ntempden=1 3493 if self.Defaults['flux'] != 'energy': 3494 intensity = 4.*const.pi*(const.planck*const.light*1.e+8/wvl)*ab*thisIoneq*em 3495 else: 3496 intensity = 4.*const.pi*ab*thisIoneq*em 3497 loss = intensity.sum() 3498 self.BoundBoundLoss = {'rate':loss, 'wvlRange':wvlRange, 'temperature':self.Temperature, 'eDensity':self.EDensity}

3499 # 3500 # ------------------------------------------------------------------------------------- 3501 #

3502 - def intensityRatioInterpolate(self,data, scale = 'lin', plot=0, verbose=0):

3503 ''' 3504 to take a set of date and interpolate against the IntensityRatio 3505 the scale can be one of 'lin'/'linear' [default], 'loglog', 'logx', 'logy', 3506 ''' 3507 # first, what variable to use 3508 if self.IntensityRatio['temperature'].max() > self.IntensityRatio['temperature'].min(): 3509 x = self.IntensityRatio['ratio'] 3510 y = self.IntensityRatio['temperature'] 3511 if verbose: 3512 print('using temperature with %i5 values'%(len(x))) 3513 print(' number of values') 3514 else: 3515 x = self.IntensityRatio['ratio'] 3516 y = self.IntensityRatio['eDensity'] 3517 # 3518 if x[0] > x[-1]: 3519 x = sorted(x) 3520 sy = [] 3521 for idx in range(len(y) -1, -1, -1): 3522 sy.append(y[idx]) 3523 else: 3524 sy = y 3525 # 3526 if 'lin' in scale: 3527 y2 = interpolate.splrep(x, sy, s=0) 3528 interpolatedData = interpolate.splev(data,y2) 3529 if plot: 3530 pl.plot(sy, x) 3531 pl.plot(interpolatedData, data, 'bD') 3532 elif scale == 'loglog': 3533 y2 = interpolate.splrep(np.log(x), np.log(sy), s=0) 3534 interpolatedData = np.exp(interpolate.splev(np.log(data),y2)) 3535 if plot: 3536 pl.loglog(sy, x) 3537 pl.loglog(interpolatedData, data, 'bD') 3538 elif scale == 'logx': 3539 y2 = interpolate.splrep(x, np.log(sy), s=0) 3540 interpolatedData = np.exp(interpolate.splev(data,y2)) 3541 if plot: 3542 pl.semilogx(sy, x) 3543 pl.semilogx(interpolatedData, data, 'bD') 3544 elif scale == 'logy': 3545 y2 = interpolate.splrep(np.log(x), sy, s=0) 3546 interpolatedData = interpolate.splev(np.log(data),y2) 3547 if plot: 3548 pl.semilogy(sy, x) 3549 pl.semilogy(interpolatedData, data, 'bD') 3550 else: 3551 print(' scale not understood = %s'%(scale)) 3552 for i, avalue in enumerate(interpolatedData): 3553 print(' data, value = %12.3e %12.3e'%(data[i], avalue)) 3554 self.IntensityRatioInterpolated = {'data':data, 'value':interpolatedData}

3555 # 3556 # ------------------------------------------------------------------------------------- 3557 #

3558 - def ioneqOne(self):

3559 ''' 3560 Provide the ionization equilibrium for the selected ion as a function of temperature. 3561 returned in self.IoneqOne 3562 ''' 3563 # 3564 if hasattr(self, 'Temperature'): 3565 temperature = self.Temperature 3566 else: 3567 return 3568 # 3569 if hasattr(self, 'IoneqAll'): 3570 ioneqAll = self.IoneqAll 3571 else: 3572 ioneqAll = io.ioneqRead(ioneqname = self.Defaults['ioneqfile']) 3573 # if chInteractive: 3574 self.ioneqAll=self.IoneqAll 3575 # 3576 ioneqTemperature = ioneqAll['ioneqTemperature'] 3577 Z=self.Z 3578 Ion=self.Ion 3579 Dielectronic=self.Dielectronic 3580 ioneqOne = np.zeros_like(temperature) 3581 # 3582 thisIoneq=ioneqAll['ioneqAll'][Z-1,Ion-1 + Dielectronic].squeeze() 3583 del ioneqAll 3584 # thisIoneq = self.Ioneq 3585 gioneq=thisIoneq > 0. 3586 goodt1=self.Temperature >= ioneqTemperature[gioneq].min() 3587 goodt2=self.Temperature <= ioneqTemperature[gioneq].max() 3588 goodt=np.logical_and(goodt1,goodt2) 3589 y2=interpolate.splrep(np.log(ioneqTemperature[gioneq]),np.log(thisIoneq[gioneq]),s=0) 3590 # 3591 if goodt.sum() > 0: 3592 if self.Temperature.size > 1: 3593 gIoneq=interpolate.splev(np.log(self.Temperature[goodt]),y2) #,der=0) 3594 ioneqOne[goodt] = np.exp(gIoneq) 3595 else: 3596 gIoneq=interpolate.splev(np.log(self.Temperature),y2) 3597 ioneqOne = np.exp(gIoneq) 3598 self.IoneqOne = ioneqOne

3599 # 3600 # ------------------------------------------------------------------------------------- 3601 #

3602 - def gofnt(self,wvlRange=0,top=10, verbose=0):

3603 """ 3604 Calculate the 'so-called' G(T) function. 3605 3606 Given as a function of both temperature and eDensity. 3607 3608 Only the top( set by 'top') brightest lines are plotted. 3609 the G(T) function is returned in a dictionary self.Gofnt 3610 """ 3611 # 3612 #self.emiss={"wvl":wvl,"emiss":em,"units":units,"plotLabels":plotLabels} 3613 # 3614 # 3615 if hasattr(self, 'Emiss'): 3616 em=copy.copy(self.Emiss) 3617 else: 3618 self.emiss() 3619 em=copy.copy(self.Emiss) 3620 # 3621 # 3622 # if hasattr(self, 'Abundance'): 3623 # ab=self.Abundance 3624 # else: 3625 # self.Abundance = io.abundanceRead() 3626 # ab=self.Abundance 3627 if not hasattr(self, 'Abundance'): 3628 self.Abundance = io.abundanceRead() 3629 # 3630 fontsize=12 3631 # 3632 emiss=em["emiss"] 3633 wvl=em["wvl"] 3634 pretty1 = em['pretty1'] 3635 pretty2 = em['pretty2'] 3636 lvl1 = em['lvl1'] 3637 lvl2 = em['lvl2'] 3638 # 3639 # temperature=self.Temperature 3640 # eDensity=self.EDensity 3641 plotLabels = em["plotLabels"] 3642 xLabel = plotLabels["xLabel"] 3643 yLabel = plotLabels["yLabel"] 3644 # 3645 # find which lines are in the wavelength range if it is set 3646 # 3647 # 3648 if type(wvlRange) != type(1): 3649 igvl=util.between(wvl,wvlRange) 3650 else: 3651 igvl=range(len(wvl)) 3652 nlines=len(igvl) 3653 if nlines ==0: 3654 print(' no lines in selected interval') 3655 return 3656 # find the top most intense lines 3657 # 3658 if top > nlines: 3659 top=nlines 3660 maxEmiss=np.zeros(nlines,'Float64') 3661 for iline in range(nlines): 3662 maxEmiss[iline]=emiss[igvl[iline]].max() 3663 for iline in range(nlines): 3664 if maxEmiss[iline]>=maxEmiss.max(): 3665 maxAll=emiss[igvl[iline]] 3666 # maxIndex = igvl[iline] 3667 # print ' maxIndex, maxAll = ', maxIndex, maxAll 3668 # line=range(nlines) 3669 igvlsort=np.take(igvl,np.argsort(maxEmiss)) 3670 topLines=igvlsort[-top:] 3671 maxWvl='%5.3f' % wvl[topLines[-1]] 3672 # maxline=topLines[-1] 3673 # 3674 # need to make sure there are no negative values before plotting 3675 good = np.where(emiss > 0.) 3676 emissMin=emiss[good].min() 3677 bad=np.where(emiss <= 0.) 3678 emiss[bad]=emissMin 3679 # 3680 topLines=topLines[wvl[topLines].argsort()] 3681 # 3682 # 3683 eDensity = self.EDensity 3684 temperature = self.Temperature 3685 # 3686 # temp=np.asarray(temperature,'Float32') 3687 ntemp = temperature.size 3688 if ntemp > 0: 3689 if temperature[0] == temperature[-1]: 3690 ntemp = 1 3691 # 3692 ndens = eDensity.size 3693 if ndens > 0: 3694 if eDensity[0] == eDensity[-1]: 3695 ndens = 1 3696 # 3697 print(' ndens = %5i ntemp = %5i'%(ndens, ntemp)) 3698 # 3699 ylabel = 'Emissivity relative to '+maxWvl 3700 title = self.Spectroscopic 3701 # 3702 # 3703 if ndens==1 and ntemp==1: 3704 print(' only a single temperature and eDensity') 3705 return 3706 elif ndens == 1: 3707 xlabel='Temperature (K)' 3708 ngofnt = temperature.size 3709 xvalues=temperature 3710 outTemperature=temperature 3711 outDensity = eDensity 3712 desc_str=' at Density = %10.2e' % eDensity[0] 3713 elif ntemp == 1: 3714 xvalues=eDensity 3715 ngofnt = eDensity.size 3716 outTemperature=np.zeros(ndens,'Float64') 3717 outTemperature.fill(temperature) 3718 outDensity=eDensity 3719 xlabel=r'$\rm{Electron Density (cm}^{-3}\rm{)}$' 3720 desc_str=' at Temperature = %10.2e' % temperature[0] 3721 else: 3722 outTemperature=temperature 3723 outDensity=eDensity 3724 xlabel='Temperature (K)' 3725 xvalues=temperature 3726 ngofnt = ntemp 3727 desc_str=' for variable Density' 3728 # 3729 # 3730 # put all actual plotting here 3731 # 3732 pl.ion() 3733 # if chInteractive: 3734 # pl.ion() 3735 # else: 3736 # pl.ioff() 3737 # 3738 pl.figure() 3739 ax = pl.subplot(111) 3740 nxvalues=len(xvalues) 3741 # maxAll is an array 3742 # print ' emiss = ', np.max(emiss[top-1]), np.max(emiss[0]) 3743 # print ' maxAll = ', maxAll 3744 # ymax = np.max(1.2*emiss[top-1]/maxAll) 3745 ymax = 1.2 3746 # print ' ymax = ', ymax 3747 ymin = ymax # np.min(emiss[0]/maxAll) # was originally = ymax 3748 for iline in range(top): 3749 tline=topLines[iline] 3750 pl.loglog(xvalues,emiss[tline]/maxAll) 3751 if np.min(emiss[tline]/maxAll) < ymin: 3752 ymin = np.min(emiss[tline]/maxAll) 3753 skip=2 3754 start=divmod(iline,nxvalues)[1] 3755 for ixvalue in range(start,nxvalues,nxvalues//skip): 3756 pl.text(xvalues[ixvalue],emiss[tline,ixvalue]/maxAll[ixvalue],str(wvl[tline])) 3757 pl.xlim(xvalues.min(),xvalues.max()) 3758 pl.ylim(ymin, ymax) 3759 # yl=pl.ylim() 3760 # pl.ylim(yl[0],1.2) 3761 pl.xlabel(xlabel,fontsize=fontsize) 3762 pl.ylabel(ylabel,fontsize=fontsize) 3763 if ndens == ntemp and ntemp > 1: 3764 pl.text(0.07, 0.5,title, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 3765 # 3766 ax2 = pl.twiny() 3767 xlabelDen=r'Electron Density (cm$^{-3}$)' 3768 pl.xlabel(xlabelDen, fontsize=fontsize) 3769 pl.loglog(eDensity,emiss[topLines[top-1]]/maxAll, visible=False) 3770 ax2.xaxis.tick_top() 3771 else: 3772 pl.ylim(ymin, ymax) 3773 pl.title(title+desc_str,fontsize=fontsize) 3774 pl.draw() 3775 # 3776 time.sleep(0.5) 3777 # 3778 # print ' topLInes = ', wvl[topLines] 3779 wvlChoices = [] 3780 for iline in range(top): 3781 tline = topLines[iline] 3782 wvlChoices.append('%12.4f %4i %4i %s - %s'%(wvl[tline], lvl1[tline], lvl2[tline], pretty1[tline], pretty2[tline])) 3783 gline = chGui.gui.selectorDialog(wvlChoices,label='Select line(s)') 3784 gline_idx=gline.selectedIndex 3785 # 3786 # 3787 gAbund=self.Abundance 3788 # 3789 try: 3790 thisIoneq=self.IoneqOne 3791 except: 3792 self.ioneqOne() 3793 thisIoneq=self.IoneqOne 3794 if verbose: 3795 print(' abundance = %10.2e '%(gAbund)) 3796 print(' index temperature ion fraction') 3797 for it, anioneq in enumerate(thisIoneq): 3798 print (' %5i %10.2e %10.2e '%(it, outTemperature[it], anioneq)) 3799 # gioneq=np.where(thisIoneq > 0.) 3800 # y2=interpolate.splrep(np.log(self.IoneqAll['ioneqTemperature'][gioneq]),np.log(thisIoneq[gioneq]),s=0) #allow smoothing,s=0) 3801 # gIoneq=interpolate.splev(np.log(temperature),y2) #,der=0) 3802 gIoneq=self.IoneqOne/eDensity 3803 # 3804 # 3805 # 3806 # plot the desired ratio 3807 pl.figure() 3808 g_line = topLines[gline_idx]# [0] 3809 ## print ' g_line = ',g_line 3810 # 3811 gofnt=np.zeros(ngofnt,'float64') 3812 for aline in g_line: 3813 gofnt+=gAbund*gIoneq*emiss[aline].squeeze() 3814 self.Gofnt={'temperature':outTemperature,'eDensity':outDensity,'gofnt':gofnt, 'index':g_line, 'wvl':wvl[g_line]} 3815 # 3816 pl.loglog(xvalues,gofnt) 3817 pl.xlim(xvalues.min(),xvalues.max()) 3818 pl.xlabel(xlabel,fontsize=fontsize) 3819 pl.ylabel('Gofnt',fontsize=fontsize) 3820 newTitle = '%9s'%(self.Spectroscopic) + '%12.3f %4i %4i %s - %s'%(wvl[g_line[0]], lvl1[g_line[0]], lvl2[g_line[0]], pretty1[g_line[0]], pretty2[g_line[0]]) 3821 if len(g_line) > 1: 3822 newTitle +='\n' 3823 for igl in g_line[1:]: 3824 newTitle += ' ' + '%12.3f %4i %4i %s - %s'%(wvl[igl], lvl1[igl], lvl2[igl], pretty1[igl], pretty2[igl]) 3825 if igl != g_line[-1]: 3826 newTitle +='\n' 3827 # pl.annotate(newTitle, xytext=(0.3, 0.3), textcoords='figure_fraction') 3828 pl.annotate(newTitle, xy=(-10, 10), 3829 xycoords='axes points', 3830 horizontalalignment='right', verticalalignment='bottom') #,fontsize=20) 3831 if ndens == ntemp and ntemp > 1: 3832 # newTitle +' '+str(wvl[g_line])+' '+desc_str 3833 pl.text(0.07, 0.5,newTitle, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 3834 # 3835 ax2 = pl.twiny() 3836 # xlabel=r'Electron Density (cm$^{-3}$)' 3837 pl.xlabel(xlabelDen, fontsize=fontsize) 3838 pl.loglog(eDensity,gofnt, visible=False) 3839 ax2.xaxis.tick_top() 3840 else: 3841 pl.title(newTitle, fontsize=fontsize)

3842 #pl.ioff() 3843 #pl.show() 3844 # return 3845 # 3846 # - - - - - - - - - - - - - - - - - - - - - - - 3847 #

3848 - def twoPhotonEmiss(self, wvl):

3849 ''' to calculate the two-photon continuum rate coefficient - only for hydrogen- and helium-like ions''' 3850 wvl = np.array(wvl, 'float64') 3851 nWvl = wvl.size 3852 if self.Z -self.Ion > 1 or self.Dielectronic: 3853 # this is not a hydrogen-like or helium-like ion 3854 self.TwoPhoton = {'emiss':np.zeros(nWvl, 'float4'), 'wvl':wvl} 3855 return 3856 else: 3857 try: 3858 pop = self.Population['population'] 3859 nTempDens = max(self.Temperature.size, self.EDensity.size) 3860 except: 3861 self.populate() 3862 pop = self.Population['population'] 3863 nTempDens = max(self.Temperature.size, self.EDensity.size) 3864 if nTempDens > 1: 3865 emiss = np.zeros((nTempDens, nWvl), 'float64') 3866 if self.EDensity.size == 1: 3867 eDensity = np.repeat(self.EDensity, nTempDens) 3868 else: 3869 eDensity = self.EDensity 3870 else: 3871 emiss = np.zeros(nWvl, 'float64') 3872 eDensity = self.EDensity 3873 if self.Z == self.Ion: 3874 # H seq 3875 l1 = 1-1 3876 l2 = 2 - 1 3877 wvl0 = 1.e+8/(self.Elvlc['ecm'][l2] - self.Elvlc['ecm'][l1]) 3878 goodWvl = wvl > wvl0 3879 y = wvl0/wvl[goodWvl] 3880 dist = io.twophotonHRead() 3881 avalue = dist['avalue'][self.Z-1] 3882 asum = dist['asum'][self.Z-1] 3883 distr1 = interpolate.splrep(dist['y0'], dist['psi0'][self.Z-1], s=0) 3884 distr = avalue*y*interpolate.splev(y, distr1)/(asum*wvl[goodWvl]) 3885 if self.Defaults['flux'] == 'energy': 3886 f = (const.light*const.planck*1.e+8)/wvl[goodWvl] 3887 else: 3888 f=1. 3889 if nTempDens == 1: 3890 emiss[goodWvl] = f*pop[l2]*distr/self.EDensity 3891 else: 3892 for it in range(nTempDens): 3893 emiss[it, goodWvl] = f*pop[it, l2]*distr/self.EDensity[it] 3894 self.TwoPhotonEmiss = {'wvl':wvl, 'emiss':emiss} 3895 else: 3896 # He seq 3897 l1 = 1-1 3898 l2 = heseqLvl2[self.Z -1] - 1 3899 wvl0 = 1.e+8/(self.Elvlc['ecm'][l2] - self.Elvlc['ecm'][l1]) 3900 goodWvl = wvl > wvl0 3901 y = wvl0/wvl[goodWvl] 3902 dist = io.twophotonHeRead() 3903 avalue = dist['avalue'][self.Z-1] 3904 distr1 = interpolate.splrep(dist['y0'], dist['psi0'][self.Z-1], s=0) 3905 distr = avalue*y*interpolate.splev(y, distr1)/wvl[goodWvl] 3906 if self.Defaults['flux'] == 'energy': 3907 f = (const.light*const.planck*1.e+8)/wvl[goodWvl] 3908 else: 3909 f=1. 3910 if nTempDens == 1: 3911 emiss[goodWvl] = f*pop[l2]*distr/self.EDensity 3912 else: 3913 for it in range(nTempDens): 3914 emiss[it, goodWvl] = f*pop[it, l2]*distr/self.EDensity[it] 3915 self.TwoPhotonEmiss = {'wvl':wvl, 'emiss':emiss}

3916 # 3917 #----------------------------------------------------------------- 3918 #

3919 - def twoPhoton(self, wvl, em=0):

3920 ''' 3921 to calculate the two-photon continuum - only for hydrogen- and helium-like ions 3922 includes the elemental abundance and the ionization equilibrium 3923 includes the emission measure if specified 3924 ''' 3925 wvl = np.array(wvl, 'float64') 3926 # 3927 if type(em) == int and em == 0: 3928 if hasattr(self, 'Em'): 3929 em = self.Em 3930 else: 3931 em = np.ones(self.NTempDen, 'float64') 3932 self.Em = em 3933 elif type(em) == float and em > 0.: 3934 em = np.ones(self.NTempDen, 'float64')*em 3935 self.Em = em 3936 elif type(em) == list or type(em) == tuple or type(em) == np.ndarray: 3937 em = np.asarray(em, 'float64') 3938 self.Em = em 3939 # so we know that it has been applied 3940 # 3941 nWvl = wvl.size 3942 if self.Z - self.Ion > 1 or self.Dielectronic: 3943 # this is not a hydrogen-like or helium-like ion 3944 print(' not doing 2 photon for %s'%(self.IonStr)) 3945 self.TwoPhoton = {'emiss':np.zeros(nWvl, 'float64'), 'wvl':wvl} 3946 return 3947 else: 3948 try: 3949 ab=self.Abundance 3950 except: 3951 self.Abundance = io.abundanceRead() 3952 ab=self.Abundance 3953 try: 3954 thisIoneq=self.IoneqOne 3955 except: 3956 self.ioneqOne() 3957 thisIoneq=self.IoneqOne 3958 try: 3959 pop = self.Population['population'] 3960 nTempDens = max(self.Temperature.size, self.EDensity.size) 3961 except: 3962 self.populate() 3963 pop = self.Population['population'] 3964 nTempDens = max(self.Temperature.size, self.EDensity.size) 3965 if nTempDens > 1: 3966 rate = np.zeros((nTempDens, nWvl), 'float64') 3967 if em.size == 1.: 3968 em = np.ones(nTempDens, 'float64')*em 3969 if self.EDensity.size == 1: 3970 eDensity = np.repeat(self.EDensity, nTempDens) 3971 else: 3972 eDensity = self.EDensity 3973 if self.Temperature.size == 1: 3974 temperature = np.repeat(self.Temperature, nTempDens) 3975 thisIoneq = np.repeat(thisIoneq, nTempDens) 3976 else: 3977 temperature = self.Temperature 3978 else: 3979 rate = np.zeros(nWvl, 'float64') 3980 eDensity = self.EDensity 3981 temperature = self.Temperature 3982 if self.Z == self.Ion: 3983 # H seq 3984 l1 = 1-1 3985 l2 = 2 - 1 3986 wvl0 = 1.e+8/(self.Elvlc['ecm'][l2] - self.Elvlc['ecm'][l1]) 3987 goodWvl = wvl > wvl0 3988 if goodWvl.sum() > 0: 3989 y = wvl0/wvl[goodWvl] 3990 dist = io.twophotonHRead() 3991 avalue = dist['avalue'][self.Z-1] 3992 asum = dist['asum'][self.Z-1] 3993 distr1 = interpolate.splrep(dist['y0'], dist['psi0'][self.Z-1], s=0) 3994 distr = avalue*y*interpolate.splev(y, distr1)/(asum*wvl[goodWvl]) 3995 if self.Defaults['flux'] == 'energy': 3996 f = (const.light*const.planck*1.e+8)/(4.*const.pi*wvl[goodWvl]) 3997 else: 3998 f=1./(4.*const.pi) 3999 if nTempDens == 1: 4000 rate[goodWvl] = f*pop[l2]*distr*ab*thisIoneq*em/eDensity 4001 else: 4002 for it in range(nTempDens): 4003 rate[it, goodWvl] = f*pop[it, l2]*distr*ab*thisIoneq[it]*em[it]/eDensity[it] 4004 self.TwoPhoton = {'wvl':wvl, 'rate':rate} 4005 4006 else: 4007 # He seq 4008 l1 = 1-1 4009 l2 = heseqLvl2[self.Z -1] - 1 4010 wvl0 = 1.e+8/(self.Elvlc['ecm'][l2] - self.Elvlc['ecm'][l1]) 4011 goodWvl = wvl > wvl0 4012 if goodWvl.sum() > 0: 4013 y = wvl0/wvl[goodWvl] 4014 dist = io.twophotonHeRead() 4015 avalue = dist['avalue'][self.Z-1] 4016 distr1 = interpolate.splrep(dist['y0'], dist['psi0'][self.Z-1], s=0) 4017 distr = avalue*y*interpolate.splev(y, distr1)/wvl[goodWvl] 4018 if self.Defaults['flux'] == 'energy': 4019 f = (const.light*const.planck*1.e+8)/(4.*const.pi*wvl[goodWvl]) 4020 else: 4021 f=1./(4.*const.pi) 4022 if nTempDens == 1: 4023 rate[goodWvl] = f*pop[l2]*distr*ab*thisIoneq*em/eDensity 4024 else: 4025 for it in range(nTempDens): 4026 rate[it, goodWvl] = f*pop[it, l2]*distr*ab*thisIoneq[it]*em[it]/eDensity[it] 4027 self.TwoPhoton = {'wvl':wvl, 'rate':rate, 'em':em}

4028 # 4029 #----------------------------------------------------------------- 4030 #

4031 - def twoPhotonLoss(self):

4032 ''' 4033 to calculate the two-photon energy loss rate - only for hydrogen- and helium-like ions 4034 includes the elemental abundance and the ionization equilibrium 4035 does not include the emission measure 4036 ''' 4037 if self.Z -self.Ion > 1 or self.Dielectronic: 4038 # this is not a hydrogen-like or helium-like ion 4039 nTempDens = max(self.Temperature.size, self.EDensity.size) 4040 # if nTempDens > 1: 4041 rate = np.zeros((nTempDens), 'float64') 4042 self.TwoPhotonLoss = {'rate':rate} 4043 else: 4044 if hasattr(self, 'Abundance'): 4045 ab=self.Abundance 4046 else: 4047 self.Abundance = io.abundanceRead() 4048 ab=self.Abundance 4049 if hasattr(self, 'IoneqOne'): 4050 thisIoneq=self.IoneqOne 4051 else: 4052 self.ioneqOne() 4053 thisIoneq=self.IoneqOne 4054 if hasattr(self, 'Population'): 4055 pop = self.Population['population'] 4056 nTempDens = max(self.Temperature.size, self.EDensity.size) 4057 else: 4058 self.populate() 4059 pop = self.Population['population'] 4060 nTempDens = max(self.Temperature.size, self.EDensity.size) 4061 if nTempDens > 1: 4062 rate = np.zeros((nTempDens), 'float64') 4063 if self.EDensity.size == 1: 4064 eDensity = np.repeat(self.EDensity, nTempDens) 4065 else: 4066 eDensity = self.EDensity 4067 else: 4068 eDensity = self.EDensity 4069 if self.Z == self.Ion: 4070 # H seq 4071 l1 = 1 - 1 4072 l2 = 2 - 1 4073 wvl0 = 1.e+8/(self.Elvlc['ecm'][l2] - self.Elvlc['ecm'][l1]) 4074 dist = io.twophotonHRead() 4075 avalue = dist['avalue'][self.Z-1] 4076 f = (avalue*const.light*const.planck*1.e+8)/wvl0 4077 if nTempDens == 1: 4078 rate = f*pop[l2]*ab*thisIoneq/eDensity 4079 else: 4080 for it in range(nTempDens): 4081 rate[it] = f*pop[it, l2]*ab*thisIoneq[it]/eDensity[it] 4082 self.TwoPhotonLoss = {'temperature':self.Temperature,'eDensity':self.EDensity,'rate':rate} 4083 else: 4084 # He seq 4085 l1 = 1 - 1 4086 l2 = heseqLvl2[self.Z -1] -1 4087 wvl0 = 1.e+8/(self.Elvlc['ecm'][l2] - self.Elvlc['ecm'][l1]) 4088 dist = io.twophotonHeRead() 4089 avalue = dist['avalue'][self.Z-1] 4090 f = (avalue*const.light*const.planck*1.e+8)/wvl0 4091 if nTempDens == 1: 4092 rate = f*pop[l2]*ab*thisIoneq/eDensity 4093 else: 4094 for it in range(nTempDens): 4095 rate[it] = f*pop[it, l2]*ab*thisIoneq[it]/eDensity[it] 4096 self.TwoPhotonLoss = {'temperature':self.Temperature,'eDensity':self.EDensity,'rate':rate}

Source Code for Module chianti.core.Ion