# encoding: utf-8
#
# @Author: Jon Holtzman
# @Date: March 2018
# @Filename: synth.py
# @License: BSD 3-Clause
# @Copyright: Jon Holtzman

# routines related to creation of synthetic spectra

from __future__ import division
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals

import os
import errno
import glob
import pdb
import shutil
import subprocess
import sys
import tempfile
import multiprocessing as mp
import numpy as np
import time
import struct
import yaml
from apogee.aspcap import aspcap
from apogee.aspcap import ferre
from apogee.speclib import atmos
from apogee.speclib import lsf
from apogee.speclib import sample
from apogee.utils import atomic
from apogee.utils import spectra
from sdss import yanny
from astropy.io import ascii
from astropy.io import fits
from astropy.table import Table
from sklearn.decomposition import PCA
from sklearn.decomposition import IncrementalPCA
import matplotlib
import matplotlib.pyplot as plt
from tools import plots
from tools import match
from tools import html
from synple import synple

colors=['r','g','b','c','m','y']

def showtime(string) :
    """ Utility routine to print a string and clock time, with flush to stdout

    Args:
        string(str) : string to print
    """
    print(string+' {:8.2f}'.format(time.time()))
    sys.stdout.flush()

def kurucz2turbo(infile,outfile,trim=0) :
    """ Convert Kurucz model atmosphere for use by Turbospectrum 
        Allows for trimming of layers

    Args:
        infile (str)  : name of input Kurucz atmosphere file
        outfile (str) : name of output file
        trim  (int)   : number of layers to remove (default=0)
    """
    try :
        fp=open(infile,'r')
    except :
        fp=open(infile+'.filled','r')
    fout=open(outfile,'w')
    layer=-1
    nlayers=-1
    for iline,line in enumerate(fp) :
        if iline == 0 : 
            # get parameters
            teff=line.split()[1]
            logg=float(line.split()[3])
        if line.split()[0] == 'READ' :
            # write header line once we know how many layers
            nlayers=int(line.split()[2])
            layer=0
            fout.write('KURUCZ {:d} 5000. {:8.2f} 0 0.\n'.format(nlayers-trim,logg))
        elif layer >= 0 and layer < nlayers:
            # write the layers
            layer+=1
            if layer > trim: fout.write(line)
    fp.close()
    fout.close()

def marcs2turbo(infile,outfile,trim=0,fill=True) :
    """ Prepare MARCS input model for Turbospectrum, allowing for trimming of layers

    Args :
        infile (str)  : name of input MARCS file
        outfile (str) : name of output MARCS fie
        trim (int)    : number of layers to trim (default=0)

    Returns :
        0 if successful, -1 if not
    """
    try:
        fp=open(infile,'r')
    except :
        if not fill : return -1
        try :
            fp=open(infile+'.filled','r')
        except :
            raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), filename)
            return -1

    fout=open(outfile,'w')
    layer = 0
    for iline,line in enumerate(fp) :
        if line.find('Number of depth') >= 0 :
            nlayers = int(line.split()[0])
            layer = 1
            fout.write(str(nlayers-trim)+' Number of depth points\n')
        else :
            try :
                i = int(line.split()[0])
                if i > trim : fout.write(line)
            except :
                fout.write(line)
    return 0

def get_atmod_file(teff,logg,mh,am,cm,nm,atmos_type='marcs',atmosroot=None,nskip=0,fill=True,workdir='./',atmosdir=None) :
    """ Copy atmosphere file to working directory in Turbospec format,  trimming as requested
    """
    # directory setup
    if atmosroot is None : atmosroot=os.environ['APOGEE_SPECLIB']+'/atmos/'
 
    # atmosphere and filename set up 
    geo = 'p'
    if atmos_type == 'kurucz' : 
        atmoscode= 'k'
        model = 'Kurucz'
        if atmosdir is None : atmosdir = '/kurucz/'
    elif atmos_type == 'marcs' : 
        atmoscode = 'm'
        model = 'MARCS'
        if atmosdir is None : atmosdir = '/marcs/MARCS_v3_2016/'
        if logg <= 3.001 : geo = 's'
    elif atmos_type == 'synspec' :
        atmoscode = 's'
        model = 'SYNSPEC'
        if atmosdir is None : atmosdir = '/marcs/MARCS_v3_2016/nltepop/'
        if logg <= 3.001 : geo = 's'
    else :
        print('unknown atmos_type: ', atmos_type)
        pdb.set_trace()
    atmosdir=atmosroot+'/'+atmosdir+'/'
    atmod=atmosdir+atmos.filename(teff,logg,mh,cm,am,model=model)

    # atmosphere: make local copy in Turbospectrum input format, allowing for trimmed layers
    # note that [N/M] is solar for atmospheres (since it doesn't have a big effect)
    outmod = workdir+'/'+os.path.basename(atmod)
    if atmos_type == 'kurucz' :
        if nskip == 0 : trim=0
        if nskip == 1 : trim=7
        if nskip == 2 : trim=15
        if nskip > 2 : return 0.,0.
        kurucz2turbo(atmod,outmod,trim=trim )

    elif atmos_type == 'marcs' :
        try :        
            ret = marcs2turbo(atmod,outmod,trim=nskip,fill=fill )
            if not fill and ret<0 : return ret
        except:
            return -2

    elif atmos_type == 'synspec' :
        if os.path.exists(atmod) :
            shutil.copy(atmod,outmod)
            shutil.copy(atmod.replace('.22','.5'),outmod.replace('.22','.5'))
            return outmod
        elif cm<-1 :
            atmod=atmosdir+atmos.filename(teff,logg,mh,-1,am,model=model)
            if os.path.exists(atmod) :
                shutil.copy(atmod,outmod)
                shutil.copy(atmod.replace('.22','.5'),outmod.replace('.22','.5'))
                return outmod
        else :
            tmp = np.round(cm/0.25)*0.25
            atmod=atmosdir+atmos.filename(teff,logg,mh,tmp,am,model=model)
        return -1

    return outmod

def get_workdir(teff,logg,mh,am,cm,nm,atmos_type='marcs',solarisotopes=False,save=False,elemgrid=None,vmicro=1.) :
    """ Create work directory name for synthesis
    """
    if atmos_type == 'kurucz' : 
        atmoscode= 'k'
    elif atmos_type == 'marcs' : 
        atmoscode = 'm'
    elif atmos_type == 'synspec' : 
        atmoscode = 's'
    else :
        print('unknown atmosphere type: ', atmos_type)
        pdb.set_trace()

    if abs(solarisotopes)==1: prefix=atmoscode+'d' 
    else : prefix=atmoscode+'g'

    if save :
        workdir=(prefix+'m{:s}a{:s}c{:s}n{:s}v{:s}'+elemgrid).format(
                 atmos.cval(mh),atmos.cval(am),atmos.cval(cm),atmos.cval(nm),atmos.cval(vmicro))
        workdir=os.environ['APOGEE_LOCALDIR']+'/'+workdir
        try: os.mkdir(workdir)
        except: pass
    else :
        workdir=tempfile.mkdtemp(dir=os.environ['APOGEE_LOCALDIR'],suffix=os.environ['HOSTNAME'])

    return workdir

def mk_synthesis(code,teff,logg,mh,am,cm,nm,wrange=[15100.,17000],dw=0.05,vmicro=2.0,solarisotopes=False,elemgrid='',welem=None,
    nlte=False,els=None,atmod=None,atmos_type='marcs',atmosroot=None,atmosdir=None,nskip=0,fill=True,cmnear=False,dospherical=True,
    linelist='20180901',h2o=0,linelistdir=None,atoms=True,molec=True, msuffix='', save=False,run=True, verbose=False) :
    """ Do synthesis

    Args:
        code (str) : 'turbospec' or 'synspec'
        teff (int) : Effective temperature
        logg (float) : log(surface gravity)
        mh (float ) : [M/H]
        am (float ) : [alpha/M]
        cm (float ) : [C/M]
        nm (float ) : [N/M]
        wrange (list) : wavelength range (default [15100,17000.])
        dw (float)  : wavelength spacing (default 0.05 A)
        vmicro (float ) : microturbulent velocity (default 2.0 km/s)
        solarisotopes (bool) : use solar isotope ratios, else "giant" isotope ratios ( default False )
        elemgrid (str) :  name of element for minigrid calculation
        welem (str) :  NOT IMPLEMENTED set of wavelength ranges
        els (list of pairs) : list of [element name, abundance] pairs
        atmod (str) : name of atmosphere model (default=None, model is determined from input parameters)
        atmosroot (str) : root atmosphere directory (default=None --> $APOGEE_SPECLIB/atmos)
        atmosdir (str) :  atmosphere directory under atmosroot (default=None --> kurucz or marcs/MARCS_v3_2016, depending on kurucz boolean)
        nskip (int)  : number of layers to strip (for kurucz, nskip*7), default=0
        fill (bool)  : use filled atmosphere hole if model requests it ( default = True), else skip synthesis
        linelist (str) : name of linelist (default='20150714')
        linelistdir (str) : directory for linelist (default=None --> $APOGEE_SPECLIB/linelists)
        h2o (int or None) : H2O linelist to use: 0=chooses by Teff, [M/H], [a/M], 1=8.5V, 2=9.5V, None= none
        save (bool ) : save temporary directory and files for synthesis (default=False)
        run (bool) :  actually do the synthesis (default=True)

    Returns:
        spec (np.array) : synthetic spectrum
        specnorm (np.array) : synthetic spectrum, normalized
    """
 
    # output directory and filename
    workdir = get_workdir(teff,logg,mh,am,cm,nm,atmos_type='marcs',solarisotopes=solarisotopes,save=save,elemgrid=elemgrid,vmicro=vmicro)
    root=(atmos_type+'_t{:04d}g{:s}m{:s}a{:s}c{:s}n{:s}v{:s}'+elemgrid).format(int(teff), atmos.cval(logg), 
                      atmos.cval(mh), atmos.cval(am), atmos.cval(cm), atmos.cval(nm),atmos.cval(vmicro))

    # atmosphere: make local copy in Turbospectrum input format, allowing for trimmed layers
    # note that [N/M] is solar for atmospheres (since it doesn't have a big effect)
    if nlte : atmos_type='synspec'
    if cmnear : 
        #if cm < am+0.125 or cm>am+0.375 : cmtmp = np.round(cm/0.25)*0.25 
        #else : cmtmp = np.floor(cm/0.25)*0.25 
        #if np.isclose(cmtmp,0.) : cmtmp=0. 
        #amtmp = np.round(am/0.25)*0.25
        #if np.isclose(amtmp,0.) : amtmp=0. 

        #new
        cmtmp, amtmp = get_cmnear(cm, am)

    else : 
        cmtmp = cm
        amtmp = am
    if atmod is None : atmod = get_atmod_file(teff,logg,mh,amtmp,cmtmp,nm,atmos_type=atmos_type,nskip=nskip, atmosroot=atmosroot,atmosdir=atmosdir,workdir=workdir)
    elif atmod == 'interp' and atmos_type == 'marcs' : 
        atmod= interpol_marcs(teff,logg,mh,amtmp,cmtmp,nm,workdir=workdir)
        atmos_type = 'marcs-interp'
    else : shutil.copy(atmod,workdir+'/'+os.path.basename(atmod))

    if type(atmod) is int :
        if atmod == -1 :        
            print('HOLE NOT SYNTHESIZED: ', atmod)
            return 0.,0.
        else :
            print('PROBLEM: ',atmod)
            fail('mkturbospec problem: ',atmod)
            return 0.,0.

    # linelists
    if linelistdir is None : linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/' 
    if code == 'turbospec' :
        linelists = []
        n_HI = len(open(linelistdir+'/turbospec.'+linelist+'.Hlinedata').readlines())
        if n_HI > 2 : linelists.append(linelistdir+'/turbospec.'+linelist+'.Hlinedata')
        if atoms : linelists.append(linelistdir+'/turbospec.'+linelist+'.atoms')
        if molec : 
            linelists.append(linelistdir+'/turbospec.'+linelist+msuffix+'.molec')
            tfactor=1
            if h2o == 0 :
              if teff <= 4000 :
                # DR16
                #if mh+am < -1.5 or teff > 3250 :
                #    linelists.append(linelistdir+'/turbospec.h2o-BC8.5V'+'.molec')
                #    tfactor=2
                #else  :
                #    linelists.append(linelistdir+'/turbospec.h2o-BC9.5V'+'.molec')
                #    tfactor=5
                if teff <= 4000 : 
                    if mh+am < -1.5 or teff >= 3500 :
                        linelists.append(linelistdir+'/turbospec.h2o_POKAZATEL-BC8.5V'+'.molec')
                        tfactor=2
                    elif teff > 3250 :
                        linelists.append(linelistdir+'/turbospec.h2o_POKAZATEL-BC9.0V'+'.molec')
                        tfactor=4
                    else  :
                        linelists.append(linelistdir+'/turbospec.h2o_POKAZATEL-BC9.5V'+'.molec')
                        tfactor=5
            elif h2o == 1 :
                linelists.append(linelistdir+'/turbospec.h2o-BC8.5V'+'.molec')
                tfactor=2
            elif h2o == 2 :
                linelists.append(linelistdir+'/turbospec.h2o-BC9.5V'+'.molec')
                tfactor=5
            elif h2o == 3 :
                linelists.append(linelistdir+'/turbospec.h2o'+'.molec')
                tfactor=5
    elif code == 'synspec' :
        linelists = [linelistdir+'/synspec/synspec.'+linelist+'.atoms']
        if molec :
            #if solarisotopes : linelists.append(linelistdir+'/synspec/synspec.'+linelist+'sun_nofeh_noh2o_noc2.molec')
            #else : linelists.append(linelistdir+'/synspec/synspec.'+linelist+'giant_nofeh_noh2o_noc2.molec')
            if solarisotopes : linelists.append(linelistdir+'/synspec/synspec.'+linelist+'sun'+msuffix+'.molec')
            else : linelists.append(linelistdir+'/synspec/synspec.'+linelist+'giant'+msuffix+'.molec')
            if h2o == 0 :
                if teff <= 4000 : 
                    if mh+am < -1.5 or teff >= 3500 :
                        linelists.append(linelistdir+'/synspec/synspec.h2o-BC8.5V'+'.molec')
                    elif teff > 3250 :
                        linelists.append(linelistdir+'/synspec/synspec.h2o-BC9.0V'+'.molec')
                    else  :
                        linelists.append(linelistdir+'/synspec/synspec.h2o-BC9.5V'+'.molec')
            elif h2o == 1 :
                linelists.append(linelistdir+'/synspec/synspec.h2o-BC8.5V.molec')
            elif h2o == 2 :
                linelists.append(linelistdir+'/synspec/synspec.h2o-BC9.5V.molec')
            elif h2o == 3 :
                linelists.append(linelistdir+'/synspec/synspec.h2o.molec')
            elif h2o == 4 :
                linelists.append(linelistdir+'/synspec/synspec.h2o-BC9.0V.molec')
            elif h2o == 5 :
                # DR16 match
                if teff < 4000 : 
                    if mh+am < -1.5 or teff > 3250 :
                        linelists.append(linelistdir+'/synspec/synspec.h2o_barber-BC8.5V'+'.molec')
                    else  :
                        linelists.append(linelistdir+'/synspec/synspec.h2o_barber-BC9.5V'+'.molec')
        #linelists = ['apogeeDR16.20180901.19','apogeeDR16_arc.20']
    else :
        print('unknown code!')
        pdb.set_trace()
    if verbose : print('linelists: ',linelists)

    # default abundances
    abundances = atomic.solar()
    abundances[2:] += mh
    abundances[6-1] += cm
    abundances[7-1] += nm
    for i in [8,10,12,14,16,18,20,22] : abundances[i-1] += am

    # abundance overrides from els, given as [X/M]
    if els is not None :
        for el in els :
            atomic_num = atomic.periodic(el[0])
            abundances[atomic_num-1] = atomic.solar(el[0]) + mh + el[1] 

    if elemgrid != '' :
        elemnum=atomic.periodic(elemgrid)[0]
        elem0=abundances[elemnum-1]
        eabun=np.arange(-0.75,1.20,0.25)
    else :
        eabun=np.array([0.])

    cwd = os.getcwd()
    os.chdir(workdir)

    # if turbospec do the opacity calculations first, then synthesis for all eabun
    if code == 'turbospec' :
        out = do_turbospec(root,atmod,linelists,mh,am,abundances,wrange,dw,save=save,run=run,
                           solarisotopes=solarisotopes,bsyn=False,atmos_type=atmos_type,vmicro=vmicro)
    elif code == 'synspec' :
        abundances = 10.**(np.array(abundances)-abundances[0])
        model_read = synple.read_model(atmod)
        diff = np.where(~np.isclose(abundances,model_read[4]))[0]
        for d in diff : 
            print('synthesis - model mismatch: ',atomic.periodic(d+1),np.log10(abundances[d]/model_read[4][d]))

    # if this is an elemgrid, loop over abundances
    for ielem,abun in enumerate(eabun) :
        file = root+'{:02d}'.format(ielem)
        if elemgrid != '' :
            abundances[elemnum-1] = elem0+abun
     
        if code == 'turbospec' :
            if dospherical and ('marcs' in atmos_type) and logg <= 3.001 :  spherical= True
            else : spherical = False
            if ielem == 0 : 
                wave,flux,fluxnorm = do_turbospec(file,atmod,linelists,mh,am,abundances,wrange,dw,
                                                  save=save,run=run,solarisotopes=solarisotopes,
                                                  babsma=root+'opac',atmos_type=atmos_type,spherical=spherical,tfactor=tfactor)
            else :
                out = do_turbospec(file,atmod,linelists,mh,am,abundances,wrange,dw,
                                   save=save,run=run,solarisotopes=solarisotopes,
                                   babsma=root+'opac',atmos_type=atmos_type,spherical=spherical,tfactor=tfactor)
        elif code == 'synspec' :
            print('synple.syn: ',atmod,wrange,linelists,dw,vmicro,save)
            # add 10*dw to range because of synple's spline interpolation
            wtmp=[wrange[0]-10*dw,wrange[1]+10*dw]
            try:
                wave,flux,cont = synple.syn(atmod,wtmp,linelist=linelists,dw=dw,vmicro=vmicro,abu=abundances,save=save,clean=not save)
            except:
                return 0., 0.
            wave=wave[10:-10]
            flux=flux[10:-10]
            cont=cont[10:-10]
            fluxnorm = flux/cont
        else :
            print('unknown synth code!')
            pdb.set_trace()

        # load into final arrays
        if ielem == 0 :
            spec = flux
            specnorm = fluxnorm
        else :
            spec=np.vstack([spec,flux])
            specnorm=np.vstack([specnorm,fluxnorm])

    os.chdir(cwd)
    if run :
        if not save : shutil.rmtree(workdir)
        return spec, specnorm

