# -*- coding: utf-8 -*-
"""
Created on Sat Oct 30 00:06:30 2021

@author: Joaquin de Jesus
"""

import numpy as np
import matplotlib.pyplot as plt
import os
from IPython import get_ipython
from scipy.optimize import curve_fit

#selecciono si las figuras aparecen en la terminal (inline) o en ventana emergente (qt5)
#get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('matplotlib', 'qt5')


#%% Cargo los datos experimentales

#Elijo el directorio donde se encuentran los archivos que voy a analizar
os.chdir (r'C:\Users\Desktop\PracticasLabo3\Filtros')

#Leo el archivo txt tiene que estar en la misma carpeta que el programa
file_name = "Salidatransferencia.txt"

Misdatos = np.loadtxt(file_name, delimiter=",",skiprows=0) 

# Ordeno los datos con la frecuencia:
Misdatos1 = Misdatos[np.argsort(Misdatos[:, 4])]

# Desempaquetamos los datos:
T = Misdatos1[:,0] # Transferencia
errorT = Misdatos1[:,1] # incertidumbre de T
fase = Misdatos1[:,2] # valores de dif fase en radianes
errorfase = Misdatos1[:,3] # valores de incertidumbre de la fase
fase_grad = 180*(Misdatos1[:,2])/(np.pi) # lee y convierte valores de dif fase en grados
errorfase_grad = 180*(Misdatos1[:,3])/(np.pi) # lee y convierte valores de incertidumbre de la fase
frec = Misdatos1[:,4] # v alores de la frecuencia (Hz)
errorfrec = Misdatos[:,5] # incertidumbre de la frecuencia
T_db = 20*np.log10(T) # Atenuación en dB
errorT_db= 20*(Misdatos1[:,1])/(Misdatos1[:,0])

#%% Grafico los datos experimentales para chequear que este todo bien


plt.close("all")
fs = "large" # tamaño de la leyenda con las labels ('small', 'medium', 'large', 'x-large', etc)
axesLabelSize = 15 # tamaño de las letras de los ejes
figSize = (10, 6) # tamaño de las figuras (ancho, alto)

## Transferencia
plt.figure(1, figsize=figSize)
plt.errorbar(frec, T, xerr=errorfrec, yerr=errorT, fmt=".b", linestyle='dotted')
plt.axhline(1, color='black', linestyle='dashed', linewidth=0.5)
plt.axhline(0, color='black', linestyle='dashed', linewidth=0.5)
plt.axhline(1/np.sqrt(2), color='green', linestyle='dashed', linewidth=1.5, label="1 / sqrt(2)")
#plt.grid('on');
#plt.axis([0,2,0,3.5]) #elijo el rango para el eje x del gráfico
plt.xscale('log')
plt.xlabel('Frecuencia (Hz)', size=axesLabelSize)
plt.ylabel('Transferencia', size=axesLabelSize)
#plt.legend(loc=3,prop={'size':12}, ncol=1, shadow=True, fancybox=True, title = "$\omega$ (rad/s)") #coloca las etiquetas en la mejor posición posible
plt.legend(fontsize=fs)

plt.figure(2, figsize=figSize)
plt.errorbar(frec, T_db, xerr=errorfrec, yerr=errorT_db,fmt=".b", linestyle='dotted')
plt.axhline(-3, color='green', linestyle='dashed', linewidth=1.5, label="-3 dB")
plt.axhline(0, color='black', linestyle='dashed', linewidth=0.5)
#plt.grid('on');
#plt.axis([0,2,0,3.5]) #elijo el rango para el eje x del gráfico
#plt.xscale('log')
plt.xlabel('Frecuencia (Hz)', size=axesLabelSize)
plt.ylabel('Atenuacion (dB)', size=axesLabelSize)
plt.legend(fontsize=fs)


plt.figure(3, figsize=figSize)
plt.errorbar(frec,fase_grad,xerr=errorfrec,yerr=errorfase_grad,fmt=".b")
plt.axhline(45, color='green', linestyle='dashed', linewidth=1.5, label="$\pm 45$ grados")
plt.axhline(0, color='black', linestyle='dashed', linewidth=0.5)
plt.axhline(90, color='grey', linestyle='dashed', linewidth=0.5)
plt.axhline(-45, color='green', linestyle='dashed', linewidth=1.5)
plt.axhline(-90, color='grey', linestyle='dashed', linewidth=0.5, label="$\pm 90$ grados")
#plt.grid('on');
#plt.axis([0,2,0,3.5]) #elijo el rango para el eje x del gráfico
plt.xscale('log')
plt.ylim(-360, 360)
plt.xlabel('Frecuencia (Hz)', size=axesLabelSize)
plt.ylabel('$\Delta$ Fase (grados)', size=axesLabelSize)
plt.legend(fontsize=fs)


plt.show()
#%% Defino funciones relevantes para el ajuste del filtro PB 

## TRANSFERENCIA
def Tr(x, f0):
    """
    x (float) - frecuencia (Hz)
    f0 (float) frecuencia de corte (Hz)
    """
    y = 1/(np.sqrt(1+np.square(x/np.abs(f0)))) #pasabajos
    return y

