<jupyter_start><jupyter_text>1. IPT System Frequency Splitting Model  Equations Derivation<jupyter_code># 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)
# 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))<jupyter_output><empty_output><jupyter_text>- Equivalent Circuit's Matrix Representation: $[Z][I] = [V]$<jupyter_code># 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]])
# Solve the System of Equations
I_m = Z_m.inv()*V_m
# 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])
# Operate separately each polynomial
display(factor(I_1_n))
display(factor(I_1_d + (I*(M_i**2)*(w_e**2))))
# 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))<jupyter_output><empty_output><jupyter_text>- Input Impedance Expression $Z_{in}$:<jupyter_code>Z_in = V_1*(1/simplify(I_1))
display(Z_in)<jupyter_output><empty_output><jupyter_text>- Freq. Splitting occurs when $\Im(Z_{in}) = 0$<jupyter_code># 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))
# Operate separately
factor(-(R_1*R_2) - (R_1*R_L) + ((R_2 + R_L)*(R_1 + R_2 + R_L)))<jupyter_output><empty_output><jupyter_text>a) 1st Condition to only have one non-trivial root: $\Delta < 0$<jupyter_code># 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
# 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)
#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()
#fig1.savefig ("delta_ipt.eps", format = 'eps')
fig1.savefig ("delta_ipt.svg", format = "svg", dpi = 1200)<jupyter_output><empty_output><jupyter_text>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$<jupyter_code># 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]}]')<jupyter_output>(620,)
Load Resistance interval: (381.0, 1000.0]<jupyter_text>- With the previous $R_L$ values, determine: Plot of $\Gamma(R_L, k_i) = -q(R_L) + \sqrt{\Delta(R_L, k_i)} \leq 0$. Only $\Delta > 0$ are importat, so the rest will be nan.<jupyter_code># 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')
# 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])
#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()<jupyter_output><empty_output><jupyter_text>- Another way yo visualize only the regions:<jupyter_code>"""
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
# 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)
# 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)
#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()
fig2.savefig ("gamma_ipt.svg", format = "svg", dpi = 1200)<jupyter_output><empty_output><jupyter_text>c) 3nd Condition to only have one non-trivial root: $(q0) \wedge (-q - \sqrt{\Delta} \leq 0)$ This is not a condition to consider since it will provide at least one non-trivial root. - The condition $q<0$ implies $R_L < \sqrt{2}\omega L$<jupyter_code># 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]})')<jupyter_output>(380,)
Load Resistance interval: [1.0, 380.0)<jupyter_text>- With the previous $R_L$ values, determine: Plot of $\Psi(R_L, k_i) = -q(R_L) - \sqrt{\Delta(R_L, k_i)} \leq 0$. Only $\Delta > 0$ are importat, so the rest will be nan.<jupyter_code># 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')
# 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])
#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()<jupyter_output><empty_output><jupyter_text>2. Freq. Splitting Model for Active Compensation Network  Equations Derivation<jupyter_code># 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)<jupyter_output><empty_output><jupyter_text>- Equivalent Circuit's Matrix Representation: $[Z][I] = [V]$<jupyter_code># 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]])
# Solve the System of Equations
I_m = Z_m.inv()*V_m
# 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])
# Operate separately each polynomial
display(factor(I_1_n))
display(factor(I_1_d + (I*(X_m**2))))
# 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))<jupyter_output><empty_output><jupyter_text>- Input Impedance Expression $Z_{in}$:<jupyter_code>Z_in = V_1*(1/simplify(I_1))
display(Z_in)<jupyter_output><empty_output><jupyter_text>- Freq. Splitting occurs when $\Im(Z_{in}) = 0$<jupyter_code># 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))
# Operate separately
factor(-(R_1*R_2) - (R_1*R_L) + ((R_2 + R_L)*(R_1 + R_2 + R_L)))<jupyter_output><empty_output><jupyter_text>- Analyze the simple non-trivial root $\chi = 0$<jupyter_code># 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)
# 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)
# Operation for the main resonance frequency
# Defining the equivalent parallel capacitance
C_pp = ((1/C_ep) + (1/C_Mp))**(-1)
simplify(C_pp)
# C_Mp also appears in the formula
simplify(C_Mp)<jupyter_output><empty_output><jupyter_text>a) 1st Condition to only have one non-trivial root: $\Delta < 0$ $L$ and $k_i$ are fixed from the IPT Model. External capacitance $C_e$ is defined according to standard values. $0 < k_c < 1$ Frequency dependence is visualized separetely.<jupyter_code># 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])<jupyter_output><empty_output><jupyter_text>- Define methods to calculate the condition $\Delta$<jupyter_code>##################################
## 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
# 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)
#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()
#fig1.savefig ("delta_ipt.eps", format = 'eps')
fig1.savefig ("delta_acn.svg", format = "svg", dpi = 1200)<jupyter_output><empty_output><jupyter_text>- Plot the behavior of resonance frequency (central) w.r.t. $k_c$ $w_r = \sqrt{\frac{1}{L {C^{\prime}_{p}}}}$<jupyter_code># 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
# Calculate the resonance frequency
f_res = (1e-3)*f_operation(L_c, C_e, k_c)
print(f_res.shape)
#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()


fig2.savefig ("fresonant_acn.svg", format = "svg", dpi = 1200)<jupyter_output><empty_output><jupyter_text>b) 2nd Condition to only have one non-trivial root: $(q>0) \wedge (\Delta>0) \wedge (-q + \sqrt{\Delta} \leq 0)$<jupyter_code># 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])
# 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 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)
# 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])
#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()
fig1.savefig ("fresonant_acn.svg", format = "svg", dpi = 1200)
# 3. $\pi$-Model and $T$-Model equivalent equations

- Equivalent capacitive model for 4 plates<jupyter_output><empty_output><jupyter_text>$\pi$-Model:<jupyter_code># 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)
# Common denominator
com_D = (1/(w*I))*(1/c_1 + 1/c_2 + 1/c_M)
# T-model equivalent c'_1
Zp_1 = ((1/(w*I))**2)*((1/c_1)*(1/c_M))/com_D
simplify(Zp_1)
# T-model equivalent c'_2
Zp_2 = ((1/(w*I))**2)*((1/c_2)*(1/c_M))/com_D
simplify(Zp_2)
# T-model equivalent c'_M
Zp_M = ((1/(w*I))**2)*((1/c_1)*(1/c_2))/com_D
simplify(Zp_M)<jupyter_output><empty_output><jupyter_text>$T$-Model:<jupyter_code># Define the symbolic variables

# Capacitance
cp_1, cp_2, cp_M = symbols('cp_1 cp_2 cp_M', real = True)
# Common numerator
com_N = ((1/(w*I))**2)*((1/(cp_1*cp_M)) + (1/(cp_1*cp_2)) + (1/(cp_M*cp_2)))
# Pi-model equivalent c_1
Z_1 = (com_N)/(1/(w*cp_2*I))
simplify(Z_1)
# Pi-model equivalent c_2
Z_2 = (com_N)/(1/(w*cp_1*I))
simplify(Z_2)
# Pi-model equivalent c_M
Z_M = (com_N)/(1/(w*cp_M*I))
simplify(Z_M)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1677845423588897, 'fraction_numerical': 0.029271145048670587, 'mean_word_length': 4.072226099092812, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 113, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '177772', 'license': 'MIT', 'n_tokens_mistral': 12206, 'n_tokens_neox': 11476, 'n_words': 2905, 'path': 'Active_CN_Model.ipynb', 'repo_name': 'SanTT19/Symbolic_Python', 'size': '29073'}