diff --git a/@BearImp/BearImp.m b/@BearImp/BearImp.m
index 9e51a353cb7eb2da5ce48e480e807908bd0a232c..50490b3fa2f4c4f61f3d441f4a7571077ca29d33 100644
--- a/@BearImp/BearImp.m
+++ b/@BearImp/BearImp.m
@@ -13,8 +13,8 @@ classdef BearImp < handle & matlab.mixin.Copyable
         resolutionPerRotation (1,1) double {mustBeNonnegative, mustBeInteger}  % Anzahl an Zwischenschritten des Winkels psi innerhalb einer 360°-Drehung
         L (1,1) struct        % Lagerparameter (gegeben)
         S (1,1) struct        % Schmierstoffparameter (gegeben)
-        method  (1,1) struct  % struct der Rechenmethoden mit möglichen Feldern T,R,G,B,H,Z
         AddOn (1,1) struct    % Installierte AddOns
+        method  (1,1) BearImpOptions  % BearImpOptions-Objekt mit gewählten Optionen der Rechenmethoden
         
         
     % Structs mit Berechnungsparamtern
@@ -28,12 +28,13 @@ classdef BearImp < handle & matlab.mixin.Copyable
     end
     
     
-    properties (Access = private)
+    properties (Access = protected)
     % Eingangsparameter Speicher
         privateInputParameters = struct('F_r',nan,'F_a',nan,'omega',nan,'T_Oil',nan,'psi',nan,'psi_calc',nan,'psi_all',nan,'ind_psi_all',nan,'f',nan,'L',struct,'S',struct,'resolutionPerRotation',nan)
-        privateMethod  = struct
         privateResults = struct('T',struct,'R',struct,'G',struct,'B',struct,'H',struct,'C',struct,'Z',struct)
         privateAddOn = struct('Parallel_Computing_Toolbox', 'auto', 'Optimization_Toolbox', 'auto')
+        privateMethod = BearImpOptions
+        optionChangedListenerHandle
     end
     
     properties (Access = public)
@@ -42,6 +43,7 @@ classdef BearImp < handle & matlab.mixin.Copyable
         materials
         inputDataTableFile = 'InputData.xlsx';
         PlotLabelTableFile = 'PlotLabel.csv';
+        displayWarnings = true;
     end
     
     properties (Dependent, SetAccess = protected, GetAccess = public)
@@ -67,20 +69,21 @@ classdef BearImp < handle & matlab.mixin.Copyable
             if nargin > 0
                 switch setup
                     case 'default'
-                        obj.setBearing('6205-C3');
-                        obj.setLube('FVA Referenz-Öl III');
+                        obj.setBearing('DGBB_6205_C3');
+                        obj.setLube('FVA Referenz-Öl III barus');
                         obj.F_r = 1000;
                         obj.F_a = 0;
                         obj.omega = 3000/60*2*pi;
-                        obj.T_Oil = 50;
-                        obj.resolutionPerRotation = 360;
+                        obj.T_Oil = 60;
+                        obj.resolutionPerRotation = 90;
                 end
             end
-            
+            obj.method = BearImpOptions;
+            obj.optionChangedListenerHandle = addlistener(obj.method,'optionChanged',@(src,event) obj.methodChangedFcn(src,event));
         end
     end
     
-    methods (Access = public)               % Später protected?  #######################
+    methods (Access = public)
     % Berechnungsmethoden in einzelnen Dateien
         calcLub  (obj)
         calcClear(obj)
@@ -89,8 +92,6 @@ classdef BearImp < handle & matlab.mixin.Copyable
         varargout = calcFilm(obj, options)
         calcCap  (obj)
         calcImp  (obj)
-    end
-    methods (Access = public)
         setBearing(obj,name)
         setFitting(obj,shaft,housing)
         setLube   (obj,name)
@@ -100,6 +101,11 @@ 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,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
     % set-Methoden
@@ -148,8 +154,6 @@ classdef BearImp < handle & matlab.mixin.Copyable
             obj.privateInputParameters.resolutionPerRotation = nan;
         end
         function set.ind_psi_all(obj,ind_psi_all)
-%             obj.up2date.R = false;
-%             obj.up2date.B = false;
             obj.privateInputParameters.ind_psi_all = ind_psi_all;
         end
         function set.AddOn(obj,AddOn)
@@ -200,25 +204,15 @@ classdef BearImp < handle & matlab.mixin.Copyable
             obj.up2date.C = false;
             obj.privateInputParameters.S = S;
         end
-        function set.method(obj,method)
-            possibleMethods.check(method,true);
-            somethingChanged = false;
-            for ii = {'T','R','G','B','H','C','Z'}
-                if isfield(obj.method,ii) && isfield (method,ii)
-                    if ~strcmp(obj.method.(ii{1}), method.(ii{1}))
-                        obj.up2date.(ii{1}) = false;
-                        somethingChanged = true;
-                    end
-                elseif isfield(obj.method,ii{1}) || isfield (method,ii{1})
-                    obj.up2date.(ii{1}) = false;
-                    somethingChanged = true;
+        function set.method(obj,newMethod)
+            oldMethod = obj.privateMethod.copy;
+            obj.privateMethod = newMethod;
+            [isEqual,methodsChanged] = eq(oldMethod,newMethod);
+            if ~isEqual
+                for mc = methodsChanged
+                    obj.up2date.(mc{:}) = false;
                 end
             end
-            if somethingChanged
-                obj.privateMethod = method;
-            else
-                warning('Keine Methode geändert')
-            end
         end
         
         function set.up2date(obj,up2date)
@@ -346,6 +340,9 @@ classdef BearImp < handle & matlab.mixin.Copyable
             end
             
    
+        end
+        function method = get.method(obj)
+            method = obj.privateMethod;
         end
         function L = get.L(obj)
             L = obj.privateInputParameters.L;
@@ -353,9 +350,6 @@ classdef BearImp < handle & matlab.mixin.Copyable
         function S = get.S(obj)
             S = obj.privateInputParameters.S;
         end
-        function method = get.method(obj)
-            method = obj.privateMethod;
-        end
         function T = get.T(obj)
             T = obj.privateResults.T;
         end
@@ -400,6 +394,12 @@ classdef BearImp < handle & matlab.mixin.Copyable
     end
     
     methods (Access = protected)
+        function methodChangedFcn(obj,src,~)
+            for mc = src.optionsLastChanged
+                obj.up2date.(mc{:}) = false;
+            end
+            src.optionsLastChanged = {};
+        end
         function executeAllObj(obj,functionHandle)
         % Führt eine Funktion für alle Einträge von obj aus
         % Funktioniert aktuell nur für Funktionen ohne Übergabeparameter
@@ -407,5 +407,36 @@ classdef BearImp < handle & matlab.mixin.Copyable
                 feval(functionHandle,obj(ii))
             end
         end
+        function cp = copyElement(obj)
+            cp = copyElement@matlab.mixin.Copyable(obj);
+            cp.method = obj.method.copy; % creat a deep copy of the method
+            cp.optionChangedListenerHandle = addlistener(cp.method,'optionChanged',@(src,event) cp.methodChangedFcn(src,event));
+        end
+        function s = saveobj(obj)
+            obj.optionChangedListenerHandle = [];
+            s = obj;
+        end
+    end
+    methods (Static, Hidden)
+        function obj = loadobj(s)
+            s.optionChangedListenerHandle = addlistener(s.method,'optionChanged',@(src,event) s.methodChangedFcn(src,event));
+            s.method.addSetMethods
+            obj = s;
+        end
+    end
+    methods (Static)
+        function obj = empty(n,m)
+            % Creates a nxm array of empty BearImp objects
+            arguments
+                n {mustBeInteger, mustBeScalarOrEmpty}
+                m {mustBeScalarOrEmpty} = n
+            end
+            obj(n,m) = BearImp;
+            for ii = 1:n
+                for jj = 1:m
+                    obj(ii,jj) = BearImp;
+                end
+            end
+        end
     end
 end
\ No newline at end of file
diff --git a/@BearImp/calcCap.m b/@BearImp/calcCap.m
index d391d9afcd16330bdad8f6057326d9a9b71996ac..06caa7b01c7b50899fb0af54194d89688d13f457 100644
--- a/@BearImp/calcCap.m
+++ b/@BearImp/calcCap.m
@@ -15,9 +15,9 @@ assert(obj.up2date.T,'Schmierstoffparameter noch nicht berechnet')
 assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht berechnet')
 assert(obj.up2date.B,'Belastungsverteilung noch nicht berechnet')
 assert(obj.up2date.H,'Schmierfilmdicke noch nicht berechnet')
-L = obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; H = obj.H; AddOn = obj.AddOn; psi = obj.psi;
+L = obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; H = obj.H; AddOn = obj.AddOn; psi = obj.psi; T_Oil = obj.T_Oil;
 C = obj.C;
-method = possibleMethods.addDefault(obj.method).C;
+C.method = obj.method.C;
 if ~isfield(C,'k_C')
     C.k_C = 1;
 end
@@ -62,22 +62,53 @@ for posBall_conductive=1:L.numberOfConductiveBalls
     tempR_li = G.R_li(ballMaterialInd(posBall_conductive));
     tempR_la = G.R_la(ballMaterialInd(posBall_conductive));
     tempR_RE = G.R_RE(ballMaterialInd(posBall_conductive));
-    if any(strcmp(method.outsideArea,{'Schneider_k_c','Schneider_k_vh'}))
-        tempG = H.G (:,:,posBall_conductive);
-        tempW = H.W (:,:,posBall_conductive);
+    if any([C.method.outsideArea == {'Schneider_k_c','Schneider_k_vh'} C.method.pressureDistribution == 'on'])
+        tempG = H.G(:,:,posBall_conductive);
+        tempW = H.W(:,:,posBall_conductive);
     end
         
     contactInd=B.contactInd(posBall_conductive,:);
     
