diff --git a/@BearImp/BearImp.m b/@BearImp/BearImp.m
index 29adea5230071a2e16cef084d78ff95b21667186..6d406fd472e4a74b0e9c0f315055eae99b8f7177 100644
--- a/@BearImp/BearImp.m
+++ b/@BearImp/BearImp.m
@@ -1,45 +1,54 @@
 classdef BearImp < handle & matlab.mixin.Copyable
 % pmd Berechnungstool Lagerimpedanz
 % Fachgebiet für Produktentwicklung und Maschinenelemente, TU Darmstadt
-% Version: 2.0.0
-% Stand: 25.06.2021
-% Autoren: Steffen Puchtler, Tobias Schirra
+% Autoren: Steffen Puchtler, Tobias Schirra, Julius van der Kuip
 
     properties (Dependent, Access = public)
     % Eingangsparameter Zugriff
-        F_r     (1,1) double {mustBeNonnegative}       % Radialkraft  in N
-        F_a     (1,1) double {mustBeNonnegative}       % Axialkraft   in N
-        n       (1,1) double {mustBePositive}          % Drehzahl     in 1/s !!!
-        T_Oil   (1,1) double {mustBeReal}              % Öltemperatur in °C
-        psi     (1,:) double {mustBeReal}              % Winkel, für die die Kapazität berechnet wird, in rad
-        L (1,1) struct  %     Lagerparameter (gegeben)
-        S (1,1) struct  %     Schmierstoffparameter (gegeben)
+        F_r      (1,1) double {mustBeNonnegative}  % Radialkraft  in N
+        F_a      (1,1) double {mustBeNonnegative}  % Axialkraft   in N
+        omega    (1,1) double {mustBePositive}     % Winkelgeschwindigkeit  in 1/s
+        T_Oil    (1,1) double {mustBeReal}         % Öltemperatur in °C
+        psi      (1,:) double {mustBeReal}         % Winkel, die vom Benutzer für die Kapazitätsberechnung verwendet wird, in rad
+        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
+        
         
     % Structs mit Berechnungsparamtern
-        T (1,1) struct  % 2.1 Tribologie (berechnete Schmierstoffparameter)         f(  S          ,T_Oil)
-        R (1,1) struct  % 2.2 Radiallagerluft                                       f(L                  )
-        G (1,1) struct  % 2.3 Geometrie (berechnet in Abhängigkeit des Lagerspiels) f(L,    R            )
-        B (1,1) struct  % 2.4 Belastungsverteilung                                  f(L,    R,G          )
-        H (1,1) struct  % 2.5 Schmierfilmdicke                                      f(L,S,T,  G,B        )
-        Z (1,1) struct  % 2.6 Impedanz                                              f(L,S,T,  G,B,H      )
+        T (1,1) struct  % 2.1 Tribologie (berechnete Schmierstoffparameter)         f(  S             T_Oil                                     )
+        R (1,1) struct  % 2.2 Radiallagerluft                                       f(L                     F_r                                 )
+        G (1,1) struct  % 2.3 Geometrie (berechnet in Abhängigkeit des Lagerspiels) f(L,    R         T_Oil                                     )
+        B (1,1) struct  % 2.4 Belastungsverteilung                                  f(L,      G             F_r,psi,AddOn       psi_calc psi_all)
+        H (1,1) struct  % 2.5 Schmierfilmdicke                                      f(L,S,T,  G,B     T_Oil,              omega                 )
+        C (1,1) struct  % 2.6 Kapazität und Widerstand                              f(L,S,T,  G,B,H             psi,AddOn       psi_calc psi_all)
+        Z (1,1) struct  % 2.7 Impedanz                                              f(L             C                                           )
     end
     
     
     properties (Access = private)
     % Eingangsparameter Speicher
-        privateInputParameters = struct('F_r',nan,'F_a',nan,'n',nan,'T_Oil',nan,'psi',nan,'f',nan,'L',struct,'S',struct)
+        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,'Z',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')
     end
     
     properties (Access = public)
         % gibt an, ob die Werte vorhanden und aktuell sind
-        up2date = struct('L',false,'S',false,'T',false,'R',false,'G',false,'B',false,'H',false,'Z',false)
+        up2date = struct('L',false,'S',false,'T',false,'R',false,'G',false,'B',false,'H',false,'C',false,'Z',false)
         materials
         inputDataTableFile = 'InputData.xlsx';
     end
     
+    properties (Dependent, SetAccess = protected, GetAccess = public)
+        psi_calc (1,:) double {mustBeReal}         % Winkel, die benötigt werden, um die Kapazität berechnen zu können, (unter ausnutzung der Symmetrie und ohne Dopplungen) in rad
+        psi_all  (1,:) double {mustBeReal}         % alle Winkel, die unter Ausnutzung der Symmetrie mit psi_calc berechnet wurden
+        ind_psi_all  (:,:) double {mustBeNonnegative, mustBeInteger}  % Indizes überführen psi_calc zu psi_all (damit kann zB. h_min(ind_psi_all) zu einem maximalen Vektor formiert werden) 
+    end
+    
     properties (Dependent, SetAccess = immutable, GetAccess = public)
         allUp2date (1,1) logical % true, wenn alle berechneten Größen auf dem aktuellen Stand sind
         ready2calc (1,1) logical % true, wenn alle Prameter gegeben zur Berechnung
@@ -58,13 +67,15 @@ classdef BearImp < handle & matlab.mixin.Copyable
                 switch setup
                     case 'default'
                         obj.setBearing('6205-C3');
-                        obj.setLube('FVA Referenz-Öl II');
+                        obj.setLube('FVA Referenz-Öl III');
                         obj.F_r = 1000;
                         obj.F_a = 0;
-                        obj.n = 1000/60;
+                        obj.omega = 3000/60*2*pi;
                         obj.T_Oil = 50;
+                        obj.resolutionPerRotation = 360;
                 end
             end
+            
         end
     end
     
@@ -75,6 +86,7 @@ classdef BearImp < handle & matlab.mixin.Copyable
         calcGeo  (obj)
         calcLoad (obj)
         calcFilm (obj)
+        calcCap  (obj)
         calcImp  (obj)
     end
     methods (Access = public)
@@ -100,26 +112,83 @@ classdef BearImp < handle & matlab.mixin.Copyable
             obj.up2date.G = false;
             obj.privateInputParameters.F_a = F_a;
         end
-        function set.n(obj,n)
+        function set.omega(obj,omega)
             obj.up2date.H = false;
-            obj.privateInputParameters.n = n;
+            obj.privateInputParameters.omega = omega;
         end
         function set.T_Oil(obj,T_Oil)
             obj.up2date.H = false;
             obj.up2date.T = false;
+            obj.up2date.G = false;
             obj.privateInputParameters.T_Oil = T_Oil;
         end
+        function set.resolutionPerRotation(obj,resolutionPerRotation)
+            obj.up2date.B = false;
+            obj.up2date.C = false;
+            obj.privateInputParameters.resolutionPerRotation = resolutionPerRotation;
+            obj.privateInputParameters.psi = (0:(2*pi)/(resolutionPerRotation):pi*2*(resolutionPerRotation-1)/resolutionPerRotation);
+        end
         function set.psi(obj,psi)
-            obj.up2date.R = false;
+            obj.up2date.C = false;
             obj.up2date.B = false;
             obj.privateInputParameters.psi = psi;
+            obj.privateInputParameters.resolutionPerRotation = nan;
+        end
+        function set.psi_calc(obj,psi_calc)
+            obj.up2date.C = false;
+            obj.up2date.B = false;
+            obj.privateInputParameters.psi_calc = psi_calc;
+            obj.privateInputParameters.resolutionPerRotation = nan;
+        end
+        function set.psi_all(obj,psi_all)
+            obj.up2date.C = false;
+            obj.up2date.B = false;
+            obj.privateInputParameters.psi_all = psi_all;
+            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)
+            obj.up2date.B = false;
+            obj.up2date.C = false;
+
+            pct = AddOn.Parallel_Computing_Toolbox;   
+            switch pct
+                case true
+                    pct  = 'on';
+                case false 
+                    pct = 'off';
+                otherwise
+            end
+
+            ot  = AddOn.Optimization_Toolbox;
+            switch ot
+                case true
+                    ot  = 'on';
+                case false 
+                    ot = 'off';
+                otherwise
+            end
+
+            check = {'auto', 'on', 'off'};
+            if ~ischar(ot)||~ischar(pct)
+                error("Wrong Input for AddOn. You have to decide between 'on', 'off' and 'auto'.")
+            elseif ~any(strcmp(pct,check)==1) || ~any(strcmp(ot,check)== 1)
+               error("Check if you spelled the AddOn Input the right way! ('on' 'off' 'auto')")
+            end
+           
+            AddOn.Optimization_Toolbox = ot; 
+            AddOn.Parallel_Computing_Toolbox = pct;
+            obj.privateAddOn = AddOn;
         end
         function set.L(obj,L)
             obj.up2date.L = true;
             obj.up2date.R = false;
             obj.up2date.G = false;
             obj.up2date.B = false;
-            obj.up2date.H = false;
             obj.up2date.Z = false;
             obj.privateInputParameters.L = L;
         end
@@ -127,13 +196,13 @@ classdef BearImp < handle & matlab.mixin.Copyable
             obj.up2date.S = true;
             obj.up2date.T = false;
             obj.up2date.H = false;
-            obj.up2date.Z = false;
+            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','Z'}
+            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;
@@ -163,39 +232,43 @@ classdef BearImp < handle & matlab.mixin.Copyable
             end
             
             % Abhängige Variablen auf nicht up2date setzen
-            if ~up2date.T,                    up2date.H = false;                    end
-            if ~up2date.R, up2date.B = false; up2date.G = false;                    end
-            if ~up2date.G, up2date.B = false; up2date.H = false; up2date.Z = false; end
-            if ~up2date.B,                    up2date.H = false; up2date.Z = false; end
-            if ~up2date.H,                    up2date.Z = false;                    end
+            if ~up2date.T,                    up2date.H = false; up2date.C = false; end
+            if ~up2date.R, up2date.G = false;                                       end
+            if ~up2date.G, up2date.B = false; up2date.H = false; up2date.C = false; end
+            if ~up2date.B,                    up2date.H = false; up2date.C = false; end
+            if ~up2date.H,                                       up2date.C = false; end
+            if ~up2date.C, up2date.Z = false;                                       end
             
             obj.up2date = up2date;
         end
         function set.T(obj,T)
             obj.up2date.H = false;
-            obj.up2date.Z = false;
+            obj.up2date.C = false;
             obj.privateResults.T = T;
         end
         function set.R(obj,R)
             obj.up2date.G = false;
-            obj.up2date.B = false;
             obj.privateResults.R = R;
         end
         function set.G(obj,G)
             obj.up2date.B = false;
             obj.up2date.H = false;
-            obj.up2date.Z = false;
+            obj.up2date.C = false;
             obj.privateResults.G = G;
         end
         function set.B(obj,B)
             obj.up2date.H = false;
-            obj.up2date.Z = false;
+            obj.up2date.C = false;
             obj.privateResults.B = B;
         end
         function set.H(obj,H)
-            obj.up2date.Z = false;
+            obj.up2date.C = false;
             obj.privateResults.H = H;
         end
+        function set.C(obj,C)
+            obj.up2date.Z = false;
+            obj.privateResults.C = C;
+        end
         function set.Z(obj,Z)
             obj.privateResults.Z = Z;
         end
@@ -209,14 +282,69 @@ classdef BearImp < handle & matlab.mixin.Copyable
         function F_a = get.F_a(obj)
             F_a = obj.privateInputParameters.F_a;
         end
-        function n = get.n(obj)
-            n = obj.privateInputParameters.n;
+        function omega = get.omega(obj)
+            omega = obj.privateInputParameters.omega;
         end
         function T_Oil = get.T_Oil(obj)
             T_Oil = obj.privateInputParameters.T_Oil;
         end
         function psi = get.psi(obj)
             psi = obj.privateInputParameters.psi;
+        end
+        function psi_calc = get.psi_calc(obj)
+            psi_calc = obj.privateInputParameters.psi_calc;
+        end
+        function psi_all = get.psi_all(obj)
+            psi_all = obj.privateInputParameters.psi_all;
+        end
+        function ind_psi_all = get.ind_psi_all(obj)
+            ind_psi_all = obj.privateInputParameters.ind_psi_all;
+        end
+        function resolutionPerRotation = get.resolutionPerRotation(obj)
+            resolutionPerRotation = obj.privateInputParameters.resolutionPerRotation;
+        end
+        function AddOn = get.AddOn(obj)
+            
+            pct = obj.privateAddOn.Parallel_Computing_Toolbox;          
+            ot  = obj.privateAddOn.Optimization_Toolbox;
+
+            switch pct
+                case 'on'
+                    AddOn.Parallel_Computing_Toolbox = true;
+                case 'off'
+                    AddOn.Parallel_Computing_Toolbox = false;
+                case 'auto'
+                    table_AddOns = matlab.addons.installedAddons;
+                    table_AddOns = table2array(table_AddOns(:,1));
+                    
+                    if any(strcmp(table_AddOns, "Parallel Computing Toolbox"),1) && matlab.addons.isAddonEnabled("Parallel Computing Toolbox")
+                        AddOn.Parallel_Computing_Toolbox = true;
+                    else
+                        AddOn.Parallel_Computing_Toolbox = false;
+                        disp("To speed up the calculation you can install the Parallel_Computing_Toolbox")
+                    end
+               
+            end
+            
+            switch ot
+                case 'on'
+                    AddOn.Optimization_Toolbox = true;
+                case 'off'
+                    AddOn.Optimization_Toolbox = false;
+                case 'auto'
+                    table_AddOns = matlab.addons.installedAddons;
+                    table_AddOns = table2array(table_AddOns(:,1));
+                    
+                    if any(strcmp(table_AddOns, "Optimization Toolbox"),1) && matlab.addons.isAddonEnabled("Optimization Toolbox")
+                        AddOn.Optimization_Toolbox = true;
+                    else
+                        AddOn.Optimization_Toolbox = false;
+                        disp("To get more precise and faster calculations you can install the Optimization_Toolbox")
+                    end
+               
+            end
+            
+   
         end
         function L = get.L(obj)
             L = obj.privateInputParameters.L;
