load('Voc.mat');
load('power_udds.mat');
socindex=0:0.001:1;
current_index=1:1:1370;
power_index=0:1:1369;
Nc=1000;
power=power_udds./Nc;
ocv=@(t,k)interp1(socindex,Voc,(1+k));
pwr=@(t)interp1(power_index,power,t);
Ro= 0.0048*exp(31.0494/(25+15.3253));
R1=@(t,k)((7.1135e-4) + (-4.3865e-4.*(1+k)) + (2.3788e-4.*(1+k).^2))*exp(347.4707/(25+79.5819));
C1=@(t,k)(335.4518) + (3.1712e3.*(1+k)) + (-1.3214e3.*(1+k).^2) + (53.2138 + (-65.4786.*(1+k)) +(44.3761.*(1+k).^2)).*25;
R2=@(t,k)(0.0288 + (-0.073).*(1+k) + 0.0605.*(1+k).^2)*exp(16.6712./25);
C2=@(t,k)(3.1887e4) + ((-1.1593e5).*(1+k)) + (1.0493e5.*(1+k).^2)+ (60.3114 + 1.0175e4.*(1+k) + (-9.5924e3.*(1+k).^2).*25);
vf=@(t,k,m,n) ocv(t,k)-m-n;
crrnt=@(t,k,m,n) (vf(t,k,m,n)-sqrt(vf(t,k,m,n).^2-(4*Ro).*pwr(t)))/2*Ro;
tspan=[1 1370];
x0=[0;0;0];
[t,y]=ode45(@(t,y)soc(t,y,Q,crrnt,R1,R2,C1,C2),tspan,x0);
v= @(t,k,m,n) ocv(t,k)- Ro.*crrnt(t,k,m,n)-m-n;
function dydt=soc(t,y,Q,crrnt,R1,R2,C1,C2)
ssoc=1+y(1);
dydt(1,1)= crrnt(t,y(1),y(2),y(3))./(-3600*Q);
dydt(2,1)=(-1/(R1(t,ssoc).*C1(t,ssoc))).*y(2) + (1./C1(t,ssoc)).*crrnt(t,y(1),y(2),y(3));
dydt(3,1)=(-1/(R2(t,ssoc).*C2(t,ssoc))).*y(3) + (1./C2(t,ssoc)).*crrnt(t,y(1),y(2),y(3));
end
