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,k,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
fun3=@(w)k(1)+k(2)*(1+(1i*w*k(2)*k(3))).^-1+k(4)*(1+(1i*w*k(4)*k(5))).^-1+k(6)*(1-1i)./sqrt(w);

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

grad=@(x)m*x+b;

sp2=@(x,w)(m*x+b)-imag(fun3(w));

[x]=lsqnonlin(@sp1,[5,x1,y1],[0.01,x1,y2],[1000,x2,y1])
tol=0.1;

hold on
plot(x(2),x(3),'x')

if(abs(grad(x)-1)<tol)
     k_x=real(fun3(x(3)))
     k_y=imag(fun3(x(3)))    
     plot(k_x,k_y,'o')
end

x=x1:0.01:x2;
y=grad(x);
plot(x,y)
end

function y=sp1(x)  %Berechnung der Impedanz im ungedämpften Fall
global m
global b
global p

y=p(1)+p(2)*(1+(1i*x(3)*p(2)*p(3))).^-1+p(4)*(1+(1i*x(3)*p(4)*p(5))).^-1+p(6)*(1-1i)./sqrt(x(3));
y2=m*x(2)+b;
y=[x(1)-x(1),x(2)-real(y),y2-imag(y)];

end