from __future__ import print_function

from . import absetup
from . import kurucz
import math
import scipy.interpolate as interpolate
import numpy
import os
import string
import subprocess
import stat
import logging
import struct
# import syntmoog

# atomdir = absetup.catpref+'/cats/Castelli/atoms_SLr/'
# atomdir = '${HOME}/cats/Castelli/atoms_modr/'
# atomdir = '${HOME}/cats/Kurucz/linelists/gfhyper100/'
# atomdir = '${HOME}/cats/Kurucz/linelists/gf18feb16mod/'
# atomdir = '${HOME}/cats/Kurucz/linelists/gf05jun16/'
# atomdir = '${HOME}/cats/Kurucz/linelists/gf18feb16modr/'
# atomdir = '${HOME}/cats/Kurucz/linelists/gf08oct17sl/'
# atomdir = absetup.catpref+'/cats/Kurucz/linelists/gf08oct17slr/'
atomdir = absetup.catpref+'/cats/Kurucz/linelists/gf08oct17slr_e/'
# atomdir = absetup.catpref+'/cats/Kurucz/linelists/gf08oct17sl_abo/'
moldir  = absetup.catpref+'/cats/Kurucz/molecules/20apr2016/'
preddir = absetup.catpref+'/cats/Castelli2015/linelists/'

atomlines = [
  ['gf0150', 1000., 1500.],
  ['gf0200', 1500., 2000.],
  ['gf0300', 2000., 3000.],
  ['gf0400', 3000., 4000.],
  ['gf0500', 4000., 5000.],
  ['gf0600', 5000., 6000.],
  ['gf0800', 6000., 8000.],
  ['gf1200', 8000., 12000.],
  ['gf3000', 12000., 30000.] 
]


def do_synbeg_sel(lam1, lam2, step, tio=True, predicted=False, nlte=False, initdir='.', elements=['All'], molecules=['All']):

# Same as do_synbeg, except that individual selected elements and molecules
# can be specified.

# The NLTE switch in principle allows different initialization procedures
# for NLTE and LTE atmospheres. The IFNLTE switch in the input to
# SYNBEG can be set to 0 or 1. However, not clear this makes any
# difference. So for now, the nlte option is ignored.

    idir = os.path.abspath(initdir)

    print('SYNBEG: Initializing line list for %0.2f < LAMBDA < %0.2f\n' % (lam1, lam2))
    logging.debug('SYNBEG: Initializing line list for %0.2f < LAMBDA < %0.2f' % (lam1, lam2) )

# Select the relevant atomic line data files
    
    for j in range(len(atomlines)):
        if (lam1 >= atomlines[j][1] and lam1 <= atomlines[j][2]):
            j1 = j
            
        if (lam2 >= atomlines[j][1] and lam2 <= atomlines[j][2]):
            j2 = j
            
# Select only the specified elements (default: All) and write to gfsel.synthe
            
    zsel = []
    for elem in [('H', 1, 0.92156), ('He', 2, 0.078438)] + absetup.stdabun:
        if elem[0] in elements or 'All' in elements:
            zsel.append(elem[1])
            print('Including %s' % elem[0])
            logging.debug('Including %s' % elem[0])
        else:
            print('Skipping %s' % elem[0])
            logging.debug('Skipping %s' % elem[0])

            
    fgf = open(idir+'/gfsel.synthe','w')
    for j in range(j1, j2+1):
        p = os.path.expandvars('%s%s.100' % (atomdir, atomlines[j][0]))
        with open(p,'r') as f:
            for l in f:
                szrd  = l[19:21]
                zrd  = int(szrd)
                if (zrd in zsel): fgf.write(l)

    fgf.close()
        

#    for j in range(j1, j2+1):
#        print atomlines[j][0]
    
    lam1nm = lam1/10. 
    lam2nm = lam2/10
    R = lam1/step
    
    if (tio):
        extension = 'tio'
    else:
        extension = 'notio'

    ifnlte = 0    
#    if (nlte):
#        ifnlte = 1
#        print('Initializing for NLTE')
#        logging.debug('Initializing for NLTE')
#        extension = extension + '-nlte'
#    else:
#        ifnlte = 0
#        print('Initializing for LTE')
#        logging.debug('Initializing for LTE')
#        extension = extension + '-lte'

# Set up the .com file to read the line data
    
    f = open(idir+'/synbeg.com','w')
#    f.write('#!/bin/csh\n')
    f.write('#!/bin/bash\n')
    f.write('cd '+idir+'\n')
    f.write('\\rm -f fort.*\n')
    f.write(absetup.binpath+'synbeg.exe<<EOF>'+idir+'/synbeg.out\n')
#    f.write('AIR       520.0     540.0     500000.   0.     0          10 .001         0   00\n')
    f.write('AIR       %-9.2f %-9.2f %-9.1f 0.     %i          10 .001         0   00\n' % (lam1nm, lam2nm, R, ifnlte))
    f.write('AIRorVAC  WLBEG     WLEND     RESOLU    TURBV  IFNLTE LINOUT CUTOFF        NREAD\n')
    f.write('EOF\n')
    
#    f.write('\\rm -f fort.11\n')
#    for j in range(j1, j2+1):
#        sav = False
#        for elem in absetup.stdabun + [('H',1,0.92156), ('He', 2, 0.078438)]:
#            if elem[0] in elements or 'All' in elements:
#                if (j==j1): print('Including %s' % elem[0])
#                ielem = elem[1]
#                f.write("awk '{if ($3>=%5.2f && $3<%5.2f) print $0}' %s%s.100 >> fort.11.TMP\n" % (ielem-0.01, ielem+0.5, atomdir, atomlines[j][0]))
#                sav = True
#        if (sav):
#            f.write("sort -n fort.11.TMP > fort.11\n")
#            f.write('rgfalllinesnew.exe>'+idir+'/%s.out\n' % (atomlines[j][0]) )
#            f.write('\\rm -f fort.11\n')

    f.write('ln -fs '+idir+'/gfsel.synthe fort.11\n')
    f.write(absetup.binpath+'rgfalllinesnew.exe>'+idir+'/rgfalllinesnew.out\n')

