diff --git a/@BearImp/BearImp.m b/@BearImp/BearImp.m
index d478834aebb9ba9b14a351feee7285c8afe82c80..4119c6503bf8e0ef8ade124d214f8c13a0ac5703 100644
--- a/@BearImp/BearImp.m
+++ b/@BearImp/BearImp.m
@@ -100,6 +100,10 @@ classdef BearImp < handle & matlab.mixin.Copyable
         value = get(obj,name)
         obj = copyToArray(obj,varargin)
     end
+
+    methods (Access = public, Static)
+        C_out = calcCap_puchtler2025(s,R_WK,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha,a,b)
+    end
     
     methods
     % set-Methoden
diff --git a/@BearImp/calcCap.m b/@BearImp/calcCap.m
index 94b55a6acefe7e2a98385bfa8f4ff0f67a50f9e6..833f17b659ff5be4ff8a3c5f1648f58328f74c1e 100644
--- a/@BearImp/calcCap.m
+++ b/@BearImp/calcCap.m
@@ -127,6 +127,8 @@ for posBall_conductive=1:L.numberOfConductiveBalls
                 for recentVec=length(vecStruct)
                     runSemianalytisch3D(vecStruct{recentVec})
                 end
+            case 'Puchtler2025'
+                runPuchtler2025(find(~contactInd))
         end
     end
 
@@ -152,6 +154,8 @@ for posBall_conductive=1:L.numberOfConductiveBalls
             for recentVec=1:length(vecStruct)
                 runSemianalytisch3D(vecStruct{recentVec})
             end
+        case 'Puchtler2025'
+            runPuchtler2025(find(contactInd))
     end
 
     if C.method.k_vh_factor == 'on'
@@ -271,6 +275,39 @@ function tobias_steffen(indices)
     end
 end
 
+function runPuchtler2025(indices)
+    if isfield(C,'C_out')
+        C_semi_i = C.C_out(1,:,posBall_conductive);
+        C_semi_o = C.C_out(2,:,posBall_conductive);
+    else
+        C_semi_i = nan( numel(temp_s)/2, 1);
+        C_semi_o = nan( numel(temp_s)/2, 1);
+    end
+    temp_alpha = B.alpha(:,:,posBall_conductive);
+    %                                   calcCap_puchtler2025(     s           ,     R_RE,     R_R ,   R_ZL   ,   R_ZR   ,   B,    R_L ,   epsilon_r  ,     alpha           ,     h_0           ,     a           ,     b           )
+    tmp_puchtler_i = @(runSemi) BearImp.calcCap_puchtler2025(temp_s(1,runSemi), tempR_RE, tempR_li,-G.d_bLi/2,-G.d_bRi/2, L.B, -G.R_bi, S.epsilon_Oel,temp_alpha(runSemi),temp_h_0(1,runSemi),temp_a(1,runSemi),temp_b(1,runSemi));
+    tmp_puchtler_o = @(runSemi) BearImp.calcCap_puchtler2025(temp_s(2,runSemi), tempR_RE, tempR_la, G.D_bLa/2, G.D_bRa/2, L.B,  G.R_ba, S.epsilon_Oel,temp_alpha(runSemi),temp_h_0(2,runSemi),temp_a(2,runSemi),temp_b(2,runSemi));
+
+    if AddOn.Parallel_Computing_Toolbox && (numel(indices) > 20 || ~isempty(gcp('nocreate'))) % Pool nur starten, wenn mehr als 20 Berechnungen nötig
+        tempC_i = nan(1,numel(indices));
+        tempC_o = nan(1,numel(indices));
+        parfor ii = 1:numel(indices)
+            tempC_i(ii) = tmp_puchtler_i(indices(ii));
+            tempC_o(ii) = tmp_puchtler_o(indices(ii));
+        end
+        C_semi_i(indices) = tempC_i;
+        C_semi_o(indices) = tempC_o;
+    else
+        for runSemi = indices
+            C_semi_i(runSemi) = tmp_puchtler_i(runSemi);
+            C_semi_o(runSemi) = tmp_puchtler_o(runSemi);
+        end
+    end
+
+    C.C_out(1,:,posBall_conductive) = C_semi_i;
+    C.C_out(2,:,posBall_conductive) = C_semi_o;
+end
+
 function runSemianalytisch3D(indices)
     
     if isfield(C,'C_out')
diff --git a/@BearImp/calcCap_puchtler2025.m b/@BearImp/calcCap_puchtler2025.m
new file mode 100644
index 0000000000000000000000000000000000000000..65ab72ce719468064c5d3227e29e9ca98a17b321
--- /dev/null
+++ b/@BearImp/calcCap_puchtler2025.m
@@ -0,0 +1,68 @@
+function C_out = calcCap_puchtler2025(s,R_RE,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha,h_0,a,b)
+    
+    arguments
+        s
+        R_RE
+        R_R
+        R_ZL
+        R_ZR
+        B
+        R_L
+        epsilon_r
+        alpha
+        h_0
+        a = 0
+        b = 0
+    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);
+    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 = 2.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)
+    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
+    r(r<h_0+R_RE) = h_0+R_RE;   %% TODO: oder = inf? %%
+    %r(r<R_RE*1.00125) = inf;
+end
\ No newline at end of file
diff --git a/@BearImp/calcGeo.m b/@BearImp/calcGeo.m
index a067fb8a551e7913a3f13f970b58b9f2ae31723b..4a97acca66cf2e2e3659758f7d32d9148dd8e4f9 100644
--- a/@BearImp/calcGeo.m
+++ b/@BearImp/calcGeo.m
@@ -38,6 +38,9 @@ G.R_RE = G.D_RE/2;
 
 G.d_m = (G.d + G.D)/2;
 
