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