-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathumist2krome_custm.py
53 lines (42 loc) · 1.41 KB
/
umist2krome_custm.py
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
import sys
# read UMIST 2022 rates and convert them to KROME format
# print(sys.argv)
network = sys.argv[1]
fname_umist = network+'.rates'
fname = "network_umist.dat" # output file
skip = ["PHOTON", "CRPHOT", "CRP"]
body = "@format:idx,R,R,P,P,P,P,tmin,tmax,rate\n"
body += "@common:user_Auv,user_alb,user_xi,user_AuvAv\n"
count = 0
for row in open(fname_umist):
srow = row.strip()
if srow == "" or srow.startswith("#"):
continue
arow = srow.split(":")
rtype = arow[1]
rr = arow[2:4]
pp = arow[4:8]
ka, kb, kc = [float(x) for x in arow[9:12]]
tmin, tmax = [float(x) for x in arow[12:14]]
# print(rr, pp, ka, kb, kc)
rate = None
if rtype == "CR":
rate = "%.2e * (Tgas / 3.0e2)**(%.2f) * (1./(1.-user_alb)) * (%.2f)" % (ka, kb, kc)
elif rtype == 'CP':
rate = "%.2e " % ka
elif rtype == "PH":
rate = "%.2e * user_xi * exp(-%.2f * user_Auv / user_AuvAv)" % (ka, kc)
else:
rate = "%.2e" % ka
if kb != 0e0:
rate += " * (Tgas / 3.0e2)**(%.2f)" % kb
if kc != 0e0:
rate += " * exp(-%.2f / Tgas)" % kc
body += "%d,%s,%s,%.2e,%.2e,%s\n" % (count, ",".join(rr), ",".join(pp), tmin, tmax, rate)
count += 1
body = body.replace(",e-,", ",E,")
for s in skip:
body = body.replace(",%s," % s, ",,")
with open(fname, "w") as f:
f.write(body)
print("Wrote %d reactions to %s" % (count, fname))