From bbce02aa2fe47f7b2910fad7cf07db158cb1f083 Mon Sep 17 00:00:00 2001
From: Jan H <Jan.Habscheid@web.de>
Date: Mon, 26 Aug 2024 10:13:20 +0200
Subject: [PATCH] DLKap Update for compressible case, solvatin not working, as
 y_n not coming from expressions for y_alpha

---
 .../CompressibleValidationAnalyticalMethod.py |  45 ++++---
 ...ncompressibleValidationAnalyticalMethod.py |  10 +-
 .../ValidationSimplifiedSystem.py             | 125 +++++++++---------
 src/Helper_DoubleLayerCapacity.py             |  92 ++++++-------
 src/__pycache__/Eq02.cpython-312.pyc          | Bin 16390 -> 16369 bytes
 src/__pycache__/Eq04.cpython-312.pyc          | Bin 20414 -> 20393 bytes
 ...Helper_DoubleLayerCapacity.cpython-312.pyc | Bin 9831 -> 9644 bytes
 7 files changed, 129 insertions(+), 143 deletions(-)

diff --git a/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py b/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
index db4b3e5..992b84d 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/CompressibleValidationAnalyticalMethod.py
@@ -75,12 +75,13 @@ nR_m = nR_mol * NA * 1/(1e-3)# [1/m^3]
 pR = 1.01325 * 1e+5 # [Pa]
 LR = 20e-8
 chi = 80 # [-]
-
+K_vec = ['incompressible', 50_000, 15_000, 1_500, 500] # Working
+# ! Not working: 100_000, 20_000
 # Parameter and bcs for the electrolyte
 Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
 a2 = (pR)/(nR_m * k * T)
-K = 1_000#'incompressible'
 kappa = 0
+K = 50_000
 Molarity = 0.01
 y_R = Molarity / nR_mol
 z_A, z_C = -1.0, 1.0
@@ -103,26 +104,26 @@ phi_left_dimless = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
 
 # Numerical calculations
 # Solution vectors
-# y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
-# for i, phi_bcs in enumerate(phi_left_dimless):
-#     y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=1e-4, max_iter=2500, return_type='Vector', relax_param=relax_param)
-#     y_S_ = 1 - y_A_ - y_C_
-#     y_A_num.append(y_A_)
-#     y_C_num.append(y_C_)
-#     y_S_num.append(y_S_)
-#     phi_num.append(phi_)
-#     p_num.append(p_)
-#     x_num.append(x_)
+y_A_num, y_C_num, y_S_num, phi_num, p_num, x_num = [], [], [], [], [], []
+for i, phi_bcs in enumerate(phi_left_dimless):
+    y_A_, y_C_, phi_, p_, x_ = solve_System_4eq(phi_bcs, phi_right, p_right, z_A, z_C, y_R, y_R, K, Lambda2, a2, number_cells, solvation=kappa, refinement_style='hard_log', rtol=1e-4, max_iter=2500, return_type='Vector', relax_param=relax_param)
+    y_S_ = 1 - y_A_ - y_C_
+    y_A_num.append(y_A_)
+    y_C_num.append(y_C_)
+    y_S_num.append(y_S_)
+    phi_num.append(phi_)
+    p_num.append(p_)
+    x_num.append(x_)
     
-# Q_num = []
-# for j in range(len(phi_left_dimless)):
-#     Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
-# Q_num = np.array(Q_num)
+Q_num = []
+for j in range(len(phi_left_dimless)):
+    Q_num.append(Q_num_(y_A_num[j], y_C_num[j], n(p_num[j], K), x_num[j]))
+Q_num = np.array(Q_num)
 
-# dx_ = phi_left_dimless[1] - phi_left_dimless[0] # [1/V], Assumption: phi^L is uniformly distributed
-# C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
-# C_DL_num = np.array(C_DL_num)
-# C_dl_num = C_dl(Q_num, phi_left_dimless)
+dx_ = phi_left_dimless[1] - phi_left_dimless[0] # [1/V], Assumption: phi^L is uniformly distributed
+C_DL_num = (Q_num[1:] - Q_num[:-1])/dx_ # [µAs/cm³]
+C_DL_num = np.array(C_DL_num)
+C_dl_num = C_dl(Q_num, phi_left_dimless)
 
 # Analytical calculations
 # K = K_vec[0]
@@ -137,7 +138,7 @@ plt.figure()
 # plt.plot(phi_left, Q_num[1] - Q_ana, label='Difference')
 # plt.plot(phi_left_dimless, Q_num - Q_ana, label='Difference')
 # plt.plot(phi_left_dimless, Q_num[1], label='Numerical')
