Nbody6
difsy4.f
Go to the documentation of this file.
1  SUBROUTINE difsy4(N,EPS,H,X,Y)
2 *
3 *
4 * Bulirsch-Stoer integrator.
5 * --------------------------
6 *
7 * Works if Gamma = (H - E)/L. For other time transformations 'eps'
8 * must be scaled appropriately such that the test is esentially
9 * of the form (H - E)/L < EPS. This is achieved by using TFAC.
10 * Convergence test: ABS(Q'*DP) < EPS*TFAC & ABS(P'*DQ) < EPS*TFAC.
11 * Reference: Mikkola (1987). In 'The Few Body Problem' p. 311.
12 * This works only if eqs. are canonical and we have P:s & Q:s.
13 * One additional eq. is allowed (e.g. for t'=...??) but not checked.
14 *
15 *
16  parameter(nmx=25,nmx7=7*nmx)
17  IMPLICIT REAL*8 (a-h,o-z)
18  common/bssave/ ep(4),dsc,facm,tfac,itfac,jc
19  REAL*8 y(n),ya(nmx),yl(nmx),ym(nmx),dy(nmx),dz(nmx),dt(nmx,7),
20  & d(7),x,xn,h,g,b,b1,u,v,c,ta,w,dabs
21 *
23  DATA dt /nmx7*0.0d0/
24 * DIMENSION EP(4)
25 * DATA EP/0.4D-1,0.16D-2,0.64D-4,0.256D-5/
26 *
27 *
28  SAVE
29  nhalf2=(n/2)*2
30  jti=0
31  fy=1.
32  DO i=1,n
33  ya(i)=y(i)
34  END DO
35  CALL derqp4(y(1),y(13),dz(1),dz(13))
36  IF (jc.GT.0)h=dsc
37  10 xn=x+h
38  bo=.false.
39  m=1
40  jr=2
41  js=3
42  DO j=1,10
43  IF(bo)THEN
44  d(2)=16d0/9.d0
45  d(4)=64.d0/9.d0
46  d(6)=256.d0/9.d0
47  ELSE
48  d(2)=2.25d0
49  d(4)=9.d0
50  d(6)=3.6d1
51  END IF
52  IF(j.GT.7)THEN
53  l=7
54  d(7)=6.4d1
55  ELSE
56  l=j
57  d(l)=m*m
58  END IF
59  konv=l.GT.3
60  m=m+m
61  g=h/m
62  b=g+g
63  m=m-1
64  DO i=1,n
65  yl(i)=ya(i)
66  ym(i)=ya(i)+g*dz(i)
67  END DO
68  DO k=1,m
69  CALL derqp4(ym(1),ym(13),dy(1),dy(13))
70  DO i=1,n
71  u=yl(i)+b*dy(i)
72  yl(i)=ym(i)
73  ym(i)=u
74  END DO
75  END DO
76 * Switch on tolerance scaling indicator.
77  itfac = 1
78  CALL derqp4(ym(1),ym(13),dy(1),dy(13))
79  kl=l.LT.2
80  gr=l.GT.5
81  fs=0.
82  DO i=1,n
83  v=dt(i,1)
84  c=(ym(i)+yl(i)+g*dy(i))*0.5d0
85  dt(i,1)=c
86  ta=c
87  IF(.NOT.kl)THEN
88  DO k=2,l
89  b1=d(k)*v
90  b=b1-c
91  w=c-v
92  u=v
93  IF(b.NE.0.d0)THEN
94  b=w/b
95  u=c*b
96  c=b1*b
97  END IF
98  v=dt(i,k)
99  dt(i,k)=u
100  ta=u+ta
101  END DO
102  is=i+n/2
103  IF(is.GT.nhalf2)is=i-(n/2)
104 * Use modified EPS for convergence test.
105  dyis=dabs(dy(is))/tfac
106  IF(i.GT.nhalf2)dyis=0.0
107  IF(konv)THEN
108  test=dabs( (y(i)-ta)*dyis )
109  IF(test.GT.eps) konv=.false.
110  END IF
111  IF(.NOT.gr)THEN
112  fv=dabs(w)*dyis
113  IF(fs.LT.fv) fs=fv
114  END IF
115  END IF
116  y(i)=ta
117  END DO
118  IF(fs.NE.0.d0)THEN
119  fa=fy
120  k=l-1
121  fy=(ep(k)/fs)**(1./(l+k))
122  fa7=0.7*fa
123  IF(l.EQ.2)fa7=0.0
126  h=h*fy
127  jti=jti+1
128  IF(jti.GT.5)THEN
129  h=0.0
130  DO i=1,n
131  y(i)=ya(i)
132  END DO
133  RETURN
134  END IF
135  goto 10
136  END IF
137  END IF
138  IF(konv)THEN
139  x=xn
140  h=h*fy
141  RETURN
142  END IF
143  d(3)=4.d0
144  d(5)=1.6d1
145  bo=.NOT.bo
146  m=jr
147  jr=js
148  js=m+m
149  END DO
150  h=0.5*h
151  goto 10
152  END