1 '''
2 a collection of reading and writing functions
3 seems to be Python 2.7 and Python 3.4 compatible
4 '''
5 import os, fnmatch, sys
6 import pickle
7 if sys.version_info[0] == 3:
8
9 import configparser
10 elif sys.version_info[0] == 2:
11
12 import ConfigParser as configparser
13
14 import numpy as np
15 import chianti.util as util
16 import chianti.constants as const
17 import chianti.Gui as chgui
18 from chianti.fortranformat import FortranRecordReader
19
20
21
23 """ read an abundanc file and returns the abundance values relative to hydrogen"""
24 abundance=np.zeros((50),'Float64')
25 xuvtop=os.environ["XUVTOP"]
26 if abundancename!='':
27
28 abundancefile=os.path.join(xuvtop,'abundance',abundancename+'.abund')
29 else:
30
31 abundir=os.path.join(xuvtop,'abundance')
32 abundlabel = 'ChiantiPy - Select an abundance file'
33
34 fname = chgui.gui.chpicker(abundir, filter='*.abund', label=abundlabel)
35 if fname == None:
36 print((' no abundance file selected'))
37 return 0
38 else:
39 abundancefile = os.path.join(abundir, fname)
40 abundancefilename=os.path.basename(fname)
41 abundancename,ext=os.path.splitext(abundancefilename)
42
43
44
45
46 input=open(abundancefile,'r')
47 s1=input.readlines()
48 input.close()
49 nlines=0
50 idx=-1
51 while idx <= 0:
52 minChar = min([5, len(s1[nlines])])
53 aline=s1[nlines][0:minChar]
54 idx=aline.find('-1')
55 nlines+=1
56 nlines-=1
57 for line in range(nlines):
58 z,ab,element=s1[line].split()
59 abundance[int(z)-1]=float(ab)
60 gz=np.nonzero(abundance)
61 abs=10.**(abundance[gz]-abundance[0])
62 abundance.put(gz,abs)
63 abundanceRef=s1[nlines+1:]
64 return {'abundancename':abundancename,'abundance':abundance,'abundanceRef':abundanceRef}
65
66
67
69 """
70 convert Z, ion to generic name 26, 13 -> fe_13
71 a duplicate of teh routine in util but needed by masterList Info
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
83 thisone = 0
84 return thisone
85
86
87
89 """
90 convert ion name string to Z and Ion
91 a duplicate of teh routine in util but needed by masterList Info
92 """
93 s2=name.split('_')
94 els=s2[0].strip()
95 i1=const.El.index(els)+1
96 ions=s2[1].strip()
97 d=ions.find('d')
98 if d >0 :
99 dielectronic=True
100 ions=ions.replace('d','')
101 else: dielectronic=False
102 higher = zion2name(int(i1), int(ions)+1)
103 lower = zion2name(int(i1), int(ions)-1)
104 return {'Z':int(i1),'Ion':int(ions),'Dielectronic':dielectronic, 'Element':els, 'higher':higher, 'lower':lower}
105
106
107
108 -def cireclvlRead(ions, filename=0, cilvl=0, reclvl=0, rrlvl=0):
109 '''
110 to read Chianti cilvl and reclvl files and return data
111 must specify type as either cilvl, reclvl or rrlvl
112 '''
113 if filename:
114 fname = filename
115 else:
116 fname = util.ion2filename(ions)
117 if cilvl:
118 paramname=fname+'.cilvl'
119 elif reclvl:
120 paramname = fname + '.reclvl'
121 elif rrlvl:
122 paramname = fname + '.rrlvl'
123 else:
124 print('either "cilvl", "reclvl" ir "rrlvl" must be specified')
125 return {}
126 if os.path.exists(paramname):
127 input=open(paramname,'r')
128 lines = input.readlines()
129 input.close()
130 else:
131 print(('file does not exist: ', paramname))
132 return {'error':'file does not exist: ' + paramname}
133
134 iline = 0
135 idx = -1
136 while idx < 0:
137 aline=lines[iline][0:5]
138 idx=aline.find('-1')
139 iline += 1
140 ndata = iline - 1
141 ntrans = ndata/2
142
143
144
145
146
147
148
149
150
151
152
153
154 ntemp = np.zeros(ntrans, 'int32')
155 iline = 0
156 for jline in range(0, ndata, 2):
157 dummy = lines[jline].replace(os.linesep, '').split()
158 ntemp[iline] = len(dummy[4:])
159 iline += 1
160 maxNtemp = ntemp.max()
161
162 temp = np.zeros((ntrans,maxNtemp), 'float64')
163 iline = 0
164 for jline in range(0, ndata, 2):
165 recdat = lines[jline].replace(os.linesep, '').split()
166 shortT = np.asarray(recdat[4:], 'float64')
167
168 t = np.resize(shortT, maxNtemp)
169 if rrlvl:
170 temp[iline] = t
171 else:
172 temp[iline] = 10.**t
173 iline += 1
174
175 lvl1 = np.zeros(ntrans, 'int64')
176 lvl2 = np.zeros(ntrans, 'int64')
177 ci = np.zeros((ntrans, maxNtemp), 'float64')
178
179 idat = 0
180 for jline in range(1, ndata, 2):
181 cidat = lines[jline].replace(os.linesep, '').split()
182 shortCi = np.asarray(cidat[4:], 'float64')
183 lvl1[idat] = int(cidat[2])
184 lvl2[idat] = int(cidat[3])
185 ci[idat] = np.resize(shortCi, maxNtemp)
186 idat += 1
187 return {'temperature':temp, 'ntemp':ntemp,'lvl1':lvl1, 'lvl2':lvl2, 'rate':ci,'ref':lines[ndata+1:], 'ionS':ions}
188
189
190
192 '''
193 possibleDefaults = {'wavelength':['angstrom', 'kev', 'nm']}
194 symbolDefaults = {'wavelength':['A', 'keV', 'nm']}
195 '''
196 initDefaults={'abundfile': 'sun_photospheric_1998_grevesse','ioneqfile': 'chianti', 'wavelength': 'angstrom', 'flux': 'energy','gui':False}
197 rcfile=os.path.join(os.environ['HOME'],'.chianti/chiantirc')
198 if os.path.isfile(rcfile):
199 print((' reading chiantirc file'))
200 config = configparser.RawConfigParser(initDefaults)
201 config.read(rcfile)
202 defaults = {}
203 for anitem in config.items('chianti'):
204 defaults[anitem[0]] = anitem[1]
205 if defaults['gui'].lower() in ('t', 'y', 'yes', 'on', 'true', '1', 1, True):
206 defaults['gui'] = True
207 elif defaults['gui'].lower() in ('f', 'n', 'no', 'off', 'false', '0', 0, False):
208 defaults['gui'] = False
209 else:
210 defaults = initDefaults
211 if verbose:
212 print((' chiantirc file (/HOME/.chianti/chiantirc) does not exist'))
213 print((' using the following defaults'))
214 for akey in list(defaults.keys()):
215 print((' %s = %s'%(akey, defaults[akey])))
216 return defaults
217
218
219
221 """
222 read chianti direct ionization .params files and return
223 {"info":info,"btf":btf,"ev1":ev1,"xsplom":xsplom,"ysplom":ysplom,"ref":hdr}
224 info={"iz":iz,"ion":ion,"nspl":nspl,"neaev":neaev}
225 cannot read dilvlparams files
226 """
227
228 if filename:
229 paramname = filename
230 else:
231 zion = util.convertName(ions)
232 if zion['Z'] < zion['Ion']:
233 print(' this is a bare nucleus that has no ionization rate')
234 return
235
236 fname = util.ion2filename(ions)
237 paramname=fname+'.diparams'
238
239 input=open(paramname,'r')
240
241 line1=input.readline()
242 indices=line1.split()
243 iz=int(indices[0])
244 ion=int(indices[1])
245 nspl=indices[2]
246 nfac=int(indices[3])
247 neaev=int(indices[4])
248 nspl=int(nspl)
249
250 header_line = FortranRecordReader(str(nspl+1)+'E10.2')
251
252
253
254 ev1=np.zeros(nfac,'Float64')
255 btf=np.zeros(nfac,'Float64')
256 xsplom=np.zeros([nfac, nspl],'Float64')
257 ysplom=np.zeros([nfac, nspl],'Float64')
258
259 for ifac in range(nfac):
260 line=input.readline()
261
262 paramdat = header_line.read(line)
263 btf[ifac]=paramdat[0]
264 xsplom[ifac]=paramdat[1:]
265 line=input.readline()
266
267 paramdat = header_line.read(line)
268 ev1[ifac]=paramdat[0]
269 ysplom[ifac]=paramdat[1:]
270 if neaev:
271 line=input.readline()
272 eacoef=line.split()
273
274 eaev=[float(avalue) for avalue in eacoef]
275
276
277
278
279
280
281
282 else:
283 eaev=0.
284 hdr=input.readlines()
285 input.close()
286 info={"iz":iz,"ion":ion,"nspl":nspl,"neaev":neaev, 'nfac':nfac}
287 if neaev:
288 info['eaev'] = eaev
289 DiParams={"info":info,"btf":btf,"ev1":ev1,"xsplom":xsplom,"ysplom":ysplom, 'eaev':eaev,"ref":hdr}
290 return DiParams
291
292
293
295 """
296 read chianti dielectronic recombination .drparams files and return
297 {'rrtype','params','ref'}
298 """
299
300
301 fname = util.ion2filename(ions)
302 paramname=fname+'.drparams'
303 if os.path.isfile(paramname):
304 input=open(paramname,'r')
305
306 lines=input.readlines()
307 input.close()
308 drtype=int(lines[0])
309 ref=lines[4:-1]
310
311 if drtype == 1:
312
313
314 header_line = FortranRecordReader('2i5,8e12.4')
315 inpt1 = header_line.read(lines[1])
316
317 eparams = np.asarray(inpt1[2:], 'float64')
318 inpt2 = header_line.read(lines[2])
319
320 cparams = np.asarray(inpt2[2:], 'float64')
321 DrParams={'drtype':drtype, 'eparams':eparams,'cparams':cparams, 'ref':ref}
322 elif drtype == 2:
323
324
325 header_line = FortranRecordReader('2i5,4e12.4')
326 inpt1 = header_line.read(lines[1])
327
328 params = np.asarray(inpt1[2:], 'float64')
329 DrParams={'drtype':drtype, 'params':params, 'ref':ref}
330 else:
331 DrParams = None
332 print((' for ion %5s unknown DR type = %5i' %(ions, drtype)))
333 else:
334 DrParams=None
335 return DrParams
336
337
338
340 '''
341 read a chianti excitation-autoionization file and return the EA ionization rate data
342 derived from splupsRead
343 {"lvl1":lvl1,"lvl2":lvl2,"ttype":ttype,"gf":gf,"de":de,"cups":cups,"bsplups":bsplups,"ref":ref}
344 '''
345 if filename:
346 splupsname = filename
347 else:
348 zion = util.convertName(ions)
349 if zion['Z'] < zion['Ion']:
350 print(' this is a bare nucleus that has no ionization rate')
351 return
352
353 fname = util.ion2filename(ions)
354 splupsname = fname+'.easplups'
355 if not os.path.exists(splupsname):
356 print((' could not find file: ', splupsname))
357 return {"lvl1":-1}
358
359 else:
360 input=open(splupsname,'r')
361 s1=input.readlines()
362 input.close()
363 nsplups=0
364 ndata=2
365 while ndata > 1:
366 s1a=s1[nsplups][:]
367 s2=s1a.split()
368 ndata=len(s2)
369 nsplups=nsplups+1
370 nsplups=nsplups-1
371 lvl1=[0]*nsplups
372 lvl2=[0]*nsplups
373 ttype=[0]*nsplups
374 gf=[0.]*nsplups
375 de=[0.]*nsplups
376 cups=[0.]*nsplups
377 nspl=[0]*nsplups
378 splups=np.zeros((nsplups,9),'Float64')
379
380
381
382 header_line1 = FortranRecordReader('6x,3i3,8e10.0')
383 header_line2 = FortranRecordReader('6x,3i3,12e10.0')
384
385 for i in range(0,nsplups):
386 try:
387 inpt = header_line1.read(s1[i])
388
389 except:
390 inpt = header_line2.read(s1[i])
391
392 lvl1[i]=inpt[0]
393 lvl2[i]=inpt[1]
394 ttype[i]=inpt[2]
395 gf[i]=inpt[3]
396 de[i]=inpt[4]
397 cups[i]=inpt[5]
398 if len(inpt) > 13:
399 nspl[i]=9
400 splups[i].put(list(range(9)),inpt[6:])
401 else:
402 nspl[i]=5
403 splups[i].put(list(range(5)),inpt[6:])
404
405 ref=[]
406 for i in range(nsplups+1,len(s1)):
407 s1a=s1[i][:-1]
408 ref.append(s1a.strip())
409
410
411 return {"lvl1":lvl1,"lvl2":lvl2,"ttype":ttype,"gf":gf,"de":de,"cups":cups
412 ,"nspl":nspl,"splups":splups,"ref":ref}
413
414
415
417 """
418 read chianti splom files and returns
419 {"lvl1":lvl1,"lvl2":lvl2,"deryd":de,"gf":gf,"eryd":eout,"omega":omout}
420 currently only works for 5 point spline fit files
421 splomRead probably does just as good a job - this function may be redundant
422 """
423
424
425 if filename:
426 input = open(filename)
427 else:
428 fname=ion2filename(ions)
429 omname=fname+extension
430 input=open(omname,'r')
431 lines=input.readlines()
432 input.close()
433 format=FortranFormat('5i3,8e10.3')
434 data=5
435 iline=0
436 lvl1=[]
437 lvl2=[]
438 ttype=[]
439 gf=[]
440 de=[]
441 om=[]
442 z=1
443 while z > 0:
444 omdat1=FortranLine(lines[iline],format)
445 z=omdat1[0]
446 if z > 0:
447 l1=omdat1[2]
448 l2=omdat1[3]
449 ttype1=omdat1[4]
450 gf1=omdat1[5]
451 de1=omdat1[6]
452 btf1=omdat1[7]
453 om1=omdat1[8:]
454
455 lvl1.append(l1)
456 lvl2.append(l2)
457 ttype.append(ttype1)
458 gf.append(gf1)
459 de.append(de1)
460 om.append(om1)
461 iline=iline+1
462 omout=np.asarray(om,'Float64')
463 ref=lines[iline:-1]
464
465
466
467
468
469
470
471 Splom={"lvl1":lvl1,"lvl2":lvl2,'ttype':ttype,"gf":gf, "deryd":de,"omega":omout, 'ref':ref}
472 return Splom
473
474
475
476 -def elvlcRead(ions, filename=0, getExtended=0, verbose=0, useTh=1):
477 """
478 reads the new format elvlc files
479 read a chianti energy level file that has 6 energy columns
480 and returns
481 {"lvl":lvl,"conf":conf,"term":term,"spin":spin,"l":l,"spd":spd,"j":j
482 ,"mult":mult,"ecm":ecm,"eryd":eryd,"ecmth":ecmth,"erydth":erydth,
483 "ecmx":ecmx,"erydx":erydx,"ref":ref,"pretty":pretty, 'ionS':ions}
484 if a energy value for ecm or eryd is zero(=unknown), the theoretical values
485 (ecmth and erydth) are inserted if useTh is true
486 Python 2.7 and 3 compliant
487 """
488
489
490 '%7i%30s%5s%5i%5s%5.1f%15.3f%15.3f \n'
491
492 header_line = FortranRecordReader('i7,a30,a5,i5,a5,f5.1,2f15.3')
493
494
495
496 if filename:
497 elvlname = filename
498 bname = os.path.basename(filename)
499 ions = bname.split('.')[0]
500 else:
501 fname = util.ion2filename(ions)
502 elvlname=fname+'.elvlc'
503 if not os.path.isfile(elvlname):
504 print((' elvlc file does not exist: %s'%(elvlname)))
505 return {'status':0}
506 status = 1
507 input=open(elvlname,'r')
508 s1=input.readlines()
509 input.close()
510 nlvls=0
511 ndata=2
512 while ndata > 1:
513 s1a=s1[nlvls][:-1]
514 s2=s1a.split()
515 ndata=len(s2)
516 nlvls=nlvls+1
517 nlvls-=1
518 if verbose:
519 print((' nlvls = %i'%(nlvls)))
520 lvl=[0]*nlvls
521 conf = [0]*nlvls
522 term=[0]*nlvls
523 label = [0]*nlvls
524 spin=[0]*nlvls
525 spd=[0]*nlvls
526 l = ['']*nlvls
527 j = [0.]*nlvls
528 mult = [0.]*nlvls
529 ecm=[0]*nlvls
530 ecmth=[0]*nlvls
531 pretty=[0]*nlvls
532 if getExtended:
533 extended = [' ']*nlvls
534 for i in range(0,nlvls):
535 if verbose:
536 print((s1[i][0:115]))
537
538 inpt =header_line.read(s1[i][0:115])
539 lvl[i]=inpt[0]
540 term[i]=inpt[1].strip()
541 label[i] = inpt[2]
542 spin[i]=inpt[3]
543 spd[i]=inpt[4].strip()
544 l[i] = const.Spd.index(spd[i])
545 j[i]=inpt[5]
546 mult[i] = 2.*inpt[5] + 1.
547 ecm[i]=inpt[6]
548 ecmth[i]=inpt[7]
549 if ecm[i] == 0.:
550 if useTh:
551 ecm[i] = ecmth[i]
552 stuff = term[i].strip() + ' %1i%1s%3.1f'%( spin[i], spd[i], j[i])
553 pretty[i] = stuff.strip()
554 if getExtended:
555 cnt = s1[i].count(',')
556 if cnt > 0:
557 idx = s1[i].index(',')
558 extended[i] = s1[i][idx+1:]
559 eryd = [x*const.invCm2ryd for x in ecm]
560 erydth = [x*const.invCm2ryd for x in ecmth]
561 ref=[]
562
563 for i in range(nlvls+1,len(s1)):
564 s1a=s1[i]
565 ref.append(s1a.strip())
566
567
568 info = {"lvl":lvl,"conf":conf, "term":term,'label':label, "spin":spin, "spd":spd, "l":l, "j":j,
569 'mult':mult, "ecm":ecm, 'eryd':eryd,'erydth':erydth, "ecmth":ecmth, "ref":ref,
570 "pretty":pretty, 'status':status, 'filename':elvlname}
571 if getExtended:
572 info['extended'] = extended
573 return info
574
575
576
577 -def elvlcWrite(info, outfile=0, addLvl=0, includeRyd=0):
578 '''
579 for files created after elvlc format change in November 2012
580 creates a .elvlc in the current directory
581 info is a dictionary that must contain the following keys
582 ionS, the Chianti style name of the ion such as c_4
583 term, a string showing the configuration
584 spin, an integer of the spin of the state in LS coupling
585 l, an integer of the angular momentum quantum number
586 spd, an string for the alphabetic symbol of the angular momemtum, S, P, D, etc
587 j, a floating point number, the total angular momentum
588 ecm, the observed energy in inverse cm, if unknown, the value is 0.
589 eryd, the observed energy in Rydbergs, if unknown, the value is 0.
590 ecmth, the calculated energy from the scattering calculation, in inverse cm
591 erydth, the calculated energy from the scattering calculation in Rydbergs
592 ref, the references in the literature to the data in the input info
593
594 the output filename will be ionS+'.elvlc' unless outfile is specified
595 addLvl is to add a constant value to the index of all levels
596 setting includeRyd will also write the Rydberg energies in the extended area, demarked by a comma
597 '''
598 if outfile:
599 elvlcName = outfile
600 else:
601 try:
602 gname = info['ionS']
603 except:
604 print(' ''ionS'' not included in input dict')
605 return
606 elvlcName = gname + '.elvlc'
607 print((' elvlc file name = ', elvlcName))
608
609
610
611
612
613 if 'label' not in info:
614 nlvl = len(info['ecm'])
615 info['label'] = [' ']*nlvl
616 if 'eryd' not in info:
617 info['eryd'] = [x*const.invCm2ryd for x in info['ecm']]
618 if 'erydth 'not in info:
619 info['erydth'] = [x*const.invCm2ryd for x in info['ecmth']]
620
621 out = open(elvlcName, 'w')
622 for i, aterm in enumerate(info['term']):
623 thisTerm = aterm.ljust(29)
624 thisLabel = info['label'][i].ljust(4)
625
626 if includeRyd:
627 pstring = '%7i%30s%5s%5i%5s%5.1f%15.3f%15.3f , %15.8f , %15.8f \n'%(i+1+addLvl, thisTerm, thisLabel, info['spin'][i], info['spd'][i],info['j'][i], info['ecm'][i], info['ecmth'][i], info['eryd'][i], info['erydth'][i])
628 else:
629 pstring = '%7i%30s%5s%5i%5s%5.1f%15.3f%15.3f \n'%(i+1+addLvl, thisTerm, thisLabel, info['spin'][i], info['spd'][i],info['j'][i], info['ecm'][i], info['ecmth'][i])
630 out.write(pstring)
631 out.write(' -1\n')
632 out.write('%filename: ' + os.path.split(elvlcName)[1] + '\n')
633
634
635
636
637
638 for aref in info['ref']:
639 out.write(aref + '\n')
640 out.write(' -1 \n')
641 out.close()
642 return
643
644
645
647 """
648 read a chianti energy level file and returns
649 {"lvl":lvl,"conf":conf,"term":term,"spin":spin,"l":l,"spd":spd,"j":j
650 ,"mult":mult,"ecm":ecm,"eryd":eryd,"ref":ref}
651 """
652
653 fstring='i5,a20,2i5,a3,i5,2f20.3'
654
655 header_line = FortranRecordReader(fstring)
656
657 if os.path.exists(filename):
658 input=open(filename,'r')
659 s1=input.readlines()
660 input.close()
661 nlvls=0
662 ndata=2
663 while ndata > 1:
664 s1a=s1[nlvls][:-1]
665 s2=s1a.split()
666 ndata=len(s2)
667 nlvls=nlvls+1
668 nlvls-=1
669 if verbose:
670 print((' nlvls = %5i'%(nlvls)))
671 lvl=[0]*nlvls
672 conf=[0]*nlvls
673 pqn=[0]*nlvls
674 l=[0]*nlvls
675 spd=[0]*nlvls
676 mult=[0]*nlvls
677 ecm=[0]*nlvls
678 ecmth=[0]*nlvls
679 for i in range(0,nlvls):
680 if verbose:
681 print((s1[i]))
682
683 inpt = header_line.read(s1[i])
684 lvl[i]=inpt[0]
685 conf[i]=inpt[1].strip()
686 pqn[i]=inpt[2]
687 l[i]=inpt[3]
688 spd[i]=inpt[4].strip()
689 mult[i]=inpt[5]
690 if inpt[6] == 0.:
691 ecm[i]=inpt[7]
692 else:
693 ecm[i]=inpt[6]
694 ecmth[i]=inpt[7]
695 ref=[]
696 for i in range(nlvls+1,len(s1)-1):
697 s1a=s1[i][:-1]
698 ref.append(s1a.strip())
699 return {"lvl":lvl,"conf":conf,'pqn':pqn,"l":l,"spd":spd,"mult":mult,
700 "ecm":ecm,'ecmth':ecmth, 'ref':ref}
701 else:
702 return {'errorMessage':' fblvl file does not exist'}
703
704
705
707 '''
708 to read the free-free gaunt factors of Sutherland, 1998, MNRAS, 300, 321.
709 this function reads the file and reverses the values of g2 and u
710 '''
711 xuvtop = os.environ['XUVTOP']
712 fileName = os.path.join(xuvtop, 'continuum','gffgu.dat' )
713 input = open(fileName)
714 lines = input.readlines()
715 input.close()
716
717
718 ngamma=41
719 nu=81
720 nvalues = ngamma*nu
721 g2 = np.zeros(ngamma, 'float64')
722 g21d = np.zeros(nvalues, 'float64')
723 u = np.zeros(nu, 'float64')
724 u1d = np.zeros(nvalues, 'float64')
725 gff = np.zeros((ngamma, nu), 'float64')
726 gff1d = np.zeros(nvalues, 'float64')
727
728 iline = 5
729 ivalue = 0
730 for ig2 in range(ngamma):
731 for iu in range(nu):
732 values = lines[iline].split()
733 u[iu] = float(values[1])
734 u1d[ivalue] = float(values[1])
735 g2[ig2] = float(values[0])
736 g21d[ivalue] = float(values[0])
737 gff[ig2, iu] = float(values[2])
738 gff1d[ivalue] = float(values[2])
739 iline+=1
740 ivalue += 1
741
742 return {'g2':g2, 'g21d':g21d, 'u':u, 'u1d':u1d, 'gff':gff, 'gff1d':gff1d}
743
744
745
747 '''
748 to read the integrated free-free gaunt factors of Sutherland, 1998, MNRAS, 300, 321.
749 '''
750 xuvtop = os.environ['XUVTOP']
751 fileName = os.path.join(xuvtop, 'continuum','gffint.dat' )
752 input = open(fileName)
753 lines = input.readlines()
754 input.close()
755
756 ngamma=41
757 g2 = np.zeros(ngamma, 'float64')
758 gffint = np.zeros(ngamma, 'float64')
759 s1 = np.zeros(ngamma, 'float64')
760 s2 = np.zeros(ngamma, 'float64')
761 s3 = np.zeros(ngamma, 'float64')
762
763 ivalue = 0
764 start = 4
765 for iline in range(start,start+ngamma):
766 values = lines[iline].split()
767 g2[ivalue] = float(values[0])
768 gffint[ivalue] = float(values[1])
769 s1[ivalue] = float(values[2])
770 s2[ivalue] = float(values[3])
771 s3[ivalue] = float(values[4])
772 ivalue += 1
773
774 return {'g2':g2, 'gffint':gffint, 's1':s1, 's2':s2, 's3':s3}
775
776
777
779 '''
780 to read in the free-free gaunt factors of Itoh et al. (ApJS 128, 125, 2000)
781 '''
782 xuvtop = os.environ['XUVTOP']
783 itohName = os.path.join(xuvtop, 'continuum', 'itoh.dat')
784 input = open(itohName)
785 lines = input.readlines()
786 input.close()
787 gff = np.zeros((30, 121), 'float64')
788 for iline in range(30):
789 gff[iline]= np.asarray(lines[iline].split(), 'float64')
790 return {'itohCoef':gff}
791
792
793
794
796 '''
797 to read CHIANTI files file containing the free-bound gaunt factors for n=1-6 from Karzas and Latter, 1961, ApJSS, 6, 167
798 returns {pe, klgfb}, the photon energy and the free-bound gaunt factors
799 '''
800 xuvtop = os.environ['XUVTOP']
801 fname = os.path.join(xuvtop, 'continuum', 'klgfb.dat')
802 input = open(fname)
803 lines = input.readlines()
804 input.close()
805
806 ngfb = int(lines[0].split()[0])
807 nume = int(lines[0].split()[1])
808
809 gfb = np.zeros((ngfb, ngfb, nume), 'float64')
810 nlines = len(lines)
811
812 pe = np.asarray(lines[1].split(), 'float64')
813 for iline in range(2, nlines):
814 data = lines[iline].split()
815 n = int(data[0])
816 l = int(data[1])
817 gfb[n-1, l] = np.array(data[2:], 'float64')
818 return {'pe':pe, 'klgfb':gfb}
819
820
821
823 """
824 reads an ioneq file and stores temperatures and ionization
825 equilibrium values in self.IoneqTemperature and self.Ioneq and returns
826 a dictionary containing these value and the reference to the literature
827 """
828 dir=os.environ["XUVTOP"]
829 ioneqdir = os.path.join(dir,'ioneq')
830 if ioneqname == '':
831
832 fname1 = chgui.gui.chpicker(ioneqdir,filter='*.ioneq',label = 'Select an Ionization Equilibrium file')
833 fname = os.path.join(ioneqdir, fname1)
834 if fname == None:
835 print(' no ioneq file selected')
836 return False
837 else:
838 ioneqfilename=os.path.basename(fname)
839 ioneqname,ext=os.path.splitext(ioneqfilename)
840 else:
841 filelist = util.listFiles(ioneqdir)
842 fname=os.path.join(dir,'ioneq',ioneqname+'.ioneq')
843 newlist = fnmatch.filter(filelist, '*.ioneq')
844 baselist = []
845 for one in newlist:
846 baselist.append(os.path.basename(one))
847 cnt = baselist.count(ioneqname+'.ioneq')
848 if cnt == 0:
849 print((' ioneq file not found: ', fname))
850 print(' the following files do exist: ')
851 for one in newlist:
852 print((os.path.basename(one)))
853 return
854 elif cnt == 1:
855 idx = baselist.index(ioneqname+'.ioneq')
856 if verbose:
857 print((' file exists: ', newlist[idx]))
858 fname = newlist[idx]
859 elif cnt > 1:
860 print((' found more than one ioneq file', fname))
861 return
862
863 input=open(fname,'r')
864 s1=input.readlines()
865 input.close()
866 ntemp,nele=s1[0].split()
867 if verbose:
868 print((' ntemp, nele = %5i %5i'%(ntemp, nele)))
869 nTemperature=int(ntemp)
870 nElement=int(nele)
871
872
873 header_linet = FortranRecordReader(str(nTemperature)+'f6.2')
874
875 ioneqTemperature = header_linet.read(s1[1])
876 ioneqTemperature = np.asarray(ioneqTemperature[:],'Float64')
877 ioneqTemperature = 10.**ioneqTemperature
878 nlines=0
879 idx=-1
880 while idx < 0:
881 aline=s1[nlines][0:5]
882 idx=aline.find('-1')
883 nlines+=1
884 nlines -= 1
885
886
887
888 header_lineq = FortranRecordReader('2i3,'+str(nTemperature)+'e10.2')
889
890 ioneqAll=np.zeros((nElement,nElement+1,nTemperature),'Float64')
891 for iline in range(2,nlines):
892
893 out = header_lineq.read(s1[iline])
894 iz=out[0]
895 ion=out[1]
896 ioneqAll[iz-1,ion-1].put(list(range(nTemperature)),np.asarray(out[2:],'Float64'))
897 ioneqRef = []
898 for one in s1[nlines+1:]:
899 ioneqRef.append(one[:-1])
900 del s1
901 return {'ioneqname':ioneqname,'ioneqAll':ioneqAll,'ioneqTemperature':ioneqTemperature,'ioneqRef':ioneqRef}
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
954 """
955 reads the ionization potential file, returns ip array in eV
956 """
957 topdir=os.environ["XUVTOP"]
958 ipname=os.path.join(topdir, 'ip','chianti.ip')
959 ipfile=open(ipname)
960 data=ipfile.readlines()
961 ipfile.close()
962 nip=0
963 ndata=2
964 maxz=0
965 while ndata > 1:
966 s1=data[nip]
967 s2=s1.split()
968 ndata=len(s2)
969 nip=nip+1
970 if int(s2[0]) > maxz:
971 maxz=int(s2[0])
972 if verbose:
973 print((' maxz = %5i'%(maxz)))
974 nip=nip-1
975 ip=np.zeros((maxz, maxz), 'Float64')
976 for aline in data[0:nip]:
977 s2=aline.split()
978 iz=int(s2[0])
979 ion=int(s2[1])
980 ip[iz-1, ion-1]=float(s2[2])
981 return ip*const.invCm2Ev
982
983
984
986 """
987 read a chianti masterlist file and return a list of ions
988 """
989 dir=os.environ["XUVTOP"]
990 fname=os.path.join(dir,'masterlist','masterlist.ions')
991 input=open(fname,'r')
992 s1=input.readlines()
993 input.close()
994 masterlist=[]
995 for i in range(0,len(s1)):
996 s1a=s1[i][:-1]
997 s2=s1a.split(';')
998 masterlist.append(s2[0].strip())
999 return masterlist
1000
1001
1002
1004 """
1005 returns information about ions in masterlist
1006 the reason for this file is to speed up multi-ion spectral calculations
1007 the information is stored in a pickled file 'masterlist_ions.pkl'
1008 if the file is not found, one will be created and the following information
1009 returned for each ion
1010 wmin, wmax : the minimum and maximum wavelengths in the wgfa file
1011 tmin, tmax : the minimum and maximum temperatures for which the ionization balance is nonzero
1012 for pickle compatibility between Python 2 and 3, have used floats
1013 """
1014 dir=os.environ["XUVTOP"]
1015 infoPath = os.path.join(dir, 'masterlist')
1016 infoName=os.path.join(dir,'masterlist','masterlist_ions.pkl')
1017
1018
1019 makeNew = force == 1 or not os.path.isfile(infoName)
1020
1021 if not makeNew:
1022
1023 pfile = open(infoName, 'rb')
1024 masterListInfo = pickle.load(pfile)
1025 pfile.close
1026 elif os.access(infoPath, os.W_OK):
1027
1028 defaults = defaultsRead()
1029 print((' defaults = %s'%(str(defaults))))
1030 ioneqName = defaults['ioneqfile']
1031 ioneq = ioneqRead(ioneqname = ioneqName)
1032 masterList = masterListRead()
1033 masterListInfo = {}
1034 haveZ = [0]*31
1035 haveStage = np.zeros((31, 31), 'Int32')
1036 haveDielectronic = np.zeros((31, 31), 'Int32')
1037 for one in masterList:
1038 if verbose:
1039 print((' ion = %s'%(one)))
1040 ionInfo = convertName(one)
1041 z = ionInfo['Z']
1042 stage = ionInfo['Ion']
1043 haveZ[z] = 1
1044 dielectronic = ionInfo['Dielectronic']
1045 if dielectronic:
1046 haveDielectronic[z, stage] = 1
1047 else:
1048 haveStage[z, stage] = 1
1049 thisIoneq = ioneq['ioneqAll'][z- 1, stage - 1 + dielectronic]
1050 good = thisIoneq > 0.
1051 goodTemp = ioneq['ioneqTemperature'][good]
1052 tmin = float(goodTemp.min())
1053 tmax = float(goodTemp.max())
1054 vgood = thisIoneq == thisIoneq.max()
1055 vgoodTemp = float(ioneq['ioneqTemperature'][vgood][0])
1056 wgfa = wgfaRead(one)
1057 nZeros = wgfa['wvl'].count(0.)
1058
1059 while nZeros > 0:
1060 wgfa['wvl'].remove(0.)
1061 nZeros = wgfa['wvl'].count(0.)
1062
1063 wvl = np.abs(np.asarray(wgfa['wvl'], 'float64'))
1064 wmin = float(wvl.min())
1065 wmax = float(wvl.max())
1066 masterListInfo[one] = {'wmin':wmin, 'wmax':wmax, 'tmin':tmin, 'tmax':tmax, 'tIoneqMax':vgoodTemp}
1067 masterListInfo['haveZ'] = haveZ
1068 masterListInfo['haveStage'] = haveStage
1069 masterListInfo['haveDielectronic'] = haveDielectronic
1070
1071
1072 for iz in range(1, 31):
1073 ions = zion2name(iz, iz+1)
1074 thisIoneq = ioneq['ioneqAll'][iz-1, iz]
1075 good = thisIoneq > 0.
1076 goodTemp = ioneq['ioneqTemperature'][good]
1077 tmin = float(goodTemp.min())
1078 tmax = float(goodTemp.max())
1079 wmin=0.
1080 wmax = 1.e+30
1081 masterListInfo[ions] = {'wmin':wmin, 'wmax':wmax, 'tmin':tmin, 'tmax':tmax}
1082 pfile = open(infoName, 'wb')
1083 pickle.dump(masterListInfo, pfile)
1084 pfile.close
1085 else:
1086
1087
1088 masterListInfo = {}
1089 for one in masterList:
1090 ionInfo = convertName(one)
1091 z = ionInfo['Z']
1092 stage = ionInfo['Ion']
1093 dielectronic = ionInfo['Dielectronic']
1094 wmin=0.
1095 wmax = 1.e+30
1096 masterListInfo[one] = {'wmin':wmin, 'wmax':wmax, 'tmin':1.e+4, 'tmax':1.e+9}
1097
1098
1099 for iz in range(1, 31):
1100 ions = zion2name(iz, iz+1)
1101 wmin=0.
1102 wmax = 1.e+30
1103 masterListInfo[ions] = {'wmin':wmin, 'wmax':wmax, 'tmin':1.e+4, 'tmax':1.e+9}
1104 pfile = open(infoName, 'wb')
1105 pickle.dump(masterListInfo, pfile)
1106 pfile.close
1107 masterListInfo = {'noInfo':'none'}
1108 return masterListInfo
1109
1110
1111
1113 """
1114 read chianti photoionization .photox files and return
1115 {"energy", "cross"} where energy is in Rydbergs and the
1116 cross section is in cm^2
1117 the photox files are not in any released version of the CHIANTI database
1118 """
1119
1120 zion = util.convertName(ions)
1121 if zion['Z'] < zion['Ion']:
1122 print((' this is a bare nucleus that has no ionization rate'))
1123 return
1124
1125 fname = util.ion2filename(ions)
1126 paramname=fname+'.photox'
1127 input=open(paramname,'r')
1128 lines = input.readlines()
1129 input.close
1130
1131
1132 dataEnd = 0
1133 lvl1 = []
1134 lvl2 = []
1135 energy = []
1136 cross = []
1137 icounter = 0
1138 while not dataEnd:
1139 lvl11 = int(lines[icounter][:8])
1140 lvl21 = int(lines[icounter][8:15])
1141 ener = lines[icounter][15:].split()
1142 energy1 = np.asarray(ener, 'float64')
1143
1144 icounter += 1
1145 irsl = int(lines[icounter][:8])
1146 ind0 = int(lines[icounter][8:15])
1147 if irsl != lvl11 or ind0 != lvl21:
1148
1149 print((' lvl1, lvl2 = %7i %7i'%(lvl11, lvl21)))
1150 print((' irsl, indo = %7i %7i'%(irsl, ind0)))
1151 return
1152 crs = lines[icounter][15:].split()
1153 cross1 = np.asarray(crs, 'float64')
1154 lvl1.append(lvl11)
1155 lvl2.append(lvl21)
1156 energy.append(energy1)
1157 cross.append(cross1)
1158 icounter += 1
1159 dataEnd = lines[icounter].count('-1')
1160 ref = lines[icounter+1:-1]
1161 cross = np.asarray(cross, 'float64')
1162 energy = np.asarray(energy, 'float64')
1163 return {'lvl1':lvl1, 'lvl2':lvl2,'energy':energy, 'cross':cross, 'ref':ref}
1164
1165
1166
1168 """
1169 read chianti radiative recombination .rrparams files and return
1170 {'rrtype','params','ref'}
1171 """
1172
1173
1174 fname = util.ion2filename(ions)
1175 paramname=fname+'.rrparams'
1176 if os.path.isfile(paramname):
1177 input=open(paramname,'r')
1178
1179 lines=input.readlines()
1180 input.close()
1181 rrtype=int(lines[0])
1182 ref=lines[3:-2]
1183
1184 if rrtype == 1:
1185
1186
1187 header_line = FortranRecordReader('3i5,e12.4,f10.5,2e12.4')
1188
1189 params = header_line.read(lines[1])
1190 RrParams={'rrtype':rrtype, 'params':params, 'ref':ref}
1191 elif rrtype == 2:
1192
1193
1194 header_line = FortranRecordReader('3i5,e12.4,f10.5,2e11.4,f10.5,e12.4')
1195
1196 params = header_line.read(lines[1])
1197 RrParams={'rrtype':rrtype, 'params':params, 'ref':ref}
1198 elif rrtype == 3:
1199
1200
1201 header_line = FortranRecordReader('2i5,2e12.4')
1202
1203 params = header_line.read(lines[1])
1204 RrParams={'rrtype':rrtype, 'params':params, 'ref':ref}
1205 else:
1206 RrParams=None
1207 print((' for ion %5s unknown RR type = %5i' %(ions, rrtype)))
1208 return RrParams
1209 else:
1210 return {'rrtype':-1}
1211
1212
1213
1215 '''
1216 to read the new format ~ version 8 scups file containing the Burgess and Tully scaled temperature and upsilons.
1217 Python 2.7/3 compliant
1218 '''
1219
1220 if filename:
1221 scupsFileName = filename
1222 bname = os.path.basename(scupsFileName)
1223 ions = bname.split('.')[0]
1224 else:
1225 fname = util.ion2filename(ions)
1226 scupsFileName = fname+'.scups'
1227 if not os.path.isfile(scupsFileName):
1228 print((' elvlc file does not exist: %s'%(scupsFileName)))
1229 return {'status':0}
1230
1231
1232 if os.path.isfile(scupsFileName):
1233 if verbose:
1234 print(' scupsFileName = %s'%(scupsFileName))
1235 inpt = open(scupsFileName)
1236 lines = inpt.readlines()
1237 inpt.close()
1238 else:
1239 print(('file does not exist: '+str(scupsFileName)))
1240 return {'errorMessage':'file does not exist' +str(scupsFileName)}
1241 return
1242
1243
1244 minusOne = 0
1245 counter = 0
1246 while not minusOne:
1247 if '-1' in lines[counter][:4]:
1248 minusOne = 1
1249 else:
1250 counter += 1
1251 ntrans = (counter)/3
1252
1253 lvl1 = []
1254 lvl2 = []
1255 de = []
1256 gf = []
1257 lim = []
1258 ttype = []
1259 cups = []
1260 ntemp = []
1261 btemp = []
1262 bscups = []
1263 counter = 0
1264
1265
1266 for itrans in range(int(ntrans)):
1267 if verbose:
1268 print((lines[counter]))
1269 print((lines[counter+1]))
1270 print((lines[counter+2]))
1271 ll1 = lines[counter].split()
1272 lvl1.append(int(ll1[0]))
1273 lvl2.append(int(ll1[1]))
1274 de.append(float(ll1[2]))
1275 gf.append(float(ll1[3]))
1276 lim.append(float(ll1[4]))
1277 ntemp.append(int(ll1[5]))
1278 ttype.append(int(ll1[6]))
1279 cups.append(float(ll1[7]))
1280 ll2 = lines[counter+1].split()
1281 ll3 = lines[counter+2].split()
1282
1283
1284 btemp.append(np.asarray(ll2, 'float64'))
1285 bscups.append(np.asarray(ll3, 'float64'))
1286 counter += 3
1287 counter += 1
1288 ref = []
1289 for aline in lines[counter:]:
1290 ref.append(aline.strip('\n'))
1291 return {'lvl1':lvl1, 'lvl2':lvl2, 'de':de, 'gf':gf, 'lim':lim, 'ttype':ttype,'cups':cups,'ntemp':ntemp, 'btemp':btemp, 'bscups':bscups, 'ntrans':ntrans, 'ref':ref}
1292
1293
1294
1296 """
1297 read chianti .splom files and return
1298 {"lvl1":lvl1,"lvl2":lvl2,"ttype":ttype,"gf":gf,"deryd":de,"c":c,"splom":splomout,"ref":hdr}
1299 still needed for ionization cross sections
1300 """
1301
1302 if type(filename) == type(None):
1303 fname = util.ion2filename(ions)
1304 if ea:
1305 splomname=fname+'.easplom'
1306 else:
1307 splomname=fname+'.splom'
1308 else:
1309 splomname = filename
1310 input=open(splomname,'r')
1311
1312 line1=input.readline()
1313
1314 remainder=line1[16:]
1315 nom=remainder.split(' ')
1316
1317 header_line = FortranRecordReader('5i3,'+str(len(nom))+'E10.2')
1318
1319 input.seek(0)
1320 lines=input.readlines()
1321 data=5
1322 iline=0
1323 lvl1=[]
1324 lvl2=[]
1325 ttype=[]
1326 gf=[]
1327 de=[]
1328 f=[]
1329 splom=[]
1330
1331 while data > 1:
1332
1333 splomdat = header_line.read(lines[iline])
1334 l1=splomdat[2]
1335 l2=splomdat[3]
1336 tt1=splomdat[4]
1337 gf1=splomdat[5]
1338 de1=splomdat[6]
1339 f1=splomdat[7]
1340 splom1=splomdat[8:]
1341 lvl1.append(int(l1))
1342 lvl2.append(int(l2))
1343 ttype.append(int(tt1))
1344 gf.append(float(gf1))
1345 de.append(float(de1))
1346 f.append(float(f1))
1347 splom.append(splom1)
1348 iline=iline+1
1349 data=len(lines[iline].split(' ',2))
1350 hdr=lines[iline+1:-1]
1351 de=np.asarray(de,'Float64')
1352 splomout=np.asarray(splom,'Float64')
1353 splomout=np.transpose(splomout)
1354 input.close()
1355
1356 splom={"lvl1":lvl1,"lvl2":lvl2,"ttype":ttype,"gf":gf,"deryd":de,"c":f
1357 ,"splom":splomout,"ref":hdr}
1358 return splom
1359
1360
1361
1362 -def splupsRead(ions, filename=0, prot=0, ci=0, diel=0):
1363 """
1364 read a chianti splups file and return
1365 {"lvl1":lvl1,"lvl2":lvl2,"ttype":ttype,"gf":gf,"de":de,"cups":cups,"bsplups":bsplups,"ref":ref}
1366 if prot >0, then reads the psplups file
1367 if ci > 0, then reads cisplups file
1368 if diel > 0, then reads dielsplups file
1369 """
1370
1371 if filename:
1372 splupsname = filename
1373 else:
1374 fname = util.ion2filename(ions)
1375 if prot:
1376 splupsname=fname+'.psplups'
1377 elif ci:
1378 splupsname=fname+'.cisplups'
1379 elif diel:
1380 splupsname=fname+'.splups'
1381 else:
1382 splupsname=fname+'.splups'
1383 if not os.path.exists(splupsname):
1384 if prot:
1385 return {'file not found':splupsname}
1386 elif ci:
1387 return {'file not found':splupsname}
1388 else:
1389 return {'file not found':splupsname}
1390
1391 else:
1392 input=open(splupsname,'r')
1393 s1=input.readlines()
1394 input.close()
1395 nsplups=0
1396 ndata=2
1397 while ndata > 1:
1398 s1a=s1[nsplups][:]
1399 s2=s1a.split()
1400 ndata=len(s2)
1401 nsplups=nsplups+1
1402 nsplups=nsplups-1
1403 lvl1=[0]*nsplups
1404 lvl2=[0]*nsplups
1405 ttype=[0]*nsplups
1406 gf = np.zeros(nsplups, 'float64')
1407 de= np.zeros(nsplups, 'float64')
1408 cups= np.zeros(nsplups, 'float64')
1409 nspl=[0]*nsplups
1410
1411 splups = [0.]*nsplups
1412 if prot:
1413
1414
1415 header_line = FortranRecordReader('3i3,3e10.3')
1416 else:
1417
1418
1419 header_line = FortranRecordReader('6x,3i3,3e10.3')
1420 for i in range(0,nsplups):
1421
1422 inpt = header_line.read(s1[i])
1423 lvl1[i]=inpt[0]
1424 lvl2[i]=inpt[1]
1425 ttype[i]=inpt[2]
1426 gf[i]=inpt[3]
1427 de[i]=inpt[4]
1428 cups[i]=inpt[5]
1429 if prot:
1430 as1 = s1[i][39:].rstrip()
1431 else:
1432 as1 = s1[i][45:].rstrip()
1433 nspl[i] = len(as1)//10
1434
1435
1436 header_line3 = FortranRecordReader(str(nspl[i])+'e10.3' )
1437
1438 inpt = header_line3.read(as1)
1439 spl1 = np.asarray(inpt[:], 'float64')
1440 splups[i] = spl1
1441
1442 ref=[]
1443 for i in range(nsplups+1,len(s1)):
1444 s1a=s1[i][:-1]
1445 ref.append(s1a.strip())
1446 if prot:
1447
1448
1449
1450 return {"lvl1":lvl1,"lvl2":lvl2,"ttype":ttype,"gf":gf,"de":de,"cups":cups
1451 ,"nspl":nspl,"splups":splups,"ref":ref, 'filename':splupsname}
1452 else:
1453
1454
1455 return {"lvl1":lvl1,"lvl2":lvl2,"ttype":ttype,"gf":gf,"de":de,"cups":cups
1456 ,"nspl":nspl,"splups":splups,"ref":ref, 'filename':splupsname}
1457
1458
1459
1461 '''
1462 to read the two-photon A values and distribution function for the H seq
1463 '''
1464 xuvtop = os.environ['XUVTOP']
1465 fName = os.path.join(xuvtop, 'continuum', 'hseq_2photon.dat')
1466 dFile = open(fName, 'r')
1467 a = dFile.readline()
1468 y0 = np.asarray(a.split())
1469 a = dFile.readline()
1470 z0 = np.asarray(a.split())
1471 nz = 30
1472 avalue = np.zeros(nz, 'float64')
1473 asum = np.zeros(nz, 'float64')
1474 psi0 = np.zeros((nz, 17), 'float64')
1475 for iz in range(nz):
1476 a=dFile.readline().split()
1477 avalue[iz] = float(a[1])
1478 asum[iz] = float(a[2])
1479 psi = np.asarray(a[3:])
1480 psi0[iz] = psi
1481 dFile.close()
1482 return {'y0':y0, 'z0':z0, 'avalue':avalue, 'asum':asum, 'psi0':psi0.reshape(30, 17)}
1483
1484
1485
1487 '''
1488 to read the two-photon A values and distribution function for the He seq
1489 '''
1490 xuvtop = os.environ['XUVTOP']
1491 fName = os.path.join(xuvtop, 'continuum', 'heseq_2photon.dat')
1492 dFile = open(fName, 'r')
1493 a = dFile.readline()
1494 y0 = np.asarray(a.split())
1495 nz = 30
1496 avalue = np.zeros(nz, 'float64')
1497 psi0 = np.zeros((nz, 41), 'float64')
1498 for iz in range(1, nz):
1499 a=dFile.readline().split()
1500 avalue[iz] = float(a[1])
1501 psi = np.asarray(a[2:])
1502 psi0[iz] = psi
1503 dFile.close()
1504 return {'y0':y0, 'avalue':avalue, 'psi0':psi0.reshape(30, 41)}
1505
1506
1507
1509 '''
1510 Reads the Verner & Yakovlev (A&AS 109, 125, 1995) photoionization cross-section data
1511 '''
1512 xuvtop = os.environ['XUVTOP']
1513 fname = os.path.join(xuvtop, 'continuum', 'verner_short.txt')
1514 input = open(fname)
1515 lines = input.readlines()
1516 input.close()
1517
1518 nlines=465
1519 maxZ = 30+1
1520 maxNel = 30 +1
1521
1522
1523
1524 pqn = np.zeros((maxZ,maxNel),'int32')
1525 l = np.zeros((maxZ,maxNel),'int32')
1526 eth = np.zeros((maxZ,maxNel),'float64')
1527 e0 = np.zeros((maxZ,maxNel),'float64')
1528 sig0 = np.zeros((maxZ,maxNel),'float64')
1529 ya = np.zeros((maxZ,maxNel),'float64')
1530 p = np.zeros((maxZ,maxNel),'float64')
1531 yw = np.zeros((maxZ,maxNel),'float64')
1532
1533 fstring='i2,i3,i2,i2,6f11.3'
1534 header_line = FortranRecordReader(fstring)
1535
1536
1537 for iline in range(nlines):
1538
1539 out = header_line.read(lines[iline])
1540 z = out[0]
1541 nel = out[1]
1542 stage = z - nel + 1
1543 pqn[z,stage] = out[2]
1544 l[z,stage] = out[3]
1545 eth[z,stage] = out[4]
1546 e0[z,stage] = out[5]
1547 sig0[z,stage] = out[6]
1548 ya[z,stage] = out[7]
1549 p[z,stage] = out[8]
1550 yw[z,stage] = out[9]
1551
1552 return {'pqn':pqn, 'l':l, 'eth':eth, 'e0':e0, 'sig0':sig0, 'ya':ya, 'p':p, 'yw':yw}
1553
1554
1555
1557 """
1558 read the version number of the CHIANTI database
1559 """
1560 xuvtop = os.environ['XUVTOP']
1561 vFileName = os.path.join(xuvtop, 'VERSION')
1562 vFile = open(vFileName)
1563 versionStr = vFile.readline()
1564 vFile.close()
1565 return versionStr.strip()
1566
1567
1568
1569 -def wgfaRead(ions, filename=0, elvlcname=-1, total=0, verbose=0):
1570 """
1571 reads chianti wgfa file and returns
1572 {"lvl1":lvl1,"lvl2":lvl2,"wvl":wvl,"gf":gf,"avalue":avalue,"ref":ref}
1573 if elvlcname is specified, the lsj term labels are returned as 'pretty1' and 'pretty2'
1574 """
1575
1576 if filename:
1577 wgfaname = filename
1578 if elvlcname < 0:
1579 elvlcname = 0
1580 elvlc = 0
1581 elif not elvlcname:
1582 elvlcname = os.path.splitext(wgfaname)[0] + '.elvlc'
1583 if os.path.isfile(elvlcname):
1584 elvlc = elvlcRead('', elvlcname)
1585 else:
1586 elvlc = 0
1587 else:
1588 elvlc = elvlcRead('',elvlcname)
1589
1590 else:
1591 fname = util.ion2filename(ions)
1592 wgfaname=fname+'.wgfa'
1593 elvlcname = fname + '.elvlc'
1594 if os.path.isfile(elvlcname):
1595 elvlc = elvlcRead('', elvlcname)
1596 else:
1597 elvlc = 0
1598 if verbose:
1599 if elvlc:
1600 print(' have elvlc data')
1601 else:
1602 print(' do not have elvlc data')
1603
1604 input=open(wgfaname,'r')
1605 s1=input.readlines()
1606 input.close()
1607 nwvl=0
1608 ndata=2
1609 while ndata > 1:
1610 s1a=s1[nwvl]
1611 s2=s1a.split()
1612 ndata=len(s2)
1613 nwvl += 1
1614 nwvl -= 1
1615 if verbose:
1616 print((' nwvl = %10i ndata = %4i'%(nwvl, ndata)))
1617 lvl1=[0]*nwvl
1618 lvl2=[0]*nwvl
1619 wvl=[0.]*nwvl
1620 gf=[0.]*nwvl
1621 avalue=[0.]*nwvl
1622 if elvlc:
1623 pretty1 = ['']*nwvl
1624 pretty2 = ['']*nwvl
1625
1626 if verbose:
1627 print((' nwvl = %10i'%(nwvl)))
1628
1629 wgfaFormat='(2i5,f15.3,2e15.3)'
1630 header_line = FortranRecordReader(wgfaFormat)
1631 for ivl in range(nwvl):
1632 if verbose:
1633 print(' index %5i %s'%(ivl, s1[ivl]))
1634
1635 inpt = header_line.read(s1[ivl])
1636 lvl1[ivl]=inpt[0]
1637 lvl2[ivl]=inpt[1]
1638 wvl[ivl]=inpt[2]
1639 gf[ivl]=inpt[3]
1640 avalue[ivl]=inpt[4]
1641 if elvlc:
1642 idx1 = elvlc['lvl'].index(inpt[0])
1643 idx2 = elvlc['lvl'].index(inpt[1])
1644 pretty1[ivl] = elvlc['pretty'][idx1]
1645 pretty2[ivl] = elvlc['pretty'][idx2]
1646
1647 ref=[]
1648
1649 for i in range(nwvl+1,len(s1)):
1650 s1a=s1[i]
1651 ref.append(s1a.strip())
1652 Wgfa={"lvl1":lvl1,"lvl2":lvl2,"wvl":wvl,"gf":gf,"avalue":avalue,"ref":ref, 'ionS':ions, 'filename':wgfaname}
1653 if total:
1654 avalueLvl = [0.]*max(lvl2)
1655 for iwvl in range(nwvl):
1656 avalueLvl[lvl2[iwvl] -1] += avalue[iwvl]
1657 Wgfa['avalueLvl'] = avalueLvl
1658
1659 if elvlc:
1660 Wgfa['pretty1'] = pretty1
1661 Wgfa['pretty2'] = pretty2
1662
1663 return Wgfa
1664
1665
1666
1667 -def wgfaWrite(info, outfile = 0, minBranch = 0.):
1668 '''
1669 to write a wgfa file
1670 info is a dictionary the contains the following elements
1671 ionS, the Chianti style name of the ion such as c_4 for C IV
1672 lvl1 - the lower level, the ground level is 1
1673 lvl2 - the upper level
1674 wvl - the wavelength in Angstroms
1675 gf - the weighted oscillator strength
1676 avalue - the A value
1677 pretty1 - descriptive text of the lower level (optional)
1678 pretty2 - descriptive text of the upper level (optiona)
1679 ref - reference text, a list of strings
1680 minBranch: the transition must have a branching ratio greater than the specified to be written to the file
1681 '''
1682
1683
1684 if outfile:
1685 wgfaname = outfile
1686 else:
1687 print(' output filename not specified, no file will be created')
1688 return
1689
1690 print((' wgfa file name = ', wgfaname))
1691 if minBranch > 0.:
1692 info['ref'].append(' minimum branching ratio = %10.2e'%(minBranch))
1693 out = open(wgfaname, 'w')
1694
1695 nlvl = max(info['lvl2'])
1696 totalAvalue = np.zeros(nlvl, 'float64')
1697 if 'pretty1' in info:
1698 pformat = '%5i%5i%15.4f%15.3e%15.3e%30s - %30s'
1699 else:
1700 pformat = '%5i%5i%15.4f%15.3e%15.3e'
1701 for itrans, avalue in enumerate(info['avalue']):
1702
1703 if abs(info['lvl1'][itrans]) > 0 and info['lvl2'][itrans] > 0:
1704 totalAvalue[info['lvl2'][itrans] -1] += avalue
1705
1706 for itrans, avalue in enumerate(info['avalue']):
1707 if avalue > 0.:
1708 branch = avalue/totalAvalue[info['lvl2'][itrans] -1]
1709 else:
1710 branch = 0.
1711 if branch > minBranch and abs(info['lvl1'][itrans]) > 0 and info['lvl2'][itrans] > 0:
1712 if 'pretty1' in info:
1713
1714 if 'transType' in info:
1715 if info['transType'][itrans] != '':
1716 lbl2 = info['pretty2']+' ' + info['transType'][itrans]
1717 else:
1718 lbl2= info['pretty2'][itrans]
1719 pstring= pformat%(info['lvl1'][itrans], info['lvl2'][itrans], info['wvl'][itrans], info['gf'][itrans], avalue, info['pretty1'][itrans].rjust(30), lbl2.ljust(30))
1720 out.write(pstring+'\n')
1721 else:
1722 pstring= pformat%(info['lvl1'][itrans], info['lvl2'][itrans], info['wvl'][itrans], info['gf'][itrans], avalue)
1723 out.write(pstring+'\n')
1724 out.write(' -1\n')
1725 out.write('%filename: ' + wgfaname + '\n')
1726 for one in info['ref']:
1727 out.write(one+'\n')
1728 out.write(' -1 \n')
1729 out.close()
1730