#!/usr/bin/env python

from __future__ import print_function

import argparse
import math
import os
import re
import sys

import numpy

import lfa

# Stop pyplot trying to use X.
import matplotlib
matplotlib.use('Agg')

import fftrv

from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt

from medsig import *

from eqwidth import *
from lsd import *
from makesky import *
from prepord import *
from read_spec import *
from vsini import *

# Constants used by the script go here.

# Figure size based on pgplot.
figsize = (10.5, 7.8)  # inches

# Velocity range to plot
velrange = 200  # km/s, plots +/- this much

# Regularization constant for LSDs.  Depends on template and possibly
# also on SNR, so may need changing.
kreg = 100

# Limb darkening coefficients for rotational broadening kernel.
# These are from Phoenix mod. atm. for Teff=3224K, logg=5, Rc band,
# least-squares method, Claret et al. (2012) 2012A&A...546A..14C
ald = 0.4629
bld = 0.3212

def do_vrad(pdf, tmplname, filename,
            tmpl_mbjd, tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk, tmpl_vbcv, tmpl_vrad,
            mbjd, wave, flux, e_flux, msk, vbcv,
            order, dovsini=0, vsini=0, savefile=False, emchop=True):
  # Extract order and clean.
  thistmpl_wave, thistmpl_flux, thistmpl_e_flux = prepord(order, tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk)
  thiswave, thisflux, thise_flux = prepord(order, wave, flux, e_flux, msk)

  # Take off sky.
  ss = makesky(thistmpl_wave, thistmpl_flux, 4)

  thistmpl_flux -= ss

  thisrawflux = numpy.copy(thisflux)

  ss = makesky(thiswave, thisflux, 4)

  thisflux -= ss

  # Mask NaNs.
  emmask = numpy.isfinite(thisflux)

  if emchop:
    medflux, sigflux = medsig(thisflux[emmask])
    emmask = numpy.logical_and(emmask,
                               thisflux < medflux + 5.0*sigflux)

  # Number of measurements (pixels).
  npix = len(thisflux)

  if dovsini:
    # Figure out vsini first.
    frv = fftrv.fftrv(nbin=32*npix, pkfit=0,
                      t_emchop=emchop, s_emchop=emchop)

    initvsini, l_vsini, l_p = do_vsini_grid(frv,
                                            thistmpl_wave, thistmpl_flux,
                                            thiswave, thisflux,
                                            ald, bld,
                                            0, 100, 101)

    avsini = initvsini - 1.0
    if avsini < 0.0:
      avsini = 0.0

    bvsini = avsini + 2.0

    newvsini, r_vsini, r_p = do_vsini_grid(frv,
                                           thistmpl_wave, thistmpl_flux,
                                           thiswave, thisflux,
                                           ald, bld,
                                           avsini, bvsini, 21)

    vsini = newvsini
    
  usevsini = vsini

  # Correlate.
  frv = fftrv.fftrv(nbin=32*npix, t_emchop=emchop, s_emchop=emchop)

  z, corr, zbest, hbest, sigt = frv.correlate(thistmpl_wave, thistmpl_flux,
                                              thiswave, thisflux,
                                              usevsini*1000 / lfa.LIGHT,
                                              ald, bld)

  vels = z * lfa.LIGHT / 1000
  vbest = zbest * lfa.LIGHT / 1000

  ww = numpy.abs(vels-vbest) < velrange

  # Optionally save to file.
  if savefile:
    basefile = stripname(filename)
  
    if dovsini:
      rotfile = basefile + "_rot_det.txt"
  
      with open(rotfile, "w") as rfp:
        for imeas, thisvsini in enumerate(l_vsini):
          rfp.write(str(thisvsini) + " " + str(l_p[imeas]) + "\n")
  
    ccffile = basefile
    if dovsini:
      ccffile += "_rot_ccf.txt"
    else:
      ccffile += "_nonrot_ccf.txt"

    with open(ccffile, "w") as cfp:
      for imeas, thisvel in enumerate(vels):
        cfp.write(str(thisvel) + " " + str(corr[imeas]) + "\n")

  # Plot.
  fig = plt.figure(figsize=figsize)

  plt.subplot(2, 1, 1)

  if dovsini:
    plt.plot(l_vsini, l_p)
    plt.axvline(vsini)
    plt.xlim(l_vsini[0], l_vsini[-1])
    plt.xlabel("vsini (km/s)")
    plt.ylabel("Correlation^2")
    plt.title("{0:s}".format(filename))
  else:
    plt.plot(thiswave, thisrawflux)
    plt.xlim(thiswave[0], thiswave[-1])
    yl = numpy.min(thisrawflux[emmask])
    yh = numpy.max(thisrawflux[emmask])
    plt.ylim(yl-0.05*(yh-yl), yh+0.05*(yh-yl))
    plt.xlabel("Wavelength (A)")
    plt.ylabel("Counts")
    plt.title("{0:s} aperture {1:d}".format(filename, order))
  
  plt.subplot(2, 1, 2)
  plt.plot(vels[ww], corr[ww])
  plt.xlim(vels[ww][0], vels[ww][-1])
  plt.xlabel("Barycentric velocity (km/s)")

  if dovsini:
    plt.ylabel("Rot. correlation")
  else:
    plt.ylabel("Nonrot. correlation")
  
  xl, xh = plt.xlim()
  yl, yh = plt.ylim()

  plt.text(xh-0.02*(xh-xl), yh-0.1*(yh-yl),
           "tmpl = {0:s}".format(tmplname),
           horizontalalignment='right')
  plt.text(xh-0.02*(xh-xl), yh-0.2*(yh-yl),
           "Bary. RV = {0:.3f} km/s".format(vbest),
           horizontalalignment='right')
  plt.text(xh-0.02*(xh-xl), yh-0.3*(yh-yl),
           "h = {0:.3f}".format(hbest),
           horizontalalignment='right')

  if usevsini > 0:
    plt.text(xh-0.02*(xh-xl), yh-0.4*(yh-yl),
             "vsini = {0:.2f} km/s".format(usevsini),
             horizontalalignment='right')

  # Indicate where junk appears.
  tellvel = vbcv + tmpl_vrad - tmpl_vbcv

  # Tellurics in target spectrum x Tellurics in template.
  plt.axvline(tellvel, color="black", linestyle="--")
  plt.text(tellvel, yh+0.01*(yh-yl), "t",
           fontsize=10, horizontalalignment="center")

  # Sky in target x template spectrum.
  plt.axvline(vbcv, color="black", linestyle="--")
  plt.text(vbcv, yh+0.01*(yh-yl), "s",
           fontsize=10, horizontalalignment="center")

  plt.axvline(vbest)
  
  pdf.savefig(fig)
  plt.close()

  return vbest, hbest, usevsini

