-
Notifications
You must be signed in to change notification settings - Fork 11
/
boris_c_function.c
109 lines (91 loc) · 2.52 KB
/
boris_c_function.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#include "boris_c_function.h"
void boris_c(int N_sub_steps, double Dtt,
double* B_field, double* B_skew,
double* xn1, double* yn1, double* zn1,
double* vxn1, double* vyn1, double* vzn1,
double* Ex_n, double* Ey_n,
double* Bx_n_custom, double* By_n_custom, double* Bz_n_custom,
int custom_B,
int N_mp, int N_multipoles,
double charge, double mass)
{
int p, isub, order;
double Ex_np, Ey_np;
double Bx_n, By_n, Bz_n;
double rexy, imxy, rexy_0;
double tBx, tBy, tBz, tBsq;
double sBx, sBy, sBz;
double vx_prime, vy_prime, vz_prime;
double vx_min, vy_min, vz_min;
double vx_plus, vy_plus, vz_plus;
double vxn1p, vyn1p, vzn1p;
double xn1p, yn1p, zn1p;
const double qm=charge/mass;
for(p=0; p<N_mp; p++)
{
Ex_np = Ex_n[p];
Ey_np = Ey_n[p];
vxn1p = vxn1[p];
vyn1p = vyn1[p];
vzn1p = vzn1[p];
xn1p = xn1[p];
yn1p = yn1[p];
zn1p = zn1[p];
for (isub=0; isub<N_sub_steps; isub++)
{
rexy = 1.;
imxy = 0.;
By_n = B_field[0];
Bx_n = B_skew[0];
Bz_n = 0.;
for(order = 1; order < N_multipoles; order++)
{
/* rexy, imxy correspond to real, imaginary part of (x+iy)^(n-1) */
rexy_0 = rexy;
rexy = rexy_0*xn1p - imxy*yn1p;
imxy = imxy*xn1p + rexy_0*yn1p;
/*
* Bx +iBy = sum[ (k + ik')(x + iy)^(n-1) ]
* where k, k' are the strengths and skew strengths of the magnet
*/
By_n += (B_field[order]*rexy - B_skew[order]*imxy);
Bx_n += (B_field[order]*imxy + B_skew[order]*rexy);
}
if (custom_B){
Bx_n += Bx_n_custom[p];
By_n += By_n_custom[p];
Bz_n += Bz_n_custom[p];
}
tBx = 0.5*qm*Dtt*Bx_n;
tBy = 0.5*qm*Dtt*By_n;
tBz = 0.5*qm*Dtt*Bz_n;
tBsq = tBx*tBx + tBy*tBy + tBz*tBz;
sBx = 2.*tBx/(1.+tBsq);
sBy = 2.*tBy/(1.+tBsq);
sBz = 2.*tBz/(1.+tBsq);
vx_min = vxn1p + 0.5*qm*Ex_np*Dtt;
vy_min = vyn1p + 0.5*qm*Ey_np*Dtt;
vz_min = vzn1p;
//v_prime = v_min + cross(v_min, tB)
vx_prime = vy_min*tBz - vz_min*tBy + vx_min;
vy_prime = vz_min*tBx - vx_min*tBz + vy_min;
vz_prime = vx_min*tBy - vy_min*tBx + vz_min;
//v_plus = v_min + cross(v_prime, sB)
vx_plus = vy_prime*sBz - vz_prime*sBy + vx_min;
vy_plus = vz_prime*sBx - vx_prime*sBz+ vy_min;
vz_plus = vx_prime*sBy - vy_prime*sBx + vz_min;
vxn1p = vx_plus + 0.5*qm*Ex_np*Dtt;
vyn1p = vy_plus + 0.5*qm*Ey_np*Dtt;
vzn1p = vz_plus;
xn1p = xn1p + vxn1p * Dtt;
yn1p = yn1p + vyn1p * Dtt;
zn1p = zn1p + vzn1p * Dtt;
}
xn1[p] = xn1p;
yn1[p] = yn1p;
zn1[p] = zn1p;
vxn1[p] = vxn1p;
vyn1[p] = vyn1p;
vzn1[p] = vzn1p;
}
}