      subroutine kepler(xm1,xm2,r,v) ! Bernd Berg Feb 10 2009.
C Exact solution of the gravitational 3D 2-body problem.
C Arguments:
C xm1,xm2    the two masses.
C r          vectors of initial positions; 
C v          vectors of initial velocities; 
C            overwritten with final positions and velocities.
      include 'implicit.sta'
      include 'constants.par'
      include 'kepler.par'
C Parameters of kepler.par, which are used here:
C G          Newton's constant.
      include 'common_kepler.f'
C common_kepler.f is initialized here and not changed by any routine. 
      dimension r(3,2),v(3,2) ! Arguments. 
      ltest=.true.
      ltest=.false.
      if(ltest) stop "subroutine kepler: stop ltest."
      xmsum=xm1+xm2       ! Total Mass.
      xmu=xm1*xm2/xmsum   ! Reduced mass.
      alpha=G*xm1*xm2     ! Gravitational coupling (Landau-Lifschiz).
      do id=1,3	
        r12(id)=r(id,1)-r(id,2)
        v12(id)=v(id,1)-v(id,2)
        pmu(id)=xmu*v12(id)
      enddo
      r12(4)=vec_norm(nd,r12)
      v12(4)=vec_norm(nd,v12)
      call vec3d_cross(r12,pmu,xl) ! Calculate angular momentum xl.
      xl(4)=vec_norm(nd,xl)
      Ekcm=fe_kin(nd,ione,xmu,v12)         ! Kinetic energy.
      U=-alpha/r12(4)                      ! Potential energy.
      Ecm=Ekcm+U                           ! Energy in CM frame.
      Ueffmin=-alpha**2*xmu/(two*xl(4)**2) ! Min. eff. Pot.
C 
      xmual=xmu*alpha
      xl2=xl(4)**2
      if(xl2.lt.em8) stop "kepler: |L|=0 not allowed."
      p=xl2/xmual ! Landau-Lifschiz parameter, 2p latus rectum.
      e=sqrt(one+two*Ecm*xl2/(xmual*alpha)) ! Eccentricity.
      cosphi0=(p/r12(4)-one)/e ! Initial angle:
      phi0=acos(cosphi0)
      r12dv12=vec_dot(nd,r12,v12)
      if(r12dv12.lt.zero) phi0=-phi0
C
      return
      end