def do_lsd(pdf, filename,
           tmpl_mbjd, tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk, tmpl_vbcv, tmpl_vrad,
           mbjd, wave, flux, e_flux, msk, vbcv,
           vrad, orders, savefile=False, emchop=True):
  # Compute at 0.5 km/s intervals.
  vl = vrad-velrange
  vh = vrad+velrange
  nv =  801

  vels, prof = lsd_multiorder(tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk,
                              wave, flux, e_flux, msk,
                              orders,
                              vl, vh, nv,
                              kreg, emchop=emchop)

  # Optionally save LSD to file.
  if savefile:
    basefile = stripname(filename)
    lsdfile = basefile + "_lsd.txt"

    with open(lsdfile, "w") as lfp:
      for imeas, thisvel in enumerate(vels):
        lfp.write(str(thisvel) + " " + str(prof[imeas]) + "\n")

  # Plot.
  fig = plt.figure(figsize=figsize)

  plt.plot(vels, prof)

  plt.xlim(vels[0], vels[-1])

  plt.xlabel("Barycentric velocity (km/s)".format(vrad))
  plt.ylabel("Normalized flux")
  plt.title("{0:s} LSD".format(filename))

  # Indicate where junk appears.
  yl, yh = plt.ylim()

  tellvel = vbcv + tmpl_vrad - tmpl_vbcv

  # Tellurics in target spectrum x Tellurics in template.
  plt.axvline(tellvel, color="black", linestyle="--")

  # Sky in target x template spectrum.
  plt.axvline(vbcv, color="black", linestyle="--")

  pdf.savefig(fig)
  plt.close()

