// To Compile     gcc -lglut -lGL -fopenmp sim.c
//

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>
#include <sys/time.h>

#define N 512
#define ES 100.0 //空間の大きさ x,y,zそれぞれ-ES~ESの大きさを持つ空間になる

GLubyte image[N][N][4];
float q_odd[N][N],q_even[N][N];
unsigned char KEY = '1';
double all_t1, all_t2;

const float Q = 1.0;
const float ka = 8987.7424380; //クーロン定数
float theta = 0.0;
int  nz, ny, nx;
float q_max, q_min;

struct color{
  unsigned char R;
  unsigned char G;
  unsigned char B;
  float A;
} color[256];

struct C{
  float R;
  float G;
  float B;
  float A;
} C[N][N];

int step=0;

float prop[N][N][5];      // 熱伝導係数
float x=0.25;             //係数  delta_t/delta_x  
unsigned char kabe[N][N]; //障害物がある格子では PROP = 0 にしたい！ .... 未実装 20120427
unsigned char log_mode=1; //表示モード normal scale or log scale
int Px =511; //東尋坊　120
int Py =0; //東尋坊　230              ..... x,y 逆だったかな？


double gettimeofday_sec()
{
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return tv.tv_sec + (double)tv.tv_usec*1e-6;
}

void simulation()
{
  //double t1, t2;
  //t1 = gettimeofday_sec();
  
  float Qx, Qy, Qz;
  float rx, ry, rz, r, rzy;
  int i, j, k;
  float intensity,log_intensity;

  float h=0.75; // 境界補正

  /*  q_odd[0][0]=100.0;*/
  /*  q_even[0][0]=100.0;*/

  /*********************************************************************************************/
  /*********************************************************************************************/
  q_even[Px][Py]=255.0;/* FIXED Heat Source*/

  i=0;                         /******************************************/
  j=0;
  //q_odd[i][j]=0.0;

  for(j=1;j<N-1;j++){
    q_odd[i][j]=q_even[i][j]+
      x*(  (prop[i][j][0]*q_even[i][j-1]  +prop[i][j][1]*q_even[i][j+1] 
	    +prop[i][j][2]*q_even[i][j]*h   +prop[i][j][3]*q_even[i+1][j])-prop[i][j][4]*q_even[i][j]);
  }

  j=N-1;
  q_odd[i][j]=q_even[i][j]+
    x*(  (prop[i][j][0]*q_even[i][j-1]  +prop[i][j][1]*q_even[i][j]*h   
	  +prop[i][j][2]*q_even[i][j]*h   +prop[i][j][3]*q_even[i+1][j])-prop[i][j][4]*q_even[i][j]); 
  
#pragma omp parallel for private(j) 
  for(i = 1; i < N-1; i++){    /******************************************/
    j=0;
    q_odd[i][j]=q_even[i][j]+
      x*(  (prop[i][j][0]*q_even[i][j]*h    +prop[i][j][1]*q_even[i][j+1] 
	    +prop[i][j][2]*q_even[i-1][j] +prop[i][j][3]*q_even[i+1][j])-prop[i][j][4]*q_even[i][j]);

    for(j=1;j<N-1;j++){
      q_odd[i][j]=q_even[i][j]+
	x*( (prop[i][j][0]*q_even[i][j-1] +prop[i][j][1]*q_even[i][j+1] 
	     +prop[i][j][2]*q_even[i-1][j] +prop[i][j][3]*q_even[i+1][j])-prop[i][j][4]*q_even[i][j]);
    }

    j=N-1;
    q_odd[i][j]=q_even[j][j]+
      x*( (prop[i][j][0]*q_even[i][j-1] +prop[i][j][1]*q_even[i][j]*h   
	   +prop[i][j][2]*q_even[i-1][j] +prop[i][j][3]*q_even[i+1][j])-prop[i][j][4]*q_even[i][j]);
  }

  i=N-1;                        /******************************************/
  j=0;
  q_odd[i][j]=q_even[i][j]+
    x*( (prop[i][j][0]*q_even[i][j]*h   +prop[i][j][1]*q_even[i][j+1] 
	 +prop[i][j][2]*q_even[i-1][j] +prop[i][j][3]*q_even[i][j]*h  )-prop[i][j][4]*q_even[i][j]);

  for(j=1;j<N-1;j++){
      q_odd[i][j]=q_even[i][j]+
	x*( (prop[i][j][0]*q_even[i][j-1] +prop[i][j][1]*q_even[i][j+1] 
	     +prop[i][j][2]*q_even[i-1][j] +prop[i][j][3]*q_even[i][j]*h  )-prop[i][j][4]*q_even[i][j]);
  }

  j=N-1;
  q_odd[i][N-1]=q_even[i][j]+
    x*( (prop[i][j][0]*q_even[i][j-1] +prop[i][j][1]*q_even[i][j]*h   
	 +prop[i][j][2]*q_even[i-1][j] +prop[i][j][3]*q_even[i][j]*h  )-prop[i][j][4]*q_even[i][j]);

  q_odd[Px][Py]=255.0;


  /*********************************************************************************************/
  /*********************************************************************************************/
  q_odd[Px][Py]=255.0; /* FIXED Heat Source*/

  i=0;                          /******************************************/
  j=0;
  //q_even[i] for(j=1;j<N-1;j++){
   q_even[i][j]=q_odd[i][j]+
     x*( (prop[i][j][0]*q_odd[i][j-1] +prop[i][j][1]*q_odd[i][j+1] 
	  +prop[i][j][2]*q_odd[i][j]*h   +prop[i][j][3]*q_odd[i+1][j] )-prop[i][j][4]*q_odd[i][j]);
}

