Basically the same as the previous post.
Previous post: [C language] Calculate the motion of charged particles in a uniform magnetic field by the Runge-Kutta method
It becomes as follows.
\frac{du_x}{dτ} = \frac{E_x}{Bv_0} + u_y
In the code, the electric field is appropriately set to $ a = 1 $.
The motion of a positively charged charged particle emitted at a velocity v0 in the x direction from the origin in a uniform magnetic field in the presence of an electric field is calculated by the Lungequtta method and illustrated.
sample.c
#include <stdio.h>
#include <math.h>
int main(void) {
FILE *data;
double x1, x2, y1, y2, vx1, vx2, vy1, vy2, t, dt, x, y, kx1, kx2, kx3, kx4, kx;
double ky1, ky2, ky3, ky4, ky, hx1, hx2, hx3, hx4, hx, hy1, hy2, hy3, hy4, hy, x0, y0;
data = fopen("data-sample.csv", "w");
x1 = 0;
y1 = 0;
vx1 = 1;
vy1 = 0;
dt = 0.01;
int a = 1;
for (int i=0; i<=5000; i++) {
t = i*dt;
kx1 = dt*vx1;
ky1 = dt*vy1;
hx1 = dt*(vy1 + a);
hy1 = -dt*vx1;
kx2 = dt*(vx1 + hx1/2.);
ky2 = dt*(vy1 + hy1/2.);
hx2 = dt*(vy1 + hy1/2. + a);
hy2 = dt*( -(vx1 + hx1/2.) );
kx3 = dt*(vx1 + hx2/2.);
ky3 = dt*(vy1 + hy2/2.);
hx3 = dt*(vy1 + hy2/2. + a);
hy3 = dt*( -(vx1 + hx2/2.) );
kx4 = dt*(vx1 + hx3);
ky4 = dt*(vy1 + hy3);
hx4 = dt*(vy1 + hy3 + a);
hy4 = dt*( -(vx1 + hx3) );
kx = (kx1 + 2*kx2 + 2*kx3 + kx4)/6.;
ky = (ky1 + 2*ky2 + 2*ky3 + ky4)/6.;
x2 = x1 + kx;
y2 = y1 + ky;
hx = (hx1 + 2*hx2 + 2*hx3 + hx4)/6.;
hy = (hy1 + 2*hy2 + 2*hy3 + hy4)/6.;
vx2 = vx1 + hx;
vy2 = vy1 + hy;
if (fmod(i, 10)<0.01) {
printf("t = %5.1f,Calculated value-> x=%f, y=%f\n", t, x1, y1);
fprintf(data,"%f, %f\n", x1, y1);
}
x1 = x2;
y1 = y2;
vx1 = vx2;
vy1 = vy2;
}
fclose(data);
return 0;
}
The motion of a negatively charged charged particle emitted at a velocity v0 in the y direction from the origin in a uniform magnetic field when there is an electric field is calculated by the Lungekutter method and illustrated.
problem.c
#include <stdio.h>
#include <math.h>
int main(void) {
FILE *data;
double x1, x2, y1, y2, vx1, vx2, vy1, vy2, t, dt, x, y, kx1, kx2, kx3, kx4, kx;
double ky1, ky2, ky3, ky4, ky, hx1, hx2, hx3, hx4, hx, hy1, hy2, hy3, hy4, hy, x0, y0;
data = fopen("data-problem.csv", "w");
x1 = 0;
y1 = 0;
vx1 = 0;
vy1 = 1;
dt = 0.01;
int a = 1;
for (int i=0; i<=5000; i++) {
t = i*dt;
kx1 = dt*vx1;
ky1 = dt*vy1;
hx1 = -dt*vy1;
hy1 = dt*(vx1 + a);
kx2 = dt*(vx1 + hx1/2.);
ky2 = dt*(vy1 + hy1/2.);
hx2 = dt*( -(vy1 + hy1/2.) );
hy2 = dt*(vx1 + hx1/2. + a);
kx3 = dt*(vx1 + hx2/2.);
ky3 = dt*(vy1 + hy2/2.);
hx3 = dt*(-(vy1 + hy2/2.) );
hy3 = dt*(vx1 + hx2/2. + a);
kx4 = dt*(vx1 + hx3);
ky4 = dt*(vy1 + hy3);
hx4 = dt*( -(vy1 + hy3) );
hy4 = dt*(vx1 + hx3 + a);
kx = (kx1 + 2*kx2 + 2*kx3 + kx4)/6.;
ky = (ky1 + 2*ky2 + 2*ky3 + ky4)/6.;
x2 = x1 + kx;
y2 = y1 + ky;
hx = (hx1 + 2*hx2 + 2*hx3 + hx4)/6.;
hy = (hy1 + 2*hy2 + 2*hy3 + hy4)/6.;
vx2 = vx1 + hx;
vy2 = vy1 + hy;
if (fmod(i, 10)<0.01) {
printf("t = %5.1f,Calculated value-> x=%f, y=%f\n", t, x1, y1);
fprintf(data,"%f, %f\n", x1, y1);
}
x1 = x2;
y1 = y2;
vx1 = vx2;
vy1 = vy2;
}
fclose(data);
return 0;
}
Recommended Posts