*****************************************************************
(1) c_agents_immune_oned.m
*****************************************************************
theta_1 = 2 * Du * tau / delta^2;
theta_T = 2 * DT * tau / delta^2;
nu = 2 * chiT * tau * cstar / delta^2;
Dc_ag = Dc * tau / delta^2;

T = T_pde / tau;
griglia = L / delta;
griglia_uninf = cell_uninf / delta;
n_uninf = h_uninf * delta;

U = zeros(2 * griglia + 1, T_pde + 1);
Tcell = zeros(2 * griglia + 1, T_pde + 1);
c = zeros(2 * griglia + 1, T_pde + 1);
vessel = ones(2 * griglia + 1, 1);

% Inizializzazione cellule non infette
U(:, 1) = [zeros(griglia - griglia_uninf, 1);
           n_uninf * ones(2 * griglia_uninf + 1, 1);
           zeros(griglia - griglia_uninf, 1)];

U_temp = U(:, 1);
Tc_temp = Tcell(:, 1);
c_temp = c(:, 1);

for i = 2:T + 1
    P = U_temp;  % Pressione basata solo sulle U

    % Diffusione cellule U e T con chemotassi
    U_temp = passeggiata_unod_smorzata(U_temp, theta_1);
    Tc_temp = passeggiata_unod_chemotaxisdiffusion(Tc_temp, c_temp, theta_T, nu, cstar);

    % Nascita/morte per pressione
    U_temp = nascitamorte_pressione(U_temp, P, K * delta, p * tau, 0);
    Tc_temp = nascitamorte_pressione(Tc_temp, P, K * delta, 0, qT * tau);

    % Inflow cellule T
    Tc_temp = inflow(Tc_temp, alphaz * sum(U_temp) * delta * tau, vessel);

    % Uccisione da parte del sistema immunitario
    U_temp = T_kill(U_temp, Tc_temp, min(1, zU * tau / delta));

    % Dinamica della sostanza chimica
    c_temp = dinamicac_unod(c_temp, Dc_ag, zeros(size(U_temp)), alphac * tau / delta, ...
                             U_temp, gammac * tau / delta, cstar, qc * tau);

    % Controllo NaN
    if any(isnan(U_temp)) || any(isnan(Tc_temp))
        i
    end

    % Salvataggio delle variabili ogni ora intera
    t = (i - 1) * tau;
    if t == floor(t)
        U(:, t + 1) = U_temp;
        Tcell(:, t + 1) = Tc_temp;
        c(:, t + 1) = c_temp;

        figure(1)
        clf
        subplot(1, 2, 1)
        plot(U_temp / delta)
        title("Cellule non infette")

        subplot(1, 2, 2)
        plot(Tc_temp / delta)
        title("Cellule T")
    end
end
% Cosa è stato rimosso:
% (1) variabili I, griglia_inf, n_inf
% (2) tutti i riferimenti a infezione_dued, T_kill su infette, beta
% (3) la diffusione e aggiornamento di cellule infette
% (4) le funzioni che coinvolgevano I (in dinamicac_unod, abbiamo passato uno zeros(...) per I)
% (5) plot relativi a I





*****************************************************************
(2) c_main_immune_oned.m
*****************************************************************
clearvars
c_parameters_immune_oned

L = 10;
T_pde = 1000;
maxit = 1;
delta = 0.05;
tau = 0.01;
dx = 0.05;
dt = 1e-3;
cell_uninf = Ru;
h_uninf = 0.9 * K;

space = [-L L];

chiT = chiT / 10;
zU = zU / 2;
S_z = 0.05;

tic
c_pde_immune_oned
time_pde = toc;

tic
seed = 1;
rng(seed)
c_agents_immune_oned
time_agents = toc;

U_plot = U;
Tcell_plot = Tcell;
c_plot = c;

h = figure(1); clf
vid = VideoWriter('video');
vid.FrameRate = 4;
c_video_oned_immune