j=N-1;
  q_even[i][j]=q_odd[i][j]+
    x*( (prop[i][j][0]*q_odd[i][j-1] +prop[i][j][1]*q_odd[i][j]*h   
	 +prop[i][j][2]*q_odd[i][j]*h   +prop[i][j][3]*q_odd[i+1][j] )-prop[i][j][4]*q_odd[i][j]);

  
#pragma omp parallel for private(j) 
for(i = 1; i < N-1; i++){     /******************************************/
  j=0;
    q_even[i][j]=q_odd[i][j]+
      x*( (prop[i][j][0]*q_odd[i][j]*h   +prop[i][j][1]*q_odd[i][j+1] 
	   +prop[i][j][2]*q_odd[i-1][j] +prop[i][j][3]*q_odd[i+1][j] )-prop[i][j][4]*q_odd[i][j]);

    for(j=1;j<N-1;j++){
      q_even[i][j]=q_odd[i][j]+
	x*( (prop[i][j][0]*q_odd[i][j-1] +prop[i][j][1]*q_odd[i][j+1] 
	     +prop[i][j][2]*q_odd[i-1][j] +prop[i][j][3]*q_odd[i+1][j] )-prop[i][j][4]*q_odd[i][j]);
      //          if(i==N/2) q_even[i][j]=255.0;
    }

    j=N-1;
    q_even[i][j]=q_odd[i][j]+
      x*( (prop[i][j][0]*q_odd[i][j-1] +prop[i][j][1]*q_odd[i][j]*h   
	   +prop[i][j][2]*q_odd[i-1][j] +prop[i][j][3]*q_odd[i+1][j] )-prop[i][j][4]*q_odd[i][j]);
 }


i=N-1;                           /******************************************/
j=0;
  q_even[i][j]=q_odd[i][j]+
    x*( (prop[i][j][0]*q_odd[i][j]*h   +prop[i][j][1]*q_odd[i][j+1] 
	 +prop[i][j][2]*q_odd[i-1][j] +prop[i][j][3]*q_odd[i][j]*h   )-prop[i][j][4]*q_odd[i][j]);

