/* FILE: h041proto.c   begins                        */
/* gcc -g h041.c -o a -L../lib -lnr -lm  */
#define NRANSI
#define NP 100
#include <stdio.h>
#include <stdlib.h>
#include <stdlib.h>
#include <limits.h>

#include <time.h>
#include <math.h>
#include "../util/nr.h"
#include "../util/nrutil.h"
#include "../util/matutl6.c"
void hist3(float *v, long m, long nrbins, int p, char *title)

/* Voit kayttaa tata rutiinia histogramman piirtamiseen */

{
  FILE *fp;
  float **xy;
  float min=v[1], max=v[1];
  int i,j;
  char *fname="hist3.dat";
  xy=matrix(1,nrbins+2,1,2);
  for(i=1;i<(nrbins+3);i++) { xy[i][1]=i-1; xy[i][2]=0; } 
  for(i=1;i<=m;i++) { 
    if(v[i]<min) min=v[i];
    if(v[i]>max) max=v[i];
  }
  for(i=1;i<=m;i++)
    for(j=0;j<nrbins;j++)
      if((v[i]-min)>=(((max-min)/nrbins)*j)&&
	 (v[i]-min)<(((max-min)/nrbins)*(j+1))) xy[j+2][2]++;
  /* If you omit the next line, on the x-axis we have indices of bins: */
  for(i=1;i<(nrbins+3);i++) { xy[i][1]=min+(i-1.)*(max-min)/(nrbins+2.); } 
  fp = fopen(fname,"w");
  for (i=1;i<(nrbins+3);i++)
    fprintf(fp,"%f   %f\n",xy[i][1],xy[i][2]); 
  fclose(fp);
  fp=fopen("gnuplt.cmd","w");
  if (p==1)
    {
      fprintf(fp,"set terminal postscript\n");
      fprintf(fp,"set output \"a.ps\"\n");
      fprintf(fp,"set title \"%s %s \"  \n",title,fname);
      fprintf(fp,"set timestamp\n");
      fprintf(fp,"set grid\n");
      fprintf(fp,"plot \"%s\" w boxes lw 3 \n",fname);
    }
  if (p==0)
    {
      fprintf(fp,"set title \"%s %s\"  \n",title,fname);
      fprintf(fp,"set timestamp\n");
      fprintf(fp,"set grid\n");
      fprintf(fp,"plot \"%s\" w boxes lw 3 \n",fname);
      fprintf(fp,"pause -1");
    } 
  fclose(fp);
  system("gnuplot gnuplt.cmd");
  free_matrix(xy,1,nrbins+2,1,2);
}


void NextXYD(float *x, float *y, float *dist)
{
  float r=rdm(0.1,0.6), phi=2*M_PI*rdm(0.0,1.0), dx,dy;
  dx= r*cos(phi); dy= r*sin(phi); *dist=r;
  if (fabs(*x+dx-0.5) <0.5) *x+=dx;
  else *x-=dx;
  if (fabs(*y+dy-0.5) <0.5) *y+=dy;
  else *y-=dy; 
}
float Dist(int i, int j, float **Persons,int m)
{
  float x=Persons[i][2], y=Persons[i][3],
         u=Persons[j][2],v=Persons[j][3];
  return pow((x-u)*(x-u)+(y-v)*(y-v),0.5);
}

