chianti.util

1 '''Utility functions, many for reading the CHIANTI database files. 2 3 Copyright 2009, 2010 Kenneth P. Dere 4 5 This software is distributed under the terms of the GNU General Public License 6 that is found in the LICENSE file 7 8 9 ''' 10 import os, fnmatch 11 #from types import * 12 #try: 13 # # for Python 3 import 14 # import configparser 15 #except ImportError: 16 # # for Python 2 import 17 # import ConfigParser as configparser 18 ##from ConfigParser import * 19 import pickle 20 from datetime import date 21 import numpy as np 22 from scipy import interpolate 23 import chianti 24 import chianti.constants as const 25 #import chianti.io as chio 26 #import chianti.gui as gui 27 # 28 #

29 -def between(array,limits):

30 ''' 31 returns an index array of elements of array where the values lie 32 between the limits given as a 2 element list or tuple 33 ''' 34 array=np.asarray(array) 35 nlines=len(array) 36 # hi=np.where(array >= limits[0],range(1,nlines+1),0) 37 # lo=np.where(array <= limits[1],range(1,nlines+1),0) 38 hi=np.where(array >= limits[0],list(range(1,nlines+1)),0) 39 lo=np.where(array <= limits[1],list(range(1,nlines+1)),0) 40 41 hilo=hi&lo 42 out=[a -1 for a in hilo if a > 0] 43 return out

44 # 45 # -------------------------------------------------- 46 #

47 -def z2element(z):

48 """ convert Z to element string """ 49 if z-1 < len(const.El): 50 thisel=const.El[z-1] 51 else: 52 thisel='' 53 return thisel

54 # 55 # ------------------------------------------------------------------------------------- 56 #

57 -def spectroscopic2name(el,roman, dielectronic=False):

58 """ 59 convert Z and ion to spectroscopic notation string 60 """ 61 elu = el.lower() 62 romanu = roman.upper() 63 idx =const.Ionstage.index(romanu) 64 gname = elu+'_'+str(idx+1) 65 return gname

66 # 67 # ------------------------------------------------------------------------------------- 68 #

69 -def zion2name(z,ion, dielectronic=False):

70 """ 71 convert Z, ion to generic name 26, 13 -> fe_13 72 """ 73 if ion == 0: 74 thisone = 0 75 elif ion == z+2: 76 thisone = 0 77 elif (z-1 < len(const.El)) and (ion <= z+1): 78 thisone=const.El[z-1]+'_'+str(ion) 79 if dielectronic: 80 thisone+='d' 81 else: 82 # this should not actually happen 83 thisone = 0 84 return thisone

85 # 86 # ------------------------------------------------------------------------------------- 87 #

88 -def zion2dir(z,ion, dielectronic=False, xuvtop=0):

89 """ convert Z to generic file name string """ 90 if xuvtop: 91 dir = xuvtop 92 else: 93 dir=os.environ["XUVTOP"] 94 if (z-1 < len(const.El)) and (ion <= z+1): 95 thisel=const.El[z-1] 96 else: 97 thisel='' 98 if z-1 < len(const.El): 99 thisone=const.El[z-1]+'_'+str(ion) 100 if dielectronic: 101 thisone+='d' 102 else: 103 thisone='' 104 if thisel != '' : 105 fname=os.path.join(dir,thisel,thisone) 106 return fname

107 # 108 # ------------------------------------------------------------------------------------- 109 #

110 -def zion2filename(z,ion, dielectronic=False, xuvtop=0):

111 """ convert Z to generic file name string """ 112 if xuvtop: 113 dir = xuvtop 114 else: 115 dir=os.environ["XUVTOP"] 116 if (z-1 < len(const.El)) and (ion <= z+1): 117 thisel=const.El[z-1] 118 else: 119 thisel='' 120 if z-1 < len(const.El): 121 thisone=const.El[z-1]+'_'+str(ion) 122 if dielectronic: 123 thisone+='d' 124 else: 125 thisone='' 126 if thisel != '' : 127 fname=os.path.join(dir,thisel,thisone,thisone) 128 return fname

129 # 130 # ------------------------------------------------------------------------------------- 131 #

132 -def zion2localFilename(z,ion, dielectronic=False):

