3次元のイジングモデル

3次元のイジングモデルのシュミレーションを行ないました。
前回同様、OpenGLで可視化しています。
macなら、

$ g++ ising.cpp -o ising -O2 -famework GLUT -framework OpenGL

でコンパイル、

$ ./ising

で実行できます。
Linuxの場合は、freeglutがインストールされている必要があります。
MyOpengl.hの #include を　#includeに書き換えて

$ g++ ising.cpp -o ising -lglut -lGLU

でコンパイル、

$ ./ising

で実行できます。
ising.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <ctime>
#include "MyOpengl.h"
 
using namespace std;
 
typedef vector <int> Sv;
typedef vector <Sv> Sm;
typedef vector <Sm> St;
const int L = 20;
const int M = 100;
const double J = 1.0;
const double B = 0.0;
double beta = 0;
St spin(L, Sm(L, Sv(L, 1)));
double ran() {
  return (double)rand()/((double)RAND_MAX+1);
}
double dE(int i, int j, int k, int Pre) {
  return -J *( Pre - spin[i][j][k])* ( (i == 0 ? 0 : spin[i-1][j][k])
				       + (i == L-1 ? 0 : spin[i+1][j][k])
				       +(j == 0 ? 0 : spin[i][j-1][k])
				       +(j == L-1 ? : spin[i][j+1][k])
				       +(k == 0 ? 0 : spin[i][j][k-1])
				       +(k == L-1 ? 0 : spin[i][j][k+1])
				       )
    -B*(Pre - spin[i][j][k]);
}    
void sweep_1(int i, int j, int k) {
  double r = ran();
  int Pre = -spin[i][j][k];
  double difE = dE(i, j, k, Pre);
  if(difE < 0) {
    spin[i][j][k] = Pre;
  } else if(r < exp(-beta*difE) ) {
    spin[i][j][k] = Pre;
  }
}
void sweep_all() {
  for(int i = 0; i < L ; i++) {
    for(int j = 0; j < L; j++) {
      for(int k =0; k < L; k++) {
	sweep_1(i, j, k);
      }
    }
  }
}
double cl(int n) {
  double sum = 0;
  for(int i = 0; i < L; i++ ) {
    for(int j = 0; j < L; j++ ) {
      for(int k = 0; k < L; k++) {
	sum += (double) spin[i][j][k] * spin[i][j][k];
      }
    }
  }
  return sum;
}
double magnet() {
  double sum = 0;
  for(int n = 0 ; n < M; n++) {
    for(int i = 0; i < L; i++) {
      for(int j = 0; j < L; j++) {
	for(int k = 0; k < L; k++) {
	  sum += spin[i][j][k];
	}
      }
    }
    sweep_all();
  }
  return sum / (M * L * L * L);
}
void sphere(int x, int y,int z, bool direction) {


  glColor4f(direction, 1-direction, 1, 0.2);
  glPushMatrix();
  glTranslated(x, y, z);
  glutSolidSphere(0.3, 10, 10);
  glPopMatrix();

}
void action() {
  glDisable(GL_DEPTH_TEST);
  glEnable(GL_BLEND);
  for(int i = 0; i < L; i++) {
    for(int j = 0; j < L; j++) {
      for(int k = 0; k < L; k++) {
	sphere(i, j, k, 0.5 + 0.5 *spin[i][j][k]);
      }
    }
  }
  glDisable(GL_BLEND);
  glEnable(GL_DEPTH_TEST);
  beta = V2;
  cout << beta << endl;
  sweep_all();
}
int main(int argc, char **argv) {
  
  srand( (unsigned)time(NULL));
  InitOpengl(argc, argv, action);

  return 0;
}

MyOpengl.h

//仮想実験室
#include <stdlib.h>
#include <time.h>
#include<stdio.h>        
#include<math.h>
#include <GLUT/glut.h>
#define GRID_NUMBER 100


typedef void (*Action)(void);
Action glutMainAction;

//画面が再描写されるごとにGLUT_T,にGLUT_DTだけ加えられる
double GLUT_T=0,GLUT_DT=0.01;

//これが１のとき自動的に画面を再描写する
int START_STOP=0;          

//パラメーターに用いることができる変数
double V1=0,V2=0,V3=0.1,dV=0.01;

//InitOpenglにコマンドラインの変数を入力することで初期化を行う.

void resize(int w,int h)
{
  glViewport(0,0,w,h);
  glLoadIdentity();
  gluPerspective(60,(double)w/h,1,100);
  glTranslated(-5, -5, -40);
}

void field(void)
{  
  glColor3d(0.8,0.8,0.8);
  int i,j;
  glBegin(GL_LINES);
  for(i=-GRID_NUMBER;i<GRID_NUMBER;i++){
    glVertex2d(-GRID_NUMBER,i);
    glVertex2d(GRID_NUMBER,i);
    glVertex2d(i,GRID_NUMBER);
    glVertex2d(i,-GRID_NUMBER);
  }
  glEnd();
}
void disp(void)
{   
  glClearColor(0, 0, 0, 0);
  glClear(GL_COLOR_BUFFER_BIT);
  //field();
  glutMainAction();
  glFlush();

}
void mouse(int button,int state,int x,int y)
{
  switch(button){
  case GLUT_LEFT_BUTTON:
    if(state==GLUT_DOWN){
      glRotated(-10,1,0,0);
      break;
    }
  case GLUT_RIGHT_BUTTON:
    if(state==GLUT_DOWN){
      glRotated(10,1,0,0);
    }
    break;

  default:break;
  }
  glutPostRedisplay();

}
void special(int key,int x,int y)
{
  switch (key){
  case GLUT_KEY_UP:
    glTranslated(0,-0.1,0);
    break;

  case GLUT_KEY_DOWN:
    glTranslated(0,0.1,0);
    break;
  case GLUT_KEY_LEFT:
    glTranslated(0.1,0,0);
    break;

  case GLUT_KEY_RIGHT:
    glTranslated(-0.1,0,0);
    break;

  default:break;
  }
  glutPostRedisplay();
} 
void key(unsigned char key,int x,int y)
{   int i;
  switch (key){
  case 'z':glRotated(1,0,0,1);
    break;
  case 'x':glRotated(-1,0,0,1);
    break;
  case 'q':exit(0);
    break;
  case 't':if(START_STOP==1){
      START_STOP=0;
    } else if(START_STOP==0){
      START_STOP=1;
    }
    break;
  case 'a':glTranslated(0,0,0.2);
    break;
  case 's':glTranslated(0,0,-0.2);
    break;
  case 'c':glRotated(1,1,0,0);
    break;
  case 'd':glRotated(-1,1,0,0);
    break;
  case 'v':V2+=dV;break;
  case 'f':V2-=dV;break;
  case 'b':V3+=dV;break;
  case 'g':V3-=dV;break;
  default:break;
  }
  glutPostRedisplay();
}
void idle(void)
{
  if(START_STOP==1){
    GLUT_T+=GLUT_DT;
    glutPostRedisplay();
  }
}
void InitOpengl(int argc, char *argv[],Action action)
{
  
  glutMainAction = action;
  glutInit(&argc,argv);
  glutInitDisplayMode(GLUT_RGBA);
  glutCreateWindow("GLUT_LABORATORY");
  glutDisplayFunc(disp);
  glutReshapeFunc(resize);
  glutMouseFunc(mouse);
  glutSpecialFunc(special);
  glutKeyboardFunc(key);
  glutIdleFunc(idle);
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  glutMainLoop();
}