% gravitational forces 

echo on;

% masses and positions, cut-off for singularities 
m = [1, 2, 3]; 
p = [0, 0; 2, 1; -1 3];
rmin = 1/5;
pause;

% define grid
box = [min(p); max(p)] + [-1 -1; 1 1];
gx = linspace(box(1,1),box(2,1),40);
gy = linspace(box(1,2),box(2,2),40);
[x,y] = meshgrid(gx,gy);
pause;

% plot potential
f = 0;
for k=1:length(m); 
   r{k} = max(rmin,sqrt((x-p(k,1)).^2+(y-p(k,2)).^2));
   f = f + m(k) ./ r{k};
end;
surf(x,y,f); shading interp; 
pause;

contourf(x,y,f);
pause;

% plot forcefield
[fx, fy] = gradient(f,gx,gy);
quiver(x,y,fx,fy); 
pause;

hold on

% show particle path for 3 start positions
h=streamline(x,y,fx,fy, [1.5,-1.5,1.5],[3.25,1.4,-.5]);
set(h,'color',[0,1,0],'linewidth',2);
pause;

hold off; echo off;

