chianti.core.Continuum

14 ''' 15 The top level class for continuum calculations. 16 17 includes methods for the calculation of the free-free and free-bound continua. 18 '''

19 - def __init__(self, ionStr, temperature, abundance=0, abundanceName=0, em=0, verbose=0):

20 nameDict = util.convertName(ionStr) 21 self.Z = nameDict['Z'] 22 self.Ion = nameDict['Ion'] 23 self.IonStr = ionStr 24 self.Dielectronic = 0 25 self.Defaults = chdata.Defaults 26 self.AbundanceName = self.Defaults['abundfile'] 27 self.IoneqName = self.Defaults['ioneqfile'] 28 # 29 # ip in eV, reading Ip of next lower level, needed for freeBound 30 31 if self.Ion == 1: 32 if verbose: 33 print(' for %s this is the neutral ions an does not produce a continuum'%(ionStr)) 34 return 35 # the Ip is only relevant to the free-free methods 36 # self.Ip = ip[self.Z-1, self.Ion-1] 37 self.Ipr = ip[self.Z-1, self.Ion-2] 38 # 39 if type(temperature) == float and temperature > 0.: 40 self.Temperature = np.asarray(temperature,'float64') 41 elif type(temperature) == list or type(temperature) == tuple or type(temperature) == np.ndarray: 42 temperature = np.asarray(temperature, 'float64') 43 self.Temperature = temperature 44 # 45 if abundance: 46 self.Abundance = abundance 47 elif abundanceName: 48 if abundanceName in sorted(chdata.Abundance.keys()): 49 self.AbundanceName = abundanceName 50 self.Abundance = chdata.Abundance[self.AbundanceName]['abundance'][self.Z-1] 51 else: 52 abundChoices = sorted(chdata.Abundance.keys()) 53 # for one in wvl[topLines]: 54 # wvlChoices.append('%12.3f'%(one)) 55 abundChoice = chgui.gui.selectorDialog(abundChoices,label='Select Abundance name') 56 abundChoice_idx = abundChoice.selectedIndex 57 self.AbundanceName = abundChoices[abundChoice_idx[0]] 58 abundanceName = self.AbundanceName 59 print(' Abundance chosen: %s '%(self.AbundanceName)) 60 self.Abundance = chdata.Abundance[self.AbundanceName]['abundance'][self.Z-1] 61 else: 62 self.AbundanceName = self.Defaults['abundfile'] 63 self.Abundance = chdata.Abundance[self.AbundanceName]['abundance'][self.Z-1] 64 # 65 self.ioneqOne() 66 # 67 if type(em) == int and em == 0: 68 em = np.ones_like(self.Temperature, 'float64') 69 self.Em = em 70 elif type(em) == float and em > 0.: 71 em = np.ones_like(self.Temperature, 'float64')*em 72 self.Em = em 73 elif type(em) == list or type(em) == tuple or type(em) == np.ndarray: 74 em = np.asarray(em, 'float64') 75 self.Em = em

76 # 77 #---------------------------------------------------------------------------------------- 78 #

79 - def freeBoundEmiss(self, wvl, verner=1):

