Nbody6
 All Files Functions Variables
mloss.f
Go to the documentation of this file.
1  SUBROUTINE mloss
2 *
3 *
4 * Mass loss from evolving stars.
5 * ------------------------------
6 *
7 * Original scheme of Elena Terlevich 1983.
8 * ----------------------------------------
9 *
10  include 'common6.h'
11  REAL*8 a(9)
12 *
13 *
14 * Find the heaviest body (exclude any merged binary components).
15  99 imax = 0
16  1 bmax = 0.0d0
17  DO 2 i = 1,n
18  IF (body(i).LT.bmax.OR.i.LE.imax) go to 2
19  bmax = body(i)
20  imax = i
21  2 CONTINUE
22 *
23  zmstar = zmbar*bmax
24 * Obtain total evolution time for ZMSTAR in solar masses.
25  tms = (2.55d+03 + 6.69d2*zmstar**2.5d0 + zmstar**4.5d0)/
26  & (3.27d-02*zmstar**1.5d0 + 3.46d-01*zmstar**4.5d0)
27  tg = 0.15*tms
28  the = tms*1.37*zmstar**(-0.881d0)
29  tev1 = tms + tg + the
30 * Scale the evolution time to model units.
31  tmdot = tev1/tscale
32 *
33  IF (tmdot.GT.time) THEN
34 * Update the maximum mass.
35  body1 = bmax
36 * Set phase indicator = -1 for new time-step list in INTGRT.
37  iphase = -1
38  RETURN
39  END IF
40 *
41  ibody = imax
42 * Include special treatment for regularized components.
43  IF (imax.GT.2*npairs) go to 4
44  ipair = kvec(imax)
45  ibody = n + ipair
46 * Continue search if merged binary component is identified.
47  IF (name(ibody).LT.0) go to 1
48 *
49 * Obtain current coordinates & velocities (unperturbed KS pair OK).
50  CALL resolv(ipair,2)
51 *
52 * Copy X0DOT since routine RESOLV skipped in KSTERM (IPHASE = -1).
53  DO 3 k = 1,3
54  x0dot(k,2*ipair-1) = xdot(k,2*ipair-1)
55  x0dot(k,2*ipair) = xdot(k,2*ipair)
56  3 CONTINUE
57 *
58 * Set variable white dwarf mass (Iben & Renzini, Ann. Rev. 21, 298).
59  4 zmstar = 0.38 + 0.15*zmstar
60 * Assume neutron star instead if mass > 6 solar masses.
61  IF (zmbar*bmax.GT.6.0) zmstar = 1.5
62  dm = body(imax) - zmstar/zmbar
63  zmass = zmass - dm
64  zmdot = zmdot + dm*zmbar
65  nmdot = nmdot + 1
66  body(imax) = zmstar/zmbar
67  IF (tidal(1).NE.0.0d0) rtide = (zmass/tidal(1))**0.3333
68 *
69 * Save old velocity for FDOT correction.
70  vi2 = 0.0
71  DO 5 k = 1,3
72  a(k+6) = xdot(k,imax)
73  vi2 = vi2 + xdot(k,imax)**2
74  5 CONTINUE
75 *
76  de = 0.5*dm*(vi2 - tidal(1)*x(1,imax)**2 - tidal(3)*x(3,imax)**2)
77  a0 = 1.0
78  IF (zmbar*bmax.LE.6.0) go to 8
79 *
80 * Assign a high velocity to neutron star (at least 4*rms velocity).
81  a0 = 4.0*sqrt(0.5d0*zmass/(rscale*vi2))
82  IF (imax.GT.2*npairs) go to 6
83  a2 = 2.0*(body(2*ipair-1) + body(2*ipair))/r(ipair)
84 * Add escape velocity from the regularized pair.
85  a0 = sqrt(a0**2 + a2/vi2)
86  6 DO 7 k = 1,3
87  xdot(k,imax) = a0*xdot(k,imax)
88  x0dot(k,imax) = xdot(k,imax)
89  7 CONTINUE
90 *
91 * Correct total energy, forces & first derivatives.
92  8 poti = 0.0d0
93  DO 20 j = 1,ntot
94  IF (j.EQ.imax) go to 20
95  rijdot = 0.0d0
96  rdvdot = 0.0d0
97 *
98  DO 10 k = 1,3
99  a(k) = x(k,imax) - x(k,j)
100  a(k+3) = a(k+6) - xdot(k,j)
101  rijdot = rijdot + a(k)*a(k+3)
102  rdvdot = rdvdot + a(k)*(xdot(k,imax) - a(k+6))
103  10 CONTINUE
104 *
105  rij2 = a(1)**2 + a(2)**2 + a(3)**2
106  IF (j.LE.n) poti = poti + body(j)/sqrt(rij2)
107  IF (j.LE.2*npairs.OR.j.EQ.ibody) go to 20
108  a3 = 1.0/(rij2*sqrt(rij2))
109  a4 = body(imax)*a3
110  a5 = dm*a3
111  a6 = 3.0*rijdot/rij2
112  a7 = 3.0*rdvdot/rij2
113 *
114  DO 11 k = 1,3
115  a(k+3) = (a(k+3) - a(k)*a6)*a5
116  IF (a0.GT.1.0) THEN
117 * Include FDOT corrections due to increased velocity.
118  a(k+3) = a(k+3) + (xdot(k,imax) - a(k+6))*a4
119  a(k+3) = a(k+3) - a7*a(k)*a4
120  END IF
121  11 CONTINUE
122 *
123 * Use neighbour list to distinguish irregular & regular terms.
124  nnb = list(1,j) + 1
125  DO 12 l = 2,nnb
126  IF (list(l,j).EQ.ibody) go to 16
127  IF (list(l,j).GT.ibody) go to 13
128  12 CONTINUE
129 *
130  13 DO 14 k = 1,3
131  f(k,j) = f(k,j) - 0.5*a(k)*a5
132  fr(k,j) = fr(k,j) - a(k)*a5
133  fdot(k,j) = fdot(k,j) - one6*a(k+3)
134  d1r(k,j) = d1r(k,j) - a(k+3)
135  14 CONTINUE
136  go to 20
137 *
138  16 DO 18 k = 1,3
139  f(k,j) = f(k,j) - 0.5*a(k)*a5
140  fi(k,j) = fi(k,j) - a(k)*a5
141  fdot(k,j) = fdot(k,j) - one6*a(k+3)
142  d1(k,j) = d1(k,j) - a(k+3)
143  18 CONTINUE
144  20 CONTINUE
145 *
146 * Correct total energy for mass loss effect.
147  de = de - dm*poti
148  be(3) = be(3) - de
149  ri = sqrt((x(1,imax) - rdens(1))**2 + (x(2,imax) - rdens(2))**2 +
150  & (x(3,imax) - rdens(3))**2)
151  WRITE (6,30) name(imax), bmax, bmax*zmbar, zmdot, de, be(3),
152  & tmdot, ri, list(1,imax)
153  30 FORMAT (/,' MASS LOSS NAME =',i5,' MI =',f8.4,' M* =',f5.1,
154  & ' ZM* =',f6.1,' DE =',f9.5,' E =',f10.6,
155  & ' TMDOT =',f6.1,' RI =',f5.2,' NNB =',i3)
156 *
157 * Set option = 2 to skip next energy check (reduced by 1 in ADJUST).
158  kz(19) = 2
159 *
160  IF (zmbar*bmax.LE.6.0) go to 50
161 * Include correction to total energy from the increased velocity.
162  de = 0.5*body(imax)*vi2*(a0**2 - 1.0)
163  be(3) = be(3) + de
164  WRITE (6,45) imax, de, tev1, a0, sqrt(vi2), be(3)
165  45 FORMAT (' RECOIL I =',i5,' DE =',f9.5,' T* =',f5.1,
166  & ' V/V0 =',f6.2,' V0 =',f5.2,' E =',f10.6)
167 *
168  50 IF (imax.LE.2*npairs) THEN
169 * Save pair index in KSPAIR and terminate regularization.
170  kspair = ipair
171  iphase = -1
172 * Indicator for skipping routine RESOLV to preserve velocity X0DOT.
173  CALL ksterm
174  ELSE
175 * Reduce the steps by velocity factor.
176  step(imax) = max(step(imax)/a0,time - t0(imax))
177  stepr(imax) = max(stepr(imax)/a0,step(imax))
178  END IF
179 *
180 * Set next mass loss time & current maximum mass before returning.
181  go to 99
182 *
183  END