Nbody6
 All Files Functions Variables
hipop.f
Go to the documentation of this file.
1  SUBROUTINE hipop
2 *
3 *
4 * Initial hierarchical population.
5 * --------------------------------
6 *
7  include 'common6.h'
8  REAL*8 xorb(2),vorb(2),xrel(3),vrel(3),px(3),qx(3),
9  & bs(mmax),xs(3,mmax),vs(3,mmax)
10  REAL*8 ran2
11 *
12 *
13 * Read input parameters (same usage as routine BINPOP).
14  READ (5,*) semi0, ecc0, ratio, range, icirc
15  nhi = nhi0
16  WRITE (6,1) nhi, semi0, ecc0, ratio, range, icirc
17  1 FORMAT (/,12x,'HIERARCHIES: NHI =',i4,' A =',f9.6,
18  & ' E =',f6.2,' RATIO =',f4.1,' RANGE =',f6.1,
19  & ' ICIRC =',i2,/)
20 *
21  IF (nhi.GT.nbin0) THEN
22  WRITE (6,2) nbin0, nhi
23  2 FORMAT (5x,'FATAL ERROR! NBIN0 =',i4,' NHI =',i4)
24  stop
25  END IF
26 *
27 * Introduce binary components by splitting the primary.
28  DO 40 i = 1,nhi
29 *
30 * Randomize perihelion, node & inclination.
31  pi = twopi*ran2(idum1)
32  omega = twopi*ran2(idum1)
33  zi = 0.5*twopi*ran2(idum1)
34 *
35 * Set transformation elements (Brouwer & Clemence p. 35).
36  px(1) = cos(pi)*cos(omega) - sin(pi)*sin(omega)*cos(zi)
37  qx(1) =-sin(pi)*cos(omega) - cos(pi)*sin(omega)*cos(zi)
38  px(2) = cos(pi)*sin(omega) + sin(pi)*cos(omega)*cos(zi)
39  qx(2) =-sin(pi)*sin(omega) + cos(pi)*cos(omega)*cos(zi)
40  px(3) = sin(pi)*sin(zi)
41  qx(3) = cos(pi)*sin(zi)
42 *
43 * Determine two-body elements for original binary.
44  i1 = 2*i - 1
45  i2 = 2*i
46  rij2 = 0.0
47  vij2 = 0.0
48  rdot = 0.0
49  DO 5 k = 1,3
50  xrel(k) = x(k,i1) - x(k,i2)
51  vrel(k) = xdot(k,i1) - xdot(k,i2)
52  rij2 = rij2 + xrel(k)**2
53  vij2 = vij2 + vrel(k)**2
54  rdot = rdot + xrel(k)*vrel(k)
55  5 CONTINUE
56  rij = sqrt(rij2)
57  zmb1 = body(i1) + body(i2)
58  a1 = 2.0/rij - vij2/zmb1
59  semi1 = 1.0/a1
60  ecc2 = (1.0 - rij/semi1)**2 + rdot**2/(semi1*zmb1)
61  ecc1 = sqrt(ecc2)
62  pmin = semi1*(1.0 - ecc1)
63 *
64 * Specify component masses (primary fraction range 0.5 - 0.9).
65  q0 = 0.5 + 0.4*ran2(idum1)
66  IF (ratio.EQ.1.0) q0 = 0.5
67  bs(i) = body(i2)
68  body(i1) = q0*body(i1)
69  body(i2) = body(i1)*(1.0 - q0)/q0
70  zmb = body(i1) + body(i2)
71 *
72 * Choose random (thermalized) or fixed eccentricity.
73  IF (ecc0.LT.0.0) THEN
74  ecc2 = ran2(idum1)
75  ecc = sqrt(ecc2)
76  ELSE
77  ecc = ecc0
78  END IF
79 *
80 * Select semi-major axis from uniform distribution in log(A) or SEMI0.
81  iter = 0
82  10 IF (range.GT.0.0) THEN
83  IF (iter.LE.5) THEN
84  exp = ran2(idum1)*log10(range)
85  semi = semi1/10.0**exp
86  ELSE
87 * Shrink by factor 2 if no success after 5 iterations.
88  semi = 0.5*semi
89  END IF
90  ELSE
91  semi = semi0
92  END IF
93 *
94 * Check stability criterion (maximum 12 tries with inclination effect).
95  pcrit = stability(body(i1),body(i2),bs(i),ecc,ecc1,zi)
96  pcrit = pcrit*semi
97  iter = iter + 1
98  IF (pmin.LT.pcrit.AND.iter.LE.12.AND.semi.LT.semi0) go to 10
99 *
100  p0 = days*semi*sqrt(semi/zmb)
101  p1 = days*semi1*sqrt(semi1/zmb1)
102  WRITE (6,20) iter, ecc, ecc1, pmin, pcrit, p0, p1
103  20 FORMAT (' HIERARCHY: IT E E1 PMIN PCRIT P0 P1 ',
104  & i4,2f7.3,1p,4e9.1)
105 *
106 * Specify relative motion at apocentre and sum binding energy.
107  xorb(1) = semi*(1.0 + ecc)
108  xorb(2) = 0.0
109  vorb(1) = 0.0
110  vorb(2) = sqrt(zmb*(1.0d0 - ecc)/(semi*(1.0d0 + ecc)))
111  ebin0 = ebin0 - 0.5*body(i1)*body(i2)/semi
112 *
113 * Transform to relative variables.
114  DO 25 k = 1,3
115  xrel(k) = px(k)*xorb(1) + qx(k)*xorb(2)
116  vrel(k) = px(k)*vorb(1) + qx(k)*vorb(2)
117  25 CONTINUE
118 *
119 * Save old secondary and set global variables for each component.
120  DO 30 k = 1,3
121  xs(k,i) = x(k,i2)
122  vs(k,i) = xdot(k,i2)
123  x(k,i1) = x(k,i1) + body(i2)*xrel(k)/zmb
124  x(k,i2) = x(k,i1) - xrel(k)
125  xdot(k,i1) = xdot(k,i1) + body(i2)*vrel(k)/zmb
126  xdot(k,i2) = xdot(k,i1) - vrel(k)
127  30 CONTINUE
128  40 CONTINUE
129 *
130 * Move single particles down by NHI to make room for outer components.
131  l = n
132  DO 50 j = 2*nbin0+1,n
133  i = nhi + l
134  body(i) = body(l)
135  DO 45 k = 1,3
136  x(k,i) = x(k,l)
137  xdot(k,i) = xdot(k,l)
138  45 CONTINUE
139  l = l - 1
140  50 CONTINUE
141 *
142 * Place hierarchical components immediately after the binaries.
143  DO 60 l = 1,nhi
144  i = 2*nbin0 + l
145  body(i) = bs(l)
146  DO 55 k = 1,3
147  x(k,i) = xs(k,l)
148  xdot(k,i) = vs(k,l)
149  55 CONTINUE
150  60 CONTINUE
151 *
152 * Update the particle number and reset NHI (might be used elsewhere).
153  n = n + nhi
154  nzero = n
155  ntot = n
156  nhi = 0
157 *
158  RETURN
159 *
160  END