#    f.write('rgfalllinesnew.exe>'+idir+'/%s.out\n' % (atomlines[j][0]) )
    f.write('\\rm -f fort.11\n')

    if (predicted):
        print('Including predicted lines')
        logging.debug('Including predicted lines')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/lowlines.bin fort.11\n')
#        f.write(absetup.binpath+'rpredict.exe>'+idir+'/predictedlines.out\n')
#        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+preddir+'/fclowlines.bin fort.11\n')
        f.write(absetup.binpath+'rpredict.exe>'+idir+'/predictedlow.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+preddir+'/fchighlines.bin fort.11\n')
        f.write(absetup.binpath+'rpredict.exe>'+idir+'/predictedhigh.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping predicted lines')
        logging.debug('Skipping predicted lines')

    if ('CH' in molecules or 'All' in molecules):
        print('Including CH')
        logging.debug('Including CH')
#        f.write('ln -fs '+absetup.catpref+'/cats/Kurucz/molecules/20apr2016/chmasseron.asc fort.11\n')
        f.write('ln -fs '+moldir+'/chmasseron.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/chmasseron.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping CH')
        logging.debug('Skipping CH')

    if ('MgH' in molecules or 'All' in molecules):    
        print('Including MgH')
        logging.debug('Including MgH')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/mghax.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/mghax.out\n')
#        f.write('\\rm -f fort.11\n')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/mghbx.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/mghbx.out\n')
        f.write('ln -fs '+moldir+'/mgh15082018.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/mgh.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping MgH')
        logging.debug('Skipping MgH')

    if ('NH' in molecules or 'All' in molecules):    
        print('Including NH')
        logging.debug('Including NH')
        f.write('ln -fs '+moldir+'/nhaxCas.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/nhax.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/nhcaCas.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/nhca.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping NH')
        logging.debug('Skipping NH')
   
    if ('OH' in molecules or 'All' in molecules):    
        print('Including OH')
        logging.debug('Including OH')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/ohax.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/ohax.out\n')
#        f.write('\\rm -f fort.11\n')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/ohxx.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/ohxx.out\n')
        f.write('ln -fs '+moldir+'/ohupdate.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/oh.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping OH')
        logging.debug('Skipping OH')
  
    if ('SiH' in molecules or 'All' in molecules):
        print('Including SiH')
        logging.debug('Including SiH')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/sihax.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sihax.out\n')
#        f.write('ln -fs '+moldir+'/sihax.asc fort.11\n')
        # New SiH data from Yurchenko et al. 2018 (via Kurucz)
        f.write('ln -fs '+moldir+'/sihaxsightly12012018.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sihax.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/sihxxsightly12012018.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sihxx.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping SiH')
        logging.debug('Skipping SiH')
   
    if ('H2' in molecules or 'All' in molecules):    
        print('Including H2')
        logging.debug('Including H2')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/h2bx.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/h2bx.out\n')
#        f.write('\\rm -f fort.11\n')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/h2cx.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/h2cx.out\n')
        f.write('ln -fs '+moldir+'/h2.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/h2.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/h2xx.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/h2xx.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping H2')
        logging.debug('Skipping H2')
    
    if ('C2' in molecules or 'All' in molecules):    
        print('Including C2')
        logging.debug('Including C2')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2ax.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ax.out\n')
#        f.write('\\rm -f fort.11\n')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2ba.asc fort.11\n') 
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ba.out\n')
#        f.write('\\rm -f fort.11\n')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2da.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2da.out\n')
#        f.write('\\rm -f fort.11\n')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2ea.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ea.out\n')
        f.write('ln -fs '+moldir+'/c2ax.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ax.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/c2ba.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ba.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/c2dabrookek.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2da.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/c2ea.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ea.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping C2')
        logging.debug('Skipping C2')
        
    if ('CN' in molecules or 'All' in molecules):    
        print('Including CN')
        logging.debug('Including CN')
        f.write('ln -fs '+moldir+'/cnaxbrookek.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/cnaxbrookek.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/cnbxbrookek.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/cnbxbrookek.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/cnxx12brooke.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/cnxx12brooke.out\n')
        f.write('\\rm -f fort.11\n')         
    else:
        print('Skipping CN')
        logging.debug('Skipping CN')

    if ('CO' in molecules or 'All' in molecules):        
        print('Including CO')
        logging.debug('Including CO')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/coax.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/coax.out\n')
#        f.write('\\rm -f fort.11\n')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/coxx.asc fort.11\n')
#        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/coxx.out\n')
        f.write('ln -fs '+moldir+'/coax.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/coax.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/coxx.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/coxx.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping CO')
        logging.debug('Skipping CO')
   
    if ('SiO' in molecules or 'All' in molecules):        
        print('Including SiO')
        logging.debug('Including SiO')
        f.write('ln -fs '+moldir+'/sioax.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sioax.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/sioex.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sioex.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('ln -fs '+moldir+'/sioxx.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sioxx.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping SiO')
        logging.debug('Skipping SiO')
        
    if ('CaH' in molecules):   
        print('Including CaH')
        logging.debug('Including CaH')
        f.write('ln -fs '+moldir+'/cah.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/cah.out\n')
        f.write('\\rm -f fort.11\n')         
    else:
        print('Skipping CaH')
        logging.debug('Skipping CaH')

    if ('FeH' in molecules):
        print('Including FeH')
        logging.debug('Including FeH')
        f.write('ln -fs '+moldir+'/fehfx.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/fehfx.out\n')
        f.write('\\rm -f fort.11\n')         
    else:
        print('Skipping FeH')
        logging.debug('Skipping FeH')

    if ('VO' in molecules):
        print('Including VO')   
        logging.debug('Including VO')   
        f.write('ln -fs '+moldir+'/voax.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/voax.out\n')
        f.write('\\rm -f fort.11\n')         
        f.write('ln -fs '+moldir+'/vobx.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/vobx.out\n')
        f.write('\\rm -f fort.11\n')         
        f.write('ln -fs '+moldir+'/vocx.asc fort.11\n')
        f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/vocx.out\n')
        f.write('\\rm -f fort.11\n')         
    else:
        print('Skipping VO')
        logging.debug('Skipping VO')
    
    if ('TiO' in molecules or 'All' in molecules):   
        print('Including TiO')    
        logging.debug('Including TiO')    
        if (tio):
            print('Including Schwenke TiO lines')
            logging.debug('Including Schwenke TiO lines')
#            f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/schwenke.bin fort.11\n')
#            f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/eschwenke.bin fort.48\n')
            f.write('ln -fs '+moldir+'/tioschwenke.bin fort.11\n')
            f.write('ln -fs '+moldir+'/eschwenke.bin fort.48\n')
            f.write(absetup.binpath+'rschwenk.exe>'+idir+'/rschwenk.out\n')
            f.write('\\rm -f fort.11\n')
            f.write('\\rm -f fort.48\n')
        else:
            print('Skipping Schwenke TiO lines')
            logging.debug('Skipping Schwenke TiO lines')
    else:
        print('Skipping TiO')
        logging.debug('Skipping TiO')
    
    if ('H2O' in molecules or 'All' in molecules):
        print('Including H2O')
        logging.debug('Including H2O')
#        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/h2ofast.bin fort.11\n')
        f.write('ln -fs '+moldir+'/h2ofastfix.bin fort.11\n')
        f.write(absetup.binpath+'rh2ofast.exe>'+idir+'/h2ofast.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping H2O')
        logging.debug('Skipping H2O')

    f.write('mv fort.20 '+idir+'/fort20.%s\n' % (extension))
    f.write('mv fort.19 '+idir+'/fort19.%s\n' % (extension))
    f.write('mv fort.93 '+idir+'/fort93.%s\n' % (extension))
    f.write('mv fort.14 '+idir+'/fort14.%s\n' % (extension))
    f.write('mv fort.12 '+idir+'/fort12.%s\n' % (extension))

    f.close()
    
#    os.system('source synbeg.com')

    os.chmod(idir+'/synbeg.com',stat.S_IXUSR | stat.S_IRUSR | stat.S_IWUSR)
    p = subprocess.Popen(idir+'/synbeg.com')
    os.waitpid(p.pid,0)

    
# def synbeg_sel(lam1, lam2, step, predicted=False, initdir='.', elements=['All'], molecules=['All']):
def synbeg(lam1, lam2, step, predicted=False, initdir='.', elements=['All'], molecules=['All']):

# Set up the input files for SYNTHE. They will be placed in the
# directory 'initdir'.
# Two sets of files are generated, with and without TiO, respectively.
# The appropriate set will then be selected by 'synthespec', depending
# on the temperature.

    do_synbeg_sel(lam1, lam2, step, tio=True, predicted=predicted, nlte=False, initdir=initdir, elements=elements, molecules=molecules)
    do_synbeg_sel(lam1, lam2, step, tio=False, predicted=predicted, nlte=False, initdir=initdir, elements=elements, molecules=molecules)
    
    
    
def do_synbeg_OLD(lam1, lam2, step, tio=True, predicted=False, nlte=False, initdir='.'):

# The NLTE switch in principle allows different initialization procedures
# for NLTE and LTE atmospheres. The IFNLTE switch in the input to
# SYNBEG can be set to 0 or 1. However, not clear this makes any
# difference. So for now, the nlte option is ignored.

    idir = os.path.abspath(initdir)

    print('SYNBEG: Initializing line list for %0.2f < LAMBDA < %0.2f\n' % (lam1, lam2))
    logging.debug('SYNBEG: Initializing line list for %0.2f < LAMBDA < %0.2f' % (lam1, lam2) )

    
    for j in range(len(atomlines)):
        if (lam1 >= atomlines[j][1] and lam1 <= atomlines[j][2]):
            j1 = j
            
        if (lam2 >= atomlines[j][1] and lam2 <= atomlines[j][2]):
            j2 = j

#    for j in range(j1, j2+1):
#        print atomlines[j][0]
    
    lam1nm = lam1/10. 
    lam2nm = lam2/10
    R = lam1/step
    
    if (tio):
        extension = 'tio'
    else:
        extension = 'notio'

    ifnlte = 0    
#    if (nlte):
#        ifnlte = 1
#        print('Initializing for NLTE')
#        logging.debug('Initializing for NLTE')
#        extension = extension + '-nlte'
#    else:
#        ifnlte = 0
#        print('Initializing for LTE')
#        logging.debug('Initializing for LTE')
#        extension = extension + '-lte'
    
    f = open(idir+'/synbeg.com','w')
#    f.write('#!/bin/csh\n')
    f.write('#!/bin/bash\n')
    f.write('cd '+idir+'\n')
    f.write('\\rm -f fort.*\n')
    f.write(absetup.binpath+'synbeg.exe<<EOF>'+idir+'/synbeg.out\n')
#    f.write('AIR       520.0     540.0     500000.   0.     0          10 .001         0   00\n')
    f.write('AIR       %-9.2f %-9.2f %-9.1f 0.     %i          10 .001         0   00\n' % (lam1nm, lam2nm, R, ifnlte))
    f.write('AIRorVAC  WLBEG     WLEND     RESOLU    TURBV  IFNLTE LINOUT CUTOFF        NREAD\n')
    f.write('EOF\n')
    
    for j in range(j1, j2+1):
#        print atomlines[j][0]
# Castelli line lists:
#        f.write('ln -fs ${HOME}/cats/Castelli/atoms/%s.100 fort.11\n' % (atomlines[j][0]))
# Kurucz lists with HFS:        
#        f.write('ln -fs ${HOME}/cats/Kurucz/linelists/gfhyper100/%s.100 fort.11\n' % (atomlines[j][0]))

        f.write('ln -fs %s%s.100 fort.11\n' % (atomdir, atomlines[j][0]))
        f.write(absetup.binpath+'rgfalllinesnew.exe>'+idir+'/%s.out\n' % (atomlines[j][0]) )
        f.write('\\rm -f fort.11\n')
        
# ln -fs ${HOME}/cats/Castelli/atoms/gf1200.100 fort.11
# rgfalllinesnew.exe>gf0600.out
# \rm -f fort.11

    if (predicted):
        print('Including predicted lines')
        logging.debug('Including predicted lines')
        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/lowlines.bin fort.11\n')
        f.write(absetup.binpath+'rpredict.exe>'+idir+'/predictedlines.out\n')
        f.write('\\rm -f fort.11\n')
    else:
        print('Skipping predicted lines')
        logging.debug('Skipping predicted lines')

#    f.write('ln -fs ${HOME}/cats/Castelli/molecules/chax.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/chax.out\n')
#    f.write('\\rm -f fort.11\n')
#    f.write('ln -fs ${HOME}/cats/Castelli/molecules/chbx.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/chbx.out\n')
#    f.write('\\rm -f fort.11\n')
#    f.write('ln -fs ${HOME}/cats/Castelli/molecules/chcx.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/chcx.out\n')
#    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Kurucz/molecules/20apr2016/chmasseron.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/chmasseron.out\n')
    f.write('\\rm -f fort.11\n')
    
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/mghax.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/mghax.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/mghbx.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/mghbx.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/nhax.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/nhax.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/nhca.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/nhca.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/ohax.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/ohax.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/ohxx.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/ohxx.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/sihax.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sihax.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/h2bx.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/h2bx.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/h2cx.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/h2cx.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2ax.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ax.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2ba.asc fort.11\n') 
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ba.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2da.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2da.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/c2ea.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/c2ea.out\n')
    f.write('\\rm -f fort.11\n')
#    f.write('ln -fs ${HOME}/cats/Castelli/molecules/cnax.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/cnax.out\n')
#    f.write('\\rm -f fort.11\n')
#    f.write('ln -fs ${HOME}/cats/Castelli/molecules/cnbx.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/cnbx.out\n')
#    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Kurucz/molecules/20apr2016/cnaxbrookek.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/cnaxbrookek.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Kurucz/molecules/20apr2016/cnbxbrookek.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/cnbxbrookek.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Kurucz/molecules/20apr2016/cnxx12brooke.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/cnxx12brooke.out\n')
    f.write('\\rm -f fort.11\n')         
     
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/coax.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/coax.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/coxx.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/coxx.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/sioax.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sioax.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/sioex.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sioex.out\n')
    f.write('\\rm -f fort.11\n')
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules/sioxx.asc fort.11\n')
    f.write(absetup.binpath+'rmolecasc.exe>'+idir+'/sioxx.out\n')
    f.write('\\rm -f fort.11\n')

#    f.write('ln -fs ${HOME}/cats/Kurucz/molecules/20apr2016/cah.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/cah.out\n')
#    f.write('\\rm -f fort.11\n')         
#
#    f.write('ln -fs ${HOME}/cats/Kurucz/molecules/20apr2016/fehfx.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/fehfx.out\n')
#    f.write('\\rm -f fort.11\n')         
#
#    f.write('ln -fs ${HOME}/cats/Kurucz/molecules/20apr2016/voax.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/voax.out\n')
#    f.write('\\rm -f fort.11\n')         
#    f.write('ln -fs ${HOME}/cats/Kurucz/molecules/20apr2016/vobx.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/vobx.out\n')
#    f.write('\\rm -f fort.11\n')         
#    f.write('ln -fs ${HOME}/cats/Kurucz/molecules/20apr2016/vocx.asc fort.11\n')
#    f.write('rmolecasc.exe>'+idir+'/vocx.out\n')
#    f.write('\\rm -f fort.11\n')         

    
    if (tio):
        print('Including Schwenke TiO lines')
        logging.debug('Including Schwenke TiO lines')
        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/schwenke.bin fort.11\n')
        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/eschwenke.bin fort.48\n')
        f.write(absetup.binpath+'rschwenk.exe>'+idir+'/rschwenk.out\n')
        f.write('\\rm -f fort.11\n')
        f.write('\\rm -f fort.48\n')
    else:
        print('Skipping Schwenke TiO lines')
        logging.debug('Skipping Schwenke TiO lines')
    
    f.write('ln -fs '+absetup.catpref+'/cats/Castelli/linelists/h2ofast.bin fort.11\n')
    f.write(absetup.binpath+'rh2ofast.exe>'+idir+'/h2ofast.out\n')
    f.write('\\rm -f fort.11\n')

    f.write('mv fort.20 '+idir+'/fort20.%s\n' % (extension))
    f.write('mv fort.19 '+idir+'/fort19.%s\n' % (extension))
    f.write('mv fort.93 '+idir+'/fort93.%s\n' % (extension))
    f.write('mv fort.14 '+idir+'/fort14.%s\n' % (extension))
    f.write('mv fort.12 '+idir+'/fort12.%s\n' % (extension))

    f.close()
    
#    os.system('source synbeg.com')

    os.chmod(idir+'/synbeg.com',stat.S_IXUSR | stat.S_IRUSR | stat.S_IWUSR)
    p = subprocess.Popen(idir+'/synbeg.com')
    os.waitpid(p.pid,0)
    
    
    

def synbeg_OLD(lam1, lam2, step, predicted=False, initdir='.'):

# Set up the input files for SYNTHE. They will be placed in the
# directory 'initdir'.
# Two sets of files are generated, with and without TiO, respectively.
# The appropriate set will then be selected by 'synthespec', depending
# on the temperature.

    do_synbeg_OLD(lam1, lam2, step, tio=True, predicted=predicted, nlte=False, initdir=initdir)
    do_synbeg_OLD(lam1, lam2, step, tio=False, predicted=predicted, nlte=False, initdir=initdir)



def mkmod(m, atmname, modname, atoms=[], abun=[], vturb=None, surfaceflux=False):

# Set up a model atmosphere input file for SYNTHE.
# This entails adding a special header to a normal ATLAS9 or ATLAS12
# model atmoshere and modifying the abundances according to the 'atoms' and
# 'abun' parameters.
# 

#    fh = 0.92156
    fhe = 0.07834      # Corresponds to Y=0.248; Grevesse & Sauval 1998
#    if ('H' in atoms):  fh = abun[atoms.index('H')]
    if ('He' in atoms): fhe = abun[atoms.index('He')]
    
# First find the H fraction

    fh, mcorr = kurucz.fhcorr(m, atoms, abun)

#    fz = 0
#    for aa in absetup.stdabun:
#        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)]
# 
#        iabun = aa[2] + ida       
#        fz = fz + 10**(iabun+m)
#
#    fh = 1 - fhe - fz

# If the user insists:        

    if ('H' in atoms):  fh = abun[atoms.index('H')]

# Write the header for SYNTHE    

    fmod = open(modname,'w')
    if (surfaceflux):
        fmod.write('SURFACE FLUX\n')
    else:
        fmod.write('SURFACE INTENSI 17 1.,.9,.8,.7,.6,.5,.4,.3,.25,.2,.15,.125,.1,.075,.05,.025,.01\n')
        
    fmod.write('ITERATIONS 1 PRINT 2 PUNCH 2\n')
    fmod.write('CORRECTION OFF\n')
    fmod.write('PRESSURE OFF\n')
    fmod.write('READ MOLECULES\n')
    fmod.write('MOLECULES ON\n')

    fatm = open(atmname,'r')
    s = fatm.readline()
    atlas12 = False
    
#    for s in file(atmname):
    while (s != ''):
        ss = s.split()
        if (ss[0] == 'TEFF'): 
            teff = float(ss[1])
            slte = ss[-1]
            nlte = (slte == 'NLTE')
            
        if (ss[0] == 'TITLE'):
            if ('ATLAS12' in s): atlas12 = True    
                            
        if ('ABUNDANCE TABLE' in s):   # Skip these ATLAS12 model cards
            while not ('READ DECK' in s): s = fatm.readline()    
                         
        if ('ABUNDANCE SCALE' in s):
            mmod = math.log10(float(ss[2]))
            dm = mmod - m
            if (abs(dm) > 0.01) and not atlas12:
                print('WARNING: Using [m/H]=%0.3f atmosphere for [m/H]=%0.3f spectrum' % (mmod, m))
                logging.debug('Using [m/H]=%0.3f atmosphere for [m/H]=%0.3f spectrum' % (mmod, m))

#            fmod.write('ABUNDANCE SCALE %9.5f %s' % (10**m, string.rstrip(s[26:],'\n')))
            fmod.write('ABUNDANCE SCALE %9.5f ABUNDANCE CHANGE 1 %7.5f 2 %8.6f' % (10**(m+mcorr), fh, fhe))
           
           
            for j in range(len(absetup.stdabun)):
                if (j % 5 == 0):
                    fmod.write('\n ABUNDANCE CHANGE')
               
                ida = 0.
                ielem = absetup.stdabun[j][1]
                selem = absetup.stdabun[j][0]
               
                if (ielem in atoms):
                    ida = abun[atoms.index(ielem)]
                   
                if (selem in atoms):
                    ida = abun[atoms.index(selem)]
    
                iabun = absetup.stdabun[j][2] + ida
                fmod.write('%3d %7.3f' % (ielem, iabun))
           
            fmod.write('\n')    
        elif ('READ DECK' in s):
                    
            fmod.write(s)
            ss = s.split()
            nrd = int(ss[2])
            logging.debug('Now read %d layers:' % nrd)
            for j in range(nrd):
                s = fatm.readline()
                logging.debug(s)
#                ss = s.split()
#                print len(s)
                if len(s) > 100:
                    ss = list(struct.unpack_from('15s9s10s10s10s10s10s10s10s10s',s.encode()))
                    fstr = '%15.8E%9.1F%10.3E%10.3E%10.3E%10.3E%10.3E%10.3E%10.3E%10.3E\n'
                else:
                    ss = list(struct.unpack_from('15s9s10s10s10s10s10s10s10s',s.encode()))
                    fstr = '%15.8E%9.1F%10.3E%10.3E%10.3E%10.3E%10.3E%10.3E%10.3E\n'
                     
                if (vturb is not None):
                    ss[6] = vturb*1e5
                fmod.write(fstr % tuple([float(x) for x in ss])) 
        else:
            if not ('ABUNDANCE CHANGE' in s):
                fmod.write(s)
        logging.debug('Reading next line:')
        s = fatm.readline()  
        logging.debug(s)  
          
           
    fmod.close()
    fatm.close()
    return (teff, nlte)


def synthespec(m, atmname, spname, vturb=2.0, vrot=0.01, 
               atoms=[], abun=[], initdir='.', workdir='.', 
               wait=True, surfaceflux=False):

# Spectral synthesis using SYNTHE.
# The input files must have been initialized with 'synbeg' and should
# be present in the 'initdir' directory. They will be copied to
# 'workdir' and renamed so the original files will remain intact.

# atoms[] is a list of elements and abun[] the corresponding offsets
# in their abundances with respect to the reference (absetup.stdabun).
# atoms[] can be given either as atomic numbers or the element
# abbreviations (i.e., ['O', 'Mg'] and [8, 12] are equally valid). 
# If surfaceflux=True, set the 'SURFACE FLUX' card and do not use
# ROTATE. Otherwise, set 'SURFACE INTENSITY' to the 17 angles and use
# ROTATE to calculate fluxes. The result should be the same if there
# is no rotation.

    cwd = os.getcwd()
    wdir = os.path.abspath(workdir)
    idir = os.path.abspath(initdir)

# First set up the model atmosphere file, using output from an ATLAS9 run

    teff, nlte = mkmod(m, atmname, wdir+'/SPEC.mod', atoms=atoms, 
                       abun=abun, vturb=vturb, surfaceflux=surfaceflux)

# The TiO lines can be skipped for hotter stars, this will speed up things
# very significantly

    if (teff > 4500):
        extension = 'notio'
        print('Teff = %0.0f K for %s: Excluding TiO lines' % (teff, atmname))
        logging.debug('Teff = %0.0f K for %s: Excluding TiO lines' % (teff, atmname))
    else:
        extension = 'tio'
        print('Teff = %0.0f K for %s: Including TiO lines' % (teff, atmname))
        logging.debug('Teff = %0.0f K for %s: Including TiO lines' % (teff, atmname))
        
# Was the model atmosphere computed with NLTE?        
#    if (nlte):
#        extension = extension+'-nlte'
#        print('%s is NLTE model' % (atmname))
#        logging.debug('%s is NLTE model' % (atmname))
#    else:
#        extension = extension+'-lte'
#        print('%s is LTE model' % (atmname))
#        logging.debug('%s is LTE model' % (atmname))

        
# Then set up the COM file to run SYNTHE. This will be put in 'workdir'    
    
    fc = open(wdir+'/synthe.com','w')
#    fc.write('#!/bin/csh\n')
    fc.write('#!/bin/bash\n')
    fc.write('cd '+wdir+'\n')
    fc.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules.dat fort.2\n')
    fc.write('ln -fs '+absetup.catpref+'/cats/Castelli/continua.dat fort.17\n')
    fc.write(absetup.binpath+'xnfpelsyn.exe< "SPEC.mod" >xnfpelsyn.out\n')
    fc.write('mv fort.10 xnftSPEC.dat\n')
    fc.write('\\rm -f fort.*\n')
    fc.write('ln -fs xnftSPEC.dat fort.10\n')
    
    fc.write('cp  '+idir+'/fort20.%s fort.20\n' % (extension))
    fc.write('cp  '+idir+'/fort19.%s fort.19\n' % (extension))
    fc.write('cp  '+idir+'/fort93.%s fort.93\n' % (extension))
    fc.write('cp  '+idir+'/fort14.%s fort.14\n' % (extension))
    fc.write('cp  '+idir+'/fort12.%s fort.12\n' % (extension))

    fc.write('ln -fs '+absetup.catpref+'/cats/Castelli/he1tables.dat fort.18\n')
    fc.write(absetup.binpath+'synthe.exe>synthe.out\n')
    fc.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules.dat fort.2\n')
    fc.write('cat <<EOF >fort.25\n')
    fc.write('0.0       0.        1.        0.        0.        0.        0.        0.\n')
    fc.write('0.\n')
    fc.write('RHOXJ     R1        R101      PH1       PC1       PSI1      PRDDOP    PRDPOW\n')
    fc.write('EOF\n')
    fc.write(absetup.binpath+'spectrv.exe<SPEC.mod>spectrv.out\n')
    fc.write('mv -f fort.7 SPEC.dat\n')
    
    fc.write('\\rm -f fort.*\n')
    fc.write('ln -fs SPEC.dat fort.1\n')

    if (not surfaceflux):
        fc.write(absetup.binpath+'rotate.exe<<EOF>rotate.out\n')
        fc.write('    1\n')
        fc.write('%0.2f\n' % (vrot))
        fc.write('EOF\n')
        fc.write('mv -f ROT1 SPECvr2.dat\n')
#    fc.write('ln -fs SPECvr2.dat fort.21\n')
#    fc.write('broaden.exe<<EOF>broaden.out\n')
#    fc.write('GAUSSIAN  100000.   RESOLUTION\n')
#    fc.write('EOF\n')
#    fc.write('mv -f fort.22 SPECvr2br.bin\n')
        fc.write('\\rm -f fort.*\n')
        fc.write('ln -fs SPECvr2.dat fort.1\n')
    
    fc.write(absetup.binpath+'converfsynnmtoa.exe\n')
    fc.write('mv -f fort.2 SPECvr2.asc\n')
    fc.write('\\rm -f fort.*\n')
    
    fc.write('echo \'# SYNTHE\' > %s\n' % (cwd+'/'+spname))
    fc.write('echo \'# atomdir = %s\' >> %s\n' % (atomdir, cwd+'/'+spname))
    fc.write('echo \'# moldir  = %s\' >> %s\n' % (moldir,  cwd+'/'+spname))
    fc.write('echo \'# atmname = %s\' >> %s\n' % (atmname, cwd+'/'+spname))
    fc.write('''awk '{if (l++ < 2) printf("# %s\\n",$0); else printf("%0.4f  %0.4f  %0.4e\\n",$1,$4,$2)}' SPECvr2.asc >>'''+cwd+'/'+spname+'\n')

    fc.close()
    
    print('SYNTHE: Computing synthetic spectrum '+spname)
    logging.debug('SYNTHE: Computing synthetic spectrum %s' % spname)
    os.chmod(wdir+'/synthe.com',stat.S_IXUSR | stat.S_IRUSR | stat.S_IWUSR)

    p = subprocess.Popen(wdir+'/synthe.com')

    if (wait):
        os.waitpid(p.pid,0)
    else:
        return p
    
## Use MOOG to compute the continuum fluxes    
#
#    fsb = open('synbeg.out','r')
#    s = fsb.readline()
#    fsb.close()
#    ss = s.split()
#    lam1 = float(ss[2]) * 10.
#    lam2 = float(ss[4]) * 10.
#    
#    syntmoog.kur2moog(atmname,'moog.atm', vturb, m)
#    fm = open('batch.par','w')
#        
#    fm.write('doflux\n')
#    fm.write('atmosphere  1\n')
#    fm.write('molecules   1\n')
#    fm.write('trudamp     1\n')
#    fm.write('lines       1\n')
#    fm.write('units       0\n')
#    fm.write('fluxlimits\n')
#    fm.write(' {0:0.3f}   {1:0.3f}   1.0\n'.format(lam1-10.,lam2+10.))
#    fm.write('plot        0\n')
#    fm.write('freeform    0\n')
#    fm.write('iraf        0\n')
#    fm.write('standard_out   \'moog.out\'\n')
#    fm.write('summary_out    \'moog.summ\'\n')
#    fm.write('model_in       \'moog.atm\'\n')
#    fm.close()
#    
#    os.system('moogsilent >& /dev/null')
#
## Read continuum fluxes from moog.out
#
#    lamc, cont = [], []
#	
#    for s in file('moog.out'):
#        if (s.find('CONTINUUM FLUX/INTENSITY') > 0):
#            ss = s.split()
#            scont = string.replace(ss[7],'D','E')
#            lamc = lamc + [float(ss[3])]
#            cont = cont + [float(scont)]
#    
#    fcont = interpolate.interp1d(lamc, cont, bounds_error=False)
#    

# Now format to our "standard" format and write to output file

#    fsp = open(spname,'w')
#    l = 0
#    for s in file(workdir+'/SPECvr2.asc'):
#        if (l < 2):
#            fsp.write('# %s' % (s))
#        else:
#            ss = s.split()
#            lam = float(ss[0])
#            flux = float(ss[1])
#            fnrm = float(ss[3])
##            fscl = fnrm * fcont(lam)
#            
#            fsp.write('{0:0.4f}  {1:0.4f}  {2:0.4e}\n'.format(lam, fnrm, flux))
#
#        l = l+1
#                    
#    fsp.close()


def idlines(m, atmname, idname, lam1, lam2, R=100000., atoms=[], abun=[], decrowd=True):

    teff, nlte = mkmod(m, atmname, 'SPEC.mod', atoms=atoms, abun=abun)
        
    lam1nm = lam1/10.
    lam2nm = lam2/10.

    for j in range(len(atomlines)):
        if (lam1 >= atomlines[j][1] and lam1 <= atomlines[j][2]):
            j1 = j

        if (lam2 >= atomlines[j][1] and lam2 <= atomlines[j][2]):
            j2 = j
            
    if (nlte):
        ifnlte = 1
    else:
        ifnlte = 0        

    fl = open('gftmp.txt','w')
    for j in range(j1, j2+1):
#        for s in file(absetup.catpref+'/cats/Castelli/atoms/'+atomlines[j][0]+'.100'):
        with open(atomdir+atomlines[j][0]+'.100','r') as f:
            for s in f:
                lam = float(s[0:11]) * 10.
                if (lam > lam1 and lam < lam2): fl.write(s)
    fl.close()
    
    with open('gftmp.txt','r') as f:
        elem = list(set([float(s[19:24]) for s in f]))
    elem.sort()
#    print elem

    extlst = ['I','II','III','IV','V','VI','VII','VIII','IX','X']
    
    fout = open(idname,'w')
    
    for ee in elem:
        eid = int(ee)
        eion = int((ee - eid)*100.1)
        if (eid == 1):
            name = 'H'
        elif (eid == 2):
            name = 'He'
        else:
            for kk in absetup.stdabun:
                if (kk[1] == eid): name = kk[0]
        
        if (eion < 10):
            ext = extlst[eion]
        else:
            ext = '*'
            
        print(ee, eid, eion, name, ext)
        
        elam = []
        fl2 = open('gfe.txt','w')

        with open('gftmp.txt','r') as f:
            for s in f:
                ef = float(s[19:24])
                if (ef == ee): 
                    fl2.write(s)
                    elam.append(float(s[0:11])*10.)
                
        fl2.close()    
        
# Run synthe for this element
        
        f = open('idsynthe.com','w')

        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules.dat fort.2\n')
        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/continua.dat fort.17\n')
        f.write(absetup.binpath+'xnfpelsyn.exe< "SPEC.mod" >xnfpelsyn.out\n')
        f.write('mv fort.10 xnftSPEC.dat\n')
        f.write('\\rm -f fort.*\n')
        f.write('ln -fs xnftSPEC.dat fort.10\n')
        f.write(absetup.binpath+'synbeg.exe<<EOF>synbeg.out\n')
#    f.write('AIR       520.0     540.0     500000.   0.     0          10 .001         0   00\n')
        f.write('AIR       %-9.2f %-9.2f %-9.1f 0.     %i          10 .001         0   00\n' % (lam1nm, lam2nm, R, ifnlte))
        f.write('AIRorVAC  WLBEG     WLEND     RESOLU    TURBV  IFNLTE LINOUT CUTOFF        NREAD\n')
        f.write('EOF\n')
    
        f.write('ln -fs gfe.txt fort.11\n')
        f.write(absetup.binpath+'rgfalllinesnew.exe>gfe.out\n')
        f.write('\\rm -f fort.11\n')   

        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/he1tables.dat fort.18\n')
        f.write(absetup.binpath+'synthe.exe>synthe.out\n')
        f.write('ln -fs '+absetup.catpref+'/cats/Castelli/molecules.dat fort.2\n')
        f.write('cat <<EOF >fort.25\n')
        f.write('0.0       0.        1.        0.        0.        0.        0.        0.\n')
        f.write('0.\n')
        f.write('RHOXJ     R1        R101      PH1       PC1       PSI1      PRDDOP    PRDPOW\n')
        f.write('EOF\n')
        f.write(absetup.binpath+'spectrv.exe<SPEC.mod>spectrv.out\n')
        f.write('mv -f fort.7 SPEC.dat\n')
        f.write('ln -fs SPEC.dat fort.1\n')
        f.write(absetup.binpath+'rotate.exe<<EOF>rotate.out\n')
        f.write('    1\n')
        f.write('2.\n')
        f.write('EOF\n')
        f.write('mv -f ROT1 SPECvr2.dat\n')
        f.write('\\rm -f fort.*\n')

        f.write('ln -fs SPECvr2.dat fort.1\n')
        f.write(absetup.binpath+'converfsynnmtoa.exe\n')
        f.write('mv -f fort.2 SPECvr2.asc\n')
        f.write('\\rm -f fort.*\n')
    
        f.write('''awk '{if (l++ > 1) printf("%0.4f  %0.4f  %0.4e\\n",$1,$4,$2)}' SPECvr2.asc > specid.txt\n''')
        f.close()
        
        os.system('source idsynthe.com')
        
# Now find out which entries in the line list produce a significant line in the spectrum

        with open('specid.txt','r') as ff:
            slam = numpy.array([float(s.split()[0]) for s in ff])

        with open('specid.txt','r') as ff:
            sflx = numpy.array([float(s.split()[1]) for s in ff])
        
        if (decrowd):
           for j in range(1,len(slam)-1):
               if (sflx[j] < sflx[j-1]) and (sflx[j] < sflx[j+1]):
                   dlam = list(abs(numpy.array(elam) - slam[j]))
                   dmin = min(dlam)
                   wdmin = dlam.index(dmin)
                   fmin = sflx[j]
                   fout.write('%10.3f %6.2f %6.4f %2s %-5s\n' % (elam[wdmin], ee, fmin, name, ext))
        else:        
            for ll in elam:
                dlam = list(abs(slam - ll))
                dmin = min(dlam)
                wdmin = dlam.index(dmin)
                fmin = sflx[wdmin]
                fout.write('%10.3f %6.2f %6.4f %2s %-5s\n' % (ll, ee, fmin, name, ext))
            
    fout.close()
        


def calcew(atmname, elem, lam, nablog, minlog, dablog, logvt=0.300, lamtol=0.001, lines=None):
# Input:
#  atmname = name of ATLAS model atmosphere
#  elem    = element code, e.g. 11.00 = Na I
#  lam     = wavelength of line, in AA
#  nablog  = Number of abundance offsets (with respect to value in atmosphere)
#  minlog  = Minimum relative abundance
#  dablog  = Log abundance step
#  logvt   = Log of microturbulent velocity (in km/s)
#  lamtol  = Tolerance when matching lambda to line list, in AA
#  lines   = File with line data (Kurucz format). Default: Castelli list
# Returns:
#  Array of equivalent widths (in mAA).

    if (nablog > 99):
        print('ERROR: NABLOG must be less than 100')
        raise Exception('Error: NABLOG too large')
        
    f = open('inpew.com','w')
    f.write('\\rm -f fort.*\n')
    if (lines is None):
        for j in range(len(atomlines)):
            if (lam >= atomlines[j][1] and lam <= atomlines[j][2]):
                f.write('ln -s '+absetup.catpref+'/cats/Castelli/atoms_SLr/'+atomlines[j][0]+'.100 fort.10\n')
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0100.100 fort.10\n')
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0150.100 fort.11\n')
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0200.100 fort.12\n')
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0300.100 fort.13\n')
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0400.100 fort.14\n')
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0500.100 fort.15\n')
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0600.100 fort.16\n') 
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf0800.100 fort.17\n') 
#        f.write('ln -s '+os.getenv('HOME')+'/cats/Castelli/atoms_SLr/gf1200.100 fort.18\n')
#        f.write('inpwidth.exe<<EOF>inpwidth.out\n')
    else:
        f.write('ln -s '+lines+' fort.10\n')
 
    f.write(absetup.binpath+'inpw_sl.exe<<EOF>inpwidth.out\n')
        
    f.write('NoName\n')           # Star name
    f.write('1\n')                # VTUR(bulence)? (yes=1/no=0)
    f.write('1\n')                # Number N of the microturb. velocities (NMAX=3)
    vturb = 10**logvt
    f.write('%0.2f\n' % vturb)    # Values of the N microturb. velocities
    f.write('0\n')                # Is some line file already existing?(yes=1/no=0)
    f.write('0\n')                # LINE ? (yes=1/no=0)
    f.write('0\n')                # PROF(ile)? (yes=1/no=0)
    f.write('0\n')                # AVER(age)? (yes=1/no=0)
    f.write('1\n')                # CURV(e)? (yes=1/no=0)
    f.write('%d %0.3f %0.3f\n' % (nablog, minlog, dablog))
    f.write('1\n')                # LINE ? (yes=1/no=0)
    f.write('%0.4f\n' % (lam/10)) # wavelength (in Nm: es. 4000A=400.0 Nm)
#    if not (lines==None):
    f.write('%0.4f\n' % (lamtol/10)) # matching tolerance (in Nm)
    f.write('10\n')               # equivalent width (in Pm: es. 200 mA=20.0 Pm)
    f.write('%0.2f\n' % elem)     # Code
    f.write('0\n')                # LINE ? (yes=1/no=0)
    f.write('1\n')                # end of data for the lines? (yes=1/no=0)
    f.write('%s\n' % atmname)     # Name of the model file
    f.write('EOF\n')
    
    f.write('mv fort.1 inpwidth.dat\n')
    f.write('\\rm -f fort.*\n')
    f.close()
    
    os.system("bash -c 'source inpew.com'")
    
    f2 = open('width9.com','w')
    f2.write('\\rm -f fort.*\n')
    f2.write('ln -s '+absetup.catpref+'/cats/Castelli2015/molecules.dat fort.2\n')
    f2.write(absetup.binpath+'width9sl.exe < inpwidth.dat > width9.out\n')
    f2.write('\\rm -f fort.*\n')
    f2.close()

    os.system("bash -c 'source width9.com'")
    
    fw = open('width9.out','r')
    s = ''
#    while (not 'NoName' in s): s = fw.readline()
    while not ('VTURB  ABUND' in s): s = fw.readline()
#    s = fw.readline()

    ss = s.split()
    abund = [float(sss) for sss in ss[2:]]

    s = fw.readline()
    ss = s.split()
    logew = [float(sss) for sss in ss[2:]]

    s = fw.readline()
    ss = s.split()
    ew = [10*float(sss) for sss in ss]

    s = fw.readline()
    ss = s.split()
    depth = [float(sss) for sss in ss[1:]]

    s = fw.readline()
    ss = s.split()
    resid = [float(sss) for sss in ss[1:]]

    s = fw.readline()
    ss = s.split()
    contin = float(ss[1])
  
    fw.close()
    
    output = {'abund': abund, 'logew': logew, 'ew': ew, 
            'depth': depth, 'resid': resid, 'contin': contin}
    
    
    return output
    
    
    
def ewlst(atmname, linelst, logvt=0.300):
# List equivalent widths for lines in 'linelst'. 
# Also includes information on residual depth and HFS*ISO weights.
# Calculations are done with WIDTH9
#
    extlst = ['I','II','III','IV','V','VI','VII','VIII','IX','X']
    
    print("# Lambda     loggf      E       Code  f(HFS*ISO) Resid     EW   Species ")
    
    with open(linelst,'r') as f:
        for l in f:
            ftmp = open('linestmp.lst','w')
            ftmp.write(l)
            ftmp.close()
            lam  = float(l[0:11])*10
            loggf = float(l[11:18])
            E1 = float(l[24:36])
            E2 = float(l[52:64])
            E = min(E1, E2)
            
            elem = float(l[19:24])
            loghfs = float(l[109:115])
            logiso = float(l[118:124])
            fhfsiso = 10**(loghfs+logiso)
            
            eid = int(elem)
            eion = int((elem - eid)*100.1)
            if (eid == 1):
                name = 'H'
            elif (eid == 2):
                name = 'He'
            else:
                for kk in absetup.stdabun:
                    if (kk[1] == eid): name = kk[0]
            
            if (eion < 10):
                ext = extlst[eion]
            else:
                ext = '*'
    
            
            ewdat = calcew(atmname, elem, lam, nablog=1, minlog=0., dablog=0., logvt=logvt, lines='linestmp.lst')
    #        print ewdat
            print("%10.3f %7.3f %12.3f %5.2f   %5.3f   %6.4f %7.1f  %2s %-5s" % (lam, loggf, E, elem, fhfsiso, ewdat['resid'][0], ewdat['ew'][0], name, ext))