Nbody6
 All Files Functions Variables
decide.f
Go to the documentation of this file.
1  SUBROUTINE decide(IPAIR,SEMI,ECC,EMAX,EMIN,TC,TG,EDAV,IQ)
2 *
3 *
4 * Merger decision.
5 * ----------------
6 *
7  include 'common6.h'
8  common/binary/ cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
9  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
10  & namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
11 *
12 *
13 * Initialize the merger skip indicator and define KS indices.
14  iq = 0
15  i1 = 2*ipair - 1
16  i2 = i1 + 1
17 * Prevent over-writing in NMERGE+1.
18  IF(nmerge.GE.mmax-1)THEN
19  iq = 1
20  goto 40
21  ENDIF
22 *
23 * Accept active spiral but ensure some updating to increase pericentre.
24  dmr = 0.d0
25  IF (kstar(n+ipair).EQ.-2) THEN
26  pmin = semi*(1.0 - ecc)
27  icirc = -1
28  CALL tcirc(pmin,ecc,i1,i2,icirc,tc1)
29 * Restrict circularization time to account for stellar evolution.
30  tc1 = min(tc1,500.0d0)
31 * Adopt termination interval 1/2 of predicted t_{circ} (at const R*).
32  dt1 = 0.5*(tc1 + 0.1)/tstar
33  IF (ecc.LT.0.1.AND.max(kstar(i1),kstar(i2)).LE.1) THEN
34  dt1 = 10.0*dt1
35  END IF
36 * Rectify SPIRAL just in case.
37  IF(list(1,i1).GT.0)THEN
38 * IQ = 1
39 * GO TO 30
40 * Accept perturbed spiral case for now!
41  dmr = -1.d0
42  CALL chrect(ipair,dmr)
43  ELSE
44  CALL chrect(ipair,dmr)
45  ENDIF
46  go to 30
47  END IF
48 *
49 * Obtain eccentricity derivative due to #JCOMP (KS resolved in INDUCE).
50  CALL edot(i1,i2,jcomp,semi,ecc,eccdot)
51 *
52 * Define c.m. index and set maximum stellar radius and peri.
53  i = n + ipair
54  rm = max(radius(i1),radius(i2))
55  pm = semi*(1.0 - ecc)/rm
56  pm = min(pm,99.9d0)
57 *
58 * WRITE (75,1) KSTAR(I1), KSTAR(I2), KSTAR(I), LIST(1,I1), NAME(I),
59 * & TTOT, ECC, EMIN, EMAX, TG, EDAV, PM
60 * 1 FORMAT (' DECIDE: K* NP NM T E EM EX TG ED PM ',
61 * & 4I4,I6,F10.3,F8.4,2F6.2,2F7.2,F6.1)
62 * CALL FLUSH(75)
63 *
64 * See whether to deform binary orbit (skip chaos case).
65  IF (kstar(i).NE.-1.AND.tg.LT.-10.0) THEN
66 * Change eccentricity slowly by a fraction 0.1*(1 - E) of decay time.
67  ecc0 = ecc
68  dt1 = 0.1*(1.0 - ecc)*tg/tstar
69  ecc = ecc + edav*dt1
70  ecc = max(ecc,emin)
71 * Impose temporary limit of 0.99*EMAX to allow for apsidal motion.
72  ecc = min(ecc,0.99*emax)
73 *
74 * Evaluate t_{circ} and (de/dt)_{circ} except for near-collision.
75  pmin = semi*(1.0 - ecc)
76  IF (pmin.GT.rm) THEN
77  icirc = -1
78  CALL tcirc(pmin,ecc,i1,i2,icirc,tc1)
79 * CALL EDOT(IPAIR,JCOMP,SEMI,ECC,ECCDOT)
80  CALL ecirc(pmin,ecc,i1,i2,icirc,tg,tc2,ecc2,edt)
81 * WRITE (6,10) ECCDOT, R(IPAIR), SEMI, TDOT2(IPAIR)
82 * 10 FORMAT (' ECIRC: ECDOT R A TD2 ',1P,4E10.2)
83  ELSE
84  tc1 = tg
85  tc2 = tc1
86  edt = 2.0*edav
87  END IF
88 *
89 * Compare eccentricity derivatives and check t_{circ} (standard case).
90  IF (abs(edt).GT.abs(edav).AND.tc1.LT.50.0) THEN
91 * Enforce new spiral on short circularization time (IQ > 0: no merger).
92  IF (kstar(i).GE.0) iq = 1
93  ecc = ecc0
94  IF (kstar(i).EQ.-2) go to 40
95  go to 30
96  END IF
97 *
98 * Rectify spiral orbit before changing eccentricity.
99  IF (kstar(i).EQ.-2) THEN
100  CALL chrect(ipair,dmr)
101  END IF
102 *
103 * Transform to exact apocentre unless done already.
104  IF (abs(tdot2(ipair)).GT.1.0d-12) THEN
105  IF (r(ipair).GT.semi.AND.tdot2(ipair).LT.0.0) THEN
106  CALL ksapo(ipair)
107  END IF
108  CALL ksperi(ipair)
109  CALL ksapo(ipair)
110 * Check spiral case already at peri.
111  ELSE IF (kstar(i).EQ.-2.AND.r(ipair).LT.semi.AND.
112  & abs(tdot2(ipair)).LT.1.0d-12) THEN
113  CALL ksapo(ipair)
114  END IF
115 *
116 * Deform orbit (H = const at apo), rectify and transform to X & XDOT.
117  CALL deform(ipair,ecc0,ecc)
118  CALL ksrect(ipair)
119  CALL resolv(ipair,1)
120 *
121 * Rectify spiral orbit (terminate on collision).
122  IF (kstar(i).EQ.-2) THEN
123  CALL chrect(ipair,dmr)
124  IF (iphase.LT.0) THEN
125  iq = 2
126  go to 40
127  END IF
128  END IF
129 *
130 * Re-initialize KS polynomials for perturbed motion (also in DEFORM).
131  IF (list(1,i1).GT.0) THEN
132  IF (r(ipair).LT.semi.AND.
133  & abs(tdot2(ipair)).LT.1.0d-12) THEN
134  CALL ksapo(ipair)
135  END IF
136  t0(i1) = time
137  imod = kslow(ipair)
138  CALL kspoly(ipair,imod)
139  tdot2(ipair) = -1.0d-20
140  END IF
141 *
142 * Produce diagnostic output for large eccentricity.
143  IF (ecc.GT.0.9) THEN
144  WRITE (75,20) name(i1), ttot, ecc0, ecc, emin, emax,
145  & eccdot, edt, tg, tc1, edav, pm
146  20 FORMAT (' DECIDE: NM T E0 E1 EM EX ED EDT TG TC EDAV ',
147  & 'PM ', i5,f9.2,4f7.3,1p,6e9.1)
148  CALL flush(75)
149  END IF
150  ELSE
151  dt1 = 10.0/tstar
152  END IF
153 *
154 * See whether growth time exceeds merger disruption time from TSTAB.
155  30 tg1 = time + dt1
156  IF (iq.EQ.0) THEN
157  tmdis(nmerge+1) = min(tg1,tmdis(nmerge+1))
158 * WRITE (6,55) TIME, NAME(I1), NAME(JCOMP), ECC, EMAX,DT1, TG, SEMI
159 * 55 FORMAT (' WATCH! T NM NMJ E EX DT1 TG A ',
160 * & F10.4,2I7,2F8.4,1P,3E10.2)
161  END IF
162 * WRITE (6,35) LIST(1,I1), ECC, EDAV*DT1, EMAX, TG1+TOFF, GAMMA(IPAIR)
163 * 35 FORMAT (' DECIDE: NP E DE EMAX TG1 G ',I4,3F7.3,F10.3,1P,E9.1)
164 *
165  40 RETURN
166 *
167  END