<reponame>EduardoCristo/CaRM
import numpy as np
from astropy.time import Time
from math import floor
import glob
from astropy.io import fits
from .bouchy_error import bouchy_err
import os
from .constants import storbitpar, pathlist, spect,ccfext
#from .constants import tepoch,P,r1, parguess, spect, pathlist
import subprocess
import scipy.optimize as op
from copy import deepcopy as dpcy
import matplotlib.pyplot as plt
import sys
import scipy.optimize as opt
import string
from .orders import orderl
#from joblib import Parallel, delayed
import multiprocessing
import emcee
from .priors import Gprior, Uprior


def get_data(path, st_type, orders):
    print(path)
    rvc = []
    rvc_err = []
    phase = []
    if spect == "HARPS":
        datafits = "/*ccf_"+str(st_type)+"_A.fits"
        dataspec = "/*e2ds_A.fits"

        data_path = sorted(glob.glob(path + datafits))
        dataspec_path = sorted(glob.glob(path + dataspec))
        dataspec_path = dpcy(dataspec_path[1:])

        # Calibration file
        cal_wav_ext = "/*wave_A.fits"
        cal_wav_file = sorted(glob.glob(path + cal_wav_ext))
        hdulist_cal_wav = fits.open(cal_wav_file[0])
        hdu_cal_wav = hdulist_cal_wav[0]
        datacal = hdu_cal_wav.data

    elif spect == "ESPRESSO" and len(pathlist) == 0:
        datafits = "/"+ccfext
        data_path = []
        data_path = np.append(data_path, glob.glob(path+datafits))
        data_path = sorted(data_path)
        print(data_path)
        datacal = None
    elif spect == "ESPRESSO" and len(pathlist) != 0:
        data_path = sorted(pathlist)
        print(data_path)
        datacal = None
    else:
        sys.exit("Spectrograph not recognized!")
    ldata = len(data_path)
    lorder = len(orders)
    rv = np.zeros((lorder, ldata))
    sigrv = np.zeros((lorder, ldata))
    phase = np.zeros((lorder, ldata))

    # Computes the wl ccfs from the individual orders
    if spect == "ESPRESSO":
        # Ver o valor do dataspec_path
        wl = np.matrix([readccf(data_path[k], None, [None for x in range(
            ldata)], st_type, [0, 170]) for k in range(ldata)]).T
        # plt.errorbar(np.array(wl[0])[0],np.array(wl[1])[0],yerr=np.array(wl[2])[0],fmt=".")
        # plt.show()
        # sys.exit()
        pool = multiprocessing.Pool()
        results = pool.map(getccf, ((order, data_path, None, [
                           None for x in range(ldata)], st_type) for order in orders))
        c = 0
    elif spect == "HARPS":
        wl = np.matrix([readccf(data_path[k], datacal, dataspec_path[k], st_type, [
                       0, 70]) for k in range(ldata)]).T
        pool = multiprocessing.Pool()
        results = pool.map(
            getccf, ((order, data_path, datacal, dataspec_path, st_type) for order in orders))
        c = 0
    for order in orders:
        rv[c] = np.array(results[c][1])[0]
        sigrv[c] = np.array(results[c][2])[0]
        phase[c] = np.array(results[c][0])[0]
        c += 1
    return(phase, rv, sigrv, np.array(wl[1])[0], np.array(wl[0])[0], np.array(wl[2])[0])


def RVerror(rv, ccf, eccf=1.):
    """
    Calculate the uncertainty on the radial velocity, following the same steps
    as the ESPRESSO DRS pipeline.

    Parameters
    ----------
    rv : array
    The velocity values where the CCF is defined.
    ccf : array
    The values of the CCF profile.
    eccf : array
    The errors on each value of the CCF profile.
    """

    ccf_slope = np.gradient(ccf, rv)
    ccf_sum = np.sum((ccf_slope / eccf)**2)
    return 1.0 / np.sqrt(ccf_sum)


def fgauss(x, a, mu, sigma2, c):
    return(c+a*np.exp(-(x-mu)**2./(2.*sigma2**2.)))


def getccf(k):
    order, data_path, datacal, dataspec_path, st_type = k
    ldata = len(data_path)
    f = [readccf(data_path[k], datacal, dataspec_path[k], st_type, order)
         for k in range(ldata)]
    f = np.matrix(f)
    f = f.T
    # plt.plot(f[1],f[2])
    # plt.show()
    return(f[0], f[1], f[2])


