# -*- coding: utf-8 -*-
"""
Created on Wed Jul 22 18:49:57 2020

@author: User
"""

#%%
#Loading modules
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
import os
from IPython import get_ipython

#%%


#selecciono el grafico en Terminal (inline) o en ventana emergente (qt5)
#get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('matplotlib', 'qt5')

#Elijo el directorio donde se encuentran los archivos que voy a analiza
os.chdir (r'G:\Mi unidad\Labo 3 1er cuat 2022\Clase 04 19\Mediciones')


file = '1000.txt'

#Importing our data
data = np.loadtxt(file,dtype=float,delimiter = ' ',skiprows= 1)



x1i = data[:,0]#tiempo
y1i = data[:,1]#señal de tensión  Vin
y2i = data[:,2]#señal de tensión  Vc


#Graficando nuestros datos en función del número de dato
plt.ion()
plt.close("all")
plt.figure(1)
plt.plot(x1i, y1i,'.-r')
plt.plot(x1i, y2i,'.-g');
plt.grid('on');
plt.xlabel('tiempo (s)');
plt.ylabel('Tensión(V)');

#%%
#Seleccionar el ciclo ON
#limitar los datos al rango de interés
j=5
datoinicial1=0
datofinal1=0
dif=0
n=5
while dif<0.5:   
    j=j+1
    dif=y1i[j]-y1i[j-n]
else:
    datoinicial1=j-n #dato inicial del ciclo ON

dif=y1i[j-n]-y1i[j]    

while dif<0.5:   
    j=j+1
    dif=y1i[j-n]-y1i[j]
else:
    datofinal1=j-2*n #dato final del ciclo ON

datoinicial2=j #inicio de ciclo OFF
datofinal2=datoinicial2+datofinal1-datoinicial1 #fin de ciclo OFF





x1=x1i[datoinicial1:datofinal1]-x1i[datoinicial1]
y1=y1i[datoinicial1:datofinal1]
y2=y2i[datoinicial1:datofinal1]




#Graficamos los datos en la zona elegida
plt.close()
plt.figure(2)
plt.plot(x1, y2,'.-g');
plt.grid('on');
plt.xlabel('tiempo(s)');
plt.ylabel('Vcapacitor(V)');

print("\033[4;35m"+"OBSERVAR LAS FIGURAS Y CHEQUEAR QUE SON LAS ESPERADAS") 


#%%
#definimos la función que vamos a usar para el ajuste, una función senoidal amortigauda + una constante

def cos_fit_fun_damped(parameters,time):
  a = parameters[0]
  omega = parameters[1]
  offset = parameters[2]
  phi = parameters[3]
  alpha = parameters[4]
  y = a * np.cos(omega * time + phi) * np.exp(-alpha*time) + offset
  return y

def get_residuals(parameters, position_data, time_data):
  theoretical_function = cos_fit_fun_damped(parameters,time_data )
  residuals = np.abs(theoretical_function - position_data)
  return residuals

#Buscar los parámetros característicos de la señal Vc
#Busca los dos primeros máximos
j1max=np.argmax(y2)
y22=y2[j1max:-1]
j1min=np.argmin(y22)
y23=y2[j1max+j1min:-1]
j2max=np.argmax(y23)+j1max+j1min

T1=x1[j1max]
T2=x1[j2max]


Amp1=y2[j1max]
Amp2=y2[j2max] 
vm = y2[-10]


#Guess parameters
guess_amplitude = Amp1-vm #se estima la amplitud
guess_omega = 2*np.pi/(T2-T1) #a partir de los datos se estima la frecuencia angular
guess_offset = vm
guess_phi = np.pi 
guess_alpha = -np.log((Amp2-vm)/(Amp1-vm))/(T2-T1)#se estima la fase inicial
guess_parameters = [guess_amplitude, guess_omega, guess_offset, guess_phi, guess_alpha]

print(guess_parameters)

#Performing the fit
res_lsq = least_squares(get_residuals, guess_parameters, args=(y2,x1))

#Fit results
best_parameters = res_lsq['x']

N=len(x1)
x1f=np.arange(x1[0],x1[N-1],step=(x1[N-1]-x1[0])/1000)
fitted_function = cos_fit_fun_damped(best_parameters, x1f)



# Calculamos la matriz de covarianza "pcov"
def calcular_cov(res,y_datos):
    U, S, V = np.linalg.svd(res.jac, full_matrices=False)
    threshold = np.finfo(float).eps * max(res.jac.shape) * S[0]
    S = S[S > threshold]
    V = V[:S.size]
    pcov = np.dot(V.T / S**2, V)

    s_sq = 2 * res.cost / (y_datos.size - res.x.size)
    pcov = pcov * s_sq
    return pcov

pcov = calcular_cov(res_lsq,y2)

# De la matriz de covarinza podemos obtener los valores de desviación estándar
# de los parametros hallados
pstd =np.sqrt(np.diag(pcov))

print('Parámetros hallados (con incertezas):')
   
print('Best Amplitude: ',round(best_parameters[0],3),'  ± ',round( pstd[0],3))
print('Best Omega: ' ,round(best_parameters[1],3),'  ± ',round( pstd[1],3))
print('Best offset: ' ,round(best_parameters[2],3),'  ± ',round( pstd[2],3))
print('Best Phi: ' ,round(best_parameters[3],3),'  ± ',round( pstd[3],3))
print('Best Alpha: ' ,round(best_parameters[4],3),'  ± ',round( pstd[4],3))

#Plotting results
plt.figure(3)
plt.scatter(x1, y2,  label='data')
plt.plot(x1f, fitted_function, color = 'red', linewidth = 2.0, label='fit')
plt.xlim(x1[0],x1[N-1])
plt.xlabel("Time (s)")
plt.ylabel("Vcapacitor(V)")
plt.title("RLC Oscillation")
plt.legend()
plt.show()

C=10**(-7) #capacidad en faradios
W02=(best_parameters[1]**2+best_parameters[4]**2)
L=1/(C*W02)
deltaW02=2*best_parameters[1]*pstd[1]+2*best_parameters[4]*pstd[4]
deltaL=(deltaW02/W02)*L
R=best_parameters[4]*2*L
deltaR=(pstd[4]/best_parameters[4]+deltaL/L)*R

print("resistencia:(", round(R,2),'  ± ',round( deltaR,2),')ohm')

print("inductancia: (", round(L,5),'  ± ',round( deltaL,5),')H')