%             phiKepler.m
% Evolution map for the Kepler problem for general in. conds.
function [q p] = phiKepler(q0,p0,h)
tol=1.e-10;
 r0=norm(q0);     v02=dot(p0,p0);    u=dot(q0,p0);
 a=-1/(v02-2/r0); En=1/sqrt(a*a*a);
 Ec=1-r0/a;       Es=u/(En*a*a);
%Initial guess for the iteration
 x=En*h*a/r0;   x1=x;
 while 1
   c=cos(x);    s=sin(x);
   x=x-(x-Ec*s+Es*(1-c)-En*h)/(1-Ec*c+Es*s);
   if abs(x-x1) < tol, break, end 
   x1=x;
 end
 c=cos(x);     s=sin(x);
 FP=1 - Ec*c+Es*s;
 ff=(c-1)*a/r0+1.d0;    g=h+(s-x)/En;
 ffp=-a*En*s/(FP*r0);   gp=(c-1)/FP+1;

 for i=1:2
   q(i)=ff*q0(i)+g*p0(i);
   p(i)=ffp*q0(i)+gp*p0(i); 
 end             