def readccf(filename, filecal, filespec, st_type, interval):
    ordi = orderl(interval)
    # CCF
    ccffile = fits.open(filename)
    header = ccffile[0].header
    # Leio RV do header:
    if spect == "ESPRESSO":
        rvheader = header["HIERARCH ESO QC CCF RV"]
        srvheader = header["HIERARCH ESO QC CCF RV ERROR"]
    else:
        None

    # Spectrum
    try:
        hdulistspec = fits.open(filespec)
    except:
        None
    if spect == "HARPS":
        hduspec = hdulistspec[0]
    elif spect == "ESPRESSO":
        None
    else:
        None
    try:
        data_spec = hduspec.data
    except:
        None

    # Time of observation MJD
    obs_time = header["MJD-OBS"]

    # Number of orbits since transit center epoch
    # tepoch=storbitpar["tepoch"]
    norb = (obs_time-storbitpar["tepoch"])/storbitpar["P"]

    # Number of full orbits
    nforb = round(norb)
    phase = norb-nforb
    if spect == "HARPS":
        dim = len(ccffile[0].data[:])
        ldim = len(ccffile[0].data[0])
    elif spect == "ESPRESSO":
        dim = len(ccffile[1].data[:])
        ldim = len(ccffile[1].data[0])
    else:
        None
    # Arrays for best fit radial velocity, error in the fit, order number
    rvp = np.zeros(dim)
    rv_err = np.zeros(dim)
    order = np.zeros(dim)
    ccfout = np.zeros(ldim)
    rvf = []
    ordsigrv = 0.
    ccferr = 0.
    if spect == "HARPS":
        ccc = 0
        for k in ordi:
            k = int(k)
            if int(sum(ccffile[0].data[k])) == 0:
                None
                #print("Warning: Some ccfs are not here!")
                # rvp[k],rv_err[k],order[k]=[0,0,0]
            else:
                ccf = ccffile[0].data[k]
                vstart, vstep = header['CRVAL1'], header['CDELT1']
                rv = np.linspace(vstart, vstart+vstep*(ccf.size-1), ccf.size)
                ccf = -ccf + abs(min(-ccf))*np.ones(len(ccf))
                ordsigrv += 1./RVerror(rv, ccf)**2.
                ccfout += ccf
                ccc += 1
    elif spect == "ESPRESSO":
        ccc = 0
        for k in ordi:
            k = int(k)
            if int(sum(ccffile[1].data[k])) == 0:
                None
                # Warns the user when some of the chosen orders are not present.
                print(("Warning: Some ccfs are not here!"+str(k)))
            else:
                ccf = ccffile[1].data[k]
                # print(k)
                # print(ccf)
                eccf = ccffile[2].data[k]
                vstart, vstep = header['HIERARCH ESO RV START'], header['HIERARCH ESO RV STEP']
                rv = np.linspace(vstart, vstart+vstep*(ccf.size-1), ccf.size)
                ccf = -ccf + abs(min(-ccf))*np.ones(len(ccf))
                """
				plt.plot(rv,ccf)
				plt.show(block=False)
				plt.pause(0.2)
				plt.close()
				"""
                ordsigrv += 1./RVerror(rv, ccf, eccf)**2.
                ccfout += ccf
                ccc += 1
    else:
        None
    #import sys
    # sys.exit("Done")
    rvout = rv
    sigrvccf = 1./np.sqrt(ordsigrv)
    mean = np.mean(rvout)
    sigma = np.sqrt(sum((rvout-mean)**2.)/(len(rvout)))
    guess = [max(ccfout), mean, sigma, 0.]
    fitpar, fitpar_cov = opt.curve_fit(fgauss, rvout, ccfout, guess)
    fitpar_err = np.sqrt(np.diag(fitpar_cov))
    """	
	print(fitpar,interval)
	print(phase,sigrvccf)
	plt.plot(rvout,ccfout,".")
	#print(fitpar)
	plt.plot(rvout,fgauss(rvout, fitpar[0], fitpar[1], fitpar[2],fitpar[3]))
	plt.show()
	"""
    return(phase, fitpar[1], sigrvccf, fitpar[2])

import numpy as np
from astropy.time import Time
from math import floor
import glob
from astropy.io import fits
from .bouchy_error import bouchy_err
import os
from .constants import storbitpar, pathlist, spect,ccfext
#from .constants import tepoch,P,r1, parguess, spect, pathlist
import subprocess
import scipy.optimize as op
from copy import deepcopy as dpcy
import matplotlib.pyplot as plt
import sys
import scipy.optimize as opt
import string
from .orders import orderl
#from joblib import Parallel, delayed
import multiprocessing
import emcee
from .priors import Gprior, Uprior


