2012夏合宿立嶋 - P.S.P.T. 計算物理研究会

簡単な二体問題

発表者　立嶋 ( 30期 )
日時　2012年度夏合宿プロ発
言語　C
目的　インターネット百科事典「Wikipedia」の記事「重心」(2013/3/9 アクセス)に添付された動画の再現

参考動画1 参考動画2 参考動画3

---

注意

• DrawPS2を利用する。
• 万有引力を受ける2つの物体をシミュレートする。
• 万有引力の式はである。
• 今回は「オイラー法」でなく「ルンゲ＝クッタ法」を用いる。
• 質量と初期条件を調整して動画の再現を図る。

参考動画1のソースコード

//二体問題(3.9)
#include "drawps2.h"
 
#define GN 160
 
DRAWPS2 drps2, graph;
 
// 変数
double t;                       //時間
double x, v, y, w;		//物体1の座標と速度
double a, b, c, d;              //物体2の座標と速度
double Gx,Gy;			//重心座標
 
// 定数             
const double dt = 0.01;         //時間変化
const double m  = 1.0;          //物体1の質量
const double M  = 1.0;          //物体2の質量
const double cy = 320.0;	//基準座標             
const double cx = 200.0;	//基準座標
const double G  = 1.0;		//万有引力定数G（ここでは簡単の為G＝1とする）       
 
// プロトタイプ宣言
void Init();
void Calculate();
void Display();
void subDisplay();
 
// 連立微分方程式
double fx(double v){ return v;}             
double fy(double w){ return w;}             
double fv(double x){ return -G*M*(x-a)/sqrt(((x-a)*(x-a)+(y-b)*(y-b))*((x-a)*(x-a)+(y-b)*(y-b))*((x-a)*(x-a)+(y-b)*(y-b)));}         
double fw(double y){ return -G*M*(y-b)/sqrt(((x-a)*(x-a)+(y-b)*(y-b))*((x-a)*(x-a)+(y-b)*(y-b))*((x-a)*(x-a)+(y-b)*(y-b)));}        
double fa(double b){ return b;}
double fb(double d){ return d;}
double fc(double a){ return -G*m*(a-x)/sqrt(((a-x)*(a-x)+(b-y)*(b-y))*((a-x)*(a-x)+(b-y)*(b-y))*((a-x)*(a-x)+(b-y)*(b-y)));}
double fd(double c){ return -G*m*(b-y)/sqrt(((a-x)*(a-x)+(b-y)*(b-y))*((a-x)*(a-x)+(b-y)*(b-y))*((a-x)*(a-x)+(b-y)*(b-y)));}
 
 
void psMain()
{
     drps2.window.Size(640, 480);      
     drps2.window.Position(50, 100);
     drps2.event.Init(Init);           
     drps2.event.Calculate(Calculate);
     drps2.event.Display(Display);
     drps2.psCreateWindow();
 
graph.window.Size(320, 240);
     graph.window.Position(680, 100);           
     graph.event.Display(subDisplay);
     graph.psCreateWindow();
}
 
void Init()
{
     int i;
 
//　初期条件
     t = 0.0;
     x = 0.5;                               
     y = 0.0;  
     v = 0.0; 
     w = 0.7;
     a = -0.5;
     b = 0.0;
     c = 0.0;
     d = -0.7;
     Gx = (m*x+M*a)/(m+M);
     Gy = (m*y+M*b)/(m+M);
 
}
 
