// To Compile     gcc -lglut -lGL -fopenmp diffuse2d.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*/
  /*********************************************************************************************/

  
#pragma omp parallel for private(j) 
  for(i = 0; i < N; i++){

    if(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]);
  
    }else if(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]);

    }else{/*********************************************************************************************/

      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]);

      /*************************************************************************************************/
	}
  }


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


 
  
#pragma omp parallel for private(j) 
  for(i = 0; i < N; i++){     

    if(i==0){ /****************************************************************************************/

      j=0;
      //q_even[i][j]=0.0;

      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]);

    }else if(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]);

    }else{/**********************************************************************************************/

      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]);

    }


  }




  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);
}


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

  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 == 'k'){
    Px = (int)(511-y);
    Py = (int)x;
    //    printf(" at (x,y)=(%d, %d) (Px,Py)=(%d %d)\n", x, y,Px,Py); 
    printf("Put a Kabe at (Px,Py)=(%d %d)\n",Px,Py);
    
    int i,j;
    for(i = 0;i<20;i++){
      for(j= 0;j<20;j++){
	kabe[Px+i][Py+j] = 0;
	kabe[Px-i][Py-j] = 0;
	kabe[Px+i][Py-j] = 0;
	kabe[Px-i][Py+j] = 0;
      }
    }
  }
  */

  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(state){
	case GLUT_DOWN:
  	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;
	  }
	  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;
}