C----------------------------------------------------------------------------------------
C	Description:  A new BPC is simulated with EGS4. The BPC is of a 
C       -----------   W-SCI-type. (W: 3.5mm, SCI: 2.6mm, x= 120mm, y= 128mm).
C	              It has 26 layers of W. The following energies 
C                     are stored in a ntupel: energy deposit in each 
C                     layer, total energy deposit
C                     in SCI-layers, total energy deposit in W-layers,
C                     total energy deposit in the BPC, 
C                     sampling fraction,
C                     energy deposit behind the last
C                     W-layer (check for long. leakage) and energy 
C                     deposit in the transverse space
C                     (check for transverse leakage) (for each event as well as 
C           	      for N events).
C
C
C	Author     : Bernd Surrow
C	--------------------------
C       ADDITIONS  : Ulrich Fricke
C
C	Revision History :
C	------------------
C	Version    Initials	 Date	       Brief Description
C	-------    --------	 ----	       -----------------
C       1.0        B.S.          10/12/94      Simulation of new BPC with EGS4 
C	2.0        B.S.          02/20/96      Include W-bar layers
C       3.0        U.F.          20/07/96      Include search for backscattered particles
C                                              into last plane of Si-Tracer
C       4.0        U.F.          16/08/96      Include absorber between BPC & Tracer 
C
C-------------------------------------------------------------------------------------------
C
C      REGION SCHEMATICS: NUMBERS INDICATE REGION NUMBERS. SYMBOLS IN BRACKETS THE ELEMENT,
C                         COMPOSITE OR MIXTURE IN THIS REGION
C
C-------------------------------------------------------------------------------------------
C
C          TRACER --><-ABSORBER-><--     AIR     --><--      BPC     --><--  AIR  -->
C
C                                              56 (AIR)                         Y
C                  --------------- - - - - - - - - - -------------------         ^
C                 |  |           |                  |   |   |      |    |        |
C                 |  |   61      |                  |   |   | ...  |    |        |
C                 |  |  (AL)     |                  | 1 |  2| ...  | 52 |        |
C   -->       63   - |           |      54 (AIR)    |   |   |      |    |  58    |
C  INCIDENT  (AIR) ^ |           |                  |   |   |      |    | (AIR)  |
C  PARTICLE        | |           |                  |-------------------          --------> Z
C   (PHOTON,       | |           |                  |       60 (W)      |
C   ELECTRON)      |  -----------  - - - - - - - - - - -------------------
C                  |                           57 (AIR)
C                  | 
C 
C             LAST TRACER PLANE: 62 (SI)
C
C                    <--------- SIDISTANCE --------->
C-------------------------------------------------------------------------------------------
C
C                           56 (AIR)                                      Y
C
C                --------------------------------                         ^
C                |        |                     |                         |
C                |        |                     |                         |
C                |        |                     |                         |
C                |        |                     |                         |
C  55 (AIR)      | 59 (W) |      1-52           |   53 (AIR)              |
C                |        |  (W,SZINTILLATOR)   |             X  <--------
C                |        |                     |
C                |------- |---------------------|
C                         |      60 (W)         |
C                         -----------------------
C
C                           57 (AIR)
C
C-------------------------------------------------------------------------------------------
C
C DIRECTIONS:   1  +X     CENTER OF SYSTEM (X,Y,Z)=(0,0,0)
C               2  -X
C               3  +Y     IS IN THE MIDDLE (X=Y=0) OF THE FIRST W-LAYER FRONTFACE OF THE BPC
C               4  -Y
C               5  +Z     REGION NUMBER 1
C               6  -Z
C-------------------------------------------------------------------------------------------
C IMPORTANT **** IMPORTANT **** IMPORTANT **** IMPORTANT *** IMPORTANT **** IMPORTANT 
C IF THE FORTRAN VERSION IS USED THE MAXIMUM NUMBER OF MEDIA AND REGIONS IS SET TO 9 AND 
C 70. ALL COMMON BLOCKS ARE DEFINED USING THIS VALUES. IF MORE ARE REQUIRED, ALL
C BLOCKS HAVE TO BE CHANGED. BY SETTING NMED AND NREG IN THE PROGRAMM ONE DEFINES THE
C NUMBERS OF MEDIA AND REGIONS WHICH ARE ACTUALLY USED. 
C IF ONE USES THE MORTRAN CODE THE MAXIMUM NUMBERS ARE DEFINED BY SETTING 
C $MXMED AND $MXREG FOR ALL COMMOM BLOCKS.
C IN BOTH CASES NREG AND NMED DEFINE THE USED NUMBERS OF REGIONS AND MEDIA
C-------------------------------------------------------------------------------------------

      integer IQUEST
      COMMON/QUEST/IQUEST(100)

      integer IDNEVT
      real VIDN1,VIDN2,VIDN3,X(512),LIMIT
      COMMON/PAWIDN/IDNEVT,VIDN1,VIDN2,VIDN3,X(512)

      INTEGER NPAWC,icycle,ISTAT,ERROR,LEN
      PARAMETER (NPAWC=500000)
      COMMON /PAWC/HMEMOR(NPAWC)

      REAL*4 ESCI,SCI1,SCI2,SCI3,SCI4,SCI5,SCI6,SCI7,SCI8,SCI9,SCI10,SCI
     *11, SCI12,SCI13,SCI14,SCI15, SCI16,SCI17,SCI18,SCI19,SCI20,SCI21,S
     *CI22,SCI23,SCI24,SCI25,SCI26

      REAL*4 ex1,ex2,ex3,ex4,ex5,ex6,ex7,ex8,ex9,ex10,ex11,ex12, ex13,ex
     *14,ex15,ey1,ey2,ey3,ey4,ey5,ey6,ey7,ey8,ey9,ey10,ey11,ey12,ey13, e
     *y14,ey15,ey16

      CHARACTER*500 SSCI
      REAL*4 EW,W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14, W15,W16,
     *W17,W18,W19,W20,W21,W22,W23,W24,W25,W26
      CHARACTER*500 SW
      CHARACTER*500 FIN

      REAL*4 ELONG,ETRANS,R,ETOT,xtrue,ytrue,ztrue,utrue,vtrue,wtrue,Etr
     *ue, xs,ys,zs,us,vs,ws,es,qs,rs
      Real Nph,Eph,Nel,Eel
      CHARACTER*500 SADD

      COMMON/BOUNDS/ECUT(60),PCUT(60),VACDST

      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP

      COMMON/MEDIA/ RLC(3),RLDU(3),RHO(3),MSGE(3),MGE(3),MSEKE(3),MEKE(3
     *),MLEKE(3),MCMFP(3),MRANGE(3),IRAYLM(3),NMED

      COMMON/MEDIAC/MEDIA(24,3)
      CHARACTER*4 MEDIA

      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(60),RHOR(60),IRAYLR(60)

      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)

      COMMON/SCI/ESCI,SCI1,SCI2,SCI3,SCI4,SCI5, SCI6,SCI7,SCI8,SCI9,SCI1
     *0, SCI11,SCI12,SCI13,SCI14,SCI15,SCI16,SCI17,SCI18,SCI19,SCI20,SCI
     *21,SCI22, SCI23,SCI24,SCI25,SCI26

      COMMON/W/EW,W1,W2,W3,W4,W5,W6,W7,W8, W9,W10,W11,W12,W13,W14,W15,W1
     *6,W17,W18,W19,W20,W21,W22,W23,W24,W25,W26

      COMMON/ENERGY/EFILL(60),ESUM(60), EFINN(31),EFIN1(15),EFIN2(16),EF
     *(40)

      COMMON/BOXES/BOXPOS(53),BOXSIZ(54,3)

      COMMON/ADD/ELONG,ETRANS,R,ETOT,xtrue,ytrue,ztrue, utrue,vtrue,wtru
     *e,Etrue,xs,ys,zs,us,vs,ws,es,qs,rs,Nph,Eph,Nel,Eel

      COMMON/FE/ex1,ex2,ex3,ex4,ex5,ex6,ex7,ex8,ex9,ex10,ex11,ex12, ex13
     *,ex14,ex15,ey1,ey2,ey3,ey4,ey5,ey6,ey7,ey8,ey9,ey10,ey11,ey12,ey13
     *, ey14,ey15,ey16

      COMMON/FILL/SH(10)

      DOUBLE PRECISION EF,EFILL,ESUM,ESCIF,ESCIN,EWF,EWN,EFIN1,EFIN2,EFI
     *NN, ETOTF,ETOTN,ETRANSN,ELONGN,RF,RN,TOTKE

	integer iseed, v, irand
	double precision xv, rand

      CHARACTER*4 TEMP(24,3)/'A','I','R',' ','A','T',' ','N','T','P',14*
     *' ','P','O','L','Y','S','T','Y','R','E','N','E',13*' ','W',23*' '/

      CALL HLIMIT(NPAWC)
      IQUEST(10) = 65000

      SSCI='ESCI,SCI1,SCI2,SCI3,SCI4,SC5,SCI6,SCI7,SCI8,SCI9,SCI10,SCI11
     *,             SCI12,SCI13,SCI14,SCI15,SCI16,SCI17,SCI18,SCI19,SCI2
     *0,SCI21,SCI22,SCI23,        SCI24,SCI25,SCI26'
      SW='EW,W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14,
     *              W15,W16,W17,W18,W19,W20,W21,W22,W23,W24,W25,W26'
      SADD='ELONG,ETRANS,R,ETOT,xtrue,ytrue,ztrue,utrue,vtrue,
     *              wtrue,Etrue,xs,ys,zs,us,vs,ws,es,qs,rs,Nph,Eph,Nel,E
     *el'
      FIN='ex1,ex2,ex3,ex4,ex5,ex6,ex7,ex8,ex9,ex10,ex11,ex12,
     *              ex13,ex14,ex15,ey1,ey2,ey3,ey4,ey5,ey6,ey7,ey8,ey9,e
     *y10,ey11,ey12,ey13,         ey14,ey15,ey16'

      CPU0=ETIME(TARRAY)
      NREG=60
      NMED=3
        DO 1011 I=1,NREG
        EFILL(I)=0.D0
1011  CONTINUE
1012  CONTINUE
        DO 1021 I=1,NREG
        ESUM(I)=0.D0
1021  CONTINUE
1022  CONTINUE
        DO 1031 I=1,15
        EFIN1(I)=0.D0
1031  CONTINUE
1032  CONTINUE
        DO 1041 I=1,16
        EFIN2(I)=0.D0
1041  CONTINUE
1042  CONTINUE
        DO 1051 I=1,31
        EFINN(I)=0.D0
1051  CONTINUE
1052  CONTINUE
        DO 1061 I=1,40
        EF(I)=0.D0
1061  CONTINUE
1062  CONTINUE
        DO 1071 I=1,40
        EF(I)=0.D0
1071  CONTINUE
1072  CONTINUE
        DO 1081 I=1,10
        SH(I)=0.0
1081  CONTINUE
1082  CONTINUE
        DO 1091 J=1,NMED
          DO 1101 I=1,24
          MEDIA(I,J)=TEMP(I,J)
1101    CONTINUE
1102    CONTINUE
1091  CONTINUE
1092  CONTINUE
      MED(1)=3
      MED(3)=3
      MED(5)=3
      MED(7)=3
      MED(9)=3
      MED(11)=3
      MED(13)=3
      MED(15)=3
      MED(17)=3
      MED(19)=3
      MED(21)=3
      MED(23)=3
      MED(25)=3
      MED(27)=3
      MED(29)=3
      MED(31)=3
      MED(33)=3
      MED(35)=3
      MED(37)=3
      MED(39)=3
      MED(41)=3
      MED(43)=3
      MED(45)=3
      MED(47)=3
      MED(49)=3
      MED(51)=3
      MED(59)=3
      MED(60)=3
      MED(2)=2
      MED(4)=2
      MED(6)=2
      MED(8)=2
      MED(10)=2
      MED(12)=2
      MED(14)=2
      MED(16)=2
      MED(18)=2
      MED(20)=2
      MED(22)=2
      MED(24)=2
      MED(26)=2
      MED(28)=2
      MED(30)=2
      MED(32)=2
      MED(34)=2
      MED(36)=2
      MED(38)=2
      MED(40)=2
      MED(42)=2
      MED(44)=2
      MED(46)=2
      MED(48)=2
      MED(50)=2
      MED(52)=2
      MED(53)=1
      MED(54)=1
      MED(55)=1
      MED(56)=1
      MED(57)=1
      MED(58)=1
        DO 1111 I=1,NREG
        ECUT(I)=0.01
        PCUT(I)=0.01
1111  CONTINUE
1112  CONTINUE
      CALL HATCH
        DO 1121 J=1,53
        BOXPOS(J)=0.0
1121  CONTINUE
1122  CONTINUE
        DO 1131 J=2,53
          DO 1141 I=1,(J-1)
          BOXPOS(J)=BOXPOS(J)+BOXSIZ(I,3)
1141    CONTINUE
1142    CONTINUE
1131  CONTINUE
1132  CONTINUE

	iseed=43980342
	call srand(iseed)

      CALL HBOOK1(10,'Number of Photons per event',50,0.,50.,0.)
      CALL HBOOK1(20,'Number of Electrons per event',50,0.,50.,0.)
      CALL HBOOK1(30,'Energy (MeV) of Photons per event',200,0.,50.,0.)
      CALL HBOOK1(40,'Energy (MeV) of Electrons per event',200,0.,50.,0.
     *)

      OPEN(UNIT=12,FILE='egs.dat',STATUS='NEW')

      CALL HROPEN(118,'sel','egs4.rz','NQ',4096,ISTAT)
      IF (ISTAT.NE.0) THEN
         write (*,*) "Error opening file : >>",'egs4.rz',"<<"
         STOP
      ENDIF 
      CALL HBNT(1,'sel',' ')
      CALL HCDIR('//sel',' ')

C      CALL HROPEN(30,'MYFILE','egs4.nt','N',4096,ISTAT)
C      CALL HBSET('BSIZE',64,ERROR)
C      CALL HBNT(1,'egs4.nt','D')

      CALL HBNAME(1,'SCI',ESCI,SSCI)
      CALL HBNAME(1,'FE',ex1,FIN)
      CALL HBNAME(1,'W',EW,SW)
      CALL HBNAME(1,'ADD',ELONG,SADD)
      IQIN=-1
      EIN=1000
      XIN=0.0
      YIN=0.0
      ZIN=0.0
      UIN=0.0
      VIN=0.0
      WIN=1.0
      IRIN=1
      WTIN=1.0
      TIMUSR=10000000000
      NCASES=10
      LIMIT=0
      CPUT1=ETIME(TARRAY)
        DO 1151 I=1,NCASES
0000    CONTINUE

C...Scatter electron arround a certain region:

  	xv = rand( )

        XIN=-0.8+1.6*xv

        CPUT2=ETIME(TARRAY)
        TIMCPU=(CPUT2-CPUT0)*1.0
        AVRCPU=((CPUT2-CPUT1)*1.0)/FLOAT(I)
        WRITE(6,1160)I,TIMCPU
1160    FORMAT(T5,' START HISTORY:',I4,T30,'TIME CPU: ',F10.3)
        IF (((TIMCPU+AVRCPU).GE.TIMUSR)) THEN
          WRITE(6,1170)
1170      FORMAT(/,T15,'SIMULATION STOPS DUE TO THE TIME LIMIT')
          LIMIT=1
          WRITE(6,1180)I-1
1180      FORMAT(T15,'NUMBER OF EVENTS CALCULATED ',I4)
          GO TO1152
        END IF
        xtrue=XIN
        ytrue=YIN
        ztrue=ZIN
        utrue=UIN
        vtrue=VIN
        wtrue=WIN
        Etrue=EIN
        CALL SHOWER(IQIN,EIN,XIN,YIN,ZIN,UIN,VIN,WIN,IRIN,WTIN)
        EWF=0.D0
          DO 1191 J=0,25
          EWF=EWF+EFILL(2*J+1)
1191    CONTINUE
1192    CONTINUE
        EW=EWF
        ESCIF=0.D0
          DO 1201 K=1,25
          ESCIF=ESCIF+EFILL(2*K)
1201    CONTINUE
1202    CONTINUE
        ESCI=ESCIF
        ETOTF=0.D0
          DO 1211 L=1,51
          ETOTF=ETOTF+EFILL(L)
1211    CONTINUE
1212    CONTINUE
        ETOT=ETOTF
        RF=0.D0
        RF=ESCIF/(ESCIF+EWF)
        R=RF
        Nph=SH(1)
        Eph=EF(32)
        Nel=SH(2)
        Eel=EF(33)
        CALL HF1(10,Nph,1.)
        CALL HF1(20,Nel,1.)
        CALL HF1(30,Eph,1.)
        CALL HF1(40,Eel,1.)
        W1=EFILL(1)
        W2=EFILL(3)
        W3=EFILL(5)
        W4=EFILL(7)
        W5=EFILL(9)
        W6=EFILL(11)
        W7=EFILL(13)
        W8=EFILL(15)
        W9=EFILL(17)
        W10=EFILL(19)
        W11=EFILL(21)
        W12=EFILL(23)
        W13=EFILL(25)
        W14=EFILL(27)
        W15=EFILL(29)
        W16=EFILL(31)
        W17=EFILL(33)
        W18=EFILL(35)
        W19=EFILL(37)
        W20=EFILL(39)
        W21=EFILL(41)
        W22=EFILL(43)
        W23=EFILL(45)
        W24=EFILL(47)
        W25=EFILL(49)
        W26=EFILL(51)
        SCI1=EFILL(2)
        SCI2=EFILL(4)
        SCI3=EFILL(6)
        SCI4=EFILL(8)
        SCI5=EFILL(10)
        SCI6=EFILL(12)
        SCI7=EFILL(14)
        SCI8=EFILL(16)
        SCI9=EFILL(18)
        SCI10=EFILL(20)
        SCI11=EFILL(22)
        SCI12=EFILL(24)
        SCI13=EFILL(26)
        SCI14=EFILL(28)
        SCI15=EFILL(30)
        SCI16=EFILL(32)
        SCI17=EFILL(34)
        SCI18=EFILL(36)
        SCI19=EFILL(38)
        SCI20=EFILL(40)
        SCI21=EFILL(42)
        SCI22=EFILL(44)
        SCI23=EFILL(46)
        SCI24=EFILL(48)
        SCI25=EFILL(50)
        SCI26=EFILL(52)
        ex1=EF(1)
        ex2=EF(2)
        ex3=EF(3)
        ex4=EF(4)
        ex5=EF(5)
        ex6=EF(6)
        ex7=EF(7)
        ex8=EF(8)
        ex9=EF(9)
        ex10=EF(10)
        ex11=EF(11)
        ex12=EF(12)
        ex13=EF(13)
        ex14=EF(14)
        ex15=EF(15)
        ey1=EF(16)
        ey2=EF(17)
        ey3=EF(18)
        ey4=EF(19)
        ey5=EF(20)
        ey6=EF(21)
        ey7=EF(22)
        ey8=EF(23)
        ey9=EF(24)
        ey10=EF(25)
        ey11=EF(26)
        ey12=EF(27)
        ey13=EF(28)
        ey14=EF(29)
        ey15=EF(30)
        ey16=EF(31)
        ELONG=EFILL(58)
        ETRANS=EFILL(53)+EFILL(55)+EFILL(56)+EFILL(57)
        CALL HFNT(1)
          DO 1221 K=1,NREG
          EFILL(K)=0.D0
1221    CONTINUE
1222    CONTINUE
          DO 1231 K=1,40
          EF(K)=0.D0
1231    CONTINUE
1232    CONTINUE
          DO 1241 K=1,10
          SH(K)=0.0
1241    CONTINUE
1242    CONTINUE
          DO 1251 K=1,15
          EFIN1(K)=0.D0
1251    CONTINUE
1252    CONTINUE
          DO 1261 K=1,16
          EFIN2(K)=0.D0
