221
222
223
224#include "implicit_f.inc"
225
226
227
228 INTEGER NSN, NMN, NIR, I0, I2SIZE, IDEL2
229 INTEGER IRECT(4,*), MSR(*), NSV(*), IRTL(*), WEIGHT(*),
230 . IADI2(NIR,*),IRUPT(*)
231
233 . a(3,*), ar(3,*),crst(2,*), ms(*),
234 . x(3,*),in(*),stifn(*),stifr(*), fskyi2(i2size,*), csts_bis(2,*)
235
236
237
238 INTEGER I, J, II, L, NN
239
241 . ss, tt, xmsi, fxi, fyi, fzi, mxi, myi, mzi,ins,
242 . x0,x1,x2,x3,x4,y0,y1,y2,y3,y4,z0,z1,z2,z3,z4,aa,
243 . xc0,yc0,zc0,sp,sm,tp,tm,xc,yc,zc,h1, h2,h3,h4,
244 . stf,h21,h22,h23,h24
245
246#include "vectorize.inc"
247 DO ii=1,nmn
248 j=msr(ii)
249 in(j)=
max(em20,in(j))
250 ENDDO
251
252 DO ii=1,nsn
253 i=nsv(ii)
254 IF (i > 0) THEN
255
256 IF (irupt(ii) == 0) THEN
257
258
259 IF (weight(i) == 1) THEN
260 l=irtl(ii)
261
262 ss=crst(1,ii)
263 tt=crst(2,ii)
264 sp=one + ss
265 sm=one - ss
266 tp=fourth*(one + tt)
267 tm=fourth*(one - tt)
268 h1=tm*sm
269 h2=tm*sp
270 h3=tp*sp
271 h4=tp*sm
272
273
274 ss=csts_bis(1,ii)
275 tt=csts_bis(2,ii)
276 sp=one + ss
277 sm=one - ss
278 tp=fourth*(one + tt)
279 tm=fourth*(one - tt)
280 h21=tm*sm
281 h22=tm*sp
282 h23=tp*sp
283 h24=tp*sm
284
285 x0 = x(1,i)
286 y0 = x(2,i)
287 z0 = x(3,i)
288
289 x1 = x(1,irect(1,l))
290 y1 = x(2,irect(1,l))
291 z1 = x(3,irect(1,l))
292 x2 = x(1,irect(2,l))
293 y2 = x(2,irect(2,l))
294 z2 = x(3,irect(2,l))
295 x3 = x(1,irect(3,l))
296 y3 = x(2,irect(3,l))
297 z3 = x(3,irect(3,l))
298 x4 = x(1,irect(4,l))
299 y4 = x(2,irect(4,l))
300 z4 = x(3,irect(4,l))
301
302 xc = x1 * h1 + x2 * h2 + x3 * h3 + x4 * h4
303 yc = y1 * h1 + y2 * h2 + y3 * h3 + y4 * h4
304 zc = z1 * h1 + z2 * h2 + z3 * h3 + z4 * h4
305
306 xc0=x0-xc
307 yc0=y0-yc
308 zc0=z0-zc
309
310 aa = xc0*xc0 + yc0*yc0 + zc0*zc0
311 ins = in(i) + aa * ms(i)
312 stf = stifr(i) + aa * stifn(i)
313
314 fxi=a(1,i)
315 fyi=a(2,i)
316 fzi=a(3,i)
317
318 mxi = ar(1,i) + yc0 * fzi - zc0 * fyi
319 myi = ar(2,i) + zc0 * fxi - xc0 * fzi
320 mzi = ar(3,i) + xc0 * fyi - yc0 * fxi
321
322 i0 = i0 + 1
323 nn = iadi2(1,i0)
324 fskyi2(6,nn) = mxi*h1
325 fskyi2(7,nn) = myi*h1
326 fskyi2(8,nn) = mzi*h1
327 fskyi2(9,nn) = ins*h21
328 fskyi2(10,nn)= stf*h1
329
330 nn = iadi2(2,i0)
331 fskyi2(6,nn) = mxi*h2
332 fskyi2(7,nn) = myi*h2
333 fskyi2(8,nn) = mzi*h2
334 fskyi2(9,nn) = ins*h22
335 fskyi2(10,nn)= stf*h2
336
337 nn = iadi2(3,i0)
338 fskyi2(6,nn) = mxi*h3
339 fskyi2(7,nn) = myi*h3
340 fskyi2(8,nn) = mzi*h3
341 fskyi2(9,nn) = ins*h23
342 fskyi2(10,nn)= stf*h3
343
344 nn = iadi2(4,i0)
345 fskyi2(6,nn) = mxi*h4
346 fskyi2(7,nn) = myi*h4
347 fskyi2(8,nn) = mzi*h4
348 fskyi2(9,nn) = ins*h24
349 fskyi2(10,nn)= stf*h4
350 ENDIF
351
352 stifn(i)=em20
353 stifr(i)=em20
354 in(i) =zero
355 ms(i) =zero
356 a(1,i)=zero
357 a(2,i)=zero
358 a(3,i)=zero
359
360 ELSEIF (weight(i) == 1) THEN
361
362
363 i0 = i0 + 1
364 nn = iadi2(1,i0)
365 fskyi2(6,nn) = zero
366 fskyi2(7,nn) = zero
367 fskyi2(8,nn) = zero
368 fskyi2(9,nn) = zero
369 fskyi2(10,nn)= zero
370
371 nn = iadi2(2,i0)
372 fskyi2(6,nn) = zero
373 fskyi2(7,nn) = zero
374 fskyi2(8,nn) = zero
375 fskyi2(9,nn) = zero
376 fskyi2(10,nn)= zero
377
378 nn = iadi2(3,i0)
379 fskyi2(6,nn) = zero
380 fskyi2(7,nn) = zero
381 fskyi2(8,nn) = zero
382 fskyi2(9,nn) = zero
383 fskyi2(10,nn)= zero
384
385 nn = iadi2(4,i0)
386 fskyi2(6,nn) = zero
387 fskyi2(7,nn) = zero
388 fskyi2(8,nn) = zero
389 fskyi2(9,nn) = zero
390 fskyi2(10,nn)= zero
391 ENDIF
392
393 ENDIF
394 ENDDO
395
396 RETURN