for(j=1;j<N-1;j++){
    q_even[i][j]=q_odd[i][j]+
      x*( (prop[i][j][0]*q_odd[i][j-1] +prop[i][j][1]*q_odd[i][j+1] 
	   +prop[i][j][2]*q_odd[i-1][j] +prop[i][j][3]*q_odd[i][j]*h   )-prop[i][j][4]*q_odd[i][j]);
 }

j=N-1;
  q_even[i][N-1]=q_odd[i][j]+
    x*( (prop[i][j][0]*q_odd[i][j-1] +prop[i][j][1]*q_odd[i][j]*h   
	 +prop[i][j][2]*q_odd[i-1][j] +prop[i][j][3]*q_odd[i][j]*h  )-prop[i][j][4]*q_odd[i][j]);


q_even[Px][Py]=255.0;


/*********************************************************************************************/
/*********************************************************************************************/
 
#pragma omp parallel for private(j,intensity,log_intensity) 
for(i = 0; i < N; i++){
  // for(i = 1; i < N-1; i++){
  for(j = 0; j < N; j++){
    intensity= q_even[i][j]>255.0? 255.0: q_even[i][j];
    log_intensity=(log10(intensity)+6)*32;
    log_intensity= (log_intensity <0.0) ? 0.0:log_intensity;
    log_intensity= (log_intensity >255.0) ? 255.0:log_intensity;

    if(log_mode==0){
      C[i][j].R=intensity;
      C[i][j].G= (intensity >1e-9)? intensity:(128 - intensity*1e11);
      C[i][j].B= (intensity >1e-9)? intensity:(255 - 2.5*intensity*1e11);
    }else{
      C[i][j].R=log_intensity;
      C[i][j].G= (intensity >1e-9)? log_intensity:(128 - intensity*1e11);
      C[i][j].B= (intensity >1e-9)? log_intensity:(255 - 2.5*intensity*1e11);
      //C[i][j].B= log_intensity;
    }

    if(kabe[i][j]==0){
      C[i][j].R=255.0;
      C[i][j].G=128.0;
      C[i][j].B=64.0;
    }

    C[i][j].A=1.0;
  }
 }
 

step++;
if(step%5000==0) printf("Step=%6d W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",step,prop[1][1][0],prop[1][1][1],prop[1][1][2],prop[1][1][3]); 
//if(step%2000==0) printf("Step=%6d  q[N/2-1][N/2]=%f  q[N/2-1][N/2-1]=%f %f\n",step,q_even[N/2-1][N/2],q_even[N/2-1][N/2-1],log10(q_even[N/2-1][N/2-1]));
//printf("Step=%d  q[N-1][N-1]=%f\n",step,q_even[N-1][N-1]);
//t2 = gettimeofday_sec();
//printf("simulation : %f\n", t2 - t1);
}