80 '''Calculates the free-bound (radiative recombination) continuum emissivity of an ion. 81 82 Uses the Gaunt factors of Karzas, W.J, Latter, R, 1961, ApJS, 6, 167 83 for recombination to the ground level, the photoionization cross sections of 84 Verner and Yakovlev, 1995, A&ASS, 109, 125 85 are used to develop the free-bound cross section 86 provides emissivity = ergs cm^-2 s^-1 str^-1 Angstrom ^-1 for an individual ion 87 does not include the elemental abundance or ionization fraction 88 the specified ion is the target ion''' 89 # 90 wvl = np.asarray(wvl, 'float64') 91 temperature = self.Temperature 92 # 93 # data for the recombined ion 94 # 95 if hasattr(self, 'Fblvl'): 96 fblvl = self.Fblvl 97 else: 98 fblvlname = util.zion2filename(self.Z,self.Ion-1)+'.fblvl' 99 self.Fblvl = chio.fblvlRead(fblvlname) 100 fblvl = self.Fblvl 101 if 'errorMessage' in fblvl.keys(): 102 self.FreeBound = fblvl 103 return 104 # 105 # need some data for the recombining ion 106 # 107 if hasattr(self, 'rFblvl'): 108 rFblvl = self.rFblvl 109 else: 110 if self.Ion == self.Z+1: 111 # then we looking for the bare ion 112 rFblvl = {'mult':[1., 1.]} 113 else: 114 rfblvlname = util.zion2filename(self.Z,self.Ion)+'.fblvl' # previously self.Ion) 115 self.rFblvl = chio.fblvlRead(rfblvlname) 116 rFblvl = self.rFblvl 117 if 'errorMessage' in rFblvl.keys(): 118 self.FreeBound = rFblvl 119 return 120 # 121 # 6/9/2010 the recombining ion is the present ion 122 # 123 nlvls = len(fblvl['lvl']) 124 # pqn = principle quantum no. n 125 pqn = fblvl['pqn'] 126 # l is angular moment quantum no. L 127 l = fblvl['l'] 128 # energy level in inverse cm 129 ecm = fblvl['ecm'] 130 # statistical weigths/multiplicities 131 mult = fblvl['mult'] 132 multr = rFblvl['mult'] 133 # 134 # Ipr 135 # 136 # for the ionization potential, must use that of the recombined ion 137 # Ipr = ip[self.Z-1, self.Ion-2] 138 # ipcm = self.Ip/const.invCm2Ev 139 iprcm = self.Ipr/const.invCm2Ev 140 #for i in range(nlvls): 141 #print(' lvl %5i ecm %12.3e Ip-wvl %12.4f Ipr-wvl %12.4f '%(i, ecm[i], 1.e+8/(ipcm-ecm[i]), 1.e+8/(iprcm-ecm[i]))) 142 # 143 # get karzas-latter Gaunt factors 144 if hasattr(self,'Klgfb'): 145 klgfb = self.Klgfb 146 else: 147 self.Klgfb = chio.klgfbRead() 148 klgfb = self.Klgfb 149 # 150 # 151 nWvl = wvl.size 152 nTemp = temperature.size 153 # 154 if verner: 155 lvl1 = 1 156 else: 157 lvl1 = 0 158 # 159 if (nTemp > 1) and (nWvl > 1): 160 mask = np.zeros((nlvls,nTemp,nWvl),'Bool') 161 fbrate = np.zeros((nlvls,nTemp,nWvl),'float64') 162 expf = np.zeros((nlvls,nTemp,nWvl),'float64') 163 ratg = np.zeros((nlvls),'float64') 164 if verner: 165 self.vernerCross(wvl) 166 vCross = self.VernerCross 167 # 168 ratg[0] = float(mult[0])/float(multr[0]) 169 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[0] 170 iprLvlErg = const.ev2Erg*iprLvlEv 171 iprLvlCm = (iprcm - ecm[0]) 172 for itemp in range(nTemp): 173 mask[0,itemp] = 1.e+8/wvl < (iprcm - ecm[0]) 174 expf[0,itemp] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature[itemp])) 175 fbrate[0,itemp] = (const.planck*const.light/(1.e-8*wvl))**5*const.verner*ratg[0]*expf[0,itemp]*vCross/temperature[itemp]**1.5 176 for ilvl in range(lvl1,nlvls): 177 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[ilvl] 178 iprLvlErg = const.ev2Erg*iprLvlEv 179 scaledE = np.log(const.ev2Ang/(iprLvlEv*wvl)) 180 thisGf = klgfb['klgfb'][pqn[ilvl]-1, l[ilvl]] 181 spl = interpolate.splrep(klgfb['pe'], thisGf) 182 gf = np.exp(interpolate.splev(scaledE, spl)) 183 ratg[ilvl] = float(mult[ilvl])/float(multr[0]) # ratio of statistical weights 184 # 185 for itemp in range(nTemp): 186 expf[ilvl,itemp] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature[itemp])) 187 mask[ilvl,itemp] = 1.e+8/wvl < (iprcm - ecm[ilvl]) 188 fbrate[ilvl,itemp] = const.freeBound*ratg[ilvl]*(iprLvlErg**2/float(pqn[ilvl]))*gf*expf[ilvl,itemp]/(temperature[itemp]**1.5*(wvl)**2) 189 fbrma = np.ma.array(fbrate) 190 fbrma.mask = mask 191 fbrma.fill_value = 0. 192 fbRate = (fbrma).sum(axis=0) 193 fbRate.fill_value = 0. 194 self.FreeBoundEmiss = {'emiss':fbRate.data, 'temperature':temperature,'wvl':wvl} 195 # 196 elif (nTemp == 1) and (nWvl > 1): 197 mask = np.zeros((nlvls,nWvl),'Bool') 198 fbrate = np.zeros((nlvls,nWvl),'float64') 199 expf = np.zeros((nlvls,nWvl),'float64') 200 #ratg = np.zeros((nlvls),'float64') 201 if verner: 202 self.vernerCross(wvl) 203 vCross = self.VernerCross 204 # 205 ratg = float(mult[0])/float(multr[0]) 206 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[0] 207 iprLvlErg = const.ev2Erg*iprLvlEv 208 iprLvlCm = (iprcm - ecm[0]) 209 # 210 mask[0] = 1.e+8/wvl < (iprcm - ecm[0]) 211 expf[0] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 212 fbrate[0] = (const.planck*const.light/(1.e-8*wvl))**5*const.verner*ratg*expf[0]*vCross/temperature**1.5 213 for ilvl in range(lvl1,nlvls): 214 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[ilvl] 215 iprLvlErg = const.ev2Erg*iprLvlEv 216 iprLvlCm = (iprcm - ecm[ilvl]) 217 mask[ilvl] = 1.e+8/wvl < iprLvlCm 218 # # print ' ilvl, ipLvlEv = ', ilvl, ipLvlEv 219 scaledE = np.log(const.ev2Ang/(iprLvlEv*wvl)) 220 thisGf = klgfb['klgfb'][pqn[ilvl]-1, l[ilvl]] 221 spl = interpolate.splrep(klgfb['pe'], thisGf) 222 gf = np.exp(interpolate.splev(scaledE, spl)) 223 ratg = float(mult[ilvl])/float(multr[0]) # ratio of statistical weights 224 expf[ilvl] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 225 fbrate[ilvl] = const.freeBound*ratg*expf[ilvl]*(iprLvlErg**2/float(pqn[ilvl]))*gf/(temperature**1.5*(wvl)**2) 226 fbrma = np.ma.array(fbrate) 227 fbrma.mask = mask 228 fbrma.fill_value = 0. 229 fbRate = fbrma.sum(axis=0) 230 fbRate.fill_value = 0. 231 # factor of 1.e-8 converts to Angstrom^-1, otherwise it would be cm^-1 232 # fbRate = (expf*fbrma).sum(axis=0) 233 #fbRma = np.ma.array(fbrate.sum(axis=0), mask=mask, fill_value=0.) 234 # fbRate.fill_value = 0. 235 self.FreeBoundEmiss = {'fbRate':fbRate, 'emiss':fbRate.data, 'temperature':temperature,'wvl':wvl, 'mask':mask, 'expf':expf,'vCross':vCross} 236 else: 237 # this is the case of a single wavelength, not going to do this anymore 238 #mask = np.zeros((nlvls,nTemp),'Bool') 239 #fbrate = np.zeros((nlvls,nTemp),'float64') 240 #expf = np.zeros((nlvls,nTemp),'float64') 241 #ratg = np.zeros((nlvls),'float64') 242 #if verner: 243 #self.vernerCross(wvl) 244 #vCross = self.VernerCross 245 ## 246 #mask[0] = 1.e+8/wvl < (ipcm - ecm[0]) 247 #ratg[0] = float(mult[0])/float(multr[0]) 248 #ipLvlEv = self.Ip - const.invCm2Ev*ecm[0] 249 #ipLvlErg = const.ev2Erg*ipLvlEv 250 #expf[0] = np.exp((ipLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 251 #fbrate[0] = (const.planck*const.light/(1.e-8*wvl))**5*const.verner*ratg[0]*expf[0]*vCross/temperature**1.5 252 #for ilvl in range(lvl1,nlvls): 253 ## scaled energy is relative to the ionization potential of each individual level 254 #ipLvlEv = self.Ip - const.invCm2Ev*ecm[ilvl] 255 #scaledE = np.log(const.ev2Ang/(ipLvlEv*wvl)) 256 #thisGf = klgfb['klgfb'][pqn[ilvl]-1, l[ilvl]] 257 #spl = interpolate.splrep(klgfb['pe'], thisGf) 258 #gf = np.exp(interpolate.splev(scaledE, spl)) 259 #mask[ilvl] = 1.e+8/wvl < (ipcm - ecm[ilvl]) 260 #ratg[ilvl] = float(mult[ilvl])/float(multr[0]) # ratio of statistical weights 261 #ipLvlErg = const.ev2Erg*ipLvlEv 262 #expf[ilvl] = np.exp((ipLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 263 #fbrate[ilvl] = const.freeBound*ratg[ilvl]*(ipLvlErg**2/float(pqn[ilvl]))*expf[ilvl]*gf/(temperature**1.5*(wvl)**2) 264 #fbrma = np.ma.array(fbrate) 265 #fbrma.mask = mask 266 #fbrma.fill_value = 0. 267 ## factor of 1.e-8 converts to Angstrom^-1, otherwise it would be cm^-1 268 #fbRate = (fbrma).sum(axis=0) 269 #fbRate.fill_value = 0. 270 #self.FreeBoundEmiss = {'emiss':fbRate.data, 'temperature':temperature,'wvl':wvl} 271 self.FreeBoundEmiss = {'errorMessage':' this is the case of a single wavelength'}

272 # 273 # ---------------------------------------------------------------------------- 274 #

275 - def freeBound(self, wvl, verner=1):

