diff --git a/@BearImp/BearImp.m b/@BearImp/BearImp.m
index 7fdfc4a8a6f33cf59a0f614a9048664912246382..50490b3fa2f4c4f61f3d441f4a7571077ca29d33 100644
--- a/@BearImp/BearImp.m
+++ b/@BearImp/BearImp.m
@@ -103,8 +103,8 @@ classdef BearImp < handle & matlab.mixin.Copyable
     end
 
     methods (Access = public, Static)
-        C_out = calcCap_puchtler2025    (s,R_WK,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha,h_0,a,b      ,options)
-        C_out = calcCap_semianalytisch3D(s,R_WK,R_R,   B_R   ,B,R_L,epsilon_r,      h_0,a,b,nt,np,options)
+        C_out = calcCap_puchtler2025    (s,R_WK,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha,h_0,a,b      ,method)
+        C_out = calcCap_semianalytisch3D(s,R_WK,R_R,   B_R   ,B,R_L,epsilon_r,      h_0,a,b,nt,np,method)
     end
     
     methods
diff --git a/@BearImp/calcCap.m b/@BearImp/calcCap.m
index f89d241e336e0c57ab865ca8a7511e077a15c13a..06caa7b01c7b50899fb0af54194d89688d13f457 100644
--- a/@BearImp/calcCap.m
+++ b/@BearImp/calcCap.m
@@ -298,8 +298,8 @@ function runPuchtler2025(indices)
     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),'CmethodDeformedArea',C.method.deformedArea,'r2smallCriterion',C.method.r2smallCriterion);
-    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),'CmethodDeformedArea',C.method.deformedArea,'r2smallCriterion',C.method.r2smallCriterion);
+    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),C.method);
+    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),C.method);
 
     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));
diff --git a/@BearImp/calcCap_puchtler2025.m b/@BearImp/calcCap_puchtler2025.m
index 511b75812d9daaf1284d37090e101b146f58efce..33c2328419354f155d0e2a03a645dadfea160c3f 100644
--- a/@BearImp/calcCap_puchtler2025.m
+++ b/@BearImp/calcCap_puchtler2025.m
@@ -1,4 +1,4 @@
-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)
+function C_out = calcCap_puchtler2025(s,R_RE,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha,h_0,a,b,method)
     
     arguments
         s
@@ -13,8 +13,7 @@ function C_out = calcCap_puchtler2025(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'
+        method = struct('deformedArea','neglect','outlet','oil','r2smallCriterion','r<1.00125*R_RE')
     end
 
     Delta_x =              -(R_R - R_RE - s)*sin(alpha);
@@ -22,16 +21,26 @@ function C_out = calcCap_puchtler2025(s,R_RE,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,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_groove = @(phi,theta) r_easy(phi,theta,R_RE,R_R,R_L,Delta_x,Delta_z,h_0,method);
     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)));
+    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)));
+    if method.outlet == 'air'
+        Theta_0lIn  = @(phi) max(atan((R_RE*sin(alpha))./(-Delta_z*cos(phi) + sign(R_L)*sqrt(R_L^2 -Delta_z^2*sin(phi).^2))) - max(a*sqrt(max(1/R_RE^2-phi.^2/b^2,0)),phi.*(phi>0)),Theta_1l(phi));
+        Theta_0rIn  = @(phi) min(atan((R_RE*sin(alpha))./(-Delta_z*cos(phi) + sign(R_L)*sqrt(R_L^2 -Delta_z^2*sin(phi).^2))) + max(a*sqrt(max(1/R_RE^2-phi.^2/b^2,0)),phi.*(phi>0)),Theta_1r(phi));
+        Theta_0lOut = @(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_0rOut = @(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));
+        phi_0 = a/R_RE/sqrt(1+a^2/b^2);
+        phi_1l = fzero(@(phi) Theta_0l(phi) - phi - Theta_1l(phi),pi/4);
+        phi_1r = fzero(@(phi) Theta_0r(phi) + phi - Theta_1r(phi),pi/4);
+        % fprintf('phi_0 = %d | phi_1l = %d | phi_1r = %d \n',phi_0,phi_1l,phi_1r)
+    end
+    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;
@@ -42,16 +51,28 @@ function C_out = calcCap_puchtler2025(s,R_RE,R_R,R_ZL,R_ZR,B,R_L,epsilon_r,alpha
     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);
