from . import absetup import os import stat import shutil import subprocess import math import logging # atomdir = absetup.catpref+'/cats/Castelli/atoms_SLrTS/' # moldir = absetup.catpref+'/cats/VALD/mol' # atomdir = os.getenv('HOME')+'/cats/VALD/atom' # atomdir = os.getenv('HOME')+'/cats/Kurucz/linelists/gf08oct17TS/' atomdir = absetup.catpref+'/cats/Kurucz/linelists/gf08oct17slr_eTS/' # atomdir = absetup.catpref+'/cats/Kurucz/linelists/gf08oct17slr_eTSa/' moldir = absetup.catpref+'/cats/Kurucz/molecules/20apr2016TS/' # moldir = absetup.catpref+'/cats/Plez/molecules/' # VALD + CN plezdir = absetup.catpref+'/cats/Plez/' # CH, TiO atomlines = [ ['vald-3000-3500a.list', 3000., 3500.], ['vald-3500-4000a.list', 3500., 4000.], ['vald-4000-4500a.list', 4000., 4500.], ['vald-4500-5000a.list', 4500., 5000.], ['vald-5000-5500a.list', 5000., 5500.], ['vald-5500-6000a.list', 5500., 6000.], ['vald-6000-6500a.list', 6000., 6500.], ['vald-6500-7000a.list', 6500., 7000.], ['vald-7000-8000a.list', 7000., 8000.], ['vald-8000-9000a.list', 8000., 9000.], ['vald-9000-10000a.list', 9000., 10000.], ['vald-10000-12000a.list', 10000., 12000.], ['vald-12000-14000a.list', 12000., 14000.], ['vald-14000-16000a.list', 14000., 16000.], ['vald-16000-18000a.list', 16000., 18000.], ['vald-18000-20000a.list', 18000., 20000.], ['vald-20000-25000a.list', 20000., 25000.] ] mollines = [ ['vald-3000-3500m.list', 3000., 3500.], ['vald-3500-4000m.list', 3500., 4000.], ['vald-4000-4500m.list', 4000., 4500.], ['vald-4500-5000m.list', 4500., 5000.], ['vald-5000-5500m.list', 5000., 5500.], ['vald-5500-6000m.list', 5500., 6000.], ['vald-6000-6500m.list', 6000., 6500.], ['vald-6500-7000m.list', 6500., 7000.], ['vald-7000-8000m.list', 7000., 8000.], ['vald-8000-9000m.list', 8000., 9000.], ['vald-9000-10000m.list', 9000., 10000.], ['vald-10000-12000m.list', 10000., 12000.], ['vald-12000-14000m.list', 12000., 14000.], ['vald-14000-16000m.list', 14000., 16000.], ['vald-16000-18000m.list', 16000., 18000.], ['vald-18000-20000m.list', 18000., 20000.], ['vald-20000-25000m.list', 20000., 25000.] ] TiOlines = [ ['TiOVALD-bsyn-3000-4000.list', 3000., 4000.], ['TiOVALD-bsyn-4000-4200.list', 4000., 4200.], ['TiOVALD-bsyn-4200-4400.list', 4200., 4400.], ['TiOVALD-bsyn-4400-4600.list', 4400., 4600.], ['TiOVALD-bsyn-4600-4800.list', 4600., 4800.], ['TiOVALD-bsyn-4800-5000.list', 4800., 5000.], ['TiOVALD-bsyn-5000-5200.list', 5000., 5200.], ['TiOVALD-bsyn-5200-5400.list', 5200., 5400.], ['TiOVALD-bsyn-5400-5600.list', 5400., 5600.], ['TiOVALD-bsyn-5600-5800.list', 5600., 5800.], ['TiOVALD-bsyn-5800-6000.list', 5800., 6000.], ['TiOVALD-bsyn-6000-6200.list', 6000., 6200.], ['TiOVALD-bsyn-6200-6400.list', 6200., 6400.], ['TiOVALD-bsyn-6400-6600.list', 6400., 6600.], ['TiOVALD-bsyn-6600-6800.list', 6600., 6800.], ['TiOVALD-bsyn-6800-7000.list', 6800., 7000.], ['TiOVALD-bsyn-7000-7200.list', 7000., 7200.], ['TiOVALD-bsyn-7200-7400.list', 7200., 7400.], ['TiOVALD-bsyn-7400-7600.list', 7400., 7600.], ['TiOVALD-bsyn-7600-7800.list', 7600., 7800.], ['TiOVALD-bsyn-7800-8000.list', 7800., 8000.], ['TiOVALD-bsyn-8000-8200.list', 8000., 8200.], ['TiOVALD-bsyn-8200-8400.list', 8200., 8400.], ['TiOVALD-bsyn-8400-8600.list', 8400., 8600.], ['TiOVALD-bsyn-8600-8800.list', 8600., 8800.], ['TiOVALD-bsyn-8800-9000.list', 8800., 9000.], ['TiOVALD-bsyn-9000-9500.list', 9000., 9500.], ['TiOVALD-bsyn-9500-10000.list', 9500., 10000.], ['TiOVALD-bsyn-10000-10500.list', 10000., 10500.], ['TiOVALD-bsyn-10500-11000.list', 10500., 11000.], ['TiOVALD-bsyn-11000-11500.list', 11000., 11500.], ['TiOVALD-bsyn-11500-12000.list', 11500., 12000.], ['TiOVALD-bsyn-12000-12500.list', 12000., 12500.], ['TiOVALD-bsyn-12500-13000.list', 12500., 13000.], ['TiOVALD-bsyn-13000-14000.list', 13000., 14000.], ['TiOVALD-bsyn-14000-15000.list', 14000., 15000.], ['TiOVALD-bsyn-15000-16000.list', 15000., 16000.], ['TiOVALD-bsyn-16000-17000.list', 16000., 17000.], ['TiOVALD-bsyn-17000-18000.list', 17000., 18000.], ['TiOVALD-bsyn-18000-19000.list', 18000., 19000.], ['TiOVALD-bsyn-19000-20000.list', 19000., 20000.], ['TiOVALD-bsyn-20000-22000.list', 20000., 22000.], ['TiOVALD-bsyn-22000-25000.list', 22000., 25000.] ] # We now take the CN lines from Kurucz and include them # in the 'vald' files. # # C12N14lines = [ # ['C12N14_Brooke_bsyn-3000-4000.list', 3000., 4000.], # ['C12N14_Brooke_bsyn-4000-5000.list', 4000., 5000.], # ['C12N14_Brooke_bsyn-5000-6000.list', 5000., 6000.], # ['C12N14_Brooke_bsyn-6000-7000.list', 6000., 7000.], # ['C12N14_Brooke_bsyn-7000-8000.list', 7000., 8000.], # ['C12N14_Brooke_bsyn-8000-10000.list', 8000., 10000.], # ['C12N14_Brooke_bsyn-10000-12000.list', 10000., 12000.], # ['C12N14_Brooke_bsyn-12000-14000.list', 12000., 14000.], # ['C12N14_Brooke_bsyn-14000-20000.list', 14000., 20000.], # ['C12N14_Brooke_bsyn-20000-25000.list', 20000., 25000.] #] # C13N14lines = [ # ['C13N14_Brooke_bsyn-3000-4000.list', 3000., 4000.], # ['C13N14_Brooke_bsyn-4000-5000.list', 4000., 5000.], # ['C13N14_Brooke_bsyn-5000-6000.list', 5000., 6000.], # ['C13N14_Brooke_bsyn-6000-7000.list', 6000., 7000.], # ['C13N14_Brooke_bsyn-7000-8000.list', 7000., 8000.], # ['C13N14_Brooke_bsyn-8000-10000.list', 8000., 10000.], # ['C13N14_Brooke_bsyn-10000-12000.list', 10000., 12000.], # ['C13N14_Brooke_bsyn-12000-14000.list', 12000., 14000.], # ['C13N14_Brooke_bsyn-14000-20000.list', 14000., 20000.], # ['C13N14_Brooke_bsyn-20000-25000.list', 20000., 25000.] # ] # C12N15lines = [ # ['C12N15_Brooke_bsyn-3000-4000.list', 3000., 4000.], # ['C12N15_Brooke_bsyn-4000-5000.list', 4000., 5000.], # ['C12N15_Brooke_bsyn-5000-6000.list', 5000., 6000.], # ['C12N15_Brooke_bsyn-6000-7000.list', 6000., 7000.], # ['C12N15_Brooke_bsyn-7000-8000.list', 7000., 8000.], # ['C12N15_Brooke_bsyn-8000-10000.list', 8000., 10000.], # ['C12N15_Brooke_bsyn-10000-12000.list', 10000., 12000.], # ['C12N15_Brooke_bsyn-12000-14000.list', 12000., 14000.], # ['C12N15_Brooke_bsyn-14000-20000.list', 14000., 20000.], # ['C12N15_Brooke_bsyn-20000-25000.list', 20000., 25000.] # ] def turboinit(lam1, lam2, step, initdir='.'): # Initialize for a set of Turbospectrum synthetic atmosphere calculations. # Turbospectrum does not require much initialization, we just need # to write lam1, lam2, and step to a file which will then be read # by 'turbospec'. We can also copy the appropriate line list(s) here. for j in range(len(atomlines)): if (lam1 >= atomlines[j][1] and lam1 <= atomlines[j][2]): j1a = j if (lam2 >= atomlines[j][1] and lam2 <= atomlines[j][2]): j2a = j for j in range(len(mollines)): if (lam1 >= mollines[j][1] and lam1 <= mollines[j][2]): j1m = j if (lam2 >= mollines[j][1] and lam2 <= mollines[j][2]): j2m = j # for j in range(len(C12N14lines)): # if (lam1 >= C12N14lines[j][1] and lam1 <= C12N14lines[j][2]): # j1C12N14 = j # # if (lam2 >= C12N14lines[j][1] and lam2 <= C12N14lines[j][2]): # j2C12N14 = j # # for j in range(len(C13N14lines)): # if (lam1 >= C13N14lines[j][1] and lam1 <= C13N14lines[j][2]): # j1C13N14 = j # # if (lam2 >= C13N14lines[j][1] and lam2 <= C13N14lines[j][2]): # j2C13N14 = j # # for j in range(len(C12N15lines)): # if (lam1 >= C12N15lines[j][1] and lam1 <= C12N15lines[j][2]): # j1C12N15 = j # # if (lam2 >= C12N15lines[j][1] and lam2 <= C12N15lines[j][2]): # j2C12N15 = j for j in range(len(TiOlines)): if (lam1 >= TiOlines[j][1] and lam1 <= TiOlines[j][2]): j1TiO = j if (lam2 >= TiOlines[j][1] and lam2 <= TiOlines[j][2]): j2TiO = j finit = open(initdir+'/tspec.init','w') finit.write("lam1= %0.3f lam2= %0.3f step= %0.4e\n" % (lam1, lam2, step)) finit.write("natom= %d\n" % (j2a-j1a+1)) for j in range(j1a,j2a+1): os.system("\cp -f %s/%s %s" % (atomdir, atomlines[j][0], initdir)) finit.write("%s\n" % (atomlines[j][0])) finit.write("nmol= %d\n" % (j2m-j1m+1)) for j in range(j1m,j2m+1): os.system("\cp -f %s/%s %s" % (moldir, mollines[j][0], initdir)) finit.write("%s\n" % (mollines[j][0])) # finit.write("nC12N14= %d\n" % (j2C12N14-j1C12N14+1)) # for j in range(j1C12N14,j2C12N14+1): # os.system("\cp -f %s/%s %s" % (plezdir, C12N14lines[j][0], initdir)) # finit.write("%s\n" % (C12N14lines[j][0])) # # finit.write("nC13N14= %d\n" % (j2C13N14-j1C13N14+1)) # for j in range(j1C13N14,j2C13N14+1): # os.system("\cp -f %s/%s %s" % (plezdir, C13N14lines[j][0], initdir)) # finit.write("%s\n" % (C13N14lines[j][0])) # # finit.write("nC12N15= %d\n" % (j2C12N15-j1C12N15+1)) # for j in range(j1C12N15,j2C12N15+1): # os.system("\cp -f %s/%s %s" % (plezdir, C12N15lines[j][0], initdir)) # finit.write("%s\n" % (C12N15lines[j][0])) finit.write("nTiO= %d\n" % (5*(j2TiO-j1TiO+1))) for iTiO in [46,47,48,49,50]: for j in range(j1TiO,j2TiO+1): os.system("\cp -f %s/TiOVALD/%d%s %s" % (plezdir, iTiO, TiOlines[j][0], initdir)) finit.write("%d%s\n" % (iTiO, TiOlines[j][0])) os.system("\cp -f %s/12CH-Masseron_bsyn.list %s" % (plezdir, initdir)) os.system("\cp -f %s/13CH-Masseron_bsyn.list %s" % (plezdir, initdir)) finit.close() def turbospec(m, atmname, spname, vturb=2.0, vrot=0.01, atoms=[], abun=[], initdir='.', workdir='.', wait=True): cwd = os.getcwd() wdir = os.path.abspath(workdir) idir = os.path.abspath(initdir) fa = open(atmname,'r') l = fa.readline() # Read the first line of the atmosphere file and try # to figure out which kind of model it is fa.close() interpol = ('INTERPOL' in l) # Interpolated MARCS model? atlas = ('TEFF' in l) and ('GRAVITY' in l) # ATLAS model? os.system("\cp -f %s %s/%s.mod" % (atmname, workdir, atmname)) shutil.rmtree('%s/DATA' % (workdir), ignore_errors=True) os.system("mkdir %s/DATA" % (workdir)) os.system("\cp -f %s/DATA/Hlinedata %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/atomicweights.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/jonabs_vac_v07.3.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/LISTE_molecules_all_v12.1.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/IRWIN_molecules_v1?.1.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/IRWIN_atoms_v07.3.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/tableU2.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/tableU3.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/CIA.H2-H2.Yi.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/CIA.H2-He.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/CIA.H2-H.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/CIA.H-He.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/Anstee-1802-2014.dat %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/stehle.dat.SUN %s/DATA\n" % (plezdir, workdir)) os.system("\cp -f %s/DATA/vcsbalmer.dat %s/DATA\n" % (plezdir, workdir)) fi = open(idir+'/tspec.init','r') l = fi.readline() ll = l.split() lam1, lam2, step = float(ll[1]), float(ll[3]), float(ll[5]) l = fi.readline() ll = l.split() natom = int(ll[1]) atomfile = [] for iatom in range(natom): l = fi.readline() atomfile.append(l[:-1]) l = fi.readline() ll = l.split() nmol = int(ll[1]) molfile = [] for imol in range(nmol): l = fi.readline() molfile.append(l[:-1]) # l = fi.readline() # ll = l.split() # nC12N14 = int(ll[1]) # C12N14file = [] # for iC12N14 in range(nC12N14): # l = fi.readline() # C12N14file.append(l[:-1]) # # l = fi.readline() # ll = l.split() # nC13N14 = int(ll[1]) # C13N14file = [] # for iC13N14 in range(nC13N14): # l = fi.readline() # C13N14file.append(l[:-1]) # # l = fi.readline() # ll = l.split() # nC12N15 = int(ll[1]) # C12N15file = [] # for iC12N15 in range(nC12N15): # l = fi.readline() # C12N15file.append(l[:-1]) l = fi.readline() ll = l.split() nTiO = int(ll[1]) TiOfile = [] for iTiO in range(nTiO): l = fi.readline() TiOfile.append(l[:-1]) fi.close() # Find the H fraction (and number of elements to print) # Note - Turbospectrum normalises abundances to H, not total - so # we just need to figure out the Solar H abundance fhe = 0.07834 # Corresponds to Y=0.248; Grevesse & Sauval 1998 fz0 = 0 nelem = 0 for aa in absetup.stdabun: ielem = aa[1] iabun = aa[2] fz0 = fz0 + 10**iabun if (ielem < 93): nelem += 1 fh0 = 1 - fhe - fz0 # fh = 0.92156 # Hydrogen abundance (by number fraction) logfh0 = math.log10(fh0) ft = open(workdir+'/turbospec.com','w') # ft.write("#!/bin/csh -f\n") ft.write("#!/bin/bash\n") ft.write("cd %s\n" % (wdir)) # ft.write("set lam_min = '%0.0f'\n" % (lam1)) # ft.write("set lam_max = '%0.0f'\n" % (lam2)) # ft.write("set deltalam = '%0.3f'\n\n" % (step)) # ft.write("set METALLIC = '%10.3f'\n" % (m)) # ft.write("set TURBVEL = '%0.1f'\n\n" % (vturb)) ft.write(absetup.binpath+"babsma_lu << EOF\n") ft.write("'LAMBDA_MIN:' '%0.0f'\n" % (lam1)) ft.write("'LAMBDA_MAX:' '%0.0f'\n" % (lam2)) ft.write("'LAMBDA_STEP:' '%0.3f'\n" % (step)) ft.write("'MODELINPUT:' '%s.mod'\n" % (atmname)) if (interpol or atlas): ft.write("'MARCS-FILE:' '.false.'\n") else: ft.write("'MARCS-FILE:' '.true.'\n") ft.write("'MODELOPAC:' '%sopac'\n" % (atmname)) ft.write("'METALLICITY:' '%0.3f'\n" % (m)) ft.write("'ALPHA/Fe :' '0.00'\n") # ft.write("'HELIUM :' '0.00'\n") # ft.write("'R-PROCESS :' '0.00'\n") # ft.write("'S-PROCESS :' '0.00'\n") # ft.write("'INDIVIDUAL ABUNDANCES:' '0'\n") # ft.write("'INDIVIDUAL ABUNDANCES:' '%d'\n" % (len(atoms))) ft.write("'INDIVIDUAL ABUNDANCES:' '%d'\n" % (nelem)) for aa in absetup.stdabun: # ida = None ida = 0. selem = aa[0] ielem = aa[1] if (ielem in atoms): ida = abun[atoms.index(ielem)] if (selem in atoms): ida = abun[atoms.index(selem)] # if (ida != None): if (aa[2] < -19.99): iabun = -99.0 else: iabun = 12 + m + aa[2] + ida - logfh0 if (ielem < 93): ft.write('%d %6.3f\n' % (ielem, iabun)) ft.write("'XIFIX:' 'T'\n") ft.write("%0.1f\n" % (vturb)) ft.write("EOF\n\n") ft.write("########################################################################\n\n") ft.write(absetup.binpath+"bsyn_lu <'''+cwd+'/'+spname+'\n') ft.write('echo \'# Turbospectrum\' > %s\n' % (cwd+'/'+spname)) ft.write('echo \'# atomdir = %s\' >> %s\n' % (atomdir, cwd+'/'+spname)) ft.write('echo \'# moldir = %s\' >> %s\n' % (moldir, cwd+'/'+spname)) ft.write('echo \'# atmname = %s\' >> %s\n' % (atmname, cwd+'/'+spname)) ft.write('''awk '{printf("%0.4f %0.4f %0.4e\\n", $1, $2, $3/3.14159265)}' resultTMP.spec >>'''+cwd+'/'+spname+'\n') ft.close() print('TURBOSPECTRUM: Computing synthetic spectrum '+spname) logging.debug('TURBOSPECTRUM: Computing synthetic spectrum %s' % spname) os.chmod(workdir+'/turbospec.com',stat.S_IXUSR | stat.S_IRUSR | stat.S_IWUSR) p = subprocess.Popen(workdir+'/turbospec.com') if (wait): os.waitpid(p.pid,0) else: return p # Notes: # SYNTHE fluxes are flux moments at the surface of the star in ergs/cm**2/s/ster/nm # Turbospectrum seems to produce total Fluxes. The difference is a factor of Pi. # So Turbospectrum output should be divided by Pi to be compatible with SYNTHE.