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

@author: Adriana Márquez

Estudio de la corriente de circuito RLC sobreamotriguado
"""

#%%
#Loading modules
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
from scipy.optimize import curve_fit
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')

print("nombre del archivo completo con terminación .txt incluida")
file = input()

#Completar con los parámetros del circuito
print("Cargar valor de la Resistencia del cirucito en ohm")
R= float(input())
print("Cargar valor de la Inductancia del cirucito en H")
L= float(input())
#print("Cargar valor de la Capacidad del cirucito en F")
#C= float(input())

#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 todos nuestros datos en función del tiempo

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)')

#%%
#limitar los datos al rango de interés
j=5
datoinicial1=0
datofinal1=0
dif=0
n=5
tolerancia=1.0
while dif<tolerancia:   
    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<tolerancia:   
    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 en función del tiempo
plt.ion()
fig, ax1 = plt.subplots()
ax1.plot(x1, y1,'.-g')
ax1.set_ylabel('Tensión entrada(V)')
ax1.set_xlabel('tiempo s')
ax2=ax1.twinx()
ax2.plot(x1, y2,'.-r')
ax2.set_ylabel('Caída de tensión en C (V)', color='r')




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

#Modelo aproximado para  VC
Vfinal=y2[-10]
def aproximacionVC(time, a,b):
    y= a*np.exp(-b*time)+Vfinal
    return y
"""
#Modelo aproximado para la corriente
def aproximacionI(time, a,b):
    y= a*np.exp(-b*time)
    return y
"""
# Ajuste de los datos
popt1, pcov01 = curve_fit(aproximacionVC, x1,y2)
#popt2, pcov02 = curve_fit(aproximacionI, x1,y2)

funcredVC=aproximacionVC(x1,*popt1)
#funcredI=aproximacionI(x1,*popt2)
# ESTIMO LOS ERRORES ( 1 sigma )
perr1 = np.sqrt(np.diag(pcov01))
#perr2 = np.sqrt(np.diag(pcov02))
N=len(x1)

"""
#modelo completo para la corriente
def modeloI(parameters1,time):
  A2= parameters1[0]
  landa12 = parameters1[1]
  landa22= parameters1[2]
  y = A2*(np.exp(landa12*time)-np.exp(landa22*time))
  return y

#Guess parameters
guess_landa12 = -popt2[1]
guess_landa22 = (-R/L- popt2[1]) #taken from theory
guess_A2 = popt2[0]

guess_parameters2 = [guess_A2, guess_landa12, guess_landa22]
"""



#modelo completo para la caida sobre el capacitor
def modeloVC(parameters1,time):
  A1= parameters1[0]
  landa11 = parameters1[1]
  landa21= parameters1[2]
  y = A1*(landa21*np.exp(landa11*time)-landa11*np.exp(landa21*time))+Vfinal
  return y

#Guess parameters

guess_landa11 = -popt1[1]
guess_landa21 = (-R/L- popt1[1]) #tomados de la teoría
guess_A1 = popt1[0]/(guess_landa21)


guess_parameters1 = [guess_A1, guess_landa11, guess_landa21]





def get_residuals1(parameters1, y_data, x_data):
  theoretical_function = modeloVC(parameters1,x_data)
  residuals1 = np.abs(theoretical_function - y_data)
  return residuals1
  
"""
def get_residuals2(parameters2, y_data, x_data):
  theoretical_function = modeloI(parameters2,x_data)
  residuals2 = np.abs(theoretical_function - y_data)
  return residuals2
