```python
# Libraries Sympy and Numpy
from sympy import*
import numpy as np

# To define automatic printing mode (Not needed anymore)
# init_printing()

# For priting with text
from IPython.display import display, Latex

# Plotting Libraries
from mpl_toolkits import mplot3d

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.ticker as ticker

#interactive plotting in separate window
#%matplotlib qt 

#normal charts inside notebooks
#%matplotlib inline 

# Use Latex in plots
plt.rcParams['text.usetex'] = False
```

# 1. IPT System Frequency Splitting Model

## Equations Derivation


```python
# Definition of the symbolic variables

# Resistances --> Coils and Load
R_1, R_2, R_L = symbols('R_1, R_2, R_L',
                       real = True)

# Inductance and Capacitance
C_1, C_2, L_1, L_2, M_i = symbols('C_1, C_2, L_1,L_2, M_i', 
                                 real = True)

# Reactance
X_1, X_2 = symbols('X_1, X_2',
                       real = True)

# Angular frequency
w_e = symbols('w_e', real = True)

# Input Voltage and Current
V_1, I_1 = symbols('V_1, I_1', real = True)
```


```python
# Reactance in each mesh
#X_1 = (w_e*L_1) - (1/(w_e*C_1))
#X_2 = (w_e*L_2) - (1/(w_e*C_2))
```

- Equivalent Circuit's Matrix Representation: $[Z][I] = [V]$


```python
# Impedance Matriz [Z]
Z_m = Matrix([[R_1 + (I*X_1), I*w_e*M_i],
            [I*w_e*M_i, R_2 + R_L + (I*X_2)]])

# Voltage Vector
V_m = Matrix([[V_1], [0]])
```


```python
# Solve the System of Equations
I_m = Z_m.inv()*V_m
```


```python
# Input Current
result_I = "Input Current: \n $${} = {}$$".format(latex(I_1), 
                                                 latex(I_m[0]))
display(Latex(result_I))

# Separate Numerator and denominator of I1
I_1_n, I_1_d = fraction(I_m[0])
```


Input Current: 
 $$I_{1} = \frac{V_{1} \left(- i R_{2} - i R_{L} + X_{2}\right)}{- i M_{i}^{2} w_{e}^{2} - i R_{1} R_{2} - i R_{1} R_{L} + R_{1} X_{2} + R_{2} X_{1} + R_{L} X_{1} + i X_{1} X_{2}}$$



```python
# Operate separately each polynomial
display(factor(I_1_n))
display(factor(I_1_d + (I*(M_i**2)*(w_e**2))))
```


$\displaystyle - i V_{1} \left(R_{2} + R_{L} + i X_{2}\right)$



$\displaystyle - i \left(R_{1} + i X_{1}\right) \left(R_{2} + R_{L} + i X_{2}\right)$



```python
# Reconstruct the factored result
I_1 = (factor(I_1_n))/(factor(I_1_d + (I*(M_i**2)*(w_e**2))) 
                       - (I*(M_i**2)*(w_e**2)))
display(simplify(I_1))
```


$\displaystyle \frac{V_{1} \left(R_{2} + R_{L} + i X_{2}\right)}{M_{i}^{2} w_{e}^{2} + \left(R_{1} + i X_{1}\right) \left(R_{2} + R_{L} + i X_{2}\right)}$


- Input Impedance Expression $Z_{in}$:


```python
Z_in = V_1*(1/simplify(I_1))
display(Z_in)
```


$\displaystyle \frac{M_{i}^{2} w_{e}^{2} + \left(R_{1} + i X_{1}\right) \left(R_{2} + R_{L} + i X_{2}\right)}{R_{2} + R_{L} + i X_{2}}$


- Freq. Splitting occurs when $\Im(Z_{in}) = 0$


```python
# Assume both reactance are equal for simplicity
# Common reactance X_e
X_e = symbols('X_e', real = True)

# Substitute the variable
Z_in = Z_in.subs([(X_1, X_e), (X_2, X_e)])

# Imaginary part of Zin
simplify(im(Z_in))
```




$\displaystyle \frac{X_{e} \left(- M_{i}^{2} w_{e}^{2} - R_{1} R_{2} - R_{1} R_{L} + X_{e}^{2} + \left(R_{2} + R_{L}\right) \left(R_{1} + R_{2} + R_{L}\right)\right)}{X_{e}^{2} + \left(R_{2} + R_{L}\right)^{2}}$




```python
# Operate separately
factor(-(R_1*R_2) - (R_1*R_L) + ((R_2 + R_L)*(R_1 + R_2 + R_L)))
```




$\displaystyle \left(R_{2} + R_{L}\right)^{2}$



## a) 1st Condition to only have one non-trivial root: $\Delta < 0$


```python
# Angular Resonance frequency (central)
w_c = 2*np.pi*(800e3)  # 800 kHz

# Assume coils' self-inductance equal [H]
L_c = 53.49e-6

# Frequency-dependent parameter
beta_c = (w_c*L_c)**2
```


```python
# Formula for condition Delta
Delta_c = lambda R, beta, k: (((R**2) - (2*beta))**2) - (4*(beta**2)*(1 - (k**2)))

# Inductive Coupling Coefficient Range k_i
k_i = np.arange(0, 1 + 0.01, 0.01)
print(k_i.shape)
print(k_i[-1])

# Load Resistance [ohms]
R_L = np.arange(1., 1000. + 1., 1.)
print(R_L.shape)
print(R_L[-1])

# Create a mesh
k_i_m, R_L_m = np.meshgrid(k_i, R_L)

# Calculate the function Delta
f_delta = Delta_c(R_L_m, beta_c, k_i_m)
print(f_delta.shape)
```

    (101,)
    1.0
    (1000,)
    1000.0
    (1000, 101)