oid Keyboard(unsigned char key, int x, int y)
{
  int i,j, k;

  if(key == 't')
    {
      Px = 120;
      Py = 230; 
    }
  if(key == 'a')
    {
      for(i = 0 ; i < N ; i++){
	for(j = 0 ;j < N ; j++){
	  prop[i][j][0] = 0;
	  prop[i][j][1] = 0;
	  prop[i][j][2] = 0;
	  prop[i][j][3] = 0;
	  prop[i][j][4] = 0;
	}
      }
      printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][1],prop[1][1][0],prop[1][1][2],prop[1][1][3]);
  
    }
  if (key == 27) exit(0);
  else if (key >= '1' && key <= '6')
    KEY = key; 

  if (key == 'l') log_mode=(log_mode==0)?1:0 ;
  if (key == 'w') {
    for(i=0;i<N;i++)
      for(j=0;j<N;j++){
	if(prop[i][j][0] >= 0.2 ){
	  prop[i][j][0] -= 0.2;
	  prop[i][j][4] -= 0.2;
	}
      }  
    //    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][0],prop[1][1][1],prop[1][1][2],prop[1][1][3]);
    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][1],prop[1][1][0],prop[1][1][2],prop[1][1][3]);
  }

  if (key == 'q') {
    for(i=0;i<N;i++)
      for(j=0;j<N;j++){
	if(prop[i][j][0] < 0.9 ){
	  prop[i][j][0] += 0.2;
	  prop[i][j][4] += 0.2;
	}
      }
    //    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][0],prop[1][1][1],prop[1][1][2],prop[1][1][3]);
    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][1],prop[1][1][0],prop[1][1][2],prop[1][1][3]);
  }
  if (key == 'e') {
    for(i=0;i<N;i++)
      for(j=0;j<N;j++){
	if(prop[i][j][1] >= 0.2 ){
	  prop[i][j][1] -= 0.2;
	  prop[i][j][4] -= 0.2;
	}
      }  
    //    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][0],prop[1][1][1],prop[1][1][2],prop[1][1][3]);
    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][1],prop[1][1][0],prop[1][1][2],prop[1][1][3]);
  }

  if (key == 'r') {
    for(i=0;i<N;i++)
      for(j=0;j<N;j++){
	if(prop[i][j][1] < 0.9 ){
	  prop[i][j][1] += 0.2;
	  prop[i][j][4] += 0.2;
	}
      }
    //    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][0],prop[1][1][1],prop[1][1][2],prop[1][1][3]);
    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][1],prop[1][1][0],prop[1][1][2],prop[1][1][3]);
  }

  if(key =='h'){
    printf("Type L key to toggle normal- and log-scale modes\n");
    printf("Use  Q/W keys to increase/decrease diffusion constant for E-direction\n");
    printf("Use  R/E keys to increase/decrease diffusion constant for W-direction\n");
  }
  //  simulation();
  glutPostRedisplay();
}

void Display(void)
{ 
  int i, j;
  float vx,vy;
  
  for(i = 0; i < N; i++){
    for(j = 0; j < N; j++){
      image[i][j][0] = (GLubyte)C[i][j].R;
      image[i][j][1] = (GLubyte)C[i][j].G;
      image[i][j][2] = (GLubyte)C[i][j].B;
      image[i][j][3] = (GLubyte)C[i][j].A;
    }
  }
  
  glClear(GL_COLOR_BUFFER_BIT);

  glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, N, N, GL_RGBA, GL_UNSIGNED_BYTE, image);
  
  glEnable(GL_TEXTURE_2D);
  glBegin(GL_QUADS);
  glTexCoord2f(0.0, 0.0);
  glVertex2f(-1.0, -1.0);
  glTexCoord2f(0.0, 1.0);
  glVertex2f(-1.0, 1.0);
  glTexCoord2f( 1.0,  1.0);
  glVertex2f( 1.0,  1.0);
  glTexCoord2f(1.0, 0.0);
  glVertex2f(1.0, -1.0);
  glEnd();
  glDisable(GL_TEXTURE_2D);

  
  vx=0.1*(prop[1][1][0]-prop[1][1][1]);
  vy=0.1*(prop[1][1][2]-prop[1][1][3]);

  glLineWidth(2.0);
  glBegin(GL_LINE_LOOP);
  glVertex2f(-0.9,0.9);
  glVertex2f(-0.7,0.9);
  glVertex2f(-0.7,0.7);
  glVertex2f(-0.9,0.7);
  glEnd();

  glBegin(GL_LINES);
  //     glColor3d(0.0, 0.0, 0.0);
  glVertex2f(-0.8, 0.8);
  glVertex2f(-0.8+vx, 0.8+vy);
  glEnd();

  glutSwapBuffers();

  //all_t2 = gettimeofday_sec();
  //printf("all : %f\n", all_t2 - all_t1);
}

void Idle()
{
  int i,j;
  theta += 0.005;

  for(i=0;i<N;i++)
    for(j=0;j<N;j++){
      if(kabe[i][j]==1){
	prop[i][j][2] = 0.5*(1-sin(theta));
	prop[i][j][3] = 0.1;//0.5*(1+sin(theta));
	prop[i][j][4] = prop[i][j][0]+prop[i][j][1]+prop[i][j][2]+prop[i][j][3];
      }
    }

  simulation();
  glutPostRedisplay();
}



