HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 仮想物理実験室

VisualC++ と OpenGL を利用した仮想物理実験室
【0-2-3】仮想物理実験室の構築 （ver1.2）

文責：遠藤 理平 （2010年9月14日） カテゴリ：仮想物理実験室(325)

「【0-2-2】仮想物理実験室の構築 （ver1.1）」をベースに、仮想物理実験室の構築 （ver1.2）を構築します。

仮想物理実験室(ver 1.2)

「【０日目】仮想物理実験室の構築 （ver1.0）」と同様に、「gl_screenshot.h」と「gl_material.h」を予め準備してください。 

//////////////////////////////////////////////////////////////////////////
// 仮想物理実験室(ver 1.2)
// 1.0+ばね関数+床
//////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <direct.h>
#include <time.h>
#include <GL/glut.h>
#include <GL/gl_material.h>
#include <GL/gl_screenshot.h>
using namespace std;
double PI = acos(-1.0);

//////////////////////////////////////////////////////////////////////////
// 変数の定義
//////////////////////////////////////////////////////////////////////////
#define _BITMAP 1//アニメーション作成用ビットマップの保存 0:しない　1:する

//--------------------------------------------------------
// 仮想物理実験室変数の定義
//--------------------------------------------------------
double t = 0.0;  //時刻
double dt= 0.01; //時間刻み
int tn = 0;      //ステップ数

//箱の位置と一辺の長さ
double box_x = 30.0, box_y = 0.0, box_z = 20.0, box_l = 10.0;
//ボールの位置と半径
double  ball1_x =-20.0, ball1_y = 15.0, ball1_z = 30.0, ball1_r =8.0;
double  ball2_x =-10.0, ball2_y = -10.0,  ball2_z = 10.0, ball2_r =4.0;

//--------------------------------------------------------
// OpenGL用変数定義
//--------------------------------------------------------
//////////////////////////////////////////
// ウィンドウ生成用
int WindowPositionX = 200;                         //生成するウィンドウ位置のX座標
int WindowPositionY = 200;                         //生成するウィンドウ位置のY座標
int WindowWidth     = 512;                         //生成するウィンドウの幅
int WindowHeight    = 512;                         //生成するウィンドウの高さ
char WindowTitle[]  = "仮想物理実験室(ver 1.1)";   //ウィンドウのタイトル
//////////////////////////////////////////
// アニメーション用　ビットマップ保存
gl_screenshot gs; //「gl_screenshot」クラスのインスタンス「gs」を宣言
//////////////////////////////////////////
// 床描画用
static GLfloat floor_planar[4];
static GLfloat floor_s = 50.0f; //床の広さ（1辺=2.0*floor_s）
static GLfloat pM[16];
static GLfloat LightPosition[4] = { -10, -20, 70, 1 }; //光源の位置
typedef struct _QUADS_VERTEX{
 GLfloat v0[3];
 GLfloat v1[3];
 GLfloat v2[3];
 GLfloat v3[3];
}QUADS_VERTEX;
static QUADS_VERTEX floor_v = {
 {  floor_s,  floor_s, 0.0f },
 { -floor_s,  floor_s, 0.0f },
 { -floor_s, -floor_s, 0.0f },
 {  floor_s, -floor_s, 0.0f },
};
void findPlane(GLfloat plane[4], GLfloat v0[3], GLfloat v1[3], GLfloat v2[3]);
void DrawFloor(bool bTexture);

//////////////////////////////////////////
// マウスドラッグ用
int cx, cy;                                             // ドラッグ開始位置
double sx, sy;                                          // マウスの絶対位置→ウィンドウ内での相対位置の換算係数
double cq[4] = { 1.0, 0.0, 0.0, 0.0 };                  // 回転の初期値 (クォータニオン)
double tq[4];                                           // ドラッグ中の回転 (クォータニオン)
double rt[16];                                          // 回転の変換行列
unsigned int listNumber;
float camera_z_pos =200.0;

//////////////////////////////////////////
// 文字描画用
int text_list;
char t_char[20];
char t_char2[20];
void DRAW_STRING(int x, int y, char *string, void *font = GLUT_BITMAP_TIMES_ROMAN_24);
void DISPLAY_TEXT(int x, int y, char *string);

