      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