```python
#interactive plotting in separate window
#%matplotlib qt 

#normal charts inside notebooks
%matplotlib inline 

# Set the plot --> size(width, height)
fig1 = plt.figure(figsize=(15,10))

################
# 1st Subplot
################

# Set the subplot for the two graphs
axs = fig1.add_subplot(1, 2, 1, projection = '3d')

# Complete figure information
fig1.suptitle('Frequency Splitting: $\Delta$ < 0', fontsize = 20)

# For a simple 3D plot
#ax = fig.gca(projection = '3d')

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R^{\'}_{L}$) [$\Omega$]')
axs.set_zlabel('Function $\Delta$')

# Title
axs.set_title('3D Plot of $\Delta$ Function', fontsize = 15)

# Plot the surface
surf = axs.plot_surface(k_i_m, R_L_m, f_delta, cmap = cm.coolwarm,
                      linewidth = 0, antialiased = False)

# Figure size (single plot)
#plt.rcParams["figure.figsize"] = (30,30)

# Add a color bar which maps values to colors.
plt.colorbar(surf, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.05)

################
# 2nd Subplot
################

# Set the subplot for the two graphs
axs = fig1.add_subplot(1, 2, 2)

# Plot the contour plot
cp = plt.contour(k_i_m, R_L_m, f_delta, cmap = cm.Spectral)

# Define the levels to label
#c_labels = [0.e0]
# To specify scientific notation
fmt = ticker.LogFormatterSciNotation()

axs.clabel(cp, inline = 1, colors = 'k', fmt = fmt, fontsize = 10)

# Color bar
plt.colorbar(cp, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.1)

# Information
axs.set_title('Contour Plot of $\Delta$ Function', fontsize = 15)

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R^{\'}_{L}$) [$\Omega$]')

axs.xaxis.grid(True)
axs.yaxis.grid(True)

# Distance between subplots
fig1.subplots_adjust(wspace=0.3)
plt.show()
```


```python
#fig1.savefig ("delta_ipt.eps", format = 'eps')
fig1.savefig ("delta_ipt.svg", format = "svg", dpi = 1200)
```

## b) 2nd Condition to only have one non-trivial root: 
$(q>0) \wedge (\Delta>0) \wedge (-q + \sqrt{\Delta} \leq 0)$

- The condition $q>0$ implies $R_{L}> \sqrt{2}\omega L$


```python
# We use the previously defined w_c, L_c, and beta_c

# RL is re-defined to comply with the condition:

# Load Resistance [ohms]
R_L = np.arange(np.ceil(np.sqrt(2)*w_c*L_c), 
                1000. + 1., 1.)   # Ceil function --> Up the inferior limit
print(R_L.shape)
print(f'Load Resistance interval: ({R_L[0]}, {R_L[-1]}]')
```

    (620,)
    Load Resistance interval: (381.0, 1000.0]


- With the previous $R_L$ values, determine:
<ol>
<li> Plot of $\Gamma(R_L, k_i) = -q(R_L) + \sqrt{\Delta(R_L, k_i)} \leq 0$.
<li> Only $\Delta > 0$ are importat, so the rest will be nan.
</ol>


```python
# Formula for condition Gamma
Gamma_c = lambda R, beta, k: -((R**2) - (2*beta)) + (np.sqrt(Delta_c(R, beta, k)))

# Create a new mesh
k_i_m, R_L_m = np.meshgrid(k_i, R_L)

print(k_i_m.shape)
print(k_i_m)
print(R_L_m.shape)
print(R_L_m)
print()

# Calculate the Gamma function
f_gamma= Gamma_c(R_L_m, beta_c, k_i_m)
print(f'With sqrt(-1): \n {f_gamma} \n')

# Since 1 is not part of the result
# Invalid values (nan) are set to 1
f_gamma = np.nan_to_num(f_gamma, nan = 1)
print(f'After setting to 1: \n {f_gamma} \n')
```

    (620, 101)
    [[0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     ...
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]]
    (620, 101)
    [[ 381.  381.  381. ...  381.  381.  381.]
     [ 382.  382.  382. ...  382.  382.  382.]
     [ 383.  383.  383. ...  383.  383.  383.]
     ...
     [ 998.  998.  998. ...  998.  998.  998.]
     [ 999.  999.  999. ...  999.  999.  999.]
     [1000. 1000. 1000. ... 1000. 1000. 1000.]]
    
    With sqrt(-1): 
     [[            nan             nan             nan ...             nan
                  nan      0.        ]
     [            nan             nan             nan ...             nan
                  nan      0.        ]
     [            nan             nan             nan ...             nan
                  nan      0.        ]
     ...
     [-12365.74627493 -12364.50058911 -12360.76354272 ...   -486.26641676
        -244.32643007      0.        ]
     [-12336.38550525 -12335.14282056 -12331.41477751 ...   -485.12790536
        -243.75454267      0.        ]
     [-12307.13565211 -12305.89595691 -12302.17688225 ...   -483.99358204
        -243.18475785      0.        ]] 
    
    After setting to 1: 
     [[ 1.00000000e+00  1.00000000e+00  1.00000000e+00 ...  1.00000000e+00
       1.00000000e+00  0.00000000e+00]
     [ 1.00000000e+00  1.00000000e+00  1.00000000e+00 ...  1.00000000e+00
       1.00000000e+00  0.00000000e+00]
     [ 1.00000000e+00  1.00000000e+00  1.00000000e+00 ...  1.00000000e+00
       1.00000000e+00  0.00000000e+00]
     ...
     [-1.23657463e+04 -1.23645006e+04 -1.23607635e+04 ... -4.86266417e+02
      -2.44326430e+02  0.00000000e+00]
     [-1.23363855e+04 -1.23351428e+04 -1.23314148e+04 ... -4.85127905e+02
      -2.43754543e+02  0.00000000e+00]
     [-1.23071357e+04 -1.23058960e+04 -1.23021769e+04 ... -4.83993582e+02
      -2.43184758e+02  0.00000000e+00]] 
    


    <ipython-input-13-0713b495cec5>:2: RuntimeWarning: invalid value encountered in sqrt
      Gamma_c = lambda R, beta, k: -((R**2) - (2*beta)) + (np.sqrt(Delta_c(R, beta, k)))