@@ -242,6 +370,9 @@ classdef BearImp < handle & matlab.mixin.Copyable
         function H = get.H(obj)
             H = obj.privateResults.H;
         end
+        function C = get.C(obj)
+            C = obj.privateResults.C;
+        end
         function Z = get.Z(obj)
             Z = obj.privateResults.Z;
         end
@@ -249,7 +380,7 @@ classdef BearImp < handle & matlab.mixin.Copyable
             allUp2date = all(struct2array(obj.up2date));
         end
         function ready2calc = get.ready2calc(obj)
-            ready2calc = all(~isnan([obj.F_r obj.F_a obj.n obj.T_Oil])) &&...
+            ready2calc = all(~isnan([obj.F_r obj.F_a obj.omega obj.T_Oil obj.psi])) &&...
                          obj.up2date.L       &&...
                          obj.up2date.S;
         end
diff --git a/@BearImp/calcCap.m b/@BearImp/calcCap.m
new file mode 100644
index 0000000000000000000000000000000000000000..dca5cd4885d4c768f9d4d5e65208747e534a19b0
--- /dev/null
+++ b/@BearImp/calcCap.m
@@ -0,0 +1,462 @@
+function calcCap(obj)
+%2.6 Kapazität und Widerstand C = C(L,S,T,G,B,H)
+%Berechnet die Kapazität der Hertz'schen Flächen, sowie der Randbereiche
+%und unbelasteter Wälzkörper je nach gewählter Methode. Zusätzlich wird der
+%Widerstand der Hertz'schen Flächen bestimmt und anschließend die Impedanz
+%berechnet.
+
+%    pmd Berechnungstool Lagerimpedanz
+%    Autor: Steffen Puchtler, Tobias Schirra, Leander, Julius van der Kuip
+
+%% Parameter prüfen
+assert(obj.up2date.L,'Lager nicht gesetzt')
+assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
+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;
+C = struct;
+method = possibleMethods.addDefault(obj.method).C;
+if ~isfield(C,'k_C')
+    C.k_C = 1;
+end
+if ~isfield(C,'k_R')
+    C.k_R = 1;
+end
+
+%% Berechnung
+
+assert(~L.allBallsSameMaterial || L.RE_A.electric_behaviour == "conductor",'Die elektrischen Eigenschaften des Lagers können nur von RE berechnet werden, die ein leitendes Material besitzen. Bitte Material der RE prüfen.')
+   
+if ~isnan(obj.psi_calc) % Falls psi_calc gegeben ist
+   psi = obj.psi_calc;  % psi_calc in calcCap als psi verwenden
+end
+
+B.numOfCagePositions = length(psi);
+   
+if L.allBallsSameMaterial %|| strcmp(L.RE_order,'AAAAAAAAA')
+   ballMaterialInd=1;
+else
+   ballMaterialInd=L.RE_order_num(L.IndBall_conductive);
+end
+
+C.C_Hertz     = zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.phi_1i      = nan  (1,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.phi_1a      = nan  (1,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.C_out       = zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.c_i         = nan  (1,B.numOfCagePositions);
+C.c_a         = nan  (1,B.numOfCagePositions);
+C.g_1i        = zeros(  B.numOfCagePositions,L.numberOfConductiveBalls);
+C.g_1a        = zeros(  B.numOfCagePositions,L.numberOfConductiveBalls);
+R_el_single = inf  (2,B.numOfCagePositions,L.numberOfConductiveBalls);
+
+for posBall_conductive=1:L.numberOfConductiveBalls
+
+    temp_s=      B.s(:,:,posBall_conductive);
+    temp_p=      H.p(:,:,posBall_conductive);
+    temp_A=      H.A(:,:,posBall_conductive);
+    temp_a=      H.a(:,:,posBall_conductive);
+    temp_b=      H.b(:,:,posBall_conductive);
+    temp_h_0th = H.h_0th(:,:,posBall_conductive);
+    tempR_li   = G.R_li(ballMaterialInd(posBall_conductive));
+    tempR_la   = G.R_la(ballMaterialInd(posBall_conductive));
+    tempR_RE   = G.R_RE(ballMaterialInd(posBall_conductive));
+        
+    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_0th(:,contactInd);
+ 
+    R_el_single(:,contactInd,posBall_conductive) = 1./C.k_R .* S.rho_el .* temp_h_0th(:,contactInd) ./ temp_A(:,contactInd);
+    C.DeltaR_i(contactInd)    = tempR_li - tempR_RE - temp_h_0th(1,contactInd);
+    C.DeltaR_a(contactInd)    = tempR_la - tempR_RE - temp_h_0th(2,contactInd);
+
+    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,:)',posBall_conductive);
+        C.DeltaR_a(B.noContactInd(posBall_conductive,:)') = tempR_la - tempR_RE - temp_s(2,B.noContactInd(posBall_conductive,:)',posBall_conductive);
+    end
+    
+    C.phi_1i(:,:,posBall_conductive) = temp_a(1,:)./G.D_RE(ballMaterialInd(posBall_conductive));
+    C.phi_1a(:,:,posBall_conductive) = temp_a(2,:)./G.D_RE(ballMaterialInd(posBall_conductive));
+    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'})])
+        calcZg
+    end
+    
+    if size(B.noContactInd,2) % Wenn Kraftfreie Wälzkörper vorhanden
+        switch method.unloadedRE
+            case 'neglect'
+            case 'Leander_Parallel'
+                leander_parallel(B.noContactInd)
+            case 'Leander_Radial'
+                leander_radial(find(contactInd==0))
+            case 'LeanderSteffen'
+                leander_steffen(B.noContactInd)
+            case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
+                tobias_steffen(B.noContactInd)
+            case 'semianalytisch3D'
+
+                    vecStruct=splitVec(contactInd,false);
+            
+                for recentVec=length(vecStruct)
+                    runSemianalytisch3D(vecStruct{recentVec})
+                end
+                
+        end
+    end
+
+    switch method.outsideArea
+        case 'neglect'
+        case 'k-factor'
+            C.C_Hertz(:,contactInd,posBall_conductive) = C.k_C .* obj.epsilon_0 .* C.epsilon_primeOil(:,contactInd) .* temp_A(:,contactInd) ./ H.h_minth(:,contactInd);
+        case 'Leander_Parallel'
+            leander_parallel(B.contactInd)
+        case 'Leander_Radial'
+            leander_radial(find(contactInd==1))
+        case 'LeanderSteffen'
+            leander_steffen(B.contactInd)
+        case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
+            tobias_steffen(B.contactInd)
+        case 'semianalytisch3D'
+
+                   vecStruct=splitVec(contactInd,true);
+                
+            for recentVec=1:length(vecStruct)
+                runSemianalytisch3D(vecStruct{recentVec})
+            end
+    end
+
+    C_el_single(:,:,posBall_conductive) = C.C_Hertz(:,:,posBall_conductive) + C.C_out(:,:,posBall_conductive);
+    
+end
+
+
+if L.allBallsSameMaterial
+
+    C_ind = check_convert_psi_calc2matrx(psi,L.Z,'useSymmetry');
+
+    C.C_el_single = convert_vek2matr(C_el_single,C_ind);
+    C.R_el_single = convert_vek2matr(R_el_single,C_ind);
+    
+    C_el_single = C.C_el_single;
+    R_el_single = C.R_el_single;
+else
+    check_convert_psi_calc2matrx(psi,L.Z,'noSymmetry');
+    C.C_el_single = C_el_single;
+    C.R_el_single = R_el_single;
+end
+
+
+
+C.C_el_aprx = 1 ./    (1./ sum(C_el_single(1,:,:),3) + 1./ sum(    C_el_single(2,:,:) ,3));  
+C.R_el_aprx = 1 ./ sum(1./     R_el_single(1,:,:),3) + 1./ sum(1./ R_el_single(2,:,:) ,3) ;
+
+%% Attribute ändern
+obj.C = C;
+obj.up2date.C = true;
+
+%% Helper-Functions
+function calcZg
+              
+        C.g_1i(contactInd,posBall_conductive) = asin( sin( C.phi_1i(1,contactInd,posBall_conductive))  .*  (-C.DeltaR_i(:,contactInd)./tempR_li .* cos(C.phi_1i(1,contactInd,posBall_conductive)) + ...
+                                                        sqrt(1-C.DeltaR_i(:,contactInd).^2./tempR_li.^2 .* sin(C.phi_1i(1,contactInd,posBall_conductive)).^2)));
+       
+        C.g_1a(contactInd,posBall_conductive) = asin( sin( C.phi_1a(1,contactInd,posBall_conductive))  .*  (-C.DeltaR_a(:,contactInd)./tempR_la .* cos(C.phi_1a(1,contactInd,posBall_conductive)) + ... %hier ein - statt +?
+                                                        sqrt(1-C.DeltaR_a(:,contactInd).^2./tempR_la.^2 .* sin(C.phi_1a(1,contactInd,posBall_conductive)).^2)));
+        
+        C.g_2i = acos(L.B_i/(2*tempR_li));
+        C.g_2a = acos(L.B_a/(2*tempR_la));
+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);
+    
+    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.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? ###
+    
+    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.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;
+    for ii = indices
+        x1=double(C.g_1i(ii));
+        x2=double(C.g_2i);
+
+        C.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) C.dA_i(p)./C.h_i(p,t,ii), x1            , x2        , 0, pi/9);
+        C.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) C.dA_a(p)./C.h_a(p,t,ii), C.g_1a(ii)+pi, C.g_2a+pi, 0, pi/9);
+%             Z.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_i(p)./Z.h_i(p,t,ii), Z.g_1i(ii)   , Z.g_2i   , 0, pi/9);
+%             Z.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_a(p)./Z.h_a(p,t,ii), Z.g_1a(ii)+pi, Z.g_2a+pi, 0, pi/9);
+    end
+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')
+        C.dA_i = @(~,beta,~) tempR_RE^2 * cos(beta);
+        C.dA_a = C.dA_i;
+    else
+        C.dA_i = @(alpha,beta,ii) tempR_li./cos(beta) .* (1 - C.DeltaR_i(ii)*cos(alpha)./sqrt(tempR_li^2-C.DeltaR_i(ii)^2*sin(alpha).^2)) .* (tempR_RE + C.h_i(alpha,beta,ii));
+        C.dA_a = @(alpha,beta,ii) tempR_la./cos(beta) .* (1 - C.DeltaR_a(ii)*cos(alpha)./sqrt(tempR_la^2-C.DeltaR_a(ii)^2*sin(alpha).^2)) .* (tempR_RE + C.h_a(alpha,beta,ii));
+    end
+    for ii = indices
+        C.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) C.dA_i(alpha,beta,ii)./C.h_i(alpha,beta,ii), C.phi_1i(ii), C.phi_2i, 0, pi/3);
+        C.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) C.dA_a(alpha,beta,ii)./C.h_a(alpha,beta,ii), C.phi_1a(ii), C.phi_2a, 0, pi/3);
+    end
+end
+
+function runSemianalytisch3D(indices)
+    
+    if isfield(C,'C_out')
+        C_semi_i = C.C_out(1,:);
+        C_semi_a = C.C_out(2,:);
+    else
+        C_semi_i = nan( numel(temp_s)/2, 1);
+        C_semi_a = nan( numel(temp_s)/2, 1);
+    end
+    
+    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);
+
+        parfor runSemi = indices
+            C_semi_i(runSemi) = feval(tmp_semi_i, runSemi); %da parfor loop nicht auf nested functions (semianalytisch3D) direkt zugreifen kann
+            C_semi_a(runSemi )= feval(tmp_semi_a, runSemi);
+        end
+
+    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);
+        end
+    end
+
+    C.C_out(1,:,posBall_conductive)=C_semi_i;
+    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('Was ist hier los')
+    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
+    %mehrere Vektoren aufgeteilt, die diese Bedingungen erfüllen.
+    
+    findInds = find(vec==status);                                  
+    startInds = unique([findInds(1) findInds(diff([0 findInds])>1)]);      % Segment Start Indices
+    
+    if length(startInds)==1
+        seperateVecs{1}=findInds;
+    else
+        findNonZero = [find(diff([0 findInds])>1)-1 length(findInds)];
+        findNonZero(findNonZero == 0) = [];
+        endInds = findInds(findNonZero);   % Segment End Indices
+        
+        for runSection = 1:length(startInds)
+            seperateVecs{runSection} = (startInds(runSection):endInds(runSection));
+        end
+    end
+end
+
+function ind_psi = check_convert_psi_calc2matrx(psi,Z,symCase)
+    % Diese Funktion überprüft, ob der Input-Vektor psi alle Winkel
+    % enthält, die benötigt werden, um die Gesamtkapazität des Wälzlagers
+    % zu berechnen.
+    % Außerdem wird eine Matrix ind_psi erstellt, die die benötigten
+    % Indizes ausgibt, um dieberechneten Einzelkapazitäten in die richtige
+    % Reihenfolge und Ordnung zu bringen.
+    
+    switch symCase
+        case 'useSymmetry'
+            symmetry = true;
+        case 'noSymmetry'
+            symmetry = false;
+        otherwise
+            error('Falscher Eingabewert für symCase! Mögliche Eingabewerte: ''useSymmetry'' oder ''noSymmetry''')
+    end
+    ind_psi = [];
+
+    if symmetry
+        psi_calc = [psi, 2*pi - psi];
+    else
+        psi_calc = psi;
+    end
+
+    for run_C = 1 : length(obj.psi)
+
+       psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj.psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
+
+        psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) = psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
+
+        psiNeedForIndCapa = round(psiNeedForIndCapa,8);
+
+        [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-8);
+
+        %Check, if all nessecary C_el_single are calculated
+        assert(~any(ind_single_psi==0),'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!') 
+
+        ind_psi(:,run_C) = ind_single_psi';
+    end
+    ind_psi(ind_psi> length(psi)) = ind_psi(ind_psi> length(psi)) - length(psi);
+end
+function matrx_out = convert_vek2matr(A,B)
+    matrx_out = nan(size(A,1),size(B,2),size(B,1));
+    for run_1A = 1:size(A,1)
+        for run_1B = 1:size(B,1)
+            matrx_out(run_1A,:,run_1B) = A(run_1A,B(run_1B,:));
+        end
+    end
+end
+end
\ No newline at end of file
diff --git a/@BearImp/calcClear.m b/@BearImp/calcClear.m
index 00c156840cd6d7e152e7edaf59b200c937bd1eeb..f0e13550af4abdd02d6b5573c15e2a928563ac08 100644
--- a/@BearImp/calcClear.m
+++ b/@BearImp/calcClear.m
@@ -5,22 +5,16 @@ function calcClear(obj)
 %Setzbeträgen, etc.
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 06.04.2021
-%    Autor: Steffen Puchtler
+%    Autor: Steffen Puchtler, Julius van der Kuip
 
 %% Parameter prüfen