133 """ convert Z to generic file name string with current directory at top""" 134 dir='.' 135 if (z-1 < len(const.El)) and (ion <= z+1): 136 thisel=const.El[z-1] 137 else: 138 thisel='' 139 if z-1 < len(const.El): 140 thisone=const.El[z-1]+'_'+str(ion) 141 if dielectronic: 142 thisone+='d' 143 else: 144 thisone='' 145 if thisel != '' : 146 fname=os.path.join(dir,thisel,thisone,thisone) 147 return fname

148 # 149 # ------------------------------------------------------------------------------------- 150 #

151 -def zion2spectroscopic(z,ion, dielectronic=False):

152 """ 153 convert Z and ion to spectroscopic notation string 154 """ 155 if (z-1 < len(const.El)) and (ion <= z+1): 156 spect=const.El[z-1].capitalize()+' '+const.Ionstage[ion-1] 157 if dielectronic: 158 spect+=' d' 159 else: spect = '' 160 return spect

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

164 -def convertName(name):

165 """ convert ion name string to Z and Ion """ 166 s2=name.split('_') 167 els=s2[0].strip() 168 i1=const.El.index(els)+1 169 ions=s2[1].strip() 170 d=ions.find('d') 171 if d >0 : 172 dielectronic=True 173 ions=ions.replace('d','') 174 else: dielectronic=False 175 higher = zion2name(int(i1), int(ions)+1) 176 lower = zion2name(int(i1), int(ions)-1) 177 return {'Z':int(i1),'Ion':int(ions),'Dielectronic':dielectronic, 'Element':els, 'higher':higher, 'lower':lower}

178 # 179 # ------------------------------------------------------------------------------------- 180 #

181 -def ion2filename(ions):

182 """ 183 convert ion string to generic file name string 184 """ 185 dir=os.environ["XUVTOP"] 186 zion=convertName(ions) 187 el=z2element(zion['Z']) 188 fname=os.path.join(dir,el,ions,ions) 189 return fname

190 # 191 # ------------------------------------------------------------------------------------- 192 #

193 -def el2z(els):

194 """ 195 from an the name of the element (1-2 letter) return Z 196 """ 197 z=const.El.index(els.lower())+1 198 return z

199 # 200 # ------------------------------------------------------------------------------------- 201 #

202 -def qrp(z,u):

203 ''' qrp(Z,u) u = E/IP 204 calculate Qr-prime (equ. 2.12) of Fontes, Sampson and Zhang 1999''' 205 # 206 aa=1.13 # aa stands for A in equ 2.12 207 # 208 if z >= 16 : 209 # use Fontes Z=20, N=1 parameters 210 dd=3.70590 211 c=-0.28394 212 d=1.95270 213 cc=0.20594 214 else: 215 # use Fontes Z=10, N=2 parameters 216 dd=3.82652 217 c=-0.80414 218 d=2.32431 219 cc=0.14424 220 # 221 if z > 20: 222 cc+=((z-20.)/50.5)**1.11 223 # 224 bu=u <= 1. 225 q=np.ma.array(u, 'Float64', mask=bu, fill_value=0.) 226 # 227 # 228 q=(aa*np.ma.log(u) + dd*(1.-1./u)**2 + cc*u*(1.-1./u)**4 + (c/u+d/u**2)*(1.-1/u))/u 229 # 230 q.set_fill_value(0.) # I don't know why this is necessary 231 return q # .set_fill_value(0.)

232 # 233 #----------------------------------------------------------- 234 #

235 -def splomDescale(splom, energy):

