function acc = orbit( t,mu1,mu2,d1,d2,r,v )
%
% computes an acceleration vector acc for a local orbit geometry that
% follows from the restricted three body problem
%
omega = sqrt((mu1+mu2)./((d1+d2).^3));
theta = omega * t;
p = [-d1*cos(theta)  d1*sin(theta) 0];
q = [+d2*cos(theta) -d2*sin(theta) 0];
term1 = sqrt(sum((r - p).^2))^3;
acc = [ 0 0 0 ];
if (term1 ~= 0),
   acc = acc - mu1/term1*(r-p);
end
term2 = sqrt(sum((r - q).^2))^3;
if (term2 ~= 0),
   acc = acc - mu2/term2*(r-q);
end
acc = acc + 2.0*omega*[ v(2) -v(1) 0 ] + omega*omega*r;