def plotreg(speclist, tresorder, xmark, xwidth, name, legend, a, b, c):
  xcen = numpy.average(xmark)

  wl = xcen*(1.0-xwidth)
  wh = xcen*(1.0+xwidth)

  ll = [None] * len(speclist)

  have = False
  ymin = None
  totyr = 0

  for ispec, spec in enumerate(speclist):
     mbjd, wave, flux, e_flux, msk = spec

     # Where are the pixels within the desired wavelength range?
     ww = numpy.logical_and(wave > wl, wave < wh)
     ww = numpy.logical_and(ww, numpy.isfinite(flux))

     # Figure out which order to plot.  Take the one with the most
     # median counts if more than one contains the desired range.
     nord, npix = wave.shape

     best_iord = None
     best_dl = None
     best_med = None

     for iord in range(nord):
       thisww = ww[iord,:]

       if numpy.any(thisww):
         # Wavelength range in order.
         lmin = numpy.min(wave[iord,thisww])
         lmax = numpy.max(wave[iord,thisww])

         # Do we have the feature?
         if xcen >= lmin and xcen <= lmax:
           dl = lmax-lmin
           thismed = numpy.median(flux[iord,thisww])
           
           # Largest amount of window covered, then best s/n.
           if best_med is None or dl > best_dl or thismed > best_med:
             best_iord = iord
             best_dl = dl
             best_med = thismed

     if best_iord is None:  # no pixels in desired range
       continue

     thiswave = wave[best_iord,:]
     thisflux = flux[best_iord,:]
     thisww = ww[best_iord,:]

     med = numpy.median(thisflux[thisww])

     if med == 0:  # can't normalize
       continue

     if not have:
       have = True
       plt.subplot(a, b, c)  # init plot

     normflux = thisflux[thisww] / med

     # Compute statistics to make mask of emission features in whole order.
     ss = makesky(thiswave, thisflux, 4)
     tt = thisflux - ss

     emmask = numpy.isfinite(thisflux)

     medtt, sigtt = medsig(tt[emmask])
     emmask = numpy.logical_and(emmask,
                                tt < medtt + 5.0*sigtt)

     # Make sure the present feature is excluded in case
     # we're dealing with an emission feature.
     featmask = numpy.logical_and(thiswave >= xcen-1,
                                  thiswave <= xcen+1)

     emmask[featmask] = 1

     # Extract the relevant part.
     thisemmask = emmask[thisww]

     # Compute plot range.
     thisymin = numpy.min(normflux[thisemmask])
     thisymax = numpy.max(normflux[thisemmask])

     if ispec == 0:
       ymin = thisymin

     thisyr = thisymax-thisymin

     plt.plot(thiswave[thisww], (normflux-thisymin) + ymin + totyr)

     totyr += thisyr

  if have:
    plt.xlim(wl, wh)
    plt.ylim(ymin-0.05*totyr, ymin+1.2*totyr)
    
    plt.locator_params(nbins=5)
    plt.tick_params(axis="y", which="both", labelleft=False)
    
    ax = plt.gca()
    ax.ticklabel_format(useOffset=False)
    
    if hasattr(xmark, "__iter__"):
      for x in xmark:
        plt.axvline(x, 0.9, 0.95, color="black")
    else:
      plt.axvline(xmark, 0.9, 0.95, color="black")    

    if legend is not None:
      plt.text(wh-0.05*(wh-wl),
               ymin-0.05*totyr+0.95*1.25*totyr,
               legend,
               horizontalalignment="right",
               verticalalignment="top",
               fontsize=10)

      # BODGE
      plt.axvspan(xcen-3, xcen+3, alpha=0.05, color="black")

    #  plt.xlabel("Wavelength (A)")
    #  plt.ylabel("Counts")
    plt.title(name)

  return have

def do_plotfeatures(pdf, filename,
                    mbjd, wave, flux, e_flux, msk,
                    vrad, ewha, e_ewha, ewcai, e_ewcai):
  npx = 4
  npy = 2
  width = 0.001
  
  fig = plt.figure(figsize=figsize)

  restwave = wave / (1.0 + vrad * 1000 / lfa.LIGHT)

  speclist = [(mbjd, restwave, flux, e_flux, msk)]

  featlist = [(2, 3933.66, width, "Ca II K", None),
              (3, 3968.47, width, "Ca II H", None),
              (37, 6562.808, width, "Halpha", "{0:.1f}$\\pm${1:.1f}A".format(ewha, e_ewha) if ewha is not None else None),
              (37, 6572.795, width, "Ca I", "{0:.1f}$\\pm${1:.1f}A".format(ewcai, e_ewcai) if ewcai is not None else None),
              (38, 6707.79, width, "Li I", None),
              (47, [8183.25, 8194.81] , 2*width, "Na I", None),
              (49, 8498.02, width, "Ca II", None),
              (49, 8542.09, width, "Ca II", None),
              (49, 8662.14, width, "Ca II", None)]

  panel = 1
  for feat in featlist:
    tresord, cenwave, xwidth, name, legend = feat

    h = plotreg(speclist, tresord, cenwave, xwidth, name, legend, npy, npx, panel)
    if h:
      panel += 1
      if panel > npx*npy:
        break

  fig.text(0.5, 0.05, "Wavelength (A)", ha="center", va="center")
  fig.text(0.05, 0.5, "Normalized counts + constant", ha="center", va="center", rotation="vertical")
  fig.text(0.5, 0.95, "{0:s}".format(filename), ha="center", va="center")
  
  if panel > 1:
    plt.tight_layout(rect=(0.05, 0.05, 0.95, 0.95))

  pdf.savefig(fig)
  plt.close()

