difsy3.f

Go to the documentation of this file.

      SUBROUTINE difsy3(N,EPS,H,X,Y)
*
*
*       Bulirsch-Stoer integrator.
*       --------------------------
*
*       Works if Gamma = (H - E)/L. For other time transformations 'eps'
*       must be scaled appropriately such that the test is esentially
*       of the form (H - E)/L < EPS. This is achieved by using TFAC.
*       Convergence test: ABS(Q'*DP) < EPS*TFAC & ABS(P'*DQ) < EPS*TFAC.
*       Reference: Mikkola (1987). In 'The Few Body Problem' p. 311.
*       This works only if eqs. are canonical and we have P:s & Q:s.
*       One additional eq. is allowed (e.g. for t'=...??) but not checked.
*
*
      parameter(nmx=17,nmx7=7*nmx)
      IMPLICIT REAL*8 (a-h,o-z)
      common/bssave/  ep(4),dsc,facm,tfac,itfac,jc
      REAL*8  y(n),ya(nmx),yl(nmx),ym(nmx),dy(nmx),dz(nmx),dt(nmx,7),
     &        d(7),x,xn,h,g,b,b1,u,v,c,ta,w,dabs
*
      LOGICAL  konv,bo,kl,gr,fybad
      DATA dt /nmx7*0.0d0/
*     DIMENSION  EP(4)
*     DATA  EP/0.4D-1,0.16D-2,0.64D-4,0.256D-5/
*
*
      SAVE
      nhalf2=(n/2)*2
      jti=0
      fy=1.
      DO i=1,n
      ya(i)=y(i)
      END DO
      CALL derqp3(y(1),y(9),dz(1),dz(9),dz(17))
      IF (jc.GT.0)h=dsc
   10 xn=x+h
      bo=.false.
      m=1
      jr=2
      js=3
      DO j=1,10
      IF(bo)THEN
      d(2)=16d0/9.d0
      d(4)=64.d0/9.d0
      d(6)=256.d0/9.d0
      ELSE
      d(2)=2.25d0
      d(4)=9.d0
      d(6)=3.6d1
      END IF
      IF(j.GT.7)THEN
      l=7
      d(7)=6.4d1
      ELSE
      l=j
      d(l)=m*m
      END IF
      konv=l.GT.3
      m=m+m
      g=h/m
      b=g+g
      m=m-1
      DO i=1,n
      yl(i)=ya(i)
      ym(i)=ya(i)+g*dz(i)
      END DO
      DO k=1,m
      CALL derqp3(ym(1),ym(9),dy(1),dy(9),dy(17))
      DO i=1,n
      u=yl(i)+b*dy(i)
      yl(i)=ym(i)
      ym(i)=u
      END DO
      END DO
*       Switch on tolerance scaling indicator.
      itfac = 1
      CALL derqp3(ym(1),ym(9),dy(1),dy(9),dy(17))
      kl=l.LT.2
      gr=l.GT.5
      fs=0.
      DO i=1,n
      v=dt(i,1)
      c=(ym(i)+yl(i)+g*dy(i))*0.5d0
      dt(i,1)=c
      ta=c
      IF(.NOT.kl)THEN
      DO k=2,l
      b1=d(k)*v
      b=b1-c
      w=c-v
      u=v
      IF(b.NE.0.d0)THEN
      b=w/b
      u=c*b
      c=b1*b
      END IF
      v=dt(i,k)
      dt(i,k)=u
      ta=u+ta
      END DO
      is=i+n/2
      IF(is.GT.nhalf2)is=i-(n/2)
*       Use modified EPS for convergence test.
      dyis=dabs(dy(is))/tfac
      IF(i.GT.nhalf2)dyis=0.0
      IF(konv)THEN
      test=dabs( (y(i)-ta)*dyis )
      IF(test.GT.eps) konv=.false.
      END IF
      IF(.NOT.gr)THEN
      fv=dabs(w)*dyis
      IF(fs.LT.fv) fs=fv
      END IF
      END IF
      y(i)=ta
      END DO
      IF(fs.NE.0.d0)THEN
      fa=fy
      k=l-1
      fy=(ep(k)/fs)**(1./(l+k))
      fa7=0.7*fa
      IF(l.EQ.2)fa7=0.0
      fybad=.NOT.((fa7.GT.fy).OR.(fy.GT.0.7))
      IF(fybad)THEN
      h=h*fy
      jti=jti+1
      IF(jti.GT.5)THEN
      h=0.0
      DO i=1,n
      y(i)=ya(i)
      END DO
      RETURN
      END IF
      goto 10
      END IF
      END IF
      IF(konv)THEN
      x=xn
      h=h*fy
      RETURN
      END IF
      d(3)=4.d0
      d(5)=1.6d1
      bo=.NOT.bo
      m=jr
      jr=js
      js=m+m
      END DO
      h=0.5*h
      goto 10
      END