      subroutine rkqcs(r,drdt,t,dttry,eps,rscal,dtdid,dtnext,
     &                    n,nStts)
      integer n, nStts
      double precision r(n),drdt(n),t,dttry,eps,rscal(n),dtdid,dtnext
c
c     Integrate ODE through one timestep and adjust timestep 
c     according to required accuracy (quality controlled step).
c     Routine based closely on Numerical Recipes's rkqs.
c
c         
c     INPUT/OUTPUT ARGUMENTS:
c
c        r -- The vector of particle 6-tuples at time t+dt, where
c             dt = dtdid
c
c        t -- On input the time t. On output, the time t+dt (which is
c             t+dtdid)
c
c     INPUT ARGUMENTS:
c
c        drdt -- The vector of particle derivatives at time t
c
c        dttry -- The timestep to try.
c        
c        eps -- A relative error. True error is given by rscal*eps
c
c        rscal -- A vector of values which provide absolute error by
c                 err(i) = rscal(i)*eps
c
c        derivs -- The name of the subroutine (function pointer) to
c                  evaluate the derivatives
c
c        n     -- The dimension of the system
c
c     OUTPUT ARGUMENTS:
c
c
c        dtdid -- The actual timestep achieved within the accuracy
c                 constraints
c
c        dtnext -- An estimated timestep to try for the next step.
c                  The estimate is based on local error behaviour.
c
c      include 'bparm.h'
c
c
      logical debug
      common /debug/ debug

      double precision SAFETY,POWGRW,POWSHR,ERRCON
      parameter( SAFETY = 0.9d0 )
      parameter( POWGRW = -0.2d0 )
      parameter( POWSHR = -0.25d0 )
c      parameter( ERRCON = (5.d0/SAFETY)**(1.D0/POWGRW) )
c      parameter( ERRCON = (5.d0/0.9d0)**(1.d0/(-0.2d0)))
c
      integer NMAX
      parameter( NMAX = 20 )
c
      integer i
      double precision rtmp(NMAX),rerr(NMAX),err,dt,dtnew
c
c      write(*,*) 'Inside rkcs'
      ERRCON = (5.d0/SAFETY)**(1.D0/POWGRW) 

      dt = dttry
c
1     call rkck(t, r, drdt, dt, rtmp, rerr, n, nStts)
      if(nStts .eq. 0)return
c
c     Error processing
c
      err = 0.0d0
      do i=1,n
        err = max(err , dabs( rerr(i) / rscal(i) ) )
      enddo
c
c     Scale error by required error
c
      err = err/eps
c
      if (err .gt. 1.d0) then
        dtnew = SAFETY*dt*(err**POWSHR)

        dt = max(dabs(dtnew),0.1*dabs(dt))*nint(dt/dabs(dt))
c
c       Check for stepsize underflow!
c

        if ( t+dt .eq. t ) then
c           call error_printf('rkqcstep: Stepsize underflow!') 
        endif
c
        goto 1 ! ugggghhh!!!

      else
        if (err .gt. ERRCON) then
          dtnext = SAFETY*dt*(err**POWGRW)
        else
          dtnext = 5.d0*dt
        endif
        dtdid = dt
        t = t+dt
c
c       Copy the solution vector into r for return
c
        do i=1,n
          r(i) = rtmp(i)
        enddo

        return
      endif
      end

c******************************************************************

      subroutine rkck(x, y, dydx, h, yout, yerr, n, nStts)

      integer n, nStts
      double precision x, y(n), dydx(n), h, yout(n), yerr(n)