//--------------------------------------------------------
// 回転用
//--------------------------------------------------------
void qmul(double r[], const double p[], const double q[]);
void qrot(double r[], double q[]);

//--------------------------------------------------------
// 関数のプロトタイプ
//--------------------------------------------------------
//////////////////////////////////////////
// メイン関数用
void Initialize(void);
void Display(void);
void Resize(int w, int h); 
void Idle(void);
void Keyboard(unsigned char key, int x, int y);
void mouse_motion(int x, int y);
void mouse_on(int button, int state, int x, int y);
void mouse_wheel(float);
//////////////////////////////////////////
// 描画用
void Calculate(void);      //計算
void DrawStructure(void);  //描画
void Ground(void);         //大地の描画
void Cuboid(double, double, double); //直方体(tate, yoko, takasa)
//////////////////////////////////////////
// 影描画用
void shadowMatrix(GLfloat *m, GLfloat plane[4],  GLfloat light[4]);
void DrawShadow(void);
//////////////////////////////////////////
// ばね描画用
void drowSolidSpring(double x0, double y0, double z0, double x1, double y1, double z1);//
void skSolidSpring(int Nm, int Ns, int Np, double radius, double ratio, double length);

//--------------------------------------------------------
// メイン関数
//--------------------------------------------------------
int main(int argc, char *argv[]){
 //SetUp();
 srand((unsigned)time(NULL));
#if _BITMAP
 _mkdir("bitmap"); //bmpファイル保存用のフォルダの作成
#endif
 glutInit(&argc, argv);                                     //環境の初期化
 glutInitWindowPosition(WindowPositionX, WindowPositionY);  //ウィンドウの位置の指定
 glutInitWindowSize(WindowWidth, WindowHeight);             //ウィンドウサイズの指定
 glutInitDisplayMode(GLUT_RGBA | GLUT_DEPTH | GLUT_DOUBLE | GLUT_STENCIL); //ディスプレイモードの指定
 glutCreateWindow(WindowTitle);                             //ウィンドウの作成
 glutDisplayFunc(Display);                                  //描画時に呼び出される関数を指定する（関数名：Display）
 glutReshapeFunc(Resize);                                   //リサイズにより呼び出される関数
 glutMouseFunc(mouse_on);                                   //マウスクリック時に呼び出される関数
 glutMotionFunc(mouse_motion);                              //マウスドラッグ解除時に呼び出される関数
 glutKeyboardFunc(Keyboard);                                //キーボード入力時に呼び出される関数を指定する（関数名：Keyboard）
 glutIdleFunc(Idle);                                        //プログラムアイドル状態時に呼び出される関数
 Initialize();                                              //初期設定の関数を呼び出す
 glutMainLoop();
 return 0;
}
//--------------------------------------------------------
// 初期設定の関数
//--------------------------------------------------------
void Initialize(void){
 glClearColor(1.0, 1.0, 1.0, 1.0); //背景色
 glEnable( GL_DEPTH_TEST ); //デプスバッファを使用：glutInitDisplayMode() で GLUT_DEPTH を指定する
 glDepthFunc( GL_LEQUAL );
 glClearDepth( 1.0 );
 //////////////////////////////////////////
 // ディスプレイリストを作成
 listNumber = glGenLists(1);
 glNewList( listNumber, GL_COMPILE );
 glEndList();
 //////////////////////////////////////////
 // マウスポインタ位置のウィンドウ内の相対的位置への換算用
 sx = 1.0 / (double)WindowWidth;
 sy = 1.0 / (double)WindowHeight;
 // 回転行列の初期化
 qrot(rt, cq);
 //////////////////////////////////////////
 // 床
 findPlane( floor_planar,
       floor_v.v0,
       floor_v.v1,
       floor_v.v2 );
 //////////////////////////////////////////
 //透視変換行列の設定
 glMatrixMode(GL_PROJECTION);//行列モードの設定（GL_PROJECTION : 透視変換行列の設定、GL_MODELVIEW：モデルビュー変換行列）
 glLoadIdentity();//行列の初期化
 gluPerspective(30.0, (double)WindowWidth/(double)WindowHeight, 0.1, 1000.0); //透視投影法の視体積gluPerspactive(th, w/h, near, far);

}
//--------------------------------------------------------
// 描画の関数
//--------------------------------------------------------
void Display(void) {
 glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT);

 //////////////////////////////////////////
 //モデルビュー変換行列の設定
 glMatrixMode(GL_MODELVIEW);//行列モードの設定（GL_PROJECTION : 透視変換行列の設定、GL_MODELVIEW：モデルビュー変換行列）
 glLoadIdentity();//行列の初期化
 glViewport(0, 0, WindowWidth, WindowHeight);

 //////////////////////////////////////////
 //視点の設定
 gluLookAt(
      0.0,  0.0,  camera_z_pos, // 視点の位置x,y,z;
      0.0,  0.0,  0.0,   // 視界の中心位置の参照点座標x,y,z
      0.0,  1.0,  0.0 ) ;  //視界の上方向のベクトルx,y,z

 //////////////////////////////////////////
 //ステンシルバッファクリア値の設定
 glClearStencil( 0 );
 glCullFace( GL_BACK );
 glEnable( GL_CULL_FACE );
 glEnable( GL_AUTO_NORMAL );
 glEnable( GL_NORMALIZE );

 //////////////////////////////////////////
 // 平面射影行列の算出
 shadowMatrix(pM,floor_planar,LightPosition);

 //////////////////////////////////////////
 // 回転
 glMultMatrixd(rt);

 //////////////////////////////////////////
 // 光源ON
 glEnable( GL_LIGHTING );
 glEnable( GL_LIGHT0 );
 glLightfv( GL_LIGHT0,GL_POSITION,LightPosition );

 //////////////////////////////////////////
 // 描画
 glPushMatrix();
   Calculate();     //計算
   DrawStructure(); //物体
   DrawShadow();    //影
 glPopMatrix();

 //////////////////////////////////////////
 //陰影OFF
 glDisable(GL_AUTO_NORMAL);
 glDisable(GL_NORMALIZE);
 glDisable(GL_LIGHTING);

 //////////////////////////////////////////
 //文字の描画
 strcpy_s(t_char2, "t = ");
 sprintf_s(t_char, "%5.2f", t);
 strcat_s(t_char2, t_char);
 DISPLAY_TEXT(5, 95, t_char2  );

 //////////////////////////////////////////
 //ビットマップの保存
 #if _BITMAP
   ostringstream fname;
   int tt = tn +10000;
   fname  << "bitmap/" << tt << ".bmp" ;//出力ファイル名
   string name = fname.str();
   gs.screenshot(name.c_str(), 24);
 #endif

 glutSwapBuffers(); //glutInitDisplayMode(GLUT_DOUBLE)でダブルバッファリングを利用可
}
//--------------------------------------------------------
// 計算と物体の描画
//--------------------------------------------------------
void Calculate(){
  t = dt * double(tn);
  tn++;
}
void DrawStructure(){
  glPushMatrix();
  glMaterialfv(GL_FRONT, GL_AMBIENT, ms_ruby.ambient);
  glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_ruby.diffuse);
  glMaterialfv(GL_FRONT, GL_SPECULAR, ms_ruby.specular);
  glMaterialfv(GL_FRONT, GL_SHININESS, &ms_ruby.shininess);
  glTranslated(0.0, 0.0, 20.0);//平行移動値の設定
  glutSolidCube(10.0);//引数：(一辺の長さ)
  glPopMatrix();
