program survey
       implicit none
       real*8 h 
       integer j

       write(*,*) "      for r=1,        "
       write(*,*) " nsteps      dt=1/nsteps   error r-r0  rel error%"
       h=0.001
       do while (h .le. 0.1)
       
       call rk(h)
       h=h+0.001
       enddo
       end
      

        
      subroutine rk(dt) 
      implicit none
      integer i,j,k,nstep
      real*8 time,px,py,pxx,pyy,vx,vy,pi,dt,vxx,vyy
      real*8 vxg,vyg,e,r1,r2,pnx,pny,ex,t
c       Define pi
      pi = 4.0d0*atan(1.0d0)
      

c===========================================
c  three ways to calculate the velocity 
c  selected by vtype
c============================================
c     angular vel=2*pi
c        dt=0.04
        nstep=1/dt+0.5d0
        px=0.0d0
        r1=1.0d0
        py=r1
        i=0
        t=0.0d0
c        write(*,*)"i,px,py,dt,r1,r2,e,ex,dt^3"
c        write(*,*)i,px,py,dt,r1,r2,e,ex,dt**3.0d0
c        write(*,*)
100    format(I4,(F25.16),E25.16,F15.8)
c        write(*,*)"i,r1,r2,e,rel"
c        write(*,100)i,r1,r2,e,(e/r1)*100.0d0
      do i=1,nstep
      t=t+dt
c     first RK step
      vx=2.0d0*pi*(-py)
      vy=2.0d0*pi*(px)

      pxx=px+vx*dt
      pyy=py+vy*dt

c     Second RK step
      vxx=2.0d0*pi*(-pyy)
      vyy=2.0d0*pi*(pxx)
      pnx=px
      pny=py
      px=0.5d0*(px+pxx+vxx*dt)
      py=0.5d0*(py+pyy+vyy*dt)
      r2=sqrt(px*px+py*py)
      e=abs(r1-r2)
      pnx=-sin(2.0d0*pi*t)
      ex=abs(px-pnx)
c      write(*,*)i,px,py,dt,r1,r2,e,ex,dt**3.0d0
c      write(*,*),px,pnx,r1,r2,ex,e
c       write(*,100)i,r1,r2,e,(e/r1)*100
        enddo
       write(*,100)nstep,dt,e,(e/r1)*100
       return
        end