void UpdateExposure(float **Persons, int m)
{
  int i,j;
  for (i=1;i<=m;i++){
    for (j=1;j<=m;j++)
      if ((j !=i) && (Persons[j][4]>0.5) && ( Dist(i,j,Persons,m)< 0.01))
	{Persons[i][4]+= 0.1;// j=m;
	}
  }
  /* Only one exposure per step is counted.*/
}
void InitPersons(float **Persons, int m)
{
  int i;
  for (i=1;i<=m;i++)
    {Persons[i][1]=i;
    Persons[i][2]=rdm(0.01, 0.99); /* Initial x -coord*/
    Persons[i][3]=rdm(0.01, 0.99); /* Initial y -coord*/
    Persons[i][4]=rdm(0.1, 0.6);   /* Initial rate of exposure */
    Persons[i][5]=0.0;             /* Initial traversed distance */
    }
}
void OneStep(float **Persons, int m)
{
  int i;
  float x,y,r;
  for (i=1;i<=m;i++)
    {x=Persons[i][2]; y=Persons[i][3];
    // printf("x= %12.4f, y= %12.4f\n",x,y);
    NextXYD(&x,&y,&r);
    Persons[i][2]=x; Persons[i][3]=y; /* Update locations */

    Persons[i][5]+=r; /* Update the distance travelled*/
    //     printf("x= %12.4f, y= %12.4f, dist=%12.4f\n",x,y,Persons[i][5]);

    UpdateExposure(Persons,m);
    //getchar();
    }
}

void DoGraph(float **Persons,int m, int iter)
{
  char *mystr= (char *)malloc(20);  
  float *x, *y;
  int i;
  x=vector(1,m);// y=vector(1,m);
  for (i=1;i<=m;i++)
    {x[i]=Persons[i][4]; 

    }
  sprintf(mystr,"Step %i",iter);
  hist3(x,m,m,0,mystr);
  free_vector(x,1,m);
  
}
void PlotDE(float **Persons,int m)
{
  int i;
  float x, y;
  FILE *fp;
  char *fname="de.dat";
  fp =fopen(fname,"w");
  for (i=1;i<=m;i++)
    {
      x=Persons[i][5];
      y=Persons[i][4];
      fprintf(fp,"%12.4f  %12.4f \n",x,y);
    }
  fclose(fp);
  fname="gnuplot.cmd";
  fp =fopen(fname,"w");
  fprintf(fp,"set title \"Exposure as a function of distance \" \n");
  fprintf(fp,"set timestamp \n set grid \n");
  fprintf(fp,"plot \"%s\" w p ps 3\n","de.dat");
  fprintf(fp,"pause -1\n");
  fclose(fp);
  system("gnuplot gnuplot.cmd");
}

void ShowLocations(float **Persons, int m)
{

/* Tee tahan aliohjelma, joka piirtää henkiloiden paikat */
  
}

void ShowInf(float **Persons,int m)
{

/* Tee tahan histogramman piirto infektion asteesta */

}

void DistDist(float **Persons,int m)
{
  int i;
  float *vec;
  float sum=0,var=0,mean=0;
  vec=vector(1,m);
  for(i=1;i<=m;i++) {
    vec[i]=Persons[i][5];
    sum+=Persons[i][5];
  }
  mean=(sum/(float)m);
  for(i=1;i<=m;i++)
    var+=((mean-Persons[i][5])*(mean-Persons[i][5]));
  var=(var/(float)(m-1));
  printf("\nAverage of distance travelled: %12.4f\n",mean);
  printf("Variance of distance travelled:  %12.4f\n",var);
  hist3(vec,m,10,0,"Distance travelled");
  free_vector(vec,1,m);
}

int main(void)
{
  int i,m=20;
  float **Persons;
  init_srand();
  Persons=matrix(1,m,1,5);
  InitPersons(Persons,m);
  for (i=1;i<=100;i++)
    {
      if (i==1) {DoGraph(Persons,m,1);
                 showmat(Persons,m,5); }
      OneStep(Persons,m);
      if ((i%20)==0) {
        printf("%d \n",i);
        DoGraph(Persons,m,i);
        ShowLocations(Persons,m);
        showmat(Persons,m,5);        
        ShowInf(Persons,m);
      }
    }
  PlotDE(Persons,m);
   showmat(Persons,m,5); 
  DistDist(Persons,5);
  free_matrix(Persons,1,m,1,5);
  return 0;
}
#undef NRANSI
/* FILE: h041proto.c ends                                       */