*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:58  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
*-- Author :
      FUNCTION FISSIO(EK1)
C
C *** GENERATION OF PHOTONS AND NEUTRONS BY FISSION ***
C *** NVE 04-MAR-1988 CERN GENEVA ***
C
C ORIGIN : H.FESEFELDT (21-MAR-1987)
C
C THE PHYSICS IS BASED ON U(238)
C FOR OTHER MATERIALS EXTRAPOLATIONS ARE USED
C
#include "geant321/s_defcom.inc"
#include "geant321/s_kginit.inc"
C
      DIMENSION SPNEUT(10)
      DIMENSION RNDM(2)
      SAVE SPNEUT
      DATA SPNEUT/10*0./
C
C --- INITIALIZATION INDICATED BY KGINIT(15) ---
      IF (KGINIT(15) .NE. 0) GO TO 10
      KGINIT(15)=1
C
      XX=1.-0.5
      XXX=SQRT(2.29*XX)
      SPNEUT(1)=EXP(-XX/0.965)*(EXP(XXX)-EXP(-XXX))/2.
      DO 1 I=2,10
      XX=I*1.-0.5
      XXX=SQRT(2.29*XX)
    1 SPNEUT(I)=SPNEUT(I-1)+EXP(-XX/0.965)*(EXP(XXX)-EXP(-XXX))/2.
      DO 2 I=1,10
    2 SPNEUT(I)=SPNEUT(I)/SPNEUT(10)
C** IN THIS ROUTINE WE USE MEV AS UNIT FOR ENERGY AND MOMENTUM
   10 NT=0
      IER(82)=IER(82)+1
      ND=IND+1
      PV( 1,MXGKPV)=PX*P
      PV( 2,MXGKPV)=PY*P
      PV( 3,MXGKPV)=PZ*P
      PV( 4,MXGKPV)=EN
      PV( 5,MXGKPV)=ABS(AMAS)
      PV( 6,MXGKPV)=NCH
      PV( 7,MXGKPV)=TOF
      PV( 8,MXGKPV)=IPART
      PV( 9,MXGKPV)=0.
      PV(10,MXGKPV)=USERW
      PV( 1,MXGKPV-1)=0.
      PV( 2,MXGKPV-1)=0.
      PV( 3,MXGKPV-1)=0.
      PV( 4,MXGKPV-1)=ATOMAS(ATNO(ND),ZNO(ND))
      PV( 5,MXGKPV-1)=PV(4,MXGKPV-1)
      PV( 6,MXGKPV-1)=ZNO(ND)
      PV( 7,MXGKPV-1)=TOF
      PV( 8,MXGKPV-1)=0.
      PV( 9,MXGKPV-1)=0.
      PV(10,MXGKPV-1)=0.
      CALL ADD(MXGKPV,MXGKPV-1,MXGKPV-2)
      PV(1,MXGKPV-2)=-PV(1,MXGKPV-2)
      PV(2,MXGKPV-2)=-PV(2,MXGKPV-2)
      PV(3,MXGKPV-2)=-PV(3,MXGKPV-2)
C** NUMBER OF NEUTRONS AND PHOTONS
      FISSIO=0.
      E1=EK1*1000.
      IF(E1.LT.1.0) E1=1.0
      AVERN=2.569+0.559*LOG(E1)
C**   TAKE THE FOLLOWING VALUE IF PHOTOFISSION IS NOT INCLUDED
      IF(IFIX(PARMAT(IND+1,8)).EQ.0)
     *AVERN=2.569+0.900*LOG(E1)
      AVERG=9.500+0.600*LOG(E1)
      CALL NORMAL(RAN)
      NN=IFIX(AVERN+RAN*1.23+0.5)
      CALL NORMAL(RAN)
      NG=IFIX(AVERG+RAN*3.+0.5)
      IF(NN.LT.1) NN=1
      IF(NG.LT.1) NG=1
      EXN=0.
      EXG=0.
C** DISTRIBUTE KINETIC ENERGY
      DO 15 I=1,NN
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      DO 11 J=1,10
      IF(RAN.LT.SPNEUT(J)) GOTO 12
   11 CONTINUE
      J=10
   12 CALL GRNDM(RNDM,1)
      EKIN=(J-1)*1.+RNDM(1)
      EXN=EXN+EKIN
      PV(4,I)=EKIN+RMASS(16)*1000.
      PV(5,I)=RMASS(16)*1000.
      PV(6,I)=0.
C** EMISSION TIME FOR NEUTRONS =0.
      PV(7,I)=TOF
      PV(8,I)=16.
      PV(9,I)=0.
      PV(10,I)=0.
   15 CONTINUE
      NT=NN
      DO 20 I=1,NG
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      NT=NT+1
      PV(4,NT)=-0.87*LOG(RAN)
      EXG=EXG+PV(4,NT)
      PV(5,NT)=0.
      PV(6,NT)=0.
C     RAN=RNDM(1)
C** EMISSION TIME FOR PHOTONS= 2.5 E-8 SEC
C     PV(7,NT)=TOF-500.*LOG(RAN)
C** CHANGED 30.7.85
      PV(7,NT)=TOF
      PV(8,NT)=1.
      PV(9,NT)=0.
      PV(10,NT)=0.
   20 CONTINUE
      IF(NT.EQ.0) GO TO 9999
      EX=EXN+EXG
      IF(NPRT(4))
     *WRITE(NEWBCD,2000) ATNO(IND+1),NN,NG,EX
      FISSIO=EX/1000.
      DO 49 I=1,NT
      PV(5,I)=PV(5,I)/1000.
      PV(4,I)=PV(4,I)/1000.
      CALL GRNDM(RNDM,2)
      COST=-1.+2.*RNDM(1)
      SINT=SQRT(ABS(1.-COST*COST))
      PHI=RNDM(2)*TWPI
      PP=SQRT(ABS(PV(4,I)**2-PV(5,I)**2))
      PV(1,I)=PP*SINT*SIN(PHI)
      PV(2,I)=PP*SINT*COS(PHI)
      PV(3,I)=PP*COST
      CALL LOR(I,MXGKPV-2,I)
   49 CONTINUE
      INTCT=INTCT+1.
      DO 50 I=1,NT
      IF(NTOT.LT.NSIZE/12) GOTO 43
      IER(39)=IER(39)+1
      GO TO 9999
   43 CALL SETTRK(I)
   50 CONTINUE
C
 2000 FORMAT(1H ,'NUCLEAR FISSION ON MATERIAL ',F6.1,
     *', NEUTRONS, PHOTONS PRODUCED= ',2I3,' WITH ',F8.4,' MEV TOTAL ENE
     *RGY')
C
 9999 CONTINUE
      END