```python
# Index of an specific value
result = np.where(f_gamma == 1)
print('Tuple of arrays returned : ', result)

# zip the 2 arrays to get the exact coordinates
listOfCoordinates= list(zip(result[0], result[1]))

# iterate over the list of coordinates
#for cord in listOfCoordinates:
#    print(cord)

# Just the first and last coordinate
print(f'1st: {listOfCoordinates[0]},\nlast: {listOfCoordinates[-1]} \n')

# Minimum value in sq_de
#ind = np.unravel_index(np.argmin(sq_delta, axis=None), sq_delta.shape)
#print(ind)
#print(sq_delta[ind])
```

    Tuple of arrays returned :  (array([  0,   0,   0, ..., 156, 156, 156]), array([ 0,  1,  2, ...,  8,  9, 10]))
    1st: (0, 0),
    last: (156, 10) 
    



```python
#interactive plotting in separate window
#%matplotlib qt 

#normal charts inside notebooks
%matplotlib inline 

# Set the plot --> size(width, height)
fig2 = plt.figure(figsize=(15,10))

################
# 1st Subplot
################

# Set the subplot for the two graphs
axs = fig2.add_subplot(1, 2, 1, projection = '3d')

# Complete figure information
fig2.suptitle('Frequency Splitting: $(q>0) \wedge (\Delta>0) \wedge (-q + \sqrt{\Delta} \leq 0)$', 
             fontsize = 20)

# For a simple 3D plot
#ax = fig.gca(projection = '3d')

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R_{L}$) [$\Omega$]')
axs.set_zlabel('Function $\Gamma$')

# Title
axs.set_title('3D Plot of $\Gamma$ Function', fontsize = 15)

# Plot the surface
surf = axs.plot_surface(k_i_m, R_L_m, f_gamma, cmap = cm.coolwarm,
                      linewidth = 0, antialiased = False)

# Figure size (single plot)
#plt.rcParams["figure.figsize"] = (30,30)

# Add a color bar which maps values to colors.
plt.colorbar(surf, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.05)

################
# 2nd Subplot
################

# Set the subplot for the two graphs
axs = fig2.add_subplot(1, 2, 2)

# Plot the contour plot

cp = plt.contour(k_i_m, R_L_m,f_gamma, 
                 levels = [-4e4, -2e4, -1e4, 0], 
                 cmap = cm.Spectral)

# Define the levels to label
#c_labels = [0.e0]
# To specify scientific notation
fmt = ticker.LogFormatterSciNotation()

axs.clabel(cp, inline = 1, colors = 'k', fmt = fmt, fontsize = 10)

# Color bar
plt.colorbar(cp, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.1)

# Information
axs.set_title('Contour Plot of $\Gamma$ Function', fontsize = 15)

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R_{L}$) [$\Omega$]')

axs.xaxis.grid(True)
axs.yaxis.grid(True)

# Distance between subplots
fig2.subplots_adjust(wspace=0.3)
plt.show()
```

### - Another way yo visualize only the regions:


```python
"""
Terms definitions for the IPT problem
"""

# Definition of term p
# Return the value of p term
def p_termi(L, k_i):
    return (L**2)*(1 - (k_i**2))

# Definition of term q
# Return the value of q term
def q_termi(R, L, C):
    return (R**2) - ((2*L)/(C))

# Definition of term r
# Return the value of r term
def r_termi(C):
    return (1/(C**(2)))


# Definition of condition Delta
def Delta_IPT(R, L, C, k_i):
    
    p_aux = p_termi(L, k_i)
    #print(f'p {p_aux}')
    
    q_aux = q_termi(R, L, C)
    #print(f'q {q_aux}')
    
    r_aux = r_termi(C)
    #print(f'r {r_aux}')
    
    # Delta condition
    delta = (q_aux**2) - (4*p_aux*r_aux)
    #print(f'delta {delta}')
    
    return delta
```


```python
# Assume coils' self-inductance equal [H]
L_c = 53.49e-6

C_r = 739.92e-12

# Inductive Coupling Coefficient Range k_i
k_i = np.arange(0, 1 + 0.01, 0.01)
print(k_i.shape)
print(k_i[-1])

# Load Resistance [ohms]
R_L = np.arange(1., 1000. + 1., 1.)
print(R_L.shape)
print(R_L[-1])

# Create a mesh
k_i_m, R_L_m = np.meshgrid(k_i, R_L)
```

    (101,)
    1.0
    (1000,)
    1000.0



```python
# Calculate the function gamma (Condition 2)
# Only check for the region, the value of
# the function is not important

f_gamma = ((q_termi(R_L_m, L_c, C_r) > 0) 
           & (Delta_IPT(R_L_m, L_c, C_r, k_i_m) > 0) 
           & (-q_termi(R_L_m, L_c, C_r) 
              + np.sqrt(Delta_IPT(R_L_m, L_c, C_r, k_i_m)) <= 0)).astype(int)
print(f_gamma.shape)
print(f_gamma)
```

    (1000, 101)
    [[0 0 0 ... 0 0 0]
     [0 0 0 ... 0 0 0]
     [0 0 0 ... 0 0 0]
     ...
     [1 1 1 ... 1 1 1]
     [1 1 1 ... 1 1 1]
     [1 1 1 ... 1 1 1]]


    <ipython-input-25-ae9987a86952>:8: RuntimeWarning: invalid value encountered in sqrt
      + np.sqrt(Delta_IPT(R_L_m, L_c, C_r, k_i_m)) <= 0)).astype(int)



