%This function solves an ODE for particles with mass and charge.

%Henry Peterson 2016

function dAdt = ode2(A,MS,QS)

%split up A
XX = A(1,:);
YY = A(2,:);
VX = A(3,:);
VY = A(4,:);

%number of particles
nParts = length(XX);

%globals for E field, Length of side
global E L

%array for accelerations
AX = zeros(1,nParts);
AY = zeros(1,nParts);

%small step (for derivative)
dx = L * 1e-6;
dy = L * 1e-6;
for j=1:nParts
    if QS(j) > 0
        %ignore protons
        AX(j) = 0;
        AY(j) = 0;
        continue;
    end
    %start a force variable (x and y comps)
    F = [0, 0];
    for k=1:nParts
        if j ~= k
            %pbc distance in the Y dir
            r = [XX(j)-XX(k), YY(j)-YY(k)];
            if r(2) > L/2
                r(2) = r(2) - L;
            elseif r(2) < -L/2
                r(2) = r(2) + L;
            end
            %get potential (see potential.m)
            V = potential(QS(j),QS(k),r);
            %now calculate the partial derivatives in each direction
            %to get the gradient 
            rdx = [r(1)+dx, r(2)];
            rdy = [r(1), r(2)+dy];
            dVdx = potential(QS(j),QS(k),rdx) - V;
            dVdx = dVdx / dx;
            dVdy = potential(QS(j),QS(k),rdy) - V;
            dVdy = dVdy / dy;
            %F = -dV/dr
            F = F + [-dVdx, -dVdy];
        end
    end
    %add in the E field and resistance
    F = F + E*QS(j);
    %F = ma
    AX(j) = F(1)/MS(j);
    AY(j) = F(2)/MS(j);
end

%return derivative of A w.r.t. time
dAdt = [VX; VY; AX; AY];