-if obj.up2date.R
-    warning('Lagerluft bereits berechnet')
-    return
-end
 assert(obj.up2date.L,'Lager nicht gesetzt')
 L = obj.L; F_r = obj.F_r;
 R = struct;
 
 %% Berechnung
 R.ue_d = 1/3.*L.ei_d + 2/3.*L.es_d - 2/3.*L.EI_d - 1/3.*L.ES_d;
-R.ue_D = 1/3.*L.ei_D + 2/3.*L.es_d - 2/3.*L.EI_D - 1/3.*L.EI_D;
+R.ue_D = 1/3.*L.ei_D + 2/3.*L.es_D - 2/3.*L.EI_D - 1/3.*L.ES_D;
 R.Deltad_ue = R.ue_d.*L.d_i./L.d;
 R.DeltaD_ue = R.ue_D.*L.D  ./L.D_G;
 if L.d <= 0.1
@@ -59,8 +53,6 @@ else
 end
 R.e_ir = 2000 .* L.d./L.ring.E .* (L.d./L.d_f)./(1-(L.d./L.d_f).^2) .* R.sigma_rir .* 1e-3;
 R.e_ar = 2000 .* L.D./L.ring.E .* (L.D_f./L.D)./(1-(L.D_f./L.D).^2) .* R.sigma_rar .* 1e-3;
-R.e_t = L.ring.alpha_t .* L.d_m .* 5; % aktuell mit fix Detla_t = 5 implementiert
-R.P_d = L.c - R.e_ir + R.e_ar - R.e_t + 0; % aktuelle mit fix e_L = 0 implementiert
 
 %% Attribute ändern
 obj.R = R;
diff --git a/@BearImp/calcFilm.m b/@BearImp/calcFilm.m
index f8f4906249da0de98d0b35525543cace95ed8174..fee783c09e2ff72917bf3558a5eb2aa1bbd9e8e7 100644
--- a/@BearImp/calcFilm.m
+++ b/@BearImp/calcFilm.m
@@ -1,54 +1,66 @@
 function calcFilm(obj)
-%2.5 Schmierfilmdicke H = H(L,S,T,G,B,T_Oil,n)
+%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)
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 07.04.2021
-%    Autor: Steffen Puchtler, Tobias Schirra
+%    Autor: Steffen Puchtler, Tobias Schirra, Julius van der Kuip
 
 %% Parameter prüfen
-if obj.up2date.H
-    warning('Belastungsverteilung betreits berechnet')
-    return
-end
 assert(obj.up2date.L,'Lager nicht gesetzt')
 assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
 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')
-L = obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; T_Oil = obj.T_Oil; n = obj.n;
+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;
 
 %% Berechnung