236 """ 237 Calculates the collision strength 238 for excitation-autoionization as a function of energy. 239 energy in eV 240 """ 241 # 242 # 243 nenergy=energy.size 244 nsplom=len(splom['deryd']) 245 # for these files, there are 5 spline points 246 nspl = 5 247 if nenergy > 1: 248 omega = np.zeros((nsplom,nenergy),"float64") 249 else: 250 omega = np.zeros(nsplom,"float64") 251 # 252 dx = 1./(float(nspl)-1.) 253 sxint = dx*np.arange(nspl) 254 for isplom in range(0,nsplom): 255 # 256 sx1 = energy/(splom['deryd'][isplom]*const.ryd2Ev) 257 good = sx1 >= 1. 258 # make sure there are some valid energies above the threshold 259 if good.sum(): 260 nbad = nenergy - good.sum() 261 c_curr = splom['c'][isplom] 262 # 263 if splom['ttype'][isplom] == 1: 264 sx = 1. - np.log(c_curr)/np.log(sx1[good] - 1. + c_curr) 265 y2 = interpolate.splrep(sxint,splom['splom'][:, isplom],s=0) #allow smoothing,s=0) 266 som = interpolate.splev(sx,y2,der=0) 267 omega[isplom, nbad:] = som*np.log(sx -1. + np.exp(1.)) 268 # 269 elif splom['ttype'][isplom] == 2: 270 sx =(sx1[good] - 1.)/(sx1[good] -1. + c_curr) 271 y2 = interpolate.splrep(sxint,splom['splom'][:, isplom],s=0) #allow smoothing,s=0) 272 som=interpolate.splev(sx,y2,der=0) 273 omega[isplom, nbad:] = som 274 # 275 elif splom['ttype'][isplom] == 3: 276 sx = (sx1[good] - 1.)/(sx1[good] -1. + c_curr) 277 y2 = interpolate.splrep(sxint,splom['splom'][:, isplom],s=0) #allow smoothing,s=0) 278 som = interpolate.splev(sx,y2,der=0) 279 omega[isplom, nbad:] = som/sx1[good]**2 280 # 281 elif splom['ttype'][isplom] == 4: 282 sx = 1. - np.log(c_curr)/np.log(sx1[good] -1. + c_curr) 283 y2 = interpolate.splrep(sxint,splom['splom'][:, isplom],s=0) #allow smoothing,s=0) 284 som=interpolate.splev(sx,y2,der=0) 285 omega[isplom, nbad:] = som*np.log(sx1[good] -1. + c_curr) 286 # 287 # 288 # 289 elif ttype > 4: 290 print((' splom t_type ne 1,2,3,4 = %4i %4i %4i'%(ttype,l1,l2))) 291 else: 292 # there are no energies above the threshold 293 pass 294 # 295 # 296 omega=np.where(omega > 0.,omega,0.) 297 # 298 return omega

299 # 300 # ------------------------------------------------------------------------------ 301 #

302 -def dilute(radius):

303 ''' 304 to calculate the dilution factor as a function distance from 305 the center of a star in units of the stellar radius 306 a radius of less than 1.0 (incorrect) results in a dilution factor of 0. 307 ''' 308 if radius >= 1.: 309 d = 0.5*(1. - np.sqrt(1. - 1./radius**2)) 310 else: 311 d = 0. 312 return d

313 # 314 # 315 # ------------------------------------------------------------------------------ 316 #

317 -def diCross(diParams, energy=0, verbose=0):

