%   EPS 109 Movie Project  Function file for wdot ode.
%   Nathan Haese  December 2, 2009
%   

function wdot = wdotHaese2(t,w)
global e1 qd q1 q2 md m1 m2
%   interparticle separations cubed: r3d1 reads as "r cubed of dust 1"
    r3d1    = (((w(1,1)-w(5,1))^2 + (w(2,1)-w(6,1))^2))^1.5;
    r3d2    = (((w(1,1)-w(9,1))^2 + (w(2,1)-w(10,1))^2))^1.5;
    r312    = (((w(5,1)-w(9,1))^2 + (w(6,1)-w(10,1))^2))^1.5;
%   force on dust by ion-1: Fd1x reads as x component of force on dust by 1
    Fd1x    = e1*qd*q1*(w(1,1)-w(5,1))/r3d1;
    Fd1y    = e1*qd*q1*(w(2,1)-w(6,1))/r3d1;
%   force on dust by ion-2
    Fd2x    = e1*qd*q2*(w(1,1)-w(9,1))/r3d2;
    Fd2y    = e1*qd*q2*(w(2,1)-w(10,1))/r3d2;
%   force on ion-1 by dust
    F1dx    = -Fd1x;
    F1dy    = -Fd1y;
%   force on ion-1 by ion-2
    F12x    = e1*q1*q2*(w(5,1)-w(9,1))/r312;
    F12y    = e1*q1*q2*(w(6,1)-w(10,1))/r312;
%   force on ion-2 by dust
    F2dx    = - Fd2x;
    F2dy    = - Fd2y;
%   force on ion-2 by ion-1
    F21x    = - F12x;
    F21y    = - F12y;
% wdot for the dust
    wdot(1,1)   =   w(3,1);                 % the old vx
    wdot(2,1)   =   w(4,1);                 % the old vy
    wdot(3,1)   =   (Fd1x + Fd2x)/md;       % x acceleration on dust
    wdot(4,1)   =   (Fd1y + Fd2y)/md;       % y acceleration on dust
% wdot for ion-1
    wdot(5,1)   =   w(7,1);
    wdot(6,1)   =   w(8,1);
    wdot(7,1)   =   (F1dx + F12x)/m1;
    wdot(8,1)   =   (F1dy + F12y)/m1;
% wdot for ion-2
    wdot(9,1)   =   w(11,1);
    wdot(10,1)  =   w(12,1);
    wdot(11,1)  =   (F2dx + F21x)/m2;
    wdot(12,1)  =   (F2dy + F21y)/m2;
