シュレディンガー方程式に従う粒子の時間発展
無限に深い２次元井戸型ポテンシャル３

シュレディンガー方程式に従う粒子の時間発展：無限に深い２次元井戸型ポテンシャルの時間発展において、 (k_x,k_y)=(0,0)、 (k_x,k_y)=(10,0)の場合のシミュレーションを行いました。 今回は(k_x,k_y)=(10,10)の場合のシミュレーションを行います。

(k_0x, k_0y) = (10,10) の場合の波数空間におけるガウス分布

(k_0x, k_0y) = (10,10) の場合の時間発展

縦軸：実空間 x 0と10[nm]がポテンシャル障壁
横軸: 実空間 y 0と10[nm]がポテンシャル障壁
色：黒色→存在確率0
パルスの中心座標：r_0 = (L/2,L/2)
時間間隔： 10^{-15}[s]

VisualC++ ＋ OpenGL

あらかじめ計算しておいたデータを読み込んで、OpenGL を用いて描画します。

//////////////////////////////////////////////////////////////////////////
// 仮想物理実験室(ver 1.2)
// データを読み込んで波動関数の時間発展を描画
//////////////////////////////////////////////////////////////////////////
#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:する
//ofstream file1( "ball1.data" );

//--------------------------------------------------------
// 仮想物理実験室変数の定義
//--------------------------------------------------------
double t = 0.0;    //時刻
double dt= 1.0;  //時間刻み
int tn = 0;        //ステップ数
int skipBMP = 1; //BMP書出し回数の間引き数
int skipCal = 0;  //描画回数の間引き
int skipDat = 1000;//データ書出し回数の間引き回数

const int N=101;           //系の長さ
const double c=1.0;           //波の速度
const double dd=2.0 ;         //系の刻み幅
double phi[N+1][N+1];

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, 300, 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[20];
char t_char2[20];
//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;    //ラインの表示

// 表面の描写用 --------------------------------------------------------
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  | 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(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);
     //////////////////////////////////////////

 //////////////////////////////////////////
 //透視変換行列の設定
 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, -140.0, 150.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();

  t = double(tn) * dt;
 //////////////////////////////////////////
 //文字の描画
 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(){
  
  string InputFolder = "C:/gnuplot/bin/3_10/", ff="";
  ostringstream fname;
  ifstream fin_t;

	fname  << InputFolder << tn << ".data" ;
  string fname_c = fname.str();
	fin_t.open(fname_c.c_str() , ios::in);	fname.str("");

  int ii=0, jj=0;
	while(!fin_t.eof()){

		double x = 0.0, y=0.0;
    double z = 0.0;
		fin_t >> x; //x座標
		fin_t >> y; //y座標
		fin_t >> z;  //phiの値
    if(x!=0.0 || y!=0.0 || z!=0.0 ) phi[ii][jj] = z * 10000.0;//
    jj++;
    if(jj==N){
      jj=0;
      ii++;
    }
  }
  fin_t.close();
  
  //cin.get();//入力待ち
  
	//// 描画 --------------------------------------------------------
	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] = phi[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,phi[i][j]);
				glVertex3d(i,j+1,phi[i][j+1]);
				glEnd();
				glBegin(GL_LINES);
				glVertex3d(i,j,phi[i][j]);
				glVertex3d(i+1,j,phi[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++;
}


//--------------------------------------------------------
// アイドル時に呼び出される関数
//--------------------------------------------------------
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;
}

波数空間のガウス分布　gnuplotテンプレート

set terminal jpeg  enhanced font "Times" 20 size 600, 480
set tics font 'Times,18'
set nokey
set size square
set rmargin 0
set lmargin 0

set pm3d
set pm3d map
set ticslevel 0
set cbrange[0:1]
set palette defined (0 "black",  1 "white")
set xrange[-30:30]
set yrange[-30:30]
#set xtics -5, 1, 5
#set ytics -5, 1, 5

set output "K.jpg"
splot "K.data" with pm3d