chianti.core._IonTrails

14 ''' 15 a collection of methods for use in ion and spectrum calculations 16 '''

17 - def __init__(self):

18 pass 19 return

20 # 21 # ------------------------------------------------------------------------------------- 22 # 23 # def gofnt(self,wvlRange=0,top=10, verbose=0): 24 # """ 25 # Calculate the 'so-called' G(T) function. 26 # 27 # Given as a function of both temperature and eDensity. 28 # 29 # Only the top( set by 'top') brightest lines are plotted. 30 # the G(T) function is returned in a dictionary self.Gofnt 31 # """ 32 # # 33 # #self.emiss={"wvl":wvl,"emiss":em,"units":units,"plotLabels":plotLabels} 34 # # 35 # # 36 # # 37 # if not hasattr(self, 'Intensity'): 38 # try: 39 # self.intensity() 40 # except: 41 # print(' intensities not calculated and emiss() is unable to calculate them') 42 # print(' perhaps the temperature and/or eDensity are not set') 43 # return 44 # # 45 # # everything in self.Intensity should be a numpy array 46 # # 47 # # 48 # # 49 # # 50 # # 51 ## if hasattr(self, 'Abundance'): 52 ## ab=self.Abundance 53 ## else: 54 ## self.Abundance = io.abundanceRead() 55 ## ab=self.Abundance 56 ## if not hasattr(self, 'Abundance'): 57 ## self.Abundance = io.abundanceRead() 58 # # 59 # fontsize=12 60 # # 61 ## emiss=em["emiss"] 62 ## wvl=em["wvl"] 63 ## pretty1 = em['pretty1'] 64 ## pretty2 = em['pretty2'] 65 ## lvl1 = em['lvl1'] 66 ## lvl2 = em['lvl2'] 67 # # 68 # intens = copy.deepcopy(self.Intensity) 69 # intensity = intens['intensity'] 70 # ionS = intens['ionS'] 71 # ionSet = set(ionS) 72 # ionNum = len(ionSet) 73 # wvl = intens['wvl'] 74 # lvl1 = intens['lvl1'] 75 # lvl2 = intens['lvl2'] 76 # pretty1 = intens['pretty1'] 77 # pretty2 = intens['pretty2'] 78 # obs = intens['obs'] 79 # avalue = intens['avalue'] 80 # # 81 # temperature=self.Temperature 82 # eDensity=self.EDensity 83 ## plotLabels=em["plotLabels"] 84 ## xLabel=plotLabels["xLabel"] 85 ## yLabel=plotLabels["yLabel"] 86 # # 87 # # find which lines are in the wavelength range if it is set 88 # # 89 # # 90 # if type(wvlRange) != type(1): 91 # igvl=util.between(wvl,wvlRange) 92 # else: 93 # igvl=range(len(wvl)) 94 # nlines=len(igvl) 95 # if nlines ==0: 96 # print(' no lines in selected interval') 97 # return 98 # # find the top most intense lines 99 # # 100 # if top > nlines: 101 # top=nlines 102 # maxEmiss=np.zeros(nlines,'Float64') 103 # for iline in range(nlines): 104 # maxEmiss[iline] = intensity[:, igvl[iline]].max() 105 # for iline in range(nlines): 106 # if maxEmiss[iline] >= maxEmiss.max(): 107 # maxAll = intensity[:, igvl[iline]] 108 ## maxIndex = igvl[iline] 109 ## print ' maxIndex, maxAll = ', maxIndex, maxAll 110 ## line=range(nlines) 111 # igvlsort=np.take(igvl,np.argsort(maxEmiss)) 112 # topLines=igvlsort[-top:] 113 # maxWvl='%5.3f' % wvl[topLines[-1]] 114 ## maxline=topLines[-1] 115 # # 116 # # need to make sure there are no negative values before plotting 117 # good = np.where(intensity > 0.) 118 # emissMin = intensity[good].min() 119 # bad = np.where(intensity <= 0.) 120 # emiss[bad] = emissMin 121 # # 122 # topLines=topLines[wvl[topLines].argsort()] 123 # # 124 # # 125 # ntemp=self.Temperature.size 126 # # 127 # ndens=self.EDensity.size 128 # # 129 # ylabel = 'Emissivity relative to '+maxWvl 130 # title = self.Spectroscopic 131 # # 132 # # 133 # if ndens==1 and ntemp==1: 134 # print(' only a single temperature and eDensity') 135 # return 136 # elif ndens == 1: 137 # xlabel='Temperature (K)' 138 # ngofnt = temperature.size 139 # xvalues=temperature 140 # outTemperature=temperature 141 # outDensity=np.zeros(ntemp,'Float64') 142 # outDensity.fill(eDensity) 143 # desc_str=' at Density = %10.2e' % eDensity 144 # elif ntemp == 1: 145 # xvalues=eDensity 146 # ngofnt = eDensity.size 147 # outTemperature=np.zeros(ndens,'Float64') 148 # outTemperature.fill(temperature) 149 # outDensity=eDensity 150 # xlabel=r'$\rm{Electron Density (cm}^{-3}\rm{)}$' 151 # desc_str=' at Temperature = %10.2e' % temperature 152 # else: 153 # outTemperature=temperature 154 # outDensity=eDensity 155 # xlabel='Temperature (K)' 156 # xvalues=temperature 157 # ngofnt = temperature.size 158 # desc_str=' for variable Density' 159 # # 160 # # 161 # # put all actual plotting here 162 # # 163 # pl.ion() 164 ## if chInteractive: 165 ## pl.ion() 166 ## else: 167 ## pl.ioff() 168 # # 169 # pl.figure() 170 # ax = pl.subplot(111) 171 # nxvalues=len(xvalues) 172 # # maxAll is an array 173 ## print ' emiss = ', np.max(emiss[top-1]), np.max(emiss[0]) 174 ## print ' maxAll = ', maxAll 175 ## ymax = np.max(1.2*emiss[top-1]/maxAll) 176 # ymax = 1.2 177 ## print ' ymax = ', ymax 178 # ymin = ymax # np.min(emiss[0]/maxAll) # was originally = ymax 179 # for iline in range(top): 180 # tline=topLines[iline] 181 # pl.loglog(xvalues,intensity[tline]/maxAll) 182 # if np.min(intensity[:, tline]/maxAll) < ymin: 183 # ymin = np.min(intensity[:, tline]/maxAll) 184 # skip=2 185 # start=divmod(iline,nxvalues)[1] 186 # for ixvalue in range(start,nxvalues,nxvalues/skip): 187 # pl.text(xvalues[ixvalue],intensity[:, tline,ixvalue]/maxAll[ixvalue],str(wvl[tline])) 188 # pl.xlim(xvalues.min(),xvalues.max()) 189 # pl.ylim(ymin, ymax) 190 ## yl=pl.ylim() 191 ## pl.ylim(yl[0],1.2) 192 # pl.xlabel(xlabel,fontsize=fontsize) 193 # pl.ylabel(ylabel,fontsize=fontsize) 194 # if ndens == ntemp and ntemp > 1: 195 # pl.text(0.07, 0.5,title, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 196 # # 197 # ax2 = pl.twiny() 198 # xlabelDen=r'Electron Density (cm$^{-3}$)' 199 # pl.xlabel(xlabelDen, fontsize=fontsize) 200 # pl.loglog(eDensity,intensity[:, topLines[top-1]]/maxAll, visible=False) 201 # ax2.xaxis.tick_top() 202 # else: 203 # pl.ylim(ymin, ymax) 204 # pl.title(title+desc_str,fontsize=fontsize) 205 # pl.draw() 206 # # 207 # time.sleep(0.5) 208 # # 209 ## print ' topLInes = ', wvl[topLines] 210 ## wvlChoices = [] 211 ## for iline in range(top): 212 ## tline = topLines[iline] 213 ## wvlChoices.append('%12.4f %4i %4i %s - %s'%(wvl[tline], lvl1[tline], lvl2[tline], pretty1[tline], pretty2[tline])) 214 ## gline = chGui.gui.selectorDialog(wvlChoices,label='Select line(s)') 215 # # 216 # selectTags = [] 217 # for itop in topLines: 218 # if ionNum == 1: 219 # selectTags.append(str(wvl[itop])) 220 # else: 221 # selectTags.append(ionS[itop]+ ' '+ str(wvl[itop])) 222 # gline = chGui.gui.selectorDialog(selectTags,label='Select line(s)') 223 # 224 # gline_idx=gline.selectedIndex 225 # # 226 # # 227 ## gAbund=self.Abundance 228 # # 229 ## try: 230 ## thisIoneq=self.IoneqOne 231 ## except: 232 ## self.ioneqOne() 233 ## thisIoneq=self.IoneqOne 234 ## if verbose: 235 ## print(' abundance = %10.2e '%(gAbund)) 236 ## print(' index temperature ion fraction') 237 ## for it, anioneq in enumerate(thisIoneq): 238 ## print (' %5i %10.2e %10.2e '%(it, outTemperature[it], anioneq)) 239 # # gioneq=np.where(thisIoneq > 0.) 240 # # y2=interpolate.splrep(np.log(self.IoneqAll['ioneqTemperature'][gioneq]),np.log(thisIoneq[gioneq]),s=0) #allow smoothing,s=0) 241 # # gIoneq=interpolate.splev(np.log(temperature),y2) #,der=0) 242 ## gIoneq=self.IoneqOne/eDensity 243 # # 244 # # 245 # # 246 # # plot the desired ratio 247 # pl.figure() 248 # g_line = topLines[gline_idx]# [0] 249 # ## print ' g_line = ',g_line 250 # # 251 # gofnt=np.zeros(ngofnt,'float64') 252 # for aline in g_line: 253 # gofnt += intensity[:, aline].squeeze() 254 # gofnt += intensity[:, aline].squeeze() 255 # self.Gofnt={'temperature':outTemperature,'eDensity':outDensity,'gofnt':gofnt, 'index':g_line, 'wvl':wvl[g_line]} 256 # # 257 # pl.loglog(xvalues,gofnt) 258 # pl.xlim(xvalues.min(),xvalues.max()) 259 # pl.xlabel(xlabel,fontsize=fontsize) 260 # pl.ylabel('Gofnt',fontsize=fontsize) 261 # if ionNum == 1: 262 # 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]]) 263 # if len(g_line) > 1: 264 # newTitle +='\n' 265 # for igl in g_line[1:]: 266 # newTitle += ' ' + '%12.3f %4i %4i %s - %s'%(wvl[igl], lvl1[igl], lvl2[igl], pretty1[igl], pretty2[igl]) 267 # if igl != g_line[-1]: 268 # newTitle +='\n' 269 # # pl.annotate(newTitle, xytext=(0.3, 0.3), textcoords='figure_fraction') 270 # else: 271 # newTitle = 'newTitle' 272 # pl.annotate(newTitle, xy=(-10, 10), 273 # xycoords='axes points', 274 # horizontalalignment='right', verticalalignment='bottom') #,fontsize=20) 275 # if ndens == ntemp and ntemp > 1: 276 ## newTitle +' '+str(wvl[g_line])+' '+desc_str 277 # pl.text(0.07, 0.5,newTitle, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 278 # # 279 # ax2 = pl.twiny() 280 ## xlabel=r'Electron Density (cm$^{-3}$)' 281 # pl.xlabel(xlabelDen, fontsize=fontsize) 282 # pl.loglog(eDensity,gofnt, visible=False) 283 # ax2.xaxis.tick_top() 284 # else: 285 # pl.title(newTitle, fontsize=fontsize) 286 # #pl.ioff() 287 # #pl.show() 288 ## return 289 # 290 # --------------------------------------------------------------------------- 291 #

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