-H.b(1,:) = (6*G.fraktE_i.*B.Q.*G.R_si  ./  (pi*G.k_i.*G.E_red)).^(1/3);
-H.b(2,:) = (6*G.fraktE_a.*B.Q.*G.R_sa  ./  (pi*G.k_a.*G.E_red)).^(1/3);
-H.a = bsxfun(@times,[G.k_i;G.k_a],H.b);
+  
+Q_temp(1,:,:)=B.Q_contInd';
+Q_temp(2,:,:)=B.Q_contInd';
+
+if L.allBallsSameMaterial
+    ballMaterialInd=1;
+else
+    ballMaterialInd=L.RE_order_num(L.IndBall_conductive);
+end    
+
+H.b(1,:,:) = (6*G.fraktE_i.*B.Q_contInd'.*G.R_si(ballMaterialInd)  ./  (pi*G.k_i.*G.E_red(ballMaterialInd))).^(1/3);
+H.b(2,:,:) = (6*G.fraktE_a.*B.Q_contInd'.*G.R_sa(ballMaterialInd)  ./  (pi*G.k_a.*G.E_red(ballMaterialInd))).^(1/3);
+H.a = [G.k_i;G.k_a].*H.b;
 H.A = pi*H.a.*H.b;
-H.p = bsxfun(@rdivide,B.Q,H.A);
+H.p =Q_temp./H.A;
 H.alpha_p = 1  ./  (S.a_1 + S.a_2*T_Oil + (S.b_1 + S.b_2*T_Oil).*H.p);
-K.u_ir = (L.d_m.^2 - L.D_wk.^2 .* cosd(G.alpha).^2) ./ (4*L.d_m) .* 2 .* pi .* n;
-K.u_ar = K.u_ir;
-K.u = [K.u_ir;K.u_ar];
-U_tmp = T.eta_0 .* K.u ./ (G.E_red .* G.R_x);
-H.U = repmat(U_tmp,1,size(B.Q,2));
-EtimesR2_tmp = G.E_red .* G.R_x.^2;
-H.W = bsxfun(@rdivide,B.Q,EtimesR2_tmp);
-H.G = bsxfun(@times,H.alpha_p,G.E_red);
-k_tmp = repmat([G.k_i;G.k_a],1,size(B.Q,2));
+H.u_ir= (G.d_m.^2 - G.D_RE(ballMaterialInd).^2 .* cosd(G.alpha).^2) ./ (4*G.d_m) .* omega;
+H.u_ar = H.u_ir;
+H.u(:,1,:) = [H.u_ir; H.u_ar];
+E_red_tmp(1,1,:)=G.E_red(ballMaterialInd);
+R_x_temp(:,1,:)=G.R_x(:,ballMaterialInd);
+H.U = T.eta_0 .* H.u ./ (E_red_tmp .* R_x_temp);
+H.W =  Q_temp./(E_red_tmp .* R_x_temp.^2);
+H.G = H.alpha_p.*E_red_tmp;
+k_tmp = [G.k_i;G.k_a];
 H.H_0   = 2.69 * H.G.^0.53 .* H.U.^0.67 .* H.W.^-0.067 .* (1-0.61*exp(-0.73.*k_tmp));
 H.H_min = 3.63 * H.G.^0.49 .* H.U.^0.68 .* H.W.^-0.073 .* (1-     exp(-0.68.*k_tmp));
-H.h_0   = bsxfun(@times,H.H_0  ,G.R_x);
-H.h_min = bsxfun(@times,H.H_min,G.R_x);
-H.L = T.eta_0 .* T.alpha_etaT .* K.u.^2  ./  (4*S.lambda);
+H.h_0   = H.H_0 .*R_x_temp;
+H.h_min = H.H_min.*R_x_temp;
+H.L = T.eta_0 .* T.alpha_etaT .* H.u.^2  ./  (4*S.lambda);
 H.C_korr = 3.94  ./  (3.94 + H.L.^0.62);
-H.h_0th   = bsxfun(@times,H.h_0  ,H.C_korr);
-H.h_minth = bsxfun(@times,H.h_min,H.C_korr);
-clear U_tmp EtimesR2_tmp k_tmp
+H.h_0th   = H.h_0  .*H.C_korr;
+H.h_minth = H.h_min.*H.C_korr;
+
+B.s=zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
+
+B.s(1,B.noContactInd') = B.intersectionBall(B.noContactInd)';
+B.s(2,B.noContactInd') = B.intersectionBall(B.noContactInd)';
+
+B.s(1,B.  contactInd') = B.intersectionBall(B.contactInd)' + H.h_0th(1,B.contactInd)';
+B.s(2,B.  contactInd') = B.intersectionBall(B.contactInd)' + H.h_0th(2,B.contactInd)'; % s beschreibt den Abstand von WK zu Laufbahnen, unter Berücksichtigung des Schmierfilms in der Lastzone
 
 %% Attribute ändern
 obj.H = H;
+obj.B = B;
 obj.up2date.H = true;
 
 end
\ No newline at end of file
diff --git a/@BearImp/calcGeo.m b/@BearImp/calcGeo.m
index ccdc75dcd6711247a2745a00a17f39c34e8e58ca..d8c107caa64cc53ad286d4064abe25202e603989 100644
--- a/@BearImp/calcGeo.m
+++ b/@BearImp/calcGeo.m
@@ -1,132 +1,92 @@
 function calcGeo(obj)
-%2.3 LagerGeometrie G = G(L,R,F_a)
+%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
 %Fläche, Winkelpositionen der Wälzkörper, etc.
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 06.04.2021
-%    Autor: Steffen Puchtler
+%    Autor: Steffen Puchtler, Julius van der Kuip
 
 %% Parameter prüfen
-if obj.up2date.G
-    warning('Lagergeometrie in Abhängigkeit des Lagerspiels bereits berechnet')
-    return
-end
 assert(obj.up2date.L,'Lager nicht gesetzt')
 assert(obj.up2date.R,'Lagerspiel noch nicht berechnet')
-L = obj.L; R = obj.R; F_a = obj.F_a;
-G = struct;
+L = obj.L; R = obj.R; T_Oil = obj.T_Oil;
+G = struct; 
 
 %% Berechnung
-G.E_red = 2  ./  ((1-L.RE_A.nu.^2)./L.RE_A.E + (1-L.ring.nu.^2)/L.ring.E);
-G.B_m = L.f_a + L.f_i - 1;
-if F_a
-    if R.P_d > 0
-        G.alpha_0 = acosd(1 - R.P_d./(2.*G.B_m.*L.D_wk));
-    else
-        G.alpha_0 = 0;
-    end
-    alpha_old = min(G.alpha_0+10, 45);
-    G.K_a = spline([0 0.02 0.08 0.14 0.2],[0 350 1700 3400 5300].*1e6,G.B_m);
-    G.iterationen_alpha = 0;
-    tmp = @(alpha_alt) cosd(G.alpha_0)./cosd(alpha_alt) - 1;
-    alpha = @(alpha_alt,K_a) ((F_a./(L.Z.*L.D_wk.^2.*K_a)) - sind(alpha_old).*tmp(alpha_alt).^1.5)  ./  ...
-                         (cosd(alpha_alt).*tmp(alpha_alt).^1.5  +  1.5*tand(alpha_alt).^2 .* tmp(alpha_alt).^0.5 .* cosd(G.alpha_0));
-else
-    G.alpha = 0;
-end
+G.alpha = 0; 
 
-while true
-    if F_a 
-        G.alpha = fzero(@(alpha_alt) alpha(alpha_alt,G.K_a),alpha_old);
-    end
-    G.gamma = L.D_wk .* cosd(G.alpha) ./ L.d_m;
-    G.sigmarho_i = 1./L.D_wk .* (4 - 1./L.f_i + 2*G.gamma./(1-G.gamma));
-    G.sigmarho_a = 1./L.D_wk .* (4 - 1./L.f_a - 2*G.gamma./(1+G.gamma));
-    G.R_xi = (L.D_wk.*(L.d_m-L.D_wk.*cosd(G.alpha))) ./ (2*L.d_m);
-    G.R_xa = (L.D_wk.*(L.d_m+L.D_wk.*cosd(G.alpha))) ./ (2*L.d_m);
-    G.R_x  = [G.R_xi;G.R_xa];
-    G.R_yi = L.f_i.*L.D_wk ./ (2*L.f_i - 1);
-    G.R_ya = L.f_a.*L.D_wk ./ (2*L.f_a - 1);
-    G.R_si = G.R_xi.*G.R_yi  ./ (G.R_xi+G.R_yi);
-    G.R_sa = G.R_xa.*G.R_ya  ./ (G.R_xa+G.R_ya);
-    G.Gamma_i = G.R_si  .*  (1./G.R_xi - 1./G.R_yi);
-    G.Gamma_a = G.R_sa  .*  (1./G.R_xa - 1./G.R_ya);
-
-    % Iteration zur Bestimmung des Elliptizitätsparameters
-    k_alt = 0;
-    G.k_i = 2;
-    G.k_a = 2;
-    integralE = @(phi,k) (1 - (1-1/k.^2).*sin(phi).^2).^0.5;
-    integralF = @(phi,k) (1 - (1-1/k.^2).*sin(phi).^2).^-0.5;
-    G.iterationen_k = [0 0];
-    while abs(k_alt - G.k_i) > 1e-7
-        G.fraktE_i = integral(@(phi)integralE(phi,G.k_i),0,pi/2);
-        G.fraktF_i = integral(@(phi)integralF(phi,G.k_i),0,pi/2);
-        k_alt = G.k_i;
-        G.k_i = ((2*G.fraktF_i-G.fraktE_i.*(1+G.Gamma_i)) ./ (G.fraktE_i.*(1-G.Gamma_i))).^0.5;
-        G.iterationen_k(1) = G.iterationen_k(1) + 1;
-        if G.iterationen_k(1) > 1000
-            error('Elliptizitätsparameter innen konvergiert nicht');
-        end
-    end
-    k_alt = 0;
-    while abs(k_alt - G.k_a) > 1e-7
-        G.fraktE_a = integral(@(phi)integralE(phi,G.k_a),0,pi/2);
-        G.fraktF_a = integral(@(phi)integralF(phi,G.k_a),0,pi/2);
-        k_alt = G.k_a;
-        G.k_a = ((2*G.fraktF_a-G.fraktE_a.*(1+G.Gamma_a)) ./ (G.fraktE_a.*(1-G.Gamma_a))).^0.5;
-        G.iterationen_k(2) = G.iterationen_k(2) + 1;
-        if G.iterationen_k(2) > 1000
-            error('Elliptizitätsparameter außen konvergiert nicht');
-        end
-    end
-    clear k_alt integralE integralF
-
-    G.delta_sterni = 2*G.fraktF_i/pi .* (pi/(2*G.k_i.^2.*G.fraktE_i))^(1/3);
-    G.delta_sterna = 2*G.fraktF_a/pi .* (pi/(2*G.k_a.^2.*G.fraktE_a))^(1/3);
-    G.K_pi = 2.15e5 * (G.sigmarho_i*1e-3).^-0.5 .* (G.delta_sterni)^-1.5 * 10^4.5;
-    G.K_pa = 2.15e5 * (G.sigmarho_a*1e-3).^-0.5 .* (G.delta_sterna)^-1.5 * 10^4.5;
-    G.K_p = (1  /  (G.K_pi.^(-2/3) + G.K_pa.^(-2/3))).^1.5;
-        
-    if F_a
-        G.K_a = G.K_p .* G.B_m.^1.5 ./ L.D_wk.^0.5;
-    
-        G.iterationen_alpha = G.iterationen_alpha + 1;
-        if G.iterationen_alpha >= 1000
-            error('Druckwinkel konvergiert nicht');
-        end
-    
-        if abs(G.alpha - alpha_old) < 1e-6
-            clear alpha_old
-            break
-        end
-        alpha_old = G.alpha;
-    else
-        break
-    end
-end
-clear tmp alpha
-G.R_wk=L.D_wk/2; %Radius der Kugel
-G.R_bi=(L.d_m-L.D_wk-R.P_d/2)/2; %30.508/2; %Minimaler innerer Rillenradius
-G.R_ba=(L.d_m+L.D_wk+R.P_d/2)/2; %46.516/2; %Maximaler äußerer Rillenradius
-G.R_li=L.D_wk*L.f_i; %Radius des inneren Rillenmittelpunkts
-G.R_la=L.D_wk*L.f_a; %Radius des äußeren Rillenmittelpunkts
-
-G.Z_calc = L.Z + 1;
-if ~mod(L.Z,2) % gerade Anzahl WK
-    G.realInd = [1:L.Z/2+1  L.Z/2:-1:2  L.Z/2+2:L.Z+1  L.Z+1:-1:L.Z/2+2];
-    G.calcInd = [1:L.Z/2+1  L.Z+1:3/2*L.Z];
-else % ungerade Anzahl WK
-    G.realInd = [1:ceil(L.Z/2)  ceil(L.Z/2):-1:2 ceil(L.Z/2)+1:L.Z+1  L.Z:-1:ceil(L.Z/2)+1];
-    G.calcInd = [1:ceil(L.Z/2) L.Z+1:L.Z+ceil(L.Z/2)];
+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.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_wk * (1 + L.RE_A.alpha_t * ((G.Delta_Ti_norm + G.Delta_Ta_norm) / 2));
+
+if ~L.allBallsSameMaterial % Materialparameter für alle WK Materialien berechnet
+    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_wk * (1 + L.RE_B.alpha_t * ((G.Delta_Ti_norm + G.Delta_Ta_norm) / 2));
 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.d_m = (G.d + G.D)/2;
+G.R_RE = G.D_RE/2;
+
+b_li = (L.D-L.d-2*L.D_wk-L.c)/2; %Breite von Innenring (bezogen auf Durchmesser)
+b_la = b_li; 
+G.b_li = b_li*(1+L.ring.alpha_t*G.Delta_Ti_norm);
+G.b_la = b_la*(1+L.ring.alpha_t*G.Delta_Ta_norm);
+
+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.R_bi = (G.d + G.b_li) / 2; %Minimaler innerer Rillenradius
+G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius
+
+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.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.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_si = 1 ./ G.sigmarho_i;
+G.R_sa = 1 ./ G.sigmarho_a;
+
+G.gamma = G.D_RE  *cos(G.alpha) ./ G.d_m;
+
+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;
+
+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)));
+
+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);
+
+G.k_i = fzero(k_i_fkt,5);
+G.k_a = fzero(k_a_fkt,5);
+
+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.K_p = sqrt(2)/3.*pi.*G.E_red.*(G.fraktF_i.*(G.sigmarho_i/(G.k_i^2.*G.fraktE_i)).^(1/3)+ ...
+        G.fraktF_a.*(G.sigmarho_a/(G.k_a^2.*G.fraktE_a)).^(1/3)).^(-3/2);
+
 G.Deltapsi = 2*pi./L.Z;
 G.psi_ = (0:L.Z-1) .* G.Deltapsi;
-G.psi_(L.Z+1:2*L.Z) = ((0:L.Z-1) + 0.5) .* G.Deltapsi;
-G.psi = G.psi_(G.calcInd);
 
 %% Attribute ändern
 obj.G = G;
diff --git a/@BearImp/calcImp.m b/@BearImp/calcImp.m
index 869ee591aa4872aef7404feeb8f48e76f5a08c49..5d85bb3e2b035175b30cf34606cf93ce31b606c4 100644
--- a/@BearImp/calcImp.m
+++ b/@BearImp/calcImp.m
@@ -1,157 +1,78 @@
 function calcImp(obj)
-%2.6 Impedanz Z = Z(L,S,T,G,B,H,
+%2.7 Impedanz Z = Z(L,S,T,G,B,H,C)
 %Berechnet die Kapazität der Hertz'schen Flächen, sowie der Randbereiche
 %und unbelasteter Wälzkörper je nach gewählter Methode. Zusätzlich wir der
 %Widerstand der Hertz'schen Flächen bestimmt und anschließend die Impedanz
 %berechnet.
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 07.04.2021
-%    Autor: Steffen Puchtler, Tobias Schirra, Leander
+%    Autor: Steffen Puchtler, Tobias Schirra, Leander, Julius van der Kuip
 
 %% Parameter prüfen
-if obj.up2date.Z
-    warning('Impedanz bereits berechnet')
-    return
-end
 assert(obj.up2date.L,'Lager nicht gesetzt')
-assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
-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; Z = obj.Z;
-method = possibleMethods.addDefault(obj.method).Z;
-if ~isfield(Z,'k_C')
-    Z.k_C = 1;
-end
-if ~isfield(Z,'k_R')
-    Z.k_R = 1;
-end
-
+assert(obj.up2date.C,'Kapazitäten noch nicht berechnet')
+L = obj.L; C = obj.C; psi = obj.psi;
+Z = struct;
 %% Berechnung
-Z.rhoRatio = 1  +  (9.73e-3 .* (H.p.*0.101972e-6).^0.75)  ./  (T.nu_38.*1e6).^0.0385;
-Z.epsilon_primeOil = (S.epsilon_Oel + 2 + 2*(S.epsilon_Oel-1).*Z.rhoRatio)  ./  (S.epsilon_Oel + 2 - (S.epsilon_Oel-1).*Z.rhoRatio);
-Z.C_Hertz = zeros(2,G.Z_calc);
-Z.C_Hertz(:,B.contactInd) = obj.epsilon_0 .* Z.epsilon_primeOil .* H.A ./ H.h_0th;
-Z.R_el_single = inf(2,G.Z_calc);
-Z.R_el_single(:,B.contactInd) = (Z.k_R./(1+Z.k_R)) .* S.rho_el .* H.h_0th ./ H.A;
-Z.DeltaR_i(B.contactInd)    = G.R_li - G.R_wk - H.h_0th(1,:);
-Z.DeltaR_a(B.contactInd)    = G.R_la - G.R_wk - H.h_0th(2,:);
-if size(B.noContactInd,2) % Wenn Kraftfreie Wälzkörper vorhanden
-    Z.DeltaR_i(B.noContactInd) = G.R_li - G.R_wk - B.s(B.noContactInd);
-    Z.DeltaR_a(B.noContactInd) = G.R_la - G.R_wk - B.s(B.noContactInd);
-end
-Z.phi_1i = zeros(1,G.Z_calc);
-Z.phi_1i(B.contactInd) = H.a(1,:)./L.D_wk;
-Z.phi_1a = zeros(1,G.Z_calc);
-Z.phi_1a(B.contactInd) = H.a(2,:)./L.D_wk;
-Z.phi_2i = asin(L.B_i/L.D_wk);
-Z.phi_2a = asin(L.B_a/L.D_wk);
-if any([strcmp(method.unloadedRE,{'Leander_Radial','LeanderSteffen'}) strcmp(method.outsideArea,{'Leander_Radial','LeanderSteffen'})])
-    calcZg
-end
-Z.C_out = zeros(2,G.Z_calc);
 
-if size(B.noContactInd,2) % Wenn Kraftfreie Wälzkörper vorhanden
-    switch method.unloadedRE
-        case 'neglect'
-        case 'Leander_Parallel'
-            leander_parallel(B.noContactInd)
-        case 'Leander_Radial'
-            leander_radial(B.noContactInd)
-        case 'LeanderSteffen'
-            leander_steffen(B.noContactInd)
-        case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
-            tobias_steffen(B.noContactInd)
-    end
+% Überprüfen, ob alle C_el_single vorhanden
+if ~isnan(obj.psi_calc)
+   psi = obj.psi_calc;
 end
+check_all_nessecary_psi(psi,L.Z,L.allBallsSameMaterial)
 
-switch method.outsideArea
-    case 'neglect'
-    case 'k-factor'
-        Z.C_Hertz(:,B.contactInd) = (1+Z.k_C) .* obj.epsilon_0 .* Z.epsilon_primeOil .* H.A ./ H.h_minth;
-    case 'Leander_Parallel'
-        leander_parallel(B.contactInd)
-    case 'Leander_Radial'
-        leander_radial(B.contactInd)
-    case 'LeanderSteffen'
-        leander_steffen(B.contactInd)
-    case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
-        tobias_steffen(B.contactInd)
-end
+%Berechnung der approximierten Impedanz
+Z.omega =@(f) 2 * pi * f; %f = frequency of current 
+Z.phi_aprx =@(f) atand( -Z.omega(f) .* Z.C_el_aprx .* Z.R_el_aprx);
+Z.Z_el_aprx =@(f) abs(1./(1./ C.R_el_aprx + 1j * Z.omega(f) .* C.C_el_aprx));
 
-Z.C_el_single = Z.C_Hertz(:,G.realInd) + Z.C_out(:,G.realInd);
-Z.C_el(1) = 1./sum(1./sum(Z.C_el_single(:,    1:L.Z),2));
-Z.C_el(2) = 1./sum(1./sum(Z.C_el_single(:,L.Z+1:end),2));
-Z.C_el(3) = (Z.C_el(1)+Z.C_el(2)) / 2;
-Z.R_el_single = Z.R_el_single(:,G.realInd);
-Z.R_el(1) = sum(1./sum(1./Z.R_el_single(:,    1:L.Z),2));
-Z.R_el(2) = sum(1./sum(1./Z.R_el_single(:,L.Z+1:end),2));
-Z.R_el(3) = (Z.R_el(1)+Z.R_el(2)) / 2;
+%Berechung der exakten Impedanz 
+if strcmp(L.cage_electr,'conductor')
+    Z.Z_el_single_iext =@(f) 1./(1./ C.R_el_single(1,:,:) + 1j * 2*pi* f .* C.C_el_single(1,:,:));
+    Z.Z_el_single_aext =@(f) 1./(1./ C.R_el_single(2,:,:) + 1j * 2*pi* f .* C.C_el_single(2,:,:));
+    Z.Z_el_iext        =@(f) 1./ sum( 1 ./ Z.Z_el_single_iext(f),3);
+    Z.Z_el_aext        =@(f) 1./ sum( 1 ./ Z.Z_el_single_aext(f),3);
+    Z.Z_el             =@(f) Z.Z_el_iext(f) + Z.Z_el_aext(f);
+elseif strcmp(L.cage_electr,'isolator')
+    Z.Z_el_single_iext =@(f) squeeze( 1./(1./ C.R_el_single(1,:,:) + 1j * 2*pi* f .* C.C_el_single(1,:,:)) );
+    Z.Z_el_single_aext =@(f) squeeze( 1./(1./ C.R_el_single(2,:,:) + 1j * 2*pi* f .* C.C_el_single(2,:,:)) );
+    Z.Z_el_REext       =@(f) squeeze( Z.Z_el_single_iext(f) + Z.Z_el_single_aext(f) );
+    Z.Z_el             =@(f) squeeze( sum(Z.Z_el_REext(f),1) );
+    
+else
+    warning("Die elektrischen Eigenschaften des Käfigs sind falsch definiert! Es kann keine Gesamtinpedanz berechnet werden. ['conductor','isolator']")
+end
 
-% Z.omega = 2*pi*f;
-% Z.Z_el = R_el ./ sqrt(1+(omega.*C_el.*R_el).^2);
-% Z.phi = atand(-omega.*C_el.*R_el);
+reZ =@(f) real(Z.Z_el(f));
+imZ =@(f) imag(Z.Z_el(f));
 
+Z.R_el =@(f) (reZ(f) + imZ(f).^2) ./ reZ(f);
+Z.C_el =@(f) -imZ(f) ./ (Z.omega(f) .* (reZ(f).^2 + imZ(f).^2));
 %% Attribute ändern
 obj.Z = Z;
 obj.up2date.Z = true;
 
-    function calcZg
-        Z.g_1i = zeros(1,G.Z_calc);
-        Z.g_1a = zeros(1,G.Z_calc);
-        Z.g_1i(B.contactInd) = asin(sin(Z.phi_1i(B.contactInd))  .*  (-Z.DeltaR_i(B.contactInd)./G.R_li.*cos(Z.phi_1i(B.contactInd)) + sqrt(1-Z.DeltaR_i(B.contactInd).^2./G.R_li.^2.*sin(Z.phi_1i(B.contactInd)).^2)));
-        Z.g_1a(B.contactInd) = asin(sin(Z.phi_1a(B.contactInd))  .*  (-Z.DeltaR_a(B.contactInd)./G.R_la.*cos(Z.phi_1a(B.contactInd)) + sqrt(1-Z.DeltaR_a(B.contactInd).^2./G.R_la.^2.*sin(Z.phi_1a(B.contactInd)).^2)));
-        Z.g_2i = acos(L.B_i/(2*G.R_li));
-        Z.g_2a = acos(L.B_a/(2*G.R_la));
-    end
-    function leander_parallel(indices)
-        Z.R_hi = @(v) G.R_bi.*sqrt(1-(G.R_wk/G.R_bi.*cos(v)).^2);
-        Z.R_ha = @(v) G.R_ba.*sqrt(1-(G.R_wk/G.R_ba.*cos(v)).^2);
-        for ii = indices
-            Z.I_i = @(v,p) abs(G.R_wk^2.*sin(v))./(-(Z.R_hi(v)-G.R_li).*sqrt(1-(G.R_wk.*sin(p).*sin(v)./(Z.R_hi(v)+G.R_li)).^2)+G.R_wk.*(1+sin(v).*cos(p))-G.R_li+G.R_bi+B.s(ii));
-            Z.I_a = @(v,p) abs(G.R_wk^2.*sin(v))./( (Z.R_ha(v)+G.R_la).*sqrt(1-(G.R_wk.*sin(p).*sin(v)./(Z.R_ha(v)+G.R_la)).^2)+G.R_wk.*(1-sin(v).*cos(p))-G.R_la-G.R_ba+B.s(ii));
-            Z.C_out(1,ii) = 4 * obj.epsilon_0 * S.epsilon_Oel * integral2(Z.I_i,Z.phi_1i(ii)+pi/2,Z.phi_2i+pi/2, pi,1.5*pi);
-            Z.C_out(2,ii) = 4 * obj.epsilon_0 * S.epsilon_Oel * integral2(Z.I_a,Z.phi_1a(ii)+pi/2,Z.phi_2a+pi/2,  0,0.5*pi);
-        end
-    end
-    function leander_radial(indices)
-        Z.c_i =  G.R_li-G.R_bi-G.R_wk-B.s;          % TODO: B.s an der Stelle sinnvoll? ###
-        Z.c_a = -G.R_la-G.R_ba+G.R_wk+B.s;          % TODO: B.s an der Stelle sinnvoll? ###
-        for ii = indices
-            Z.I_i = @(g,t) abs(G.R_li*(G.R_li.*sin(g)+G.R_bi))./(sqrt(G.R_li.^2+G.R_bi.^2+Z.c_i(ii).^2+2*G.R_bi*G.R_li.*sin(g)+2*Z.c_i(ii)*(G.R_li.*sin(g)+G.R_bi).*cos(t))-G.R_wk);
-            Z.I_a = @(g,t) abs(G.R_la*(G.R_la.*sin(g)+G.R_ba))./(sqrt(G.R_la.^2+G.R_ba.^2+Z.c_a(ii).^2+2*G.R_ba*G.R_la.*sin(g)+2*Z.c_a(ii)*(G.R_la.*sin(g)+G.R_ba).*cos(t))-G.R_wk);
-            Z.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(Z.I_i,-pi/2+Z.g_1i(ii),  -Z.g_2i, 0,pi/9);
-            Z.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(Z.I_a, pi/2+Z.g_1a(ii),pi-Z.g_2a, 0,pi/9);
-        end
-    end
-    function leander_steffen(indices)
-        Z.dA_i = @(p) G.R_li * (G.R_bi + G.R_li*cos(p));
-        Z.dA_a = @(p) G.R_la * (G.R_ba + G.R_la*cos(p));
-        Z.x_kmi = G.R_bi + Z.DeltaR_i;
-        Z.x_kma = G.R_ba - Z.DeltaR_a;
-        Z.h_i = @(p,t,ii) sqrt(G.R_bi^2 + G.R_li^2 + Z.x_kmi(ii)^2 + 2*G.R_bi*G.R_li*cos(p) - 2*Z.x_kmi(ii)*cos(t).*(G.R_bi+G.R_li*cos(p)))  -  G.R_wk;
-        Z.h_a = @(p,t,ii) sqrt(G.R_ba^2 + G.R_la^2 + Z.x_kma(ii)^2 + 2*G.R_ba*G.R_la*cos(p) - 2*Z.x_kma(ii)*cos(t).*(G.R_ba+G.R_la*cos(p)))  -  G.R_wk;
-        for ii = indices
-            Z.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_i(p)./Z.h_i(p,t,ii), Z.g_1i(ii)   , Z.g_2i   , 0, pi/9);
-            Z.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_a(p)./Z.h_a(p,t,ii), Z.g_1a(ii)+pi, Z.g_2a+pi, 0, pi/9);
-        end
-    end
-    function tobias_steffen(indices)
-        Z.h_i = @(alpha,beta,ii) (sqrt(G.R_li^2 - Z.DeltaR_i(ii)^2*sin(alpha).^2) - Z.DeltaR_i(ii)*cos(alpha))  ./  cos(beta) - G.R_wk;
-        Z.h_a = @(alpha,beta,ii) (sqrt(G.R_la^2 - Z.DeltaR_a(ii)^2*sin(alpha).^2) - Z.DeltaR_a(ii)*cos(alpha))  ./  cos(beta) - G.R_wk;
-        if strcmp(method.unloadedRE,'TobiasSteffen_Kugelfläche')
-            Z.dA_i = @(~,beta,~) G.R_wk^2 * cos(beta);
-            Z.dA_a = Z.dA_i;
+%% Helper Functions
+function check_all_nessecary_psi(psi,Z,symCase)
+        
+        if symCase
+            psi_calc = [psi, 2*pi - psi];
         else
-            Z.dA_i = @(alpha,beta,ii) G.R_li./cos(beta) .* (1 - Z.DeltaR_i(ii)*cos(alpha)./sqrt(G.R_li^2-Z.DeltaR_i(ii)^2*sin(alpha).^2)) .* (G.R_wk + Z.h_i(alpha,beta,ii));
-            Z.dA_a = @(alpha,beta,ii) G.R_la./cos(beta) .* (1 - Z.DeltaR_a(ii)*cos(alpha)./sqrt(G.R_la^2-Z.DeltaR_a(ii)^2*sin(alpha).^2)) .* (G.R_wk + Z.h_a(alpha,beta,ii));
+            psi_calc = psi;
         end
-        for ii = indices
-            Z.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) Z.dA_i(alpha,beta,ii)./Z.h_i(alpha,beta,ii), Z.phi_1i(ii), Z.phi_2i, 0, pi/3);
-            Z.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) Z.dA_a(alpha,beta,ii)./Z.h_a(alpha,beta,ii), Z.phi_1a(ii), Z.phi_2a, 0, pi/3);
+        
+        for run_C = 1 : length(obj.psi)
+            
+           psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj.psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
+
+            psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) = psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
+           
+            psiNeedForIndCapa = round(psiNeedForIndCapa,8);
+
+            [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-8);
+
+            %Check, if all nessecary C_el_single are calculated
+            assert(~any(ind_single_psi==0),'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!') 
         end
