clear all;
close all;

rIn = 5; 
rOut = 100;
t1 = [2.0,4.0,8.0,12.0,20.0];
c1 = [0.0243291,0.0368201,0.0465258,0.0480842,0.0499363];
t2 = [1.0,3.0,5.0,10.0,15.0];
c2 = [0.0141734,0.031606,0.0405562,0.0482163,0.0496632];

d1 = log(1-rOut*c1/rIn)/rOut;
d2 = log(1-rOut*c2/rIn)/rOut;

% least square approximation for the first data set
c = polyfit(-d1,t1,1)
V = c(1);
fprintf('The volume of a heart V = %10.5f.\n',V);
E = sum((t1 - c(1)*(-d1) - c(2)).^2);
fprintf('The total square error E = %10.5f.\n',E);
n = length(d1);
dApr = linspace(d1(1),d1(n),100);
tApr = polyval(c,-dApr);
figure(1);
plot(d1,t1,'g*',dApr,tApr,'b');
xlabel('log(1-r_{out}*c/r_{in})/r_{out}');
ylabel('t');
figure(2);
m = length(t1);
tApr = linspace(t1(1),t1(m),100);
cApr = rIn*(1-exp(-rOut*tApr/V))/rOut; 
plot(t1,c1,'*g',tApr,cApr,'b');
xlabel('t');
ylabel('c');

% least square approximation for the second data set
c = polyfit(-d2,t2,1)
V = c(1);
fprintf('The volume of a heart V = %10.5f.\n',V);
E = sum((t2 - c(1)*(-d2) - c(2)).^2);
fprintf('The total square error E = %10.5f.\n',E);
n = length(d2);
dApr = linspace(d2(1),d2(n),100);
tApr = polyval(c,-dApr);
figure(3);
plot(d2,t2,'g*',dApr,tApr,'b');
xlabel('log(1-r_{out}*c/r_{in})/r_{out}');
ylabel('t');
figure(4);
m = length(t2);
tApr = linspace(t2(1),t2(m),100);
cApr = rIn*(1-exp(-rOut*tApr/V))/rOut; 
plot(t2,c2,'*g',tApr,cApr,'b');
xlabel('t');
ylabel('c');