-# plt.plot(phi_left_dimless, Q_num, label='Numerical')
+plt.plot(phi_left_dimless, Q_num, label='Numerical')
 plt.plot(phi_left_dimless, Q_ana, label='Analytical')
 plt.grid()
 plt.legend()
@@ -150,7 +151,7 @@ plt.figure()
 # plt.plot(Phi_pot_center(phi_left), C_DL_num[1] - C_DL_ana, label='Difference')
 # plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num - C_DL_ana, label='Difference')
 # plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num[1], label='Numerical')
-# plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num, label='Numerical')
+plt.plot(Phi_pot_center(phi_left_dimless), C_DL_num, label='Numerical')
 plt.plot(Phi_pot_center(phi_left_dimless), C_DL_ana, label='Analytical')
 plt.grid()
 plt.legend()
diff --git a/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
index f23debe..cf633d7 100644
--- a/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
+++ b/examples/ReproducableCode/DoubleLayerCapacity/IncompressibleValidationAnalyticalMethod.py
@@ -44,7 +44,7 @@ chi = 80 # [-]
 Lambda2 = (k*T*epsilon0*(1+chi))/(e0**2 * nR_m * (LR)**2)
 a2 = (pR)/(nR_m * k * T)
 K = 'incompressible'
-kappa = 0
+kappa = 5
 Molarity = 0.01
 y_R = Molarity / nR_mol
 z_A, z_C = -1.0, 1.0
@@ -61,7 +61,7 @@ rtol = 1e-4 # ! Change back to 1e-8
 # phi^L domain
 Vol_start = 0.1 # ! Change back to 0
 Volt_end = 0.75
-n_Volts = 5#0
+n_Volts = 10#0
 
 phi_left = np.linspace(Vol_start, Volt_end, n_Volts) * e0/(k*T)
 
@@ -99,7 +99,7 @@ C_DL_ana = C_dl(Q_ana, phi_left)
 # Plotting
 plt.figure()
 # plt.plot(phi_left, Q_num - Q_ana, label='Difference')
-# plt.plot(phi_left, Q_num, label='Numerical')
+plt.plot(phi_left, Q_num, label='Numerical')
 plt.plot(phi_left, Q_ana, label='Analytical')
 plt.grid()
 plt.legend()
@@ -110,7 +110,7 @@ plt.show()
 
 plt.figure()
 # plt.plot(Phi_pot_center(phi_left), C_DL_num - C_DL_ana, label='Difference')
-# plt.plot(Phi_pot_center(phi_left), C_dl_num, label='Numerical')
+plt.plot(Phi_pot_center(phi_left), C_dl_num, label='Numerical')
 plt.plot(Phi_pot_center(phi_left), C_DL_ana, label='Analytical')
 plt.grid()
 plt.legend()
@@ -120,5 +120,5 @@ plt.tight_layout()
 plt.show()
 
 print('phi_left:', phi_left)
-print('Q_num:', Q_num)
+# print('Q_num:', Q_num)
 print('Q_ana:', Q_ana)
\ No newline at end of file
diff --git a/examples/ReproducableCode/ValidationSimplifiedSystem.py b/examples/ReproducableCode/ValidationSimplifiedSystem.py
index 9d5c47d..dc92552 100644
--- a/examples/ReproducableCode/ValidationSimplifiedSystem.py
+++ b/examples/ReproducableCode/ValidationSimplifiedSystem.py
@@ -22,7 +22,6 @@ del sys.path[0]
 # Further imports
 import matplotlib.pyplot as plt
 import numpy as np
-import pandas as pd
 
 # Define the parameters and boundary conditions
 phi_left = 4.0
@@ -36,59 +35,11 @@ Lambda2 = 8.553e-6
 a2 = 7.5412e-4
 solvation = 15
 number_cells = 1024#*4
-refinement_style = 'log'
+refinement_style = 'hard_log'
 rtol = 1e-8
 
-# # Incompressible
-# K = 'incompressible'
-
-# # solve the complete system
-# y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
-
-# # solve the simplified system
-# y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
-
-# # Evaluate the difference
-# plt.figure()
-# # plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
-# # plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
-# plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-# plt.figure()
-# # plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
-# # plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
-# plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-# plt.figure()
-# # plt.plot(x_4eq, phi_4eq, label='phi_4eq')
-# # plt.plot(x_2eq, phi_2eq, label='phi_2eq')
-# plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-# plt.figure()
-# # plt.plot(x_4eq, p_4eq, label='p_4eq')
-# # plt.plot(x_2eq, p_2eq, label='p_2eq')
-# plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
-# plt.grid()
-# plt.xlim(0, 0.05)
-# plt.legend()
-# plt.show()
-
-
-
-# Compressible
-K = 10_000
+# Incompressible
+K = 'incompressible'
 
 # solve the complete system
 y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