276 ''' 277 to calculate the free-bound (radiative recombination) continuum rate coefficient of an ion, where 278 the ion is taken to be the recombined ion, 279 including the elemental abundance and the ionization equilibrium population 280 uses the Gaunt factors of Karzas, W.J, Latter, R, 1961, ApJS, 6, 167 281 for recombination to the ground level, the photoionization cross sections of 282 Verner and Yakovlev, 1995, A&ASS, 109, 125 283 are used to develop the free-bound cross section 284 includes the elemental abundance and the ionization fraction 285 provides emissivity = ergs cm^-2 s^-1 str^-1 Angstrom ^-1 286 ''' 287 wvl = np.asarray(wvl, 'float64') 288 temperature = self.Temperature 289 # 290 if hasattr(self, 'IoneqOne'): 291 gIoneq = self.IoneqOne 292 else: 293 self.ioneqOne() 294 gIoneq = self.IoneqOne 295 # 296 # put in freefree to go through ipymspectrum 297 if not np.any(gIoneq) > 0: 298 self.FreeBound = {'errorMessage':' no non-zero values of ioneq'} 299 return 300 # 301 em = self.Em 302 # 303 # the recombined ion contains that data for fblvl 304 # 305 if hasattr(self,'Fblvl'): 306 fblvl = self.Fblvl 307 else: 308 fblvlname = util.zion2filename(self.Z,self.Ion-1)+'.fblvl' 309 self.Fblvl = chio.fblvlRead(fblvlname) 310 fblvl = self.Fblvl 311 # in case there is no fblvl file 312 if 'errorMessage' in fblvl.keys(): 313 self.FreeBound = fblvl 314 return 315 # 316 # need data for the current/recombining ion 317 # 318 if hasattr(self,'rFblvl'): 319 rFblvl = self.rFblvl 320 else: 321 if self.Ion == self.Z + 1: 322 # this is a bare ion 323 rFblvl = {'mult':[1., 1.]} 324 else: 325 rfblvlname = util.zion2filename(self.Z,self.Ion)+'.fblvl' 326 self.rFblvl = chio.fblvlRead(rfblvlname) 327 rFblvl = self.rFblvl 328 if 'errorMessage' in rFblvl.keys(): 329 self.FreeBound = rFblvl 330 return 331 # 332 #gIoneq = self.IoneqOne 333 # 334 # 335 abund = self.Abundance 336 # 337 # 338 nlvls = len(fblvl['lvl']) 339 # pqn = principle quantum no. n 340 pqn = np.asarray(fblvl['pqn'], 'float64') 341 # l is angular moment quantum no. L 342 l = fblvl['l'] 343 # energy level in inverse cm 344 ecm = fblvl['ecm'] 345 # statistical weigths/multiplicities 346 mult = fblvl['mult'] 347 multr = rFblvl['mult'] 348 # ipcm = self.Ip/const.invCm2Ev 349 # 350 # 351 # for the ionization potential, must use that of the recombined ion 352 # 353 iprcm = self.Ipr/const.invCm2Ev 354 # 355 # get karzas-latter Gaunt factors 356 if hasattr(self,'Klgfb'): 357 klgfb = self.Klgfb 358 else: 359 self.Klgfb = chio.klgfbRead() 360 klgfb = self.Klgfb 361 # 362 nWvl = wvl.size 363 nTemp = temperature.size 364 # 365 if verner: 366 lvl1 = 1 367 else: 368 lvl1 = 0 369 # 370 nWvl = wvl.size 371 nTemp = temperature.size 372 # 373 # 374 if (nTemp > 1) and (nWvl > 1): 375 mask = np.zeros((nlvls,nTemp,nWvl),'Bool') 376 fbrate = np.zeros((nlvls,nTemp,nWvl),'float64') 377 expf = np.zeros((nlvls,nTemp,nWvl),'float64') 378 ratg = np.zeros((nlvls),'float64') 379 #fbRate = np.zeros((nTemp,nWvl),'float64') 380 if verner: 381 self.vernerCross(wvl) 382 vCross = self.VernerCross 383 ratg[0] = float(mult[0])/float(multr[0]) 384 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[0] 385 iprLvlErg = const.ev2Erg*iprLvlEv 386 iprLvlCm = (iprcm - ecm[0]) 387 for itemp in range(nTemp): 388 mask[0,itemp] = 1.e+8/wvl < (iprcm - ecm[0]) 389 expf[0,itemp] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature[itemp])) 390 fbrate[0,itemp] = em[itemp]*(const.planck*const.light/(1.e-8*wvl))**5*const.verner*gIoneq[itemp]*ratg[0]*expf[0,itemp]*vCross/temperature[itemp]**1.5 391 for ilvl in range(lvl1,nlvls): 392 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[ilvl] 393 iprLvlErg = const.ev2Erg*iprLvlEv 394 scaledE = np.log(const.ev2Ang/(iprLvlEv*wvl)) 395 thisGf = klgfb['klgfb'][pqn[ilvl]-1, l[ilvl]] 396 spl = interpolate.splrep(klgfb['pe'], thisGf) 397 gf = np.exp(interpolate.splev(scaledE, spl)) 398 ratg[ilvl] = float(mult[ilvl])/float(multr[0]) # ratio of statistical weights 399 # 400 for itemp in range(nTemp): 401 expf[ilvl] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature[itemp])) 402 expf[ilvl,itemp] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature[itemp])) 403 mask[ilvl,itemp] = 1.e+8/wvl < (iprcm - ecm[ilvl]) 404 fbrate[ilvl,itemp] = em[itemp]*const.freeBound*gIoneq[itemp]*ratg[ilvl]*(iprLvlErg**2/float(pqn[ilvl]))*gf*expf[ilvl,itemp]/(temperature[itemp]**1.5*(wvl)**2) 405 fbrma = np.ma.array(fbrate) 406 fbrma.mask = mask 407 fbrma.fill_value = 0. 408 fbRate = (fbrma).sum(axis=0) 409 fbRate.fill_value = 0. 410 self.FreeBound = {'rate':abund*fbRate.data, 'temperature':temperature,'wvl':wvl,'em':em} 411 # 412 elif (nTemp == 1) and (nWvl > 1): 413 mask = np.zeros((nlvls,nWvl),'Bool') 414 fbrate = np.zeros((nlvls,nWvl),'float64') 415 expf = np.zeros((nlvls,nWvl),'float64') 416 ratg = np.zeros((nlvls),'float64') 417 if verner: 418 self.vernerCross(wvl) 419 vCross = self.VernerCross 420 # mask is true for bad values 421 ratg[0] = float(mult[0])/float(multr[0]) 422 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[0] 423 iprLvlErg = const.ev2Erg*iprLvlEv 424 iprLvlCm = (iprcm - ecm[0]) 425 # 426 mask[0] = 1.e+8/wvl < (iprcm - ecm[0]) 427 expf[0] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 428 fbrate[0] = em*(const.planck*const.light/(1.e-8*wvl))**5*const.verner*ratg[0]*expf[0]*vCross/temperature**1.5 429 # 430 for ilvl in range(lvl1,nlvls): 431 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[ilvl] 432 iprLvlErg = const.ev2Erg*iprLvlEv 433 iprLvlCm = (iprcm - ecm[ilvl]) 434 # scaled energy is relative to the ionization potential of each individual level 435 scaledE = np.log(const.ev2Ang/(iprLvlEv*wvl)) 436 thisGf = klgfb['klgfb'][pqn[ilvl]-1, l[ilvl]] 437 spl = interpolate.splrep(klgfb['pe'], thisGf) 438 gf = np.exp(interpolate.splev(scaledE, spl)) 439 mask[ilvl] = 1.e+8/wvl < iprLvlCm 440 ratg[ilvl] = float(mult[ilvl])/float(multr[0]) # ratio of statistical weights 441 expf[ilvl] = np.exp((iprLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 442 fbrate[ilvl] = em*const.freeBound*ratg[ilvl]*(iprLvlErg**2/float(pqn[ilvl]))*expf[ilvl]*gf/(temperature**1.5*(wvl)**2) 443 fbrma = np.ma.array(fbrate) 444 fbrma.mask = mask 445 fbrma.fill_value = 0. 446 # factor of 1.e-8 converts to Angstrom^-1, otherwise it would be cm^-1 447 fbRate = abund*gIoneq*fbrma.sum(axis=0) 448 fbRate.fill_value = 0. 449 self.FreeBound = {'rate':fbRate.data, 'temperature':temperature,'wvl':wvl} 450 #elif (nTemp > 1) and (nWvl == 1): 451 else: 452 #mask = np.zeros((nlvls,nTemp),'Bool') 453 #fbrate = np.zeros((nlvls,nTemp),'float64') 454 #expf = np.zeros((nlvls,nTemp),'float64') 455 #ratg = np.zeros((nlvls),'float64') 456 #if verner: 457 #self.vernerCross(wvl) 458 #vCross = self.VernerCross 459 #mask[0] = 1.e+8/wvl < (ipcm - ecm[0]) 460 #ratg[0] = float(mult[0])/float(multr[0]) 461 #ipLvlEv = self.Ip - const.invCm2Ev*ecm[0] 462 #ipLvlErg = const.ev2Erg*ipLvlEv 463 #expf[0] = np.exp((ipLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 464 #fbrate[0] = self.Em[itemp]*(const.planck*const.light/(1.e-8*wvl))**5*const.verner*ratg[0]*expf[0]*vCross/temperature**1.5 465 #for ilvl in range(lvl1,nlvls): 466 ## scaled energy is relative to the ionization potential of each individual level 467 #ipLvlEv = self.Ip - const.invCm2Ev*ecm[ilvl] 468 #scaledE = np.log(const.ev2Ang/(ipLvlEv*wvl)) 469 #thisGf = klgfb['klgfb'][pqn[ilvl]-1, l[ilvl]] 470 #spl = interpolate.splrep(klgfb['pe'], thisGf) 471 #gf = np.exp(interpolate.splev(scaledE, spl)) 472 #mask[ilvl] = 1.e+8/wvl < (ipcm - ecm[ilvl]) 473 #ratg[ilvl] = float(mult[ilvl])/float(multr[0]) # ratio of statistical weights 474 #ipLvlErg = const.ev2Erg*ipLvlEv 475 #expf[ilvl] = np.exp((ipLvlErg - 1.e+8*const.planck*const.light/wvl)/(const.boltzmann*temperature)) 476 #fbrate[ilvl] = self.Em[itemp]*const.freeBound*ratg[ilvl]*(ipLvlErg**2/float(pqn[ilvl]))*expf[ilvl]*gf/(temperature**1.5*(wvl)**2) 477 #fbrma = np.ma.array(fbrate) 478 #fbrma.mask = mask 479 #fbrma.fill_value = 0. 480 ## factor of 1.e-8 converts to Angstrom^-1, otherwise it would be cm^-1 481 #fbRate = abund*gIoneq*(fbrma).sum(axis=0) 482 #fbRate.fill_value = 0. 483 #self.FreeBound = {'rate':fbRate.data, 'temperature':temperature,'wvl':wvl, 'em':self.Em} 484 self.FreeBoundEmiss = {'errorMessage':' this is the case of a single wavelength'}