-    C.rhoRatio = 1  +  (9.73e-3 .* (temp_p.*0.101972e-6).^0.75)  ./  (T.nu_38.*1e6).^0.0385;
-    C.epsilon_primeOil = (S.epsilon_Oel + 2 + 2*(S.epsilon_Oel-1).*C.rhoRatio)  ./  (S.epsilon_Oel + 2 - (S.epsilon_Oel-1).*C.rhoRatio);
-    C.C_Hertz(:,contactInd,posBall_conductive) = obj.epsilon_0 .* C.epsilon_primeOil(:,contactInd) .* temp_A(:,contactInd) ./ temp_h_0(:,contactInd);
+    switch C.method.rhoRatio
+        case 'Hartung' % Brüser 1972 (4.33, p. 43) (default)
+            C.rhoRatio = @(p) 1  +  (9.73e-3 .* (p.*0.101972e-6).^0.75)  ./  (T.nu_38.*1e6).^0.0385;
+        case 'Dowson'  % Brüser 1972 (4.27, p. 41)
+            C.rhoRatio = @(p) 1 + 0.0057.*(p.*0.101972e-6)./(1+0.0165.*(p.*0.101972e-6)); % assuming [p] = kP/mm²
+        case 'Bradbury'% Brüser 1972 (4.30, p. 43)
+            C.rhoRatio = @(p) 1 ./ (1 - 0.087.*log(0.0684.*(p.*0.101972e-6)+1));          % assuming [p] = kP/mm²
+        case 'Tait'    % Bair 2007 (Table 6.4, p. 123)
+            K_0prime = 11; K_00 = 9e9; beta_K = 6.5e-3; a_V = 8e-4;
+            K_0 = K_00 * exp(-beta_K*(T_Oil+273.15));    % by Fakhreddine and Zoller in Bair 2007 (4.12, p. 69)
+            VbyV_0 = @(p) 1 - 1/(1+K_0prime)*log(1+p./K_0*(1+K_0prime)); % Bair 2007 (4.6 , p. 68)
+            V_0byV_R = 1./(1 - S.alpha_rho.*(T_Oil-15));                 % Bair 2007 (4.13, p. 70)
+            %V_0byV_R = 1 + a_V*(T_Oil-15);                              % Bair 2007 (4.14, p.70) depreciated
+            C.rhoRatio = @(p) 1./(VbyV_0(p) .* V_0byV_R);
+    end
+    C.epsilon_primeOil = @(p) (S.epsilon_Oel + 2 + 2*(S.epsilon_Oel-1).*C.rhoRatio(p))  ./  (S.epsilon_Oel + 2 - (S.epsilon_Oel-1).*C.rhoRatio(p));
+    if C.method.pressureDistribution == 'on'
+        for io = 1:size(temp_a,1)
+            for singleContInd = find(contactInd)
+                p_max = 3/2*temp_p(io,singleContInd); % Hamrock/Dowson (3.6), S. 71
+                p_fcn = @(x,y) p_max .* sqrt(1-(x./temp_a(io,singleContInd)).^2-(y./temp_b(io,singleContInd)).^2);
+                epsByH_fcn = @(x,y) 4 * obj.epsilon_0 .* C.epsilon_primeOil(p_fcn(x,y)) ./ temp_h_0(io,singleContInd);
+                C.C_Hertz(io,singleContInd,posBall_conductive) = integral2(epsByH_fcn,0,temp_a(io,singleContInd),0,@(x) temp_b(io,singleContInd).*sqrt(1-(x./temp_a(io,singleContInd)).^2));
+            end
+        end
+    else
+        C.C_Hertz(:,contactInd,posBall_conductive) = obj.epsilon_0 .* C.epsilon_primeOil(temp_p(:,contactInd)) .* temp_A(:,contactInd) ./ temp_h_0(:,contactInd);
+    end
+
+    if C.method.roughnessCorrection == 'Napel'
+        A = (L.Rz_RE + [L.Rz_IR;L.Rz_OR])/4;
+        C.napelFactor = 1./sqrt(1 - (A./temp_h_0(:,contactInd)).^2);
+        C.C_Hertz(:,contactInd,posBall_conductive) = C.C_Hertz(:,contactInd,posBall_conductive) .* C.napelFactor;
+    end
  
     R_el_single(:,contactInd,posBall_conductive) = 1./C.k_R .* S.rho_el .* temp_h_0(:,contactInd) ./ temp_A(:,contactInd);
     C.DeltaR_i(contactInd)    = tempR_li - tempR_RE - temp_h_0(1,contactInd);
     C.DeltaR_a(contactInd)    = tempR_la - tempR_RE - temp_h_0(2,contactInd);
 
-    if size(B.noContactInd,2) % Wenn Kraftfreie Wälzkörper vorhanden
+    if size(B.noContactInd,2) % Wenn kraftfreie Wälzkörper vorhanden
         C.DeltaR_i(B.noContactInd(posBall_conductive,:)') = tempR_li - tempR_RE - temp_s(1,B.noContactInd(posBall_conductive,:)');
         C.DeltaR_a(B.noContactInd(posBall_conductive,:)') = tempR_la - tempR_RE - temp_s(2,B.noContactInd(posBall_conductive,:)');
     end
@@ -87,12 +118,12 @@ for posBall_conductive=1:L.numberOfConductiveBalls
     C.phi_2i = asin(L.B_i/G.D_RE(ballMaterialInd(posBall_conductive)));
     C.phi_2a = asin(L.B_a/G.D_RE(ballMaterialInd(posBall_conductive)));
     
-    if any([strcmp(method.unloadedRE,{'Leander_Radial','LeanderSteffen'}) strcmp(method.outsideArea,{'Leander_Radial','LeanderSteffen'})])
+    if any([C.method.unloadedRE=={'Leander_Radial','LeanderSteffen'} C.method.outsideArea=={'Leander_Radial','LeanderSteffen'}])
         calcZg
     end
     
     if any(B.noContactInd,'all') % Wenn Kraftfreie Wälzkörper vorhanden
-        switch method.unloadedRE
+        switch C.method.unloadedRE
             case 'neglect'
             case 'stateOfTheArt'
                 stateOfTheArt(find(contactInd==0))
@@ -109,10 +140,12 @@ for posBall_conductive=1:L.numberOfConductiveBalls
                 for recentVec=length(vecStruct)
                     runSemianalytisch3D(vecStruct{recentVec})
                 end
+            case 'Puchtler2025'
+                runPuchtler2025(find(~contactInd))
         end
     end
 
-    switch method.outsideArea
+    switch C.method.outsideArea
         case 'neglect'
         case 'k-factor'
             C.C_out(:,contactInd,posBall_conductive) = (C.k_C-1) .* C.C_Hertz(:,contactInd,posBall_conductive);
@@ -120,11 +153,7 @@ for posBall_conductive=1:L.numberOfConductiveBalls
             stateOfTheArt(find(contactInd==1))
         case 'Schneider_k_c'
             C.k_c  = 9.5923 * H.U.^0.3305 .* tempG.^0.3413 .* tempW.^(-0.3342) .* G.k.^0.1391;
-            %C.C_Hertz(:,contactInd,posBall_conductive) = k_c(:,contactInd,posBall_conductive) .* C.C_Hertz(:,contactInd,posBall_conductive);
             C.C_out(:,contactInd,posBall_conductive) = (C.k_c (:,contactInd,posBall_conductive)-1) .* C.C_Hertz(:,contactInd,posBall_conductive);
-        case 'Schneider_k_vh'
-            C.k_vh = 0.7772 * H.U.^0.0451 .* tempG.^0.1522 .* tempW.^(-0.0245) .* G.k.^(-0.0921);
-            C.C_out(:,contactInd,posBall_conductive) = (C.k_vh(:,contactInd,posBall_conductive)-1) .* C.C_Hertz(:,contactInd,posBall_conductive);
         case 'Leander_Parallel'
             leander_parallel(B.contactInd)
         case 'Leander_Radial'
@@ -138,9 +167,16 @@ 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'
+        C.k_vh = 0.7772 * H.U.^0.0451 .* tempG.^0.1522 .* tempW.^-0.0245 .* G.k.^-0.0921;
+        C.C_Hertz(:,:,posBall_conductive) = C.C_Hertz(:,:,posBall_conductive) .* C.k_vh(:,:,posBall_conductive);
     end
 
-    C_el_single(:,:,posBall_conductive) = C.C_Hertz(:,:,posBall_conductive) + C.C_out(:,:,posBall_conductive);
+    C_el_single(:,:,posBall_conductive) = C.C_Hertz(:,:,posBall_conductive) + C.C_out(:,:,posBall_conductive); %#ok<AGROW>
     
 end
 
@@ -195,36 +231,36 @@ end
 
 
 function leander_parallel(indices)
-    C.R_hi = @(v) G.R_bi.*sqrt(1-(tempR_RE/G.R_bi.*cos(v)).^2);
-    C.R_ha = @(v) G.R_ba.*sqrt(1-(tempR_RE/G.R_ba.*cos(v)).^2);
+    C.R_hi = @(v) G.R_ri.*sqrt(1-(tempR_RE/G.R_ri.*cos(v)).^2);
+    C.R_ha = @(v) G.R_ra.*sqrt(1-(tempR_RE/G.R_ra.*cos(v)).^2);
     
     for ii = indices
-        C.I_i = @(v,p) abs(tempR_RE^2.*sin(v))./(-(C.R_hi(v)-tempR_li).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_hi(v)+tempR_li)).^2)+tempR_RE.*(1+sin(v).*cos(p))-tempR_li+G.R_bi+B.s(ii));
-        C.I_a = @(v,p) abs(tempR_RE^2.*sin(v))./( (C.R_ha(v)+tempR_la).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_ha(v)+tempR_la)).^2)+tempR_RE.*(1-sin(v).*cos(p))-tempR_la-G.R_ba+B.s(ii));
+        C.I_i = @(v,p) abs(tempR_RE^2.*sin(v))./(-(C.R_hi(v)-tempR_li).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_hi(v)+tempR_li)).^2)+tempR_RE.*(1+sin(v).*cos(p))-tempR_li+G.R_ri+B.s(ii));
+        C.I_a = @(v,p) abs(tempR_RE^2.*sin(v))./( (C.R_ha(v)+tempR_la).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_ha(v)+tempR_la)).^2)+tempR_RE.*(1-sin(v).*cos(p))-tempR_la-G.R_ra+B.s(ii));
         C.C_out(1,ii) = 4 * obj.epsilon_0 * S.epsilon_Oel * integral2(C.I_i,C.phi_1i(ii)+pi/2,C.phi_2i+pi/2, pi,1.5*pi);
         C.C_out(2,ii) = 4 * obj.epsilon_0 * S.epsilon_Oel * integral2(C.I_a,C.phi_1a(ii)+pi/2,C.phi_2a+pi/2,  0,0.5*pi);
     end
 end
 
 function leander_radial(indices)
-    C.c_i(1,indices) =  tempR_li - G.R_bi - tempR_RE - temp_s(1,indices);          % TODO: B.s an der Stelle sinnvoll? ###
-    C.c_a(1,indices) = -tempR_la - G.R_ba + tempR_RE + temp_s(2,indices);          % TODO: B.s an der Stelle sinnvoll? ###
+    C.c_i(1,indices) =  tempR_li - G.R_ri - tempR_RE - temp_s(1,indices);          % TODO: B.s an der Stelle sinnvoll? ###
+    C.c_a(1,indices) = -tempR_la - G.R_ra + tempR_RE + temp_s(2,indices);          % TODO: B.s an der Stelle sinnvoll? ###
     
     for ii = indices
-        C.I_i = @(g,t) abs(tempR_li*(tempR_li.*sin(g)+G.R_bi))./(sqrt(tempR_li.^2+G.R_bi.^2+C.c_i(ii').^2+2*G.R_bi*tempR_li.*sin(g)+2*C.c_i(ii')*(tempR_li.*sin(g)+G.R_bi).*cos(t))-tempR_RE);
-        C.I_a = @(g,t) abs(tempR_la*(tempR_la.*sin(g)+G.R_ba))./(sqrt(tempR_la.^2+G.R_ba.^2+C.c_a(ii').^2+2*G.R_ba*tempR_la.*sin(g)+2*C.c_a(ii')*(tempR_la.*sin(g)+G.R_ba).*cos(t))-tempR_RE);
+        C.I_i = @(g,t) abs(tempR_li*(tempR_li.*sin(g)+G.R_ri))./(sqrt(tempR_li.^2+G.R_ri.^2+C.c_i(ii').^2+2*G.R_ri*tempR_li.*sin(g)+2*C.c_i(ii')*(tempR_li.*sin(g)+G.R_ri).*cos(t))-tempR_RE);
+        C.I_a = @(g,t) abs(tempR_la*(tempR_la.*sin(g)+G.R_ra))./(sqrt(tempR_la.^2+G.R_ra.^2+C.c_a(ii').^2+2*G.R_ra*tempR_la.*sin(g)+2*C.c_a(ii')*(tempR_la.*sin(g)+G.R_ra).*cos(t))-tempR_RE);
         C.C_out(1,ii,posBall_conductive) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(C.I_i,-pi/2 + C.g_1i(ii,posBall_conductive),   -C.g_2i, 0, pi/9);
         C.C_out(2,ii,posBall_conductive) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(C.I_a, pi/2 + C.g_1a(ii,posBall_conductive), pi-C.g_2a, 0, pi/9);
     end 
 end
 
 function leander_steffen(indices)
