clear;clc;

load time_Hadji
load Tcellcount

rho = 0.2; %4.85e-2*1.75; %4.85e-2*1.75; %1e-2
C = 55.6;  %6e2; %3e2

t_final=100;

%Not actual parameter values in simulation - must
%copy and paste into MM_3CM.m file

% lamb = 10;
% dT = 0.01;
% k = 3.87e-3;    %4.5e-5; %4.57e-5
% dI = 0.28;       %0.4;
% NdI = 5.37e-1;   %32; %38
% dV = 2.39;       %2.9;

lamb = 10;
dT = 0.01;
k = 4.72e-3;    %4.5e-5; %4.57e-5
dI = 0.69;       %0.4;
NdI = 5.37e-1;   %32; %38
dV = 2.39;       %2.9;


R_0 = (lamb*k*NdI)/(dT*dI*dV)
R_m = rho/(k*C)

%Eu_T = lamb/dT
Ei_T= (dI*dV)/(k*NdI)

b = C*k*(1 - R_m) + dT*(1 - R_0);
c = C*dT*(1 - R_0);

Ei_V = (1/(2*k))*(-b + sqrt(b^2 - 4*k*c))

p =[lamb,k,dT,rho,C,dI,NdI,dV];

%y0 = [Ei_T; 0; Ei_V]; 
y0 = [1e3; 0; 1e-3];

f = @(t,x) MM_3CM_allp(p,t,x);  %function to call ode15s


[ts,ys] = ode15s(f, [0,t_final], y0);

%y1 = [Ei_T; 0; 1e-1]; 
%[ts1,ys1] = ode15s(f, [0,t_final], y1);

figure;
plot(time_Hadji, Tcellcount, '-r', ts, ys(:, 1),'--b', 'linewidth', 1), axis([0,t_final, 300, 1100]) 
%ts1, ys1(:, 1),'--b'), axis([0,t_final, 0, 1200])
%hold on

legend('Full Model', 'Reduced Model')
title('T-cell Population vs Time','interpreter','latex','fontweight','bold','fontsize',14)
xlabel('Time (days)','interpreter','latex','fontsize',12)
ylabel('T-cell Population (mm$^{-3}$)','interpreter','latex','fontsize',12)

% td = (0:5)*21;
% yd = [1e3, 940, 502, 600, 660, 680];
% plot(td, yd, 'ob')

%figure;
%plot(ts, log10(ys(:, 3)), 'r')
%hold on
%plot(ts1, log10(ys1(:, 3)), '--r')