485 # 486 # ---------------------------------------------------------------------------------------- 487 # 488 #

489 - def freeBoundLoss(self):

490 ''' 491 to calculate the free-bound (radiative recombination) energy loss rate coefficient of an ion, 492 the ion is taken to be the recombined iion, 493 including the elemental abundance and the ionization equilibrium population 494 uses the Gaunt factors of Karzas, W.J, Latter, R, 1961, ApJS, 6, 167 495 provides rate = ergs cm^-2 s^-1 496 ''' 497 # 498 temperature = self.Temperature 499 # 500 if hasattr(self, 'Fblvl'): 501 fblvl = self.Fblvl 502 else: 503 fblvlname = util.zion2filename(self.Z,self.Ion-1)+'.fblvl' 504 if os.path.isfile(fblvlname): 505 self.Fblvl = chio.fblvlRead(fblvlname) 506 fblvl = self.Fblvl 507 else: 508 self.FreeBoundLoss = {'errorMessage':' file does not exist %s .fblvl'%(fblvlname)} 509 return 510 # need some data for the recombining/target ion 511 # 512 if hasattr(self, 'rFblvl'): 513 rFblvl = self.rFblvl 514 else: 515 if self.Ion == self.Z + 1: 516 # then we looking for the bare ion 517 rFblvl = {'mult':[1., 1.]} 518 else: 519 rfblvlname = util.zion2filename(self.Z,self.Ion)+'.fblvl' 520 if os.path.isfile(rfblvlname): 521 self.rFblvl = chio.fblvlRead(rfblvlname) 522 rFblvl = self.rFblvl 523 else: 524 self.FreeBoundLoss = {'errorMessage':' file does not exist %s .fblvl'%(rfblvlname)} 525 return 526 # 527 gIoneq = self.IoneqOne 528 # 529 abund = self.Abundance 530 # 531 # 532 nlvls = len(fblvl['lvl']) 533 # pqn = principle quantum no. n 534 pqn = np.asarray(fblvl['pqn'], 'float64') 535 # l is angular moment quantum no. L 536 l = fblvl['l'] 537 # energy level in inverse cm 538 ecm = fblvl['ecm'] 539 # statistical weigths/multiplicities 540 mult = fblvl['mult'] 541 multr = rFblvl['mult'] 542 # 543 # 544 # for the ionization potential, must use that of the recombined ion 545 # 546 iprcm = self.Ipr/const.invCm2Ev 547 # 548 # get karzas-latter Gaunt factors 549 if hasattr(self, 'Klgfb'): 550 klgfb = self.Klgfb 551 else: 552 self.Klgfb = chio.klgfbRead() 553 klgfb = self.Klgfb 554 # 555 nTemp = temperature.size 556 # statistical weigths/multiplicities 557 # 558 # 559 #wecm=1.e+8/(ipcm-ecm) 560 # 561 # sometime the rFblvl file does not exist 562 if 'mult' in fblvl.keys() and 'mult' in rFblvl.keys(): 563 # 564 nlvls = len(fblvl['lvl']) 565 # pqn = principle quantum no. n 566 pqn = fblvl['pqn'] 567 # l is angular moment quantum no. L 568 l = fblvl['l'] 569 # energy level in inverse cm 570 ecm = fblvl['ecm'] 571 mult = fblvl['mult'] 572 multr = rFblvl['mult'] 573 fbrate = np.zeros((nlvls,nTemp),'float64') 574 ratg = np.zeros((nlvls),'float64') 575 for ilvl in range(nlvls): 576 # scaled energy is relative to the ionization potential of each individual level 577 # will add the average energy of a free electron to this to get typical photon energy to 578 # evaluate the gaunt factor 579 hnuEv = 1.5*const.boltzmann*temperature/const.ev2Erg 580 iprLvlEv = self.Ipr - const.invCm2Ev*ecm[ilvl] 581 scaledE = np.log(hnuEv/iprLvlEv) 582 thisGf = klgfb['klgfb'][pqn[ilvl]-1, l[ilvl]] 583 spl = interpolate.splrep(klgfb['pe'], thisGf) 584 gf = np.exp(interpolate.splev(scaledE, spl)) 585 ratg[ilvl] = float(mult[ilvl])/float(multr[0]) # ratio of statistical weights 586 iprLvlErg = const.ev2Erg*iprLvlEv 587 fbrate[ilvl] = ratg[ilvl]*(iprLvlErg**2/float(pqn[ilvl]))*gf/np.sqrt(temperature) 588 fbRate = abund*gIoneq*const.freeBoundLoss*(fbrate.sum(axis=0)) 589 else: 590 fbRate = np.zeros((nTemp),'float64') 591 self.FreeBoundLoss = {'rate':fbRate, 'temperature':temperature}