def get_data(path, st_type, orders):
    print(path)
    rvc = []
    rvc_err = []
    phase = []
    if spect == "HARPS":
        datafits = "/*ccf_"+str(st_type)+"_A.fits"
        dataspec = "/*e2ds_A.fits"

        data_path = sorted(glob.glob(path + datafits))
        dataspec_path = sorted(glob.glob(path + dataspec))
        dataspec_path = dpcy(dataspec_path[1:])

        # Calibration file
        cal_wav_ext = "/*wave_A.fits"
        cal_wav_file = sorted(glob.glob(path + cal_wav_ext))
        hdulist_cal_wav = fits.open(cal_wav_file[0])
        hdu_cal_wav = hdulist_cal_wav[0]
        datacal = hdu_cal_wav.data

    elif spect == "ESPRESSO" and len(pathlist) == 0:
        datafits = "/"+ccfext
        data_path = []
        data_path = np.append(data_path, glob.glob(path+datafits))
        data_path = sorted(data_path)
        print(data_path)
        datacal = None
    elif spect == "ESPRESSO" and len(pathlist) != 0:
        data_path = sorted(pathlist)
        print(data_path)
        datacal = None
    else:
        sys.exit("Spectrograph not recognized!")
    ldata = len(data_path)
    lorder = len(orders)
    rv = np.zeros((lorder, ldata))
    sigrv = np.zeros((lorder, ldata))
    phase = np.zeros((lorder, ldata))

    # Computes the wl ccfs from the individual orders
    if spect == "ESPRESSO":
        # Ver o valor do dataspec_path
        wl = np.matrix([readccf(data_path[k], None, [None for x in range(
            ldata)], st_type, [0, 170]) for k in range(ldata)]).T
        # plt.errorbar(np.array(wl[0])[0],np.array(wl[1])[0],yerr=np.array(wl[2])[0],fmt=".")
        # plt.show()
        # sys.exit()
        pool = multiprocessing.Pool()
        results = pool.map(getccf, ((order, data_path, None, [
                           None for x in range(ldata)], st_type) for order in orders))
        c = 0
    elif spect == "HARPS":
        wl = np.matrix([readccf(data_path[k], datacal, dataspec_path[k], st_type, [
                       0, 70]) for k in range(ldata)]).T
        pool = multiprocessing.Pool()
        results = pool.map(
            getccf, ((order, data_path, datacal, dataspec_path, st_type) for order in orders))
        c = 0
    for order in orders:
        rv[c] = np.array(results[c][1])[0]
        sigrv[c] = np.array(results[c][2])[0]
        phase[c] = np.array(results[c][0])[0]
        c += 1
    return(phase, rv, sigrv, np.array(wl[1])[0], np.array(wl[0])[0], np.array(wl[2])[0])


def RVerror(rv, ccf, eccf=1.):
    """
    Calculate the uncertainty on the radial velocity, following the same steps
    as the ESPRESSO DRS pipeline.

    Parameters
    ----------
    rv : array
    The velocity values where the CCF is defined.
    ccf : array
    The values of the CCF profile.
    eccf : array
    The errors on each value of the CCF profile.
    """

    ccf_slope = np.gradient(ccf, rv)
    ccf_sum = np.sum((ccf_slope / eccf)**2)
    return 1.0 / np.sqrt(ccf_sum)


def fgauss(x, a, mu, sigma2, c):
    return(c+a*np.exp(-(x-mu)**2./(2.*sigma2**2.)))


def getccf(k):
    order, data_path, datacal, dataspec_path, st_type = k
    ldata = len(data_path)
    f = [readccf(data_path[k], datacal, dataspec_path[k], st_type, order)
         for k in range(ldata)]
    f = np.matrix(f)
    f = f.T
    # plt.plot(f[1],f[2])
    # plt.show()
    return(f[0], f[1], f[2])