//ボール
  glPushMatrix();
  glMaterialfv(GL_FRONT, GL_AMBIENT, ms_silver.ambient);
  glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_silver.diffuse);
  glMaterialfv(GL_FRONT, GL_SPECULAR, ms_silver.specular);
  glMaterialfv(GL_FRONT, GL_SHININESS, &ms_silver.shininess);
  glTranslated(ball1_x, ball1_y, ball1_z); //平行移動値の設定
  glutSolidSphere(ball1_r, 20, 20);        //引数：(半径, Z軸まわりの分割数, Z軸に沿った分割数)
  glPopMatrix();
  glPushMatrix();
  glTranslated(ball2_x, ball2_y, ball2_z); //平行移動値の設定
  glutSolidSphere(ball2_r, 20, 20);        //引数：(半径, Z軸まわりの分割数, Z軸に沿った分割数)
  glPopMatrix();
  //バネ
  glPushMatrix();
  glMaterialfv(GL_FRONT, GL_AMBIENT, ms_yellow_rubber.ambient);
  glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_yellow_rubber.diffuse);
  glMaterialfv(GL_FRONT, GL_SPECULAR, ms_yellow_rubber.specular);
  glMaterialfv(GL_FRONT, GL_SHININESS, &ms_yellow_rubber.shininess);  
  drowSolidSpring(ball1_x, ball1_y, ball1_z, ball2_x, ball2_y, ball2_z);
  glPopMatrix();
  //直方体
  glPushMatrix();
    glMaterialfv(GL_FRONT, GL_AMBIENT, ms_ruby.ambient);
    glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_ruby.diffuse);
    glMaterialfv(GL_FRONT, GL_SPECULAR, ms_ruby.specular);
    glMaterialfv(GL_FRONT, GL_SHININESS, &ms_ruby.shininess);
    glTranslated( 30.0, 30.0, 0.0);//平行移動値の設定
    Cuboid(20.0, 10.0, 5.0);
  glPopMatrix();
}

