SUBROUTINE EVENT(NERR)
C ****************** COMMON COMMON *************************** 0030
COMMON MAP(2000),PARS(1000),MISC(27),KLIST(500),MTABLE(2)
DIMENSION ZMAP(2000)
DIMENSION HTABLE(7, 100), KTABLE(7, 100), OTABLE(7, 50)
DIMENSION RTABLE(9, 20, 2), ITABLE(6, 20), LTABLE(9, 20, 2)
DIMENSION JTABLE(7, 50), RBANK(10), DCTGT(3), DCTRA(3)
DIMENSION TITLE(48), NC(48), DIRINC(3)
DIMENSION PARA(1000),NPARA(1000),SNAME(1000),NAME(1000),TABLE(100)
DIMENSION TBLMS (30), ADIS (40), ZEROS (10) , AD(50)
DIMENSION NTABLE(100), HEAD(11), NBRNCH(10)
DIMENSION OBANK(90), IBANK(90), OB(90),VM(3), VN(3)
DIMENSION PPLANE(3),FMASS(8),PTBL(8),ETBL(8),DCSO(3,8),DCSINC(3)
EQUIVALENCE (KTABLE,HTABLE,MAP)
EQUIVALENCE (OTABLE,JTABLE,MAP(701)), (RTABLE,LTABLE,MAP(1051)),
1 (ITABLE,MAP(1411)), (VAL,IVAL,MAP(1531)),
2 (WGT,MAP(1631))
EQUIVALENCE ( MAP (1731), NC(1) ),( MAP (1870), TITLE(1) )
EQUIVALENCE (NCFLAG,MAP(1869)), (WEIGHT,MAP(1978)),
1 (NTAPE,MAP(1988)), (EINC,MAP(1998)),
2 (PINC,MAP(1999)), (BINC,MAP(2000))
EQUIVALENCE (IPS,MAP(1973)), (WPS,MAP(1974))
EQUIVALENCE (MAP(1975), GMINP), (MAP(1976), GMAXP)
EQUIVALENCE (MAP(1977), GSCALE), (WSCALE,MAP(1972))
EQUIVALENCE (MTOT,MAP(1987)) 7/13/68
EQUIVALENCE (PARA,NPARA,PARS),(SNAME,NAME,MAP(1))
EQUIVALENCE (ZMAP,MAP),(TABLE,PARS(101),NTABLE)
EQUIVALENCE ( TBLMS (1), PARA (1) ),( NTOT, NPARA (99) )
EQUIVALENCE ( AD (1), PARA (40) ),( ZEROS (1), PARA (40) )
EQUIVALENCE (RSCALE, PARA(98)), (DIRINC(1), PARA(95))
EQUIVALENCE (ADIS(1),AD(11)), (OBANK(1),OB(1))
EQUIVALENCE (OBANK(1), MAP(1779)),(OBANK(1), IBANK(1))
EQUIVALENCE (OBANK(1), PPLANE(1)), (OBANK(4), FMASS(1))
EQUIVALENCE (OBANK(12), PTBL(1)), (OBANK(20), ETBL(1) )
EQUIVALENCE (OBANK(28), DCSINC(1)),(OBANK(31), DCSO(1))
EQUIVALENCE (OBANK(55), VM(1)), (OBANK(58), VN(1))
EQUIVALENCE ( OBANK(61), DCTGT(1)), (OBANK(64), RBANK(1 ))
EQUIVALENCE (OBANK(74), DCTRA(1)) , (BENUC, PARA(94))
EQUIVALENCE (PI, MISC), (RADIAN, MISC(2)), (NIT, MISC(3)),
1 (NOT, MISC(4)), (HEAD, MISC(5)), (NBRNCH, MISC(16)),
2 (NPAGE, MISC(26)), (NORD, MISC(27))
C *************** END COMMON COMMON *************************
NTRY = 0
10 NTRY = NTRY+1
C TRY 50 TIMES BEFORE GIVING UP --- THE FAILURE MAY HAVE BEEN
C THE RESULT OF AN UNPHYSICAL THROW FOR A RESONANCE
IF (NTRY-50) 12, 12, 500
12 CONTINUE
NERR = 0 0610
WPS = 1.0 0620
IPS = IPS 0630
C FILL UP RESONANCE BANK 0640
DO 50 K = 20, 29 0650
KK = K + 10 0660
LL = K - 19 0670
CALL RANRES (TBLMS(K), TBLMS(KK), RBANK(LL)) 0680
50 CONTINUE 0690
DO 150 K = 1,20 0700
NPROD = ITABLE(1,K) 0710
IF (NPROD -1) 150,150,75 0720
75 TMASS = RTABLE(9,K,1) 0730
IF (TMASS) 800, 76, 76
C A NEGATIVE MASS
800 NERR = 23
GO TO 10
76 CALL FERMI (RTABLE(9, K, 2), BENUC, TP, DCTGT(1)) 0740
DO 80 KK = 1, NPROD 0750
CALL GETMS(LTABLE(KK,K,1), FMASS(KK)) 0760
IF(FMASS(KK))800,80,80 0770
80 CONTINUE 0800
C TEST FOR VERTEX A 0810
IF (K-1) 90,85,90 0820
85 DO 86 KK = 1,3 0830
86 DCSINC(KK) = DIRINC(KK) 0840
BMASS = BINC 0850
EB = EINC 0860
IF (BMASS) 892, 87, 892
87 IF (GMINP) 892, 892, 89
C INVOKE BREMMSTRAHLUNG SPECTRUM 0910
89 CALL RANDOM (RAN)
PINC = GMINP * EXP( ABS( GSCALE * RAN) ) 0930
EINC = PINC 0940
EB = PINC 0950
C FILL OTABLE(I,1) WITH BEAM QUANTITIES
892 DO 894 KK=1,3
894 OTABLE(KK,1) = DCSINC(KK)
OTABLE(4,1) = PINC
OTABLE(5,1) = EINC
GO TO 100 0960
90 II = ITABLE(3,K) 0970
KK = ITABLE(4,K) 0980
CALL GETMS (LTABLE(KK, II, 1), BMASS) 0990
JJ = ITABLE(2,II) + KK 1000
DO 92 KK = 1,3 1010
92 DCSINC(KK) = OTABLE(KK,JJ) 1020
EB = OTABLE(5,JJ) 1030
100 IF (BMASS) 800, 1002, 1002
1002 IF (EB) 108, 1004, 1004
1004 CONTINUE
C GET ADIS ENTRY 1050
LL = ITABLE(5, K) + 1 1060
C SET PPLANE 1070
IF (ITABLE(6,K)) 103,103,105
103 DO 104 KK = 1,3 1120
104 PPLANE(KK) = 0.0 1130
GO TO 107 1140
105 MM = ITABLE(6,K) 1150
DO 106 KK= 1,3 1160
106 PPLANE(KK) = OTABLE(KK, MM) 1170
107 CONTINUE 1180
CALL PSEGEN (TMASS, TP, DCTGT(1), BMASS, EB, DCSINC(1), NPROD, FMA 1200
1SS(1), PTBL(1), ETBL(1), DCSO(1,1), AD(LL), PPLANE(1), WATE) 1210
IF (ETBL(1) + 2.0) 108,108,109 1220
C NOT ENOUGH ENERGY
108 NERR = 21 1230
GO TO 10
109 IF (ETBL(1) ) 110,110,115 1250
C REPEATED ZERO MOMENTA (OR DIVIDE CHECKS)
110 NERR = 22 1260
GO TO 500 1270
115 GO TO (118, 116), IPS 1280
116 WPS = WPS * WATE 1290
118 CONTINUE 1300
DO 120 NN = 1, NPROD 1310
II = ITABLE(2, K) + NN 1320
OTABLE(5, II) = ETBL(NN) 1330
OTABLE(4,II) = PTBL(NN) 1340
DO 120 KK = 1,3 1350
120 OTABLE(KK,II) = DCSO(KK,NN) 1360
150 CONTINUE 1370
500 CONTINUE 1380
WPS = WPS*WSCALE 1390
RETURN 1400
END 1410
SUBROUTINE PSEGEN ( TMASS, TP, DCTGT, BMASS, EB, DCSIN, NPROD, 0010
1 FMASS, PTBL, ETBL, DCSO, A, PPLANE, W ) 0040
C VERSION OF 1/69--COMPATIBLE WITH INVARIANT RDECAY 1/69
C ALSO RETURNS WEIGHT TO CONVERT TO NON-INVARIANT PHASE SPASE 1/69
DIMENSION DCSINC(3) ,FMASS(8) ,PTBL(8), DCSIN(3) 0050
DIMENSION ETBL(8) ,DCSO (3,8) ,R(3,8) , DCTGT(3) 0060
DIMENSION P(8,4) ,X(8) 2/69
DIMENSION PP(3,8) ,ROT(3,3) 1/69
DIMENSION A(10) , PPLANE(3) 0090
CALL DVCHK (K000FX) 0140
NDVCHK = -1 0150
C CALCULATE ENERGY IN CENTER OF MASS 0190
PB = SQRT( EB**2 - BMASS**2) 0200
IF (TP) 11, 3, 11 0210
3 DO 4 K = 1, 3 0220
4 DCSINC(K) = DCSIN(K) 0230
DEN = PB 0240
TE = TMASS 0250
IF (TMASS) 5, 5, 10 0260
5 OMEGA = BMASS 0270
GO TO 12 0280
10 OMEGA = SQRT( TMASS ** 2 + BMASS ** 2 + 2.0 * TMASS * EB ) 0290
GO TO 12 0300
11 DEN = 0 0310
DO 121 K = 1, 3 0320
DCSINC(K) = DCSIN(K) * PB + DCTGT(K) * TP 0330
121 DEN = DEN + DCSINC(K)**2 0340
TE = SQRT( TMASS**2 + TP**2 ) 0350
OMEGA = SQRT(( EB + TE ) **2 - DEN ) 0360
DEN = SQRT(DEN) 0370
DO 131 K = 1, 3 0380
131 DCSINC(K) = DCSINC(K) / DEN 0390
12 ESUM = OMEGA 0400
DO 13 K = 1, NPROD 0410
ESUM = ESUM - FMASS(K) 0420
I = NPROD-K+1 2/69
X(I) = FMASS(K) 2/69
13 CONTINUE 0430
BETA = DEN / ( TE + EB) 0440
GAMMA = (EB + TE ) / OMEGA 0450
IF (ESUM) 14, 15, 15 0460
C NOT ENOUGH ENERGY * EXIT AND SET ETBL(1) = -10.0 0470
14 ETBL(1) = -10.0 0480
GO TO 401 0490
15 NDVCHK = NDVCHK + 1 0500
IF (NDVCHK - 50) 17, 16, 16 0510
C EXIT IF DIVIDE CHECK 50 TIMES AND SET ETBL(1) = -1.0 0520
16 ETBL(1) = -1.0 0530
GO TO 401 0540
17 CONTINUE 0550
CALL RDECAY (NPROD, OMEGA, X, P, FLAG, A, TMASS, BMASS) 2/69
IF (FLAG) 14, 18, 16 1/69
18 CONTINUE 0580
C CALCULATE NON-INVARIANT WEIGHT 1/69
W = 1.0 0600
DO 100 I = 1, NPROD 0610
W = W*P(I,4) 1/69
100 CONTINUE 0640
IF (NPROD - 2) 401, 2001, 201 0650
2001 CONTINUE 0660
DO 2002 K = 1, 3 0670
IF (PPLANE(K) ) 2003, 2002, 2003 0680
2002 CONTINUE 0690
GO TO 201 0700
2003 CONTINUE 0710
C CALCULATE ROTATION MATRIX FOR TRANSVERSE POLARIZATION 0720
ALFC = PPLANE(2) * DCSINC(3) - PPLANE(3) * DCSINC(2) 0730
BETC = PPLANE(3) * DCSINC(1) - PPLANE(1) * DCSINC(3) 0740
GAMC = PPLANE(1) * DCSINC(2) - PPLANE(2) * DCSINC(1) 0750
C NORMALOIZE NORMAL VECTOR 0760
DEN = SQRT( ALFC**2 + BETC**2 + GAMC**2 ) 0770
IF (DEN) 2008, 2008, 2009 0780
C SET RANDOM NORMAL 0790
2008 CONTINUE 0800
CALL RANDOM ( ALFC ) 6/68
CALL RANDOM ( RSN ) 6/68
BETC = SIGN( SQRT( 1.0 - ALFC **2), RSN ) 0820
GAMC = 0.0 0830
GO TO 2007 0840
2009 ALFC = ALFC / DEN 0850
BETC = BETC / DEN 0860
GAMC = GAMC / DEN 0870
2007 CONTINUE 0880
C SET UP X-AXIS ALONG BEAM DIRECTION AND Z-AXIS ALONG POLARIZATION D 0890
COST = GAMC 0900
SINT = SQRT( 1.0 - GAMC**2) 0910
IF (SINT) 2012, 2011, 2012 0920
2011 COSP = 1.0 0930
SINP = 0.0 0940
GO TO 2014 0950
2012 COSP = ALFC/SINT 0960
SINP = BETC / SINT 0970
2014 CONTINUE 0980
ALFCP = DCSINC(1) * COST * COSP + DCSINC(2) * COST * SINP 0990
1 - DCSINC(3) * SINT 1000
BETCP = - DCSINC(1) * SINP + DCSINC(2) * COSP 1010
GAMCP = DCSINC(1) * SINT * COSP + DCSINC(2) * SINT * SINP 1020
1 + DCSINC(3) * COST 1030
C LORENTZ TRANSFORMATION ALONG THE X AXIS 1040
DO 210 K = 1, NPROD 1050
I = NPROD-K+1 2/69
PP(1, K) = GAMMA * ( P(I,1) + BETA * P(I,4) ) 2/69
PP(2, K) = P(I,2) 2/69
PP(3, K) = P(I,3) 2/69
ETBL(K) = GAMMA * ( P(I,4) + BETA * P(I,1) ) 2/69
PTBL(K) = SQRT( ETBL(K) ** 2 - FMASS(K) ** 2) 1100
210 CONTINUE 1110
ROT (1,1) =-BETCP * SINP + ALFCP * COST * COSP 1120
ROT (1, 2) = - ALFCP * SINP - BETCP * COST * COSP 1130
ROT(1, 3) = SINT * COSP 1140
ROT (2, 1) = + BETCP * COSP + ALFCP * COST * SINP 1150
ROT (2, 2) = ALFCP * COSP - BETCP * COST * SINP 1160
ROT (2,3) = SINT * SINP 1170
ROT (3,1) =-ALFCP * SINT 1180
ROT (3,2) = BETCP * SINT 1190
ROT (3,3) = COST 1200
GO TO 220 1210
201 CONTINUE 1220
C LONGITUDINAL POLARIZATION * PROD ANGULAR DIST * NO ASSYMMETRY 1230
C LORENTZ TRANSFORM ALONG THE Z AXIS 1240
DO 200 K = 1, NPROD 1250
I = NPROD-K+1 2/69
PP(1, K) = P(I,1) 2/69
PP(2, K) = P(I,2) 2/69
PP(3, K) = GAMMA * ( P(I,3) + BETA * P(I,4) ) 2/69
ETBL(K) = GAMMA * ( P(I,4) + BETA * P(I,3) ) 2/69
PTBL(K) = SQRT( ETBL(K) ** 2 - FMASS(K) ** 2) 1300
200 CONTINUE 1310
COST = DCSINC(3) 1320
SINT = SQRT( 1.0 - COST ** 2) 1330
IF (SINT) 203, 202, 203 1340
202 COSP = 1.0 1350
SINP = 0.0 1360
GO TO 204 1370
203 COSP = DCSINC(1) / SINT 1380
SINP = DCSINC(2) / SINT 1390
204 CONTINUE 1400
C NEXT ROTATE Z AXIS TO COINCIDE WITH DIRINC 1410
C SET UP ROTATION MATRIX ROT 1420
ROT(1,1) = COST * COSP 1430
ROT(1,2) = -SINP 1440
ROT(1,3) = SINT * COSP 1450
ROT(2,1) = COST * SINP 1460
ROT(2,2) = COSP 1470
ROT(2,3) = SINT * SINP 1480
ROT(3,1) = -SINT 1490
ROT(3,2) = 0.0 1500
ROT(3,3) = COST 1510
C CALCULATE MOMENTUM COMPONENTS 1520
220 CONTINUE 1530
DO 300 L = 1, NPROD 1540
DO 300 M = 1, 3 1550
R(M, L) = 0.0 1560
DO 300 N = 1, 3 1570
R(M, L) = R(M, L) + ROT(M, N) * PP(N, L) 1580
300 CONTINUE 1590
C CALCULATE DIRECTIONS FOR PRODUCT PARTICLES 1600
DO 400 I = 1, NPROD 1610
IF (PTBL(I)) 312, 15, 312 1640
312 CONTINUE 1650
DO 400 K = 1, 3 1690
DCSO(K,I) = R(K,I) / PTBL(I) 1700
400 CONTINUE 1710
CALL DVCHK (K000FX) 1712
GO TO(15,401,15),K000FX 1714
401 CONTINUE 1720
RETURN 1730
END 1740
SUBROUTINE RDECAY (N, E, X, P, FLAG, A, TMASS, BMASS)
C THIS VERSION OF 1/69 GENERATES INVARIANT PHASE SPACE EVENTS
C IT IS A MODIFICATION OF THE J.H.U. ROUTINES PSEVNT AND ISODK
C REFERENCE-- D.T. GILLESPIE, DISSERTATION(1968), APPENDICES C AND D
C JOHNS HOPKINS UNIVERSITY (HIGH ENERGY PHYSICS GROUP)
C * * * * * * * * * * * * * *
DIMENSION P(8,4), Q(8,4), X(8), Z(8) 1/69
DIMENSION ZL(8), POWER(7), A(10) 1/69
DIMENSION NSAVE(5), XSAVE(8,5), ESAVE(5), WMSAVE(5) 1/69
DOUBLE PRECISION P10(3),BETA(3),TLOR(4,4),GAMMA,DELTA,PJ3
PS(A,B,C) = ((A+B+C)*(A-B-C)*(A+B-C)*(A-B+C))/(A*A)
C * * * * * * * * * * * * * *
IF (N .GT. 0) GO TO 10 1/69
IMAX = 0 1/69
DO 4 I = 1,5 1/69
4 NSAVE(I) = 0 1/69
GO TO 400 1/69
C
C INITIALIZE.
10 FLAG = 0. 1/69
NP1 = N + 1
NM1 = N - 1
NM2 = N - 2
C DETERMINE DIST TYPE--ITYP = 0 (UNIFORM), 1 (E**AT), 2 (UNIFORM 1/69
C WITHIN LIMITS), 3 (COSINE SERIES) 1/69
ITYP = 0 1/69
IF (A(1) .NE. 0) ITYP = 1 1/69
IF (A(10) .NE. 0) ITYP = ITYP+2 1/69
DO 14 I = 2,9 1/69
14 IF (A(I) .NE. 0) ITYP = 3 1/69
DO 20 I = 1,3 1/69
20 Q(N,I) = 0. 1/69
Q(N,4) = E 1/69
Z(N) = E
Z(1) = X(1)
ZL(1) = X(1)
DO 22 I=2,N
22 ZL(I) = ZL(I-1) + X(I)
T = Z(N) - ZL(N) 1/69
IF (T .GT. 0.) GO TO 26 1/69
FLAG = -1. 1/69
GO TO 400 1/69
26 IF (NM2 .EQ. 0) GO TO 300 1/69
DO 28 I=2,NM1
FIM1 = I - 1
28 POWER(I) = 1./FIM1
C FIND WMAXS.
C FIRST SEE IF IT IS BEING SAVED 1/69
DO 86 I = 1,IMAX 1/69
IF (NSAVE(I) .NE. N) GO TO 86 1/69
DO 82 J = 1,N 1/69
IF (XSAVE(J,I) .NE. X(J)) GO TO 86 1/69
82 CONTINUE 1/69
IF (ESAVE(I) .LT. E) GO TO 84 1/69
WMAXS = WMSAVE(I) 1/69
GO TO 200 1/69
84 ISAV = I 1/69
GO TO 100 1/69
86 CONTINUE 1/69
ISAV = IMAX+1 1/69
IF (ISAV .LE. 5) IMAX = ISAV 1/69
C CALCULATE WMAXS 1/69
100 DO 104 J=1,NM2
I = N - J
104 Z(I) = ZL(I) + (Z(I+1) - ZL(I+1)) * .5**POWER(I)
WMAXS = 1. 2/69
DO 106 K=2,N 2/69
106 WMAXS = WMAXS*PS(Z(K),Z(K-1),X(K))
DZ = T/3. 1/69
108 DZ = DZ/5.
IF (DZ .LT. T/10000.) GO TO 122 1/69
110 KDZ = 0
DO 120 J=1,NM2
I = N - J
112 Z(I) = Z(I) + DZ
IF (Z(I) .GE. (Z(I+1)-X(I))) GO TO 115
WS = 1. 2/69
DO 113 K=2,N 2/69
113 WS = WS*PS(Z(K),Z(K-1),X(K))
IF (WS .LE. WMAXS) GO TO 115
WMAXS = WS
KDZ = 1
GO TO 112
115 Z(I) = Z(I) - DZ
116 Z(I) = Z(I) - DZ
IF (Z(I) .LE. (Z(I-1)+X(I))) GO TO 119
WS = 1. 2/69
DO 117 K=2,N 2/69
117 WS = WS*PS(Z(K),Z(K-1),X(K))
IF (WS .LE. WMAXS) GO TO 119
WMAXS = WS
KDZ = 1
GO TO 116
119 Z(I) = Z(I) + DZ
120 CONTINUE
IF (KDZ .EQ. 1) GO TO 110
GO TO 108
122 WMAXS = 1.001*WMAXS
C SAVE WMAXS 1/69
IF (ISAV .GT. 5) GO TO 200 1/69
NSAVE(ISAV) = N 1/69
DO 124 J = 1,N 1/69
124 XSAVE(J,ISAV) = X(J) 1/69
ESAVE(ISAV) = E 1/69
WMSAVE(ISAV) = WMAXS 1/69
C END INITIALIZATION
C
C GENERATE EVENTS. FIRST GENERATE INTERMEDIATE MASSES.
200 DO 206 J=1,NM2 1/69
IF (J.NE.1 .OR. ITYP.NE.1) GO TO 204 1/69
CALL EXPDIS (A, Z(N-1), E, TMASS, BMASS, X(N), ZL(N-1), N) 2/69
GO TO 206 1/69
204 I = N - J 1/69
CALL RANDOM (R) 1/69
Z(I) = ZL(I) + (Z(I+1) - ZL(I+1)) * ABS(R)**POWER(I) 1/69
206 CONTINUE 1/69
WS = 1. 2/69
DO 208 K=2,N 2/69
208 WS = WS*PS(Z(K),Z(K-1),X(K))
210 IF (WS .LE. WMAXS) GO TO 214
C WEIGHTING FAULT
FLAG = 1. 1/69
GO TO 400 1/69
214 CALL RANDOM (R) 1/69
IF ((WS/WMAXS) .LT. R**2) GO TO 200
C HAVE INTERMEDIATE MASSES. NOW GENERATE VERTEX DECAYS.
300 DO 303 J=1,NM1 1/69
I = NP1 - J
IM1 = I - 1
U02 = Z(I)/2.
V = (Z(IM1) - X(I))*(Z(IM1) + X(I))/(2.*Z(I))
E10 = U02 + V
E20 = U02 - V
P120 = SQRT((E10 - Z(IM1))*(E10 + Z(IM1)))
IF (J.NE.1 .OR. ITYP.EQ.0) GO TO 301 1/69
CALL RANDIS (ITYP, A, CPOL, E, TMASS, BMASS, X(N), Z(N-1)) 1/69
GO TO 302 1/69
301 CALL RANDOM (CPOL) 1/69
302 CALL RANDOM (R) 1/69
AZM = 3.141593*R 1/69
P120XY = P120*SQRT((1.0-CPOL)*(1.0+CPOL))
P10(1) = P120XY*COS(AZM)
P10(2) = P120XY*SIN(AZM)
P10(3) = P120*CPOL
BETA(1) = Q(I,1)/Q(I,4)
BETA(2) = Q(I,2)/Q(I,4)
BETA(3) = Q(I,3)/Q(I,4)
GAMMA = Q(I,4)/Z(I)
DELTA = (GAMMA*GAMMA)/(1.0+GAMMA)
TLOR(1,1) = 1.0 + DELTA*BETA(1)*BETA(1)
TLOR(2,2) = 1.0 + DELTA*BETA(2)*BETA(2)
TLOR(3,3) = 1.0 + DELTA*BETA(3)*BETA(3)
TLOR(1,2) = DELTA*BETA(1)*BETA(2)
TLOR(1,3) = DELTA*BETA(1)*BETA(3)
TLOR(2,3) = DELTA*BETA(2)*BETA(3)
TLOR(1,4) = GAMMA*BETA(1)
TLOR(2,4) = GAMMA*BETA(2)
TLOR(3,4) = GAMMA*BETA(3)
TLOR(4,4) = GAMMA
TLOR(2,1) = TLOR(1,2)
TLOR(3,1) = TLOR(1,3)
TLOR(4,1) = TLOR(1,4)
TLOR(3,2) = TLOR(2,3)
TLOR(4,2) = TLOR(2,4)
TLOR(4,3) = TLOR(3,4)
DO 303 K = 1,4 1/69
PJ3 = TLOR(K,1)*P10(1) + TLOR(K,2)*P10(2) + TLOR(K,3)*P10(3) 1/69
Q(IM1,K) = +PJ3 + TLOR(K,4)*E10 1/69
303 P(I,K) = -PJ3 + TLOR(K,4)*E20 1/69
DO 304 I=1,4
304 P(1,I) = Q(1,I) 1/69
400 RETURN 1/69
END