hold off
clc
clear all

tic

load('test messung.mat')    %Laden des Test-Impedanzspektrum
Frequenz(1)=[];
Z_re(1)=[];
Z_im(1)=[];
w=Frequenz*2*pi; 
load('imp_Soc.mat')     %Laden der Impedanzspektren bei unterschiedlichem SoC
Soc=Soc.';              %Matrix transponieren
werte=['100%';'78% ';'66% ';'47% ';'23% ';'0%  '];  %gemessene SoC

x3=[36,16,0.00005,92,0.005,1.2];  
x4=[36,16,0.00009,0.75,97,0.004,0.9,1.2];
fun=@(x,w)x(1)+x(2)*(1+(1i*w*x(2)*x(3))).^-1+x(4)*(1+(1i*w*x(4)*x(5))).^-1+x(6)*(1-1i)./sqrt(w); 
fun2=@(x2,w)x2(1)+x2(2)*(1+(1i*w*x2(2)*x2(3)).^x2(4)).^-1+x2(5)*(1+(1i*w*x2(5)*x2(6)).^x2(7)).^-1+x2(8)*(1-1i)./sqrt(w);  %Impedanzberechnung mit 2*gedämpften RC-Glied + Rs und Warburg-Impedanz

Fl=Flaeche(x4,2)

n=-2:0.01:3;
w1=10.^n;
w1=w1*2*pi;

Flaeche2=abs(trapz(Soc(:,1),Soc(:,2)))

100/Flaeche2*Fl

plot(fun(x3,w1),'-.')
plot(fun2(x4,w1),'-.')
%hold on
%plot(Soc(:,1),Soc(:,2))
%set(gca,'Ydir','reverse')
disp('Fertig')

toc

function [A]=Flaeche(x,n)
fun=@(x,w)x(1)+x(2)*(1+(1i*w*x(2)*x(3))).^-1+x(4)*(1+(1i*w*x(4)*x(5))).^-1+x(6)*(1-1i)./sqrt(w);
fun2=@(x2,w)x2(1)+x2(2)*(1+(1i*w*x2(2)*x2(3)).^x2(4)).^-1+x2(5)*(1+(1i*w*x2(5)*x2(6)).^x2(7)).^-1+x2(8)*(1-1i)./sqrt(w);  %Impedanzberechnung mit 2*gedämpften RC-Glied + Rs und Warburg-Impedanz