592 # 593 # ---------------------------------------------------------------------------------------- 594 #

595 - def vernerCross(self,wvl):

596 '''Calculates the photoionization cross section. 597 598 The data are from Verner and Yakovlev 599 the cross section refers to the next lower ionization stage''' 600 try: 601 vernerDat = self.VernerDat 602 except: 603 self.VernerDat = chio.vernerRead() 604 vernerDat = self.VernerDat 605 z = self.Z 606 stage = self.Ion 607 # ip = self.Ip 608 iprcm = self.Ipr/const.invCm2Ev 609 ecm = self.Fblvl['ecm'] 610 # 611 en = const.ev2Ang/wvl 612 y = en/vernerDat['e0'][z,stage-1] 613 fy= vernerDat['sig0'][z,stage-1]*((y - 1.)**2 + vernerDat['yw'][z,stage-1]**2) * y**(-5.5 - vernerDat['l'][z,stage-1] + 0.5*vernerDat['p'][z,stage-1]) * (1. + np.sqrt(y/vernerDat['ya'][z,stage-1]))**(-vernerDat['p'][z,stage-1]) 614 # mask = en < vernerDat['eth'][z,stage] 615 # will use Chianti values for energy of ground level 616 mask = (1.e+8/wvl) < (iprcm - ecm[0]) 617 vCross = np.ma.array(fy) 618 vCross.mask = mask 619 vCross.fill_value = 0. 620 # cross-section will be output in cm^2 621 self.VernerCross = vCross*1.e-18

622 # 623 # ---------------------------------------------------------------------------------------- 624 #

625 - def freeFree(self, wvl):

626 ''' 627 Calculates the free-free emission for a single ion. 628 629 Includes elemental abundance and ionization equilibrium population 630 and the emission measure (em) if specified 631 Uses Itoh where valid and Sutherland elsewhere 632 ''' 633 em=self.Em 634 temperature = self.Temperature 635 nTemp = temperature.size 636 # 637 if self.Ion == 1: 638 self.FreeFree = {'errorMessage':' freefree is not produced by neutrals'} 639 return 640 else: 641 wvl = np.asarray(wvl, 'float64') 642 ffs = self.sutherland(wvl) 643 # 644 ffi = self.itoh(wvl) 645 ff = ffs['suthFf'] 646 if not 'errorMessage' in ffi.keys(): 647 iff = ffi['itohFf'] 648 itohMask = np.logical_not(iff.mask) 649 ff[itohMask] = iff[itohMask] 650 ff = np.asarray(ff,'float64') 651 # 652 if hasattr(self, 'IoneqOne'): 653 gIoneq = self.IoneqOne 654 else: 655 self.ioneqOne() 656 gIoneq = self.IoneqOne 657 # 658 # put in freefree to go through ipymspectrum 659 if not np.any(gIoneq) > 0: 660 self.FreeFree = {'errorMessage':' no non-zero values of ioneq'} 661 return 662 # 663 if type(gIoneq) == float: 664 # only one temperature specified 665 if gIoneq == 0.: 666 ffRate = np.zeros(wvl.size) 667 self.FreeFree = {'rate':ffRate, 'temperature':self.Temperature,'wvl':wvl} 668 return 669 else: 670 if gIoneq.sum() == 0.: 671 ffRate = np.zeros((self.Temperature.size, wvl.size), 'float64') 672 self.FreeFree = {'rate':ffRate, 'temperature':self.Temperature,'wvl':wvl} 673 return 674 # 675 #for one in gIoneq.shape: 676 #print('gIoneq shape = %10i'%(one)) 677 #for one in ff.shape: 678 #print('ff shape = %10i'%(one)) 679 # 680 if wvl.size > 1: 681 #gIoneq = gIoneq.repeat(wvl.size).reshape(self.Temperature.size,wvl.size) 682 #em = em.repeat(wvl.size).reshape(self.Temperature.size,wvl.size) 683 # 684 try: 685 abund = self.Abundance 686 except: 687 self.AbundanceAll = util.abundanceRead(abundancename = self.AbundanceName) 688 self.Abundance = self.AbundanceAll['abundance'][self.Z-1] 689 abund = self.Abundance 690 # 691 if nTemp == 1: 692 ffRate = em*const.freeFree*(self.Z)**2*abund*gIoneq*ff 693 elif nTemp > 1: 694 ffRate = np.zeros((nTemp,wvl.size),'float64') 695 for it in range(nTemp): 696 ffRate[it] = (em[it]*const.freeFree*(self.Z)**2*abund*gIoneq[it]*ff[it]) # .squeeze() 697 self.FreeFree = {'rate':ffRate, 'temperature':self.Temperature,'wvl':wvl,'ff':ff}

698 # 699 # ---------------------------------------------------------------------------------------- 700 #

701 - def freeFreeEmiss(self, wvl):