318 ''' 319 Calculate the direct ionization cross section. 320 diParams obtained by util.diRead with the following keys: 321 ['info', 'ysplom', 'xsplom', 'btf', 'ev1', 'ref', 'eaev'] 322 Given as a function of the incident electron energy in eV 323 returns a dictionary - {'energy':energy, 'cross':cross} 324 ''' 325 iso=diParams['info']['iz'] - diParams['info']['ion'] + 1 326 energy = np.array(energy, 'float64') 327 if not energy.any(): 328 btenergy=0.1*np.arange(10) 329 btenergy[0]=0.01 330 dum=np.ones(len(btenergy)) 331 [energy, dum] = descale_bti(btenergy, dum, 2., diParams['ev1'][0]) 332 energy=np.asarray(energy, 'float64') 333 # 334 if iso == 1 and self.Z >= 6: 335 # hydrogenic sequence 336 ryd=27.2113845/2. 337 u=energy/self.Ip 338 ev1ryd=self.Ip/ryd 339 a0=0.5291772108e-8 340 a_bohr=const.pi*a0**2 # area of bohr orbit 341 if self.Z >= 20: 342 ff = (140.+(self.Z/20.)**3.2)/141. 343 else: 344 ff = 1. 345 # qr = util.qrp(self.Z,u)*ff 346 qr = qrp(self.Z,u)*ff 347 bb = 1. # hydrogenic 348 qh = bb*a_bohr*qr/ev1ryd**2 349 diCross = {'energy':energy, 'cross':qh} 350 elif iso == 2 and self.Z >= 10: 351 # use 352 ryd=27.2113845/2. 353 u=energy/self.Ip 354 ev1ryd=self.Ip/ryd 355 a0=0.5291772108e-8 356 a_bohr=const.pi*a0**2 # area of bohr orbit 357 if self.Z >= 20: 358 ff=(140.+(self.Z/20.)**3.2)/141. 359 else: 360 ff=1. 361 # qr=util.qrp(self.Z,u)*ff 362 qr = qrp(self.Z,u)*ff 363 bb = 2. # helium-like 364 qh=bb*a_bohr*qr/ev1ryd**2 365 diCross={'energy':energy, 'cross':qh} 366 else: 367 cross=np.zeros(len(energy), 'Float64') 368 369 for ifac in range(diParams['info']['nfac']): 370 # prob. better to do this with masked arrays 371 goode=energy > diParams['ev1'][ifac] 372 if goode.sum() > 0: 373 dum=np.ones(len(energy)) 374 btenergy, btdum = scale_bti(energy[goode],dum[goode], diParams['btf'][ifac], diParams['ev1'][ifac]) 375 # these interpolations were made with the scipy routine used here 376 y2=interpolate.splrep(diParams['xsplom'][ifac], diParams['ysplom'][ifac], s=0) 377 btcross=interpolate.splev(btenergy, y2, der=0) 378 energy1, cross1 = descale_bti(btenergy, btcross, diParams['btf'][ifac], diParams['ev1'][ifac] ) 379 offset=len(energy)-goode.sum() 380 if verbose: 381 pl.plot(diParams['xsplom'][ifac], diParams['ysplom'][ifac]) 382 pl.plot(btenergy, btcross) 383 if offset > 0: 384 seq=[np.zeros(offset, 'Float64'), cross1] 385 cross1=np.hstack(seq) 386 cross+=cross1*1.e-14 387 return {'energy':energy, 'cross':cross}

388 # 389 # ------------------------------------------------------------------------------ 390 #

391 -def diCross1(diParams, energy=0, verbose=0):

392 ''' 393 Calculate the direct ionization cross section. 394 diParams obtained by util.diRead with the following keys: 395 ['info', 'ysplom', 'xsplom', 'btf', 'ev1', 'ref', 'eaev'] 396 Given as a function of the incident electron energy in eV 397 returns a dictionary - {'energy':energy, 'cross':cross} 398 this version tests whether using the seq and hstack works 399 so using a different approach 400 ''' 401 iso=diParams['info']['iz'] - diParams['info']['ion'] + 1 402 energy = np.array(energy, 'float64') 403 if not energy.any(): 404 btenergy=0.1*np.arange(10) 405 btenergy[0]=0.01 406 dum=np.ones(len(btenergy)) 407 [energy, dum] = descale_bti(btenergy, dum, 2., diParams['ev1'][0]) 408 energy=np.asarray(energy, 'float64') 409 # 410 if iso == 1 and self.Z >= 6: 411 # hydrogenic sequence 412 ryd=27.2113845/2. 413 u=energy/self.Ip 414 ev1ryd=self.Ip/ryd 415 a0=0.5291772108e-8 416 a_bohr=const.pi*a0**2 # area of bohr orbit 417 if self.Z >= 20: 418 ff = (140.+(self.Z/20.)**3.2)/141. 419 else: 420 ff = 1. 421 # qr = util.qrp(self.Z,u)*ff 422 qr = qrp(self.Z,u)*ff 423 bb = 1. # hydrogenic 424 qh = bb*a_bohr*qr/ev1ryd**2 425 diCross = {'energy':energy, 'cross':qh} 426 elif iso == 2 and self.Z >= 10: 427 # use 428 ryd=27.2113845/2. 429 u=energy/self.Ip 430 ev1ryd=self.Ip/ryd 431 a0=0.5291772108e-8 432 a_bohr=const.pi*a0**2 # area of bohr orbit 433 if self.Z >= 20: 434 ff=(140.+(self.Z/20.)**3.2)/141. 435 else: 436 ff=1. 437 # qr=util.qrp(self.Z,u)*ff 438 qr = qrp(self.Z,u)*ff 439 bb=2. # helium-like 440 qh=bb*a_bohr*qr/ev1ryd**2 441 diCross={'energy':energy, 'cross':qh} 442 else: 443 cross=np.zeros(len(energy), 'Float64') 444 445 for ifac in range(diParams['info']['nfac']): 446 # prob. better to do this with masked arrays 447 goode=energy > diParams['ev1'][ifac] 448 if goode.sum() > 0: 449 dum=np.ones(len(energy)) 450 btenergy, btdum = scale_bti(energy[goode],dum[goode], diParams['btf'][ifac], diParams['ev1'][ifac]) 451 # these interpolations were made with the scipy routine used here 452 y2=interpolate.splrep(diParams['xsplom'][ifac], diParams['ysplom'][ifac], s=0) 453 btcross=interpolate.splev(btenergy, y2, der=0) 454 energy1, cross1 = descale_bti(btenergy, btcross, diParams['btf'][ifac], diParams['ev1'][ifac] ) 455 offset=len(energy)-goode.sum() 456 if verbose: 457 pl.plot(diParams['xsplom'][ifac], diParams['ysplom'][ifac]) 458 pl.plot(btenergy, btcross) 459 # if offset > 0: 460 # seq=[np.zeros(offset, 'Float64'), cross1] 461 # cross1=np.hstack(seq) 462 cross[offset:]+=cross1*1.e-14 463 return {'energy':energy, 'cross':cross}

