Nbody6
lagr.f
Go to the documentation of this file.
1  SUBROUTINE lagr(C)
2 *
3 *
5 * -----------------
6 *
7  include 'common6.h'
8  REAL*8 r2,mrt,rt
9  common/work1/ r2(nmax)
10  parameter(lx=11)
11  REAL*8 c(3),flagr(lx),rlagr(lx),rm(lx),dens(lx),vr(lx),avm(lx)
12 * DATA FLAGR/-1.9,-1.7,-1.5,-1.3,-1.1,-.9,-.7,-.5,-.3,-.1/
13 * DATA FLAGR/0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,
14 * & 0.75,0.9/
15  DATA flagr/0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.625,0.75,0.9/
16 *
17 *
18 * Set square radii of single particles & c.m. bodies.
19  np = 0
20  DO 10 i = ifirst,ntot
21  np = np + 1
22  r2(np) = (x(1,i) - c(1))**2 + (x(2,i) - c(2))**2 +
23  & (x(3,i) - c(3))**2
24  jlist(np) = i
25  10 CONTINUE
26 *
27 * Sort square distances with respect to the centre C.
28  CALL sort1(np,r2,jlist)
29
30 * Determine Lagrangian radii for specified mass fractions.
31  zm = 0.0
32  DO 20 il = 1,lx
33  zm1 = zm
34  zm = 0.0
35  zmh = flagr(il)*zmass
36  vm2 = 0.0
37  IF (il.GT.1) THEN
38  r1 = sqrt(r2(i))
39  iprev = i
40  ELSE
41  r1 = 0.0
42  iprev = 0
43  END IF
44  i = 0
45  15 i = i + 1
46  im = jlist(i)
47  zm = zm + body(im)
48 * Sum mass-weighted square velocity in current shell.
49  IF (zm.GT.zm1) THEN
50  vm2 = vm2 + body(im)*(xdot(1,im)**2 + xdot(2,im)**2 +
51  & xdot(3,im)**2)
52  END IF
53  IF (zm.LT.zmh) go to 15
54  rlagr(il) = sqrt(r2(i))
55  IF (kz(7).LT.5) go to 20
56 * Form mean square velocity and density (print outer shell radii).
57  dm = zm - zm1
58  IF (abs(dm).LT.1.0d-10) THEN
59  vr(il) = 0.0
60  ELSE
61  vr(il) = sqrt(vm2/dm)
62  END IF
63  dv = 2.0*twopi/3.0*(r2(i)**1.5 - r1**3)
64  dens(il) = dm/dv
65  IF (dens(il).LE.0.0d0) dens(il) = 1.0
66  rm(il) = 0.5*(r1 + rlagr(il))
67 * Obtain the average mass in solar units.
68  IF (i.GT.iprev) THEN
69  avm(il) = dm/float(i - iprev)*smu
70  ELSE
71  avm(il) = 0.0
72  END IF
73  20 CONTINUE
74 *
75 * Obtain half-mass radius separately.
76  zm = 0.0
77  zmh = 0.5*zmass
78  i = 0
79  25 i = i + 1
80  im = jlist(i)
81  zm = zm + body(im)
82  IF (zm.LT.zmh) go to 25
83 *
84 * Replace approximate half-mass radius by actual value.
85  rscale = sqrt(r2(i))
86 *
87 * Determine tidal radius and mass if kz(14) = 2.
88  if (kz(14).eq.2.and.tidal(1).ne.0.) then
89  i = np
90  mrt = zmass
91  30 continue
92  if (mrt.gt.0.) then
93  rt = (mrt/tidal(1))**(1./3.)
94  else
95  rt = 0.0
96  endif
97  IF (sqrt(r2(i)).gt.rt.and.rt.gt.0.and.i.gt.1) then
98  im = jlist(i)
99  mrt = mrt - body(im)
100  i = i - 1
101  go to 30
102  endif
103  endif
104
105  IF ((kz(7).EQ.2.OR.kz(7).EQ.4).AND.time.GE.tnext) THEN
106  IF (kz(14).EQ.2) WRITE (6,*) tphys,mrt,rt
107  END IF
108 *
109  IF (kz(7).GE.3.AND.time.GE.tnext) THEN
110  IF (kz(14).EQ.2) WRITE (14,35) tphys, mrt, rt
111  35 FORMAT (3x,'TIDAL RADIUS TPHYS MRT RT ',f8.1,1p,2e10.2)
112  END IF
113 *
114 * Check output options (line printer or unit 14 or both).
115  IF (kz(7).EQ.2.OR.kz(7).EQ.4.AND.time.GE.tnext) THEN
116  WRITE (6,40) (log10(rlagr(k)),k=1,lx)
117  40 FORMAT (/,' LAGR: ',13f7.3)
118  END IF
119 *
120  IF (kz(7).GE.3.AND.time.GE.tnext) THEN
121  WRITE (14,50) ttot, (log10(rlagr(k)),k=1,lx)
122  50 FORMAT (' LAGR: ',f7.1,13f7.3)
123  CALL flush(14)
124  END IF
125 *
126  IF (kz(7).EQ.5.AND.time.GE.tnext) THEN
127  WRITE (26,60) ttot, (dens(k),k=1,lx)
128  60 FORMAT (' DENSITY (T =',f7.1,'): ',1p,13e10.2)
129  WRITE (26,65) (rm(k),k=1,lx)
130  65 FORMAT (' DISTANCE: ',1p,13e10.2)
131  CALL flush(26)
132  WRITE (27,70) ttot, (vr(k),k=1,lx)
133  70 FORMAT (' VELOCITY (T =',f7.1,'): ',1p,13e10.2)
134  CALL flush(27)
135  WRITE (36,75) ttot, (avm(k),k=1,lx)
136  75 FORMAT (' AVERAGE MASS (T =',f7.1,'): ',13f7.3)
137  CALL flush(36)
138  END IF
139  IF (kz(7).EQ.6.AND.time.GE.tnext) THEN
140  WRITE (28,80) ttot, (avm(k), rm(k),k=1,lx)
141  80 FORMAT (' PAIRWISE <M> <RM> (T =',f7.1,'): ',
142  & 13(0p,f7.3,1p,e10.2))
143  CALL flush(28)
144  END IF
145 *
146  RETURN
147 *
148  END