-    end
+end
 end
\ No newline at end of file
diff --git a/@BearImp/calcLoad.m b/@BearImp/calcLoad.m
index af85f86062688c3d1384ac65a1104f89c45a5ab0..1148c8219064abe938b80df053d793b9ac0a2544 100644
--- a/@BearImp/calcLoad.m
+++ b/@BearImp/calcLoad.m
@@ -1,74 +1,197 @@
 function calcLoad(obj)
-%2.4 Belastungsverteilung B = B(L,R,G,F_r,F_a)
+%2.4 Belastungsverteilung B = B(L,G,F_r,psi)
 % Berechnet die Belastungsverteilung auf die einzelnen Wälzkörper.
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 07.04.2021
-%    Autor: Steffen Puchtler
+%    Autor: Steffen Puchtler, Julius van der Kuip
 
 %% Parameter prüfen
-if obj.up2date.B
-    warning('Belastungsverteilung bereits berechnet')
-    return
-end
 assert(obj.up2date.L,'Lager nicht gesetzt')
-assert(obj.up2date.R,'Lagerspiel noch nicht berechnet')
 assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht berechnet')
-L = obj.L; R = obj.R; G = obj.G; F_r = obj.F_r; F_a = obj.F_a;
-B = struct;
+L = obj.L; G = obj.G; F_r = obj.F_r; AddOn = obj.AddOn; psi=obj.psi;
+B = struct; 
+
 
 %% Berechnung
-if F_a
-    Fratio = F_r/F_a*tand(G.alpha);
-    if Fratio >= 1
-        error('Axialkraft zu gering');
+    
+if ~isnan(obj.psi_calc) 
+    psi = obj.psi_calc; % alle folgenden Rechnungen werden mit dem reduzierten Psi berechnet, um alle benötigten Werte zu berechnen
+end
+
+B.numOfCagePositions = length(psi);
+B.delta_s=nan(B.numOfCagePositions,2);
+
+if L.allBallsSameMaterial
+
+    B.delta_s=deltaShaft(psi,F_r,B,G,L,B.delta_s,AddOn)';
+     
+    for IndPosBearing=1:B.numOfCagePositions
+        B.Q(IndPosBearing) = totalReactionForcePerBallOnlyY(G,L.IndBall_conductive,B.delta_s(:,IndPosBearing),psi(IndPosBearing));
     end
-    Table_Fratio   = [1 0.9318 0.8964 0.8601 0.8225 0.7835 0.7427 0.6995 0.6529 0.6 0.4538 0.3088 0.185 0.0831   0.0394   0.0192  0.0076];
-    Table_epsilon  = [0 0.2    0.3    0.4    0.5    0.6    0.7    0.8    0.9    1   1.25   1.67   2.5   5       10       20      50     ];
-    B.epsilon_0 = interp1(Table_Fratio,Table_epsilon,Fratio,'linear','extrap');
-
-    psi_1 = @(eps) acos(1-2*min(max(eps,0),1));
-    J_r = @(eps) 1/2/pi*integral(@(psi) (1 - 1./2./eps.*(1-cos(psi))).^1.5 .* cos(psi),-psi_1(eps),psi_1(eps));
-    J_a = @(eps) 1/2/pi*integral(@(psi) (1 - 1./2./eps.*(1-cos(psi))).^1.5            ,-psi_1(eps),psi_1(eps));
-    B.epsilon = fzero(@(eps) J_r(eps) - J_a(eps).*Fratio,B.epsilon_0);
-    B.J_r = J_r(B.epsilon);
-    B.psi_1 = psi_1(B.epsilon);
-    B.Q_max = F_r ./ (L.Z.*B.J_r.*cosd(G.alpha));
-    clear Fratio Table_Fratio Table_epsilon psi_1 J_r J_a
+    
 else
-    B.delta_r = abs(R.P_d);
-    delta_r_old = 0;
-    integralJ = @(psi,epsilon) (1 - 1./2./epsilon.*(1-cos(psi))).^1.5 .* cos(psi);
-    B.iterations_delta = 0;
-    while abs(delta_r_old - B.delta_r) > 1e-13
-        B.psi_1 = acos(R.P_d./2./B.delta_r);
-        B.epsilon = 1/2  *  (1 - (R.P_d./2./B.delta_r));
-        B.J_repsilon = 1/2/pi*integral(@(psi) integralJ(psi,B.epsilon),-B.psi_1,B.psi_1);
-        delta_r_old = B.delta_r;
-        B.delta_r = (F_r ./ (L.Z.*G.K_p.*B.J_repsilon)).^(2/3) + 0.5*R.P_d;
-        B.iterations_delta = B.iterations_delta + 1;
+
+    B.delta_s=deltaShaft(psi,F_r,B,G,L,B.delta_s,AddOn)';
+
+    B.Q=zeros(L.numberOfConductiveBalls,B.numOfCagePositions);
+    
+    for posBall_conductive=1:L.numberOfConductiveBalls
+
+        B.psi_conductive(posBall_conductive,:)=[psi(1+round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z):B.numOfCagePositions) psi(1: ...
+                                    (round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z)))];
+        B.extendedTotalReactionForce=true;
+        
+        for IndPosBearing=1:B.numOfCagePositions
+            B.Q(posBall_conductive,IndPosBearing)=totalReactionForcePerBallForYandX(B,L,G,B.delta_s(:,IndPosBearing),posBall_conductive,B.psi_conductive(posBall_conductive,IndPosBearing));
+        end
     end
-    clear delta_r_old integralJ
-    B.Q_max = G.K_p .* (B.delta_r-R.P_d/2).^1.5;
 end
 
+B.  contactInd =  B.Q~=0;  % Indizes der Wälzkörper, die belastet sind (Q ungleich 0)
+B.noContactInd =  B.Q==0;  % Indizes der Wälzkörper, die nicht belastet sind
+B.Q_contInd=B.Q; %##um leichter B.Q_contInd anzupassen, falls doch nur ~=0 berechnet werden sollen
+
+    for posBall_conductive=1:L.numberOfConductiveBalls
+        B.psi_conductive(posBall_conductive,:)=[psi(1+round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z):B.numOfCagePositions) psi(1: ...
+                                    (round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z)))];
+        for IndPosBall=1:B.numOfCagePositions
+            B.intersectionBall(posBall_conductive,IndPosBall) = deflectionForYandX(B,L,G,posBall_conductive,B.delta_s(:,IndPosBall),B.psi_conductive(posBall_conductive,IndPosBall))./2;
+        end %Beschreibt die Überschneidung der WK mit den Laufbahnen (<0 -> Kraft wird übertragen; >0 -> keine Kraft wird übertragen)
+    end
 
-B.Q = (B.Q_max .* (1 - 1/2./B.epsilon .* (1 - cos(G.psi))).^1.5) .* ((abs(G.psi)<=B.psi_1)|(abs(G.psi-2*pi)<=B.psi_1));
-B.  contactInd = find( B.Q);  % Indizes der Wälzkörper, die belastet sind (Q ungleich 0)      ################## wäre logical(B.Q) nicht sinnvoller?
-B.noContactInd = find(~B.Q);  % Indizes der Wälzkörper, die nicht belastet sind
-if size(B.noContactInd,2)     % Wenn Kraftfreie Wälzkörper vorhanden
-    B.Q = B.Q(B.contactInd);  % Unbelastete Kontakte vorläufig eliminieren
-    assert(~F_a,'Unbelastete Wälzkörper aktuell ausschließlich für reine Radiallasten implementiert') % da B.delta_r für kombinierte Lasten unbekannt
-    B.s(B.noContactInd) = 0.5*(-B.delta_r.*cos(G.psi(B.noContactInd)) + sqrt(G.R_ba^2-B.delta_r^2.*sin(G.psi(B.noContactInd)).^2)-G.R_bi-L.D_wk);
-end
-% Visualisierung der Lastverteilung
-% psi_kont = linspace(0,2*pi,73);
-% Q_kont = (B.Q_max .* (1 - 1/2./B.epsilon .* (1 - cos(psi_kont))).^1.5) .* ((abs(psi_kont)<=B.psi_1)|(abs(psi_kont-2*pi)<=B.psi_1));
-% polar(psi_kont,Q_kont);
 
 %% Attribute ändern
 obj.B = B;
 obj.up2date.B = true;
 