1261    CONTINUE
1262    CONTINUE
1151  CONTINUE
1152  CONTINUE
      CALL HRPUT(0,'egs.rz','N')
      CPUT2=ETIME(TARRAY)
      IF ((LIMIT.EQ.0)) THEN
        AVRCPU=((CPUT2-CPUT1)*1.0)/FLOAT(NCASES)
      END IF
      TIMCPU=(CPUT2-CPUT0)*1.0
      WRITE(6,1270)TIMCPU
1270  FORMAT(T15,'TOTAL CPU TIME : ',F7.3)
      WRITE(6,1280)AVRCPU
1280  FORMAT(T15,'AVERAGE TIME PER HISTORY: ' F7.3)
      TOTKE=0.D0
      IF ((LIMIT.EQ.0)) THEN
        TOTKE=EIN*NCASES
      ELSE
        TOTKE=EIN*(I-1)
      END IF
      ETOTN=0.D0
      ESCIN=0.D0
      EWN=0.D0
      ELONGN=0.D0
      ETRANSN=0.D0
      RN=0.D0
        DO 1291 I=1,52
        ETOTN=ETOTN+ESUM(I)
1291  CONTINUE
1292  CONTINUE
        DO 1301 I=1,26
        ESCIN=ESCIN+ESUM(2*I)
1301  CONTINUE
1302  CONTINUE
        DO 1311 I=0,25
        EWN=EWN+ESUM(2*I+1)
1311  CONTINUE
1312  CONTINUE
      ELONGN=ESUM(58)
      ETRANSN=ESUM(53)+ESUM(55)+ESUM(56)+ESUM(57)
      RN=ESCIN/(ESCIN+EWN)
      ETOTN=ETOTN/TOTKE
      ESCIN=ESCIN/TOTKE
      EWN=EWN/TOTKE
      ELONGN=ELONGN/TOTKE
      ETRANSN=ETRANSN/TOTKE
      WRITE(6,1320)