load('imp_Soc.mat')     %Laden der Impedanzspektren bei unterschiedlichem SoC
Soc=Soc.';  
l=size(Soc);
kor_x=0;
kor_y=0;
w_s=0;
fl=0;
fl_buf=0;
A=0;
k=1;
for i=1:l(1)-1
    disp(i)
    [kor_x(end+1,1),kor_y(end+1,1),w_s(end+1,1)]=schnittpunkt(Soc(i,1),Soc(i+1,1),Soc(i,2),Soc(i+1,2),x,n);
    if(kor_x(end)==0)
        X=[Soc(i,1),Soc(i+1,1)];
        Y=[Soc(i,2),Soc(i+1,2)];
        fl=fl+trapz(X,Y);
        if(i==l(1)-1)
            if(n==1)
            fl_buf=trapz(real(fun(x,0.01*2*pi:0.01:w_sp)),imag(fun(x,0.01*2*pi:0.01:w_sp)));
            plot(real(fun(x,0.01*2*pi:0.01:w_sp)),imag(fun(x,0.01*2*pi:0.01:w_sp)))
            A=abs(abs(fl)-fl_buf)+A;
            else
            fl_buf=trapz(real(fun2(x,0.01*2*pi:0.01:w_sp)),imag(fun2(x,0.01*2*pi:0.01:w_sp)));
            plot(real(fun2(x,0.01*2*pi:0.01:w_sp)),imag(fun2(x,0.01*2*pi:0.01:w_sp)))
            A=abs(abs(fl)-fl_buf)+A;
            end
        end
    else
        X=[Soc(i,1),kor_x(end)];
        Y=[Soc(i,2),kor_y(end)];
        fl=fl+trapz(X,Y);
        hold on
        if(k==1)
            if(n==1)
        fl_buf=trapz(real(fun(x,w_s(end):0.01:1000*2*pi)),imag(fun(x,w_s(end):0.01:1000*2*pi)));
        plot(real(fun(x,w_s(end):0.01:1000*2*pi)),imag(fun(x,w_s(end):0.01:1000*2*pi)))
        w_sp=w_s(end);
            else
        fl_buf=trapz(real(fun2(x,w_s(end):0.01:1000*2*pi)),imag(fun2(x,w_s(end):0.01:1000*2*pi)));
        plot(real(fun2(x,w_s(end):0.01:1000*2*pi)),imag(fun2(x,w_s(end):0.01:1000*2*pi)))
        w_sp=w_s(end);
            end
        else
            if(n==1)
            fl_buf=trapz(real(fun(x,w_s(end):0.01:w_sp)),imag(fun(x,w_s(end):0.01:w_sp)));
            plot(real(fun(x,w_s(end):0.01:w_sp)),imag(fun(x,w_s(end):0.01:w_sp)))
            w_sp=w_s(end);
            else
            fl_buf=trapz(real(fun2(x,w_s(end):0.01:w_sp)),imag(fun2(x,w_s(end):0.01:w_sp)));
            plot(real(fun2(x,w_s(end):0.01:w_sp)),imag(fun2(x,w_s(end):0.01:w_sp)))
            w_sp=w_s(end);
            end                
        end
        A=abs(abs(fl)-fl_buf)+A;
        X=[kor_x(end),Soc(i+1,1)];
        Y=[kor_y(end),Soc(i+1,2)];
        fl=trapz(X,Y);
        k=2;
    end
end
sp=[kor_x,kor_y,w_s]

end

function [k_x,k_y,w_sp]=schnittpunkt(x1,x2,y1,y2,p,q)
fun=@(x,w)x(1)+x(2)*(1+(1i*w*x(2)*x(3))).^-1+x(4)*(1+(1i*w*x(4)*x(5))).^-1+x(6)*(1-1i)./sqrt(w);
fun2=@(x2,w)x2(1)+x2(2)*(1+(1i*w*x2(2)*x2(3)).^x2(4)).^-1+x2(5)*(1+(1i*w*x2(5)*x2(6)).^x2(7)).^-1+x2(8)*(1-1i)./sqrt(w);  %Impedanzberechnung mit 2*gedämpften RC-Glied + Rs und Warburg-Impedanz

k_x=0;
k_y=0;
w_sp=0;
m=(y2-y1)/(x2-x1);
b=((y2-y1)/(x2-x1))*(-x1)+y1;

grad=@(x)m*x+b;
n=-2:0.001:3;
w1=10.^n;
i=1;
l=0;
y=grad(x);
tol=0.01;
while(l==0) 
    k=0;
    j=1;
while(k==0)
    if(q==1)
if(abs(real(fun(p,w1(j)))-x(i))<tol)
if(abs(imag(fun(p,w1(j)))-grad(x(i)))<tol)
    k_x=x(i);
    k_y=grad(x(i));
    w_sp=w1(j);
    plot(x(i),grad(x(i)),'rx')
    k=1;
    l=1;
    hold on
end
end
    else
        if(abs(real(fun2(p,w1(j)))-x(i))<tol)
if(abs(imag(fun2(p,w1(j)))-grad(x(i)))<tol)
    k_x=x(i);
    k_y=grad(x(i));
    w_sp=w1(j);
    plot(x(i),grad(x(i)),'rx')
    k=1;
    l=1;
    hold on
end
        end
    end
j=j+1;
if(j>(length(w1)-1)) 
    k=1;
end
    end
    i=i+1;
    if(i>(length(x)-1)) 
        l=1;
    end
end
plot(x,y)
end