-end
\ No newline at end of file
+end
+
+
+function delta_s=deltaShaft(psi,F_r,B,G,L,delta_s,AddOn)
+
+    for I=1:B.numOfCagePositions
+        B.psi_run=G.psi_+psi(I);
+        delta_s=solveEquilibrOfForces(I,F_r,B,G,L,delta_s,AddOn);
+    end
+end
+
+function delta_s=solveEquilibrOfForces(I,F_r,B,G,L,delta_s,AddOn) 
+
+    start= [0.0001;0.00001];
+   
+    if AddOn.Optimization_Toolbox
+        options = optimoptions('fsolve', 'StepTolerance', 1e-20, 'FunctionTolerance', 1e-20, 'FunValCheck', 'on', 'OptimalityTolerance',1e-20,'Display','off');%, 'SpecifyObjectiveGradient', true );%,'Algorithm', 'levenberg-marquardt' );%,"PlotFcns",@optimplotfval);
+        delta_s(I,:) = fsolve(@(delta_s) EoF(delta_s,B,G,L,F_r,AddOn), start,options);
+        
+    else
+         
+        options = optimset('MaxFunEvals',10000,"MaxIter", 10000,"TolFun", 1e-8,"TolX", 1e-8);%,"PlotFcns",@optimplotfval);
+        delta_s(I,:) = fminsearch(@(delta_s) EoF(delta_s,B,G,L,F_r,AddOn), start,options);
+    end
+    
+end
+
+function EoF = EoF(delta_s,B,G,L,F_r,AddOn)
+
+    if delta_s(2) == 0
+        
+        delta_s = delta_s(1);
+        
+        EoF = abs(sum(yReactionForceOnlyY(B,G,L,delta_s))-F_r);
+
+    elseif AddOn.Optimization_Toolbox
+         EoF= [F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s)); ...
+            sum(xReactionForcePerBallForYandX(B,G,L,delta_s))];
+    else
+        EoF = abs(F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s)))+ ...
+                abs(sum(xReactionForcePerBallForYandX(B,G,L,delta_s)));
+    end
+end
+
+
+
+function Q_y_y=yReactionForceOnlyY(B,G,L,delta_s)
+    Q_y_y=zeros(1,L.Z);
+    
+    for IndBall=1:L.Z
+            Q_y_y(IndBall)=yReactionForcePerBallOnlyY(B,G,IndBall,delta_s);
+    end   
+end
+
+function Q_y_y_perball_cont=yReactionForcePerBallOnlyY(B,G,IndBall,delta_s)
+   Q_y_y_perball_cont = totalReactionForcePerBallOnlyY(G,IndBall,delta_s,B.psi_run)* ...
+                       cos(deflectionAngleOnlyY(B,G,IndBall,delta_s));
+end
+
+function Q_total_y_perball_cont=totalReactionForcePerBallOnlyY(G,IndBall,delta_s,psi_run)
+   Q_total_y_perball_cont = abs((deflectionOnlyY(G,IndBall,delta_s,psi_run)<0)* ...
+                            (deflectionOnlyY(G,IndBall,delta_s,psi_run)))^(3/2)*G.K_p;
+end
+
+function delta_nury=deflectionOnlyY(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;
+end
+
+function psi_k_y=deflectionAngleOnlyY(B,G,IndBall,delta_s)
+    psi_k_y = B.psi_run(IndBall)+atan(sin(B.psi_run(IndBall))*delta_s/(G.R_ba-G.R_RE)-cos(B.psi_run(IndBall))*delta_s);
+end
+
+
+
+
+function Q_y_yandx_perball=yReactionForcePerBallForYandX(B,G,L,delta_s)
+    Q_y_yandx_perball=zeros(1,L.Z);
+    
+    for IndBall=1:L.Z
+         Q_y_yandx_perball(IndBall)=totalReactionForcePerBallForYandX(B,L,G,delta_s,IndBall,B.psi_run)* ...
+         cos(deflectionAngleForYandX(B,L,G,IndBall,delta_s));
+    end
+end
+
+function Q_total_yandx_perball=totalReactionForcePerBallForYandX(B,L,G,delta_s,IndBall,psi_run)
+
+    if isfield(B,'extendedTotalReactionForce')
+        Q_total_yandx_perball=abs(((deflectionForYandX(B,L,G,IndBall,delta_s,psi_run)<0)).* ...
+                               deflectionForYandX(B,L,G,IndBall,delta_s,psi_run)).^(3/2).*G.K_p(L.RE_order_num(L.IndBall_conductive(IndBall)));
+    else
+        Q_total_yandx_perball=abs(((deflectionForYandX(B,L,G,IndBall,delta_s,psi_run)<0)).* ...
+                               deflectionForYandX(B,L,G,IndBall,delta_s,psi_run)).^(3/2).*G.K_p(L.RE_order_num(IndBall));
+
+    end
+end
+
+function delta_yandx=deflectionForYandX(B,L,G,IndBall,delta_s,psi_run)
+
+    if isfield(B,'extendedTotalReactionForce')
+        delta_yandx = sin(psi_run-atan(delta_s(2,:)./delta_s(1,:))).*delta_s(2,:)./...
+                            sin(atan(delta_s(2,:)./delta_s(1,:)))./(sin(atan(sin(psi_run-atan(delta_s(2,:)./delta_s(1,:))).*delta_s(2,:)./...
+                            sin(atan(delta_s(2,:)./delta_s(1,:)))./(G.R_ba-G.R_RE(L.RE_order_num(L.IndBall_conductive(IndBall)))-cos(psi_run-atan(delta_s(2,:)./delta_s(1,:))).* ...
+                            delta_s(2,:)./sin(atan(delta_s(2,:)./delta_s(1,:)))))))-G.R_bi-G.R_RE(L.RE_order_num(L.IndBall_conductive(IndBall)));
+
+    else
+        delta_yandx = sin(psi_run(IndBall)-atan(delta_s(2,:)./delta_s(1,:))).*delta_s(2,:)./...
+                        sin(atan(delta_s(2,:)./delta_s(1,:)))./(sin(atan(sin(psi_run(IndBall)-atan(delta_s(2,:)./delta_s(1,:))).*delta_s(2,:)./...
+                        sin(atan(delta_s(2,:)./delta_s(1,:)))./(G.R_ba-G.R_RE(L.RE_order_num(IndBall))-cos(psi_run(IndBall)-atan(delta_s(2,:)./delta_s(1,:))).* ...
+                        delta_s(2,:)./sin(atan(delta_s(2,:)./delta_s(1,:)))))))-G.R_bi-G.R_RE(L.RE_order_num(IndBall));
+
+    end        
+end
+
+function psi_k_yandx=deflectionAngleForYandX(B,L,G,IndBall,delta_s)
+
+    psi_k_yandx = B.psi_run(IndBall)+atan(sin(B.psi_run(IndBall)-atan(delta_s(2)/delta_s(1)))*delta_s(2)/...
+                    sin(atan(delta_s(2)/delta_s(1)))/(G.R_ba-G.R_RE(L.RE_order_num(IndBall))-cos(B.psi_run(IndBall)-atan(delta_s(2)/ ...
+                    delta_s(1)))*delta_s(2)/sin(atan(delta_s(2)/delta_s(1))))); %(bei radialer und axialer Verschiebung)
+
+end
+
+
+function Q_x_yandx_perball=xReactionForcePerBallForYandX(B,G,L,delta_s)    
+    Q_x_yandx_perball=zeros(1,L.Z);
+
+    for IndBall=1:L.Z
+        Q_x_yandx_perball(IndBall)=totalReactionForcePerBallForYandX(B,L,G,delta_s,IndBall,B.psi_run)*sin(deflectionAngleForYandX(B,L,G,IndBall,delta_s));
+    end
+end
diff --git a/@BearImp/calcLub.m b/@BearImp/calcLub.m
index 8f0cc3bc388a39415cbd227aecbabc5bb962698e..58c25ebbc36ae73892e7256c309c4a703ee8e0d8 100644
--- a/@BearImp/calcLub.m
+++ b/@BearImp/calcLub.m
@@ -3,15 +3,9 @@ function calcLub(obj)
 %Berechnet die temperaturabhängigen Schmierstoffparameter
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 06.04.2021
 %    Autor: Steffen Puchtler
 
 %% Parameter prüfen
-if obj.up2date.T
-    warning('Schmierstoffparameter sind bereits berechnet')
-    return
-end
 assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
 S = obj.S; T_Oil = obj.T_Oil;
 T = struct;
diff --git a/@BearImp/calculate.m b/@BearImp/calculate.m
index fc1f29e74a82a82f0265defcbf81ccd0b02f9c0c..e4a66902fd87cb9dbaeb0f3b47861d1b39d2348e 100644
--- a/@BearImp/calculate.m
+++ b/@BearImp/calculate.m
@@ -1,38 +1,101 @@
-function calculate (obj)
+function calculate (obj, forceMode)
 %Berechnet alle Rechnungsabschnitte, die noch nicht berechnet wurden,
 %sofern alle Eingangsgrößen gesetzt sind.
 
+%forcemode
+%   - 'nonForced' -> Bereits berechnete und noch aktuelle größen werden nicht erneut berechnet (default)
+%   - 'forced'    -> Alle Größen werden erneut errechnet
+
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 28.04.2021
 %    Autor: Steffen Puchtler
 
+    arguments
+       obj
+       forceMode = 'nonForced'
+    end
+
     if ~isscalar(obj)
         executeAllObj(obj,@calculate)
         return
     end
     
-    if obj.allUp2date
+    switch forceMode
+    case 'forced'
+        force = true;
+    case 'nonForced'
+        force = false;
+    otherwise
+        error('forceMode not defined. Please select between ''forced'' and ''nonForced''.')
+    end
+    
+    if obj.allUp2date && ~force
         warning('Alle Werte bereits berechnet')
         return
     end
     assert(obj.ready2calc,'Noch nicht alle Eingangsparameter gesetzt')
-    if ~obj.up2date.T
+    
+    obj.psi_calc = find_nessecary_psi(obj.psi,obj.L.Z,obj.L.allBallsSameMaterial);
+    
+    if ~obj.up2date.T || force
         obj.calcLub
     end
-    if ~obj.up2date.R
+    if ~obj.up2date.R || force
         obj.calcClear
     end
-    if ~obj.up2date.G
+    if ~obj.up2date.G || force
         obj.calcGeo
     end
-    if ~obj.up2date.B
+    if ~obj.up2date.B || force
         obj.calcLoad
     end
-    if ~obj.up2date.H
+    if ~obj.up2date.H || force
         obj.calcFilm
     end
-    if ~obj.up2date.Z
+    if ~obj.up2date.C || force
+        obj.calcCap
+    end
+    if ~obj.up2date.Z || force
         obj.calcImp
     end
+    
+    function psi_calc = find_nessecary_psi(psi,Z, symCase) %(angle, number of balls)
+        %this function creates a vector with all positions necessary to
+        %calculate the capacity of the bearing
+        
+        switch symCase
+            case 'useSymmetry'
+                symmetry = true;
+            case 'noSymmetry'
+                symmetry = false;
+            case obj.L.allBallsSameMaterial % hier wird durch die Eigenschaft, ob alle RE das selbe Material besitzen, die Symmetrie-Eigenschaft des Lagers festgelegt
+                symmetry = obj.L.allBallsSameMaterial;
+            otherwise
+                error('Falscher Eingabewert für symCase! Mögliche Eingabewerte: ''useSymmetry'' oder ''noSymmetry''')
+        end
+        
+        
+        psi_calc=[];
+
+        for indPsi = 1 : length(psi)
+
+            psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + psi(indPsi); % Bsp. Ergebnis: [2,42,82,122,162,202,242,282,322]/360*2*pi (Offset 2 Grad) im Abstand der WK und dem gewünschten Offsetz von psi
+
+            psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) = psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
+
+            psi_calc = [psi_calc; psiNeedForIndCapa];
+        end
+        psi_calc=round(psi_calc,8);
+        psi_calc = uniquetol(psi_calc, 1e-8); %doppelte Werte löschen
+        
+        
+        if symmetry %Es entsteht eine Symmetrie in der vertikalen Achse, wodurch die Berechnung von zB. der Kapazität nur für eine Hälfte der Umdrehung berechnet werden muss
+                psi_calc(psi_calc >  pi) = 2*pi - psi_calc(psi_calc > pi);
+                tempPsiComp              = 2*pi - psi_calc(psi_calc<pi - eps(1e8)); % um Rechenungenauigkeiten durch irrationale Zahl zu vermeiden, wird der Filten erweitert
+                
+                psi_calc = uniquetol(psi_calc, 1e-8);
+                psi_calc(ismembertol(psi_calc,tempPsiComp,1e-8)) = []; %gespiegelte Punkte löschen 
+        end
+    end
+
+    
 end
\ No newline at end of file
diff --git a/@BearImp/copyToArray.m b/@BearImp/copyToArray.m
index 84a71e69b4b0fe369fbfb6f33374745107625326..0c3938f578738de581a3a9dc02a9b27b9894bd6b 100644
--- a/@BearImp/copyToArray.m
+++ b/@BearImp/copyToArray.m
@@ -4,8 +4,6 @@ function objArray = copyToArray(obj,varargin)
 %objArray = obj.copyToArray(m,n)    > m x n Matrix
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 12.05.2021
 %    Autor: Steffen Puchtler
 
     switch nargin
