# -*- coding: utf-8 -*-
"""
Created on Thu Aug 13 14:48:19 2020

@author: Adriana Marquez
"""
#%%
import numpy as np
import matplotlib.pyplot as plt
import os
from IPython import get_ipython
from scipy.optimize import least_squares


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

#Analisis de datos a partir de varios archivos
#Poner la ruta de la carpeta donde se encuentran los archivos

os.chdir (r'C:\Users\User\Desktop\Adriana\Laboratorio 3\RLC\RLC_paralelo')

#Leo el archivo txt tiene que estar en la misma carpeta que el programa

file1 = 'Salidatransferencia20ohm.txt'
#en el archivo que leo las columnas están separadas por coma, la primera fila tiene el título entonces la salteo

Misdatos1 = np.loadtxt(file1, delimiter=",",skiprows=0) 
#print(Misdatos1)



"""
#leer los datos desde una archivo csv delimitado por tabs
print("nombre del archivo completo con terminación .csv incluida")
file = input()
Misdatos1=np.genfromtxt(file,delimiter='\t',skip_header=1)
"""
"""
#Datos ordenados col0 :transmitancia, col1: incertidumbre transmitancia
col2: diferencia de fase, col3: incertidumbre de diferencia de fase
col4: frecuencia, col5: incertidumbre de frecuencia
"""
frec= Misdatos1[:,0] #valores de frecuencia
tran= Misdatos1[:,3] #valores transferencia
errort = Misdatos1[:,4] #incertidumbre de transferencia
diffase= Misdatos1[:,5]
#errordf= Misdatos1[:,3]



        
#graficamos la transferencia en función de la frecuencia
plt.ion()
plt.close("all")
plt.figure(1)
plt.errorbar(frec, tran,yerr=errort,fmt=".b" );
plt.grid('on');
plt.xlabel('frecuencia(1/s)');
plt.ylabel('Transferencia');
plt.title('transferencia vs frecuencia')

#graficamos la diferencia de fase en función de la frecuencia
plt.figure(2)
plt.scatter(frec, diffase )
plt.grid('on')
plt.xlabel('frecuencia(1/s)')
plt.ylabel('Diferencia de fase')
plt.title('Diferencia de fase vs frecuencia')


#%%

#Ajuste para filtro pasa banda

def transferencia(parameters, x):
  f0 = parameters[0]
  Q = parameters[1]
  A= parameters[2]
  y=A*np.sqrt(1/(1+(1-(f0/x)**2)**2*(x*Q/f0)**2))
  return y

def get_residuals(parameters, y_data, x_data):
  theoretical_function = transferencia(parameters, x_data)
  residuals = np.abs(theoretical_function - y_data)
  return residuals

#Cargar los parámetros extraidos del gráfico
guess_f0 = 1250.0  #frecuencia de resonancia
guess_fs =  1840.0 # frecuencia superior correspondiente a Tmax/(2)^0.5
guess_fi = 875.0   # frecuencia inferior correspondiente a Tmax/(2)^0.5
guess_Tmax= 0.35  # valor de Tmax
guess_Q = guess_f0/(guess_fs-guess_fi)
guess_parameters = [guess_f0, guess_Q,guess_Tmax]

res_lsq = least_squares(get_residuals, guess_parameters, args=(tran,frec),method='trf')


best_parameters = res_lsq['x']

# 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,tran)


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




frecfit=np.arange(frec[0],frec[-1], (frec[-1]-frec[0])/100)
fitted_function = transferencia(best_parameters, frecfit)
plt.figure(3)
plt.scatter(frec, tran,label='data1')
plt.plot(frecfit, fitted_function, color = 'red', linewidth = 2.0, label='fit')
#N=len(frec)
#plt.xlim(x1[0],x1[N-1])
plt.xlabel("Frecuencia (1/s)")
plt.ylabel("Transferencia(V)")
fresonancia=int(best_parameters[0])
print(fresonancia)
plt.title('frecuencia de resonancia='+ str(fresonancia))


print("PARAMETROS DEL AJUSTE  de la transferencia")
print('Frecuencia de resonancia: (' ,round(best_parameters[0],3),'  ± ',round( pstd[0],3), ') 1/s')
print('Factor de mérito: ',round(best_parameters[1],3),'  ± ',round( pstd[1],3))
print('Máxima Transferencia: ',round(best_parameters[2],3),'  ± ',round( pstd[2],3))

Resistencia=20.0/best_parameters[2]
Inductancia= best_parameters[1]*Resistencia/(best_parameters[0]*2*np.pi)
Capacidad= 1/(best_parameters[0]*2*np.pi)**2/Inductancia

print (Resistencia, 'ohm')
print(Inductancia, 'H')
print(Capacidad, 'F')