@@ -98,20 +49,20 @@ y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right,
 
 # Evaluate the difference
 plt.figure()
-plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
-plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
-# plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
+# plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
+# plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
+plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
 plt.show()
 
 plt.figure()
-plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
-plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
-# plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
+# plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
+# plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
+plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
 plt.show()
 
@@ -120,7 +71,7 @@ plt.figure()
 # plt.plot(x_2eq, phi_2eq, label='phi_2eq')
 plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
 plt.show()
 
@@ -129,6 +80,54 @@ plt.figure()
 # plt.plot(x_2eq, p_2eq, label='p_2eq')
 plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
 plt.grid()
-plt.xlim(0, 0.1)
+plt.xlim(0, 0.05)
 plt.legend()
-plt.show()
\ No newline at end of file
+plt.show()
+
+
+
+# # Compressible
+# K = 10_000
+
+# # solve the complete system
+# y_A_4eq, y_C_4eq, phi_4eq, p_4eq, x_4eq = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# # solve the simplified system
+# y_A_2eq, y_C_2eq, phi_2eq, p_2eq, x_2eq = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=0.05, refinement_style='hard_log', return_type='Vector', max_iter=1_000, rtol=rtol)
+
+# # Evaluate the difference
+# plt.figure()
+# plt.plot(x_4eq, y_A_4eq, label='y_A_4eq')
+# plt.plot(x_2eq, y_A_2eq, label='y_A_2eq')
+# # plt.plot(x_4eq, y_A_4eq - y_A_2eq, label='y_A_4eq - y_A_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# plt.plot(x_4eq, y_C_4eq, label='y_C_4eq')
+# plt.plot(x_2eq, y_C_2eq, label='y_C_2eq')
+# # plt.plot(x_4eq, y_C_4eq - y_C_2eq, label='y_C_4eq - y_C_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# # plt.plot(x_4eq, phi_4eq, label='phi_4eq')
+# # plt.plot(x_2eq, phi_2eq, label='phi_2eq')
+# plt.plot(x_4eq, phi_4eq - phi_2eq, label='phi_4eq - phi_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
+
+# plt.figure()
+# # plt.plot(x_4eq, p_4eq, label='p_4eq')
+# # plt.plot(x_2eq, p_2eq, label='p_2eq')
+# plt.plot(x_4eq, p_4eq - p_2eq, label='p_4eq - p_2eq')
+# plt.grid()
+# plt.xlim(0, 0.1)
+# plt.legend()
+# plt.show()
\ No newline at end of file
diff --git a/src/Helper_DoubleLayerCapacity.py b/src/Helper_DoubleLayerCapacity.py
index 7714945..6fa7ab4 100644
--- a/src/Helper_DoubleLayerCapacity.py
+++ b/src/Helper_DoubleLayerCapacity.py
@@ -4,6 +4,7 @@ Jan.Habscheid@rwth-aachen.de
 '''
 
 import numpy as np
+from scipy.optimize import fixed_point, fsolve
 
 
 # Helper functions
@@ -157,66 +158,51 @@ def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float
         # Assume E_p = p_R = 0        
         D_A = y_A_R / (np.exp(-(solvation+1)*a2*p_R-z_A*phi_R))
         D_C = y_C_R / (np.exp(-(solvation+1)*a2*p_R-z_C*phi_R))
-        D_N = y_N_R / (np.exp(-(solvation+1)*a2*p_R-z_N*phi_R))
-
-        D_tilde_A = np.log(D_A)
-        D_tilde_C = np.log(D_C)
-        D_tilde_N = np.log(D_N)
-
-        CLambda_L = - z_A * phi_L + D_tilde_A - z_C * phi_L + D_tilde_C - z_N * phi_L + D_tilde_N
-        CLambda_R = - z_A * phi_R + D_tilde_A - z_C * phi_R + D_tilde_C - z_N * phi_R + D_tilde_N
-
-        Lambda = np.sqrt(Lambda2)
-        a = np.sqrt(a2)
-
-        N = 3
-
-        Q_DL = Lambda * a * (np.sqrt(-CLambda_L/(-(N-1)*(solvation+1)-1)) - np.sqrt(-CLambda_R/(-(N-1)*(solvation+1)-1)))
-
-        # CLambda_L = np.log(D_A * np.exp(-z_A * phi_L) + D_C * np.exp(-z_C * phi_L) + D_N * np.exp(-z_N * phi_L))
-        # CLambda_R = np.log(D_A * np.exp(-z_A * phi_R) + D_C * np.exp(-z_C * phi_R) + D_N * np.exp(-z_N * phi_R))
-
+        D_N = y_N_R / (np.exp(-a2*p_R))
+
+        Q_DL = []
+        dx_phi_L = 0
+        for phi_L_ in phi_L:
+            def func(p_L):
+                # A_Term = D_A * np.exp(-z_A*phi_L)
+                A_Term = D_A * np.exp(-(solvation+1)*a2*p_L-z_A*phi_L_)
+                C_Term = D_C * np.exp(-(solvation+1)*a2*p_L-z_C*phi_L_)
+                N_Term = D_N * np.exp(-a2*p_L)
+
+                return (A_Term + C_Term + N_Term) - 1
+            p_L = fsolve(func, dx_phi_L)
+            dx_phi_L = np.sqrt(p_L - E_p) * np.sqrt(2*a2/Lambda2)
+            Q_DL_ = Lambda2 * dx_phi_L
+            Q_DL.append(Q_DL_)
+        Q_DL = np.array(Q_DL)
+
+        # def func_L(p_L):
+        #     CLambda = D_A * np.exp(-z_A * phi_L - solvation*a2*p_L) + D_C * np.exp(-z_C * phi_L - solvation*a2*p_L) + D_N
+        #     return np.log(CLambda)/a2
+        
+        # p_L = fixed_point(func_L, 0, maxiter=5_000)
         # Lambda = np.sqrt(Lambda2)
-        # Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2) * (np.sqrt(CLambda_L/(solvation+1) - E_p * a2) - np.sqrt(CLambda_R/(solvation+1) - E_p * a2))
-    else:
-        C_A = y_A_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_A*phi_R))
-        C_C = y_C_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_C*phi_R))
-        C_N = y_N_R / ((K + p_R  - 1)**(-(solvation+1)*a2*K)*np.exp(-z_N*phi_R))
-
-        # # Lambda_tilda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
-        # # Lambda_tilda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
-
-        # # # Lambda_hat_L = np.exp(np.log(Lambda_tilda_L) / (a2 * K))
-        # # # Lambda_hat_R = np.exp(np.log(Lambda_tilda_R) / (a2 * K))
-
-        # # dx_Phi_L = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_L))
-        # # dx_Phi_R = np.sqrt(2 * a2/Lambda2 * (1 - K - E_p + Lambda_hat_R))
-
-        # # Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
-
-        # Lambda_tilde_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
-        # Lambda_tilde_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
-
-        # Lambda_hat_L = (1/Lambda_tilde_L)**(-1/((solvation+1)*a2*K))
-        # Lambda_hat_R = (1/Lambda_tilde_R)**(-1/((solvation+1)*a2*K))
-
         # a = np.sqrt(a2)
-        # Lambda = np.sqrt(Lambda2)
-        # K_tilde = 1 - K + E_p
+        # Q_DL = np.sqrt(2) * Lambda * a * np.sqrt(p_L - E_p)
 
-        # Q_DL = Lambda * a * np.sqrt(2) * (np.sqrt(Lambda_hat_L - K_tilde) - np.sqrt(Lambda_hat_R - K_tilde))
+        # ! ToDo: Implement full analytical for solvation == 0
+        
+    else:
+        C_A = y_A_R / ((K + p_R  - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
+        C_C = y_C_R / ((K + p_R  - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
+        C_N = y_N_R / ((K + p_R  - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
 
-        Lambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
-        Lambda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N * np.exp(-z_N*phi_R)
+        CLambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N
+        # CLambda_R = C_A * np.exp(-z_A*phi_R) + C_C * np.exp(-z_C*phi_R) + C_N
 
-        B = -a2 * K
-        A = 1/2 * Lambda2/a2
-        K_tilde = K + E_p - 1
+        Left = np.sqrt((1/CLambda_L)**(-1/(a2*K)) + 1 - K - E_p)
+        # Right = np.sqrt((1/CLambda_R)**(-1/(a2*K)) + 1 - K - E_p)
 
-        dx_Phi_L = (np.power(1/Lambda_L, B) - K_tilde) / A 
-        dx_Phi_R = (np.power(1/Lambda_R, B) - K_tilde) / A
+        Lambda = np.sqrt(Lambda2)
+        a = np.sqrt(a2)
 
-        Q_DL = Lambda2 * (dx_Phi_L - dx_Phi_R)
+        Q_DL = Lambda * a * np.sqrt(2) * (Left)# - Right)
+                     
     return Q_DL
 
 def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, solvation:float) -> float:
diff --git a/src/__pycache__/Eq02.cpython-312.pyc b/src/__pycache__/Eq02.cpython-312.pyc
index db7936e73aa933adfe164ec1a59346d55a066edb..22914ed1a5ac4beb0cd95fb993ee210cac42c52d 100644
GIT binary patch
delta 46
zcmZo`VEkCmb()u#i-CcGA!++YE>B*O0{x8q+*JLn#5{eM{N&Qy)Vz}7%{9DTMgUK^
B4>$k-

delta 67
zcmexZ-`2p$b()u#i-CcG;o;+rT%Np|t@^onCHl$wp~b01#rj!^c^QfNF8Rr&xv6<2
W#rgq7`RPT8xw)Bn>6_DeyNmz}3K#4E

diff --git a/src/__pycache__/Eq04.cpython-312.pyc b/src/__pycache__/Eq04.cpython-312.pyc
index 62a3b60dcbc1854782221fd1d152844cbb07b5d6..4fcbdc9745e29e883aed5afc77b94d4200e0387f 100644
GIT binary patch
delta 479
zcmdltpK;}UM()$Ryj%<n3=HLB$2W2d@ro4aXXNLm>Srb9>AU17m*%GCl@w2Q<24m_
z$uCXHN%cvrOf7OwEJ#ewEUDbw&-<K-kz;eYz$p&KjLlAxX^fLMNr_CpCTlWTL{4O~
zft=~)Sm`n;SvCd+24MyUhR<OP3=C75r!$l=g6st066VRzOl>4rLj{5vRx<e&i8C-T
ztYj+UVPIgGtY?;#XUM?7@S}m@j<EO~rWwYw*%v6SU|iw2l69fZ2F)Fb8+9)#Ib4)*
zIADB2_@L!Q37^Ztz883Wi$WO~7?N2)mOwFHXF%03)G$u|U~0q)cF<-y^IRsDm5fFF
zlP_9nZr)}o#mHE<`J&Z5M#jd;8*FqKn>OFI$!2D3ne6HCit8i;1A_$v14D7hW?RRZ
z%<O)Q-`E&A{U+aa>9?DxJRy08@dD-rftQ8U*NCi0UM+c<$9PN9MIMVge8Ll*XE4tQ
zUZA|7@Un>J8j}sot1U0{neJh_$Y=GDnUU9e^H$dn%v^#(>E_JO%~>-TC$IN3U~HTG
hz*CR0V=}+jX2#CRhrRwVo}7HZdmH1^&E-Ds%mAUhj)DLH

delta 402
zcmZ2EpK;%OM()$Ryj%<n3=9t+AKl0;#H-n+pPN^rpR6BRoLW?@pOu)Gk*M#IpIn-o
znpaY+A5fH^UX+-do0*qB*^t*%#LqcDuec<$q!c9Un^{p(T9mrElJ_|iqt50;fm0le
zRhxAr(-<f3k>Q%$AiH|AvUHggKMMl`gD?XF!)F%;28OB3lld&gC!aLc<Oyb2$>dig
z&cML1l5sM(SyHGD0|Ud428KJrq7#xQo6iVZV7Q=gq3O(+6_GmxSH@nHwZ159y(9U6
z^UlnR!tNJ%+>1gP7#Na4#)EtW!XP%1Rnr+t7$;vev)L?ap2fseR6O~(mFDI(mQsw2
z#hZ^?-D6~FU|?XF{Lxy4jg66kp{Q~5Wt(hf#^%ZP4zIXQFfcG!FfcF_$8R=voXO1P
z!T61hk<)|m3)f^VH<ii1Trwu_atqkJ-1Q|hm$Fa>C-VzV)(qjvvpfwLTPL6P)MIOB
lU|`T>pZw8IYO<u)6vmFp>%9Ijo|wGadmH1E%}GA)%m9+6e@*}Z

diff --git a/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc b/src/__pycache__/Helper_DoubleLayerCapacity.cpython-312.pyc
index 82d88e6411714334b2d2d3d89e81ae6b20e02a4c..7fcfd455a306d8842cf22799aa1f0ea7210cd3a3 100644
GIT binary patch
delta 3228
zcmaFvv&Nh6G%qg~0|NuYJI2##;<6L@B<e+&7#OBAq%cG=q%fv1<uFDur7)&2=P>88
zM6u+uMzJz7q_CtiW$}P?q2d%)?BZ!GDQqn)QS5m1aip=Nu(z;8ai*|xfDD;vCtgpA
z4Qwf#Ei6%7DO@1UsICI>Qn*vOQ#q>G85nT7ku8O%g(Zq7MJa_BVulz<BQj3m14&Gd
zW5n(>0gx=RdI&p3FqI=kD3v3HHI*ZcFGaY8HHtSyB$XpYG?gPoER_Q;FP_SgB9Y3G
zBALp8<Xov#juh!sjuaWF2~sJtEi55Xd}J6Qm&%bMpURP<fMlFPied{(6n`a?rt)M1
zCKGW@rd!-;nH8xi@df#rc_p{l(u(tQ%2GF{FtswuvobI+6f-g~Fev<L)6dAyP1Vmz
z%+q(tPcF?(%_}LMEX&fu%E!RKz`c1H%N0gX4h9B>A|?>Q4I)6%3DI4|%fP^Jiz%g|
zNRWYnK^kPLI0FMi1H%Ue1|H!FhMiV3*cUL)<hsbMbb&>wh=228HZMj~ka2!RoD2*M
znv7tLnoLDPAcMiq0@<UW05O@xIX)$)NQ8lbL3gqfhb|-k<^qnnjEo$U`MAW{L_u1_
zCTnxa@`9bgm{%mjz`$TXIh4y<fFI;A=E``-Tg;X5&P6ht7jRWFX-I-ZSxbr%3#v3Z
ziey0s$bkruz9M-LO94bEPIlzg76Kc}78svbnj2rF#=yW3GP#b|j*)-zUfuvhkn5Ql
z7#KcRurM&RGq<y}bG37~OSMb4%e2c*XGmeJ<*8w><*ngJVanEEU?}=9*?>>Mpp!d=
zxe64n44qt6tPBh_>>#p}xrPlyS8*{gtOmsj0|P@9Hv>bL)Z_--npkjXlAe5lPa&$4
ziyLY>iitAm47GeIjJ5nVtS$_(DGUr<vK>+#(j77#Tph9<+#SpvENL7mEG-<h0wti>
z0H^#Kfg1kRplD%WV5k+G%q1YtDp<o?#Wi`AfJQ_LE7U8TASc5uV`X5d<*(ta;RkuF
zR;Y$A3lw5t;}ArR5KKKILr)qDLkZOL3=COZlXduICHQLif*EQUYB*}xLBW>73<{zY
zw$91#1$CJ^nJ3E$Ns0*9FyfGB0n3ZjFkzGTU>2<ro?IX#2Qp`ZkO)(a2#CknIk{g*
zib=F;vc7<_SqgKF7|2^GY@J*+93ZcNcvTz>43K!N;Q%=Z6c&h>s}-NjC7{VDF_}w1
zPNb6y=EGXa8i^W7kP#pWnG{x-1W2ex97NWFy{pN-*-ylqMN<Nl5<yud86GMO3=C`x
z3=G1cboYsYfnh4gbcRkAkh6_CIa3%xLB!3_2~JZr3?LIbxjG@PMow9vsA246#ixxM
zhc>28HjuXZPIhhvhz-cTsAaBUtYHGVq6VDGQy4*^6U?B=<W~gph$d4JI0G`LRupLR
z++r?>_qoO97$1^alzWTK8A|&>Xl55s)^&+@zQybk?^mP;%D8L=8JY1u@kRReprC=|
zR8Wz_11iWA6cieu09tmWmF6W^$pyx{_{67V=H{dp7sn^&CF<Ga<R>TQ6x->6<d=c;
zw=*y>{Agf!ASpAaXnxtuvK5>g6jyRD$lDNfQO$e@^F=l59nJ^Xc6wb@bGt0>eo@lB
zegfMKq7@^lT~M?jd4=-A^x5?nCG{q-J&=%FAU2bKh1x|4gAI-s#Z4~=m==Kwt|CxH
z4k;knz(H1I1oANpC=ghR3ybPYZm}g66r|>*++s~EDoU)>lq*sPsn7rsnjk^}L@0s?
zEfAp%B6L6ms4^_l1+hS3TBHv$lesG1?-sLbd;z!^zQvJJ5f2XPB5RO3V~`S7a5&y#
zb_PXDtTQNHoa6m&aXR}X<|d^i#{1l2@kvc9xy1(I-(pNGG6!k401=iT!U`1Bpwd=B
zK_MYQ0R*5FhzZ6H4h{}QyptbDim7;kw19N@GBGgxZ(#Tk#K^$S)8KrAi?_k~20O<E
z9+^v%C8d(<XB5pyUZ6ZPeRBOpKJ^af2mFE^EH{KjIyi2KNq=JC=9KSn`_9D1YxbFe
zjo0i8AA_*u43QZXmxVMsobQNBPA{KWKBIU+@XX4~;+kum*1NBC-{80*X`{;y@r^#0
zbsR5>J9c>75RsmdaZyCA!}SJ_K)-jV_YBsnJaQkzxjE&()Jrf32=!NVR?J}OuIpg=
zV8y^A(BS<+o`GLv2ICCpnd}P;XY#CIT4A)3bw%Pz_KV^MTZFcYZ4}!fyn|&&=uWQ7
zhRzqooiFmabTHrG6DD0tJ;{0=h=|Wony)rfZGrLvkp-y>r7lb9Toln+QM@H+d)UUX
z9l<+X4v6mbx@_uxQQ!Tdi2DiUiz2=qoL|^L9xd-IpTRgob0+tKq7_Cf5?7i;EUdgJ
zZU|1jI~XsDTkdeaC~kk5-=Tx$3mb!w$aL|E;xoiA3o3T7-QnTuckgtc;MDEi!T5k*
zbOvK*)f}h!?lavNIL`LE$gj3S<s!d+2g?I){s!+N{>|00d91;pWS}WPZsPU^B_baX
z;R_-_X}w5jvW>E&1h{VJ0jKGd%v?~pP~^|Rz;I!5u5xreC<|&br9<ke5Bp&Zdrg5`
z?D6p_`N{F|w^;K^a|<fLne!H3adKusrCxqPNoH<lRcaBaT~Gwhtn7Iy;QUzxYE0c?
zP0PtoECIK8z-^ja%*7=|;5JJUsCiKY@@$bmgEa#~5xAWNatF99`^90Co1apelWJEK
z#=yV;st=0!nHU&8Ff%eTeq>^1Wckd+z$kp5L2q)Z$_bYkMn-`drXLtU^a@r8y&)7r
zAAs^tK;<t)#6l!4q-H@G1*IPtVkQ@=idlmU_`ty6%M`=}k^8{F5Y7~X;AAl6A~;1%
U<xHQ!mVT&VntVugAq&U@0QQ6VBLDyZ

delta 3354
zcmZ4E{oIG|G%qg~0|Ns?*o-4-@`@ArB<fk17#OBAq%cG=q%fv1<uFDur7)*5WwC)&
zq2d%4?BZ!GDXc9lQ7m}$v8J)4u(hy6v8Ax1nNbh51CM5`j!9!n;b>urVo%{jGZ1DT
zS1Ly;YZW^K15RhNrEs^fC`ECmD5mh>FoYK*KG}{@ydGP4@PlNb_90Oz0;#Mif~l-2
zEUBz%d?`XLtWjJk!l|q&BB`t?qN%KKd9hU16!BEn6p2*U6ohjnQ(04_Qdv`^p(aSB
z$h5G8L~)a0glsBnid-sdiae5W3MmRLEKxj_OqxoQ9hgiucQdsznlLjkFcdQ~Ffb_m
znysIkSE8S+A6lGRRIHzsn3s{L?~<Qfnwy$eQmh|Pl%HOdn46oKmp)marG=G;fq{W*
z^A46Pj2fU2Dq><_U?}1Q5ujuSF|vr8fq~%`Q%Xe<KLZ1U)Z{H}3XHs)ud;bDnu7HB
z6|plgFlaJ@Wi^?K1VFmM4gwjbpa9X%;vAonQzXQ|z@RfZghQ8+cXJcRTt-H=$qHQJ
zY{DQdB9m>oWO>1MG3FIXF)%RLO-|>s7T^Wh&RiMqc#F9*-nmF>^CqrJCJk|rC~HYk
zVnLNAN0Bth02vSg(pMx4V#$FB`N_e&+CpGM*#hJ9N^|3jR2Uc-f+tVlwPWO+e4aNz
znS+6W0p#4zMr`1m#8}Hy!(Pi<!;!+2t--)Bfw3rWasZ!#K@E2ba}@&v12;nrR~0J*
zLk&BaWUgTY(N$aw468vgz`(#z#m&G_%QtxeZcQvWH1SXVz^7n?!zLaangmLq_A)SJ
zu}tO@kgXR)7C}ncP-YE3%%!zL1XK#ZRSMTIb2HSig8W(|2nvT976=dI^BN(r2oFTW
z2`nTG5~>veYXb)nTpO4NRs++95P|3toh%@v$tX5iKuBJk3+@E*6qZ_v8nGI2kRw2H
ziXz+$H4-VzAhWm`V8I5Gg@{#gPwwN_2(1w%qzozUxEL5f2GmN{NCq?1Fw}6=u!B4T
zF%V`2H$%<j`2xC(HIoYkq-3OO7zrsmD<CB!UBiS|*<=ktX=a%ksmU8T<i+8E4f8$7
zQ8m&axyiEyr5I&EEJmKma|9(=YGomsCKn2bGC|a^vT!rhuuWbdC<;+ISyNC9R8-f>
z)o^h$)Np`ItCde-t>FZ*YPi5jR-uNYh8N-4T1BM90m>!N<O8Epm}_J~`fB7?gK|0p
z14E4>C?;wYK%%uu=(^ot8X<HFbB!EGAr9@zHJmj{HOe3ZK~99khbCL8!)86%V5ZH#
zMZPiCi-EE@s0?}u%E0Lmb064ivKE1h2j<j@f?Lcv`RTV<iVKTMZm|~Rm!}qKsuZb$
z^r?Xebr1n5DL^a@5TOYov_OOih|mTRIv_$9ROT>O#rxf2c8xE%#a!<KDwJH}oo_L_
z#QWXibcrv?%t=X&cSJCq5lp{ZoX$Rpxk)LB@jeJ<&@DEI;4Q|)Ta13-((@LxGss|P
zkipLJez!OvD&l?Mte{(rPPZ669dEIFL+rf8kx~&KkdX;e17!sjgGxFDB0v%MWJ4)2
zWp4%shC)z9QpU`{(7^B^h>?Mtr@{FK7jJ{}4R($TJTjLjXG<kdE|XHOe<C6_-EE@V
z4CBe(9h`S?avq3DPH>&!I?;E5=S30S4a^%-H*#MTvEGq+QN;a(@`ZrVliC+WLN7!l
zUKB~};JiV!>Uxm!CxXJ$#V3l-5Wg&_*ui#(hp*qg(|v+dw|58Q10LQE#s~bOGZ;Io
z<^;_TpBcU&b3xPw!wrEOjW>wxkUU_xQ~IKz>w(0JhVB<tJub_7o^d+we$xFyXxxR6
zxC=?i7gAC$CZ%3TNV||xc%i8HVn*@Bgpv!TWfv35E{2p}uJ@_9C|hxnzp{hnLkuGW
zk3fU>2VX`8evuiBGn{9#FEE_RvqE!&<4Wxnd>e``YFO@2yr^M)QQT%v(EhNUVF!W_
zxSS9@=yloF|Dw46MZSOz<{Ny%q-d?DKrgb^4@L|k;&X)Ni_H{UAiRKOLFhuR%aYm`
zg|$~CZ*ki0w$W{e^A3{(tUIkPo48%nbGs<)c0%|}(fP9alVumeqb`I+U5HM&7?yC^
zKk=ex;zi-a4vrhbVq}<y)3^t`{2k5{6g%BJ7;o_KPf+jlouPD*M<%g@@rID_gyf0p
zGaM&scd&z!qsK&#8G$n@W=5?Lyey)<BK@L>=?=k*!d4v|p!A7b_5r`x42=a57o~JB
z^6Pc5+~DSK@Gj!ryjUiWwH}luGzG{lFML3ShA)Wl0}-I&1DwB$K}{s|d<khL@dU=Z
z_{67V=Ef)HB^I#;FfcIOn><x1nvr+%M<t0U-dpVP@hSPq@$t7<^Gb6IDj@|OdtORn
zQBh)L5vVzLi#07LKd}Vd90NDFZZQ{^6oDI5;I>Z@$WuiilOYW~Q2VNgdvdL^xG)bB
z1H%VqMn=YuOw5cdpSc(qg(t64K2e{|$S5$w^aBHkUcm~XH-tjy15o}6sQiV9Sct@h
z)GR2YpcKletogveP{-)Z^pSy~iqV+~BK(1Y!Ivor!3k%IL2xpda+$>7>WY}knLdM^
P@ga(NGP~+R7Le-!?FKA9

-- 
GitLab