   program kepler
    use const
    implicit none
    real(16) r(3),th(3),x(3),y(3)
    real(16) cth(3),sth(3)
    real(16) nr(3),nth(3),nx(3),ny(3)    
    real(16) ncth(3),nsth(3)
   !!!!!!!!!!!!!!!!!!!!!    
    kang=1.0
    e=0.5 
    p=1.0 
    th(1)=.7854Q0
    cth(1)=cos(th(1))
    sth(1)=sin(th(1))
    r(1) =p/( 1.0 + e*cos(th(1)) )
    x(1)=r(1)*cos(th(1))
    y(1)=r(1)*sin(th(1))
    th(2)=kang/r(1)**2 ! theta dot 
    r(2) =sth(1)*e*p*th(2) ! r dot
    r(2) =r(2)/( 1.0 + e*cth(1))**2 
    th(3)=-2.0*kang*r(2)/r(1)**3 ! theta dbl dot    
! 1 st  2nd derivs of sin, cos   
    call sine(th,sth)
    call cosn(th,cth)
! derivs of x and y
    call prod(r,cth,x)
    call prod(r,sth,y)
    call numer(th,nr,nx,ny)
    print *,"e,p=",e,p
    print 109,"r =",r
    print 109,"nr=",nr
109 format(a,3g20.12)
    print 109,"th =",th
    print 109,"nth=",nth
    print 109,"x =",x
    print 109,"nx=",nx
    print 109,"y =",y
    print 109,"ny=",ny
    read(*,*)
  end
  subroutine eval(th,r,x,y)
  use const
  real(16) th,cth,sth,r,x,y
  cth=cos(th)
  sth=sin(th)
  r=p/(1.q0+e*cth)
  x=r*cth
  y=r*sth
  end subroutine
  subroutine numer(th,r,x,y)
  real(16) th(3),r(3),x(3),y(3)
  real(16) xth(3),xr(3),xx(3),xy(3)
  real(16) dt/1q-6/,dt2
  xth(2)=th(1)
  xth(1)=th(1)-th(2)*dt
  xth(3)=th(1)+th(2)*dt
  do i=1,3
    call eval(xth(i),xr(i),xx(i),xy(i))
    print 109,i,xth(i),xr(i),xx(i),xy(i)
109 format(I2,4G20.12)
  enddo
  dt2=dt*dt
  r(1)=xr(2)
  r(2)=(xr(3)-xr(1))/(2.0*dt)
  r(3)=(xr(3)+xr(1)-2.0*xr(2))/dt2
  x(1)=xx(2)
  x(2)=(xx(3)-xx(1))/(2.0*dt)
  x(3)=(xx(3)+xx(1)-2.0*xx(2))/dt2
  y(1)=xy(2)
  y(2)=(xy(3)-xy(1))/(2.0*dt)
  y(3)=(xy(3)+xy(1)-2.0*xy(2))/dt2
  end
  subroutine prod(A,B,C)
  real(16) a(3),b(3),c(3)
!  1 st deriv.
  c(1)=a(1)*b(1)
  c(2)=a(1)*b(2)+a(2)*b(1)
! 2nd deriv
  c(3)=a(1)*b(3)+a(3)*b(1)+2.0*a(2)*b(2)
  end
  subroutine cosn(a,b)
  real(16) a(3),b(3),sina1,cosa1
  cosa1=cos(a(1))
  b(1)=cosa1
  sina1=sin(a(1))
  b(2)=-a(2)*sina1
  b(3)=-sina1*a(3)-cosa1*a(2)**2
  end
  subroutine sine(a,b)
  real(16) a(3),b(3),sina1,cosa1
  cosa1=cos(a(1))
  sina1=sin(a(1))
  b(1)=cosa1  
  b(2)=a(2)*sina1
  b(3)=sina1*a(3)+cosa1*a(2)**2
  end
  end
  
  
  