2012後期プロ発立嶋 - P.S.P.T. 計算物理研究会

荷電粒子のドリフト運動

発表者　立嶋 ( 30期 )
日時　2012年度後期プロ発
言語　C
目的　電磁気学の講義内容の再現

---

解説

ラーモア運動

磁場中の荷電粒子の運動方程式は

である。簡単のため磁場がz軸方向を向くようとして、この微分方程式を解くと

となり磁場中の荷電粒子は螺旋運動することがわかる。

ドリフト運動

電磁場中の荷電粒子の運動方程式は

である。簡単のため電場がxy平面に平行であるようとして、この微分方程式を解くと

となり電磁場中の荷電粒子はドリフト運動することがわかる。

注意

• DrawPS2を利用する。
• 「ルンゲ＝クッタ法」を用いる。

ソースコード

//荷電粒子のドリフト運動(3.9)
#include "drawps2.h"
#include <math.h>
 
#define GN 160
 
DRAWPS2 drps2, graph;
 
// 変数
double t;			//時間
double x, v, y, w;      	//荷電粒子の座標と速度
 
//定数
const double dt = 0.01;         //時間変化
const double cy = 320.0;        //基準座標    
const double cx = 200.0;        //基準座標        
const double m  = 1.0;          //荷電粒子の質量
const double q  = 1.0;          //荷電粒子の電荷
const double B  = 1.0;          //磁場
const double E  = -0.05;        //電場
 
 
// プロトタイプ宣言
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 q*w*B+E;}        
double fw(double y){ return -q*v*B;}       
 
 
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.5;                   
}
 
void Calculate()
{
     int i;
 
//ルンゲ＝クッタ法
     double kx[4],kv[4];
     double ky[4],kw[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;
 
     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;
}
 
void Display()           
{
//　描写設定
     //基準線 
     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
     );
 
     //荷電粒子
     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
     );    
}
 
void subDisplay()
{
     int i;
     char str[64];
 
     graph.EraseOrbit();
 
     sprintf(str, "x : %f", x);
     graph.DrawText(2, 2, str, RGB(0, 0, 0));
     sprintf(str, "y : %f", y);
     graph.DrawText(2, 22, str, RGB(0, 0, 0));
}