2011夏合宿田中 - P.S.P.T. 計算物理研究会

タイトル　剛体の自由回転

• 発表者　田中
• 夏合宿
• 使用言語 c //動作を確認したコンパイラ名も
• 対象OS 　windows
• ライブラリ等 opengl glut
• 添付資料
• 添付ソースコード

概要 結果 考察 感想 参考文献 使用方法 ソースコード

---

概要

力が加わってない時の剛体の運動方程式は

であるが剛体系に座標を移すと

だがは剛体系の座標に付随している量なので時間的に変化はしないので

である。　よって

で　 となるような剛体系の座標を選ぶと、運動方程式は

となる。　これを展開すると

長さがa,b,cの直方体の場合、この式は

となる。

結果

夏合宿には間に合わなかったが、その後できた。

考察

感想

単純な形の剛体については実現できたと思う　今後はコマ、地球の摂動なども実現できるようにしたい。

参考文献

ゴールドステイン　古典力学
opengl入門　など

---

使用方法

ソースコード

#include <stdlib.h>
#include <GL/glut.h>
#include <math.h>
/* 角速度　*/
double wx = 0.0; /* x軸 */
double wy = 3.0; /* y軸 */
double wz = 1.0; /* z軸 */
const double dt=0.01;
const double pi=3.1415926535 ;
//aは長方体のx軸方向の長さ　bはy軸　cはz軸である
const double a=0.3,b=2.3,c=0.3;
double t=0, v,r,s;
 
double f1(double t,double wy,double wz);
double f2(double t,double wx,double wz);
double f3(double t,double wx,double wy);
 
 
 GLdouble vertex[][3] = {
  { -0.5, -0.5, -0.5 },
  { 0.5, -0.5, -0.5 },
  { 0.5, 0.5, -0.5 },
  { -0.5, 0.5, -0.5 },
  { -0.5, -0.5, 0.5 },
  { 0.5, -0.5, 0.5 },
  { 0.5, 0.5, 0.5 },
  { -0.5, 0.5, 0.5 }
 
};
 
int face[][4] = {
  { 0, 1, 2, 3 },
  { 1, 5, 6, 2 },
  { 5, 4, 7, 6 },
  { 4, 0, 3, 7 },
  { 4, 5, 1, 0 },
  { 3, 2, 6, 7 }
};
 
GLdouble normal[][3] = {
  { 0.0, 0.0,-1.0 },
  { 1.0, 0.0, 0.0 },
  { 0.0, 0.0, 1.0 },
  {-1.0, 0.0, 0.0 },
  { 0.0,-1.0, 0.0 },
  { 0.0, 1.0, 0.0 }
};
 
GLfloat light0pos[] = { 0.0, 3.0, 5.0, 1.0 };
GLfloat light1pos[] = { 5.0, 3.0, 0.0, 1.0 };
 
GLfloat green[] = { 0.0, 1.0, 0.0, 1.0 };
GLfloat yerow[] = { 0.0, 2.0, 0.2, 1.5 };
 
  void drawaxis3D(double w){
  glBegin(GL_LINES);
  glColor3ub(192,0,0);
  glVertex3d(-w, 0, 0);
  glVertex3d( w, 0, 0);
  glColor3ub(192,0,0);
  glVertex3d( 0,-w, 0);
  glVertex3d( 0, w, 0);
  glColor3ub(192,0,0);
  glVertex3d( 0, 0,-w);
  glVertex3d( 0, 0, w);  
  glEnd();
  }
 
 
void idle(void)
{ 
 
  glutPostRedisplay();
}
 
