crm/testsaocrm_adf04_S.py at main · AtomDB/crm · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
import saocrm, pylab, numpy, pickle
import pyatomdb
import sys


import file_readers.read_adf04 as read_adf04

Z=16 # Sulphur
z1 = int(sys.argv[1]) # ioncharge +1
elsymb=pyatomdb.atomic.Ztoelsymb(Z)
ionsymb = elsymb+"_%i"%(z1)
ionchsymb = elsymb+"$^{%i+}$"%(z1-1)


C = saocrm.Ion(Z,z1)

# ADF04 only, since all the data is in one file, pre-read the file
a04datfile = 'data/adas/%s%i_adf04.dat'%(elsymb.lower(), z1-1)
a04dat =  read_adf04.read_adf04(a04datfile)

print('reading levels')
C.read_lv(saocrm.readAdf04Lev,[a04dat])
print('reading transitions')
C.read_trn(saocrm.readAdf04Trn,[a04dat])
print('reading cross section')
C.read_ec(saocrm.readAdf04Ce,[a04dat])
#print('reading AI data')
#C.read_aut(saocrm.readFacAI,['/export1/projects/EBIT/test_Si/Ar/Ar02a.ai'])

print('reading done')
Angstroms= True


Elist = numpy.arange(2000,5001,100)
eDensity = 1e0
BROADEN=0.02 #eV

wvbins = numpy.arange(1,300, 0.1)
if Angstroms:
  ebins = 12398.425/wvbins[::-1]
#ebins=numpy.arange(2000.,10001.)
C.set_spectral_bins(wvbins, units='A')
spec = numpy.zeros([len(Elist), len(wvbins)-1])
llist = {}
dE = 0.0
fig = pylab.figure()
fig.show()
ax = fig.add_subplot(111)

#Elist=[2000,3000,4000,5000,6000,7000,8000]
Elist = [2000,2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000,5000,6000,8000,10000,14000]

Elist = [2000,3000,4000]

for iE, E in enumerate(Elist):
  C.gather_rates(E,eDensity, dE=dE)
  print('rates gathered')
  C.solve_levels(E,eDensity)
  print('levels solved')
  print(C.levpop[:20])
  C.calc_lines()
  print('lines calculated')

  llist[E]=C.return_lines()
#  spec=C.return_spectrum()
#s = numpy.append(spec[0],spec)

  specb=C.return_spectrum(BROADEN, broadenunits='A')
#  if Angstroms:
    #specb=specb[::-1]
  spec[iE,:]=specb
  sb = numpy.append(specb[0],specb)
  #if Angstroms:
    #ax.plot(wvbins,  sb, drawstyle='steps')
  #else:
  ax.plot(C.spectral_bin_edges, sb, drawstyle='steps')
  pylab.draw()
#  zzz=input()
ret = {}
ret['Elist'] = Elist
ret['ebins'] =wvbins
ret['spec'] = spec
ret['lines']=llist
pickle.dump(ret, open('spec_v1_ne1cm3_%s.pkl'%(ionsymb),'wb'))
#fig = pylab.figure()
#fig.show()
#ax = fig.add_subplot(111)
#ax.plot(C.spectral_bin_edges, s, drawstyle='steps')
#ax.plot(C.spectral_bin_edges, sb, drawstyle='steps')
#pylab.draw()
#zzz=input()