C*****This subroutine fills the common NeutronsGR: Number_N, Neut_En(3),
C     teta_n(3), fi_n(3),
C*****where
C*****Number_N - number of neutrons in case of Gigantic resonance
C     (Number_N=0, if no);
C*****NeutE_n(i) - energy of each generated neutron, MeV
C*****teta_n(i) - teta (polar) angle in reference coordinate system
C*****fi_n(i) - fi (azimuth) angle in reference coordinate system Input
C     parameters from GEANT: Step - GEANT step, cm En-gamma - gamma
C     energy in tungsten, GeV px_gamma, py_gamma, pz_gamma - gamma
C     momentum projections, GeV
      subroutine GPGRES(Step, En_g, px_g, py_g, pz_g)
      
      integer i
      real En_exc
      real cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n
      
C     Input parameters
*     Energy of gamma, in GeV
      real En_gamma                     
*     Gamma momentum projections, multiplied by c, in GeV
      real px_gamma, py_gamma, pz_gamma 
*     GEANT step, in cm   
      real Step                         

      common/bind/En_bind(3)
#include "gpgneut.inc"
c     
*     
*     Let's trasnform energy and momentum in MeV. Remember, GEANT is GeV
*     
      open(77,file='testgig.dat',access='append')
      En_gamma=En_g/1000.
      px_gamma=px_g/1000.
      py_gamma=py_g/1000.
      pz_gamma=pz_g/1000.
*     
*     Exciting energy of nucleus in Gigantic resonance, MeV
      En_exc=En_gamma                   
      Number_N=0
      do i=1,3
         Neut_En(i)=0.
         teta_n(i)=0.
         fi_n(i)=0.
      enddo
      
C     Simulate number of neutrons in Gygantic resonance process
C     (Number_n=0 if no interaction)
      call GRESsimulate(En_gamma, Step, Number_N) 
      
      if (Number_N.gt.0) then
         cos_teta=pz_gamma/sqrt(px_gamma**2+py_gamma**2+pz_gamma**2)
         sin_teta=sqrt(1.-cos_teta**2)
         cos_fi=px_gamma/sqrt(px_gamma**2+py_gamma**2)
         sin_fi=py_gamma/sqrt(px_gamma**2+py_gamma**2)

         do i=1, Number_N
            if (En_exc.gt.En_bind(i)) then
               call SpectrumSpace(Number_N, i, En_exc,
     &              cos_teta, sin_teta, cos_fi, sin_fi, 
     &              En_n, alfa_n, beta_n)
               Neut_En(i)=En_n
               teta_n(i)=alfa_n
               fi_n(i)=beta_n
               En_exc=En_exc-En_bind(i)-En_n
            endif
         enddo

      endif
      write(77,*) Number_N, step, en_gamma
      close(77)      
      RETURN
      end

c************************************************************************
*************************************************************************
C*****Simulate generated number of
C*****neutrons
C******************************
      subroutine GRESsimulate(En_gamma, Step, Number_N)
c     Output parameter: Number_N (as default Number_N=0) 
      real En_gamma, Step
      integer Number_N
      
      real Sigma
      real*4 tmp
      real Lambda1, Lambda2, Lambda3         

#include "gpgsigma.inc"
      
      if ((En_gamma.gt.9.097).and.(En_gamma.lt.28.608)) then
CCCCC Calculate Lambda1
         Sigma=-1.      
         do i=1,63
            if ((En_gamma.ge.energy1(i)).and.(En_gamma.lt.energy1(i+1))) 
     &           Sigma=cross1(i)+(cross1(i+1)-cross1(i))*
     &           (En_gamma-energy1(i))/
     &           (energy1(i+1)-energy1(i))
         enddo
         if (Sigma.le.0.) then
            Lambda1=1000.
         else
 1          call GRNDM(tmp,1)
            if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 1
            Lambda1=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22)
         endif
CCCCCCCalculate Lambda2
         Sigma=-1.
         do i=1,61
            if ((En_gamma.ge.energy2(i)).and.(En_gamma.lt.energy2(i+1)))
     &           Sigma=cross2(i)+(cross2(i+1)-cross2(i))*
     &           (En_gamma-energy2(i))/
     &           (energy2(i+1)-energy2(i))
         enddo
         if (Sigma.le.0.) then
            Lambda2=1000.
         else
 2          call GRNDM(tmp,1)
            if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 2      
            Lambda2=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22)      
         endif              
         