-    C.dA_i = @(p) tempR_li * (G.R_bi + tempR_li*cos(p));
-    C.dA_a = @(p) tempR_la * (G.R_ba + tempR_la*cos(p));
-    C.x_kmi = G.R_bi + C.DeltaR_i;
-    C.x_kma = G.R_ba - C.DeltaR_a;
-    C.h_i = @(p,t,ii) sqrt(G.R_bi^2 + tempR_li^2 + C.x_kmi(ii)^2 + 2*G.R_bi*tempR_li*cos(p) - 2*C.x_kmi(ii)*cos(t).*(G.R_bi+tempR_li*cos(p)))  -  tempR_RE;
-    C.h_a = @(p,t,ii) sqrt(G.R_ba^2 + tempR_la^2 + C.x_kma(ii)^2 + 2*G.R_ba*tempR_la*cos(p) - 2*C.x_kma(ii)*cos(t).*(G.R_ba+tempR_la*cos(p)))  -  tempR_RE;
+    C.dA_i = @(p) tempR_li * (G.R_ri + tempR_li*cos(p));
+    C.dA_a = @(p) tempR_la * (G.R_ra + tempR_la*cos(p));
+    C.x_kmi = G.R_ri + C.DeltaR_i;
+    C.x_kma = G.R_ra - C.DeltaR_a;
+    C.h_i = @(p,t,ii) sqrt(G.R_ri^2 + tempR_li^2 + C.x_kmi(ii)^2 + 2*G.R_ri*tempR_li*cos(p) - 2*C.x_kmi(ii)*cos(t).*(G.R_ri+tempR_li*cos(p)))  -  tempR_RE;
+    C.h_a = @(p,t,ii) sqrt(G.R_ra^2 + tempR_la^2 + C.x_kma(ii)^2 + 2*G.R_ra*tempR_la*cos(p) - 2*C.x_kma(ii)*cos(t).*(G.R_ra+tempR_la*cos(p)))  -  tempR_RE;
     for ii = indices
         x1=double(C.g_1i(ii));
         x2=double(C.g_2i);
@@ -239,7 +275,7 @@ end
 function tobias_steffen(indices)
     C.h_i = @(alpha,beta,ii) (sqrt(tempR_li^2 - C.DeltaR_i(ii)^2*sin(alpha).^2) - C.DeltaR_i(ii)*cos(alpha))  ./  cos(beta) - tempR_RE;
     C.h_a = @(alpha,beta,ii) (sqrt(tempR_la^2 - C.DeltaR_a(ii)^2*sin(alpha).^2) - C.DeltaR_a(ii)*cos(alpha))  ./  cos(beta) - tempR_RE;
-    if strcmp(method.unloadedRE,'TobiasSteffen_Kugelfläche')
+    if C.method.unloadedRE == 'TobiasSteffen_Kugelfläche'
         C.dA_i = @(~,beta,~) tempR_RE^2 * cos(beta);
         C.dA_a = C.dA_i;
     else
@@ -252,6 +288,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),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));
+        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')
@@ -264,8 +333,8 @@ function runSemianalytisch3D(indices)
     
     if AddOn.Parallel_Computing_Toolbox
 
-        tmp_semi_i = @(runSemi) semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_bi, S.epsilon_Oel,temp_a(1,runSemi),temp_b(1,runSemi),110,110);
-        tmp_semi_a = @(runSemi) semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ba, S.epsilon_Oel,temp_a(2,runSemi),temp_b(2,runSemi),110,110);
+        tmp_semi_i = @(runSemi) BearImp.calcCap_semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_ri, S.epsilon_Oel,temp_h_0(1,runSemi),temp_a(1,runSemi),temp_b(1,runSemi),110,110,'CmethodDeformedArea',C.method.deformedArea,'r2smallCriterion',C.method.r2smallCriterion);
+        tmp_semi_a = @(runSemi) BearImp.calcCap_semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ra, S.epsilon_Oel,temp_h_0(2,runSemi),temp_a(2,runSemi),temp_b(2,runSemi),110,110,'CmethodDeformedArea',C.method.deformedArea,'r2smallCriterion',C.method.r2smallCriterion);
 
         parfor runSemi = indices
             C_semi_i(runSemi) = feval(tmp_semi_i, runSemi); %da parfor loop nicht auf nested functions (semianalytisch3D) direkt zugreifen kann
@@ -275,8 +344,9 @@ function runSemianalytisch3D(indices)
     else
 
         for runSemi = 1 : numel(temp_s)/2
-           C_semi_i(runSemi) = semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_bi, S.epsilon_Oel,temp_a(1,runSemi),temp_b(1,runSemi),110,110);
-           C_semi_a(runSemi) = semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ba, S.epsilon_Oel,temp_a(2,runSemi),temp_b(2,runSemi),110,110);
+            %                           calcCap_semianalytisch3D(     s           ,     R_WK,     R_R ,   B_R,   B,    R_L ,   epsilon_r  ,     h_0           ,     a           ,     b           , nt, np)
+            C_semi_i(runSemi) = BearImp.calcCap_semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_ri, S.epsilon_Oel,temp_h_0(1,runSemi),temp_a(1,runSemi),temp_b(1,runSemi),110,110,'CmethodDeformedArea',C.method.deformedArea,'r2smallCriterion',C.method.r2smallCriterion);
+            C_semi_a(runSemi) = BearImp.calcCap_semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ra, S.epsilon_Oel,temp_h_0(2,runSemi),temp_a(2,runSemi),temp_b(2,runSemi),110,110,'CmethodDeformedArea',C.method.deformedArea,'r2smallCriterion',C.method.r2smallCriterion);
         end
     end
 
@@ -284,130 +354,6 @@ function runSemianalytisch3D(indices)
     C.C_out(2,:,posBall_conductive)=C_semi_a;
 end
 