464 # 465 # ------------------------------------------------------------------------------------- 466 #

467 -def eaCross(diparams, easplom, elvlc, energy=None, verbose=False):

468 ''' 469 Provide the excitation-autoionization cross section. 470 471 Energy is given in eV. 472 ''' 473 energy = np.asarray(energy, 'float64') 474 if not energy.any(): 475 btenergy=0.1*np.arange(10) 476 btenergy[0]=0.01 477 dum=np.ones(len(btenergy)) 478 [energy, dum] = descale_bti(btenergy, dum, 2., min(easplom['deryd'])) 479 energy=np.asarray(energy, 'float64') 480 # 481 omega = splomDescale(easplom, energy) 482 # 483 # need to replicate neaev 484 485 if diparams['info']['neaev'] > 0: 486 f1 = np.ones(omega.shape[0]) 487 else: 488 f1 = diparams['info']['eaev'] 489 490 totalCross = np.zeros_like(energy) 491 ntrans = omega.shape[0] 492 for itrans in range(ntrans): 493 lvl1 = easplom['lvl1'][itrans] 494 mult = 2.*elvlc['j'][lvl1 - 1] + 1. 495 cross = f1[itrans]*const.bohrCross*omega[itrans]/(mult.energy/const.ryd2Ev) 496 totalCross += cross 497 return {'energy':energy, 'cross':totalCross}

498 # 499 #----------------------------------------------------------- 500 #

501 -def eaDescale(easplups, temperature):