diff --git a/@BearImp/get.m b/@BearImp/get.m
index aba3ffe55816ce7fe3de2a73e9a71e6112f70a44..8534fd1d3fc6d1bb2a6e0dd205a30ebaf54e056e 100644
--- a/@BearImp/get.m
+++ b/@BearImp/get.m
@@ -3,8 +3,6 @@ function value = get(obj, name)
 %auch mit BearImp-Arrays bzw. -Matrizen
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 12.05.2021
 %    Autor: Steffen Puchtler
 
     arguments
diff --git a/@BearImp/set.m b/@BearImp/set.m
index 655c3808975d9456372b85f45a67a1c9d3be07e3..474b8684afa74c36c3dfe76a2d250419351806ef 100644
--- a/@BearImp/set.m
+++ b/@BearImp/set.m
@@ -3,8 +3,6 @@ function set(obj, name, value)
 %auch mit BearImp-Arrays bzw. -Matrizen
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 12.05.2021
 %    Autor: Steffen Puchtler
 
     arguments
diff --git a/@BearImp/setBearing.m b/@BearImp/setBearing.m
index 7babfa974b82372284a5f42930b90733ff22a61e..566568beba1114cba62a14eb66067931d0f5c9f6 100644
--- a/@BearImp/setBearing.m
+++ b/@BearImp/setBearing.m
@@ -3,9 +3,7 @@ function setBearing(obj,name)
 %setBearing(name) -> lädt das angegebene Lager
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 17.06.2021
-%    Autor: Steffen Puchtler
+%    Autor: Steffen Puchtler, Julius van der Kuip
 
     arguments
         obj
@@ -49,6 +47,31 @@ function setBearing(obj,name)
 
     L.d_m = (L.d + L.D)/2;
 
+     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
+                                       % kann umgangen werden, indem RE_order mit nur einem Buchstaben befüllt wird
+        L.IndBall_conductive=1;        % da bei psi=0 sich ein RE befindet und alle RE leitend sind
+        L.numberOfConductiveBalls=1;
+     else
+        L.allBallsSameMaterial = false;
+        L.RE_order_num = double(L.RE_order)-64;   
+        
+        if L.RE_A.electric_behaviour == "conductor"
+            OrderBall_conductive(find(L.RE_order_num==1))=find(L.RE_order_num==1);
+        end
+        if isfield(L,'RE_B') && L.RE_B.electric_behaviour == "conductor"
+            OrderBall_conductive(find(L.RE_order_num==2))=find(L.RE_order_num==2);
+        end
+        if isfield(L,'RE_C') && L.RE_C.electric_behaviour == "conductor"
+            OrderBall_conductive(find(L.RE_order_num==3))=find(L.RE_order_num==3);
+        end
+        
+        L.IndBall_conductive = nonzeros(OrderBall_conductive)';
+        L.numberOfConductiveBalls = length(L.IndBall_conductive);
+    
+    end
+    
     L = orderfields(L);
     obj.set('L',L)
     obj.setFitting('drawing')
diff --git a/@BearImp/setFitting.m b/@BearImp/setFitting.m
index 01bfbaf3cd04b0eeca357d3dea6057cdb186f356..8085b66eb3efa15d661b7a0bc5f7169b857c87b0 100644
--- a/@BearImp/setFitting.m
+++ b/@BearImp/setFitting.m
@@ -4,8 +4,6 @@ function setFitting(obj,shaft,housing)
 %setFitting(shaft,housing) -> lädt die angegebene Welle & Gehäuse
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 17.06.2021
 %    Autor: Steffen Puchtler
 
     arguments
@@ -70,7 +68,7 @@ end
 
 function setHousingTolerances(obj,fittingTable,index)
     obj.L.EI_D = fittingTable(index,:).Delta_lower_housing;
-    obj.L.ES_D = fittingTable(index,:).Delta_lower_housing;
+    obj.L.ES_D = fittingTable(index,:).Delta_upper_housing;
     obj.L.D_G  = fittingTable(index,:).D_G;
 end
     
\ No newline at end of file
diff --git a/@BearImp/setLube.m b/@BearImp/setLube.m
index 62c59ecfac9fc6414909865b521d6172e2542a1d..c3daf9d5ea0b459c864e45d99ee81e671e0e8906 100644
--- a/@BearImp/setLube.m
+++ b/@BearImp/setLube.m
@@ -3,8 +3,6 @@ function setLube(obj,name)
 %setLube(name) -> lädt den angegebenen Schmierstoff
 
 %    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 16.06.2021
 %    Autor: Steffen Puchtler
 
     arguments
diff --git a/Compare2OnlineCalc.mlx b/Compare2OnlineCalc.mlx
new file mode 100644
index 0000000000000000000000000000000000000000..599851ce2d991408f64d22ebc0293cd6fd6b4f11
Binary files /dev/null and b/Compare2OnlineCalc.mlx differ
diff --git a/InputData.xlsx b/InputData.xlsx
index 84f4701b50666be93a4d290855fc8c5099309528..039bf4744b57d3239b440a857bace73ba5a70007 100644
Binary files a/InputData.xlsx and b/InputData.xlsx differ
diff --git a/possibleMethods.m b/possibleMethods.m
index 05b23cc5183d17bacd83cd11e00beb4f8308dd9d..972b21550285c58edbfabdf70b1e2079db8d03ce 100644
--- a/possibleMethods.m
+++ b/possibleMethods.m
@@ -4,9 +4,9 @@ classdef possibleMethods
 % überprüfen des Methodenstructs
 
 % pmd Berechnungstool Lagerimpedanz
-% Version: 1.0.0
-% Stand: 22.04.2021
-% Autor: Steffen Puchtler
+% Version: 2.0.0
+% Stand: 02.01.2022
+% 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
@@ -24,16 +24,16 @@ classdef possibleMethods
         function s = H
             s = {'Hamrock/Dowson','Dowson/Higginson','Moes'};  % Beispiel, noch nicht implementiert
         end