def do_ew(tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk, tmpl_blaze,
          restwave, flux, e_flux, msk, blaze,
          feat, width, mskwidth,
          cont1_awave, cont1_bwave,
          cont2_awave, cont2_bwave):

  feat_awave = feat - 0.5*width
  feat_bwave = feat + 0.5*width

  featcorr1_awave = feat_awave

  if mskwidth is not None:
    featcorr1_bwave = feat - 0.5*mskwidth
    featcorr2_awave = feat + 0.5*mskwidth

  featcorr2_bwave = feat_bwave

  wl = min(feat_awave, cont1_awave, cont2_awave)
  wh = max(feat_bwave, cont1_bwave, cont2_bwave)

  # Where are the pixels within the desired wavelength range?
  ww = numpy.logical_and(restwave > wl, restwave < wh)

  # Figure out which order to use.  Take the one with the most
  # median counts if more than one contains the desired range.
  nord, npix = restwave.shape

  best_iord = None
  best_med = None
  
  for iord in range(nord):
    thisww = ww[iord,:]
    
    if numpy.any(thisww):
      # Wavelength range in order.
      lmin = numpy.min(restwave[iord,:])
      lmax = numpy.max(restwave[iord,:])
      
      # Do we have the whole feature?
      if wl >= lmin and wh <= lmax:
        thismed = numpy.median(flux[iord,thisww])
           
        # Largest amount of window covered, then best s/n.
        if best_med is None or thismed > best_med:
          best_iord = iord
          best_med = thismed

  if best_iord is not None:
    thistmpl_wave = tmpl_wave[best_iord,:]
    thistmpl_rawflux = tmpl_flux[best_iord,:]
    thistmpl_rawe_flux = tmpl_e_flux[best_iord,:]
    thistmpl_blaze = tmpl_blaze[best_iord,:]

    thiswave = restwave[best_iord,:]
    thisrawflux = flux[best_iord,:]
    thisrawe_flux = e_flux[best_iord,:]
    thisblaze = blaze[best_iord,:]

    # Divide blaze, carefully.
    tmpl_divb = numpy.where(thistmpl_blaze > 0, 1.0/thistmpl_blaze, 0.0)

    thistmpl_flux = thistmpl_rawflux * tmpl_divb
    thistmpl_e_flux = thistmpl_rawe_flux * tmpl_divb

    divb = numpy.where(thisblaze > 0, 1.0/thisblaze, 0.0)

    thisflux = thisrawflux * divb
    thise_flux = thisrawe_flux * divb

    # Equivalent width correction for continuum.
    tmpl_awave, tmpl_bwave, tmpl_width = subpixel(thistmpl_wave)

    if mskwidth is not None:
      # Continuum calculation for template (same as target).
      tmpl_cont1, tmpl_e_cont1 = average(tmpl_awave, tmpl_bwave, tmpl_width,
                                         thistmpl_flux, thistmpl_e_flux,
                                         cont1_awave, cont1_bwave)
      
      tmpl_cont2, tmpl_e_cont2 = average(tmpl_awave, tmpl_bwave, tmpl_width,
                                         thistmpl_flux, thistmpl_e_flux,
                                         cont2_awave, cont2_bwave)
      
      tmpl_cont = 0.5*(tmpl_cont1+tmpl_cont2)
      
      # Sum "false EW" contribution at sides of feature window.
      ss1, sv1, sl1 = eqwidth_sum(tmpl_awave, tmpl_bwave, tmpl_width,
                                  thistmpl_flux, thistmpl_e_flux,
                                  0.0, 0.0,
                                  featcorr1_awave, featcorr1_bwave)
      
      ss2, sv2, sl2 = eqwidth_sum(tmpl_awave, tmpl_bwave, tmpl_width,
                                  thistmpl_flux, thistmpl_e_flux,
                                  0.0, 0.0,
                                  featcorr2_awave, featcorr2_bwave)
      
      tmpl_false = (ss1+ss2) / (sl1+sl2)
      
      # Adjustment to continuum to null "false EW".
      contcorr = tmpl_false / tmpl_cont
    else:
      contcorr = 1.0

    awave, bwave, width = subpixel(thiswave)

    # Continuum calculation for target with correction applied.
    cont1, e_cont1 = average(awave, bwave, width, thisflux, thise_flux,
                             cont1_awave, cont1_bwave)
    
    cont2, e_cont2 = average(awave, bwave, width, thisflux, thise_flux,
                             cont2_awave, cont2_bwave)
    
    cont = 0.5*(cont1+cont2) * contcorr
    e_cont = 0.5*math.hypot(e_cont1, e_cont2) * contcorr

    if cont != 0:
      ewha, e_ewha = eqwidth(awave, bwave, width, thisflux, thise_flux,
                             cont, e_cont,
                             feat_awave, feat_bwave)

      # Report as absorption EW for historical reasons.
      return -ewha, e_ewha
    else:
      return None, None

  return None, None