"""

#Performing the fit
res_lsq1 = least_squares(get_residuals1, guess_parameters1, args=(y2,x1))
#res_lsq2 = least_squares(get_residuals2, guess_parameters2, args=(y2,x1))
#Fit results
best_parameters1 = res_lsq1['x']
#best_parameters2 = res_lsq2['x']

x1f=np.arange(x1[0],x1[N-1],step=(x1[N-1]-x1[0])/500)
fitted_functionVC = modeloVC(best_parameters1, x1f)
#fitted_functionI = modeloI(best_parameters2, 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

pcov1 = calcular_cov(res_lsq1,y2)
#pcov2 = calcular_cov(res_lsq2,y2)

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

print()
print('Modelo aproximado para el Vcpactior =a*np.exp(-b*time)+Vfinal')#
print('a = (', round(popt1[0],6), '+/-',  round(perr1[0],6), ')V')
print('b = (', round(popt1[1],3), '+/-',  round(perr1[1],3), ') 1/s')
"""
print('Modelo aproximado para la corriente =a*np.exp(-b*time)')#
print('a = (', round(popt2[0],6), '+/-',  round(perr2[0],6), ')V')
print('b = (', round(popt2[1],3), '+/-',  round(perr2[1],3), ') 1/s')
#best=popt1[1]
#res1=L/best*(1/(L*C)-best**2)
#print('valor aproximado de resistencia (ohm): ', int(res1))

print()

print('parametros iniciales para el ajuste de la corriente')
print(guess_parameters2)
print()
print('Parámetros hallados con el modelo completo de la corriente (con incertezas):')
   
print('Best A2: ',(best_parameters2[0]),'  ± ',( pstd2[0]))
print('Best landa1: ' ,best_parameters2[1],'  ± ', pstd2[1])
print('Best landa2: ' ,best_parameters2[2],'  ± ',( pstd2[2]))
#print('Best V0: ' ,round(best_parameters[3],3),'  ± ',round( pstd[3],3))
"""
print()

print('parametros iniciales para el ajuste de la V capacitor')
print(guess_parameters1)
print()
print('Parámetros hallados con el modelo completo de la Vcapacitor (con incertezas):')
   
print('Best A1: ',best_parameters1[0],'  ± ', pstd1[0])
print('Best landa1: ' ,(best_parameters1[1]),'  ± ',( pstd1[1]))
print('Best landa2: ' ,(best_parameters1[2]),'  ± ',( pstd1[2]))

fitted_functionVC=modeloVC(best_parameters1,x1)
#fitted_functionI=modeloI(best_parameters2,x1)

#Plotting results
plt.figure(3)
plt.scatter(x1, y2,  label='data')
plt.plot(x1,  fitted_functionVC, color = 'red', linewidth = 2.0, label='modelo completo')
plt.plot(x1, funcredVC, color = 'green', linewidth = 2.0, label='modelo aproximado')
plt.xlim(x1[0],x1[-1])
plt.xlabel("Time (s)")
plt.ylabel("V Capacitor (V)")
plt.title("RLC Sobreamortiguado")
plt.legend()
plt.show()

"""
plt.figure(4)
plt.scatter(x1, y2,  label='data')
plt.plot(x1, fitted_functionI, color = 'red', linewidth = 2.0, label='modelo completo')
plt.plot(x1, funcredI, color = 'green', linewidth = 2.0, label='modelo aproximado')
plt.xlim(x1[0],x1[-1])
plt.ylim(y2[-1],y2[0])
plt.xlabel("Time (s)")
plt.ylabel("Corriente (A)")
plt.title("RLC Sobreamortiguado")
plt.legend()
plt.show()


print()
print('A partir del ajuste de la corriente')
ResInduci=-(best_parameters2[1]+best_parameters2[2])
Induci= Vfinal/best_parameters2[0]/(best_parameters2[1]-best_parameters2[2])
Capi=1/Induci*1/(((best_parameters2[1]+best_parameters2[2])/2)**2-((best_parameters2[1]-best_parameters2[2])/2)**2)
Resi=ResInduci*Induci

print("resistencia(ohm):", round(Resi,3))
print("Capacidad(C):", (Capi))
print ('Inductancia (H)',Induci)
"""
print()
 
print('A partir del ajuste de la Vcapacitor')
ResInduc=-(best_parameters1[1]+best_parameters1[2])
Induc= Vfinal/best_parameters1[0]/(best_parameters1[1]-best_parameters1[2])
Cap=1/Induc*1/(((best_parameters1[1]+best_parameters1[2])/2)**2-((best_parameters1[1]-best_parameters1[2])/2)**2)
Res=ResInduc*Induc

print("resistencia(ohm):", round(Res,3))
print("Capacidad(C):", (Cap))
print ('Inductancia (H)',Induc)