-        function s = Z
-            s.unloadedRE  = {           'neglect','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'};
-            s.outsideArea = {'k-factor','neglect','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'};
+        function s = C
+            s.unloadedRE  = {           'neglect','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
+            s.outsideArea = {'k-factor','neglect','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.H = 'Hamrock/Dowson';
-            s.Z.unloadedRE  = 'Leander_Radial';
-            s.Z.outsideArea = 'Leander_Radial';
+            s.C.unloadedRE  = 'semianalytisch3D';
+            s.C.outsideArea = 'semianalytisch3D';
         end
         
         function method = addDefault(method)
diff --git a/testAll.m b/testAll.m
deleted file mode 100644
index f9d504501a55cbd1c1dbf32dd27c18b74cadd0e6..0000000000000000000000000000000000000000
--- a/testAll.m
+++ /dev/null
@@ -1,52 +0,0 @@
-% Skript, um die Funktionsfähigkeit von BearImp zu testen. Muss vor jedem
-% Merge ausgeführt werden. Neue Berechnungsmethoden etc. müssen hier
-% ergänzt werden.
-
-%    pmd Berechnungstool Lagerimpedanz
-%    Version: 1.0.0
-%    Stand: 07.05.2021
-%    Autor: Steffen Puchtler
-clear b
-
-numTests = 4;
-b(1,numTests) = BearImp;   % leeren Vektor erzeugen
-b.setBearing('6205-C3')
-b.setLube('FVA Referenz-Öl II')
-b.set('F_r'  , [1000  4000 7000 1000])
-b.set('F_a'  , [   0     0    0 1000])
-b.set('n'    , [  20    20   20  150])
-b.set('T_Oil', [  20    60   80   60])
-
-method(1).Z.unloadedRE  = 'neglect';
-method(1).Z.outsideArea = 'Leander_Parallel';
-method(2).Z.unloadedRE  = 'Leander_Radial';
-method(2).Z.outsideArea = 'LeanderSteffen';
-method(3).Z.unloadedRE  = 'TobiasSteffen_Kugelfläche';
-method(3).Z.outsideArea = 'TobiasSteffen_Laufbahnfläche';
-method(4).Z.unloadedRE  = 'neglect';
-method(4).Z.outsideArea = 'k-factor';
-b.set('method',method)
-
-waitbarHandle = waitbar(0,'','CreateCancelBtn','assignin(''base'',''cancel'',true)');
-cancel = false;
-
-for ii = 1:numTests
-    waitbar((ii-1)/numTests,waitbarHandle,['Test ' num2str(ii) ' of ' num2str(numTests)]);
-    b(ii).calculate
-    if cancel
-        if ii < numTests %#ok<UNRCH>
-            warning('Test unterbrochen nach Durchgang %i',ii)
-        end
-        break
-    end
-end
-delete(waitbarHandle)
-
-C_el = cell2mat(b.get('Z.C_el')');
-R_el = cell2mat(b.get('Z.R_el')');
-
-assert(all([isreal(C_el) isreal(R_el)],'all'),'Imaginäre Werte vorhanden')
-assert(all([C_el<1e-8    C_el>1e-11  ],'all'),'Ungewöhnliche Werte für C_el')
-assert(all([R_el<1e9     R_el>1e7    ],'all'),'Ungewöhnliche Werte für R_el')
-
-msgbox('Test erfolgreich')
\ No newline at end of file
diff --git a/testProgram.m b/testProgram.m
new file mode 100644
index 0000000000000000000000000000000000000000..46a73ca59b43cfaa6c0b302c56e594ad86bb7dbf
--- /dev/null
+++ b/testProgram.m
@@ -0,0 +1,389 @@
+% Skript, um die Funktionsfähigkeit von BearImp zu testen. Muss vor jedem
+% Merge ausgeführt werden. Neue Berechnungsmethoden etc. müssen hier
+% ergänzt werden.
+
+%    pmd Berechnungstool Lagerimpedanz
+%    Autor: Julius van der Kuip
+
+clear all
+
+%% 1. Testlauf 'default' (a) [Standartfall Stahllager mit Voreinstellungen]
+
+% Berechnung
+a = BearImp('default');
+a.calculate
+
+% Plotvariablen definieren
+a_delta_s_1 = halfSym(a.B.delta_s(1,:),'axSym');
+a_delta_s_2 = halfSym(a.B.delta_s(2,:),'pointSym');
+a_B_Q       = halfSym(a.B.Q, 'axSym');
+a_B_s       = halfSym(a.B.s, 'axSym');
+a_H_h_0th   = halfSym(a.H.h_0th, 'axSym');
+a_C_C_el_single = change2mid(a.C.C_el_single(:,:,1));
+a_C_R_el_single = change2mid(a.C.R_el_single(:,:,1));
+a_C_C_el(1,:,:) = change2mid(a.C.C_el_aprx);
+a_C_C_el(2,:,:) = change2mid(a.Z.C_el(30e3));
+a_C_C_el = squeeze(a_C_C_el);
+a_C_R_el(1,:,:) = change2mid(a.C.R_el_aprx);
+a_C_R_el(2,:,:) = change2mid(a.Z.R_el(30e3));
+a_C_R_el = squeeze(a_C_R_el);
+a_Z_Z_el = change2mid(a.Z.Z_el(30e3));
+
+% Plotten
+figure('name', ['Test1_Stahl_default: ' a.L.Labelling ' bei F_r=' num2str(a.F_r), ' N und omega= ' num2str(a.omega), ' 1/s und T_Oil= ' num2str(a.T_Oil) ' °C'])
+
+subplot(3,4,1)
+plotBearImp(a_delta_s_1*10^6, 'Verschiebung der Welle in y-Richtung', "\delta_y in \mu m")
+
+subplot(3,4,2)
+plotBearImp(a_delta_s_2*10^6, 'Verschiebung der Welle in x-Richtung',"\delta_x in \mu m") 
+
+subplot(3,4,3)
+plotBearImp(a_B_Q*10^-3, 'Reaktionskraft der Wälzkörper',"Q in kN") 
+
+subplot(3,4,4)
+plotBearImp(a_B_s*10^6, 'Überschneidung WK mit Laufbahn',"s in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,5)
+plotBearImp(a_H_h_0th*10^6, 'Schmierfilmdicke',"h_{0,th} in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,6)
+plotBearImp(a_C_C_el_single*10^12, 'Einzelkapazität Kontaktfläche',"C_{el,single} in pF")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,7)
+plotBearImp(a_C_R_el_single*10^-3, 'Einzelwiderstand Kontaktfläche',"R_{el,single} in k\Omega")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,8)
+plotBearImp(a_C_C_el*10^12, 'Gesamtkapazität Lager',"C_{el} in pF")
+legend('approximiert','exakt')
+
+subplot(3,4,9)
+plotBearImp(a_C_R_el*10^-3, 'Gesamtwiderstand Lager',"R_{el} in k\Omega")
+legend('approximiert','exakt')
+
+subplot(3,4,10)
+plotBearImp(abs(a_Z_Z_el), 'Betrag Gesamtimpedanz Lager bei Frequenz 30 kHz',"|Z_{el}| in \Omega")
+
+subplot(3,4,11)
+plotBearImp(angle(a_Z_Z_el)/pi*180, 'Phase Gesamtimpedanz Lager bei Frequenz 30 kHz',"\Phi_{el} in °")
+
+msgText = '1. Testlauf erfolgreich abgeschlosen';
+msg = msgbox(msgText,'Programm Update');
+
+%% 2. Testlauf (b) [werden ohne obj.calculate trotzdem bei korrektem psi die richtigen Werte berechnet?']
+
+% Berechnung
+b = BearImp('default');
+
+b.calcLub
+b.calcClear
+b.calcGeo
+b.calcLoad
+b.calcFilm
+b.calcCap
+b.calcImp
+
+% Plotvariablen definieren
+b_delta_s_1 = change2mid(b.B.delta_s(1,:));
+b_delta_s_2 = change2mid(b.B.delta_s(2,:));
+b_B_Q       = change2mid(b.B.Q);
+b_B_s       = change2mid(b.B.s);
+b_H_h_0th   = change2mid(b.H.h_0th);
+b_C_C_el_single = change2mid(b.C.C_el_single(:,:,1));
+b_C_R_el_single = change2mid(b.C.R_el_single(:,:,1));
+b_C_C_el(1,:,:) = change2mid(b.C.C_el_aprx);
+b_C_C_el(2,:,:) = change2mid(b.Z.C_el(30e3));
+b_C_C_el = squeeze(b_C_C_el);
+b_C_R_el(1,:,:) = change2mid(b.C.R_el_aprx);
+b_C_R_el(2,:,:) = change2mid(b.Z.R_el(30e3));
+b_C_R_el = squeeze(b_C_R_el);
+b_Z_Z_el = change2mid(b.Z.Z_el(30e3));
+
+% Plotten
+figure('name', ['Test2_Stahl_ohneCalculate: ' b.L.Labelling ' bei F_r=' num2str(b.F_r), ' N und omega= ' num2str(b.omega), ' 1/s und T_Oil= ' num2str(b.T_Oil) ' °C'])
+
+subplot(3,4,1)
+plotBearImp(b_delta_s_1*10^6, 'Verschiebung der Welle in y-Richtung', "\delta_y in \mu m")
+
+subplot(3,4,2)
+plotBearImp(b_delta_s_2*10^6, 'Verschiebung der Welle in x-Richtung',"\delta_x in \mu m") 
+
+subplot(3,4,3)
+plotBearImp(b_B_Q*10^-3, 'Reaktionskraft der Wälzkörper',"Q in kN") 
+
+subplot(3,4,4)
+plotBearImp(b_B_s*10^6, 'Überschneidung WK mit Laufbahn',"s in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,5)
+plotBearImp(b_H_h_0th*10^6, 'Schmierfilmdicke',"h_{0,th} in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,6)
+plotBearImp(b_C_C_el_single*10^12, 'Einzelkapazität Kontaktfläche',"C_{el,single} in pF")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,7)
+plotBearImp(b_C_R_el_single*10^-3, 'Einzelwiderstand Kontaktfläche',"R_{el,single} in k\Omega")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,8)
+plotBearImp(b_C_C_el*10^12, 'Gesamtkapazität Lager',"C_{el} in pF")
+legend('approximiert','exakt')
+
+subplot(3,4,9)
+plotBearImp(b_C_R_el*10^-3, 'Gesamtwiderstand Lager',"R_{el} in k\Omega")
+legend('approximiert','exakt')
+
+subplot(3,4,10)
+plotBearImp(abs(b_Z_Z_el), 'Betrag Gesamtimpedanz Lager bei Frequenz 30 kHz',"|Z_{el}| in \Omega")
+
+subplot(3,4,11)
+plotBearImp(angle(b_Z_Z_el)/pi*180, 'Phase Gesamtimpedanz Lager bei Frequenz 30 kHz',"\Phi_{el} in °")
+
+msgText = {msgText;''; '2. Testlauf erfolgreich abgeschlosen!';'Vergleichen Sie die Verläufe des 1. Testlaufs mit den Verläufen des 2. Testlaufs.';'-> Diese sollten identisch aussehen!'};
+msgbox(msgText,'Programm Update','replace'  ); 
+
+%% 3. Testlauf (c) [Mixlager mit default-Einstellungen]
+% Berechnung
+c = BearImp('default');
+c.setBearing("6205-C3 mix")
+c.calculate
+
+assert(length(c.psi) == length(c.psi_calc), 'psi_calc entspricht nicht dem optimalen Winkelvektor! Es ist die Funktion ''find_nessecary_psi'' in caculate zu überprüfen.')
+
+% Plotvariablen definieren
+c_delta_s_1 = change2mid(c.B.delta_s(1,:));
+c_delta_s_2 = change2mid(c.B.delta_s(2,:));
+c_B_Q       = change2mid(c.B.Q);
+c_B_s       = change2mid(c.B.s);
+c_H_h_0th   = change2mid(c.H.h_0th);
+c_C_C_el_single = change2mid(c.C.C_el_single(:,:,1));
+c_C_R_el_single = change2mid(c.C.R_el_single(:,:,1));
+c_C_C_el(1,:,:) = change2mid(c.C.C_el_aprx);
+c_C_C_el(2,:,:) = change2mid(c.Z.C_el(30e3));
+c_C_C_el = squeeze(c_C_C_el);
+c_C_R_el(1,:,:) = change2mid(c.C.R_el_aprx);
+c_C_R_el(2,:,:) = change2mid(c.Z.R_el(30e3));
+c_C_R_el = squeeze(c_C_R_el);
+c_Z_Z_el = change2mid(c.Z.Z_el(30e3));
+
+% Plotten
+figure('name', ['Test3_Mix_default: ' c.L.Labelling ' bei F_r=' num2str(c.F_r), ' N und omega= ' num2str(c.omega), ' 1/s und T_Oil= ' num2str(c.T_Oil) ' °C'])
+
+subplot(3,4,1)
+plotBearImp(c_delta_s_1*10^6, 'Verschiebung der Welle in y-Richtung', "\delta_y in \mu m")
+
+subplot(3,4,2)
+plotBearImp(c_delta_s_2*10^6, 'Verschiebung der Welle in x-Richtung',"\delta_x in \mu m") 
+
+subplot(3,4,3)
+plotBearImp(c_B_Q*10^-3, 'Reaktionskraft der Wälzkörper',"Q in kN") 
+
+subplot(3,4,4)
+plotBearImp(c_B_s*10^6, 'Überschneidung WK mit Laufbahn',"s in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,5)
+plotBearImp(c_H_h_0th*10^6, 'Schmierfilmdicke',"h_{0,th} in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,6)
+plotBearImp(c_C_C_el_single*10^12, 'Einzelkapazität Kontaktfläche',"C_{el,single} in pF")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,7)
+plotBearImp(c_C_R_el_single*10^-3, 'Einzelwiderstand Kontaktfläche',"R_{el,single} in k\Omega")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,8)
+plotBearImp(c_C_C_el*10^12, 'Gesamtkapazität Lager',"C_{el} in pF")
+legend('approximiert','exakt')
+
+subplot(3,4,9)
+plotBearImp(c_C_R_el*10^-3, 'Gesamtwiderstand Lager',"R_{el} in k\Omega")
+legend('approximiert','exakt')
+
+subplot(3,4,10)
+plotBearImp(abs(c_Z_Z_el), 'Betrag Gesamtimpedanz Lager bei Frequenz 30 kHz',"|Z_{el}| in \Omega")
+
+subplot(3,4,11)
+plotBearImp(angle(c_Z_Z_el)/pi*180, 'Phase Gesamtimpedanz Lager bei Frequenz 30 kHz',"\Phi_{el} in °")
+
+msgText{6} = ''; 
+msgText{7} = '3. Testlauf erfolgreich abgeschlosen!';
+msgbox(msgText,'Programm Update', 'replace'); 
+
+%% 4. Testlauf (d) [Mix-Lager mit wenigeren Abtastpunkten]
+% Berechnung
+d = BearImp;
+
+d.F_r = 3000; %Bereich zwischen 920-1020
+d.F_a = 0;
+d.omega=2000*2*pi/60;
+d.resolutionPerRotation=10;
+d.T_Oil = 50;
+d.setBearing("6205-C3 mix")
+d.setLube("FVA Referenz-Öl III")
+
+d.calculate
+
+% Plotvariablen definieren
+d_delta_s_1 = change2mid(d.B.delta_s(1,:));
+d_delta_s_2 = change2mid(d.B.delta_s(2,:));
+d_B_Q       = change2mid(d.B.Q);
+d_B_s       = change2mid(d.B.s);
+d_H_h_0th   = change2mid(d.H.h_0th);
+d_C_C_el_single = change2mid(d.C.C_el_single(:,:,1));
+d_C_R_el_single = change2mid(d.C.R_el_single(:,:,1));
+d_C_C_el(1,:,:) = change2mid(d.C.C_el_aprx);
+d_C_C_el(2,:,:) = change2mid(d.Z.C_el(30e3));
+d_C_C_el = squeeze(d_C_C_el);
+d_C_R_el(1,:,:) = change2mid(d.C.R_el_aprx);
+d_C_R_el(2,:,:) = change2mid(d.Z.R_el(30e3));
+d_C_R_el = squeeze(d_C_R_el);
+d_Z_Z_el = change2mid(d.Z.Z_el(30e3));
+
+% Plotten
+figure('name', ['Test4_Mix_wenigerPunkte: ' d.L.Labelling ' bei F_r=' num2str(d.F_r), ' N und omega= ' num2str(d.omega), ' 1/s und T_Oil= ' num2str(d.T_Oil) ' °C'])
+
+subplot(3,4,1)
+plotBearImp(d_delta_s_1*10^6, 'Verschiebung der Welle in y-Richtung', "\delta_y in \mu m")
+
+subplot(3,4,2)
+plotBearImp(d_delta_s_2*10^6, 'Verschiebung der Welle in x-Richtung',"\delta_x in \mu m") 
+
+subplot(3,4,3)
+plotBearImp(d_B_Q*10^-3, 'Reaktionskraft der Wälzkörper',"Q in kN") 
+
+subplot(3,4,4)
+plotBearImp(d_B_s*10^6, 'Überschneidung WK mit Laufbahn',"s in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,5)
+plotBearImp(d_H_h_0th*10^6, 'Schmierfilmdicke',"h_{0,th} in \mu m")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,6)
+plotBearImp(d_C_C_el_single*10^12, 'Einzelkapazität Kontaktfläche',"C_{el,single} in pF")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,7)
+plotBearImp(d_C_R_el_single*10^-3, 'Einzelwiderstand Kontaktfläche',"R_{el,single} in k\Omega")
+legend('Innenkontakt','Außenkontakt')
+
+subplot(3,4,8)
+plotBearImp(d_C_C_el*10^12, 'Gesamtkapazität Lager',"C_{el} in pF")
+legend('approximiert','exakt')
+
+subplot(3,4,9)
+plotBearImp(d_C_R_el*10^-3, 'Gesamtwiderstand Lager',"R_{el} in k\Omega")
+legend('approximiert','exakt')
+
+subplot(3,4,10)
+plotBearImp(abs(d_Z_Z_el), 'Betrag Gesamtimpedanz Lager bei Frequenz 30 kHz',"|Z_{el}| in \Omega")
+
+subplot(3,4,11)
+plotBearImp(angle(d_Z_Z_el)/pi*180, 'Phase Gesamtimpedanz Lager bei Frequenz 30 kHz',"\Phi_{el} in °")
+
+msgText{8} = ''; 
+msgText{9} = '4. Testlauf erfolgreich abgeschlosen!';
+msgbox(msgText,'Programm Update','replace'); 
+
+%% 5. Testlauf (e,f) [werden Fehlermeldungen richtig ausgeführt?]
+
+e = BearImp('default');
+e.psi = [0, 20, 30, 90, 130, 180, 201, 200, 270, 360]/360*2*pi;
+try
+    e.calculate
+ catch ME1
+    error('Es wurden fälschlicher Weise nicht alle benötigten Werte für die Gesamtkapazitätsberechnung berechnet! (siehe ME1)')
+end
+try
+check_all_necessary_psi(e.psi,e.psi,e.L.Z,true)
+catch ME2
+    strMsg2 = 'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!';
+    assert(strcmp(strMsg2,ME2.message), 'Obwohl nicht alle Werte für Gesamtkapazität berechnet wurden, wird kein Error angezeigt! (siehe ME2)')
+end
+
+f = BearImp('default');
+f.calculate
+try
+check_all_necessary_psi(f.psi_calc,f.psi,f.L.Z,true)
+catch ME3
+    error('psi_calc ermittelt nicht alle benötigten Winkel zur Berechnung der Gesamtkapazität! (siehe ME3)')
+end
+
+try
+    f.F_r = 200;
+    f.calcImp
+    
+catch ME4
+    strMsg4 = 'Kapazitäten noch nicht berechnet';
+    assert(strcmp(strMsg4,ME4.message), 'up2date beinhaltet fehler, da nach Veränderung der radialen Kraft nicht alle benötigten Werte berechnet wurden und das Programm es trotzdem zulässt! (siehe ME4) ')
+end
+
+msgText{10} = ''; 
+msgText{11} = '5. Testlauf erfolgreich abgeschlosen!';
+msgText{12} = 'Fehlermeldungen werden richtig ausgeführt und schützen den Benutzer vor Fehlberechnungen.';
+msgText{13} = ''; 
+msgText{14} = 'Das TestProgram ist durchgeführt!';
+msgbox(msgText,'Programm Update', 'replace'); 
+
+%% Helper functions
+
+function X = change2mid(vek,Ind)
+    lv=length(vek);
+    if ~exist('Ind','var')
+        X = [vek(  :, lv/2 + 1:lv), vek(:,   1:lv / 2)]; %, vek(:,   lv/2 + 1:lv / 2 + 2)
+    else
+        X = [vek(Ind, lv/2 + 1:lv), vek(Ind, 1:lv / 2), vek(Ind, lv/2 + 1:lv / 2 + 2)];
+    end
+end
+
+function X = halfSym(vek,symCase)
+    switch symCase
+        case 'axSym'
+            X = [fliplr(vek),vek(:,2:end-1)];
+        case 'pointSym'
+            X = [-fliplr(vek),vek(:,2:end-1)];
+    end
+end
+
+function plotBearImp(vek, pltTitle, pltYlabel)
+    len = length(vek);
+    plot(1:len,vek,'linewidth',1) %length(a.B.delta_s(1,:)
+    title(pltTitle)
+    xlabel("Drehung Käfig in °")
+    ylabel(pltYlabel)
+    set(gca,'XTick',(1:floor(len/12):len+1));
+    set(gca,'XTickLabel',(-180:30:181)); 
+    xlim([0 len+1])
+end
+
+function check_all_necessary_psi(psi,obj_psi,Z,symCase)
+        
+        if symCase
+            psi_calc = [psi, 2*pi - psi];
+        else
+            psi_calc = psi;
+        end
+        
+        for run_C = 1 : length(obj_psi)
+            
+           psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj_psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
+
+            psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) = psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch %*0.9999999
+           
+            psiNeedForIndCapa = round(psiNeedForIndCapa,8);
+
+            [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-6);
+
+            %Check, if all nessecary C_el_single are calculated
+            assert(~any(ind_single_psi==0),'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!') 
+        end
+end
\ No newline at end of file