702 '''Calculates the free-free emissivity for a single ion. 703 does not include element abundance or ionization fraction 704 Uses Itoh where valid and Sutherland elsewhere''' 705 if self.Ion == 1: 706 self.FreeFreeEmiss = {'errorMessage':' freefree is not produced by neutrals'} 707 else: 708 wvl = np.asarray(wvl, 'float64') 709 ffs = self.sutherland(wvl) 710 ff = ffs['suthFf'] 711 ffi = self.itoh(wvl) 712 if 'errorMessage' not in ffi.keys(): 713 iff = ffi['itohFf'] 714 itohMask = np.logical_not(iff.mask) 715 ff[itohMask] = iff[itohMask] 716 # 717 ffRate = (const.freeFree*(self.Z)**2*ff).squeeze() 718 self.FreeFree = {'rate':ffRate, 'temperature':self.Temperature,'wvl':wvl}

719 # 720 # ---------------------------------------------------------------------------------------- 721 #

722 - def freeFreeLoss(self):

723 ''' 724 Calculates the total free-free emission for a single ion. 725 726 Includes elemental abundance and ionization equilibrium population. 727 728 Uses Itoh where valid and Sutherland elsewhere 729 ''' 730 temperature = self.Temperature 731 #print(' doing freeFreeLoss for %s'%(self.IonStr)) 732 # 733 if self.Ion == 1: 734 self.FreeFree = {'errorMessage':' freefree is not produced by neutrals'} 735 else: 736 # 737 # should be done with masked arrays 738 # 739 #if hasattr(self, 'Gffint'): 740 #gffint = self.Gffint['gffint'] 741 #g2 = self.Gffint['g2'] 742 #else: 743 #self.Gffint = chio.gffintRead() 744 #gffint = self.Gffint['gffint'] 745 #g2 = self.Gffint['g2'] 746 ## 747 #gamma2 = self.Ip*const.ev2Erg/(const.boltzmann*temperature) 748 ## 749 #try: 750 #spl = interpolate.splrep(g2, gffint) 751 #gff = interpolate.splev(np.log(gamma2), spl) 752 #except: 753 #print(' for %s problem with interpolate.splev(np.log(gamma2), spl)'%(self.IonStr)) 754 #self.FreeFreeLoss = {'rate':np.zeros_like(temperature),'temperature':temperature,'errorMessage':'problems with interpolate'} 755 #return 756 # 757 # this is reallly close enough 758 # 759 gff = 1.2 760 # 761 if hasattr(self, 'IoneqOne'): 762 gIoneq = self.IoneqOne 763 else: 764 self.ioneqOne() 765 gIoneq = self.IoneqOne 766 # 767 if hasattr(self, 'Abundance'): 768 abund = self.Abundance 769 else: 770 self.AbundanceAll = util.abundanceRead(abundancename = self.AbundanceName) 771 self.Abundance = self.AbundanceAll['abundance'][self.Z-1] 772 abund = self.Abundance 773 # 774 #print('gff min %12.e max %12.e'%(gff.min(),gff.max())) 775 #if np.isnan(gff.any()): 776 #ffRate = np.zeros_like(temperature) 777 #else: 778 ffRate = const.freeFreeLoss*(self.Z)**2*abund*gIoneq*gff*np.sqrt(temperature) 779 self.FreeFreeLoss = {'rate':ffRate, 'temperature':temperature}

780 # 781 # ---------------------------------------------------------------------------------------- 782 #

783 - def klgfbInterp(self, wvl, n, l):

784 '''A Python version of the CHIANTI IDL procedure karzas_xs. 785 786 Interpolates free-bound gaunt factor of Karzas and Latter, (1961, Astrophysical Journal 787 Supplement Series, 6, 167) as a function of wavelength (wvl).''' 788 try: 789 klgfb = self.Klgfb 790 except: 791 self.Klgfb = util.klgfbRead() 792 klgfb = self.Klgfb 793 # get log of photon energy relative to the ionization potential 794 sclE = np.log(self.Ip/(wvl*const.ev2ang)) 795 thisGf = klgfb['klgfb'][n-1, l] 796 spl = interpolate.splrep(klgfb['pe'], thisGf) 797 gf = interpolate.splev(sclE, spl) 798 return gf

799 800 # 801 # ---------------------------------------------------------------------------------------- 802 #

803 - def itoh(self, wvl):

