40
41
42
46
47
48
49#include "implicit_f.inc"
50
51
52
53#include "mvsiz_p.inc"
54
55
56
57#include "units_c.inc"
58#include "param_c.inc"
59#include "vect01_c.inc"
60#include "scr03_c.inc"
61#include "scr11_c.inc"
62#include "com04_c.inc"
63
64
65
66 INTEGER MAT(*),NIX,IX(NIX,*)
67 INTEGER,INTENT(IN) :: M151_ID
68 my_real pm(npropm,nummat), tb(*),x(3,*),vdet2
69 TYPE(DETONATORS_STRUCT_)::DETONATORS
70
71
72
73 INTEGER I, N, , MT, IOPT
74 INTEGER NDETPS,NDETSG,NECRAN,NDETPL,NDETCORD,NPE,NPE2,NDETPS_NO_SHADOW,NDETPS_SHADOW
75 LOGICAL HAS_DETONATOR(MVSIZ)
77 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
78 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
79 . xc(mvsiz), yc(mvsiz), zc(mvsiz), bt(mvsiz),
80 . dl(mvsiz), alt, xlp, ylp, zlp, xlp1,
81 . ylp1, zlp1, xlp2, ylp2, zlp2, xl0, yl0, zl0, xl1, yl1, zl1,
82 . xl2, yl2, zl2, ps1, ps2, dl1, dl2, s1, s2, s3,
83 . nx, ny, nz , nn, vdet
84 INTEGER :: NODE1, NODE2, NODE3, NODE4, II, GRNOD_ID, INOD, NNOD, NODE_ID
85 INTEGER :: I_SHADOW_FLAG
86
87
88 ndetps = detonators%N_DET_POINT
89 ndetsg = detonators%N_DET_LINE
90 necran = detonators%N_DET_WAVE_SHAPER
91 ndetpl = detonators%N_DET_PLANE
92 ndetcord = detonators%N_DET_CORD
93
94
95 ndetps_no_shadow = 0
96 ndetps_shadow = 0
97 DO i = 1, detonators%N_DET_POINT
98 i_shadow_flag = detonators%POINT(i)%SHADOW
99 IF(i_shadow_flag == 0)THEN
100 ndetps_no_shadow = ndetps_no_shadow + 1
101 ELSE
102 ndetps_shadow = ndetps_shadow + 1
103 ENDIF
104 ENDDO
105
106
107
108 IF(detonators%N_DET - ndetps_shadow > 0) THEN
109
110 DO i = lft, llt
111 ii = i + nft
112 node1 = ix(2, ii)
113 node2 = ix(3, ii)
114 node3 = ix(4, ii)
115 node4 = ix(5, ii)
116 y1(i) = x(2, node1)
117 y2(i) = x(2, node2)
118 y3(i) = x(2, node3)
119 y4(i) = x(2, node4)
120 z1(i) = x(3, node1)
121 z2(i) = x(3, node2)
122 z3(i) = x(3, node3)
123 z4(i) = x(3, node4)
124 ENDDO
125
126
127
128
129 DO i=lft,llt
130 tb(i)=ep21
131 has_detonator(i)=.false.
132 xc(i)=zero
133 yc(i)=fourth*(y1(i)+y2(i)+y3(i)+y4(i))
134 zc(i)=fourth*(z1(i)+z2(i)+z3(i)+z4(i))
135 END DO
136
137
138
139
140 IF(ndetps /= 0) THEN
141 DO i=lft,llt
142 DO n=1,ndetps
143 i_shadow_flag = detonators%POINT(n
144 IF(i_shadow_flag /= 0)cycle
145 mtl=detonators%POINT(n)%MAT
146 grnod_id=detonators%POINT(n)%GRNOD_ID
147 IF(mtl == 0 .OR. mtl == mat(i) .OR. mtl == m151_id) THEN
148
149 IF(grnod_id == 0)THEN
150 alt=detonators%POINT(n)%TDET
151 xlp=detonators%POINT(n)%XDET
152 ylp=detonators%POINT(n)%YDET
153 zlp=detonators%POINT(n)%ZDET
154 dl(i) =(xc(i)-xlp)**2+(yc(i)-ylp)**2+(zc(i)-zlp)**2
155 dl(i)=sqrt(dl(i))
156 bt(i) =alt+dl(i)/pm(38,mat(i))
157 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
158 has_detonator(i)=.true.
159
160 ELSE
161 nnod = detonators%POINT(n)%NNOD
162 alt=detonators%POINT(n)%TDET
163 has_detonator(i)=.true.
164 DO inod=1,nnod
165 node_id=detonators%POINT(n)%NODLIST(inod)
166 xlp=x(1,node_id)
167 ylp=x(2,node_id)
168 zlp=x(3,node_id)
169 dl(i) =(xc(i)-xlp)**2+(yc(i)-ylp)**2+(zc(i)-zlp)**2
170 dl(i)=sqrt(dl(i))
171 bt(i) =alt+dl(i)/pm(38,mat(i))
172 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
173 ENDDO
174 ENDIF
175 ENDIF
176 END DO
177 END DO
178 ENDIF
179
180
181
182
183 IF(ndetsg /= 0) THEN
184 DO n=1,ndetsg
185 alt=detonators%LINE(n)%TDET
186 mtl=detonators%LINE(n)%MAT
187 xlp1=detonators%LINE(n)%XDET_1
188 ylp1=detonators%LINE(n)%YDET_1
189 zlp1=detonators%LINE(n)%ZDET_1
190 xlp2=detonators%LINE(n)%XDET_2
191 ylp2=detonators%LINE(n)%YDET_2
192 zlp2=detonators%LINE(n)%ZDET_2
193 DO i=lft,llt
194 IF(mtl == 0 .OR. mtl == mat(i) .OR. mtl == m151_id) THEN
195 xl0 =(xlp1-xlp2)
196 yl0 =(ylp1-ylp2)
197 zl0 =(zlp1-zlp2)
198 xl1 =(xc(i)-xlp1)
199 yl1 =(yc(i)-ylp1)
200 zl1 =(zc(i)-zlp1)
201 xl2 =(xc(i)-xlp2)
202 yl2 =(yc(i)-ylp2)
203 zl2 =(zc(i)-zlp2)
204 ps1 =xl1*xl0+yl1*yl0+zl1*zl0
205 ps2 =xl2*xl0+yl2*yl0+zl2*zl0
206 IF(ps1*ps2 > zero)THEN
207 dl1 =sqrt(xl1**2+yl1**2+zl1**2)
208 dl2 =sqrt(xl2**2+yl2**2+zl2**2)
210 ELSE
211 s1 =yl1*zl0 - zl1*yl0
212 s2 =zl1*xl0 - xl1*zl0
213 s3 =xl1*yl0 - yl1*xl0
214 dl(i)=sqrt((s1**2+s2**2+s3**2)/(xl0**2+yl0**2+zl0**2))
215 ENDIF
216 bt(i)=alt+dl(i)/pm(38,mat(i))
217 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
218 has_detonator(i)= .true.
219 END IF
220 END DO
221 END DO
222 ENDIF
223
224
225
226
227 IF(necranTHEN
228 DO n=1,necran
229 alt=detonators%WAVE_SHAPER(n)%TDET
230 mtl=detonators%WAVE_SHAPER(n)%MAT
231 vdet =detonators%WAVE_SHAPER(n)%VDET
232 yd =detonators%WAVE_SHAPER(n)%YDET
233 zd =detonators%WAVE_SHAPER(n)%ZDET
234 npe=detonators%WAVE_SHAPER(n)%NUMNOD
235 dto0=alt
236 vdto=pm(38,mat(1))
237 IF(vdet == zero)vdet=pm(38,mat(1))
238
239 CALL ecran1(detonators%WAVE_SHAPER(n),x,vdet)
240
241 DO i=lft,llt
242 IF(mtl /= mat(i) .AND. mtl /= 0 .AND. mtl /= m151_id) cycle
243 yd =detonators%WAVE_SHAPER(n)%YDET
244 zd =detonators%WAVE_SHAPER(n)%ZDET
245 dto0=alt
246 yl=yc(i)
247 zl=zc(i)
248 CALL ecran2(detonators%WAVE_SHAPER(n),x,vdet)
249 bt(i) =dto
250 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
251 has_detonator(i)= .true.
252 END DO
253 END DO
254 ENDIF
255
256
257
258
259 IF(ndetpl /= 0) THEN
260 DO n=1,ndetpl
261 alt=detonators%PLANE(n)%TDET
262 mtl=detonators%PLANE(n)%MAT
263 xlp=detonators%PLANE(n)%XDET
264 ylp=detonators%PLANE(n)%YDET
265 zlp=detonators%PLANE(n)%ZDET
266 nx=detonators%PLANE(n)%NX
267 ny=detonators%PLANE(n)%NY
268 nz=detonators%PLANE(n)%NZ
269 nn=sqrt(nx**2+ny**2+nz**2)
271 dl1=xlp*nx + ylp*ny + zlp*nz
272 dl1 = dl1/nn
273 DO i=lft,llt
274 IF(mtl == 0 .OR. mtl == mat(i) .OR. mtl == m151_id) THEN
275
276
277
278
279 dl(i) = (xc(i)-xlp)*nx + (yc(i)-ylp)*ny + (zc(i)-zlp)*nz
280 dl(i) = abs(dl(i))
281 dl(i) = dl(i)/nn
282 bt(i) =alt+dl(i)/pm(38,mat(i))
283 IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
284 has_detonator(i)=.true.
285 END IF
286 END DO
287 END DO
288 ENDIF
289
290
291
292
293 IF(ndetcord /= 0) THEN
294 DO n=1,ndetcord
295 alt = detonators%CORD(n)%TDET
296 mtl = detonators%CORD(n)%MAT
297 vdet2 = detonators%CORD(n)%VDET
298 iopt = detonators%CORD(n)%IOPT
299 npe2 = detonators%CORD(n)%NUMNOD
300 dto0 = alt
301 mt = mtl
302 IF(mt == 0)mt=mat(1)
303 vdto = pm(38,mt)
304 IF(mtl /= mat(1) .AND. mtl /= 0 .AND. mtl /= m151_id) cycle
305 dto0 = alt
306 CALL detcord(detonators%CORD(n),x,xc,yc,zc,vdto,vdet2,alt,bt,tb,has_detonator,iopt)
307 END do
308 ENDIF
309
310
311 END IF
312
313
314 RETURN
subroutine detcord(detonator_cord, x, xc, yc, zc, vdet, vdet2, alt, bt, tb, has_detonator, iopt)
subroutine ecran1(detonator, x, vdet)
subroutine ecran2(detonator_wave_shaper, x, vdet)