def get_cmnear(cm,am) :
    """ given input [C/M] and [alpha/M], return nearest model
        grid point, but try to ensure that it is on the same
        side of the C/O boundary

    """
    amtmp = np.round(am/0.25)*0.25
    cmtmp=np.round(cm/0.25)*0.25
    # if synthesis has C/O<-0.025 ([C/O]<0.245), make sure adoped model does 
    if (cm-am < 0.245) and (cmtmp-amtmp > 0.245) : cmtmp -= 0.25
    # if synthesis has C/O>-0.025 ([C/O]>0.245), make sure adoped model does
    if (cm-am > 0.245) and (cmtmp-amtmp < 0.245) : cmtmp += 0.25
    if np.isclose(cmtmp,0.) : cmtmp=0.
    if np.isclose(amtmp,0.) : amtmp=0.
    #print('cmnear: ',cm,cmtmp,am,amtmp)
    return cmtmp, amtmp


def do_turbospec(file,atmod,linelists,mh,am,abundances,wrange,dw,save=False,run=True,solarisotopes=False,babsma=None,bsyn=True,
                 atmos_type='marcs',spherical=True,vmicro=1.0,tfactor=1.,verbose=False ) :
    """ Runs Turbospectrum for specified input parameters

    Args:

    Returns:
    """

    # Turbospectrum setup
    try: os.symlink(os.environ['APOGEE_DIR']+'/src/turbospec/DATA','./DATA')
    except: pass

    if verbose : stdout = None
    else : stdout = open(os.devnull, 'w')

    # individual element grid?
    nels = len(abundances)

    welem=np.array(wrange)
    # only compute opacities for a single nominal abundance
    if babsma is None :
        fout=open(file+'_babsma.csh','w')
        fout.write("#!/bin/csh -f\n")
        fout.write("{:s}/bin/babsma_lu << EOF\n".format(os.environ['APOGEE_DIR']))
        fout.write("'LAMBDA_MIN:'   '{:12.3f}'\n".format(welem.min()-dw))
        fout.write("'LAMBDA_MAX:'   '{:12.3f}'\n".format(welem.max()+dw))
        fout.write("'LAMBDA_STEP:'  '{:8.3f}'\n".format(dw))
        fout.write("'MODELINPUT:'  '{:s}'\n".format(os.path.basename(atmod)))
        if atmos_type != 'marcs' : fout.write("'MARCS-FILE:'  '.false.'\n")
        fout.write("'MODELOPAC:'  '{:s}opac'\n".format(os.path.basename(file)))
        fout.write("'METALLICITY:'  '{:8.3f}'\n".format(mh))
        fout.write("'ALPHA/Fe:'  '{:8.3f}'\n".format(am))
        fout.write("'HELIUM:'  '{:8.3f}'\n".format(0.00))
        fout.write("'R-PROCESS:'  '{:8.3f}'\n".format(0.00))
        fout.write("'S-PROCESS:'  '{:8.3f}'\n".format(0.00))
        fout.write("'INDIVIDUAL ABUNDANCES:'  '{:2d}'\n".format(nels))
        for iel,abun in enumerate(abundances) :
            fout.write("{:5d}  {:8.3f}\n".format(iel+1,abun))
        if not solarisotopes :
          fout.write("'ISOTOPES:'  '2'\n")
          # adopt ratio of 12C/13C=15
          fout.write("   6.012 0.9375\n")
          fout.write("   6.013 0.0625\n")
        fout.write("'XIFIX:'  'T'\n")
        fout.write("{:8.3f}\n".format(vmicro))
        fout.write("EOF\n")
        fout.close()
        if run :
            os.chmod(file+'_babsma.csh', 0o777)
            subprocess.call(['time','./'+os.path.basename(file)+'_babsma.csh'],stdout=stdout,stderr=stdout)
        babsma = os.path.basename(file)+'opac'

    if not bsyn : return

    # create bsyn control file
    bsynfile = file
    fout = open(bsynfile+'.inp','w')
    fout.write("'LAMBDA_STEP:'  '{:8.3f}'\n".format(dw))
    fout.write("'LAMBDA_MIN:'   '{:12.3f}'\n".format(welem.min()))
    fout.write("'LAMBDA_MAX:'   '{:12.3f}'\n".format(welem.max()))
    fout.write("'INTENSITY/FLUX:'  'Flux'\n")
    fout.write("'COS(THETA):'  '1.00'\n")
    fout.write("'ABFIND:'  '.false.'\n")
    fout.write("'MODELINPUT:'  '{:s}'\n".format(os.path.basename(atmod)))
    if atmos_type != 'marcs' : fout.write("'MARCS-FILE:'  '.false.'\n")
    fout.write("'MODELOPAC:'  '{:s}'\n".format(babsma))
    fout.write("'RESULTFILE:'  '{:s}'\n".format(os.path.basename(file)))
    fout.write("'METALLICITY:'  '{:8.3f}'\n".format(mh))
    fout.write("'ALPHA/Fe:'  '{:8.3f}'\n".format(am))
    fout.write("'HELIUM:'  '{:8.3f}'\n".format(0.00))
    fout.write("'R-PROCESS:'  '{:8.3f}'\n".format(0.00))
    fout.write("'S-PROCESS:'  '{:8.3f}'\n".format(0.00))
    fout.write("'INDIVIDUAL ABUNDANCES:'  '{:2d}'\n".format(len(abundances)))
    for iel,abun in enumerate(abundances) :
        fout.write("{:5d}  {:8.3f}\n".format(iel+1,abun))
    if  not solarisotopes :
        fout.write("'ISOTOPES:'  '2'\n")
        # adopt ratio of 12C/13C=15
        fout.write("   6.012 0.9375\n")
        fout.write("   6.013 0.0625\n")
    fout.write("'NFILES:'  '{:4d}'\n".format(len(linelists)))
    for linelist in linelists: 
        fout.write(linelist+"\n")
    if spherical: fout.write("'SPHERICAL:'  'T'\n")
    else : fout.write("'SPHERICAL:'  'F'\n")
    fout.write("30\n")
    fout.write("300.00\n")
    fout.write("15\n")
    fout.write("1.3\n")
    fout.close()

    # control file, with special handling in case bsyn goes into infinite loop ...
    fout = open(file+"_bsyn.csh",'w')
    fout.write("#!/bin/csh -f\n")
    fout.write("{:s}/bin/bsyn_lu < {:s} &\n".format(os.environ['APOGEE_DIR'],os.path.basename(bsynfile)+'.inp'))
    fout.write('set bsynjob = $!\n')
    fout.write("set ok = 0\n")
    fout.write("set runtime = `ps -p $bsynjob -o cputime | tail -1 | awk -F: '{print ($1*3600)+($2*60)+$3}'`\n")
    tmax=120*int(.05/min([0.05,dw]))*tfactor
    fout.write('while ( $runtime < {:d} && $ok == 0 )\n'.format(tmax))
    fout.write('  usleep 200000\n')
    fout.write("  set runtime = `ps -p $bsynjob -o cputime | tail -1 | awk -F: '{print ($1*3600)+($2*60)+$3}'`\n")
    fout.write('  if ( `ps -p $bsynjob -o comm=` == "" ) then\n')
    fout.write('    echo process done, exiting!\n')
    fout.write('    set ok = 1\n')
    fout.write('  endif\n')
    fout.write('end\n')
    fout.write('if ( $ok == 0 ) then\n')
    fout.write('  echo expired, killing job\n')
    fout.write('  kill $bsynjob\n')
    fout.write('endif\n')
    fout.close()
    if run :
        os.chmod(file+'_bsyn.csh', 0o777)
        subprocess.call(['time','./'+os.path.basename(file)+'_bsyn.csh'],stdout=stdout,stderr=stdout)
        try:
            out=np.loadtxt(file)
            wave=out[:,1]
            specnorm=out[:,1]
            spec=out[:,2]
        except :
            print('failed...',file,atmod,mh,am)
            return 0.,0.,0.
        return wave,spec,specnorm