def readccf(filename, filecal, filespec, st_type, interval):
    ordi = orderl(interval)
    # CCF
    ccffile = fits.open(filename)
    header = ccffile[0].header
    # Leio RV do header:
    if spect == "ESPRESSO":
        rvheader = header["HIERARCH ESO QC CCF RV"]
        srvheader = header["HIERARCH ESO QC CCF RV ERROR"]
    else:
        None

    # Spectrum
    try:
        hdulistspec = fits.open(filespec)
    except:
        None
    if spect == "HARPS":
        hduspec = hdulistspec[0]
    elif spect == "ESPRESSO":
        None
    else:
        None
    try:
        data_spec = hduspec.data
    except:
        None

    # Time of observation MJD
    obs_time = header["MJD-OBS"]

    # Number of orbits since transit center epoch
    # tepoch=storbitpar["tepoch"]
    norb = (obs_time-storbitpar["tepoch"])/storbitpar["P"]

    # Number of full orbits
    nforb = round(norb)
    phase = norb-nforb
    if spect == "HARPS":
        dim = len(ccffile[0].data[:])
        ldim = len(ccffile[0].data[0])
    elif spect == "ESPRESSO":
        dim = len(ccffile[1].data[:])
        ldim = len(ccffile[1].data[0])
    else:
        None
    # Arrays for best fit radial velocity, error in the fit, order number
    rvp = np.zeros(dim)
    rv_err = np.zeros(dim)
    order = np.zeros(dim)
    ccfout = np.zeros(ldim)
    rvf = []
    ordsigrv = 0.
    ccferr = 0.
    if spect == "HARPS":
        ccc = 0
        for k in ordi:
            k = int(k)
            if int(sum(ccffile[0].data[k])) == 0:
                None
                #print("Warning: Some ccfs are not here!")
                # rvp[k],rv_err[k],order[k]=[0,0,0]
            else:
                ccf = ccffile[0].data[k]
                vstart, vstep = header['CRVAL1'], header['CDELT1']
                rv = np.linspace(vstart, vstart+vstep*(ccf.size-1), ccf.size)
                ccf = -ccf + abs(min(-ccf))*np.ones(len(ccf))
                ordsigrv += 1./RVerror(rv, ccf)**2.
                ccfout += ccf
                ccc += 1
    elif spect == "ESPRESSO":
        ccc = 0
        for k in ordi:
            k = int(k)
            if int(sum(ccffile[1].data[k])) == 0:
                None
                # Warns the user when some of the chosen orders are not present.
                print(("Warning: Some ccfs are not here!"+str(k)))
            else:
                ccf = ccffile[1].data[k]
                # print(k)
                # print(ccf)
                eccf = ccffile[2].data[k]
                vstart, vstep = header['HIERARCH ESO RV START'], header['HIERARCH ESO RV STEP']
                rv = np.linspace(vstart, vstart+vstep*(ccf.size-1), ccf.size)
                ccf = -ccf + abs(min(-ccf))*np.ones(len(ccf))
                """
				plt.plot(rv,ccf)
				plt.show(block=False)
				plt.pause(0.2)
				plt.close()
				"""
                ordsigrv += 1./RVerror(rv, ccf, eccf)**2.
                ccfout += ccf
                ccc += 1
    else:
        None
    #import sys
    # sys.exit("Done")
    rvout = rv
    sigrvccf = 1./np.sqrt(ordsigrv)
    mean = np.mean(rvout)
    sigma = np.sqrt(sum((rvout-mean)**2.)/(len(rvout)))
    guess = [max(ccfout), mean, sigma, 0.]
    fitpar, fitpar_cov = opt.curve_fit(fgauss, rvout, ccfout, guess)
    fitpar_err = np.sqrt(np.diag(fitpar_cov))
    """	
	print(fitpar,interval)
	print(phase,sigrvccf)
	plt.plot(rvout,ccfout,".")
	#print(fitpar)
	plt.plot(rvout,fgauss(rvout, fitpar[0], fitpar[1], fitpar[2],fitpar[3]))
	plt.show()
	"""
    return(phase, fitpar[1], sigrvccf, fitpar[2])
#from .constants import tepoch,P,r1, parguess, spect, pathlist #from joblib import Parallel, delayed # Calibration file # Computes the wl ccfs from the individual orders # Ver o valor do dataspec_path # plt.errorbar(np.array(wl[0])[0],np.array(wl[1])[0],yerr=np.array(wl[2])[0],fmt=".") # plt.show() # sys.exit() Calculate the uncertainty on the radial velocity, following the same steps     as the ESPRESSO DRS pipeline.      Parameters     ----------     rv : array     The velocity values where the CCF is defined.     ccf : array     The values of the CCF profile.     eccf : array     The errors on each value of the CCF profile. # plt.plot(f[1],f[2]) # plt.show() # CCF # Leio RV do header: # Spectrum # Time of observation MJD # Number of orbits since transit center epoch # tepoch=storbitpar["tepoch"] # Number of full orbits # Arrays for best fit radial velocity, error in the fit, order number #print("Warning: Some ccfs are not here!") # rvp[k],rv_err[k],order[k]=[0,0,0] # Warns the user when some of the chosen orders are not present. # print(k) # print(ccf) plt.plot(rv,ccf) 				plt.show(block=False) 				plt.pause(0.2) 				plt.close() #import sys # sys.exit("Done") print(fitpar,interval) 	print(phase,sigrvccf) 	plt.plot(rvout,ccfout,".") 	#print(fitpar) 	plt.plot(rvout,fgauss(rvout, fitpar[0], fitpar[1], fitpar[2],fitpar[3])) 	plt.show()