// ■分子動力学法(Molecular Dynamics Method)
// 右ボタンで停止、左ボタンで再開します。
#include "myWin.h"
#include "stdlib.h"
#include "time.h"
#include "math.h"
#define frand()((double)rand()/(RAND_MAX+1))
#define Xmax 10 // X最大位置
#define Ymax 10 // Y最大位置
#define Nmax 100 // 分子個数の最大数
#define EPSI 1.0 //フィッティングパラメータε(ポテンシャルの深さ)
#define SIGMA 1 //フィッティングパラメータσ(衝突直径)
typedef struct{ double X, Y;} point;
static point F[Nmax]; //分子に働く力
static point R[2][Nmax]; //分子の位置
static point V[2][Nmax]; //分子の速度
double E1; //計算用（E1 = -48ε・σ^13 の値)
double E2; //計算用（E2 = 24ε・σ^7 の値)
double MG; //質量（分子の質量を1とする）
double DT; //時間刻み
double TM; //時間(今はカウントのみ。将来拡張用)
int N; //分子の数
int ID1, ID2; //現時刻と次時刻切替え用添字,
static HDC hBuff; static HBITMAP hBM; //ビットマップ表示用
void drawEllipse(int X1, int Y1, int X2, int Y2, int lColor, int bColor){//■端部の円形
HPEN pen=CreatePen(PS_SOLID, 1, lColor) ; SelectObject(hBuff, pen) ;
HBRUSH brush=CreateSolidBrush(bColor) ; SelectObject(hBuff, brush);
Ellipse(hBuff,X1,Y1, X2,Y2) ; DeleteObject(brush) ;
DeleteObject(pen);
}
void drawBox(int X1, int Y1, int X2, int Y2, int lColor, int bColor){//■端部の円形
HPEN pen=CreatePen(PS_SOLID, 1, lColor) ; SelectObject(hBuff, pen) ;
HBRUSH brush=CreateSolidBrush(bColor) ; SelectObject(hBuff, brush);
Rectangle(hBuff,X1-10,Y1-10, X2+10,Y2+10); DeleteObject(brush) ;
brush=CreateSolidBrush(0xFFFFFF) ; SelectObject(hBuff, brush);
Rectangle(hBuff,X1,Y1, X2,Y2) ; DeleteObject(brush);
DeleteObject(pen);
}
void crearBitmap(){ //■ビットマップクリア
SelectObject(hBuff,GetStockObject(NULL_PEN));
PatBlt(hBuff,0, 0,1000,400,WHITENESS);
}
void compForce(){//■力の計算(Lennard-Jones Potentialによる力)
for(int i=0; i<N;i++){ // Φ(r)=4ε{(σ/r)^12 - (σ/r)^6
F[i].X = F[i].Y = 0; // F(r) =4ε{-12(σ/r)^13 + 6(σ/r)^7}
for(int j=0; j<N;j++){ // E1 = -48ε・σ^13 , E2 = 24 *ε・σ^7
if(i != j){
double dX = R[ID1][j].X - R[ID1][i].X;
double dY = R[ID1][j].Y - R[ID1][i].Y;
double DR = sqrt(dX * dX + dY * dY), FF ;
if(DR < 1) DR =1;
FF = E1 / pow(DR,13) + E2 / pow(DR,7);
F[i].X += FF * dX / DR; F[i].Y += FF * dY / DR;
}
}
}
}
void compVel(){//■速度計算
for(int i=0;i<N;i++){
V[ID2][i].X = V[ID1][i].X + DT * F[i].X / MG;
V[ID2][i].Y = V[ID1][i].Y + DT * F[i].Y / MG;
}
double T = 0; // 温度境界条件(温度一定とするための処理)
for(int i=0;i<N;i++){
double VX = V[ID2][i].X, VY = V[ID2][i].Y;
T += VX * VX + VY * VY;
}
T = sqrt(T) *0.1;
if(T > 0.00001){
for(int i=0;i<N;i++){
V[ID2][i].X /= T; V[ID2][i].Y /= T;
}
}
}
void compPlace(){//■位置計算
for(int i=0;i<N;i++){
R[ID2][i].X = R[ID1][i].X + DT * V[ID2][i].X;
R[ID2][i].Y = R[ID1][i].Y + DT * V[ID2][i].Y;
}
}
void boundWall(){//■壁衝突
double HS = 0.5 * SIGMA, XM=Xmax-HS, YM=Ymax-HS;
for(int i=0;i<N;i++){
double IDX=R[ID2][i].X, IDY=R[ID2][i].Y;
if(IDX < HS){ R[ID2][i].X = HS; V[ID2][i].X = abs(V[ID2][i].X);}
if(IDX > XM){ R[ID2][i].X = XM; V[ID2][i].X = -abs(V[ID2][i].X);}
if(IDY < HS){ R[ID2][i].Y = HS; V[ID2][i].Y = abs(V[ID2][i].Y);}
if(IDY > YM){ R[ID2][i].Y = YM; V[ID2][i].Y = -abs(V[ID2][i].Y);}
}
}
void initData(){//初期設定
N=10; DT=0.05; MG = 1; ID1 = 0; ID2 = 1;TM=0;
E1 = -48 * EPSI * pow((double)SIGMA ,13);
E2 = 24 * EPSI * pow((double)SIGMA , 7);
for(int i=0;i<N;i++) F[i].X = F[i].Y = V[0][i].X =V[0][i].Y = 0;
srand((int)time(NULL));int i=0;
for(int i=0;i<N;i++){
R[0][i].X = frand() * (Xmax - 2)+1 ;
R[0][i].Y = frand() * (Ymax - 2)+1 ;
V[0][i].X=frand(); V[0][i].Y=frand();
int iflag=true, NN=0;
while(iflag && NN<100){//近すぎる分子があったらずらすが、
NN++; //100回を限度とする。
for(int j=0; j<2;j++) for(int k=0; k<i;k++){
iflag=false;
double DX=R[0][i].X-R[0][k].X, DY=R[0][i].Y-R[0][k].Y;
if(DX*DX+DY*DY < 5 ){
iflag=true;
R[0][i].X+=2;if(R[0][i].X> Xmax-1)R[0][i].X=1;
R[0][i].Y+=2;if(R[0][i].Y> Ymax-1)R[0][i].Y=1;
}
}
}
}
}
void display(){//■表示
crearBitmap(); drawBox(30,20, 30+Xmax*20, 20+Ymax*20, 0,0xFF0000);
for (int i = 0; i<N; i++){
drawEllipse((int)(R[ID1][i].X*20+20), (int)(R[ID1][i].Y*20+10),
(int)(R[ID1][i].X*20+40), (int)(R[ID1][i].Y*20+30),0, 0xFF);
}
}
void exec(){
TM=TM+DT;compForce(); compVel(); compPlace(); boundWall();
ID1=ID2; ID2=1-ID2; display();
}
void procTimer(HWND hw, WPARAM wp,LPARAM lp){ //■タイマー処理
exec(); InvalidateRect(hw,NULL,TRUE);
}
void procLButtonDown(HWND hw, WPARAM wp,LPARAM lp){ //■左ボタン処理
SetTimer(hw,1,10,NULL); // 画面キャプチャのときここをコメントにする
//exec(); InvalidateRect(hw,NULL,TRUE); // アニメーションのときここをコメントにする
}
void procRButtonDown(HWND hw, WPARAM wp,LPARAM lp){ //■左ボタン処理
KillTimer(hw,1);
}
void initBitmap(HWND hw){ //■ビットマップ初期設定
HDC hdc=GetDC(hw); hBM=CreateCompatibleBitmap(hdc,1000,400);
hBuff=CreateCompatibleDC(hdc); SelectObject(hBuff,hBM);
crearBitmap(); ReleaseDC(hw,hdc);
}
void procCreate(HWND hw, WPARAM wp,LPARAM lp){ //■Create時処理
SetWindowText(hw,TEXT("分子運動 左ボタンで開始, 右ボタンで停止"));
MoveWindow(hw,0,0,280,280,TRUE);
initBitmap(hw); initData();display(); InvalidateRect(hw,NULL,TRUE);
procLButtonDown(hw, wp,lp);
}
void procPaint(HWND hw, WPARAM wp,LPARAM lp){ //■Paint時処理
PAINTSTRUCT ps; HDC hdc=BeginPaint(hw, &ps);
BitBlt(hdc, 0,0,1000,400,hBuff,0,0,SRCCOPY);
EndPaint(hw,&ps);
}
LRESULT CALLBACK WndProc(HWND hw, UINT msg, WPARAM wp,LPARAM lp){//■Window Procedure
switch(msg){
case WM_DESTROY : PostQuitMessage(0) ; return 0;
case WM_CREATE : procCreate(hw,wp,lp) ; return 0;
case WM_PAINT : procPaint(hw,wp,lp) ; return 0;
case WM_TIMER : procTimer(hw,wp,lp) ; return 0;
case WM_LBUTTONDOWN : procLButtonDown(hw,wp,lp) ; break;
case WM_RBUTTONDOWN : procRButtonDown(hw,wp,lp) ; break;
}
return DefWindowProc(hw,msg,wp,lp);
}