CCCCCCCalculate Lambda3
         Sigma=-1.
         do i=1,26
            if ((En_gamma.ge.energy3(i)).and.(En_gamma.lt.energy3(i+1)))
     &           Sigma=cross3(i)+(cross3(i+1)-cross3(i))*
     &           (En_gamma-energy3(i))/
     &           (energy3(i+1)-energy3(i))
         enddo
         if (Sigma.le.0.) then
            Lambda3=1000.
         else
 3          call GRNDM(tmp,1)
            if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 3      
            Lambda3=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22)         
         endif
         
         if (Lambda1.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=1
         if (Lambda2.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=2      
         if (Lambda3.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=3
      endif
      RETURN
      end

c************************************************************************
c************************************************************************
******Simulate neutron energy and flight direction*********************** 

      subroutine SpectrumSpace(Number_N, Num, En_exc,
     &     cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n)
c     Output parameters: En_n, alfa_n, beta_n
      common/bind/En_bind(3)

      real tmp(2)     
      real Pi, Temperature, lev_dens, cos_alfa, sin_alfa, beta
      real cos_alfa_n, cos_beta_n, sin_beta_n
      integer Number_N, Num 
      
      real En_exc
      real cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n

      Pi=4.*atan(1.)	
      if (Number_N.eq.1) then
         En_n=En_exc-En_bind(1)

      else if ((Number_N.eq.2).and.(Num.eq.1)) then
C************************************************************************
C     SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book
C     is: V.C. Barashenkov, V.D. Toneev, "The interactions of high
C     energy particles and atomic nuclei with nuclei", in Russian
C     Calculation of density levels parameter (see (6.31) from reference
C     book)
         lev_dens=0.085*184.
C     Calculation of nucleus temperature (see (6.30) from reference book)
C     (see (6.64) and comments after)       
         Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) 
         
C     Neutron evaporation energy spectrum simulation by exception method
 4       call GRNDM(tmp,2)
C     The average binding energy of neutron
         x_tmp=(En_exc-En_bind(1)-En_bind(2))*tmp(1) 
         y_tmp=(exp(-1.)/Temperature)*tmp(2)
         if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature)))
     &        goto 4
         En_n=x_tmp
      else if ((Number_N.eq.2).and.(Num.eq.2)) then
         En_n=En_exc-En_bind(2)
      else if ((Number_N.eq.3).and.(Num.eq.1)) then
C************************************************************************
C     SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book is:
C     V.C. Barashenkov, V.D. Toneev, "The interactions of high energy
C     particles and atomic nuclei with nuclei", in Russian Calculation
C     of density levels parameter (see (6.31) from reference book)
         lev_dens=0.085*184.
C     Calculation of nucleus temperature (see (6.30) from reference book)
C     (see (6.64) and comments after)       
         Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) 

C     Neutron evaporation energy spectrum simulation by exception method
 5       call GRNDM(tmp,2)
C     Average binding energy of neutron
         x_tmp=(En_exc-En_bind(1)-En_bind(2)-En_bind(3))*tmp(1)
         y_tmp=(exp(-1.)/Temperature)*tmp(2)
         if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature)))
     &        goto 5
         En_n=x_tmp
         
      else if ((Number_N.eq.3).and.(Num.eq.2)) then
C************************************************************************
C     SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book
C     is: V.C. Barashenkov, V.D. Toneev,  "The interactions of high
C     energy particles and atomic nuclei with nuclei", in Russian
C     Calculation of density levels parameter (see (6.31) from reference
C     book) 
         lev_dens=0.085*183.
C     Calculation of nucleus temperature (see (6.30) from reference book)
C     (see (6.64) and comments after)       
         Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) 
         
C     Neutron evaporation energy spectrum simulation by exception method
 6       call GRNDM(tmp,2)
C     Average binding energy of neutron
         x_tmp=(En_exc-En_bind(2)-En_bind(3))*tmp(1)
         y_tmp=(exp(-1.)/Temperature)*tmp(2)
         if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature)))
     &        goto 6
         En_n=x_tmp
         
      else if ((Number_N.eq.3).and.(Num.eq.3)) then
         En_n=En_exc-En_bind(3)
         
      endif
C************************************************************************
C     Neutron direction simulation in gamma coordinate system
      call GRNDM(tmp,2)
      cos_alfa=-1.+2.*tmp(1)
      beta=2.*Pi*tmp(2)
      sin_alfa=sqrt(1.-cos_alfa**2)
C     Return to the GEANT reference coordinate system
      cos_alfa_n=-sin_alfa*cos(beta)*(sin_teta)+cos_alfa*(cos_teta)
      cos_beta_n=sin_alfa*cos(beta)*(cos_teta*cos_fi)-
     &     sin_alfa*sin(beta)*(sin_fi)+cos_alfa*(sin_teta*cos_fi)
      sin_beta_n=sin_alfa*cos(beta)*(cos_teta*sin_fi)+
     &     sin_alfa*sin(beta)*(cos_fi)+cos_alfa*(sin_teta*sin_fi)
      if (abs(cos_beta_n).lt.1.E-6) cos_beta_n=0.
      if (abs(sin_beta_n).lt.1.E-6) sin_beta_n=0.
      if ((cos_beta_n.eq.0.).and.(sin_beta_n.eq.0.)) beta_n=0.
      
      alfa_n=acos(cos_alfa_n)
      if ((cos_beta_n.ge.0.).and.(sin_beta_n.gt.0.)) 
     &     beta_n=atan(sin_beta_n/cos_beta_n)
      if ((cos_beta_n.ge.0.).and.(sin_beta_n.lt.0.))
     &     beta_n=2.*Pi+atan(sin_beta_n/cos_beta_n)     
      if ((cos_beta_n.lt.0.).and.(sin_beta_n.le.0.))
     &     beta_n=Pi+atan(sin_beta_n/cos_beta_n)      
      if ((cos_beta_n.lt.0.).and.(sin_beta_n.ge.0.))
     &     beta_n=Pi+atan(sin_beta_n/cos_beta_n)             
      
      RETURN
      END
