%% Modele SIR 
N = 1000000;    % Population totale
 beta = 0.2;     % Taux de contamination
 gamma = 0.1;    % Taux de guérison
 delta = 0.01;   % Taux de mortalité
 % 
 sird(N, beta, gamma, delta);
 %% on augmente la mortalité
 sird(N,beta,gamma,0.1);
 % Lorsqu'on augmente le paramètre $\delta$ 
 % on constate que l'épidémie s'arrête! 
 % Le virus est pourtant beaucoup plus virulent.
 % Mais les personnes infectées décèdent trop vite 
 % 
 % et n'ont pas le temps de contaminer d'autres personnes.
 %% on s'arrange pour que R_0 =1
 sird(N, 0.11, 0.1, 0.01);
 %Lorsque $R_0 = 1$ on retrouve le même phénomène, 
 % l'épidémie ne se déclenche pas. 
 % Le nombre d'infectés n'augmente pas suffisamment.

% Définition des équations différentielles du modèle SIR
    function dydt = deriv(t, y, N, beta, gamma, delta)
        % N : population totale
        % beta : taux de contamination
        % gamma : taux de guérison
        % delta : taux de mortalité

        S = y(1);
        I = y(2);
        R = y(3);
        D = y(4);

        dSdt = -beta * S * I / N;
        dIdt = beta * S * I / N - gamma * I - delta * I;
        dRdt = gamma * I;
        dDdt = delta * I;
        
        dydt = [dSdt; dIdt; dRdt; dDdt];
    end

    function SIRD = sird(N, beta, gamma, delta)
        % Nombre initial d'individus infectés, rétablis et décédés, le reste est sain
        I0 = 1;
        R0 = 0;
        D0 = 0;
        S0 = N - I0 - R0 - D0;

        % Vecteur des conditions initiales
        y0 = [S0; I0; R0; D0];

        % Vecteur des temps (2 ans)
        t = linspace(0, 730, 730);

        % Résolution des équations différentielles du modèle SIR
        sol = ode45(@(t, y) deriv(t, y, N, beta, gamma, delta), [0, 730], y0);

        SIRD = deval(sol, t);

        % Traçage des résultats
        plot(t, SIRD(1,:), 'b', 'LineWidth', 1.5); hold on;
        plot(t, SIRD(2,:), 'r', 'LineWidth', 1.5);
        plot(t, SIRD(3,:), 'g', 'LineWidth', 1.5);
        plot(t, SIRD(4,:), 'k', 'LineWidth', 1.5);
        xlabel('temps (jours)');
        ylabel('Nombre d''individus (millions)');
        
        % Calcul du R_0
        R_0 = beta * S0 / (N * (gamma + delta));
        title(['R_0 = ' num2str(R_0)]);
        
        % Affichage des valeurs finales
        text(300, 0.95e6, ['Infectés = ' num2str(SIRD(2,end))]);
        text(300, 0.85e6, ['Rétablis = ' num2str(SIRD(3,end))]);
        text(300, 0.75e6, ['Décédés = ' num2str(SIRD(4,end))]);
        
        legend('Sains', 'Infectés', 'Rétablis', 'Décédés');
        grid on;
        hold off;
        %saveas(gcf, 'sir.png'); pour sauvegarder
        disp(['Nombre maximum de personnes infectées: ', num2str(round(max(SIRD(2,:))))]);
        % Vérification de la conservation de la population
        disp(['Conservation de la population: ', num2str(max(abs(N - sum(SIRD, 1))))]);
    end