<file_sep>import sys
import os
import numpy as np
import matplotlib.pyplot as plt

alpha1,alpha2,alpha3 = np.loadtxt('ALPHA_AMP_values_V.txt',usecols=(1,2,3)).T
kalpha1,kalpha2,kalpha3 = np.loadtxt('Kate_ALPHA_AMP_values_V.txt',usecols=(1,2,3)).T
print '#Inf in Alpha 1: {}'.format(len(np.isnan(alpha1)==0))
print '#Inf in Alpha 2: {}'.format(len(np.isnan(alpha2)==0))
print '#Inf in Alpha 3: {}'.format(len(np.isnan(alpha3==0)))
print '#Inf in Kate Alpha 1: {}'.format(len(np.isnan(kalpha1)==0))
print '#Inf in Kate Alpha 2: {}'.format(len(np.isnan(kalpha2)==0))
print '#Inf in Kate Alpha 3: {}'.format(len(np.isnan(kalpha3)==0))
print len(np.isnan(alpha2) )
#alpha1 = alpha1[~np.isnan(alpha2)]
#alpha2 = alpha2[~np.isnan(alpha2)]
#alpha3 = alpha3[~np.isnan(alpha2)]
#
xx=np.where(np.isnan(alpha2))
alpha1=np.delete(alpha1,xx)
alpha2=np.delete(alpha2,xx)
alpha3=np.delete(alpha3,xx)
kalpha1=np.delete(kalpha1,xx)
kalpha2=np.delete(kalpha2,xx)
kalpha3=np.delete(kalpha3,xx)

#kalpha1 = alpha1[~np.isnan(alpha2)]
#kalpha2 = alpha2[~np.isnan(alpha2)]
#kalpha3 = alpha3[~np.isnan(alpha2)]

n_bins = 50
fig,(ax0,ax1,ax2) = plt.subplots(1,3,figsize=(15,8))
n1, bins1, patches1 = ax0.hist(kalpha1, n_bins,color='red', normed=0, histtype='step', label='Kate ALPHA')
n, bins, patches = ax0.hist(alpha1, bins1,color='blue', normed=0, histtype='step', label='My ALPHA')

nn1, nbins1, npatches1 = ax1.hist(kalpha2, n_bins,color='red', normed=0, histtype='step', label='Kate ALPHA')
nn, nbins, npatches = ax1.hist(alpha2, nbins1,color='blue', normed=0, histtype='step', label='My ALPHA')

nnn1, nnbins1, nnpatches1 = ax2.hist(kalpha3, n_bins,color='red', normed=0, histtype='step', label='Kate ALPHA')
nnn, nnbins, nnpatches = ax2.hist(alpha3, nnbins1,color='blue', normed=0, histtype='step', label='My ALPHA')

ax0.set_xlabel('ALPHA 1')
ax1.set_xlabel('ALPHA 2')
ax2.set_xlabel('ALPHA 3')
ax0.set_ylabel('Histogram Density')
ax1.set_ylabel('Histogram Density')
ax2.set_ylabel('Histogram Density')
ax0.legend(loc=2)
ax0.set_title('Comparison of Kate\'s and Vivek\'s Continuum Fitting ')
fig.tight_layout()
fig.savefig('Continuumfitting_comparison_V.png')
plt.show()


<file_sep>import sys
import os
import numpy as np
from astropy import io
from astropy.io import fits 
from astropy.io import ascii
from matplotlib import pyplot as plt
from astropy import constants as const
from astropy import units as U
from astropy.coordinates import SkyCoord
from dustmaps.sfd import SFDQuery
from specutils import extinction  
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit 
from matplotlib.backends.backend_pdf import PdfPages
from lmfit import minimize, Parameters
from astropy.modeling.models import Voigt1D

def dered_flux(Av,wave,flux):
    dered_flux = np.zeros(len(wave))
    for i in range(len(wave)):
        ext_lambda= extinction.extinction_ccm89(wave[i] * U.angstrom, a_v= Av, r_v= 3.1)
        tau_lambda= ext_lambda/(1*1.086) 
        dered_flux[i]= flux[i] * np.exp(tau_lambda)
    return dered_flux

def myVoigt(x, amp, center, fwhm_l, fwhm_g, scale, alpha):
    v1= Voigt1D(x_0=center, amplitude_L=amp, fwhm_L=fwhm_l, fwhm_G=fwhm_g)
    powerlaw = scale*x**alpha
    voigt = v1(x)
    voigt_tot = (voigt+powerlaw)
    return voigt_tot

def myGaussHermite(x, amp, center, sig, skew, kurt, scale, alpha):
    c1=-np.sqrt(3); c2=-np.sqrt(6); c3=2/np.sqrt(3); c4=np.sqrt(6)/3; c5=np.sqrt(6)/4
    gausshermite = amp*np.exp(-.5*((x-center)/sig)**2)*(1+skew*(c1*((x-center)/sig)+c3*((x-center)/sig)**3)+kurt*(c5+c2*((x-center)/sig)**2+c4*((x-center)/sig)**4))
    powerlaw = scale*x**alpha
    gaustot_gh = (gausshermite+powerlaw)
    return gaustot_gh


def myDoubleGauss(x, amp1, center1, sig1, amp2, center2, sig2, scale, alpha):
    gaus1=amp1*np.exp(-.5*((x-center1)/sig1)**2)
    gaus2=amp2*np.exp(-.5*((x-center2)/sig2)**2)
    powerlaw = scale*x**alpha
    gaustot_2g= (gaus1+gaus2+powerlaw)
    
    return gaustot_2g

def ext_coeff(lamb):
    inv_lamba=[0.45,0.61,0.8,1.82,2.27,2.7,3.22,3.34,3.46,3.6,3.75,3.92,4.09,4.28,4.50,4.73,5.00,5.24,5.38,5.52,5.70,5.88,6.07,6.27,6.48,6.72,6.98,7.23,7.52,7.84]
    smc_ext=[-2.61,-2.47,-2.12,0.0,1.0,1.67,2.29,2.65,3.0,3.15,3.49,3.91,4.24,4.53,5.3,5.85,6.38,6.76,6.9,7.17,7.71,8.01,8.49,9.06,9.28,9.84,10.8,11.51,12.52,13.54]
    xy=np.interp(1.0/(lamb*10**(-4)),inv_lamba,smc_ext)
    ext_lamb=(xy+2.98)/3.98 # Rv=2.98
    #return(z,ext_lamb)
    return(ext_lamb)

def myfunct(wave, a , b):
    return a*ext_coeff(wave)*wave**(b)


def waveRange(wave):
    #wav_range=[[1300.,1350.],[1420.,1470.],[1700.,1800.],[2080.,2215],[2480.,2655],[3225.,3900.],[4200.,4230.],[4435.,4700.],[5200.,5700.]]
    wav_range= [(1400,1850)]


    mwav_range = []
    for rr in wav_range:
        if ((np.min(wave) < rr[0]) & (np.max(wave) > rr[1])):
            nrr0 = rr[0] ; nrr1 = rr[1]
        elif ((np.min(wave) > rr[0]) & (np.min(wave) < rr[1])& (np.max(wave) < rr[1])):
            nrr0 = np.min(wave) ; nrr1 = rr[1]
        elif ((np.min(wave) < rr[0]) & (np.max(wave) > rr[0]) &(np.max(wave) < rr[1])):
            nrr0 = rr[0] ; nrr1 = np.max(wave)
        else :
            continue
        mwav_range.append([nrr0,nrr1])
    return mwav_range

def compute_alpha_voigt(wl, spec, ivar, wav_range, per_value=[1,99.5]):
    print ' Voigt Routine begins'
    spec[np.isnan(spec)] = 0
    ivar[np.isnan(ivar)] = 0
    print ivar
    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])
    #plt.plot(wl,spec)
    for j in range(len(wav_range)):
        #print wav_range[j]
        #print min(wl),max(wl)
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        #print wl[temp],len(spec),len(wl),spec[temp],ivar[temp]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
        #print len(tempspec)
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut[0]) & (tempspec < cut[1])  & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))
    print 'voigt Debug',len(wavelength)
    pv0=[1.0, 1545., 50., 8.0,1.0,1.0]
    param_boundsv = ([0,1540,0.0,0.0,0.0,-np.inf],[np.inf,1550,np.inf,20,np.inf,np.inf])

    try:
        #plt.plot(wavelength,spectra)
        #plt.show()
        poptv, pcovv = curve_fit(myVoigt, wavelength, spectra, pv0, sigma=1.0/np.sqrt(invar),bounds=param_boundsv)
    except (RuntimeError, TypeError):
        AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv, CHISQv, DOFv = np.nan, np.nan, np.nan, np.nan,np.nan, np.nan, np.nan,  np.nan
    else:
        AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv = poptv[0], poptv[1], poptv[2], poptv[3], poptv[4], poptv[5]
        CHISQv = np.sum(invar * (spectra - myVoigt(wavelength, poptv[0], poptv[1], poptv[2], poptv[3], poptv[4], poptv[5]))**2)
        # DOF = N - n  , n = 2
        DOFv = len(spectra) - 6
    print 'Voight Routine ends' 
    print 'Compute Alpha:',  AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv, CHISQv, DOFv
    return  AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv, CHISQv, DOFv



def compute_alpha_gh(wl, spec, ivar, wav_range, per_value=[1,99.5]):
    print ' GH Routine begins'
    spec[np.isnan(spec)] = 0
    ivar[np.isnan(ivar)] = 0
    print ivar
    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])
    #plt.plot(wl,spec)
    for j in range(len(wav_range)):
        #print wav_range[j]
        #print min(wl),max(wl)
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        #print wl[temp],len(spec),len(wl),spec[temp],ivar[temp]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
        print len(tempspec)
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut[0]) & (tempspec < cut[1])  & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))
    print '2H Debug',len(wavelength)
    pgh0=[1.0, 1545., 8., 0.2, 0.4,1.0,1.0]
    param_boundsgh = ([0,1540,-20,-np.inf,-np.inf,0,-np.inf],[np.inf,1555,20,np.inf,np.inf,np.inf,np.inf])
    try:
        #plt.plot(wavelength,spectra)
        #plt.show()
        poptgh, pcovgh = curve_fit(myGaussHermite, wavelength, spectra, pgh0, sigma=1.0/np.sqrt(invar),bounds=param_boundsgh)
    except (RuntimeError, TypeError):
        AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh, SCALEgh, ALPHAgh, CHISQgh, DOFgh = np.nan, np.nan, np.nan, np.nan,np.nan, np.nan, np.nan,  np.nan, np.nan
    else:
        AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh, SCALEgh, ALPHAgh = poptgh[0], poptgh[1], poptgh[2], poptgh[3], poptgh[4], poptgh[5], poptgh[6]
        CHISQgh = np.sum(invar * (spectra - myGaussHermite(wavelength, poptgh[0], poptgh[1], poptgh[2], poptgh[3], poptgh[4], poptgh[5], poptgh[6]))**2)
        # DOF = N - n  , n = 2
        DOFgh = len(spectra) - 7
    print 'GH Routine ends' 
    print 'Compute Alpha:', AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh,SCALEgh, ALPHAgh, CHISQgh, DOFgh
    return AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh,SCALEgh, ALPHAgh, CHISQgh, DOFgh


def compute_alpha_2g(wl, spec, ivar, wav_range, per_value=[1,99.5]):
    print '2G Routine begins'
    spec[np.isnan(spec)] = 0
    ivar[np.isnan(ivar)] = 0

    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])
    #plt.plot(wl,spec)
    for j in range(len(wav_range)):
        #print wav_range[j]
        #print min(wl),max(wl)
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        #print wl[temp],len(spec),len(wl),spec[temp],ivar[temp]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
        print len(tempspec)
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut[0]) & (tempspec < cut[1])  & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))
    print '2G Debug',len(wavelength)
    pg20=[1.0, 1545., 8., 1.0, 1550., 54.,1.0,1.0]
    param_bounds2g = ([0,1540,0,0,1540,-np.inf,-np.inf,-np.inf],[np.inf,1555,np.inf,np.inf,1555,np.inf,np.inf,np.inf])
    try:
        #plt.plot(wavelength,spectra)
        #plt.show()
        popt2g, pcov2g = curve_fit(myDoubleGauss, wavelength, spectra, pg20, sigma=1.0/np.sqrt(invar),bounds=param_bounds2g)
    except (RuntimeError, TypeError):
        AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g, CHISQ2g, DOF2g = np.nan, np.nan, np.nan, np.nan,np.nan, np.nan, np.nan, np.nan, np.nan, np.nan
    else:
        AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g = popt2g[0], popt2g[1], popt2g[2], popt2g[3], popt2g[4], popt2g[5], popt2g[6], popt2g[7]
        CHISQ2g = np.sum(invar * (spectra - myDoubleGauss(wavelength, popt2g[0], popt2g[1], popt2g[2], popt2g[3], popt2g[4], popt2g[5], popt2g[6], popt2g[7]))**2)
        # DOF = N - n  , n = 2
        DOF2g = len(spectra) - 8 
    print '2G Routine ends' 
    print 'Compute Alpha:' ,AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g, CHISQ2g, DOF2g
    return AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g, CHISQ2g, DOF2g 


def maskOutliers_v(wave, flux, weight, amp, center, sigmal, sigmag, scale, alpha):
    model = myVoigt(wave,amp, center, sigmal, sigmag,  scale, alpha)
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    print 'Weights Before mask',weight
    #ww = np.where (np.abs(fluxdiff) > 3*std)
    ww = np.where (((np.abs(fluxdiff) > 2*std) & (flux <= np.median(flux))) | ((np.abs(fluxdiff) > 4*std) & (flux >= np.median(flux))))
    #nwave = np.delete(wave,ww)
    #nflux = np.delete(flux,ww)
    weight[ww] = 0#np.delete(weight,ww)
    print 'Weights After mask',weight
    return wave,flux,weight


def maskOutliers_gh(wave, flux, weight, amp, center, sigma, skew, kurt, scale, alpha):
    model = myGaussHermite(wave,amp, center, sigma, skew, kurt, scale, alpha)
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    print 'Weights Before mask',weight
    #ww = np.where (np.abs(fluxdiff) > 3*std)
    ww = np.where (((np.abs(fluxdiff) > 2*std) & (flux <= np.median(flux))) | ((np.abs(fluxdiff) > 4*std) & (flux >= np.median(flux))))
    #nwave = np.delete(wave,ww)
    #nflux = np.delete(flux,ww)
    weight[ww] = 0#np.delete(weight,ww)
    print 'Weights After mask',weight
    return wave,flux,weight


def maskOutliers_2g(wave,flux,weight,amp1,center1,sigma1,amp2,center2,sigma2,scale,alpha):
    model = myDoubleGauss(wave,amp1,center1,sigma1,amp2,center2,sigma2,scale,alpha)
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    #ww = np.where (np.abs(fluxdiff) > 3*std)
    ww = np.where (((np.abs(fluxdiff) > 2*std) & (flux <= np.median(flux))) | ((np.abs(fluxdiff) > 4*std) & (flux >= np.median(flux))))
    #nwave = np.delete(wave,ww)
    #nflux = np.delete(flux,ww)
    weight[ww] = 0#np.delete(weight,ww)
    
    return wave,flux,weight


def minmax(vec):
    return min(vec),max(vec)