void Calculate()
{
     int i;
 
//ルンゲ＝クッタ法
     double kx[4],kv[4];
     double ky[4],kw[4];
     double ka[4],kc[4];
     double kb[4],kd[4];            
 
     kx[0] = fx(v)*dt;
     kv[0] = fv(x)*dt;
     kx[1] = fx(v+kv[0]/2.0)*dt;                     
     kv[1] = fv(x+kx[0]/2.0)*dt;                     
     kx[2] = fx(v+kv[1]/2.0)*dt;
     kv[2] = fv(x+kx[1]/2.0)*dt;
     kx[3] = fx(v+kv[2])*dt;
     kv[3] = fv(x+kx[2])*dt;
 
     ky[0] = fy(w)*dt;       
     kw[0] = fw(y)*dt;
     ky[1] = fy(w+kw[0]/2.0)*dt;
     kw[1] = fw(y+ky[0]/2.0)*dt;
     ky[2] = fy(w+kw[1]/2.0)*dt;
     kw[2] = fw(y+ky[1]/2.0)*dt;
     ky[3] = fy(w+kw[2])*dt;
     kw[3] = fw(y+ky[2])*dt;
 
     ka[0] = fa(c)*dt;
     kc[0] = fc(a)*dt;
     ka[1] = fa(c+kc[0]/2.0)*dt;
     kc[1] = fc(a+ka[0]/2.0)*dt;
     ka[2] = fa(c+kc[1]/2.0)*dt;
     kc[2] = fc(a+ka[1]/2.0)*dt;
     ka[3] = fa(c+kc[2])*dt;
     kc[3] = fc(a+ka[2])*dt;
 
     kb[0] = fb(d)*dt;
     kd[0] = fd(d)*dt;
     kb[1] = fb(d+kd[0]/2.0)*dt;
     kd[1] = fd(b+kb[0]/2.0)*dt;
     kb[2] = fb(d+kd[1]/2.0)*dt;
     kd[2] = fd(b+kb[1]/2.0)*dt;
     kb[3] = fb(d+kd[2])*dt;
     kd[3] = fd(b+kb[2])*dt;      
 
     t += dt;
     x += (kx[0]+2.0*kx[1]+2.0*kx[2]+kx[3])/6.0;
     v += (kv[0]+2.0*kv[1]+2.0*kv[2]+kv[3])/6.0;
     y += (ky[0]+2.0*ky[1]+2.0*ky[2]+ky[3])/6.0;  
     w += (kw[0]+2.0*kw[1]+2.0*kw[2]+kw[3])/6.0;
 
     a += (ka[0]+2.0*ka[1]+2.0*ka[2]+ka[3])/6.0;
     c += (kc[0]+2.0*kc[1]+2.0*kc[2]+kc[3])/6.0;
     b += (kb[0]+2.0*kb[1]+2.0*kb[2]+kb[3])/6.0;
     d += (kd[0]+2.0*kd[1]+2.0*kd[2]+kd[3])/6.0;
 
 
     Gx = (m*x+M*a)/(m+M);
     Gy = (m*y+M*b)/(m+M);
}
 
void Display()               
{
//　描写設定
 
     //基準線
     drps2.EraseOrbit();                
          drps2.SetPen(PS_SOLID, 1, RGB(0,250, 0));      
 
          drps2.DrawEllipse(
          cy-5.0,
          cx-5.0,           
          cy+5.0,
          cx+5.0
     );
 
     drps2.DrawLine(
          cy-30.0,
          cx,           
          cy+30.0,
          cx
     );
 
     drps2.DrawLine(
          cy,
          cx-30.0,           
          cy,
          cx+30.0
     );
 
     //物体1と物体2を結ぶ補助線
     drps2.SetPen(PS_SOLID, 1, RGB(150, 150, 150));
 
      drps2.DrawLine(
          cy+y*200,
          cx+x*200,           
          cy+b*200,
          cx+a*200
     );
 
     //1個目
     drps2.SetPen(PS_SOLID, 1, RGB(0, 0,250));
 
          drps2.DrawEllipse(
          cy+y*200.0-10.0,
          cx+x*200.0-10.0,           
          cy+y*200.0+10.0,
          cx+x*200.0+10.0
     );
 
     drps2.SetPen(PS_SOLID, 1, RGB(250, 0, 0));
 
     //2個目
     drps2.DrawEllipse(
          cy+b*200.0-10.0,
          cx+a*200.0-10.0,           
          cy+b*200.0+10.0,
          cx+a*200.0+10.0
     );
 
     //重心
     drps2.SetPen(PS_SOLID, 1, RGB(0, 0, 0));
 
     drps2.DrawEllipse(
          cy+Gy*200.0-5.0,
          cx+Gx*200.0-5.0,           
          cy+Gy*200.0+5.0,
          cx+Gx*200.0+5.0
     );
 
}
 
void subDisplay()
{
     int i;
     char str[64];
 
     graph.EraseOrbit();
 
     sprintf(str, "Gx : %f", Gx);
     graph.DrawText(2, 2, str, RGB(0, 0, 0));
     sprintf(str, "Gy : %f", Gy);
     graph.DrawText(2, 22, str, RGB(0, 0, 0));
}