35
36
37
38#include "implicit_f.inc"
39
40
41
42 INTEGER INOD(4)
43
45 . s,t,hsl
47 . x(3,*),v(3,*)
49 . e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
50 . ls1,ls2,lt1,lt2
51
52
53
54#include "com08_c.inc"
55
56
57
58 INTEGER I,J,J1,J2,J3,J4,II,JJ,L,NN,NIR
59
61 . x21,x32,x34,x41,y21,y32,y34,y41,z21,z32,z34,z41,
62 . sm,sp,tm,tp,sm2,sp2,tm2,tp2,sm3,sp3,tm3,tp3,suma,s1,s2
64 . xm(4),ym(4),zm(4)
65
66 nir=4
67 IF (inod(3) == inod(4)) THEN
68 nir=3
69 xm(4)=zero
70 ym(4)=zero
71 zm(4)=zero
72 ENDIF
73 DO jj=1,nir
74 j=inod(jj)
75 xm(jj)=x(1,j)+ v(1,j)*dt12*half
76 ym(jj)=x(2,j)+ v(2,j)*dt12*half
77 zm(jj)=x(3,j)+ v(3,j)*dt12*half
78 ENDDO
79
80 x21 = xm(2) - xm(1)
81 x32 = xm(3) - xm(2)
82 x34 = xm(3) - xm(4)
83 x41 = xm(4) - xm(1)
84 y21 = ym(2) - ym(1)
85 y32 = ym(3) - ym(2)
86 y34 = ym(3) - ym(4)
87 y41 = ym(4) - ym(1)
88 z21 = zm(2) - zm(1)
89 z32 = zm(3) - zm(2)
90 z34 = zm(3) - zm(4)
91 z41 = zm(4) - zm(1)
92
93 e1x = (x21+x34 )
94 e1y = (y21+y34 )
95 e1z = (z21+z34 )
96
97 e2x = (x32+x41 )
98 e2y = (y32+y41 )
99 e2z = (z32+z41 )
100
101 e3x = e1y*e2z-e1z*e2y
102 e3y = e1z*e2x-e1x*e2z
103 e3z = e1x*e2y-e1y*e2x
104
105 suma = e3x*e3x+e3y*e3y+e3z*e3z
106
107 e3x = e3x * suma
108 e3y = e3y * suma
109 e3z = e3z * suma
110
111 s1 = e1x*e1x+e1y*e1y+e1z*e1z
112 s2 = e2x*e2x+e2y*e2y+e2z*e2z
113 suma = sqrt(s1/s2)
114 e1x = e1x + (e2y*e3z-e2z*e3y)*suma
115 e1y = e1y + (e2z*e3x-e2x*e3z)*suma
116 e1z = e1z + (e2x*e3y-e2y*e3x)*suma
117
118 suma = e1x*e1x+e1y*e1y+e1z*e1z
119 suma = one /
max(sqrt(suma),em20)
120 e1x = e1x * suma
121 e1y = e1y * suma
122 e1z = e1z * suma
123
124 e2x = e3y * e1z - e3z * e1y
125 e2y = e3z * e1x - e3x * e1z
126 e2z = e3x * e1y - e3y * e1x
127
128 ls1 = sqrt(x21**2 + y21**2 + z21**2)
129 ls2 = sqrt(x34**2 + y34**2 + z34**2)
130 lt1 = sqrt(x41**2 + y41**2 + z41**2)
131 lt2 = sqrt(x32**2 + y32**2 + z32**2)
132
133 RETURN