def mkturbospec(teff,logg,mh,am,cm,nm,wrange=[15100.,17000],dw=0.05,vmicro=2.0,solarisotopes=False,elemgrid='',welem=None,
    els=None,atmod=None,kurucz=True,atmosroot=None,atmosdir=None,nskip=0,fill=True,
    linelist='20150714',h2o=0,linelistdir=None,atoms=True,molec=True,
    save=False,run=True) :
    """ Runs Turbospectrum for specified input parameters

    Args:
        teff (int) : Effective temperature
        logg (float) : log(surface gravity)
        mh (float ) : [M/H]
        am (float ) : [alpha/M]
        cm (float ) : [C/M]
        nm (float ) : [N/M]
        wrange (list) : wavelength range (default [15100,17000.])
        dw (float)  : wavelength spacing (default 0.05 A)
        vmicro (float ) : microturbulent velocity (default 2.0 km/s)
        solarisotopes (bool) : use solar isotope ratios, else "giant" isotope ratios ( default False )
        elemgrid (str) :  name of element for minigrid calculation
        welem (str) :  NOT IMPLEMENTED set of wavelength ranges
        els (list of pairs) : list of [element name, abundance] pairs
        atmod (str) : name of atmosphere model (default=None, model is determined from input parameters)
        atmosroot (str) : root atmosphere directory (default=None --> $APOGEE_SPECLIB/atmos)
        atmosdir (str) :  atmosphere directory under atmosroot (default=None --> kurucz or marcs/MARCS_v3_2016, depending on kurucz boolean)
        nskip (int)  : number of layers to strip (for kurucz, nskip*7), default=0
        fill (bool)  : use filled atmosphere hole if model requests it ( default = True), else skip synthesis
        linelist (str) : name of linelist (default='20150714')
        linelistdir (str) : directory for linelist (default=None --> $APOGEE_SPECLIB/linelists)
        h2o (int or None) : H2O linelist to use: None=none, 1=8.5V, 2=9.5V, 0=chooses by Teff and [alpha/H]
        save (bool ) : save temporary directory and files for synthesis (default=False)
        run (bool) :  actually do the synthesis (default=True)

    Returns:
        spec (np.array) : synthetic spectrum
        specnorm (np.array) : synthetic spectrum, normalized
    """

    # directory setup
    if atmosroot is None : atmosroot=os.environ['APOGEE_SPECLIB']+'/atmos/'
    if linelistdir is None : linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/' 
 
    # atmosphere and filename set up 
    geo = 'p'
    if kurucz : 
        atmoscode= 'k'
        model = 'Kurucz'
        if atmosdir is None : atmosdir = '/kurucz/'
    else : 
        atmoscode = 'm'
        model = 'MARCS'
        if atmosdir is None : atmosdir = '/marcs/MARCS_v3_2016/'
        if logg <= 3.001 : geo = 's'
    atmosdir=atmosroot+'/'+atmosdir+'/'

    if abs(solarisotopes)==1: prefix=atmoscode+'d' 
    else : prefix=atmoscode+'g'

    # output directory and filename
    if save :
        workdir=(prefix+'m{:s}a{:s}c{:s}n{:s}v{:s}'+elemgrid).format(
                 atmos.cval(mh),atmos.cval(am),atmos.cval(cm),atmos.cval(nm),atmos.cval(vmicro))
        workdir=os.environ['APOGEE_LOCALDIR']+'/'+workdir
        try: os.mkdir(workdir)
        except: pass
    else :
        workdir=tempfile.mkdtemp(dir=os.environ['APOGEE_LOCALDIR'],suffix=os.environ['HOSTNAME'])

    root=workdir+'/'+(prefix+'t{:04d}g{:s}m{:s}a{:s}c{:s}n{:s}v{:s}'+elemgrid).format(int(teff), atmos.cval(logg), 
                      atmos.cval(mh), atmos.cval(am), atmos.cval(cm), atmos.cval(nm),atmos.cval(vmicro))
    print(root)
    if save :
        stdout = None
    else :
        stdout = open(os.devnull, 'w')

    # atmosphere: make local copy in Turbospectrum input format, allowing for trimmed layers
    # note that [N/M] is solar for atmospheres (since it doesn't have a big effect)
    if atmod is None :
        atmod=atmosdir+atmos.filename(teff,logg,mh,cm,am,model=model)
    if kurucz :
        if nskip == 0 : trim=0
        if nskip == 1 : trim=7
        if nskip == 2 : trim=15
        if nskip > 2 : return 0.,0.
        kurucz2turbo(atmod,workdir+'/'+os.path.basename(atmod),trim=trim )
    else :
        try :        
            ret = marcs2turbo(atmod,workdir+'/'+os.path.basename(atmod),trim=nskip,fill=fill )
            if not fill and ret<0 :
                print('HOLE NOT SYNTHESIZED: ', atmod)
                return 0.,0.
        except:
            print('PROBLEM: ',atmod)
            fail('mkturbospec problem: '+atmod)
            return 0.,0.

    # Turbospectrum setup
    try: os.symlink(os.environ['APOGEE_DIR']+'/src/turbospec/DATA',workdir+'/DATA')
    except: pass

    # individual element grid?
    nels = 3
    if elemgrid != '' :
        elemnum=atomic.periodic(elemgrid)[0]
        eabun=atomic.solar(elemgrid)[0]+np.arange(-0.75,1.20,0.25)
    else :
        eabun=np.array([0.])
    if els is not None :
        nels += len(els)
        if 'C' in np.array(els)[:,0] : nels-=1
        if 'N' in np.array(els)[:,0] : nels-=1
        if 'O' in np.array(els)[:,0] : nels-=1

    # welem only computes in windows, but it is easier/faster to compute the whole range with a linelist that only 
    #   has lines in windows!
    wrange=np.array(wrange)
    if welem is not None :
        print('welem not yet implemented!')
        pdb.set_trace()
    else :
        welem=wrange

    sz=welem.shape
    if len(welem) > 1 : 
        # not tested!
        nrange = welem[0]
    else :
        drange=1

    # loop over individual element abundances
    for ielem,abun in enumerate(eabun) :
        file = root+'{:02d}'.format(ielem)
        # only compute opacities for a single nominal abundance
        if ielem == 0 :
            fout=open(root+'_babsma.csh','w')
            fout.write("#!/bin/csh -f\n")
            fout.write("{:s}/bin/babsma_lu << EOF\n".format(os.environ['APOGEE_DIR']))
            fout.write("'LAMBDA_MIN:'   '{:12.3f}'\n".format(welem.min()-dw))
            fout.write("'LAMBDA_MAX:'   '{:12.3f}'\n".format(welem.max()+dw))
            fout.write("'LAMBDA_STEP:'  '{:8.3f}'\n".format(dw))
            fout.write("'MODELINPUT:'  '{:s}'\n".format(os.path.basename(atmod)))
            if kurucz : fout.write("'MARCS-FILE:'  '.false.'\n")
            fout.write("'MODELOPAC:'  '{:s}opac'\n".format(os.path.basename(root)))
            fout.write("'METALLICITY:'  '{:8.3f}'\n".format(mh))
            fout.write("'ALPHA/Fe:'  '{:8.3f}'\n".format(am))
            fout.write("'HELIUM:'  '{:8.3f}'\n".format(0.00))
            fout.write("'R-PROCESS:'  '{:8.3f}'\n".format(0.00))
            fout.write("'S-PROCESS:'  '{:8.3f}'\n".format(0.00))
            fout.write("'INDIVIDUAL ABUNDANCES:'  '{:2d}'\n".format(nels))
            if els is None or not 'C' in np.array(els)[:,0] :
                fout.write("    6  {:8.3f}\n".format(8.39+mh+cm))
            if els is None or not 'N' in np.array(els)[:,0] :
                fout.write("    7  {:8.3f}\n".format(7.78+mh+nm))
            if els is None or not 'O' in np.array(els)[:,0] :
                fout.write("    8  {:8.3f}\n".format(8.66+mh+am))
            if els is not None :
                for el in els :
                  num=atomic.periodic(el[0])[0]
                  ab=atomic.solar(el[0])[0]+el[1]
                  fout.write("{:5d}  {:8.3f}\n".format(num,ab+mh))
            if not solarisotopes :
              fout.write("'ISOTOPES:'  '2'\n")
              # adopt ratio of 12C/13C=15
              fout.write("   6.012 0.9375\n")
              fout.write("   6.013 0.0625\n")
            fout.write("'XIFIX:'  'T'\n")
            fout.write("{:8.3f}\n".format(vmicro))
            fout.write("EOF\n")
            fout.close()
            if run :
                os.chmod(root+'_babsma.csh', 0o777)
                cwd = os.getcwd()
                os.chdir(workdir)
                subprocess.call(['time','./'+os.path.basename(root)+'_babsma.csh'],stdout=stdout)
            if elemgrid != '' : nels+=1

        # create bsyn control file
        bsynfile = file+'bsyn{:02d}.inp'.format(ielem)
        fout = open(bsynfile,'w')
        fout.write("'LAMBDA_STEP:'  '{:8.3f}'\n".format(dw))
        fout.write("'LAMBDA_MIN:'   '{:12.3f}'\n".format(welem.min()))
        fout.write("'LAMBDA_MAX:'   '{:12.3f}'\n".format(welem.max()))
        fout.write("'INTENSITY/FLUX:'  'Flux'\n")
        fout.write("'COS(THETA):'  '1.00'\n")
        fout.write("'ABFIND:'  '.false.'\n")
        fout.write("'MODELINPUT:'  '{:s}'\n".format(os.path.basename(atmod)))
        if kurucz : fout.write("'MARCS-FILE:'  '.false.'\n")
        fout.write("'MODELOPAC:'  '{:s}'\n".format(os.path.basename(root)+'opac'))
        fout.write("'RESULTFILE:'  '{:s}'\n".format(os.path.basename(file)))
        fout.write("'METALLICITY:'  '{:8.3f}'\n".format(mh))
        fout.write("'ALPHA/Fe:'  '{:8.3f}'\n".format(am))
        fout.write("'HELIUM:'  '{:8.3f}'\n".format(0.00))
        fout.write("'R-PROCESS:'  '{:8.3f}'\n".format(0.00))
        fout.write("'S-PROCESS:'  '{:8.3f}'\n".format(0.00))
        fout.write("'INDIVIDUAL ABUNDANCES:'  '{:2d}'\n".format(nels))
        if els is None or not 'C' in np.array(els)[:,0] :
            fout.write("    6  {:8.3f}\n".format(8.39+mh+cm))
        if els is None or not 'N' in np.array(els)[:,0] :
            fout.write("    7  {:8.3f}\n".format(7.78+mh+nm))
        if els is None or not 'O' in np.array(els)[:,0] :
            fout.write("    8  {:8.3f}\n".format(8.66+mh+am))
        if elemgrid != '' : fout.write("{:6d}  {:8.3f}\n".format(elemnum,abun+mh))
        if els is not None :
            for el in els :
              num=atomic.periodic(el[0])[0]
              ab=atomic.solar(el[0])[0]+el[1]
              fout.write("{:5d}  {:8.3f}\n".format(num,ab+mh))
        if  not solarisotopes :
          fout.write("'ISOTOPES:'  '2'\n")
          # adopt ratio of 12C/13C=15
          fout.write("   6.012 0.9375\n")
          fout.write("   6.013 0.0625\n")
        # do we need to add the H2O linelist?
        if h2o is None : 
          h2o=0
        elif h2o==0 :
          if teff < 4000 :
            if mh+am < -1.5 or teff > 3250 :
              h2o=1    
            else  :
              h2o=2
          else :
              h2o=0
        nlists=0
        # if no HI lines, don't use that list: it takes a while to read/process
        n_HI = len(open(linelistdir+'/turbospec.'+linelist+'.Hlinedata').readlines())
        if n_HI > 2 : nlists+=1
        if atoms : nlists+=1
        if molec : nlists+=1
        # if we are using H2O, add that list
        if h2o > 0 : nlists+=1
        fout.write("'NFILES:'  '{:4d}'\n".format(nlists))
        if n_HI >= 3 : fout.write(linelistdir+'/turbospec.'+linelist+'.Hlinedata\n')
        if atoms: fout.write(linelistdir+'/turbospec.'+linelist+'.atoms\n')
        if molec: fout.write(linelistdir+'/turbospec.'+linelist+'.molec\n')
        if h2o == 1 :
            fout.write(linelistdir+'/turbospec.h2o-BC8.5V'+'.molec\n')
        elif h2o == 2 :
            fout.write(linelistdir+'/turbospec.h2o-BC9.5V'+'.molec\n')
        if geo == 's' :
            fout.write("'SPHERICAL:'  'T'\n")
        else :
            fout.write("'SPHERICAL:'  'F'\n")
        fout.write("30\n")
        fout.write("300.00\n")
        fout.write("15\n")
        fout.write("1.3\n")
        fout.close()

        # control file, with special handling in case bsyn goes into infinite loop ...
        fout = open(root+"_bsyn.csh",'w')
        fout.write("#!/bin/csh -f\n")
        fout.write("{:s}/bin/bsyn_lu < {:s} &\n".format(os.environ['APOGEE_DIR'],os.path.basename(bsynfile)))
        fout.write('set bsynjob = $!\n')
        fout.write("set ok = 0\n")
        fout.write("set runtime = `ps -p $bsynjob -o cputime | tail -1 | awk -F: '{print ($1*3600)+($2*60)+$3}'`\n")
        tmax=120*int(.05/min([0.05,dw]))
        if h2o == 1 : tmax*=2
        if h2o == 2 : tmax*=5
        fout.write('while ( $runtime < {:d} && $ok == 0 )\n'.format(tmax))
        fout.write('  usleep 200000\n')
        fout.write("  set runtime = `ps -p $bsynjob -o cputime | tail -1 | awk -F: '{print ($1*3600)+($2*60)+$3}'`\n")
        fout.write('  if ( `ps -p $bsynjob -o comm=` == "" ) then\n')
        fout.write('    echo process done, exiting!\n')
        fout.write('    set ok = 1\n')
        fout.write('  endif\n')
        fout.write('end\n')
        fout.write('if ( $ok == 0 ) then\n')
        fout.write('  echo expired, killing job\n')
        fout.write('  kill $bsynjob\n')
        fout.write('endif\n')
        fout.close()
        if run :
            os.chmod(root+'_bsyn.csh', 0o777)
            os.chdir(workdir)
            subprocess.call(['time','./'+os.path.basename(root)+'_bsyn.csh'],stdout=stdout)
            os.chdir(cwd)
            try:
                out=np.loadtxt(file)
                if ielem == 0 : 
                    spec=out[:,2]
                    specnorm=out[:,1]
                else :
                    spec=np.vstack([spec,out[:,2]])
                    specnorm=np.vstack([specnorm,out[:,1]])
            except :
                print('failed...',file)
                fail('mkturbospec array problem: {:8d} {:8.2f} {:8.2f} {:8.2f}  {:8.2f} {:8.2f} {:8.2f}'.format(
                               teff,logg,mh,am,cm,nm,vmicro))
                return 0.,0.

    if run :
        if not save : shutil.rmtree(workdir)
        return spec, specnorm

