VisualC++ と OpenGL を利用した仮想物理実験室
波動方程式のテスト
今後、マックスウェル方程式とシュレディンガー方程式のアニメーション化の準備のため、
簡単な波動方程式のテストソースをつくってみました。
参考ページ
■VisualC++ と OpenGL を利用した仮想物理実験室
仮想物理実験室の構築 （ver1.0）
■OpenGL（Akita National College of Technology Yamamoto's Laboratory ）
VisualC++ ＋ OpenGL のソース
//////////////////////////////////////////////////////////////////////////
// 波動方程式のテスト
//////////////////////////////////////////////////////////////////////////
#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 0//アニメーション作成用ビットマップの保存 0:しない 1:する
//ofstream file1( "ball1.data" );
//--------------------------------------------------------
// 仮想物理実験室変数の定義
//--------------------------------------------------------
double t = 0.0; //時刻
double dt= 0.1; //時間刻み
int tn = 0; //ステップ数
int skipBMP = 1; //BMP書出し回数の間引き数
int skipCal = 0; //描画回数の間引き
int skipDat = 1000;//データ書出し回数の間引き回数
const int N=51; //系の長さ
const double c=1.0; //波の速度
const double dd=2.0 ; //系の刻み幅
double z[3][N][N];
void SetUp(void){
}
void calculate();
//--------------------------------------------------------
// OpenGL用変数定義
//--------------------------------------------------------
//////////////////////////////////////////
// ウィンドウ生成用
int WindowPositionX = 200; //生成するウィンドウ位置のX座標
int WindowPositionY = 200; //生成するウィンドウ位置のY座標
int WindowWidth = 512; //生成するウィンドウの幅
int WindowHeight = 512; //生成するウィンドウの高さ
char WindowTitle[] = "仮想物理実験室(ver 1.2)"; //ウィンドウのタイトル
//////////////////////////////////////////
// アニメーション用 ビットマップ保存
gl_screenshot gs; //「gl_screenshot」クラスのインスタンス「gs」を宣言
static GLfloat LightPosition[4] = { 0, 0, 70, 1 }; //光源の位置
//////////////////////////////////////////
// マウスドラッグ用
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[10];
char t_char2[10];
//string t_string;
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);
// モードの選択 --------------------------------------------------------
int _Surface = 1; //表面の表示
int _Line = 0; //ラインの表示
//// 「_Surface」「_Line」はどちらかが「1」にする
int _Bitmap = 0; // ビットマップ保存（gifアニメーション作成用）
// 表面の描写用 --------------------------------------------------------
GLfloat ctlpoints[N][N][3];
GLfloat knots[2*N+8] = {0.0}; //ノット数 = 制御点数 + 次数 + 1(今回は3次) N=51 -> 110
GLUnurbsObj *theNurb;
//--------------------------------------------------------
// メイン関数
//--------------------------------------------------------
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); //ディスプレイモードの指定
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(0.0, 0.0, 0.0, 1.0); //背景色
glEnable( GL_DEPTH_TEST ); //デプスバッファを使用：glutInitDisplayMode() で GLUT_DEPTH を指定する
glDepthFunc( GL_LEQUAL );
glClearDepth( 1.0 );
//// デプスバッファとは：描画対象の物体の座標を考慮して、ピクセルあたりの描画内容を決める
theNurb = gluNewNurbsRenderer();
gluNurbsProperty(theNurb, GLU_SAMPLING_TOLERANCE, 25.0);
gluNurbsProperty(theNurb, GLU_DISPLAY_MODE, GLU_FILL);
for(int i=0; i<2*N+8 ; i++){
knots[i] = float(i); //ノットの定義
}
// マウスポインタ位置のウィンドウ内の相対的位置への換算用
sx = 1.0 / (double)WindowWidth;
sy = 1.0 / (double)WindowHeight;
// 回転行列の初期化
qrot(rt, cq);
//初期値---------------------------------------------
for(int i=0;i<N;i++) {
for(int j=0;j<N;j++) {
if(i>20 && i<30 && j>20 && j<30) {
z[0][i][j]=20.0;
}
else {
z[0][i][j]=0.0;
}
}
}
//t=dt*1---------------------------------------------
for(int i=1;i<N-1;i++) {
for(int j=1;j<N-1;j++) {
z[1][i][j]=z[0][i][j]+c*c/2.0*dt*dt/(dd*dd)*(z[0][i+1][j]+z[0][i-1][j]+z[0][i][j+1]+z[0][i][j-1]-4.0*z[0][i][j]);
}
}
for(int i=0;i<N;i++) {
z[1][i][0]=z[1][i][N-1]=z[1][0][i]=z[1][N-1][i]=0.0;//境界条件
}
//////////////////////////////////////////
//透視変換行列の設定
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, -90.0, 60.0, // 視点の位置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 );
//////////////////////////////////////////
// 回転
glMultMatrixd(rt);
//////////////////////////////////////////
// 光源ON
glEnable( GL_LIGHTING );
glEnable( GL_LIGHT0 );
glLightfv( GL_LIGHT0,GL_POSITION,LightPosition );
//////////////////////////////////////////
// 描画
glPushMatrix();//階層化
calculate();
//zahyou();
glPopMatrix();
//////////////////////////////////////////
//文字の描画
glColor3d(0.0, 0.0, 0.0);
strcpy_s(t_char2, "t = ");
sprintf_s(t_char, "%5.2f", t);
strcat_s(t_char2, t_char);
DISPLAY_TEXT(5, 95, t_char2 );
//////////////////////////////////////////
//陰影OFF
glDisable(GL_AUTO_NORMAL);
glDisable(GL_NORMALIZE);
glDisable(GL_LIGHTING);
//////////////////////////////////////////
//ビットマップの保存
#if _BITMAP
if(tn%skipBMP == skipCal){
ostringstream fname;
int tt = tn/skipBMP +10000;
fname << "bitmap/" << tt << ".bmp" ;//出力ファイル名
string name = fname.str();
gs.screenshot(name.c_str(), 24);
cout << tn << " " << t << endl;
}
#endif
glutSwapBuffers(); //glutInitDisplayMode(GLUT_DOUBLE)でダブルバッファリングを利用可
}
//////////////////////////////////////////////////////////////////////////
// 計算＋描画
//////////////////////////////////////////////////////////////////////////
void calculate(){
for(int i=1;i<N-1;i++) {
for(int j=1;j<N-1;j++) {
z[2][i][j]=2.0*z[1][i][j]-z[0][i][j]+c*c*dt*dt/(dd*dd)*(z[1][i+1][j]+z[1][i-1][j]+z[1][i][j+1]+z[1][i][j-1]-4.0*z[1][i][j]);
}
}
for(int i=0;i<N;i++) {
z[2][i][0]=z[2][i][N-1]=z[2][0][i]=z[2][N-1][i]=0.0; //境界条件
}
/*次の計算のために配列の数値を入れかえる。ここで過去の情報は失われる。*/
for(int i=0;i<N;i++) {
for(int j=0;j<N;j++) {
z[0][i][j]=z[1][i][j];
z[1][i][j]=z[2][i][j];
}
}
//// 描画 --------------------------------------------------------
for(int i=0;i<N-1;i++){
for(int j=0;j<N-1;j++){
if(_Surface){
ctlpoints[i][j][0] = float(i-(N-1)/2);
ctlpoints[i][j][1] = float(j-(N-1)/2);
ctlpoints[i][j][2] = z[2][i][j];
}
if(_Line){
glDisable(GL_LIGHTING);
glColor3d(1.0,1.0,1.0);
//glColor3d(abs(ez[i][j]), abs(ez[i][j]), abs(ez[i][j]));
glBegin(GL_LINES);
glVertex3d(i,j,z[2][i][j]);
glVertex3d(i,j+1,z[2][i][j+1]);
glEnd();
glBegin(GL_LINES);
glVertex3d(i,j,z[2][i][j]);
glVertex3d(i+1,j,z[2][i+1][j]);
glEnd();
glEnable(GL_LIGHTING);
}
}
}
if(_Surface){
// 表面表示時の効果 ----------------------------------------------------------
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);
gluBeginSurface(theNurb);
gluNurbsSurface(theNurb,
N+3, knots,
N+3, knots,
N * 3,
3,
&ctlpoints[0][0][0],
4, 4,
GL_MAP2_VERTEX_3);
gluEndSurface(theNurb);
}
tn++;
t = double(tn) * dt;
}
//--------------------------------------------------------
// アイドル時に呼び出される関数
//--------------------------------------------------------
void Idle(){
glutPostRedisplay(); //glutDisplayFunc()を１回実行する
}
//--------------------------------------------------------
// サイズ変更
//--------------------------------------------------------
void Resize(int x, int y){
WindowWidth = x;
WindowHeight = y;
glMatrixMode(GL_MODELVIEW);//行列モードの設定（GL_PROJECTION : 透視変換行列の設定、GL_MODELVIEW：モデルビュー変換行列）
glViewport(0, 0, WindowWidth, WindowHeight);
glMatrixMode(GL_PROJECTION); //行列モードの設定（GL_PROJECTION : 射影行列）
glLoadIdentity();//行列の初期化
gluPerspective(30.0, (double)WindowWidth/(double)WindowHeight, 0.1, 1000.0); //透視投影法の視体積gluPerspactive(th, w/h, near, far);
}
//--------------------------------------------------------
// 文字描画
//--------------------------------------------------------
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(1.0, 1.0, 1.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 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;
}