```python
#interactive plotting in separate window
#%matplotlib qt 

#normal charts inside notebooks
%matplotlib inline 

# Set the plot --> size(width, height)
fig2 = plt.figure(figsize=(15,10))

################
# 1st Subplot
################

# Set the subplot for the two graphs
axs = fig2.add_subplot(1, 2, 1, projection = '3d')

# Complete figure information
fig2.suptitle('Frequency Splitting: $(q>0) \wedge (\Delta>0) \wedge (-q + \sqrt{\Delta} \leq 0)$', 
             fontsize = 20)

# For a simple 3D plot
#ax = fig.gca(projection = '3d')

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R^{\'}_{L}$) [$\Omega$]')
axs.set_zlabel('Function $\Gamma$')

# Title
axs.set_title('3D Plot of $\Gamma$ Function', fontsize = 15)

# Plot the surface
surf = axs.plot_surface(k_i_m, R_L_m, f_gamma, cmap = cm.coolwarm,
                      linewidth = 0, antialiased = False)

# Figure size (single plot)
#plt.rcParams["figure.figsize"] = (30,30)

# Add a color bar which maps values to colors.
plt.colorbar(surf, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.05)

################
# 2nd Subplot
################

# Set the subplot for the two graphs
axs = fig2.add_subplot(1, 2, 2)

# Plot the contour plot

cp = plt.contourf(k_i_m, R_L_m,f_gamma, 
                 levels = [-1, 0, 1],
                 cmap = cm.Spectral)

# Define the levels to label
#c_labels = [0.e0]
# To specify scientific notation
fmt = ticker.LogFormatterSciNotation()

axs.clabel(cp, inline = 1, colors = 'k', fmt = fmt, fontsize = 10)

# Color bar
plt.colorbar(cp, orientation = 'horizontal', shrink=0.6, aspect=20, pad = 0.1)

# Information
axs.set_title('Contour Plot of $\Gamma$ Function', fontsize = 15)

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R^{\'}_{L}$) [$\Omega$]')

axs.xaxis.grid(True)
axs.yaxis.grid(True)

# Distance between subplots
fig2.subplots_adjust(wspace=0.5)
plt.show()
```


```python
fig2.savefig ("gamma_ipt.svg", format = "svg", dpi = 1200)
```

## c) 3nd Condition to only have one non-trivial root: 
$(q<0) \wedge (\Delta>0) \wedge (-q - \sqrt{\Delta} \leq 0)$ <br>
This is not a condition to consider since it will provide at least <br>
one non-trivial root.

- The condition $q<0$ implies $R_L < \sqrt{2}\omega L$


```python
# We use the previously defined w_c, L_c, and beta_c

# RL is re-defined to comply with the condition:

# Load Resistance [ohms]
R_L = np.arange(1., (np.sqrt(2)*w_c*L_c), 1.)
print(R_L.shape)
print(f'Load Resistance interval: [{R_L[0]}, {R_L[-1]})')
```

    (380,)
    Load Resistance interval: [1.0, 380.0)


- With the previous $R_L$ values, determine:
<ol>
<li> Plot of $\Psi(R_L, k_i) = -q(R_L) - \sqrt{\Delta(R_L, k_i)} \leq 0$.
<li> Only $\Delta > 0$ are importat, so the rest will be nan.
</ol>


```python
# Formula for condition Psi
Psi_c = lambda R, beta, k: -((R**2) - (2*beta)) - (np.sqrt(Delta_c(R, beta, k)))

# Create a new mesh
k_i_m, R_L_m = np.meshgrid(k_i, R_L)

print(k_i_m.shape)
print(k_i_m)
print(R_L_m.shape)
print(R_L_m)
print()

# Calculate the Psi function
f_psi= Psi_c(R_L_m, beta_c, k_i_m)
print(f'With sqrt(-1): \n {f_psi} \n')

# Since 1 is not part of the result
# Invalid values (nan) are set to 1
f_psi = np.nan_to_num(f_psi, nan = 1.)
print(f'After setting to 1: \n {f_psi} \n')
```

    (380, 101)
    [[0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     ...
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]
     [0.   0.01 0.02 ... 0.98 0.99 1.  ]]
    (380, 101)
    [[  1.   1.   1. ...   1.   1.   1.]
     [  2.   2.   2. ...   2.   2.   2.]
     [  3.   3.   3. ...   3.   3.   3.]
     ...
     [378. 378. 378. ... 378. 378. 378.]
     [379. 379. 379. ... 379. 379. 379.]
     [380. 380. 380. ... 380. 380. 380.]]
    
    With sqrt(-1): 
     [[            nan 143239.12214754 141740.019677   ...   2891.66489673
        1445.83234529      0.        ]
     [            nan 143611.91175624 141894.01763856 ...   2891.72612341
        1445.86264939      0.        ]
     [            nan             nan 142173.38925858 ...   2891.82817368
        1445.91315905      0.        ]
     ...
     [            nan             nan             nan ...             nan
                  nan      0.        ]
     [            nan             nan             nan ...             nan
                  nan      0.        ]
     [            nan             nan             nan ...             nan
                  nan      0.        ]] 
    
    After setting to 1: 
     [[1.00000000e+00 1.43239122e+05 1.41740020e+05 ... 2.89166490e+03
      1.44583235e+03 0.00000000e+00]
     [1.00000000e+00 1.43611912e+05 1.41894018e+05 ... 2.89172612e+03
      1.44586265e+03 0.00000000e+00]
     [1.00000000e+00 1.00000000e+00 1.42173389e+05 ... 2.89182817e+03
      1.44591316e+03 0.00000000e+00]
     ...
     [1.00000000e+00 1.00000000e+00 1.00000000e+00 ... 1.00000000e+00
      1.00000000e+00 0.00000000e+00]
     [1.00000000e+00 1.00000000e+00 1.00000000e+00 ... 1.00000000e+00
      1.00000000e+00 0.00000000e+00]
     [1.00000000e+00 1.00000000e+00 1.00000000e+00 ... 1.00000000e+00
      1.00000000e+00 0.00000000e+00]] 
    


    <ipython-input-12-f4cf8d1a4a28>:2: RuntimeWarning: invalid value encountered in sqrt
      Psi_c = lambda R, beta, k: -((R**2) - (2*beta)) - (np.sqrt(Delta_c(R, beta, k)))



```python
# Index of an specific value
result = np.where(f_psi == 1)
print('Tuple of arrays returned : ', result)

# zip the 2 arrays to get the exact coordinates
listOfCoordinates= list(zip(result[0], result[1]))

# iterate over the list of coordinates
#for cord in listOfCoordinates:
#    print(cord)

# Just the first and last coordinate
print(f'1st: {listOfCoordinates[0]},\nlast: {listOfCoordinates[-1]} \n')

# Minimum value in sq_de
#ind = np.unravel_index(np.argmin(sq_delta, axis=None), sq_delta.shape)
#print(ind)
#print(sq_delta[ind])
```

    Tuple of arrays returned :  (array([  0,   1,   2, ..., 379, 379, 379]), array([ 0,  0,  0, ..., 97, 98, 99]))
    1st: (0, 0),
    last: (379, 99) 
    



