#include <stdio.h>
#include <math.h>
#define cities 200    //number or cities/airports/superorganisms
#define totaltime 90   //maximum simulated time

double S[cities], I[cities], R[cities], Saux, Iaux, Raux, population, It[cities];
int departure, arrived, i, j, time, voos[cities][cities], flag[cities];
float alpha, beta, gama, I_zero;
FILE *f1;

main()
{
    //open ANAC domestic flighs file
    f1=fopen("domestic-flighs-file.txt", "r");
    
    //main variables initialization
    alpha=0.0;                 //traveling infection parameter (modified SIR)
    beta=0.3035;           //contact infection parameter (original SIR)
    gama=0.0;                //recovering parameter (original SIR)
    
    population=1.0e6;  //initial population
     for (i=0; i<cities; i++) 
     {
         S[i]=population;     //suscectible population of city i
         I[i]=0.0;                    //infected population of citi i
         R[i]=0.0;                   //recovered population of city i
         flag[i]=0;
        for (j=0; j<cities; j++) voos[i][j]=0; //flights from citi i to j
     }
     
    //read ANAC domestic flights 
    while (!feof(f1))
   {
      fscanf (f1, "%d  %d\n", &departure, &arrived);
      if (arrived && departure) voos[arrived][departure]++;
   }
   fclose(f1);  //close ANAC file
   
    // initial infected condition (only abroad)
     I[0]=I_zero;
     
     //flights from Italy and Spain 
     voos[0][40]=196;   //Guarulhos
     voos[0][37]=60;     //Galeão
     voos[0][36]= 9;     //Fortaleza
     voos[0][82]= 9;    //Recife
     voos[0][94]=13;      //Salvador
     
    //print initial condition
    printf ("%d  ", 0);
    for (i=1; i<155; i++)
        printf ("%lf         ", 100*I[i]/population);
    printf ("\n");
     
     
     for (time=1; time<totaltime; time++) //time loop
     {
        
       //calculate travel factor It
      for (i=1; i<cities; i++)    It[i]=0.0; 
       for (i=1; i<cities; i++)
           for (j=0; j<cities; j++) 
           {
             
                  It[i]+=alpha*(1.0*voos[j][i])*I[j];
                  if (I[i]>=It[i]) flag[i]=1;
                  else  flag[i]=0;
           }
           

               
        //communitary infection
       Saux=0; Iaux=0; Raux=0;
       for (i=1; i<cities; i++)
       {
           //SIR model
       Saux= S[i] - (beta/population) * S[i] *(I[i]+It[i]);
       Iaux= I[i] + (beta/population )* S[i]* (I[i]+It[i]) - gama*(I[i]+It[i]);
       Raux= R[i] + gama*(I[i]+It[i]);
       if  (Saux<=0.0) 
       {
           Saux=0.0;
           Iaux=population;
       }
       S[i]=Saux;
       I[i]=Iaux;
       R[i]=Raux;
    }
    
    //print % of infected
     printf ("%d  ", time);
    for (i=1; i<155; i++)
        printf ("%lf        ", 100*I[i]/population);
    printf ("\n");
  }
          
}