def do_ews(tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk, tmpl_blaze,
           wave, flux, e_flux, msk, blaze,
           vrad):

  restwave = wave / (1.0 + vrad * 1000 / lfa.LIGHT)

  if tmpl_blaze is None:
    tmpl_blaze = numpy.ones_like(flux)

  if blaze is None:
    blaze = numpy.ones_like(flux)

  haline = 6562.808
  hawid = 6.0

  ewha, e_ewha = do_ew(tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk, tmpl_blaze,
                       restwave, flux, e_flux, msk, blaze,
                       haline, hawid, None,
                       haline-1.5*hawid, haline-0.5*hawid,
                       haline+0.5*hawid, haline+1.5*hawid)

  caline = 6572.795
  cawid = 6.0

  ewcai, e_ewcai = do_ew(tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk, tmpl_blaze,
                         restwave, flux, e_flux, msk, blaze,
                         caline, cawid, None,
                         caline-1.5*cawid, caline-0.5*cawid,
                         caline+0.5*cawid, caline+1.5*cawid)

  # Attempts to find a way to measure Li seem to be epic fail, so far.
#  ewli, e_ewli = do_ew(tmpl_wave, tmpl_flux, tmpl_e_flux, tmpl_msk, tmpl_blaze,
#                       restwave, flux, e_flux, msk, blaze,
#                       6707.79, 8.0, 4.0,
#                       6692.00, 6700.00,
#                       6692.00, 6700.00)

# 1.0614
# 1.0708

  return ewha, e_ewha, ewcai, e_ewcai

def get_range(s):
  ll = s.split(":")
  if len(ll) == 2:
    return map(float, ll)
  else:
    raise argparse.ArgumentError

ap = argparse.ArgumentParser()
ap.add_argument("template", help="template spectrum file or @list of files to be stacked")
ap.add_argument("filelist", metavar="file", nargs="+", help="spectrum file or @list of files to be stacked")
ap.add_argument("-c", metavar="ra dec pmra pmdec plx vrad epoch", type=str, help="specify catalogue information for target")
ap.add_argument("-E", action="store_true", help="don't remove emission lines from spectrum")
ap.add_argument("-f", action="store_true", help="save correlations and LSDs to ascii files")
ap.add_argument("-o", type=int, help="override order number used for analysis")
ap.add_argument("-p", help="set plot filename suffix (default inone)")
ap.add_argument("-R", action="store_true", help="don't use cosmic rejection when stacking")
ap.add_argument("-r", metavar="vsini", type=float, help="use fixed rotation for final velocity analysis")
ap.add_argument("-w", type=get_range, default=None, help="set wavelength range -w='min:max' in A")

if len(sys.argv) == 1:
  ap.print_help(sys.stderr)
  sys.exit(1)

args = ap.parse_args()

# New read_spec structure.
rs = read_spec()

if args.o is not None:
  rs.overrideorder = args.o

if args.p is not None:
  plotsuffix = args.p
else:
  plotsuffix = "inone"
  
emchop = not args.E

# Read template.
tmplsp = rs.read_spec(args.template, istmpl=1, wantstruct=True, doreject=not args.R)

