Nbody6
 All Files Functions Variables
chinit.f
Go to the documentation of this file.
1  SUBROUTINE chinit(ISUB)
2 *
3 *
4 * Initialization of chain system.
5 * -------------------------------
6 *
7  include 'common6.h'
8  parameter(nmx=10,nmx3=3*nmx,nmx4=4*nmx,nmxm=nmx*(nmx-1)/2)
9  REAL*8 m,mass,mc,mij,mkk,ang(3),firr(3),fd(3)
10  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
11  & zz(nmx3),wc(nmx3),mc(nmx),
12  & xi(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
13  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
14  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
15  & listc(lmax)
16  common/cpert/ rgrav,gpert,ipert,npert
17  common/chreg/ timec,tmax,rmaxc,cm(10),namec(6),nstep1,kz27,kz30
18  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
19  & names(ncmax,5),isys(5)
20  common/ccoll2/ qk(nmx4),pk(nmx4),rik(nmx,nmx),SIZE(nmx),vstar1,
21  & ecoll1,rcoll,qperi,istar(nmx),icoll,isync,ndiss1
22  common/incond/ x4(3,nmx),xdot4(3,nmx)
23  common/echain/ ech
24  common/slow3/ gcrit,kz26
25 *
26 *
27 * Define chain membership.
28  CALL setsys
29 *
30 * Initialize c.m. variables.
31  DO 2 k = 1,7
32  cm(k) = 0.0d0
33  2 CONTINUE
34 *
35 * Transform to the local c.m. reference frame.
36  DO 4 l = 1,nch
37  j = jlist(l)
38  SIZE(l) = radius(j)
39  istar(l) = kstar(j)
40 * Place the system in first single particle locations.
41  cm(7) = cm(7) + m(l)
42  DO 3 k = 1,3
43  x4(k,l) = x(k,j)
44  xdot4(k,l) = xdot(k,j)
45  cm(k) = cm(k) + m(l)*x4(k,l)
46  cm(k+3) = cm(k+3) + m(l)*xdot4(k,l)
47  3 CONTINUE
48  4 CONTINUE
49 *
50 * Set c.m. coordinates & velocities of subsystem.
51  DO 5 k = 1,6
52  cm(k) = cm(k)/cm(7)
53  5 CONTINUE
54 *
55 * Specify initial conditions for chain regularization.
56  lk = 0
57  DO 8 l = 1,nch
58  DO 7 k = 1,3
59  lk = lk + 1
60  x4(k,l) = x4(k,l) - cm(k)
61  xdot4(k,l) = xdot4(k,l) - cm(k+3)
62  xch(lk) = x4(k,l)
63  vch(lk) = xdot4(k,l)
64  7 CONTINUE
65  8 CONTINUE
66 *
67 * Calculate internal energy and and save in chain energy.
68  CALL const(xch,vch,m,nch,energy,ang,gam)
69  ech = energy
70 *
71 * Set sum of mass products and save separations & RINV for CHLIST.
72  sum = 0.0d0
73  rsum = 0.0d0
74  DO 10 l = 1,nch-1
75  DO 9 k = l+1,nch
76  sum = sum + m(l)*m(k)
77  rlk2 = (x4(1,l) - x4(1,k))**2 + (x4(2,l) - x4(2,k))**2 +
78  & (x4(3,l) - x4(3,k))**2
79  rsum = rsum + sqrt(rlk2)
80  rinv(l) = 1.0/sqrt(rlk2)
81  9 CONTINUE
82  10 CONTINUE
83 *
84 * Reduce RSUM by geometrical factor and check upper limit from IMPACT.
85  IF (nch.EQ.4) rsum = 0.5*rsum
86  rsum = min(float(nch-1)*rsum/float(nch),rmin)
87 *
88 * Define gravitational radius for initial perturber list.
89  rgrav = sum/abs(energy)
90 *
91 * Avoid small value after collision (CHTERM improves perturbers).
92  IF (nch.GT.2) THEN
93  rgrav = min(rgrav,0.5*rsum)
94  END IF
95 *
96 * Set global index of c.m. body and save name (SUBSYS sets NAME = 0).
97  IF (timec.GT.0.0d0) ich0 = ich
98  ich = jlist(1)
99  name0 = name(ich)
100 *
101 * Define subsystem indicator (ISYS = 1, 2, 3 for triple, quad, chain).
102  isys(nsub+1) = 3
103 *
104 * Form ghosts and initialize c.m. motion in ICOMP (= JLIST(1)).
105  CALL subsys(nch,cm)
106 *
107 * Copy neighbour list for ghost removal.
108  nnb = list(1,ich)
109  DO 20 l = 2,nnb+1
110  jpert(l-1) = list(l,ich)
111  20 CONTINUE
112 *
113 * Check possible switch of reference body on second call from CHAIN.
114  IF (timec.GT.0.0d0.AND.ich.NE.ich0) THEN
115 * Add #ICH to neighbour & perturber lists before removing all ghosts.
116  CALL nbrest(ich0,1,nnb)
117  END IF
118 *
119 * Remove ghosts (saved in JLIST) from neighbour lists of #ICH.
120  CALL nbrem(ich,nch,nnb)
121 *
122 * Remove ghosts from list of ICOMP (use NTOT as dummy here).
123  jpert(1) = icomp
124  CALL nbrem(ntot,nch,1)
125 *
126 * Initialize perturber list for integration of chain c.m.
127  CALL chlist(ich)
128 *
129 * Perform differential F & FDOT corrections due to perturbers.
130  DO 25 k = 1,3
131  firr(k) = 0.0d0
132  fd(k) = 0.0
133  25 CONTINUE
134  CALL chfirr(ich,0,x(1,ich),xdot(1,ich),firr,fd)
135  DO 30 k = 1,3
136  f(k,ich) = f(k,ich) + 0.5*firr(k)
137  fdot(k,ich) = fdot(k,ich) + one6*fd(k)
138  d1(k,ich) = d1(k,ich) + fd(k)
139  30 CONTINUE
140 *
141 * Take maximum integration interval equal to c.m. step.
142  tmax = step(ich)
143 *
144 * Check next treatment time of perturbers.
145  CALL tchain(nsub,tsmin)
146  tmax = min(tmax,tsmin)
147 *
148 * Copy binding energy and output & capture option for routine CHAIN.
149  cm(8) = e(3)
150  kz26 = kz(26)
151  kz27 = kz(27)
152  kz30 = kz(30)
153 * Copy velocity scale factor to VSTAR1.
154  vstar1 = vstar
155 *
156 * Assign new subsystem index and begin chain regularization.
157  isub = nsub
158  nchain = nchain + 1
159 *
160 * Set phase indicator < 0 to ensure new time-step list in INTGRT.
161  iphase = -1
162 *
163  RETURN
164 *
165  END