#Version _ Tried small sigma values for both Gaussians in GausDouble and GausHermite
#Version II Tried making one of the sigmas high, one 1 and two 54 GausDouble; Same as before for GausHermite
#Version III Changing both the central wavelength to 1550 for GausDouble
df = np.genfromtxt('tdss_allmatches_crop_edit.dat',names=['ra','dec','z','pmf1','pmf2','pmf3'],dtype=(float,float,float,'|S15','|S15','|S15'))
wav_range= [(1280,1350),(1700,1800),(1950,2200),(2650,2710),(3010,3700),(3950,4050),(4140,4270)]
pp = PdfPages('ContinuumNormalization_plus_EmissionLineFits_IV.pdf')
#fx=open('ALPHA_AMP_values_V.txt','w')
#fx1=open('Kate_ALPHA_AMP_values_V.txt','w')
#for i in range(len(df['pmf1'])):
for i in range(len(df)):
#for i in range(12):

    print 'Norm_Spectra/Normspec_'+df['pmf1'][i]+'.txt' 
    print 'Norm_Spectra/Normspec_'+df['pmf2'][i]+'.txt' 
    print 'Norm_Spectra/Normspec_'+df['pmf3'][i]+'.txt' 
    if ((not os.path.isfile('Norm_Spectra/Normspec_'+df['pmf1'][i]+'.txt')) | (not os.path.isfile('Norm_Spectra/Normspec_'+df['pmf2'][i]+'.txt')) | ( not os.path.isfile('Norm_Spectra/Normspec_'+df['pmf3'][i]+'.txt'))):
        continue
    else:
        data1 = np.loadtxt('Norm_Spectra/Normspec_'+df['pmf1'][i]+'.txt')
        data2 = np.loadtxt('Norm_Spectra/Normspec_'+df['pmf2'][i]+'.txt')
        data3 = np.loadtxt('Norm_Spectra/Normspec_'+df['pmf3'][i]+'.txt')
    wave1 = data1.T[0] ;    flux1 = data1.T[1] ; weight1 =  data1.T[2]; mask1 = data1.T[3]
    wave2 = data2.T[0] ;    flux2 = data2.T[1] ; weight2 =  data2.T[2]; mask2 = data2.T[3]
    wave3 = data3.T[0] ;    flux3 = data3.T[1] ; weight3 =  data3.T[2]; mask3 = data3.T[3]
    print len(wave1),len(flux1),len(weight1)
    #print weight1
    #Initial Mask for broad absorptions
    cut1 = np.nanpercentile(flux1, 20)
    cut2 = np.nanpercentile(flux2, 20)
    cut3 = np.nanpercentile(flux3, 20)
    xx1 = np.where(flux1 < cut1)[0]
    xx2 = np.where(flux2 < cut2)[0]
    xx3 = np.where(flux3 < cut3)[0]
    print cut1
    print flux1[xx1]
    weight1[xx1] = 0
    weight2[xx2] = 0
    weight3[xx3] = 0
    
        
    #Fit powerlaw iterate and mask outliers
    iteration = 3
    for j in range(iteration):
        if j == 0:
            print 'iteration 1 Begins'
            AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1, CHISQgh1, DOFgh1   = compute_alpha_gh(wave1, flux1, weight1,  waveRange(wave1))
            AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 ,CHISQ2g1, DOF2g1   = compute_alpha_2g(wave1, flux1, weight1,  waveRange(wave1))
            AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1, CHISQv1, DOFv1   = compute_alpha_voigt(wave1, flux1, weight1,  waveRange(wave1))
            print 'Object 1 Done'
            AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2, CHISQgh2, DOFgh2   = compute_alpha_gh(wave2, flux2, weight2,  waveRange(wave2))
            AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 ,CHISQ2g2, DOF2g2   = compute_alpha_2g(wave2, flux2, weight2,  waveRange(wave2))
            AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2, CHISQv2, DOFv2   = compute_alpha_voigt(wave2, flux2, weight2,  waveRange(wave2))
            print 'Object 2 Done'
            AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3, CHISQgh3, DOFgh3   = compute_alpha_gh(wave3, flux3, weight3,  waveRange(wave3))
            AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 ,CHISQ2g3, DOF2g3   = compute_alpha_2g(wave3, flux3, weight3,  waveRange(wave3))
            AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3, CHISQv3, DOFv3   = compute_alpha_voigt(wave3, flux3, weight3,  waveRange(wave3))
            nwavegh1 = wave1; nfluxgh1=flux1;nweightgh1 = weight1
            nwavegh2 = wave2; nfluxgh2=flux2;nweightgh2 = weight2
            nwavegh3 = wave3; nfluxgh3=flux3;nweightgh3 = weight3

            nwave2g1 = wave1; nflux2g1=flux1;nweight2g1 = weight1
            nwave2g2 = wave2; nflux2g2=flux2;nweight2g2 = weight2
            nwave2g3 = wave3; nflux2g3=flux3;nweight2g3 = weight3


            nwavev1 = wave1; nfluxv1=flux1;nweightv1 = weight1
            nwavev2 = wave2; nfluxv2=flux2;nweightv2 = weight2
            nwavev3 = wave3; nfluxv3=flux3;nweightv3 = weight3
            print 'iteration 1 Ends'
            continue
        else:
            nwavegh1,nfluxgh1,nweightgh1 = maskOutliers_gh(nwavegh1, nfluxgh1, nweightgh1, AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1)
            nwavegh2,nfluxgh2,nweightgh2 = maskOutliers_gh(nwavegh2, nfluxgh2, nweightgh2, AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2)
            nwavegh3,nfluxgh3,nweightgh3 = maskOutliers_gh(nwavegh3, nfluxgh3, nweightgh3, AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3)
            
            nwave2g1,nflux2g1,nweight2g1 = maskOutliers_2g(nwave2g1, nflux2g1, nweight2g1, AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1)
            nwave2g2,nflux2g2,nweight2g2 = maskOutliers_2g(nwave2g2, nflux2g2, nweight2g2, AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2)
            nwave2g3,nflux2g3,nweight2g3 = maskOutliers_2g(nwave2g3, nflux2g3, nweight2g3, AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3)

            nwavev1,nfluxv1,nweightv1 = maskOutliers_v(nwavev1, nfluxv1, nweightv1, AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1 )
            nwavev2,nfluxv2,nweightv2 = maskOutliers_v(nwavev2, nfluxv2, nweightv2, AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2 )
            nwavev3,nfluxv3,nweightv3 = maskOutliers_v(nwavev3, nfluxv3, nweightv3, AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3 )

            
            AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1, CHISQgh1, DOFgh1   = compute_alpha_gh(nwavegh1, nfluxgh1, nweightgh1,  waveRange(nwavegh1))
            AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 ,CHISQ2g1, DOF2g1   = compute_alpha_2g(nwave2g1, nflux2g1, nweight2g1,  waveRange(nwave2g1))
            AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1, CHISQv1, DOFv1   = compute_alpha_voigt(nwavev1, nfluxv1, nweightv1,  waveRange(nwavev1))
            AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2, CHISQgh2, DOFgh2   = compute_alpha_gh(nwavegh2, nfluxgh2, nweightgh2,  waveRange(nwavegh2))
            AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 ,CHISQ2g2, DOF2g2   = compute_alpha_2g(nwave2g2, nflux2g2, nweight2g2,  waveRange(nwave2g2))
            AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2, CHISQv2, DOFv2   = compute_alpha_voigt(nwavev2, nfluxv2, nweightv2,  waveRange(nwavev2))
            AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3, CHISQgh3, DOFgh3   = compute_alpha_gh(nwavegh3, nfluxgh3, nweightgh3,  waveRange(nwavegh3))
            AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 ,CHISQ2g3, DOF2g3   = compute_alpha_2g(nwave2g3, nflux2g3, nweight2g3,  waveRange(nwave2g3))
            AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3, CHISQv3, DOFv3   = compute_alpha_voigt(nwavev3, nfluxv3, nweightv3,  waveRange(nwavev3))



    fig,((ax1,rax1),(ax2,rax2),(ax3,rax3))=plt.subplots(3,2,figsize=(20,10))
    ax1.plot(wave1,flux1,label=str(df['pmf1'][i]))
    ax1.plot(wave1,weight1,alpha=0.2)
    ax1.plot(wave1[weight1 > 0],flux1[weight1>0],'.')
    ax2.plot(wave2,flux2,label=str(df['pmf2'][i]))
    ax2.plot(wave2,weight2,alpha=0.2)
    ax2.plot(wave2[weight2 > 0],flux2[weight2>0],'.')

    ax3.plot(wave3,flux3,label=str(df['pmf3'][i]))
    ax3.plot(wave3,weight3,alpha=0.2)
    ax3.plot(wave3[weight3 > 0],flux3[weight3>0],'.')

    #ax1.plot(cwave1,cflux1,'--')
    pgh1=ax1.plot(wave1,myGaussHermite(wave1,AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1),':',color='red',lw=3,label='GaussHermite+Power law')
    pgh2=ax2.plot(wave2,myGaussHermite(wave2,AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2),':',color='red',lw=3,label='GaussHermite+Power law')
    pgh3=ax3.plot(wave3,myGaussHermite(wave3,AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3),':',color='red',lw=3,label='GaussHermite+Power law')
    
    p2g1=ax1.plot(wave1,myDoubleGauss(wave1,AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 ),':',color='blue',lw=3,label='Double Gaussian + Power law')
    p2g2=ax2.plot(wave2,myDoubleGauss(wave2,AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 ),':',color='blue',lw=3,label='Double Gaussian + Power law')
    p2g3=ax3.plot(wave1,myDoubleGauss(wave3,AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 ),':',color='blue',lw=3,label='Double Gaussian + Power law')
    
    p2v1 = ax1.plot(wave1,myVoigt(wave1,AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1),':',color='orange',lw=3,label='Voigt+ Power law')
    p2v2 = ax2.plot(wave2,myVoigt(wave2,AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2),':',color='orange',lw=3,label='Voigt+ Power law')
    p2v3 = ax3.plot(wave3,myVoigt(wave3,AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3),':',color='orange',lw=3,label='Voigt+ Power law')
    ax1.axhline(cut1,ls=':')
    ax2.axhline(cut1,ls=':')
    ax3.axhline(cut1,ls=':')
    string1gh = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPgh1,CENTERgh1,(CHISQgh1/DOFgh1))
    string12g = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPa2g1,CENTERa2g1,(CHISQ2g1/ DOF2g1))
    string2gh = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPgh2,CENTERgh2,(CHISQgh2/ DOFgh2))
    string22g = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPa2g2,CENTERa2g2,(CHISQ2g2/ DOF2g2))
    string3gh = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPgh3,CENTERgh3,(CHISQgh3/ DOFgh3))
    string32g = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPa2g3,CENTERa2g3,(CHISQ2g3/ DOF2g3))
    ax1.set_xlim(np.min(wave1),np.max(wave1))
    ax2.set_xlim(np.min(wave1),np.max(wave1))
    ax3.set_xlim(np.min(wave1),np.max(wave1))
    xlim=ax1.get_xlim()
    ylim1 = ax1.get_ylim()
    ylim2 = ax1.get_ylim()
    ylim3 = ax1.get_ylim()
    ax1.set_ylim(-0.05,4)
    ax2.set_ylim(-0.05,4)
    ax3.set_ylim(-0.05,4)
    ax1.text(xlim[1]-0.15*(xlim[1] - xlim[0]),ylim1[0] + 0.1*(ylim1[1] - ylim1[0]),str(df['pmf1'][i]))
    ax2.text(xlim[1]-0.15*(xlim[1] - xlim[0]),ylim1[0] + 0.1*(ylim1[1] - ylim1[0]),str(df['pmf2'][i]))
    ax3.text(xlim[1]-0.15*(xlim[1] - xlim[0]),ylim1[0] + 0.1*(ylim1[1] - ylim1[0]),str(df['pmf3'][i]))
    #ax1.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim1[1] - 0.2*(ylim1[1] - ylim1[0]),string12g)
    #ax2.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim2[1] - 0.1*(ylim1[1] - ylim2[0]),string2gh)
    #ax2.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim2[1] - 0.2*(ylim1[1] - ylim2[0]),string22g)
    #ax3.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim3[1] - 0.1*(ylim1[1] - ylim3[0]),string3gh)
    #ax3.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim3[1] - 0.2*(ylim1[1] - ylim2[0]),string32g)
    ax1.legend([pgh1,p2g1,p2v1],['Gaus-Hermite','2-Gaus','Voigt'],loc=3)
    ax1.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f' \
        %(AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1),xy=(.05,.95), \
        xycoords='axes fraction',ha="left", va="top", \
    bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax1.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f' \
        %(AMPa2g1, AMPb2g1, CENTERa2g1, CENTERb2g1, SIGMAa2g1, SIGMAb2g1),xy=(.95,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax1.annotate('Voigt:\nAmp$_1$ = %.2f\nCenter = %.2f\n$\sigma_l$ = %.2f\n$\sigma_g$ = %.2f' \
        %(AMPv1, CENTERv1, SIGMALv1, SIGMAGv1),xy=(.35,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)

    
    ax2.legend([pgh2,p2g2,p2v2],['Gaus-Hermite','2-Gaus','Voigt'],prop={'size':8},loc='center left')
    ax2.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f' \
        %(AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2),xy=(.05,.95), \
        xycoords='axes fraction',ha="left", va="top", \
    bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax2.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f' \
        %(AMPa2g2, AMPb2g2, CENTERa2g2, CENTERb2g2, SIGMAa2g2, SIGMAb2g2),xy=(.95,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax2.annotate('Voigt:\nAmp$_1$ = %.2f\nCenter = %.2f\n$\sigma_l$ = %.2f\n$\sigma_g$ = %.2f' \
        %(AMPv2, CENTERv2, SIGMALv2, SIGMAGv2),xy=(.35,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)


    ax3.legend([pgh3,p2g3,p2v3],['Gaus-Hermite','2-Gaus','Voigt'],prop={'size':8},loc='center left')
    ax3.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f' \
        %(AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3),xy=(.05,.95), \
        xycoords='axes fraction',ha="left", va="top", \
    bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax3.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f' \
        %(AMPa2g3, AMPb2g3, CENTERa2g3, CENTERb2g3, SIGMAa2g3, SIGMAb2g3),xy=(.95,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax3.annotate('Voigt:\nAmp$_1$ = %.2f\nCenter = %.2f\n$\sigma_l$ = %.2f\n$\sigma_g$ = %.2f' \
        %(AMPv3, CENTERv3, SIGMALv3, SIGMAGv3),xy=(.35,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)

    #Save the spectra
    filenamegh1 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf1'][i]+'_Em_gh.txt'
    filenamegh2 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf2'][i]+'_Em_gh.txt'
    filenamegh3 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf3'][i]+'_Em_gh.txt'
    
    filename2g1 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf1'][i]+'_Em_2g.txt'
    filename2g2 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf2'][i]+'_Em_2g.txt'
    filename2g3 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf3'][i]+'_Em_2g.txt'

    filenamev1 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf1'][i]+'_Em_v.txt'
    filenamev2 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf2'][i]+'_Em_v.txt'
    filenamev3 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf3'][i]+'_Em_v.txt'


    # Residuals & Weights
    resgh1 = flux1 / myGaussHermite(wave1,AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1)
    resgh2 = flux2 / myGaussHermite(wave2,AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2)
    resgh3 = flux3 / myGaussHermite(wave3,AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3)
    
    wresgh1 = weight1 / myGaussHermite(wave1,AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1)
    wresgh2 = weight2 / myGaussHermite(wave2,AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2)
    wresgh3 = weight3 / myGaussHermite(wave3,AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3)

    res2g1 = flux1 / myDoubleGauss(wave1,AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 )
    res2g2 = flux2 / myDoubleGauss(wave2,AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 )
    res2g3 = flux3 / myDoubleGauss(wave3,AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 )

    wres2g1 = weight1 / myDoubleGauss(wave1,AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 )
    wres2g2 = weight2 / myDoubleGauss(wave2,AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 )
    wres2g3 = weight3 / myDoubleGauss(wave3,AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 )

    resv1 = flux1 / myVoigt(wave1,AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1)
    resv2 = flux2 / myVoigt(wave2,AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2)
    resv3 = flux3 / myVoigt(wave3,AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3)

    wresv1 = weight1 / myVoigt(wave1,AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1)
    wresv2 = weight2 / myVoigt(wave2,AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2)
    wresv3 = weight3 / myVoigt(wave3,AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3)

    #Save Begins
    np.savetxt(filenamegh1,zip(wave1,resgh1,wresgh1,mask1), fmt='%10.5f')
    np.savetxt(filenamegh2,zip(wave2,resgh2,wresgh2,mask2), fmt='%10.5f')
    np.savetxt(filenamegh3,zip(wave3,resgh3,wresgh3,mask3), fmt='%10.5f')

    np.savetxt(filename2g1,zip(wave1,res2g1,wres2g1,mask1), fmt='%10.5f')
    np.savetxt(filename2g2,zip(wave2,res2g2,wres2g2,mask2), fmt='%10.5f')
    np.savetxt(filename2g3,zip(wave3,res2g3,wres2g3,mask3), fmt='%10.5f')

    np.savetxt(filenamev1,zip(wave1,resv1,wresv1,mask1), fmt='%10.5f')
    np.savetxt(filenamev2,zip(wave2,resv2,wresv2,mask2), fmt='%10.5f')
    np.savetxt(filenamev3,zip(wave3,resv3,wresv3,mask3), fmt='%10.5f')

    #Plot Begins
    presgh1,pres2g1,presv1,=rax1.plot(wave1,resgh1,'k--',wave1,res2g1,'r-',wave1,resv1,'b:',alpha=0.7)
    presgh2,pres2g2,presv2,=rax2.plot(wave2,resgh2,'k--',wave2,res2g2,'r-',wave2,resv2,'b:',alpha=0.7)
    presgh3,pres2g3,presv3,=rax3.plot(wave3,resgh3,'k--',wave3,res2g3,'r-',wave3,resv3,'b:',alpha=0.7)

    rax1.set_ylabel('Normalized Flux')
    rax1.set_xlabel('Wavelength ($\AA$)')
    rax1.legend([presgh1,pres2g1,presv1],['Gaus-Hermite','2-Gaus','Voigt'],numpoints=4,prop={'size':8},loc='lower right')
    rax2.set_ylabel('Normalized Flux')
    rax2.set_xlabel('Wavelength ($\AA$)')
    rax2.legend([presgh2,pres2g2,presv2],['Gaus-Hermite','2-Gaus','Voigt'],numpoints=4,prop={'size':8},loc='lower right')
    rax3.set_ylabel('Normalized Flux')
    rax3.set_xlabel('Wavelength ($\AA$)')
    rax3.legend([presgh3,pres2g3,presv3],['Gaus-Hermite','2-Gaus', 'Voigt'],numpoints=4,prop={'size':8},loc='lower right')
    rax1.axhline(1.0,ls=':',color='blue',alpha=0.5)
    rax2.axhline(1.0,ls=':',color='blue',alpha=0.5)
    rax3.axhline(1.0,ls=':',color='blue',alpha=0.5)
    rax1.set_ylim(0,2)
    rax2.set_ylim(0,2)
    rax3.set_ylim(0,2)
    fig.tight_layout()
    fig.savefig(pp,format='pdf')
    #plt.show()
    plt.clf()
pp.close()

<file_sep>import sys
import os
import numpy as np
from astropy import io
from astropy.io import fits 
from astropy.io import ascii
from matplotlib import pyplot as plt
from astropy import constants as const
from astropy import units as U
from astropy.coordinates import SkyCoord
#from dustmaps.sfd import SFDQuery
#from specutils import extinction  
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit 
from matplotlib.backends.backend_pdf import PdfPages
#from lmfit import minimize, Parameters
from astropy.modeling.models import Voigt1D
from collections import Iterable

def flatten(items):
    """Yield items from any nested iterable; see REF."""
    for x in items:
        if isinstance(x, Iterable) and not isinstance(x, (str, bytes)):
            return  flatten(x)
        else:
            return x

def velocity(wave,ion_wave,z):
    #print wave
    #print z,ion_wave
    c = 299792.0# in km/s
    vel =np.zeros(len(wave))
    zabs = np.zeros(len(wave))
    for i in range(len(wave)):
        #print i,wave[i],zabs[i],vel[i]
        zabs[i] = (wave[i]/ion_wave) -1.0
        vel[i] = -((((1+z)/(1+zabs[i]))**2-1.0)/(((1+z)/(1+zabs[i]))**2+1))*c
    return vel



def balBounds(vel,flux,threshold=0.9,min_vel=2000):
  start = -1
  contigous_segments = []
  for idx,x in enumerate(zip(vel,flux)):
    if start < 0 and x[1] <= threshold:
      start = idx
      #print 'start', start
    elif start >= 0 and x[1] >= threshold:
      dur = idx-start
      veldiff = vel[idx] - vel[start]
      #print 'velocity difference: ',dur, veldiff
      #Multiply by -1 to take care of blueshifted velocities
      if veldiff >= min_vel:
        contigous_segments.append((vel[start],vel[start+dur]))
      start = -1
  return contigous_segments


def mergeVelocityBounds(bounds1,bounds2,bounds3):

    masterbounds = bounds1+bounds2+bounds3
    sortedbounds = sorted(masterbounds, key=lambda t: t[0])
    print 'MasterBounds',masterbounds
    print 'SortedBounds',sortedbounds
    merged = []

    for higher in sortedbounds:
        print 'Higher',higher
        print 'M',merged
        if not merged:
            merged.append(higher)
        else:
            lower = merged[-1]
            #print 'To set the merging of adjacent BALs',higher[0] - lower[1]
            # this condition to merge close by BALs
            if higher[0] > lower[1]:
                if higher[0] - lower[1] <= 50:
                
                    merged[-1] = (lower[0], higher[1])  # replace by merged interval
                else:
                    merged.append(higher)
            elif higher[0] <= lower[1]:
                upper_bound = max(lower[1], higher[1])
                merged[-1] = (lower[0], upper_bound)  # replace by merged interval
            else:
                merged.append(higher)
    print 'Merged',merged
    return merged
def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
    The Savitzky-Golay filter removes high frequency noise from data.
    It has the advantage of preserving the original shape and
    features of the signal better than other types of filtering
    approaches, such as moving averages techniques.
    Parameters
    ----------
    y : array_like, shape (N,)
        the values of the time history of the signal.
    window_size : int
        the length of the window. Must be an odd integer number.
    order : int
        the order of the polynomial used in the filtering.
        Must be less then `window_size` - 1.
    deriv: int
        the order of the derivative to compute (default = 0 means only smoothing)
    Returns
    -------
    ys : ndarray, shape (N)
        the smoothed signal (or it's n-th derivative).
    Notes
    -----
    The Savitzky-Golay is a type of low-pass filter, particularly
    suited for smoothing noisy data. The main idea behind this
    approach is to make for each point a least-square fit with a
    polynomial of high order over a odd-sized window centered at
    the point.
    Examples
    --------
    t = np.linspace(-4, 4, 500)
    y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
    ysg = savitzky_golay(y, window_size=31, order=4)
    import matplotlib.pyplot as plt
    plt.plot(t, y, label='Noisy signal')
    plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
    plt.plot(t, ysg, 'r', label='Filtered signal')
    plt.legend()
    plt.show()
    References
    ----------
    .. [1] <NAME>, <NAME>, Smoothing and Differentiation of
       Data by Simplified Least Squares Procedures. Analytical
       Chemistry, 1964, 36 (8), pp 1627-1639.
    .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing
       W.H. Press, <NAME>, <NAME>, <NAME>
       Cambridge University Press ISBN-13: 9780521880688
    """
    import numpy as np
    from math import factorial
    
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except (ValueError, msg):
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    # precompute coefficients
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    # pad the signal at the extremes with
    # values taken from the signal itself
    firstvals = np.asarray(y[0]) - np.abs( y[1:half_window+1][::-1] - np.asarray(y[0]) )
    lastvals = np.asarray(y[-1]) + np.abs(y[-half_window-1:-1][::-1] - np.asarray(y[-1]))
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')

def topHat(mvel1,merged):
    #Plot a top hat function
    mask1 = np.ones(len(mvel1))*(1==0)
    #print 'VelMerged',mvel1
    indices = list()
    for m in merged:
        #print m[0],m[1]
        ix = np.where((mvel1 >= m[0]) & (mvel1 <=m[1]))[0]
        #print len(ix)
        if len(ix)>0:
            indices.append(ix)
    #print len(indices)
    indices = np.array(indices)
    findices = [y for x in indices for y in x]
    #print findices,
    #print len(findices)
    if len(findices) >0:
        mask1[findices] = True

    return mask1

c_light=299792.0 
cw=1549.0
df = np.genfromtxt('tdss_allmatches_crop_edit.dat',names=['ra','dec','z','pmf1','pmf2','pmf3'],dtype=(float,float,float,'|S15','|S15','|S15'))
pp = PdfPages('BALcomponents_search_I.pdf')
fx=open('BALcomponents_output.txt','w')
#for i in range(len(df['pmf1'])):
for i in range(len(df)):
#for i in range(5):

    filename1 = 'Normspec_'+df['pmf1'][i]+'_cEm_2g.txt' 
    filename2 = 'Normspec_'+df['pmf2'][i]+'_cEm_2g.txt' 
    filename3 = 'Normspec_'+df['pmf3'][i]+'_cEm_2g.txt' 
    if ((not os.path.isfile('EmissionLines_Norm_Spectra/Normspec_'+df['pmf1'][i]+'_cEm_2g.txt')) | (not os.path.isfile('EmissionLines_Norm_Spectra/Normspec_'+df['pmf2'][i]+'_cEm_2g.txt')) | ( not os.path.isfile('EmissionLines_Norm_Spectra/Normspec_'+df['pmf3'][i]+'_cEm_2g.txt'))):
        continue
    else:
        data1 = np.loadtxt('EmissionLines_Norm_Spectra/Normspec_'+df['pmf1'][i]+'_cEm_2g.txt')
        data2 = np.loadtxt('EmissionLines_Norm_Spectra/Normspec_'+df['pmf2'][i]+'_cEm_2g.txt')
        data3 = np.loadtxt('EmissionLines_Norm_Spectra/Normspec_'+df['pmf3'][i]+'_cEm_2g.txt')
    wave1 = data1.T[0] ;    bflux1 = data1.T[1] ; weight1 =  data1.T[2]#; mask1 = data1.T[3]
    wave2 = data2.T[0] ;    bflux2 = data2.T[1] ; weight2 =  data2.T[2]#; mask2 = data2.T[3]
    wave3 = data3.T[0] ;    bflux3 = data3.T[1] ; weight3 =  data3.T[2]#; mask3 = data3.T[3]

    if ((len(bflux1[~np.isnan(bflux1)]) < 1) | (len(bflux2[~np.isnan(bflux2)]) < 1) | (len(bflux3[~np.isnan(bflux3)]) < 1)  ) :
        filename1 = 'Normspec_'+df['pmf1'][i]+'_cEm_v.txt' 
        filename2 = 'Normspec_'+df['pmf2'][i]+'_cEm_v.txt' 
        filename3 = 'Normspec_'+df['pmf3'][i]+'_cEm_v.txt' 

        data1 = np.loadtxt('EmissionLines_Norm_Spectra/Normspec_'+df['pmf1'][i]+'_cEm_v.txt')
        data2 = np.loadtxt('EmissionLines_Norm_Spectra/Normspec_'+df['pmf2'][i]+'_cEm_v.txt')
        data3 = np.loadtxt('EmissionLines_Norm_Spectra/Normspec_'+df['pmf3'][i]+'_cEm_v.txt')
        wave1 = data1.T[0] ;    bflux1 = data1.T[1] ; weight1 =  data1.T[2]#; mask1 = data1.T[3]
        wave2 = data2.T[0] ;    bflux2 = data2.T[1] ; weight2 =  data2.T[2]#; mask2 = data2.T[3]
        wave3 = data3.T[0] ;    bflux3 = data3.T[1] ; weight3 =  data3.T[2]#; mask3 = data3.T[3]

    
    #Work on smoothed spectrum to get rid og high frequency noise
    flux1 = savitzky_golay(bflux1, 5, 2)
    flux2 = savitzky_golay(bflux2, 5, 2)
    flux3 = savitzky_golay(bflux3, 5, 2)
    print len(wave1),len(flux1),len(weight1)    
    
    #Convert wavelengths to velocities
    vel1 = velocity(wave1, cw, 0.0) ; vel2 = velocity(wave2, cw, 0.0) ; vel3 = velocity(wave3, cw, 0.0)
    bounds1 = balBounds(vel1,flux1) ; bounds2 = balBounds(vel2,flux2); bounds3 = balBounds(vel3, flux3)    
    print 'bounds1', bounds1
    print 'bounds2', bounds2
    print 'bounds3', bounds3
    
    #Merge the overlapping regions
    merged = mergeVelocityBounds(bounds1,bounds2,bounds3)
    ntroughs = len(merged)
    
    #Check if each component is present in atleast 2 epochs
    #topHat function returns a mask with ones inside the absorption
    mask1 = topHat(vel1,bounds1) ; mask2 = topHat(vel1,bounds2) ; mask3 = topHat(vel1,bounds3)
    for kk,mm in enumerate(merged):
        #print mm
        #print topHat(vel1,[mm])
        #print mask1
        track_absorption = 0
        if np.sum( mask1.any() and topHat(vel1,[mm])) > 0:
            track_absorption += 1
        if np.sum(mask2.any() and topHat(vel1,[mm])) > 0:
            track_absorption += 1
        if np.sum(mask3.any() and topHat(vel1,[mm])) > 0:
            track_absorption += 1
        if track_absorption < 2:
            merged.pop(kk)
    
    maskall = topHat(vel1,merged)
    #Get the min max of each bounds
    ubounds = []
    if (len(bounds1) == len(bounds2) == len(bounds3)):
      for k in range(len(bounds1)):

        uboundsl,uboundsu = min(bounds1[k][0],bounds2[k][0],bounds3[k][0]),max(bounds1[k][1],bounds2[k][1],bounds3[k][1])
        ubounds.append((uboundsl,uboundsu))
    else:
        print 'Manually adjust ',df[i]
    
    #Plot the spectra and check the bounds
    fig, (ax1,ax2,ax3) = plt.subplots(3,1,sharex=True, sharey=True,figsize=(15,8))
    ax1.plot(vel1,bflux1,color='black',alpha=0.3,label=str(df['pmf1'][i]))
    ax1.plot(vel1,flux1,color='magenta',alpha=0.3)
    ax1.plot(vel1,maskall*0.5,color='green')
    for bb in bounds1:
        ax1.axvline(bb[0],ls='--',color='blue',lw=1)
        ax1.axvline(bb[1],ls='--',color='red',lw=1)
    ax2.plot(vel2,bflux2,color='black',alpha=0.3,label=str(df['pmf2'][i]))
    ax2.plot(vel2,flux2,color='magenta',alpha=0.3)
    ax2.plot(vel1,maskall*0.5,color='green')
    for bb in bounds2:
        ax2.axvline(bb[0],ls='--',color='blue',lw=1)
        ax2.axvline(bb[1],ls='--',color='red',lw=1)
    ax3.plot(vel3,bflux3,color='black',alpha=0.3,label=str(df['pmf3'][i]))
    ax3.plot(vel3,flux3,color='magenta',alpha=0.3)
    ax3.plot(vel1,maskall*0.5,color='green')
    for bb in bounds3:
        ax3.axvline(bb[0],ls='--',color='blue',lw=1)
        ax3.axvline(bb[1],ls='--',color='red',lw=1)
    for ubb in ubounds:
        #print 'final bounds',ubounds
        ax1.axvline(ubb[0],ls='-',color='blue',lw=1)
        ax2.axvline(ubb[0],ls='-',color='blue',lw=1)
        ax3.axvline(ubb[0],ls='-',color='blue',lw=1)
        ax1.axvline(ubb[1],ls='-',color='red',lw=1)
        ax2.axvline(ubb[1],ls='-',color='red',lw=1)
        ax3.axvline(ubb[1],ls='-',color='red',lw=1)

    #ax1.set_xlabel(r'Velocity km s$^{-1}$')
    #ax2.set_xlabel(r'Velocity km s$^{-1}$')
    ax3.set_xlabel(r'Velocity km s$^{-1}$')
    
    ax1.axhline(1.0,ls='--',color='orange')
    ax2.axhline(1.0,ls='--',color='orange')
    ax3.axhline(1.0,ls='--',color='orange')

    ax1.axhline(0.9,ls='--',color='gold')
    ax2.axhline(0.9,ls='--',color='gold')
    ax3.axhline(0.9,ls='--',color='gold')

    ax1.axvspan(-30000,1000,color='grey',alpha=0.1)
    ax2.axvspan(-30000,1000,color='grey',alpha=0.1)
    ax3.axvspan(-30000,1000,color='grey',alpha=0.1)
    ax1.set_ylim(0,2)
    ax2.set_ylim(0,2)
    ax3.set_ylim(0,2)
    #ax1.set_ylabel(r'Normalized Flux')
    #ax2.set_ylabel(r'Normalized Flux')
    ax3.set_ylabel(r'Normalized Flux')
    ax1.legend(loc=1)
    ax2.legend(loc=1)
    ax3.legend(loc=1)
    
    for j in range(len(merged)):
        merge = merged[j]
        print>>fx,'{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}'.format(i,filename1,merge[1],merge[0],0,ntroughs, 2000.0000,2000.0000,10.0,10.0)
        print>>fx,'{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}'.format(i,filename2,merge[1],merge[0],0,ntroughs, 2000.0000,2000.0000,10.0,10.0)
        print>>fx,'{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}'.format(i,filename3,merge[1],merge[0],0,ntroughs, 2000.0000,2000.0000,10.0,10.0)
    #plt.show()
    #plt.clf()
    # Some book keeping
    
    fig.tight_layout()
    fig.savefig(pp,format='pdf')
   # sinp = raw_input('Type Some Key to Continue: ')
pp.close()
fx.close()
<file_sep>#!/usr/bin/python
#-*- coding: utf-8 -*-

'''
Code for continuum normalization of SDSS spectra
@author   : <NAME>
@date     : 21-Jan-2018
@version  : 1.0
'''



import sys
import os
import numpy as np
from astropy import io
from astropy.io import fits 
from astropy.io import ascii
from matplotlib import pyplot as plt
from astropy import constants as const
from astropy import units as U
from astropy.coordinates import SkyCoord
from dustmaps.sfd import SFDQuery
from specutils import extinction  
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
import scipy.optimize as optimization
from scipy import optimize
from scipy.optimize import curve_fit

def maskOutliers(wave,flux,weight,popt):
    model = powerlawFunc(wave,popt[0],popt[1],popt[2])
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    ww = np.where (np.abs(fluxdiff) > 1*std)
    nwave = np.delete(wave,ww)
    nflux = np.delete(flux,ww)
    nweight = np.delete(weight,ww)
    
    return nwave,nflux,nweight

def contWave(wave):
    linefree_regions = [(1250,1350),(1700,1800),(1950,2200),(2650,2710),(2950,3700),(3950,4050)]
    finalcond = False
    for lfr in linefree_regions:
        cond = ((wave >= lfr[0]) & (wave <= lfr[1]))
        finalcond = finalcond | cond
    indices = np.where(finalcond)
    return indices


def powerlawFunc(xdata, amp,index):
    return  amp*xdata**index


def fitPowerlaw(wave,flux,sigma,p):
    p0 = p
    popt,pcov =  optimization.curve_fit(powerlawFunc,wave,flux,p0,sigma)
    return popt,pcov

def myfunct(wave, a , b):
    return a*wave**(b)


def compute_alpha(wl, spec, ivar, wav_range, per_value):
    # print 'Routine begins'
    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])

    for j in range(len(wav_range)):
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
            
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut) & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))

    try:
        popt, pcov = curve_fit(myfunct, wavelength, spectra, sigma=1.0/np.sqrt(invar))
    except (RuntimeError, TypeError):
        AMP, ALPHA, CHISQ, DOF = np.nan, np.nan, np.nan, np.nan
    else:
        AMP, ALPHA = popt[0], popt[1]
        CHISQ = np.sum(invar * (spectra - myfunct(wavelength, popt[0], popt[1]))**2)
        # DOF = N - n  , n = 2
        DOF = len(spectra) - 2
    # print 'Routine ends' 

    return AMP, ALPHA, CHISQ, DOF

per_value= 0.2
wav_range=[[1300.,1350.],[1420.,1470.],[1700.,1800.],[2080.,2215],[2480.,2655],[3225.,3900.],[4200.,4230.],[4435.,4700.],[5200.,5700.]]
mwav_range = []
        
data = fits.open('Kate_Sources/spec-0389-51795-0332.fits')[1].data
flux =data.flux
wave = (10**data.loglam)/3.85
weight = (data.ivar*(data.and_mask == 0))

for rr in wav_range:
    if ((np.min(wave) < rr[0]) & (np.max(wave) > rr[1])):
        nrr0 = rr[0] ; nrr1 = rr[1]
    elif ((np.min(wave) > rr[0]) & (np.min(wave) < rr[1])& (np.max(wave) < rr[1])):
        nrr0 = np.min(wave) ; nrr1 = rr[1]
    elif ((np.min(wave) < rr[0]) & (np.max(wave) > rr[0]) &(np.max(wave) < rr[1])):
        nrr0 = rr[0] ; nrr1 = np.max(wave)
    else :
        continue
    mwav_range.append([nrr0,nrr1])
        
        #print mwav_range


AMP, ALPHA, CHISQ, DOF = compute_alpha(wave, flux, weight, mwav_range, per_value)
cwave = wave[contWave(wave)]
cflux = flux[contWave(wave)]
cweight = weight[contWave(wave)]
#sigma = cweight
#p0=[1,1]
#popt,pcov =  optimization.curve_fit(powerlawFunc,cwave,cflux,p0,sigma)
#print popt
plt.plot(cwave,cflux)
plt.plot(wave,AMP*wave**ALPHA,'--')
plt.show()
 
<file_sep>import numpy as np
import os


def velocity(wave,ion_wave,z):
    #print wave
    #print z,ion_wave
    c = 299792.0# in km/s
    vel =np.zeros(len(wave))
    zabs = np.zeros(len(wave))
    for i in range(len(wave)):
        #print i,wave[i],zabs[i],vel[i]
        zabs[i] = (wave[i]/ion_wave) -1.0
        vel[i] = -((((1+z)/(1+zabs[i]))**2-1.0)/(((1+z)/(1+zabs[i]))**2+1))*c
    return vel

cw=1549.0

pf = np.genfromtxt('BALcomponents_output.txt',names=('index','filename','vmin','vmax','nt','ntroughs','paddmin','paddmax','ew','ewerr'),dtype=(int,'|S35',float,float,float,float,float,float,float))
uindex = np.unique(pf['index'])
for i in  range(len(uindex)):
    match = np.where(pf['index'] == uindex[i])[0]
    mpf = pf[match]
    for k in range(len(mpf)):
        wave,flux,err = np.loadtxt('EmissionLines_Norm_Spectra/'+str(mpf['filename'][k])).T
        vel = velocity(wave, cw, 0.0)/1000.
        print 'Match', match, uindex[i]
        print mpf
        ntroughs = mpf['ntroughs'][0]
        print 'Ntroughs: ',ntroughs
        for ii in range(int(ntroughs)):
            print uindex[i]
            print 'index',ii,(mpf['vmax'][ii] - mpf['paddmin'][ii])/1000, (mpf['vmin'][ii] + mpf['paddmax'][ii])/1000.
            xx = np.where((vel >= (mpf['vmax'][ii] - mpf['paddmin'][ii])/1000.) & ( vel <= (mpf['vmin'][ii]+mpf['paddmax'][ii])/1000.))[0]
            nvel=vel[xx] ; nflux = flux[xx] ; nerr = err[xx]
            print 'Final',min(nvel),max(nvel)
            sfilename = 'spec-'+mpf['filename'][k].split('_')[1]+'_'+mpf['filename'][k].split('_')[2]+'_'+mpf['filename'][k].split('_')[3]
            savefilename = 'fivecol/'+sfilename+'.crop.5col'+'_'+str(ii)
            print savefilename
            np.savetxt(savefilename,zip(nvel,nflux,nerr,nflux,nerr), fmt='%10.5f')
            #fsdf=raw_input()
<file_sep>#!/usr/bin/python
#-*- coding: utf-8 -*-

import sys
import argparse
import os
import numpy as np
from astropy import io
from astropy.io import fits as pyfits
from astropy.io import ascii
from astropy.table import Table
from matplotlib import pyplot as plt
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy import optimize 
import matplotlib.cm as cm
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib import rcParams
#import crop_spectrum #not used in this version 
from Observation import Observation
import CCF_interp as myccf
from scipy import stats 

#Wavelengths of each species
nv_0 = 1238.821
siiv_0 = 1393.755
civ_0 = 1548.202
aliii_0 = 1854.7164
mgii_0 = 2796.352
lightspeed = 2.99792458e5
nu_0 = civ_0

#crop = True #(False if I've already done cropping and just want to redo the lags). This is from an old version; not used here. 

#Convert wavelength to velocity
def lam2vel(wavelength):
    zlambda = (wavelength-ion_0)/ion_0
    R = 1./(1+zlambda)
    vel_ion = -lightspeed*(R**2-1)/(R**2+1)
    return vel_ion

###################################################################################################
## Directories and list of files to use
###################################################################################################
tdssdatadir = 'EmissionLines_Norm_Spectra/'
bossdatadir = 'EmissionLines_Norm_Spectra/'
sdssdatadir = 'EmissionLines_Norm_Spectra/'

#bi_file = '/Users/cjg235/balqsos/varbal/measurements/BI_out_multitrough.dat'
bi_file = '/Users/vzm83/Acc_BALQSO/BALcomponents_output.txt' #'../../measurements/pyccf_tries_69kms_feb/trough_complexes.dat'
datafile = '/Users/vzm83/Acc_BALQSO/tdss_allmatches_crop_edit.dat'#'../../sample_selection/tdss_allmatches_crop_edit.dat' 

#Read in file containing all plate-mjd-fibers of pairs of spectra I want to use 
targ_obs = list()
with open(datafile) as file:
    for line in file:
        if line.startswith('#'):
            continue
        else:
            targ_obs.append(line.rstrip("\n"))

num_obj = np.size(targ_obs)
print 'Number of different targets: ', num_obj

target_observations = [[] for _ in range(num_obj)]
print target_observations
#exit()

#bi_info = np.genfromtxt(bi_file, dtype = None, names = ['files', 'num_troughs', 'BI', 'BI_err', 'BIround', 'Vmax', 'Verr', 'Vmaxround', 'Vmin', 'Verr2', 'Vminround', 'Chi2'], skiprows = 1)
#num_bi = np.size(bi_info['num_troughs'])
bi_info = np.genfromtxt(bi_file, dtype = None, names = ['speci', 'files', 'Vmin', 'Vmax', 'numt', 'num_troughs', 'highpad', 'lowpad', 'ew', 'ew_err'])
num_bi = np.size(bi_info['num_troughs'])
x = np.genfromtxt('linefits.txt', dtype = None, names = ['numobj', 'spec1', 'spec2', 'spec3', 'linefit'])
linefit = x['linefit']
num_obj=119
for i in xrange(0, num_obj):
    values = targ_obs[i].split()
    print values
    ra = values[0]
    dec = values[1]
    redshift = values[2]
    valuesn = sorted(values[3:])
    valuesize = np.size(valuesn)
    line_fit = linefit[i]
    if line_fit == 'gauss-hermite':
        file_suffix = '_cEm_gh.txt'
    elif line_fit == 'voigt':
        file_suffix = '_cEm_v.txt'
    else:
        file_suffix = '_cEm_2g.txt'
    
    #print valuesize, valuesn
    for j in xrange(0, valuesize):
        obs = Observation() 
        platemjdfiber = valuesn[j]
        plate = platemjdfiber.split('-')[0] 
        mjd = int(platemjdfiber.split('-')[1])
        fiber = platemjdfiber.split('-')[2]
        if mjd < 55050:
            #print 'sdss'
            filename = sdssdatadir+'spec-'+valuesn[j]+file_suffix
            survey = 'SDSS'
        if np.logical_and(mjd >= 55100., mjd <= 56850.):
            #print 'boss'
            filename = bossdatadir+'spec-'+valuesn[j]+file_suffix        
            survey = 'BOSS'
        if mjd > 56850.:
            #print 'tdss'
            filename = tdssdatadir+'spec-'+valuesn[j]+file_suffix   
            survey = 'TDSS'
        
        spec_filename = 'NormSpec_'+valuesn[j]+file_suffix
        match = np.where(bi_info['files'] == spec_filename)[0]
        print 'Match' ,match, spec_filename
        if np.size(match) == 0:
            obs.filename = filename
            obs.spec_name = spec_filename
            obs.mjd = mjd
            obs.pmf = platemjdfiber
            obs.survey = survey
            obs.redshift = redshift 
            obs.num_troughs = [0]
            obs.v_max = 0
            obs.v_min = 0
            print i, spec_filename
        else:
            obs.filename = filename
            obs.spec_name = spec_filename
            obs.mjd = mjd
            obs.pmf = platemjdfiber
            obs.survey = survey
            obs.redshift = redshift 
            obs.num_troughs = bi_info['num_troughs'][match]
            obs.v_max = bi_info['Vmax'][match]
            obs.v_min = bi_info['Vmin'][match]

        target_observations[i].append(obs)


print 'Carrying out CCF analysis' 

#Output files
lagsout = open('shift_measurements.dat', 'w')
lagsout.write('#Plot, Target-i, Spec1       Spec2       delta_t (years)  Centroid       UpErr     Lowerr       Peak       UpErr     Lowerr       R_val     Accel (cm/s^2)    UpErr     Lowerr       Shift      UpErr     Lowerr     vmin vmax  \n') 
lagsout.close()
lagsout3 = open('shift_measurements_3sigma.dat', 'w')
lagsout3.write('#Plot, Target-i, Spec1       Spec2       delta_t (years)  Centroid       UpErr     Lowerr       Peak       UpErr     Lowerr       R_val     Accel (cm/s^2)    UpErr     Lowerr       Shift      UpErr     Lowerr     vmin vmax  \n') 
lagsout3.close()

pdf_pages = PdfPages('all_ccfs.pdf')
numplot = 0
for i in xrange(0, num_obj):   #num_obj
    target = target_observations[i]
    numrep = np.size(target)
    n_troughs = target[0].num_troughs[0]
    if n_troughs > 0:
        print n_troughs
        for z in xrange(0, n_troughs):
            suffix = str(z)
            for j in xrange(0, numrep-1):
                k = 0
                while (j+k) < (numrep-1):
                    k+=1
                    print k
                    print 'Targetfilename',target[j].filename
                    spec1 = 'fivecol/'+(target[j].filename).split('/')[1]+'.crop.5col_'+suffix
                    spec2 = 'fivecol/'+(target[j+k].filename).split('/')[1]+'.crop.5col_'+suffix
                    print spec1,spec2
                    prefix1 = ((target[j].filename).split('/')[1]).split('-')[1]+'-'+ ((target[j].filename).split('/')[1]).split('-')[2]+'-'+ ((target[j].filename).split('/')[1]).split('-')[3]
                    prefix2 = ((target[j+k].filename).split('/')[1]).split('-')[1]+'-'+ ((target[j+k].filename).split('/')[1]).split('-')[2]+'-'+ ((target[j+k].filename).split('/')[1]).split('-')[3]
                    prefix = prefix1+'-vs-'+prefix2+'-'+suffix
                    print prefix
                    print 'i = ', i, ' Cross correlating: ', prefix
                    deltamjd = abs(int(target[j+k].mjd) - int(target[j].mjd))
                    deltamjd_rest = (deltamjd/(1+float(target[j].redshift)))/365.
            
                    #Calculate lag with python CCF program
                    nsim = 10000
                    wave1, flux1, err1 = np.loadtxt(spec1, unpack = True, usecols = [0, 1, 2])
                    wave2, flux2, err2 = np.loadtxt(spec2, unpack = True, usecols = [0, 1, 2])
                    print z
                    print target[j].v_min
                    vmin = target[j].v_min[z]
                    vmax = target[j].v_max[z] 
                
                    #Run the CCF on the two spectra
                    #Run the CCF from -2000 km/s to +2000 km/s and at 69 km/s interpolation
                    tlag_peak, status_peak, tlag_centroid, status_centroid, ccf_pack, max_r, status_r,peak_r = myccf.peakcent(wave1, flux1, wave2, flux2, -2.001, 2.001, 0.069, imode = 0, sigmode = 0.2)
                    tlags_peak, tlags_centroid, nsuccess_peak, nfail_peak, nsuccess_centroid, nfail_centroid, max_rvals, nfail_maxrvals,pvale = myccf.xcor_mc(wave1, flux1, abs(err1), wave2, flux2, abs(err2), -2.001, 2.001, 0.069, nsim = nsim, mcmode=2, sigmode = 0.2)
                    #plt.plot(ccf_pack[1], ccf_pack[0])
                    #plt.show()
                    #plt.hist(tlags_centroid, color = 'b')
                    #plt.hist(tlags_peak, color = 'r')
                    #plt.show()
                    
                    #Write out to file
                    centfile = open('ccf_files/centdist_'+prefix, 'w')
                    peakfile = open('ccf_files/peakdist_'+prefix, 'w')
                    ccf_file = open('ccf_files/ccf_'+prefix, 'w')
                    lag = ccf_pack[1]
                    r = ccf_pack[0]
                    for m in xrange(0, np.size(tlags_centroid)):
                        centfile.write('%5.5f \n'%(tlags_centroid[m]))
                    for m in xrange(0, np.size(tlags_peak)):
                        peakfile.write('%5.5f \n'%(tlags_peak[m]))
                    centfile.close()
                    peakfile.close()
                    for m in xrange(0, np.size(lag)):
                        ccf_file.write('%5.5f    %5.5f  \n'%(lag[m], r[m]))
                    ccf_file.close() 

                                #Calculate lag with python CCF program
                    centfile = 'ccf_files/centdist_'+prefix
                    peakfile = 'ccf_files/peakdist_'+prefix
                    ccf_file = 'ccf_files/ccf_'+prefix
                    lag, r = np.loadtxt(ccf_file, unpack = True)
                    tlags_centroid = np.loadtxt(centfile, unpack = True)
                    tlags_peak = np.loadtxt(peakfile, unpack = True) 

                    perclim = 84.1344746
                    perclim3sig = 99.8650102
                    
                    #For non-Gaussian Distribution
                    centau = stats.scoreatpercentile(tlags_centroid, 50)
                    centau_uperr = (stats.scoreatpercentile(tlags_centroid, perclim))-centau
                    centau_loerr = centau-(stats.scoreatpercentile(tlags_centroid, (100.-perclim)))
                    centau_uperr3 = (stats.scoreatpercentile(tlags_centroid, perclim3sig))-centau
                    centau_loerr3 = centau-(stats.scoreatpercentile(tlags_centroid, (100.-perclim3sig)))
                    print 'Centroid, error: ', centau, centau_loerr, centau_uperr
                    
                    peaktau = stats.scoreatpercentile(tlags_peak, 50)
                    peaktau_uperr = (stats.scoreatpercentile(tlags_peak, perclim))-centau
                    peaktau_loerr = centau-(stats.scoreatpercentile(tlags_peak, (100.-perclim)))
                    peaktau_uperr3 = (stats.scoreatpercentile(tlags_peak, perclim3sig))-centau
                    peaktau_loerr3 = centau-(stats.scoreatpercentile(tlags_peak, (100.-perclim3sig)))
                    print 'Peak, errors: ', peaktau, peaktau_uperr, peaktau_loerr

                    #Convert the measured shift in velocity space to acceleration.
                    rval = np.max(r) 
                    accel_meas = (centau*1000/deltamjd_rest)*(1.0e5)/(365.*24.*60.*60.)   #in cm/s/s
                    accel_uperr = (centau_uperr*1000./deltamjd_rest)*(1.0e5)/(365.*24.*60.*60.)
                    accel_loerr = (centau_loerr*1000./deltamjd_rest)*(1.0e5)/(365.*24.*60.*60.)
                    accel_uperr3 = (centau_uperr3*1000./deltamjd_rest)*(1.0e5)/(365.*24.*60.*60.)
                    accel_loerr3 = (centau_loerr3*1000./deltamjd_rest)*(1.0e5)/(365.*24.*60.*60.)

                    shift = centau
                    shift_uperr = centau_uperr
                    shift_loerr = centau_loerr
                    shift_uperr3 = centau_uperr3
                    shift_loerr3 = centau_loerr3

                    print 'Shift Error ratio: ', shift_uperr/shift_loerr
                    print 'Accel Error ratio: ', accel_uperr/accel_loerr
                    print 'Shift3 Error ratio: ', shift_uperr3/shift_loerr3
                    print 'Accel3 Error ratio: ', accel_uperr3/accel_loerr3
                    
                    #write all of this out to a file.
                    lagsout = open('shift_measurements.dat', 'a')
                    lagsout3 = open('shift_measurements_3sigma.dat', 'a')
                    lagsout.write('%i     %i     %s     %s     %5.5f     %5.5f     %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f  \n'%(numplot, \
                                i, spec1,spec2, deltamjd_rest, centau, centau_uperr, centau_loerr, peaktau, peaktau_uperr, peaktau_loerr, rval, accel_meas, accel_uperr, accel_loerr, shift, shift_uperr, shift_loerr, vmin/1000., vmax/1000.))
                    lagsout.close() 
                    lagsout3.write('%i     %i     %s     %s     %5.5f     %5.5f     %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f      %5.5f  \n'%(numplot, \
                                i, spec1,spec2, deltamjd_rest, centau, centau_uperr3, centau_loerr3, peaktau, peaktau_uperr3, peaktau_loerr3, rval, accel_meas, accel_uperr3, accel_loerr3, shift, shift_uperr3, shift_loerr3, vmin/1000., vmax/1000.))
                    lagsout3.close() 

                    #Make plots of the spectra, CCF, CCCD, and CCPD, etc 
                    spec1w, spec1f, spec1er = np.loadtxt(spec1, unpack = True, usecols = [0, 1, 2])
                    spec2w, spec2f, spec2er = np.loadtxt(spec2, unpack = True, usecols = [0, 1, 2])
                    spec1w = -spec1w ; spec2w = - spec2w
                    spec2wn = spec2w-shift

                    fig = plt.figure()
                    fig.subplots_adjust(hspace=0.2, wspace = 0.1)
            
                    ax1 = fig.add_subplot(2, 1, 1)
                    ax1.plot(spec1w, spec1f, color = 'k', label = 'First epoch')
                    ax1.plot(spec2w, spec2f, color = 'r', label = 'Second epoch')
                    try:
                        ax1.plot(spec2wn, spec2f, color = 'b', linestyle = '--', label = 'Shifted')
                    except:
                        print 'No shift measured!'
                    ax1.legend(loc = 'lower left', fontsize = 8)
                    ax1.set_title('Spectrum i = %i, Trough %s/%i, $\Delta$t$_{\\rm rest}$ = %5.3f years'%(i, suffix, n_troughs-1, deltamjd_rest), fontsize = 20)
                    ax1.text(0.05, 0.9, prefix, fontsize = 15, transform = ax1.transAxes)
                    ax1.set_xlim(30, -2)
                    ax1.set_ylabel('Normalized Flux')
                    ax1.set_xlabel('Velocity (10$^3$) km/s)')
            
                    xmin, xmax = -2,2
                    ax2 = fig.add_subplot(2, 3, 4)
                    ax2.set_ylabel('Correlation coefficient')
                    ax2.text(0.2, 0.85, 'CCF ', horizontalalignment = 'center', verticalalignment = 'center', transform = ax2.transAxes, fontsize = 16)
                    ax2.set_xlim(xmin, xmax)
                    ax2.set_ylim(0, 1.0)
                    try:
                        ax2.plot(lag, r, color = 'k')
                    except:
                        print 'No CCF data'
                
                    ax3 = fig.add_subplot(2, 3, 5, sharex = ax2)
                    #ax3.set_ylabel('N')
                    ax3.set_xlim(xmin, xmax)
                    ax3.axes.get_yaxis().set_ticks([])
                    ax3.set_xlabel('Shift: %5.1f + %5.1f - %5.1f km/s, Accel: %5.3f+ %5.3f - %5.3f cm s$^{-2}$'%(shift*1000, shift_uperr*1000, shift_loerr*1000, accel_meas, accel_uperr, accel_loerr), fontsize = 15) 
                    ax3.text(0.2, 0.85, 'CCCD ', horizontalalignment = 'center', verticalalignment = 'center', transform = ax3.transAxes, fontsize = 16)
                    try:
                        ax3.hist(tlags_centroid, color = 'k')
                    except:
                        print 'No centroid data'
                
                    ax4 = fig.add_subplot(2, 3, 6, sharex = ax2)
                    ax4.set_ylabel('N')
                    ax4.yaxis.tick_right()
                    ax4.yaxis.set_label_position('right') 
                    #ax4.set_xlabel('Lag (\AA)')
                    ax4.set_xlim(xmin, xmax)
                    ax4.text(0.2, 0.85, 'CCPD ', horizontalalignment = 'center', verticalalignment = 'center', transform = ax4.transAxes, fontsize = 16)
                    try:
                        ax4.hist(tlags_peak, color = 'k')
                    except:
                        print 'No peak data'
                
                    pdf_pages.savefig(fig)
                    plt.savefig('ccf_plots/ccfplot_'+prefix+'.png', format = 'png', bbox_inches = 'tight')
                    
                    plt.close(fig)
                    #ttt = raw_input('After plot')
            
pdf_pages.close()
<file_sep>import sys
import os
import numpy as np
from astropy import io
from astropy.io import fits 
from astropy.io import ascii
from matplotlib import pyplot as plt
from astropy import constants as const
from astropy import units as U
from astropy.coordinates import SkyCoord
from dustmaps.sfd import SFDQuery
from specutils import extinction  
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit 
from matplotlib.backends.backend_pdf import PdfPages
from lmfit import minimize, Parameters

def gaussfunc_gh(paramsin,x):
    print paramsin
    amp=paramsin['amp'].value
    center=paramsin['center'].value
    sig=paramsin['sig'].value
    c1=-np.sqrt(3); c2=-np.sqrt(6); c3=2/np.sqrt(3); c4=np.sqrt(6)/3; c5=np.sqrt(6)/4
    skew=paramsin['skew'].value
    kurt=paramsin['kurt'].value
    scale = paramsin['scale'].value ; alpha = paramsin['alpha'].value 
 
    gaustot_gh=amp*np.exp(-.5*((x-center)/sig)**2)*(1+skew*(c1*((x-center)/sig)+c3*((x-center)/sig)**3)+kurt*(c5+c2*((x-center)/sig)**2+c4*((x-center)/sig)**4))+scale*x**alpha
    return gaustot_gh


def gaussfunc_2g(paramsin,x):
    amp1=paramsin['amp1'].value; amp2=paramsin['amp2'].value;
    center1=paramsin['center1'].value; center2=paramsin['center2'].value;
    sig1=paramsin['sig1'].value; sig2=paramsin['sig2'].value;
    scale = paramsin['scale'].value ; alpha = paramsin['alpha'].value 
    gaus1=amp1*np.exp(-.5*((x-center1)/sig1)**2)
    gaus2=amp2*np.exp(-.5*((x-center2)/sig2)**2)
    gaustot_2g=(gaus1+gaus2+scale*x**alpha)
    return gaustot_2g


def PLot(vels,stackspec , fit_gh, fit_2g, pars_gh, pars_2g, resid_gh, resid_2g):
    fig3=plt.figure(3)
    f1=fig3.add_axes((.1,.3,.8,.6))
        #xstart, ystart, xwidth, yheight --> units are fraction of the image from bottom left
    
    plt.plot(vels,stackspec,'k.')
    pgh,=plt.plot(vels,fit_gh,'b')
    p2g,=plt.plot(vels,fit_2g,'r')
    f1.set_xticklabels([]) #We will plot the residuals below, so no x-ticks on this plot
    plt.title('Multiple Gaussian Fit Example')
    plt.ylabel('Amplitude (Some Units)')
    f1.legend([pgh,p2g],['Gaus-Hermite','2-Gaus'],prop={'size':10},loc='center left')
     
    from matplotlib.ticker import MaxNLocator
    plt.gca().yaxis.set_major_locator(MaxNLocator(prune='lower')) #Removes lowest ytick label
    
    f1.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f' \
        %(pars_gh[0],pars_gh[1],pars_gh[2],pars_gh[3],pars_gh[4]),xy=(.05,.95), \
        xycoords='axes fraction',ha="left", va="top", \
    bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    f1.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f' \
        %(pars_2g[0],pars_2g[3],pars_2g[1],pars_2g[4],pars_2g[2],pars_2g[5]),xy=(.95,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    
    f2=fig3.add_axes((.1,.1,.8,.2))
     
    resgh,res2g,=plt.plot(vels,resid_gh,'k--',vels,resid_2g,'k')
     
    plt.ylabel('Residuals')
    plt.xlabel('Velocity (km s$^{-1}$)')
    f2.legend([resgh,res2g],['Gaus-Hermite','2-Gaus'],numpoints=4,prop={'size':9},loc='upper left')
     
    plt.savefig(pp,format='pdf')
    plt.show()
     
    plt.clf()


df = np.genfromtxt('tdss_allmatches_crop_edit.dat',names=['ra','dec','z','pmf1','pmf2','pmf3'],dtype=(float,float,float,'|S15','|S15','|S15'))
wav_range= [(1280,1350),(1700,1800),(1950,2200),(2650,2710),(3010,3700),(3950,4050),(4140,4270)]
pp = PdfPages('ContinuumNormalization_plus_EmissionLineFits.pdf')
#fx=open('ALPHA_AMP_values_V.txt','w')
#fx1=open('Kate_ALPHA_AMP_values_V.txt','w')
#for i in range(len(df['pmf1'])):
#for i in range(len(df)):
for i in range(4):

    print 'Norm_Spectra/Normspec_'+df['pmf1'][i]+'.txt' 
    print 'Norm_Spectra/Normspec_'+df['pmf2'][i]+'.txt' 
    print 'Norm_Spectra/Normspec_'+df['pmf3'][i]+'.txt' 
    if ((not os.path.isfile('Norm_Spectra/Normspec_'+df['pmf1'][i]+'.txt')) | (not os.path.isfile('Norm_Spectra/Normspec_'+df['pmf2'][i]+'.txt')) | ( not os.path.isfile('Norm_Spectra/Normspec_'+df['pmf3'][i]+'.txt'))):
        continue
    else:
        data1 = np.loadtxt('Norm_Spectra/Normspec_'+df['pmf1'][i]+'.txt')
        data2 = np.loadtxt('Norm_Spectra/Normspec_'+df['pmf2'][i]+'.txt')
        data3 = np.loadtxt('Norm_Spectra/Normspec_'+df['pmf3'][i]+'.txt')
    wave1 = data1.T[0] ;    flux1 = data1.T[1] ; weight1 =  data1.T[2]
    wave2 = data2.T[0] ;    flux2 = data2.T[1] ; weight2 =  data2.T[2]
    wave3 = data3.T[0] ;    flux3 = data3.T[1] ; weight3 =  data3.T[2]
    print len(wave1),len(flux1),len(weight1)
    
    print weight1
    
    p_gh=Parameters()
    p_gh.add('amp',value=np.max(flux1),vary=True);
    p_gh.add('center',value=1550,min=1540,max=1560,vary=True);
    p_gh.add('sig',value=15,min=3,max=10,vary=True);
    p_gh.add('skew',value=0,vary=True,min=None,max=None);
    p_gh.add('kurt',value=0,vary=True,min=None,max=None);
    p_gh.add('scale',value=1,vary=True,min=None,max=None);
    p_gh.add('alpha',value=-1.5,vary=True,min=-3,max=3);
 
    p_2g=Parameters()
    p_2g.add('amp1',value=np.max(flux1)/2.,min=.1*np.max(flux1),max=np.max(flux1),vary=True);
    p_2g.add('center1',value=1550,min=1540,max=1560,vary=True);
    p_2g.add('sig1',value=15,min=3,max=10,vary=True);
    p_2g.add('amp2',value=np.max(flux1)/2.,min=.1*np.max(flux1),max=np.max(flux1),vary=True);
    p_2g.add('center2',value=1545,min=1540,max=1560,vary=True);
    p_2g.add('sig2',value=15,min=3,max=10,vary=True);
    p_2g.add('scale',value=1,vary=True,min=None,max=None);
    p_2g.add('alpha',value=-1.5,vary=True,min=-3,max=3);

    gausserr_gh = lambda p,x,y: gaussfunc_gh(p,x)-y
    gausserr_2g = lambda p,x,y: gaussfunc_2g(p,x)-y

    fitout_gh=minimize(gausserr_gh,p_gh,args=(wave1,flux1))
    fitout_2g=minimize(gausserr_2g,p_2g,args=(wave1,flux1))
    
    pars_gh=[p_gh['amp'].value,p_gh['center'].value,p_gh['sig'].value,p_gh['skew'].value,p_gh['kurt'].value]
    pars_2g=[p_2g['amp1'].value,p_2g['center1'].value,p_2g['sig1'].value,p_2g['amp2'].value,p_2g['center2'].value,p_2g['sig2'].value]

    fit_gh=gaussfunc_gh(p_gh,wave1)
    fit_2g=gaussfunc_2g(p_2g,wave1)
    resid_gh=fit_gh- flux1
    resid_2g=fit_2g- flux1

    print('Fitted Parameters (Gaus+Hermite):\nAmp = %.2f , Center = %.2f , Disp = %.2f\nSkew = %.2f , Kurt = %.2f' \
%(pars_gh[0],pars_gh[1],pars_gh[2],pars_gh[3],pars_gh[4]))
 
    print('Fitted Parameters (Double Gaussian):\nAmp1 = %.2f , Center1 = %.2f , Sig1 = %.2f\nAmp2 = %.2f , Center2 = %.2f , Sig2 = %.2f' \
%(pars_2g[0],pars_2g[1],pars_2g[2],pars_2g[3],pars_2g[4],pars_2g[5]))
    PLot(wave1,flux1,fit_gh, fit_2g, pars_gh, pars_2g, resid_gh, resid_2g)
<file_sep>import numpy as np
import os


def velocity(wave,ion_wave,z):
    #print wave
    #print z,ion_wave
    c = 299792.0# in km/s
    vel =np.zeros(len(wave))
    zabs = np.zeros(len(wave))
    for i in range(len(wave)):
        #print i,wave[i],zabs[i],vel[i]
        zabs[i] = (wave[i]/ion_wave) -1.0
        vel[i] = -((((1+z)/(1+zabs[i]))**2-1.0)/(((1+z)/(1+zabs[i]))**2+1))*c
    return vel

cw=1549.0

pf = np.genfromtxt('BALcomponents_output.txt',names=('index','filename','vmin','vmax','nt','ntroughs','paddmin','paddmax','ew','ewerr'),dtype=(int,'|S35',float,float,float,float,float,float,float))


for i in  range(len(pf)):
    wave,flux,err = np.loadtxt(pf['filename'][i]).T
    vel = velocity(wave1, cw, 0.0)
    xx = np.where((vel >= pf['vmin'][i] - pf['paddmin'][i]) & ( vel <= pf['vmax'][i]+pf['paddmax'][i]))[]
    nvel=vel[xx]/1000 ; nflux = flux[xx] ; nerr = 1.0/np.sqrt(err[xx])
    sfilename = 'spec-'+pf['filename'][i].split(_)[1]+pf['filename'][i].split(_)[2]pf['filename'][i].split(_)[3]
    savefilename = 'fivecol/'+sfilename+'.crop.5col_0'
    print savefilename
    #np.savetxt(savefilename,zip(nvel,nflux,nerr,nflux,nerr))
<file_sep>import numpy as np
import os 


df = np.genfromtxt('tdss_allmatches_crop_edit.dat',names=['ra','dec','z','pmf1','pmf2','pmf3'],dtype=(float,float,float,'|S15','|S15','|S15'))

pf = np.genfromtxt('BALcomponents_output.txt',names=('index','filename','vmin','vmax','nt','ntroughs','paddmin','paddmax','ew','ewerr'),dtype=(int,'|S35',float,float,float,float,float,float,float))
print pf['filename']
fpmf = []
fittype =[]
for k in range(len(pf['filename'])):
    fpmf.append(pf['filename'][k].strip().split('_')[1])
    if (pf['filename'][k].strip().split('_')[3].split('.')[0] == 'v'):
        fittype.append('voigt')
    else:
        fittype.append('2guass')
ab=open('linefits.txt','w')
fpmf = np.array(fpmf)
#print fpmf
#print fittype[0]
for i in range(len(df)):
    #print type(df['pmf1'][i]),type(fpmf)
    xx = np.where(fpmf == str(df['pmf1'][i]).strip())[0]
    print xx
    if len(xx)> 0:
        print>>ab,'{}\t{}\t{}\t{}\t{}'.format( i,df['pmf1'][i],df['pmf2'][i],df['pmf3'][i],fittype[xx[0]])

    #dksh=raw_input()
ab.close()
<file_sep>#!/usr/bin/python
#-*- coding: utf-8 -*-

import scipy
import numpy as np

class Spectrum(object):
    '''
    The Spectrum object
    '''
    def __init__(self):
        '''
        Initialize spectrum object
        '''
        self._wavelengths= None
        self._flux = None
        self._flux_error = None
        self._c = None
        self._units = None
        
    @property
    def wavelengths(self):
        return self._wavelengths
        
    @wavelengths.setter
    def wavelengths(self, new_wave):
        self._wavelengths = new_wave
        
    @property
    def flux(self):
        return self._flux
        
    @flux.setter
    def flux(self, new_flux):
        self._flux = new_flux

    @property
    def flux_error(self):
        return self._flux_error
        
    @flux_error.setter
    def flux_error(self, new_ferr):
        self._flux_error = new_ferr

    @property
    def c(self):
        return self._c
        
    @c.setter
    def c(self, new_c):
        self._c = new_c

    @property
    def units(self):
        return self._units
        
    @units.setter
    def units(self, new_units):
        self._units = new_units

    @property
    def redshift(self):
        return self._redshift
        
    @units.setter
    def redshift(self, redshift):
        self._redshift = new_redshift

    @property
    def objname(self):
        return self._objname
        
    @units.setter
    def objname(self, objname):
        self._objname = new_objname










<file_sep>#!/usr/bin/python
#-*- coding: utf-8 -*-

import sys
import numpy as np
from astropy import io
from astropy.io import fits as pyfits
from astropy.io import ascii
from astropy.table import Table
from matplotlib import pyplot as plt
from Spectrum import Spectrum 
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
from pylab import *
from scipy import optimize 

def runningMeanFast(x,N):
    return np.convolve(x, np.ones((N,))/N)[(N-1):]

def SmoothBoxCar(x,N):
    boxcar=np.ones(N)
    return convolve(x, boxcar/boxcar.sum())

#Define the continuum power law function
#reddening is E(B - V)

def redpowerlaw(x, amp, index, reddening):
  w = x/10000  #convert wavelengths to microns
  a = [185, 27, 0.005, 0.010, 0.012, 0.030]
  l = [0.042, 0.08, 0.22, 9.7, 18, 25]
  b = [90, 5.50, -1.95, -1.95, -1.80, 0]
  n = [2, 4, 2, 2, 2, 2]
  K = [2.89, 0.91, 0.02, 1.55, 1.72, 1.89]
  R_V = 2.93
  
  xi = 0
  for i in xrange(0,6):
      xi+= a[i] / (np.power((w/l[i]), n[i]) + np.power((l[i]/w), n[i]) + b[i])

  ElB = (xi*(1+R_V)-R_V)*reddening  ##E(lambda - V)
  extinction = np.power(10, -0.4*ElB)    ##10^(-0.4*E(lambda - V))

  return (amp*(x/2000.)**(index))*np.power(10, (-0.4*(xi*(1+R_V)-R_V)*reddening))

####Write a fitting program ################
def fit_three_params(x, y, errs):

    fittingx = x
    fittingy = y
    weights = errs
    #print x[0:10], y[0:10], weights[0:10]
    p_init = [8.,  -1.0,  0.05]
       
    out, var_matrix = optimize.curve_fit(redpowerlaw, fittingx, fittingy, p_init, maxfev = 10000)
    #print out
    #print var_matrix
    pfin = out
    #cov = out[1]
    redco = pfin[2] 
    index = pfin[1]
    amp = pfin[0]
    #print amp, index, redco
    #errors = np.sqrt(np.diag(cov))
    #indexErr = errors[0]
    #ampErr = errors[1]
    #redcoErr = errors[2]
    modelval = redpowerlaw(fittingx, amp, index, redco)
   
    return amp, index, redco

def continuum_fit(wave, flux, error, redshift): ###########################
    spectrum = Spectrum()
    spectrum.wavelengths = wave
    spectrum.flux = flux
    spectrum.flux_error = error
    
    smooth=False
    N_kern=5 #pxl
    plotfit=True 
    fitcon=True 

    ############################################
    ###########################################
    ####### Define line-free regions and calculate weighting scheme. This will need to 
    ####### be updated to deal with different wavelength regions. 
    maxwav = np.max(spectrum.wavelengths)
    minwav = np.min(spectrum.wavelengths)

    #RLF regions from Filiz ak et al., updated a bit, and also using a different one for low-z targets. 
    if np.logical_and(redshift > 1.65, redshift < 1.85):
        index1=np.where(np.logical_and(spectrum.wavelengths>=1280.,spectrum.wavelengths<=1350.))[0]
        index2=np.where(np.logical_and(spectrum.wavelengths>=1425.,spectrum.wavelengths<=1450.))[0]
        index3=np.where(np.logical_and(spectrum.wavelengths>=1700.,spectrum.wavelengths<=1800.))[0]
        index4=np.where(np.logical_and(spectrum.wavelengths>=1950.,spectrum.wavelengths<=2200.))[0]
        index5=np.where(np.logical_and(spectrum.wavelengths>=2650.,spectrum.wavelengths<=2710.))[0]
        index6=np.where(np.logical_and(spectrum.wavelengths>=3010.,spectrum.wavelengths<=3700.))[0]
        index7=np.where(np.logical_and(spectrum.wavelengths>=3950.,spectrum.wavelengths<=4050.))[0]
        index8=np.where(np.logical_and(spectrum.wavelengths>=4140.,spectrum.wavelengths<=4270.))[0]
        index9=np.where(np.logical_and(spectrum.wavelengths>=4400.,spectrum.wavelengths<=4770.))[0]
        index10=np.where(np.logical_and(spectrum.wavelengths>=5100.,spectrum.wavelengths<=6400.))[0]
    if np.logical_or(redshift >= 1.85, redshift <= 1.65):
        index1=np.where(np.logical_and(spectrum.wavelengths>=1280.,spectrum.wavelengths<=1350.))[0]
        index2=np.where(np.logical_and(spectrum.wavelengths>=1700.,spectrum.wavelengths<=1800.))[0]
        index3=np.where(np.logical_and(spectrum.wavelengths>=1950.,spectrum.wavelengths<=2200.))[0]
        index4=np.where(np.logical_and(spectrum.wavelengths>=2650.,spectrum.wavelengths<=2710.))[0]
        index5=np.where(np.logical_and(spectrum.wavelengths>=3010.,spectrum.wavelengths<=3700.))[0]
        index6=np.where(np.logical_and(spectrum.wavelengths>=3950.,spectrum.wavelengths<=4050.))[0]
        index7=np.where(np.logical_and(spectrum.wavelengths>=4140.,spectrum.wavelengths<=4270.))[0]
        index8=np.where(np.logical_and(spectrum.wavelengths>=4400.,spectrum.wavelengths<=4770.))[0]
        index9=np.where(np.logical_and(spectrum.wavelengths>=5100.,spectrum.wavelengths<=6400.))[0]
        index10=np.where(spectrum.wavelengths>=6900.)[0]
        #print 'made it to here!'
    indexes=np.concatenate((index1, index2, index3, index4, index5, index6, index7, index8, index9, index10))
    #Calculate weights for each pixel 
    n1, n2, n3, n4, n5, n6, n7, n8, n9, n10 = np.size(index1), np.size(index2), np.size(index3), \
        np.size(index4), np.size(index5), np.size(index6), np.size(index7), np.size(index8), np.size(index9), np.size(index10) 
    ntot = n1+n2+n3+n4+n5+n6+n7+n8+n9+n10
    w1 = list() 
    w2 = list() 
    w3 = list() 
    w4 = list() 
    w5 = list() 
    w6 = list()
    w7 = list()
    w8 = list()
    w9 = list()
    w10 = list() 
    for j in xrange(0, n1):
        w1.append(ntot/float(n1))
    for j in xrange(0, n2):
        w2.append(ntot/float(n2))
    for j in xrange(0, n3):
        w3.append(ntot/float(n3))
    for j in xrange(0, n4):
         w4.append(ntot/float(n4))
    for j in xrange(0, n5):
        w5.append(ntot/float(n5))
    for j in xrange(0, n6):
        w6.append(ntot/float(n6))
    for j in xrange(0, n7):
        w7.append(ntot/float(n7))
    for j in xrange(0, n8):
        w8.append(ntot/float(n8))
    for j in xrange(0, n9):
        w9.append(ntot/float(n9))
    for j in xrange(0, n10):
        w10.append(ntot/float(n10))

    allweights = concatenate((w1, w2, w3, w4, w5, w6, w7, w8, w9, w10))
    ############################################
    ###########################################
    ##Fit the spectrum to get the actual fit first. 
    spec_to_fit1 = Spectrum()
    spec_to_fit1.wavelengths = spectrum.wavelengths[indexes]
    spec_to_fit1.flux =spectrum.flux[indexes]
    spec_to_fit1.flux_error = spectrum.flux_error[indexes]
    num_wavs_fit = np.size(spec_to_fit1.wavelengths)

    toweight = True
    if toweight: 
        weights_specfit = allweights
    else:
        weights_specfit = spec_to_fit1.flux_error
        #weights_specfit = ones(num_wavs_fit) 
    
    #Iteratively fit continuum and drop 3sigma outliers with 10 iterations. 
    newregions = spec_to_fit1.wavelengths
    newflux = spec_to_fit1.flux
    newweights = weights_specfit
    newfluxerr = spec_to_fit1.flux_error
    checksize = 10 
    while checksize > 0:
        amp, index, redco = fit_three_params(newregions, newflux, newweights)
        first_continuum_fit = redpowerlaw(spectrum.wavelengths, amp, index, redco)

    ###iterate and remove pixels that deviate from too 
        con_fit_fitreg = redpowerlaw(newregions, amp, index, redco)
        diff = con_fit_fitreg - newflux
        sigma = np.std(diff)
        toolarge = np.where(diff >= 3*sigma)
        justright = np.where(diff <= 3*sigma)
        #print np.size(newregions)
        #print np.size(toolarge)
        checksize = np.size(toolarge) 

        newregions = newregions[justright]
        newflux = newflux[justright] 
        newweights = newweights[justright]
        newfluxerr = newfluxerr[justright] 
   
    finalspecfitwave, finalspecfitflux, finalweights_specfit = newregions, newflux, newweights

    #final continuum fit
    amp, index, redco = fit_three_params(finalspecfitwave, finalspecfitflux, finalweights_specfit)
    continuum_fit = redpowerlaw(spectrum.wavelengths, amp, index, redco)

    #print 'Amplitude, Index, Reddening', amp, index, redco 
    spec_to_fit = Spectrum()
    spec_to_fit.wavelengths = finalspecfitwave
    spec_to_fit.flux = finalspecfitflux
    spec_to_fit.flux_error = newfluxerr 
    num_wavs_fit = np.size(spec_to_fit.wavelengths)
            
    toweight = True
    if toweight: 
        weights_specfit = allweights
    else:
        weights_specfit = spec_to_fit.flux_error
        #weights_specfit = ones(num_wavs_fit) 
    
    ############################################
    ###########################################
    ##### Now do Monte Carlo simulations to get the uncertainties. 
    num_samples = 100
    num_wavelengths_fit = np.size(spec_to_fit.wavelengths)
    num_wavelengths_tot = np.size(spectrum.wavelengths)
            
    if toweight:
        weights = allweights
    else:
        #weights = ones(num_wavelengths_fit)
        weights = ones(num_wavelengths_fit)
            
    consamples = zeros((num_samples, num_wavelengths_tot))
    ampsamples = zeros(num_samples)
    indexsamples = zeros(num_samples)
        
    meancon = zeros(num_wavelengths_tot) 
    errcon = zeros(num_wavelengths_tot)
      
    #Monte carlo iterations for uncertainties
    for j in xrange(0, num_samples):
        specfit = zeros(num_wavelengths_fit)
        for k in xrange(0, num_wavelengths_fit):
            specfit[k] = spec_to_fit.flux[k]+randn()*spec_to_fit.flux_error[k]

        plot2 = False
        if plot2:
            figure2, ax3 = plt.subplots(1, figsize=(12,4))
            ax3.errorbar(spectrum.wavelengths, spectrum.flux, yerr = spectrum.flux_error, color = 'b', label = 'Spec')
            ax3.plot(spectrum.wavelengths, spectrum.flux, color = 'k', label = 'Spec')
            ax3.plot(spec_to_fit.wavelengths, specfit, color = 'r', label = 'Fitting flux')
            ymin = min(0, np.percentile(spectrum.flux,  1))
            ymax = 1.25*np.percentile(spectrum.flux, 99)
            ax3.set_xlim(min(spectrum.wavelengths), max(spectrum.wavelengths))
            ax3.set_ylim(ymin, ymax)
            ax3.set_xlabel('Rest-frame Wavelength (\AA)', fontsize = 20)
            ax3.set_ylabel('Flux $(10^{-17} erg/cm^2/s/\AA)$', fontsize = 20)
            #ax2.set_title('SDSS '+str(object_info['name'][i])+ ', $z$ = '+str(object_info['z'][i]), fontsize = 22)
            ax3.tick_params(axis='x', labelsize=19)
            ax3.tick_params(axis='y', labelsize=19)
            ax3.legend()
            ax3.grid()
            plt.show()
            plt.close(figure2)

        ampout, indexout, redout = fit_three_params(spec_to_fit.wavelengths, specfit, weights)
        contflux = redpowerlaw(spectrum.wavelengths, ampout, indexout, redout)
        ampsamples[j] = ampout
        indexsamples[j] = indexout 
        for m in xrange(0, num_wavelengths_tot):
            consamples[j,m] += contflux[m]
    for j in xrange(0, num_wavelengths_tot):
        meancon[j] = mean(consamples[:,j])
        errcon[j] = std(consamples[:,j])
        
    conflux = continuum_fit
    confluxerr = errcon 
    mean_amp = amp
    err_amp = np.std(ampsamples)
    mean_index = index
    err_index = np.std(indexsamples)

    plotfit = False
    if plotfit:
        figure, ax2 = plt.subplots(1, figsize=(12,4))
        ax2.plot(spectrum.wavelengths, spectrum.flux, color = 'k', label = 'Spec')
        ax2.plot(spectrum.wavelengths, conflux, color = 'b', label = 'Continuum', linewidth = 2)
        ymin = min(0, np.percentile(spectrum.flux,  1))
        ymax = 1.25*np.percentile(spectrum.flux, 99)
        ax2.set_xlim(min(spectrum.wavelengths), max(spectrum.wavelengths))
        ax2.set_ylim(ymin, ymax)
        ax2.set_xlabel('Rest-frame Wavelength (\AA)', fontsize = 20)
        ax2.set_ylabel('Flux $(10^{-17} erg/cm^2/s/\AA)$', fontsize = 20)
        #ax2.set_title('SDSS '+str(object_info['name'][i])+ ', $z$ = '+str(object_info['z'][i]), fontsize = 22)
        ax2.tick_params(axis='x', labelsize=19)
        ax2.tick_params(axis='y', labelsize=19)
        ax2.legend()
        ax2.grid()
        plt.show()

    #Return the continuum fit, flux, uncertainty associated with the Monte Carlo iterations, and fitting regions. 
    return spectrum.wavelengths, conflux, confluxerr, mean_amp, mean_index, redco, spec_to_fit.wavelengths, spec_to_fit.flux
    
<file_sep>import numpy as np
import os

df = np.genfromtxt('tdss_allmatches_crop_edit.dat',names=['ra','dec','z','pmf1','pmf2','pmf3'],dtype=(float,float,float,'|S15','|S15','|S15'))



def download_spectra(plate, mjd, fiber, dirname='.'):
    '''  Downloads SDSS spectra from DR14 and puts it in dirname
         Change the SDSS URL to download from a different location
    '''
    FITS_FILENAME = 'spec-%(plate)04i-%(mjd)05i-%(fiber)04i.fits'
    SDSS_URL = ('https://data.sdss.org/sas/dr14/eboss/spectro/redux/v5_10_0/spectra/%(plate)04i/'
            'spec-%(plate)04i-%(mjd)05i-%(fiber)04i.fits')
    SDSS_URL = ('https://data.sdss.org/sas/dr8/sdss/spectro/redux/26/spectra/%(plate)04i/'
            'spec-%(plate)04i-%(mjd)05i-%(fiber)04i.fits')
    # print SDSS_URL % dict(plate=plate,mjd=mjd,fiber=fiber)
    download_url = 'wget  '+SDSS_URL % dict(plate=plate,mjd=mjd,fiber=fiber)
    print download_url
    os.system(download_url)
    mv_cmd='mv '+FITS_FILENAME % dict(plate=plate,mjd=mjd,fiber=fiber) + ' '+dirname+'/.'
    #print mv_cmd
    os.system(mv_cmd)



for i in range(len(df['pmf1'])):
    print df['pmf1'][i]
    plate = int(df['pmf1'][i].split('-')[0])
    mjd = int(df['pmf1'][i].split('-')[1])
    fiber = int(df['pmf1'][i].split('-')[2])
    print plate,mjd,fiber
    download_spectra(int(df['pmf1'][i].split('-')[0]),int(df['pmf1'][i].split('-')[1]),int(df['pmf1'][i].split('-')[2]),'Kate_Sources')
    download_spectra(int(df['pmf2'][i].split('-')[0]),int(df['pmf2'][i].split('-')[1]),int(df['pmf2'][i].split('-')[2]),'Kate_Sources')
    download_spectra(int(df['pmf3'][i].split('-')[0]),int(df['pmf3'][i].split('-')[1]),int(df['pmf3'][i].split('-')[2]),'Kate_Sources')

<file_sep>import sys
import os
import numpy as np
from astropy import io
from astropy.io import fits 
from astropy.io import ascii
from matplotlib import pyplot as plt
from astropy import constants as const
from astropy import units as U
from astropy.coordinates import SkyCoord
from dustmaps.sfd import SFDQuery
from specutils import extinction  
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit 
from matplotlib.backends.backend_pdf import PdfPages
from lmfit import minimize, Parameters
from astropy.modeling.models import Voigt1D

def dered_flux(Av,wave,flux):
    dered_flux = np.zeros(len(wave))
    for i in range(len(wave)):
        ext_lambda= extinction.extinction_ccm89(wave[i] * U.angstrom, a_v= Av, r_v= 3.1)
        tau_lambda= ext_lambda/(1*1.086) 
        dered_flux[i]= flux[i] * np.exp(tau_lambda)
    return dered_flux

def myVoigt(x, amp, center, fwhm_l, fwhm_g, scale, alpha):
    v1= Voigt1D(x_0=center, amplitude_L=amp, fwhm_L=fwhm_l, fwhm_G=fwhm_g)
    powerlaw = scale*x**alpha
    voigt = v1(x)
    voigt_tot = (voigt+powerlaw)
    return voigt_tot

def myGaussHermite(x, amp, center, sig, skew, kurt, scale, alpha):
    c1=-np.sqrt(3); c2=-np.sqrt(6); c3=2/np.sqrt(3); c4=np.sqrt(6)/3; c5=np.sqrt(6)/4
    gausshermite = amp*np.exp(-.5*((x-center)/sig)**2)*(1+skew*(c1*((x-center)/sig)+c3*((x-center)/sig)**3)+kurt*(c5+c2*((x-center)/sig)**2+c4*((x-center)/sig)**4))
    powerlaw = scale*x**alpha
    gaustot_gh = (gausshermite+powerlaw)
    return gaustot_gh


def myDoubleGauss(x, amp1, center1, sig1, amp2, center2, sig2, scale, alpha):
    gaus1=amp1*np.exp(-.5*((x-center1)/sig1)**2)
    gaus2=amp2*np.exp(-.5*((x-center2)/sig2)**2)
    powerlaw = scale*x**alpha
    gaustot_2g= (gaus1+gaus2+powerlaw)
    
    return gaustot_2g

def ext_coeff(lamb):
    inv_lamba=[0.45,0.61,0.8,1.82,2.27,2.7,3.22,3.34,3.46,3.6,3.75,3.92,4.09,4.28,4.50,4.73,5.00,5.24,5.38,5.52,5.70,5.88,6.07,6.27,6.48,6.72,6.98,7.23,7.52,7.84]
    smc_ext=[-2.61,-2.47,-2.12,0.0,1.0,1.67,2.29,2.65,3.0,3.15,3.49,3.91,4.24,4.53,5.3,5.85,6.38,6.76,6.9,7.17,7.71,8.01,8.49,9.06,9.28,9.84,10.8,11.51,12.52,13.54]
    xy=np.interp(1.0/(lamb*10**(-4)),inv_lamba,smc_ext)
    ext_lamb=(xy+2.98)/3.98 # Rv=2.98
    #return(z,ext_lamb)
    return(ext_lamb)

def myfunct(wave, a , b):
    return a*ext_coeff(wave)*wave**(b)


def waveRange(wave):
    #wav_range=[[1300.,1350.],[1420.,1470.],[1700.,1800.],[2080.,2215],[2480.,2655],[3225.,3900.],[4200.,4230.],[4435.,4700.],[5200.,5700.]]
    wav_range= [(1400,1850)]


    mwav_range = []
    for rr in wav_range:
        if ((np.min(wave) < rr[0]) & (np.max(wave) > rr[1])):
            nrr0 = rr[0] ; nrr1 = rr[1]
        elif ((np.min(wave) > rr[0]) & (np.min(wave) < rr[1])& (np.max(wave) < rr[1])):
            nrr0 = np.min(wave) ; nrr1 = rr[1]
        elif ((np.min(wave) < rr[0]) & (np.max(wave) > rr[0]) &(np.max(wave) < rr[1])):
            nrr0 = rr[0] ; nrr1 = np.max(wave)
        else :
            continue
        mwav_range.append([nrr0,nrr1])
    return mwav_range

def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
    The Savitzky-Golay filter removes high frequency noise from data.
    It has the advantage of preserving the original shape and
    features of the signal better than other types of filtering
    approaches, such as moving averages techniques.
    Parameters
    ----------
    y : array_like, shape (N,)
        the values of the time history of the signal.
    window_size : int
        the length of the window. Must be an odd integer number.
    order : int
        the order of the polynomial used in the filtering.
        Must be less then `window_size` - 1.
    deriv: int
        the order of the derivative to compute (default = 0 means only smoothing)
    Returns
    -------
    ys : ndarray, shape (N)
        the smoothed signal (or it's n-th derivative).
    Notes
    -----
    The Savitzky-Golay is a type of low-pass filter, particularly
    suited for smoothing noisy data. The main idea behind this
    approach is to make for each point a least-square fit with a
    polynomial of high order over a odd-sized window centered at
    the point.
    Examples
    --------
    t = np.linspace(-4, 4, 500)
    y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
    ysg = savitzky_golay(y, window_size=31, order=4)
    import matplotlib.pyplot as plt
    plt.plot(t, y, label='Noisy signal')
    plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
    plt.plot(t, ysg, 'r', label='Filtered signal')
    plt.legend()
    plt.show()
    References
    ----------
    .. [1] <NAME>, <NAME>, Smoothing and Differentiation of
       Data by Simplified Least Squares Procedures. Analytical
       Chemistry, 1964, 36 (8), pp 1627-1639.
    .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing
       W.H. Press, <NAME>, <NAME>, <NAME>
       Cambridge University Press ISBN-13: 9780521880688
    """
    import numpy as np
    from math import factorial
    
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except (ValueError, msg):
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    # precompute coefficients
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    # pad the signal at the extremes with
    # values taken from the signal itself
    firstvals = np.asarray(y[0]) - np.abs( y[1:half_window+1][::-1] - np.asarray(y[0]) )
    lastvals = np.asarray(y[-1]) + np.abs(y[-half_window-1:-1][::-1] - np.asarray(y[-1]))
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')



def compute_alpha_voigt(wl, spec, ivar, wav_range, per_value=[1,99.5]):
    print ' Voigt Routine begins'
    spec[np.isnan(spec)] = 0
    ivar[np.isnan(ivar)] = 0
    print ivar
    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])
    #plt.plot(wl,spec)
    for j in range(len(wav_range)):
        #print wav_range[j]
        #print min(wl),max(wl)
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        #print wl[temp],len(spec),len(wl),spec[temp],ivar[temp]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
        #print len(tempspec)
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut[0]) & (tempspec < cut[1])  & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))
    print 'voigt Debug',len(wavelength)
    pv0=[20.0, 1545., 50., 8.0,1.0,-1.0]
    param_boundsv = ([0,1540,0.0,0.0,0.0,-np.inf],[np.inf,1550,np.inf,20,np.inf,np.inf])

    try:
        #plt.plot(wavelength,spectra)
        #plt.show()
        poptv, pcovv = curve_fit(myVoigt, wavelength, spectra, pv0, sigma=1.0/np.sqrt(invar),bounds=param_boundsv)
    except (RuntimeError,TypeError):
        AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv, CHISQv, DOFv = np.nan, np.nan, np.nan, np.nan,np.nan, np.nan, np.nan,  np.nan
    else:
        AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv = poptv[0], poptv[1], poptv[2], poptv[3], poptv[4], poptv[5]
        CHISQv = np.sum(invar * (spectra - myVoigt(wavelength, poptv[0], poptv[1], poptv[2], poptv[3], poptv[4], poptv[5]))**2)
        # DOF = N - n  , n = 2
        DOFv = len(spectra) - 6
    print 'Voight Routine ends' 
    print 'Compute Alpha:',  AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv, CHISQv, DOFv
    return  AMPv, CENTERv, SIGMALv, SIGMAGv,  SCALEv, ALPHAv, CHISQv, DOFv



def compute_alpha_gh(wl, spec, ivar, wav_range, per_value=[1,99.5]):
    print ' GH Routine begins'
    print ivar
    spec[np.isnan(spec)] = 0
    ivar[np.isnan(ivar)] = 0
    print ivar
    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])
    #plt.plot(wl,spec)
    for j in range(len(wav_range)):
        #print wav_range[j]
        #print min(wl),max(wl)
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        #print wl[temp],len(spec),len(wl),spec[temp],ivar[temp]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
        print len(tempspec)
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut[0]) & (tempspec < cut[1])  & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))
    print '2H Debug',len(wavelength)
    pgh0=[20.0, 1545., 8., 0.2, 0.4,21.0,-1.0]
    param_boundsgh = ([0,1540,-20,-1,-np.inf,0,-np.inf],[np.inf,1555,20,1,np.inf,np.inf,np.inf])
    try:
        #plt.plot(wavelength,spectra)
        #plt.show()
        poptgh, pcovgh = curve_fit(myGaussHermite, wavelength, spectra, pgh0, sigma=1.0/np.sqrt(invar),bounds=param_boundsgh)
    except (RuntimeError, TypeError):
        AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh, SCALEgh, ALPHAgh, CHISQgh, DOFgh = np.nan, np.nan, np.nan, np.nan,np.nan, np.nan, np.nan,  np.nan, np.nan
    else:
        AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh, SCALEgh, ALPHAgh = poptgh[0], poptgh[1], poptgh[2], poptgh[3], poptgh[4], poptgh[5], poptgh[6]
        CHISQgh = np.sum(invar * (spectra - myGaussHermite(wavelength, poptgh[0], poptgh[1], poptgh[2], poptgh[3], poptgh[4], poptgh[5], poptgh[6]))**2)
        # DOF = N - n  , n = 2
        DOFgh = len(spectra) - 7
    print 'GH Routine ends' 
    print 'Compute Alpha:', AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh,SCALEgh, ALPHAgh, CHISQgh, DOFgh
    return AMPgh, CENTERgh, SIGMAgh, SKEWgh, KURTgh,SCALEgh, ALPHAgh, CHISQgh, DOFgh