void Calculate()
{     
  double k0[3],k1[3],k2[3],k3[3],k4[3];
  t=t+dt;
  k0[0]=dt*f1(t,wy,wz);
  k0[1]=dt*f2(t,wx,wz);
  k0[2]=dt*f3(t,wx,wy);
 
  k1[0]=dt*f1(t+dt/2.0,wy+k0[1]/2.0,wz+k0[2]/2.0);
  k1[1]=dt*f2(t+dt/2.0,wx+k0[0]/2.0,wz+k0[2]/2.0);
  k1[2]=dt*f3(t+dt/2.0,wx+k0[0]/2.0,wy+k0[1]/2.0);
 
  k2[0]=dt*f1(t+dt/2.0,wy+k1[1]/2.0,wz+k1[2]/2.0);
  k2[1]=dt*f2(t+dt/2.0,wx+k1[0]/2.0,wz+k1[2]/2.0);
  k2[2]=dt*f3(t+dt/2.0,wx+k1[0]/2.0,wy+k1[1]/2.0);
 
  k3[0]=dt*f1(t+dt,wy+k2[1],wz+k2[2]);
  k3[1]=dt*f2(t+dt,wx+k2[0],wz+k2[2]);
  k3[2]=dt*f3(t+dt,wx+k2[0],wy+k2[1]);
 
  wx=wx+(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
  wy=wy+(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
  wz=wz+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
  v = sqrt(pow(wx,2)+pow(wy,2)+pow(wz,2));
  r=t*v;
  s=r*180/pi;
}
 
void display(void)
{
  int i,j;
 
  glClear(GL_COLOR_BUFFER_BIT| GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();
  /* 視点位置と視線方向 */
  gluLookAt(3.0, 4.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
 
  drawaxis3D(10);
 
  /* 光源の位置設定 */
  glLightfv(GL_LIGHT0, GL_POSITION, light0pos);
  glLightfv(GL_LIGHT1, GL_POSITION, light1pos);
 
  Calculate();
 
  /* 図形の回転 */
  glRotated((double)s, wx, wy, -wz);
 
 
  glScaled(a,b,c);
 
  /* 図形の描画 */
  glColor3d(0.0, 0.0, 0.0);
  glBegin(GL_QUADS);
  for (j = 0; j < 6; ++j) {
  glNormal3dv(normal[j]);
    for (i = 0; i < 4; ++i) {
      glVertex3dv(vertex[face[j][i]]);
    }
  }
 
  glEnd();
  glPopMatrix();
 
  glutSwapBuffers();
 
  /* 一周回ったら回転角を 0 に戻す */
  if (r >= 360) r = 0;
}
 
 
void resize(int w, int h)
{
  glViewport(0, 0, w, h);
 
  /* 透視変換行列の設定 */
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(30.0, (double)w / (double)h, 1.0, 100.0);
 
  /* モデルビュー変換行列の設定 */
  glMatrixMode(GL_MODELVIEW);
}
 
 void mouse(int button, int state, int x, int y)
{
  switch (button) {
  case GLUT_LEFT_BUTTON:  
    if (state == GLUT_DOWN) {
      /* アニメーション開始 */
      glutIdleFunc(idle);
    }
     else  {
      /* アニメーション停止 */
      glutIdleFunc(0);
    }
    break;
  case GLUT_RIGHT_BUTTON:
    if (state == GLUT_DOWN) {
      /* コマ送り (1ステップだけ進める) */
        glutPostRedisplay();
    }
    break;
  default:
    break;
  }
}  
 
void keyboard(unsigned char key, int x, int y)
{
  switch (key) {
  case 'q':
  case 'Q':
  case '\033':  /* '\033' は ESC の ASCII コード */
    exit(0);
  case 'z' :
    glutIdleFunc(idle);
    break;
    case 'x' :
    glutIdleFunc(0);
    break;
  default:
    break;
  }
}
 
void init(void)
{
  glClearColor(1.0, 1.0, 1.0, 1.0);
  glEnable(GL_DEPTH_TEST);
  glEnable(GL_CULL_FACE);
  glCullFace(GL_FRONT);
 
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);
  glEnable(GL_LIGHT1);
  glLightfv(GL_LIGHT1, GL_DIFFUSE, green);
  glLightfv(GL_LIGHT1, GL_SPECULAR, green);
  glLightfv(GL_LIGHT0, GL_DIFFUSE, yerow);
  glLightfv(GL_LIGHT0, GL_SPECULAR, yerow);
}
 
int main(int argc, char *argv[])
{
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_RGBA| GLUT_DOUBLE| GLUT_DEPTH);
  glutCreateWindow(argv[0]);
  glutDisplayFunc(display);
  glutReshapeFunc(resize);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keyboard);
  init();
  glutMainLoop();
  return 0;
}
 
double f1(double t,double wy,double wz)
{
        double r;
 
        r=((pow(c,2)-pow(b,2))/(pow(b,2)+pow(c,2)))*wy*wz;
 
        return r;
}
 
double f2(double t,double wx,double wz)
{
        double r;
 
        r=((pow(a,2)-pow(c,2))/(pow(a,2)+pow(c,2)))*wz*wx;
 
        return r;
}
 
double f3(double t, double wx,double wy)
{
        double r;
 
        r=((pow(b,2)-pow(a,2))/(pow(a,2)+pow(b,2)))*wx*wy;
 
        return r;
}