// ■分子動力学法(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);
}