274
275
276
277#include "implicit_f.inc"
278#include "comlock.inc"
279
280
281
282 INTEGER NRTMFT, NRTMLT , IRECT(4,*), NRADM
283
285 . x(3,*), nod_normal(3,*), rcurv(*), anglm(*), anglt
286
287
288
289 INTEGER I ,, N2, N3, N4
291 . x1, x2, x3, x4,
292 . y1, y2, y3, y4,
293 . z1, z2, z3, z4,
294 . nnx1, nnx2, nnx3, nnx4,
295 . nny1, nny2, nny3, nny4,
296 . nnz1, nnz2, nnz3, nnz4,
297 . surfx, surfy, surfz,
298 . erx, ery, erz, dnx, dny, dnz, dnt, ll, aaa, rr,
299 . x13, y13, z13, x24, y24, z24, nx, ny, nz, cc
300
301 rcurv(nrtmft:nrtmlt) = ep30
302 anglm(nrtmft:nrtmlt) = ep30
303
304 DO i=nrtmft, nrtmlt
305 n1=irect(1,i)
306 n2=irect(2,i)
307 n3=irect(3,i)
308 n4=irect(4,i)
309
310 x1=x(1,n1)
311 y1=x(2,n1)
312 z1=x(3,n1)
313
314 x2=x(1,n2)
315 y2=x(2,n2)
316 z2=x(3,n2)
317
318 x3=x(1,n3)
319 y3=x(2,n3)
320 z3=x(3,n3)
321
322 x4=x(1,n4)
323 y4=x(2,n4)
324 z4=x(3,n4)
325
326 nnx1=nod_normal(1,n1)
327 nny1=nod_normal(2,n1)
328 nnz1=nod_normal(3,n1)
329
330 nnx2=nod_normal(1,n2)
331 nny2=nod_normal(2,n2)
332 nnz2=nod_normal(3,n2)
333
334 nnx3=nod_normal(1,n3)
335 nny3=nod_normal(2,n3)
336 nnz3=nod_normal(3,n3)
337
338 nnx4=nod_normal(1,n4)
339 nny4=nod_normal(2,n4)
340 nnz4=nod_normal(3,n4)
341
342
343 erx = (x2+x3)-(x1+x4)
344 ery = (y2+y3)-(y1+y4)
345 erz = (z2+z3)-(z1+z4)
346
347
348 ll = sqrt(erx*erx+ery*ery+erz*erz)
349 aaa = one / ll
350 erx = erx*aaa
351 ery = ery*aaa
352 erz = erz*aaa
353
354 dnx= (nnx2+nnx3)-(nnx1+nnx4)
355 dny= (nny2+nny3)-(nny1+nny4)
356 dnz= (nnz2+nnz3)-(nnz1+nnz4)
357
358
359 dnt=(dnx*erx+dny*ery+dnz*erz)
360
361 rr=ll/
max(em20,abs(dnt))
362 rcurv(i)=
min(rcurv(i),rr)
363
364 erx = (x4+x3)-(x1+x2)
365 ery = (y4+y3)-(y1+y2)
366 erz = (z4+z3)-(z1+z2)
367
368
369 ll = sqrt(erx*erx+ery*ery+erz*erz)
370 aaa = one / ll
371 erx = erx*aaa
372 ery = ery*aaa
373 erz = erz*aaa
374
375 dnx= (nnx4+nnx3)-(nnx1+nnx2)
376 dny= (nny4+nny3)-(nny1+nny2)
377 dnz= (nnz4+nnz3)-(nnz1+nnz2)
378
379
380 dnt=(dnx*erx+dny*ery+dnz*erz)
381
382 rr=ll/(nradm*
max(em20,abs(dnt)))
383 rcurv(i)=
min(rcurv(i),rr)
384
385
386
387
388
389 x13 = x3 - x1
390 y13 = y3 - y1
391 z13 = z3 - z1
392
393 x24 = x4 - x2
394 y24 = y4 - y2
395 z24 = z4 - z2
396
397 surfx = y13*z24 - z13*y24
398 surfy = z13*x24 - x13*z24
399 surfz = x13*y24 - y13*x24
400
401 aaa=one/
max(em30,sqrt(surfx*surfx+surfy*surfy+surfz*surfz))
402 surfx = surfx * aaa
403 surfy = surfy * aaa
404 surfz = surfz * aaa
405
406 cc=(surfx*nnx1+surfy*nny1+surfz*nnz1)/
max(em20,anglt)
407 anglm(i)=
min(anglm(i),cc)
408
409 cc=(surfx*nnx2+surfy*nny2+surfz*nnz2)/
max(em20,anglt)
410 anglm(i)=
min(anglm(i),cc)
411
412 cc=(surfx*nnx3+surfy*nny3+surfz*nnz3)/
max(em20,anglt)
413 anglm(i)=
min(anglm(i),cc)
414
415 cc=(surfx*nnx4+surfy*nny4+surfz*nnz4)/
max(em20,anglt)
416 anglm(i)=
min(anglm(i),cc)
417 ENDDO
418
419 RETURN