// Simulacao de fluxo de pessoas 
// Autor: Luiz Duczmal e leandro Alves Pereira
// 2004/04/26
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#include <gl/glut.h>
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876
long  nx,ny, pixsize;
int depth;
long n2, peoples, auxi, auxj, trocak, u,u2, nsimsmax, nsims;
long g=0,replicacoes , lance,lancemax,lancelen,lancelenmax;
long aux,pause1,pause2;
int printflag,exibesimulacao, r1, r2;
int somastatus, t1,t2,i,j,k,a,b1,b3,i00,j00,i000,j000,ntries;
char layout[50], trajetoria[50],campo[50];
long idum; //Semente da geracao de numeros aleatorios
double maximo = -1, pEntra, temp, somatemp, somataux , somataux2, mataux[3][3], mataux2[5][5], matw2[5][5], kd=100,ks=0.0, // pEntra=PInicial*pSai
       pSai,media,    //
       pInicial; //
double auxr;
int  ocup[100][100], tabn[100][100],tipo[100][100], traj[100][100], dur[100][100],mark[100][100], auxw[100][100];
int ni,nj;
long tab[100][100];
long orderi[10000],orderj[10000],smat;
double D[200][200], matw[3][3], mat3[3][3],mat5[5][5],mat7[7][7],soma,acc3[9],acc5[25],r,S[200][200];
float rnd(long idum0);
double logpow( double x1, double x2 );
double power( double x1, double x2 );
double sqr( double x );
int movever(int i);
int movehor(int i);
int mov(int range, int i0, int j0, int *ni, int *nj);
void plote(int x, int y, int t ,int c);
void finner(int i000, int j000, int *i00, int *j00);
void fbound(int i000, int j000, int *i00, int *j00);
FILE *par;
FILE *sai;
FILE *leia;
FILE *leiadir;
FILE *leiacampo;
// Função callback chamada para fazer o desenho
void redefine (int pa)
{
 int t;
  soma = 0.0; 
  smat=9;
  switch (pa){
   case 1:
    if ( ( tab[orderi[k]][orderj[k+1]] + tab[orderi[k]][orderj[k+2]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.03; matw[1][0]=0.07; matw[2][0]=0.10;
    matw[0][1]=0.05; matw[1][1]=0.15; matw[2][1]=0.40;
    matw[0][2]=0.03; matw[1][2]=0.07; matw[2][2]=0.10;       
    matw2[0][0]=0.00; matw2[1][0]=0.00; matw2[2][0]=0.00; matw2[3][0]=0.06; matw2[4][0]=0.10;
    matw2[0][1]=0.00; matw2[1][1]=0.00; matw2[2][1]=0.00; matw2[3][1]=0.05; matw2[4][1]=0.15;
    matw2[0][2]=0.00; matw2[1][2]=0.00; matw2[2][2]=0.01; matw2[3][2]=0.07; matw2[4][2]=0.20;      
    matw2[0][3]=0.00; matw2[1][3]=0.00; matw2[2][3]=0.00; matw2[3][3]=0.05; matw2[4][3]=0.15;      
    matw2[0][4]=0.00; matw2[1][4]=0.00; matw2[2][4]=0.00; matw2[3][4]=0.06; matw2[4][4]=0.10;         
    break;
   case 2:
    if ( ( tab[orderi[k-1]][orderj[k-1]] + tab[orderi[k-2]][orderj[k]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.10; matw[1][0]=0.40; matw[2][0]=0.10;
    matw[0][1]=0.07; matw[1][1]=0.15; matw[2][1]=0.07;
    matw[0][2]=0.03; matw[1][2]=0.05; matw[2][2]=0.03;
    matw2[0][0]=0.10; matw2[1][0]=0.15; matw2[2][0]=0.20; matw2[3][0]=0.15; matw2[4][0]=0.10;
    matw2[0][1]=0.06; matw2[1][1]=0.05; matw2[2][1]=0.07; matw2[3][1]=0.05; matw2[4][1]=0.06;
    matw2[0][2]=0.00; matw2[1][2]=0.00; matw2[2][2]=0.01; matw2[3][2]=0.00; matw2[4][2]=0.00;      
    matw2[0][3]=0.00; matw2[1][3]=0.00; matw2[2][3]=0.00; matw2[3][3]=0.00; matw2[4][3]=0.00;      
    matw2[0][4]=0.00; matw2[1][4]=0.00; matw2[2][4]=0.00; matw2[3][4]=0.00; matw2[4][4]=0.00;         
    break;
   case 3: 
      if ( ( tab[orderi[k+1]][orderj[k]] + tab[orderi[k+2]][orderj[k]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.03; matw[1][0]=0.05; matw[2][0]=0.03;
    matw[0][1]=0.07; matw[1][1]=0.15; matw[2][1]=0.07;
    matw[0][2]=0.10; matw[1][2]=0.40; matw[2][2]=0.10;
    matw2[0][0]=0.00; matw2[1][0]=0.00; matw2[2][0]=0.00; matw2[3][0]=0.00; matw2[4][0]=0.00;
    matw2[0][1]=0.00; matw2[1][1]=0.00; matw2[2][1]=0.00; matw2[3][1]=0.00; matw2[4][1]=0.00;
    matw2[0][2]=0.00; matw2[1][2]=0.00; matw2[2][2]=0.01; matw2[3][2]=0.00; matw2[4][2]=0.00;      
    matw2[0][3]=0.06; matw2[1][3]=0.05; matw2[2][3]=0.07; matw2[3][3]=0.05; matw2[4][3]=0.06;      
    matw2[0][4]=0.10; matw2[1][4]=0.15; matw2[2][4]=0.20; matw2[3][4]=0.15; matw2[4][4]=0.10;         
    break;
   case 4:
     if ( ( tab[orderi[k]][orderj[k-1]] + tab[orderi[k]][orderj[k-2]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.10; matw[1][0]=0.07; matw[2][0]=0.03;
    matw[0][1]=0.40; matw[1][1]=0.15; matw[2][1]=0.05;
    matw[0][2]=0.10; matw[1][2]=0.07; matw[2][2]=0.03;   
    matw2[0][0]=0.10; matw2[1][0]=0.06; matw2[2][0]=0.00; matw2[3][0]=0.00; matw2[4][0]=0.00;
    matw2[0][1]=0.15; matw2[1][1]=0.05; matw2[2][1]=0.00; matw2[3][1]=0.00; matw2[4][1]=0.00;
    matw2[0][2]=0.20; matw2[1][2]=0.07; matw2[2][2]=0.01; matw2[3][2]=0.00; matw2[4][2]=0.00;      
    matw2[0][3]=0.15; matw2[1][3]=0.05; matw2[2][3]=0.00; matw2[3][3]=0.00; matw2[4][3]=0.00;      
    matw2[0][4]=0.10; matw2[1][4]=0.06; matw2[2][4]=0.00; matw2[3][4]=0.00; matw2[4][4]=0.00; 
    break;    
    case 5:
   if ( ( tab[orderi[k-1]][orderj[k+1]] + tab[orderi[k-2]][orderj[k+2]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.07; matw[1][0]=0.10; matw[2][0]=0.40;
    matw[0][1]=0.03; matw[1][1]=0.15; matw[2][1]=0.10;
    matw[0][2]=0.05; matw[1][2]=0.03; matw[2][2]=0.07;
    matw2[0][0]=0.00; matw2[1][0]=0.06; matw2[2][0]=0.10; matw2[3][0]=0.15; matw2[4][0]=0.20;
    matw2[0][1]=0.00; matw2[1][1]=0.00; matw2[2][1]=0.05; matw2[3][1]=0.07; matw2[4][1]=0.15;
    matw2[0][2]=0.00; matw2[1][2]=0.00; matw2[2][2]=0.01; matw2[3][2]=0.05; matw2[4][2]=0.10;      
    matw2[0][3]=0.00; matw2[1][3]=0.00; matw2[2][3]=0.00; matw2[3][3]=0.00; matw2[4][3]=0.06;      
    matw2[0][4]=0.00; matw2[1][4]=0.00; matw2[2][4]=0.00; matw2[3][4]=0.00; matw2[4][4]=0.00;  
    break;
  
   case 6:
        if ( ( tab[orderi[k+1]][orderj[k+1]] + tab[orderi[k+2]][orderj[k+2]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.05; matw[1][0]=0.03; matw[2][0]=0.07;
    matw[0][1]=0.03; matw[1][1]=0.15; matw[2][1]=0.10;
    matw[0][2]=0.07; matw[1][2]=0.10; matw[2][2]=0.40;
    matw2[0][0]=0.00; matw2[1][0]=0.00; matw2[2][0]=0.00; matw2[3][0]=0.00; matw2[4][0]=0.00;
    matw2[0][1]=0.00; matw2[1][1]=0.00; matw2[2][1]=0.00; matw2[3][1]=0.00; matw2[4][1]=0.06;
    matw2[0][2]=0.00; matw2[1][2]=0.00; matw2[2][2]=0.01; matw2[3][2]=0.05; matw2[4][2]=0.10;      
    matw2[0][3]=0.00; matw2[1][3]=0.00; matw2[2][3]=0.05; matw2[3][3]=0.07; matw2[4][3]=0.15;      
    matw2[0][4]=0.00; matw2[1][4]=0.06; matw2[2][4]=0.10; matw2[3][4]=0.15; matw2[4][4]=0.20;       
    break;
  case 7:  
       if ( ( tab[orderi[k-1]][orderj[k-1]] + tab[orderi[k-2]][orderj[k-2]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.40; matw[1][0]=0.10; matw[2][0]=0.07;
    matw[0][1]=0.10; matw[1][1]=0.15; matw[2][1]=0.03;
    matw[0][2]=0.07; matw[1][2]=0.03; matw[2][2]=0.05;
    matw2[0][0]=0.20; matw2[1][0]=0.15; matw2[2][0]=0.10; matw2[3][0]=0.06; matw2[4][0]=0.00;
    matw2[0][1]=0.15; matw2[1][1]=0.07; matw2[2][1]=0.05; matw2[3][1]=0.00; matw2[4][1]=0.00;
    matw2[0][2]=0.10; matw2[1][2]=0.05; matw2[2][2]=0.01; matw2[3][2]=0.00; matw2[4][2]=0.00;      
    matw2[0][3]=0.06; matw2[1][3]=0.00; matw2[2][3]=0.00; matw2[3][3]=0.00; matw2[4][3]=0.00;      
    matw2[0][4]=0.00; matw2[1][4]=0.00; matw2[2][4]=0.00; matw2[3][4]=0.00; matw2[4][4]=0.00;  
    break;
    case 8:
         
    if ( ( tab[orderi[k+1]][orderj[k-1]] + tab[orderi[k+2]][orderj[k-2]]) > 0 ){  
    smat=9;
    }else
    {smat=25;}
    matw[0][0]=0.07; matw[1][0]=0.03; matw[2][0]=0.05;
    matw[0][1]=0.10; matw[1][1]=0.15; matw[2][1]=0.03;
    matw[0][2]=0.40; matw[1][2]=0.10; matw[2][2]=0.07;
    matw2[0][0]=0.00; matw2[1][0]=0.00; matw2[2][0]=0.00; matw2[3][0]=0.00; matw2[4][0]=0.00;
    matw2[0][1]=0.06; matw2[1][1]=0.00; matw2[2][1]=0.00; matw2[3][1]=0.00; matw2[4][1]=0.00;
    matw2[0][2]=0.10; matw2[1][2]=0.05; matw2[2][2]=0.01; matw2[3][2]=0.00; matw2[4][2]=0.00;      
    matw2[0][3]=0.15; matw2[1][3]=0.07; matw2[2][3]=0.05; matw2[3][3]=0.00; matw2[4][3]=0.00;      
    matw2[0][4]=0.20; matw2[1][4]=0.15; matw2[2][4]=0.10; matw2[3][4]=0.06; matw2[4][4]=0.00;  
    break;
  }   
   kd = kd*0.999;
   ks = ks*0.999;
   if (smat==9) {
   somataux = 0.0;
   mataux[0][0]= matw[0][0] * exp(kd*D[orderi[k-1]][orderj[k-1]]) * exp(ks*S[orderi[k-1]][orderj[k-1]]);
   mataux[0][1]= matw[0][1] * exp(kd*D[orderi[k]][orderj[k-1]]) * exp(ks*S[orderi[k]][orderj[k-1]]);
   mataux[0][2]= matw[0][2] * exp(kd*D[orderi[k+1]][orderj[k-1]]) * exp(ks*S[orderi[k+1]][orderj[k-1]]);
   mataux[1][0]= matw[1][0] * exp(kd*D[orderi[k-1]][orderj[k]]) * exp(ks*S[orderi[k-1]][orderj[k]]);
   mataux[1][1]= matw[1][1] * exp(kd*D[orderi[k]][orderj[k]]) * exp(ks*S[orderi[k]][orderj[k]]);
   mataux[1][2]= matw[1][2] * exp(kd*D[orderi[k+1]][orderj[k]]) * exp(ks*S[orderi[k+1]][orderj[k]]);
   mataux[2][0]= matw[2][0] * exp(kd*D[orderi[k-1]][orderj[k+1]]) * exp(ks*S[orderi[k-1]][orderj[k+1]]);
   mataux[2][1]= matw[2][1] * exp(kd*D[orderi[k]][orderj[k+1]]) * exp(ks*S[orderi[k]][orderj[k+1]]);
   mataux[2][2]= matw[2][2] * exp(kd*D[orderi[k+1]][orderj[k+1]]) * exp(ks*S[orderi[k+1]][orderj[k+1]]);
    
   }
   else{
   somataux2= 0.0;
   mataux2[0][0]= matw2[0][0]* exp(kd*D[orderi[k-2]][orderj[k-2]]) * exp(ks*S[orderi[k-2]][orderj[k-2]]);
   mataux2[0][1]= matw2[0][1]* exp(kd*D[orderi[k-1]][orderj[k-2]]) * exp(ks*S[orderi[k-1]][orderj[k-2]]);//
   mataux2[0][2]= matw2[0][2]* exp(kd*D[orderi[k]][orderj[k-2]]) * exp(ks*S[orderi[k]][orderj[k-2]]);
   mataux2[0][3]= matw2[0][3]* exp(kd*D[orderi[k+1]][orderj[k-2]]) * exp(ks*S[orderi[k+1]][orderj[k-2]]);
   mataux2[0][4]= matw2[0][4]* exp(kd*D[orderi[k+2]][orderj[k-2]]) * exp(ks*S[orderi[k+2]][orderj[k-2]]);
   mataux2[1][0]= matw2[1][0]* exp(kd*D[orderi[k-2]][orderj[k-1]]) * exp(ks*S[orderi[k-2]][orderj[k-1]]);
   mataux2[1][1]= matw2[1][1]* exp(kd*D[orderi[k-1]][orderj[k-1]]) * exp(ks*S[orderi[k-1]][orderj[k-1]]);
   mataux2[1][2]= matw2[1][2]* exp(kd*D[orderi[k]][orderj[k-1]]) * exp(ks*S[orderi[k]][orderj[k-1]]);
   mataux2[1][3]= matw2[1][3]* exp(kd*D[orderi[k+1]][orderj[k-1]]) * exp(ks*S[orderi[k+1]][orderj[k-1]]);
   mataux2[1][4]= matw2[1][4]* exp(kd*D[orderi[k+2]][orderj[k-1]]) * exp(ks*S[orderi[k+2]][orderj[k-1]]); 
   mataux2[2][0]= matw2[2][0]* exp(kd*D[orderi[k-2]][orderj[k]]) * exp(ks*S[orderi[k-2]][orderj[k]]);
   mataux2[2][1]= matw2[2][1]* exp(kd*D[orderi[k-1]][orderj[k]]) * exp(ks*S[orderi[k-1]][orderj[k]]);//
   mataux2[2][2]= matw2[2][2]* exp(kd*D[orderi[k]][orderj[k]]) * exp(ks*S[orderi[k]][orderj[k]]);
   mataux2[2][3]= matw2[2][3]* exp(kd*D[orderi[k+1]][orderj[k]]) * exp(ks*S[orderi[k+1]][orderj[k]]);
   mataux2[2][4]= matw2[2][4]* exp(kd*D[orderi[k+2]][orderj[k]]) * exp(ks*S[orderi[k+2]][orderj[k]]); 
   mataux2[3][0]= matw2[3][0]* exp(kd*D[orderi[k-2]][orderj[k+1]]) * exp(ks*S[orderi[k-2]][orderj[k+1]]);//
   mataux2[3][1]= matw2[3][1]* exp(kd*D[orderi[k-1]][orderj[k+1]]) * exp(ks*S[orderi[k-1]][orderj[k+1]]);
   mataux2[3][2]= matw2[3][2]* exp(kd*D[orderi[k]][orderj[k+1]]) * exp(ks*S[orderi[k]][orderj[k+1]]);
   mataux2[3][3]= matw2[3][3]* exp(kd*D[orderi[k+1]][orderj[k+1]]) * exp(ks*S[orderi[k+1]][orderj[k+1]]); 
   mataux2[3][4]= matw2[3][4]* exp(kd*D[orderi[k+2]][orderj[k+1]]) * exp(ks*S[orderi[k+2]][orderj[k+1]]);//
   mataux2[4][0]= matw2[4][0]* exp(kd*D[orderi[k-2]][orderj[k+2]]) * exp(ks*S[orderi[k-2]][orderj[k+2]]);
   mataux2[4][1]= matw2[4][1]* exp(kd*D[orderi[k-1]][orderj[k+2]]) * exp(ks*S[orderi[k-1]][orderj[k+2]]);
   mataux2[4][2]= matw2[4][2]* exp(kd*D[orderi[k]][orderj[k+2]]) * exp(ks*S[orderi[k]][orderj[k+2]]); 
   mataux2[4][3]= matw2[4][3]* exp(kd*D[orderi[k+1]][orderj[k+2]]) * exp(ks*S[orderi[k+1]][orderj[k+2]]);//
   mataux2[4][4]= matw2[4][4]* exp(kd*D[orderi[k+2]][orderj[k+2]]) * exp(ks*S[orderi[k+2]][orderj[k+2]]);
 }
    if (smat == 9){
   for (t1=0;t1<=2;t1++)
    for (t2=0;t2<=2;t2++){
      somataux += mataux[t1][t2];
    }
   for (t1=0;t1<=2;t1++)
    for (t2=0;t2<=2;t2++){
      mat3[t1][t2]= mataux[t1][t2]/somataux;
      }
     for(t=0;t<9;t++){soma+=mat3[t/3][t%3];acc3[t]=soma;}         
    }
    else{
     for (t1=0;t1<=4;t1++)
    for (t2=0;t2<=4;t2++){
      somataux2 += mataux2[t1][t2];
    }
   for (t1=0;t1<=4;t1++)
    for (t2=0;t2<=4;t2++){
      mat5[t1][t2]= mataux2[t1][t2]/somataux2;
      }
   for(t=0;t<25;t++){soma+=mat5[t/5][t%5];acc5[t]=soma;}               
    }
} 
void Desenha(void)
{
     glMatrixMode(GL_MODELVIEW);
     glLoadIdentity();
   // Limpa a janela de visualização com a cor de fundo especificada
     glClear(GL_COLOR_BUFFER_BIT);
    // Especifica que a cor corrente é vermelha
     //         R     G     B
     glColor3f(0.7f, 0.7f, 0.8f);
     // Desenha n pontos com a cor corrente
g=0;
refaz:;
temp=0;
kd=100;
for(r1=0;r1<=200;r1++){
        for(r2=0;r2<=200;r2++){
         D[r1][r2] = 0;
       //  S[r1][r2] = 0;
        }    
     }    
g++;
n2=0; //Numero total de celulas ativas
for(i=0;i<nx;i++){ //Inicializa os estados das celulas
  for(j=0;j<ny;j++)
    if(rnd(idum)<pInicial && (tipo[i][j] & 1)==1){
       tab[i][j]=1;
       orderi[n2]=i;//Inicializa a ordem de processamento
       orderj[n2]=j;
       n2++;
    }
    else tab[i][j]=0;
}
peoples = n2;
nsims=0;
iniciodes:;
nsims++;
for(k=0;k<n2;k++){ //embaralha a ordem de processamento
   auxi=orderi[k];
   auxj=orderj[k];
   trocak=(int)(rnd(idum)*n2);
   orderi[k]=orderi[trocak];
   orderj[k]=orderj[trocak];
   orderi[trocak]=auxi;
   orderj[trocak]=auxj;
}
for(k=0;k<n2;k++){//varredura
   iniciovarr:;
   i00=orderi[k];
   j00=orderj[k];
   //if(tipo[i00][j00] & 4){}// in
   if(tipo[i00][j00] & 8){ // out
      if(rnd(idum)<pSai){
         tab[i00][j00]=0;
         n2--;
         if(n2>k){
            orderi[k]=orderi[n2];
            orderj[k]=orderj[n2];
            goto iniciovarr;
         }
      }
   }
   else
   {
      redefine(traj[orderi[k]][orderj[k]]); 
      D[orderi[k]][orderj[k]]+=0.1;
     if (D[orderi[k]][orderj[k]]> maximo){
               maximo = D[orderi[k]][orderj[k]];
      }                              
      mov(2, i00, j00, &ni, &nj);
      orderi[k]=ni;
      orderj[k]=nj;
   }
   if (exibesimulacao) plote(i00,j00,pixsize,8);
}
/*********************************/
if (exibesimulacao){ 
 for(i=0;i<ny;i++){
   for(j=0;j<nx;j++)
    plote(i,j,pixsize,tab[i][j]);
    glFlush();
   for(i=0;i<nx;i++){
   for(j=0;j<ny;j++){
    aux = tipo[i][j];
    if (aux ==3) aux =1;
    plote(i,j,pixsize,aux);
   }
 }
 }
} 
/**********************************/
 if (n2>0) goto iniciodes;
 temp = (nsims*0.298142)/60;
 //fprintf(sai,"numero de simulações %d\n ",nsims);
 fprintf(sai," %d tempo total de evacuação e numero de simulações: %f minutos / %d\n",g,temp, nsims);
 //fprintf(sai,"          numero de pessoas no dominio: %d pessoas\n",peoples);
  somatemp += temp;
 if (g<replicacoes) goto refaz;
 media = somatemp/replicacoes;
 fprintf(sai,"tempo medio de evacuação total: %f minutos\n",media);
 fprintf(sai," ");
/*  for(i=0;i<=256;i++){
        for(j=0;j<=256;j++){
         fprintf(sai," %d ",D[j][i]);
        }    
     fprintf(sai,"\n");      
     }    */
fim:;

}
// Inicializa parâmetros de rendering
void Inicializa (void)
{   
    // Define a cor de fundo da janela de visualização como preta
    glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
    //glClearColor(100, 20, 180, 1.0f);
}

// Função callback chamada quando o tamanho da janela é alterado 
void AlteraTamanhoJanela(GLsizei w, GLsizei h)
{
                   // Evita a divisao por zero
                   if(h == 0) h = 1;
                   // Especifica as dimensões da Viewport
                   glViewport(0, 0, w, h);
                   // Inicializa o sistema de coordenadas
                   glMatrixMode(GL_PROJECTION);
                   glLoadIdentity();
                   // Estabelece a janela de seleção (left, right, bottom, top)
                   if (w <= h) 
                           gluOrtho2D (0.0f, 250.0f, 0.0f, 250.0f*h/w);
                   else 
                           gluOrtho2D (0.0f, 250.0f*w/h, 0.0f, 250.0f);
}
// Programa Principal 
int main(void)
{
     par=fopen("parametros.txt","r");
     fscanf(par,"%s",layout);
     fscanf(par,"%s",trajetoria);
     fscanf(par,"%s",campo);
     fscanf(par,"%d",&nx);
     fscanf(par,"%d",&ny);
     fscanf(par,"%d",&pixsize);
     fscanf(par,"%d",&depth);
     fscanf(par,"%d",&nsimsmax);
     fscanf(par,"%d",&lancemax);
     fscanf(par,"%d",&lancelenmax);
     fscanf(par,"%d",&pause1);
     fscanf(par,"%d",&pause2);
     fscanf(par,"%d",&printflag);
     fscanf(par,"%d",&ntries);
     fscanf(par,"%d",&idum);
     fscanf(par,"%lg",&pEntra);
     fscanf(par,"%lg",&pSai);
     fscanf(par,"%lg",&pInicial);
     fscanf(par,"%d",&smat);
     fscanf(par,"%d",&replicacoes);
     fscanf(par,"%d",&exibesimulacao);
     fclose(par);
    leia=fopen(layout,"r");
     for(i=0;i<ny;i++)
        for(j=0;j<nx;j++)
         fscanf(leia,"%d",&tipo[j][i]);
     fclose(leia);    
     leiadir=fopen(trajetoria,"r");
     for(i=0;i<ny;i++)
        for(j=0;j<nx;j++)
         fscanf(leiadir,"%d",&traj[j][i]);
     fclose(leiadir);
     leiacampo=fopen(campo,"r");
    for(i=0;i<ny;i++)
       for(j=0;j<nx;j++)
         fscanf(leiacampo,"%lg",&S[i][j]);
     fclose(leiacampo);
     for(j=0;j<ny;j++)
        for(i=0;i<nx;i++)
           ocup[i][j]=1;
           
     //inicialização aqui (coisas que só são acionadas uma única vez)
     sai=fopen("saida.txt","w");

     glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
     glutInitWindowSize(700,700);
     glutInitWindowPosition(10,10);
     glutCreateWindow("Simulação de tráfego de pessoas ");
      glutDisplayFunc(Desenha);
      glutReshapeFunc(AlteraTamanhoJanela);
      Inicializa();
      glutMainLoop();
      fclose(sai);
      printf("concluido");
}
double logpow(double x1, double x2){ /* power(2,3)=8 */
  return( (log(x1)*x2) );
}
double power(double x1, double x2){ /* power(2,3)=8 */
  return( exp(log(x1)*x2) );
}
double sqr( double x ){
  return (x*x);
}
int movever(int j){
  double auxr;
  auxr=rnd(idum);
  if(auxr<0.25)return(j-1);
  else if(auxr<0.75)return(j);
  else return(j+1);
}
int movehor(int i){
  double auxr;
  auxr=rnd(idum);
  if(auxr<0.10)return(i-1);
  else if(auxr<0.20)return(i);
  else return(i+1);
}
int mov(int range, int i0, int j0, int *ni, int *nj){
  double aux1,aux2;
  int p,q,i,j,k,cont=0,i00,j00,i0min,i0max,j0min,j0max;
  i0min=i0-2*range; if(i0min<0)i0min=0;
  i0max=i0+2*range; if(i0max>=nx)i0max=nx-1;
  j0min=j0-2*range; if(j0min<0)j0min=0;
  j0max=j0+2*range; if(j0max>=ny)j0max=ny-1;
  for(i=i0min;i<=i0max;i++)// Inicializa mark[][]
     for(j=j0min;j<=j0max;j++)
        mark[i][j]=0;
  i0min=i0-range; if(i0min<0)i0min=0;
  i0max=i0+range; if(i0max>=nx)i0max=nx-1;
  j0min=j0-range; if(j0min<0)j0min=0;
  j0max=j0+range; if(j0max>=ny)j0max=ny-1;
  for(i=i0min;i<=i0max;i++){//Localiza as ativas proximas
     for(j=j0min;j<=j0max;j++){
        if(tab[i][j]==1 && i!=i0 && j!=j0)
        for(lance=0;lance<lancemax;lance++){
           i000=i00=i;
           j000=j00=j;
           for(lancelen=0;lancelen<lancelenmax;lancelen++){
              if(tipo[i00][j00] & 2 ){//fronteira do dominio
                 fbound(i000,j000,&i00,&j00);
                 mark[i00][j00]++;
              }
              else{//dentro do dominio
                 finner(i000,j000,&i00,&j00);
                 mark[i00][j00]++;
              }
           }
        }
     }
  }
  if(printflag){
  fprintf(sai,"%d %d:\n",i0,j0);
  for(i=i0min;i<=i0max;i++){
     for(j=j0min;j<=j0max;j++)fprintf(sai,"%d ",mark[i][j]);
     fprintf(sai,"\n");
  }
  fprintf(sai,"\n");
  }
  r=rnd(idum);
  k=0;
  if(smat==9){
    while(k<smat && r>acc3[k])k++;
    *ni=k/3-1+ i0;
    *nj=k%3-1+ j0;
    if(*ni!=i0 || *nj!=j0){
      tab[i0][j0]=0;
      if(*ni<0 || *ni>=nx || *nj<0 || *nj>=ny || mark[*ni][*nj]>0 ||
         tab[*ni][*nj]>0 || (tipo[*ni][*nj]&1)==0){
        *ni=i0;
        *nj=j0;
        tab[i0][j0]=1;
      }
      else{
        tab[*ni][*nj]=1;
      }
    }
  }
  else{//smat==25
    while(k<smat && r>acc5[k])k++;
    *ni=k/5-2+ i0;
    *nj=k%5-2+ j0;
    if(*ni!=i0 || *nj!=j0){
      tab[i0][j0]=0;
      if(*ni<0 || *ni>=nx || *nj<0 || *nj>=ny || mark[*ni][*nj]>0 ||
         tab[*ni][*nj]>0 || (tipo[*ni][*nj]&1)==0){
        *ni=i0;
        *nj=j0;
        tab[i0][j0]=1;
      }
      else{
        tab[*ni][*nj]=1;
      }
    }
  }
}
void finner(int i000, int j000, int *i00, int *j00){
  int i,j,k;
  double r;
  r=rnd(idum);
  k=0;
  if(smat==9){
    while(k<smat && r>acc3[k])k++;
    *i00=k/3-1+ *i00;
    *j00=k%3-1+ *j00;
  }
  else{//smat==25
    while(k<smat && r>acc5[k])k++;
    *i00=k/5-2+ *i00;
    *j00=k%5-2+ *j00;
  }
}
void fbound(int i000, int j000, int *i00, int *j00){
 int i,j,k,tent=0,tentmax=5;
  double r;
  do{
     r=rnd(idum);
     k=0;
     if(smat==9){
       while(k<smat && r>acc3[k])k++;
       i=k/3-1+ *i00;
       j=k%3-1+ *j00;
     }
     else{//smat==25
       while(k<smat && r>acc5[k])k++;
       i=k/5-2+ *i00;
       j=k%5-2+ *j00;
     }
     tent++;
  }while((tab[i][j] & 1)!=1 && tent<tentmax);
  if(tent<tentmax){
     *i00=i;     *j00=j;
  }
}
float rnd(long idum0){
   long k;
   float ans;
   do{
     idum0 ^= MASK;
     k=(idum0)/IQ;
     idum0=IA*(idum0-k*IQ)-IR*k;
     if(idum0<0) idum0+=IM;
     ans=AM*(idum0);
     idum0^=MASK;
     idum=idum0;
   }while(ans>0.999999);
   return(ans);
}
void plote(int x, int y, int t ,int c){
  int i,j,a,b;
  double R=0.0,G=0.0,B=0.0;
  if(c==8)glColor3f(1.0,1.0,1.0);
  else{
     if(c>0){R=G=B=0.3;}
     c/=2; if(c%2)R+=0.7;
     c/=2; if(c%2)G+=0.7;
     c/=2; if(c%2)B+=0.7;
     glColor3f(R,G,B);
  }
  a=x*t;
  b=y*t;
  glBegin(GL_QUADS);
               glVertex2i(a+10,b+10);
               glVertex2i(a+t+10,b+10);
               glVertex2i(a+t+10,b+t+10);
               glVertex2i(a+10,b+t+10);
  glEnd();
  }