%% =========================================================================
%  RITRATTI DI FASE 3D (u,i,z) per il modello ODE di viroterapia oncolitica
%  - Tre equilibri considerati:
%      E1: assenza tumore
%      E2: solo tumore non infetto
%      E3: infezione attiva 
%
%  - Vettori e traiettorie 3D; proiezioni 2D; piani nulcline; classificazione stabilita' 
%  - Solver adattivo (ode45/ode15s) con NonNegative e controllo blow-up
%  
%% =========================================================================

clear; close all; clc;

%% -------------------- PARAMETRI DI BASE (personalizzare per la tesi) ----
par.p     = 1.87e-2;    % proliferazione tumore (1/h)
par.q     = 8.34e-3;    % clearance cellule infette (1/h)
par.q_z   = 7.50e-3;    % decadimento cellule immunitarie (1/h)
par.S_z   = 5.00e-2;    % sorgente baseline immunitaria (cell/h)
par.K     = 1.00e4;     % capacita'  portante (cell)
par.zeta  = 5.00e-1;    % forza citotossica immunitaria (1/h)
par.alpha = 5.00e-2;    % reclutamento immunitario indotto da i (1/h)
par.beta  = 3.00e-2;    % efficienza d'infezione virale (1/h)

% Tolleranze e limiti di sicurezza
tol.zero    = 1e-12;
tol.exists  = 1e-12;
MAXPOP      = 2e2*par.K;   % soglia massima assoluta per evento di arresto

% Controllo sanita'  parametri
check_param_consistency(par, tol);

%% -------------------- Calcolo equilibri e classificazione ----------------
[E1, ok1, info1] = equilibrium_E1(par, tol);
[E2, ok2, info2] = equilibrium_E2(par, tol);
[E3, ok3, info3] = equilibrium_E3(par, tol);  

% Classificazione stabilita'  (nodo/fuoco/sella)
if ok1, [~, lab1, col1] = classify_equilibrium(E1, par);
else,   lab1 = info1; col1=[.5 .5 .5]; end
if ok2, [~, lab2, col2] = classify_equilibrium(E2, par);
else,   lab2 = info2; col2=[.5 .5 .5]; end
if ok3, [~, lab3, col3] = classify_equilibrium(E3, par);
else,   lab3 = info3; col3=[.5 .5 .5]; end

fprintf('\n================= EQUILIBRI E CLASSIFICAZIONE =================\n');
disp_equilibrium("E1 (assenza di tumore)", E1, ok1, lab1);
disp_equilibrium("E2 (solo tumore non infetto)", E2, ok2, lab2);
disp_equilibrium("E3 (infezione attiva)", E3, ok3, lab3);
fprintf('================================================================\n\n');