def compute_alpha_2g(wl, spec, ivar, wav_range, per_value=[1,99.5]):
    print '2G Routine begins'
    spec[np.isnan(spec)] = 0
    ivar[np.isnan(ivar)] = 0

    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])
    #plt.plot(wl,spec)
    for j in range(len(wav_range)):
        #print wav_range[j]
        #print min(wl),max(wl)
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        #print wl[temp],len(spec),len(wl),spec[temp],ivar[temp]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
        print len(tempspec)
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut[0]) & (tempspec < cut[1])  & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))
    print '2G Debug',len(wavelength)
    pg20=[10.0, 1545., 9., 10.0, 1550., 18.,20.0,1]
    param_bounds2g = ([0,1540,5,0,1540,0,0,-np.inf],[np.inf,1555,np.inf,np.inf,1555,1000,np.inf,np.inf])
    try:
        #plt.plot(wavelength,spectra)
        #plt.show()
        popt2g, pcov2g = curve_fit(myDoubleGauss, wavelength, spectra, pg20, sigma=1.0/np.sqrt(invar),bounds=param_bounds2g)
    except (RuntimeError,  TypeError):
        AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g, CHISQ2g, DOF2g = np.nan, np.nan, np.nan, np.nan,np.nan, np.nan, np.nan, np.nan, np.nan, np.nan
    else:
        AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g = popt2g[0], popt2g[1], popt2g[2], popt2g[3], popt2g[4], popt2g[5], popt2g[6], popt2g[7]
        CHISQ2g = np.sum(invar * (spectra - myDoubleGauss(wavelength, popt2g[0], popt2g[1], popt2g[2], popt2g[3], popt2g[4], popt2g[5], popt2g[6], popt2g[7]))**2)
        # DOF = N - n  , n = 2
        DOF2g = len(spectra) - 8 
    print '2G Routine ends' 
    print 'Compute Alpha:' ,AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g, CHISQ2g, DOF2g
    return AMPa2g, CENTERa2g, SIGMAa2g, AMPb2g, CENTERb2g, SIGMAb2g, SCALE2g, ALPHA2g, CHISQ2g, DOF2g 