```python
#interactive plotting in separate window
#%matplotlib qt 

#normal charts inside notebooks
%matplotlib inline 

# To specify scientific notation
fmt = ticker.LogFormatterSciNotation()

# Set the plot --> size(width, height)
fig = plt.figure(figsize=(15,10))

################
# 1st Subplot
################

# Set the subplot for the two graphs
axs = fig.add_subplot(1, 2, 1, projection = '3d')

# Complete figure information
fig.suptitle('Frequency Splitting: $(q<0) \wedge (\Delta>0) \wedge (-q - \sqrt{\Delta} \leq 0)$', 
             fontsize = 20)

# For a simple 3D plot
#ax = fig.gca(projection = '3d')

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R_{L}$) [$\Omega$]')
axs.set_zlabel('Function $\Psi$')

# Title
axs.set_title('3D Plot of $\Psi$ Function', fontsize = 15)

# Plot the surface
surf = axs.plot_surface(k_i_m, R_L_m, f_psi, cmap = cm.coolwarm,
                      linewidth = 0, antialiased = False)

# Figure size (single plot)
#plt.rcParams["figure.figsize"] = (30,30)

# Add a color bar which maps values to colors.
plt.colorbar(surf, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.05)

################
# 2nd Subplot
################

# Set the subplot for the two graphs
axs = fig.add_subplot(1, 2, 2)

# Plot the contour plot

cp = plt.contour(k_i_m, R_L_m,f_psi,
                 levels = [0, 2e4, 6e4, 1e5],
                 cmap = cm.Spectral)

# Define the levels to label
#c_labels = [-2e10, 0, 1e10]

axs.clabel(cp, inline = 1, colors = 'k', fmt = fmt, fontsize = 10)

# Color bar
plt.colorbar(cp, format = fmt, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.1)

# Information
axs.set_title('Contour Plot of $\Psi$ Function', fontsize = 15)

# Axis labels
axs.set_xlabel('Inductive Coupling Coefficient ($k_{i}$)')
axs.set_ylabel('Load Resistance ($R_{L}$) [$\Omega$]')

axs.xaxis.grid(True)
axs.yaxis.grid(True)

# Distance between subplots
fig.subplots_adjust(wspace=0.3)
plt.show()
```

# 2. Freq. Splitting Model for Active Compensation Network


## Equations Derivation


```python
# Definition of the symbolic variables

# Resistances --> Coils and Load
R_1, R_2, R_L = symbols('R_1, R_2, R_L',
                       real = True)

# Inductance and Capacitance
C_1, C_2, L_1, L_2, M_i = symbols('C_1, C_2, L_1,L_2, M_i', 
                                 real = True)

# Reactance
X_1, X_2, X_m = symbols('X_1, X_2, X_m',
                       real = True)

# Angular frequency
w_e = symbols('w_e', real = True)

# Input Voltage and Current
V_1, I_1 = symbols('V_1, I_1', real = True)
```

- Equivalent Circuit's Matrix Representation: $[Z][I] = [V]$


```python
# Impedance Matriz [Z]
Z_m = Matrix([[R_1 + (I*X_1), I*X_m],
            [I*X_m, R_2 + R_L + (I*X_2)]])

# Voltage Vector
V_m = Matrix([[V_1], [0]])
```


```python
# Solve the System of Equations
I_m = Z_m.inv()*V_m
```


```python
# Input Current
result_I = "Input Current: \n $${} = {}$$".format(latex(I_1), 
                                                 latex(I_m[0]))
display(Latex(result_I))

# Separate Numerator and denominator of I1
I_1_n, I_1_d = fraction(I_m[0])
```


Input Current: 
 $$I_{1} = \frac{V_{1} \left(- i R_{2} - i R_{L} + X_{2}\right)}{- i R_{1} R_{2} - i R_{1} R_{L} + R_{1} X_{2} + R_{2} X_{1} + R_{L} X_{1} + i X_{1} X_{2} - i X_{m}^{2}}$$



```python
# Operate separately each polynomial
display(factor(I_1_n))
display(factor(I_1_d + (I*(X_m**2))))
```


$\displaystyle - i V_{1} \left(R_{2} + R_{L} + i X_{2}\right)$



$\displaystyle - i \left(R_{1} + i X_{1}\right) \left(R_{2} + R_{L} + i X_{2}\right)$



```python
# Reconstruct the factored result
I_1 = (factor(I_1_n))/(factor(I_1_d + (I*(X_m**2))) 
                       - (I*(X_m**2)))
display(simplify(I_1))
```


$\displaystyle \frac{V_{1} \left(R_{2} + R_{L} + i X_{2}\right)}{X_{m}^{2} + \left(R_{1} + i X_{1}\right) \left(R_{2} + R_{L} + i X_{2}\right)}$


- Input Impedance Expression $Z_{in}$:


```python
Z_in = V_1*(1/simplify(I_1))
display(Z_in)
```


$\displaystyle \frac{X_{m}^{2} + \left(R_{1} + i X_{1}\right) \left(R_{2} + R_{L} + i X_{2}\right)}{R_{2} + R_{L} + i X_{2}}$


- Freq. Splitting occurs when $\Im(Z_{in}) = 0$


```python
# Assume both reactance are equal for simplicity
# Common reactance X_e
X_e = symbols('X_e', real = True)

# Substitute the variable
Z_in = Z_in.subs([(X_1, X_e), (X_2, X_e)])

# Imaginary part of Zin
simplify(im(Z_in))
```




$\displaystyle \frac{X_{e} \left(- R_{1} R_{2} - R_{1} R_{L} + X_{e}^{2} - X_{m}^{2} + \left(R_{2} + R_{L}\right) \left(R_{1} + R_{2} + R_{L}\right)\right)}{X_{e}^{2} + \left(R_{2} + R_{L}\right)^{2}}$




```python
# Operate separately
factor(-(R_1*R_2) - (R_1*R_L) + ((R_2 + R_L)*(R_1 + R_2 + R_L)))
```