def fail(out,file='FAILURE') :
    """ Routine to log to FAILURE file

    Args:
        out (str) : string to write to file
        file (str) : name of failure file (default='FAILURE')
    """
    ferr = open(file,'a+')
    ferr.write(out+'\n')
    ferr.close()

def prange(start,delta,n) :
    """ Routine to return vector of values given start, delta, n as strings

    Args :
        start (str) : starting vlue
        delta (str) : delta
        n (str) : number of values

    Returns :
        np.array : array of n values, from start incremented by delta
    """
    return float(start)+np.arange(int(n))*float(delta)

def get_vmicro(vmicrofit,vmicro,teff=None,logg=None,mh=None) :
    """ Return vmicro from input functional type and coefficients
   
    Args :
        vmicrofit (int) : input vmicro code to set fitting function
        vmicro (float)  : input vmicro fitting function coefficients

    Returns:
        vmicro (float) : microturbulent velocity

    """ 
    if vmicrofit == 0 :
        vm = vmicro[0]
    elif vmicrofit == 1 :
        #cubic in logg
        vm = 10.**(vmicro[0]+vmicro[1]*logg+vmicro[2]*logg**2+vmicro[3]*logg**3)
    elif vmicrofit == 2 :
        #linear with terms in  teff, logg, mh
        vm = 10.**(vmicro[0]+vmicro[1]*teff+vmicro[2]*logg+vmicro[3]*mh)
    else :
        print('need to implement vmicrofit: ', vmicrofit)
        pdb.set_trace()
    return vm

def mkgrid(planfile,code=None,clobber=False,save=False,run=True,atoms=True,molec=True) :
    """ Create a grid of synthetic spectra using Turbospectrum given specifications in  input parameter file
        Outputs results in FITS file

    Args :
        planfile (str) : name of input Yanny planfile with grid specifications
        clobber ( bool ) : skip trying to use previously saved syntheses (by [M/H]) (default = False)
        save (bool) :  save all temporary files for syntheses (default = False)
        run (bool)  : actually run syntheses (default = True)
    """

    # Read planfile
    if not os.path.isfile(planfile): 
        print('{:s} does not exist'.format(planfile))
        return
    p=yaml.safe_load(open(planfile,'r'))

    # grid specifications from input planfile
    if code == None : code = p['synthcode'] if p.get('synthcode') else None
    wrange=p['wrange']
    dw=float(p['dw'])
    h2o = int(p['h2o']) if p.get('h2o') else 0
    vacuum = int(p['vacuum']) if p.get('vacuum') else 0
    kurucz = True if p['atmos'] == 'kurucz' else False
    marcsdir = p['marcsdir'] if p.get('marcsdir') else None
    solarisotopes = int(p['solarisotopes']) if p.get('solarisotopes') else 0
    solarisotopes = True if abs(solarisotopes) == 1 else False
    nlte = p['nlte'] if p.get('nlte') else False
    cmnear = p['cmnear'] if p.get('cmnear') else False
    enhanced_o = p['enhanced_o'] if p.get('enhanced_o') else 0
    elem = p['elem'] if p.get('elem') else ''
    maskdir = p['maskdir'] if p.get('maskdir') else None
    vmicrofit = int(p['vmicrofit']) if p.get('vmicrofit') else 0
    vmicro = np.array(p['vmicro']) if p.get('vmicro') else 0.
    vmacrofit = int(p['vmacrofit']) if p.get('vmacrofit') else 0
    vmacro = p['vmacro'] if p.get('vmacro') else 0
    try: spherical = p['spherical'] 
    except KeyError : spherical = True
    specdir = os.environ['APOGEE_SPECLIB']+'/synth/'+p['specdir'] if p.get('specdir') else './'
    linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/' 
    linelist = p['linelist'] if p.get('linelist') else None
    oa0 = p['oa0'] if p.get('oa0') else 0.
    doa = p['doa'] if p.get('doa') else 0.
    noa = p['noa'] if p.get('noa') else 1

    # wavelength array
    nspec=int((wrange[1]-wrange[0])/dw)+1
    rawwave=wrange[0]+np.arange(nspec)*dw

    # if element minigrid, create mini linelist and get wavelengths to store
    if elem == '' :
        nelem=1
        gdspec=range(nspec)
        nwind = 1
        pixels=[[0,nspec]]
    else :
        # number of minigrid abundances
        nelem=8
        wvac,wair = mini_linelist(elem,linelist,maskdir=maskdir)
        nwind=wair.shape[0]
        gdspec=[]
        pixels=[]
        for iwind in range(nwind) :
            pix=np.where( (rawwave >= wair[iwind,0]) & (rawwave <= wair[iwind,1]) )[0]
            gdspec.extend(pix)
            pixels.append([pix[0],pix[-1]+1])
        nspec=len(gdspec)
        # for elem with all waves, use next line and comment out previous 8
        #gdspec=range(nspec)
        #nwind = 1
        #pixels=[[0,nspec]]

    # make the grid(s)
    for oa in prange(oa0,doa,noa) :
     for am in prange(p['am0'],p['dam'],p['nam']) :
      if enhanced_o : oa = [('O',2*am)]
      else : oa = [('O',oa+am)]
      for cm in prange(p['cm0'],p['dcm'],p['ncm']) :
        for nm in prange(p['nm0'],p['dnm'],p['nnm']) :
          specdata=np.zeros([nelem,int(p['nmh']),int(p['nlogg']),int(p['nteff']),nspec],dtype=np.float32)
          #specnormdata=np.zeros([nelem,int(p['nmh']),int(p['nlogg']),int(p['nteff']),nspec],dtype=np.int16)
          specnormdata=np.zeros([nelem,int(p['nmh']),int(p['nlogg']),int(p['nteff']),nspec],dtype=np.float32)
          # allow for restart after certain number of [M/H] have been completed
          if clobber :
              nmh = 0
          else :
              try :
                  # does output file exist?
                  old=fits.open(specdir+'/'+p['name']+elem+'.fits')
                  specdata=old[0].data
                  specnormdata=old[1].data
                  if len(old[0].shape) < 5 : specdata=np.expand_dims(specdata,axis=0)
                  if len(old[1].shape) < 5 : specnormdata=np.expand_dims(specnormdata,axis=0)
                  # is it a partially completed file with nmh card, or a completed file?
                  try:
                      nmh =old[0].header['nmh']
                  except :
                      # file is completed
                      nmh = p['nmh']
                      print('file already done!')
                      return
              except :
                  nmh = 0
          for imh,mh in enumerate(prange(p['mh0'],p['dmh'],p['nmh'])) :
           if imh >= nmh :
            for ilogg,logg in enumerate(prange(p['logg0'],p['dlogg'],p['nlogg'])) :
              for iteff,teff in enumerate(prange(p['teff0'],p['dteff'],p['nteff'])) :

                print(teff, logg, mh)
                sys.stdout.flush()
                nskip=0 
                dskip = 1 if kurucz else 2
                vout = get_vmicro(vmicrofit,vmicro,teff=teff,logg=logg,mh=mh)
                print(teff,logg,mh,am,cm,nm,oa,vout)
                while nskip >= 0 and nskip < 10 :
                  if code is None :
                      spec,specnorm=mkturbospec(int(teff),logg,mh,am,cm,nm,els=oa,
                        wrange=wrange,dw=dw,atmosdir=marcsdir,
                        elemgrid=elem,linelistdir=linelistdir+'/'+elem+'/',linelist=linelist,vmicro=vout,
                        solarisotopes=solarisotopes,
                        nskip=nskip,kurucz=kurucz,run=run,save=save) 
                  else :
                      spec,specnorm=mk_synthesis(code,int(teff),logg,mh,am,cm,nm,els=oa,cmnear=cmnear,
                        wrange=wrange,dw=dw,atmosdir=marcsdir,
                        elemgrid=elem,linelistdir=linelistdir+'/'+elem+'/',linelist=linelist,vmicro=vout,
                        solarisotopes=solarisotopes,nlte=nlte,dospherical=spherical,
                        nskip=nskip,atmos_type=p['atmos'],run=run,save=save,atoms=atoms,molec=molec,h2o=h2o) 
                  nskip = nskip+dskip if isinstance(spec,float) else -1
                if nskip > 0 : 
                    print('FAILED synthesis',nskip)
                    fail('failed synthesis: {:8d} {:8.2f} {:8.2f} {:8.2f}  {:8.2f} {:8.2f} {:8.2f} {:d}'.format(
                               int(teff),logg,mh,am,cm,nm,vout,nskip))
                try:
                    if elem == '' :
                        specdata[0,imh,ilogg,iteff,:]=spec
                        #specnormdata[0,imh,ilogg,iteff,:]=np.round((specnorm-0.5)*65534.).astype(int)
                        specnormdata[0,imh,ilogg,iteff,:]=specnorm
                    else :
                        specdata[:,imh,ilogg,iteff,:]=spec[:,gdspec]
                        #specnormdata[:,imh,ilogg,iteff,:]=np.round((specnorm[:,gdspec]-0.5)*65534).astype(int)
                        specnormdata[:,imh,ilogg,iteff,:]=specnorm
                except :
                    print(specdata.shape)
                    specdata[:,imh,ilogg,iteff,:]=0.
                    specnormdata[:,imh,ilogg,iteff,:]=-32767
                    try: lspec=len(spec)
                    except: lspec=1
                    fail('error loading specdata: {:8d} {:8.2f} {:8.2f} {:8.2f}  {:8.2f} {:8.2f} {:8.2f} {:d}'.format(
                               int(teff),logg,mh,am,cm,nm,vout,lspec))

            # FITS header and output after each metallicity subgrid
            # for minigrids, each section is output in a separate HDU
            p1=0
            hdulist=fits.HDUList()
            for iwind in range(nwind) :
                p2 = p1 + pixels[iwind][1]-pixels[iwind][0]
                print(iwind,p1,p2)
                if iwind == 0 :
                    hdu=fits.PrimaryHDU(np.squeeze(specdata[:,:,:,:,p1:p2]))
                else :
                    hdu=fits.ImageHDU(np.squeeze(specdata[:,:,:,:,p1:p2]))
                p1 += pixels[iwind][1]-pixels[iwind][0]
                idim=1
                # for elem with all waves, use next line and comment out following 7
                #spectra.add_dim(hdu.header,rawwave[0],rawwave[1]-rawwave[0],1,'WAVELENGTH',idim)
                if elem == '' :
                    spectra.add_dim(hdu.header,rawwave[0],rawwave[1]-rawwave[0],1,'WAVELENGTH',idim)
                else :
                    hdu.header.append(('ELEM',elem))
                    gd = np.where( (rawwave >= wair[iwind,0]) & (rawwave <= wair[iwind,1]) )[0]
                    spectra.add_dim(hdu.header,rawwave[gd[0]],rawwave[1]-rawwave[0],1,'WAVELENGTH',idim)
                    hdu.header.append(('CDELT1',rawwave[1]-rawwave[0]))
                if int(p['nteff']) > 1 :
                    idim+=1
                    spectra.add_dim(hdu.header,float(p['teff0']),float(p['dteff']),1,'TEFF',idim)
                if int(p['nlogg']) > 1 :
                    idim+=1
                    spectra.add_dim(hdu.header,float(p['logg0']),float(p['dlogg']),1,'LOGG',idim)
                if int(p['nmh']) > 1 :
                    idim+=1
                    spectra.add_dim(hdu.header,float(p['mh0']),float(p['dmh']),1,'M_H',idim)
                if elem != '' :
                    idim+=1
                    spectra.add_dim(hdu.header,-0.75,0.25,1,elem,idim)
                hdu.header['LOGW'] = 0
                if p.get('width') : hdu.header['width'] = p['width']
                if p.get('linelist') : hdu.header['linelist'] = p['linelist']
                if p.get('solarisotopes') : hdu.header['isotopes'] = int(p['solarisotopes'])
                if p['synthcode'] == 'asset'  : hdu.header.add_comment('ASSET generated synthetic spectra')
                if p['synthcode'] == 'turbospec' : hdu.header.add_comment('Turbospec generated synthetic spectra')
                if p['synthcode'] == 'moog ' : hdu.header.add_comment('MOOG generated synthetic spectra')
                hdu.header.add_comment('APOGEE_VER:'+os.environ['APOGEE_VER'])
                # following card for partial completion output
                if imh+1 < int(p['nmh']) : hdu.header['nmh'] = imh+1
                try : os.mkdir(specdir)
                except: pass
                hdulist.append(hdu)
            hdunorm=fits.ImageHDU(np.squeeze(specnormdata))
            hdunorm.header.extend(hdu.header.copy(strip=True))
            #hdunorm.header['BZERO'] = 0.5
            #hdunorm.header['BSCALE'] = 1./65534.
            hdulist.append(hdunorm)
            hdulist.writeto(specdir+'/'+p['name']+elem+'.fits',overwrite=True)

          return hdulist