//--------------------------------------------------------
// アイドル時に呼び出される関数
//--------------------------------------------------------
void Idle(){
 glutPostRedisplay(); //glutDisplayFunc()を１回実行する
}
//---------------------------------------------------------------------------------------------------
// サイズ変更
//---------------------------------------------------------------------------------------------------
void Resize(int x, int y)
{
	WindowWidth = x;
	WindowHeight = y;
	if ( WindowWidth < 1 ) WindowWidth = 1;
	if ( WindowHeight < 1 ) WindowHeight = 1;
}

//--------------------------------------------------------
// キーボード入力時に呼び出される関数
//--------------------------------------------------------
void Keyboard(unsigned char key, int x, int y){
 switch ( key )
 {
 case 'i':
    camera_z_pos -= 10.0;
    break;
 case 'o':
    camera_z_pos += 10.0;
    break;

 default:
   break;
 }
 cout << camera_z_pos << endl;
}
//--------------------------------------------------------
// マウスドラッグ時
//--------------------------------------------------------
void mouse_motion(int x, int y){
 double dx, dy, a;

 // マウスポインタの位置のドラッグ開始位置からの変位
 dx = (x - cx) * sx;
 dy = (y - cy) * sy;

 // マウスポインタの位置のドラッグ開始位置からの距離
 a = sqrt(dx * dx + dy * dy);

 if( a != 0.0 )
 {
   // マウスのドラッグに伴う回転のクォータニオン dq を求める
   double ar = a * 2.0 * PI * 0.5;
   double as = sin(ar) / a;
   double dq[4] = { cos(ar), dy * as, dx * as, 0.0 };

   // 回転の初期値 cq に dq を掛けて回転を合成
   qmul(tq, dq, cq);

   // クォータニオンから回転の変換行列を求める
   qrot(rt, tq);
 }
}
//--------------------------------------------------------
// マウスクリック時
//--------------------------------------------------------
void mouse_on(int button, int state, int x, int y){
 switch (button) {
 case 0:
   switch (state) {
   case 0:
     // ドラッグ開始点を記録
     cx = x;
     cy = y;
     break;
   case 1:
     // 回転の保存
     cq[0] = tq[0];
     cq[1] = tq[1];
     cq[2] = tq[2];
     cq[3] = tq[3];
     break;
   default:
     break;
   }
   break;
 default:
   break;
 }
 cout << x << " " << y<<endl;
}
//--------------------------------------------------------
// 床平面の方程式と行列の計算
//--------------------------------------------------------
void findPlane(
 GLfloat plane[4],  // 作成する平面方程式の係数
 GLfloat v0[3],    // 頂点１
 GLfloat v1[3],    // 頂点２
 GLfloat v2[3])    // 頂点３
{
 GLfloat vec0[3], vec1[3];

 // Need 2 vectors to find cross product.
 vec0[0] = v1[0] - v0[0];
 vec0[1] = v1[1] - v0[1];
 vec0[2] = v1[2] - v0[2];

 vec1[0] = v2[0] - v0[0];
 vec1[1] = v2[1] - v0[1];
 vec1[2] = v2[2] - v0[2];

 // find cross product to get A, B, and C of plane equation
 plane[0] = vec0[1] * vec1[2] - vec0[2] * vec1[1];
 plane[1] = -(vec0[0] * vec1[2] - vec0[2] * vec1[0]);
 plane[2] = vec0[0] * vec1[1] - vec0[1] * vec1[0];

 plane[3] = -(plane[0] * v0[0] + plane[1] * v0[1] + plane[2] * v0[2]);
}
void shadowMatrix(
 GLfloat *m,      // 作成する行列のポインタ
 GLfloat plane[4],  // 射影する表面の平面方程式の係数
 GLfloat light[4])  // 光源の同時座標値
{
 GLfloat dot;

 // Find dot product between light position vector and ground plane normal.
 dot = plane[0] * light[0] +
     plane[1] * light[1] +
     plane[2] * light[2] +
     plane[3] * light[3];

 m[0]  = dot - light[0] * plane[0];
 m[4]  = 0.f - light[0] * plane[1];
 m[8]  = 0.f - light[0] * plane[2];
 m[12] = 0.f - light[0] * plane[3];

 m[1]  = 0.f - light[1] * plane[0];
 m[5]  = dot - light[1] * plane[1];
 m[9]  = 0.f - light[1] * plane[2];
 m[13] = 0.f - light[1] * plane[3];

 m[2]  = 0.f - light[2] * plane[0];
 m[6]  = 0.f - light[2] * plane[1];
 m[10] = dot - light[2] * plane[2];
 m[14] = 0.f - light[2] * plane[3];

 m[3]  = 0.f - light[3] * plane[0];
 m[7]  = 0.f - light[3] * plane[1];
 m[11] = 0.f - light[3] * plane[2];
 m[15] = dot - light[3] * plane[3];
}
//--------------------------------------------------------
// 床の描画と影の描画
//--------------------------------------------------------
void DrawFloor(bool bTexture){
  double pich =10.0;
 if( bTexture ){
 // 床にテクスチャを使う時はココで設定する
 //  glBindTexture( GL_TEXTURE_2D, );

   glDisable(GL_LIGHTING);
    glBegin(GL_QUADS);
    for (int j = -20; j < 20; j++) {
      for (int i = -20; i < 20; i++) {
        if((i+j)%2==0) glColor4f( 0.6, 0.6, 0.6, 1.0);
        else glColor4f( 0.3, 0.3, 0.3, 1.0);
        glVertex3d((GLdouble)(i-0.5)*pich, (GLdouble)(j-0.5)*pich, 0);
        glVertex3d((GLdouble)(i+0.5)*pich, (GLdouble)(j-0.5)*pich, 0);
        glVertex3d((GLdouble)(i+0.5)*pich, (GLdouble)(j+0.5)*pich, 0);
        glVertex3d((GLdouble)(i-0.5)*pich, (GLdouble)(j+0.5)*pich, 0);
      }
    }
    glEnd();
   glEnable(GL_LIGHTING);
 }else{
   glDisable(GL_LIGHTING);
    glBegin(GL_QUADS);
    for (int j = -20; j < 20; j++) {
      for (int i = -20; i < 20; i++) {
        glVertex3d((GLdouble)(i-0.5)*pich, (GLdouble)(j-0.5)*pich, 0);
        glVertex3d((GLdouble)(i+0.5)*pich, (GLdouble)(j-0.5)*pich, 0);
        glVertex3d((GLdouble)(i+0.5)*pich, (GLdouble)(j+0.5)*pich, 0);
        glVertex3d((GLdouble)(i-0.5)*pich, (GLdouble)(j+0.5)*pich, 0);
      }
    }
    glEnd();
   glEnable(GL_LIGHTING);
 }
}
void DrawShadow(void){
 /////////////////////////////////////////////
 //床のステンシルを付ける
 glEnable(GL_STENCIL_TEST);
 glStencilFunc( GL_ALWAYS, 1, ~0);
 //これから描画するもののステンシル値にすべて１タグをつける
 glStencilOp(GL_KEEP,GL_KEEP ,GL_REPLACE);
 DrawFloor( true );//床の描画

 /////////////////////////////////////////////
 //カラー・デプスバッファマスクをセットする
 //これで以下の内容のピクセルの色の値は、書き込まれない。
 glColorMask(0,0,0,0);
 glDepthMask(0);
 /////////////////////////////////////////////
 //床にオブジェクトの影のステンシルを付ける
 glEnable(GL_STENCIL_TEST);
 glStencilFunc( GL_EQUAL, 1, ~0);
 //これから描画するもののステンシル値にすべて１タグをつける
 glStencilOp(GL_KEEP,GL_KEEP ,GL_INCR);
 glDisable(GL_DEPTH_TEST);
 glPushMatrix();
   glMultMatrixf(pM);
   DrawStructure();
 glPopMatrix();
 glEnable(GL_DEPTH_TEST);

 /////////////////////////////////////////////
 //ビットマスクを解除
 glColorMask(1,1,1,1);
 glDepthMask(1);

 /////////////////////////////////////////////
 //影をつける
 glStencilFunc( GL_EQUAL, 2, ~0 );
 glEnable(GL_BLEND);
   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
   glColor4f(0.1f, 0.1f, 0.1f, 0.5f);
   glDisable(GL_DEPTH_TEST);
     DrawFloor( false );//床の描画
   glEnable(GL_DEPTH_TEST);
   glDisable(GL_BLEND);
 glDisable(GL_STENCIL_TEST);

}