$\displaystyle \left(R_{2} + R_{L}\right)^{2}$



- Analyze the simple non-trivial root $\chi = 0$


```python
# Definition of the symbolic variables

# Capacitances T-Model
C_ep, C_Mp = symbols('C_ep, C_Mp',
                       real = True)

# Capacitances Pi-Model
C_e, C_M = symbols('C_e, C_M',
                       real = True)

# Capacitive Coupling Coefficient
k_c = symbols('k_c', real = True)
```


```python
# Replace the k_c
C_M = k_c*C_e

# Equivalent Pi --> T Model
C_ep = C_e + (2*C_M)
C_Mp = (C_e/C_M)*(C_ep)
```


```python
# Operation for the main resonance frequency
# Defining the equivalent parallel capacitance
C_pp = ((1/C_ep) + (1/C_Mp))**(-1)
simplify(C_pp)
```




$\displaystyle \frac{C_{e} \left(2 k_{c} + 1\right)}{k_{c} + 1}$




```python
# C_Mp also appears in the formula
simplify(C_Mp)
```




$\displaystyle 2 C_{e} + \frac{C_{e}}{k_{c}}$



## a) 1st Condition to only have one non-trivial root: $\Delta < 0$
<ol>
<li> $L$ and $k_i$ are fixed from the IPT Model.
<li> External capacitance $C_e$ is defined according to standard values.
<li> $0 < k_c < 1$
<li> Frequency dependence is visualized separetely.
</ol>


```python
# Defining the Inductive Coupling Coefficient k_i
k_i = np.array([0.214])

# Assume coils' self-inductance equal [H]
L_c = 53.49e-6

# External Capacitance [F]
C_e = np.array([0.668e-9])
```

- Define methods to calculate the condition $\Delta$


```python
##################################
## Define Methods for Calculations
##################################

# Definition of Equivalent Parallel Capacitance C_pp
def C_ppri(C, k_c):
    return (C*(1 + (2*k_c)))/(1 + k_c)

# Definition of the Mutual Capacitance in T-Model (C_Mp)
def C_Mpri(C, k_c):
    return C*((1/k_c) + 2)

# Definition of term p
# Return the value of p term
def p_term(L, k_i):
    return (L**2)*(1 - (k_i**2))

# Definition of term q
# Return the value of q term
def q_term(R, L, C, k_c, k_i):
    return (R**2) + (2*L*((k_i/C_Mpri(C, k_c)) 
                          - (1/C_ppri(C, k_c))))

# Definition of term r
# Return the value of r term
def r_term(C, k_c):
    return (1/((C_ppri(C, k_c))**2)) - (1/((C_Mpri(C, k_c))**2))


# Definition of condition Delta
def Delta_c1(R, L, C, k_c, k_i):
    
    # Auxiliar variables
    C_Mp = C_Mpri(C, k_c)
    #print(f'C_Mp {C_Mp}')
    
    C_pp = C_ppri(C, k_c)
    #print(f'C_pp {C_pp}')
    
    p_aux = p_term(L, k_i)
    #print(f'p {p_aux}')
    
    q_aux = q_term(R, L, C, k_c, k_i)
    #print(f'q {q_aux}')
    
    r_aux = r_term(C, k_c)
    #print(f'r {r_aux}')
    
    # Delta condition
    delta = (q_aux**2) - (4*p_aux*r_aux)
    #print(f'delta {delta}')
    
    return delta
```


```python
# Capacitive Coupling Coefficient Range k_c
k_c = np.arange(0.01, 1 + 0.01, 0.01)
print(k_c.shape)
print(k_c)

# Load Resistance [ohms]
R_L = np.arange(1., 1000. + 1., 1.)
print(R_L.shape)
print(R_L[-1])

# Create a mesh
k_c_m, R_L_m = np.meshgrid(k_c, R_L)

# Calculate the function Delta
f_delta = Delta_c1(R_L_m, L_c, C_e, k_c_m, k_i)
print(f_delta.shape)
```

    (100,)
    [0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13 0.14
     0.15 0.16 0.17 0.18 0.19 0.2  0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28
     0.29 0.3  0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4  0.41 0.42
     0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5  0.51 0.52 0.53 0.54 0.55 0.56
     0.57 0.58 0.59 0.6  0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7
     0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8  0.81 0.82 0.83 0.84
     0.85 0.86 0.87 0.88 0.89 0.9  0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98
     0.99 1.  ]
    (1000,)
    1000.0
    (1000, 100)



```python
#interactive plotting in separate window
#%matplotlib qt 

#normal charts inside notebooks
%matplotlib inline 

# Set the plot --> size(width, height)
fig1 = plt.figure(figsize=(15,10))

################
# 1st Subplot
################

# Set the subplot for the two graphs
axs = fig1.add_subplot(1, 2, 1, projection = '3d')

# Complete figure information
fig1.suptitle('Frequency Splitting Active CN: $\Delta$ Function', fontsize = 20)

# For a simple 3D plot
#ax = fig.gca(projection = '3d')

# Axis labels
axs.set_xlabel('Capacitive Coupling Coefficient ($k_c$)')
axs.set_ylabel('Load Resistance ($R_L$) [$\Omega$]')
axs.set_zlabel('Function $\Delta$')

# Title
axs.set_title('3D Plot of $\Delta$ Function', fontsize = 15)

# Plot the surface
surf = axs.plot_surface(k_c_m, R_L_m, f_delta, cmap = cm.coolwarm,
                      linewidth = 0, antialiased = False)

# Figure size (single plot)
#plt.rcParams["figure.figsize"] = (30,30)

# Add a color bar which maps values to colors.
plt.colorbar(surf, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.05)

################
# 2nd Subplot
################

# Set the subplot for the two graphs
axs = fig1.add_subplot(1, 2, 2)

# Plot the contour plot
cp = plt.contour(k_c_m, R_L_m, f_delta,
                  levels = [-0.5e11, 0, 1.5e11, 3e11, 4.5e11],
                  cmap = cm.Spectral)

# Define the levels to label
#c_labels = [0.e0]
# To specify scientific notation
fmt = ticker.LogFormatterSciNotation()

axs.clabel(cp, inline = 1, colors = 'k', fmt = fmt, fontsize = 10)

# Color bar
plt.colorbar(cp, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.1)

# Information
axs.set_title('Contour Plot of $\Delta$ Function', fontsize = 15)

# Axis labels
axs.set_xlabel('Capacitive Coupling Coefficient ($k_c$)')
axs.set_ylabel('Load Resistance ($R_L$) [$\Omega$]')

axs.xaxis.grid(True)
axs.yaxis.grid(True)

# Distance between subplots
fig1.subplots_adjust(wspace=0.3)
plt.show()
```


