ドリフト運動（電場×磁場;垂直）（重力あり）

本日のソースコード

#include <>

main(){

double a,B,E,g,m,q,C,x,y,z,Vx,Vy,Vz,Vx1,Vy1,dt;

//C=q*B/m（定数）

int i,n;

a=0.0; //空気抵抗
B=2.0; //磁束密度//z方向
E=1.0; //電場//x方向
g=9.8; //重力加速度
m=1.0; //粒子質量
q=5.0; //粒子電荷
x=0.0; //
y=0.0; //
z=0.0; //
Vx=2.0; //
Vy=0.0; //
Vz=0.0; //


dt=0.0001; //時間（パラメータ）
for(i=0;i<100000;i++)
{
C=q*B/m;
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt+q*E*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
Vz=Vz-g*dt-a*Vz*dt;
printf("%e %e %e \n", x, y, z );
}
}

x-y平面への出力結果
x方向へ少しずつ広がっているのは微小時間tの値が大きいことによる誤差と考えられます。

x軸+方向に電場をかけた結果、ドリフト運動してy軸-方向へドリフト運動します。
さらに、振動の中心がすこしx軸+方向へずれているのが分かります。

z軸-方向へ重力場がはたらいているのでだんだん上昇速度が小さくなっています。

三次元plotその2 空気抵抗あり

右の図は下のプログラムの場合です


#include

main(){

double a,x,y,z,Vx,Vy,Vz,Vx1,Vy1,C,dt;

int i,n;

a=0.1;
x=1.0;
y=0.0;
z=0.0;
Vx=2.0;
Vy=0.0;
Vz=0.1;
C=10;
dt=0.0001;
for(i=0;i<100000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e %e \n", x, y,z );
}
}

三次元のplot

先週のグラフを三次元で表してみました
（z軸方向に等速度運動）

gnuplot > plot "re"
→gnuplot > splot "re"
で三次元のplotが表示される
splotはファイルreに保存されたデータから一行目をx座標、二行目をy座標、三行目をz座標としてplotする。


#include

main(){

double a,x,y,z,Vx,Vy,Vz,Vx1,Vy1,C,dt;

int i,n;

a=0.0;
x=1.0;
y=0.0;
z=0.0;
Vx=2.0;
Vy=0.0;
Vz=0.1;
C=10;
dt=0.0001;
for(i=0;i<100000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e %e \n", x, y,z );
}
}

来週までの宿題

・プロットする点の減らし方
・磁場が変化する場合
・"三"次元プロット

磁場中の荷電粒子の運動（空気抵抗あり）

一様磁場：画面手前方向
電荷：＋
初速度：＋x軸方向

#include

main(){

double a,x,y,Vx,Vy,Vx1,Vy1,C,dt;

int i,n;

a=0.1;
x=1.0;
y=0.0;
Vx=2.0;
Vy=0.0;
C=10;
dt=0.0001;
for(i=0;i<1000000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;

Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e \n", x, y);
}