*
* $Id$
*
* $Log$
* Revision 1.1  1995/10/24 10:21:05  cernlib
* Initial revision
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.40  by  S.Giani
*-- Author :
      SUBROUTINE PHPNUC
C
C *** DOUBLE PRECISION VERSION OF THE PHASE SPACE ROUTINE "PHASP"
C *** THIS ROUTINE MUST BE CALLED BY THE NUCLEAR INTERACTION ROUTINE
C *** "NUCREC" (SEE ALSO COMMENTS THEREIN). THE REASON IS SIMPLY THAT
C *** ENERGY-MOMENTUM CALCULATIONS ARE NOT POSSIBLE WITHIN ONLY
C *** 6 DIGITS OF ACCURACY FOR TOTAL ENERGIES
C *** IN THE ORDER OF HUNDREDS OF GEV (URANIUM NUCLEUS), COMPARED WITH
C *** KINETIC ENERGIES IN THE ORDER OF MEV (NEUTRONS, PROTONS AND
C *** PHOTONS IN THE REACTIONS A(X,Y(GAMMA,GAMMA))A'). IN THE ORIGINAL
C *** GHEISHA8 CODE ALL THESE CALCULATIONS ARE DONE IN DOUBLE PRECISION
C *** HMF 29-AUG-1989 RWTH AACHEN
C
C CALLED BY : NUCREC
C ORIGIN    : H.FESEFELDT
C
#if !defined(CERNLIB_SINGLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#endif
      REAL RNDM(1)
C
#include "geant321/s_prntfl.inc"
#include "geant321/s_nucio.inc"
C
      SAVE  KNT, TWOPI, FFQ
      DIMENSION EMM(18)
      DIMENSION RNO(50)
      DIMENSION EM(18),PD(18),EMS(18),SM(18),FFQ(18),PCM1(90)
      EQUIVALENCE (NT,NPG),(AMASS(1),EM(1)),(PCM1(1),PCM(1,1))
      DATA  FFQ/0.,3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
     2                       256.3704, 268.4705, 240.9780, 189.2637,
     3                       132.1308,  83.0202,  47.4210,  24.8295,
     4                        12.0006,   5.3858,   2.2560,   0.8859/
      DATA  KNT , TWOPI /  1 , 6.2831853073 /
C
C --- Initialise local arrays and the result array PCM ---
      CALL VZERO(PCM,90)
      DO 80 JZERO=1,18
         EMM(JZERO)=0.
         PD(JZERO) =0.
         EMS(JZERO)=0.
         SM(JZERO) =0.
  80  CONTINUE
C
      KNT = KNT + 1
      IF (.NOT.NPRT(3).AND..NOT.NPRT(4)) GOTO 100
      WRITE(NEWBCD,1200) NPG,TECM,(AMASS(JK),JK=1,NPG)
  100 CONTINUE
  150 IF (NT .LT. 2)  GO TO 1001
      IF (NT .GT. 18)  GO TO 1002
      NTM1=NT-1
      NTM2=NT-2
      NTP1=NT+1
      NTNM4 = 3*NT - 4
      EMM(1)=EM(1)
      TM=0.0
      DO 200 I=1,NT
      EMS(I)=EM(I)**2
      TM=TM+EM(I)
 200  SM(I)=TM
      WGT=1.
 210  TECMTM=TECM-TM
      IF (TECMTM .LE. 0.0)  GO TO 1000
      EMM(NT)=TECM
      IF (KGENEV.GT.1) GO TO 400
      EMMAX=TECMTM+EM(1)
      EMMIN=0.0
      WTMAX=1.0
      DO 350 I=2,NT
      EMMIN=EMMIN+EM(I-1)
      EMMAX=EMMAX+EM(I)
  350 WTMAX=WTMAX*DPDNUC(EMMAX,EMMIN,EM(I))
      WTMAXQ=1.0/WTMAX
      GO TO 455
  400 WTMAXQ=TECMTM**NTM2*FFQ(NT) / TECM
  455 CONTINUE
      DO 457 I= 1, NTNM4
      CALL GRNDM(RNDM,1)
#if defined(CERNLIB_SINGLE)
  457 RNO(I) = RNDM(1)
#endif
#if !defined(CERNLIB_SINGLE)
  457 RNO(I) = DBLE(RNDM(1))
#endif
      IF(NTM2) 900,509,460
  460 CONTINUE
      CALL DLPNUC(RNO,NTM2)
      DO 508 J=2,NTM1
  508 EMM(J)=RNO(J-1)*(TECMTM)+SM(J)
  509 WGT=WTMAXQ
      IR=NTM2
      DO 530 I=1,NTM1
      PD(I)=DPDNUC(EMM(I+1),EMM(I),EM(I+1))
  530 WGT=WGT*PD(I)
      PCM(1,1)=0.0
      PCM(2,1)=PD(1)
      PCM(3,1)=0.0
      DO 570 I=2,NT
      PCM(1,I)=0.0
      PCM(2,I) = -PD(I-1)
      PCM(3,I)=0.0
      IR=IR+1
      BANG=TWOPI*RNO(IR)
      CB=COS(BANG)
      SB=SIN(BANG)
      IR=IR+1
      C=2.0*RNO(IR)-1.0
      S=SQRT(1.0-C*C)
      IF(I.EQ.NT) GO TO 1567
      ESYS=SQRT(PD(I)**2+EMM(I)**2)
      BETA=PD(I)/ESYS
      GAMA=ESYS/EMM(I)
      DO 568 J=1,I
      NDX = 5*J - 5
      AA= PCM1(NDX+1)**2 + PCM1(NDX+2)**2 + PCM1(NDX+3)**2
      PCM1(NDX+5) = SQRT(AA)
      PCM1(NDX+4) = SQRT(AA+EMS(J))
      CALL DOTNUC(C,S,CB,SB,PCM,J)
      PSAVE = GAMA*(PCM(2,J)+BETA*PCM(4,J))
  568 PCM(2,J)=PSAVE
      GO TO 570
 1567 DO 1568 J=1,I
      AA=PCM(1,J)**2 + PCM(2,J)**2 + PCM(3,J)**2
      PCM(5,J)=SQRT(AA)
      PCM(4,J)=SQRT(AA+EMS(J))
      CALL DOTNUC(C,S,CB,SB,PCM,J)
 1568 CONTINUE
  570 CONTINUE
  900 CONTINUE
      RETURN
 1000 DO 212 I=1,NPG
      PCM(1,I)=0.
      PCM(2,I)=0.
      PCM(3,I)=0.
      PCM(4,I)=AMASS(I)
  212 PCM(5,I)=AMASS(I)
      WGT=0.
      RETURN
 1001 IF(NPRT(3).OR.NPRT(4)) WRITE(NEWBCD,1101)
      GO TO 1050
 1002 WRITE(NEWBCD,1102)
 1050 WRITE(NEWBCD,1150) KNT
      WRITE(NEWBCD,1200) NPG,TECM,(AMASS(JK),JK=1,NPG)
      STOP
 1100 FORMAT ('0*PHPNUC* AVAILABLE ENERGY NEGATIVE')
 1101 FORMAT ('0*PHPNUC* LESS THAN 2 OUTGOING PARTICLES')
 1102 FORMAT ('0*PHPNUC* MORE THAN 18 OUTGOING PARTICLES')
 1150 FORMAT ('0*PHPNUC* ABOVE ERROR DETECTED IN PHASP',
     $ ' AT CALL NUMBER ',I7)
 1200 FORMAT ('0*PHPNUC* INPUT DATA TO PHPNUC. NPG = ',I6/
     $ ' TECM = ',E15.7,' PARTICLE MASSES = ',5E15.7/(42X,5E15.7))
      END
