224
225
226
227#include "implicit_f.inc"
228
229
230
231 INTEGER NIR
232
234 . rs(3),rx(4),ry(4),rz(4),fmx(4),fmy(4),fmz(4),h(4),stifm,mxs,mys,mzs,stifmr,
235 . betax,betay
236
237
238
239#include "units_c.inc"
240#include "param_c.inc"
241#include "com01_c.inc"
242
243
244
245 INTEGER J
246 my_real mx,my,mz,d1,d2,dd,rax,rbx,ray,rby,mxy,myx,myz,mzy,mxz,mzx,mxx,myy,mzz,rxx,ryy,fxx,fyy,fa,fb
247 my_real rsx(4),rsy(4),rsz(4),df(nir),mat(3,3),mati(3,3),bx,by,bz,
alpha
248
249
250
251 rsx(1) = rx(1) - rs(1)
252 rsy(1) = ry(1) - rs(2)
253 rsz(1) = rz(1) - rs(3)
254 rsx(2) = rx(2) - rs(1)
255 rsy(2) = ry(2) - rs(2)
256 rsz(2) = rz(2) - rs(3)
257 rsx(3) = rx(3) - rs(1)
258 rsy(3) = ry(3) - rs(2)
259 rsz(3) = rz(3) - rs(3)
260 rsx(4) = rx(4) - rs(1)
261 rsy(4) = ry(4) - rs(2)
262 rsz(4) = rz(4) - rs(3)
263
264 mx = zero
265 my = zero
266 mz = zero
267 DO j=1,nir
268 mx = mx + rsy(j)*fmz(j) - rsz(j)*fmy(j)
269 my = my + rsz(j)*fmx(j) - rsx(j)*fmz(j)
270 mz = mz + rsx(j)*fmy(j) - rsy(j)*fmx(j)
271 ENDDO
272
273 bx= rsx(1)*h(1)+rsx(2)*h(2)+rsx(3)*h(3)+rsx(4)*h(4)
274 by= rsy(1)*h(1)+rsy(2)*h(2)+rsy(3)*h(3)+rsy(4)*h(4)
275 bz= rsz(1)*h(1)+rsz(2)*h(2)+rsz(3)*h(3)+rsz(4)*h(4)
276
277
278 IF (nir == 4) THEN
279 rxx = rx(2)+rx(3)-rx(1)-rx(4)
280 ryy = ry(3)+ry(4)-ry(1)-ry(2)
281
282
283
284
285 alpha = abs(rxx)/(abs(rxx)+abs(ryy))
286
287
288
289 mxy = rsx(1)*fmy(1)+rsx(2)*fmy(2)+rsx(3)*fmy(3)+rsx(4)*fmy(4)-
alpha*mzs
290 myx = rsy(1)*fmx(1)+rsy(2)*fmx(2)+rsy(3)*fmx(3)+rsy(4)*fmx(4)+(one-
alpha)*mzs
291 mzz = mxy - myx
292
293
294
295
296
297
298 fyy = (betay*mxy-(one-betax)*myx)/rxx
299 fxx = ((one-betay)*mxy-betax*myx)/ryy
300
301 stifm=
max(sqrt((one-betax)*by*by+ betay*bx*bx)/abs(rxx),sqrt(betax*by*by+ (one-betay)*bx*bx)/abs(ryy))
302 stifmr=one/(abs(rxx)+abs(ryy))
303
304 fmx(1) = fmx(1) - fxx
305 fmx(2) = fmx(2) - fxx
306 fmx(3) = fmx(3) + fxx
307 fmx(4) = fmx(4) + fxx
308
309 fmy(1) = fmy(1) + fyy
310 fmy(2) = fmy(2) - fyy
311 fmy(3) = fmy(3) - fyy
312 fmy(4) = fmy(4) + fyy
313
314
315
316
317
318 myz = rsy(1)*fmz(1)+rsy(2)*fmz(2)+rsy(3)*fmz(3)+rsy(4)*fmz(4)
319 mzy = rsz(1)*fmy(1)+rsz(2)*fmy(2)+rsz(3)*fmy(3)+rsz(4)*fmy(4)
320 mxx = betax*(myz - mzy - mxs)
321 mzx = rsz(1)*fmx(1)+rsz(2)*fmx(2)+rsz(3)*fmx(3)+rsz(4)*fmx(4)
322 mxz = rsx(1)*fmz(1)+rsx(2)*fmz(2)+rsx(3)*fmz(3)+rsx(4)*fmz(4)
323 myy = betay*(mzx - mxz - mys)
324
325 rax = rx(1) - rx(3)
326 rbx = rx(4) - rx(2)
327 ray = ry(1) - ry(3)
328 rby = ry(4) - ry(2)
329 d1 = -mxx*rbx - myy*rby
330 d2 = mxx*rax + myy*ray
331 dd = ray*rbx - rax*rby
332 fa = d1 / dd
333 fb = d2 / dd
334
336 . sqrt((by*by+bz*bz)*rbx*rbx*betax+(bz*bz+bx*bx)*rby*rby*betay)/abs(dd),
337 . sqrt((by*by+bz*bz)*rax*rax*betax+(bz*bz+bx*bx)*ray*ray*betay)/abs(dd))
338
339 stifmr=
max(stifmr,(abs(rbx*betax)+abs(rby*betay))/abs(dd),(abs(rax*betax)+abs(ray*betay))/abs(dd))
340
341 fmz(1) = fmz(1) + fa
342 fmz(2) = fmz(2) - fb
343 fmz(3) = fmz(3) - fa
344 fmz(4) = fmz(4) + fb
345
346
347 ELSEIF (nir == 3) THEN
348 rxx = rx(2)+rx(3)-two*rx(1)
349 ryy = two*ry(3)-ry(1)-ry(2)
350
351
352
353 alpha = abs(rxx)/(abs(rxx)+abs(ryy))
354
355
356
357
358 mxy = rsx(1)*fmy(1)+rsx(2)*fmy(2)+rsx(3)*fmy(3)-
alpha*mzs
359 myx = rsy(1)*fmx(1)+rsy(2)*fmx(2)+rsy(3)*fmx(3)+(one-
alpha)*mzs
360 mzz = mxy - myx
361
362
363
364
365
366
367 fyy = -(betay*mxy-(one-betax)*myx)/rxx
368 fxx = ((one-betay)*mxy-betax*myx)/ryy
369
370 stifm=
max(sqrt((one-betax)*by*by+ betay*bx*bx)/abs(rxx),sqrt(betax*by*by+ (one-betay)*bx*bx)/abs(ryy))
371 stifmr=one/(abs(rxx)+abs(ryy))
372
373 fmx(1) = fmx(1) - fxx
374 fmx(2) = fmx(2) - fxx
375 fmx(3) = fmx(3) + two * fxx
376
377 fmy(1) = fmy(1) - two * fyy
378 fmy(2) = fmy(2) + fyy
379 fmy(3) = fmy(3) + fyy
380
381
382
383
384
385 mx=betax*(mx+(two*rsz(1)-rsz(2)-rsz(3))*fyy - mxs)
386 my=betay*(my+(two*rsz(3)-rsz(1)-rsz(2))*fxx - mys)
387
388 mat(1,1) = rsy(1)
389 mat(1,2) = rsy(2)
390 mat(1,3) = rsy(3)
391 mat(2,1) = rsx(1)
392 mat(2,2) = rsx(2)
393 mat(2,3) = rsx(3)
394 mat(3,1) = one
395 mat(3,2) = one
396 mat(3,3) = one
398 df(1) = -mati(1,1)*mx + mati(1,2)*my
399 df(2) = -mati(2,1)*mx + mati(2,2)*my
400 df(3) = -mati(3,1)*mx + mati(3,2)*my
401 DO j=1,nir
402 fmz(j) = fmz(j) + df(j)
403 ENDDO
404
405 bx= rsx(1)*h(1)+rsx(2)*h(2)+rsx(3)*h(3)
406 by= rsy(1)*h(1)+rsy(2)*h(2)+rsy(3)*h(3)
407 bz= rsz(1)*h(1)+rsz(2)*h(2)+rsz(3)*h(3)
408 stifm=
max(stifm,sqrt(bx*bx+by*by+bz*bz)*sqrt(
max(
409 . betax*(mati(1,1)*mati(1,1)+mati(2,1)*mati(2,1)+mati(3,1)*mati(3,1)),
410 . betay*(mati(1,2)*mati(1,2)+mati(2,2)*mati(2,2)+mati(3,2)*mati(3,2)),
411 . mati(1,3)*mati(1,3)+mati(2,3)*mati(2,3)+mati(3,3)*mati(3,3))))
412
413 stifmr=
max(stifmr,abs(mati(1,1))+abs(mati(1,2)),abs(mati(2,1))+abs(mati(2,2)),abs(mati(3,1))+abs(mati(3,2)))
414
415
416 ENDIF
417
418 RETURN