7 IMPLICIT DOUBLE PRECISION(
d)
10 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
12 common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
14 dimension dps(5),kfl(3),pmq(3),
px(3),
py(3),gam(3),ie(2),pr(2),
15 &
in(9),dhm(4),dhg(4),dp(5,5),irank(2),mju(4),iju(3),pju(5,5),
16 &tju(5),kfjh(2),njs(2),kfjs(2),pjs(4,5)
19 four(i,j)=
p(i,4)*
p(j,4)-
p(i,1)*
p(j,1)-
p(i,2)*
p(j,2)-
p(i,3)*
p(j,3)
20 dfour(i,j)=dp(i,4)*dp(j,4)-dp(i,1)*dp(j,1)-dp(i,2)*dp(j,2)-
34 IF(i.GT.
min(
n,mstu(4)-mstu(32)))
THEN
35 CALL
luerrm(12,
'(LUSTRF:) failed to reconstruct jet system')
36 IF(mstu(21).GE.1)
RETURN
38 IF(
k(i,1).NE.1.AND.
k(i,1).NE.2.AND.
k(i,1).NE.41) goto 110
41 kq=kchg(kc,2)*isign(1,
k(i,2))
43 IF(
n+5*
np+11.GT.mstu(4)-mstu(32)-5)
THEN
44 CALL
luerrm(11,
'(LUSTRF:) no more memory left in LUJETS')
45 IF(mstu(21).GE.1)
RETURN
53 120 dps(j)=dps(j)+
p(i,j)
58 dps(4)=dps(4)+
max(0.,
p(
n+
np,4)-
p(i,4))
60 IF(kq.NE.2) kqsum=kqsum+kq
63 IF(kqsum.EQ.kq) mju(1)=
n+
np
64 IF(kqsum.NE.kq) mju(2)=
n+
np
66 IF(
k(i,1).EQ.2.OR.
k(i,1).EQ.41) goto 110
68 CALL
luerrm(12,
'(LUSTRF:) unphysical flavour combination')
69 IF(mstu(21).GE.1)
RETURN
73 CALL ludbrb(
n+1,
n+
np,0.,0.,-dps(1)/dps(4),-dps(2)/dps(4),
86 IF(i.EQ.
n+nr.AND.iabs(
k(
n+1,2)).NE.21) goto 140
89 IF(
k(i,1).EQ.41.OR.
k(i1,1).EQ.41) goto 140
90 IF(mju(1).NE.0.AND.i1.LT.mju(1).AND.iabs(
k(i1,2)).NE.21)
92 IF(mju(2).NE.0.AND.i.GT.mju(2).AND.iabs(
k(i,2)).NE.21) goto 140
93 pap=
sqrt((
p(i,1)**2+
p(i,2)**2+
p(i,3)**2)*(
p(i1,1)**2+
94 &
p(i1,2)**2+
p(i1,3)**2))
95 pvp=
p(i,1)*
p(i1,1)+
p(i,2)*
p(i1,2)+
p(i,3)*
p(i1,3)
96 pdr=4.*(pap-pvp)**2/(paru13**2*pap+2.*(pap-pvp))
97 IF(
pdr.LT.pdrmin)
THEN
104 IF(pdrmin.LT.paru12.AND.ir.EQ.
n+nr)
THEN
106 150
p(
n+1,j)=
p(
n+1,j)+
p(
n+nr,j)
111 ELSEIF(pdrmin.LT.paru12)
THEN
113 160
p(ir,j)=
p(ir,j)+
p(ir+1,j)
114 p(ir,5)=
sqrt(
max(0.,
p(ir,4)**2-
p(ir,1)**2-
p(ir,2)**2-
120 IF(ir.EQ.
n+nr-1)
k(ir,2)=
k(
n+nr,2)
122 IF(mju(1).GT.ir) mju(1)=mju(1)-1
123 IF(mju(2).GT.ir) mju(2)=mju(2)-1
134 IF(ntry.GT.100.AND.ntryr.LE.4)
THEN
138 ELSEIF(ntry.GT.100)
THEN
139 CALL
luerrm(14,
'(LUSTRF:) caught in infinite loop')
140 IF(mstu(21).GE.1)
RETURN
143 IF(mju(1).EQ.0.AND.mju(2).EQ.0) goto 500
146 IF(mju(jt).EQ.0) goto 490
155 DO 200 i1=
n+1+(jt-1)*(nr-1),
n+nr+(jt-1)*(1-nr),js
156 IF(
k(i1,2).NE.21.AND.iu.LE.2)
THEN
161 200 pju(iu,j)=pju(iu,j)+
p(i1,j)
163 210 pju(iu,5)=
sqrt(pju(iu,1)**2+pju(iu,2)**2+pju(iu,3)**2)
164 IF(
k(iju(3),2)/100.NE.10*
k(iju(1),2)+
k(iju(2),2).AND.
165 &
k(iju(3),2)/100.NE.10*
k(iju(2),2)+
k(iju(1),2))
THEN
166 CALL
luerrm(12,
'(LUSTRF:) unphysical flavour combination')
167 IF(mstu(21).GE.1)
RETURN
171 t12=(pju(1,1)*pju(2,1)+pju(1,2)*pju(2,2)+pju(1,3)*pju(2,3))/
173 t13=(pju(1,1)*pju(3,1)+pju(1,2)*pju(3,2)+pju(1,3)*pju(3,3))/
175 t23=(pju(2,1)*pju(3,1)+pju(2,2)*pju(3,2)+pju(2,3)*pju(3,3))/
177 t11=
sqrt((2./3.)*(1.-t12)*(1.-t13)/(1.-t23))
178 t22=
sqrt((2./3.)*(1.-t12)*(1.-t23)/(1.-t13))
179 tsq=
sqrt((2.*t11*t22+t12-1.)*(1.+t12))
180 t1f=(tsq-t22*(1.+t12))/(1.-t12**2)
181 t2f=(tsq-t11*(1.+t12))/(1.-t12**2)
183 220 tju(j)=-(t1f*pju(1,j)/pju(1,5)+t2f*pju(2,j)/pju(2,5))
184 tju(4)=
sqrt(1.+tju(1)**2+tju(2)**2+tju(3)**2)
186 230 pju(iu,5)=tju(4)*pju(iu,4)-tju(1)*pju(iu,1)-tju(2)*pju(iu,2)-
190 IF(pju(1,5)+pju(2,5).GT.pju(1,4)+pju(2,4))
THEN
210 IF(
is.EQ.1) dp(1,j)=
p(is1,j)
212 250
IF(
is.EQ.
ns) dp(2,j)=-pju(iu,j)
213 IF(
is.EQ.
ns) dp(2,4)=
sqrt(pju(iu,1)**2+pju(iu,2)**2+pju(iu,3)**2)
214 IF(
is.EQ.
ns) dp(2,5)=0.
218 IF(dp(3,5)+2.*dhkc+dp(4,5).LE.0.)
THEN
219 dp(1,4)=
sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
220 dp(2,4)=
sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
225 dhks=
sqrt(dhkc**2-dp(3,5)*dp(4,5))
226 dhk1=0.5*((dp(4,5)+dhkc)/dhks-1.)
227 dhk2=0.5*((dp(3,5)+dhkc)/dhks-1.)
229 p(in1,5)=
sqrt(dp(3,5)+2.*dhkc+dp(4,5))
231 p(in1,j)=(1.+dhk1)*dp(1,j)-dhk2*dp(2,j)
232 260
p(in1+1,j)=(1.+dhk2)*dp(2,j)-dhk1*dp(1,j)
237 IF(ntry.GT.100.AND.ntryr.LE.4)
THEN
241 ELSEIF(ntry.GT.100)
THEN
242 CALL
luerrm(14,
'(LUSTRF:) caught in infinite loop')
243 IF(mstu(21).GE.1)
RETURN
247 ie(1)=
k(
n+1+(jt/2)*(
np-1),3)
252 DO 280 in1=
n+nr+2+jq,
n+nr+4*
ns-2+jq,4
269 dp(1,4)=
sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
270 dp(2,4)=
sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
271 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
272 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
273 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
274 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
275 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
276 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
277 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
279 dhcx1=dfour(3,1)/dhc12
280 dhcx2=dfour(3,2)/dhc12
281 dhcxx=1d0/
sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
282 dhcy1=dfour(4,1)/dhc12
283 dhcy2=dfour(4,2)/dhc12
284 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
285 dhcyy=1d0/
sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
287 dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
289 310
p(
in(6)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
294 IF(2*i-nsav.GE.mstu(4)-mstu(32)-5)
THEN
295 CALL
luerrm(11,
'(LUSTRF:) no more memory left in LUJETS')
296 IF(mstu(21).GE.1)
RETURN
305 330 CALL
lukfdi(kfl(1),0,kfl(3),
k(i,2))
306 IF(
k(i,2).EQ.0) goto 270
307 IF(mstj(12).GE.3.AND.irankj.EQ.1.AND.iabs(kfl(1)).LE.10.AND.
308 &iabs(kfl(3)).GT.10)
THEN
309 IF(
rlu(0).GT.parj(19)) goto 330
313 pr(1)=
p(i,5)**2+(
px(1)+
px(3))**2+(
py(1)+
py(3))**2
314 CALL
luzdis(kfl(1),kfl(3),pr(1),
z)
315 gam(3)=(1.-
z)*(gam(1)+pr(1)/
z)
320 IF(
in(1)+1.EQ.
in(2).AND.
z*
p(
in(1)+2,3)*
p(
in(2)+2,3)*
321 &
p(
in(1),5)**2.GE.pr(1))
THEN
323 p(
in(2)+2,4)=pr(1)/(
p(
in(1)+2,4)*
p(
in(1),5)**2)
327 ELSEIF(
in(1)+1.EQ.
in(2))
THEN
331 IF(
in(2).GT.
n+nr+4*
ns) goto 270
340 360
IF(
in(1).GT.
n+nr+4*
ns.OR.
in(2).GT.
n+nr+4*
ns.OR.
341 &
in(1).GT.
in(2)) goto 270
342 IF(
in(1).NE.
in(4).OR.
in(2).NE.
in(5))
THEN
348 dp(1,4)=
sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
349 dp(2,4)=
sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
351 IF(dhc12.LE.1
e-2)
THEN
358 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
359 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
360 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
361 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
362 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
363 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
364 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
365 dhcx1=dfour(3,1)/dhc12
366 dhcx2=dfour(3,2)/dhc12
367 dhcxx=1d0/
sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
368 dhcy1=dfour(4,1)/dhc12
369 dhcy2=dfour(4,2)/dhc12
370 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
371 dhcyy=1d0/
sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
373 dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
375 380
p(
in(3)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
380 IF(abs(pxp**2+
pyp**2-
px(3)**2-
py(3)**2).LT.0.01)
THEN
391 DO 390 in1=
in(4),
in(1)-4,4
392 390
p(i,j)=
p(i,j)+
p(in1+2,3)*
p(in1,j)
393 DO 400 in2=
in(5),
in(2)-4,4
394 400
p(i,j)=
p(i,j)+
p(in2+2,3)*
p(in2,j)
401 DO 410 in2=
in(1)+1,
in(2),4
402 DO 410 in1=
in(1),in2-1,4
404 dhg(1)=dhg(1)+
p(in1+2,1)*
p(in2+2,1)*dhc
405 IF(in1.EQ.
in(1)) dhg(2)=dhg(2)-
p(in2+2,1)*dhc
406 IF(in2.EQ.
in(2)) dhg(3)=dhg(3)+
p(in1+2,1)*dhc
407 410
IF(in1.EQ.
in(1).AND.in2.EQ.
in(2)) dhg(4)=dhg(4)-dhc
410 dhs1=dhm(3)*dhg(4)-dhm(4)*dhg(3)
411 IF(abs(dhs1).LT.1
e-4) goto 270
412 dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(2)*dhg(3)-dhg(4)*
413 &(
p(i,5)**2-dhm(1))+dhg(2)*dhm(3)
414 dhs3=dhm(2)*(gam(3)-dhg(1))-dhg(2)*(
p(i,5)**2-dhm(1))
415 p(
in(2)+2,4)=0.5*(
sqrt(
max(0d0,dhs2**2-4.*dhs1*dhs3))/abs(dhs1)-
417 IF(dhm(2)+dhm(4)*
p(
in(2)+2,4).LE.0.) goto 270
418 p(
in(1)+2,4)=(
p(i,5)**2-dhm(1)-dhm(3)*
p(
in(2)+2,4))/
419 &(dhm(2)+dhm(4)*
p(
in(2)+2,4))
422 IF(
p(
in(2)+2,4).GT.
p(
in(2)+2,3))
THEN
426 IF(
in(2).GT.
n+nr+4*
ns) goto 270
433 ELSEIF(
p(
in(1)+2,4).GT.
p(
in(1)+2,3))
THEN
443 430 pju(iu+3,j)=pju(iu+3,j)+
p(i,j)
444 IF(
p(i,4).LE.0.) goto 270
445 pju(iu+3,5)=tju(4)*pju(iu+3,4)-tju(1)*pju(iu+3,1)-
446 &tju(2)*pju(iu+3,2)-tju(3)*pju(iu+3,3)
447 IF(pju(iu+3,5).LT.pju(iu,5))
THEN
452 IF(
in(3).NE.
in(6))
THEN
455 440
p(
in(6)+1,j)=
p(
in(3)+1,j)
459 p(
in(jq)+2,3)=
p(
in(jq)+2,3)-
p(
in(jq)+2,4)
460 450
p(
in(jq)+2,1)=
p(
in(jq)+2,1)-(3-2*jq)*
p(
in(jq)+2,4)
465 IF(iabs(kfl(1)).GT.10) goto 270
469 460 pju(iu+3,j)=pju(iu+3,j)-
p(i+1,j)
474 kfjs(jt)=
k(
k(mju(jt+2),3),2)
475 kfls=2*
int(
rlu(0)+3.*parj(4)/(1.+3.*parj(4)))+1
476 IF(kfjh(1).EQ.kfjh(2)) kfls=3
477 IF(ista.NE.i) kfjs(jt)=isign(1000*
max(iabs(kfjh(1)),
478 &iabs(kfjh(2)))+100*
min(iabs(kfjh(1)),iabs(kfjh(2)))+
481 pjs(jt,j)=pju(1,j)+pju(2,j)+
p(mju(jt),j)
482 480 pjs(jt+2,j)=pju(4,j)+pju(5,j)
483 pjs(jt,5)=
sqrt(
max(0.,pjs(jt,4)**2-pjs(jt,1)**2-pjs(jt,2)**2-
488 500
IF(mju(1).NE.0.AND.mju(2).NE.0)
THEN
491 ELSEIF(mju(1).NE.0)
THEN
494 ELSEIF(mju(2).NE.0)
THEN
497 ELSEIF(iabs(
k(
n+1,2)).NE.21)
THEN
505 510 w2sum=w2sum+
p(
n+nr+
is,1)
509 w2sum=w2sum-
p(
n+nr+nb,1)
510 IF(w2sum.GT.w2ran.AND.nb.LT.nr) goto 520
515 is1=
n+
is+nb-1-nr*((
is+nb-2)/nr)
516 is2=
n+
is+nb-nr*((
is+nb-1)/nr)
519 IF(iabs(
k(is1,2)).EQ.21) dp(1,j)=0.5*dp(1,j)
520 IF(is1.EQ.mju(1)) dp(1,j)=pjs(1,j)-pjs(3,j)
522 IF(iabs(
k(is2,2)).EQ.21) dp(2,j)=0.5*dp(2,j)
523 530
IF(is2.EQ.mju(2)) dp(2,j)=pjs(2,j)-pjs(4,j)
527 IF(dp(3,5)+2.*dhkc+dp(4,5).LE.0.)
THEN
530 dp(1,4)=
sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2+dp(1,5)**2)
531 dp(2,4)=
sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2+dp(2,5)**2)
534 dhks=
sqrt(dhkc**2-dp(3,5)*dp(4,5))
535 dhk1=0.5*((dp(4,5)+dhkc)/dhks-1.)
536 dhk2=0.5*((dp(3,5)+dhkc)/dhks-1.)
538 p(in1,5)=
sqrt(dp(3,5)+2.*dhkc+dp(4,5))
540 p(in1,j)=(1.+dhk1)*dp(1,j)-dhk2*dp(2,j)
541 540
p(in1+1,j)=(1.+dhk2)*dp(2,j)-dhk1*dp(1,j)
546 IF(ntry.GT.100.AND.ntryr.LE.4)
THEN
550 ELSEIF(ntry.GT.100)
THEN
551 CALL
luerrm(14,
'(LUSTRF:) caught in infinite loop')
552 IF(mstu(21).GE.1)
RETURN
561 IF(mju(jt).NE.0) irank(jt)=njs(jt)
562 IF(
ns.GT.nr) irank(jt)=1
563 ie(jt)=
k(
n+1+(jt/2)*(
np-1),3)
564 in(3*jt+1)=
n+nr+1+4*(jt/2)*(
ns-1)
565 in(3*jt+2)=
in(3*jt+1)+1
566 in(3*jt+3)=
n+nr+4*
ns+2*jt-1
567 DO 570 in1=
n+nr+2+jt,
n+nr+4*
ns-2+jt,4
576 IF(
ns.EQ.1.AND.mju(1)+mju(2).EQ.0) CALL
luptdi(0,
px(1),
py(1))
581 IF(mju(jt).NE.0) kfl(jt)=kfjs(jt)
588 kfl(3)=
int(1.+(2.+parj(2))*
rlu(0))*(-1)**
int(
rlu(0)+0.5)
589 CALL
lukfdi(kfl(3),0,kfl(1),kdump)
591 IF(iabs(kfl(1)).GT.10.AND.
rlu(0).GT.0.5)
THEN
592 kfl(2)=-(kfl(1)+isign(10000,kfl(1)))
593 ELSEIF(iabs(kfl(1)).GT.10)
THEN
594 kfl(1)=-(kfl(2)+isign(10000,kfl(2)))
599 pr3=
min(25.,0.1*
p(
n+nr+1,5)**2)
600 590 CALL
luzdis(kfl(1),kfl(2),pr3,
z)
601 zr=pr3/(
z*
p(
n+nr+1,5)**2)
602 IF(zr.GE.1.) goto 590
607 in1=
n+nr+3+4*(jt/2)*(
ns-1)
610 p(in1,3)=(2-jt)*(1.-
z)+(jt-1)*
z
613 600
p(in1+1,3)=(2-jt)*(1.-zr)+(jt-1)*zr
618 IF(jt.EQ.1.OR.
ns.EQ.nr-1)
THEN
626 dp(1,4)=
sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
627 dp(2,4)=
sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
628 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
629 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
630 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
631 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
632 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
633 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
634 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
636 dhcx1=dfour(3,1)/dhc12
637 dhcx2=dfour(3,2)/dhc12
638 dhcxx=1d0/
sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
639 dhcy1=dfour(4,1)/dhc12
640 dhcy2=dfour(4,2)/dhc12
641 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
642 dhcyy=1d0/
sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
644 dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
646 620
p(in3+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
651 630
p(in3+3,j)=
p(in3+1,j)
656 IF(mju(1)+mju(2).GT.0)
THEN
658 IF(njs(jt).EQ.0) goto 660
660 650
p(
n+nrs,j)=
p(
n+nrs,j)-pjs(jt+2,j)
666 IF(2*i-nsav.GE.mstu(4)-mstu(32)-5)
THEN
667 CALL
luerrm(11,
'(LUSTRF:) no more memory left in LUJETS')
668 IF(mstu(21).GE.1)
RETURN
671 IF(iabs(kfl(3-jt)).GT.10) jt=3-jt
674 irank(jt)=irank(jt)+1
681 680 CALL
lukfdi(kfl(jt),0,kfl(3),
k(i,2))
682 IF(
k(i,2).EQ.0) goto 550
683 IF(mstj(12).GE.3.AND.irank(jt).EQ.1.AND.iabs(kfl(jt)).LE.10.AND.
684 &iabs(kfl(3)).GT.10)
THEN
685 IF(
rlu(0).GT.parj(19)) goto 680
689 pr(jt)=
p(i,5)**2+(
px(jt)+
px(3))**2+(
py(jt)+
py(3))**2
694 wmin=parj(32+mstj(11))+pmq(1)+pmq(2)+parj(36)*pmq(3)
695 IF(iabs(kfl(jt)).GT.10.AND.iabs(kfl(3)).GT.10) wmin=
696 &wmin-0.5*parj(36)*pmq(3)
698 IF(wrem2.LT.0.10) goto 550
699 IF(wrem2.LT.
max(wmin*(1.+(2.*
rlu(0)-1.)*parj(37)),
700 &parj(32)+pmq(1)+pmq(2))**2) goto 810
703 CALL
luzdis(kfl(jt),kfl(3),pr(jt),
z)
707 &
mod(kfl2a/1000,10)).GE.4)
THEN
708 pr(jr)=(pmq(jr)+pmq(3))**2+(
px(jr)-
px(3))**2+(
py(jr)-
py(3))**2
709 pw12=
sqrt(
max(0.,(wrem2-pr(1)-pr(2))**2-4.*pr(1)*pr(2)))
710 z=(wrem2+pr(jt)-pr(jr)+pw12*(2.*
z-1.))/(2.*wrem2)
711 pr(jr)=(pmq(jr)+parj(32+mstj(11)))**2+(
px(jr)-
px(3))**2+
713 IF((1.-
z)*(wrem2-pr(jt)/
z).LT.pr(jr)) goto 810
715 gam(3)=(1.-
z)*(gam(jt)+pr(jt)/
z)
720 IF(
in(1)+1.EQ.
in(2).AND.
z*
p(
in(1)+2,3)*
p(
in(2)+2,3)*
721 &
p(
in(1),5)**2.GE.pr(jt))
THEN
723 p(
in(jr)+2,4)=pr(jt)/(
p(
in(jt)+2,4)*
p(
in(1),5)**2)
727 ELSEIF(
in(1)+1.EQ.
in(2))
THEN
728 p(
in(jr)+2,4)=
p(
in(jr)+2,3)
731 IF(js*
in(jr).GT.js*
in(4*jr)) goto 550
733 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
740 710
IF(js*
in(1).GT.js*
in(3*jr+1).OR.js*
in(2).GT.js*
in(3*jr+2).OR.
741 &
in(1).GT.
in(2)) goto 550
742 IF(
in(1).NE.
in(3*jt+1).OR.
in(2).NE.
in(3*jt+2))
THEN
748 dp(1,4)=
sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
749 dp(2,4)=
sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
751 IF(dhc12.LE.1
e-2)
THEN
752 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
758 dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
759 dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
760 dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
761 IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
762 IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
763 IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
764 IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
765 dhcx1=dfour(3,1)/dhc12
766 dhcx2=dfour(3,2)/dhc12
767 dhcxx=1d0/
sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
768 dhcy1=dfour(4,1)/dhc12
769 dhcy2=dfour(4,2)/dhc12
770 dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
771 dhcyy=1d0/
sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
773 dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
775 730
p(
in(3)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
782 IF(abs(pxp**2+
pyp**2-
px(3)**2-
py(3)**2).LT.0.01)
THEN
791 p(i,j)=
px(jt)*
p(
in(3*jt+3),j)+
py(jt)*
p(
in(3*jt+3)+1,j)+
793 DO 740 in1=
in(3*jt+1),
in(1)-4*js,4*js
794 740
p(i,j)=
p(i,j)+
p(in1+2,3)*
p(in1,j)
795 DO 750 in2=
in(3*jt+2),
in(2)-4*js,4*js
796 750
p(i,j)=
p(i,j)+
p(in2+2,3)*
p(in2,j)
803 DO 760 in2=
in(1)+1,
in(2),4
804 DO 760 in1=
in(1),in2-1,4
806 dhg(1)=dhg(1)+
p(in1+2,jt)*
p(in2+2,jt)*dhc
807 IF(in1.EQ.
in(1)) dhg(2)=dhg(2)-js*
p(in2+2,jt)*dhc
808 IF(in2.EQ.
in(2)) dhg(3)=dhg(3)+js*
p(in1+2,jt)*dhc
809 760
IF(in1.EQ.
in(1).AND.in2.EQ.
in(2)) dhg(4)=dhg(4)-dhc
812 dhs1=dhm(jr+1)*dhg(4)-dhm(4)*dhg(jr+1)
813 IF(abs(dhs1).LT.1
e-4) goto 550
814 dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(jt+1)*dhg(jr+1)-dhg(4)*
815 &(
p(i,5)**2-dhm(1))+dhg(jt+1)*dhm(jr+1)
816 dhs3=dhm(jt+1)*(gam(3)-dhg(1))-dhg(jt+1)*(
p(i,5)**2-dhm(1))
817 p(
in(jr)+2,4)=0.5*(
sqrt(
max(0d0,dhs2**2-4.*dhs1*dhs3))/abs(dhs1)-
819 IF(dhm(jt+1)+dhm(4)*
p(
in(jr)+2,4).LE.0.) goto 550
820 p(
in(jt)+2,4)=(
p(i,5)**2-dhm(1)-dhm(jr+1)*
p(
in(jr)+2,4))/
821 &(dhm(jt+1)+dhm(4)*
p(
in(jr)+2,4))
824 IF(
p(
in(jr)+2,4).GT.
p(
in(jr)+2,3))
THEN
825 p(
in(jr)+2,4)=
p(
in(jr)+2,3)
828 IF(js*
in(jr).GT.js*
in(4*jr)) goto 550
830 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
835 ELSEIF(
p(
in(jt)+2,4).GT.
p(
in(jt)+2,3))
THEN
836 p(
in(jt)+2,4)=
p(
in(jt)+2,3)
845 780
p(
n+nrs,j)=
p(
n+nrs,j)-
p(i,j)
846 IF(
p(i,4).LE.0.) goto 550
852 IF(
in(3).NE.
in(3*jt+3))
THEN
854 p(
in(3*jt+3),j)=
p(
in(3),j)
855 790
p(
in(3*jt+3)+1,j)=
p(
in(3)+1,j)
859 p(
in(jq)+2,3)=
p(
in(jq)+2,3)-
p(
in(jq)+2,4)
860 800
p(
in(jq)+2,jt)=
p(
in(jq)+2,jt)-js*(3-2*jq)*
p(
in(jq)+2,4)
869 CALL
lukfdi(kfl(jr),-kfl(3),kfldmp,
k(i,2))
870 IF(
k(i,2).EQ.0) goto 550
872 pr(jr)=
p(i,5)**2+(
px(jr)-
px(3))**2+(
py(jr)-
py(3))**2
879 dhr1=
four(
n+nrs,
in(3*jq+2))/dhc12
880 dhr2=
four(
n+nrs,
in(3*jq+1))/dhc12
881 IF(
in(4).NE.
in(7).OR.
in(5).NE.
in(8))
THEN
884 pr(3-jq)=
p(i+(jt+jq-3)**2-1,5)**2+(
px(3-jq)+(2*jq-3)*js*
885 &
px(3))**2+(
py(3-jq)+(2*jq-3)*js*
py(3))**2
889 wrem2=wrem2+(
px(1)+
px(2))**2+(
py(1)+
py(2))**2
891 IF(mju(1)+mju(2).NE.0.AND.i.EQ.isav+2.AND.
fd.GE.1.) goto 180
892 IF(
fd.GE.1.) goto 550
893 fa=wrem2+pr(jt)-pr(jr)
894 IF(mstj(11).EQ.2) prev=0.5*
fd**parj(37+mstj(11))
896 &parj(37+mstj(11))*(pr(1)+pr(2))**2))
902 &4.*wrem2*pr(jt))),float(js))
905 &
p(
in(3*jq+3)+1,j)+0.5*(dhr1*(
fa+
fb)*
p(
in(3*jq+1),j)+
906 &dhr2*(
fa-
fb)*
p(
in(3*jq+2),j))/wrem2
907 820
p(i,j)=
p(
n+nrs,j)-
p(i-1,j)
911 DO 830 i=nsav+1,nsav+
np
914 IF(mstu(16).NE.2)
THEN
932 840
v(nsav,j)=
v(ip,j)
933 p(nsav,5)=
sqrt(
max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))
945 860
p(i-nsav+
n,j)=
p(i,j)
947 DO 880 i=
n+1,2*
n-nsav
948 IF(
k(i,3).NE.ie(1)) goto 880
953 IF(mstu(16).NE.2)
k(i1,3)=nsav
955 DO 900 i=2*
n-nsav,
n+1,-1
956 IF(
k(i,3).EQ.ie(1)) goto 900
961 IF(mstu(16).NE.2)
k(i1,3)=nsav
965 CALL ludbrb(nsav+1,
n,0.,0.,dps(1)/dps(4),dps(2)/dps(4),