def mkgridlink(planfile,suffix=None) :
    """  DEVELOPMENT : create coarse grid by merging syntheses from multiple grids
    """

    # Read planfile
    if not os.path.isfile(planfile): 
        print('{:s} does not exist'.format(planfile))
        return
    p=yanny.yanny(planfile,np=True)

    linelist=p['linelist'][2:]
    if suffix is None :
        suffix='_'+p['smooth']

    for ivm,vm in enumerate(spectra.vector(p['vt0'],p['dvt'],p['nvt'])) :
      for icm,cm in enumerate(spectra.vector(p['cm0'],p['dcm'],p['ncm'])) :
        for inm,nm in enumerate(spectra.vector(p['nm0'],p['dnm'],p['nnm'])) :
          for iam,am in enumerate(spectra.vector(p['am0'],p['dam'],p['nam'])) :
            file=('a{:s}c{:s}n{:s}v{:s}.fits').format(
                   atmos.cval(am),atmos.cval(cm),atmos.cval(nm),atmos.cval(10**vm))
            print(file)

            GKg= fits.open('../../giantisotopes/tgGK_'+linelist+suffix+'/rbf_'+file)[0]
            Mg= fits.open('../../giantisotopes/tgM_'+linelist+suffix+'/rbf_'+file)[0]
            Fg= fits.open('../../giantisotopes/tgF_'+linelist+suffix+'/'+file)[0]
            Fd= fits.open('../../solarisotopes/tdF_'+linelist+suffix+'/rbf_'+file)[0]
            GKd= fits.open('../../solarisotopes/tdGK_'+linelist+suffix+'/rbf_'+file)[0]
            Md= fits.open('../../solarisotopes/tdM_'+linelist+suffix+'/rbf_'+file)[0]
            nwave = GKg.shape[-1]
            s=np.zeros([int(p['nmh']),int(p['nlogg']),int(p['nteff']),nwave],dtype=np.float32)
            grids = [Fg, GKg, Mg, Fd, GKd, Md]
            for imh,mh in enumerate(spectra.vector(p['mh0'],p['dmh'],p['nmh'])) :
              for ilogg,logg in enumerate(spectra.vector(p['logg0'],p['dlogg'],p['nlogg'])) :
                for iteff,teff in enumerate(spectra.vector(p['teff0'],p['dteff'],p['nteff'])) :
                    igrid=0
                    i,j,k=getindex(Fg.header,(2,3,4),(teff,logg,mh))
                    if i<0 or j<0 or k<0 :
                        igrid=1
                        i,j,k=getindex(GKg.header,(2,3,4),(teff,logg,mh))
                    if i<0 or j<0 or k<0 :
                        igrid=2
                        i,j,k=getindex(Mg.header,(2,3,4),(teff,logg,mh))
                    if i<0 or j<0 or k<0 :
                        igrid=3
                        i,j,k=getindex(Fd.header,(2,3,4),(teff,logg,mh))
                    if i<0 or j<0 or k<0 :
                        igrid=4
                        i,j,k=getindex(GKd.header,(2,3,4),(teff,logg,mh))
                    if i<0 or j<0 or k<0 :
                        igrid=5
                        i,j,k=getindex(Md.header,(2,3,4),(teff,logg,mh))
                    if i<0 or j<0 or k<0 :
                        print("can't find model to fill!",mh,logg,teff)
                    else :
                        #print(mh,logg,teff,igrid,i,j,k)
                        try:
                            #print(grids[igrid].data.shape)
                            s[imh,ilogg,iteff,:]=grids[igrid].data[k,j,i,:]
                        except:
                            pdb.set_trace()
            hdulist=fits.HDUList()
            hdu=fits.ImageHDU(s)
            idim=1
            hdu.header['CRVAL1'] = GKg.header['CRVAL1']
            hdu.header['CDELT1'] = GKg.header['CDELT1']
            hdu.header['CTYPE1'] = GKg.header['CTYPE1']
            if int(p['nteff']) > 1 :
                idim+=1
                spectra.add_dim(hdu.header,float(p['teff0']),float(p['dteff']),1,'TEFF',idim)
            if int(p['nlogg']) > 1 :
                idim+=1
                spectra.add_dim(hdu.header,float(p['logg0']),float(p['dlogg']),1,'LOGG',idim)
            if int(p['nmh']) > 1 :
                idim+=1
                spectra.add_dim(hdu.header,float(p['mh0']),float(p['dmh']),1,'M_H',idim)
            hdulist.append(hdu)
            hdulist.writeto(file,overwrite=True)

def comp(testdir,libdir,name='ap00cp00np00vp12.fits',thresh=0.005,yr=None,pause=True,outdir='./',mh=None,names=None) :
    """ compare test grid with library grid
    """
    test = fits.open(testdir+'/'+name)[1]
    lib = fits.open(libdir+'/'+name)[1]
    wave=spectra.fits2vector(test.header,1)
    nmh,nlogg,nte,nwave = test.shape

    sfig,sax=plots.multi(nte,nlogg,hspace=0.001,wspace=0.001,figsize=(10,4))
    sfig.suptitle(testdir+' / '+libdir)
    fig,ax=plots.multi(nte,nlogg,hspace=0.001,wspace=0.001,figsize=(10,4))
    fig.suptitle(testdir+' / '+libdir)
    color=['m','b','g','r']
    if mh is None : 
        mhlist =  spectra.fits2vector(test.header,4)
        mhlab= False
    else : 
        mhlist = [mh]
        mhlab= True
        
    for ite,te in enumerate(spectra.fits2vector(test.header,2)) :
        for ilogg,logg in enumerate(spectra.fits2vector(test.header,3)) :
            if mhlab :
                sax[ilogg,ite].text(0.1,0.9,'Te: {:.0f} logg: {:6.2f} [M/H]: {:5.1f}'.format(te,logg,mh),
                               transform=ax[ilogg,ite].transAxes)
                ax[ilogg,ite].text(0.1,0.9,'Te: {:.0f} logg: {:6.2f} [M/H]: {:5.1f}'.format(te,logg,mh),
                               transform=ax[ilogg,ite].transAxes)
            else :
                sax[ilogg,ite].text(0.1,0.9,'Te: {:.0f} logg: {:6.2f}'.format(te,logg),
                               transform=ax[ilogg,ite].transAxes)
                ax[ilogg,ite].text(0.1,0.9,'Te: {:.0f} logg: {:6.2f}'.format(te,logg),
                               transform=ax[ilogg,ite].transAxes)
            for mh in reversed(mhlist) :
                ind=getindex(test.header,[2,3,4],[te,logg,mh])
                imh=ind[2]
                print('test',mhlab)
                print(te,logg,mh)
                print(ind)
                ind=getindex(lib.header,[2,3,4],[te,logg,mh])
                print('lib')
                print(ind)
                print('')
                ax[ilogg,ite].plot(wave,test.data[imh,ilogg,ite,:]-lib.data[ind[2],ind[1],ind[0],:],
                    color=color[imh],label='{:.2}'.format(mh))
                plt.draw()
                if names is not None : label=names[1]
                else : label='{:.2}'.format(mh)
                sax[ilogg,ite].plot(wave,lib.data[ind[2],ind[1],ind[0],:],label=label)
                if names is not None : label=names[0]
                else : label='{:.2}'.format(mh)
                sax[ilogg,ite].plot(wave,test.data[imh,ilogg,ite,:],label=label)
                sax[ilogg,ite].set_ylim([0,1.2])
                if yr is not None : ax[ilogg,ite].set_ylim(yr)
                plt.draw()
                if np.abs(test.data[imh,ilogg,ite,:]-lib.data[ind[2],ind[1],ind[0],:]).max() > thresh:
                    if pause: pdb.set_trace()
    ax[-1,-1].legend(fontsize='xx-small')
    sax[-1,-1].legend(fontsize='xx-small')
    if mhlab :
      outfile = os.path.basename(testdir)+'_'+os.path.basename(libdir)+'_{:.1f}'.format(mh)
    else :
      outfile = os.path.basename(testdir)+'_'+os.path.basename(libdir)
    fig.savefig(outdir+'/'+outfile+'.png')
    sfig.savefig(outdir+'/'+outfile+'_spec.png')
    plt.close()
    plt.close()
    return [outfile+'_spec.png',outfile+'.png']

def allcomp() :
    """ create grid of synthesis comparisons, with web page
    """

    matplotlib.use('TkAgg')
    #grids = ['solarisotopes/testMd']
    #tab=[]
    #for grid in grids :
    #    try: out=comp('synspec/marcs/'+grid+'_h2o2','synspec/marcs/'+grid+'_h2o3',outdir='comp',pause=False,yr=[0.98,1.02])
    #    except : out=['','']
    #    tab.append(out)
    #    try: out=comp('synspec/marcs/'+grid+'_h2o4','synspec/marcs/'+grid+'_h2o3',outdir='comp',pause=False,yr=[0.98,1.02])
    #    except : out=['','']
    #    tab.append(out)
    #    try: out=comp('synspec/marcs/'+grid+'_h2o1','synspec/marcs/'+grid+'_h2o3',outdir='comp',pause=False,yr=[0.98,1.02])
    #    except : out=['','']
    #    tab.append(out)
    #    try: out=comp('synspec/marcs/'+grid+'_noh2o','synspec/marcs/'+grid+'_h2o3',outdir='comp',pause=False,yr=[0.98,1.02])
    #    except : out=['','']
    #    tab.append(out)
    #    try: out=comp('synspec/marcs/'+grid+'_h2o0','synspec/marcs/'+grid+'_h2o3',outdir='comp',pause=False,yr=[0.98,1.02])
    #    except : out=['','']
    #    tab.append(out)
    #yt=['-9.5EV','-9.0EV','-8.5EV','No H2O','DR17 strawman']
    #html.htmltab(tab,file='comp/comph2o.html',ytitle=yt)
    grids = ['solarisotopes/testFd','solarisotopes/testGKd','solarisotopes/testMd',
             'giantisotopes/testGKg','giantisotopes/testMg']
    libs = ['solarisotopes/tdF_180901','solarisotopes/tdGK_180901','solarisotopes/tdM_180901',
             'giantisotopes/tgGK_180901','giantisotopes/tgM_180901']
    grids = ['giantisotopes/testMg']
    libs = ['giantisotopes/tgM_180901']
    tab=[]
    ytit=[]
    for grid,lib in zip(grids,libs):
      for mh in [-2.5,-1.5,-0.5,0.5] :
        try: out=comp('synspec/marcs/'+grid,'turbospec/marcs/'+lib,outdir='comp',pause=False,yr=[-.5,0.5],mh=mh,names=['syn','turbo'])
        except KeyboardInterrupt: 
            raise
        except: out=['','']
        tab.append(out)
        ytit.append('{:s} {:6.1f}'.format(grid,mh))
    html.htmltab(tab,file='comp/comp_syn_turbo.html',ytitle=ytit)

def getindex(header,axes,vals) :
    """ DEVELOPMENT : get index of requested model from input grid header

    """
    out=[]
    for ax,val in zip(axes,vals) :
        i = (val-header['CRVAL'+str(ax)])/header['CDELT'+str(ax)]
        if i+1 > header['NAXIS'+str(ax)] : i=-1.
        out.append(int(round(i)))
    return out

