###############################################################################
# CARGAMOS LAS LIBRERIAS DE FUNCIONES HECHAS QUE VAMOS A USAR
###############################################################################
import numpy as np #PARA HACER CUENTAS
import xlrd #PARA CARGAR ARCHIVOS DE EXCEL
import matplotlib.pyplot as plt #PARA GRAFICAR
import scipy.optimize as scop #PARA ANALISIS DE DATOS
plt.close('all') #CIERRA CUALQUIER GRAFICO QUE HAYA ABIERTO

###############################################################################
# DEFINIMOS NOSOTROS OTRAS DE LAS FUNCIONES QUE VAYAMOS A USAR
###############################################################################
def abrir_columna_excel(nombre_archivo,columna):
  wb = xlrd.open_workbook(nombre_archivo) 
  Hoja1 = wb.sheet_by_index(0) 
  N_datos=len(Hoja1.col_values(columna))
  datos=[]
  for i in range(1,N_datos):
    if type(Hoja1.cell_value(i, columna))==float:
      datos.append(Hoja1.cell_value(i, columna))
  return datos

def residuo(parametros,x,y_medida):
    y_teorica=gaussiana(x,parametros)
    diferencia=y_medida-y_teorica
    return diferencia
  
def gaussiana(x,parametros):
    mu,sigma,A=parametros
    y_teorica=A*np.exp(-(x-mu)**2/(2*sigma**2))
    return y_teorica

def graficar_ajuste_intervalos(periodos_obturacion,parametros_ajuste):
  mu=parametros_ajuste[0][0]
  sigma=parametros_ajuste[0][1]
  x_teorico=np.linspace(mu-4*sigma,mu+4*sigma,1000)
  y_teorica=gaussiana(x_teorico,parametros_ajuste[0])
  # GRAFICO EL AJUSTE Y LOS INTERVALOS
  plt.plot(x_teorico,y_teorica,label='ajuste')
  plt.axvspan(mu-sigma,mu+sigma, alpha=0.25, color='orange',label='1 sigma (68.2 %)')
  plt.axvspan(mu-2*sigma,mu+2*sigma, alpha=0.125, color='orange',label='2 sigma (95.4 %)')
  plt.axvspan(mu-3*sigma,mu+3*sigma, alpha=0.0625, color='orange',label='3 sigma (99.6 %)')
  plt.axvspan(mu-4*sigma,mu+4*sigma, alpha=0.03125, color='orange',label='4 sigma (99.8 %)')
  plt.legend()
###############################################################################
# ESCRIBIMOS LA SECUENCIA DE PROCESOS DE NUESTRO PROGRAMA
###############################################################################
# CARGAMOS LAS COLUMNAS DE TIEMPOS Y VOLTAJES
nombre_archivo = ("C:\\Users\\Usuario\\Dropbox\\Clase Estadística\\datos_carrito.xlsx")
columna=0
tiempos=abrir_columna_excel(nombre_archivo,columna)
columna=1
voltajes=abrir_columna_excel(nombre_archivo,columna)

plt.figure()
plt.plot(tiempos,voltajes)
plt.xlabel('Tiempo (s)');plt.ylabel('Voltaje')

# Los paso de listas a arrays
tiempos=np.asarray(tiempos)
voltajes=np.asarray(voltajes)

# Los voltajes que van de 0 a 5 y tienen ruido los hago binarios (True o False)
obturador=voltajes>2.5

# Obtengo los tiempos de obturacion o desobturacion
tiempos_de_obturacion=[]
tiempos_de_desobturacion=[]
for i in range(len(obturador)-1):
  if obturador[i]==True and obturador[i+1]==False:
    tiempos_de_obturacion.append(tiempos[i+1])
  if obturador[i]==False and obturador[i+1]==True:
    tiempos_de_desobturacion.append(tiempos[i+1])

# Verifico graficamente que se hayan calculado bien
plt.figure()
plt.plot(tiempos,obturador)
plt.plot(tiempos_de_obturacion,np.ones(len(tiempos_de_obturacion)),'o',label='Obturacion')
plt.plot(tiempos_de_desobturacion,np.ones(len(tiempos_de_desobturacion)),'o',label='Desobturacion')
plt.xlabel('Tiempo (s)');plt.ylabel('Señal Binaria');plt.legend()

# Calculo los periodos haciendo la diferencia de tiempos de obturacion
periodos_obturacion=np.diff(tiempos_de_obturacion)
periodos_desobturacion=np.diff(tiempos_de_desobturacion)

# Hago un histograma y ajusto con una gaussiana para obtener la media y el sigma
plt.figure()
histo=plt.hist(periodos_obturacion)

plt.xlabel('Periodo(s)');plt.ylabel('Frecuencia')
# DEFINO X e Y PARA AJUSTAR
x=[]
for i in range(len(histo[1])-1):
  x.append((histo[1][i+1]+histo[1][i])/2)
x=np.asarray(x)
y=histo[0]

# AJUSTO
parametros_iniciales=[np.mean(periodos_obturacion),np.std(periodos_obturacion),1]  
parametros_ajuste=scop.leastsq(residuo, parametros_iniciales, args=(x,y), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)
mu=parametros_ajuste[0][0]
sigma=parametros_ajuste[0][1]

graficar_ajuste_intervalos(periodos_obturacion,parametros_ajuste)

cant_cifras=7
print('El periodo es ( '+str(mu)[:cant_cifras]+' ± '+str(sigma)[:cant_cifras]+' )s')