293 ''' 294 List the line intensities 295 296 wvlRange, a 2 element tuple, list or array determines the wavelength range 297 298 Top specifies to plot only the top strongest lines, default = 10 299 300 normalize = 1 specifies whether to normalize to strongest line, default = 0 301 rewrite of emissList 302 ''' 303 # 304 # 305 if not hasattr(self, 'Intensity'): 306 try: 307 self.intensity() 308 except: 309 print(' intensities not calculated and emiss() is unable to calculate them') 310 print(' perhaps the temperature and/or eDensity are not set') 311 return 312 # 313 # everything in self.Intensity should be a numpy array 314 # 315 # 316 # 317 temperature = self.Temperature 318 eDensity = self.EDensity 319 em = self.Em 320 # 321 ndens = eDensity.size 322 ntemp = temperature.size 323 # 324 # 325 intens = copy.deepcopy(self.Intensity) 326 intensity = intens['intensity'] 327 ionS = intens['ionS'] 328 wvl = intens['wvl'] 329 lvl1 = intens['lvl1'] 330 lvl2 = intens['lvl2'] 331 pretty1 = intens['pretty1'] 332 pretty2 = intens['pretty2'] 333 obs = intens['obs'] 334 avalue = intens['avalue'] 335 # 336 if ndens == 1 and ntemp == 1: 337 dstr = ' - Density = %10.2e (cm$^{-3}$)' %(eDensity) 338 tstr = ' - T = %10.2e (K)' %(temperature) 339 elif ndens == 1 and ntemp > 1: 340 if index < 0: 341 index = ntemp/2 342 print('using index = %5i specifying temperature = %10.2e'%(index, temperature[index])) 343 self.Message = 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 344 345 intensity=intensity[index] 346 elif ndens > 1 and ntemp == 1: 347 if index < 0: 348 index = ntemp/2 349 print('using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index])) 350 self.Message = 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 351 intensity=intensity[index] 352 elif ndens > 1 and ntemp > 1: 353 if index < 0: 354 index = ntemp/2 355 print('using index = %5i specifying temperature = %10.2e, eDensity = %10.2e em = %10.2e'%(index, temperature[index], eDensity[index], em[index])) 356 self.Message = 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2eem = %10.2e'%(index, temperature[index], eDensity[index], em[index]) 357 intensity=intensity[index] 358 # 359 # print('shpae of intensity 0 = %5i %5i'%(intensity.shape)) 360 # 361 if wvlRange: 362 wvlIndex=util.between(wvl,wvlRange) 363 elif wvlRanges: 364 wvlIndex = [] 365 for awvlRange in wvlRanges: 366 wvlIndex.extend(util.between(wvl,awvlRange)) 367 else: 368 wvlIndex = range(wvl.size) 369 # 370 # get lines in the specified wavelength range 371 # 372 # print('shpae of intensity 1 = %5i %5i'%(intensity.shape)) 373 # 374 intensity = intensity[wvlIndex] 375 ionS = ionS[wvlIndex] 376 wvl = wvl[wvlIndex] 377 lvl1 = lvl1[wvlIndex] 378 lvl2 = lvl2[wvlIndex] 379 avalue = avalue[wvlIndex] 380 pretty1 = pretty1[wvlIndex] 381 pretty2 = pretty2[wvlIndex] 382 obs = obs[wvlIndex] 383 # 384 self.Error = 0 385 if wvl.size == 0: 386 print('No lines in this wavelength interval') 387 self.Error = 1 388 self.Message = 'No lines in this wavelength interval' 389 return 390 # 391 elif top == 0: 392 top = wvl.size 393 elif top > wvl.size: 394 top = wvl.size 395 # 396 # sort by intensity 397 # 398 isrt = np.argsort(intensity) 399 ionS = ionS[isrt[-top:]] 400 wvl = wvl[isrt[-top:]] 401 lvl1 = lvl1[isrt[-top:]] 402 lvl2 = lvl2[isrt[-top:]] 403 obs = obs[isrt[-top:]] 404 intensity = intensity[isrt[-top:]] 405 avalue = avalue[isrt[-top:]] 406 pretty1 = pretty1[isrt[-top:]] 407 pretty2 = pretty2[isrt[-top:]] 408 # 409 # print('shpae of intensity 2= %5i %5i'%(intensity.shape)) 410 # 411 # must follow setting top 412 # 413 if relative: 414 intensity = intensity/intensity[:top].max() 415 # 416 # 417 idx = np.argsort(wvl) 418 fmt1 = '%5s %5s %5s %25s - %-25s %12s %12s %12s %3s' 419 fmt = '%5s %5i %5i %25s - %-25s %12.4f %12.3e %12.2e %1s' 420 print(' ') 421 print(' ------------------------------------------') 422 print(' ') 423 print(fmt1%('Ion','lvl1','lvl2','lower','upper','Wvl(A)','Intensity','A value','Obs')) 424 for kdx in idx: 425 print(fmt%(ionS[kdx], lvl1[kdx], lvl2[kdx], pretty1[kdx], pretty2[kdx], wvl[kdx], intensity[kdx], avalue[kdx], obs[kdx])) 426 print(' ') 427 print(' ------------------------------------------') 428 print(' ') 429 # 430 self.Intensity['wvlTop'] = wvl[idx] 431 self.Intensity['intensityTop'] = intensity[idx] 432 if outFile: 433 fmt1a = '%5s %5s %5s %25s - %-25s %12s %12s %12s %3s \n' 434 fmt = '%5s %5i %5i %25s - %-25s %12.4f %12.3e %12.2e %1s \n' 435 outpt = open(outFile, 'w') 436 outpt.write(fmt1a%('Ion','lvl1','lvl2','lower','upper','Wvl(A)','Intensity','A value','Obs')) 437 for kdx in idx: 438 outpt.write(fmt%(ionS[kdx], lvl1[kdx], lvl2[kdx], pretty1[kdx], pretty2[kdx], wvl[kdx], intensity[kdx], avalue[kdx], obs[kdx])) 439 outpt.close()

