# Diskretni SIR model

S[n]:=S[n-1]-β/N *S[n-1]*I[n-1];
I[n]:=I[n-1]+β/N *S[n-1]*I[n-1]-γ*I[n-1];
R[n]:=R[n-1]+γ*I[n-1];

## Modeliranje epidemije gripe prema: J. M. Mahaffy, "Discrete SIR Models: Influenza" http://jmahaffy.sdsu.edu/courses/f17/math636/beamer/sir.pdf

S[0]: 157756;
I[0]: 3;
R[0]: 0;
N: S[0]+I[0]+R[0];
β: 3.9928;
γ: 3.5170;

s: makelist([n,S[n]],n,0,48);
i:  makelist([n,I[n]],n,0,48);
r:  makelist([n,R[n]],n,0,48);

crtaj(niz):=wxplot2d ([discrete, niz],[style,points],[xlabel, "n"],[ylabel,""],[color,red]);

crtaj(i);

wxplot2d([discrete,i],[style,lines],[xlabel,"n"],[ylabel,""],[legend,"I"],[color,red]);

wxplot2d([discrete,i],[style,lines],[xlabel,"n"],[ylabel,""],[legend,"I"],[color,red],[y,0,N]);

wxplot2d([discrete,s],[style,lines],[xlabel,"n"],[ylabel,""],[legend,"S"],[color,blue],[y,0,N]);

wxplot2d([discrete,r],[style,lines],[xlabel,"n"],[ylabel,""],[legend,"R"],[color,black],[y,0,N]);

wxplot2d([[discrete,s],[discrete,i],[discrete,r]],[style,lines],[xlabel,"n"],[ylabel,""],[legend,"S","I","R"],[color,blue,red,black]);

## Modeliranje zaraze Corona virusom u RH

N: 4.5*10^6;
i0: 10;

crtaj(b,g,m):=block([], 
     kill(S,I,R), 
     S[n]:=S[n-1]-β/N *S[n-1]*I[n-1], 
     I[n]:=I[n-1]+β/N *S[n-1]*I[n-1]-γ*I[n-1],
     R[n]:=R[n-1]+γ*I[n-1], 
     S[0]: N-i0, I[0]: i0, R[0]: 0, β: b, γ: g,
     s: makelist([n,S[n]],n,0,m),
     i:  makelist([n,I[n]],n,0,m),
     r:  makelist([n,R[n]],n,0,m),
     wxplot2d([[discrete,s],[discrete,i],[discrete,r]],[style,lines],[xlabel,"n"],[ylabel,""],[legend,"S","I","R"],[color,blue,red,black],[y,0,1.2*N])
);

crtaj(0.2,1/14,200);

crtaj(0.1,1/14,800);

crtaj(0.05,1/14,5000);

crtaj(0.08,1/14,2000);

crtaj(1/14+0.001,1/14,20000);

## Modificirani model: nema imuniteta

S[n]:=S[n-1]-β/N *S[n-1]*I[n-1]+γ*(1-τ)*I[n-1];
I[n]:=I[n-1]+β/N *S[n-1]*I[n-1]-γ*I[n-1];
R[n]:=R[n-1]+γ*τ*I[n-1];

N: 4.5*10^6;
i0: 10;
τ: 0.03;

crtaj2(b,g,m):=block([], 
     kill(S,I,R), 
     S[n]:=S[n-1]-β/N *S[n-1]*I[n-1]+γ*(1-τ)*I[n-1],
     I[n]:=I[n-1]+β/N *S[n-1]*I[n-1]-γ*I[n-1],
     R[n]:=R[n-1]+γ*τ*I[n-1], 
     S[0]: N-i0, I[0]: i0, R[0]: 0, β: b, γ: g,
     s: makelist([n,S[n]],n,0,m),
     i:  makelist([n,I[n]],n,0,m),
     r:  makelist([n,R[n]],n,0,m),
     wxplot2d([[discrete,s],[discrete,i],[discrete,r]],[style,lines],[xlabel,"n"],[ylabel,""],[legend,"S","I","R"],[color,blue,red,black],[y,0,1.2*N])
);

crtaj2(0.2,1/14,500);

crtaj2(0.1,1/14,2000);

crtaj2(0.08,1/14,3000);

crtaj2(1.5,1/14,30);

crtaj2(2,1/14,30);

crtaj2(3,1/14,30);

crtaj(1.8,1/14,30);

crtaj(2.5,1/14,30);

crtaj(3,1/14,30);