```python
#fig1.savefig ("delta_ipt.eps", format = 'eps')
fig1.savefig ("delta_acn.svg", format = "svg", dpi = 1200)
```

- Plot the behavior of resonance frequency (central) w.r.t. $k_c$ <br>
$w_r = \sqrt{\frac{1}{L {C^{\prime}_{p}}}}$


```python
# Define central resonance frequency
def f_operation(L, C, k_c):
    freq = (1/(2*np.pi))*np.sqrt((1)/(L*C_ppri(C, k_c)))
    return freq
```


```python
# Calculate the resonance frequency
f_res = (1e-3)*f_operation(L_c, C_e, k_c)
print(f_res.shape)
```

    (100,)



```python
#interactive plotting in separate window
%matplotlib qt 

#normal charts inside notebooks
#%matplotlib inline 

# Set the plot --> size(width, height)
fig2 = plt.figure(figsize=(10,10))
ax = fig2.gca()

# Axis labels
ax.set_xlabel('Capacitive Coupling Coefficient $k_{c}$')
ax.set_ylabel('Resonance Frequency ($f_{r}$) [kHz]')

# Title
ax.set_title('System\'s Resonance Frequency', fontsize = 25)

# Plot the surface
lines = ax.plot(k_c, f_res, 'g', linewidth = 2.5)

ax.xaxis.grid(True)
ax.yaxis.grid(True)

# Axis limits
ax.set_xlim([0, 1])

# Legend
#plt.legend(('Maxwell'), loc = 'upper right')

plt.show()


```


```python
fig2.savefig ("fresonant_acn.svg", format = "svg", dpi = 1200)
```

## b) 2nd Condition to only have one non-trivial root: 
$(q>0) \wedge (\Delta>0) \wedge (-q + \sqrt{\Delta} \leq 0)$


```python
# Defining the Inductive Coupling Coefficient k_i
k_i = np.array([0.214])

# Assume coils' self-inductance equal [H]
L_c = np.array([53.49e-6])

# External Capacitance [F]
C_e = np.array([0.66e-9])
```


```python
# Capacitive Coupling Coefficient Range k_c
k_c = np.arange(0.01, 1 + 0.01, 0.01)
print(k_c.shape)
print(k_c)

# Load Resistance [ohms]
R_L = np.arange(1., 1000. + 1., 1.)
print(R_L.shape)
print(R_L[-1])

# Create a mesh
k_c_m, R_L_m = np.meshgrid(k_c, R_L)
```

    (100,)
    [0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13 0.14
     0.15 0.16 0.17 0.18 0.19 0.2  0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28
     0.29 0.3  0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4  0.41 0.42
     0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5  0.51 0.52 0.53 0.54 0.55 0.56
     0.57 0.58 0.59 0.6  0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7
     0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8  0.81 0.82 0.83 0.84
     0.85 0.86 0.87 0.88 0.89 0.9  0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98
     0.99 1.  ]
    (1000,)
    1000.0



```python
# Calculate the function gamma (Condition 2)
# Only check for the region, the value of
# the function is not important

f_gamma = ((q_term(R_L_m, L_c, C_e, k_c_m, k_i) > 0) 
           & (Delta_c1(R_L_m, L_c, C_e, k_c_m, k_i) > 0) 
           & (-q_term(R_L_m, L_c, C_e, k_c_m, k_i) 
              + np.sqrt(Delta_c1(R_L_m, L_c, C_e, k_c_m, k_i)) <= 0)).astype(int)
print(f_gamma.shape)
print(f_gamma)
```

    (1000, 100)
    [[0 0 0 ... 0 0 0]
     [0 0 0 ... 0 0 0]
     [0 0 0 ... 0 0 0]
     ...
     [1 1 1 ... 1 1 1]
     [1 1 1 ... 1 1 1]
     [1 1 1 ... 1 1 1]]


    <ipython-input-17-fd1090a33310>:8: RuntimeWarning: invalid value encountered in sqrt
      + np.sqrt(Delta_c1(R_L_m, L_c, C_e, k_c_m, k_i)) <= 0)).astype(int)



```python
# Index of an specific value
result = np.where(np.isnan(f_gamma))
print('Tuple of arrays returned : ', result)

# zip the 2 arrays to get the exact coordinates
listOfCoordinates= list(zip(result[0], result[1]))

# iterate over the list of coordinates
#for cord in listOfCoordinates:
#    print(cord)

# Just the first and last coordinate
print(f'1st: {listOfCoordinates[0]},\nlast: {listOfCoordinates[-1]} \n')

# Minimum value in sq_de
#ind = np.unravel_index(np.argmin(sq_delta, axis=None), sq_delta.shape)
#print(ind)
#print(sq_delta[ind])
```


