%% Full Code
[odefunc, IC] = electricSet(10, 1e-12);
tspan = [0 1e-12];
[time, data] = ode45(odefunc,tspan,IC);
Local functions
function [odefunc, ICs] = electricSet(N, spacing)
pc = physConst();
c = pc.c;
q = [pc.q -pc.q];
m = [pc.mp pc.me];
elecprots = round(rand(1,N)+1,0); % Generates N 1's or 2's. 1's are protons and 2's are electrons.
q = q(elecprots);
m = m(elecprots);
ICs = zeros([1,6*N]);
for i = 1:N
n = i-1;
x = rand(1,1)*spacing;
y = rand(1,1)*spacing;
z = rand(1,1)*spacing;
v = (.2*randn([1,1])+.4)*c;
direction = rand(1,3)-.5;
direction = direction / norm(direction);
vx = v*direction(1);
vy = v*direction(2);
vz = v*direction(3);
ICs(6*n+1:6*n+6) = [x y z vx vy vz];
end
[Ex, Ey, Ez] = electricField(N,q);
odefunc = @(t,vals) 0;
for i = 1:N
n = i-1;
odefunc = @(t,vals) [odefunc(t,vals);... % Repeating pattern of [x y z vx vy vz]
vals(6*n+4);...
vals(6*n+5);...
vals(6*n+6);...
Ex(vals(6*n+1),vals(6*n+2),6*n+3,vals)*q(i)/m(i);...
Ey(vals(6*n+1),vals(6*n+2),6*n+3,vals)*q(i)/m(i);...
Ez(vals(6*n+1),vals(6*n+2),6*n+3,vals)*q(i)/m(i)];
end
end
function [Ex, Ey, Ez] = electricField(N,q)
pc = physConst;
k = pc.kC;
Ex = @(x,y,z,vals) 0;
Ey = @(x,y,z,vals) 0;
Ez = @(x,y,z,vals) 0;
for i = 1:N
n = i-1;
r = norm([x-vals(3*n+1), y-vals(3*n+2), z-vals(3*n+3)]);
direction = [x-vals(3*n+1), y-vals(3*n+2), z-vals(3*n+3)] / r;
magnitude = k*q(i)/(r^2);
Ex = @(x,y,z,vals) Ex(x,y,z,vals) + magnitude*direction(1);
Ey = @(x,y,z,vals) Ey(x,y,z,vals) + magnitude*direction(2);
Ez = @(x,y,z,vals) Ez(x,y,z,vals) + magnitude*direction(3);
end
end