def maskOutliers_v(wave, flux, weight, amp, center, sigmal, sigmag, scale, alpha):
    model = myVoigt(wave,amp, center, sigmal, sigmag,  scale, alpha)
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    print 'Weights Before mask',weight
    #ww = np.where (np.abs(fluxdiff) > 3*std)
    ww = np.where (((np.abs(fluxdiff) > 2*std) & (flux <= np.median(flux))) | ((np.abs(fluxdiff) > 4*std) & (flux >= np.median(flux))))
    #nwave = np.delete(wave,ww)
    #nflux = np.delete(flux,ww)
    weight[ww] = 0#np.delete(weight,ww)
    print 'Weights After mask',weight
    return wave,flux,weight


def maskOutliers_gh(wave, flux, weight, amp, center, sigma, skew, kurt, scale, alpha):
    model = myGaussHermite(wave,amp, center, sigma, skew, kurt, scale, alpha)
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    print 'Weights Before mask',weight
    #ww = np.where (np.abs(fluxdiff) > 3*std)
    ww = np.where (((np.abs(fluxdiff) > 2*std) & (flux <= np.median(flux))) | ((np.abs(fluxdiff) > 4*std) & (flux >= np.median(flux))))
    #nwave = np.delete(wave,ww)
    #nflux = np.delete(flux,ww)
    weight[ww] = 0#np.delete(weight,ww)
    print 'Weights After mask',weight
    return wave,flux,weight