//--------------------------------------------------------
// 文字描画
//--------------------------------------------------------
void DISPLAY_TEXT(int x, int y, char *string){
glDisable(GL_LIGHTING);
glDisable(GL_LIGHT0);

 glPushAttrib(GL_ENABLE_BIT);
 glMatrixMode(GL_PROJECTION);
 glPushMatrix();
 glLoadIdentity();
 gluOrtho2D(0, 100, 0, 100);
 glMatrixMode(GL_MODELVIEW);
 glPushMatrix();
 glLoadIdentity();
 glColor3f(0.0, 0.0, 0.0);
 glCallList(text_list);
 glPopMatrix();
 glMatrixMode(GL_PROJECTION);
 glPopMatrix();
 glPopAttrib();
 glMatrixMode(GL_MODELVIEW);
 text_list=glGenLists(1);
 glNewList(text_list,GL_COMPILE);

 DRAW_STRING(x, y, string , GLUT_BITMAP_TIMES_ROMAN_24);
 glEndList();

glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
}
void DRAW_STRING(int x, int y, char *string, void *font){
 int len, i;
 glRasterPos2f(x, y);
 len = (int) strlen(string);
 for (i = 0; i < len; i++){
   glutBitmapCharacter(font, string[i]);
 }
}
//--------------------------------------------------------
// マウスドラッグによる回転
//--------------------------------------------------------
//////////////////////////////////////////////
// クォータニオンの積 r <- p x q
static void qmul(double r[], const double p[], const double q[])
{
 r[0] = p[0] * q[0] - p[1] * q[1] - p[2] * q[2] - p[3] * q[3];
 r[1] = p[0] * q[1] + p[1] * q[0] + p[2] * q[3] - p[3] * q[2];
 r[2] = p[0] * q[2] - p[1] * q[3] + p[2] * q[0] + p[3] * q[1];
 r[3] = p[0] * q[3] + p[1] * q[2] - p[2] * q[1] + p[3] * q[0];
}
/////////////////////////////////////////////
// 回転の変換行列 r <- クォータニオン q
static void qrot(double r[], double q[]){
 double x2 = q[1] * q[1] * 2.0;
 double y2 = q[2] * q[2] * 2.0;
 double z2 = q[3] * q[3] * 2.0;
 double xy = q[1] * q[2] * 2.0;
 double yz = q[2] * q[3] * 2.0;
 double zx = q[3] * q[1] * 2.0;
 double xw = q[1] * q[0] * 2.0;
 double yw = q[2] * q[0] * 2.0;
 double zw = q[3] * q[0] * 2.0;

 r[ 0] = 1.0 - y2 - z2;
 r[ 1] = xy + zw;
 r[ 2] = zx - yw;
 r[ 4] = xy - zw;
 r[ 5] = 1.0 - z2 - x2;
 r[ 6] = yz + xw;
 r[ 8] = zx + yw;
 r[ 9] = yz - xw;
 r[10] = 1.0 - x2 - y2;
 r[ 3] = r[ 7] = r[11] = r[12] = r[13] = r[14] = 0.0;
 r[15] = 1.0;
}
/////////////////////////////////////////////
// ばねの描画
void drowSolidSpring(double x0, double y0, double z0, double x1, double y1, double z1){//
  double l = sqrt(pow(x1-x0,2)+pow(y1-y0,2)+pow(z1-z0,2));
  double rx = (x1-x0)/l;
  double ry = (y1-y0)/l;
  double radian_theta, theta;
  radian_theta = acos((z1-z0)/l);
  theta = radian_theta/(2.0*PI) * 360.0;
  glTranslated(x0, y0, z0);
  if(rx==0.0 && ry==0.0)  glRotated( theta , 1.0, 0.0, 0.0);
  else  glRotated( theta , -ry, rx, 0.0);
  skSolidSpring(10, 10, 5, 1.0, 0.8, l);
}
void skSolidSpring(int Nm, int Ns, int Np, double radius, double ratio, double length){
  //
  //Nm:主円周分割点数,Ns:断面円周分割点数,Np:ピッチ数
  //radius:主円周軸半径,ratio:断面半径/主円周軸半径,length:全長
  int i, ii, j, jj, k, k1, k2;
  double phai, ph1, ph2; //主円周分割点のx軸に対する偏角
  double theta, th1, th2; //断面円周分割点のx-z平面に対する偏角
  float p[21][21][101][3];//p[44100][3]; //各点の座標
  float p1[3], p2[3], p3[3], p4[3];
  double rr[21], zz[21];
  double pitch, dp, hh;
  double r1, r2;

  r1 = radius; //主円周軸半径
  r2 = ratio * r1;//断面半径
  if( Nm > 20 ) Nm = 20;
  if( Ns > 20 ) Ns = 20;
  if( Np > 100 ) Np = 100;
  pitch = length / (double)Np;
    //縮んだときの制限
  if(pitch < 2 * r2) pitch = 2.0f * r2 ;

  dp = pitch / (double)Nm;

  //基本断面(x-z)の座標
  for(j = 0; j < Ns; j++)
  { theta = PI - 2.0 * PI*(double)j/(double)Ns;
    rr[j] = r1 + r2 * cos(theta); //原点からの距離
    zz[j] = r2 * sin(theta);//ｚ
  }

  //他の断面の座標
  hh = 0;
  for(k = 0; k < Np; k++)
    for(i = 0; i < Nm; i++)
    { phai = 2.0 * PI * (double)i/(double)Nm;
      for(j = 0; j < Ns; j++)
      {
        p[i][j][k][0] = (float)(rr[j] * cos(phai)); //x座標
        p[i][j][k][1] = (float)(rr[j] * sin(phai)); //y
        p[i][j][k][2] = (float)(zz[j] + hh) ;              //z
      }
      hh += dp;//中心軸の高さがdpずつ増加
    }

  //最終端(k=Np-1,i=Nm)
  k = Np - 1; i = Nm;
  for(j = 0; j < Ns; j++){
    phai = 0.0;
    p[i][j][k][0] = (float)(rr[j] * cos(phai)); //x座標
    p[i][j][k][1] = (float)(rr[j] * sin(phai)); //y
    p[i][j][k][2] = (float)(zz[j] + hh) ;            //z
  }

  //頂点列を定義し描画
  for(k = 0; k < Np; k++)
  for(i = 0; i < Nm; i++){
    ii = i+1;
    k1 = k; k2 = k;
    if(ii == Nm) {
      if(k < Np-1) { ii = 0; k2 = k + 1; }
    }
    ph1 = 2.0*PI*(double)i / (double)Nm;
    ph2 = 2.0*PI*(double)ii / (double)Nm;
    for(j = 0;j < Ns; j++)
    {
      jj = j+1;
      if(jj == Ns) jj = 0;
      th1 = PI - 2.0 * PI * (double)j / (double)Ns;
      th2 = PI - 2.0 * PI * (double)jj / (double)Ns;
      //面の頂点   x座標                       y座標                       z座標
      p1[0] = p[i][j][k1][0]   ; p1[1] = p[i][j][k1][1]   ; p1[2] = p[i][j][k1][2];
      p2[0] = p[i][jj][k1][0]  ; p2[1] = p[i][jj][k1][1]  ; p2[2] = p[i][jj][k1][2];
      p3[0] = p[ii][jj][k2][0] ; p3[1] = p[ii][jj][k2][1] ; p3[2] = p[ii][jj][k2][2];
      p4[0] = p[ii][j][k2][0]  ; p4[1] = p[ii][j][k2][1]  ; p4[2] = p[ii][j][k2][2];

      glBegin(GL_QUADS);
        glNormal3d(cos(th1)*cos(ph1),cos(th1)*sin(ph1),sin(th1));glVertex3fv(p1);
        glNormal3d(cos(th2)*cos(ph1),cos(th2)*sin(ph1),sin(th2));glVertex3fv(p2);
        glNormal3d(cos(th2)*cos(ph2),cos(th2)*sin(ph2),sin(th2));glVertex3fv(p3);
        glNormal3d(cos(th1)*cos(ph2),cos(th1)*sin(ph2),sin(th1));glVertex3fv(p4);
      glEnd();
    }
  }
  //始端（k=0)
  glBegin(GL_POLYGON);
    glNormal3d(0.0,-1.0,0.0);
    for(j = Ns-1; j >= 0; j--) glVertex3fv(p[0][j][0]);
  glEnd();
  //終端（k=Np)
  glBegin(GL_POLYGON);
    glNormal3d(0.0,1.0,0.0);
    for(j = 0; j < Ns; j++) glVertex3fv(p[Nm][j][Np-1]);
  glEnd();

}
/////////////////////////////////////////////
// 直方体
void Cuboid(double a , double b, double c){
  GLdouble vertex[][3] = {
      { -a/2.0, -b/2.0, -c/2.0 },
      {  a/2.0, -b/2.0, -c/2.0 },
      {  a/2.0,  b/2.0, -c/2.0 },
      { -a/2.0,  b/2.0, -c/2.0 },
      { -a/2.0, -b/2.0,  c/2.0 },
      {  a/2.0, -b/2.0,  c/2.0 },
      {  a/2.0,  b/2.0,  c/2.0 },
      { -a/2.0,  b/2.0,  c/2.0 }
    };
  int face[][4] = {//面の定義
      { 0, 1, 2, 3 },
      { 1, 2, 6, 5 },
      { 4, 5, 6, 7 },
      { 0, 4, 7, 3 },
      { 0, 1, 5, 4 },
      { 2, 3, 7, 6 }
    };
  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 }
  };
  glBegin(GL_QUADS);
  for (int j = 0; j < 6; ++j) {
    glNormal3dv(normal[j]); //法線ベクトルの指定
    for (int i = 0; i < 4; ++i) {
      glVertex3dv(vertex[face[j][i]]);
    }
  }
  glEnd();
}