804 '''Calculates free-free emission with the free-free gaunt factors of Itoh et al. (ApJS 128, 125, 2000). 805 806 the relativistic values are valid for 6. < log10(T) < 8.5 and -4. < log10(u) < 1.''' 807 wvl = np.array(wvl, 'float64') 808 try: 809 itohCoef = self.ItohCoef 810 except: 811 self.ItohCoef = chio.itohRead()['itohCoef'][self.Z-1].reshape(11, 11) 812 itohCoef = self.ItohCoef 813 try: 814 t = (np.log10(self.Temperature) -7.25)/1.25 815 except: 816 errorMessage = ' temperature undefined in continuum.itoh' 817 print(errorMessage) 818 return {'errorMessage':errorMessage} 819 if type(self.Temperature) == float: 820 nTemp = 1 821 else: 822 nTemp = self.Temperature.size 823 # 824 nWvl = wvl.size 825 # 826 if (nTemp > 1) and (nWvl > 1): 827 u = (const.planck*const.light*1.e+8/const.boltzmann)*np.outer(1./self.Temperature, 1./wvl ) 828 lU = (np.log10(u) + 1.5)/2.5 829 lT = (np.log10(self.Temperature) -7.25)/1.25 830 g = np.zeros((nTemp, nWvl), 'float64') 831 rad = np.ma.zeros((nTemp, nWvl), 'float64') 832 for itemp in range(nTemp): 833 for j in range(11): 834 for i in range(11): 835 g[itemp]+= itohCoef[i,j]*(lT[itemp]**i)*(lU[itemp]**j) 836 rad[itemp] = (1./wvl)**2*g[itemp]*np.exp(-u[itemp])/np.sqrt(self.Temperature[itemp]) 837 tArray = np.zeros((1, len(self.Temperature)), 'float64') 838 tArray[0] = self.Temperature 839 t2Array = np.repeat(tArray.transpose(), len(wvl), axis=1) 840 nonValidT1 = np.log10(t2Array) < 6. 841 nonValidT2 = np.log10(t2Array) > 8.5 842 nonValidT = np.logical_or(nonValidT1, nonValidT2) 843 nonValidU1 = np.log10(u) < -4. 844 nonValidU2 = np.log10(u) > 1. 845 nonValidU = np.logical_or(nonValidU1, nonValidU2) 846 nonValid = np.logical_or(nonValidT, nonValidU) 847 rad.mask = nonValid 848 rad.set_fill_value(0.) 849 g=np.ma.array(g, mask=nonValid, fill_value=0.) 850 return {'itohGff':g, 'itohFf':rad} 851 elif (nTemp > 1) and (nWvl == 1): 852 u = (const.planck*const.light*1.e+8/const.boltzmann)/(self.Temperature*wvl ) 853 lU = (np.log10(u) + 1.5)/2.5 854 lT = (np.log10(self.Temperature) -7.25)/1.25 855 g = np.zeros((nTemp), 'float64') 856 rad = np.ma.zeros((nTemp), 'float64') 857 for itemp in range(nTemp): 858 for j in range(11): 859 for i in range(11): 860 g[itemp]+= itohCoef[i,j]*(lT[itemp]**i)*(lU[itemp]**j) 861 rad[itemp] = (1./wvl)**2*g[itemp]*np.exp(-u[itemp])/np.sqrt(self.Temperature[itemp]) 862 nonValidT1 = np.log10(self.Temperature) < 6. 863 nonValidT2 = np.log10(self.Temperature) > 8.5 864 nonValidT = np.logical_or(nonValidT1, nonValidT2) 865 nonValidU1 = np.log10(u) < -4. 866 nonValidU2 = np.log10(u) > 1. 867 nonValidU = np.logical_or(nonValidU1, nonValidU2) 868 nonValid = np.logical_or(nonValidT, nonValidU) 869 rad.mask = nonValid 870 rad.set_fill_value(0.) 871 g=np.ma.array(g, mask=nonValid, fill_value=0.) 872 return {'itohGff':g, 'itohFf':rad} 873 elif (nTemp == 1) and (nWvl > 1): 874 if (np.log10(self.Temperature) < 6.) or (np.log10(self.Temperature > 8.5)): 875 errorMessage ='invalid temperature in continuum.itoh()' 876 return {'errorMessage':errorMessage} 877 else: 878 u = (const.planck*const.light*1.e+8/const.boltzmann)/(self.Temperature*wvl ) 879 lU = (np.log10(u) + 1.5)/2.5 880 lT = (np.log10(self.Temperature) -7.25)/1.25 881 g = np.zeros(nWvl, 'float64') 882 rad = np.ma.zeros((nWvl), 'float64') 883 for j in range(11): 884 for i in range(11): 885 g+= itohCoef[i,j]*(lT**i)*(lU**j) 886 rad = np.ma.array((1./wvl)**2*g*np.exp(-u)/np.sqrt(self.Temperature), 'Float64') 887 nonValidU1 = np.log10(u) < -4. 888 nonValidU2 = np.log10(u) > 1. 889 nonValidU = np.logical_or(nonValidU1, nonValidU2) 890 nonValid = nonValidU 891 rad.mask = nonValid 892 rad.set_fill_value(0.) 893 g=np.ma.array(g, mask=nonValid, fill_value=0.) 894 return {'itohGff':g, 'itohFf':rad} 895 elif (nTemp == 1) and (nWvl == 1): 896 u = (const.planck*const.light*1.e+8/const.boltzmann)/(self.Temperature*wvl ) 897 if (np.log10(self.Temperature) < 6.) or (np.log10(self.Temperature > 8.5)): 898 errorMessage ='invalid temperature in continuum.itoh()' 899 return {'errorMessage':errorMessage} 900 elif (np.log10(u) < -4.) or (np.log10(u) > 8.5): 901 errorMessage ='invalid temperature/wavelength in continuum.itoh()' 902 return {'errorMessage':errorMessage} 903 else: 904 u = (const.planck*const.light*1.e+8/const.boltzmann)/(self.Temperature*wvl ) 905 lU = (np.log10(u) + 1.5)/2.5 906 lT = (np.log10(self.Temperature) -7.25)/1.25 907 g = np.zeros(nWvl, 'float64') 908 rad = np.ma.zeros((nWvl), 'float64') 909 for j in range(11): 910 for i in range(11): 911 g+= itohCoef[i,j]*(lT**i)*(lU**j) 912 rad = np.ma.array((1./wvl)**2*g*np.exp(-u)/np.sqrt(self.Temperature), 'Float64') 913 nonValidU1 = np.log10(u) < -4. 914 nonValidU2 = np.log10(u) > 1. 915 nonValidU = np.logical_or(nonValidU1, nonValidU2) 916 nonValid = nonValidU 917 rad.mask = nonValid 918 rad.set_fill_value(0.) 919 return {'itohGff':g, 'itohFf':rad}

920 # 921 # - - - - - - - - - - - - - - - - - - - - - - - 922 #

923 - def sutherland(self, wvl):