1320  FORMAT(/,T5,'ENERGY DEPOSIT MeV',T28,'ENERGY DEPOSIT %', T50,'REGI
     *ON',T60,'MEDIUM')
        DO 1331 L=1,NREG
        WRITE(6,1340)ESUM(L),100*ESUM(L)/TOTKE,L,MED(L)
1340    FORMAT(T10,F15.3,T30,F10.3,T52,I2,T61,I2)
1331  CONTINUE
1332  CONTINUE
        DO 1351 J=1,31
        WRITE(6,1360)EFINN(J),100*EFINN(J)/TOTKE,J,MED(J)
1360    FORMAT(T10,F15.3,T30,F10.3,T52,I2,T61,I2)
1351  CONTINUE
1352  CONTINUE
      ETOTN=100*ETOTN
      WRITE(6,1370)
1370  FORMAT(/,T20,'Fraction of energy deposited in BPC in %:')
      WRITE(6,1380)ETOTN
1380  FORMAT(T20,G15.10)
      WRITE(6,1390)
1390  FORMAT(/,T20,'Longitudinal energy leakage in %:')
      WRITE(6,1400)100*ELONGN
1400  FORMAT(T20,G15.10)
      WRITE(6,1410)
1410  FORMAT(/,T20,'Transversal energy leakage in %:')
      WRITE(6,1420)100*ETRANSN
1420  FORMAT(T20,G15.10)
      WRITE(6,1430)
1430  FORMAT(/,T20,'Sampling Fraction in %:')
      WRITE(6,1440)100*RN
1440  FORMAT(T20,G15.10)
      WRITE(12,*) ''
      WRITE(12,*) 'Fraction of energy deposited in BPC in %:'
      WRITE(12,83) ETOTN

      WRITE(12,*) ''
      WRITE(12,*) 'Longitudinal energy leakage in %:'
      WRITE(12,83) 100*ELONGN

      WRITE(12,*) ''
      WRITE(12,*) 'Transversal energy leakage in %:'
      WRITE(12,83) 100*ETRANSN

      WRITE(12,*) ''
      WRITE(12,*) 'Sampling Fraction in %:'
      WRITE(12,83) 100*RN

83    FORMAT(20X,G15.13)

      ETOTN=ETOTN*0.01
      WRITE(6,1450)
1450  FORMAT(/,T20,'This quantity should be close to unity:')
      WRITE(6,1460)ETOTN
1460  FORMAT(T20,G15.13)
      WRITE(12,*) ''
      WRITE(12,*) 'This quantity should be close to unity:'
      WRITE(12,84) ETOTN
 84   FORMAT(20X,G15.13)

 85   FORMAT(1X,F15.3,2X,F10.3,10X,I2,10X,2A)
      WRITE(12,*) 'ENERGY DEPOSIT MeV  ENERGY DEPOSIT %'//
     + ' REGION    MEDIUM'

      DO 86 L=1,NREG
      WRITE(12,85) ESUM(L),100*ESUM(L)/TOTKE,L,MED(L)
 86   CONTINUE

      DO 87 L=1,31
      WRITE(12,85) EFINN(L),100*EFINN(L)/TOTKE,L,MED(L)  		
 87   CONTINUE

      CLOSE(12)

      CALL HCDIR('//sel',' ')
      CALL HROUT(1,icycle,' ')
      CALL HREND('sel')
      CLOSE(118)

C      CALL HROUT(1,icycle,' ')
C      CALL HREND('MYFILE')

      WRITE(6,1470)
1470  FORMAT(/,T25,'END OF SIMULATION',//)
      STOP
      END

      SUBROUTINE AUSGAB(IARG)

      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/ENERGY/EFILL(60),ESUM(60), EFINN(31),EFIN1(15),EFIN2(16),EF
     *(40)
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/FILL/SH(10)
      DOUBLE PRECISION EF,EFILL,ESUM,EFIN1,EFIN2,EFINN
      IRK=IR(NP)
        DO 1481 I=0,12
        IF ((IRK.EQ.(2+4*I))) THEN
          IF ((X(NP).LE.-5.2.AND.X(NP).GE.-6.0)) THEN
            EF(1)=EF(1)+EDEP
            EFINN(1)=EFINN(1)+EDEP
          END IF
          IF ((X(NP).LE.-4.4.AND.X(NP).GE.-5.2)) THEN
            EF(2)=EF(2)+EDEP
            EFINN(2)=EFINN(2)+EDEP
          END IF
          IF ((X(NP).LE.-3.6.AND.X(NP).GE.-4.4)) THEN
            EF(3)=EF(3)+EDEP
            EFINN(3)=EFINN(3)+EDEP
          END IF
          IF ((X(NP).LE.-2.8.AND.X(NP).GE.-3.6)) THEN
            EF(4)=EF(4)+EDEP
            EFINN(4)=EFINN(4)+EDEP
          END IF
          IF ((X(NP).LE.-2.0.AND.X(NP).GE.-2.8)) THEN
            EF(5)=EF(5)+EDEP
            EFINN(5)=EFINN(5)+EDEP
          END IF
          IF ((X(NP).LE.-1.2.AND.X(NP).GE.-2.0)) THEN
            EF(6)=EF(6)+EDEP
            EFINN(6)=EFINN(6)+EDEP
          END IF
          IF ((X(NP).LE.-0.4.AND.X(NP).GE.-1.2)) THEN
            EF(7)=EF(7)+EDEP
            EFINN(7)=EFINN(7)+EDEP
          END IF
          IF ((X(NP).LE.0.4.AND.X(NP).GE.-0.4)) THEN
            EF(8)=EF(8)+EDEP
            EFINN(8)=EFINN(8)+EDEP
          END IF
          IF ((X(NP).LE.1.2.AND.X(NP).GE.0.4)) THEN
            EF(9)=EF(9)+EDEP
            EFINN(9)=EFINN(9)+EDEP
          END IF
          IF ((X(NP).LE.2.0.AND.X(NP).GE.1.2)) THEN
            EF(10)=EF(10)+EDEP
            EFINN(10)=EFINN(10)+EDEP
          END IF
          IF ((X(NP).LE.2.8.AND.X(NP).GE.2.0)) THEN
            EF(11)=EF(11)+EDEP
            EFINN(11)=EFINN(11)+EDEP
          END IF
          IF ((X(NP).LE.3.6.AND.X(NP).GE.2.8)) THEN
            EF(12)=EF(12)+EDEP
            EFINN(12)=EFINN(12)+EDEP
          END IF
          IF ((X(NP).LE.4.4.AND.X(NP).GE.3.6)) THEN
            EF(13)=EF(13)+EDEP
            EFINN(13)=EFINN(13)+EDEP
          END IF
          IF ((X(NP).LE.5.2.AND.X(NP).GE.4.4)) THEN
            EF(14)=EF(14)+EDEP
            EFINN(14)=EFINN(14)+EDEP
          END IF
          IF ((X(NP).LE.6.0.AND.X(NP).GE.5.2)) THEN
            EF(15)=EF(15)+EDEP
            EFINN(15)=EFINN(15)+EDEP
          END IF
        END IF
1481  CONTINUE
1482  CONTINUE
        DO 1491 I=0,12
        IF ((IRK.EQ.(4+4*I))) THEN
          IF ((Y(NP).LE.-5.6.AND.Y(NP).GE.-6.4)) THEN
            EF(16)=EF(16)+EDEP
            EFINN(16)=EFINN(16)+EDEP
          END IF
          IF ((Y(NP).LE.-4.8.AND.Y(NP).GE.-5.6)) THEN
            EF(17)=EF(17)+EDEP
            EFINN(17)=EFINN(17)+EDEP
          END IF
          IF ((Y(NP).LE.-4.0.AND.Y(NP).GE.-4.8)) THEN
            EF(18)=EF(18)+EDEP
            EFINN(18)=EFINN(18)+EDEP
          END IF
          IF ((Y(NP).LE.-3.2.AND.Y(NP).GE.-4.0)) THEN
            EF(19)=EF(19)+EDEP
            EFINN(19)=EFINN(19)+EDEP
          END IF
          IF ((Y(NP).LE.-2.4.AND.Y(NP).GE.-3.2)) THEN
            EF(20)=EF(20)+EDEP
            EFINN(20)=EFINN(20)+EDEP
          END IF
          IF ((Y(NP).LE.-1.6.AND.Y(NP).GE.-2.4)) THEN
            EF(21)=EF(21)+EDEP
            EFINN(21)=EFINN(21)+EDEP
          END IF
          IF ((Y(NP).LE.-0.8.AND.Y(NP).GE.-1.6)) THEN
            EF(22)=EF(22)+EDEP
            EFINN(22)=EFINN(22)+EDEP
          END IF
          IF ((Y(NP).LE.0.0.AND.Y(NP).GE.-0.8)) THEN
            EF(23)=EF(23)+EDEP
            EFINN(23)=EFINN(23)+EDEP
          END IF
          IF ((Y(NP).LE.0.8.AND.Y(NP).GE.0.0)) THEN
            EF(24)=EF(24)+EDEP
            EFINN(24)=EFINN(24)+EDEP
          END IF
          IF ((Y(NP).LE.1.6.AND.Y(NP).GE.0.8)) THEN
            EF(25)=EF(25)+EDEP
            EFINN(25)=EFINN(25)+EDEP
          END IF
          IF ((Y(NP).LE.2.4.AND.Y(NP).GE.1.6)) THEN
            EF(26)=EF(26)+EDEP
            EFINN(26)=EFINN(26)+EDEP
          END IF
          IF ((Y(NP).LE.3.2.AND.Y(NP).GE.2.4)) THEN
            EF(27)=EF(27)+EDEP
            EFINN(27)=EFINN(27)+EDEP
          END IF
          IF ((Y(NP).LE.4.0.AND.Y(NP).GE.3.2)) THEN
            EF(28)=EF(28)+EDEP
            EFINN(28)=EFINN(28)+EDEP
          END IF
          IF ((Y(NP).LE.4.8.AND.Y(NP).GE.4.0)) THEN
            EF(29)=EF(29)+EDEP
            EFINN(29)=EFINN(29)+EDEP
          END IF
          IF ((Y(NP).LE.5.6.AND.Y(NP).GE.4.8)) THEN
            EF(30)=EF(30)+EDEP
            EFINN(30)=EFINN(30)+EDEP
          END IF
          IF ((Y(NP).LE.6.4.AND.Y(NP).GE.5.6)) THEN
            Ef(31)=EF(31)+EDEP
            EFINN(31)=EFINN(31)+EDEP
          END IF
        END IF
1491  CONTINUE
1492  CONTINUE
      ESUM(IRK)=ESUM(IRK)+EDEP
      EFILL(IRK)=EFILL(IRK)+EDEP
      IF ((IARG.EQ.3.AND.IR(NP).EQ.54.AND.IQ(NP).EQ.0)) THEN
        SH(1)=SH(1)+1.
        EF(32)=EF(32)+E(NP)
      END IF
      IF ((IARG.EQ.3.AND.IR(NP).EQ.54.AND.IQ(NP).EQ.+1)) THEN
        SH(2)=SH(2)+1.
        EF(33)=EF(33)+E(NP)
      END IF
      RETURN
      END

      SUBROUTINE HOWFAR

      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/ENERGY/EFILL(60),ESUM(60), EFINN(31),EFIN1(15),EFIN2(16),EF
     *(40)
      COMMON/BOXES/BOXPOS(53),BOXSIZ(54,3)
      COMMON/FILL/SH(10)
      DOUBLE PRECISION EF,EFILL,ESUM,EFIN1,EFIN2,EFINN
      U(NP)=U(NP)/SQRT(U(NP)*U(NP)+V(NP)*V(NP)+W(NP)*W(NP))
      V(NP)=V(NP)/SQRT(U(NP)*U(NP)+V(NP)*V(NP)+W(NP)*W(NP))
      W(NP)=W(NP)/SQRT(U(NP)*U(NP)+V(NP)*V(NP)+W(NP)*W(NP))
      IRL=IR(NP)
      IF ((IRL.EQ.53.OR.IRL.EQ.54.OR.IRL.EQ.55.OR.IRL.EQ.56.OR. IRL.EQ.5
     *7.OR.IRL.EQ.58)) THEN
        IDISC=1
        RETURN
      ELSE
        CALL BOX(DSTVAL,DSTNST,IDIR)
        IF ((IRL.EQ.59)) THEN
          IF ((DSTVAL.GT.USTEP)) THEN
            RETURN
          ELSE
            USTEP=DSTVAL
            IF ((IDIR.EQ.2)) THEN
              IRNEW=55
              RETURN
            ELSE IF((IDIR.EQ.3)) THEN
              IRNEW=56
              RETURN
            ELSE IF((IDIR.EQ.4)) THEN
              IRNEW=57
              RETURN
            ELSE IF((IDIR.EQ.5)) THEN
              IRNEW=58
              RETURN
            ELSE IF((IDIR.EQ.6)) THEN
              IRNEW=54
              RETURN
            ELSE IF((IDIR.EQ.1)) THEN
                DO 1501 I=1,52
                IF ((Z(NP)+USTEP*W(NP).GE.BOXPOS(I).
     &          AND.Z(NP)+USTEP*W(NP).LT.BOXPOS(I+1))) THEN
                  IRNEW=I
                  RETURN
                END IF
1501          CONTINUE
1502          CONTINUE
            END IF
          END IF
        END IF
        IF ((IRL.EQ.60)) THEN
          IF ((DSTVAL.GT.USTEP)) THEN
            RETURN
          ELSE
            USTEP=DSTVAL
            IF ((IDIR.EQ.1)) THEN
              IRNEW=53
              RETURN
            ELSE IF((IDIR.EQ.2)) THEN
              IRNEW=55
              RETURN
            ELSE IF((IDIR.EQ.4)) THEN
              IRNEW=57
              RETURN
            ELSE IF((IDIR.EQ.5)) THEN
              IRNEW=58
              RETURN
            ELSE IF((IDIR.EQ.6)) THEN
              IRNEW=54
              RETURN
            ELSE IF((IDIR.EQ.3)) THEN
                DO 1511 I=1,52
                IF ((Z(NP)+USTEP*W(NP).GE.BOXPOS(I).
     &          AND.Z(NP)+USTEP*W(NP).LT.BOXPOS(I+1))) THEN
                  IRNEW=I
                  RETURN
                END IF
1511          CONTINUE
1512          CONTINUE
            END IF
          END IF
        ELSE
          IF ((DSTVAL.GT.USTEP)) THEN
            RETURN
          ELSE
            USTEP=DSTVAL
            IF ((IDIR.EQ.1)) THEN
              IRNEW=53
              RETURN
            ELSE IF((IDIR.EQ.2)) THEN
              IRNEW=59
              RETURN
            ELSE IF((IDIR.EQ.3)) THEN
              IRNEW=56
              RETURN
            ELSE IF((IDIR.EQ.4)) THEN
              IRNEW=60
              RETURN
            ELSE IF((IDIR.EQ.5)) THEN
              IF ((IRL.EQ.52)) THEN
                IRNEW=58
                RETURN
              END IF
              IRNEW=IRL+1
              RETURN
            ELSE IF((IDIR.EQ.6)) THEN
              IF ((IRL.EQ.1)) THEN
                IRNEW=54
                RETURN
              END IF
              IRNEW=IRL-1
              RETURN
            END IF
          END IF
        END IF
      END IF
      RETURN
      END
      SUBROUTINE BOX(DSTVAL,DSTNST,IDIR)
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/BOXES/BOXPOS(53),BOXSIZ(54,3)
      DIS=1.0E04
      IF ((IR(NP).EQ.59)) THEN
        IRTEMP=IR(NP)-6
        DIS1=-X(NP)-6.0
        DIS2=6.5+X(NP)
        DIS3=BOXSIZ(IRTEMP,2)-Y(NP)
        DIS4=2*BOXSIZ(IRTEMP,2)-DIS3
        DIS6=Z(NP)
        DIS5=15.86-Z(NP)
      ELSE IF((IR(NP).EQ.60)) THEN
        IRTEMP=IR(NP)-6
        DIS1=BOXSIZ(IRTEMP,1)-X(NP)
        DIS2=2*BOXSIZ(IRTEMP,1)-DIS1
        DIS3=-Y(NP)-6.4
        DIS4=Y(NP)+6.9
        DIS6=Z(NP)
        DIS5=15.86-Z(NP)
      ELSE
        IRTEMP=IR(NP)
        DIS1=BOXSIZ(IRTEMP,1)-X(NP)
        DIS2=2*BOXSIZ(IRTEMP,1)-DIS1
        DIS3=BOXSIZ(IRTEMP,2)-Y(NP)
        DIS4=2*BOXSIZ(IRTEMP,2)-DIS3
        DIS6=Z(NP)-BOXPOS(IRTEMP)
        DIS5=BOXSIZ(IRTEMP,3)-DIS6
      END IF
      DSTNST=AMIN1(DIS1,DIS2,DIS3,DIS4,DIS5,DIS6)
      IF ((DSTNST.LT.0.0)) THEN
        DSTNST=0.0
      END IF
      DNEAR(NP)=DSTNST
      IF ((U(NP).LT.0.0)) THEN
        DIS=-DIS2/U(NP)
        IDIR=2
      ELSE IF((U(NP).GT.0.0)) THEN
        DIS=DIS1/U(NP)
        IDIR=1
      END IF
      IF ((V(NP).LT.0.0)) THEN
        DIS4=-DIS4/V(NP)
        IF ((DIS4.LT.DIS)) THEN
          DIS=DIS4
          IDIR=4
        END IF
      ELSE IF((V(NP).GT.0.0)) THEN
        DIS3=DIS3/V(NP)
        IF ((DIS3.LT.DIS)) THEN
          DIS=DIS3
          IDIR=3
        END IF
      END IF
      IF ((W(NP).LT.0.0)) THEN
        DIS6=-DIS6/W(NP)
        IF ((DIS6.LT.DIS)) THEN
          DIS=DIS6
          IDIR=6
        END IF
      ELSE IF((W(NP).GT.0.0)) THEN
        DIS5=DIS5/W(NP)
        IF ((DIS5.LT.DIS)) THEN
          DIS=DIS5
          IDIR=5
        END IF
      END IF
      DSTVAL=DIS
      IF ((DSTVAL.LT.0.0)) THEN
        DSTVAL=0.0
      END IF
      RETURN
      END
      BLOCK DATA BOXDATA
      COMMON/BOXES/BOXPOS(53),BOXSIZ(54,3)
      DATA BOXSIZ/ 52*6.0,0.25,6.0, 52*6.4,6.4,0.25, 0.35,0.26,0.35,0.26
     *,0.35,0.26,0.35,0.26,0.35,0.26,0.35, 0.26,0.35,0.26,0.35,0.26,0.35
     *,0.26,0.35,0.26,0.35,0.26,0.35, 0.26,0.35,0.26,0.35,0.26,0.35,0.26
     *,0.35,0.26,0.35,0.26,0.35, 0.26,0.35,0.26,0.35,0.26,0.35,0.26,0.35
     *,0.26,0.35,0.26,0.35, 0.26,0.35,0.26,0.35,0.26,15.86,15.86/
      END
      BLOCK DATA
      COMMON/BOUNDS/ECUT(60),PCUT(60),VACDST
      COMMON/BREMPR/ DL1(6,3),DL2(6,3),DL3(6,3),DL4(6,3),DL5(6,3),DL6(6,
     *3), ALPHI(2,3),BPAR(2,3),DELPOS(2,3), ASYM(3,50,2), WA(3,50),PZ(3,
     *50),ZELEM(3,50),RHOZ(3,50), PWR2I(50), DELCM(3),NNE(3),ZBRANG(3),
     *IBRDST,IPRDST,IBRSPL,NBRSPL,FBRSPL
      CHARACTER*4 ASYM
      INTEGER EPSTFL(3)
      COMMON/ELECIN/ EKELIM, EKE0(3),EKE1(3),CMFP0(3),CMFP1(3),RANGE0(3)
     *,RANGE1(3), XR0(3),TEFF0(3),BLCC(3),XCC(3), PICMP0(1,3),PICMP1(1,3
     *),EICMP0(1,3),EICMP1(1,3), ESIG0(150,3),ESIG1(150,3),PSIG0(150,3),
     *PSIG1(150,3),EDEDX0(150,3),EDEDX1(150,3),PDEDX0(150,3),PDEDX1(150,
     *3),EBR10(150,3),EBR11(150,3),PBR10(150,3),PBR11(150,3),PBR20(150,3
     *),PBR21(150,3),TMXS0(150,3),TMXS1(150,3), CMFPE0(1,3),CMFPE1(1,3),
     *CMFPP0(1,3),CMFPP1(1,3),ERANG0(1,3),ERANG1(1,3),PRANG0(1,3),PRANG1
     *(1,3),CXC2E0(1,3),CXC2E1(1,3),CXC2P0(1,3),CXC2P1(1,3),CLXAE0(1,3),
     *CLXAE1(1,3),CLXAP0(1,3),CLXAP1(1,3), THR0(1,1),THR1(1,1),THR2(1,1)
     *, THRI0(1,1),THRI1(1,1),THRI2(1,1), FSTEP(16),FSQR(16), VERT1(1000
     *),VERT2(100,16), BLC0,BLC1,RTHR0,RTHR1,RTHRI0,RTHRI1, ICOMP, MPEEM
     *(1,3), MSMAP(200), MSTEPS,JRMAX,MXV1, MXV2,NBLC,NRNTH,NRNTHI, IUNR
     *ST(3),EPSTFL,IAPRIM(3)
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/MEDIA/ RLC(3),RLDU(3),RHO(3),MSGE(3),MGE(3),MSEKE(3),MEKE(3
     *),MLEKE(3),MCMFP(3),MRANGE(3),IRAYLM(3),NMED
      COMMON/MEDIAC/MEDIA(24,3)
      CHARACTER*4 MEDIA
      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(60),RHOR(60),IRAYLR(60)
      COMMON/MULTS/ B0G21,B1G21,G210(7),G211(7),G212(7), B0G22,B1G22,G22
     *0(8),G221(8),G222(8), B0G31,B1G31,G310(11),G311(11),G312(11), B0G3
     *2,B1G32,G320(25),G321(25),G322(25), B0BGB,B1BGB,BGB0(8),BGB1(8),BG
     *B2(8), NG21,NG22,NG31,NG32,NBGB
      COMMON/PATHCM/B0PTH,B1PTH,PTH0(6),PTH1(6),PTH2(6),NPTH
      COMMON/PHOTIN/ EBINDA(3), GE0(3),GE1(3), GMFP0(200,3),GMFP1(200,3)
     *,GBR10(200,3),GBR11(200,3),GBR20(200,3),GBR21(200,3), RCO0(3),RCO1
     *(3), RSCT0(100,3),RSCT1(100,3), COHE0(200,3),COHE1(200,3), MPGEM(1
     *,3), NGR(3)
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIIN/SINC0,SINC1,SIN0(1002),SIN1(1002)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/EDGE/IEDGFL(60),EKALPH(3),EKBETA(3),BKR1(3),BKR2(3)
      COMMON/USERSC/ SMAX,SMAXIR(60),ESTEPE,ESTEPR(60),ESAVE(60), IEDEP,
     *NOMSCT(60),NOPLC(60)
      COMMON/USERVR/ CEXPTR,GWAIT,IFORCE,NFMIN,NFMAX,NFTIME,ISOURC,IFPB,
     *IQINC,MONOEN
      COMMON/USERXT/IPHTER(60)
      COMMON/RANDOM/IXX
      CHARACTER*4 MEDIA1(24)
      EQUIVALENCE(MEDIA1(1),MEDIA(1,1))
      DATA ECUT/60*0./,PCUT/60*0./,VACDST/1.E8/
      DATA EKELIM/0./,ICOMP/1/
      DATA IAUSFL/5*1,20*0/,RHOF/1.0/
      DATA NMED /1/, MEDIA1/'N','A','I',' ',' ',' ',' ',' ',' ',' ',' ',
     *' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '/
      DATA IRAYLM/3*0/
      DATA KMPI/12/,KMPO/8/,DUNIT/1./,NOSCAT/0/
      DATA MED/60*1/,RHOR/60*0./,IRAYLR/60*0/
      DATA NG21/ 7/,B0G21/ 2.0000E+00/,B1G21/ 5.0000E+00/
      DATA G210( 1),G211( 1),G212( 1)/-9.9140E-04, 2.7672E+00,-1.1544E+0
     *0/
      DATA G210( 2),G211( 2),G212( 2)/-9.9140E-04, 2.7672E+00,-1.1544E+0
     *0/
      DATA G210( 3),G211( 3),G212( 3)/-7.1017E-02, 3.4941E+00,-3.0773E+0
     *0/
      DATA G210( 4),G211( 4),G212( 4)/-7.3556E-02, 3.5487E+00,-3.1989E+0
     *0/
      DATA G210( 5),G211( 5),G212( 5)/ 3.6658E-01, 2.1162E+00,-2.0311E+0
     *0/
      DATA G210( 6),G211( 6),G212( 6)/ 1.4498E+00,-5.9717E-01,-3.2951E-0
     *1/
      DATA G210( 7),G211( 7),G212( 7)/ 1.4498E+00,-5.9717E-01,-3.2951E-0
     *1/
      DATA NG22/ 8/,B0G22/ 2.0000E+00/,B1G22/ 6.0000E+00/
      DATA G220( 1),G221( 1),G222( 1)/-5.2593E-04, 1.4285E+00,-1.2670E+0
     *0/
      DATA G220( 2),G221( 2),G222( 2)/-5.2593E-04, 1.4285E+00,-1.2670E+0
     *0/
      DATA G220( 3),G221( 3),G222( 3)/-6.4819E-02, 2.2033E+00,-3.6399E+0
     *0/
      DATA G220( 4),G221( 4),G222( 4)/ 3.7427E-02, 1.6630E+00,-2.9362E+0
     *0/
      DATA G220( 5),G221( 5),G222( 5)/ 6.1955E-01,-6.2713E-01,-6.7859E-0
     *1/
      DATA G220( 6),G221( 6),G222( 6)/ 1.7584E+00,-4.0390E+00, 1.8810E+0
     *0/
      DATA G220( 7),G221( 7),G222( 7)/ 2.5694E+00,-6.0484E+00, 3.1256E+0
     *0/
      DATA G220( 8),G221( 8),G222( 8)/ 2.5694E+00,-6.0484E+00, 3.1256E+0
     *0/
      DATA NG31/ 11/,B0G31/ 2.0000E+00/,B1G31/ 9.0000E+00/
      DATA G310( 1),G311( 1),G312( 1)/ 4.9437E-01, 1.9124E-02, 1.8375E+0
     *0/
      DATA G310( 2),G311( 2),G312( 2)/ 4.9437E-01, 1.9124E-02, 1.8375E+0
     *0/
      DATA G310( 3),G311( 3),G312( 3)/ 5.3251E-01,-6.1555E-01, 4.5595E+0
     *0/
      DATA G310( 4),G311( 4),G312( 4)/ 6.6810E-01,-2.2056E+00, 8.9293E+0
     *0/
      DATA G310( 5),G311( 5),G312( 5)/-3.8262E+00, 2.5528E+01,-3.3862E+0
     *1/
      DATA G310( 6),G311( 6),G312( 6)/ 4.2335E+00,-1.0604E+01, 6.6702E+0
     *0/
      DATA G310( 7),G311( 7),G312( 7)/ 5.0694E+00,-1.4208E+01, 1.0456E+0
     *1/
      DATA G310( 8),G311( 8),G312( 8)/ 1.4563E+00,-3.3275E+00, 2.2601E+0
     *0/
      DATA G310( 9),G311( 9),G312( 9)/-3.2852E-01, 1.2938E+00,-7.3254E-0
     *1/
      DATA G310(10),G311(10),G312(10)/-2.2489E-01, 1.0713E+00,-6.1358E-0
     *1/
      DATA G310(11),G311(11),G312(11)/-2.2489E-01, 1.0713E+00,-6.1358E-0
     *1/
      DATA NG32/ 25/,B0G32/ 2.0000E+00/,B1G32/ 2.3000E+01/
      DATA G320( 1),G321( 1),G322( 1)/ 2.9907E-05, 4.7318E-01, 6.5921E-0
     *1/
      DATA G320( 2),G321( 2),G322( 2)/ 2.9907E-05, 4.7318E-01, 6.5921E-0
     *1/
      DATA G320( 3),G321( 3),G322( 3)/ 2.5820E-03, 3.5853E-01, 1.9776E+0
     *0/
      DATA G320( 4),G321( 4),G322( 4)/-5.3270E-03, 4.9418E-01, 1.4528E+0
     *0/
      DATA G320( 5),G321( 5),G322( 5)/-6.6341E-02, 1.4422E+00,-2.2407E+0
     *0/
      DATA G320( 6),G321( 6),G322( 6)/-3.6027E-01, 4.7190E+00,-1.1380E+0
     *1/
      DATA G320( 7),G321( 7),G322( 7)/-2.7953E+00, 2.6694E+01,-6.0986E+0
     *1/
      DATA G320( 8),G321( 8),G322( 8)/-3.6091E+00, 3.4125E+01,-7.7512E+0
     *1/
      DATA G320( 9),G321( 9),G322( 9)/ 1.2491E+01,-7.1103E+01, 9.4496E+0
     *1/
      DATA G320(10),G321(10),G322(10)/ 1.9637E+01,-1.1371E+02, 1.5794E+0
     *2/
      DATA G320(11),G321(11),G322(11)/ 2.1692E+00,-2.5019E+01, 4.5340E+0
     *1/
      DATA G320(12),G321(12),G322(12)/-1.6682E+01, 6.2067E+01,-5.5257E+0
     *1/
      DATA G320(13),G321(13),G322(13)/-2.1539E+01, 8.2651E+01,-7.7065E+0
     *1/
      DATA G320(14),G321(14),G322(14)/-1.4344E+01, 5.5193E+01,-5.0867E+0
     *1/
      DATA G320(15),G321(15),G322(15)/-5.4990E+00, 2.3874E+01,-2.3140E+0
     *1/
      DATA G320(16),G321(16),G322(16)/ 3.1029E+00,-4.4708E+00, 2.1318E-0
     *1/
      DATA G320(17),G321(17),G322(17)/ 6.0961E+00,-1.3670E+01, 7.2823E+0
     *0/
      DATA G320(18),G321(18),G322(18)/ 8.6179E+00,-2.0950E+01, 1.2536E+0
     *1/
      DATA G320(19),G321(19),G322(19)/ 7.5064E+00,-1.7956E+01, 1.0520E+0
     *1/
      DATA G320(20),G321(20),G322(20)/ 5.9838E+00,-1.4065E+01, 8.0342E+0
     *0/
      DATA G320(21),G321(21),G322(21)/ 4.4959E+00,-1.0456E+01, 5.8462E+0
     *0/
      DATA G320(22),G321(22),G322(22)/ 3.2847E+00,-7.6709E+00, 4.2445E+0
     *0/
      DATA G320(23),G321(23),G322(23)/ 1.9514E+00,-4.7505E+00, 2.6452E+0
     *0/
      DATA G320(24),G321(24),G322(24)/ 4.8808E-01,-1.6910E+00, 1.0459E+0
     *0/
      DATA G320(25),G321(25),G322(25)/ 4.8808E-01,-1.6910E+00, 1.0459E+0
     *0/
      DATA NBGB/ 8/,B0BGB/ 1.5714E+00/,B1BGB/ 2.1429E-01/
      DATA BGB0( 1),BGB1( 1),BGB2( 1)/-1.0724E+00, 2.8203E+00,-3.5669E-0
     *1/
      DATA BGB0( 2),BGB1( 2),BGB2( 2)/ 3.7136E-01, 1.4560E+00,-2.8072E-0
     *2/
      DATA BGB0( 3),BGB1( 3),BGB2( 3)/ 1.1396E+00, 1.1910E+00,-5.2070E-0
     *3/
      DATA BGB0( 4),BGB1( 4),BGB2( 4)/ 1.4908E+00, 1.1267E+00,-2.2565E-0
     *3/
      DATA BGB0( 5),BGB1( 5),BGB2( 5)/ 1.7342E+00, 1.0958E+00,-1.2705E-0
     *3/
      DATA BGB0( 6),BGB1( 6),BGB2( 6)/ 1.9233E+00, 1.0773E+00,-8.1806E-0
     *4/
      DATA BGB0( 7),BGB1( 7),BGB2( 7)/ 2.0791E+00, 1.0649E+00,-5.7197E-0
     *4/
      DATA BGB0( 8),BGB1( 8),BGB2( 8)/ 2.0791E+00, 1.0649E+00,-5.7197E-0
     *4/
      DATA NPTH/ 6/,B0PTH/ 2.0000E+00/,B1PTH/ 1.8182E+01/
      DATA PTH0( 1),PTH1( 1),PTH2( 1)/ 1.0000E+00, 9.8875E-01, 2.5026E+0
     *0/
      DATA PTH0( 2),PTH1( 2),PTH2( 2)/ 1.0000E+00, 9.8875E-01, 2.5026E+0
     *0/
      DATA PTH0( 3),PTH1( 3),PTH2( 3)/ 1.0060E+00, 7.8657E-01, 4.2387E+0
     *0/
      DATA PTH0( 4),PTH1( 4),PTH2( 4)/ 1.0657E+00,-2.5051E-01, 8.7681E+0
     *0/
      DATA PTH0( 5),PTH1( 5),PTH2( 5)/ 1.6971E+00,-7.5600E+00, 2.9946E+0
     *1/
      DATA PTH0( 6),PTH1( 6),PTH2( 6)/ 1.6971E+00,-7.5600E+00, 2.9946E+0
     *1/
      DATA RMT2/1.022007/,RMSQ/.2611245/
      DATA PI/3.141593/,TWOPI/6.283185/,PI5D2/7.853982/
      DATA RM/.5110034/
      DATA IPHTER /60*0/
      DATA IEDGFL/60*0/, BKR1/3*0/,BKR2/3*0/, EKALPH/3*0/,EKBETA/3*0/
      END
      SUBROUTINE ANNIH
      DOUBLE PRECISION PAVIP,PESG1,PESG2
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/RANDOM/IXX
      PAVIP=E(NP)+PRM
      AVIP=PAVIP
      A=AVIP/RM
      AI=1.0/A
      G=A-1.0
      T=G-1.0
      P=SQRT(A*T)
      POT=P/T
      EP0=1.0/(A+P)
1011  CONTINUE
        IXX=IXX*663608943
        RNNO01=0.5 + IXX*0.23283064E-09
        EP=EP0*EXP(RNNO01*ALOG((1.0-EP0)/EP0))
        IXX=IXX*663608943
        RNNO02=0.5 + IXX*0.23283064E-09
        REJF=1.0-EP+AI*AI*(2.0*G-1.0/EP)
        IF((RNNO02.LE.REJF))GO TO1012
      GO TO 1011
1012  CONTINUE
      EP=AMAX1(EP,1.0-EP)
      ESG1=AVIP*EP
      PESG1=ESG1
      E(NP)=PESG1
      PESG2=PAVIP-PESG1
      E(NP+1)=PESG2
      ESG2=PESG2
      IQ(NP)=0
      COSTHE=AMIN1(1.00,(ESG1-RM)*POT/ESG1)
      SINTHE=SQRT(1.0-COSTHE*COSTHE)
      CALL UPHI(2,1)
      NP=NP+1
      IQ(NP)=0
      COSTHE=AMIN1(1.00,(ESG2-RM)*POT/ESG2)
      SINTHE=-SQRT(1.0-COSTHE*COSTHE)
      CALL UPHI(3,2)
      RETURN
      END
      SUBROUTINE BHABHA
      DOUBLE PRECISION PEIP,PEKSE2,PESE1,PESE2
      DOUBLE PRECISION PEKIN,H1,DCOSTH
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/RANDOM/IXX
      PEIP=E(NP)
      EIP=PEIP
      PEKIN=PEIP-PRM
      EKIN=PEKIN
      T0=EKIN/RM
      E0=T0+1.
      YY=1./(T0+2.)
      E02=E0*E0
      BETAI2=E02/(E02-1.)
      EP0=TE(MEDIUM)/EKIN
      EP0C=1.-EP0
      Y2=YY*YY
      YP=1.-2.*YY
      YP2=YP*YP
      B4=YP2*YP
      B3=B4+YP2
      B2=YP*(3.+Y2)
      B1=2.-Y2
1011  CONTINUE
        IXX=IXX*663608943
        RNNO03=0.5 + IXX*0.23283064E-09
        BR=EP0/(1.-EP0C*RNNO03)
        IXX=IXX*663608943
        RNNO04=0.5 + IXX*0.23283064E-09
        REJF2=EP0C*(BETAI2-BR*(B1-BR*(B2-BR*(B3-BR*B4))))
        IF((RNNO04.LE.REJF2))GO TO1012
      GO TO 1011
1012  CONTINUE
      IF ((BR.LT.0.5)) THEN
        IQ(NP+1)=-1
      ELSE
        IQ(NP)=-1
        IQ(NP+1)=1
        BR=1.-BR
      END IF
      BR=AMAX1(BR,0.0)
      PEKSE2=BR*EKIN
      PESE1=PEIP-PEKSE2
      PESE2=PEKSE2+PRM
      ESE1=PESE1
      ESE2=PESE2
      E(NP)=PESE1
      E(NP+1)=PESE2
      H1=(PEIP+PRM)/PEKIN
      DCOSTH=MIN(1.0D0,H1*(PESE1-PRM)/(PESE1+PRM))
      SINTHE=DSQRT(1.D0-DCOSTH)
      COSTHE=DSQRT(DCOSTH)
      CALL UPHI(2,1)
      NP=NP+1
      DCOSTH=H1*(PESE2-PRM)/(PESE2+PRM)
      SINTHE=-DSQRT(1.D0-DCOSTH)
      COSTHE=DSQRT(DCOSTH)
      CALL UPHI(3,2)
      RETURN
      END
      SUBROUTINE BREMS
      DOUBLE PRECISION PEIE,PESG,PESE
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/BREMPR/ DL1(6,3),DL2(6,3),DL3(6,3),DL4(6,3),DL5(6,3),DL6(6,
     *3), ALPHI(2,3),BPAR(2,3),DELPOS(2,3), ASYM(3,50,2), WA(3,50),PZ(3,
     *50),ZELEM(3,50),RHOZ(3,50), PWR2I(50), DELCM(3),NNE(3),ZBRANG(3),
     *IBRDST,IPRDST,IBRSPL,NBRSPL,FBRSPL
      CHARACTER*4 ASYM
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/RANDOM/IXX
      DATA AI2LN2/0.7213475/
      PEIE=E(NP)
      EIE=PEIE
      NP=NP+1
      IF ((EIE.LT.50.0)) THEN
        LVX=1
        LVL0=0
      ELSE
        LVX=2
        LVL0=3
      END IF
      ABREMS=FLOAT(IFIX(1.44269*ALOG(EIE/AP(MEDIUM))))
1011  CONTINUE
        IXX=IXX*663608943
        RNNO06=0.5 + IXX*0.23283064E-09
        IF ((0.5.LT.((ABREMS*ALPHI(LVX,MEDIUM)+0.5)*RNNO06))) THEN
          IXX=IXX*663608943
          RNNO07=0.5 + IXX*0.23283064E-09
          IDISTR=ABREMS*RNNO07
          P=PWR2I(IDISTR+1)
          LVL=LVL0+1
          IXX=IXX*663608943
          RNNO08=0.5 + IXX*0.23283064E-09
          IF ((RNNO08.GE.AI2LN2)) THEN
1021        CONTINUE
              IXX=IXX*663608943
              RNNO09=0.5 + IXX*0.23283064E-09
              IXX=IXX*663608943
              RNNO10=0.5 + IXX*0.23283064E-09
              IXX=IXX*663608943
              RNNO11=0.5 + IXX*0.23283064E-09
              H=AMAX1(RNNO10,RNNO11)
              BR=1.0-0.5*H
              IF((RNNO09.LE.0.5/BR))GO TO1022
            GO TO 1021
1022        CONTINUE
          ELSE
            IXX=IXX*663608943
            RNNO12=0.5 + IXX*0.23283064E-09
            BR=RNNO12*0.5
          END IF
          BR=BR*P
        ELSE
          IXX=IXX*663608943
          RNNO13=0.5 + IXX*0.23283064E-09
          IXX=IXX*663608943
          RNNO14=0.5 + IXX*0.23283064E-09
          BR=AMAX1(RNNO13,RNNO14)
          LVL=LVL0+2
        END IF
        ESG=EIE*BR
        IF((ESG.LT.AP(MEDIUM)))GO TO1011
        PESG=ESG
        PESE=PEIE-PESG
        ESE=PESE
        IF((ESE.LT.RM))GO TO1011
        DEL = BR/ESE
        IF((DEL.GE.DELPOS(LVX,MEDIUM)))GO TO1011
        DELTA = DELCM(MEDIUM)*DEL
        IF ((DELTA.LT.1.0)) THEN
          REJF=DL1(LVL,MEDIUM)+DELTA*(DL2(LVL,MEDIUM)+DELTA*DL3(LVL,MEDI
     *    UM))
        ELSE
          REJF=DL4(LVL,MEDIUM)+DL5(LVL,MEDIUM)*ALOG(DELTA+DL6(LVL,MEDIUM
     *    ))
        END IF
        IXX=IXX*663608943
        RNSCRN=0.5 + IXX*0.23283064E-09
        IF((RNSCRN.LE.REJF))GO TO1012
      GO TO 1011
1012  CONTINUE
      IF ((IBRDST.EQ.0)) THEN
        THETA=RM/EIE
      ELSE IF((IBRDST.EQ.1)) THEN
        ZTARG=ZBRANG(MEDIUM)
        TTEIE=EIE/RM
        TTESE=ESE/RM
        ESEDEI=TTESE/TTEIE
        Y2MAX=(PI*TTEIE)**2
        RJARG1=(1.0+ESEDEI**2)
        RJARG2=3.0*RJARG1-2.0*ESEDEI
        RJARG3=((1.0-ESEDEI)/(2.0*TTEIE*ESEDEI))**2
        Y2TST1=(1.+0.0E0)**2
        REJMIN= (4.+LOG(RJARG3+ZTARG/Y2TST1))*(4.*ESEDEI*0.0E0/Y2TST1-RJ
     *  ARG1)+RJARG2
        Y2TST1=(1.+1.0E0)**2
        REJMID= (4.+LOG(RJARG3+ZTARG/Y2TST1))*(4.*ESEDEI*1.0E0/Y2TST1-RJ
     *  ARG1)+RJARG2
        Y2TST1=(1.+Y2MAX)**2
        REJMAX= (4.+LOG(RJARG3+ZTARG/Y2TST1))*(4.*ESEDEI*Y2MAX/Y2TST1-RJ
     *  ARG1)+RJARG2
        REJTOP=MAX(REJMIN,REJMID,REJMAX)
1031    CONTINUE
          IXX=IXX*663608943
          Y2TST=0.5 + IXX*0.23283064E-09
          Y2TST=Y2TST/(1.0-Y2TST+1.0/Y2MAX)
          Y2TST1=(1.+Y2TST)**2
          REJTST= (4.+LOG(RJARG3+ZTARG/Y2TST1))*(4.*ESEDEI*Y2TST/Y2TST1-
     *    RJARG1)+RJARG2
          IXX=IXX*663608943
          RTEST=0.5 + IXX*0.23283064E-09
          IF(((RTEST.LE.(REJTST/REJTOP))))GO TO1032
        GO TO 1031
1032    CONTINUE
        THETA=SQRT(Y2TST)/TTEIE
      END IF
      CALL UPHI(1,3)
      IF ((ESG.LE.ESE)) THEN
        IQ(NP)=0
        E(NP)=PESG
        E(NP-1)=PESE
      ELSE
        IQ(NP)=IQ(NP-1)
        IQ(NP-1)=0
        E(NP)=PESE
        E(NP-1)=PESG
        T=U(NP)
        U(NP)=U(NP-1)
        U(NP-1)=T
        T=V(NP)
        V(NP)=V(NP-1)
        V(NP-1)=T
        T=W(NP)
        W(NP)=W(NP-1)
        W(NP-1)=T
      END IF
      RETURN
      END
      SUBROUTINE COMPT
      DOUBLE PRECISION PEIG,PESG,PESE
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/RANDOM/IXX
      PEIG=E(NP)
      EIG=PEIG
      EGP=EIG/RM
      BR0I=1.+2.*EGP
      ALPH1=ALOG(BR0I)
      ALPH2=EGP*(BR0I+1.)/(BR0I*BR0I)
      SUMALP = ALPH1+ALPH2
1011  CONTINUE
        IXX=IXX*663608943
        RNNO15=0.5 + IXX*0.23283064E-09
        IF ((ALPH1.GE.SUMALP*RNNO15)) THEN
          IXX=IXX*663608943
          RNNO16=0.5 + IXX*0.23283064E-09
          BR=EXP(ALPH1*RNNO16)/BR0I
        ELSE
          IXX=IXX*663608943
          BRP=0.5 + IXX*0.23283064E-09
          IXX=IXX*663608943
          RNNO18=0.5 + IXX*0.23283064E-09
          IF ((EGP.GE.(EGP+1.)*RNNO18)) THEN
            IXX=IXX*663608943
            RNNO19=0.5 + IXX*0.23283064E-09
            BRP=AMAX1(BRP,RNNO19)
          END IF
          BR=((BR0I-1.)*BRP+1.)/BR0I
        END IF
        ESG=BR*EIG
        A1MIBR = 1.-BR
        TEMP=RM*A1MIBR/ESG
        SINTHE=AMAX1(0.0,TEMP*(2.0-TEMP))
        IXX=IXX*663608943
        RNNO20=0.5 + IXX*0.23283064E-09
        REJF3=1.0-BR*SINTHE/(1.0+BR*BR)
        IF((RNNO20.LE.REJF3))GO TO1012
      GO TO 1011
1012  CONTINUE
      SINTHE=SQRT(SINTHE)
      COSTHE=1.0-TEMP
      PESG=ESG
      PESE=PEIG-PESG+PRM
      ESE=PESE
      CALL UPHI(2,1)
      NP=NP+1
      PSQ=ESE**2-RMSQ
      IF ((PSQ.LE.0.0)) THEN
        COSTHE=0.
        SINTHE=-1.
      ELSE
        COSTHE=(ESE+ESG)*A1MIBR/SQRT(PSQ)
        SINTHE=-SQRT(AMAX1(0.0,1.0-COSTHE*COSTHE))
      END IF
      CALL UPHI(3,2)
      IF ((ESE.LE.ESG)) THEN
        IQ(NP)=-1
        E(NP)=PESE
        E(NP-1)=PESG
      ELSE
        IQ(NP)=0
        IQ(NP-1)=-1
        E(NP)=PESG
        E(NP-1)=PESE
        T=U(NP)
        U(NP)=U(NP-1)
        U(NP-1)=T
        T=V(NP)
        V(NP)=V(NP-1)
        V(NP-1)=T
        T=W(NP)
        W(NP)=W(NP-1)
        W(NP-1)=T
      END IF
      RETURN
      END
      SUBROUTINE ELECTR(IRCODE)
      DOUBLE PRECISION PEIE
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/BOUNDS/ECUT(60),PCUT(60),VACDST
      COMMON/BREMPR/ DL1(6,3),DL2(6,3),DL3(6,3),DL4(6,3),DL5(6,3),DL6(6,
     *3), ALPHI(2,3),BPAR(2,3),DELPOS(2,3), ASYM(3,50,2), WA(3,50),PZ(3,
     *50),ZELEM(3,50),RHOZ(3,50), PWR2I(50), DELCM(3),NNE(3),ZBRANG(3),
     *IBRDST,IPRDST,IBRSPL,NBRSPL,FBRSPL
      CHARACTER*4 ASYM
      INTEGER EPSTFL(3)
      COMMON/ELECIN/ EKELIM, EKE0(3),EKE1(3),CMFP0(3),CMFP1(3),RANGE0(3)
     *,RANGE1(3), XR0(3),TEFF0(3),BLCC(3),XCC(3), PICMP0(1,3),PICMP1(1,3
     *),EICMP0(1,3),EICMP1(1,3), ESIG0(150,3),ESIG1(150,3),PSIG0(150,3),
     *PSIG1(150,3),EDEDX0(150,3),EDEDX1(150,3),PDEDX0(150,3),PDEDX1(150,
     *3),EBR10(150,3),EBR11(150,3),PBR10(150,3),PBR11(150,3),PBR20(150,3
     *),PBR21(150,3),TMXS0(150,3),TMXS1(150,3), CMFPE0(1,3),CMFPE1(1,3),
     *CMFPP0(1,3),CMFPP1(1,3),ERANG0(1,3),ERANG1(1,3),PRANG0(1,3),PRANG1
     *(1,3),CXC2E0(1,3),CXC2E1(1,3),CXC2P0(1,3),CXC2P1(1,3),CLXAE0(1,3),
     *CLXAE1(1,3),CLXAP0(1,3),CLXAP1(1,3), THR0(1,1),THR1(1,1),THR2(1,1)
     *, THRI0(1,1),THRI1(1,1),THRI2(1,1), FSTEP(16),FSQR(16), VERT1(1000
     *),VERT2(100,16), BLC0,BLC1,RTHR0,RTHR1,RTHRI0,RTHRI1, ICOMP, MPEEM
     *(1,3), MSMAP(200), MSTEPS,JRMAX,MXV1, MXV2,NBLC,NRNTH,NRNTHI, IUNR
     *ST(3),EPSTFL,IAPRIM(3)
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/MEDIA/ RLC(3),RLDU(3),RHO(3),MSGE(3),MGE(3),MSEKE(3),MEKE(3
     *),MLEKE(3),MCMFP(3),MRANGE(3),IRAYLM(3),NMED
      COMMON/MEDIAC/MEDIA(24,3)
      CHARACTER*4 MEDIA
      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(60),RHOR(60),IRAYLR(60)
      COMMON/PATHCM/B0PTH,B1PTH,PTH0(6),PTH1(6),PTH2(6),NPTH
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIIN/SINC0,SINC1,SIN0(1002),SIN1(1002)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/USERSC/ SMAX,SMAXIR(60),ESTEPE,ESTEPR(60),ESAVE(60), IEDEP,
     *NOMSCT(60),NOPLC(60)
      COMMON/USERVR/ CEXPTR,GWAIT,IFORCE,NFMIN,NFMAX,NFTIME,ISOURC,IFPB,
     *IQINC,MONOEN
      COMMON/USERXT/IPHTER(60)
      COMMON/RANDOM/IXX
      IRCODE=1
      IROLD=IR(NP)
      IRL=IR(NP)
      MEDIUM=MED(IRL)
1010  CONTINUE
1011    CONTINUE
        LELEC=IQ(NP)
        PEIE=E(NP)
        EIE=PEIE
        IF ((EIE.LE.ECUT(IRL))) THEN
          GO TO 1020
        END IF
        MEDIUM=MED(IRL)
1030    CONTINUE
1031      CONTINUE
          IF ((MEDIUM.NE.0)) THEN
            IXX=IXX*663608943
            RNNE1=0.5 + IXX*0.23283064E-09
            IF ((RNNE1.EQ.0.0)) THEN
              RNNE1=1.E-30
            END IF
            DEMFP=AMAX1(-ALOG(RNNE1),1.E-6)
            EKE=EIE-RM
            ELKE=ALOG(EKE)
            LELKE=EKE1(MEDIUM)*ELKE+EKE0(MEDIUM)
            IF ((LELEC.LT.0)) THEN
              SIG0=ESIG1(LELKE,MEDIUM)*ELKE+ESIG0(LELKE,MEDIUM)
            ELSE
              SIG0=PSIG1(LELKE,MEDIUM)*ELKE+PSIG0(LELKE,MEDIUM)
            END IF
          END IF
1040      CONTINUE
1041        CONTINUE
            IF ((MEDIUM.EQ.0)) THEN
              TSTEP=VACDST
              USTEP=TSTEP
              TUSTEP=USTEP
            ELSE
              RHOF=RHOR(IRL)/RHO(MEDIUM)
              SIG=SIG0*RHOF
              IF ((SIG.LE.0.0)) THEN
                TSTEP=VACDST
              ELSE
                TSTEP=DEMFP/SIG
              END IF
              TMXS=TMXS1(LELKE,MEDIUM)*ELKE+TMXS0(LELKE,MEDIUM)
              TP=200*TEFF0(MEDIUM)
              IF ((ESTEPR(IRL).NE.0)) THEN
                TMXS=TMXS*ESTEPR(IRL)
              END IF
              TMXS = AMIN1(TMXS,TP,SMAX,SMAXIR(IRL))
              TMXS=TMXS/RHOF
              TUSTEP=AMIN1(TSTEP,TMXS)
              IF ((LELEC.LT.0)) THEN
                DEDX0=EDEDX1(LELKE,MEDIUM)*ELKE+EDEDX0(LELKE,MEDIUM)
              ELSE
                DEDX0=PDEDX1(LELKE,MEDIUM)*ELKE+PDEDX0(LELKE,MEDIUM)
              END IF
              DEDX=RHOF*DEDX0
              RANGE=(EIE-ECUT(IRL)+0.0001)/DEDX
              BETA2=AMAX1(1.E-8,1.-RMSQ/EIE**2)
              TSCAT=RLDU(MEDIUM)*(EIE*BETA2*0.094315)**2
              TSCAT=TSCAT/RHOF
              TUSTEP=AMIN1(TUSTEP,0.3*TSCAT,RANGE)
              USTEP=TUSTEP*(1.-TUSTEP/TSCAT)
            END IF
            IRNEW=IR(NP)
            IDISC=0
            USTEP0=USTEP
            IF (((USTEP.GT.DNEAR(NP)  ))) THEN
              CALL HOWFAR
            END IF
            IF ((RANGE.LT.DNEAR(NP) .AND. E(NP).LE.ESAVE(IRL) .AND. MEDI
     *      UM.NE.0)) THEN
              IF ((LELEC.EQ.-1)) THEN
                IDISC=1
              ELSE
                IDISC=99
              END IF
            END IF
            IF ((IDISC.GT.0)) THEN
              GO TO 1050
            END IF
            IF ((USTEP.LE.0.0)) THEN
              IF ((USTEP.LT.-1.E-4)) THEN
                IERUST=IERUST+1
                WRITE(6,1060)IERUST,USTEP,IR(NP),IRNEW, IROLD,X(NP),Y(NP
     *          ),Z(NP),SQRT(X(NP)**2+Y(NP)**2)
1060            FORMAT(I6,' NEGATIVE USTEP=',E12.6,' IR,IRNEW,IROLD=',3I
     *          4,'X,Y,Z,R=', 4E10.3)
                IF ((IERUST.GT.1000)) THEN
                  WRITE(6,1070)
1070              FORMAT(///'0CALLED EXIT, TOO MANY USTEP ERRORS'///)
                  CALL EXIT
                END IF
              END IF
              USTEP=0.0
            END IF
            IF ((USTEP.EQ.0.0.OR.MEDIUM.EQ.0)) THEN
              IF ((USTEP.NE.0.0)) THEN
                VSTEP=USTEP
                TVSTEP=VSTEP
                EDEP=PZERO
                IARG=0
                IF ((IAUSFL(IARG+1).NE.0)) THEN
                  CALL AUSGAB(IARG)
                END IF
                X(NP)=X(NP)+U(NP)*VSTEP
                Y(NP)=Y(NP)+V(NP)*VSTEP
                Z(NP)=Z(NP)+W(NP)*VSTEP
                DNEAR(NP)=DNEAR(NP)-VSTEP
                IROLD=IR(NP)
              END IF
              IR(NP)=IRNEW
              IRL=IRNEW
              MEDIUM=MED(IRL)
              IF ((EIE.LE.ECUT(IRL))) THEN
                GO TO 1020
              END IF
              IF ((USTEP.NE.0.0)) THEN
                IARG=5
                IF ((IAUSFL(IARG+1).NE.0)) THEN
                  CALL AUSGAB(IARG)
                END IF
              END IF
              GO TO 1031
            END IF
            VSTEP=USTEP
            IF ((USTEP.EQ.USTEP0)) THEN
              TVSTEP=TUSTEP
            ELSE
              VSTP=VSTEP/TSCAT
              IPTH=B0PTH+B1PTH*VSTP
              IF ((IPTH.GT.NPTH)) THEN
                WRITE(6,1080)VSTP,IPTH,NPTH
1080            FORMAT(' OUT OF BOUNDS IPTH: VSTP,IPTH,NPTH=',1PG15.6,2I
     *          10)
                STOP
              END IF
              PTH=PTH0(IPTH)+VSTP*(PTH1(IPTH)+VSTP*PTH2(IPTH))
              TVSTEP=PTH*VSTEP
            END IF
            IF (( (IEDEP.LE.0).AND.(EKE.LT.1.0) )) THEN
              DDXTMP=DEDX
              EKETMP=EKE
              ELKTMP=ELKE
              EKE=EKETMP-0.5*DEDX*TVSTEP
              ELKE=ALOG(EKE)
              LELKE=EKE1(MEDIUM)*ELKE+EKE0(MEDIUM)
              IF ((LELEC.LT.0)) THEN
                DEDX01=EDEDX1(LELKE,MEDIUM)*ELKE+EDEDX0(LELKE,MEDIUM)
              ELSE
                DEDX01=PDEDX1(LELKE,MEDIUM)*ELKE+PDEDX0(LELKE,MEDIUM)
              END IF
              DEDX=RHOF*DEDX01
              DEDX=0.5*(DDXTMP+DEDX)
              EKE=EKETMP
              ELKE=ELKTMP
            END IF
            DE=DEDX*TVSTEP
            EDEP=DE
            EKEF=EKE-DE
            EOLD=EIE
            ENEW=EOLD-DE
            CALL MSCAT
            IARG=0
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
            X(NP)=X(NP)+U(NP)*VSTEP
            Y(NP)=Y(NP)+V(NP)*VSTEP
            Z(NP)=Z(NP)+W(NP)*VSTEP
            DNEAR(NP)=DNEAR(NP)-VSTEP
            IROLD=IR(NP)
            CALL UPHI(2,1)
            PEIE=PEIE-EDEP
            EIE=PEIE
            E(NP)=PEIE
            IF ((EIE.LE.ECUT(IRL))) THEN
              GO TO 1020
            END IF
            MEDOLD=MEDIUM
            IF ((MEDIUM.NE.0)) THEN
              EKEOLD=EKE
              EKE=EIE-RM
              ELKE=ALOG(EKE)
              LELKE=EKE1(MEDIUM)*ELKE+EKE0(MEDIUM)
            END IF
            IF ((IRNEW.NE.IROLD)) THEN
              IR(NP)=IRNEW
              IRL=IRNEW
              MEDIUM=MED(IRL)
            END IF
            IF ((EIE.LE.ECUT(IRL))) THEN
              GO TO 1020
            END IF
            IARG=5
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
            IF((IDISC.LT.0))GO TO 1050
            IF((MEDIUM.NE.MEDOLD))GO TO 1031
            DEMFP=AMAX1(0.,DEMFP-TVSTEP*SIG)
            IF(((DEMFP.LT.1.E-6)))GO TO1042
          GO TO 1041
1042      CONTINUE
          IF ((LELEC.LT.0)) THEN
            SIGF=ESIG1(LELKE,MEDIUM)*ELKE+ESIG0(LELKE,MEDIUM)
          ELSE
            SIGF=PSIG1(LELKE,MEDIUM)*ELKE+PSIG0(LELKE,MEDIUM)
          END IF
          IXX=IXX*663608943
          RFICT=0.5 + IXX*0.23283064E-09
          IF(((RFICT.LE.SIGF/SIG0)))GO TO1032
        GO TO 1031
1032    CONTINUE
        IF ((LELEC.LT.0)) THEN
          EBR1=EBR11(LELKE,MEDIUM)*ELKE+EBR10(LELKE,MEDIUM)
          IF((EKE.LE.AP(MED(IRL))))EBR1=0.0
          IXX=IXX*663608943
          RNNO24=0.5 + IXX*0.23283064E-09
          IF ((RNNO24.LT.EBR1)) THEN
            GO TO 1090
          ELSE
            IF ((E(NP).LE.THMOLL(MEDIUM))) THEN
              IF ((EBR1.LE.0.0)) THEN
                GO TO 1010
              END IF
              GO TO 1090
            END IF
            IARG=8
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
            CALL MOLLER
            IARG=9
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
          END IF
          GO TO 1010
        END IF
        PBR1=PBR11(LELKE,MEDIUM)*ELKE+PBR10(LELKE,MEDIUM)
        IF((EKE.LE.AP(MED(IRL))))PBR1=0.0
        IXX=IXX*663608943
        RNNO25=0.5 + IXX*0.23283064E-09
        IF ((RNNO25.LT.PBR1)) THEN
          GO TO 1090
        END IF
        PBR2=PBR21(LELKE,MEDIUM)*ELKE+PBR20(LELKE,MEDIUM)
        IF ((RNNO25.LT.PBR2)) THEN
          IARG=10
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          CALL BHABHA
          IARG=11
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
        ELSE
          IARG=12
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          CALL ANNIH
          IARG=13
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          GO TO 1012
        END IF
      GO TO 1011
1012  CONTINUE
      RETURN
1090  IARG=6
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      CALL BREMS
      IF ((IBRSPL.EQ.1)) THEN
        IF (((NBRSPL.GT.1).AND.((NP+NBRSPL).GE.125))) THEN
1101      CONTINUE
            WRITE(6,1110)125,NBRSPL,(2*NBRSPL+1)/3
1110        FORMAT('0*** WARNING ***. STACK SIZE = ',I4,' MIGHT OVERFLOW
     *'/ '                 NBRSPL BEING REDUCED, ',I4,'-->',I4/)
            NBRSPL=(2*NBRSPL+1)/3
            FBRSPL=1.0/FLOAT(NBRSPL)
            IF ((NBRSPL.EQ.1)) THEN
              WRITE(6,1120)125
1120          FORMAT('0*** WARNING ***. STACK SIZE = ',I4,' IS TOO SMALL
     *'/ '                 BREMSSTRAHLUNG SPLITTING NOW SHUT OFF'/)
              IBRSPL=0
            END IF
            IF((((NP+NBRSPL).LT.125)))GO TO1102
          GO TO 1101
1102      CONTINUE
        END IF
        IF ((IQ(NP).EQ.0)) THEN
          NPSTRT=NP-1
          FDUMMY = U(NP-1)
          U(NP-1) = U(NP)
          U(NP) = FDUMMY
          FDUMMY = V(NP-1)
          V(NP-1) = V(NP)
          V(NP) = FDUMMY
          FDUMMY = W(NP-1)
          W(NP-1) = W(NP)
          W(NP) = FDUMMY
          FDUMMY = E(NP-1)
          E(NP-1) = E(NP)
          E(NP) = FDUMMY
          FDUMMY = WT(NP-1)
          WT(NP-1) = WT(NP)
          WT(NP) = FDUMMY
          IDUMMY = IQ(NP-1)
          IQ(NP-1) = IQ(NP)
          IQ(NP) = IDUMMY
          IDUMMY = LATCH(NP-1)
          LATCH(NP-1) = LATCH(NP)
          LATCH(NP) = IDUMMY
        ELSE
          NPSTRT=NP
        END IF
        WT(NP-1)=WT(NP-1)*FBRSPL
        FRSTBR=E(NP-1)
        E(NP)=E(NP)+E(NP-1)
        IARG=7
        IF ((IAUSFL(IARG+1).NE.0)) THEN
          CALL AUSGAB(IARG)
        END IF
        ICSPLT=1
1131    IF(ICSPLT.GE.NBRSPL)GO TO 1132
          ICSPLT=ICSPLT+1
          IARG=6
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          CALL BREMS
          IF ((IQ(NP).EQ.0)) THEN
            FDUMMY = U(NP-1)
            U(NP-1) = U(NP)
            U(NP) = FDUMMY
            FDUMMY = V(NP-1)
            V(NP-1) = V(NP)
            V(NP) = FDUMMY
            FDUMMY = W(NP-1)
            W(NP-1) = W(NP)
            W(NP) = FDUMMY
            FDUMMY = E(NP-1)
            E(NP-1) = E(NP)
            E(NP) = FDUMMY
            FDUMMY = WT(NP-1)
            WT(NP-1) = WT(NP)
            WT(NP) = FDUMMY
            IDUMMY = IQ(NP-1)
            IQ(NP-1) = IQ(NP)
            IQ(NP) = IDUMMY
            IDUMMY = LATCH(NP-1)
            LATCH(NP-1) = LATCH(NP)
            LATCH(NP) = IDUMMY
          END IF
          WT(NP-1)=WT(NP-1)*FBRSPL
          E(NP)=E(NP)+E(NP-1)
          IARG=7
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
        GO TO 1131
1132    CONTINUE
        E(NP)=E(NP)-FRSTBR
        FDUMMY = U(NP)
        U(NP) = U(NPSTRT)
        U(NPSTRT) = FDUMMY
        FDUMMY = V(NP)
        V(NP) = V(NPSTRT)
        V(NPSTRT) = FDUMMY
        FDUMMY = W(NP)
        W(NP) = W(NPSTRT)
        W(NPSTRT) = FDUMMY
        FDUMMY = E(NP)
        E(NP) = E(NPSTRT)
        E(NPSTRT) = FDUMMY
        FDUMMY = WT(NP)
        WT(NP) = WT(NPSTRT)
        WT(NPSTRT) = FDUMMY
        IDUMMY = IQ(NP)
        IQ(NP) = IQ(NPSTRT)
        IQ(NPSTRT) = IDUMMY
        IDUMMY = LATCH(NP)
        LATCH(NP) = LATCH(NPSTRT)
        LATCH(NPSTRT) = IDUMMY
      END IF
      IARG=7
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      IF ((IQ(NP).EQ.0)) THEN
        RETURN
      ELSE
        GO TO 1010
      END IF
1020  IF ((EIE.GT.AE(MEDIUM))) THEN
        IDR=1
        IF ((LELEC.LT.0)) THEN
          EDEP=PEIE-PRM
        ELSE
          EDEP=PEIE-PRM
        END IF
      ELSE
        IDR=2
        EDEP=PEIE-PRM
      END IF
      IARG=IDR
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
1140  CONTINUE
      IF ((LELEC.GT.0)) THEN
        IF ((EDEP.LT.PEIE)) THEN
          IXX=IXX*663608943
          COSTHE=0.5 + IXX*0.23283064E-09
          IXX=IXX*663608943
          FLIP=0.5 + IXX*0.23283064E-09
          IF ((FLIP.LE.0.5)) THEN
            COSTHE=-COSTHE
          END IF
          SINTHE=SQRT(1.0-COSTHE**2)
          E(NP)=PRM
          IQ(NP)=0
          U(NP)=0.
          V(NP)=0.
          W(NP)=1.
          CALL UPHI(2,1)
          NP=NP+1
          E(NP)=PRM
          IQ(NP)=0
          X(NP)=X(NP-1)
          Y(NP)=Y(NP-1)
          Z(NP)=Z(NP-1)
          IR(NP)=IR(NP-1)
          WT(NP)=WT(NP-1)
          DNEAR(NP)=DNEAR(NP-1)
          LATCH(NP)=LATCH(NP-1)
          U(NP)=-U(NP-1)
          V(NP)=-V(NP-1)
          W(NP)=-W(NP-1)
          IARG=14
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          RETURN
        END IF
      END IF
      NP=NP-1
      IRCODE=2
      RETURN
1050  IDISC=ABS(IDISC)
      IF (((LELEC.LT.0).OR.(IDISC.EQ.99))) THEN
        EDEP=PEIE-PRM
      ELSE
        EDEP=PEIE+PRM
      END IF
      IARG=3
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      IF((IDISC.EQ.99))GOTO 1140
      IRCODE=2
      NP=NP-1
      RETURN
      END
      SUBROUTINE HATCH
      CHARACTER*4 MBUF(72),MDLABL(8)
      DIMENSION ZEROS(3),LOK(3)
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/BOUNDS/ECUT(60),PCUT(60),VACDST
      COMMON/BREMPR/ DL1(6,3),DL2(6,3),DL3(6,3),DL4(6,3),DL5(6,3),DL6(6,
     *3), ALPHI(2,3),BPAR(2,3),DELPOS(2,3), ASYM(3,50,2), WA(3,50),PZ(3,
     *50),ZELEM(3,50),RHOZ(3,50), PWR2I(50), DELCM(3),NNE(3),ZBRANG(3),
     *IBRDST,IPRDST,IBRSPL,NBRSPL,FBRSPL
      CHARACTER*4 ASYM
      INTEGER EPSTFL(3)
      COMMON/ELECIN/ EKELIM, EKE0(3),EKE1(3),CMFP0(3),CMFP1(3),RANGE0(3)
     *,RANGE1(3), XR0(3),TEFF0(3),BLCC(3),XCC(3), PICMP0(1,3),PICMP1(1,3
     *),EICMP0(1,3),EICMP1(1,3), ESIG0(150,3),ESIG1(150,3),PSIG0(150,3),
     *PSIG1(150,3),EDEDX0(150,3),EDEDX1(150,3),PDEDX0(150,3),PDEDX1(150,
     *3),EBR10(150,3),EBR11(150,3),PBR10(150,3),PBR11(150,3),PBR20(150,3
     *),PBR21(150,3),TMXS0(150,3),TMXS1(150,3), CMFPE0(1,3),CMFPE1(1,3),
     *CMFPP0(1,3),CMFPP1(1,3),ERANG0(1,3),ERANG1(1,3),PRANG0(1,3),PRANG1
     *(1,3),CXC2E0(1,3),CXC2E1(1,3),CXC2P0(1,3),CXC2P1(1,3),CLXAE0(1,3),
     *CLXAE1(1,3),CLXAP0(1,3),CLXAP1(1,3), THR0(1,1),THR1(1,1),THR2(1,1)
     *, THRI0(1,1),THRI1(1,1),THRI2(1,1), FSTEP(16),FSQR(16), VERT1(1000
     *),VERT2(100,16), BLC0,BLC1,RTHR0,RTHR1,RTHRI0,RTHRI1, ICOMP, MPEEM
     *(1,3), MSMAP(200), MSTEPS,JRMAX,MXV1, MXV2,NBLC,NRNTH,NRNTHI, IUNR
     *ST(3),EPSTFL,IAPRIM(3)
      COMMON/MEDIA/ RLC(3),RLDU(3),RHO(3),MSGE(3),MGE(3),MSEKE(3),MEKE(3
     *),MLEKE(3),MCMFP(3),MRANGE(3),IRAYLM(3),NMED
      COMMON/MEDIAC/MEDIA(24,3)
      CHARACTER*4 MEDIA
      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(60),RHOR(60),IRAYLR(60)
      COMMON/PHOTIN/ EBINDA(3), GE0(3),GE1(3), GMFP0(200,3),GMFP1(200,3)
     *,GBR10(200,3),GBR11(200,3),GBR20(200,3),GBR21(200,3), RCO0(3),RCO1
     *(3), RSCT0(100,3),RSCT1(100,3), COHE0(200,3),COHE1(200,3), MPGEM(1
     *,3), NGR(3)
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIIN/SINC0,SINC1,SIN0(1002),SIN1(1002)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/USERSC/ SMAX,SMAXIR(60),ESTEPE,ESTEPR(60),ESAVE(60), IEDEP,
     *NOMSCT(60),NOPLC(60)
      COMMON/USERVR/ CEXPTR,GWAIT,IFORCE,NFMIN,NFMAX,NFTIME,ISOURC,IFPB,
     *IQINC,MONOEN
      COMMON/USERXT/IPHTER(60)
      COMMON/RANDOM/IXX
      DATA MDLABL/' ','M','E','D','I','U','M','='/,LMDL/8/,LMDN/24/,DUNI
     *TO/1./
      DATA I1ST/1/,NSINSS/37/,MXSINC/1002/,ISTEST/0/,NRNA/1000/
1010  FORMAT(1X,14I5)
1020  FORMAT(1X,1PE14.5,4E14.5)
1030  FORMAT(72A1)
      IF ((I1ST.NE.0)) THEN
        I1ST=0

c       CHANGE INITIALISATION FOR RANDOM GENERATOR
c       READ VALUE FORM VALUE RSEED
        OPEN(UNIT=29,FILE='RANDOM_SEED',STATUS='OLD')
        READ(29,*) IXX
        CLOSE(29)

        IF ((IXX.EQ.0)) THEN
          IXX=987654321
        END IF
        IF ((SMAX.LE.0.0)) THEN
          SMAX=1.E10
        END IF
          DO 1041 J=1,60
          IF ((SMAXIR(J).LE.0.0)) THEN
            SMAXIR(J)=SMAX
          END IF
1041    CONTINUE
1042    CONTINUE
        LATCHI=0.0
        PRM=RM
        PRMT2=2.D0*PRM
        PZERO=0.0D0
        NISUB=MXSINC-2
        FNSSS=NSINSS
        WID=PI5D2/FLOAT(NISUB)
        WSS=WID/(FNSSS-1.0)
        ZEROS(1)=0.
        ZEROS(2)=PI
        ZEROS(3)=TWOPI
          DO 1051 ISUB=1,MXSINC
          SX=0.
          SY=0.
          SXX=0.
          SXY=0.
          XS0=WID*FLOAT(ISUB-2)
          XS1=XS0+WID
          IZ=0
            DO 1061 IZZ=1,3
            IF (((XS0.LE.ZEROS(IZZ)).AND.(ZEROS(IZZ).LE.XS1))) THEN
              IZ=IZZ
              GO TO1062
            END IF
1061      CONTINUE
1062      CONTINUE
          IF ((IZ.EQ.0)) THEN
            XSI=XS0
          ELSE
            XSI=ZEROS(IZ)
          END IF
            DO 1071 ISS=1,NSINSS
            XS=WID*FLOAT(ISUB-2)+WSS*FLOAT(ISS-1)-XSI
            YS=SIN(XS+XSI)
            SX=SX+XS
            SY=SY+YS
            SXX=SXX+XS*XS
            SXY=SXY+XS*YS
1071      CONTINUE
1072      CONTINUE
          IF ((IZ.NE.0)) THEN
            SIN1(ISUB)=SXY/SXX
            SIN0(ISUB)=-SIN1(ISUB)*XSI
          ELSE
            DEL=FNSSS*SXX-SX*SX
            SIN1(ISUB)=(FNSSS*SXY-SY*SX)/DEL
            SIN0(ISUB)=(SY*SXX-SX*SXY)/DEL - SIN1(ISUB)*XSI
          END IF
1051    CONTINUE
1052    CONTINUE
        SINC0=2.0
        SINC1=1.0/WID
        IF ((ISTEST.NE.0)) THEN
          ADEV=0.
          RDEV=0.
          S2C2MN=10.
          S2C2MX=0.
            DO 1081 ISUB=1,NISUB
              DO 1091 ISS=1,NSINSS
              THETA=WID*FLOAT(ISUB-1)+WSS*FLOAT(ISS-1)
              CTHET=PI5D2-THETA
              LTHETA=SINC1*THETA+SINC0
              LCTHET=SINC1*CTHET+SINC0
              SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
              COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
              SINT=SIN(THETA)
              COST=COS(THETA)
              ASD=ABS(SINTHE-SINT)
              ACD=ABS(COSTHE-COST)
              ADEV=AMAX1(ADEV,ASD,ACD)
              IF((SINT.NE.0.0))RDEV=AMAX1(RDEV,ASD/ABS(SINT))
              IF((COST.NE.0.0))RDEV=AMAX1(RDEV,ACD/ABS(COST))
              S2C2=SINTHE**2+COSTHE**2
              S2C2MN=AMIN1(S2C2MN,S2C2)
              S2C2MX=AMAX1(S2C2MX,S2C2)
              IF ((ISUB.LT.11)) THEN
                WRITE(6,1100)THETA,SINTHE,SINT,COSTHE,COST
1100            FORMAT(1PE20.7,4E20.7)
              END IF
1091        CONTINUE
1092        CONTINUE
1081      CONTINUE
1082      CONTINUE
          WRITE(6,1110)MXSINC,NSINSS
1110      FORMAT(' SINE TESTS,MXSINC,NSINSS=',2I5)
          WRITE(6,1120)ADEV,RDEV,S2C2MN,S2C2MX
1120      FORMAT(' ADEV,RDEV,S2C2(MN,MX) =',1PE16.8,3E16.8)
          ADEV=0.
          RDEV=0.
          S2C2MN=10.
          S2C2MX=0.
            DO 1131 IRN=1,NRNA
            IXX=IXX*663608943
            THETA=0.5 + IXX*0.23283064E-09
            THETA=THETA*PI5D2
            CTHET=PI5D2-THETA
            LTHETA=SINC1*THETA+SINC0
            LCTHET=SINC1*CTHET+SINC0
            SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
            COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
            SINT=SIN(THETA)
            COST=COS(THETA)
            ASD=ABS(SINTHE-SINT)
            ACD=ABS(COSTHE-COST)
            ADEV=AMAX1(ADEV,ASD,ACD)
            IF((SINT.NE.0.0))RDEV=AMAX1(RDEV,ASD/ABS(SINT))
            IF((COST.NE.0.0))RDEV=AMAX1(RDEV,ACD/ABS(COST))
            S2C2=SINTHE**2+COSTHE**2
            S2C2MN=AMIN1(S2C2MN,S2C2)
            S2C2MX=AMAX1(S2C2MX,S2C2)
1131      CONTINUE
1132      CONTINUE
          WRITE(6,1140)NRNA
1140      FORMAT(' TEST AT ',I7,' RANDOM ANGLES IN (0,5*PI/2)')
          WRITE(6,1150)ADEV,RDEV,S2C2MN,S2C2MX
1150      FORMAT(' ADEV,RDEV,S2C2(MN,MX) =',1PE16.8,3E16.8)
        END IF
        P=1.
          DO 1161 I=1,50
          PWR2I(I)=P
          P=P/2.
1161    CONTINUE
1162    CONTINUE
      END IF
        DO 1171 J=1,NMED
1180    CONTINUE
          DO 1181 I=1,60
          IF ((IRAYLR(I).EQ.1.AND.MED(I).EQ.J)) THEN
            IRAYLM(J)=1
            GO TO 1182
          END IF
1181    CONTINUE
1182    CONTINUE
1171  CONTINUE
1172  CONTINUE
      REWIND KMPI
      OPEN(UNIT=KMPI,STATUS='OLD')
      NM=0
        DO 1191 IM=1,NMED
        LOK(IM)=0
        IF ((IRAYLM(IM).EQ.1)) THEN
          WRITE(6,1200)IM
1200      FORMAT(' RAYLEIGH OPTION REQUESTED FOR MEDIUM NUMBER',I3,/)
        END IF
1191  CONTINUE
1192  CONTINUE
1210  CONTINUE
1211    CONTINUE
1220    CONTINUE
1221      CONTINUE
          READ(KMPI,1030,END=1230)MBUF
            DO 1241 IB=1,LMDL
            IF((MBUF(IB).NE.MDLABL(IB)))GO TO 1221
1241      CONTINUE
1242      CONTINUE
1250      CONTINUE
            DO 1251 IM=1,NMED
              DO 1261 IB=1,LMDN
              IL=LMDL+IB
              IF((MBUF(IL).NE.MEDIA(IB,IM)))GO TO 1251
              IF((IB.EQ.LMDN))GO TO 1222
1261        CONTINUE
1262        CONTINUE
1251      CONTINUE
1252      CONTINUE
        GO TO 1221
1222    CONTINUE
        IF((LOK(IM).NE.0))GO TO 1220
        LOK(IM)=1
        NM=NM+1
        WRITE(KMPO,1270)IM,MBUF
1270    FORMAT(' DATA FOR MEDIUM #',I3,', WHICH IS:',72A1)
        READ(KMPI,1,ERR=1280) (MBUF(I),I=1,5),RHO(IM),NNE(IM),IUNRST(IM)
     *  ,EPSTFL(IM),IAPRIM(IM)
1       FORMAT(5A1,5X,F11.0,4X,I2,9X,I1,9X,I1,9X,I1)
        GO TO 1290
1280    BACKSPACE(KMPI)
        READ(KMPI,2)(MBUF(I),I=1,5),RHO(IM),NNE(IM),IUNRST(IM),EPSTFL(IM
     *  ), IAPRIM(IM)
2       FORMAT(5A1,5X,F11.0,4X,I2,26X,I1,9X,I1,9X,I1)
1290    CONTINUE
        WRITE(KMPO,1300)(MBUF(I),I=1,5),RHO(IM),NNE(IM)
1300    FORMAT (5A1,',RHO=',1PG11.4,',NE=',I2,',COMPOSITION IS :')
          DO 1311 IE=1,NNE(IM)
          READ(KMPI,1320)(MBUF(I),I=1,6),(ASYM(IM,IE,I),I=1,2), ZELEM(IM
     *    ,IE),WA(IM,IE),PZ(IM,IE),RHOZ(IM,IE)
1320      FORMAT (6A1,2A1,3X,F3.0,3X,F9.0,4X,F12.0,6X,F12.0)
          WRITE(KMPO,1330)(MBUF(I),I=1,6),(ASYM(IM,IE,I),I=1,2), ZELEM(I
     *    M,IE),WA(IM,IE),PZ(IM,IE),RHOZ(IM,IE)
1330      FORMAT (6A1,2A1,',Z=',F3.0,',A=',F9.3,',PZ=',1PE12.5,',RHOZ=',
     *    1PE12.5)
1311    CONTINUE
1312    CONTINUE
        WRITE(KMPO,1340)
1340    FORMAT(' ECHO READ:$LGN(RLC,AE,AP,UE,UP(IM))')
        READ(KMPI,1020)RLC(IM),AE(IM),AP(IM),UE(IM),UP(IM)
        WRITE(KMPO,1020)RLC(IM),AE(IM),AP(IM),UE(IM),UP(IM)
        TE(IM)=AE(IM)-RM
        THMOLL(IM)=TE(IM)*2. + RM
        WRITE(KMPO,1350)
1350    FORMAT(' ECHO READ:$LGN(MSGE,MGE,MSEKE,MEKE,MLEKE,MCMFP,MRANGE(I
     *M)),IRAYL')
        READ(KMPI,1010)MSGE(IM),MGE(IM),MSEKE(IM),MEKE(IM),MLEKE(IM),MCM
     *  FP(IM),MRANGE(IM),IRAYL
        WRITE(KMPO,1010)MSGE(IM),MGE(IM),MSEKE(IM),MEKE(IM),MLEKE(IM),MC
     *  MFP(IM),MRANGE(IM),IRAYL
        NSGE=MSGE(IM)
        NGE=MGE(IM)
        NSEKE=MSEKE(IM)
        NEKE=MEKE(IM)
        NLEKE=MLEKE(IM)
        NCMFP=MCMFP(IM)
        NRANGE=MRANGE(IM)
        WRITE(KMPO,1360)
1360    FORMAT(' ECHO READ:($LGN(DL(I,IM)/1,2,3,4,5,6/),I=1,6)')
        READ(KMPI,1020)(DL1(I,IM),DL2(I,IM),DL3(I,IM),DL4(I,IM),DL5(I,IM
     *  ),DL6(I,IM),I=1,6)
        WRITE(KMPO,1020)(DL1(I,IM),DL2(I,IM),DL3(I,IM),DL4(I,IM),DL5(I,I
     *  M),DL6(I,IM),I=1,6)
        WRITE(KMPO,1370)
1370    FORMAT(' ECHO READ:DELCM(IM),($LGN(ALPHI,BPAR, DELPOS(I,IM)),I=1
     *,2)')
        READ(KMPI,1020)DELCM(IM),(ALPHI(I,IM),BPAR(I,IM),DELPOS(I,IM),I=
     *  1,2)
        WRITE(KMPO,1020)DELCM(IM),(ALPHI(I,IM),BPAR(I,IM),DELPOS(I,IM),I
     *  =1,2)
        WRITE(KMPO,1380)
1380    FORMAT(' ECHO READ:$LGN(XR0,TEFF0,BLCC,XCC(IM))')
        READ(KMPI,1020)XR0(IM),TEFF0(IM),BLCC(IM),XCC(IM)
        WRITE(KMPO,1020)XR0(IM),TEFF0(IM),BLCC(IM),XCC(IM)
        WRITE(KMPO,1390)
1390    FORMAT(' ECHO READ:$LGN(EKE(IM)/0,1/)')
        READ(KMPI,1020)EKE0(IM),EKE1(IM)
        WRITE(KMPO,1020)EKE0(IM),EKE1(IM)
        WRITE(KMPO,1400)
1400    FORMAT(' ECHO READ:($LGN(ESIG,PSIG,EDEDX,PDEDX,EBR1,PBR1,PBR2, T
     *MXS(I,IM)/0,1/),I=1,NEKE)')
        READ(KMPI,1020)(ESIG0(I,IM),ESIG1(I,IM),PSIG0(I,IM),PSIG1(I,IM),
     *  EDEDX0(I,IM),EDEDX1(I,IM),PDEDX0(I,IM),PDEDX1(I,IM),EBR10(I,IM),
     *  EBR11(I,IM),PBR10(I,IM),PBR11(I,IM),PBR20(I,IM),PBR21(I,IM),TMXS
     *  0(I,IM),TMXS1(I,IM),I=1,NEKE)
        WRITE(KMPO,1020)(ESIG0(I,IM),ESIG1(I,IM),PSIG0(I,IM),PSIG1(I,IM)
     *  ,EDEDX0(I,IM),EDEDX1(I,IM),PDEDX0(I,IM),PDEDX1(I,IM),EBR10(I,IM)
     *  ,EBR11(I,IM),PBR10(I,IM),PBR11(I,IM),PBR20(I,IM),PBR21(I,IM),TMX
     *  S0(I,IM),TMXS1(I,IM),I=1,NEKE)
        WRITE(KMPO,1410)
1410    FORMAT(' ECHO READ:EBINDA(IM),$LGN(GE(IM)/0,1/)')
        READ(KMPI,1020)EBINDA(IM),GE0(IM),GE1(IM)
        WRITE(KMPO,1020)EBINDA(IM),GE0(IM),GE1(IM)
        WRITE(KMPO,1420)
1420    FORMAT(' ECHO READ:($LGN(GMFP,GBR1,GBR2(I,IM)/0,1/),I=1,NGE)')
        READ(KMPI,1020)(GMFP0(I,IM),GMFP1(I,IM),GBR10(I,IM),GBR11(I,IM),
     *  GBR20(I,IM),GBR21(I,IM),I=1,NGE)
        WRITE(KMPO,1020)(GMFP0(I,IM),GMFP1(I,IM),GBR10(I,IM),GBR11(I,IM)
     *  ,GBR20(I,IM),GBR21(I,IM),I=1,NGE)
        IF ((IRAYLM(IM).EQ.1.AND.IRAYL.NE.1)) THEN
          WRITE(6,1430)IM
1430      FORMAT(' STOPPED IN HATCH: REQUESTED RAYLEIGH OPTION FOR MEDIU
     *M',I3, /,' BUT RAYLEIGH DATA NOT INCLUDED IN DATA CREATED BY PEGS.
     *')
          STOP
        END IF
        IF ((IRAYL.EQ.1)) THEN
          WRITE(KMPO,1440)
1440      FORMAT(' ECHO READ:NGR(IM)')
          READ(KMPI,1010)NGR(IM)
          WRITE(KMPO,1010)NGR(IM)
          NGRIM=NGR(IM)
          WRITE(KMPO,1450)
1450      FORMAT(' ECHO READ:$LGN(RCO(IM)/0,1/)')
          READ(KMPI,1020)RCO0(IM),RCO1(IM)
          WRITE(KMPO,1020)RCO0(IM),RCO1(IM)
          WRITE(KMPO,1460)
1460      FORMAT(' ECHO READ:($LGN(RSCT(I,IM)/0,1/),I=1,NGRIM)')
          READ(KMPI,1020)(RSCT0(I,IM),RSCT1(I,IM),I=1,NGRIM)
          WRITE(KMPO,1020)(RSCT0(I,IM),RSCT1(I,IM),I=1,NGRIM)
          WRITE(KMPO,1470)
1470      FORMAT(' ECHO READ:($LGN(COHE(I,IM)/0,1/),I=1,NGE)')
          READ(KMPI,1020)(COHE0(I,IM),COHE1(I,IM),I=1,NGE)
          WRITE(KMPO,1020)(COHE0(I,IM),COHE1(I,IM),I=1,NGE)
          IF ((IRAYLM(IM).NE.1)) THEN
            WRITE(6,1480)IM
1480        FORMAT(' RAYLEIGH DATA AVAILABLE FOR MEDIUM',I3, ' BUT OPTIO
     *N NOT REQUESTED.',/)
          END IF
        END IF
        IF((NM.GE.NMED))GO TO1212
      GO TO 1211
1212  CONTINUE
      DUNITR=DUNIT
      IF ((DUNIT.LT.0.0)) THEN
        ID=MAX0(1,MIN0(3,IFIX(-DUNIT)))
        DUNIT=RLC(ID)
      END IF
      IF ((DUNIT.NE.1.0)) THEN
        WRITE(6,1490)DUNITR,DUNIT
1490    FORMAT(' DUNIT REQUESTED&USED ARE:',1PE14.5,E14.5,'(CM.)')
      END IF
        DO 1501 IM=1,NMED
        DFACT=RLC(IM)/DUNIT
        DFACTI=1.0/DFACT
        I=1
          GO TO 1513
1511      I=I+1
1513      IF(I-(MEKE(IM)).GT.0)GO TO 1512
          ESIG0(I,IM)=ESIG0(I,IM)*DFACTI
          ESIG1(I,IM)=ESIG1(I,IM)*DFACTI
          PSIG0(I,IM)=PSIG0(I,IM)*DFACTI
          PSIG1(I,IM)=PSIG1(I,IM)*DFACTI
          EDEDX0(I,IM)=EDEDX0(I,IM)*DFACTI
          EDEDX1(I,IM)=EDEDX1(I,IM)*DFACTI
          PDEDX0(I,IM)=PDEDX0(I,IM)*DFACTI
          PDEDX1(I,IM)=PDEDX1(I,IM)*DFACTI
          TMXS0(I,IM)=TMXS0(I,IM)*DFACT
          TMXS1(I,IM)=TMXS1(I,IM)*DFACT
        GO TO 1511
1512    CONTINUE
        I=1
          GO TO 1523
1521      I=I+1
1523      IF(I-(MLEKE(IM)).GT.0)GO TO 1522
          ERANG0(I,IM)=ERANG0(I,IM)*DFACT
          ERANG1(I,IM)=ERANG1(I,IM)*DFACT
          PRANG0(I,IM)=PRANG0(I,IM)*DFACT
          PRANG1(I,IM)=PRANG1(I,IM)*DFACT
        GO TO 1521
1522    CONTINUE
        TEFF0(IM)=TEFF0(IM)*DFACT
        BLCC(IM)=BLCC(IM)*DFACTI
        XCC(IM)=XCC(IM)*SQRT(DFACTI)
        RLDU(IM)=RLC(IM)/DUNIT
        I=1
          GO TO 1533
1531      I=I+1
1533      IF(I-(MGE(IM)).GT.0)GO TO 1532
          GMFP0(I,IM)=GMFP0(I,IM)*DFACT
          GMFP1(I,IM)=GMFP1(I,IM)*DFACT
        GO TO 1531
1532    CONTINUE
1501  CONTINUE
1502  CONTINUE
      VACDST=VACDST*DUNITO/DUNIT
      DUNITO=DUNIT
        DO 1541 JR=1,60
        MD=MED(JR)
        IF (((MD.GE.1).AND.(MD.LE.NMED))) THEN
          ECUT(JR)=AMAX1(ECUT(JR),AE(MD))
          PCUT(JR)=AMAX1(PCUT(JR),AP(MD))
          IF ((RHOR(JR).EQ.0.0)) THEN
            RHOR(JR)=RHO(MD)
          END IF
        END IF
1541  CONTINUE
1542  CONTINUE
      IF ((IBRDST.EQ.1)) THEN
          DO 1551 IM=1,NMED
          ZBRANG(IM)=0.0
          PZNORM=0.0
            DO 1561 IE=1,NNE(IM)
            ZBRANG(IM)= ZBRANG(IM)+PZ(IM,IE)*ZELEM(IM,IE)*(ZELEM(IM,IE)+
     *      1.0)
            PZNORM=PZNORM+PZ(IM,IE)
1561      CONTINUE
1562      CONTINUE
          ZBRANG(IM)=(8.116224E-05)*(ZBRANG(IM)/PZNORM)**(1./3.)
1551    CONTINUE
1552    CONTINUE
      END IF
      IF ((IPRDST.GT.0)) THEN
          DO 1571 IM=1,NMED
          ZBRANG(IM)=0.0
          PZNORM=0.0
            DO 1581 IE=1,NNE(IM)
            ZBRANG(IM)= ZBRANG(IM)+PZ(IM,IE)*ZELEM(IM,IE)*(ZELEM(IM,IE)+
     *      1.0)
            PZNORM=PZNORM+PZ(IM,IE)
1581      CONTINUE
1582      CONTINUE
          ZBRANG(IM)=(8.116224E-05)*(ZBRANG(IM)/PZNORM)**(1./3.)
1571    CONTINUE
1572    CONTINUE
      END IF
      IF ((NMED.EQ.1)) THEN
        WRITE(6,1590)
1590    FORMAT(' EGS SUCCESSFULLY ''HATCHED'' FOR ONE MEDIUM.')
      ELSE
        WRITE(6,1600)NMED
1600    FORMAT(' EGS SUCCESSFULLY ''HATCHED'' FOR ',I5,' MEDIA.')
      END IF
      CLOSE(UNIT=KMPI,STATUS='KEEP')
      RETURN
1230  WRITE(6,1610)KMPI
1610  FORMAT(' END OF FILE ON UNIT ',I2,//, ' PROGRAM STOPPED IN HATCH B
     *ECAUSE THE',/, ' FOLLOWING NAMES WERE NOT RECOGNIZED:',/)
        DO 1621 IM=1,NMED
        IF ((LOK(IM).NE.1)) THEN
          WRITE(6,1630)(MEDIA(I,IM),I=1,LMDN)
1630      FORMAT(40X,'''',24A1,'''')
        END IF
1621  CONTINUE
1622  CONTINUE
      STOP
      END
      SUBROUTINE MOLLER
      DOUBLE PRECISION PEIE,PEKSE2,PESE1,PESE2
      DOUBLE PRECISION PEKIN,H1,DCOSTH
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/RANDOM/IXX
      IRCODE=1
      PEIE=E(NP)
      EIE=PEIE
      PEKIN=PEIE-PRM
      EKIN=PEKIN
      T0=EKIN/RM
      E0=T0+1.0
      EXTRAE = EIE - THMOLL(MEDIUM)
      E02=E0*E0
      BETAI2=E02/(E02-1.0)
      EP0=TE(MEDIUM)/EKIN
      G1=(1.-2.*EP0)*BETAI2
      G2=T0*T0/E02
      G3=(2.*T0+1.)/E02
1011  CONTINUE
        IXX=IXX*663608943
        RNNO27=0.5 + IXX*0.23283064E-09
        BR = TE(MEDIUM)/(EKIN-EXTRAE*RNNO27)
        R=BR/(1.-BR)
        IXX=IXX*663608943
        RNNO28=0.5 + IXX*0.23283064E-09
        REJF4=G1*(1.+G2*BR*BR+R*(R-G3))
        IF((RNNO28.LE.REJF4))GO TO1012
      GO TO 1011
1012  CONTINUE
      PEKSE2=BR*EKIN
      PESE1=PEIE-PEKSE2
      PESE2=PEKSE2+PRM
      ESE1=PESE1
      ESE2=PESE2
      E(NP)=PESE1
      E(NP+1)=PESE2
      H1=(PEIE+PRM)/PEKIN
      DCOSTH=H1*(PESE1-PRM)/(PESE1+PRM)
      SINTHE=DSQRT(1.D0-DCOSTH)
      COSTHE=DSQRT(DCOSTH)
      CALL UPHI(2,1)
      NP=NP+1
      IQ(NP)=-1
      DCOSTH=H1*(PESE2-PRM)/(PESE2+PRM)
      SINTHE=-DSQRT(1.D0-DCOSTH)
      COSTHE=DSQRT(DCOSTH)
      CALL UPHI(3,2)
      RETURN
      END
      SUBROUTINE MSCAT
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      INTEGER EPSTFL(3)
      COMMON/ELECIN/ EKELIM, EKE0(3),EKE1(3),CMFP0(3),CMFP1(3),RANGE0(3)
     *,RANGE1(3), XR0(3),TEFF0(3),BLCC(3),XCC(3), PICMP0(1,3),PICMP1(1,3
     *),EICMP0(1,3),EICMP1(1,3), ESIG0(150,3),ESIG1(150,3),PSIG0(150,3),
     *PSIG1(150,3),EDEDX0(150,3),EDEDX1(150,3),PDEDX0(150,3),PDEDX1(150,
     *3),EBR10(150,3),EBR11(150,3),PBR10(150,3),PBR11(150,3),PBR20(150,3
     *),PBR21(150,3),TMXS0(150,3),TMXS1(150,3), CMFPE0(1,3),CMFPE1(1,3),
     *CMFPP0(1,3),CMFPP1(1,3),ERANG0(1,3),ERANG1(1,3),PRANG0(1,3),PRANG1
     *(1,3),CXC2E0(1,3),CXC2E1(1,3),CXC2P0(1,3),CXC2P1(1,3),CLXAE0(1,3),
     *CLXAE1(1,3),CLXAP0(1,3),CLXAP1(1,3), THR0(1,1),THR1(1,1),THR2(1,1)
     *, THRI0(1,1),THRI1(1,1),THRI2(1,1), FSTEP(16),FSQR(16), VERT1(1000
     *),VERT2(100,16), BLC0,BLC1,RTHR0,RTHR1,RTHRI0,RTHRI1, ICOMP, MPEEM
     *(1,3), MSMAP(200), MSTEPS,JRMAX,MXV1, MXV2,NBLC,NRNTH,NRNTHI, IUNR
     *ST(3),EPSTFL,IAPRIM(3)
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(60),RHOR(60),IRAYLR(60)
      COMMON/MULTS/ B0G21,B1G21,G210(7),G211(7),G212(7), B0G22,B1G22,G22
     *0(8),G221(8),G222(8), B0G31,B1G31,G310(11),G311(11),G312(11), B0G3
     *2,B1G32,G320(25),G321(25),G322(25), B0BGB,B1BGB,BGB0(8),BGB1(8),BG
     *B2(8), NG21,NG22,NG31,NG32,NBGB
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIIN/SINC0,SINC1,SIN0(1002),SIN1(1002)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/RANDOM/IXX
      VSTEFF=TVSTEP*RHOF
      OMEGA0=BLCC(MEDIUM)*VSTEFF/BETA2
      IF ((OMEGA0.LE.1.0)) THEN
        SINTHE=0.0
        COSTHE=1.0
        THETA=0.0
        NOSCAT=NOSCAT+1
        RETURN
      END IF
      BLC=ALOG(OMEGA0)
      IF ((BLC.LE.1.30685)) THEN
        B=1.530394*BLC
      ELSE
        IB=B0BGB+BLC*B1BGB
        IF ((IB.GT.NBGB)) THEN
          WRITE(6,1010)IB
1010      FORMAT(' NBGB<IB=',I5)
        END IF
        B=BGB0(IB)+BLC*(BGB1(IB)+BLC*BGB2(IB))
      END IF
      XR=XCC(MEDIUM)*SQRT(VSTEFF*B)/(EOLD*BETA2)
      IF ((B.GT.2.0)) THEN
        BI=1./B
      ELSE
        BI=0.5
      END IF
      BMD=1.0+1.75*BI
      BM1=(1.0-2.0/B)/BMD
      BM2=(1.0+0.025*BI)/BMD
1021  CONTINUE
        IXX=IXX*663608943
        RMS1=0.5 + IXX*0.23283064E-09
        IF ((RMS1.LE.BM1)) THEN
          IXX=IXX*663608943
          RMS2=0.5 + IXX*0.23283064E-09
          IF ((RMS2.EQ.0.0)) THEN
            RMS2=1.E-30
          END IF
          THR=SQRT(AMAX1(0.,-ALOG(RMS2)))
        ELSE IF((RMS1.LE.BM2)) THEN
          IXX=IXX*663608943
          RMS3=0.5 + IXX*0.23283064E-09
          IXX=IXX*663608943
          RMS4=0.5 + IXX*0.23283064E-09
          ETA=AMAX1(RMS3,RMS4)
          I31=B0G31+ETA*B1G31
          G31=G310(I31)+ETA*(G311(I31)+ETA*G312(I31))
          I32=B0G32+ETA*B1G32
          G32=G320(I32)+ETA*(G321(I32)+ETA*G322(I32))
          G3=G31+G32*BI
          IXX=IXX*663608943
          RMS5=0.5 + IXX*0.23283064E-09
          IF((RMS5.GT.G3))GO TO1021
          THR=1.0/ETA
        ELSE
          IXX=IXX*663608943
          RMS6=0.5 + IXX*0.23283064E-09
          THR=RMS6
          I21=B0G21+THR*B1G21
          G21=G210(I21)+THR*(G211(I21)+THR*G212(I21))
          I22=B0G22+THR*B1G22
          G22=G220(I22)+THR*(G221(I22)+THR*G222(I22))
          G2=G21+G22*BI
          IXX=IXX*663608943
          RMS7=0.5 + IXX*0.23283064E-09
          IF((RMS7.GT.G2))GO TO1021
        END IF
        THETA=THR*XR
        IF((THETA.GE.PI))GO TO1021
        LTHETA=SINC1*THETA+SINC0
        SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
        IXX=IXX*663608943
        RMS8=0.5 + IXX*0.23283064E-09
        IF(((RMS8**2*THETA.LE.SINTHE)))GO TO1022
      GO TO 1021
1022  CONTINUE
      CTHET=PI5D2-THETA
      LCTHET=SINC1*CTHET+SINC0
      COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
      RETURN
      END
      SUBROUTINE PAIR
      DOUBLE PRECISION PEIG,PESE1,PESE2
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/BREMPR/ DL1(6,3),DL2(6,3),DL3(6,3),DL4(6,3),DL5(6,3),DL6(6,
     *3), ALPHI(2,3),BPAR(2,3),DELPOS(2,3), ASYM(3,50,2), WA(3,50),PZ(3,
     *50),ZELEM(3,50),RHOZ(3,50), PWR2I(50), DELCM(3),NNE(3),ZBRANG(3),
     *IBRDST,IPRDST,IBRSPL,NBRSPL,FBRSPL
      CHARACTER*4 ASYM
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/RANDOM/IXX
      PEIG=E(NP)
      EIG=PEIG
      IF ((EIG.LE.2.1)) THEN
        ESE2=PRM
      ELSE
        IF (EIG.LT.50.) THEN
          LVX=1
          LVL0=0
        ELSE
          LVX=2
          LVL0=3
        END IF
1011    CONTINUE
          IXX=IXX*663608943
          RNNO30=0.5 + IXX*0.23283064E-09
          IXX=IXX*663608943
          RNNO31=0.5 + IXX*0.23283064E-09
          IF ((RNNO31.GE.BPAR(LVX,MEDIUM))) THEN
            LVL=LVL0+1
            IXX=IXX*663608943
            RNNO32=0.5 + IXX*0.23283064E-09
            IXX=IXX*663608943
            RNNO33=0.5 + IXX*0.23283064E-09
            BR=0.5*(1.0-AMAX1(RNNO32,RNNO33,RNNO30))
          ELSE
            LVL=LVL0+3
            BR=RNNO30*0.5
          END IF
          IF((BR.EQ.0.0))GO TO1011
          DEL=1.0/(EIG*BR*(1.0-BR))
          IF((DEL.GE.DELPOS(LVX,MEDIUM)))GO TO1011
          DELTA=DELCM(MEDIUM)*DEL
          IF ((DELTA.LT.1.0)) THEN
            REJF=DL1(LVL,MEDIUM)+DELTA*(DL2(LVL,MEDIUM)+DELTA*DL3(LVL,ME
     *      DIUM))
          ELSE
            REJF=DL4(LVL,MEDIUM)+DL5(LVL,MEDIUM)*ALOG(DELTA+DL6(LVL,MEDI
     *      UM))
          END IF
          IXX=IXX*663608943
          RNSCRN=0.5 + IXX*0.23283064E-09
          IF((((RNSCRN.LE.REJF).AND.(BR.GE.(RM/EIG)))))GO TO1012
        GO TO 1011
1012    CONTINUE
        ESE2=BR*EIG
      END IF
      PESE2=ESE2
      PESE1=PEIG-PESE2
      E(NP)=PESE1
      E(NP+1)=PESE2
      IF (((IPRDST.EQ.1).OR.((IPRDST.EQ.2).AND.(EIG.LT.4.14)))) THEN
          DO 1021 ICHRG=1,2
          IF ((ICHRG.EQ.1)) THEN
            ESE=PESE1
          ELSE
            ESE=ESE2
          END IF
          PSE=SQRT(MAX(0.0,(ESE-RM)*(ESE+RM)))
          IXX=IXX*663608943
          COSTHE=0.5 + IXX*0.23283064E-09
          COSTHE=1.0-2.0*COSTHE
          SINTHE=RM*SQRT((1.0-COSTHE)*(1.0+COSTHE))/(PSE*COSTHE+ESE)
          COSTHE=(ESE*COSTHE+PSE)/(PSE*COSTHE+ESE)
          IF ((ICHRG.EQ.1)) THEN
            CALL UPHI(2,1)
          ELSE
            NP=NP+1
            SINTHE=-SINTHE
            CALL UPHI(3,2)
          END IF
1021    CONTINUE
1022    CONTINUE
        IF ((RNNO34.LE.0.5)) THEN
          IQ(NP)=1
          IQ(NP-1)=-1
        ELSE
          IQ(NP)=-1
          IQ(NP-1)=1
        END IF
        RETURN
      ELSE IF(((IPRDST.EQ.2).AND.(EIG.GE.4.14))) THEN
        ZTARG=ZBRANG(MEDIUM)
        TTEIG=EIG/RM
          DO 1031 ICHRG=1,2
          IF ((ICHRG.EQ.1)) THEN
            ESE=PESE1
          ELSE
            ESE=ESE2
          END IF
          TTESE=ESE/RM
          TTPSE=SQRT((TTESE-1.0)*(TTESE+1.0))
          ESEDEI=TTESE/(TTEIG-TTESE)
          ESEDER=1.0/ESEDEI
          XIMIN=1.0/(1.0+(3.141593*TTESE)**2)
          REJMIN = 2.0+3.0*(ESEDEI+ESEDER) - 4.00*(ESEDEI+ESEDER+1.0-4.0
     *    *(XIMIN-0.5)**2)*( 1.0+0.25*LOG( ((1.0+ESEDER)*(1.0+ESEDEI)/(2
     *    .*TTEIG))**2+ZTARG*XIMIN**2 ) )
          YA=(2.0/TTEIG)**2
          XITRY=MAX(0.01,MAX(XIMIN,MIN(0.5,SQRT(YA/ZTARG))))
          GALPHA=1.0+0.25*LOG(YA+ZTARG*XITRY**2)
          GBETA=0.5*ZTARG*XITRY/(YA+ZTARG*XITRY**2)
          GALPHA=GALPHA-GBETA*(XITRY-0.5)
          XIMID=GALPHA/(3.0*GBETA)
          IF ((GALPHA.GE.0.0)) THEN
            XIMID=0.5-XIMID+SQRT(XIMID**2+0.25)
          ELSE
            XIMID=0.5-XIMID-SQRT(XIMID**2+0.25)
          END IF
          XIMID=MAX(0.01,MAX(XIMIN,MIN(0.5,XIMID)))
          REJMID = 2.0+3.0*(ESEDEI+ESEDER) - 4.00*(ESEDEI+ESEDER+1.0-4.0
     *    *(XIMID-0.5)**2)*( 1.0+0.25*LOG( ((1.0+ESEDER)*(1.0+ESEDEI)/(2
     *    .*TTEIG))**2+ZTARG*XIMID**2 ) )
          REJTOP=1.02*MAX(REJMIN,REJMID)
1041      CONTINUE
            IXX=IXX*663608943
            XITST=0.5 + IXX*0.23283064E-09
            REJTST = 2.0+3.0*(ESEDEI+ESEDER) - 4.00*(ESEDEI+ESEDER+1.0-4
     *      .0*(XITST-0.5)**2)*( 1.0+0.25*LOG( ((1.0+ESEDER)*(1.0+ESEDEI
     *      )/(2.*TTEIG))**2+ZTARG*XITST**2 ) )
            IXX=IXX*663608943
            RTEST=0.5 + IXX*0.23283064E-09
            THETA=SQRT(1.0/XITST-1.0)/TTESE
            IF((((RTEST.LE.(REJTST/REJTOP)).AND.(THETA.LT.PI))))GO TO104
     *      2
          GO TO 1041
1042      CONTINUE
          SINTHE=SIN(THETA)
          COSTHE=COS(THETA)
          IF ((ICHRG.EQ.1)) THEN
            CALL UPHI(2,1)
          ELSE
            NP=NP+1
            SINTHE=-SINTHE
            CALL UPHI(3,2)
          END IF
1031    CONTINUE
1032    CONTINUE
        IF ((RNNO34.LE.0.5)) THEN
          IQ(NP)=1
          IQ(NP-1)=-1
        ELSE
          IQ(NP)=-1
          IQ(NP-1)=1
        END IF
        RETURN
      ELSE
        THETA=RM/EIG
      END IF
      CALL UPHI(1,1)
      NP=NP+1
      SINTHE=-SINTHE
      CALL UPHI(3,2)
      IXX=IXX*663608943
      RNNO34=0.5 + IXX*0.23283064E-09
      IF ((RNNO34.LE.0.5)) THEN
        IQ(NP)=1
        IQ(NP-1)=-1
      ELSE
        IQ(NP)=-1
        IQ(NP-1)=1
      END IF
      RETURN
      END
      SUBROUTINE PHOTO
      DOUBLE PRECISION PEIG
      COMMON/BOUNDS/ECUT(60),PCUT(60),VACDST
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/EDGE/IEDGFL(60),EKALPH(3),EKBETA(3),BKR1(3),BKR2(3)
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/MEDIA/ RLC(3),RLDU(3),RHO(3),MSGE(3),MGE(3),MSEKE(3),MEKE(3
     *),MLEKE(3),MCMFP(3),MRANGE(3),IRAYLM(3),NMED
      COMMON/MEDIAC/MEDIA(24,3)
      CHARACTER*4 MEDIA
      COMMON/PHOTIN/ EBINDA(3), GE0(3),GE1(3), GMFP0(200,3),GMFP1(200,3)
     *,GBR10(200,3),GBR11(200,3),GBR20(200,3),GBR21(200,3), RCO0(3),RCO1
     *(3), RSCT0(100,3),RSCT1(100,3), COHE0(200,3),COHE1(200,3), MPGEM(1
     *,3), NGR(3)
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/USERSC/ SMAX,SMAXIR(60),ESTEPE,ESTEPR(60),ESAVE(60), IEDEP,
     *NOMSCT(60),NOPLC(60)
      COMMON/USERVR/ CEXPTR,GWAIT,IFORCE,NFMIN,NFMAX,NFTIME,ISOURC,IFPB,
     *IQINC,MONOEN
      COMMON/USERXT/IPHTER(60)
      COMMON/RANDOM/IXX
      PEIG=E(NP)
      IF ((E(NP).LE.EBINDA(MEDIUM))) THEN
        EDEP=PEIG
        IBLOBE=1
      ELSE
        IF ((IEDGFL(IR(NP)).NE.0)) THEN
          IXX=IXX*663608943
          BR=0.5 + IXX*0.23283064E-09
          IF ((BR.GT.BKR1(MEDIUM))) THEN
            ENEW=0.0
          ELSE
            IF ((BR.LE.BKR2(MEDIUM))) THEN
              ENEW=EKALPH(MEDIUM)
            ELSE
              ENEW=EKBETA(MEDIUM)
            END IF
          END IF
        ELSE
          ENEW=0.0
        END IF
        EDEP=EBINDA(MEDIUM)-ENEW
        E(NP)=EDEP
        IBLOBE=0
      END IF
      IARG=4
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      IF ((IBLOBE.EQ.1)) THEN
        E(NP)=PZERO
      ELSE
        IQ(NP)=-1
        E(NP)=PEIG-EBINDA(MEDIUM)+PRM
        IF ((IPHTER(IR(NP)).EQ.1)) THEN
          EELEC=E(NP)
          IF ((EELEC.GT.ECUT(IR(NP)))) THEN
            BETA=SQRT((EELEC-RM)*(EELEC+RM))/EELEC
            GAMMA=EELEC/RM
            ALPHA=0.5*GAMMA-0.5+1./GAMMA
            RATIO=BETA/ALPHA
1011        CONTINUE
              IXX=IXX*663608943
              RNPHT=0.5 + IXX*0.23283064E-09
              RNPHT=2.*RNPHT-1.
              IF ((RATIO.LE.0.2)) THEN
                FKAPPA=RNPHT+0.5*RATIO*(1.-RNPHT)*(1.+RNPHT)
                COSTHE=(BETA+FKAPPA)/(1.+BETA*FKAPPA)
                XI=1./(1.-BETA*COSTHE)
              ELSE
                XI=GAMMA*GAMMA*(1.+ALPHA*(SQRT(1.+RATIO*(2.*RNPHT+RATIO)
     *          )-1.))
                COSTHE=(1.-1./XI)/BETA
              END IF
              SINTH2=MAX(0.,(1.-COSTHE)*(1.+COSTHE))
              IXX=IXX*663608943
              RNPHT2=0.5 + IXX*0.23283064E-09
              IF(RNPHT2.LE.0.5*(1.+GAMMA)*SINTH2*XI/GAMMA)GO TO1012
            GO TO 1011
1012        CONTINUE
            SINTHE=SQRT(SINTH2)
            CALL UPHI(2,1)
          END IF
        END IF
        IF ((IEDGFL(IR(NP)).NE.0)) THEN
          IF ((BR.GT.BKR1(MEDIUM))) THEN
            RETURN
          END IF
          NP=NP+1
          E(NP)=ENEW
          IQ(NP)=0
          IXX=IXX*663608943
          RNISO=0.5 + IXX*0.23283064E-09
          COSTHE=2.0*RNISO-1.0
          SINTHE=SQRT(1.0-COSTHE*COSTHE)
          U(NP)=0.0
          V(NP)=0.0
          W(NP)=1.0
          CALL UPHI(2,1)
          X(NP)=X(NP-1)
          Y(NP)=Y(NP-1)
          Z(NP)=Z(NP-1)
          IR(NP)=IR(NP-1)
          WT(NP)=WT(NP-1)
          DNEAR(NP)=DNEAR(NP-1)
          LATCH(NP)=LATCH(NP-1)
        END IF
      END IF
      RETURN
      END
      SUBROUTINE EDGSET(NREGLO,NREGHI)
      COMMON/EDGE/IEDGFL(60),EKALPH(3),EKBETA(3),BKR1(3),BKR2(3)
      COMMON/MEDIA/ RLC(3),RLDU(3),RHO(3),MSGE(3),MGE(3),MSEKE(3),MEKE(3
     *),MLEKE(3),MCMFP(3),MRANGE(3),IRAYLM(3),NMED
      COMMON/MEDIAC/MEDIA(24,3)
      CHARACTER*4 MEDIA
      COMMON/PHOTIN/ EBINDA(3), GE0(3),GE1(3), GMFP0(200,3),GMFP1(200,3)
     *,GBR10(200,3),GBR11(200,3),GBR20(200,3),GBR21(200,3), RCO0(3),RCO1
     *(3), RSCT0(100,3),RSCT1(100,3), COHE0(200,3),COHE1(200,3), MPGEM(1
     *,3), NGR(3)
      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(60),RHOR(60),IRAYLR(60)
      DIMENSION EALF(100),EBET(100),OMEG(100),PHOTOK(100),PKA(100)
      DATA EALF/0. , 0. , 0.054, 0.109, 0.184, 0.297, 0.393, 0.524, 0.67
     *5, 0.849, 1.041, 1.255, 1.487, 1.739, 2.014, 2.307, 2.622, 2.957,
     *3.312, 3.690, 4.088, 4.508, 4.949, 5.411, 5.895, 6.400, 6.925, 7.4
     *72, 8.041, 8.631, 9.243, 9.876, 10.532,11.210,11.907,12.630,13.375
     *,14.142,14.933,15.746, 16.584,17.443,18.327,19.235,20.167,21.122,2
     *2.103,23.108, 24.138,25.192,26.272,27.378,28.510,29.667,30.851,32.
     *062, 33.297,34.564,35.858,37.179,38.528,39.906,41.313,42.750, 44.2
     *18,45.714,47.242,48.801,50.392,52.104,53.670,55.356, 57.078,58.832
     *,60.620,62.443,64.303,66.200,68.133,70.103, 72.113,74.161,76.246,7
     *8.378,80.547,82.757,85.018,87.321, 89.662,92.050,94.491,96.977,99.
     *516,102.10,104.74,107.44, 110.20,113.03,115.93,118.89/
      DATA EBET/0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0
     *. , 0. , 1.838, 2.142, 2.468, 2.817, 3.191, 3.589, 4.012, 4.459, 4
     *.931, 5.427, 5.947, 6.492, 7.059, 7.649, 8.265, 8.907, 9.572,10.26
     *3,10.983, 11.730,12.503,13.300,14.126,14.980,15.859,16.767,17.700,
     * 18.661,19.648,20.785,21.705,22.778,23.878,25.008,26.166, 27.354,2
     *8.573,29.821,31.103,32.414,33.757,35.131,36.535, 37.966,39.431,40.
     *930,42.460,44.024,45.622,47.253,48.918, 50.618,52.352,54.123,55.93
     *0,57.772,59.652,61.572,63.531, 65.529,67.564,69.642,71.759,73.919,
     *76.121,78.367,80.656, 82.991,85.370,87.796,90.273,92.794,95.365,97
     *.989,100.66, 103.39,106.17,109.01,111.90,114.85,117.86,120.93,124.
     *08, 127.29,130.58,133.96,137.41/
      DATA OMEG/0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. , 0. ,0. ,0. ,0. ,0.0357,
     *0.0470,0.0604,0.0761, 0.0942,0.115 ,0.138 ,0.163 ,0.190 ,0.219 ,0.
     *250 ,0.282 , 0.314 ,0.347 ,0.381 ,0.414 ,0.445 ,0.479 ,0.510 ,0.54
     *0 , 0.567 ,0.596 ,0.622 ,0.646 ,0.669 ,0.691 ,0.711 ,0.730 , 0.748
     * ,0.764 ,0.779 ,0.793 ,0.807 ,0.819 ,0.830 ,0.840 , 0.850 ,0.859 ,
     *0.867 ,0.875 ,0.882 ,0.889 ,0.895 ,0.901 , 0.906 ,0.911 ,0.915 ,0.
     *920 ,0.924 ,0.928 ,0.931 ,0.934 , 0.937 ,0.940 ,0.943 ,0.945 ,0.94
     *8 ,0.950 ,0.952 ,0.954 , 0.956 ,0.957 ,0.959 ,0.961 ,0.962 ,0.963
     *,0.964 ,0.966 , 0.967 ,0.968 ,0.969 ,0.970 ,0.971 ,0.972 ,0.972 ,0
     *.973 , 0.974 ,0.975 ,0.975 ,0.976 ,0.977 ,0.977 ,0.978 ,0.978 , 0.
     *979 ,0.979 ,0.980 ,0.980 /
      DATA PHOTOK/1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 0.92
     *8, 0.922, 0.916, 0.911, 0.907, 0.903, 0.900, 0.896, 0.893, 0.890,
     *0.888, 0.885, 0.883, 0.880, 0.878, 0.876, 0.874, 0.872, 0.871, 0.8
     *70, 0.869, 0.868, 0.866, 0.865, 0.863, 0.861, 0.860, 0.858, 0.856,
     * 0.855, 0.853, 0.852, 0.850, 0.848, 0.846, 0.845, 0.843, 0.841, 0.
     *840, 0.839, 0.838, 0.836, 0.835, 0.833, 0.832, 0.830, 0.829, 0.827
     *, 0.825, 0.824, 0.823, 0.821, 0.819, 0.818, 0.816, 0.815, 0.814, 0
     *.813, 0.812, 0.811, 0.809, 0.807, 0.805, 0.803, 0.802, 0.800, 0.79
     *8, 0.796, 0.794, 0.792, 0.790, 0.788, 0.786, 0.784, 0.782, 0.780,
     *0.778, 0.776, 0.774, 0.772, 0.770, 0.767, 0.765, 0.762, 0.759, 0.7
     *56, 0.754, 0.751, 0.749, 0.746/
      DATA PKA/1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1.
     * , 1. , 1. , 1. , 0.955, 0.940, 0.924, 0.914, 0.904, 0.901, 0.898,
     * 0.897, 0.895, 0.894, 0.893, 0.8925,0.892, 0.891, 0.890, 0.887, 0.
     *884, 0.879, 0.875, 0.870, 0.866, 0.862, 0.859, 0.855, 0.852, 0.849
     *, 0.847, 0.844, 0.842, 0.840, 0.838, 0.835, 0.833, 0.831, 0.829, 0
     *.826, 0.824, 0.821, 0.819, 0.817, 0.815, 0.813, 0.812, 0.811, 0.81
     *0, 0.809, 0.808, 0.807, 0.806, 0.805, 0.804, 0.803, 0.802, 0.801,
     *0.800, 0.799, 0.798, 0.797, 0.795, 0.794, 0.794, 0.793, 0.792, 0.7
     *91, 0.790, 0.789, 0.788, 0.787, 0.786, 0.785, 0.784, 0.783, 0.782,
     * 0.781, 0.781, 0.780, 0.779, 0.778, 0.778, 0.777, 0.776, 0.775, 0.
     *774, 0.773, 0.773/
      WRITE(6,1010)
1010  FORMAT(' OUTPUT FROM SUBROUTINE EDGSET:'/ ' REGION MEDIUM  Zeff  N
     *AME                    ', '     K-EDGE    EKALPH    EKBETA      BK
     *R1      BKR2'/)
        DO 1021 JJ=NREGLO,NREGHI
        IZ=IEDGFL(JJ)
        IMED=MED(JJ)
        IF (IMED.GT.0.) THEN
          IF ((IZ.GT.0 .AND. IZ.LE.100)) THEN
            BKR1(IMED)=OMEG(IZ)*PHOTOK(IZ)
            BKR2(IMED)=BKR1(IMED)*PKA(IZ)
            EKALPH(IMED)=EALF(IZ)*1.E-3
            EKBETA(IMED)=EBET(IZ)*1.E-3
            WRITE(6,1030)JJ,IMED,IZ,(MEDIA(I,IMED),I=1,24),EBINDA(IMED),
     *      EKALPH(IMED), EKBETA(IMED),BKR1(IMED),BKR2(IMED)
1030        FORMAT(3(1X,I6),1X,24A1,1X,1P5E10.2)
          ELSE
            WRITE(6,1040)JJ,IMED,IZ,(MEDIA(I,IMED),I=1,24)
1040        FORMAT(3(1X,I6),1X,24A1,' -- K-EDGE FLUORESCENCE OPTION NOT
     *ENABLED', ' FOR THIS REGION --')
          END IF
        ELSE
          WRITE(6,1050)JJ,IMED
1050      FORMAT(2(1X,I6),' VACUUM')
        END IF
1021  CONTINUE
1022  CONTINUE
      RETURN
      END
      SUBROUTINE PHOTON(IRCODE)
      DOUBLE PRECISION PEIG
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/BOUNDS/ECUT(60),PCUT(60),VACDST
      COMMON/MEDIA/ RLC(3),RLDU(3),RHO(3),MSGE(3),MGE(3),MSEKE(3),MEKE(3
     *),MLEKE(3),MCMFP(3),MRANGE(3),IRAYLM(3),NMED
      COMMON/MEDIAC/MEDIA(24,3)
      CHARACTER*4 MEDIA
      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(60),RHOR(60),IRAYLR(60)
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/PHOTIN/ EBINDA(3), GE0(3),GE1(3), GMFP0(200,3),GMFP1(200,3)
     *,GBR10(200,3),GBR11(200,3),GBR20(200,3),GBR21(200,3), RCO0(3),RCO1
     *(3), RSCT0(100,3),RSCT1(100,3), COHE0(200,3),COHE1(200,3), MPGEM(1
     *,3), NGR(3)
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP(3),AE(3),UP(3),UE(3),TE(3),THMOLL
     *(3)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/USEFUL/PZERO,PRM,PRMT2,RM,MEDIUM,MEDOLD,IBLOBE
      DOUBLE PRECISION PZERO,PRM,PRMT2
      COMMON/USERSC/ SMAX,SMAXIR(60),ESTEPE,ESTEPR(60),ESAVE(60), IEDEP,
     *NOMSCT(60),NOPLC(60)
      COMMON/USERVR/ CEXPTR,GWAIT,IFORCE,NFMIN,NFMAX,NFTIME,ISOURC,IFPB,
     *IQINC,MONOEN
      COMMON/USERXT/IPHTER(60)
      COMMON/RANDOM/IXX
      IRCODE=1
      PEIG=E(NP)
      EIG=PEIG
      IRL=IR(NP)
      MEDIUM=MED(IRL)
      IF ((EIG.LE.PCUT(IRL))) THEN
        GO TO 1010
      END IF
1020  CONTINUE
1021    CONTINUE
        GLE=ALOG(EIG)
        IXX=IXX*663608943
        RNNO35=0.5 + IXX*0.23283064E-09
        IF ((RNNO35 .EQ. 0.0)) THEN
          RNNO35=1.E-30
        END IF
        DPMFP=-ALOG(RNNO35)
        IF ((CEXPTR.NE.0.0)) THEN
          IF ((W(NP).GT.0.0)) THEN
            TEMP=CEXPTR*W(NP)
            BEXPTR=1./(1.-TEMP)
            DPMFP=DPMFP*BEXPTR
            WT(NP)=WT(NP)*BEXPTR*EXP(-DPMFP*TEMP)
          END IF
        END IF
        IROLD=IR(NP)
1030    CONTINUE
1031      CONTINUE
          IF ((MEDIUM.NE.0)) THEN
            LGLE=GE1(MEDIUM)*GLE+GE0(MEDIUM)
            GMFPR0=GMFP1(LGLE,MEDIUM)*GLE+GMFP0(LGLE,MEDIUM)
          END IF
1040      CONTINUE
1041        CONTINUE
            IF ((MEDIUM.EQ.0)) THEN
              TSTEP=VACDST
            ELSE
              RHOF=RHOR(IRL)/RHO(MEDIUM)
              GMFP=GMFPR0/RHOF
              IF ((IRAYLR(IRL).EQ.1)) THEN
                COHFAC=COHE1(LGLE,MEDIUM)*GLE+COHE0(LGLE,MEDIUM)
                GMFP=GMFP*COHFAC
              END IF
              TSTEP=GMFP*DPMFP
            END IF
            IRNEW=IR(NP)
            IDISC=0
            USTEP=TSTEP
            TUSTEP=USTEP
            IF ((USTEP.GT.DNEAR(NP))) THEN
              CALL HOWFAR
            END IF
            IF ((IDISC.GT.0)) THEN
              GO TO 1050
            END IF
            VSTEP=USTEP
            TVSTEP=VSTEP
            EDEP=PZERO
            IARG=0
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
            X(NP)=X(NP)+U(NP)*USTEP
            Y(NP)=Y(NP)+V(NP)*USTEP
            Z(NP)=Z(NP)+W(NP)*USTEP
            DNEAR(NP)=DNEAR(NP)-USTEP
            IF ((MEDIUM.NE.0)) THEN
              DPMFP=AMAX1(0.,DPMFP-USTEP/GMFP)
            END IF
            IROLD=IR(NP)
            MEDOLD=MEDIUM
            IF ((IRNEW.NE.IROLD)) THEN
              IR(NP)=IRNEW
              IRL=IRNEW
              MEDIUM=MED(IRL)
              IF ((EIG.LE.PCUT(IRL))) THEN
                GO TO 1010
              END IF
            END IF
            IARG=5
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
            IF((IDISC.LT.0))GO TO 1050
            IF((MEDIUM.NE.MEDOLD))GO TO 1042
            IF ((MEDIUM.NE.0.AND.DPMFP.LE.1.E-6)) THEN
              GO TO 1032
            END IF
          GO TO 1041
1042      CONTINUE
        GO TO 1031
1032    CONTINUE
        IF ((IRAYLR(IRL).EQ.1)) THEN
          IXX=IXX*663608943
          RNNO37=0.5 + IXX*0.23283064E-09
          IF ((RNNO37.LE.(1.0-COHFAC))) THEN
            IARG=23
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
1060        CONTINUE
1061          CONTINUE
              IXX=IXX*663608943
              XXX=0.5 + IXX*0.23283064E-09
              LXXX=RCO1(MEDIUM)*XXX+RCO0(MEDIUM)
              X2=RSCT1(LXXX,MEDIUM)*XXX+RSCT0(LXXX,MEDIUM)
              Q2=X2*RMSQ/(20.60744*20.60744)
              COSTHE=1.-Q2/(2.*E(NP)*E(NP))
              IF((ABS(COSTHE).GT.1.0))GO TO 1060
              CSQTHE=COSTHE*COSTHE
              REJF=(1.0+CSQTHE)/2.0
              IXX=IXX*663608943
              RNNORJ=0.5 + IXX*0.23283064E-09
              IF((RNNORJ.LE.REJF))GO TO1062
            GO TO 1061
1062        CONTINUE
            SINTHE=SQRT(1.0-CSQTHE)
            CALL UPHI(2,1)
            IARG=24
            IF ((IAUSFL(IARG+1).NE.0)) THEN
              CALL AUSGAB(IARG)
            END IF
            GOTO 1020
          END IF
        END IF
        IXX=IXX*663608943
        RNNO36=0.5 + IXX*0.23283064E-09
        GBR1=GBR11(LGLE,MEDIUM)*GLE+GBR10(LGLE,MEDIUM)
        IF (((RNNO36.LE.GBR1).AND.(E(NP).GT.RMT2) )) THEN
          IARG=15
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          CALL PAIR
          IARG=16
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          GO TO 1022
        END IF
        GBR2=GBR21(LGLE,MEDIUM)*GLE+GBR20(LGLE,MEDIUM)
        IF ((RNNO36.LT.GBR2)) THEN
          IARG=17
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          CALL COMPT
          IF ((IQ(NP).EQ.0 .AND. E(NP-1).LT.ECUT(IR(NP-1)))) THEN
            IQ(NP)=IQ(NP-1)
            IQ(NP-1)=0
            T=E(NP)
            E(NP)=E(NP-1)
            E(NP-1)=T
            T=U(NP)
            U(NP)=U(NP-1)
            U(NP-1)=T
            T=V(NP)
            V(NP)=V(NP-1)
            V(NP-1)=T
            T=W(NP)
            W(NP)=W(NP-1)
            W(NP-1)=T
          END IF
          IARG=18
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          IF((IQ(NP).NE.0))GO TO 1022
        ELSE
          IARG=19
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          CALL PHOTO
          IF ((NP.EQ.0)) THEN
            IRCODE=2
            RETURN
          END IF
          IARG=20
          IF ((IAUSFL(IARG+1).NE.0)) THEN
            CALL AUSGAB(IARG)
          END IF
          IF((IQ(NP).EQ.-1))GO TO 1022
        END IF
        PEIG=E(NP)
        EIG=PEIG
        IF((EIG.LT.PCUT(IRL)))GO TO 1010
      GO TO 1021
1022  CONTINUE
      RETURN
1010  IF ((EIG.GT.AP(MEDIUM))) THEN
        IDR=1
      ELSE
        IDR=2
      END IF
      EDEP=PEIG
      IARG=IDR
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      IRCODE=2
      NP=NP-1
      RETURN
1050  EDEP=PEIG
      IARG=3
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      IRCODE=2
      NP=NP-1
      RETURN
      END
      SUBROUTINE SHOWER(IQI,EI,XI,YI,ZI,UI,VI,WI,IRI,WTI)
      DOUBLE PRECISION DEG,DPGL,DEI,DPI,DCSTH,DCOSTH,PI0MSQ
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/RANDOM/IXX
      DATA PI0MSQ/1.8215416D4/
      NP=1
      DNEARI=0.0
      IQ(1)=IQI
      E(1)=EI
      U(1)=UI
      V(1)=VI
      W(1)=WI
      X(1)=XI
      Y(1)=YI
      Z(1)=ZI
      IR(1)=IRI
      WT(1)=WTI
      DNEAR(1)=DNEARI
      LATCH(1)=LATCHI
      IF ((IQI.EQ.2)) THEN
        IF ((EI.LE.PI0MSQ)) THEN
          WRITE(6,1010)EI
1010      FORMAT(//,' STOPPED IN SUBROUTINE SHOWER---PI-ZERO OPTION INVO
     *KED', /,' BUT THE TOTAL ENERGY WAS TOO SMALL (EI=',G15.5,' MEV)')
          STOP
        END IF
        IXX=IXX*663608943
        CSTH=0.5 + IXX*0.23283064E-09
        DCSTH=CSTH
        DEI=EI
        DPI=DSQRT(DEI*DEI-PI0MSQ)
        DEG=DEI+DPI*DCSTH
        DPGL=DPI+DEI*DCSTH
        DCOSTH=DPGL/DEG
        COSTHE=DCOSTH
        SINTHE=DSQRT(1.D0-DCOSTH*DCOSTH)
        IQ(1)=0
        E(1)=DEG/2.
        CALL UPHI(2,1)
        NP=2
        DEG=DEI-DPI*DCSTH
        DPGL=DPI-DEI*DCSTH
        DCOSTH=DPGL/DEG
        COSTHE=DCOSTH
        SINTHE=-DSQRT(1.D0-DCOSTH*DCOSTH)
        IQ(2)=0
        E(2)=DEG/2.
        CALL UPHI(3,2)
      END IF
1020  CONTINUE
1021    CONTINUE
        IF((IQ(NP).EQ.0))GO TO 1030
1041    CONTINUE
1050      CALL ELECTR(IRCODE)
          IF((IRCODE.EQ.2))GO TO1042
1030      CALL PHOTON(IRCODE)
          IF((IRCODE.EQ.2))GO TO1042
        GO TO 1041
1042    CONTINUE
        IF((NP.LE.0))GO TO1022
      GO TO 1021
1022  CONTINUE
      RETURN
      END
      SUBROUTINE UPHI(IENTRY,LVL)
      COMMON/QDEBUG/QDEBUG
      LOGICAL QDEBUG
      COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP, RHOF,EOLD,ENEW
     *,EKE,ELKE,BETA2,GLE,TSCAT, IDISC,IROLD,IRNEW,IAUSFL(25)
      DOUBLE PRECISION EDEP
      COMMON/STACK/E(125),X(125),Y(125),Z(125),U(125),V(125),W(125),DNEA
     *R(125),WT(125),IQ(125),IR(125),LATCH(125),LATCHI, NP
      DOUBLE PRECISION E
      COMMON/UPHIIN/SINC0,SINC1,SIN0(1002),SIN1(1002)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      COMMON/RANDOM/IXX
      IARG=21
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      GO TO (1010,1020,1030),IENTRY
      GO TO 1040
1010  CONTINUE
      LTHETA=SINC1*THETA+SINC0
      SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
      CTHET=PI5D2-THETA
      LCTHET=SINC1*CTHET+SINC0
      COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
1020  IXX=IXX*663608943
      RNNO38=0.5 + IXX*0.23283064E-09
      PHI=RNNO38*TWOPI
      LPHI=SINC1*PHI+SINC0
      SINPHI=SIN1(LPHI)*PHI+SIN0(LPHI)
      CPHI=PI5D2-PHI
      LCPHI=SINC1*CPHI+SINC0
      COSPHI=SIN1(LCPHI)*CPHI+SIN0(LCPHI)
1030  GO TO (1050,1060,1070),LVL
      GO TO 1040
1050  A=U(NP)
      B=V(NP)
      C=W(NP)
      GO TO 1080
1070  A=U(NP-1)
      B=V(NP-1)
      C=W(NP-1)
1060  X(NP)=X(NP-1)
      Y(NP)=Y(NP-1)
      Z(NP)=Z(NP-1)
      IR(NP)=IR(NP-1)
      WT(NP)=WT(NP-1)
      DNEAR(NP)=DNEAR(NP-1)
      LATCH(NP)=LATCH(NP-1)
1080  SINPS2=A*A+B*B
      IF ((SINPS2.LT.1.0E-20)) THEN
        U(NP)=SINTHE*COSPHI
        V(NP)=SINTHE*SINPHI
        W(NP)=C*COSTHE
      ELSE
        SINPSI=SQRT(SINPS2)
        US=SINTHE*COSPHI
        VS=SINTHE*SINPHI
        SINDEL=B/SINPSI
        COSDEL=A/SINPSI
        U(NP)=C*COSDEL*US-SINDEL*VS+A*COSTHE
        V(NP)=C*SINDEL*US+COSDEL*VS+B*COSTHE
        W(NP)=-SINPSI*US+C*COSTHE
      END IF
      IARG=22
      IF ((IAUSFL(IARG+1).NE.0)) THEN
        CALL AUSGAB(IARG)
      END IF
      RETURN
1040  WRITE(6,1090)IENTRY,LVL
1090  FORMAT(' STOPPED IN UPHI WITH IENTRY,LVL=',2I6)
      STOP
      END
