# -*- coding: utf-8 -*-
"""
Created on Mon Aug 9, 2021

@author: jbobowsk
"""

# This tutorial demonstrates how to work with and manipulate complex numbers.

# Import the usual modules for number crunching (numPy) and plotting (matplotlib)
import numpy as np
import matplotlib.pyplot as plt

# We will attempt to plot the magnitude and phase of the transfer function 
# for an op-amp based low-pass filter. 

# Here is the circuit that has been analyzed.  It is part of a larger circuit
# designed by Carl Michal for an Earth's-Field NMR apparatus: 
# https://iopscience.iop.org/article/10.1088/0957-0233/21/10/105902
import matplotlib.image as mpimg
img = mpimg.imread('low-pass.png')
plt.figure(figsize = (10, 7))
imgplot = plt.imshow(img)
plt.axis('off')

# The magnitude and phase of the transfer function v_out/v_in will be plotted
# as a function of frequency.  The frequency range will be 10 Hz to 100 kHz. 
# Set up a list of frequencies evenly spaced on a logarithmic scale.
ff = np.logspace(1, 5, 1000) # Hz
w = 2*np.pi*ff # rad/s

# Enter in component values and calculate the impedance of the capacitor.
# Notice the use of 1j to represent sqrt(-1).
R2 = 47e3 # ohms
R3 = 12e3 # ohms
C2 = 1e-9 # F
Zc2 = 1/(1j*w*C2) # ohms

# Here is the transfer function.  The voltage ratio v_out/v_in is complex.
vRatio = Zc2**2*(1 + R3/R2)/((R2 + Zc2)**2 - Zc2*R3)
print(vRatio)

# We can use abs(...) to find the magnitude of the complex numbers.  If z = a+j*b,
# abs(z) is equivalent to sqrt(a^2 + b^2).
print(abs(vRatio))

# We can use np.angle(...) to find the phase of the complex numbers.  If z = a+j*b,
# np.angle(z) is equivalent to arctan(b/a).
print(np.angle(vRatio)/np.pi)

# Here's a plot of the magnitude of the voltge ratio.
plt.figure()
plt.loglog(ff, abs(vRatio), 'k-', linewidth = 2.5)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Voltage Ratio')
plt.axis((10, 1e5, 1e-3, 2));

# Here's a plot of the phase in radians scaled by Pi.
plt.figure()
plt.semilogx(ff, np.angle(vRatio)/np.pi, 'r-', linewidth = 1.5)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Phase divided by Pi')
plt.axis((10, 1e5, -1, 0.05));