440 # 441 # --------------------------------------------------------------------------- 442 #

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

444 '''Plot the line intensities. 445 446 wvlRange, a 2 element tuple, list or array determines the wavelength range 447 448 Top specifies to plot only the top strongest lines, default = 10 449 450 linLog specifies a linear or log plot, want either lin or log, default = lin 451 452 normalize = 1 specifies whether to normalize to strongest line, default = 0''' 453 # 454 title=self.Spectroscopic 455 # 456 if hasattr(self, 'Intensity'): 457 intens = self.Intensity['intensity'] 458 else: 459 try: 460 self.intensity(em=em) 461 intens = self.Intensity['intensity'] 462 except: 463 print(' emissivities not calculated and emiss() is unable to calculate them') 464 print(' perhaps the temperature and/or eDensity are not set') 465 return 466 wvl = self.Intensity['wvl'] 467 temperature = self.Temperature 468 eDensity = self.EDensity 469 # 470 ndens = eDensity.size 471 ntemp = temperature.size 472 # 473 if ndens == 1 and ntemp == 1: 474 dstr = ' - Density = %10.2e (cm$^{-3}$)' %(eDensity) 475 tstr = ' - T = %10.2e (K)' %(temperature) 476 elif ndens == 1 and ntemp > 1: 477 if index < 0: 478 index = ntemp/2 479 print('using index = %5i specifying temperature = %10.2e'%(index, temperature[index])) 480 self.Message = 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 481 # if chInteractive: 482 # print 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 483 # else: 484 # self.Message = 'using index = %5i specifying temperature = %10.2e'%(index, temperature[index]) 485 intensity=self.Intensity['intensity'][index] 486 dstr=' - Density = %10.2e (cm$^{-3}$)' % eDensity 487 tstr=' - T = %10.2e (K)' % temperature[index] 488 elif ndens > 1 and ntemp == 1: 489 if index < 0: 490 index = ntemp/2 491 print('using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index])) 492 self.Message = 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 493 # if chInteractive: 494 # print 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 495 # else: 496 # self.Message = 'using index =%5i specifying eDensity = %10.2e'%(index, eDensity[index]) 497 intensity = self.Intensity['intensity'][index] 498 dstr=' - Density = %10.2e (cm$^{-3}$)' % eDensity[index] 499 tstr=' - T = %10.2e (K)' % temperature 500 elif ndens > 1 and ntemp > 1: 501 if index < 0: 502 index = ntemp/2 503 print('using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index])) 504 self.Message = 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index]) 505 # if chInteractive: 506 # print 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index]) 507 # else: 508 # self.Message = 'using index = %5i specifying temperature = %10.2e, eDensity = %10.2e'%(index, temperature[index], eDensity[index]) 509 intensity = self.Intensity['intensity'][index] 510 dstr=' - Density = %10.2e (cm$^{-3}$)' % eDensity[index] 511 tstr=' - T = %10.2e (K)' % temperature[index] 512 if type(wvlRange) != type(None): 513 wvlIndex = util.between(wvl, wvlRange) 514 else: 515 wvlIndex = range(wvl.size) 516 intensity = intensity[wvlIndex] 517 wvl = wvl[wvlIndex] 518 # 519 self.Error = 0 520 if wvl.size == 0: 521 print('No lines in this wavelength interval') 522 self.Error = 1 523 self.Message = 'No lines in this wavelength interval' 524 # if chInteractive: 525 # print 'No lines in this wavelength interval' 526 # else: 527 # self.Error = 1 528 # self.Message = 'No lines in this wavelength interval' 529 return 530 elif top == 0: 531 top = wvl.size 532 elif wvl.size > top: 533 intsrt = np.argsort(intensity) 534 wvl = wvl[intsrt[-top:]] 535 intensity = intensity[intsrt[-top:]] 536 else: 537 top = wvl.size 538 # must follow setting top 539 # 540 pl.figure() 541 ylabel = 'intensity' 542 if relative: 543 intensity = intensity/intensity[-1] 544 ylabel += ' (Relative)' 545 # 546 xlabel = 'Wavelength ('+self.Defaults['wavelength'] +')' 547 # 548 ymin = 10.**(np.log10(intensity[0].min()).round(0)-0.5 ) 549 # 550 pl.ion() 551 # if chInteractive: 552 # pl.ion() 553 # else: 554 # pl.ioff() 555 # 556 for idx in range(top): 557 xx=[wvl[idx], wvl[idx]] 558 if linLog == 'lin': 559 yy=[0., intensity[idx]] 560 pl.plot(xx, yy) 561 else: 562 yy=[ymin/10., intensity[idx]] 563 pl.semilogy(xx, yy) 564 pl.xlabel(xlabel) 565 pl.ylabel(ylabel) 566 pl.title(title+tstr+dstr) 567 if wvlRange: 568 pl.axis([wvlRange[0], wvlRange[1], ymin, intensity.max()]) 569 if plotFile: 570 pl.savefig(plotFile)