502 """ 503 Calculates the effective collision strengths (upsilon) 504 for excitation-autoionization as a function of temperature. 505 a duplicate of ion.eaDescale() 506 """ 507 # 508 # xt=kt/de 509 # 510 # need to make sure elvl is >0, except for ground level 511 # 512 ntemp=temperature.size 513 nsplups=len(easplups['de']) 514 if ntemp > 1: 515 ups=np.zeros((nsplups,ntemp),"Float64") 516 else: 517 ups=np.zeros(nsplups,"Float64") 518 # 519 for isplups in range(0,nsplups): 520 l1=easplups["lvl1"][isplups]-1 521 l2=easplups["lvl2"][isplups]-1 522 ttype=easplups["ttype"][isplups] 523 cups=easplups["cups"][isplups] 524 nspl=easplups["nspl"][isplups] 525 de=easplups["de"][isplups] 526 dx=1./(float(nspl)-1.) 527 ## print easplups["easplups"][l1,l2] 528 splups=easplups["splups"][isplups,0:nspl] 529 kte=const.boltzmannEv*temperature/(const.ryd2Ev*de) 530 # 531 if ttype ==1: 532 st=1.-np.log(cups)/np.log(kte+cups) 533 xs=dx*np.arange(nspl) 534 y2=interpolate.splrep(xs,splups,s=0) #allow smoothing,s=0) 535 sups=interpolate.splev(st,y2,der=0) 536 ups[isplups]=sups*np.log(kte+np.exp(1.)) 537 # 538 if ttype == 2: 539 st=kte/(kte+cups) 540 xs=dx*np.arange(nspl) 541 y2=interpolate.splrep(xs,splups,s=0) 542 sups=interpolate.splev(st,y2,der=0) 543 ups[isplups]=sups 544 # 545 if ttype == 3: 546 st=kte/(kte+cups) 547 xs=dx*np.arange(nspl) 548 y2=interpolate.splrep(xs,splups,s=0) 549 sups=interpolate.splev(st,y2,der=0) 550 ups[isplups]=sups/(kte+1.) 551 # 552 if ttype == 4: 553 st=1.-np.log(cups)/np.log(kte+cups) 554 xs=dx*np.arange(nspl) 555 y2=interpolate.splrep(xs,splups,s=0) 556 sups=interpolate.splev(st,y2,der=0) 557 ups[isplups]=sups*np.log(kte+cups) 558 # 559 if ttype == 5: 560 st=kte/(kte+cups) 561 xs=dx*np.arange(nspl) 562 y2=interpolate.splrep(xs,splups,s=0) #allow smoothing,s=0) 563 sups=interpolate.splev(st,y2,der=0) 564 ups[isplups]=sups/(kte+0.) 565 # 566 # 567 elif ttype > 5: print((' t_type ne 1,2,3,4,5 = %5i %5i %5i'%(ttype,l1,l2))) 568 # 569 # 570 ups=np.where(ups > 0.,ups,0.) 571 # 572 # easplups['ups'] = ups 573 return ups

574 # 575 # ------------------------------------------------------------------------------------- 576 #

577 -def listFiles(path):

578 ''' 579 walks the path and subdirectories to return a list of files 580 ''' 581 alist=os.walk(path) 582 # print(' getting file list') 583 listname=[] 584 for (dirpath,dirnames,filenames) in alist: 585 if len(dirnames) == 0: 586 for f in filenames: 587 file=os.path.join(dirpath,f) 588 if os.path.isfile(file): 589 listname.append(file) 590 else: 591 for f in filenames: 592 file=os.path.join(dirpath,f) 593 if os.path.isfile(file): 594 listname.append(file) 595 return listname

596 # 597 #----------------------------------------------------------- 598 #

599 -def scale_bti(evin,crossin,f,ev1):

600 """ 601 apply BT ionization scaling to (energy, cross-section) 602 returns [bte,btx] 603 """ 604 u=evin/ev1 605 bte=1.-np.log(f)/np.log(u-1.+f) 606 btx=u*crossin*(ev1**2)/(np.log(u)+1.) 607 return [bte,btx]

608 # 609 #----------------------------------------------------------- 610 #

611 -def descale_bti(bte,btx,f,ev1):

612 """ 613 descale BT ionization scaling 614 returns [energy,cross-section] 615 """ 616 u=1.-f+np.exp(np.log(f)/(1.-bte)) 617 energy=u*ev1 618 cross=(np.log(u)+1.)*btx/(u*ev1**2) 619 return [energy,cross]

620 # 621 #----------------------------------------------------------- 622 #

623 -def descale_bt(bte,btomega,f,ev1):

624 """ 625 descale BT excitation scaling 626 returns [energy,collision strength] 627 """ 628 u=1.-f+np.exp(np.log(f)/(1.-bte)) 629 energy=u*ev1 630 omega=(np.log(u)-1.+np.exp(1.))*btomega 631 return [energy,omega]

632 # 633 #----------------------------------------------------------- 634 #

635 -def scale_bt(evin,omega,f,ev1):

636 """ 637 apply BT excitation scaling to (energy, collision strength) 638 returns [bte,btomega] 639 """ 640 u=evin/ev1 641 bte=1.-np.log(f)/np.log(u-1.+f) 642 btomega=omega/(np.log(u)-1.+np.exp(1.)) 643 return [bte,btomega]

644

Source Code for Module chianti.util