Edited at

More than 3 years have passed since last update.

[参考]

# 一様磁場中の荷電粒子の運動

m\frac{d \mathbf{v} }{dt} = q \mathbf{v} × \mathbf{B}


$\mathbf{B}$は一様で$z$方向、q>0とする。

$x$成分

\frac{dv_x}{dt} = \frac{qB}{m}v_y


$y成分$

\frac{dv_y}{dt} = -\frac{qB}{m}v_x


u_x = \frac{v_x}{v_0}
\\
u_y = \frac{v_y}{v_0}
\\
τ = \frac{qBt}{m}
\\
X = \frac{qBx}{mv_0}
\\
Y = \frac{qBy}{mv_0}


B = B_0(1 + α(\frac{y}{R_L})^2) = B_0(1 + αY^2)
\\
( Y = \frac{qB}{mv_0}y )


$R_L$ はラーモア半径（Larmor radius）

# ドリフトを考慮する

yの正方向に行くほど、$B$が強くなる場合を考える。

\frac{dv_x}{dt} = \frac{qB}{m}v_y
\\
\hspace{125px} = \frac{q}{m}B_0(1 + αY^2)v_y
\\
⇔ \frac{d(v_0u_x)}{dt}= \frac{q}{m}B_0(1 + αY^2)v_y
\\
⇔ \frac{du_x}{dt}= \frac{q}{m}B_0(1 + αY^2)\frac{v_y}{v_0}
\\
⇔ \frac{du_x}{dt}= \frac{q}{m}B_0(1 + αY^2)u_y \ … ①


ここで

\frac{du_x}{dt} = \frac{du_x}{dτ}\frac{dτ}{dt} = \frac{du_x}{dτ}\frac{qB_0}{m} \ …②


①と②より

\frac{du_x}{dτ}\frac{qB_0}{m} = \frac{qB_0}{m}(1 + αY^2)u_y
\\
∴ \frac{du_x}{dτ} = (1 + αY^2)u_y


\frac{du_y}{dτ} = -(1 + αY^2)u_x


# サンプル

sample.c

#include <stdio.h>
#include <math.h>

int main(void) {

FILE *fp;

double x1, x2, x, y1, y2, y, vx1, vx2, vy1, vy2, t, dt, kx1, kx2, kx3, kx4, kx;
double ky1, ky2, ky3, ky4, ky, ux1, ux2, ux3, ux4, ux, uy1, uy2, uy3, uy4, uy, al;

fp = fopen("data-sample.csv", "w");

x1 = 0;
y1 = 0;
vx1 = 1;
vy1 = 0;
dt = 0.01;
al = 0.1;

for (int i=0; i<=5000; i++) {
t = i*dt;

kx1 = dt*vx1;
ky1 = dt*vy1;
ux1 = dt*(1 + al*y1*y1)*vy1;
uy1 = -dt*(1 + al*y1*y1)*vx1;

kx2 = dt*(vx1 + ux1/2.);
ky2 = dt*(vy1 + uy1/2.);
ux2 = dt*( (1 + al*y1*y1)*vy1 + uy1/2. );
uy2 = dt*( -( (1 + al*y1*y1)*vx1 + ux1/2. ) );

kx3 = dt*(vx1 + ux2/2.);
ky3 = dt*(vy1 + uy2/2.);
ux3 = dt*( (1 + al*y1*y1)*vy1 + uy2/2. );
uy3 = dt*( -( (1 + al*y1*y1)*vx1 + ux2/2. ) );

kx4 = dt*(vx1 + ux3);
ky4 = dt*(vy1 + uy3);
ux4 = dt*( (1 + al*y1*y1)*vy1 + uy3 );
uy4 = dt*( -( (1 + al*y1*y1)*vx1 + ux3 ) );

kx = (kx1 + 2*kx2 + 2*kx3 + kx4)/6.;
ky = (ky1 + 2*ky2 + 2*ky3 + ky4)/6.;

x2 = x1 + kx;
y2 = y1 + ky;

ux = (ux1 + 2*ux2 + 2*ux3 + ux4)/6.;
uy = (uy1 + 2*uy2 + 2*uy3 + uy4)/6.;
vx2 = vx1 + ux;
vy2 = vy1 + uy;

if (fmod(i, 10)<0.01) {
printf("t = %5.1f, 計算値 -> x=%f, y=%f\n", t, x1, y1);
fprintf(fp,"%f, %f \n", x1, y1);
}

x1 = x2;
y1 = y2;
vx1 = vx2;
vy1 = vy2;

}
fclose(fp);

return 0;
}


こんな感じにグルグル回る。