###############################################################################
# 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
from scipy.signal import find_peaks_cwt
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

###############################################################################
# ESCRIBIMOS LA SECUENCIA DE PROCESOS DE NUESTRO PROGRAMA
###############################################################################

nombre_archivo='C:\\Users\\Usuario\\Dropbox\\Clase Estadística\\datos_osilaciones.xlsx'
columna=0
t=abrir_columna_excel(nombre_archivo,columna)
columna=1
x=abrir_columna_excel(nombre_archivo,columna)

t=np.asarray(t)
x=np.asarray(x)

plt.plot(t,x)

nventana = 20
indices_picos = find_peaks_cwt(x,np.arange(30,100)) # Existen otros metodos, por ejemplo en https://docs.scipy.org/doc/scipy/reference/signal.html


plt.plot(t[indices_picos], x[indices_picos],'o')

#sacamos el primer y ultimo pico
indices_picos=indices_picos[1:-1]
plt.plot(t[indices_picos], x[indices_picos],'o')
plt.xlabel('Tiempo')
plt.ylabel('Posición')


tiempos_picos = t[indices_picos]
dt_picos = tiempos_picos[1:len(tiempos_picos)]-tiempos_picos[0:(len(tiempos_picos)-1)]
w = 2*np.pi/dt_picos
w_media = np.mean(w)
w_error = np.std(w,ddof=1)/np.sqrt(len(w))
print('Frecuencia angular = '+str(w_media)+' ± '+str(w_error))


amplitudes_picos = x[indices_picos] - np.mean(x)

amplitudes_picos_error = 0.1
tiempos_picos_error = 0.1

plt.figure()
plt.errorbar(tiempos_picos,amplitudes_picos,xerr=tiempos_picos_error,yerr=amplitudes_picos_error,color='b',fmt='o')#,10000,fmt='-o',color='g')
plt.xlabel('Tiempo')
plt.ylabel('Amplitud')

# Linealizo
# Modelo: x ~ A*exp(-b*t)
# ln(x) ~ ln(A*exp(-b*t)) = ln(A) + ln(exp(-b*t)) = ln(A) - b*t
# ln(x) ~ ln(A) - b*t
log_amplitudes_picos = np.log(amplitudes_picos)
log_amplitudes_picos_error = amplitudes_picos_error/amplitudes_picos
plt.figure()
plt.errorbar(tiempos_picos,log_amplitudes_picos,xerr=tiempos_picos_error,yerr=log_amplitudes_picos_error,color='b',fmt='o')#,10000,fmt='-o',color='g')
plt.xlabel('Tiempo')
plt.ylabel('log(Amplitud)')