+G.R_li = L.f_i*G.D_RE; %Laufbahnradius des Innenrings (Dowson S.48) %der im Lagerquerschnitt gemessene Laufbahnrillenradius des Außenrings
+G.R_la = L.f_a*G.D_RE; %Laufbahnradius des Außenrings
+
 b_li = (L.D-L.d-2*L.D_RE-L.c)/2; %Breite von Innenring (bezogen auf Durchmesser)
 b_la = b_li; 
 G.b_li = b_li*(1+L.ring.alpha_t*G.Delta_Ti_norm);
@@ -48,8 +51,6 @@ G.P_d = G.D-G.d-2*G.D_RE-G.b_li-G.b_la;
 G.B_i = L.B_i*(1+L.ring.alpha_t*G.Delta_Ti_norm);
 G.B_a = L.B_a*(1+L.ring.alpha_t*G.Delta_Ta_norm);
 
-G.R_li = L.f_i*G.D_RE; %Laufbahnradius des Innenrings (Dowson S.48) %der im Lagerquerschnitt gemessene Laufbahnrillenradius des Außenrings
-G.R_la = L.f_a*G.D_RE; %Laufbahnradius des Außenrings
 
 G.b = L.f_i + L.f_a - 1; % Harris B
 if ~isempty(options.alpha)
@@ -87,6 +88,10 @@ G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius                 %
 
 G.R_ri = G.R_bi + G.R_RE*(1-cos(G.alpha)); % Race radius at point of contact
 G.R_ra = G.R_ba - G.R_RE*(1-cos(G.alpha));
+G.d_bLi = L.d_bLi * (1 + L.ring.alpha_t * G.Delta_Ti_norm);
+G.d_bRi = L.d_bRi * (1 + L.ring.alpha_t * G.Delta_Ti_norm);
+G.D_bLa = L.D_bLa * (1 + L.ring.alpha_t * G.Delta_Ta_norm);
+G.D_bRa = L.D_bRa * (1 + L.ring.alpha_t * G.Delta_Ta_norm);
 
 G.sigmarho_i = (2./G.R_RE + 1./G.R_ri - 1./G.R_li);
 G.sigmarho_a = (2./G.R_RE - 1./G.R_ra - 1./G.R_la);                         % abh. von alpha
diff --git a/@BearImp/setBearing.m b/@BearImp/setBearing.m
index 4a28309b1e1624201b73c1887aa45405116ff29b..9374607efad88e293029b4e82d0bab19c27f4cf6 100644
--- a/@BearImp/setBearing.m
+++ b/@BearImp/setBearing.m
@@ -47,6 +47,13 @@ function setBearing(obj,name)
 
     L.d_m = (L.d + L.D)/2;
 
+    switch L.Type
+        case 'ACBB'
+            if L.c > 1 % L.c doubles as input of free contact angle for ACBB. A value less than 1 is interpreted as a clearance
+                L.c = 2*L.D_RE*(L.f_i+L.f_a-1)*(1-cosd(L.c));
+            end
+    end
+
      if ~isfield(L,'RE_order')
         L.RE_order_num = ones(1,L.Z);
         L.allBallsSameMaterial = true; % wenn die Variable RE_order nicht erzeugt wird, ist nur ein RE Material vorhanden
diff --git a/BearImpOptions.m b/BearImpOptions.m
index 23104846f1c7ee10fa6831679aad26bedc57008d..cb28646f544cb9b46707e5d1a9461353ac7df190 100644
--- a/BearImpOptions.m
+++ b/BearImpOptions.m
@@ -35,8 +35,8 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
             possibleOption.H.thermalCorrection = {'on','off'};
             possibleOption.H.starvationCorrection = {'on','off'};
             possibleOption.H.roughnessCorrection = {'on','off'};
-            possibleOption.C.unloadedRE  = {           'neglect','stateOfTheArt',                'Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
-            possibleOption.C.outsideArea = {'k-factor','neglect','stateOfTheArt','Schneider_k_c','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
+            possibleOption.C.unloadedRE  = {           'neglect','stateOfTheArt',                'Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D','Puchtler2025'};
+            possibleOption.C.outsideArea = {'k-factor','neglect','stateOfTheArt','Schneider_k_c','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D','Puchtler2025'};
             possibleOption.C.k_vh_factor = {'on','off'};
             possibleOption.C.pressureDistribution = {'on','off'};
             possibleOption.C.rhoRatio = {'Hartung','Dowson','Bradbury'};
diff --git a/InputData.xlsx b/InputData.xlsx
index 7491cb4460984b4e0e7b0e3db642a0e459ee3229..08e6ac3ea019368fc5993b662dab839618ea8974 100644
Binary files a/InputData.xlsx and b/InputData.xlsx differ