figure(2)
clf
U_sum = sum(U);
T_sum = sum(Tcell);
u1_sum = sum(u1') * dx;
u3_sum = sum(u3') * dx;

plot(0:T_pde, U_sum, 'LineWidth', 2)
hold on
plot(0:T_pde, T_sum, 'LineWidth', 2)
hold on
plot(0:T_pde, u1_sum, 'k:', 'LineWidth', 2)
hold on
plot(0:T_pde, u3_sum, 'k:', 'LineWidth', 2)

axis([0 T_pde 0 max(u1_sum)])

legend('Uninfected', 'Immune', 'Location', 'northwest')







*****************************************************************
(3) c_parameters_immune_oned.m                                   
*****************************************************************
p = log(2)/37;
alpha = 3500; % inutilizzato senza infezione
Du = ((6.2 - 2.6)/2 / 40 / 24)^2 / p; % diffusione cellule non infette
Di = Du; % rimosso dal modello
K = 1e4;
Ru = 2.6; % raggio iniziale non infette in mm
% Ri = 1; % rimosso

Dc = 0.8/24;
qc = 2/24;
cstar = 5^(1/3);
alphac = 2.5 / K;
gammac = alphac / 50;

qT = 0.18 / 24;
alphaz = qT / 50;

zU = 10 / K / 10;
DT = 0.0023;
chiT = 1.65;
% Rimossi parametri legati all'infezione (beta, q, Ri, zI, ecc.)





*****************************************************************
(4) c_pde_immune_oned.m
*****************************************************************
T = T_pde / dt;
x = -L:dx:L;
Nx = 2 * L / dx;

D = Du / K;

Nx_u = 2 * Ru / dx; % punti griglia non infette

% Inizializzazione popolazioni
u1 = zeros(Nx + 1, T_pde + 1); % cellule non infette
u3 = zeros(Nx + 1, T_pde + 1); % cellule T
u4 = zeros(Nx + 1, T_pde + 1); % sostanza chimica

u = [zeros((Nx - Nx_u)/2, 1);
     h_uninf * ones(Nx_u + 1, 1);
     zeros((Nx - Nx_u)/2, 1)];

z = zeros(Nx + 1, 1);  % cellule T
c = zeros(Nx + 1, 1);  % concentrazione chimica

u1(:, 1) = u;
u3(:, 1) = z;
u4(:, 1) = c;

for k = 1:T
    rho = u;

    % DIFFUSIONE cellule non infette
    vicino_destra = [u(2:Nx+1); 0];
    vicino_sinistra = [0; u(1:Nx)];
    u = u + Du * (vicino_sinistra + vicino_destra - 2*u) * dt / dx^2;
    u(1) = u(2);
    u(Nx+1) = u(Nx);

    % DIFFUSIONE sostanza chimica
    vicino_destra = [c(2:Nx+1); 0];
    vicino_sinistra = [0; c(1:Nx)];
    c = c + Dc * (vicino_sinistra + vicino_destra - 2*c) * dt / dx^2;
    c(1) = c(2);
    c(Nx+1) = c(Nx);

    % CHEMOTASSI cellule T
    c_destra = [c(2:Nx+1); 0];
    z_destra = [z(2:Nx+1); 0];
    w = chiT * (c_destra - c) / dx;
    Fz = max(w, 0) .* z - max(-w, 0) .* z_destra - DT * (z_destra - z) / dx;
    z = z + (dt / dx) * ([0; Fz(1:Nx)] - [Fz(1:Nx); 0]);

    % REAZIONI LOCALI
    u = u + dt * (p * (1 - rho / K) .* u - zU * u .* z);
    z = z + dt * (-qT * z + alphaz * sum(u) * dx); % afflusso in funzione del danno

    % sostanza chimica prodotta da u
    c = c + dt * ((gammac * u) .* (cstar - c) - qc * c);

    % Salvataggio dati ogni 5 ore
    if mod(k, 1/dt) == 0
        u1(:, k*dt + 1) = u;
        u3(:, k*dt + 1) = z;
        u4(:, k*dt + 1) = c;

        figure(1)
        clf
        sgtitle('t=' + string(k*dt) + ' h')
        subplot(1, 2, 1)
        plot(x, u, 'LineWidth', 1)
        hold on
        plot(x, c * 200, 'LineWidth', 1)
        legend('u (non infette)', 'c (x200)')

        subplot(1, 2, 2)
        plot(x, z, 'LineWidth', 1)
        hold on
        plot(x, c * 10, 'LineWidth', 1)
        legend('T cell', 'c (x10)')
    end
end

t = 0:T_pde;

u1 = u1';
u3 = u3';
u4 = u4';

% Modifiche principali:
%(1)rimosso completamente i, u2, Nx_i, e la diffusione/termine reattivo associato
%(2)eliminato beta, q, zI, alphac*i, ecc
%(3)la produzione della sostanza chimica è ora solo da cellule non infette (u)
%(4)il salvataggio è ora in u1, u3, u4 (coerente con i plot)




*****************************************************************
(5) c_video.m
*****************************************************************
open(vid)

theta_1 = 2 * Du * tau / delta;
minAx = 0;
maxAx = K;
griglia = L / delta;

vu = sqrt(Du * p / 2);

for i = 1:10:T_pde+1
    clf

    % Cellule non infette + confronto
    plot((0:griglia)*delta, U_plot(:, i)/delta, 'LineWidth', 1)
    hold on
    plot(x, u1(i, :), '--')
    hold on
    xline(L/2 + Ru + (i - 1) * vu, '--')
    xline(L/2 - Ru - (i - 1) * vu, '--')
    axis([0 griglia * delta minAx maxAx])
    title('Cellule non infette â t = ' + string(i - 1) + ' h')
    legend('Agenti', 'PDE', 'Fronte dx', 'Fronte sx')

    frame = getframe(h);
    writeVideo(vid, frame);
end

close(vid)

% Cosa è stato rimosso:
%(1)tracce di I_plot, u2 (infette)
%(2)linee teoriche legate all'infezione (q/beta, p/beta, ecc.)





*****************************************************************
(6) c_video_oned_immune.m
*****************************************************************
open(vid)

minAx = 0;
maxAx = K;

for i = 1:10:T_pde+1
    clf
    sgtitle('t = ' + string(i - 1) + ' h')

    % Subplot: cellule non infette + confronto PDE
    subplot(1, 2, 1)
    plot(-L:delta:L, U_plot(:, i)/delta, 'LineWidth', 1)
    hold on
    plot(x, u1(i, :), 'k:', 'LineWidth', 1)
    plot(-L:delta:L, c_plot(:, i)/10, 'g--', 'LineWidth', 1)
    legend('U (Agenti)', 'u1 (PDE)', 'c/10')
    title('Cellule non infette')
    axis([0 L minAx maxAx])

    % Subplot: cellule T e chimica
    subplot(1, 2, 2)
    plot(-L:delta:L, Tcell_plot(:, i)/delta, 'LineWidth', 1)
    hold on
    plot(-L:delta:L, c_plot(:, i), 'LineWidth', 1)
    plot(x, u3(i, :), 'k:', 'LineWidth', 1)
    plot(x, u4(i, :), 'k--', 'LineWidth', 1)
    legend('Tcell (Agenti)', 'c (Agenti)', 'u3 (PDE)', 'u4 (PDE)')
    title('Cellule T e sostanza chimica')
    axis([0 L minAx maxAx])

    frame = getframe(h);
    writeVideo(vid, frame);
end

close(vid)

% Cosa è stato rimosso:
%(1)tutte le tracce di I_plot, u2 (cellule infette)
%(2)le linee di soglia virale/infezione




*****************************************************************
(7) dinamicac_unod.m
*****************************************************************
function c = dinamicac_unod(c, D, ~, alphac, U, gammac, cstar, qc)
% Aggiornamento della sostanza chimica (solo da U)
% ~: placeholder per compatibilitÃ  (era I, ora ignorato)

c_right = [c(2:end); 0];
c_left = [0; c(1:end-1)];

% Produzione chimica solo da cellule U
c = c + D * (c_left + c_right - 2 * c) + (alphac * 0 + gammac * U) .* (cstar - c) - qc * c;
c(1) = c(2);
c(end) = c(end - 1);




*****************************************************************
(8) inflow.m
*****************************************************************
function v = inflow(v,p_inflow,vessel)
%Funzione che prende in input un vettore colonna a valori interi positivi,
%una probabilitÃ  di ingresso dai vasi sanguigni e la posizione dei vasi
%stessi e restituisce il vettore ottenuto a seguito dell'ingresso. Vessel
%vale 1 nei punti del vaso e 0 altrimenti.

a=v(:);
b=vessel(:);

a(b==1)=a(b==1)+(rand(sum(b),1)<p_inflow);

v(:)=a;





*****************************************************************
(9) nascitamorte_pressione.m
*****************************************************************
function v = nascitamorte_pressione(v, P, K, p, q)
%Funzione che prende in input una matrice v a valori interi positivi, una
%matrice pressione P delle stesse dimensioni, una carrying capacity K, una 
%probabilitÃ  di nascita p e una di morte q e restituisce la matrice 
%ottenuta dopo riproduzione e morte. Se la pressione Ã¨ superiore a K, le
%cellule muoiono prima a tasso p per riequilibrare.
%NB: le nascite precedono le morti, muoiono anche i neonati


v=v+binornd(v,p*(1-P/K).*(P<K));
v=v-binornd(v,p*(P/K-1).*(P>=K)); %muoiono le eventuali cellule in eccesso
v=binornd(v,1-q);





*****************************************************************
(10) passeggiata_unod_chemotaxisdiffusion.m
*****************************************************************
function w = passeggiata_unod_chemotaxisdiffusion(v,c,theta,chi,cmax)
%Funzione che prende in input un vettore colonna a valori interi positivi 
%di lunghezza almeno 3, un vettore di concentrazione chimica della stessa 
%dimensione, la probabilitÃ  massima di moto e il valore che normalizza la 
%concentrazione e restituisce il vettore ottenuto compiendo un passo di 
%un passegiata secondo chemotassi

% c=min(c,cmax);

% P_destra=[P(2:n);0];
% P_sinistra=[0;P(1:n-1)];

% a=v(2:n-1);

%NB: per usare mnrnd Ã¨ fondamentale che i vettori siano colonne
prob_s=theta/2+chi*subplus(c(1:end-2)-c(2:end-1))/(2*cmax);
prob_d=theta/2+chi*subplus(c(3:end)-c(2:end-1))/(2*cmax);

b=mnrnd(v(2:end-1),[prob_s 1-prob_s-prob_d prob_d]);

w=[b(:,1);0;0]+[0;b(:,2);0]+[0;0;b(:,3)];

%estremi
a=binornd(v(1),theta/2+chi*subplus(c(2)-c(1))/(2*cmax));
w(1)=w(1)+v(1)-a;
w(2)=w(2)+a;

a=binornd(v(end),theta/2+chi*subplus(c(end-1)-c(end))/(2*cmax));
w(end)=w(end)+v(end)-a;
w(end-1)=w(end-1)+a;





*****************************************************************
(11) passeggiata_unod_smorzata.m
*****************************************************************
function w = passeggiata_unod_smorzata(v,theta)
%Funzione che prende in input un vettore colonna a valori interi positivi 
%di lunghezza almeno 3 e restituisce il vettore
%ottenuto compiendo un passo di una passegiata aleatoria simmetrica con
%probabilitÃ  di restare sul posto pari a 1-theta.

n=length(v);

% P_destra=[P(2:n);0];
% P_sinistra=[0;P(1:n-1)];

% a=v(2:n-1);

%NB: per usare mnrnd Ã¨ fondamentale che i vettori siano colonne
b=mnrnd(v(2:n-1),[theta/2 1-theta theta/2]);

w=[b(:,1);0;0]+[0;b(:,2);0]+[0;0;b(:,3)];

%estremi
a=binornd(v(1),theta/2);
w(1)=w(1)+v(1)-a;
w(2)=w(2)+a;

a=binornd(v(n),theta/2);
w(n)=w(n)+v(n)-a;
w(n-1)=w(n-1)+a;




*****************************************************************
(12) T_kill.m
*****************************************************************
function v = T_kill(v, T_cell, z)
% Uccisione stocastica delle cellule (non infette) da parte delle T

prob_max = 1;
v = binornd(v, 1 - min(prob_max, z * T_cell));

% Questa funzione Ã¨ giÃ  generica e funziona anche se applicata solo a cellule U








