File indexing completed on 2024-04-06 12:13:52
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122 SUBROUTINE KTCLUS(IMODE,PP,NN,ECUT,Y,*)
0123 IMPLICIT NONE
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 INTEGER IMODE,NN
0139 DOUBLE PRECISION PP(4,*)
0140 DOUBLE PRECISION ECUT,Y(*),ONE
0141 ONE=1
0142 CALL KTCLUR(IMODE,PP,NN,ONE,ECUT,Y,*999)
0143 RETURN
0144 999 RETURN 1
0145 END
0146
0147 SUBROUTINE KTCLUR(IMODE,PP,NN,R,ECUT,Y,*)
0148 IMPLICIT NONE
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164 INTEGER NMAX,IM,IMODE,TYPE,ANGL,MONO,RECO,N,I,J,NN,
0165 & IMIN,JMIN,KMIN,NUM,HIST,INJET,IABBR,NABBR
0166 PARAMETER (NMAX=512,NABBR=7)
0167 DOUBLE PRECISION PP(4,*)
0168 DOUBLE PRECISION R,ECUT,Y(*),P,KT,ETOT,RSQ,KTP,KTS,KTPAIR,KTSING,
0169 & KTMIN,ETSQ,KTLAST,KTMAX,KTTMP
0170 LOGICAL FIRST
0171 CHARACTER TITLE(4,4)*10
0172
0173
0174
0175
0176
0177
0178
0179
0180 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0181 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0182 DIMENSION INJET(NMAX),IABBR(NABBR)
0183 DATA FIRST,TITLE,IABBR/.TRUE.,
0184 & 'e+e- ','ep ','pe ','pp ',
0185 & 'angle ','DeltaR ','f(DeltaR) ','**********',
0186 & 'no ','yes ','**********','**********',
0187 & 'E ','Pt ','Pt**2 ','**********',
0188 & 1111,2111,3111,4111,4211,4212,4223/
0189
0190 IM=IMODE
0191 IF (IM.GE.1.AND.IM.LE.NABBR) IM=IABBR(IM)
0192 TYPE=MOD(IM/1000,10)
0193 ANGL=MOD(IM/100 ,10)
0194 MONO=MOD(IM/10 ,10)
0195 RECO=MOD(IM ,10)
0196 IF (NN.GT.NMAX.OR.NN.LT.1.OR.(NN.LT.2.AND.TYPE.EQ.1))
0197 & CALL KTWARN('KTCLUS',100,*999)
0198 IF (TYPE.LT.1.OR.TYPE.GT.4.OR.ANGL.LT.1.OR.ANGL.GT.4.OR.
0199 & MONO.LT.1.OR.MONO.GT.2.OR.RECO.LT.1.OR.RECO.GT.3)
0200 & CALL KTWARN('KTCLUS',101,*999)
0201 IF (FIRST) THEN
0202 WRITE (6,'(/,1X,54(''*'')/A)')
0203 & ' KTCLUS: written by Mike Seymour, July 1992.'
0204 WRITE (6,'(A)')
0205 & ' Last modified November 2000.'
0206 WRITE (6,'(A)')
0207 & ' Please send comments or suggestions to Mike.Seymour@rl.ac.uk'
0208 WRITE (6,'(/A,I2,2A)')
0209 & ' Collision type =',TYPE,' = ',TITLE(TYPE,1)
0210 WRITE (6,'(A,I2,2A)')
0211 & ' Angular variable =',ANGL,' = ',TITLE(ANGL,2)
0212 WRITE (6,'(A,I2,2A)')
0213 & ' Monotonic definition =',MONO,' = ',TITLE(MONO,3)
0214 WRITE (6,'(A,I2,2A)')
0215 & ' Recombination scheme =',RECO,' = ',TITLE(RECO,4)
0216 IF (R.NE.1) THEN
0217 WRITE (6,'(A,F5.2)')
0218 & ' Radius parameter =',R
0219 IF (TYPE.NE.4) WRITE (6,'(A)')
0220 & ' R.NE.1 is strongly discouraged for this collision type!'
0221 ENDIF
0222 WRITE (6,'(1X,54(''*'')/)')
0223 FIRST=.FALSE.
0224 ENDIF
0225
0226 N=NN
0227 NUM=NN
0228 CALL KTCOPY(PP,N,P,(RECO.NE.1))
0229 ETOT=0
0230 DO 100 I=1,N
0231 ETOT=ETOT+P(4,I)
0232 100 CONTINUE
0233 IF (ETOT.EQ.0) CALL KTWARN('KTCLUS',102,*999)
0234 IF (ECUT.EQ.0) THEN
0235 ETSQ=1/ETOT**2
0236 ELSE
0237 ETSQ=1/ECUT**2
0238 ENDIF
0239 RSQ=R**2
0240
0241 DO 210 I=1,N-1
0242 DO 200 J=I+1,N
0243 KTP(J,I)=-1
0244 KTP(I,J)=KTPAIR(ANGL,P(1,I),P(1,J),KTP(J,I))
0245 200 CONTINUE
0246 210 CONTINUE
0247
0248 DO 230 I=1,N
0249 KTS(I)=KTSING(ANGL,TYPE,P(1,I))
0250 230 CONTINUE
0251 KTMAX=0
0252
0253 300 CONTINUE
0254
0255 CALL KTPMIN(KTP,NMAX,N,IMIN,JMIN)
0256
0257 CALL KTSMIN(KTS,NMAX,N,KMIN)
0258
0259 KTMIN=KTP(IMIN,JMIN)
0260 KTTMP=RSQ*KTS(KMIN)
0261 IF ((TYPE.GE.2.AND.TYPE.LE.4).AND.
0262 & (KTTMP.LE.KTMIN.OR.N.EQ.1))
0263 & KTMIN=KTTMP
0264 KT(N)=KTMIN
0265 Y(N)=KT(N)*ETSQ
0266
0267 IF (KTMIN.LT.KTMAX.AND.MONO.GT.1) CALL KTWARN('KTCLUS',1,*999)
0268 IF (KTMIN.GE.KTMAX) KTMAX=KTMIN
0269
0270 IF (KTMIN.EQ.KTTMP) THEN
0271 CALL KTMOVE(P,KTP,KTS,NMAX,N,KMIN,1)
0272
0273 HIST(N)=KMIN
0274 INJET(N)=KMIN
0275 DO 400 I=N,NN
0276 IF (INJET(I).EQ.KMIN) THEN
0277 KTLAST(I)=KTMAX
0278 INJET(I)=0
0279 ELSEIF (INJET(I).EQ.N) THEN
0280 INJET(I)=KMIN
0281 ENDIF
0282 400 CONTINUE
0283
0284 ELSE
0285 CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,TYPE,ANGL,MONO,RECO)
0286 CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,1)
0287
0288 HIST(N)=IMIN*NMAX+JMIN
0289 INJET(N)=IMIN
0290 DO 600 I=N,NN
0291 IF (INJET(I).EQ.JMIN) THEN
0292 INJET(I)=IMIN
0293 ELSEIF (INJET(I).EQ.N) THEN
0294 INJET(I)=JMIN
0295 ENDIF
0296 600 CONTINUE
0297 ENDIF
0298
0299 N=N-1
0300 IF (N.GT.1 .OR. N.GT.0.AND.(TYPE.GE.2.AND.TYPE.LE.4)) GOTO 300
0301 IF (N.EQ.1) THEN
0302 KT(N)=1D20
0303 Y(N)=KT(N)*ETSQ
0304 ENDIF
0305 RETURN
0306 999 RETURN 1
0307 END
0308
0309 SUBROUTINE KTYCUT(ECUT,NY,YCUT,NJET,*)
0310 IMPLICIT NONE
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323 INTEGER NY,NJET(NY),NMAX,HIST,I,J,NUM
0324 PARAMETER (NMAX=512)
0325 DOUBLE PRECISION YCUT(NY),ETOT,RSQ,P,KT,KTP,KTS,ETSQ,ECUT,KTLAST,
0326 & ROUND
0327 PARAMETER (ROUND=0.99999D0)
0328 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0329 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0330 IF (ETOT.EQ.0) CALL KTWARN('KTYCUT',100,*999)
0331 IF (ECUT.EQ.0) THEN
0332 ETSQ=1/ETOT**2
0333 ELSE
0334 ETSQ=1/ECUT**2
0335 ENDIF
0336 DO 100 I=1,NY
0337 NJET(I)=0
0338 100 CONTINUE
0339 DO 210 I=NUM,1,-1
0340 DO 200 J=1,NY
0341 IF (NJET(J).EQ.0.AND.KT(I)*ETSQ.GE.ROUND*YCUT(J)) NJET(J)=I
0342 200 CONTINUE
0343 210 CONTINUE
0344 RETURN
0345 999 RETURN 1
0346 END
0347
0348 SUBROUTINE KTYSUB(ECUT,NY,YCUT,YMAC,NSUB,*)
0349 IMPLICIT NONE
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366 INTEGER NY,NSUB(NY),NMAX,HIST,I,J,NUM
0367 PARAMETER (NMAX=512)
0368 DOUBLE PRECISION YCUT(NY),YMAC,ETOT,RSQ,P,KT,KTP,KTS,ETSQ,ECUT,
0369 & KTLAST,ROUND
0370 PARAMETER (ROUND=0.99999D0)
0371 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0372 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0373 IF (ETOT.EQ.0) CALL KTWARN('KTYSUB',100,*999)
0374 IF (ECUT.EQ.0) THEN
0375 ETSQ=1/ETOT**2
0376 ELSE
0377 ETSQ=1/ECUT**2
0378 ENDIF
0379 DO 100 I=1,NY
0380 NSUB(I)=0
0381 100 CONTINUE
0382 DO 210 I=NUM,1,-1
0383 DO 200 J=1,NY
0384 IF (NSUB(J).EQ.0.AND.KT(I)*ETSQ.GE.ROUND*YCUT(J)) NSUB(J)=I
0385 IF (NSUB(J).NE.0.AND.KTLAST(I)*ETSQ.LT.ROUND*YMAC)
0386 & NSUB(J)=NSUB(J)-1
0387 200 CONTINUE
0388 210 CONTINUE
0389 RETURN
0390 999 RETURN 1
0391 END
0392
0393 SUBROUTINE KTBEAM(ECUT,Y,*)
0394 IMPLICIT NONE
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405 INTEGER NMAX,HIST,NUM,I,J
0406 PARAMETER (NMAX=512)
0407 DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,ECUT,ETSQ,Y(*),KTLAST
0408 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0409 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0410 IF (ETOT.EQ.0) CALL KTWARN('KTBEAM',100,*999)
0411 IF (ECUT.EQ.0) THEN
0412 ETSQ=1/ETOT**2
0413 ELSE
0414 ETSQ=1/ECUT**2
0415 ENDIF
0416 J=1
0417 DO 100 I=1,NUM
0418 IF (HIST(I).LE.NMAX) THEN
0419 Y(J)=ETSQ*KT(I)
0420 J=J+1
0421 ENDIF
0422 100 CONTINUE
0423 DO 200 I=J,NUM
0424 Y(I)=0
0425 200 CONTINUE
0426 RETURN
0427 999 RETURN 1
0428 END
0429
0430 SUBROUTINE KTJOIN(ECUT,YMAC,Y,*)
0431 IMPLICIT NONE
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447 INTEGER NMAX,HIST,NUM,I,J
0448 PARAMETER (NMAX=512)
0449 DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,ECUT,ETSQ,Y(*),YMAC,KTLAST,
0450 & ROUND
0451 PARAMETER (ROUND=0.99999D0)
0452 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0453 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0454 IF (ETOT.EQ.0) CALL KTWARN('KTJOIN',100,*999)
0455 IF (ECUT.EQ.0) THEN
0456 ETSQ=1/ETOT**2
0457 ELSE
0458 ETSQ=1/ECUT**2
0459 ENDIF
0460 J=1
0461 DO 100 I=1,NUM
0462 IF (HIST(I).GT.NMAX.AND.ETSQ*KTLAST(I).GE.ROUND*YMAC) THEN
0463 Y(J)=ETSQ*KT(I)
0464 J=J+1
0465 ENDIF
0466 100 CONTINUE
0467 DO 200 I=J,NUM
0468 Y(I)=0
0469 200 CONTINUE
0470 RETURN
0471 999 RETURN 1
0472 END
0473
0474 SUBROUTINE KTRECO(RECO,PP,NN,ECUT,YCUT,YMAC,PJET,JET,NJET,NSUB,*)
0475 IMPLICIT NONE
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499 INTEGER NMAX,RECO,NUM,N,NN,NJET,NSUB,JET(*),HIST,IMIN,JMIN,I,J
0500 PARAMETER (NMAX=512)
0501 DOUBLE PRECISION PP(4,*),PJET(4,*)
0502 DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,ETSQ,YCUT,YMAC,KTLAST,
0503 & ROUND
0504 PARAMETER (ROUND=0.99999D0)
0505 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0506 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0507
0508 IF (RECO.LT.1.OR.RECO.GT.3) THEN
0509 PRINT *,'RECO=',RECO
0510 CALL KTWARN('KTRECO',100,*999)
0511 ENDIF
0512
0513 N=NN
0514 IF (NUM.NE.NN) CALL KTWARN('KTRECO',101,*999)
0515 CALL KTCOPY(PP,N,P,(RECO.NE.1))
0516 IF (ECUT.EQ.0) THEN
0517 ETSQ=1/ETOT**2
0518 ELSE
0519 ETSQ=1/ECUT**2
0520 ENDIF
0521
0522 100 IF (ETSQ*KT(N).LT.ROUND*YCUT) THEN
0523 IF (HIST(N).LE.NMAX) THEN
0524 CALL KTMOVE(P,KTP,KTS,NMAX,N,HIST(N),0)
0525 ELSE
0526 IMIN=HIST(N)/NMAX
0527 JMIN=HIST(N)-IMIN*NMAX
0528 CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO)
0529 CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0)
0530 ENDIF
0531 N=N-1
0532 IF (N.GT.0) GOTO 100
0533 ENDIF
0534
0535 NJET=N
0536 NSUB=N
0537 IF (N.EQ.0) RETURN
0538
0539 DO 210 I=1,NJET
0540 IF (RECO.EQ.1) THEN
0541 DO 200 J=1,4
0542 PJET(J,I)=P(J,I)
0543 200 CONTINUE
0544 ELSE
0545 PJET(1,I)=P(6,I)*COS(P(8,I))
0546 PJET(2,I)=P(6,I)*SIN(P(8,I))
0547 PJET(3,I)=P(6,I)*SINH(P(7,I))
0548 PJET(4,I)=P(6,I)*COSH(P(7,I))
0549 ENDIF
0550 JET(I)=I
0551 210 CONTINUE
0552
0553 300 IF (ETSQ*KT(N).LT.ROUND*YMAC) THEN
0554 IF (HIST(N).LE.NMAX) THEN
0555 IMIN=0
0556 JMIN=HIST(N)
0557 NSUB=NSUB-1
0558 ELSE
0559 IMIN=HIST(N)/NMAX
0560 JMIN=HIST(N)-IMIN*NMAX
0561 IF (ETSQ*KTLAST(N).LT.ROUND*YMAC) NSUB=NSUB-1
0562 ENDIF
0563 DO 310 I=1,NJET
0564 IF (JET(I).EQ.JMIN) JET(I)=IMIN
0565 IF (JET(I).EQ.N) JET(I)=JMIN
0566 310 CONTINUE
0567 N=N-1
0568 IF (N.GT.0) GOTO 300
0569 ENDIF
0570 RETURN
0571 999 RETURN 1
0572 END
0573
0574 SUBROUTINE KTINCL(RECO,PP,NN,PJET,JET,NJET,*)
0575 IMPLICIT NONE
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593 INTEGER NMAX,RECO,NUM,N,NN,NJET,JET(*),HIST,IMIN,JMIN,I,J
0594 PARAMETER (NMAX=512)
0595 DOUBLE PRECISION PP(4,*),PJET(4,*)
0596 DOUBLE PRECISION P,KT,KTP,KTS,ETOT,RSQ,KTLAST
0597 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0598 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0599
0600 IF (RECO.LT.1.OR.RECO.GT.3) CALL KTWARN('KTINCL',100,*999)
0601
0602 N=NN
0603 IF (NUM.NE.NN) CALL KTWARN('KTINCL',101,*999)
0604 CALL KTCOPY(PP,N,P,(RECO.NE.1))
0605
0606 DO 100 I=1,NN
0607 JET(I)=I
0608 100 CONTINUE
0609
0610 NJET=0
0611 200 IF (N.GT.0) THEN
0612 IF (HIST(N).LE.NMAX) THEN
0613 IMIN=0
0614 JMIN=HIST(N)
0615 NJET=NJET+1
0616 IF (RECO.EQ.1) THEN
0617 DO 300 J=1,4
0618 PJET(J,NJET)=P(J,JMIN)
0619 300 CONTINUE
0620 ELSE
0621 PJET(1,NJET)=P(6,JMIN)*COS(P(8,JMIN))
0622 PJET(2,NJET)=P(6,JMIN)*SIN(P(8,JMIN))
0623 PJET(3,NJET)=P(6,JMIN)*SINH(P(7,JMIN))
0624 PJET(4,NJET)=P(6,JMIN)*COSH(P(7,JMIN))
0625 ENDIF
0626 CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0)
0627 ELSE
0628 IMIN=HIST(N)/NMAX
0629 JMIN=HIST(N)-IMIN*NMAX
0630 CALL KTMERG(P,KTP,KTS,NMAX,IMIN,JMIN,N,0,0,0,RECO)
0631 CALL KTMOVE(P,KTP,KTS,NMAX,N,JMIN,0)
0632 ENDIF
0633 DO 400 I=1,NN
0634 IF (JET(I).EQ.JMIN) JET(I)=IMIN
0635 IF (JET(I).EQ.N) JET(I)=JMIN
0636 IF (JET(I).EQ.0) JET(I)=-NJET
0637 400 CONTINUE
0638 N=N-1
0639 GOTO 200
0640 ENDIF
0641
0642 DO 500 I=1,NN
0643
0644 IF (JET(I).GE.0) CALL KTWARN('KTINCL',102,*999)
0645 JET(I)=-JET(I)
0646 500 CONTINUE
0647 RETURN
0648 999 RETURN 1
0649 END
0650
0651 SUBROUTINE KTISUB(N,NY,YCUT,NSUB,*)
0652 IMPLICIT NONE
0653
0654
0655
0656
0657
0658
0659
0660
0661
0662
0663
0664
0665 INTEGER N,NY,NSUB(NY),NMAX,HIST,I,J,NUM,NM
0666 PARAMETER (NMAX=512)
0667 DOUBLE PRECISION YCUT(NY),ETOT,RSQ,P,KT,KTP,KTS,KTLAST,ROUND,EPS
0668 PARAMETER (ROUND=0.99999D0)
0669 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0670 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0671 DATA EPS/1D-6/
0672 DO 100 I=1,NY
0673 NSUB(I)=0
0674 100 CONTINUE
0675
0676 NM=0
0677 J=0
0678 DO 110 I=NUM,1,-1
0679 IF (HIST(I).LE.NMAX) J=J+1
0680 IF (J.EQ.N) THEN
0681 NM=I
0682 GOTO 120
0683 ENDIF
0684 110 CONTINUE
0685 120 CONTINUE
0686
0687 IF (NM.EQ.0) CALL KTWARN('KTISUB',100,*999)
0688 DO 210 I=NUM,1,-1
0689 DO 200 J=1,NY
0690 IF (NSUB(J).EQ.0.AND.RSQ*KT(I).GE.ROUND*YCUT(J)*KT(NM))
0691 & NSUB(J)=I
0692 IF (NSUB(J).NE.0.AND.ABS(KTLAST(I)-KTLAST(NM)).GT.EPS)
0693 & NSUB(J)=NSUB(J)-1
0694 200 CONTINUE
0695 210 CONTINUE
0696 RETURN
0697 999 RETURN 1
0698 END
0699
0700 SUBROUTINE KTIJOI(N,Y,*)
0701 IMPLICIT NONE
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713 INTEGER NMAX,HIST,NUM,I,J,N,NM
0714 PARAMETER (NMAX=512)
0715 DOUBLE PRECISION ETOT,RSQ,P,KT,KTP,KTS,Y(*),KTLAST,EPS
0716 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0717 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0718 DATA EPS/1D-6/
0719
0720 NM=0
0721 J=0
0722 DO 100 I=NUM,1,-1
0723 IF (HIST(I).LE.NMAX) J=J+1
0724 IF (J.EQ.N) THEN
0725 NM=I
0726 GOTO 105
0727 ENDIF
0728 100 CONTINUE
0729 105 CONTINUE
0730
0731 IF (NM.EQ.0) CALL KTWARN('KTIJOI',100,*999)
0732 J=1
0733 DO 110 I=1,NUM
0734 IF (HIST(I).GT.NMAX.AND.ABS(KTLAST(I)-KTLAST(NM)).LT.EPS) THEN
0735 Y(J)=RSQ*KT(I)/KT(NM)
0736 J=J+1
0737 ENDIF
0738 110 CONTINUE
0739 DO 200 I=J,NUM
0740 Y(I)=0
0741 200 CONTINUE
0742 RETURN
0743 999 RETURN 1
0744 END
0745
0746 SUBROUTINE KTIREC(RECO,PP,NN,N,YCUT,PSUB,NSUB,*)
0747 IMPLICIT NONE
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764 INTEGER NMAX,RECO,NUM,NN,NJET,NSUB,JET,HIST,I,J,N,NM
0765 PARAMETER (NMAX=512)
0766 DOUBLE PRECISION PP(4,*),PSUB(4,*)
0767 DOUBLE PRECISION ECUT,P,KT,KTP,KTS,ETOT,RSQ,YCUT,YMAC,KTLAST
0768 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0769 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0770 DIMENSION JET(NMAX)
0771
0772 NM=0
0773 J=0
0774 DO 100 I=NUM,1,-1
0775 IF (HIST(I).LE.NMAX) J=J+1
0776 IF (J.EQ.N) THEN
0777 NM=I
0778 GOTO 110
0779 ENDIF
0780 100 CONTINUE
0781 110 CONTINUE
0782
0783 IF (NM.EQ.0) CALL KTWARN('KTIREC',102,*999)
0784
0785 ECUT=SQRT(KT(NM)/RSQ)
0786 YMAC=RSQ
0787 CALL KTRECO(RECO,PP,NN,ECUT,YCUT,YMAC,PSUB,JET,NJET,NSUB,*999)
0788
0789 NSUB=0
0790 DO 210 I=1,NJET
0791 IF (JET(I).EQ.HIST(NM)) THEN
0792 NSUB=NSUB+1
0793 DO 200 J=1,4
0794 PSUB(J,NSUB)=PSUB(J,I)
0795 200 CONTINUE
0796 ENDIF
0797 210 CONTINUE
0798 RETURN
0799 999 RETURN 1
0800 END
0801
0802 SUBROUTINE KTWICH(ECUT,YCUT,JET,NJET,*)
0803 IMPLICIT NONE
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813
0814
0815
0816
0817
0818 INTEGER JET(*),NJET,NTEMP
0819 DOUBLE PRECISION ECUT,YCUT
0820 CALL KTWCHS(ECUT,YCUT,YCUT,JET,NJET,NTEMP,*999)
0821 RETURN
0822 999 RETURN 1
0823 END
0824
0825 SUBROUTINE KTWCHS(ECUT,YCUT,YMAC,JET,NJET,NSUB,*)
0826 IMPLICIT NONE
0827
0828
0829
0830
0831
0832
0833
0834
0835
0836
0837
0838
0839
0840
0841
0842
0843
0844
0845
0846 INTEGER NMAX,JET(*),NJET,NSUB,HIST,NUM,I,J,JSUB
0847 PARAMETER (NMAX=512)
0848 DOUBLE PRECISION P1(4,NMAX),P2(4,NMAX)
0849 DOUBLE PRECISION ECUT,YCUT,YMAC,ZERO,ETOT,RSQ,P,KTP,KTS,KT,KTLAST
0850 COMMON /KTCOMM/ETOT,RSQ,P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),
0851 & KT(NMAX),KTLAST(NMAX),HIST(NMAX),NUM
0852 DIMENSION JSUB(NMAX)
0853
0854
0855 DATA ((P1(J,I),I=1,NMAX),J=1,4),ZERO
0856 & /NMAX*1,NMAX*0,NMAX*0,NMAX*1,0/
0857
0858 CALL KTRECO(1,P1,NUM,ECUT,ZERO,YCUT,P2,JET,NJET,NSUB,*999)
0859
0860 CALL KTRECO(1,P1,NUM,ECUT,YCUT,YMAC,P2,JSUB,NJET,NSUB,*999)
0861
0862 DO 10 I=1,NUM
0863 IF (JET(I).NE.0) THEN
0864 IF (JSUB(JET(I)).EQ.0) JET(I)=0
0865 ENDIF
0866 10 CONTINUE
0867 RETURN
0868 999 RETURN 1
0869 END
0870
0871 SUBROUTINE KTFRAM(IOPT,CMF,SIGN,Z,XZ,N,P,Q,*)
0872 IMPLICIT NONE
0873
0874
0875
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897 INTEGER IOPT,I,N
0898 DOUBLE PRECISION CMF(4),SIGN,Z(4),XZ(4),P(4,N),Q(4,N),
0899 & R(4,4),NEW(4),OLD(4)
0900 IF (IOPT.LT.0.OR.IOPT.GT.1) CALL KTWARN('KTFRAM',200,*999)
0901
0902 CALL KTUNIT(R)
0903 CALL KTLBST(0,R,CMF,*999)
0904
0905 IF (SIGN.NE.0) THEN
0906 CALL KTVMUL(R,Z,OLD)
0907 IF (OLD(1).NE.0.OR.OLD(2).NE.0.OR.OLD(3).NE.0) THEN
0908 NEW(1)=0
0909 NEW(2)=0
0910 NEW(3)=SIGN
0911 NEW(4)=ABS(SIGN)
0912 CALL KTRROT(R,OLD,NEW,*999)
0913
0914 CALL KTVMUL(R,XZ,OLD)
0915 IF (OLD(1).NE.0.OR.OLD(2).NE.0) THEN
0916 NEW(1)=1
0917 NEW(2)=0
0918 NEW(3)=0
0919 NEW(4)=1
0920 OLD(3)=0
0921
0922
0923
0924 CALL KTRROT(R,OLD,NEW,*999)
0925 ENDIF
0926 ENDIF
0927 ENDIF
0928
0929 IF (IOPT.EQ.1) CALL KTINVT(R,R)
0930
0931 DO 30 I=1,N
0932 CALL KTVMUL(R,P(1,I),Q(1,I))
0933 30 CONTINUE
0934 RETURN
0935 999 RETURN 1
0936 END
0937
0938 SUBROUTINE KTBREI(IOPT,PLEP,PHAD,POUT,N,P,Q,*)
0939 IMPLICIT NONE
0940
0941
0942
0943
0944
0945
0946
0947
0948
0949
0950
0951
0952
0953
0954
0955
0956
0957
0958
0959 INTEGER IOPT,N
0960 DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N),
0961 & CMF(4),Z(4),XZ(4),DOT,QDQ
0962
0963 IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTBREI',200,*999)
0964
0965 DOT=ABS(PHAD)*(ABS(PLEP)-POUT(4))-PHAD*(PLEP-POUT(3))
0966 QDQ=(ABS(PLEP)-POUT(4))**2-(PLEP-POUT(3))**2-POUT(2)**2-POUT(1)**2
0967 CMF(1)=DOT*( -POUT(1))
0968 CMF(2)=DOT*( -POUT(2))
0969 CMF(3)=DOT*( PLEP -POUT(3))-QDQ* PHAD
0970 CMF(4)=DOT*(ABS(PLEP)-POUT(4))-QDQ*ABS(PHAD)
0971
0972 Z(1)=0
0973 Z(2)=0
0974 Z(3)=PHAD
0975 Z(4)=ABS(PHAD)
0976 XZ(1)=0
0977 XZ(2)=0
0978 XZ(3)=0
0979 XZ(4)=0
0980
0981 IF (IOPT.LE.1) THEN
0982 CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999)
0983 ELSE
0984 CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999)
0985 ENDIF
0986 RETURN
0987 999 RETURN 1
0988 END
0989
0990 SUBROUTINE KTHADR(IOPT,PLEP,PHAD,POUT,N,P,Q,*)
0991 IMPLICIT NONE
0992
0993
0994
0995
0996
0997
0998
0999
1000 INTEGER IOPT,N
1001 DOUBLE PRECISION PLEP,PHAD,POUT(4),P(4,N),Q(4,N),
1002 & CMF(4),Z(4),XZ(4)
1003
1004 IF (IOPT.LT.0.OR.IOPT.GT.3) CALL KTWARN('KTHADR',200,*999)
1005
1006 CMF(1)= -POUT(1)
1007 CMF(2)= -POUT(2)
1008 CMF(3)= PLEP -POUT(3)+ PHAD
1009 CMF(4)=ABS(PLEP)-POUT(4)+ABS(PHAD)
1010
1011 Z(1)=0
1012 Z(2)=0
1013 Z(3)=PHAD
1014 Z(4)=ABS(PHAD)
1015 XZ(1)=0
1016 XZ(2)=0
1017 XZ(3)=0
1018 XZ(4)=0
1019
1020 IF (IOPT.LE.1) THEN
1021 CALL KTFRAM(IOPT,CMF,PHAD,Z,XZ,N,P,Q,*999)
1022 ELSE
1023 CALL KTFRAM(IOPT-2,CMF,PHAD,Z,POUT,N,P,Q,*999)
1024 ENDIF
1025 RETURN
1026 999 RETURN 1
1027 END
1028
1029 FUNCTION KTPAIR(ANGL,P,Q,ANGLE)
1030 IMPLICIT NONE
1031
1032
1033
1034
1035
1036 INTEGER ANGL
1037 DOUBLE PRECISION P(9),Q(9),KTPAIR,R,KTMDPI,ANGLE,ETA,PHI,ESQ
1038
1039 R=ANGLE
1040 IF (ANGL.EQ.1) THEN
1041 IF (R.LE.0) R=2*(1-(P(1)*Q(1)+P(2)*Q(2)+P(3)*Q(3))*(P(5)*Q(5)))
1042 ESQ=MIN(P(4),Q(4))**2
1043 ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN
1044 IF (R.LE.0) THEN
1045 ETA=P(7)-Q(7)
1046 PHI=KTMDPI(P(8)-Q(8))
1047 IF (ANGL.EQ.2) THEN
1048 R=ETA**2+PHI**2
1049 ELSE
1050 R=2*(COSH(ETA)-COS(PHI))
1051 ENDIF
1052 ENDIF
1053 ESQ=MIN(P(9),Q(9))
1054 ELSEIF (ANGL.EQ.4) THEN
1055 ESQ=(1d0/(P(5)*Q(5))-P(1)*Q(1)-P(2)*Q(2)-
1056 &P(3)*Q(3))*2D0/(P(5)*Q(5))/(0.0001D0+1d0/P(5)+1d0/Q(5))**2
1057 R=1d0
1058 ELSE
1059 ktpair = 0D0
1060 CALL KTWARN('KTPAIR',200,*999)
1061 STOP
1062 ENDIF
1063 KTPAIR=ESQ*R
1064 IF (ANGLE.LT.0) ANGLE=R
1065 999 END
1066
1067 FUNCTION KTSING(ANGL,TYPE,P)
1068 IMPLICIT NONE
1069
1070
1071
1072
1073
1074 DOUBLE PRECISION KTSING,P(9)
1075 DOUBLE PRECISION COSTH,R,SMALL
1076 INTEGER ANGL,TYPE
1077 DATA SMALL/1D-4/
1078 IF (ANGL.EQ.1.OR.ANGL.EQ.4) THEN
1079 COSTH=P(3)*P(5)
1080 IF (TYPE.EQ.2) THEN
1081 COSTH=-COSTH
1082 ELSEIF (TYPE.EQ.4) THEN
1083 COSTH=ABS(COSTH)
1084 ELSEIF (TYPE.NE.1.AND.TYPE.NE.3) THEN
1085 ktsing = 0D0
1086 CALL KTWARN('KTSING',200,*999)
1087 STOP
1088 ENDIF
1089 R=2*(1-COSTH)
1090
1091 IF (R.LT.SMALL) R=(P(1)**2+P(2)**2)*P(5)**2
1092 KTSING=P(4)**2*R
1093 ELSEIF (ANGL.EQ.2.OR.ANGL.EQ.3) THEN
1094 KTSING=P(9)
1095 ELSE
1096 ktsing = 0D0
1097 CALL KTWARN('KTSING',201,*999)
1098 STOP
1099 ENDIF
1100 999 END
1101
1102 SUBROUTINE KTPMIN(A,NMAX,N,IMIN,JMIN)
1103 IMPLICIT NONE
1104
1105 INTEGER NMAX,N,IMIN,JMIN,KMIN,I,J,K
1106
1107
1108 DOUBLE PRECISION A(*),AMIN
1109 K=1+NMAX
1110 KMIN=K
1111 AMIN=A(KMIN)
1112 DO 110 J=0,N-2
1113 DO 100 I=K,K+J
1114 IF (A(I).LT.AMIN) THEN
1115 KMIN=I
1116 AMIN=A(KMIN)
1117 ENDIF
1118 100 CONTINUE
1119 K=K+NMAX
1120 110 CONTINUE
1121 JMIN=KMIN/NMAX+1
1122 IMIN=KMIN-(JMIN-1)*NMAX
1123 END
1124
1125 SUBROUTINE KTSMIN(A,NMAX,N,IMIN)
1126 IMPLICIT NONE
1127
1128 INTEGER N,NMAX,IMIN,I
1129 DOUBLE PRECISION A(NMAX)
1130 IMIN=1
1131 DO 100 I=1,N
1132 IF (A(I).LT.A(IMIN)) IMIN=I
1133 100 CONTINUE
1134 END
1135
1136 SUBROUTINE KTCOPY(A,N,B,ONSHLL)
1137 IMPLICIT NONE
1138
1139
1140 INTEGER I,N
1141 DOUBLE PRECISION A(4,N)
1142 LOGICAL ONSHLL
1143 DOUBLE PRECISION B(9,N),ETAMAX,SINMIN,EPS
1144 DATA ETAMAX,SINMIN,EPS/10,0,1D-6/
1145
1146 IF (SINMIN.EQ.0) SINMIN=1/COSH(ETAMAX)
1147 DO 100 I=1,N
1148 B(1,I)=A(1,I)
1149 B(2,I)=A(2,I)
1150 B(3,I)=A(3,I)
1151 B(4,I)=A(4,I)
1152 B(5,I)=SQRT(A(1,I)**2+A(2,I)**2+A(3,I)**2)
1153 IF (ONSHLL) B(4,I)=B(5,I)
1154 IF (B(5,I).EQ.0) B(5,I)=1D-10
1155 B(5,I)=1/B(5,I)
1156 B(9,I)=A(1,I)**2+A(2,I)**2
1157 B(6,I)=SQRT(B(9,I))
1158 B(7,I)=B(6,I)*B(5,I)
1159 IF (B(7,I).GT.SINMIN) THEN
1160 B(7,I)=A(4,I)**2-A(3,I)**2
1161 IF (B(7,I).LE.EPS*B(4,I)**2.OR.ONSHLL) B(7,I)=B(9,I)
1162 B(7,I)=LOG((B(4,I)+ABS(B(3,I)))**2/B(7,I))/2
1163 ELSE
1164 B(7,I)=ETAMAX+2
1165 ENDIF
1166 B(7,I)=SIGN(B(7,I),B(3,I))
1167 IF (A(1,I).EQ.0 .AND. A(2,I).EQ.0) THEN
1168 B(8,I)=0
1169 ELSE
1170 B(8,I)=ATAN2(A(2,I),A(1,I))
1171 ENDIF
1172 100 CONTINUE
1173 END
1174
1175 SUBROUTINE KTMERG(P,KTP,KTS,NMAX,I,J,N,TYPE,ANGL,MONO,RECO)
1176 IMPLICIT NONE
1177
1178
1179
1180
1181
1182 INTEGER ANGL,RECO,TYPE,I,J,K,N,NMAX,MONO
1183 DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX),PT,PTT,
1184 & KTMDPI,KTUP,PI,PJ,ANG,KTPAIR,KTSING,ETAMAX,EPS
1185 KTUP(I,J)=KTP(MAX(I,J),MIN(I,J))
1186 DATA ETAMAX,EPS/10,1D-6/
1187 IF (J.LE.I) CALL KTWARN('KTMERG',200,*999)
1188
1189 IF (MONO.GT.1) THEN
1190 DO 100 K=1,N
1191 IF (K.NE.I.AND.K.NE.J) THEN
1192 IF (RECO.EQ.1) THEN
1193 PI=P(4,I)
1194 PJ=P(4,J)
1195 ELSEIF (RECO.EQ.2) THEN
1196 PI=P(6,I)
1197 PJ=P(6,J)
1198 ELSEIF (RECO.EQ.3) THEN
1199 PI=P(9,I)
1200 PJ=P(9,J)
1201 ELSE
1202 CALL KTWARN('KTMERG',201,*999)
1203 STOP
1204 ENDIF
1205 IF (PI.EQ.0.AND.PJ.EQ.0) THEN
1206 PI=1
1207 PJ=1
1208 ENDIF
1209 KTP(MAX(I,K),MIN(I,K))=
1210 & (PI*KTUP(I,K)+PJ*KTUP(J,K))/(PI+PJ)
1211 ENDIF
1212 100 CONTINUE
1213 ENDIF
1214 IF (RECO.EQ.1) THEN
1215
1216 P(1,I)=P(1,I)+P(1,J)
1217 P(2,I)=P(2,I)+P(2,J)
1218 P(3,I)=P(3,I)+P(3,J)
1219
1220 P(5,I)=SQRT(P(1,I)**2+P(2,I)**2+P(3,I)**2)
1221 P(4,I)=P(5,I) ! JA (Massless scheme)
1222 IF (P(5,I).EQ.0) THEN
1223 P(5,I)=1
1224 ELSE
1225 P(5,I)=1/P(5,I)
1226 ENDIF
1227 ELSEIF (RECO.EQ.2) THEN
1228
1229 PT=P(6,I)+P(6,J)
1230 IF (PT.EQ.0) THEN
1231 PTT=1
1232 ELSE
1233 PTT=1/PT
1234 ENDIF
1235 P(7,I)=(P(6,I)*P(7,I)+P(6,J)*P(7,J))*PTT
1236 P(8,I)=KTMDPI(P(8,I)+P(6,J)*PTT*KTMDPI(P(8,J)-P(8,I)))
1237 P(6,I)=PT
1238 P(9,I)=PT**2
1239 ELSEIF (RECO.EQ.3) THEN
1240
1241 PT=P(9,I)+P(9,J)
1242 IF (PT.EQ.0) THEN
1243 PTT=1
1244 ELSE
1245 PTT=1/PT
1246 ENDIF
1247 P(7,I)=(P(9,I)*P(7,I)+P(9,J)*P(7,J))*PTT
1248 P(8,I)=KTMDPI(P(8,I)+P(9,J)*PTT*KTMDPI(P(8,J)-P(8,I)))
1249 P(6,I)=P(6,I)+P(6,J)
1250 P(9,I)=P(6,I)**2
1251 ELSE
1252 CALL KTWARN('KTMERG',202,*999)
1253 STOP
1254 ENDIF
1255
1256 IF (MONO.LE.0) RETURN
1257
1258 IF (ANGL.NE.1.AND.RECO.EQ.1) THEN
1259 P(9,I)=P(1,I)**2+P(2,I)**2
1260 P(7,I)=P(4,I)**2-P(3,I)**2
1261 IF (P(7,I).LE.EPS*P(4,I)**2) P(7,I)=P(9,I)
1262 IF (P(7,I).GT.0) THEN
1263 P(7,I)=LOG((P(4,I)+ABS(P(3,I)))**2/P(7,I))/2
1264 IF (P(7,I).GT.ETAMAX) P(7,I)=ETAMAX+2
1265 ELSE
1266 P(7,I)=ETAMAX+2
1267 ENDIF
1268 P(7,I)=SIGN(P(7,I),P(3,I))
1269 IF (P(1,I).NE.0.AND.P(2,I).NE.0) THEN
1270 P(8,I)=ATAN2(P(2,I),P(1,I))
1271 ELSE
1272 P(8,I)=0
1273 ENDIF
1274 ELSEIF (ANGL.EQ.1.AND.RECO.NE.1) THEN
1275 P(1,I)=P(6,I)*COS(P(8,I))
1276 P(2,I)=P(6,I)*SIN(P(8,I))
1277 P(3,I)=P(6,I)*SINH(P(7,I))
1278 P(4,I)=P(6,I)*COSH(P(7,I))
1279 IF (P(4,I).NE.0) THEN
1280 P(5,I)=1/P(4,I)
1281 ELSE
1282 P(5,I)=1
1283 ENDIF
1284 ENDIF
1285 ANG=0
1286 DO 200 K=1,N
1287 IF (K.NE.I.AND.K.NE.J) THEN
1288 IF (MONO.GT.1) ANG=KTUP(I,K)
1289 KTP(MIN(I,K),MAX(I,K))=
1290 & KTPAIR(ANGL,P(1,I),P(1,K),ANG)
1291 ENDIF
1292 200 CONTINUE
1293 KTS(I)=KTSING(ANGL,TYPE,P(1,I))
1294 999 END
1295
1296 SUBROUTINE KTMOVE(P,KTP,KTS,NMAX,N,J,IOPT)
1297 IMPLICIT NONE
1298
1299
1300 INTEGER I,J,N,NMAX,IOPT
1301 DOUBLE PRECISION P(9,NMAX),KTP(NMAX,NMAX),KTS(NMAX)
1302 DO 100 I=1,9
1303 P(I,J)=P(I,N)
1304 100 CONTINUE
1305 IF (IOPT.LE.0) RETURN
1306 DO 110 I=1,J-1
1307 KTP(I,J)=KTP(I,N)
1308 KTP(J,I)=KTP(N,I)
1309 110 CONTINUE
1310 DO 120 I=J+1,N-1
1311 KTP(J,I)=KTP(I,N)
1312 KTP(I,J)=KTP(N,I)
1313 120 CONTINUE
1314 KTS(J)=KTS(N)
1315 END
1316
1317 SUBROUTINE KTUNIT(R)
1318 IMPLICIT NONE
1319
1320 DOUBLE PRECISION R(4,4)
1321 INTEGER I,J
1322 DO 20 I=1,4
1323 DO 10 J=1,4
1324 R(I,J)=0
1325 IF (I.EQ.J) R(I,J)=1
1326 10 CONTINUE
1327 20 CONTINUE
1328 END
1329
1330 SUBROUTINE KTLBST(IOPT,R,A,*)
1331 IMPLICIT NONE
1332
1333
1334
1335
1336
1337
1338
1339 INTEGER IOPT,I,J
1340 DOUBLE PRECISION R(4,4),A(4),B(4),C(4,4),M
1341 DO 10 I=1,4
1342 B(I)=A(I)
1343 10 CONTINUE
1344 M=B(4)**2-B(1)**2-B(2)**2-B(3)**2
1345 IF (M.LE.0) CALL KTWARN('KTLBST',100,*999)
1346 M=SQRT(M)
1347 B(4)=B(4)+M
1348 M=1/(M*B(4))
1349 IF (IOPT.EQ.0) THEN
1350 B(4)=-B(4)
1351 ELSEIF (IOPT.NE.1) THEN
1352 CALL KTWARN('KTLBST',200,*999)
1353 STOP
1354 ENDIF
1355 DO 30 I=1,4
1356 DO 20 J=1,4
1357 C(I,J)=B(I)*B(J)*M
1358 IF (I.EQ.J) C(I,J)=C(I,J)+1
1359 20 CONTINUE
1360 30 CONTINUE
1361 C(4,4)=C(4,4)-2
1362 CALL KTMMUL(C,R,R)
1363 RETURN
1364 999 RETURN 1
1365 END
1366
1367 SUBROUTINE KTRROT(R,A,B,*)
1368 IMPLICIT NONE
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380 DOUBLE PRECISION R(4,4),M(4,4),A(4),B(4),C(4),D(4),AL,BL,CL,DL,EPS
1381
1382
1383
1384 PARAMETER (EPS=0.5D-6)
1385 AL=A(1)**2+A(2)**2+A(3)**2
1386 BL=B(1)**2+B(2)**2+B(3)**2
1387 IF (AL.LE.0.OR.BL.LE.0) CALL KTWARN('KTRROT',100,*999)
1388 AL=1/SQRT(AL)
1389 BL=1/SQRT(BL)
1390 CL=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*AL*BL
1391
1392 IF (CL.GE.1-EPS) THEN
1393 RETURN
1394
1395 ELSEIF (CL.LE.-1+EPS) THEN
1396 IF (ABS(B(2)).GT.EPS) THEN
1397 C(1)= 0
1398 C(2)=-B(3)
1399 C(3)= B(2)
1400
1401 ELSE
1402 C(1)= B(3)
1403 C(2)= 0
1404 C(3)=-B(1)
1405 ENDIF
1406
1407 ELSE
1408 C(1)=A(2)*B(3)-A(3)*B(2)
1409 C(2)=A(3)*B(1)-A(1)*B(3)
1410 C(3)=A(1)*B(2)-A(2)*B(1)
1411 ENDIF
1412 CL=C(1)**2+C(2)**2+C(3)**2
1413 IF (CL.LE.0) CALL KTWARN('KTRROT',101,*999)
1414 CL=1/SQRT(CL)
1415
1416 D(1)=A(2)*C(3)-A(3)*C(2)
1417 D(2)=A(3)*C(1)-A(1)*C(3)
1418 D(3)=A(1)*C(2)-A(2)*C(1)
1419 DL=AL*CL
1420 M(1,1)=A(1)*AL
1421 M(1,2)=A(2)*AL
1422 M(1,3)=A(3)*AL
1423 M(1,4)=0
1424 M(2,1)=C(1)*CL
1425 M(2,2)=C(2)*CL
1426 M(2,3)=C(3)*CL
1427 M(2,4)=0
1428 M(3,1)=D(1)*DL
1429 M(3,2)=D(2)*DL
1430 M(3,3)=D(3)*DL
1431 M(3,4)=0
1432 M(4,1)=0
1433 M(4,2)=0
1434 M(4,3)=0
1435 M(4,4)=1
1436 CALL KTMMUL(M,R,R)
1437
1438 D(1)=B(2)*C(3)-B(3)*C(2)
1439 D(2)=B(3)*C(1)-B(1)*C(3)
1440 D(3)=B(1)*C(2)-B(2)*C(1)
1441 DL=BL*CL
1442 M(1,1)=B(1)*BL
1443 M(2,1)=B(2)*BL
1444 M(3,1)=B(3)*BL
1445 M(1,2)=C(1)*CL
1446 M(2,2)=C(2)*CL
1447 M(3,2)=C(3)*CL
1448 M(1,3)=D(1)*DL
1449 M(2,3)=D(2)*DL
1450 M(3,3)=D(3)*DL
1451 CALL KTMMUL(M,R,R)
1452 RETURN
1453 999 RETURN 1
1454 END
1455
1456 SUBROUTINE KTVMUL(M,A,B)
1457 IMPLICIT NONE
1458
1459
1460
1461
1462 DOUBLE PRECISION M(4,4),A(4),B(4),C(4)
1463 INTEGER I,J
1464 DO 20 I=1,4
1465 C(I)=0
1466 DO 10 J=1,4
1467 C(I)=C(I)+M(I,J)*A(J)
1468 10 CONTINUE
1469 20 CONTINUE
1470 DO 30 I=1,4
1471 B(I)=C(I)
1472 30 CONTINUE
1473 END
1474
1475 SUBROUTINE KTMMUL(A,B,C)
1476 IMPLICIT NONE
1477
1478
1479
1480
1481 DOUBLE PRECISION A(4,4),B(4,4),C(4,4),D(4,4)
1482 INTEGER I,J,K
1483 DO 30 I=1,4
1484 DO 20 J=1,4
1485 D(I,J)=0
1486 DO 10 K=1,4
1487 D(I,J)=D(I,J)+A(I,K)*B(K,J)
1488 10 CONTINUE
1489 20 CONTINUE
1490 30 CONTINUE
1491 DO 50 I=1,4
1492 DO 40 J=1,4
1493 C(I,J)=D(I,J)
1494 40 CONTINUE
1495 50 CONTINUE
1496 END
1497
1498 SUBROUTINE KTINVT(A,B)
1499 IMPLICIT NONE
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509 DOUBLE PRECISION A(4,4),B(4,4),C(4,4)
1510 INTEGER I,J
1511
1512 DO 20 I=1,4
1513 DO 10 J=1,4
1514 C(I,J)=A(J,I)
1515 10 CONTINUE
1516 20 CONTINUE
1517
1518 DO 30 I=1,3
1519 C(4,I)=-C(4,I)
1520 C(I,4)=-C(I,4)
1521 30 CONTINUE
1522
1523 DO 50 I=1,4
1524 DO 40 J=1,4
1525 B(I,J)=C(I,J)
1526 40 CONTINUE
1527 50 CONTINUE
1528 END
1529
1530 FUNCTION KTMDPI(PHI)
1531 IMPLICIT NONE
1532
1533 DOUBLE PRECISION KTMDPI,PHI,PI,TWOPI,THRPI,EPS
1534 PARAMETER (PI=3.14159265358979324D0,TWOPI=6.28318530717958648D0,
1535 & THRPI=9.42477796076937972D0)
1536 PARAMETER (EPS=1D-15)
1537 KTMDPI=PHI
1538 IF (KTMDPI.LE.PI) THEN
1539 IF (KTMDPI.GT.-PI) THEN
1540 GOTO 100
1541 ELSEIF (KTMDPI.GT.-THRPI) THEN
1542 KTMDPI=KTMDPI+TWOPI
1543 ELSE
1544 KTMDPI=-MOD(PI-KTMDPI,TWOPI)+PI
1545 ENDIF
1546 ELSEIF (KTMDPI.LE.THRPI) THEN
1547 KTMDPI=KTMDPI-TWOPI
1548 ELSE
1549 KTMDPI=MOD(PI+KTMDPI,TWOPI)-PI
1550 ENDIF
1551 100 IF (ABS(KTMDPI).LT.EPS) KTMDPI=0
1552 END
1553
1554 SUBROUTINE KTWARN(SUBRTN,ICODE,*)
1555
1556
1557
1558
1559
1560
1561 INTEGER ICODE
1562 CHARACTER*6 SUBRTN
1563 WRITE (6,10) SUBRTN,ICODE
1564 10 FORMAT(/' KTWARN CALLED FROM SUBPROGRAM ',A6,': CODE =',I4/)
1565 IF (ICODE.LT.100) RETURN
1566 IF (ICODE.LT.200) RETURN 1
1567 STOP
1568 END
1569
1570
1571