-function C_Zusatz=semianalytisch3D(s,R_WK,R_R,B_R,B,R_L,epsilon_r,varargin)
-% Semianalytische Berechnung der Kapazität von unbelasteten Wälzkörpern
-% sowie Randbereichen belasteter WK am Innen und Außenring
-% Autor: Steffen Puchtler
-% Herleitung: Masterthesis "Optimierung des Berechnungsverfahrens für die
-% elektrische Kapazität von EHD-Kontakten unter Berücksichtigung des realen
-% elektrischen Feldes", S. 23-26
-%
-% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r)
-%   belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b)
-
-% Diskretisierung nicht Default:
-% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,0,0,nt,np)
-%   belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b,nt,np)
-
-% Inputparameter:
-% S     Minimaler Abstand zwischen Wälzkörper und Laufbahn  m
-% R_WK  Wälzkörperradius                                    m
-% f     Schmiegung (R_R / D_WK)                             -
-% B_R   Rillenbreite                                        m
-% B     Breite des Lagers                                   m
-% R_L   Radius der Laufbahn                                 m  !!! R_L < 0 am Innenring
-% epsilon_r     Relative Permittivität des Schmierstoffs    -
-% a,b   Halbachsen der Hertzschen Flächen                   m
-% nt,np Anzahl der Diskretisierngen in theta bzw. phi       -
-
-% profile on
-
-    switch nargin
-        case 7
-            nt = 60; np = 60;
-            belastet = false;
-        case 9
-            nt = 110; np = 110;
-            a  = varargin{1};
-            b  = varargin{2};
-            belastet = a~=0 && b~=0;
-        case 11
-            a  = varargin{1};
-            b  = varargin{2};
-            nt = varargin{3};
-            np = varargin{4};
-            belastet = a~=0 && b~=0;
-        otherwise
-            warning('The number of input arguments is wrong! Please check the input of semianalytisch3D.')
-    end
-    
-    DeltaR = R_R-R_WK-s;
-    R_T = R_L - R_R;
-    Dz = R_R-sqrt(R_R^2-B_R^2/4);
-    r_yz = @(phi) (R_WK-R_L+s).*cos(phi) + sqrt((R_L-Dz).^2-(R_WK-R_L+s).^2.*sin(phi).^2)*sign(R_L);
-    if belastet
-        Theta_0 = @(phi) a*sqrt(max(1/R_WK^2-phi.^2/b^2,0));
-    else
-        Theta_0 = @(~) 0;
-    end
-    Theta_1 = @(phi) atan(B_R/2./r_yz(phi));
-    Theta_2 = @(phi) atan(B  /2./r_yz(phi));
-    if R_L > 0
-        phi_1 = pi/2;
-    else
-        phi_1 = asin(R_L/(R_L-s-R_WK))*0.99;
-    end
-    p = linspace(0,phi_1,np);
-    P = repmat(p',1,nt);
-    T_t = zeros(size(P));
-    for ii = 1:np
-        T_t(ii,:) = linspace(Theta_0(p(ii)),Theta_1(p(ii)),nt);
-    end
-
-    r = zeros(size(T_t));
-    for ii = 1:size(T_t,1)
-        for jj = 1:size(T_t,2)
-            theta = T_t(ii,jj);
-            phi = P(ii,jj);
-            fun = @(r) (R_T+R_R*sqrt(1-(r.*sin(theta)/R_R).^2)).^2 - (r.*sin(phi).*cos(theta)).^2 - (DeltaR+R_T+r.*cos(phi).*cos(theta)).^2;
-%             options=optimset('Display','off','FunValCheck','off','OutputFcn',[],'PlotFcns',[],'TolX',eps);
-            r(ii,jj) = fzero(fun,R_WK);
-        end
-    end
-    clear phi theta
-    invalid = r(:)>1;
-    r2small = r(:)/R_WK<1.00125;
-    if any(invalid)
-        warning('%i Entries of r have been removed',sum(invalid))
-        C_Zusatz=nan;
-        return
-    elseif  any(r2small)
-%         fprintf('%i Entries of r have been removed \n',sum(r2small))
-        r(r2small)=inf;
-    
-    end
-    
-    integrand_t = R_WK^2./(r-R_WK).*cos(T_t);
-    firstInteg = zeros(1,np);
-    for ii = 1:np
-        firstInteg(ii) = trapz(T_t(ii,:),integrand_t(ii,:));
-    end
-    C_Rille = 4*epsilon_r * obj.epsilon_0 * trapz(P(:,1)',firstInteg);
-
-
-    T_k = zeros(size(P));
-    for ii = 1:np
-        T_k(ii,:) = linspace(Theta_1(p(ii)),Theta_2(p(ii)),nt);
-    end
-    integrand_k = R_WK^2./((r_yz(P)-s)./cos(T_k)-R_WK).*cos(T_k);
-    firstInteg = zeros(1,np);
-    for ii = 1:np
-        firstInteg(ii) = trapz(T_k(ii,:),integrand_k(ii,:));
-    end
-    C_Rand = 4*epsilon_r * obj.epsilon_0 * trapz(P(:,1)',firstInteg);
-
-    C_Zusatz=C_Rille + C_Rand;
-%     if nargout == 1
-%         varargout = {C_Rille + C_Rand};
-%     elseif nargout == 2
-%         varargout = {C_Rille,C_Rand};
-%     end
-
-% profile off
-% profile viewer
-
-end
-
 function seperateVecs=splitVec(vec,status)
     %da die parfor-Schleife nur Vektoren bearbeiten kann, die in Ganzzahl
     %Schritten auf- oder absteigen, wird in dieser Funktion jeder Vektor in
diff --git a/@BearImp/calcCap_puchtler2025.m b/@BearImp/calcCap_puchtler2025.m
new file mode 100644
index 0000000000000000000000000000000000000000..fca459c99733f135f3ab9ed800775846be5434f0
--- /dev/null
+++ b/@BearImp/calcCap_puchtler2025.m
@@ -0,0 +1,109 @@
+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
+        R_RE
+        R_R
+        R_ZL
+        R_ZR
+        B
+        R_L
+        epsilon_r
+        alpha
+        h_0
+        a = 0
+        b = 0
+        method = struct('deformedArea','neglect','outlet','oil','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,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);
+
+    if R_L > 0
+        phi_1 = pi/2;
+    else
+        phi_1 = asin(R_L/(R_L-s-R_RE))*0.9;
+    end
+    
+    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)));
+    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)),(a/R_RE+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)),(a/R_RE+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 = 0;%a/R_RE/sqrt(1+a^2/b^2);
+        if Theta_0l(phi_1) - (a/R_RE+phi_1) > Theta_1l(phi_1)
+            phi_1l = phi_1;
+        else
+            phi_1l = fzero(@(phi) real(Theta_0l(phi) - (a/R_RE+phi) - Theta_1l(phi)),0); % True value is always real. But on iterations, imaginary solutions might occur, resulting in 
+        end
+        if Theta_0r(phi_1) + (a/R_RE+phi_1) < Theta_1r(phi_1)
+            phi_1r = phi_1;
+        else
+            phi_1r = fzero(@(phi) real(Theta_0r(phi) + (a/R_RE+phi) - Theta_1r(phi)),0);
+        end
+        % 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)));
+    
+    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);
+    
+    switch method.outlet
+        case 'oil'
+            absTol = 50e-4; % 50e-4 *2 corresponds to 0.2 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 = 100e-4; % 100e-4 corresponds to 0.2 pF with a permittivity of 2.2
+            fun_grooveOut = @(phi,theta) R_RE^2./(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*2); 
+            C_grooveLOut = 8.854187e-12             * integral2(fun_grooveOut, phi_0,phi_1,Theta_0lIn ,Theta_0lOut,'AbsTol',absTol*2);
+            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,method)
+    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 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 method.deformedArea
+        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/calcCap_semianalytisch3D.m b/@BearImp/calcCap_semianalytisch3D.m
new file mode 100644
index 0000000000000000000000000000000000000000..ed6f41e560e6f7a8d7c878733ab2fa6062211742
--- /dev/null
+++ b/@BearImp/calcCap_semianalytisch3D.m
@@ -0,0 +1,129 @@
+function C_out = calcCap_semianalytisch3D(s,R_WK,R_R,B_R,B,R_L,epsilon_r,h_0,a,b,nt,np,options)
+% Semianalytische Berechnung der Kapazität von unbelasteten Wälzkörpern
+% sowie Randbereichen belasteter WK am Innen und Außenring
+% Autor: Steffen Puchtler
+% Herleitung: Masterthesis "Optimierung des Berechnungsverfahrens für die
+% elektrische Kapazität von EHD-Kontakten unter Berücksichtigung des realen
+% elektrischen Feldes", S. 23-26
+%
+% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r)
+%   belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b)
+
+% Diskretisierung nicht Default:
+% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,0,0,nt,np)
+%   belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b,nt,np)
+
+% Inputparameter:
+% S     Minimaler Abstand zwischen Wälzkörper und Laufbahn  m
+% R_WK  Wälzkörperradius                                    m
+% f     Schmiegung (R_R / D_WK)                             -
+% B_R   Rillenbreite                                        m
+% B     Breite des Lagers                                   m
+% R_L   Radius der Laufbahn                                 m  !!! R_L < 0 am Innenring
+% epsilon_r     Relative Permittivität des Schmierstoffs    -
+% a,b   Halbachsen der Hertzschen Flächen                   m
+% nt,np Anzahl der Diskretisierngen in theta bzw. phi       -
+
+% profile on
+
+    arguments
+        s
+        R_WK
+        R_R
+        B_R
+        B
+        R_L
+        epsilon_r
+        h_0
+        a = 0
+        b = 0
+        nt = 110
+        np = 110
+        options.CmethodDeformedArea = 'neglect'
+        options.r2smallCriterion = 'r<1.00125*R_RE'
+    end
+
+    belastet = a~=0 && b~=0;
+
+    if ~belastet && nargin <= 9
+        nt = 60; np = 60;
+    end
+    
+    DeltaR = R_R-R_WK-s;
+    R_T = R_L - R_R;
+    Dz = R_R-sqrt(R_R^2-B_R^2/4);
+    r_yz = @(phi) (R_WK-R_L+s).*cos(phi) + sqrt((R_L-Dz).^2-(R_WK-R_L+s).^2.*sin(phi).^2)*sign(R_L);
+    if belastet
+        Theta_0 = @(phi) a*sqrt(max(1/R_WK^2-phi.^2/b^2,0));
+    else
+        Theta_0 = @(~) 0;
+    end
+    Theta_1 = @(phi) atan(B_R/2./r_yz(phi));
+    Theta_2 = @(phi) atan(B  /2./r_yz(phi));
+    if R_L > 0
+        phi_1 = pi/2;
+    else
+        phi_1 = asin(R_L/(R_L-s-R_WK))*0.99;
+    end
+    p = linspace(0,phi_1,np);
+    P = repmat(p',1,nt);
+    T_t = zeros(size(P));
+    for ii = 1:np
+        T_t(ii,:) = linspace(Theta_0(p(ii)),Theta_1(p(ii)),nt);
+    end
+
+    r = zeros(size(T_t));
+    for ii = 1:size(T_t,1)
+        for jj = 1:size(T_t,2)
+            theta = T_t(ii,jj);
+            phi = P(ii,jj);
+            fun = @(r) (R_T+R_R*sqrt(1-(r.*sin(theta)/R_R).^2)).^2 - (r.*sin(phi).*cos(theta)).^2 - (DeltaR+R_T+r.*cos(phi).*cos(theta)).^2;
+%             options=optimset('Display','off','FunValCheck','off','OutputFcn',[],'PlotFcns',[],'TolX',eps);
+            r(ii,jj) = fzero(fun,R_WK);
+        end
+    end
+    clear phi theta
+    switch options.r2smallCriterion
+        case 'r<1.00125*R_RE'
+            r2small = r(:)/R_WK<1.00125;
+        case 'r<h_0+R_RE'
+            r2small = r(:)<h_0+R_WK;
+    end
+%   fprintf('%i Entries of r have been changed \n',sum(r2small))
+    switch options.CmethodDeformedArea
+        case 'neglect'
+            r(r2small) = inf;
+        case 'filmThickness'
+            r(r2small) = h_0+R_WK;
+    end
+    
+    integrand_t = R_WK^2./(r-R_WK).*cos(T_t);
+    firstInteg = zeros(1,np);
+    for ii = 1:np
+        firstInteg(ii) = trapz(T_t(ii,:),integrand_t(ii,:));
+    end
+    C_Rille = 4*epsilon_r * 8.854187817e-12 * trapz(P(:,1)',firstInteg);
+
+
+    T_k = zeros(size(P));
+    for ii = 1:np
+        T_k(ii,:) = linspace(Theta_1(p(ii)),Theta_2(p(ii)),nt);
+    end
+    integrand_k = R_WK^2./((r_yz(P)-s)./cos(T_k)-R_WK).*cos(T_k);
+    firstInteg = zeros(1,np);
+    for ii = 1:np
+        firstInteg(ii) = trapz(T_k(ii,:),integrand_k(ii,:));
+    end
+    C_Rand = 4*epsilon_r * 8.854187817e-12 * trapz(P(:,1)',firstInteg);
+
+    C_out = C_Rille + C_Rand;
+%     if nargout == 1
+%         varargout = {C_Rille + C_Rand};
+%     elseif nargout == 2
+%         varargout = {C_Rille,C_Rand};
+%     end
+
+% profile off
+% profile viewer
+
+end
\ No newline at end of file
diff --git a/@BearImp/calcClear.m b/@BearImp/calcClear.m
index 1a32a7589b9d8f1f3ae3189bae47906e9801f0c6..2423ccfd650edb559a2b514aff8ec6c5d749acba 100644
--- a/@BearImp/calcClear.m
+++ b/@BearImp/calcClear.m
@@ -7,6 +7,11 @@ function calcClear(obj)
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler, Julius van der Kuip
 
+if ~isscalar(obj)
+    executeAllObj(obj,@calcClear)
+    return
+end
+
 %% Parameter prüfen
 assert(obj.up2date.L,'Lager nicht gesetzt')
 L = obj.L; F_r = obj.F_r;
diff --git a/@BearImp/calcFilm.m b/@BearImp/calcFilm.m
index b734c6fea58d2c4852fa26412dc14b815189cce9..5c6e38770cd4c61aa9751b31975345aa4d6a73c9 100644
--- a/@BearImp/calcFilm.m
+++ b/@BearImp/calcFilm.m
@@ -1,4 +1,4 @@
-function varargout = calcFilm(obj, options)
+function h_0 = calcFilm(obj, options)
 %2.5 Schmierfilmdicke H = H(S,T,G,B,T_Oil,n)
 %Berechnet die Hertz'schen Flächen sowie die Schmierfilmdicken aller
 %Einzelkontakt (innen und außen für alle belasteten Wälzkörper)
@@ -35,7 +35,7 @@ end
 L=obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; T_Oil = obj.T_Oil; omega = obj.omega;
 H = struct;
 
-method = possibleMethods.addDefault(obj.method).H;
+H.method = obj.method.H;
 
 %% Berechnung
   
@@ -71,47 +71,67 @@ H.W =  Q_temp./(E_red_tmp .* R_x_temp.^2);
 H.L = H.G .* H.U.^0.25;
 H.M = H.W.*(2.*H.U).^(-3/4);
 
-switch method
-    case {'Hamrock/Dowson','Hamrock/Dowson_withoutThermCorr'}
+switch H.method.filmThicknessFormula
+    case 'Hamrock/Dowson'
         H.H_0   = 2.69 * H.G.^0.53 .* H.U.^0.67 .* H.W.^-0.067 .* (1-0.61*exp(-0.73.*G.k));
         H.H_min = 3.63 * H.G.^0.49 .* H.U.^0.68 .* H.W.^-0.073 .* (1-     exp(-0.68.*G.k));
-        H.h_0raw   = H.H_0 .*R_x_temp; % without correction factors
+        H.h_0   = H.H_0  .*R_x_temp;
         H.h_min = H.H_min.*R_x_temp;
-        H.L_therm = T.eta_0 .* T.alpha_etaT .* H.u.^2  ./  (4*S.lambda);
-        if strcmp(method,'Hamrock/Dowson')
-            H.C_korr  = 3.94  ./  (3.94 + H.L_therm.^0.62);
-        else
-            H.C_korr = 1;
-        end
-        H.h_0   = H.h_0raw  .*H.C_korr;
-        H.h_minth = H.h_min.*H.C_korr;
         
     case 'Moes'
-        
         H.lamda = R_x_temp./R_y_temp;
         H.N = H.W .* H.lamda.^0.5 .* H.U.^(-0.75);
         
-        H.C_RI = 145 * (1+0.796 * H.lamda .^(14/15)) .^(-15/7);
+        H.C_RI = 145  * (1+0.796 * H.lamda .^(14/15)) .^(-15/7);
         H.C_RP = 1.29 * (1+0.691 * H.lamda) .^(-2/3);
-        H.C_EI = 3.18 * (1+0.006 * log( H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-14/25);
-        H.C_EP = 1.48 * (1+0.006 * log( H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-7/20);
+        H.C_EI = 3.18 * (1+0.006 * log(H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-14/25);
+        H.C_EP = 1.48 * (1+0.006 * log(H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-7/20);
         
         H.H_RI = H.C_RI .* H.N.^(-2);
         H.H_RP = H.C_RP .* H.L.^(2/3);
         H.H_EI = H.C_EI .* H.N.^(-2/15);
         H.H_EP = H.C_EP .* H.N.^(-1/12) .* H.L.^(3/4);
-        H.s = 3/2.* (1+ exp(-1.2 .* H.H_EI ./ H.H_RI));
+        H.s = 3/2.* (1 + exp(-1.2 .* H.H_EI ./ H.H_RI));
         
         H.H_0    = ((H.H_RI.^(3/2) + (H.H_EI.^(-4) + 0.1 .* H.lamda.^4).^(-3/8)).^(2*H.s./3) + (H.H_RP .^(-8) + H.H_EP .^(-8) ).^(-H.s ./8)) .^(1./H.s);
         H.h_0    = H.H_0 .*R_x_temp .* H.U.^(0.5);
 end
 
+if H.method.thermalCorrection == 'on'
+    H.L_therm = T.eta_0 .* T.alpha_etaT .* H.u.^2  ./  (4*S.lambda);
+    H.C_thCorr  = 3.94  ./  (3.94 + H.L_therm.^0.62);
+    H.h_0     = H.h_0.*H.C_thCorr;
+    if isfield(H,'h_min')
+        H.h_min = H.h_min .*H.C_thCorr;
+    end
+end
+
+if H.method.starvationCorrection == 'on' % according to SKF according to Schneider 2021 !!! In my opinion a factor for the friction not for lubrication film thickness
+    H.K_rs = 3e-8; % 3e-8 for low level oil bath/oil jet lubrication, 6e−8 for grease/oil-air lubrication
+    H.K_Z = 3.1; % 3.1 for single row deep groove ball bearings, 5.1 for cylindrical roller bearings with cage
+    T.nu_0 = T.eta_0 ./ T.rho(T_Oil); %(kin. viskosity at ambient pressure at oil temperature)
+    n = omega/2/pi*60; % in rpm
+    H.C_rsCorr  = 1./exp(H.K_rs*T.nu_0*1e6*n*(L.d+L.D)*1e3*sqrt(H.K_Z/(2*(L.D-L.d)*1e3)));
+    H.h_0     = H.h_0.*H.C_rsCorr;
+    if isfield(H,'h_min')
+        H.h_min = H.h_min .*H.C_rsCorr;
+    end
+end
+
+if H.method.roughnessCorrection == 'on' % according to Kreil 2009 (S. 79)
+    H.t_r = -0.0992; % -0.0992 for longitudinal grinding, −1.0600 for transversal grinding
+    H.m_r = 0.3 * log(H.u+55) + H.t_r;
+    H.Rz = sqrt(L.Rz_RE^2 + [L.Rz_IR;L.Rz_OR].^2);
+    H.Delta_h_r = H.m_r .* H.Rz;
+    H.h_0 = H.h_0 - H.Delta_h_r;
+end
+
 if isPreCalc
-    varargout{1} = H.h_0;
+    h_0 = H.h_0;
     return
 end
 
-if strcmp(possibleMethods.addDefault(obj.method).B,'static')
+if B.method == 'static'
     assert(ndims(H.h_0) <= 2,'h_0 must be two dimensional to calculate film for static load distribution') % probably obsolete
     B.s=zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
 
@@ -123,26 +143,26 @@ end
 
 % Überprüfen, ob die Schmierfilmberechnung innerhalb des vorgegeben
 % Bereichs liegt
-switch method
-case {'Hamrock/Dowson','Hamrock/Dowson_withoutThermCorr'}
-    rangeLoadParM = [25,500]; % Quelle: Non-Dimensional Groups, Marian 2020, Fig7
-    rangeViscParL = [5,15];
-    
-    if strcmp(possibleMethods.addDefault(obj.method).B,'dynamic')
-       warning('It is not recommended to use Hamrock/Dowson as lubricant film calculation when dynamic load calculation is selected, because the loads become very small outside the load range.') 
-    end
-case 'Moes'
-    rangeLoadParM = [5,1000];
-    rangeViscParL = [1,25];
+switch H.method.filmThicknessFormula
+    case 'Hamrock/Dowson'
+        rangeLoadParM = [25,500]; % Quelle: Non-Dimensional Groups, Marian 2020, Fig7
+        rangeViscParL = [5,15];
+        
+        if B.method == 'dynamic' && obj.displayWarnings
+           warning('It is not recommended to use Hamrock/Dowson as lubricant film calculation when dynamic load calculation is selected, because the loads become very small outside the load range.') 
+        end
+    case 'Moes'
+        rangeLoadParM = [5,1000];
+        rangeViscParL = [1,25];
 end
         
-if max(rangeLoadParM) < max(H.M,[],'all') || min(rangeLoadParM) > min(H.M,[],'all')
+if (max(rangeLoadParM) < max(H.M,[],'all') || min(rangeLoadParM) > min(H.M,[],'all')) && obj.displayWarnings
 warning(['The force for the lubricant film thickness calculation (%.2f to %.2f) is outside the permissible range (%.0f to %.0f), ' ...
          'which means that the calculated lubricant film is no longer meaningful!'], ...
          min(H.M,[],'all'), max(H.M,[],'all'), min(rangeLoadParM), max(rangeLoadParM))
 end
 
-if max(rangeViscParL) < max(H.L,[],'all') || min(rangeViscParL) > min(H.L,[],'all')
+if (max(rangeViscParL) < max(H.L,[],'all') || min(rangeViscParL) > min(H.L,[],'all')) && obj.displayWarnings
 warning(['The viscosity parameter for the lubricant film thickness calculation (%.2f to %.2f) is outside the permissible range (%.0f to %.0f), ' ...
          'which means that the calculated lubricant film is no longer meaningful!'], ...
          min(H.L,[],'all'), max(H.L,[],'all'), min(rangeViscParL), max(rangeViscParL))
diff --git a/@BearImp/calcGeo.m b/@BearImp/calcGeo.m
index ba88b490a5035fb417d5d676ab1088422be533fb..2034a93178273159a13d6bdf7463798782e5d239 100644
--- a/@BearImp/calcGeo.m
+++ b/@BearImp/calcGeo.m
@@ -1,4 +1,4 @@
-function calcGeo(obj)
+function calcGeo(obj,options)
 %2.3 LagerGeometrie G = G(L,R,T_Oil)
 %Berechnet verschiedene Größen der Lagergeometrie wie z.B.
 %Berührungswinkel, Rollradien, Elliptizitätsparameter der Hertz'schen
@@ -7,31 +7,45 @@ function calcGeo(obj)
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler, Julius van der Kuip
 
+arguments
+    obj
+    options.Delta_Tia = 5 % Temperaturerhöhung des Innenrings gegenüber des Außenrings
+    options.alpha = []
+end
+
+if ~isscalar(obj)
+    executeAllObj(obj,@calcGeo)
+    return
+end
+
 %% Parameter prüfen
 assert(obj.up2date.L,'Lager nicht gesetzt')
 assert(obj.up2date.R,'Lagerspiel noch nicht berechnet')
-L = obj.L; R = obj.R; T_Oil = obj.T_Oil;
-G = struct; 
+L = obj.L; R = obj.R; T_Oil = obj.T_Oil; F_a = obj.F_a;
+G = struct; G.method = obj.method.G;
 
 %% Berechnung
-G.alpha = 0; 
-
-G.Delta_Ta_norm = T_Oil - 20; %Temperaturerhöhung der Welle bezogen auf die Normtemperatur von 20°C (Innenring heißer als Außenring)
-G.Delta_Ti_norm = G.Delta_Ta_norm + 5; %Temperaturerhöhung des Gehäuses bezogen auf die Normtemperatur von 20°C #######Idee: 5°C durch festlegbare Variable ersetzen
+G.Delta_Tia = options.Delta_Tia;                 % Temperaturerhöhung des Innenrings gegenüber des Außenrings
+G.Delta_Ta_norm = T_Oil - 20;                    % Temperaturerhöhung der Welle bezogen auf die Normtemperatur von 20°C (Innenring heißer als Außenring)
+G.Delta_Ti_norm = G.Delta_Ta_norm + G.Delta_Tia; % Temperaturerhöhung des Gehäuses bezogen auf die Normtemperatur von 20°C
+G.Delta_Tm_norm = (G.Delta_Ta_norm + G.Delta_Ti_norm)/2;
 
 G.E_red(1) = 2  ./  ((1-L.RE_A.nu.^2)./L.RE_A.E + (1-L.ring.nu.^2)/L.ring.E);
-G.D_RE(1) = L.D_RE * (1 + L.RE_A.alpha_t * ((G.Delta_Ti_norm + G.Delta_Ta_norm) / 2));
+G.D_RE(1) = L.D_RE * (1 + L.RE_A.alpha_t * G.Delta_Tm_norm);
 
-if ~L.allBallsSameMaterial % Materialparameter für alle WK Materialien berechnet
+if ~L.allBallsSameMaterial % Materialparameter für alle WK Materialien berechnen
     G.E_red(2) = 2  ./  ((1-L.RE_B.nu.^2)./L.RE_B.E + (1-L.ring.nu.^2)/L.ring.E);
-    G.D_RE(2) = L.D_RE * (1 + L.RE_B.alpha_t * ((G.Delta_Ti_norm + G.Delta_Ta_norm) / 2));
+    G.D_RE(2) = L.D_RE * (1 + L.RE_B.alpha_t * G.Delta_Tm_norm);
 end
 
 G.D = (L.D - R.e_ar) * (1 + L.ring.alpha_t * G.Delta_Ta_norm);
 G.d = (L.d + R.e_ir) * (1 + L.ring.alpha_t * G.Delta_Ti_norm);
+G.R_RE = G.D_RE/2;
 
 G.d_m = (G.d + G.D)/2;
-G.R_RE = G.D_RE/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; 
@@ -40,33 +54,65 @@ G.b_la = b_la*(1+L.ring.alpha_t*G.Delta_Ta_norm);
 
 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_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.b = L.f_i + L.f_a - 1; % Harris B
+if ~isempty(options.alpha)
+    G.alpha = options.alpha;
+elseif ~F_a || G.method.alphaForRaceRadius == 'zero'
+    G.alpha = 0;
+elseif G.method.alphaForRaceRadius == 'estimate'
+    if G.P_d > 0
+        G.alpha_0 = acos(1 - G.P_d/(2*G.D_RE*G.b));
+    else
+        G.alpha_0 = 0;
+    end
+    G.k_ax_est = interp1([0 0.02 0.08 0.14 0.2],[0 50   250  500  760]*1e3,G.b,'spline')*6894.757; % Harris.2001 Fig. 7.5, p. 247
+    alpha = pi/3; % Stable start value
+    numberOfIterations = 1;
+    while true
+        cosAlphaRatioTerm = (cos(G.alpha_0)/cos(alpha) - 1);
+        alpha_next = alpha + (            F_a/(L.Z*L.D_RE^2*G.k_ax_est) - sin(alpha)*cosAlphaRatioTerm^1.5               )  /  ... Harris.2001 (7.34), p. 247
+                             ( cos(alpha)*cosAlphaRatioTerm^1.5  +  1.5*tan(alpha)^2*cosAlphaRatioTerm^0.5*cos(G.alpha_0) );
+        if abs(alpha-alpha_next) <= 100*eps(alpha)
+            break
+        elseif numberOfIterations >= 1000
+            error('Estimation of mounted contact angle did not converge')
+            %break
+        else
+            alpha = alpha_next;
+            numberOfIterations = numberOfIterations + 1;
+        end
+    end
+    G.alpha = alpha_next;
+end
 
 G.R_bi = (G.d + G.b_li) / 2; %Minimaler innerer Rillenradius
-G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius
+G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius                 %  unabh. von alpha
 
-G.sigmarho_i = (2./G.R_RE + 1./G.R_bi - 1./G.R_li);
-G.sigmarho_a = (2./G.R_RE - 1./G.R_ba - 1./G.R_la);
+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.costau_i = (1./G.R_bi + 1./G.R_li) ./ G.sigmarho_i;
-G.costau_a = (-1./G.R_ba + 1./G.R_la) ./ G.sigmarho_a;
+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
+G.R_si = 1 ./ G.sigmarho_i;                                                                     % in calcFilm
+G.R_sa = 1 ./ G.sigmarho_a;                                                 % abh. von alpha    % in calcFilm
 
-G.R_xi = 1 ./ ( 1./G.R_RE + 1./G.R_bi );
-G.R_xa = 1 ./ ( 1./G.R_RE - 1./G.R_ba );
-G.R_x  = [G.R_xi;G.R_xa];
+G.R_xi = 1 ./ ( 1./G.R_RE + 1./G.R_ri );
+G.R_xa = 1 ./ ( 1./G.R_RE - 1./G.R_ra );
+G.R_x  = [G.R_xi;G.R_xa];                                                   % abh. von alpha    % in calcFilm & calcCap
 
 G.R_yi = 1 ./ ( 1./G.R_RE - 1./G.R_li );
 G.R_ya = 1 ./ ( 1./G.R_RE - 1./G.R_la );
 G.R_y  = [G.R_yi;G.R_ya];
 
-G.R_si = 1 ./ G.sigmarho_i;
-G.R_sa = 1 ./ G.sigmarho_a;
-
-G.gamma = G.D_RE  *cos(G.alpha) ./ G.d_m;
+G.gamma = G.D_RE  *cos(G.alpha) ./ G.d_m;                                   % abh. von alpha   % in calcFilm
 
 integralF = @(phi,k) (1 - (1-1/k.^2).*sin(phi).^2).^-0.5;
 integralE = @(phi,k) (1 - (1-1/k.^2).*sin(phi).^2).^0.5;
@@ -74,24 +120,24 @@ integralE = @(phi,k) (1 - (1-1/k.^2).*sin(phi).^2).^0.5;
 fraktF = @(k) integral( @(phi) integralF(phi,k),0,pi/2 );
 fraktE = @(k) integral( @(phi) integralE(phi,k),0,pi/2 );
 
-G.F_i = ((G.gamma/(1-G.gamma)+G.D_RE/(2*G.R_li))/(2+G.gamma/(1-G.gamma)-G.D_RE/(2*G.R_li)));
-G.F_a = ((-G.gamma/(G.gamma+1)+G.D_RE/(2*G.R_la))/(2-G.gamma/(1+G.gamma)-G.D_RE/(2*G.R_la)));
+G.F_i = (( G.gamma/(1-G.gamma)+G.D_RE/(2*G.R_li))/(2+G.gamma/(1-G.gamma)-G.D_RE/(2*G.R_li)));
+G.F_a = ((-G.gamma/(1+G.gamma)+G.D_RE/(2*G.R_la))/(2-G.gamma/(1+G.gamma)-G.D_RE/(2*G.R_la)));     % abh. von alpha
 
 k_i_fkt = @(k_i) (1-2/(k_i^(2)-1)*((fraktF(k_i)/fraktE(k_i))-1)-G.F_i);
-k_a_fkt = @(k_a) (1-2/(k_a^(2)-1)*((fraktF(k_a)/fraktE(k_a))-1)-G.F_a);
+k_a_fkt = @(k_a) (1-2/(k_a^(2)-1)*((fraktF(k_a)/fraktE(k_a))-1)-G.F_a);     % abh. von alpha
 
-G.k_i = fzero(k_i_fkt,5);
-G.k_a = fzero(k_a_fkt,5);
-G.k = [G.k_i;G.k_a];
+G.k_i = fzero(k_i_fkt,5);                                                                       % in calcFilm
+G.k_a = fzero(k_a_fkt,5);                                                                       % in calcFilm
+G.k = [G.k_i;G.k_a];                                                        % abh. von alpha    % in calcFilm & calcCap
 
 G.fraktF_i = fraktF(G.k_i);
 G.fraktF_a = fraktF(G.k_a);
-G.fraktE_i = fraktE(G.k_i);
-G.fraktE_a = fraktE(G.k_a);
+G.fraktE_i = fraktE(G.k_i);                                                                     % in calcFilm
+G.fraktE_a = fraktE(G.k_a);                                                 % abh. von alpha    % in calcFilm
 
-G.K_pi = pi.*G.k_i.*G.E_red.*(G.fraktE_i./(4.5*G.sigmarho_i.*G.fraktF_i.^3)).^0.5;
-G.K_pa = pi.*G.k_a.*G.E_red.*(G.fraktE_a./(4.5*G.sigmarho_a.*G.fraktF_a.^3)).^0.5;
-G.K_p= 1./(1./G.K_pi.^(2/3)+1./G.K_pa.^(2/3)).^1.5;
+G.K_pi = pi.*G.k_i.*G.E_red.*(G.fraktE_i./(4.5*G.sigmarho_i.*G.fraktF_i.^3)).^0.5;              % in calcLoad
+G.K_pa = pi.*G.k_a.*G.E_red.*(G.fraktE_a./(4.5*G.sigmarho_a.*G.fraktF_a.^3)).^0.5;              % in calcLoad
+G.K_p= 1./(1./G.K_pi.^(2/3)+1./G.K_pa.^(2/3)).^1.5;                         % abh. von alpha    % in calcLoad
 
 G.Deltapsi = 2*pi./L.Z;
 G.psi_ = (0:L.Z-1) .* G.Deltapsi;
diff --git a/@BearImp/calcLoad.m b/@BearImp/calcLoad.m
index fb2d85ece2c8dd94742d6b7de53af58b89de4e88..2fde9309998deaadf0dc996436229bf2c07e2adb 100644
--- a/@BearImp/calcLoad.m
+++ b/@BearImp/calcLoad.m
@@ -11,7 +11,7 @@ assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht
 L = obj.L; G = obj.G; F_r = obj.F_r; F_a = obj.F_a; AddOn = obj.AddOn; psi=obj.psi;
 B = struct; 
 
-B.method = possibleMethods.addDefault(obj.method).B;
+B.method = obj.method.B;
 %% Berechnung
     
 if ~isnan(obj.psi_calc) 
@@ -20,7 +20,7 @@ end
 
 % delta_intp1 (the distance between the raceways) needs to be calculated
 % before delta_s in the dynamic case
-if strcmp(B.method,'dynamic')
+if B.method == 'dynamic'
     n_ir = obj.omega * 60 / (2*pi);
     n_ar = 0;
     B.v_cage = pi/60*G.d_m.*((1-cos(G.alpha)./(G.d_m./G.D_RE)).*n_ir./2+(1+cos(G.alpha)./(G.d_m./G.D_RE)).*n_ar./2); %Umfangsgeschwindigkeit des Käfigs [m/s] Brändlein S.88; Dowson u_i/o/c S.58
@@ -80,7 +80,7 @@ else
     end
 end
 
-if strcmp(B.method,'dynamic')
+if B.method == 'dynamic'
    assert(any(max(B.Q_intp1) > max(B.Q) & min(B.Q_intp1) < min(B.Q)), 'The interpolated lubricant film thicknesses are extrapolated outside the calculated range! This leads to inaccurate results.')
 end
 
@@ -93,7 +93,7 @@ for posBall_conductive=1:L.numberOfConductiveBalls
                                 (round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z)))];
 
     % intersectionBall (intersection of Ball and inner/outer racerway)
-    if strcmp(B.method,'static') % static assumption: Ball in middle of gap
+    if B.method == 'static' % static assumption: Ball in middle of gap
         for IndPosBearing=1:B.numOfCagePositions
             deltaXYZ = deflectionForXYZ(L,G,L.IndBall_conductive(posBall_conductive),B.delta_s(:,IndPosBearing),B.psi_conductive(posBall_conductive,IndPosBearing));
             B.intersectionBall(1,IndPosBearing,posBall_conductive) = (deltaXYZ - G.D_RE(L.RE_order_num(posBall_conductive)))./2;
@@ -206,7 +206,7 @@ end
 
 function Q_total_y_perball_cont=totalReactionForcePerBallOnlyY(G,IndBall,delta_s,psi_run)
    delta_nury = (G.R_ba-G.R_RE-cos(psi_run(IndBall))*delta_s(1))/cos(atan(sin(psi_run(IndBall))*delta_s(1)/ ...
-                (G.R_ba-G.R_RE-cos(psi_run(IndBall))*delta_s(1))))-G.R_RE-G.R_bi - G.D_RE; %- G.D_RE richtig?
+                (G.R_ba-G.R_RE-cos(psi_run(IndBall))*delta_s(1))))-G.R_RE-G.R_bi - G.D_RE;
    Q_total_y_perball_cont = abs((delta_nury<0)*(delta_nury))^(3/2)*G.K_p;
 end
 
diff --git a/@BearImp/calcLub.m b/@BearImp/calcLub.m
index 79b871500ac4c612074cce797ea8076ce27ff8ce..663df16eb2ad4711c7fe946e3f48b41b4b4c9eae 100644
--- a/@BearImp/calcLub.m
+++ b/@BearImp/calcLub.m
@@ -5,12 +5,17 @@ function calcLub(obj)
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler
 
+if ~isscalar(obj)
+    executeAllObj(obj,@calcLub)
+    return
+end
+
 %% Parameter prüfen
 assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
 S = obj.S; T_Oil = obj.T_Oil;
 T = struct;
 
-T.method = possibleMethods.addDefault(obj.method).T;
+T.method = obj.method.T;
 
 %% Berechnung
 T.rho = @(theta) S.rho_15.*(1 - S.alpha_rho.*(theta-15));
@@ -23,17 +28,30 @@ switch T.method
         T.eta_00 = T.eta_040 - T.alpha_eta.*40;
         T.alpha_nu = (S.nu_100 - S.nu_40)/60;
         T.nu_0 = S.nu_40 - T.alpha_nu.*40;
-        T.eta_0 = T.eta_00 + T.alpha_eta.*T_Oil;
-        T.nu_38 = T.nu_0 + T.alpha_nu.*38;
+        T.eta = @(theta) T.eta_00 + T.alpha_eta.*theta;
+        T.nu  = @(theta) T.nu_0 + T.alpha_nu.*theta;
+    case 'DIN 51563'
+        assert(~any(isnan([S.nu_40 S.nu_100])),'Viscosity values not given for the selected lubricant. Choose method.T = ''Vogel''')
+        W_40  = log10(log10(S.nu_40 *1e6 + 0.8));                         % DIN 51563 (1)
+        W_100 = log10(log10(S.nu_100*1e6 + 0.8));                         % DIN 51563 (2)
+        m = (W_40-W_100)/(log10(100+273.15)-log10(40+273.15));            % DIN 51563 (3)
+        W_x = @(theta) m*(log10(40+273.15) - log10(theta+273.15)) + W_40; % DIN 51563 (4)
+        T.nu  = @(theta) (10.^(10.^W_x(theta)) - 0.8)*1e-6;
+        T.eta = @(theta) T.rho(theta).*T.nu(theta);
     case 'Vogel'
         assert(~any(isnan([S.B S.C S.K])),'Vogel-Parameters not given for the selected lubricant. Choose method.T = ''linear''')
-        T.eta = @(T) S.K .* exp(S.B ./ (T + S.C));
-        T.eta_0    = T.eta(T_Oil);
-        T.eta_040  = T.eta(  40 );
-        T.eta_0100 = T.eta( 100 );
-        T.nu_38 = T.eta(38) ./ T.rho(38);
+        T.eta = @(theta) S.K .* exp(S.B ./ (theta + S.C));
+        T.nu  = @(theta) T.eta(theta) ./ T.rho(theta);
 end
+
+if T.method == 'DIN 51563' || T.method == 'Vogel'
+    T.eta_040  = T.eta( 40);
+    T.eta_0100 = T.eta(100);
+end
+
 T.alpha_etaT = log(T.eta_040/T.eta_0100)/60;
+T.eta_0    = T.eta(T_Oil);
+T.nu_38 = T.nu(38);
 
 %% Attribute ändern
 obj.T = T;
diff --git a/@BearImp/set.m b/@BearImp/set.m
index e512d212daf3ab4c26c56055da0ccf0d8c3aac16..26c5a0633814bdc1d1ab5bf9f746a09e2c415a55 100644
--- a/@BearImp/set.m
+++ b/@BearImp/set.m
@@ -1,6 +1,6 @@
 function set(obj, name, value)
-%Weist einem Attribut eines BearImp-Objekts einen Wert zu. Funktioniert
-%auch mit BearImp-Arrays bzw. -Matrizen
+% Weist einem Attribut eines BearImp-Objekts einen Wert zu. Funktioniert
+% auch mit BearImp-Arrays bzw. -Matrizen
 
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler
@@ -24,13 +24,23 @@ function set(obj, name, value)
 end
 
 function setObjArrayPropToSingleValue(objArray, name, value)
+    propStack = strsplit(name,'.');
     for ii = 1:numel(objArray)
-        objArray(ii).(name) = value;
+        switch numel(propStack)
+            case 1
+                objArray(ii).(name) = value;
+            case 2
+                objArray(ii).(propStack{1}).(propStack{2}) = value;
+            case 3
+                objArray(ii).(propStack{1}).(propStack{2}).(propStack{3}) = value;
+            otherwise
+                error('Deeper property stacks are not yet implemented')
+        end
     end
 end
 
 function setObjArrayPropToValueArray(objArray, name, valueArray)
     for ii = 1:numel(objArray)
-        objArray(ii).(name) = valueArray(ii);
+        setObjArrayPropToSingleValue(objArray(ii), name, valueArray(ii))
     end
 end
\ No newline at end of file
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
new file mode 100644
index 0000000000000000000000000000000000000000..cbb0f3a045b7e11800ba8755031daa0573dd426c
--- /dev/null
+++ b/BearImpOptions.m
@@ -0,0 +1,377 @@
+classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & matlab.mixin.Copyable
+% Enthält die zu verwendenden Berechnungsmethoden der einzelnen
+% Berechnungsabschnitte (T, R, G, B, H, C), deren default-Methoden und 
+% eine Funktion zum überprüfen des Methodenstructs
+%
+% Limitations: Works for optionLevel 1:3 only (e.g. main.C.unloadedRE) and
+%              all options must be strings or char arrays
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Julius van der Kuip
+
+    properties (Hidden)
+        optionName
+        optionLevel = 1
+        optionChosen = 'default'
+        optionsLastChanged = {}
+        parent
+    end
+
+    properties (Dependent)
+        possibleSubOptions  (1,:) cell
+        hasSubSubOptions (1,1) logical
+    end
+
+    events
+        optionChanged
+    end
+
+    methods (Static,Hidden) 
+        function possibleOption = Possible
+            possibleOption.G.alphaForRaceRadius = {'zero','estimate'};
+            possibleOption.T = {'linear','Vogel','DIN 51563'};
+            possibleOption.B = {'static','dynamic'};
+            possibleOption.H.filmThicknessFormula = {'Hamrock/Dowson','Moes'};
+            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','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'};
+            possibleOption.C.roughnessCorrection = {'off','Napel'};
+            possibleOption.C.rhoRatio = {'Hartung','Dowson','Bradbury','Tait'};
+        end
+
+        function defaultOption = Default
+        % Gibt ein Rechenmethoden-Struct mit den Default-Methoden aus
+            defaultOption.G.alphaForRaceRadius = 'estimate';
+            defaultOption.T = 'Vogel';
+            defaultOption.B = 'dynamic';
+            defaultOption.H.filmThicknessFormula = 'Moes';
+            defaultOption.H.thermalCorrection = 'off';
+            defaultOption.H.starvationCorrection = 'off';
+            defaultOption.H.roughnessCorrection = 'off';
+            defaultOption.C.unloadedRE  = 'Puchtler2025';
+            defaultOption.C.outsideArea = 'Puchtler2025';
+            defaultOption.C.deformedArea= 'neglect';
+            defaultOption.C.outlet = 'air';
+            defaultOption.C.r2smallCriterion = 'r<h_0+R_RE';
+            defaultOption.C.k_vh_factor = 'off';
+            defaultOption.C.pressureDistribution = 'on';
+            defaultOption.C.roughnessCorrection = 'off';
+            defaultOption.C.rhoRatio = 'Tait';
+        end
+    end
+
+    methods (Static, Hidden) 
+        function subOptionsCell = subOptionsSetWith
+        % Gibt an, welche suboptions geändert werden sollen, wenn versucht
+        % wird, eine Option mit Unteroptionen direkt zu setzen.
+            subOptionsCell.H = {'filmThicknessFormula'};
+            subOptionsCell.C = {'unloadedRE','outsideArea'};
+        end
+
+        function allCategories = AllCategories
+            allCategories = [{'main'} fieldnames(BearImpOptions.Possible)'];
+            mainNames = fieldnames(BearImpOptions.Possible);
+            for ii = 1:numel(mainNames)
+                if isstruct(BearImpOptions.Possible.(mainNames{ii}))
+                    allCategories = [allCategories fieldnames(BearImpOptions.Possible.(mainNames{ii}))']; %#ok<AGROW>
+                end
+            end
+        end
+    end
+
+    methods
+        function obj = BearImpOptions(optionName,parent)
+            arguments
+                optionName = 'main'
+                parent = [];
+            end
+            %mustBeMember(optionName,BearImpOptions.AllCategories)
+
+            obj.optionName  = optionName;
+            obj.parent = parent;
+
+            if ~isempty(parent)
+                obj.optionLevel = parent.optionLevel + 1;
+            end
+
+            if isempty(obj.possibleSubOptions)
+                return
+            end
+            
+            for sub = obj.possibleSubOptions
+                if obj.hasSubSubOptions
+                    newProp = obj.addprop(sub{:});
+                    obj.(sub{:}) = BearImpOptions(sub{:},obj);
+                    newProp.SetMethod = @(obj,optionChoice) obj.set_optionChosen(sub{:},optionChoice);
+                    newProp.NonCopyable = false;
+                end
+            end
+
+            if ~obj.hasSubSubOptions
+                newProp = obj.addprop('defaultOption');
+                newProp.NonCopyable = false;
+                switch obj.optionLevel
+                    case 2
+                        obj.defaultOption = BearImpOptions.Default.(obj.optionName);
+                    case 3
+                        obj.defaultOption = BearImpOptions.Default.(obj.parent.optionName).(obj.optionName);
+                end
+            end
+        end
+    end
+
+    methods (Hidden)
+        function str = getOptionID(obj)
+            % Gerates an abbreviated string according to the options chosen
+            str = ['T' lower(obj.T.char(1)) '_B' lower(obj.B.char(1)) '_H' lower(obj.H.filmThicknessFormula.char(1)) ];
+            if obj.H.thermalCorrection    == 'on'    , str = [str 't']; end
+            if obj.H.starvationCorrection == 'on'    , str = [str 's']; end
+            if obj.H.roughnessCorrection  == 'on'    , str = [str 'r']; end
+            str = [str '_C' lower(obj.C.outsideArea.char(1:2))];
+            if obj.C.k_vh_factor          == 'on'    , str = [str 'k']; end
+            if obj.C.pressureDistribution == 'on'    , str = [str 'p']; end
+            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)
+            arguments
+                obj
+                n {mustBeInteger} = []
+            end
+            if strcmp(obj.optionChosen,'default') && ~obj.hasSubSubOptions
+                charOut = obj.defaultOption;
+            else
+                charOut = obj.optionChosen;
+            end
+            if ~isempty(n), charOut = charOut(n); end
+        end
+
+        function stringOut = string(obj)
+            stringOut = string(char(obj));
+        end
+
+        function cellOut = cell(obj)
+            cellOut = {char(obj)};
+        end
+
+        function logicalOut = ne(obj,otherObj)
+            nargoutchk(0,1)
+            logicalOut = ~eq(obj,otherObj);
+        end
+
+        function [logicalOut,diffs] = eq(obj,otherObj)
+            nargoutchk(0,2)
+            logicalOut = true;
+            diffs = {};
+            if ischar(otherObj) || isstring(otherObj) || iscellstr(otherObj)
+                logicalOut = strcmp(char(obj),otherObj);
+                return
+            end
+            if ischar(obj) || isstring(obj) || iscellstr(obj)
+                logicalOut = strcmp(char(otherObj),obj);
+                return
+            end
+            try
+                if obj.hasSubSubOptions
+                    for sub = obj.possibleSubOptions
+                        if obj.(sub{:}).hasSubSubOptions
+                            for subsub = obj.(sub{:}).possibleSubOptions
+                                if ~strcmp(obj.(sub{:}).(subsub{:}).optionChosen,otherObj.(sub{:}).(subsub{:}).optionChosen)
+                                    logicalOut = false;
+                                    if obj.optionLevel == 1
+                                        diffs = [diffs {obj.(sub{:}).optionName}]; %#ok<AGROW>
+                                    else
+                                        diffs = [diffs {obj.optionName}]; %#ok<AGROW>
+                                    end
+                                end
+                            end
+                        else
+                            if ~strcmp(obj.(sub{:}).optionChosen,otherObj.(sub{:}).optionChosen)
+                                logicalOut = false;
+                                if obj.optionLevel == 1
+                                    diffs = [diffs {obj.(sub{:}).optionName}]; %#ok<AGROW>
+                                else
+                                    diffs = [diffs {obj.optionName}]; %#ok<AGROW>
+                                end
+                            end
+                        end
+                    end
+                else
+                    logicalOut = strcmp(obj.optionChosen,otherObj.optionChosen);
+                end
+            catch
+                warning('Unsimilar BearImpOptions objects were compared')
+                logicalOut = false;
+            end
+            diffs = unique(diffs);
+        end
+
+        function set_optionChosen(obj,subOptionName,optionChoice)
+            objBefore = obj.copy;
+            if ischar(optionChoice) || isstring(optionChoice)
+                assert(~obj.(subOptionName).hasSubSubOptions || ...
+                    isfield(BearImpOptions.subOptionsSetWith,subOptionName) || ...
+                    strcmp(optionChoice,'default'),...
+                    [obj.(subOptionName).optionName ' has suboptions. Set suboptions individually or set to ''default''.'])
+                if isfield(BearImpOptions.subOptionsSetWith,subOptionName)
+                    for subsub = BearImpOptions.subOptionsSetWith.(subOptionName)
+                        obj.(subOptionName).(subsub{:}) = optionChoice;
+                    end
+                    return
+                end
+                mustBeMember(optionChoice,[obj.(subOptionName).possibleSubOptions {'default'}])
+                if ~obj.(subOptionName).hasSubSubOptions && strcmp(optionChoice,obj.(subOptionName).defaultOption)
+                    optionChoice = 'default';
+                end
+                obj.(subOptionName).optionChosen = optionChoice;
+                if strcmp(optionChoice,'default')
+                    if obj.(subOptionName).hasSubSubOptions
+                        for sub = obj.(subOptionName).possibleSubOptions
+                            obj.(subOptionName).(sub{:}).optionChosen = 'default';
+                        end
+                    end
+                    obj.optionChosen = 'default';
+                else
+                    obj.optionChosen = 'custom suboptions';
+                    if obj.optionLevel == 2
+                        obj.parent.optionChosen = 'custom suboptions';
+                    end
+                end
+            else
+                assert(isa(optionChoice,'BearImpOptions'),'Wrong input datatype for setting an option. Use a string, char or BearImpOptions object.')
+                assert(strcmp(obj.(subOptionName).optionName,optionChoice.optionName),'The name of the BearImpOptions object differs from the one delivered.')
+
+                prop = findprop(obj,subOptionName);
+                setMethod = prop.SetMethod; % Wired workaround so that the set function does not call itself recursively
+                prop.SetMethod = [];
+                obj.(subOptionName) = optionChoice;
+                prop.SetMethod = setMethod;
+            end
+            [isEqual,optionsChanged] = eq(obj,objBefore);
+            if ~isEqual
+                mainObj = obj.getMainObj;
+                mainObj.optionsLastChanged = optionsChanged;
+                notify(mainObj,'optionChanged')
+            end
+        end
+    end
+
+    methods (Access = protected)
+        function propgrp = getPropertyGroups(obj)
+            if ~isscalar(obj)
+                propgrp = getPropertyGroups@matlab.mixin.CustomDisplay(obj);
+            else
+                if obj.hasSubSubOptions
+                    for sub = obj.possibleSubOptions
+                        if obj.(sub{:}).hasSubSubOptions
+                            for subsub = obj.(sub{:}).possibleSubOptions
+                                propList.([sub{:} '_' subsub{:}]) = obj.(sub{:}).(subsub{:}).displayText;
+                            end
+                        else
+                            propList.(sub{:}) = obj.(sub{:}).displayText;
+                        end
+                    end
+                else
+                    propList = struct(obj.optionName,obj.displayText);
+                end
+                propgrp = matlab.mixin.util.PropertyGroup(propList);
+            end
+        end
+        function outText = displayText(obj)
+            outText = obj.optionChosen;
+            if strcmp(outText,'default') && ~obj.hasSubSubOptions
+                outText = [obj.defaultOption ' (default)'];
+            end
+        end
+        function cp = copyElement(obj)
+            cp = copyElement@matlab.mixin.Copyable(obj);
+            if obj.hasSubSubOptions
+                for sub = obj.possibleSubOptions
+                    prop = findprop(cp,sub{:});
+                    prop.SetMethod = []; % Wired workaround so that the set function is not called as it calls copy resulting in a recursion
+                    cp.(sub{:}) = obj.(sub{:}).copy;
+                    prop.SetMethod = @(cp,optionChoice) cp.set_optionChosen(sub{:},optionChoice);
+                    cp.(sub{:}).parent = cp;
+                end
+            end
+        end
+    end
+
+    methods (Hidden)
+        function addSetMethods(obj)
+            % Needs to be called after loading a BearImpOptions object as
+            % the set methods of dynamic properties are not loaded
+            % alongside
+            if obj.hasSubSubOptions
+                for sub = obj.possibleSubOptions
+                    prop = findprop(obj,sub{:});
+                    prop.SetMethod = @(obj,optionChoice) obj.set_optionChosen(sub{:},optionChoice);
+                    prop.NonCopyable = false;
+                    obj.(sub{:}).addSetMethods
+                end
+            end
+        end
+    end
+    
+    methods (Access = protected)
+        function mainObj = getMainObj(obj)
+            switch obj.optionLevel
+                case 1
+                    mainObj = obj;
+                case 2
+                    mainObj = obj.parent;
+                case 3
+                    mainObj = obj.parent.parent;
+            end
+        end
+    end
+
+    methods % get functions
+        function possibleSubOptions = get.possibleSubOptions (obj)
+            if obj.optionLevel == 1
+                possibleSubOptions = fieldnames(BearImpOptions.Possible)';
+            elseif obj.optionLevel == 2
+                if isstruct(BearImpOptions.Possible.(char(obj.optionName)))
+                    possibleSubOptions = fieldnames(BearImpOptions.Possible.(char(obj.optionName)))';
+                elseif iscell(BearImpOptions.Possible.(char(obj.optionName)))
+                    possibleSubOptions = BearImpOptions.Possible.(char(obj.optionName));
+                end
+            elseif obj.optionLevel == 3
+                if obj.hasSubSubOptions
+                    possibleSubOptions = fieldnames(BearImpOptions.Possible.(obj.parent.optionName))';
+                else
+                    possibleSubOptions = BearImpOptions.Possible.(obj.parent.optionName).(obj.optionName);
+                end
+            elseif isstruct(BearImpOptions.Possible.(obj.parent.optionName)) && ismember(obj.optionName,fieldnames(BearImpOptions.Possible.(obj.parent.optionName)))
+                
+            else
+                possibleSubOptions = {};
+            end
+        end
+
+        function hasSubSubOptions = get.hasSubSubOptions(obj)
+            switch obj.optionLevel
+                case 1
+                    hasSubSubOptions = true;
+                case 2
+                    hasSubSubOptions = isstruct(BearImpOptions.Possible.(obj.optionName));
+                case 3
+                    if ~isstruct(BearImpOptions.Possible.(obj.parent.optionName))
+                        hasSubSubOptions = false;
+                    else
+                        hasSubSubOptions = isstruct(BearImpOptions.Possible.(obj.parent.optionName).(obj.optionName));
+                    end
+            end
+        end
+    end
+end
\ No newline at end of file
diff --git a/InputData.xlsx b/InputData.xlsx
index b05ad322385e1b14a930a628e0de318cd9048e8f..7ec4e845e300f7f45c330e5dff0bbe7031791434 100644
Binary files a/InputData.xlsx and b/InputData.xlsx differ
diff --git a/possibleMethods.m b/possibleMethods.m
deleted file mode 100644
index f92daad6ff0cda72758ec925432bc398c5783912..0000000000000000000000000000000000000000
--- a/possibleMethods.m
+++ /dev/null
@@ -1,112 +0,0 @@
-classdef possibleMethods
-% Gibt Auskunft über mögliche Berechnungsmethoden der einzelnen
-% Berechnungsabschnitte (T, R, G, B, H, C) und enthält eine Funktion zum
-% überprüfen des Methodenstructs
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Julius van der Kuip
-
-    methods (Static) % Diese Klasse enthält ausschließlich statische Methoden, sodass kein Objekt initialisiert werden muss
-        function s = T
-            s = {'linear','Vogel'};
-        end
-        function s = R
-            s = {};
-        end
-        function s = G
-            s = {};
-        end
-        function s = B
-            s = {'static','dynamic'};
-        end
-        function s = H
-            s = {'Hamrock/Dowson','Hamrock/Dowson_withoutThermCorr','Moes'};
-        end
-        function s = C
-            s.unloadedRE  = {           'neglect','stateOfTheArt',                                 'Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
-            s.outsideArea = {'k-factor','neglect','stateOfTheArt','Schneider_k_c','Schneider_k_vh','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
-        end
-        
-        function s = Default
-        % Gibt ein Rechenmethoden-Struct mit den Default-Methoden aus
-            s.T = 'Vogel';
-            s.B = 'dynamic';
-            s.H = 'Hamrock/Dowson';
-            s.C.unloadedRE  = 'semianalytisch3D';
-            s.C.outsideArea = 'semianalytisch3D';
-        end
-        
-        function method = addDefault(method)
-        % Ergänzt ein unvollständiges Rechenmethoden-Struct mit den
-        % Default-Werten
-            arguments
-                method (1,1) struct
-            end
-            defaultMethod = possibleMethods.Default;
-            for ii = fields(defaultMethod)'
-                if ~isfield(method,ii{1}) % Falls Feld ii nicht besetzt,
-                    method.(ii{1}) = defaultMethod.(ii{1}); % überschreibe mit default
-                elseif isstruct(defaultMethod.(ii{1}))  % falls zweite Ebene existiert
-                    for jj = fields(defaultMethod.(ii{1}))'
-                        if ~isfield(method.(ii{1}),jj{1}) % Falls Feld ii in jj nicht besetzt
-                            method.(ii{1}).(jj{1}) = defaultMethod.(ii{1}).(jj{1});
-                        end
-                    end   
-                end
-            end
-        end
-        
-        function success = check(method,throwError)
-        % prüft, ob das Struct method ausschließlich valide Feld- und Methodennamen enthält
-        % optional: throwError = true gibt error-Meldungen aus
-            arguments
-                method (1,1) struct
-                throwError (1,1) logical = false
-            end
-            success = true;
-            for ii = fields(method)'
-                if ~any(strcmp(ii,methods(possibleMethods)')) % Prüft Feldname in erster Ebene
-                    success = false;
-                    if throwError
-                        error(['%s ist kein gültiger Feldname im Methoden-Struct\n',...
-                               'Möglich an dieser Stelle: %s'],...
-                               ii{1}, strjoin(methods(possibleMethods)')   )
-                    end
-                    break
-                end
-                if isstruct(method.(ii{1}))           % falls zweite Ebene existiert        
-                    for jj = fields(method.(ii{1}))'
-                        if ~any(strcmp(jj,fields(possibleMethods.(ii{1}))'))                   % Prüft Feldname in zweiter Ebene
-                            success = false;
-                            if throwError
-                                error(['%s is kein gültiger Feldname für den Berechnungsabschnitt %s\n',...
-                                       'Möglich an dieser Stelle: %s'],...
-                                       jj{1},ii{1},strjoin(fields(possibleMethods.(ii{1}))')   )
-                            end
-                            break
-                        end
-                        if ~any(strcmp(method.(ii{1}).(jj{1}),possibleMethods.(ii{1}).(jj{1}))) % Prüft Methodenname in zweiter Ebene
-                            success = false;
-                            if throwError
-                                error(['%s ist kein gültiger Methodenname zur Berechnung von %s im Berechnungsabschnitt %s\n',...
-                                       'Möglich an dieser Stelle: %s'],...
-                                       method.(ii{1}).(jj{1}), jj{1}, ii{1}, strjoin(possibleMethods.(ii{1}).(jj{1}))   )
-                            end
-                            break
-                        end
-                    end
-                else
-                    if ~any(strcmp(method.(ii{1}),possibleMethods.(ii{1}))) % Prüft Methodenname in erster Ebene
-                        success = false;
-                        if throwError
-                            error(['%s ist kein gültiger Methodenname für den Berechnungsabschnitt %s\n',...
-                                   'Möglich an dieser Stelle: %s'],...
-                                   method.(ii{1}), ii{1}, strjoin(possibleMethods.(ii{1}))   )
-                        end
-                        break
-                    end
-                end
-            end
-        end
-    end
-end
\ No newline at end of file