Nbody6
shrink.f
Go to the documentation of this file.
1  SUBROUTINE shrink(TMIN)
2 *
3 *
4 * Shrinking of large time-steps.
5 * ------------------------------
6 *
7  include 'common6.h'
8  parameter(dtm = 0.03125d0)
9  REAL*8 rij(3),vij(3)
10 *
11 *
12 * Determine impact to each high-velocity particle for large STEPR.
13  nnb = listv(1)
14  DO 20 l = 1,nnb
15  i = listv(l+1)
16  DO 10 j = ifirst,ntot
17 * Form relative quantities for large regular time-steps.
18  IF (stepr(j).LT.dtm.OR.j.EQ.i) go to 10
19 * Check whether body #I is already a neighbour (accept LIST > I).
20  nnb1 = list(1,j) + 1
21  DO 1 k = 2,nnb1
22  IF (list(k,j).GT.i) go to 2
23  IF (list(k,j).EQ.i) go to 10
24  1 CONTINUE
25  2 rv = 0.0
26  vv = 0.0
27  DO 3 k = 1,3
28  rij(k) = x(k,i) - x(k,j)
29  vij(k) = xdot(k,i) - xdot(k,j)
30  rv = rv + rij(k)*vij(k)
31  vv = vv + vij(k)**2
32  3 CONTINUE
33 *
34 * Skip treatment for receding particles.
35  IF (rv.GE.0.0d0) go to 10
36 *
37 * Evaluate time of minimum approach and truncate to next time.
38  dt = -rv/vv
39  it = 0
40  5 dt = min(t0r(j) + stepr(j) - time,dt)
41 *
42 * Obtain minimum impact parameter for straight-line orbit.
43  r2 = 0.0
44  fj2 = 0.0
45  DO 6 k = 1,3
46  r2 = r2 + (rij(k) + vij(k)*dt)**2
47  fj2 = fj2 + f(k,j)**2
48  6 CONTINUE
49 *
50 * Compare force at minimum distance with total force on body #J (> 0).
51  fi2 = (body(i)/r2)**2
52  IF (fi2.LT.0.04*fj2.OR.body(j).EQ.0.0d0) go to 10
53 *
54 * See whether the regular time-step can be shortened.
55  IF (t0r(j) + 0.5*stepr(j).GE.tmin.AND.it.LT.5) THEN
56  WRITE (29,8) j, sqrt(r2), sqrt(fi2/fj2), dt, stepr(j)
57  8 FORMAT (' SHRINK J R FI/FJ DT STEP ',
58  & i6,f7.3,f6.2,1p,2e10.2)
59  CALL flush(29)
60  stepr(j) = 0.5*stepr(j)
61  IF (step(j).GT.stepr(j)) THEN
62  step(j) = 0.5*step(j)
63  tnew(j) = t0(j) + step(j)
64  END IF
65  it = it + 1
66  go to 5
67  END IF
68  10 CONTINUE
69  20 CONTINUE
70 *
71  RETURN
72 *
73  END