Nbody6
 All Files Functions Variables
hrplot.f
Go to the documentation of this file.
1  SUBROUTINE hrplot
2 *
3 *
4 * HR diagnostics of evolving stars.
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  REAL*8 lums(10),tscls(20),gb(10)
12  REAL*8 m0,m1,m2,lum,lum2,mc,me,k2
13 *
14 *
15  WRITE (82,1) npairs, tphys
16  1 FORMAT (' ## BEGIN',i8,f9.1)
17 * Define the number of objects (rather than single stars).
18  ns = n - npairs - nmerge - (nch - 1)
19 * Choose the number of singles (alternative definition, triples only).
20 * NS = N - 2*NPAIRS - NMERGE - NCH
21  imerge = 0
22  nb = 0
23  nstar = 0
24  WRITE (83,1) ns, tphys
25 *
26  DO 20 i = 1,n
27  m0 = body0(i)*zmbar
28  m1 = body(i)*zmbar
29 * Replace ghost mass of single star with original value from merged KS.
30  IF (m1.EQ.0.0.AND.i.GE.ifirst) THEN
31  im = 0
32 * Search merger table for ghost to identify corresponding index.
33  DO 2 k = 1,nmerge
34  IF (nameg(k).EQ.name(i)) THEN
35  im = k
36  END IF
37  2 CONTINUE
38 * Skip any ghosts associated with chain regularization.
39  IF (im.EQ.0) THEN
40  WRITE (6,3) i, nch
41  3 FORMAT (' WARNING! HRPLOT I NCH ',i6,i4)
42  go to 20
43  END IF
44  m1 = cm(2,im)*zmbar
45  END IF
46 *
47 * Obtain stellar parameters at current epoch.
48  kw = kstar(i)
49  age = max(tplot,tev0(i))*tstar - epoch(i)
50  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
51  CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
52  & rm,lum,kw,mc,rcc,me,re,k2)
53  IF (i.LT.ifirst) THEN
54  jpair = kvec(i)
55  j2 = 2*jpair
56  IF (i.EQ.j2) go to 20
57  j1 = 2*jpair - 1
58  icm = n + jpair
59  m2 = body(j2)*zmbar
60  ri = (x(1,icm) - rdens(1))**2 + (x(2,icm) - rdens(2))**2 +
61  & (x(3,icm) - rdens(3))**2
62 * Check for ghost binary.
63  IF (m1.EQ.0.0) THEN
64  im = 0
65 * Search merger table to identify corresponding index of c.m. name.
66  DO 4 k = 1,nmerge
67  IF (namem(k).EQ.name(icm)) THEN
68  im = k
69  END IF
70  4 CONTINUE
71  IF (im.EQ.0) go to 20
72 * Copy masses and obtain evolution parameters for first component.
73  m1 = cm(3,im)*zmbar
74  m2 = cm(4,im)*zmbar
75  kw = kstar(j1)
76  age = max(tplot,tev0(j1))*tstar - epoch(j1)
77  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
78  CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
79  & rm,lum,kw,mc,rcc,me,re,k2)
80  END IF
81  rj = r(jpair)
82  hj = h(jpair)
83 * Determine merger & ghost index for negative c.m. name.
84  IF (name(n+jpair).LT.0) THEN
85  CALL findj(j1,j,imerge)
86 * Skip second binary of quadruple.
87  IF (name(j).GT.nzero) go to 20
88  m1 = cm(1,imerge)*zmbar
89  m2 = cm(2,imerge)*zmbar
90  hj = hm(imerge)
91  rj = sqrt(xrel(1,imerge)**2 + xrel(2,imerge)**2 +
92  & xrel(3,imerge)**2)
93 * Re-define index of second component and obtain parameters of M1.
94  j2 = j
95  age = max(tplot,tev0(j1))*tstar - epoch(j1)
96  kw = kstar(j1)
97  CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
98  CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
99  & rm,lum,kw,mc,rcc,me,re,k2)
100  END IF
101  m0 = body0(j2)*zmbar
102  kw2 = kstar(j2)
103  age = max(tplot,tev0(j2))*tstar - epoch(j2)
104  CALL star(kw2,m0,m2,tm,tn,tscls,lums,gb,zpars)
105  CALL hrdiag(m0,age,m2,tm,tn,tscls,lums,gb,zpars,
106  & rm2,lum2,kw2,mc,rcc,me,re,k2)
107  ri = sqrt(ri)/rc
108 * Specify relevant binary mass.
109  IF (body(j1).GT.0.0d0) THEN
110  bodyi = (m1 + m2)/zmbar
111  ELSE
112  bodyi = cm(3,imerge) + cm(4,imerge)
113  END IF
114  semi = -0.5*bodyi/hj
115  ecc2 = (1.0 - rj/semi)**2
116  ecc = sqrt(ecc2)
117  pb = days*semi*sqrt(abs(semi)/bodyi)
118  pb = min(pb,99999.9d0)
119  pb = log10(abs(pb))
120  semi = log10(abs(semi*su))
121  r1 = log10(rm)
122  r2 = log10(rm2)
123  zl1 = log10(lum)
124  zl2 = log10(lum2)
125  te1 = 0.25*(zl1 - 2.0*r1) + 3.7
126  te2 = 0.25*(zl2 - 2.0*r2) + 3.7
127  WRITE (82,5) name(j1), name(j2), kw, kw2, kstar(icm),
128  & ri, ecc, pb, semi, m1, m2, zl1, zl2, r1, r2, te1, te2
129  5 FORMAT (2i6,2i3,i4,f6.1,f6.3,10f7.3)
130  nb = nb + 1
131  ELSE
132 * Create output file for single stars (skip chain subsystem or ghost).
133  IF (name(i).EQ.0.OR.body(i).EQ.0.0d0) go to 20
134  ri = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
135  & (x(3,i) - rdens(3))**2
136  ri = sqrt(ri)/rc
137  ri = min(ri,99.0d0)
138  r1 = log10(rm)
139  zl1 = log10(lum)
140 * Form LOG(Te) using L = 4*pi*R**2*\sigma*T**4 and solar value 3.7.
141  te = 0.25*(zl1 - 2.0*r1) + 3.7
142  WRITE (83,10) name(i), kw, ri, m1, zl1, r1, te
143  10 FORMAT (i10,i4,5f10.3)
144  nstar = nstar + 1
145  END IF
146  20 CONTINUE
147  WRITE (82,30) nb
148  WRITE (83,30) nstar
149  30 FORMAT (' ## END',i8)
150 *
151 * Update next plotting time.
152  tplot = tplot + dtplot
153  CALL flush(82)
154  CALL flush(83)
155 *
156  RETURN
157 *
158  END