c
c     Advance the ODE system defined by y'(n)=dydx(n) from x
c     to x+h via a Cash-Karp Runge-Kutta step. Error estimates
c     for the y(n) are returned in yerr(n).
c
c     REFERENCE:
c
c        Numerical Recipes, p. 713.
c
c     INPUT ARGUMENTS:
c
c        x -- Starting value for this step
c
c        y -- The vector of y values at this value of x
c
c        dydx -- The vector of derivatives at this value of x.
c
c        h -- The step size
c
c        n -- The dimension of the system of ODEs.
c
c        derivs -- The subroutine derivs(x, y(1:n), dydx(1:n))
c                  that returns the vector dydx(1:n) for given
c                  x and y(1:n).
c
c     OUTPUT ARGUMENTS:
c
c        yout -- The updated solution vector y(x+h)
c
c        yerr -- A vector of error values representing the error
c                in each of the y(1:n). The error is estimated
c                as the difference in the fourth and fifth order
c                Runge Kutta formulae.
c
c     The maximum dimension of the system is NMAX
c
c
      logical debug
      common /debug/ debug

      integer NMAX
      parameter( NMAX = 50 )
c
c     The Cash-Carp parameters
c
      double precision A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43,B51,
     &       B52,B53,B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,
     &       DC1,DC3,DC4,DC5,DC6
      parameter ( A2=0.2D0, A3=0.3D0, A4=0.6D0, A5=1.D0, 
     &            A6=0.875D0, B21=0.2D0, B31=3.D0/40.D0,
     &            B32=9.D0/40.D0, B41=0.3D0, B42=-0.9D0,
     &            B43=1.2D0, B51=-11.D0/54.D0, B52=2.5D0,
     &            B53=-70.D0/27.D0, B54=35.D0/27.D0,
     &            B61=1631.D0/55296.D0, B62=175.D0/512.D0,
     &            B63=575.D0/13824.D0, B64=44275.D0/110592.D0,
     &            B65=253.D0/4096.D0,
     &            C1=37.D0/378.D0, C3=250.D0/621.D0,
     &            C4=125.D0/594.D0, C6=512.D0/1771.D0,
     &            DC1=C1-2825.D0/27648.D0, 
     &            DC3=C3-18575.D0/48384.D0,
     &            DC4=C4-13525.D0/55296.D0,
     &            DC5=-277.D0/14336.D0, DC6=C6-0.25D0 )
c
      integer i
      double precision ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),
     &       ak6(NMAX),ytemp(NMAX)

c
c     First step
c
      do i=1,n
        ytemp(i)=y(i)+B21*h*dydx(i)
      end do

      call derivs(x+A2*h, ytemp, ak2, nStts)
	if(nStts .eq. 0)return
    
c
c     Second step
c
      do i=1,n
        ytemp(i)=y(i)+h*(B31*dydx(i)+B32*ak2(i))
      end do

      call derivs(x+A3*h, ytemp, ak3, nStts)
	if(nStts .eq. 0)return
c
c     Third step
c
      do i=1,n
        ytemp(i)=y(i)+h*(B41*dydx(i)+B42*ak2(i)+B43*ak3(i))
      end do

      call derivs(x+A4*h,ytemp,ak4,nStts)
	if(nStts .eq. 0)return
c
c     Fourth step
c
      do i=1,n
        ytemp(i)=y(i)+h*(B51*dydx(i)+B52*ak2(i)+B53*ak3(i)+
     &           B54*ak4(i))
      end do

      call derivs(x+A5*h,ytemp,ak5,nStts)
	if(nStts .eq. 0)return
c
c     Fifth step
c
      do i=1,n
        ytemp(i)=y(i)+h*(B61*dydx(i)+B62*ak2(i)+B63*ak3(i)+
     &           B64*ak4(i)+B65*ak5(i))
      end do

      call derivs(x+A6*h,ytemp,ak6,nStts)
	if(nStts .eq. 0)return
c
c     Evaluate yout
c
      do i=1,n
        yout(i)=y(i)+h*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+
     &          C6*ak6(i))
      end do
c
c     The error is estimated as the difference between the 
c     fourth and fith order Runge-Kutta methods
c
      do i=1,n
        yerr(i)=h*(DC1*dydx(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i)
     &         +DC6*ak6(i))
      end do
c
      return
      end






