Nbody6
hcorr.f
Go to the documentation of this file.
1  SUBROUTINE hcorr(I,DM,RNEW)
2 *
3 *
4 * Mass loss correction of KS orbit.
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),iflag(mmax)
11 *
12 *
13 * Copy pair index and determine secondary component (note I1 = I here).
14  ipair = kspair
15  i1 = i
16  i2 = 2*ipair - 1
17  IF (i2.EQ.i1) THEN
18  i2 = i2 + 1
19  END IF
20 *
21 * Define c.m. index and old semi-major axis.
22  j = n + ipair
23  semi0 = -0.5d0*body(j)/h(ipair)
24 *
25 * Distinguish different cases (KS, triple or outer ghost of quadruple).
26  IF (body(i).GT.0.0d0) THEN
27  zmu0 = body(i1)*body(i2)/body(j)
28  ELSE
29  im = 0
30  jj = i
31  IF (i.LT.ifirst) jj = n + kvec(i)
32 * Determine merger index for inner and outer binary masses.
33  DO 1 k = 1,nmerge
34  IF (nameg(k).EQ.name(jj)) im = k
35  1 CONTINUE
36  IF (im.EQ.0) THEN
37  WRITE (6,5) name(i), name(jj), (nameg(k),k=1,nmerge)
38  5 FORMAT (' DANGER HCORR! NMI NMJJ NMG ',12i6)
39  stop
40  END IF
41 * Form old reduced mass from inner and outer binary.
42  zm1 = cm(1,im) + cm(2,im)
43  zm2 = cm(3,im) + cm(4,im)
44  zmu0 = zm1*zm2/body(j)
45  END IF
46 *
47 * Obtain energy change from MK*A = const and H = -M/(2*A) (DCH 8/96).
48  dh = dm/semi0*(1.0 - 0.5*dm/body(j))
49 *
50 * Reduce mass of c.m. body and subtract energy correction terms.
51  body(j) = body(j) - dm
52  IF (body(i).GT.0.0d0) THEN
53  zmu1 = (body(i1) - dm)*body(i2)/body(j)
54  ELSE
55  zmu1 = zm1*(zm2 - dm)/body(j)
56  END IF
57  emdot = emdot - zmu1*dh - (zmu1 - zmu0)*h(ipair)
58 * Retain binding energy gained in past interactions.
59  egrav = egrav - zmu1*dh - (zmu1 - zmu0)*h(ipair)
60 *
61 * Include diagnostic warning for large relative energy change.
62  IF (dh.GT.0.2*abs(h(ipair))) THEN
63  WRITE (6,10) name(i), dh, h(ipair), r(ipair)/semi0, dm*smu
64  10 FORMAT (' WARNING! LARGE CORRECTION NM DH H R/A DMS ',
65  & i6,2f8.2,2f7.2)
66  END IF
67 *
68 * Update the binding energy due to mass loss DM and set new radius.
69  h(ipair) = h(ipair) + dh