clear all
clc
% %%Segers Parameters          Definitions
% Emax=1.7; %mmHg/ml          {slope of end-systolic P-V relation}
% Vd=-15; %ml                 {intercept of end-systolic P-V relation}
% LVEDV=120; %ml              {LV end diastolic volume}
% Pv=5; %mmHg                 {venous filling pressure}
% bpm=75; %bpm                 {heart rate}
% tP=0.3; %s                  {time to reach maximal elastance}
% %Windkessel Parameters      Definitions
% Z0=0.033; %mmHg/ml*s        {impedence of aorta}
% L=0.005; %mmHg/ml*s^2       {inertia of blood in arteries}
% Rall=0.6:0.12:1.8; %mmHg/ml*s  {total peripheral resistance}
% Call=0.5:0.15:2; %ml/mmHg      {total compliance}

%%
tic
LVEDV = 311;
bpm   = 75;
Vd    = 15;

Emax  = 1.7;
Emin  = 0.037;

R    = 1.15;
C    = 0.87;
Z0   = 0.033;
L    = 0.005;
Cv   = 80;
CAoP = 0.01;

T    = 60/bpm; %s

cycles = 7;

load('elastance_segers.mat') %loads in elastance from Segers paper
etime = elastance(:,1)*(75/bpm);
dt    = etime(2)-etime(1);

elast=((elastance(:,2)-elastance(1,2))*(Emax-Emin)/(1.298-elastance(1,2))+Emin);
% elast(end-99:end)=[];
% elast=[ones(100,1)*elast(1);elast];

elast2=[elast; elast; elast];
% NN=5;
% for i=1001:2000
% elast2(i)=mean(elast2(i-NN:i+NN));
% elast(i-1000,1)=elast2(i);
% end
for i=1001:2000
    delast(i-1000,1)=(-elast2(i+2)+8*elast2(i+1)-8*elast2(i-1)+elast2(i-2))/(12*dt);
end

time=linspace(0,T*cycles,1000*cycles);

options=odeset('MaxStep',1e-4);

% R_run=[0.6:0.12:1.8];
% C_run=[0.5:0.15:2];

R_run = 1.8;
C_run = 0.5;

Ea_RCsv  = zeros(length(R_run),length(C_run));
Ea_RC    = zeros(cycles, length(R_run)*length(C_run));
Ea_EQ    = Ea_RC;
Ea_SG    = Ea_RC;
R_sv     = Ea_RCsv;
C_sv     = Ea_RCsv;
Pss_RC   = Ea_RC;
R_rel    = Ea_RCsv;
C_rel    = Ea_RCsv;
% Ea_RCrel = Ea_RC;
Pes_RC   = Ea_RC;
SV_RC    = Ea_RC;

Il   = zeros(cycles * length(elastance), length(R_run)*length(C_run));
AoP  = Il;
Pa   = Il;
Pv   = Il;
LVP  = Il;

nprout = 0;

for i=1:length(R_run)
    for j=1:length(C_run)
                
        runnum = (i-1)*length(C_run) + j;
        
        prout = ['RC Iteration ', num2str(runnum), ' of ', num2str(length(R_run)*length(C_run)) ' '];

        fprintf(repmat('\b',1,nprout));
        fprintf(prout)
        nprout = numel(prout);

        R = R_run(i);
        C = C_run(j);
        
        
        initial = [0 120 120 18 elast(1)*(LVEDV-Vd)]; % [Il AoP Pa Pv LVP]
        
        z = [[bpm;zeros(999,1)], elast, delast, etime, [Z0;zeros(999,1)], [L;zeros(999,1)], ...
            [R;zeros(999,1)], [C;zeros(999,1)],[Cv;zeros(999,1)],[CAoP;zeros(999,1)]];

        [t,y] = ode45(@wind_4el_a,time,initial,options,z);

        Il(:,runnum)   = y(:,1);
        AoP(:,runnum)  = y(:,2);
        Pa(:,runnum)   = y(:,3);
        Pv(:,runnum)   = y(:,4);
        LVP(:,runnum)  = y(:,5);
%         dAoP(:,runnum) = ezdiff(AoP(:,runnum), dt);
%         dAoP(length(AoP(:,runnum)),runnum) = dAoP(length(AoP(:,runnum))-1,runnum); % tack on last point
        
        Iz = (AoP - Pa) / Z0;
        
        thinger = floor(60/bpm/dt);
        
        for k = 1:cycles
            Pes_RC(k,runnum) = mean(AoP(1000*(k-1) + 1:(1000*k),runnum));
            SV_RC(k,runnum)  = mean(Il(1000*(k-1) + 1:(1000*k),runnum) + Iz(1000*(k-1) + 1:(1000*k),runnum)) * 60/bpm;
            Ea_RC(k,runnum)  = Pes_RC(k,runnum)/SV_RC(k,runnum);
            Pss_RC(k,runnum) = Pa(1000*(k-1) + 1,runnum) - Pa(1000*k,runnum);
            Ea_EQ(k,runnum)  = (1 - C * (AoP(1000*k,runnum) - AoP(1000*(k-1) + 1,runnum)) / SV_RC(k, runnum)) * R / T;
            Ea_SG(k,runnum)  = -0.127 + 1.023 * R/T + 0.314 / C;
        end
        
        R_sv(i,j)     = R;
        C_sv(i,j)     = C;

        elastp = [];
        
        for ll = 1:cycles
            elastp = [elastp;elast];
        end
        
        Veff(:,runnum) = LVP(:,runnum)./elastp;
        LVV(:,runnum)  = Veff(:,runnum)+Vd;
        
    end
end


figure(1)
for i=1:length(C_run)
%     figure(i)
    hold on;
    for j=1:length(R_run)
               
        runnum = (i-1)*length(C_run) + j;
        plot(time,AoP(:,runnum),time,LVP(:,runnum),time,Il(:,runnum)+(AoP(:,runnum)-Pa(:,runnum))/Z0,time,Pv(:,runnum))
%         plot(time,AoP(:,runnum), time,LVP(:,runnum), time, Pv(:,runnum))
        pause(1)
    end
end

hold off;

toc