def maskOutliers_2g(wave,flux,weight,amp1,center1,sigma1,amp2,center2,sigma2,scale,alpha):
    model = myDoubleGauss(wave,amp1,center1,sigma1,amp2,center2,sigma2,scale,alpha)
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    #ww = np.where (np.abs(fluxdiff) > 3*std)
    ww = np.where (((np.abs(fluxdiff) > 2*std) & (flux <= np.median(flux))) | ((np.abs(fluxdiff) > 4*std) & (flux >= np.median(flux))))
    #nwave = np.delete(wave,ww)
    #nflux = np.delete(flux,ww)
    weight[ww] = 0#np.delete(weight,ww)
    
    return wave,flux,weight


def minmax(vec):
    return min(vec),max(vec)

#Version _ Tried small sigma values for both Gaussians in GausDouble and GausHermite
#Version II Tried making one of the sigmas high, one 1 and two 54 GausDouble; Same as before for GausHermite
#Version III Changing both the central wavelength to 1550 for GausDouble
df = np.genfromtxt('tdss_allmatches_crop_edit.dat',names=['ra','dec','z','pmf1','pmf2','pmf3'],dtype=(float,float,float,'|S15','|S15','|S15'))
wav_range= [(1280,1350),(1700,1800),(1950,2200),(2650,2710),(3010,3700),(3950,4050),(4140,4270)]
pp = PdfPages('ContinuumNormalization_plus_EmissionLineFits_V.pdf')
#fx=open('ALPHA_AMP_values_V.txt','w')
#fx1=open('Kate_ALPHA_AMP_values_V.txt','w')
#for i in range(len(df['pmf1'])):
for i in range(len(df)):
#for i in range(12):
    print 'Kate_Sources/spec-'+df['pmf1'][i]+'.fits' 
    print 'Kate_Sources/spec-'+df['pmf2'][i]+'.fits' 
    print 'Kate_Sources/spec-'+df['pmf3'][i]+'.fits' 
    data1 = fits.open('Kate_Sources/spec-'+df['pmf1'][i]+'.fits')[1].data
    data2 = fits.open('Kate_Sources/spec-'+df['pmf2'][i]+'.fits')[1].data
    data3 = fits.open('Kate_Sources/spec-'+df['pmf3'][i]+'.fits')[1].data
    wave1 = 10**data1.loglam.copy()
    flux1 = data1.flux.copy()
    weight1 =  (data1.ivar).copy()
    sigma1 =  1.0/np.sqrt((data1.ivar).copy())
    wave2 = 10**data2.loglam.copy()
    flux2 = data2.flux.copy()
    weight2 =  (data2.ivar).copy()
    sigma2 =  1.0/np.sqrt((data2.ivar).copy())
    wave3 = 10**data3.loglam.copy()
    flux3 = data3.flux.copy()
    weight3 =  (data3.ivar).copy()
    sigma3 =  1.0/np.sqrt((data3.ivar).copy())
    print len(wave1),len(flux1),len(weight1)
    print weight1
    print data1.and_mask
    
    #clean QSO
    #sn1= flux1*np.sqrt(weight1) ; sn2= flux2*np.sqrt(weight2); sn3= flux3*np.sqrt(weight3)
    #w1 = (weight1>0)&((sn1<-10)|(sn1>80)); w2 = (weight2>0)&((sn2<-10)|(sn2>80)) ; w3 = (weight3>0)&((sn3<-10)|(sn3>80))
    #print w1,w2,w3
    #weight1[w1] = 0; flux1[w1] = 0
    #weight2[w2] = 0; flux2[w2] = 0
    #weight3[w3] = 0; flux3[w3] = 0
    
    #de-redden the flux
    info = fits.open('Kate_Sources/spec-'+df['pmf1'][i]+'.fits')[2].data
    coords = SkyCoord(info['RA'],info['DEC'],unit='degree',frame='icrs')
    sfd = SFDQuery()
    eb_v = sfd(coords)
    dered_flux1 = dered_flux(3.1*eb_v,wave1,flux1)
    dered_flux2 = dered_flux(3.1*eb_v,wave2,flux2)
    dered_flux3 = dered_flux(3.1*eb_v,wave3,flux3)
    
    # Change wavelengths to rest wavellegths
    rwave1 = wave1/(1.0+df['z'][i])
    rwave2 = wave2/(1.0+df['z'][i])
    rwave3 = wave3/(1.0+df['z'][i])
    
    #Clip all the spectra to same wavelength bounds Necessary to capture variation in alpha
    wmin1,wmax1 = minmax(rwave1) 
    wmin2,wmax2 = minmax(rwave2) 
    wmin3,wmax3 = minmax(rwave3)
    wlim1,wlim2 = max([wmin1,wmin2,wmin3]),min([wmax1,wmax2,wmax3])
    l1 = np.where((rwave1 >= 1410) & (rwave1 <= 1650))[0]
    l2 = np.where((rwave2 >= 1410) & (rwave2 <= 1650))[0]
    l3 = np.where((rwave3 >= 1410) & (rwave3 <= 1650))[0]
    crwave1 = rwave1[l1];cdered_flux1 = dered_flux1[l1];cweight1 = weight1[l1];csigma1= sigma1[l1]
    crwave2 = rwave2[l2];cdered_flux2 = dered_flux2[l2];cweight2 = weight2[l2];csigma2 = sigma2[l2]
    crwave3 = rwave3[l3];cdered_flux3 = dered_flux3[l3];cweight3 = weight3[l3];csigma3 = sigma3[l3]

    #Initial Mask for broad absorptions
    cut1 = np.nanpercentile(cdered_flux1, 20)
    cut2 = np.nanpercentile(cdered_flux2, 20)
    cut3 = np.nanpercentile(cdered_flux3, 20)
    xx1 = np.where(cdered_flux1 < cut1)[0]
    xx2 = np.where(cdered_flux2 < cut2)[0]
    xx3 = np.where(cdered_flux3 < cut3)[0]
    print cut1
    print dered_flux1[xx1]
    cweight1[xx1] = 0
    cweight2[xx2] = 0
    cweight3[xx3] = 0
    
        
    #Fit powerlaw iterate and mask outliers
    iteration = 3
    for j in range(iteration):
        if j == 0:
            print 'iteration 1 Begins'
            print j,'GH',df['pmf1'][i]
            AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1, CHISQgh1, DOFgh1   = compute_alpha_gh(crwave1, cdered_flux1, cweight1,  waveRange(crwave1))
            print j,'2G',df['pmf1'][i]
            AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 ,CHISQ2g1, DOF2g1   = compute_alpha_2g(crwave1, cdered_flux1, cweight1,  waveRange(crwave1))
            print j,'V',df['pmf1'][i]
            AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1, CHISQv1, DOFv1   = compute_alpha_voigt(crwave1, cdered_flux1, cweight1,  waveRange(crwave1))
            print 'Object 1 Done'
            print j,'GH',df['pmf2'][i]
            AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2, CHISQgh2, DOFgh2   = compute_alpha_gh(crwave2, cdered_flux2, cweight2,  waveRange(crwave2))
            print j,'2G',df['pmf2'][i]
            AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 ,CHISQ2g2, DOF2g2   = compute_alpha_2g(crwave2, cdered_flux2, cweight2,  waveRange(crwave2))
            print j,'V',df['pmf2'][i]
            AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2, CHISQv2, DOFv2   = compute_alpha_voigt(crwave2, cdered_flux2, cweight2,  waveRange(crwave2))
            print 'Object 2 Done'
            print j,'GH',df['pmf3'][i]
            AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3, CHISQgh3, DOFgh3   = compute_alpha_gh(crwave3, cdered_flux3, cweight3,  waveRange(crwave3))
            print j,'2G',df['pmf3'][i]
            AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 ,CHISQ2g3, DOF2g3   = compute_alpha_2g(crwave3, cdered_flux3, cweight3,  waveRange(crwave3))
            print j,'V',df['pmf3'][i]
            AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3, CHISQv3, DOFv3   = compute_alpha_voigt(crwave3, cdered_flux3, cweight3,  waveRange(crwave3))
            print 'Object 3 Done'
            nwavegh1 = crwave1; nfluxgh1=cdered_flux1;nweightgh1 = cweight1
            nwavegh2 = crwave2; nfluxgh2=cdered_flux2;nweightgh2 = cweight2
            nwavegh3 = crwave3; nfluxgh3=cdered_flux3;nweightgh3 = cweight3

            nwave2g1 = crwave1; nflux2g1=cdered_flux1;nweight2g1 = cweight1
            nwave2g2 = crwave2; nflux2g2=cdered_flux2;nweight2g2 = cweight2
            nwave2g3 = crwave3; nflux2g3=cdered_flux3;nweight2g3 = cweight3


            nwavev1 = crwave1; nfluxv1=cdered_flux1;nweightv1 = cweight1
            nwavev2 = crwave2; nfluxv2=cdered_flux2;nweightv2 = cweight2
            nwavev3 = crwave3; nfluxv3=cdered_flux3;nweightv3 = cweight3
            print 'iteration 1 Ends'
            continue
        else:
            nwavegh1,nfluxgh1,nweightgh1 = maskOutliers_gh(nwavegh1, nfluxgh1, nweightgh1, AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1)
            nwavegh2,nfluxgh2,nweightgh2 = maskOutliers_gh(nwavegh2, nfluxgh2, nweightgh2, AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2)
            nwavegh3,nfluxgh3,nweightgh3 = maskOutliers_gh(nwavegh3, nfluxgh3, nweightgh3, AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3)
            
            nwave2g1,nflux2g1,nweight2g1 = maskOutliers_2g(nwave2g1, nflux2g1, nweight2g1, AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1)
            nwave2g2,nflux2g2,nweight2g2 = maskOutliers_2g(nwave2g2, nflux2g2, nweight2g2, AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2)
            nwave2g3,nflux2g3,nweight2g3 = maskOutliers_2g(nwave2g3, nflux2g3, nweight2g3, AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3)

            nwavev1,nfluxv1,nweightv1 = maskOutliers_v(nwavev1, nfluxv1, nweightv1, AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1 )
            nwavev2,nfluxv2,nweightv2 = maskOutliers_v(nwavev2, nfluxv2, nweightv2, AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2 )
            nwavev3,nfluxv3,nweightv3 = maskOutliers_v(nwavev3, nfluxv3, nweightv3, AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3 )

            
            print j,'GH',df['pmf1'][i]
            AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1, CHISQgh1, DOFgh1   = compute_alpha_gh(nwavegh1, nfluxgh1, nweightgh1,  waveRange(nwavegh1))
            print j,'2G',df['pmf1'][i]
            AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 ,CHISQ2g1, DOF2g1   = compute_alpha_2g(nwave2g1, nflux2g1, nweight2g1,  waveRange(nwave2g1))
            print j,'V',df['pmf1'][i]
            AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1, CHISQv1, DOFv1   = compute_alpha_voigt(nwavev1, nfluxv1, nweightv1,  waveRange(nwavev1))
            print j,'GH',df['pmf2'][i]
            AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2, CHISQgh2, DOFgh2   = compute_alpha_gh(nwavegh2, nfluxgh2, nweightgh2,  waveRange(nwavegh2))
            print j,'2G',df['pmf2'][i]
            AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 ,CHISQ2g2, DOF2g2   = compute_alpha_2g(nwave2g2, nflux2g2, nweight2g2,  waveRange(nwave2g2))
            print j,'V',df['pmf2'][i]
            AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2, CHISQv2, DOFv2   = compute_alpha_voigt(nwavev2, nfluxv2, nweightv2,  waveRange(nwavev2))
            print j,'GH',df['pmf3'][i]
            AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3, CHISQgh3, DOFgh3   = compute_alpha_gh(nwavegh3, nfluxgh3, nweightgh3,  waveRange(nwavegh3))
            print j,'2G',df['pmf3'][i]
            AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 ,CHISQ2g3, DOF2g3   = compute_alpha_2g(nwave2g3, nflux2g3, nweight2g3,  waveRange(nwave2g3))
            print j,'V',df['pmf3'][i]
            AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3, CHISQv3, DOFv3   = compute_alpha_voigt(nwavev3, nfluxv3, nweightv3,  waveRange(nwavev3))



    fig,((ax1,rax1),(ax2,rax2),(ax3,rax3))=plt.subplots(3,2,figsize=(20,10))
    ax1.plot(crwave1,cdered_flux1)
    ax1.plot(crwave1,cweight1,alpha=0.2)
    ax1.plot(crwave1[cweight1 > 0],cdered_flux1[cweight1>0],'.')
    ax2.plot(crwave2,cdered_flux2)
    ax2.plot(crwave2,cweight2,alpha=0.2)
    ax2.plot(crwave2[cweight2 > 0],cdered_flux2[cweight2>0],'.')

    ax3.plot(crwave3,cdered_flux3)
    ax3.plot(crwave3,cweight3,alpha=0.2)
    ax3.plot(crwave3[cweight3 > 0],cdered_flux3[cweight3>0],'.')

    #ax1.plot(cwave1,cflux1,'--')
    pgh1=ax1.plot(crwave1,myGaussHermite(crwave1,AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1),':',color='red',lw=3,label='Gauss-Hermite')
    pgh2=ax2.plot(crwave2,myGaussHermite(crwave2,AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2),':',color='red',lw=3,label='Gauss-Hermite')
    pgh3=ax3.plot(crwave3,myGaussHermite(crwave3,AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3),':',color='red',lw=3,label='Gauss-Hermite')
    
    p2g1=ax1.plot(crwave1,myDoubleGauss(crwave1,AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 ),':',color='blue',lw=3,label='2-Gaussian')
    p2g2=ax2.plot(crwave2,myDoubleGauss(crwave2,AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 ),':',color='blue',lw=3,label='2-Gaussian')
    p2g3=ax3.plot(crwave1,myDoubleGauss(crwave3,AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 ),':',color='blue',lw=3,label='2-Gaussian')
    
    p2v1 = ax1.plot(crwave1,myVoigt(crwave1,AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1),':',color='orange',lw=3,label='Voigt')
    p2v2 = ax2.plot(crwave2,myVoigt(crwave2,AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2),':',color='orange',lw=3,label='Voigt')
    p2v3 = ax3.plot(crwave3,myVoigt(crwave3,AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3),':',color='orange',lw=3,label='Voigt')
    ax1.axhline(cut1,ls=':')
    ax2.axhline(cut1,ls=':')
    ax3.axhline(cut1,ls=':')
    string1gh = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPgh1,CENTERgh1,(CHISQgh1/DOFgh1))
    string12g = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPa2g1,CENTERa2g1,(CHISQ2g1/ DOF2g1))
    string2gh = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPgh2,CENTERgh2,(CHISQgh2/ DOFgh2))
    string22g = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPa2g2,CENTERa2g2,(CHISQ2g2/ DOF2g2))
    string3gh = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPgh3,CENTERgh3,(CHISQgh3/ DOFgh3))
    string32g = 'AMP: {0:4.3f}   CENTER: {1:4.3f} RCHIS:{2:4.3f}'.format(AMPa2g3,CENTERa2g3,(CHISQ2g3/ DOF2g3))
    ax1.set_xlim(np.min(crwave1),np.max(crwave1))
    ax2.set_xlim(np.min(crwave1),np.max(crwave1))
    ax3.set_xlim(np.min(crwave1),np.max(crwave1))
    xlim=ax1.get_xlim()
    ylim1 = ax1.get_ylim()
    ylim2 = ax1.get_ylim()
    ylim3 = ax1.get_ylim()
    ax1.set_ylim(min(savitzky_golay(cdered_flux1,101,2))-3*np.std(cdered_flux1),max(savitzky_golay(cdered_flux1,101,2))+5*np.std(cdered_flux1))
    ax2.set_ylim(min(savitzky_golay(cdered_flux2,101,2))-3*np.std(cdered_flux2),max(savitzky_golay(cdered_flux2,101,2))+5*np.std(cdered_flux2))
    ax3.set_ylim(min(savitzky_golay(cdered_flux3,101,2))-3*np.std(cdered_flux3),max(savitzky_golay(cdered_flux3,101,2))+5*np.std(cdered_flux3))
    ax1.annotate(str(df['pmf1'][i]) ,xy=(.95,.15), xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax2.annotate(str(df['pmf2'][i]),xy=(.95,.15), xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax3.annotate(str(df['pmf3'][i]),xy=(.95,.15), xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    #ax1.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim1[1] - 0.2*(ylim1[1] - ylim1[0]),string12g)
    #ax2.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim2[1] - 0.1*(ylim1[1] - ylim2[0]),string2gh)
    #ax2.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim2[1] - 0.2*(ylim1[1] - ylim2[0]),string22g)
    #ax3.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim3[1] - 0.1*(ylim1[1] - ylim3[0]),string3gh)
    #ax3.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim3[1] - 0.2*(ylim1[1] - ylim2[0]),string32g)
    #ax1.legend([pgh1,p2g1,p2v1],['Gaus-Hermite','2-Gaus','Voigt'],loc=3)
    ax1.legend(loc=3)
    ax1.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f\nrChi2 = %.3f' \
        %(AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1,(CHISQgh1/ DOFgh1)),xy=(.05,.95), \
        xycoords='axes fraction',ha="left", va="top", \
    bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax1.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f\nrChi2 = %.3f' \
        %(AMPa2g1, AMPb2g1, CENTERa2g1, CENTERb2g1, SIGMAa2g1, SIGMAb2g1,(CHISQ2g1/ DOF2g1)),xy=(.95,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax1.annotate('Voigt:\nAmp$_1$ = %.2f\nCenter = %.2f\n$\sigma_l$ = %.2f\n$\sigma_g$ = %.2f\nrChi2 = %.3f' \
        %(AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,(CHISQv1/ DOFv1)),xy=(.35,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)

    ax2.legend(loc=3)
    #ax2.legend([pgh2,p2g2,p2v2],['Gaus-Hermite','2-Gaus','Voigt'],prop={'size':8},loc='center left')
    ax2.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f\nrChi2 = %.3f'\
        %(AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2,(CHISQgh2/ DOFgh2)),xy=(.05,.95), \
        xycoords='axes fraction',ha="left", va="top", \
    bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax2.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f\nrChi2 = %.3f' \
        %(AMPa2g2, AMPb2g2, CENTERa2g2, CENTERb2g2, SIGMAa2g2, SIGMAb2g2,(CHISQ2g2/ DOF2g2)),xy=(.95,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax2.annotate('Voigt:\nAmp$_1$ = %.2f\nCenter = %.2f\n$\sigma_l$ = %.2f\n$\sigma_g$ = %.2f\nrChi2 = %.3f' \
        %(AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,(CHISQv2/ DOFv2)),xy=(.35,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)

    ax3.legend(loc=3)
    #ax3.legend([pgh3,p2g3,p2v3],['Gaus-Hermite','2-Gaus','Voigt'],prop={'size':8},loc='center left')
    ax3.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f\nrChi2 = %.3f' \
        %(AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3,(CHISQgh3/ DOFgh3)),xy=(.05,.95), \
        xycoords='axes fraction',ha="left", va="top", \
    bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax3.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f\nrChi2 = %.3f' \
        %(AMPa2g3, AMPb2g3, CENTERa2g3, CENTERb2g3, SIGMAa2g3, SIGMAb2g3,(CHISQ2g3/ DOF2g3)),xy=(.95,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    ax3.annotate('Voigt:\nAmp$_1$ = %.2f\nCenter = %.2f\n$\sigma_l$ = %.2f\n$\sigma_g$ = %.2f\nrChi2 = %.3f' \
        %(AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,(CHISQv3/ DOFv3)),xy=(.35,.95), \
        xycoords='axes fraction',ha="right", va="top", \
        bbox=dict(boxstyle="round", fc='1'),fontsize=10)
    
    ax1.set_ylabel('Flux')
    ax1.set_xlabel('Wavelength ($\AA$)')
    ax2.set_ylabel('Flux')
    ax2.set_xlabel('Wavelength ($\AA$)')
    ax3.set_ylabel('Flux')
    ax3.set_xlabel('Wavelength ($\AA$)')
    #Save the spectra
    filenamegh1 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf1'][i]+'_cEm_gh.txt'
    filenamegh2 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf2'][i]+'_cEm_gh.txt'
    filenamegh3 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf3'][i]+'_cEm_gh.txt'
    
    filename2g1 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf1'][i]+'_cEm_2g.txt'
    filename2g2 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf2'][i]+'_cEm_2g.txt'
    filename2g3 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf3'][i]+'_cEm_2g.txt'

    filenamev1 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf1'][i]+'_cEm_v.txt'
    filenamev2 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf2'][i]+'_cEm_v.txt'
    filenamev3 = 'EmissionLines_Norm_Spectra/NormSpec_'+df['pmf3'][i]+'_cEm_v.txt'


    # Residuals & Weights
    resgh1 = cdered_flux1 / myGaussHermite(crwave1,AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1)
    resgh2 = cdered_flux2 / myGaussHermite(crwave2,AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2)
    resgh3 = cdered_flux3 / myGaussHermite(crwave3,AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3)
    
    wresgh1 = csigma1 / myGaussHermite(crwave1,AMPgh1, CENTERgh1, SIGMAgh1, SKEWgh1, KURTgh1, SCALEgh1, ALPHAgh1)
    wresgh2 = csigma2 / myGaussHermite(crwave2,AMPgh2, CENTERgh2, SIGMAgh2, SKEWgh2, KURTgh2, SCALEgh2, ALPHAgh2)
    wresgh3 = csigma3 / myGaussHermite(crwave3,AMPgh3, CENTERgh3, SIGMAgh3, SKEWgh3, KURTgh3, SCALEgh3, ALPHAgh3)

    res2g1 = cdered_flux1 / myDoubleGauss(crwave1,AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 )
    res2g2 = cdered_flux2 / myDoubleGauss(crwave2,AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 )
    res2g3 = cdered_flux3 / myDoubleGauss(crwave3,AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 )

    wres2g1 = csigma1 / myDoubleGauss(crwave1,AMPa2g1, CENTERa2g1, SIGMAa2g1, AMPb2g1, CENTERb2g1, SIGMAb2g1, SCALE2g1, ALPHA2g1 )
    wres2g2 = csigma2 / myDoubleGauss(crwave2,AMPa2g2, CENTERa2g2, SIGMAa2g2, AMPb2g2, CENTERb2g2, SIGMAb2g2, SCALE2g2, ALPHA2g2 )
    wres2g3 = csigma3 / myDoubleGauss(crwave3,AMPa2g3, CENTERa2g3, SIGMAa2g3, AMPb2g3, CENTERb2g3, SIGMAb2g3, SCALE2g3, ALPHA2g3 )

    resv1 = cdered_flux1 / myVoigt(crwave1,AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1)
    resv2 = cdered_flux2 / myVoigt(crwave2,AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2)
    resv3 = cdered_flux3 / myVoigt(crwave3,AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3)

    wresv1 = csigma1 / myVoigt(crwave1,AMPv1, CENTERv1, SIGMALv1, SIGMAGv1,  SCALEv1, ALPHAv1)
    wresv2 = csigma2 / myVoigt(crwave2,AMPv2, CENTERv2, SIGMALv2, SIGMAGv2,  SCALEv2, ALPHAv2)
    wresv3 = csigma3 / myVoigt(crwave3,AMPv3, CENTERv3, SIGMALv3, SIGMAGv3,  SCALEv3, ALPHAv3)

    #Save Begins
    np.savetxt(filenamegh1,zip(crwave1,resgh1,wresgh1), fmt='%10.5f')
    np.savetxt(filenamegh2,zip(crwave2,resgh2,wresgh2), fmt='%10.5f')
    np.savetxt(filenamegh3,zip(crwave3,resgh3,wresgh3), fmt='%10.5f')

    np.savetxt(filename2g1,zip(crwave1,res2g1,wres2g1), fmt='%10.5f')
    np.savetxt(filename2g2,zip(crwave2,res2g2,wres2g2), fmt='%10.5f')
    np.savetxt(filename2g3,zip(crwave3,res2g3,wres2g3), fmt='%10.5f')

    np.savetxt(filenamev1,zip(crwave1,resv1,wresv1), fmt='%10.5f')
    np.savetxt(filenamev2,zip(crwave2,resv2,wresv2), fmt='%10.5f')
    np.savetxt(filenamev3,zip(crwave3,resv3,wresv3), fmt='%10.5f')

    #Plot Begins
    presgh1,pres2g1,presv1,=rax1.plot(crwave1,resgh1,'k--',crwave1,res2g1,'r-',crwave1,resv1,'b:',alpha=0.7)
    presgh2,pres2g2,presv2,=rax2.plot(crwave2,resgh2,'k--',crwave2,res2g2,'r-',crwave2,resv2,'b:',alpha=0.7)
    presgh3,pres2g3,presv3,=rax3.plot(crwave3,resgh3,'k--',crwave3,res2g3,'r-',crwave3,resv3,'b:',alpha=0.7)

    rax1.set_ylabel('Normalized Flux')
    rax1.set_xlabel('Wavelength ($\AA$)')
    rax1.legend([presgh1,pres2g1,presv1],['Gaus-Hermite','2-Gaus','Voigt'],numpoints=4,prop={'size':8},loc='lower right')
    rax2.set_ylabel('Normalized Flux')
    rax2.set_xlabel('Wavelength ($\AA$)')
    rax2.legend([presgh2,pres2g2,presv2],['Gaus-Hermite','2-Gaus','Voigt'],numpoints=4,prop={'size':8},loc='lower right')
    rax3.set_ylabel('Normalized Flux')
    rax3.set_xlabel('Wavelength ($\AA$)')
    rax3.legend([presgh3,pres2g3,presv3],['Gaus-Hermite','2-Gaus', 'Voigt'],numpoints=4,prop={'size':8},loc='lower right')
    rax1.axhline(1.0,ls=':',color='blue',alpha=0.5)
    rax2.axhline(1.0,ls=':',color='blue',alpha=0.5)
    rax3.axhline(1.0,ls=':',color='blue',alpha=0.5)
    rax1.set_ylim(0,2)
    rax2.set_ylim(0,2)
    rax3.set_ylim(0,2)
    fig.tight_layout()
    fig.savefig(pp,format='pdf')
    #plt.show()
    plt.clf()
pp.close()

<file_sep>#!/usr/bin/python
#-*- coding: utf-8 -*-

'''
Code for continuum normalization of SDSS spectra
@author   : <NAME>
@date     : 21-Jan-2018
@version  : 1.0
'''



import sys
import os
import numpy as np
from astropy import io
from astropy.io import fits 
from astropy.io import ascii
from matplotlib import pyplot as plt
from astropy import constants as const
from astropy import units as U
from astropy.coordinates import SkyCoord
from dustmaps.sfd import SFDQuery
from specutils import extinction  
import scipy
from scipy.ndimage.filters import convolve
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit 
from matplotlib.backends.backend_pdf import PdfPages
from con_fit_iterative import *

def lam2vel(wavelength,ion_w):
    zlambda = (wavelength-ion_w)/ion_w
    R = 1./(1+zlambda)
    vel_ion = -const.c.to('km/s').value*(R**2-1)/(R**2+1)
    return vel_ion



def dered_flux(Av,wave,flux):
    dered_flux = np.zeros(len(wave))
    for i in range(len(wave)):
        ext_lambda= extinction.extinction_ccm89(wave[i] * U.angstrom, a_v= Av, r_v= 3.1)
        tau_lambda= ext_lambda/(1*1.086) 
        dered_flux[i]= flux[i] * np.exp(tau_lambda)
    return dered_flux

def ext_coeff(lamb):
    
    inv_lamba=[0.45,0.61,0.8,1.82,2.27,2.7,3.22,3.34,3.46,3.6,3.75,3.92,4.09,4.28,4.50,4.73,5.00,5.24,5.38,5.52,5.70,5.88,6.07,6.27,6.48,6.72,6.98,7.23,7.52,7.84]
    smc_ext=[-2.61,-2.47,-2.12,0.0,1.0,1.67,2.29,2.65,3.0,3.15,3.49,3.91,4.24,4.53,5.3,5.85,6.38,6.76,6.9,7.17,7.71,8.01,8.49,9.06,9.28,9.84,10.8,11.51,12.52,13.54]
    xy=np.interp(1.0/(lamb*10**(-4)),inv_lamba,smc_ext)
    ext_lamb=(xy+2.98)/3.98 # Rv=2.98
    #return(z,ext_lamb)
    return(ext_lamb)

def myfunct(wave, a , b):
    return a*ext_coeff(wave)*wave**(b)


def compute_alpha(wl, spec, ivar, wav_range, per_value=[10,90]):
    # print 'Routine begins'
    spec[np.isnan(spec)] = 0
    ivar[np.isnan(ivar)] = 0

    wavelength, spectra, invar = np.array([]), np.array([]), np.array([])
    #plt.plot(wl,spec)
    for j in range(len(wav_range)):
        #print wav_range[j]
        #print min(wl),max(wl)
        temp = np.where((wl > wav_range[j][0]) & (wl < wav_range[j][1]))[0]
        #print wl[temp],len(spec),len(wl),spec[temp],ivar[temp]
        tempspec, tempivar  = spec[temp], ivar[temp]
        #print tempspec
        print len(tempspec)
        #Mask out metal absorption lines
        cut = np.percentile(tempspec, per_value)
        #print 'cut',cut
        blah = np.where((tempspec > cut[0]) & (tempspec < cut[1])  & (tempivar > 0))[0]
        wave = wl[temp][blah]

        wavelength = np.concatenate((wavelength, wave))
        spectra = np.concatenate((spectra, tempspec[blah]))
        invar = np.concatenate((invar, tempivar[blah]))
    print 'Debug',len(wavelength)
    p0=[1.0,1.0]
    try:
        #plt.plot(wavelength,spectra)
        #plt.show()
        popt, pcov = curve_fit(myfunct, wavelength, spectra, p0, sigma=1.0/np.sqrt(invar))
    except (RuntimeError, TypeError):
        AMP, ALPHA, CHISQ, DOF = np.nan, np.nan, np.nan, np.nan
    else:
        AMP, ALPHA = popt[0], popt[1]
        CHISQ = np.sum(invar * (spectra - myfunct(wavelength, popt[0], popt[1]))**2)
        # DOF = N - n  , n = 2
        DOF = len(spectra) - 2
    # print 'Routine ends' 
    print 'Compute Alpha:', AMP,ALPHA,CHISQ,DOF
    return AMP, ALPHA, CHISQ, DOF


def waveRange(wave):
    #wav_range=[[1300.,1350.],[1420.,1470.],[1700.,1800.],[2080.,2215],[2480.,2655],[3225.,3900.],[4200.,4230.],[4435.,4700.],[5200.,5700.]]
    wav_range= [(1250,1350),(1700,1800),(1950,2200),(2650,2710),(2950,3700),(3950,4050)]


    mwav_range = []
    for rr in wav_range:
        if ((np.min(wave) < rr[0]) & (np.max(wave) > rr[1])):
            nrr0 = rr[0] ; nrr1 = rr[1]
        elif ((np.min(wave) > rr[0]) & (np.min(wave) < rr[1])& (np.max(wave) < rr[1])):
            nrr0 = np.min(wave) ; nrr1 = rr[1]
        elif ((np.min(wave) < rr[0]) & (np.max(wave) > rr[0]) &(np.max(wave) < rr[1])):
            nrr0 = rr[0] ; nrr1 = np.max(wave)
        else :
            continue
        mwav_range.append([nrr0,nrr1])
    return mwav_range


def powerlawFunc(xdata, scale, amp,index):
    return scale + amp*ext_coeff(xdata)*np.power(xdata,index)


def fitPowerlaw(wave,flux,weight,scale = 1, amp=1,index=1): 
    from lmfit import minimize, Parameters, fit_report, Minimizer
    import numpy as np
    import scipy.optimize as optimization
    x0= [scale, amp,index]
    xdata=np.asarray(wave)
    ydata=np.asarray(flux)
    sigma=np.asarray(weight)
    print len(xdata),len(ydata),len(sigma) 
        
    popt, pcov = optimization.curve_fit(powerlawFunc, xdata, ydata, x0, sigma)
    print popt
    #popt, pcov = optimization.curve_fit(func, xdata, ydata, x0)
    model = powerlawFunc(wave,popt[0],popt[1],popt[2])
    chi2 = ((flux - model)*np.sqrt(weight))**2
    rchi2 = np.sum(chi2)/(len(xdata) - 2)
    print 'Reduced Chi Square : {0}  Number of points: {1}'.format(rchi2,len(xdata))
    return (popt,pcov)   

def maskOutliers(wave,flux,weight,amp,alpha):
    model = myfunct(wave,amp,alpha)
    std =np.std(flux[weight > 0])
    fluxdiff = flux - model
    ww = np.where (np.abs(fluxdiff) > 3*std)
    #nwave = np.delete(wave,ww)
    #nflux = np.delete(flux,ww)
    weight[ww] = 0#np.delete(weight,ww)
    
    return wave,flux,weight

def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
    The Savitzky-Golay filter removes high frequency noise from data.
    It has the advantage of preserving the original shape and
    features of the signal better than other types of filtering
    approaches, such as moving averages techniques.
    Parameters
    ----------
    y : array_like, shape (N,)
        the values of the time history of the signal.
    window_size : int
        the length of the window. Must be an odd integer number.
    order : int
        the order of the polynomial used in the filtering.
        Must be less then `window_size` - 1.
    deriv: int
        the order of the derivative to compute (default = 0 means only smoothing)
    Returns
    -------
    ys : ndarray, shape (N)
        the smoothed signal (or it's n-th derivative).
    Notes
    -----
    The Savitzky-Golay is a type of low-pass filter, particularly
    suited for smoothing noisy data. The main idea behind this
    approach is to make for each point a least-square fit with a
    polynomial of high order over a odd-sized window centered at
    the point.
    Examples
    --------
    t = np.linspace(-4, 4, 500)
    y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
    ysg = savitzky_golay(y, window_size=31, order=4)
    import matplotlib.pyplot as plt
    plt.plot(t, y, label='Noisy signal')
    plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
    plt.plot(t, ysg, 'r', label='Filtered signal')
    plt.legend()
    plt.show()
    References
    ----------
    .. [1] <NAME>, <NAME>, Smoothing and Differentiation of
       Data by Simplified Least Squares Procedures. Analytical
       Chemistry, 1964, 36 (8), pp 1627-1639.
    .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing
       W.H. Press, <NAME>, <NAME>, <NAME>
       Cambridge University Press ISBN-13: 9780521880688
    """
    import numpy as np
    from math import factorial
    
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except (ValueError, msg):
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    # precompute coefficients
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    # pad the signal at the extremes with
    # values taken from the signal itself
    firstvals = np.asarray(y[0]) - np.abs( y[1:half_window+1][::-1] - np.asarray(y[0]) )
    lastvals = np.asarray(y[-1]) + np.abs(y[-half_window-1:-1][::-1] - np.asarray(y[-1]))
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')





def contWave(wave):
    linefree_regions = [(1280,1350),(1700,1800),(1950,2200),(2650,2710),(2950,3700),(3950,4050)]
    finalcond = False
    for lfr in linefree_regions:
        cond = ((wave >= lfr[0]) & (wave <= lfr[1]))
        finalcond = finalcond | cond
    indices = np.where(finalcond)
    return indices


def minmax(vec):
    return min(vec),max(vec)

#version I = no maskOutliers percentile 25 75
#version II = with maskOutliers percentile 10 90
#version III = clipped wavelengths to same bounds; to deal with the difference in spectral indices
#version IV = Kate model included, ALPHA_AMP_values file is not different from III
#version V =Reddened power law

df = np.genfromtxt('tdss_allmatches_crop_edit.dat',names=['ra','dec','z','pmf1','pmf2','pmf3'],dtype=(float,float,float,'|S15','|S15','|S15'))
wav_range= [(1280,1350),(1700,1800),(1950,2200),(2650,2710),(3010,3700),(3950,4050),(4140,4270)]
pp = PdfPages('ContinuumNormalization_plots_V.pdf')
fx=open('ALPHA_AMP_values_V.txt','w')
fx1=open('Kate_ALPHA_AMP_values_V.txt','w')
#for i in range(len(df['pmf1'])):
for i in range(len(df)):
#for i in range(10):

    print 'Kate_Sources/spec-'+df['pmf1'][i]+'.fits' 
    print 'Kate_Sources/spec-'+df['pmf2'][i]+'.fits' 
    print 'Kate_Sources/spec-'+df['pmf3'][i]+'.fits' 
    data1 = fits.open('Kate_Sources/spec-'+df['pmf1'][i]+'.fits')[1].data
    data2 = fits.open('Kate_Sources/spec-'+df['pmf2'][i]+'.fits')[1].data
    data3 = fits.open('Kate_Sources/spec-'+df['pmf3'][i]+'.fits')[1].data
    wave1 = 10**data1.loglam.copy()
    flux1 = data1.flux.copy()
    weight1 =  (data1.ivar*(data1.and_mask == 0)).copy()
    wave2 = 10**data2.loglam.copy()
    flux2 = data2.flux.copy()
    weight2 =  (data2.ivar*(data2.and_mask == 0)).copy()
    wave3 = 10**data3.loglam.copy()
    flux3 = data3.flux.copy()
    weight3 =  (data3.ivar*(data3.and_mask == 0)).copy()
    print len(wave1),len(flux1),len(weight1)
    print weight1
    print data1.and_mask
    
    #clean QSO
    sn1= flux1*np.sqrt(weight1) ; sn2= flux2*np.sqrt(weight2); sn3= flux3*np.sqrt(weight3)
    w1 = (weight1>0)&((sn1<-10)|(sn1>80)); w2 = (weight2>0)&((sn2<-10)|(sn2>80)) ; w3 = (weight3>0)&((sn3<-10)|(sn3>80))
    print w1,w2,w3
    weight1[w1] = 0; flux1[w1] = 0
    weight2[w2] = 0; flux2[w2] = 0
    weight3[w3] = 0; flux3[w3] = 0
    
    #de-redden the flux
    info = fits.open('Kate_Sources/spec-'+df['pmf1'][i]+'.fits')[2].data
    coords = SkyCoord(info['RA'],info['DEC'],unit='degree',frame='icrs')
    sfd = SFDQuery()
    eb_v = sfd(coords)
    dered_flux1 = dered_flux(3.1*eb_v,wave1,flux1)
    dered_flux2 = dered_flux(3.1*eb_v,wave2,flux2)
    dered_flux3 = dered_flux(3.1*eb_v,wave3,flux3)
    
    # Change wavelengths to rest wavellegths
    rwave1 = wave1/(1.0+df['z'][i])
    rwave2 = wave2/(1.0+df['z'][i])
    rwave3 = wave3/(1.0+df['z'][i])
    
    #Clip all the spectra to same wavelength bounds Necessary to capture variation in alpha
    wmin1,wmax1 = minmax(rwave1) 
    wmin2,wmax2 = minmax(rwave2) 
    wmin3,wmax3 = minmax(rwave3)
    wlim1,wlim2 = max([wmin1,wmin2,wmin3]),min([wmax1,wmax2,wmax3])
    l1 = np.where((rwave1 >= wlim1) & (rwave1 <= wlim2))[0]
    l2 = np.where((rwave2 >= wlim1) & (rwave2 <= wlim2))[0]
    l3 = np.where((rwave3 >= wlim1) & (rwave3 <= wlim2))[0]
    crwave1 = rwave1[l1];cdered_flux1 = dered_flux1[l1];cweight1 = weight1[l1]
    crwave2 = rwave2[l2];cdered_flux2 = dered_flux2[l2];cweight2 = weight2[l2]
    crwave3 = rwave3[l3];cdered_flux3 = dered_flux3[l3];cweight3 = weight3[l3]

    # Choose line free region
    #cwave1  = rwave1[contWave(rwave1)];cflux1 = flux1[contWave(rwave1)] ; cweight1 = weight1[contWave(rwave1)]
    #cwave2  = rwave2[contWave(rwave2)];cflux2 = flux2[contWave(rwave2)] ; cweight2 = weight2[contWave(rwave2)]
    #cwave3  = rwave3[contWave(rwave3)];cflux3 = flux3[contWave(rwave3)] ; cweight3 = weight3[contWave(rwave3)]
    
    #Fit powerlaw iterate and mask outliers
    iteration = 3
    for j in range(iteration):
        if j == 0:
            AMP1, ALPHA1, CHISQ1, DOF1 = compute_alpha(crwave1, cdered_flux1, cweight1, waveRange(crwave1))
            AMP2, ALPHA2, CHISQ2, DOF2 = compute_alpha(crwave2, cdered_flux2, cweight2, waveRange(crwave2))
            AMP3, ALPHA3, CHISQ3, DOF3 = compute_alpha(crwave3, cdered_flux3, cweight3, waveRange(crwave3))
            print AMP1,ALPHA1
            nwave1 = crwave1; nflux1=cdered_flux1;nweight1 = cweight1
            nwave2 = crwave2; nflux2=cdered_flux2;nweight2 = cweight2
            nwave3 = crwave3; nflux3=cdered_flux3;nweight3 = cweight3
            print 'iteration number',j+1#,popt1[0],popt1[1],popt1[2]
            print 'iteration 1 completed'
            continue
        else:
            nwave1,nflux1,nweight1 = maskOutliers(nwave1,nflux1,nweight1,AMP1,ALPHA1)
            AMP1, ALPHA1, CHISQ1, DOF1 = compute_alpha(nwave1, nflux1, nweight1, waveRange(nwave1))
            nwave2,nflux2,nweight2 = maskOutliers(nwave2,nflux2,nweight2,AMP2,ALPHA2)
            AMP2, ALPHA2, CHISQ2, DOF2 = compute_alpha(nwave2, nflux2, nweight2, waveRange(nwave2))
            nwave3,nflux3,nweight3 = maskOutliers(nwave3,nflux3,nweight3,AMP3,ALPHA3)
            AMP3, ALPHA3, CHISQ3, DOF3 = compute_alpha(nwave3, nflux3, nweight3, waveRange(nwave3))
           # print 'iteration number',j+1,popt1[0],popt1[1],popt1[2]
           # nwave1,nflux1,nweight1 = maskOutliers(cwave1,cflux1*1e15,cweight1,popt1)
           # popt1,pcov1 = fitPowerlaw(nwave1,nflux1,nweight1,popt1[0],popt1[1],popt1[2])
           # nwave2,nflux2,nweight2 = maskOutliers(cwave2,cflux2*1e15,cweight2,popt2)
           # popt2,pcov2 = fitPowerlaw(nwave2,nflux2,nweight2,popt2[0].popt2[1],popt2[2])
           # nwave3,nflux3,nweight3 = maskOutliers(cwave3,cflux3*1e15,cweight3,popt3)
           # popt3,pcov3 = fitPowerlaw(nwave3,nflux3,nweight3,popt3[0],popt3[1],popt3[2])
        

    print>>fx,'{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}'.format(df['pmf1'][i],ALPHA1,ALPHA2,ALPHA3,AMP1,AMP2,AMP3)
    #Try Kate's continuum normalization fit
    kcwave1,kcflux1,kcfluxerr1,kmean_amp1,kmean_index1,kredco1,kwave1,kflux1=continuum_fit(wave1,flux1,weight1,df['z'][i])
    kcwave2,kcflux2,kcfluxerr2,kmean_amp2,kmean_index2,kredco2,kwave2,kflux2=continuum_fit(wave2,flux2,weight2,df['z'][i])
    kcwave3,kcflux3,kcfluxerr3,kmean_amp3,kmean_index3,kredco3,kwave3,kflux3=continuum_fit(wave3,flux3,weight3,df['z'][i])
    print>>fx1,'{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}'.format(df['pmf1'][i],kmean_index1,kmean_index2,kmean_index3,kmean_amp1,kmean_amp2,kmean_amp3)
    #
    #Plot for testing
    fig,(ax1,ax2,ax3)=plt.subplots(3,1,figsize=(15,8))
    ax1.plot(crwave1,cdered_flux1,label=str(df['pmf1'][i]))
    #ax1.plot(cwave1,cflux1,'--')
    ax1.plot(crwave1,myfunct(crwave1,AMP1,ALPHA1),':',color='red',lw=3,label='my model')
    ax1.plot(kcwave1/(1.0+df['z'][i]),kcflux1,'--',color='cyan',lw=3,label='Kate model')
    #ax1.plot(rwave1,dered_flux1,'--')
    #ax1.plot(rwave1,weight1,':',label='weight')
    #ax1.plot(rwave1,data1.flux,label='Clipped')
    string1 = 'AMP: {0:4.3f}   ALPHA: {1:4.3f}   rCHISQ: {2:4.3f}'.format(AMP1,ALPHA1,CHISQ1/DOF1)
    
    for ll in wav_range :
            ax1.axvspan(ll[0], ll[1], alpha=0.25, color='cyan')
            ax2.axvspan(ll[0], ll[1], alpha=0.25, color='cyan')
            ax3.axvspan(ll[0], ll[1], alpha=0.25, color='cyan')
    ax2.plot(crwave2,cdered_flux2,label=str(df['pmf2'][i]))
    #ax2.plot(cwave2,cflux2,'--')
    ax2.plot(crwave2,myfunct(crwave2,AMP2,ALPHA2),':',color='red',lw=3,label='my model')
    ax2.plot(kcwave2/(1.0+df['z'][i]),kcflux2,'--',color='cyan',lw=3,label='Kate model')
    #ax2.plot(rwave2,dered_flux2,'--')
    #ax2.plot(rwave2,weight2,':',label='weight')
    #ax2.plot(rwave2,data2.flux,'--',label='Clipped')
    string2 = 'AMP: {0:4.3f}   ALPHA: {1:4.3f}   rCHISQ: {2:4.3f}'.format(AMP2,ALPHA2,CHISQ2/DOF2)

    ax3.plot(crwave3,cdered_flux3,label=str(df['pmf3'][i]))
    #ax3.plot(cwave3,cflux3,'--')
    ax3.plot(crwave3,myfunct(crwave3,AMP3,ALPHA3),':',color='red',lw=3,label='my model')
    ax3.plot(kcwave3/(1.0+df['z'][i]),kcflux3,'--',color='cyan',lw=3,label='Kate model')
    #ax3.plot(rwave3,dered_flux3,'--')
    #ax3.plot(rwave3,weight3,':',label='weight')
    #ax3.plot(rwave3,data3.flux,'--',label='Clipped')
    string3 = 'AMP: {0:4.3f}   ALPHA: {1:4.3f}   rCHISQ: {2:4.3f}'.format(AMP3,ALPHA3,CHISQ3/DOF3)

    ax1.set_xlim(np.min(crwave1),np.max(crwave1))
    ax2.set_xlim(np.min(crwave1),np.max(crwave1))
    ax3.set_xlim(np.min(crwave1),np.max(crwave1))
    xlim=ax1.get_xlim()
    ylim1 = ax1.get_ylim()
    ylim2 = ax2.get_ylim()
    ylim3 = ax3.get_ylim()
    ax1.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim1[1] - 0.1*(ylim1[1] - ylim1[0]),string1)
    ax2.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim2[1] - 0.1*(ylim2[1] - ylim2[0]),string2)
    ax3.text(xlim[0]+0.1*(xlim[1] - xlim[0]),ylim3[1] - 0.1*(ylim3[1] - ylim3[0]),string3)

    ax1.legend(loc=1)
    ax2.legend(loc=1)
    ax3.legend(loc=1)
    fig.tight_layout()
    fig.savefig(pp,format='pdf')
    #plt.show()
fx.close()
fx1.close()
pp.close()
<file_sep>#!/usr/bin/python
#-*- coding: utf-8 -*-

import scipy
import numpy as np

class Observation(object):
    '''
    The Observation object
    '''
    def __init__(self):
        '''
        Initialize Observation object
        '''
        self._pmf= None
        self._snr = None
        self._survey = None
        self._filename = None
       
    @property
    def pmf(self):
        return self._pmf
        
    @pmf.setter
    def pmf(self, new_pmf):
        self._pmf = new_pmf
        
    @property
    def snr(self):
        return self._snr
        
    @snr.setter
    def snr(self, new_snr):
        self._snr = new_snr

    @property
    def survey(self):
        return self._survey
        
    @survey.setter
    def survey(self, new_survey):
        self._survey = new_survey

    @property
    def filename(self):
        return self._filename
        
    @filename.setter
    def filename(self, new_filename):
        self._filename = new_filename
 
    @property
    def redshift(self):
        return self._redshift
        
    @redshift.setter
    def redshift(self, new_redshift):
        self._redshift = new_redshift

    @property
    def bi(self):
        return self._bi
        
    @bi.setter
    def bi(self, new_bi):
        self._bi = new_bi

    @property
    def bi_err(self):
        return self._bi_err
        
    @bi_err.setter
    def bi_err(self, new_bi_err):
        self._bi_err = new_bi_err

    @property
    def v_max(self):
        return self._v_max
        
    @v_max.setter
    def v_max(self, new_v_max):
        self._v_max = new_v_max

    @property
    def v_min(self):
        return self._v_min
        
    @v_min.setter
    def v_min(self, new_v_min):
        self._v_min = new_v_min

    @property
    def mjd(self):
        return self._mjd
        
    @mjd.setter
    def mjd(self, new_mjd):
        self._mjd = new_mjd

    @property
    def spec_name(self):
        return self._spec_name
        
    @spec_name.setter
    def spec_name(self, new_spec_name):
        self._spec_name = new_spec_name

    @property
    def num_troughs(self):
        return self._num_troughs
        
    @num_troughs.setter
    def num_troughs(self, new_num_troughs):
        self._num_troughs = new_num_troughs

    @property
    def ew(self):
        return self._ew
        
    @ew.setter
    def ew(self, new_ew):
        self._ew = new_ew

    @property
    def ew_err(self):
        return self._ew_err
        
    @ew_err.setter
    def ew_err(self, new_err_ew_err):
        self._ew_err = new_err_ew_err