VisualC++ と OpenGL を利用した仮想物理実験室

第０章　仮想物理実験室の構築

第１章　様々な運動

第２章　ニュートンの運動方程式

第３章　剛体の運動（エネルギー保存則と運動量保存則）

付録

  • 【Ａ－１】参考文献
    ・（A-1-1）OpenGL について
    ・（A-1-2）VisualC++ について
    ・（A-1-3）物理シミュレーション
    ・（A-1-4）数値計算

未分類

力学

量子力学

波動論



▲このページのトップ ｜ NPO法人 natural science トップ

関連記事

仮想物理実験室

【第８話】有限の高さの障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.29 【遠藤 理平｜仮想物理実験室】

【第７話】無限に高い障壁に向けた電子パルスの照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.24 【遠藤 理平｜仮想物理実験室】


【第６話】無限に高い障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.23 【遠藤 理平｜仮想物理実験室】

【第５話】電子パルスの運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.21 【遠藤 理平｜仮想物理実験室】


【第４話】電子パルスの作り方【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.20 【遠藤 理平｜仮想物理実験室】

【第３話】自由粒子の運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.19 【遠藤 理平｜仮想物理実験室】


【第２話】自由粒子のスナップショット【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.18 【遠藤 理平｜仮想物理実験室】

【第１話】プログラムの動作確認【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.17 【遠藤 理平｜仮想物理実験室】


１０日で作る量子ドットコンピュータ～動作原理を数値シミュレーション！～【５日目】１量子ビット万能量子ゲートを設計しよう！ 2019.09.06 【遠藤 理平｜仮想物理実験室】

１０日で作る量子ドットコンピュータ～動作原理を数値シミュレーション！～【４日目】量子井戸の電子状態の基礎実験３：２重量子井戸で量子ビットを改良！ 2019.08.27 【遠藤 理平｜仮想物理実験室】



▲このページのトップ ｜ NPO法人 natural science トップ