%% -------------------- Dominio e impostazioni grafiche --------------------
% Dominio automatico (puo' essere raffinato per la tesi)
z_ref = max([eps, ok1*E1.z, ok2*E2.z, ok3*E3.z]);
uMax  = 1.05*par.K;               % u in [0, K]
iMax  = 1.05*par.K;               % i in [0, K] (scala comparabile)
zMax  = max(1.25*z_ref, 0.25*par.K*par.zeta); % z dinamico
bounds.u = [0, uMax];
bounds.i = [0, iMax];
bounds.z = [0, zMax];

% Configurazione ritratti
cfg.nSeedsEq     = 12;  % semi casuali attorno ad ogni equilibrio
cfg.nSeedsBox    = 24;  % semi nel volume (uniformi/quasi)
cfg.tFinalBase   = 400; % orizzonte base (h) per traiettorie generiche
cfg.plotPlanes   = true;
cfg.plotVectorField = true;
cfg.gridVF       = 9;   % risoluzione quiver3 (per lato)
cfg.normalizeQuiver = true;
cfg.useParallel  = license('test','Distrib_Computing_Toolbox') && ~isempty(ver('parallel'));

%% -------------------- Prepara figure e piani nulcline --------------------
fig3d = figure('Color','w','Name','Ritratto di fase 3D (u,i,z)');
hold on; grid on; box on;
xlabel('u  [cell]'); ylabel('i  [cell]'); zlabel('z  [cell]');
title('Ritratto di fase 3D - traiettorie, campi, nulcline','FontWeight','bold');

% Disegno piani nulcline (opzionale)
if cfg.plotPlanes
    draw_nullcline_planes(bounds, par);
end

% Marca gli equilibri con colore in base alla stabilita'
if ok1, plot3(E1.u,E1.i,E1.z,'o','MarkerFaceColor',col1,'MarkerEdgeColor','k','MarkerSize',8,'DisplayName',['E1: ' char(lab1)]); end
if ok2, plot3(E2.u,E2.i,E2.z,'s','MarkerFaceColor',col2,'MarkerEdgeColor','k','MarkerSize',8,'DisplayName',['E2: ' char(lab2)]); end
if ok3, plot3(E3.u,E3.i,E3.z,'^','MarkerFaceColor',col3,'MarkerEdgeColor','k','MarkerSize',8,'DisplayName',['E3: ' char(lab3)]); end
legend('Location','northeastoutside');

% Campo vettoriale 3D (quiver3) su griglia regolare
if cfg.plotVectorField
    draw_vector_field(bounds, par, cfg.gridVF, cfg.normalizeQuiver);
end

view(40,25); axis tight;
xlim(bounds.u); ylim(bounds.i); zlim(bounds.z);

%% -------------------- Generazione dei semi (condizioni iniziali) --------
seeds = [];

% Semi attorno a ciascun equilibrio (se esiste)
if ok1
    seeds = [seeds, seed_around_equilibrium(E1, bounds, cfg.nSeedsEq, [0.10, 0.10, 0.00])];
end
if ok2
    seeds = [seeds, seed_around_equilibrium(E2, bounds, cfg.nSeedsEq, [0.08, 0.08, 0.00], 1e-4)];
end
if ok3
    seeds = [seeds, seed_around_equilibrium(E3, bounds, cfg.nSeedsEq, [0.08, -0.06, 0.05])];
end

% Semi nel volume (copertura globale)
seeds = [seeds, seed_in_box(bounds, cfg.nSeedsBox)];

% Rimuove duplicati/negativi e clippa nel dominio
seeds = sanitize_seeds(seeds, bounds);

%% -------------------- Integrazione delle traiettorie ---------------------
fprintf('Integrazione di %d traiettorie...\n', size(seeds,2));
if cfg.useParallel
    parfor k = 1:size(seeds,2)
        integrate_and_plot(seeds(:,k), par, bounds, MAXPOP, cfg.tFinalBase);
    end
else
    for k = 1:size(seeds,2)
        integrate_and_plot(seeds(:,k), par, bounds, MAXPOP, cfg.tFinalBase);
    end
end
drawnow;

%% -------------------- Proiezioni 2D (u-i), (u-z), (i-z) -----------------
make_2d_projections(par, bounds, E1, ok1, E2, ok2, E3, ok3);

%% ============================ FUNZIONI LOCALI ============================

function check_param_consistency(par, tol)
    fields = {'p','q','q_z','S_z','K','zeta','alpha','beta'};
    for k=1:numel(fields)
        v = par.(fields{k});
        if ~(isnumeric(v)&&isscalar(v)&&isfinite(v))
            error('Parametro "%s" non valido (deve essere scalare finito).', fields{k});
        end
        if v < -tol.zero
            error('Parametro "%s" negativo: controllare i dati.', fields{k});
        end
    end
    if par.K <= tol.zero
        error('La capacita  portante K deve essere positiva.');
    end
    if par.q_z <= tol.zero
        warning('q_z nullo o molto piccolo: attenzione alle divisioni in z*.');
    end
end

function [E1, ok, txt] = equilibrium_E1(par, tol)
    if par.q_z <= tol.zero
        ok=false; E1=struct('u',NaN,'i',NaN,'z',NaN); txt="q_z nullo: z* non definito"; return;
    end
    E1.u=0; E1.i=0; E1.z=par.S_z/par.q_z;
    ok = (E1.z > tol.exists); txt = ternary(ok,"stabile/instabile a seconda di lambda","non fisico (z*<=0)");
end

function [E2, ok, txt] = equilibrium_E2(par, tol)
    if par.q_z<=tol.zero || par.p<=tol.zero
        ok=false; E2=struct('u',NaN,'i',NaN,'z',NaN); txt="p o q_z non validi"; return;
    end
    E2.u = par.K - (par.zeta*par.S_z)/(par.p*par.q_z);
    E2.i = 0;
    E2.z = par.S_z/par.q_z;
    ok = (E2.u > tol.exists) && (E2.z > tol.exists);
    txt = ternary(ok,"esiste (biologico)","non fisico (u*<=0 o z*<=0)");
end

function [E3, ok, txt] = equilibrium_E3(par, tol)

    E3 = struct('u',NaN,'i',NaN,'z',NaN);
    denom = (par.beta + par.p) * (par.alpha*par.zeta + par.beta*par.q_z);
    if denom <= tol.zero
        ok=false; txt="denominatore nullo o non positivo in i*"; return;
    end
    num = par.K*par.p*par.q_z*(par.beta - par.q) - par.zeta*par.S_z*(par.beta + par.p);
    i3 = num / denom;
    if par.q_z <= tol.zero
        ok=false; txt="q_z nullo: z* non definito"; return;
    end
    z3 = par.alpha*i3/par.q_z + par.S_z/par.q_z;
    if par.beta <= tol.zero
        ok=false; txt="beta nullo o troppo piccolo: u* non definito"; return;
    end
    u3 = par.K*par.q/par.beta + (par.zeta/par.beta)*z3;
    ok = (i3 > tol.exists) && (u3 > tol.exists) && (z3 > tol.exists);
    if ok, E3.u=u3; E3.i=i3; E3.z=z3; txt="esiste (biologico)";
    else, txt="almeno una componente <=0: E3 non fisico"; end
end

function [isStable, label, col] = classify_equilibrium(E, par)
% Classifica tramite autovalori della Jacobiana in E
    J = jacobian_state([E.u;E.i;E.z], par);
    ev = eig(J); re = real(ev); imv = imag(ev);
    if all(re < -1e-10)
        if any(abs(imv) > 1e-9), label="fuoco stabile"; col=[0.13 0.55 0.13];
        else,                      label="nodo stabile";  col=[0.80 0.00 0.00]*0 + [0.9 0.2 0.2]; end %#ok
        isStable = true; return;
    end
    if any(abs(re) <= 1e-10)
        label="critico (non iperbolico)"; col=[0.93 0.69 0.13]; isStable=false; return;
    end
    if sum(re>0)==1, label="sella (indice 1)"; col=[0 0 0]; isStable=false; return; end
    if sum(re>0)==2, label="sella (indice 2)"; col=[0 0 0]; isStable=false; return; end
    if sum(re>0)==3
        if any(abs(imv)>1e-9), label="fuoco instabile"; else, label="nodo instabile"; end
        col=[0 0 0]; isStable=false; return;
    end
    label="indeterminato"; col=[0 0 0]; isStable=false;
end

function J = jacobian_state(x, par)
    u=x(1); i=x(2); z=x(3);
    J11 = par.p - (2*par.p*u)/par.K - (par.p*i)/par.K - (par.beta*i)/par.K - (par.zeta*z)/par.K;
    J12 = -(par.p/par.K)*u - (par.beta/par.K)*u;
    J13 = -(par.zeta/par.K)*u;
    J21 =  (par.beta/par.K)*i;
    J22 =  (par.beta/par.K)*u - par.q - (par.zeta/par.K)*z;
    J23 = -(par.zeta/par.K)*i;
    J31 = 0;
    J32 = par.alpha;
    J33 = -par.q_z;
    J = [J11 J12 J13; J21 J22 J23; J31 J32 J33];
end

function draw_nullcline_planes(bounds, par)
% Disegna piani nulcline (traslucidi) e piani coordinate u=0, i=0
    % Griglie per superfici
    Nu=30; Ni=30;
    [U,I] = meshgrid(linspace(bounds.u(1),bounds.u(2),Nu), ...
                     linspace(bounds.i(1),bounds.i(2),Ni));

    % Piano dot{z}=0
    if par.q_z>0
        Z3 = (par.alpha/par.q_z).*I + par.S_z/par.q_z;
        s1 = surf(U,I,clip(Z3,bounds.z)); set(s1,'FaceAlpha',0.25,'EdgeColor','none','FaceColor',[0.20 0.63 0.17]); %#ok
        text(bounds.u(2),bounds.i(1),clip((par.alpha/par.q_z)*bounds.i(1)+par.S_z/par.q_z,bounds.z),'<----{dz/dt}=0','Color',[0.20 0.63 0.17]);
    end

    % Piano dot{i}=0 (parte non banale)
    if par.zeta>0
        Z2 = (par.beta/par.zeta).*U - (par.K/par.zeta)*par.q;
        s2 = surf(U,I,clip(Z2,bounds.z)); set(s2,'FaceAlpha',0.25,'EdgeColor','none','FaceColor',[0.07 0.32 0.80]); %#ok
        text(bounds.u(2),bounds.i(2),clip((par.beta/par.zeta)*bounds.u(2)-(par.K/par.zeta)*par.q,bounds.z),'<----{di/dt}=0','Color',[0.07 0.32 0.80]);
    end
    % Piano i=0
    fill3([bounds.u(1) bounds.u(2) bounds.u(2) bounds.u(1)], ...
          [0 0 0 0], ...
          [bounds.z(1) bounds.z(1) bounds.z(2) bounds.z(2)], ...
          [0.07 0.32 0.80],'FaceAlpha',0.10,'EdgeColor','none');

    % Piano dot{u}=0
    if par.zeta>0
        Z1 = (par.p*par.K - par.p.*U - (par.p+par.beta).*I)/par.zeta;
        s3 = surf(U,I,clip(Z1,bounds.z)); set(s3,'FaceAlpha',0.25,'EdgeColor','none','FaceColor',[0.85 0.33 0.10]); %#ok
        text(bounds.u(1),bounds.i(2),clip((par.p*par.K - par.p*bounds.u(1) - (par.p+par.beta)*bounds.i(2))/par.zeta,bounds.z),'<----{du/dt}=0','Color',[0.85 0.33 0.10]);
    end
    % Piano u=0
    fill3([0 0 0 0], ...
          [bounds.i(1) bounds.i(2) bounds.i(2) bounds.i(1)], ...
          [bounds.z(1) bounds.z(1) bounds.z(2) bounds.z(2)], ...
          [0.85 0.33 0.10],'FaceAlpha',0.10,'EdgeColor','none');
end

function draw_vector_field(bounds, par, N, normalize)
% Campo vettoriale 3D su griglia uniforme, con quiver3
    [U,I,Z] = meshgrid(linspace(bounds.u(1),bounds.u(2),N), ...
                       linspace(bounds.i(1),bounds.i(2),N), ...
                       linspace(bounds.z(1),bounds.z(2),N));
    F1 = par.p.*U.*(1 - (U+I)/par.K) - (par.beta/par.K).*U.*I - (par.zeta/par.K).*U.*Z;
    F2 = (par.beta/par.K).*U.*I - par.q.*I - (par.zeta/par.K).*I.*Z;
    F3 = par.alpha.*I - par.q_z.*Z + par.S_z;

    if normalize
        L = sqrt(F1.^2 + F2.^2 + F3.^2) + eps;
        F1 = F1./L; F2 = F2./L; F3 = F3./L;
    end

    step = max(1, round(N/6));
    quiver3(U(1:step:end,1:step:end,1:step:end), ...
            I(1:step:end,1:step:end,1:step:end), ...
            Z(1:step:end,1:step:end,1:step:end), ...
            F1(1:step:end,1:step:end,1:step:end), ...
            F2(1:step:end,1:step:end,1:step:end), ...
            F3(1:step:end,1:step:end,1:step:end), ...
            0.8,'Color',[0.1 0.1 0.1],'LineWidth',0.7,'AutoScale','on','AutoScaleFactor',0.9, ...
            'DisplayName','campo vettoriale (normalizzato)');
end

function S = seed_around_equilibrium(E, bounds, n, perc, extra_i)
% Semi gaussiani/percentuali attorno a un equilibrio (clippati nel dominio)
    if nargin<5, extra_i=0; end
    rng(1);  % riproducibilita'
    du = max(perc(1)*max(1,E.u), 1e-6);
    di = max(abs(perc(2))*max(1,E.i), 1e-6);
    dz = max(perc(3)*max(1,E.z), 1e-6);
    S = [E.u + du*randn(1,n);
         max(E.i + di*randn(1,n), 0) + extra_i;   % piccola semina per i (se richiesta)
         max(E.z + dz*randn(1,n), 0)];
    S = sanitize_seeds(S, bounds);
end

function S = seed_in_box(bounds, n)
% Semi nel volume (quasi-uniformi)
    rng(2);
    u = linspace(bounds.u(1),bounds.u(2),ceil(n^(1/3)));
    i = linspace(bounds.i(1),bounds.i(2),ceil(n^(1/3)));
    z = linspace(bounds.z(1),bounds.z(2),ceil(n^(1/3)));
    [U,I,Z] = ndgrid(u,i,z);
    P = [U(:)'; I(:)'; Z(:)'];
    idx = round(linspace(1,size(P,2),n));
    S = P(:,idx);
end

function S = sanitize_seeds(S, bounds)
% Clipping nel dominio e rimozione duplicati ravvicinati
    S(1,:) = min(max(S(1,:), bounds.u(1)), bounds.u(2));
    S(2,:) = min(max(S(2,:), bounds.i(1)), bounds.i(2));
    S(3,:) = min(max(S(3,:), bounds.z(1)), bounds.z(2));
    % rimozione quasi-duplicati
    if ~isempty(S)
        [~,ia] = unique(round(S',6),'rows'); S = S(:,sort(ia));
    end
end

function integrate_and_plot(x0, par, bounds, MAXPOP, tFinalBase)
% Integrazione singola traiettoria e tracciamento 3D
    J0  = jacobian_state(x0, par);
    ev0 = eig(J0);
    stiff = is_stiff(ev0);
    solver = ternary(stiff, @ode15s, @ode45);

    tspan = [0, suggest_tfinal(ev0, tFinalBase)];
    opts = odeset('RelTol',1e-8,'AbsTol',1e-10, ...
                  'NonNegative',[1 2 3], ...
                  'Events',@(t,y) stop_events(t,y,MAXPOP,bounds), ...
                  'Jacobian',@(t,y) jacobian_state(y,par));
    try
        [t,X] = solver(@(tt,yy) rhs(tt,yy,par), tspan, x0, opts);
        plot3(X(:,1),X(:,2),X(:,3),'LineWidth',1.25);
    catch ME
        warning('Integrazione fallita per x0=[%.3g,%.3g,%.3g]: %s', x0(1),x0(2),x0(3), ME.message);
    end
end

function tf = suggest_tfinal(ev, base)
% Orizzonte temporale suggerito dai tassi locali
    re = -real(ev(:)); re = re(re>1e-6);
    if isempty(re), tf = base;
    else, tf = max(base, 25/min(re)); end
    tf = min(tf, 2.0e4);
end

function stiff = is_stiff(ev)
    re = abs(real(ev)); re = re(re>1e-12);
    if isempty(re), stiff=false; return; end
    stiff = (max(re)/min(re) > 200) || (max(re) > 1);
end

function [value,isterminal,direction] = stop_events(~,y,MAXPOP,bounds)
% Eventi: (i) superamento MAXPOP, (ii) uscita dal dominio, (iii) negativita' 
    cond1 = all(y < MAXPOP);
    cond2 = (y(1)>=bounds.u(1) && y(1)<=bounds.u(2) && ...
             y(2)>=bounds.i(1) && y(2)<=bounds.i(2) && ...
             y(3)>=bounds.z(1) && y(3)<=bounds.z(2));
    cond3 = all(y >= 0);
    value = double(cond1 && cond2 && cond3); % si annulla quando violato
    isterminal = 1;
    direction  = -1;
end

function dx = rhs(~, x, par)
% Sistema ODE
    u=x(1); i=x(2); z=x(3);
    f1 = par.p*u*(1 - (u+i)/par.K) - (par.beta/par.K)*u*i - (par.zeta/par.K)*u*z;
    f2 = (par.beta/par.K)*u*i - par.q*i - (par.zeta/par.K)*i*z;
    f3 = par.alpha*i - par.q_z*z + par.S_z;
    dx = [f1; f2; f3];
end

function make_2d_projections(par, bounds, E1, ok1, E2, ok2, E3, ok3)
% Proiezioni 2D con quiver e traiettorie (integrazioni veloci) su tre piani
    planes = {'z = z_E1','z = z_E3','i = 0'};
    Zvals  = [ok1*E1.z + ~ok1*mean(bounds.z), ok3*E3.z + ~ok3*mean(bounds.z)];
    figure('Color','w','Name','Proiezioni 2D');
    tiledlayout(2,2,'TileSpacing','compact','Padding','compact');

    % (u,i) @ z circa z_E1
    nexttile; proj_quiver_ui(par, bounds, Zvals(1));
    hold on; title('(u,i) su z=z_{rif1}'); xlabel('u'); ylabel('i'); grid on; box on;

    % (u,i) @ z circa z_E3 (se disponibile, altrimenti media)
    nexttile; proj_quiver_ui(par, bounds, Zvals(2));
    hold on; title('(u,i) su z=z_{rif2}'); xlabel('u'); ylabel('i'); grid on; box on;

    % (u,z) @ i=0
    nexttile; proj_quiver_uz(par, bounds, 0);
    hold on; title('(u,z) su i=0'); xlabel('u'); ylabel('z'); grid on; box on;

    % (i,z) @ u circa K/2
    nexttile; proj_quiver_iz(par, bounds, 0.5*(bounds.u(1)+bounds.u(2)));
    hold on; title('(i,z) su u=K/2'); xlabel('i'); ylabel('z'); grid on; box on;

    % Marca proiezioni degli equilibri
    if ok1
        nexttile(1); plot(E1.u,E1.i,'ko','MarkerFaceColor',[.7 .7 .7]); 
        nexttile(2); plot(E1.u,E1.i,'ko','MarkerFaceColor',[.7 .7 .7]); 
        nexttile(3); plot(E1.u,E1.z,'ko','MarkerFaceColor',[.7 .7 .7]); 
        nexttile(4); plot(E1.i,E1.z,'ko','MarkerFaceColor',[.7 .7 .7]); 
    end
    if ok2
        nexttile(1); plot(E2.u,E2.i,'ks','MarkerFaceColor',[.5 .5 .5]); 
        nexttile(2); plot(E2.u,E2.i,'ks','MarkerFaceColor',[.5 .5 .5]); 
        nexttile(3); plot(E2.u,E2.z,'ks','MarkerFaceColor',[.5 .5 .5]); 
        nexttile(4); plot(E2.i,E2.z,'ks','MarkerFaceColor',[.5 .5 .5]); 
    end
    if ok3
        nexttile(1); plot(E3.u,E3.i,'k^','MarkerFaceColor',[.3 .3 .3]); 
        nexttile(2); plot(E3.u,E3.i,'k^','MarkerFaceColor',[.3 .3 .3]); 
        nexttile(3); plot(E3.u,E3.z,'k^','MarkerFaceColor',[.3 .3 .3]); 
        nexttile(4); plot(E3.i,E3.z,'k^','MarkerFaceColor',[.3 .3 .3]); 
    end
end

function proj_quiver_ui(par, bounds, zfix)
% Campo proiettato nel piano (u,i) a z=zfix
    Nu=24; Ni=24;
    [U,I] = meshgrid(linspace(bounds.u(1),bounds.u(2),Nu), ...
                     linspace(bounds.i(1),bounds.i(2),Ni));
    Z = zfix*ones(size(U));
    F1 = par.p.*U.*(1 - (U+I)/par.K) - (par.beta/par.K).*U.*I - (par.zeta/par.K).*U.*Z;
    F2 = (par.beta/par.K).*U.*I - par.q.*I - (par.zeta/par.K).*I.*Z;
    L  = sqrt(F1.^2 + F2.^2) + eps;
    quiver(U,I,F1./L,F2./L,0.8,'AutoScale','on','AutoScaleFactor',0.9,'Color',[0.2 0.2 0.2]);
end

function proj_quiver_uz(par, bounds, ifix)
% Campo proiettato nel piano (u,z) a i=ifix
    Nu=24; Nz=24;
    [U,Z] = meshgrid(linspace(bounds.u(1),bounds.u(2),Nu), ...
                     linspace(bounds.z(1),bounds.z(2),Nz));
    I = ifix*ones(size(U));
    F1 = par.p.*U.*(1 - (U+I)/par.K) - (par.beta/par.K).*U.*I - (par.zeta/par.K).*U.*Z;
    F3 = par.alpha.*I - par.q_z.*Z + par.S_z;
    L  = sqrt(F1.^2 + F3.^2) + eps;
    quiver(U,Z,F1./L,F3./L,0.8,'AutoScale','on','AutoScaleFactor',0.9,'Color',[0.2 0.2 0.2]);
end

function proj_quiver_iz(par, bounds, ufix)
% Campo proiettato nel piano (i,z) a u=ufix
    Ni=24; Nz=24;
    [I,Z] = meshgrid(linspace(bounds.i(1),bounds.i(2),Ni), ...
                     linspace(bounds.z(1),bounds.z(2),Nz));
    U = ufix*ones(size(I));
    F2 = (par.beta/par.K).*U.*I - par.q.*I - (par.zeta/par.K).*I.*Z;
    F3 = par.alpha.*I - par.q_z.*Z + par.S_z;
    L  = sqrt(F2.^2 + F3.^2) + eps;
    quiver(I,Z,F2./L,F3./L,0.8,'AutoScale','on','AutoScaleFactor',0.9,'Color',[0.2 0.2 0.2]);
end

function disp_equilibrium(name, E, ok, label)
    if ok
        fprintf('%s: u*=%.4g, i*=%.4g, z*=%.4g  - %s\n', name, E.u, E.i, E.z, label);
    else
        fprintf('%s: NON disponibile (%s)\n', name, label);
    end
end

function y = clip(x, rng)
% Taglia i valori fuori range per visualizzare superfici nelle box limits
    y = min(max(x, rng(1)), rng(2));
end

function val = ternary(cond,a,b)
    if cond, val=a; else, val=b; end
end
