add Kucucz Molecules

ISO.h

@@ -64,7 +64,7 @@ int Init(Molecule &m, Param &param, char (*qFilename)[160]){
 	}
 	fscanf (pFile, "%s", m.mName);
 	fgets(sp, 4, pFile);
-	//can be "Number of Isotopes =" or "Number of Isotopes ="
+	//can be "Number of Isotopes =" or "Number of Isotopologues ="
 	//scan until "="
 	fscanf(pFile, "%[^=]s", sp);
 	if(strcmp(sp, "Number of Isotopes ") != 0 && strcmp(sp, "Number of Isotopologues ") != 0){
@@ -135,6 +135,7 @@ int Init(Molecule &m, Param &param, char (*qFilename)[160]){
 		}
 		if(param.dataBase == 2){
 			//ExoMol
+			//Kurucz Molecules
 			if(m.nFiles > 1){
 				sprintf(m.dataFilename[i], "%s%s__%05d-%05d.", param.path, m.mName, m.fileLimit[i], m.fileLimit[i + 1]);
 			}

KuruczMolecules.py

@@ -0,0 +1,191 @@
+import numpy as np
+import math
+import struct
+import os
+import argparse
+
+# This script calculates the molecular line intensities
+# It writes data files to prepare the HELIOS-K binary files
+
+
+# Date: January 2022
+# Author: Simon Grimm
+
+
+#choose if the file contains wavenumber or not
+Wavenumber = 2			#1: Wavenumber, 2: vacuum wavelength, 3: air based wavelenght
+
+
+def main(name, mass, pfname, printA):
+	filename= '%s.asc' % name
+
+	outname = "%s.bin" % name
+
+	print(name, outname, mass)
+
+
+	output_file = open(outname,"wb")
+
+	numax = 0.0
+	nl = 0	
+
+	if(printA == 1):
+		Afile = open("%s_A.dat" % name, "w")
+
+	with open(filename) as f:
+		line = f.readlines()
+
+
+		for ii in range(len(line)):
+			l = line[ii]
+
+			#E in cm^-1
+			#atomic line list format
+			if(Wavenumber == 1):
+				wn = float(l[0:10])
+				wl = 1.0E7/wn		#wavelenght in nm
+			if(Wavenumber == 2):
+				wl = float(l[0:10])	#wavelenght in nm
+				wn = 1.0E7/wl
+			if(Wavenumber == 3):
+				wlAir = float(l[0:10])	#wavelenght in nm
+
+
+			loggf = float(l[10:17])
+			JLow = float(l[17:22])
+			ELow = float(l[22:32])
+			JUP = float(l[32:37])
+			EUP = float(l[37:48])
+			code = l[48:70]
+
+
+			e = 4.80320425E-10      #electron charge in cgs units [statcoulomb = cm^(3/2) g^(1/2) s^-1]
+			c = 2.99792458E10       #Speed of light cm/s
+			me = 9.1093835611E-28   #mass of electron in g
+			NA = 6.0221412927e23	#Avogadro Constant  1/mol
+
+
+			ELow = abs(ELow)
+			EUP = abs(EUP)
+
+
+			#somethimes ELOW is larger than EUP
+
+			if(ELow > EUP):
+				t = EUP
+				EUP = ELow
+				ELow = t
+
+				t = JUP
+				JUP = JLow
+				JLow = t
+
+			#convert air wavelength to vacuum wavelength
+			#http://www.astro.uu.se/valdwiki/Air-to-vacuum%20conversion
+
+			if(Wavenumber == 3):
+				if(wlAir > 200):
+					wlAir = wlAir * 10  #convert nm to Angstrom
+					s = 10000.0 / wlAir
+					n = 1.0 + 0.00008336624212083 + 0.02408926869968 / (130.1065924522 - s*s) + 0.0001599740894897 / (38.92568793293 - s*s)
+					wl = wlAir * n
+					wl = wl * 0.1  #convert Angstrom to nm
+				else:
+					wl = wlAir
+				wn = 1.0E7/wl		#wavelenght in nm
+
+
+			gUP = 2 * JUP + 1
+			gLow = 2 * JLow + 1
+
+			A = 8.0 * math.pi * wn * wn * (10.0**loggf) / gUP * math.pi * e * e / (me * c)
+
+			nl = nl + 1
+
+			numax = max(numax, wn)
+
+			S = math.pi * e * e * 10.0**loggf * NA / (c * c * me * mass) 
+
+			A = 8.0 * math.pi * wn * wn * 10.0**loggf / gUP * math.pi * e * e / (me * c)
+
+			if(printA == 1):
+				print(i, wn, A, ELow, gUP, mass, file = Afile)	
+			#print(wn, 1.0E7/wn, loggf, ELow, EUP, JLow, JUP, mass, A)
+
+			#print(wn, S, A, ELow, EUP, gLow, gUP)
+
+			s = struct.pack('d', wn)
+			output_file.write(s)
+			s = struct.pack('d', S)
+			output_file.write(s)
+			s = struct.pack('d', ELow)
+			output_file.write(s)
+			s = struct.pack('d', A)
+			output_file.write(s)
+
+
+
+	print(" Lines:",nl, end='')
+
+	output_file.close()
+	if(printA == 1):
+		Afile.close()
+
+
+	if(nl > 0):
+		f = open("%s.param" % name,'w')
+
+		print("Database = 2", file = f)
+		print("Molecule number = 1", file = f)
+		print("Name = %s" % name, file = f)
+		print("Number of Isotopologues = 1", file = f)
+		print("#Id Abundance      Q(296K)   g     Molar Mass(g)  partition file :", file = f)
+		print("0 1.0             0.0       0      %s        %s" % (mass, pfname), file = f)
+		print("Number of columns in partition File = 2", file = f)
+		print("Number of line/transition files = 1", file = f)
+		print("Number of lines per file :", file = f)
+		print("%d" % nl, file = f)
+		print("Line file limits :", file = f)
+		print("0", file = f)
+		print("%d" % (int(numax)+1), file = f)
+		print("#ExoMol :", file = f)
+		print("Number of states = 0", file = f)
+		print("Number of columns in transition files = 0", file = f)
+		print("Default value of Lorentzian half-width for all lines = 0.07", file = f)
+		print("Default value of temperature exponent for all lines = 0.5", file = f)
+		print("Version = %s" % filename, file = f)
+
+		f.close()
+
+
+
+if __name__ == '__main__':
+
+	parser = argparse.ArgumentParser()
+
+	parser.add_argument('-name', '--name', type=str,
+		#help='name', default = 'coax')
+		help='name', default = 'tiototo')
+
+	parser.add_argument('-mass', '--mass', type=float,
+		#help='Isotopologue mass g/mol', default = 28.0101)
+		help='Isotopologue mass g/mol', default = 61.947545)
+
+	parser.add_argument('-pfFile', '--pfFile', type=str,
+		#help='partition function file name', default = '12C-16O__Li2015.pf')
+		help='partition function file name', default = '46Ti-16O__Toto.pf')
+
+	parser.add_argument('-printA', '--printA', type=int,
+		help='print A to file', default = 0)
+
+
+	args = parser.parse_args()
+	name = args.name
+	mass = args.mass
+	pfname = args.pfFile
+	printA = args.printA
+
+	print("name: %s, mass: %g; pfFile: %s; printA: %d" % (name, mass, pfname, printA))
+
+	main(name, mass, pfname, printA)
+

heliosk.cu

@@ -846,23 +846,24 @@ printf("%g %g %g %g\n", param.numax, param.numin, param.dnu, (param.numax - para
 		int rbvs = 0;
 
 		if(param.dataBase == 2){
-			//Exomol		nu, S, El, A 
+			//Exomol		nu, S, El, A
+			//Kurucz Molecules
 			readBufferN = 4;
 		}
 		if(param.dataBase == 20){
 			//Exomol super lines	nu, S
 			readBufferN = 2;
 		}
 		if(param.dataBase == 30){
-			//Kurucz		nu, S, El, A, Gamma_nat
+			//Kurucz Atoms		nu, S, El, A, Gamma_nat
 			readBufferN = 5;
 		}
 		if(param.dataBase == 31){
-			//NIST
+			//NIST Atoms
 			readBufferN = 5;
 		}
 		if(param.dataBase == 32){
-			//VALD
+			//VALD Atoms
 			readBufferN = 5;
 		}
 

tools/checkBinary.py

@@ -1,13 +1,33 @@
 import struct
+import argparse
 
-filename = "/home/sigrimm/scratch/EXOMOL/1H2-16O__POKAZATEL__00000-00100.bin"
 
-file = open(filename, "rb")
+#filename = "/home/sigrimm/scratch/EXOMOL/1H2-16O__POKAZATEL__00000-00100.bin"
 
-for i in range(10):
-	nu = struct.unpack('d', file.read(8))[0]
-	S = struct.unpack('d', file.read(8))[0]
-	EL = struct.unpack('d', file.read(8))[0]
-	A = struct.unpack('d', file.read(8))[0]
 
-	print(nu, S, EL, A)
+def main(filename):
+	file = open(filename, "rb")
+
+	for i in range(100000000000):
+		try:
+			nu = struct.unpack('d', file.read(8))[0]
+			S = struct.unpack('d', file.read(8))[0]
+			EL = struct.unpack('d', file.read(8))[0]
+			A = struct.unpack('d', file.read(8))[0]
+
+			print(nu, S, EL, A)
+		except:
+			return
+
+
+if __name__ == '__main__':
+
+	parser = argparse.ArgumentParser()
+
+	parser.add_argument('-name', '--name', type=str,
+		help='name', default = '')
+
+	args = parser.parse_args()
+	name = args.name
+
+	main(name)