//Create by Wubo
//Revise on 08/25/06 03:32

#include "d:\scilab-4.0\routines\stack-c.h"
void sirs __PARAMS((int *a, int n, int m, int *b, int *nodenum,int *step,int *n0, int *tr, int *ti, double *PB,double *inf));

int intsirs(fname) 
char *fname;
{
    static int l1, m1, n1, l2, m2, n2, l3,m3,n3,l4,l5,l6,l7,l8,l9;
    static int minlhs=0, maxlhs=1, minrhs=8, maxrhs=8;   //why the minlhs must equal zero??

    /* Check number of inputs (rhs=2) and outputs (lhs=1) */
    CheckRhs(minrhs,maxrhs) ;
    CheckLhs(minlhs,maxlhs) ;
    //m2=10;
    GetRhsVar(1, "i", &m1, &n1, &l1);   //the following sentence must be wrote together!!!! 
    GetRhsVar(2, "i", &m1, &n1, &l2);  //*b 
    GetRhsVar(3, "i", &m1, &m1, &l3);  //*nodenum
    GetRhsVar(4, "i", &m1, &m1, &l4);  //*step
    GetRhsVar(5, "i", &m1, &m1, &l5);  //*n0;
    GetRhsVar(6, "i", &m1, &m1, &l6);  //*tr
    GetRhsVar(7, "i", &m1, &m1, &l7);  //*ti
    GetRhsVar(8, "d", &m1, &m1, &l8);  //*PB
    m2=*istk(l4);
    CreateVar(9, "d", &m1, &m2, &l9);  
    sirs(istk(l1), m1, n1, istk(l2),istk(l3),istk(l4),istk(l5),istk(l6),istk(l7),stk(l8),stk(l9));
    /*  Return C (3)  */
    LhsVar(1) = 9;
    return 0;
}


#include  <math.h>
#include  <stdio.h>
#include  <stdlib.h>
#include  <time.h>

int N;
int T;
int N0;
int TR;
int TI;
double P;
const double PI=3.1415926;


int      **node;
int      *state;        //state of the node {S==0,I==1,R==2}
int      *temp;
int      *nodetime;     //account the time of each node
double   *inf;          //fraction of I 
double   *sus;
double   *avedeg;       //ave deg of infected nodes
double   *rou;        
int      *deg; 



void input(int *head,int *tail,int edge,int N)
{
    int aa,bb;
    int i,j;
    int e=0;
    for (i=0;i<N;i++) {
        deg[i]=0;
        for (j=0;j<2000;j++) {
            node[i][j]=0;
        }
    }
    while (e<edge) {
        aa=head[e]-1;
        bb=tail[e]-1;
        node[aa][deg[aa]]=bb;
        node[bb][deg[bb]]=aa;
        deg[aa]++;
        deg[bb]++;
        e++;
    }
}


void initial(int T,int N)
{
    int i;
    int num;
    int temp2;
    for (i=0;i<T;i++) {
        inf[i]=0;
        sus[i]=0;
        avedeg[i]=0;
    }
    for (i=0;i<N;i++) {
        state[i]=temp[i]=0;
        nodetime[i]=0;
    }
    num=0;
    while (num<N0) {
        temp2=rand()%N;
        if (state[temp2]==0) {
            state[temp2]=temp[temp2]=1;
            num++;
        }
    }
}

void infection(int t,int N)
{
    int i,j;
    int num;
    double frand;
    for (i=0;i<N;i++) {
        state[i]=temp[i];
    }
    for (i=0;i<N;i++) {
        if (state[i]==0) {
            num=0;
            for (j=0;j<deg[i];j++) {
                if (state[node[i][j]]==1) {
                    num++;
                }
            }
            frand=(double)rand()/RAND_MAX;
            if (frand<(P*num/(double)deg[i])) {
                temp[i]=1;
            }
        } else {
            if (state[i]==1) {
                if (nodetime[i]==TI)
                    temp[i]=2;
                nodetime[i]++;
            } else {
                if (nodetime[i]==(TI+TR)) {
                    temp[i]=0;
                    nodetime[i]=0;
                } else {
                    nodetime[i]++;
                }
            }
        }
    }
}

void  tongji(int av,int t,int N)
{

    int i;
    double a,b;
    a=b=0;
    for (i=0;i<N;i++) {
        if (temp[i]==1) {
            inf[t]++;
            avedeg[t]+=deg[i];
        }
        if (temp[i]==0) {
            sus[t]++;
        }
        if (nodetime[i]!=0) {
            a+=sin(2*PI*(nodetime[i]-1)/(TI+TR));
            b+=cos(2*PI*(nodetime[i]-1)/(TI+TR));
        }
    }
    a=(double)a/N;
    b=(double)b/N;
    rou[t]=pow(a*a+b*b,0.5);
    avedeg[t]=(double)avedeg[t]/inf[t];
    inf[t]=(double)inf[t]/N;
    if (t%200==0) {
        sciprint("%i %e %i %e\n",av,P,t,inf[t]);
    }

    //      sciprint("%e\n",inf[t]);

}




//void sirs __PARAMS((int *a, int n, int m, int *b, int *node,int *step,int *n0, int *tr, int *ti, double *PB,double *inf));

//#define N 10000
//#define T 1000



void sirs(a,n,m,b,nodenum,step,n0,tr,ti,PB,infm)
int a[],b[],nodenum[],step[],n0[],tr[],ti[];
double infm[],PB[];
int n,m;
{
    int i,k,j;
    double kk;

    N=*(nodenum);
    T=*(step);
    N0=*(n0);
    TR=*(tr);
    TI=*(ti);
    P=*(PB);

    node=malloc(N*sizeof(int *));
    for (i=0;i<N;i++)
        node[i]=malloc(2000*sizeof(int));   //generate node[N][2000]

    state=(int *) malloc(N * sizeof(int));
    temp=(int *)malloc(N*sizeof(int));
    nodetime=(int *)malloc(N*sizeof(int));
    inf=(double *)malloc(T*sizeof(double));
    sus=(double *)malloc(T*sizeof(double));
    avedeg=(double *)malloc(T*sizeof(double));
    rou=(double *)malloc(T*sizeof(double));
    deg=(int *)malloc(N*sizeof(int));

    srand((unsigned)time(NULL));
    input(a,b,m,N);
    for (j=0;j<1;j++) {
        initial(T,N);
        for (i=0;i<T;i++) {
            infection(i,N);
            tongji(j,i,N);
        }
        kk=0;
        for (k=(T-10000);k<T;k++) {
            kk+=(rou[k]/10000);
        }
    }

    for (i=1;i<T;i++) {
        infm[i]=inf[i];
    }

    free(node);
    free(state);
    free(temp);
    free(nodetime);
    free(inf);
    free(sus);
    free(avedeg);
    free(rou);
    free(deg);

}