def mkgridlsf(planfile,highres=9,fiber=None,ls=None,apred=None,prefix=None,telescope=None) :
    """ Create a grid of LSF-convolved spectra given specifications in input parameter file and existing raw syntheses

    Args :
        planfile (str) : name of Yanny parameter file with specifications
        highres (int ) : number of subpixels to use for convolution (default=9)
        fiber (int or None ) : fiber number for LSF, if None use parameter file (default=None)
        ls () : lsf convolution kernel (?), if precomputed (default=None)
        apred (str) : name of reduction version to get LSF and WAVE from (default='r8')
        prefix (str) : string to prepend for input syntheses file, e.g. 'rbf_' to use rbf-filled grids (default='')
    """

    # Read planfile
    if not os.path.isfile(planfile): 
        print('{:s} does not exist'.format(planfile))
        return
    p=yaml.safe_load(open(planfile,'r'))

    # header information
    specdir = os.environ['APOGEE_SPECLIB']+'/synth/'+p['specdir'] if p.get('specdir') else './'
    if fiber is None : fiber=p['lsffiber']
    if isinstance(fiber,int): fiber= [fiber]
    if apred is None :apred = p['apred'] if p.get('apred') else 'r12'
    if telescope is None : telescope = p['telescope'] if p.get('telescope') else 'apo25m'
    lsfid=int(p.get('lsfid'))
    waveid=int(p.get('waveid'))
    vmacrofit = int(p['vmacrofit']) if p.get('vmacrofit') else 0
    vmacro_arr=p['vmacro']
    kernel=p['kernel'] if p.get('kernel') else 'rot'

    if ls is None :
        x, ls = getlsf(lsfid,waveid,apred=apred,telescope=telescope,fiber=fiber,highres=highres)
    else :
        x = ls[0]
        ls = ls[1]

    # read in raw spectra, from RBF interpolated if we have it
    if prefix is None :  
        if p.get('r0') and float(p['r0']) >= -0.001 : prefix = 'rbf_'
        else : prefix=''

    # get the synthesis
    # for regular grids, we just want to process the first extension,
    # for minigrids, we want to process all except the last (which is the normalized synthesis)
    print(specdir+'/'+prefix+p['name']+'.fits')
    speclist = fits.open(specdir+'/'+prefix+p['name']+'.fits')
    nexten = len(speclist) - 1

    hdulist=fits.HDUList()
    for exten in range(nexten) :

        specdata = speclist[exten]
        npix = specdata.data.shape[-1]
        nspec=1
        for i in range(len(specdata.data.shape)-1) :
            nspec*=specdata.data.shape[i]
        print('nspec: ', nspec)
        specdata.data=np.reshape(specdata.data,(nspec,npix))

        # synthesis is in air, we want vacuum
        ws=spectra.fits2vector(specdata.header,1)
        ws=spectra.airtovac(ws)

    # output wavelength grid
    #wa=aspcap.apStarWave()
    #nout=wa.shape[0]

        # is this a minigrid?
        try :
            if p['elem'] == '' : nelem = 1
            else : nelem=8
        except : nelem=1

        # create vmacro array
        vmacro=[]
        #dlam=np.log10(wa[1])-np.log10(wa[0])
        for l in range(nelem) :
          for k,mh in enumerate(prange(p['mh0'],p['dmh'],p['nmh'])) :
            for j,logg in enumerate(prange(p['logg0'],p['dlogg'],p['nlogg'])) :
              for i,teff in enumerate(prange(p['teff0'],p['dteff'],p['nteff'])) :
                if vmacrofit == 1 :
                    vm = 10.**(0.470794-0.254*mh)
                    vm = 10.**(vmacro_arr[0]+vmacro_arr[1]*teff+vmacro_arr[2]*logg+vmacro_arr[3]*mh)
                    vm = vm if vm<15 else 15.
                elif vmacrofit == 0 :
                    vm = 0.
                vmacro.append(vm)
        vmacro=np.array(vmacro)

        # LSF and rotation convolution all spectra at the same time
        nmh=int(p['nmh'])
        nlogg=int(p['nlogg'])
        nteff=int(p['nteff'])
        nrot = int(p['nrot'])

        if nrot == 1 :
            vrot=10.**.176
            smoothdata,waveout=lsf.convolve(ws,specdata.data,lsf=ls,xlsf=x,vrot=vrot,vmacro=vmacro)
            nout=smoothdata.shape[-1]
            smoothdata=np.reshape(smoothdata,(nelem,nmh,nlogg,nteff,nout)).astype(np.float32)
        else :
            for irot,vrot in enumerate(prange(p['rot0'],p['drot'],p['nrot'])) :
                if kernel == 'rot' :
                    smooth,waveout=lsf.convolve(ws,specdata.data,lsf=ls,xlsf=x,vrot=10.**vrot,vmacro=vmacro)
                elif kernel == 'gauss' :
                    smooth,waveout=lsf.convolve(ws,specdata.data,lsf=ls,xlsf=x,vmacro=10.**vrot)
                else :
                    print('Unknown kernel!')
                    pdb.set_trace()
                nout=smooth.shape[-1]
                if irot == 0 : smoothdata=np.zeros([nrot,nelem,nmh,nlogg,nteff,nout],dtype=np.float32)
                smoothdata[irot,:,:,:,:,:]=np.reshape(smooth,(nelem,nmh,nlogg,nteff,nout)).astype(np.float32)

        specdata.data=np.reshape(specdata.data,(nelem,nmh,nlogg,nteff,npix))
        if exten == 0 : hdu=fits.PrimaryHDU(np.squeeze(smoothdata))
        else : hdu=fits.ImageHDU(np.squeeze(smoothdata))
        hdu.header.extend(specdata.header.copy(strip=True))
        #hdu.header['CRVAL1'] = aspcap.logw0
        hdu.header['CRVAL1'] = np.log10(waveout[0])
        hdu.header['CDELT1'] = aspcap.dlogw
        hdu.header['CTYPE1'] = 'LOG(WAVELENGTH)'
        if nrot > 1 :
            hdu.header.insert('CTYPE4',('CRVAL5',float(p['rot0']),''),after=True)
            hdu.header.insert('CRVAL5',('CDELT5',float(p['drot']),''),after=True)
            hdu.header.insert('CDELT5',('CRPIX5',1,''),after=True)
            hdu.header.insert('CRPIX5',('CTYPE5','LOG(VSINI)',''),after=True)
        hdu.header['INFILE'] = p['specdir']+'/'+prefix+p['name']+'.fits'
        hdu.header['APRED'] = apred
        hdu.header['LSFID'] = lsfid
        hdu.header['WAVEID'] = waveid
        hdu.header['HIGHRES'] = highres
        hdu.header.add_comment('LSF convolved spectra')
        hdu.header.add_comment('APOGEE_VER:'+os.environ['APOGEE_VER'])
        hdulist.append(hdu)

    hdulist.writeto(p['name']+'.fits',overwrite=True)

    return smoothdata

def complsf(name) :
    """ TEST : routine to compare calculated LSFs between Python and IDL
    """
    new = fits.open(name+'.fits')[0]
    a = fits.open('old_'+name+'.fits')[1].data
    b = fits.open('old_'+name+'.fits')[2].data
    c = fits.open('old_'+name+'.fits')[3].data

    for k in range(new.header['NAXIS4']) :
      for j in range(new.header['NAXIS3']) :
        for i in range(new.header['NAXIS2']) :
          old = np.append(a[k,j,i,:],b[k,j,i,:])
          old = np.append(old,c[k,j,i,:])

          plt.clf()
          plt.plot(aspcap.aspcap2apStar(old))
          plt.plot(new.data[k,j,i,:]/600000.)
          plt.plot(aspcap.aspcap2apStar(old)/new.data[k,j,i,:]*600000.)
          plt.draw()
          pdb.set_trace()

def mkspec(input) :
    """ Makes a single spectrum given input pars
        Used by mksynth for multi-processor calculations

    Args :
        pars (list ) : list of spectrum parameters: [Teff, logg, [M/H], [alpha/M], [C/M], [N/M], vmicro, vrot, 21 abundances]

    Returns :
        pars : input parameters
        spec : synthetic spectrum
    """
    pars=input[0]
    indata=input[1]

    teff=pars[0].astype('int')
    logg=pars[1]
    mh=pars[2]
    if indata['interp'] == 'interp' :
        am=pars[3]
        cm=pars[4]
    else :
        am=round(pars[3]/0.25)*0.25
        cm=round(pars[4]/0.25)*0.25
    #nm=round(pars[5]/0.5)*0.5
    nm=pars[5]
    vmicro=pars[6]
    vrot=pars[7]
    elems=[]
    els = ['O','Na','Mg','Al','Si','P','S','K','Ca','Ti','V','Cr','Mn','Co','Fe','Ni','Cu','Ge','Rb','Ce','Nd']
    for j,el in enumerate(els) :
        elems.append([el,pars[8+j]])
    print(('{:d}'+6*'{:8.3f}').format(teff,logg,mh,vmicro,am,cm,nm))
    spec,specnorm=mk_synthesis('turbospec',teff,logg,mh,am,cm,nm,vmicro=vmicro,els=elems,atmos_type=indata['atmos_type'],fill=False,
                              linelist=indata['linelist'],linelistdir=indata['linelistdir'],atmod=indata['interp'],
                              wrange=indata['wrange'],dw=indata['dw'],h2o=indata['h2o'],atoms=indata['atoms'],save=True,verbose=False)
    return pars,specnorm
    

def mksynth(inpars,havefits=False,threads=8,highres=9,lsfid=14600018,waveid=13140000,apred='r12',telescope='apo25m',
            fiber='combo',linelist='20180901',linelistdir=None,atmos_type='marcs',interp=None,h2o=0,atoms=True,plot=False,lines=None,ls=None) :
    """ Make a series of spectra from parameters in an input file, with parallel processing for turbospec
        Outputs to FITS file {file}.fits

    Args:
        inpars (str) : name of input file with parameters of spectra to calculate
        threads (int) : number of parallel processes
        highres (int) : number of subpixels for LSF convolution (default=9)
        waveid (int) : ID for wavelength calibration file (default=14600018)
        lsfid (int) : ID for LSF calibration file (default=13140000)
        apred (str) : reduction version for waveid,lsfid (default='r12')
        fiber (int or str) : fiber for LSF (default='combo')
        plot (bool) : plot each spectrum (default=False)
        lines (list) : specify limited range of input lines to calculate (default=None, i.e. all lines)
    """
    if havefits :
        # have model parameters in a FITS table
        tab = Table.read(inpars)
        names = pars.colnames
        pars=[]
        for name in names :
            pars.append(tab[name])
        pars = np.array(pars).T
       
    elif type(inpars) == np.ndarray :
        if len(inpars.shape) == 1 : pars = np.array([inpars])
        else : pars = inpars
 
    else :
        # have model parameters in an ASCII table
        names=ascii.read(inpars).colnames
        pars=np.loadtxt(inpars)
        if lines is not None :
            pars = pars[lines[0]:min([len(pars),lines[1]])]
            suffix = '_{:d}'.format(lines[0])
        else :
            suffix = ''

    indata={}
    indata['linelist']=linelist
    indata['linelistdir']=linelistdir
    indata['atmos_type']=atmos_type
    indata['interp']=interp
    indata['h2o']=h2o
    indata['atoms']=atoms
    indata['wrange']=[15100.,17000.]
    indata['dw']=0.05
    nspec=38001
    inputs=[]
    for par in pars :
        inputs.append((par,indata))

    if threads == 0 :
        specs=[]
        for input in inputs :
            specs.append(mkspec(input))
    else :
        pool = mp.Pool(threads)
        specs = pool.map_async(mkspec, inputs).get()
        pool.close()
        pool.join()

    # convolved and bundle output spectra into output fits file
    wa=aspcap.apStarWave()
    if lsfid is not None :
        prefix='lsf_'
        if ls is None :
            x, ls = getlsf(lsfid,waveid,prefix=prefix,apred=apred,telescope=telescope,fiber=fiber,highres=highres)
        else :
            x=ls[0]
            ls=ls[1]

    out=[]
    conv=[]
    outpar=[]
    if plot : plt.clf()
    ws=np.linspace(15100.,17000., nspec)
    # synthesis is in air, we want vacuum
    ws=spectra.airtovac(ws)
    print('convolving...')
    for spec in specs :
        if isinstance(spec[1],np.ndarray) :
          if spec[1].sum() > 0.001 :
            mh=spec[0][2]
            vmacro = 10.**(0.470794-0.254*mh)
            vmacro = vmacro if vmacro<15 else 15.
            vrot=spec[0][7]
            if vrot < 0.5 : vrot=None
            if lsfid is not None :
                # convolve one at a time because we have different vrot for each
                z,waveout=lsf.convolve(ws,spec[1],lsf=ls,xlsf=x,vmacro=vmacro,vrot=vrot)
                conv.append(np.squeeze(z))
            out.append(spec[1])
            outpar.append(spec[0])
            if plot : plt.plot(wa,np.squeeze(z))

    if type(inpars) == np.ndarray : return(conv)

    # write the spectra out
    hdu=fits.HDUList()
    h=fits.ImageHDU(outpar)
    h.header['NPAR' ] = len(names)
    for ipar,par in enumerate(names) : h.header['PAR{:d}'.format(ipar)] = par
    hdu.append(h)
    hdu.append(fits.ImageHDU(out))
    if lsfid is not None :
        h=fits.ImageHDU(conv)
        h.header['CRVAL1'] = np.log10(wa[0])
        h.header['CDELT1'] = np.log10(wa[1])-np.log10(wa[0])
        h.header['CTYPE1'] = 'LOG(WAVELENGTH)'
        hdu.append(h)
    hdu.writeto(file+suffix+'.fits',overwrite=True)
    return file+suffix+'.fits' 

def getlsf(lsfid,waveid,apred='r12',telescope='apo25m',highres=9,prefix='lsf_',fiber='combo',clobber=False,fill=False) :
    """ Create LSF FITS file or read if already created
    """
    lsfile = prefix+'{:08d}_{:08d}.fits'.format(lsfid,waveid)
    while os.path.isfile(lsfile+'.lock') :
        # if another process is creating LSF wait until done
        print('waiting for lock: ',lsfile+'.lock')
        time.sleep(10)

    if os.path.isfile(lsfile) and not clobber :
        # if file exists, read it
        x=fits.open(lsfile)[1].data
        ls=fits.open(lsfile)[2].data
    else :
        fp = open(lsfile+'.lock','w')
        fp.close()
        # lsf.get does the real work
        x,ls = lsf.get(lsfid,waveid,fiber,highres=highres,apred=apred,telescope=telescope)
        hdu=fits.HDUList()
        hdu.append(fits.PrimaryHDU())
        hdu[0].header['APRED'] = apred
        hdu[0].header['LSFID'] = lsfid
        hdu[0].header['WAVEID'] = waveid
        hdu[0].header['HIGHRES'] = highres
        for i,f in enumerate(fiber) :
            hdu[0].header['FIBER{:d}'.format(i)] = f
        hdu.append(fits.ImageHDU(x))
        hdu.append(fits.ImageHDU(ls))
        hdu.writeto(lsfile,overwrite=True)
        os.remove(lsfile+'.lock')

    if fill :
        # for all non-finite pixels, fill in LSF from nearest good pixel
        gd = np.where(np.isfinite(ls[:,0]))[0]
        mask = np.zeros(ls.shape[0],dtype=bool)
        mask[gd] = True
        bd = np.where(mask == False)[0]
        for i in bd:
            j = np.argmin(np.abs(i-gd))
            ls[i,:] = ls[gd[j],:]
    return x, ls

   
