% Author: John Noonan 
% Class: EPS109 Fall 2016
% Date: 11.21.16
% Final Project
function ydot = John_Noonan_kepler_ode_three_objects(t, y)
    rs1 = [y(1); y(2); y(3)];
    vs1 = [y(4); y(5); y(6)];
    rs2 = [y(7); y(8); y(9)];
    vs2 = [y(10); y(11); y(12)];
    rp = [y(13); y(14); y(15)];
    vp = [y(16); y(17); y(18)];
    
    global G ms1 ms2 mp
    
    norm_ps1 = norm(rs1 - rp);
    norm_ps2 = norm(rs2 - rp);
    norm_s12 = norm(rs1 - rs2);
    
    Fs1px = ((-1 * G * ms1 * mp) / (norm_ps1)^3) * (rs1(1) - rp(1));
    Fs1py = ((-1 * G * ms1 * mp) / (norm_ps1)^3) * (rs1(2) - rp(2));
    Fs1pz = ((-1 * G * ms1 * mp) / (norm_ps1)^3) * (rs1(3) - rp(3));
    
    Fps1x = -Fs1px;
    Fps1y = -Fs1py;
    Fps1z = -Fs1pz;
    
    Fs2px = ((-1 * G * ms2 * mp) / (norm_ps2)^3) * (rs2(1) - rp(1));
    Fs2py = ((-1 * G * ms2 * mp) / (norm_ps2)^3) * (rs2(2) - rp(2));
    Fs2pz = ((-1 * G * ms2 * mp) / (norm_ps2)^3) * (rs2(3) - rp(3));
    
    Fps2x = -Fs2px;
    Fps2y = -Fs2py;
    Fps2z = -Fs2pz;
  
    Fs12x = ((-1 * G * ms1 * ms2) / (norm_s12)^3) * (rs1(1) - rs2(1));
    Fs12y = ((-1 * G * ms1 * ms2) / (norm_s12)^3) * (rs1(2) - rs2(2));
    Fs12z = ((-1 * G * ms1 * ms2) / (norm_s12)^3) * (rs1(3) - rs2(3));

    Fs21x = -Fs12x;
    Fs21y = -Fs12y;
    Fs21z = -Fs12z;
    
    ydot = [vs1(1); vs1(2); vs1(3); ...
           (Fs1px + Fs12x)/ms1; (Fs1py + Fs12y)/ms1; (Fs1pz + Fs12z)/ms1; ...
            vs2(1); vs2(2); vs2(3); ...
           (Fs2px + Fs21x)/ms2; (Fs2py + Fs21y)/ms2; (Fs2pz + Fs21z)/ms2; ...
            vp(1); vp(2); vp(3); ...
           (Fps1x + Fps2x)/mp; (Fps1y + Fps2y)/mp; (Fps1z + Fps2z)/mp];

    
    