void mouse(int button, int state, int x, int y)
{

  int i,j;

  // x,y は　左上が　(0,0) で →x  ↓がy 
  // Sim の座標系は　左下が (0,0) で，↑がx →がy 

  switch (button) {
  case GLUT_LEFT_BUTTON:
    Px = (int)(511-y);
    Py = (int)x;
    //    printf(" at (x,y)=(%d, %d) (Px,Py)=(%d %d)\n", x, y,Px,Py); 
    printf("Diffusion Source Placed  at (Px,Py)=(%d %d)\n",Px,Py); 
    printf("W=%4.2f E=%4.2f N=%4.2f S=%4.2f\n",prop[1][1][1],prop[1][1][0],prop[1][1][2],prop[1][1][3]);
    break;
  case GLUT_RIGHT_BUTTON:
    for(i=0;i<N;i++)
      for(j=0;j<N;j++){
	q_even[i][j]=0.0;
	q_odd[i][j]=0.0;
	Px=512;
	Py=0;
      }

    printf("initialized\n");
    break;
  case GLUT_MIDDLE_BUTTON:
    printf("exit!\n");
    exit(-1);
  default:
    break;
  }

  // printf(" button is ");
  
  /* switch (state) {
  case GLUT_UP:
    printf("up");
    break;
  case GLUT_DOWN:
    printf("down");
    break;
  default:
    break;
  }
  */  




}








void Init(char *progname)
{
  simulation();

  GLuint texName;
  int width = N, height = N;
  glGenTextures(1 , &texName);
  glBindTexture(GL_TEXTURE_2D , texName);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH | GLUT_ALPHA);
  glutInitWindowSize(width, height);
  glutInitWindowPosition(0, 0); 
  glutCreateWindow(progname);
  glClearColor(0.0, 0.0, 0.0, 0.0);

  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, N, N, 0, GL_RGBA, GL_UNSIGNED_BYTE, image);
}

int main(int argc, char** argv)
{

  int i,j,k;

  printf("Type L key to toggle normal- and log-scale modes\n");
  printf("Use  Q/W keys to increase/decrease diffusion constant for W-direction\n");
  printf("Use  R/E keys to increase/decrease diffusion constant for E-direction\n");
  printf("Diffusion Constants for N/S directions are automatically changed using sin function\n");

  unsigned char map[512][512][3];
  int fd;

  if((fd=open("./orig512.raw",0))<0){
    perror("open"); exit(-1);
  }
  if(read(fd, map, 512*512*3) != 512*512*3){
    perror("read\n");
    exit(-1);
  }
  close(fd);


  for(i=0;i<N;i++){
    //if(fgets(buf,N,stdin)==NULL) exit(-1);
    //    printf("%s\n",buf);
    for(j=0;j<N;j++){
      q_odd[i][j]=0.0;
      q_even[i][j]=0.0;
      prop[i][j][4]=0.0;

      if(map[i][j][0]>=128){
	kabe[N-i-1][j]=1;
	prop[N-i-1][j][0]=1.0;  //W
	prop[N-i-1][j][1]=0.2;  //E
	prop[N-i-1][j][2]=1.0;  //N
	prop[N-i-1][j][3]=1.0;  //S
      }else{
	kabe[N-i-1][j]=0;
	prop[N-i-1][j][0]=0.0;  //W
	prop[N-i-1][j][1]=0.0;  //E
	prop[N-i-1][j][2]=0.0;  //N
	prop[N-i-1][j][3]=0.0;  //S
      }

      prop[N-i-1][j][4]=prop[N-i-1][j][0]+prop[N-i-1][j][1]+prop[N-i-1][j][2]+prop[N-i-1][j][3];
    }
  }

  glutInit(&argc, argv); 
  Init(argv[0]);
  glutDisplayFunc(Display);
  glutIdleFunc(Idle);
  glutKeyboardFunc(Keyboard);
  glutMouseFunc(mouse);
  glutMainLoop(); 

  return 0;
}