def elemsens(files=None,outfile='elemsens',highres=9,waveid=13140000,lsfid=14600018,apred='r12',telescope='apo25m',fiber='combo',
             calc=False,plot=True,ls=None,htmlfile='elemsens.html',filt=None,filtdir=None,linelist='20180901') :
    """ Create spectra at a range of parameters with individual abundances varied independently
    """
    # get ref.dat and {el}.dat with 0.1 adjustment
    files=sample.elemsens()
    # get ref.dat and {el}_0.dat with -8 adjustment
    files_0=sample.elemsens(delta=-8,suffix='_0')
    # get ref.dat and {el}_comp.dat with "complement" of abundances from elements in same group
    files_comp=sample.elemsens(suffix='_comp',complement=True)
    hdu=fits.HDUList()
    hdumask=fits.HDUList()
    if ls is None :
        xls, ls = getlsf(lsfid,waveid,apred=apred,telescope=telescope,fiber=fiber,highres=highres)
    else :
        xls = ls[0]
        ls = ls[1]
    grid=[]
    ytit=[]
    fp=html.head('elem.html')
    for i,(name,name0,name_comp) in enumerate(zip(files[3:],files_0[3:],files_comp[3:])) :
        elem=name.replace('.dat','')
        print(name)
        if name == 'C.dat' or name == 'N.dat' : continue
        if calc : 
            mini_linelist(elem,linelist,only=True,clobber=True)
            ## default abundance, lines of element only
            #out0=mksynth('ref.dat',threads=16,ls=(xls,ls),linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/'+elem+'_only',linelist=linelist,h2o=None)
            #os.rename('ref.dat.fits',elem+'_ref.fits')
            ## enhanced abundance, lines of element only
            #out1=mksynth(name,threads=16,ls=(xls,ls),linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/'+elem+'_only',linelist=linelist,h2o=None)
            #os.rename(name+'.fits',elem+'.fits')
            ## default abundance, all lines
            #out2=mksynth('ref.dat',threads=16,ls=(xls,ls),linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/',linelist=linelist,h2o=None)
            #os.rename('ref.dat.fits',elem+'_ref_full.fits')
            ## enhanced abundance, all lines
            #out3=mksynth(name,threads=16,ls=(xls,ls),linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/',linelist=linelist,h2o=None)
            #os.rename(name+'.fits',elem+'_full.fits')
            ## "zero" abundance, all lines
            #out4=mksynth(name0,threads=16,ls=(xls,ls),linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/',linelist=linelist,h2o=None)
            #os.rename(name0+'.fits',elem+'_0_full.fits')
            # "complement" abundance, all lines
            out5=mksynth(name_comp,threads=16,ls=(xls,ls),linelistdir=os.environ['APOGEE_SPECLIB']+'/linelists/',linelist=linelist,h2o=None)
            os.rename(name_comp+'.fits',elem+'_comp_full.fits')
        #ehdu=fits.open(out)[2]
        ehdu=fits.open(elem+'.fits')[2]
        ref=fits.open(elem+'_ref.fits')[2].data
        ref_full=fits.open(elem+'_ref_full.fits')[2].data
        full=fits.open(elem+'_full.fits')[2].data
        full0=fits.open(elem+'_0_full.fits')[2].data
        comp=fits.open(elem+'_comp_full.fits')[2].data
        ehdu.header['ELEM'] = elem
        ehdu.header['CRVAL1'] = aspcap.logw0
        ehdu.header['CDELT1'] = aspcap.dlogw
        ehdu.header['CTYPE1'] = 'LOG(WAVELENGTH)'
        if plot :
            if i==-1 :
                ref=ehdu.data
            else :
                if filtdir is not None: 
                    filt=aspcap.aspcap2apStar(ascii.read(filtdir+ehdu.header['ELEM']+'.filt',format='fixed_width_no_header')['col1'])
                    mask=aspcap.aspcap2apStar(ascii.read(filtdir+ehdu.header['ELEM']+'.mask',format='fixed_width_no_header')['col1'])
                    if elem == 'Ti' :
                        filt2=aspcap.aspcap2apStar(ascii.read(filtdir+ehdu.header['ELEM']+'II.filt',format='fixed_width_no_header')['col1'])
                        mask2=aspcap.aspcap2apStar(ascii.read(filtdir+ehdu.header['ELEM']+'II.mask',format='fixed_width_no_header')['col1'])
                    wind=ascii.read(filtdir+ehdu.header['ELEM']+'.wind')
                    j=np.where(wind['col3'] > -1)[0]
                    nwind=len(j)
                    wrange=[]
                    for jj in j : wrange.append((wind['col1'][jj],wind['col2'][jj]))
                fig,ax=plots.multi(3,3,hspace=0.001,wspace=0.001,xtickrot=60,figsize=(12,8))
                nspec = ehdu.data.shape[0]
                ispec=0
                gd=[]
                wgrid=[]
                xt=[]
                for ix in range(3) :
                    te=3500+ix*1000
                    xt.append('Teff {:6.0f}'.format(te))
                    row=[]
                    yt=[]
                    for iy in range(3) :
                        logg=1.+iy*2.
                        yt.append('log g: {:5.1f}'.format(logg))
                        x=10.**spectra.fits2vector(ehdu.header,1)
                        y=ehdu.data[ispec,:]/ref[ispec,:]
                        z=ref_full[ispec,:]
                        w=full[ispec,:]
                        plots.plotl(ax[iy,ix],x,y,yr=[0.9,1.2],xr=[15100,16950],color='b')
                        plots.plotl(ax[iy,ix],x,z,yr=[0.9,1.2],xr=[15100,16950],color='orange')
                        plots.plotl(ax[iy,ix],x,w,yr=[0.9,1.2],xr=[15100,16950],color='g')
                        plots.plotl(ax[iy,ix],x,full0[ispec,:],yr=[0.9,1.2],xr=[15100,16950],color='y')
                        plots.plotl(ax[iy,ix],x,comp[ispec,:]/ref_full[ispec,:],yr=[0.9,1.2],xr=[15100,16950],color='c')
                        ax[iy,ix].text(0.05,0.9,'Teff:{:6.0f} logg:{:6.1f}'.format(te,logg),transform=ax[iy,ix].transAxes)
                        if mask is not None: plots.plotl(ax[iy,ix],x,mask*0.2+1.005)
                        j=np.where(y < 0.99)[0]
                        gd.extend(j)
                        if nwind < 12 :
                            wfig,wax=plots.multi(nwind,1,brokenx=True,wspace=0.1,figsize=(8,4),xtickrot=45)
                            if nwind == 1 : wax=np.array([wax])
                            for iw in range(nwind) :
                                plots.plotl(wax[iw],x,y,xr=wrange[iw],yr=[0.5,1.2],color='b')
                                plots.plotl(wax[iw],x,z,color='orange')   # reference abun
                                plots.plotl(wax[iw],x,w,color='g')        # with 0.1 enhanced element
                                plots.plotl(wax[iw],x,full0[ispec,:],color='y') # without any element
                                plots.plotl(wax[iw],x,comp[ispec,:]/ref_full[ispec,:],color='c')  # with 0.1 ennhanced other elements
                                if mask is not None: plots.plotl(wax[iw],x,mask*0.2+1.005,color='r')
                                if filt is not None: plots.plotl(wax[iw],x,filt*0.2+1.005,ls=':',color='r')
                                if elem == 'Ti' and iw == nwind - 1 :
                                  plots.plotl(wax[iw],x,mask2*0.2+1.005,color='r')
                                  plots.plotl(wax[iw],x,filt2*0.2+1.005,ls=':',color='r')
                                wax[iw].ticklabel_format(style='plain',useOffset=False)
                            figname='{:s}_{:d}.png'.format(ehdu.header['ELEM'].strip(),ispec)
                            wfig.savefig(figname)
                            plt.close()
                            row.append(figname)
                        ispec+=1
                    wgrid.append(row)
                html.htmltab(wgrid,file=elem+'.html',xtitle=xt,ytitle=yt)
                fp.write('<h2>{:s}</h2>'.format(elem))
                if nwind < 12 :
                    fp.write(html.table(wgrid,xtitle=xt,ytitle=yt))
                else : fp.write('<br> Number of windows: {:d}'.format(nwind))

                mask=np.zeros(ehdu.data.shape[-1])
                mask[list(set(gd))] = 1.
                figname = ehdu.header['ELEM'].strip()+'.png'
                fig.savefig(figname)
                plt.close()
                grid.append([figname])
                ytit.append(ehdu.header['ELEM'])
                mhdu=fits.ImageHDU(mask)
                mhdu.header['ELEM'] = name.replace('.dat','')
                mhdu.header['CRVAL1'] = aspcap.logw0
                mhdu.header['CDELT1'] = aspcap.dlogw
                mhdu.header['CTYPE1'] = 'LOG(WAVELENGTH)'
                hdumask.append(mhdu)
        ehdu.data -= ref
        hdu.append(ehdu)
    # output single file
    if outfile is not None: 
        hdu.writeto(outfile+'.fits',overwrite=True)
        hdumask.writeto(outfile+'_mask.fits',overwrite=True)
    #html.htmltab(grid,ytitle=ytit,file=htmlfile)
    html.tail(fp)

def mkmask(file='elemsens')  :
    """ Routine to suggest possible masks
    """
    mask=fits.open(file+'.fits')
    els=[]
    for i in range(len(mask)) : els.append(mask[i].header['ELEM'])
    els = np.array(els)
    alphas=np.array(['O','Mg','Si','S','Ca','Ti'])
    metals=np.array(['Na','Al','P','K','V','Cr','Mn','Co','Fe','Ni','Cu','Ge','Rb','Ce','Nd'])
    fig,ax=plots.multi(1,2,hspace=0.001,sharex=True)
    x=10.**spectra.fits2vector(mask[0].header,1)
    for i,el in enumerate(els) :
        print(el)
        gd=[]
        ax[0].cla()
        for k in range(9) : 
            plots.plotl(ax[0],x,mask[i].data[k,:])
            gd.extend(np.where(mask[i].data[k,:] < -0.01)[0])
        if el in alphas :
            bd=[]
            for al in alphas :
                ax[1].cla()
                if el != al :
                    print(el,al)
                    j=np.where(els == al)[0][0]
                    for k in range(9) : 
                        plots.plotl(ax[1],x,mask[j].data[k,:])
                        #plots.plotl(ax[1],x,mask[j].data[k,:]/mask[i].data[k,:] )
                        bd.extend(np.where((mask[j].data[k,:]/mask[i].data[k,:] < 0.2) & (mask[j].data[k,:]<-0.01) )[0])
                    plt.show()
                    #gd=np.where(abs(mask[j].data) > 0.)[0]
                    #mask[i].data[gd] = -1.*mask[i].data[gd]
        elif el in metals :
            for al in metals :
                ax[1].cla()
                if el != al :
                    print(el,al)
                    j=np.where(els == al)[0][0]
                    for k in range(9) : 
                        plots.plotl(ax[1],x,mask[j].data[k,:])
                        #plots.plotl(ax[1],x,mask[j].data[k,:]/mask[i].data[k,:] )
                        bd.extend(np.where((mask[j].data[k,:]/mask[i].data[k,:] < 0.2) & (mask[j].data[k,:]<-0.01) )[0])
                    #gd=np.where(abs(mask[j].data) > 0.)[0]
                    #print(el,al,j,len(gd))
                    #mask[i].data[gd] = -1.*abs(mask[i].data[gd])
                    plt.show()

        new=np.zeros(mask[i].data.shape[-1])
        new[gd] = 1.
        new[bd] = -1*new[bd]
        ax[1].cla()
        plots.plotl(ax[1],x,new)
          
        plt.draw()
        plt.show()
        pdb.set_trace()

def filter_lines(infile,outfile,wind,nskip=0) :
    """ Read from input linelist file, output comments and lines falling in windows of [w1,w2] to outfile

    Args :
        infile (str) : name of input file
        outfile (str) : name of output file
        wind (list of [w1,w2] pairs) : list of window ranges
        nskip (int) : number of header lines to skip (default=0)

    Returns :
        nout (int) : number of windows
    """
    fin=open(infile,'r')
    fout=open(outfile,'w')
    nout = 0
    nwind=wind.shape[0]
    for iline,line in enumerate(fin) :
        if iline >= nskip : 
            w=line.split()[0]
            if w[0] == '#' :
                    fout.write(line)
            else :
              for i in range(nwind) :
                if (float(w) >= wind[i,0]) and (float(w) <=wind[i,1]) :
                    fout.write(line)
                    nout+=1
    fout.close()
    return nout
 