924 '''Calculates the free-free continuum using the free-free gaunt factor calculations of Sutherland, 1998, MNRAS, 300, 321. 925 926 the wavelengths (wvl) will be sorted, first thing''' 927 # 928 wvl = np.array(wvl, 'float64') 929 nWvl = wvl.size 930 # factor = 5.44436e-39 931 try: 932 temperature = self.Temperature 933 except: 934 errorMessage = ' temperature undefined in continuum.sutherland' 935 print(errorMessage) 936 return {'errorMessage':errorMessage} 937 # read in the gaunt factors, if necessary and get interpolator 938 try: 939 gffInterpolator = self.GffInterpolator 940 except: 941 self.Gff = chio.gffRead() 942 gff = self.Gff 943 iu=(np.log10(gff['u1d']) + 4.)*10. 944 ig=(np.log10(gff['g21d']) + 4.)*5. 945 gaunt = gff['gff'] 946 #tr=Triangulation(iu.flatten(),ig.flatten()) 947 tr = Triangulation(iu,ig) 948 self.GffInterpolator = LinearTriInterpolator(tr, gaunt.flatten()) 949 gffInterpolator = self.GffInterpolator 950 # this is the old version 951 #self.Gff = chio.gffRead() 952 #gff = self.Gff 953 #iu=(np.log10(gff['u1d']) + 4.)*10. 954 #ig=(np.log10(gff['g21d']) + 4.)*5. 955 #gaunt = gff['gff'] 956 ##tr=Triangulation(iu.flatten(),ig.flatten()) 957 #tr=Triangulation(iu,ig) 958 #self.GffInterpolator = tr.nn_interpolator(gaunt.flatten()) 959 #gffInterpolator = self.GffInterpolator 960 # 961 gga = np.array((float(self.Z)**2*const.ryd2erg/const.boltzmann)*(1./temperature),'float64') 962 nonValidGg1 = np.log10(gga) < -4. 963 nonValidGg2 = np.log10(gga) > 4. 964 nonValidGg = np.logical_or(nonValidGg1, nonValidGg2) 965 ggOut = np.ma.array(gga, mask = nonValidGg, fill_value=True) 966 iGg = np.ma.array((np.log10(gga) + 4.)*5., mask=nonValidGg, fill_value=0.) 967 # 968 if nonValidGg.sum(): 969 errorMessage = 'no valid temperatures in continuum.sutherland' 970 print(errorMessage) 971 return {'errorMessage':errorMessage} 972 else: 973 #iUu = np.ma.array((np.log10(uu) + 4.)*10., mask=nonValidUu, fill_value=0.) 974 #iGg = (np.log10(gg) + 4.)*5. 975 #print ' iGg.shape = ',iGg, iGg.shape 976 # 977 nWvl = wvl.size 978 nTemp = temperature.size 979 # 980 if (nTemp > 1) and (nWvl > 1): 981 ff = np.ma.zeros((nWvl, nTemp), 'float64') 982 gffOut1 = np.ma.zeros((nWvl, nTemp), 'float64') 983 gffOutMask = np.zeros((nWvl, nTemp), 'Bool') 984 uuOut = np.zeros((nWvl, nTemp), 'float64') 985 for iwvl in range(nWvl): 986 uu = ((const.planck*const.light*1.e+8/const.boltzmann)/(wvl[iwvl]*temperature)) #.flatten() 987 nonValidUu1 = np.log10(uu) < -4. 988 nonValidUu2 = np.log10(uu) > 4. 989 nonValidUu = np.logical_or(nonValidUu1, nonValidUu2) 990 gffOutMask[iwvl] = nonValidUu 991 uuOut[iwvl] = np.ma.array(uu, mask=nonValidUu, fill_value=True) 992 iUu = np.ma.array((np.log10(uu) + 4.)*10., mask=nonValidUu, fill_value=0.) 993 gffOut1[iwvl] = gffInterpolator(iUu, iGg) 994 wvlt = 1./(wvl[iwvl]**2*np.sqrt(temperature)) # was sortedTemperature 995 ff[iwvl] = (np.exp(-uuOut[iwvl])*gffOut1[iwvl]*wvlt) 996 gffOut1.mask = gffOutMask 997 gffOut1.set_fill_value(0.) 998 gffOut = gffOut1.transpose() 999 ff.mask = gffOutMask 1000 ff.set_fill_value(0.) 1001 return {'suthFf':ff.transpose(), 'suthGff':gffOut} 1002 # 1003 if (nTemp == 1) and (nWvl > 1): 1004 uu = ((const.planck*const.light*1.e+8/const.boltzmann)/(wvl*temperature)) # .flatten() 1005 nonValidUu1 = np.log10(uu) < -4. 1006 nonValidUu2 = np.log10(uu) > 4. 1007 nonValidUu = np.logical_or(nonValidUu1, nonValidUu2) 1008 gffOutMask = nonValidUu 1009 iUu = (np.log10(uu) + 4.)*10. 1010 gffOut1 = gffInterpolator(iUu, iGg.repeat(nWvl)) 1011 wvlt = 1./(wvl**2*np.sqrt(temperature)) 1012 ff = np.ma.array(np.exp(-uu)*gffOut1*wvlt) 1013 ff.mask=gffOutMask 1014 ff.set_fill_value(0.) 1015 gffOut = np.ma.array(gffOut1, mask=gffOutMask, fill_value=0.) 1016 return {'suthFf':ff, 'suthGff':gffOut, 'iUu':iUu, 'gffOut1':gffOut1, 'wvlt':wvlt, 'iGg':iGg.repeat(nWvl), 'gffInterpolator':gffInterpolator} 1017 #elif (nTemp > 1) and (nWvl == 1): 1018 else: 1019 #print ' igg.shape = ',iGg.shape 1020 #gffOut1 = np.ma.zeros((nTemp), 'float64') 1021 #gffOutMask = np.zeros((nTemp), 'Bool') 1022 #uuOut = np.zeros((nTemp), 'float64') 1023 # 1024 uu = (const.planck*const.light*1.e+8/const.boltzmann) /(wvl*temperature).flatten() 1025 nonValidUu1 = np.log10(uu) < -4. 1026 nonValidUu2 = np.log10(uu) > 4. 1027 nonValidUu = np.logical_or(nonValidUu1, nonValidUu2) 1028 gffOutMask = nonValidUu 1029 uuOut = np.ma.array(uu, mask=nonValidUu, fill_value=True) 1030 #iUu = np.ma.array((np.log10(uu) + 4.)*10., mask=nonValidUu, fill_value=0.) 1031 iUu = (np.log10(uu) + 4.)*10. 1032 #print ' iUu.shape = ',iUu.shape 1033 gffOut1 = gffInterpolator(iUu, iGg.flatten()) 1034 # 1035 wvlt = 1./(wvl**2*np.sqrt(temperature)) 1036 ff1 = np.exp(-uuOut)*gffOut1*wvlt 1037 ff = np.ma.array(ff1, mask=gffOutMask, fill_value=0.) 1038 gffOut = np.ma.array(gffOut1, mask=gffOutMask, fill_value=0.) 1039 return {'suthFf':ff, 'suthGff':gffOut}

1040 #elif (nTemp == 1) and (nWvl == 1): 1041 ##iGg = (np.log10(gg) + 4.)*5. 1042 #uu = (const.planck*const.light*1.e+8/const.boltzmann)/(wvl*temperature) 1043 #nonValidUu1 = np.log10(uu) < -4. 1044 #nonValidUu2 = np.log10(uu) > 4. 1045 #nonValidUu = np.logical_or(nonValidUu1, nonValidUu2) 1046 #if not nonValidUu: 1047 #iUu = (np.log10(uu) + 4.)*10. 1048 #gffOut = gffInterpolator(iUu, iGg) 1049 #wvlt = 1./(wvl**2*np.sqrt(temperature)) 1050 #ff = np.exp(-uu)*gffOut*wvlt 1051 #return {'suthFf':ff, 'suthGff':gffOut} 1052 #else: 1053 #errorMessage = 'invalid Temperature/Wavlength in continuum.sutherland' 1054 #print(errorMessage) 1055 #return {'errorMessage':errorMessage} 1056 # 1057 # ---------------------------------------------------------------------------------------- 1058 #

1059 - def ioneqOne(self):

1060 '''Provide the ionization equilibrium for the selected ion as a function of temperature. 1061 returned in self.IoneqOne 1062 this is a duplicate of the method ion.ioneqOne ''' 1063 # 1064 try: 1065 temperature = self.Temperature 1066 except: 1067 return 1068 # 1069 try: 1070 ioneqAll = self.IoneqAll 1071 except: 1072 self.IoneqAll = chio.ioneqRead(ioneqname = self.Defaults['ioneqfile']) 1073 ioneqAll=self.IoneqAll 1074 # 1075 ioneqTemperature = ioneqAll['ioneqTemperature'] 1076 Z=self.Z 1077 Ion=self.Ion 1078 Dielectronic=self.Dielectronic 1079 ioneqOne = np.zeros_like(temperature) 1080 # 1081 thisIoneq=ioneqAll['ioneqAll'][Z-1,Ion-1-Dielectronic].squeeze() 1082 # thisIoneq = self.Ioneq 1083 gioneq=thisIoneq > 0. 1084 goodt1=self.Temperature >= ioneqTemperature[gioneq].min() 1085 goodt2=self.Temperature <= ioneqTemperature[gioneq].max() 1086 goodt=np.logical_and(goodt1,goodt2) 1087 y2=interpolate.splrep(np.log(ioneqTemperature[gioneq]),np.log(thisIoneq[gioneq]),s=0) 1088 # 1089 if goodt.sum() > 0: 1090 if self.Temperature.size > 1: 1091 gIoneq=interpolate.splev(np.log(self.Temperature[goodt]),y2) #,der=0) 1092 ioneqOne[goodt] = np.exp(gIoneq) 1093 else: 1094 gIoneq=interpolate.splev(np.log(self.Temperature),y2) 1095 ioneqOne = np.exp(gIoneq) 1096 else: 1097 ioneqOne = 0. 1098 # 1099 self.IoneqOne = ioneqOne

Source Code for Module chianti.core.Continuum