if args.w is not None:  # manual wavelength range
  tmplsp.msk = numpy.logical_and(tmplsp.wave >= args.w[0], tmplsp.wave < args.w[1])

tmplname = stripname(args.template)

filelist = args.filelist
nf = len(args.filelist)

donehdr = False

for ifile, filename in enumerate(filelist):
  sp = rs.read_spec(filename, src_str=args.c, wantstruct=True, doreject=not args.R)

  if args.w is not None:  # manual wavelength range
    sp.msk = numpy.logical_and(sp.wave >= args.w[0], sp.wave < args.w[1])

  basefile = stripname(filename)
  outfile = basefile + "_" + plotsuffix + ".pdf"

  pdf = PdfPages(outfile)

  # Velocity.
  vrad, hbest, vsinibest = do_vrad(pdf, tmplname, filename,
                                   tmplsp.mbjd, tmplsp.wave, tmplsp.flux, tmplsp.e_flux, tmplsp.msk, tmplsp.vbcv, tmplsp.vrad,
                                   sp.mbjd, sp.wave, sp.flux, sp.e_flux, sp.msk, sp.vbcv,
                                   order=rs.singleorder, dovsini=0, savefile=args.f, emchop=emchop)

  # LSD.
  lsdorders = rs.multiorder
  if args.o is not None:
    lsdorders = [ args.o ]

  do_lsd(pdf, filename,
         tmplsp.mbjd, tmplsp.wave, tmplsp.flux, tmplsp.e_flux, tmplsp.msk, tmplsp.vbcv, tmplsp.vrad,
         sp.mbjd, sp.wave, sp.flux, sp.e_flux, sp.msk, sp.vbcv,
         vrad,
         orders=lsdorders, savefile=args.f, emchop=emchop)

  if args.r is None:
    vrad, hbest, vsinibest = do_vrad(pdf, tmplname, filename,
                                     tmplsp.mbjd, tmplsp.wave, tmplsp.flux, tmplsp.e_flux, tmplsp.msk, tmplsp.vbcv, tmplsp.vrad,
                                     sp.mbjd, sp.wave, sp.flux, sp.e_flux, sp.msk, sp.vbcv,
                                     order=rs.singleorder, dovsini=1, savefile=args.f, emchop=emchop)
  elif args.r > 0:
    vrad, hbest, vsinibest = do_vrad(pdf, tmplname, filename,
                                     tmplsp.mbjd, tmplsp.wave, tmplsp.flux, tmplsp.e_flux, tmplsp.msk, tmplsp.vbcv, tmplsp.vrad,
                                     sp.mbjd, sp.wave, sp.flux, sp.e_flux, sp.msk, sp.vbcv,
                                     order=rs.singleorder, dovsini=0, vsini=args.r, savefile=args.f, emchop=emchop)

  if sp.blaze is not None:
    # Equivalent width.
    ewha, e_ewha, ewcai, e_ewcai = do_ews(tmplsp.wave, tmplsp.flux, tmplsp.e_flux, tmplsp.msk, tmplsp.blaze,
                                          sp.wave, sp.flux, sp.e_flux, sp.msk, sp.blaze,
                                          vrad)
  else:
    ewha = None
    e_ewha = None
    ewcai = None
    e_ewcai = None

  # Plot features.
  do_plotfeatures(pdf, filename,
                  sp.mbjd, sp.wave, sp.flux, sp.e_flux, sp.msk,
                  vrad, ewha, e_ewha, ewcai, e_ewcai)

  pdf.close()

  if not donehdr:
    ss = "#Filename"
    lf = len(ss)
    lt = len(filename)

    if lt > lf:
      ss += " " * (lt-lf)

    print(ss + " BJD-TDB      vrad     h         exptime vsini ewha     e_ewha ewca    e_ewca")
    
    donehdr = True
  
  if ewha is not None and ewcai is not None:
    print("{0:s} {1:12.4f} {2:8.3f} {3:9.6f} {4:7.2f} {5:5.2f} {6:8.4f} {7:6.4f} {8:7.4f} {9:6.4f}".format(filename, lfa.ZMJD+sp.mbjd, vrad, hbest, sp.exptime, vsinibest, ewha, e_ewha, ewcai, e_ewcai))
  else:
    print("{0:s} {1:12.4f} {2:8.3f} {3:9.6f} {4:7.2f} {5:5.2f}".format(filename, lfa.ZMJD+sp.mbjd, vrad, hbest, sp.exptime, vsinibest))

sys.exit(0)