+    
+    switch method.outlet
+        case 'oil'
+            absTol = 25e-4; % 25e-4 *2 corresponds to 0.1 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);
+        case 'air'
+            absTol = 50e-4; % 50e-4 corresponds to 0.1 pF with a permittivity of 2.2
+            fun_grooveOut = @(phi,theta) R_RE^2.*1./(h_0./epsilon_r + (r_groove(phi,theta)-R_RE-h_0)) .* cos(theta);
+            C_grooveRIn  = 8.854187e-12 * epsilon_r * integral2(fun_groove   ,-phi_1,phi_1r,Theta_0rIn,Theta_1r   ,'AbsTol',absTol); 
+            C_grooveLIn  = 8.854187e-12 * epsilon_r * integral2(fun_groove   ,-phi_1,phi_1l,Theta_1l  ,Theta_0lIn ,'AbsTol',absTol);
+            C_grooveROut = 8.854187e-12             * integral2(fun_grooveOut, phi_0,phi_1,Theta_0rOut,Theta_0rIn ,'AbsTol',absTol); 
+            C_grooveLOut = 8.854187e-12             * integral2(fun_grooveOut, phi_0,phi_1,Theta_0lIn ,Theta_0lOut,'AbsTol',absTol);
+            C_grooveR = C_grooveRIn + C_grooveROut;
+            C_grooveL = C_grooveLIn + C_grooveLOut;
+    end
     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)
+function r = r_easy(phi,theta,R_RE,R_R,R_L,Delta_x,Delta_z,h_0,method)
     assert(all(size(phi)==size(theta)))
 
     r = nan(size(phi));
@@ -65,13 +86,13 @@ function r = r_easy(phi,theta,R_RE,R_R,R_L,Delta_x,Delta_z,h_0,options)
         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
+    switch method.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
+    switch method.deformedArea
         case 'neglect'
             r(r2small) = inf;
         case 'filmThickness'
diff --git a/BearImpOptions.m b/BearImpOptions.m
index 57bb705d9718b53405ded9d3b8c7625ac7665774..0931166f19bbd7b1141d1ea8cb7454a88104d862 100644
--- a/BearImpOptions.m
+++ b/BearImpOptions.m
@@ -38,6 +38,7 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
             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.deformedArea= {'neglect','filmThickness'};
+            possibleOption.C.outlet      = {'oil','air'};
             possibleOption.C.r2smallCriterion = {'r<1.00125*R_RE','r<h_0+R_RE'};
             possibleOption.C.k_vh_factor = {'on','off'};
             possibleOption.C.pressureDistribution = {'on','off'};
@@ -57,6 +58,7 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
             defaultOption.C.unloadedRE  = 'semianalytisch3D';
             defaultOption.C.outsideArea = 'semianalytisch3D';
             defaultOption.C.deformedArea= 'neglect';
+            defaultOption.C.outlet = 'oil';
             defaultOption.C.r2smallCriterion = 'r<1.00125*R_RE';
             defaultOption.C.k_vh_factor = 'off';
             defaultOption.C.pressureDistribution = 'on';
@@ -138,6 +140,7 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
             if obj.C.roughnessCorrection  ~= 'off'   , str = [str 'r']; end
             if obj.C.deformedArea == 'filmThickness' , str = [str 'h']; end
             if obj.C.r2smallCriterion == 'r<h_0+R_RE', str = [str 'a']; end
+            if obj.C.outlet == 'air', str = [str 'o']; end
             str = [str lower(obj.C.rhoRatio.char(1))];
         end
         function charOut = char(obj,n)