571 # 572 # idx = np.argsort(wvl) 573 # self.Intens['wvlTop'] = wvl[idx] 574 # self.Intensity['intensTop'] = intensity[idx] 575 # 576 # ------------------------------------------------------------------------------------- 577 #

578 - def intensityRatio(self,wvlRange=None, wvlRanges=None,top=10):

579 """ 580 Plot the ratio of 2 lines or sums of lines. 581 Shown as a function of density and/or temperature. 582 For a single wavelength range, set wvlRange = [wMin, wMax] 583 For multiple wavelength ranges, set wvlRanges = [[wMin1,wMax1],[wMin2,wMax2], ...] 584 A plot of relative emissivities is shown and then a dialog appears for the user to 585 choose a set of lines. 586 """ 587 # 588 # self.Emiss={"temperature":temperature,"density":density,"wvl":wvl,"emiss":em, 589 # "plotLabels":plotLabels} 590 # 591 # 592 if not hasattr(self, 'Intensity'): 593 try: 594 self.intensity() 595 except: 596 print(' intensities not calculated and emiss() is unable to calculate them') 597 print(' perhaps the temperature and/or eDensity are not set') 598 return 599 # 600 # everything in self.Intensity should be a numpy array 601 # 602 # intens = copy.copy(self.Intensity) 603 # intensity = intens['intensity'] 604 # 605 # 606 fontsize=14 607 # 608 # temperature = self.Temperature 609 eDensity = self.EDensity 610 temperature = self.Temperature 611 # intensity = intens['intensity'] 612 # 613 # temp=np.asarray(temperature,'Float32') 614 ntemp = temperature.size 615 if ntemp > 0: 616 if temperature[0] == temperature[-1]: 617 ntemp = 1 618 # 619 ndens = eDensity.size 620 if ndens > 0: 621 if eDensity[0] == eDensity[-1]: 622 ndens = 1 623 # 624 print(' ndens = %5i ntemp = %5i'%(ndens, ntemp)) 625 626 ionS = self.Intensity['ionS'] 627 # see if we are dealing with more than a single ion 628 ionSet = set(ionS) 629 ionNum = len(ionSet) 630 wvl = self.Intensity["wvl"] 631 # plotLabels = intens["plotLabels"] 632 # xLabel = plotLabels["xLabel"] 633 # yLabel = plotLabels["yLabel"] 634 # 635 # find which lines are in the wavelength range if it is set 636 # 637 # 638 if wvlRange: 639 igvl=util.between(wvl,wvlRange) 640 if len(igvl) == 0: 641 print('no lines in wavelength range %12.2f - %12.2f'%(wvlRange[0], wvlRange[1])) 642 return 643 elif wvlRanges: 644 igvl = [] 645 for awvlRange in wvlRanges: 646 igvl.extend(util.between(wvl,awvlRange)) 647 if len(igvl) == 0: 648 print('no lines in wavelength ranges specified ') 649 return 650 else: 651 igvl=range(len(wvl)) 652 # 653 nlines=len(igvl) 654 # 655 # print ' nlines = ',nlines 656 # print ' iglv = ',igvl 657 igvl=np.take(igvl,wvl[igvl].argsort()) 658 # find the top most intense lines 659 # 660 if top > nlines: 661 top=nlines 662 # 663 intensity = self.Intensity['intensity'] 664 maxIntens = np.zeros(nlines,'Float64') 665 for iline in range(nlines): 666 maxIntens[iline] = intensity[:, igvl[iline]].max() 667 for iline in range(nlines): 668 if maxIntens[iline]==maxIntens.max(): 669 maxAll=intensity[:, igvl[iline]] 670 # line=range(nlines) 671 igvlsort=np.take(igvl,np.argsort(maxIntens)) 672 # print 'igvlsort = ', igvlsort 673 topLines=igvlsort[-top:] 674 # print ' topLines = ', topLines 675 maxWvl='%5.3f' % wvl[topLines[-1]] 676 # maxline=topLines[-1] 677 # 678 topLines=topLines[wvl[topLines].argsort()] 679 # 680 # 681 # need to make sure there are no negative values before plotting 682 good = intensity > 0. 683 intensMin = intensity[good].min() 684 bad = intensity <= 0. 685 intensity[bad] = intensMin 686 # 687 # 688 # ntemp=self.Temperature.size 689 # 690 # ndens=self.EDensity.size 691 # 692 ylabel='Intensity relative to '+maxWvl 693 if ionNum == 1: 694 title=self.Spectroscopic 695 else: 696 title = '' 697 # 698 # 699 if ndens==1 and ntemp==1: 700 print(' only a single temperature and eDensity') 701 return 702 elif ndens == 1: 703 xlabel='Temperature (K)' 704 xvalues=self.Temperature 705 outTemperature=self.Temperature 706 # outDensity=np.zeros(ntemp,'Float64') 707 # outDensity.fill(self.EDensity[0]) 708 outDensity = self.EDensity 709 desc_str=' Density = %10.2e (cm)$^{-3}$' % self.EDensity[0] 710 elif ntemp == 1: 711 xvalues=self.EDensity 712 # outTemperature=np.zeros(ndens,'Float64') 713 # outTemperature.fill(self.Temperature[0]) 714 outTemperature = self.Temperature 715 outDensity=self.EDensity 716 xlabel=r'$\rm{Electron Density (cm)^{-3}}$' 717 desc_str=' Temp = %10.2e (K)' % self.Temperature[0] 718 else: 719 outTemperature=self.Temperature 720 outDensity=self.EDensity 721 xlabel='Temperature (K)' 722 xvalues=self.Temperature 723 desc_str=' Variable Density' 724 # 725 # put all actual plotting here 726 # 727 pl.ion() 728 # if chInteractive: 729 # pl.ion() 730 # else: 731 # pl.ioff() 732 # 733 # maxAll is an array 734 ymax = np.max(intensity[:, topLines[0]]/maxAll) 735 ymin = ymax 736 pl.figure() 737 ax = pl.subplot(111) 738 nxvalues=len(xvalues) 739 for iline in range(top): 740 tline=topLines[iline] 741 pl.loglog(xvalues,intensity[:, tline]/maxAll) 742 if np.min(intensity[:, tline]/maxAll) < ymin: 743 ymin = np.min(intensity[:, tline]/maxAll) 744 if np.max(intensity[:, tline]/maxAll) > ymax: 745 ymax = np.max(intensity[:, tline]/maxAll) 746 skip=2 747 start=divmod(iline,nxvalues)[1] 748 for ixvalue in range(start,nxvalues,nxvalues//skip): 749 if ionNum == 1: 750 text = '%10.4f'%(wvl[tline]) 751 else: 752 text = '%s %10.4f'%(ionS[tline], wvl[tline]) 753 pl.text(xvalues[ixvalue], intensity[ixvalue, tline]/maxAll[ixvalue], text) 754 pl.xlim(xvalues.min(),xvalues.max()) 755 # pl.ylim(ymin, ymax) 756 pl.xlabel(xlabel,fontsize=fontsize) 757 pl.ylabel(ylabel,fontsize=fontsize) 758 if ndens == ntemp and ntemp > 1: 759 pl.text(0.07, 0.5,title, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 760 # 761 ax2 = pl.twiny() 762 xlabelDen=r'Electron Density (cm$^{-3}$)' 763 pl.xlabel(xlabelDen, fontsize=fontsize) 764 pl.loglog(eDensity,intensity[:, topLines[top-1]]/maxAll, visible=False) 765 ax2.xaxis.tick_top() 766 pl.ylim(ymin/1.2, 1.2*ymax) 767 else: 768 pl.ylim(ymin/1.2, 1.2*ymax) 769 pl.title(title+desc_str,fontsize=fontsize) 770 pl.draw() 771 # need time to let matplotlib finish plotting 772 time.sleep(0.5) 773 # 774 # get line selection 775 # 776 selectTags = [] 777 for itop in topLines: 778 if ionNum == 1: 779 selectTags.append(str(wvl[itop])) 780 else: 781 selectTags.append(ionS[itop]+ ' '+ str(wvl[itop])) 782 # 783 # numden = chGui.gui.choice2Dialog(wvl[topLines]) 784 numden = chGui.gui.choice2Dialog(selectTags) 785 # 786 # num_idx and den_idx are tuples 787 # 788 num_idx=numden.numIndex 789 if len(num_idx) == 0: 790 print(' no numerator lines were selected') 791 return 792 # 793 den_idx=numden.denIndex 794 if len(den_idx) == 0: 795 print(' no denominator lines were selected') 796 return 797 # 798 numIntens=np.zeros(len(xvalues),'Float64') 799 for aline in num_idx: 800 numIntens += intensity[:, topLines[aline]] 801 # 802 denIntens = np.zeros(len(xvalues),'Float64') 803 for aline in den_idx: 804 denIntens += intensity[:, topLines[aline]] 805 # 806 # plot the desired ratio 807 # maxAll is an array 808 pl.figure() 809 ax = pl.subplot(111) 810 pl.loglog(xvalues,numIntens/denIntens) 811 pl.xlim(xvalues.min(),xvalues.max()) 812 pl.xlabel(xlabel,fontsize=fontsize) 813 pl.ylabel('Ratio ('+self.Defaults['flux']+')',fontsize=fontsize) 814 if ionNum == 1: 815 desc = ionS[0] 816 else: 817 desc = '' 818 for aline in num_idx: 819 if ionNum == 1: 820 desc += ' ' + str(wvl[topLines[aline]]) 821 else: 822 desc += ' ' + ionS[topLines[aline]] + ' ' + str(wvl[topLines[aline]]) 823 desc += ' / ' 824 for aline in den_idx: 825 if ionNum == 1: 826 desc += ' ' + str(wvl[topLines[aline]]) 827 else: 828 desc += ' ' + ionS[topLines[aline]] + ' ' + str(wvl[topLines[aline]]) 829 if ndens == ntemp and ntemp > 1: 830 pl.text(0.07, 0.5,desc, horizontalalignment='left', verticalalignment='center', fontsize=fontsize, transform = ax.transAxes) 831 # 832 ax2 = pl.twiny() 833 xlabelDen=r'Electron Density (cm$^{-3}$)' 834 pl.xlabel(xlabelDen, fontsize=fontsize) 835 pl.loglog(eDensity,numIntens/denIntens, visible=False) 836 ax2.xaxis.tick_top() 837 else: 838 # pl.ylim(ymin, ymax) 839 pl.title(desc,fontsize=fontsize) 840 # desc=title+' '+str(wvl[num_line])+' / '+str(wvl[den_line])+' '+desc_str 841 # pl.title(desc, fontsize=fontsize) 842 # pl.title(title+' '+str(wvl[num_line])+' / '+str(wvl[den_line])+' '+desc_str,fontsize=fontsize) 843 # pl.draw() 844 # pl.ioff() 845 # pl.show() 846 # 847 # intensityRatioFileName=self.IonStr 848 # for aline in num_idx: 849 # intensityRatioFileName+= '_%3i'%(wvl[topLines[aline]]) 850 # intensityRatioFileName+='_2' 851 # for aline in den_idx: 852 # intensityRatioFileName+= '_%3i'%(wvl[topLines[aline]]) 853 cnt = desc.count(' ') 854 intensityRatioFileName = desc.replace(' ', '_', cnt) + '.rat' 855 intensityRatioFileName = intensityRatioFileName.lstrip('_').replace('_/_','-') 856 self.IntensityRatio={'ratio':numIntens/denIntens,'desc':desc, 857 'temperature':outTemperature,'eDensity':outDensity,'filename':intensityRatioFileName, 'numIdx':num_idx, 'denIdx':den_idx}