```python
#interactive plotting in separate window
#%matplotlib qt 

#normal charts inside notebooks
%matplotlib inline 

# Set the plot --> size(width, height)
fig1 = plt.figure(figsize=(15,10))

################
# 1st Subplot
################

# Set the subplot for the two graphs
axs = fig1.add_subplot(1, 2, 1, projection = '3d')

# Complete figure information
fig1.suptitle('Frequency Splitting: $(q>0) \wedge (\Delta>0) \wedge (-q + \sqrt{\Delta} \leq 0)$', 
             fontsize = 20)

# For a simple 3D plot
#ax = fig.gca(projection = '3d')

# Axis labels
axs.set_xlabel('Capacitive Coupling Coefficient ($k_c$)')
axs.set_ylabel('Load Resistance ($R_{L}$) [$\Omega$]')
axs.set_zlabel('Function $\Gamma$')

# Title
axs.set_title('3D Plot of $\Gamma$ Function', fontsize = 15)

# Plot the surface
surf = axs.plot_surface(k_c_m, R_L_m, f_gamma, cmap = cm.coolwarm,
                      linewidth = 0, antialiased = False)

# Figure size (single plot)
#plt.rcParams["figure.figsize"] = (30,30)

# Add a color bar which maps values to colors.
plt.colorbar(surf, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.05)

################
# 2nd Subplot
################

# Set the subplot for the two graphs
axs = fig1.add_subplot(1, 2, 2)

# Plot the contour plot
cp = plt.contourf(k_c_m, R_L_m, f_gamma,
                  levels = [-1, 0, 1],
                 cmap = cm.Spectral)

# Define the levels to label
#c_labels = [0.e0]
# To specify scientific notation
fmt = ticker.LogFormatterSciNotation()

axs.clabel(cp, inline = 1, colors = 'k', fmt = fmt, fontsize = 10)

# Color bar
plt.colorbar(cp, orientation = 'horizontal', shrink=0.6,
             aspect=20, pad = 0.09)

# Information
axs.set_title('Contour Plot of $\Gamma$ Function', fontsize = 15)

# Axis labels
axs.set_xlabel('Capacitive Coupling Coefficient ($k_c$)')
axs.set_ylabel('Load Resistance ($R_L$) [$\Omega$]')

axs.xaxis.grid(True)
axs.yaxis.grid(True)

# Distance between subplots
fig1.subplots_adjust(wspace=0.3)
plt.show()
```


```python
fig1.savefig ("fresonant_acn.svg", format = "svg", dpi = 1200)
```


```python
# 3. $\pi$-Model and $T$-Model equivalent equations

- Equivalent capacitive model for 4 plates
```

## $\pi$-Model:


```python
# Define the symbolic variables

# Capacitance
c_1, c_2, c_M = symbols('c_1 c_2 c_M', real = True)

# Angular frequency
w = symbols('w', real = True)
```


```python
# Common denominator
com_D = (1/(w*I))*(1/c_1 + 1/c_2 + 1/c_M)
```


```python
# T-model equivalent c'_1
Zp_1 = ((1/(w*I))**2)*((1/c_1)*(1/c_M))/com_D
simplify(Zp_1)
```


```python
# T-model equivalent c'_2
Zp_2 = ((1/(w*I))**2)*((1/c_2)*(1/c_M))/com_D
simplify(Zp_2)
```


```python
# T-model equivalent c'_M
Zp_M = ((1/(w*I))**2)*((1/c_1)*(1/c_2))/com_D
simplify(Zp_M)
```

## $T$-Model:


```python
# Define the symbolic variables

# Capacitance
cp_1, cp_2, cp_M = symbols('cp_1 cp_2 cp_M', real = True)
```


```python
# Common numerator
com_N = ((1/(w*I))**2)*((1/(cp_1*cp_M)) + (1/(cp_1*cp_2)) + (1/(cp_M*cp_2)))
```


```python
# Pi-model equivalent c_1
Z_1 = (com_N)/(1/(w*cp_2*I))
simplify(Z_1)
```


```python
# Pi-model equivalent c_2
Z_2 = (com_N)/(1/(w*cp_1*I))
simplify(Z_2)
```


```python
# Pi-model equivalent c_M
Z_M = (com_N)/(1/(w*cp_M*I))
simplify(Z_M)
```

# 4. Latex Formulas 

$\Im(Z_{in}) = \frac{\chi [(R^{'}_{L} + R_{2})^{2} + \chi^{2} - (\omega M)^{2}]}{(R^{'}_{L} + R_{2})^{2} + \chi^{2}} = 0$

$\chi = \omega L - \frac{1}{\omega C}$

$k_{i} = \frac{M}{L} $ 

$R^{'}_{L} \gg R_{2} \Rightarrow (R^{'}_{L} + R_{2}) \approx R^{'}_{L}$

$f(\omega^{2}) = \left[ L^{2}(1 - k_{i}^{2}) \right] (\omega^{2})^{2} + \left[ R^{'2}_{L} - \frac{2L}{C} \right] \omega^{2} + \frac{1}{C^{2}}$

$\omega = \sqrt{\frac{-q \pm \sqrt{\Delta}}{2p}}$

$\Delta =  q^{2} - 4pr$

$-q \pm \sqrt{\Delta} \leq 0$

$\Im(Z_{in}) = \frac{\chi [(R^{'}_{L} + R_{2})^{2} + \chi^{2} - \chi_{m}^{2}]}{(R^{'}_{L} + R_{2})^{2} + \chi^{2}} = 0$

$\chi = \omega L - \frac{1}{\omega C^{'}_{p}}$

$\chi_{m} = \omega M - \frac{1}{\omega C^{'}_{M}}$

$\frac{1}{C^{'}_{p}} = \frac{1}{C^{'}} + \frac{1}{C^{'}_{M}} = \frac{C + C_{M}}{C(C + 2C_{M})}$

$f(\omega^{2}) = \left[ L^{2}(1 - k_{i}^{2}) \right] (\omega^{2})^{2} + \left[ R^{'2}_{L} + 2L \left(\frac{k_{i}}{C^{'}_{M}} - \frac{1}{C^{'}_{p}} \right) \right] \omega^{2} + \left( \frac{1}{C^{'2}_{p}} - \frac{1}{C^{'2}_{M}} \right)$

$k_{c} = \frac{C_{M}}{C} $ 


```python

```

3746cff2a0ac284d2be70101aaf1f9ee848ec599
ipynb
Jupyter Notebook
Active_CN_Model.ipynb
SanTT19/Symbolic_Python
3bfaec4fc62de52e6a592cf3db333e605a35fcbc
MIT
Active_CN_Model.ipynb
SanTT19/Symbolic_Python
3bfaec4fc62de52e6a592cf3db333e605a35fcbc
MIT
Active_CN_Model.ipynb
SanTT19/Symbolic_Python
3bfaec4fc62de52e6a592cf3db333e605a35fcbc
MIT
Qwen/Qwen-72B
1. YES
2. YES