Trailing-Edge
-
PDP-10 Archives
-
-
There are no other files named in the archive.
SUBROUTINE FERMI (PFERMI, BENUC, RANP, RDCS )
CFERMI SUBROUTINE TO OBTAIN FERMI MOMENTUM DISTRIBUTION AND DIRECTION 0020
C BENUC IS 2.0 * (PROTON MASS) * (NUCLEON BINDING ENERGY) IN UNITS 0030
C OF (ENERGY)**2 0040
C PFERMI IS THE FERMI MOMENTUM INPUT AND RANP IS THE RETURNED MOMENT 0050
C RANP IS THE SET OF RETURNED RANDUM DIRECTION COSINES 0060
IF (PFERMI) 7, 5, 7 0070
5 RANP = 0.0 0080
GO TO 500 0090
7 IF (BENUC) 8, 5, 8 0100
8 PMAX = 2.0 * PFERMI 0110
WMAX = 0.7 * (PFERMI)**2 0120
10 CONTINUE 0130
CALL RANDOM (RAN)
P = ABS(PMAX * RAN) 0150
W = P**2 / (1.0 +EXP((P**2 - PFERMI**2) / BENUC )) 0160
CALL RANDOM (RAN)
IF ( W - ABS(RAN * WMAX) ) 10, 20, 20 0180
20 RANP = P 0190
CALL RANDCS(RDCS) 0200
500 RETURN 0210
END 0220
SUBROUTINE RANDCS(RDCS)
CRANDCS SUBROUTINE TO GENERATE A SET OF RANDOM DIRECTION COSINES 0020
DIMENSION RDCS(3) 0030
CALL RANDOM (RAN)
RDCS(3) = RAN 0060
SINA = SQRT( 1.0 - RAN**2 ) 0070
CALL RANDOM (RAN)
RDCS(2) = SIN( 3.14159 * RAN) * SINA 0090
RDCS(1) = COS( 3.14159 * RAN) * SINA 0100
RETURN 0110
END 0120
SUBROUTINE RANDIS (ITYP, A, RCOS, E, TMASS, BMASS, X, Z)
C VERSION OF 1/69--COMPATIBLE WITH INVARIANT RDECAY 1/69
C NOTE--RCOS IS MINUS COS THETA FOR THE FIRST PARTICLE 1/69
C IMPROVEMENTS OF 7/69 AVOID THE POSSIBILITY OF OVERFLOW 7/69
DIMENSION A(10) 1/69
PS(A,B,C) = ((A+B+C)*(A-B-C)*(A+B-C)*(A-B+C))/(A*A)
C 1/69
GO TO (300, 100, 200), ITYP 1/69
C
C UNIFORM ANG DIS BETWEEN LIMITS 1. TO A(10) OR A(10) TO -1.
C
100 SIGN=A(10)/ABS(A(10)) 1/69
DEL=SIGN-A(10) 1/69
CALL RANDOM(RCOS) 6/68
RCOS=DEL*ABS(RCOS)-SIGN 1/69
GO TO 400 1/69
C
C COSINE SERIES DISTRIBUTION
C
200 NARG=-11 1/69
SUMMAX =0.0 0140
DO 2016 K = 1, 21 0150
NARG=NARG+1 0160
ARG=FLOAT(NARG)/10. 0161
SUM = A(1) 0170
DO 2014 L = 2, 10 0180
LL = L - 1 0190
2014 SUM = SUM + A(L) *(ARG**LL) 0200
IF ( SUMMAX - SUM) 2015, 2016, 2016 0210
2015 SUMMAX = SUM 0220
ARGMAX = ARG 0230
2016 CONTINUE 0240
SUMMAX = 1.1 * SUMMAX 0250
2005 CALL RANDOM( RCOS ) 6/68
FTETA = A(1) 0270
DO 2017 KK = 2, 10 0280
LL = KK - 1 0290
2017 FTETA = FTETA + A(KK) * (-RCOS)**LL 1/69
CALL RANDOM ( RBETA ) 6/68
BETA = SUMMAX * ABS(RBETA) 0320
IF (FTETA - BETA) 2005, 400, 400 1/69
C
C EXPONENTIAL MOMENTUM TRANSFER DISTRIBUTION EXP(A(1)*T) TIMES 1/69
C PHASE SPACE--T IS BETWEEN THE TARGET AND THE FIRST PARTICLE 1/69
C
300 XLG = .5*ABS(A(1))*SQRT(PS(E,TMASS,BMASS)*PS(E,X,Z)) 7/69
CALL RANDOM (RCOS) 1/69
IF (XLG.LT.1.E-4) GO TO 400 7/69
IF (XLG.GT.10.) GO TO 320 7/69
RCOS = ALOG(1.+ABS(RCOS)*(EXP(2.*XLG)-1.))/XLG - 1. 1/69
RCOS = AMIN1(RCOS,1.) 7/69
GO TO 400 7/69
320 RCOS = ALOG(ABS(RCOS))/XLG + 1. 7/69
RCOS = AMAX1(RCOS,-1.) 7/69
400 RETURN 1/69
END
SUBROUTINE EXPDIS (A, Z, E, TMASS, BMASS, X, ZL, N)
C VERSION OF 1/69--COMPATIBLE WITH INVARIANT RDECAY
C THROWS MASS WHEN EXPONENTIAL MOM TRANSFER DIST IS REQUESTED
C IMPROVEMENTS OF 7/69 AVOID THE POSSIBILITY OF OVERFLOW 7/69
C DIMENSION F(NI+1) WHERE NI IS SOME POWER OF TWO
DIMENSION F(257), S(6)
DATA NI/256/, S/6*0./
PS(A,B,C) = ((A+B+C)*(A-B-C)*(A+B-C)*(A-B+C))/(A*A)
C
IF (A.EQ.S(1).AND.E.EQ.S(2).AND.TMASS.EQ.S(3).AND.BMASS.EQ.S(4)
1 .AND.X.EQ.S(5).AND.ZL.EQ.S(6).AND.N.EQ.NSV) GO TO 40
S(1) = A
S(2) = E
S(3) = TMASS
S(4) = BMASS
S(5) = X
S(6) = ZL
NSV = N
C
C INITIALIZE--NUMERICALLY INTEGRATE DISTRIBUTION FUNCTION
20 NM3 = N-3
PTS = .25*PS(E,TMASS,BMASS)
C = 2.*ABS(A)*SQRT(PTS) 7/69
B = -2.*ABS(A)*SQRT(PTS+TMASS**2) 7/69
XS = X**2 7/69
DZ = (E-X-ZL)/NI
Z0 = ZL-.5*DZ
F(1) = 0.
DO 22 I=1,NI
ZZ = Z0+I*DZ
PRS = .25*PS(E,X,ZZ) 7/69
CP = C*SQRT(PRS) 7/69
BE = B*SQRT(PRS+XS) 7/69
IF (CP.GT.10.) GO TO 21 7/69
XX = EXP(BE)*SINH(CP) 7/69
GO TO 22 7/69
21 XX = .5*EXP(BE+CP) 7/69
22 F(I+1) = F(I)+(XX*(ZZ-ZL)**NM3)/CP 7/69
C NORMALIZE THE INTEGRATION
DO 24 I=2,NI
24 F(I) = F(I)/F(NI+1)
F(NI+1) = 1.
C
C INTERPOLATE Z FROM INTEGRATED DISTRIBUTION FUNCTION
40 CALL RANDOM (R1)
R1 = ABS(R1)
IL = 1
IU = NI+1
42 IM = (IL+IU)/2
IF(F(IM)-R1) 44, 46, 48
44 IL = IM
GO TO 50
46 IL = IM
GO TO 52
48 IU = IM
50 IF ((IU-IL).NE.1) GO TO 42
52 Z = ZL+DZ*(IL-1.+(R1-F(IL))/(F(IU)-F(IL)))
RETURN
END
SUBROUTINE RANLTH (RADL,RHO,RANL)
CRANLTH SUBROUTINE RANLTH(RADL,RHO,RANL)*1 0020
C RADL IS RADIATION LENGTH 0030
C RHO IS DENSITY IF RADL IS GIVEN IN CM 0040
C RANL IS RANDOM LENGTH IN CENTIMETERS*RETURNED 0050
CALL RANDOM (TEMP)
ZETON = ABS(TEMP) 0070
IF (ZETON) 5,5,10 0080
5 RANL=100.0*RADL 0090
GO TO 15 0100
10 RANL=-(RADL/RHO)*ALOG(ZETON) 0110
15 CONTINUE 0120
RETURN 0130
END 0140
SUBROUTINE RANRES (CV,SIGMA,RES)
REAL LAMBDA
DIMENSION ARG(25), LAMBDA(25)
DATA ARG/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,
1 0.65,0.7,0.75,0.8,0.85,0.88,0.91,0.93,0.95,0.97,0.98,0.99/
DATA LAMBDA/0.,0.0787,0.15838,0.2401,0.3249,0.4142,0.5095,0.6128,
1 0.7265,0.8541,1.,1.1708,1.3764,1.6319,1.963,2.414,3.078,4.198
2 ,5.242,7.026,9.058,12.706,21.20,31.82,63.66/
C
IF(SIGMA) 10,5,10
5 RES=CV
GO TO 20
10 CALL RANDOM(RAN)
S=SIGN(1.,RAN)
F = FINT(1, ABS(RAN),25,ARG,LAMBDA)
RES = (S*F*SIGMA)/2. + CV
20 RETURN
END
SUBROUTINE RGAUSS(CV,SIGMA,RG)
DIMENSION GAUSS(28), STDEV(28)
DATA GAUSS/0.5,0.53983,0.57926,0.61791,0.65542,0.69146,0.72575,
1 0.75804,0.78814,0.81594,0.84134,0.86433,0.88493,0.90320,
2 0.91924,0.93319,0.94520,0.95543,0.96407,0.97128,0.97725,
3 0.98610,0.99180,0.99534,0.99744,0.99865,0.99977,0.99997/
DATA STDEV/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.1,1.2,1.3,
1 1.4,1.5,1.6,1.7,1.8,1.9,2.,2.2,2.4,2.6,2.8,3.,3.5,4./
C
CALL RANDOM(RAN)
S=SIGN(1.,RAN)
RAN = ABS(RAN)/2. + 0.5
F = FINT(1,RAN,28,GAUSS,STDEV)
RG = S*F*SIGMA + CV
RETURN
END
SUBROUTINE ITRNDM(I,J)
DATA K /999/
J=K
RETURN
END
FUNCTION FINT(NARG,ARG,NENT,ENT,TABLE)
CERN TC LIBRARY -- GENERAL SECTION -- 9-AUG-65 FINT PAGE 1
DIMENSION ARG(5),NENT(5),ENT(10),TABLE(10)
DIMENSION D(5),NCOMB(5),IENT(5)
KD=1
M=1
JA=1
DO 5 I=1,NARG
NCOMB(I)=1
JB=JA-1+NENT(I)
DO 2 J=JA,JB
IF (ARG(I).LE.ENT(J)) GO TO 3
2 CONTINUE
J=JB
3 IF (J.NE.JA) GO TO 4
J=J+1
4 JR=J-1
D(I)=(ENT(J)-ARG(I))/(ENT(J)-ENT(JR))
IENT(I)=J-JA
KD=KD+IENT(I)*M
M=M*NENT(I)
5 JA=JB+1
FINT=0.
10 FAC=1.
IADR=KD
IFADR=1
DO 15 I=1,NARG
IF (NCOMB(I).EQ.0) GO TO 12
FAC=FAC*(1.-D(I))
GO TO 15
12 FAC=FAC*D(I)
IADR=IADR-IFADR
15 IFADR=IFADR*NENT(I)
FINT=FINT+FAC*TABLE(IADR)
IL=NARG
40 IF (NCOMB(IL).EQ.0) GO TO 80
NCOMB(IL)=0
IF (IL.EQ.NARG) GO TO 10
IL=IL+1
DO 50 K=IL,NARG
50 NCOMB(K)=1
GO TO 10
80 IL=IL-1
IF(IL.NE.0) GO TO 40
RETURN
END 15