## ATENUACION
def A(x, f0):
    """
    x (float) - frecuencia (Hz)
    f0 (float) frecuencia de corte (Hz)
    """
    y = 20*np.log10(Tr(x, f0))
    return y

## DESFASAJE
def phi(x, f0):
    """
    x (float) - frecuencia (Hz)
    f0 (float) frecuencia de corte (Hz)
    """
    y = -(180.0 / np.pi) * np.arctan(x / f0)
    return y

#%% Ajuste

## Calculamos la frecuencia de corte esperada con los valores nominales de R y C.
## Esto es relevante pues va a servir como valor inicial de f0 para el ajuste

Rnominal = 2400 # (ohms)
Cnominal = 1 * 10**(-6) # (F)

w0_esperado = 1 / (float(Rnominal) * float(Cnominal))
f0_esperado = w0_esperado / (2*np.pi)
print("f0_esperado = %.3f Hz" % f0_esperado)

f0_guess = [f0_esperado] # usamos al f0_esperado como valor inicial para el ajuste

poptA, pcovA = curve_fit(f=A, xdata=frec, ydata=T_db, p0=f0_guess) ## Ajuste atenuacion
poptT, pcovT = curve_fit(f=Tr, xdata=frec, ydata=T, p0=f0_guess) ## Ajuste Transferencia
poptPhi, pcovPhi = curve_fit(f=phi, xdata=frec, ydata=fase_grad, p0=f0_guess) ## Ajuste fase

#### Desempaquetamos los parametros y sus errores obtenidos en el ajuste
f0_fitA = poptA[0]
f0_fitAerr = np.sqrt(np.diag(pcovA))
txtFitA = "fit\nf0 = (%.3f +- %.3f) Hz" % (f0_fitA, f0_fitAerr) # texto para el label del fit

f0_fitT = poptT[0]
f0_fitTerr = np.sqrt(np.diag(pcovT))
txtFitT = "fit\nf0 = (%.3f +- %.3f) Hz" % (f0_fitT, f0_fitTerr) # texto para el label del fit

f0_fitPhi = poptPhi[0]
f0_fitPhierr = np.sqrt(np.diag(pcovPhi))
txtFitPhi = "fit\nf0 = (%.3f +- %.3f) Hz" % (f0_fitPhi, f0_fitPhierr) # texto para el label del fit

################################
###### Graficamos datos y ajuste

### Definimos algunas cosas para que los graficos queden lindos y uniformes
fplot_data = np.logspace(np.log10(min(frec)*0.9), np.log10(max(frec)*1.1), 500) # vector para graficar mas lindo el ajuste. Equiespaciado logartimicamente
colorFit = "k" # color del ajuste
lwFit = 2 # linewidth del ajuste. El grosor de la linea del ajuste en el grafico
alphaDatos = 0.6 # nivel de opacidad de los markers de los datos (0 es totalmente transparente, 1 es totalmente opaco)
colorDatos = 'm' # color de los markers de los datos
markerDatos = '^' # tipo de marker para los datos ('^' triangulitos, 'o' circulitos, 's', cuadraditos, '*' estrellitas)
fs = "large" # tamaño de la leyenda con las labels ('small', 'medium', 'large', 'x-large', etc)
axesLabelSize = 15 # tamaño de las letras de los ejes
figSize = (10, 6) # tamaño de las figuras (ancho, alto)
title = "Filtro PB RC\ndatos: %s" % file_name
#title = "Filtro PB RC"

## Ahora si, graficamos
plt.close("all")

## Atenuacion
plt.figure(1, figsize=figSize)
plt.errorbar(frec, T_db, yerr=errorT_db, fmt=markerDatos, color=colorDatos, label='datos', alpha=alphaDatos)
plt.plot(fplot_data, A(fplot_data, *poptA), label=txtFitA, color=colorFit, lw=lwFit)
plt.xlabel("frecuencia (Hz)", fontsize=axesLabelSize)
plt.ylabel("Atenuacion (dB)", fontsize=axesLabelSize)
plt.xscale("log")
plt.title(title)
plt.grid()
plt.legend(fontsize=fs)

## Transferencia
plt.figure(2, figsize=figSize)
plt.errorbar(frec, T, yerr=errorT, fmt=markerDatos, color=colorDatos ,label='datos', alpha=alphaDatos)
plt.plot(fplot_data, Tr(fplot_data, *poptT), label=txtFitT, color=colorFit, lw=lwFit)
plt.xlabel("frecuencia (Hz)", fontsize=axesLabelSize)
plt.ylabel("Transferencia", fontsize=axesLabelSize)
plt.xscale("log")
plt.title(title)
plt.grid()
plt.legend(fontsize=fs)

## Fase
plt.figure(3, figsize=figSize)
plt.errorbar(frec, fase_grad, yerr=errorfase_grad, fmt=markerDatos, color=colorDatos, label='datos', alpha=alphaDatos)
plt.plot(fplot_data, phi(fplot_data, *poptPhi), label=txtFitPhi, color=colorFit, lw=lwFit)
plt.xlabel("frecuencia (Hz)", fontsize=axesLabelSize)
plt.ylabel("dif. de fase (grados)", fontsize=axesLabelSize)
plt.xscale("log")
plt.title(title)
plt.grid()
plt.legend(fontsize=fs)