重力多体系の計算

重力多体系のシュミレーションをOpenGLを使ってアニメーションにしました。
計算には、ルンゲクッタを用いました。

手持ちのmacで動くように書いたので環境によっては書き直す必要があるかもしれません。
macなら二つのファイルを同じディレクトリーに置いて、

> g++ many-body-system.c -o many-body-system -O3 -framework GLUT -framework -OpenGL

でコンパイル後、

> ./many-body-system

で実行できるはずです。

many-body-system.c

#include <stdlib.h>
#include<stdio.h>
#include"MyOpengl.h"
#include<math.h>
#define N 10
#define DIM 3
#define Random (double)rand()/RAND_MAX
static double X[N][DIM];
static double U[N][DIM];
static double R[N][N];
static double M[N];
double f1(double X[N][3],double V[N][3],int i,int j)
{
  return V[j][i];
}
double f2(double X[N][3],double V[N][3],int i,int j)
{
  double F=0;
  for(int n=0;n<N;n++){
    for(int m=n;m<N;m++){
      double r=0;
      for(int k=0;k<DIM;k++){
	r+=(X[m][k]-X[n][k])*(X[m][k]-X[n][k]);
      }
      R[n][m]=R[m][n]=r;
    }
  }
  for(int n=0;n<N;n++){
    if(n!=j){
      F+=-M[n]*(X[j][i]-X[n][i])/pow(R[j][n]+0.01,1.5);
    }
  }
  return F;
}
void calc_draw(void)
{
  double K[N][4][2][DIM],X2[N][DIM],U2[N][DIM];
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      X2[j][i]=X[j][i];
      U2[j][i]=U[j][i];
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      K[j][0][0][i]=GLUT_DT*f1(X2,U2,i,j);
      K[j][0][1][i]=GLUT_DT*f2(X2,U2,i,j);
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      X2[j][i]=X[j][i]+K[j][0][0][i]/2;
      U2[j][i]=U[j][i]+K[j][0][1][i]/2;
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      K[j][1][0][i]=GLUT_DT*f1(X2,U2,i,j);
      K[j][1][1][i]=GLUT_DT*f2(X2,U2,i,j);
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      X2[j][i]=X[j][i]+K[j][1][0][i]/2;
      U2[j][i]=U[j][i]+K[j][1][1][i]/2;
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      K[j][2][0][i]=GLUT_DT*f1(X2,U2,i,j);
      K[j][2][1][i]=GLUT_DT*f2(X2,U2,i,j);
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      X2[j][i]=X[j][i]+K[j][2][0][i];
      U2[j][i]=U[j][i]+K[j][2][1][i];
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      K[j][3][0][i]=GLUT_DT*f1(X2,U2,i,j);
      K[j][3][1][i]=GLUT_DT*f2(X2,U2,i,j);
    }
  }
  for(int j=0;j<N;j++){
    for(int i=0;i<DIM;i++){
      X[j][i]=X[j][i]+(K[j][0][0][i]+2*K[j][1][0][i]+2*K[j][2][0][i]+K[j][3][0][i])/6.0;
      U[j][i]=U[j][i]+(K[j][0][1][i]+2*K[j][1][1][i]+2*K[j][2][1][i]+K[j][3][1][i])/6.0;
    }
  }
  for(int j=0;j<N;j++){
    glPushMatrix();
    glTranslated(X[j][0],X[j][1],X[j][2]);
    glColor3d(1,0,0);
    glutSolidSphere(0.5,20,20);
    glPopMatrix();
  }
}
  
      
void glutMainAction(void){
  calc_draw();
};
void init(void){
  for(int i=0;i<N;i++){
    for(int j=0;j<DIM;j++){
      X[i][j]=Random*10-5;
      U[i][j]=0;
    }
    M[i]=1;
  }
};

int main(int argc,char *argv[]){
  init();
  InitOpengl(argc,argv);
  return 0;
}

MyOpenGL.h

//仮想実験室

#include<stdio.h>        
#include<math.h>
#include<GLUT/glut.h>
#define GRID_NUMBER 100

void glutMainAction(void);//宣言のみ、からの関数、実際の挙動を書き込めばdisp()関数ないで実行される。
void init();
double GLUT_T=0,GLUT_DT=0.01;//画面が再描写されるごとにGLUT_T,にGLUT_DTだけ加えられる
int START_STOP=0;           //これが１のとき自動的に画面を再描写する
double V1=0,V2=0,V3=0.1,dV=0.1;   //パラメーターに用いることができる変数
//InitOpenglにコマンドラインの変数を入力することで初期化を行う.

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

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(1,1,1,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':init();break;
  case 'd':V1-=dV;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[])
{
  glutInit(&argc,argv);
  glutInitDisplayMode(GLUT_RGBA);
  glutCreateWindow("GLUT_LABORATORY");
  glutDisplayFunc(disp);
  glutReshapeFunc(resize);
  glutMouseFunc(mouse);
  glutSpecialFunc(special);
  glutKeyboardFunc(key);
  glutIdleFunc(idle);
  glutMainLoop();
}