858 # 859 # ------------------------------------------------------------------------------------- 860 #

861 - def intensityRatioSave(self,outFile=0):

862 ''' 863 Save the intensity ratio to a file. 864 865 The intensity ratio as a function to temperature and eDensity is saved to an asciii file. 866 867 Descriptive information is included at the top of the file. 868 ''' 869 if not outFile: 870 outFile = self.IntensityRatio['filename'] 871 # if chInteractive: 872 print(' saving ratio to filename = %s'%(outFile)) 873 if hasattr(self, 'IntensityRatio'): 874 temperature=self.IntensityRatio['temperature'] 875 eDensity=self.IntensityRatio['eDensity'] 876 ratio=self.IntensityRatio['ratio'] 877 out = open(outFile,'w') 878 nvalues=len(ratio) 879 # 880 # need to add 7 lines to maintain IDL like files 881 # 882 out.write(outFile+'\n') #1 883 out.write(self.IntensityRatio['desc']+'\n') #2 884 out.write(' created with ChiantiPy using CHIANTI version '+ chdata.ChiantiVersion +'\n') #3 885 out.write(' columns are temperature, eDensity, ratio'+'\n') #5 886 tunit = 'K' 887 out.write(' temperature in '+tunit+', electron eDensity in cm^(-3)'+'\n') #6 888 out.write(' ratio given in '+self.Defaults['flux']+'\n') #4 889 out.write(' '+'\n') #7 890 for ivalue in range(nvalues): 891 s='%12.3e %12.3e %12.3e \n' % (temperature[ivalue],eDensity[ivalue],ratio[ivalue]) 892 out.write(s) 893 out.close() 894 else: 895 # if chInteractive: 896 print(' in .intensityRatioSave(), no IntensityRatio is found')

Source Code for Module chianti.core._IonTrails