31
32
33
34#include "implicit_f.inc"
35
36
37
38#include "param_c.inc"
39#include "com08_c.inc"
40
41
42
43 INTEGER NN, NC
44 INTEGER LLL(*),IADLL(*),NPBYL(NNPBY,*),LPBYL(*),COMNTAG(*)
45
47 . rbyl(nrby,*),x(3,*),v(3,*),vr(3,*),a(3,*),ar(3,*),
48 . mass(*),iner(*)
49
50
51
52 INTEGER I, J, N, NS, MSL, MSL2, IFX, IFR, NFIX, NFRE
53 my_real xx,yy,zz,xy,yz,xz,ixx,iyy,izz,ixy,ixz,iyz,
54 . jxx,jyy,jzz,jxy,jxz,jyz,jxy2,jxz2,jyz2,det,
55 . b1,b2,b3,c1,c2,c3,vx1,vx2,vx3,mmas,usdt,ddt,
56 . vg(3),ag(3),xm(3),vtm(3),vrm(3),atm(3),arm(3)
57
58 msl = npbyl(2,nn)
59 msl2= msl*2
60 nfix= 0
61 nfre= 0
62 ifx = 0
63 ifr = 0
64
65 ns = lpbyl(msl)
66 IF (comntag(ns)==1) THEN
67 nfix = nfix + 1
68 lpbyl(msl+nfix) = ns
69 ELSE
70 nfre = nfre + 1
71 lpbyl(msl2+nfre) = ns
72 ENDIF
73
74 DO n=1,msl-1
75 ns = lpbyl(n)
76 IF (comntag(ns)==1) THEN
77 nfix = nfix + 1
78 lpbyl(msl+nfix) = ns
79 ELSE
80 nfre = nfre + 1
81 lpbyl(msl2+nfre) = ns
82 ENDIF
83 ENDDO
84 IF (nfix<=1) nfix = 0
85
86 IF (nfix/=0) THEN
87 ifx = lpbyl(msl+1)
88 jxx = zero
89 jyy = zero
90 jzz = zero
91 jxy = zero
92 jyz = zero
93 jxz = zero
94 mmas = zero
95 DO j=1,3
96 xm(j) = zero
97 vg(j) = zero
98 ag(j) = zero
99 vtm(j) = zero
100 atm(j) = zero
101 vrm(j) = zero
102 arm(j) = zero
103 ENDDO
104
105 DO i=1,nfix
106 n = lpbyl(msl+i)
107 mmas = mmas + mass(n)
108 DO j=1,3
109 xm(j)= xm(j) + x(j,n)*mass(n)
110 vtm(j)=vtm(j) + v(j,n)*mass(n)
111 atm(j)=atm(j) + a(j,n)*mass(n)
112 ENDDO
113 ENDDO
114 DO j=1,3
115 xm(j) = xm(j) / mmas
116 vtm(j) = vtm(j) / mmas
117 atm(j) = atm(j) / mmas
118 ENDDO
119 DO i=1,nfix
120 n = lpbyl(msl+i)
121 xx = x(1,n) - xm(1)
122 yy = x(2,n) - xm(2)
123 zz = x(3,n) - xm(3)
124
125
126
127
128
129
130 ag(1)=ag(1) + ar(1,n)*iner(n)+mass(n)*(yy*a(3,n)-zz*a(2,n))
131 ag(2)=ag(2) + ar(2,n)*iner(n)+mass(n)*(zz*a(1,n)-xx*a(3,n))
132 ag(3)=ag(3) + ar(3,n)*iner(n)+mass(n)*(xx*a(2,n)-yy*a(1,n))
133
134 xy=xx*yy
135 yz=yy*zz
136 xz=xx*zz
137
138 xx=xx*xx
139 yy=yy*yy
140 zz=zz*zz
141 ixx = iner(n)+(yy+zz)*mass(n)
142 iyy = iner(n)+(xx+zz)*mass(n)
143 izz = iner(n)+(xx+yy)*mass(n)
144 ixy = xy*mass(n)
145 iyz = yz*mass(n)
146 ixz = xz*mass(n)
147 jxx = jxx + ixx
148 jyy = jyy + iyy
149 jzz = jzz + izz
150 jxy = jxy - ixy
151 jyz = jyz - iyz
152 jxz = jxz - ixz
153 ENDDO
154
155 jxy2 = jxy*jxy
156 jyz2 = jyz*jyz
157 jxz2 = jxz*jxz
158 det = jxx*jyy*jzz-jxx*jyz2-jyy*jxz2-jzz*jxy2-two*jxy*jyz*jxz
159 det = one / det
160 b1 = det*(jzz*jyy-jyz2)
161 b2 = det*(jxx*jzz-jxz2)
162 b3 = det*(jyy*jxx-jxy2)
163 c1 = det*(jxx*jyz+jxz*jxy)
164 c2 = det*(jyy*jxz+jxy*jyz)
165 c3 = det*(jzz*jxy+jyz*jxz)
166
167 vrm(1) = vr(1,ifx)
168 vrm(2) = vr(2,ifx)
169 vrm(3) = vr(3,ifx)
170
171 vg(1) = vrm(1)*jxx + vrm(2)*jxy + vrm(3)*jxz
172 vg(2) = vrm(1)*jxy + vrm(2)*jyy + vrm(3)*jyz
173 vg(3) = vrm(1)*jxz + vrm(2)*jyz + vrm(3)*jzz
174
175
176
177
178
179
180 ag(1) = ag(1) - vrm(2)*vg(3) + vrm(3)*vg(2)
181 ag(2) = ag(2) - vrm(3)*vg(1) + vrm(1)*vg(3)
182 ag(3) = ag(3) - vrm(1)*vg(2) + vrm(2)*vg(1)
183
184 arm(1)= ag(1)*b1 + ag(2)*c3 + ag(3)*c2
185 arm(2)= ag(1)*c3 + ag(2)*b2 + ag(3)*c1
186 arm(3)= ag(1)*c2 + ag(2)*c1 + ag(3)*b3
187
188 IF (nfre==0) THEN
189
190 usdt = one / dt12
191 ddt = half * dt12
192 DO j=1,3
193 vrm(j) = vrm(j) + arm(j)*dt12
194 ENDDO
195 DO n=1,msl
196 ns = lpbyl(n)
197 DO j=1,3
198 ar(j,ns) = (vrm(j)-vr(j,ns)) * usdt
199 ENDDO
200 xx = x(1,ns)-xm(1)
201 yy = x(2,ns)-xm(2)
202 zz = x(3,ns)-xm(3)
203 vx1 = vrm(2)*zz - vrm(3)*yy
204 vx2 = vrm(3)*xx - vrm(1)*zz
205 vx3 = vrm(1)*yy - vrm(2)*xx
206 a(1,ns) = atm(1) + usdt*
207 . (vtm(1)-v(1,ns)+vx1+ddt*(vrm(2)*vx3-vrm(3)*vx2))
208 a(2,ns) = atm(2) + usdt*
209 . (vtm(2)-v(2,ns)+vx2+ddt*(vrm(3)*vx1-vrm(1)*vx3))
210 a(3,ns) = atm(3) + usdt*
211 . (vtm(3)-v(3,ns)+vx3+ddt*(vrm(1)*vx2-vrm(2)*vx1))
212
213
214
215
216 ENDDO
217 ELSE
218
219 ifr = lpbyl(2*msl+1)
220 rbyl(10,nn) = mass(ifx)
221 rbyl(14,nn) = v(1,ifx)
222 rbyl(15,nn) = v(2,ifx)
223 rbyl(16,nn) = v(3,ifx)
224 rbyl(17,nn) = vr(1,ifx)
225 rbyl(18,nn) = vr(2,ifx)
226 rbyl(19,nn) = vr(3,ifx)
227 rbyl(20,nn) = a(1,ifx)
228 rbyl(21,nn) = a(2,ifx)
229 rbyl(22,nn) = a(3,ifx)
230 rbyl(23,nn) = ar(1,ifx)
231 rbyl(24,nn) = ar(2,ifx)
232 rbyl(25,nn) = ar(3,ifx)
233
234 rbyl(1,nn) = b1
235 rbyl(2,nn) = b2
236 rbyl(3,nn) = b3
237 rbyl(4,nn) = c1
238 rbyl(5,nn) = c2
239 rbyl(6,nn) = c3
240 rbyl(11,nn) = xm(1)
241 rbyl(12,nn) = xm(2)
242 rbyl(13,nn) = xm(3)
243
244 mass(ifx) = mmas
245 v(1,ifx) = vtm(1)
246 v(2,ifx) = vtm(2)
247 v(3,ifx) = vtm(3)
248 vr(1,ifx) = vrm(1)
249 vr(2,ifx) = vrm(2)
250 vr(3,ifx) = vrm(3)
251 a(1,ifx) = atm(1)
252 a(2,ifx) = atm(2)
253 a(3,ifx) = atm(3)
254 ar(1,ifx) = arm(1)
255 ar(2,ifx) = arm(2)
256 ar(3,ifx) = arm(3)
257 ENDIF
258 ENDIF
259 npbyl(4,nn) = nfix
260 npbyl(5,nn) = nfre
261 npbyl(7,nn) = ifx
262 npbyl(8,nn) = ifr
263
264 RETURN