def mini_linelist(elem,linelist,maskdir=None,only=False,clobber=False,outdir=None,delta_loggf=0., delta_vdw=0.) :
    """ Produce abbreviated Turbospec linelists, e.g. for minigrid construction, given mask file and linelist file IN AIR
        With only, produce linelist with only lines from input element 
        Return arrays of wavelength ranges wind,wair
    """

    # get window ranges in vacuum and convert to air
    if maskdir is not None :
        wind=np.loadtxt(os.environ['APOGEE_DIR']+'/data/windows/'+maskdir+'/'+elem+'.wave')
        nwind=wind.shape[0]
        wair=spectra.vactoair(wind)
    else : 
        nwind = 1
        wair=np.zeros([2,2])
        wair[0,0] = -1.
        wair[0,1] = 1.e10
        wind = wair

    # setup output directory
    if outdir == None :
        if only : outdir = os.environ['APOGEE_SPECLIB']+'/linelists/'+elem+'_only/'
        else : outdir = os.environ['APOGEE_SPECLIB']+'/linelists/'+elem+'/'
    try: os.mkdir(outdir)
    except: pass

    # if files are being create by other process, wait until done
    while os.path.isfile(outdir+elem+'.lock') : 
        print('waiting for lock: ',outdir+elem+'.lock')
        time.sleep(10)

    # if files are already created, return, otherwise open .lock file and create
    if not clobber and os.path.isfile(outdir+elem+'.done') : return wind,wair
    fp = open(outdir+elem+'.lock','w')
    fp.close()

    # convert Turbospectrum files to filtered Turbospectrum files
    # Turbospectrum files are in air wavelengths
    lists=['turbospec.'+linelist+'.atoms','turbospec.'+linelist+'.molec',
           'turbospec.'+linelist+'.Hlinedata','turbospec.h2o-BC8.5V.molec','turbospec.h2o-BC9.5V.molec']
    for i,linelist in enumerate(lists) :
        filepath=os.environ['APOGEE_SPECLIB']+'/linelists/'+linelist
        fout=open(outdir+linelist,'w')
        with open(filepath) as fp:  
            out = ''
            nelem = 0
            n = 0
            line = fp.readline()
            while line :
                if line[0] == "'" :
                    # we have a new element
                    if nelem > 0 :
                        # if it's not the first element, write out the previous one!
                        if n > 0 :
                            if 'molec' in linelist :
                                tmp = int(float(head.split("'")[1]))
                                elemcode = [tmp//100,tmp%100]
                            else :
                                elemcode = [int(float(head.split("'")[1]))]
                            if not only or atomic.periodic(elem) in elemcode :
                                j=head.split("'")[2].split()[0]
                                # for the header line, include the new number of lines
                                fout.write("'"+head.split("'")[1]+"'   "+j+'{:10d}\n'.format(n))
                                # write the accumlated data output
                                fout.write(out)
                            n=0
                    head = line
                    if 'molec' in linelist :
                        tmp = int(float(head.split("'")[1]))
                        elemcode = [tmp//100,tmp%100]
                    else :
                        elemcode = [int(float(head.split("'")[1]))]
                    # start the line data output with the comment line
                    out = fp.readline()
                    nelem += 1
                else :
                    # accumulate the linelist for this element if it's within the desired range
                    if 'atoms' in linelist :
                        pars=parse_turbo(line)
                        w = pars[0]
                    else :
                        w = line.split()[0]
                    for i in range(nwind) :
                      if (float(w) >= wair[i,0]) and (float(w) <=wair[i,1]) : 
                          if 'atoms' in linelist :
                              if atomic.periodic(elem) in elemcode :
                                  pars[2] += delta_loggf
                                  pars[3] += delta_vdw
                              out=out+write_turbo(pars)
                          else :
                              out=out+line  
                          n+=1
                line = fp.readline()          
            if n > 0 :
                # last element
                if 'molec' in linelist :
                    tmp = int(float(head.split("'")[1]))
                    elemcode = [tmp//100,tmp%100]
                else :
                    elemcode = [int(float(head.split("'")[1]))]
                if not only or atomic.periodic(elem) in elemcode :
                    j=head.split("'")[2].split()[0]
                    fout.write("'"+head.split("'")[1]+"'   "+j+'{:10d}\n'.format(n))
                    fout.write(out)
        fout.close()


    # write .done file and remove .lock
    fp = open(outdir+elem+'.done','w')
    fp.close()
    os.remove(outdir+elem+'.lock') 
    return wind,wair

def plotcross(a,val=[0,0,0],hard=None,sum=True) :
    """ plot cross sections of input 3D grid of synthetic spectra

        Args:
            a : input HDU with synthetic spectral grid
           val : 3 ([M/H], logg, Teff) grid indices to use when varying other dimensions (default=[0,0,0])
    """
    dim=a.data.shape
    # if we have a rotation dimension, take the lowest rotation
    if len(dim) == 5 : 
        data=np.squeeze(a.data[0,:,:,:,:])
        dim=data.shape
    else : data = a.data
    n=dim[3]
    x=10.**spectra.fits2vector(a.header,1)
    mh=spectra.fits2vector(a.header,4)
    logg=spectra.fits2vector(a.header,3)
    teff=spectra.fits2vector(a.header,2)
    colors=['r','g','b','c','m','y','black']*3
    for i in range(dim[0]) :
        if i == 0 : fig,ax=aspcap.plot(x,data[i,val[1],val[2],:],color=colors[i],sum=sum)
        else : aspcap.plot(x,data[i,val[1],val[2],:],ax=ax,color=colors[i],sum=sum)
    fig.suptitle('[M/H] varied from {:6.2} to {:6.2f} at logg {:6.1f}, Teff {:6.0f}'.format(mh[0],mh[-1],logg[val[1]],teff[val[2]]))
    if hard is not None : 
        fig.savefig(hard+'_mh.pdf')
        plt.close()
    for i in range(dim[1]) :
        if i == 0 : fig,ax=aspcap.plot(x,data[val[0],i,val[2],:],color=colors[i],sum=sum)
        else : aspcap.plot(x,data[val[0],i,val[2],:],ax=ax,color=colors[i],sum=sum)
    fig.suptitle('log g varied from {:6.1} to {:6.2f} at [M/H] {:6.2f}, Teff {:6.0f}'.format(logg[0],logg[-1],mh[val[0]],teff[val[2]]))
    if hard is not None : 
        fig.savefig(hard+'_logg.pdf')
        plt.close()
    for i in range(dim[2]) :
        if i == 0 : fig,ax=aspcap.plot(x,data[val[0],val[1],i,:],color=colors[i],sum=sum)
        else : aspcap.plot(x,data[val[0],val[1],i,:],ax=ax,color=colors[i],sum=sum)
    fig.suptitle('Teff varied from {:6.0f} to {:6.0f} at logg {:6.1f}, [M/H] {:6.2f}'.format(teff[0],teff[-1],logg[val[1]],mh[val[0]]))
    if hard is not None : 
        fig.savefig(hard+'_teff.pdf')
        plt.close()

    dim=a.data.shape
    if len(dim) == 5 :
        for i in range(dim[0]) :
            if i == 0 : fig,ax=aspcap.plot(x,data[i,val[0],val[1],val[2],:],color=colors[i],sum=sum)
            else : aspcap.plot(x,data[i,val[0],val[1],val[2],:],ax=ax,color=colors[i],sum=sum)
        #fig.suptitle('vsini varied from {:6.0f} to {:6.0f} at logg {:6.1f}, [M/H] {:6.2f}'.format(teff[0],teff[-1],logg[val[1]],mh[val[0]]))
        if hard is not None : 
            fig.savefig(hard+'_vsini.pdf')
            plt.close()

def co( teff=3400,logg=1,mh=0,am=0,cm=0,nm=0,vmicro=1.2,nearest=False) :
    """ Create set of syntheses with different C/O ratios
        using 1) synthesis with grid atmospheres, and 2) synthesis
        with custom self-consistent atmospheres
        Also create FERRE ipf file to use to make FERRE interpolated set
    """
    fp=open('synth.ipf','w')
    syn=[]
    syn2=[]
    #for carbon in np.arange(0,0.51,0.025) :
    for co in [0.5,0.537,0.75,0.899,0.925,0.951,0.971,0.991] :
        carbon = np.log10(co) + 0.27
        if nearest: cm = round(carbon/0.25)*0.25
        if np.isclose(cm,0.) : cm= 0.
        # self-consistent synthesis
        out = mk_synthesis('turbospec',teff,logg,mh,am,cm,nm,vmicro=vmicro,
                 els=[['C',carbon]],atmod=os.environ['APOGEE_SPECLIB']+'/atmos/marcs/CO/3400g1.0m1.0z0.00a0.00c{:5.3f}s0.00.mod'.format(co))
        syn2.append(out[0])
        # synthesis with grid atmosphere
        out = mk_synthesis('turbospec',teff,logg,mh,am,cm,nm,vmicro=vmicro,
                 els=[['C',carbon]] )
        syn.append(out[0])
        fp.write('test_{:.3f}{:9.3f}{:9.3f}{:9.3f}{:9.3f}{:9.3f}{:9.3f}{:9.3f}\n'.format(carbon,np.log10(vmicro),carbon,nm,am,mh,logg,teff))
    fp.close()

    lsfid=14600018
    waveid=13140000
    highres=9
    apred='dr17'
    prefix='lsf_'
    telescope='apo25m'
    fiber=150
    x, ls = getlsf(lsfid,waveid,prefix=prefix,apred=apred,telescope=telescope,fiber=fiber,highres=highres)
    conv=[]
    conv2=[]
    ws=np.linspace(15100.,17000., 38001)
    # synthesis is in air, we want vacuum
    ws=spectra.airtovac(ws)
    for spec,spec2 in zip(syn,syn2) :
        vmacro = 10.**(0.470794-0.254*mh)
        vmacro = vmacro if vmacro<15 else 15.
        # convolve one at a time because we have different vrot for each
        z,waveout=lsf.convolve(ws,spec,lsf=ls,xlsf=x,vmacro=vmacro,vrot=None)
        z=np.squeeze(z)
        conv.append(z/np.nanmedian(z))
        z,waveout=lsf.convolve(ws,spec2,lsf=ls,xlsf=x,vmacro=vmacro,vrot=None)
        z=np.squeeze(z)
        conv2.append(z/np.nanmedian(z))

    asp=[]
    for spec in conv2 : asp.append(aspcap.apStar2aspcap(spec)/np.nanmedian(aspcap.apStar2aspcap(spec)))
    asp=np.array(asp)
    ferre.writespec('synth.obs',asp)
    ferre.writespec('synth.err',asp*0+0.001)

    return syn,syn2,conv,conv2

def plotco(syn,syn2,ferre) :
    """ Plot the C/O results
    """
    co = [0.5,0.537,0.75,0.899,0.925,0.951,0.971,0.991] 
    fig,ax=plots.multi(2,3,hspace=0.001,sharex=True,figsize=(16,8),wspace=0.5)
    for i in range(8) :
        tit = '{:8.3f} [{:6.2f}]'.format(co[i],np.log10(co[i])+0.27)
        ax[0,0].plot(aspcap.apStarWave(),syn[i])
        ax[1,0].plot(aspcap.apStarWave(),syn2[i])
        ax[2,0].plot(aspcap.apStarWave(),aspcap.aspcap2apStar(ferre[i]))
        ax[0,1].plot(aspcap.apStarWave(),syn[i]-syn2[i],label=tit)
        ax[1,1].plot(aspcap.apStarWave(),syn2[i]-syn2[i],label=tit)
        ax[2,1].plot(aspcap.apStarWave(),aspcap.aspcap2apStar(ferre[i])-syn2[i])
    ax[1,1].legend(fontsize='x-small')
    for i in range(3) : ax[i,0].set_ylim(0.4,1.25)
    for i in range(3) : ax[i,1].set_ylim(-0.2,0.2)

def interpol_marcs(teff,logg,mh,am,cm,nm,workdir='./') :
    """ Interpolate in atmosphere grid for desired atmosphere, using Thomas Masseron's interpol_modeles
    """

    # find surrounding grid points in APOGEE grid
    if teff < 4000. :
       tlo = int(teff/100.)*100.
       thi = tlo + 100.
    else :
       tlo = int(teff/250.)*250.
       thi = tlo + 250.
    # add 10 to make sure we are positive for int
    glo = int((logg+10)/0.5)*0.5 - 10
    ghi = glo + 0.5
    zlo = int((mh+10)/0.25)*0.25 - 10
    zhi = zlo + 0.25

    cmtmp,amtmp = get_cmnear(cm, am)

    cwd = os.getcwd()
    os.chdir(workdir)
    if np.isclose(teff,tlo) : teff=tlo+0.01
    if np.isclose(teff,thi) : teff=thi-0.01
    if np.isclose(logg,glo) : logg=glo+0.001
    if np.isclose(logg,ghi) : logg=ghi-0.001
    if np.isclose(mh,zlo) : mh=zlo+0.001
    if np.isclose(mh,zhi) : mh=zhi=-0.001

    # set up input file
    test = '.false.'
    cmd=['interpol_modeles']
    models=[]
    for t in [tlo,thi] :
        for g in [glo,ghi] :
            for z in [zlo,zhi] :
                models.append(get_atmod_file(t,g,z,amtmp,cmtmp,nm,fill=True))
    model_out=os.path.basename(atmos.filename(teff,logg,mh,cmtmp,amtmp,model='marcs'))
    model_alt='interp.alt'

    fp=open(model_out.replace('.mod','.inp'),'w')
    for model in models :
        fp.write(os.path.basename(model)+'\n')
    fp.write(os.path.basename(model_out)+'\n')
    fp.write(os.path.basename(model_alt)+'\n')
    fp.write('{:20.12f}\n'.format(teff))
    fp.write('{:18.12f}\n'.format(logg))
    fp.write('{:18.12f}\n'.format(mh))
    fp.write('.false.\n')
    fp.write('.false. {:12.6f} {:12.6f} {:12.6f}\n'.format((teff-tlo)/(thi-tlo),(logg-glo)/(ghi-glo),(mh-zlo)/(zhi-zlo)))   #fracinput
    fp.write('.true.\n')    #optimize
    fp.write('.false.\n')   #binary
    fp.write('dummy\n')
    fp.close()

    fp=open(model_out.replace('.mod','.inp'))
    p=subprocess.call(cmd,stdin=fp,stdout=subprocess.DEVNULL,stderr=subprocess.DEVNULL)
    fp.close()

    os.chdir(cwd)

    return model_out

def parse_turbo(line) :

    out=[]
    out.append(float(line[0:10]))
    out.append(float(line[10:17]))
    out.append(float(line[17:25]))
    out.append(float(line[25:35]))
    out.append(float(line[35:42]))
    out.append(float(line[42:52]))
    out.append(float(line[52:59]))
    out.append(line[62])
    out.append(line[66])
    out.append(float(line[68:72]))
    out.append(float(line[72:76]))
    out.append(line[78:-2])

    return out

def write_turbo(pars) :

    return "{:10.3f}{:7.3f}{:8.3f}{:10.2f}{:7.1f}{:10.2e}{:7.2f}  '{:s}' '{:s}'{:4.1f}{:4.1f} '{:s}'\n".format(*pars)
           #pars[0],pars[1],pars[2],pars[3],pars[4],pars[5],pars[6],pars[7],pars[8],pars[9
#         1         2         3         4         5         6         7         8         9        10
#12345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
# 16549.594  4.521  -3.312     -6.21    2.0  5.01e+06  -2.39  'x' 'x' 0.0 1.0 'LI  I   1s2.4p 2P   1s2.11s 2S'
# 16814.155  4.541  -2.790     -5.54    2.0  3.02e+07   0.00  'x' 'x' 0.0 1.0 'LI  I   4d  2D      11p  2P   '
