diff --git a/@BearImp/calcCap_puchtler2025.m b/@BearImp/calcCap_puchtler2025.m
index 0fea296884379c5d79ace6e41ce5ed07944fa938..511b75812d9daaf1284d37090e101b146f58efce 100644
--- a/@BearImp/calcCap_puchtler2025.m
+++ b/@BearImp/calcCap_puchtler2025.m
@@ -1,80 +1,80 @@
-function C_out = calcCap_puchtler2025(s,R_RE,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha,h_0,a,b,options)
-
- arguments
- s
- R_RE
- R_R
- R_ZL
- R_ZR
- B
- R_L
- epsilon_r
- alpha
- h_0
- a = 0
- b = 0
- options.CmethodDeformedArea = 'neglect'
- options.r2smallCriterion = 'r<1.00125*R_RE'
- end
-
- Delta_x = -(R_R - R_RE - s)*sin(alpha);
- Delta_z = (R_L - R_R) + (R_R - R_RE - s)*cos(alpha);
- B_RhR = sqrt(R_R^2 - (R_R - R_L + R_ZR)^2);
- B_RhL = sqrt(R_R^2 - (R_R - R_L + R_ZL)^2);
-
- r_groove = @(phi,theta) r_easy(phi,theta,R_RE,R_R,R_L,Delta_x,Delta_z,h_0,options);
- r_rimL = @(phi,theta) (-Delta_z*cos(phi) + sign(R_L)*sqrt(R_ZL^2 - Delta_z^2*sin(phi).^2))./cos(theta);
- r_rimR = @(phi,theta) (-Delta_z*cos(phi) + sign(R_L)*sqrt(R_ZR^2 - Delta_z^2*sin(phi).^2))./cos(theta);
-
- Theta_0l = @(phi) atan((R_RE*sin(alpha))./(-Delta_z*cos(phi) + sign(R_L)*sqrt(R_L^2 -Delta_z^2*sin(phi).^2))) - a*sqrt(max(1/R_RE^2-phi.^2/b^2,0));
- Theta_0r = @(phi) atan((R_RE*sin(alpha))./(-Delta_z*cos(phi) + sign(R_L)*sqrt(R_L^2 -Delta_z^2*sin(phi).^2))) + a*sqrt(max(1/R_RE^2-phi.^2/b^2,0));
- Theta_1l = @(phi) atan((+B_RhL - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZL^2-Delta_z^2*sin(phi).^2)));
- Theta_1r = @(phi) atan((-B_RhR - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZR^2-Delta_z^2*sin(phi).^2)));
- Theta_2l = @(phi) atan((+B/2 - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZL^2-Delta_z^2*sin(phi).^2)));
- Theta_2r = @(phi) atan((-B/2 - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZR^2-Delta_z^2*sin(phi).^2)));
-
- if R_L > 0
- phi_1 = pi/2;
- else
- phi_1 = asin(R_L/(R_L-s-R_RE))*0.9;
- end
-
- fun_groove = @(phi,theta) R_RE^2./(r_groove(phi,theta)-R_RE) .* cos(theta);
- fun_rimR = @(phi,theta) R_RE^2./(r_rimR (phi,theta)-R_RE) .* cos(theta);
- fun_rimL = @(phi,theta) R_RE^2./(r_rimL (phi,theta)-R_RE) .* cos(theta);
- absTol = 12.5e-4; % 2.5e-4 * 2 corresponds to 0.01 pF with a permittivity of 2.2
-
- C_grooveR = 2 * 8.854187e-12 * epsilon_r * integral2(fun_groove,0,phi_1,Theta_0r,Theta_1r,'AbsTol',absTol);
- C_grooveL = 2 * 8.854187e-12 * epsilon_r * integral2(fun_groove,0,phi_1,Theta_1l,Theta_0l,'AbsTol',absTol);
- C_rimR = 2 * 8.854187e-12 * epsilon_r * integral2(fun_rimR ,0,phi_1,Theta_1r,Theta_2r,'AbsTol',absTol);
- C_rimL = 2 * 8.854187e-12 * epsilon_r * integral2(fun_rimL ,0,phi_1,Theta_2l,Theta_1l,'AbsTol',absTol);
- C_out = C_grooveR + C_grooveL + C_rimR + C_rimL;
-end
-
-function r = r_easy(phi,theta,R_RE,R_R,R_L,Delta_x,Delta_z,h_0,options)
- assert(all(size(phi)==size(theta)))
-
- r = nan(size(phi));
- R_L_R_R = R_L-R_R; % pre calculate for efficiency
- R_R_inv_sq = 1/R_R^2;
-
- for ii = 1:numel(phi)
- thisPhi = phi(ii);
- thisTheta = theta(ii);
-
- fun = @(r) real((r*sin(thisPhi)*cos(thisTheta)./(R_L_R_R+R_R*sqrt(1-R_R_inv_sq*(r*sin(thisTheta)-Delta_x).^2))).^2 + ((r*cos(thisPhi)*cos(thisTheta)+Delta_z)./(R_L_R_R+R_R*sqrt(1-R_R_inv_sq*(r*sin(thisTheta)-Delta_x).^2))).^2 - 1); % real for fzero to not interrupt at imaginary values. Results are real anyway
- r(ii) = fzero(fun,1.2*R_RE);
- end
- switch options.r2smallCriterion
- case 'r<1.00125*R_RE'
- r2small = r(:)/R_RE<1.00125;
- case 'r<h_0+R_RE'
- r2small = r(:)<h_0+R_RE;
- end
- switch options.CmethodDeformedArea
- case 'neglect'
- r(r2small) = inf;
- case 'filmThickness'
- r(r2small) = h_0+R_RE;
- end
+function C_out = calcCap_puchtler2025(s,R_RE,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha,h_0,a,b,options)
+
+ arguments
+ s
+ R_RE
+ R_R
+ R_ZL
+ R_ZR
+ B
+ R_L
+ epsilon_r
+ alpha
+ h_0
+ a = 0
+ b = 0
+ options.CmethodDeformedArea = 'neglect'
+ options.r2smallCriterion = 'r<1.00125*R_RE'
+ end
+
+ Delta_x = -(R_R - R_RE - s)*sin(alpha);
+ Delta_z = (R_L - R_R) + (R_R - R_RE - s)*cos(alpha);
+ B_RhR = sqrt(R_R^2 - (R_R - R_L + R_ZR)^2);
+ B_RhL = sqrt(R_R^2 - (R_R - R_L + R_ZL)^2);
+
+ r_groove = @(phi,theta) r_easy(phi,theta,R_RE,R_R,R_L,Delta_x,Delta_z,h_0,options);
+ r_rimL = @(phi,theta) (-Delta_z*cos(phi) + sign(R_L)*sqrt(R_ZL^2 - Delta_z^2*sin(phi).^2))./cos(theta);
+ r_rimR = @(phi,theta) (-Delta_z*cos(phi) + sign(R_L)*sqrt(R_ZR^2 - Delta_z^2*sin(phi).^2))./cos(theta);
+
+ Theta_0l = @(phi) atan((R_RE*sin(alpha))./(-Delta_z*cos(phi) + sign(R_L)*sqrt(R_L^2 -Delta_z^2*sin(phi).^2))) - a*sqrt(max(1/R_RE^2-phi.^2/b^2,0));
+ Theta_0r = @(phi) atan((R_RE*sin(alpha))./(-Delta_z*cos(phi) + sign(R_L)*sqrt(R_L^2 -Delta_z^2*sin(phi).^2))) + a*sqrt(max(1/R_RE^2-phi.^2/b^2,0));
+ Theta_1l = @(phi) atan((+B_RhL - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZL^2-Delta_z^2*sin(phi).^2)));
+ Theta_1r = @(phi) atan((-B_RhR - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZR^2-Delta_z^2*sin(phi).^2)));
+ Theta_2l = @(phi) atan((+B/2 - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZL^2-Delta_z^2*sin(phi).^2)));
+ Theta_2r = @(phi) atan((-B/2 - Delta_x)./(Delta_z*cos(phi) - sign(R_L)*sqrt(R_ZR^2-Delta_z^2*sin(phi).^2)));
+
+ if R_L > 0
+ phi_1 = pi/2;
+ else
+ phi_1 = asin(R_L/(R_L-s-R_RE))*0.9;
+ end
+
+ fun_groove = @(phi,theta) R_RE^2./(r_groove(phi,theta)-R_RE) .* cos(theta);
+ fun_rimR = @(phi,theta) R_RE^2./(r_rimR (phi,theta)-R_RE) .* cos(theta);
+ fun_rimL = @(phi,theta) R_RE^2./(r_rimL (phi,theta)-R_RE) .* cos(theta);
+ absTol = 25e-4; % 2.5e-4 * 2 corresponds to 0.01 pF with a permittivity of 2.2
+
+ C_grooveR = 2 * 8.854187e-12 * epsilon_r * integral2(fun_groove,0,phi_1,Theta_0r,Theta_1r,'AbsTol',absTol);
+ C_grooveL = 2 * 8.854187e-12 * epsilon_r * integral2(fun_groove,0,phi_1,Theta_1l,Theta_0l,'AbsTol',absTol);
+ C_rimR = 2 * 8.854187e-12 * epsilon_r * integral2(fun_rimR ,0,phi_1,Theta_1r,Theta_2r,'AbsTol',absTol);
+ C_rimL = 2 * 8.854187e-12 * epsilon_r * integral2(fun_rimL ,0,phi_1,Theta_2l,Theta_1l,'AbsTol',absTol);
+ C_out = C_grooveR + C_grooveL + C_rimR + C_rimL;
+end
+
+function r = r_easy(phi,theta,R_RE,R_R,R_L,Delta_x,Delta_z,h_0,options)
+ assert(all(size(phi)==size(theta)))
+
+ r = nan(size(phi));
+ R_L_R_R = R_L-R_R; % pre calculate for efficiency
+ R_R_inv_sq = 1/R_R^2;
+
+ for ii = 1:numel(phi)
+ thisPhi = phi(ii);
+ thisTheta = theta(ii);
+
+ fun = @(r) real((r*sin(thisPhi)*cos(thisTheta)./(R_L_R_R+R_R*sqrt(1-R_R_inv_sq*(r*sin(thisTheta)-Delta_x).^2))).^2 + ((r*cos(thisPhi)*cos(thisTheta)+Delta_z)./(R_L_R_R+R_R*sqrt(1-R_R_inv_sq*(r*sin(thisTheta)-Delta_x).^2))).^2 - 1); % real for fzero to not interrupt at imaginary values. Results are real anyway
+ r(ii) = fzero(fun,1.2*R_RE);
+ end
+ switch options.r2smallCriterion
+ case 'r<1.00125*R_RE'
+ r2small = r(:)/R_RE<1.00125;
+ case 'r<h_0+R_RE'
+ r2small = r(:)<h_0+R_RE;
+ end
+ switch options.CmethodDeformedArea
+ case 'neglect'
+ r(r2small) = inf;
+ case 'filmThickness'
+ r(r2small) = h_0+R_RE;
+ end
end
\ No newline at end of file
diff --git a/@BearImp/calcGeo.m b/@BearImp/calcGeo.m
index 412b70513026e8a660bcbdf86b4135888d3d5e84..2034a93178273159a13d6bdf7463798782e5d239 100644
--- a/@BearImp/calcGeo.m
+++ b/@BearImp/calcGeo.m
@@ -79,8 +79,8 @@ elseif G.method.alphaForRaceRadius == 'estimate'
if abs(alpha-alpha_next) <= 100*eps(alpha)
break
elseif numberOfIterations >= 1000
- warning('Estimation of mounted contact angle did not converge')
- break
+ error('Estimation of mounted contact angle did not converge')
+ %break
else
alpha = alpha_next;
numberOfIterations = numberOfIterations + 1;