% Created by Nick Fidalgo, Chris Gatesman, Phillip Less, Brett Peters, Samantha Sandwick, and Alex Shetler
%% Graphical Solution - calls the functions listed below to provide a graphical representation of Model 1
%% Solutions and plot

format long;    
beta = 0.00237254167;
sigma = 5;
xi = 1.8741*10^-4;
phi = 0.35;
gamma = 483;
B_0 = 174;
D_0 = 0;
E_0 = 45;
% A(t) parameters and dosages
P = [0 .3 .6 1.2 2.4 3.6 4.8; 0 11.29 14.29 22.44 69.44 76.54 110.66; 10^-10 .82 .57 .58 .2 .3 .17];
[m,n] = size(P);
tspan = linspace(0,21);

for i = 1:n
    %B(t)
    p = P(:,i);
    dose = p(1);
    A_0 = p(2);
    alpha = p(3);
    t = [0:21];
    plot(t, eq(alpha, A_0, beta, sigma, xi, phi, gamma, t, B_0, E_0),'Linewidth',1.3);
    hold on
end

title('$$\tilde{B(t)}$$: Healthy Heptatocytes', 'FontSize', 18, 'Interpreter','Latex')
legend('0 mg/kg','0.3 mg/kg','0.6 mg/kg','1.2 mg/kg','2.4 mg/kg','3.6 mg/kg','4.8 mg/kg', 'FontSize', 15,'Interpreter','Latex')
xlabel('Days', 'FontSize',12,'Interpreter','Latex')
figure

for i = 1:n
    %D(t)
    p = P(:,i);
    dose = p(1);
    A_0 = p(2);
    alpha = p(3);
    t = [0:21];
    plot(t, ew(alpha, A_0, beta, sigma, xi, phi, gamma, t, B_0, E_0),'Linewidth',1.3);
    hold on
end

title('$$\tilde{D(t)}$$: Damaged Hepatocytes', 'FontSize', 18,'Interpreter','Latex')
legend('0 mg/kg','0.3 mg/kg','0.6 mg/kg','1.2 mg/kg','2.4 mg/kg','3.6 mg/kg','4.8 mg/kg', 'FontSize', 15,'Interpreter','Latex')
xlabel('Days', 'FontSize',12)
figure

for i = 1:n
    %E(t)
    p = P(:,i);
    dose = p(1);
    A_0 = p(2);
    alpha = p(3);
    t = [0:21];
    plot(t, e(alpha, A_0, beta, sigma, xi, phi, gamma, t, B_0, E_0),'Linewidth',1.3);
    hold on
end

title('$$\tilde{E(t)}$$: Molecules of ALT in the Serum', 'FontSize', 18,'Interpreter','Latex')
legend('0 mg/kg','0.3 mg/kg','0.6 mg/kg','1.2 mg/kg','2.4 mg/kg','3.6 mg/kg','4.8 mg/kg', 'FontSize', 15, 'Interpreter','Latex')
xlabel('Days', 'FontSize',12, 'Interpreter','Latex')

%This is B(t)
function etq = eq(alpha, A_0, beta, sigma, xi, phi, gamma, t, B_0, E_0)
    t = t/sigma;
    beta_h = (beta*A_0)/(sigma*B_0);
    gamma_h = (gamma*B_0)/E_0;
    %xi_h = (xi*B_0)/(sigma*E_0);
    %D_h = D_0/B_0;
    phi_h = phi/sigma;
    xigamma_h = xi*gamma*B_0/(sigma*E_0);
    alpha_h = alpha/sigma;
    
    if (beta_h > (alpha/sigma))
        disp("beta_h > alpha/sigma");
    end    
    d1 = exp((-alpha/sigma)*t);

    etq = ((beta_h/alpha_h)*(d1-1))+1;
end

%This is D(t)
function etw = ew(alpha, A_0, beta, sigma, xi, phi, gamma, t, B_0, E_0)
    t = t/sigma;
    beta_h = (beta*A_0)/(sigma*B_0);
    gamma_h = (gamma*B_0)/E_0;
    %xi_h = (xi*B_0)/(sigma*E_0);
    %D_h = D_0/B_0;
    phi_h = phi/sigma;
    xigamma_h = xi*gamma*B_0/(sigma*E_0);
    
    c1 = ((sigma*beta_h)/(sigma-alpha));
    d1 = exp((-alpha/sigma)*t);
    c2 = ((-((sigma*beta_h)/(sigma-alpha))));
    d2 = exp(-t);
   
    etw = c1*d1 + c2*d2;
    
    if (beta_h > (alpha/sigma))
        disp("beta_h > alpha/sigma");
    end    
end

%This is for E(t) 
function et = e(alpha, A_0, beta, sigma, xi, phi, gamma, t, B_0, E_0)
    t = t/sigma;
    beta_h = (beta*A_0)/(sigma*B_0);
    gamma_h = (gamma*B_0)/E_0;
    xi_h = (xi*B_0)/(sigma*E_0);
    %D_h = D_0/B_0;
    phi_h = phi/sigma;
    xigamma_h = xi*gamma*B_0/(sigma*E_0);
    
    c1 = ((gamma_h*sigma*beta_h)/(sigma-alpha)) + ((xigamma_h*sigma*beta_h)/alpha);
    d1 = exp((-alpha/sigma)*t);
    c2 = ((gamma_h*sigma*beta_h/(sigma-alpha)));
    d2 = exp(-t);
    c3 = ((xigamma_h*sigma*beta_h)/alpha) - xigamma_h;
    
    c = 1 - (c1/(phi_h-(alpha/sigma))) + c2/(phi_h-1) + c3/(phi_h);
    et = (c1/(phi_h-(alpha/sigma)))*d1 - (c2/(phi_h-1))*d2 - c3/(phi_h) + c*exp(-phi_h*t);
    
         

    if (beta_h > (alpha/sigma))
        disp("beta_h > alpha/sigma");
    end    
    
end