*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2020:  CERN                                        *
*     All Rights Reserved.                                             *
*                                                                      *
*     Source routine library for FLUKA 4:                              *
*                                                                      *
*     Original source routine                                          *
*     Created on 07 January 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Created on 24 September 2020 by David Horvath & Roberto Versaci  *
*                                             ELI Beamlines            *
*                                                                      *
*     Modified on 23 May 2022      by         David Horvath            *
*                                             ELI Beamlines            *
*                                                                      *
*  Collection of subroutines and functions containing and extending    *
*  the original FLUKA source routine                                   *
*                                                                      *
*----------------------------------------------------------------------*

*======================================================================*

*  Module containing the commons necessary for the source routine

*======================================================================*

      module source_library

         include 'dblprc.inc'
         include 'dimpar.inc'
         include 'iounit.inc'

         include 'beamcm.inc'
         include 'fheavy.inc'
         include 'flkstk.inc'
         include 'ioiocm.inc'
         include 'ltclcm.inc'
         include 'paprop.inc'
         include 'sourcm.inc'
         include 'sumcou.inc'

         type phase_space
            integer :: pc
            double precision :: x, y, z
            double precision :: u, v, w
            double precision :: m_e, wei
         end type phase_space

      end module source_library

*======================================================================*

*  Basic subroutines

*======================================================================*

      subroutine initialization()

         use source_library
         use source_variables

         implicit none

         integer ionid
         logical, save :: lfirst = .true.

         if ( lfirst ) then
            TKESUM = ZERZER
            LUSSRC = .true.

            lfirst = .false.
         end if

         particle_code = IJBEAM

         heavyion_atomic_number = IPROZ
         heavyion_mass_number = IPROA
         heavyion_isomer = IPROM

         if ( particle_code .eq. -2 ) then
            ionid = heavyion_isomer * 100000 + mod ( heavyion_atomic_number, 100 ) * 1000 + heavyion_mass_number
            ionid = ionid * 100 + KXHEAV
            call DCDION ( ionid )
            call SETION ( ionid )

            momentum_energy = PBEAM / AM (ionid) * AMUC12
         else
            momentum_energy = PBEAM
         end if

         energy_logical_flag = .false.

         particle_weight = ONEONE

         divergence_x = DIVBM
         divergence_y = DIVBM
         gaussian_divergence_logical_flag = LDVGSS

         coordinate_x = XBEAM
         coordinate_y = YBEAM
         coordinate_z = ZBEAM

         direction_cosx = UBEAM
         direction_cosy = VBEAM
         direction_cosz = WBEAM
         direction_flag = 0

         polarization_cosx = -TWOTWO
         polarization_cosy = ZERZER
         polarization_cosz = ZERZER

         particle_age = ZERZER
         kshort_component = -TWOTWO
         delayed_radioactive_decay = ZERZER

         return

      end

*======================================================================*

      subroutine set_primary()

         use source_variables

         implicit none

         integer ionid

         call set_internal_flags()
         call set_beam_type( ionid )
         call set_particle_momentum_energy_weight( ionid )
         call set_particle_coordinates()
         call set_particle_direction()
         call set_particle_polarization()
         call set_particle_age()
         call search_starting_region()

      end

*======================================================================*

      subroutine set_internal_flags()

         use source_library

         implicit none

         integer counter

*     NPFLKA is the stack counter: of course any time source
*     is called it must be =0
         NPFLKA = NPFLKA + 1
*     From this point ...
*     Particle generation (1 for primaries)
         LOFLK  (NPFLKA) = 1
*     User dependent flag:
         LOUSE  (NPFLKA) = 0
*     No channeling:
         KCHFLK (NPFLKA) = 0
         ECRFLK (NPFLKA) = ZERZER
*     Extra infos:
         INFSTK (NPFLKA) = 0
         LNFSTK (NPFLKA) = 0
         ANFSTK (NPFLKA) = ZERZER
*     Parent variables:
         IPRSTK (NPFLKA) = 0
         EKPSTK (NPFLKA) = ZERZER
*     User dependent spare variables:
         do counter = 1, MKBMX1
            SPAREK (counter,NPFLKA) = ZERZER
         enddo
*     User dependent spare flags:
         do counter = 1, MKBMX2
            ISPARK (counter,NPFLKA) = 0
         enddo
*     Save the track number of the stack particle:
         ISPARK (MKBMX2,NPFLKA) = NPFLKA
         NPARMA = NPARMA + 1
         NUMPAR (NPFLKA) = NPARMA
         NEVENT (NPFLKA) = 0
         DFNEAR (NPFLKA) = +ZERZER
*     ... to this point: don't change anything

         return

      end

*======================================================================*

      subroutine set_beam_type( ionid )

         use source_library
         use source_variables

         implicit none

         logical lisnut
         logical lightion_flag
         logical tritium_flag

         integer ij, ionid

*     Statement function:
         lisnut (ij) = INDEX ( PRNAME (ij), 'NEUTRI' ) .gt. 0

         lightion_flag = heavyion_atomic_number .le. 3 .and. heavyion_mass_number .le. 4
         tritium_flag = heavyion_atomic_number .eq. 1 .and. heavyion_mass_number .eq. 3

         if ( particle_code .eq. -2 ) then

*        (Radioactive) isotope:
            if ( LRDBEA ) then
               if ( lightion_flag .and. ( .not. tritium_flag ) ) then
                  call FLABRT ('set_beam_type','Use explicit light ion particle')
               end if
               call STISBM ( heavyion_mass_number, heavyion_atomic_number, heavyion_isomer )
               IJHION = heavyion_isomer * 100000 + mod ( heavyion_atomic_number, 100 ) * 1000 + heavyion_mass_number
               IJHION = IJHION * 100 + KXHEAV
               ionid  = IJHION
               call DCDION ( ionid )
               call SETION ( ionid )
               LFRPHN (NPFLKA) = .false.

*        Heavy ion:
            else
               if ( lightion_flag ) then
                  call FLABRT ('set_beam_type','Use explicit light ion particle')
               end if
               IJHION = heavyion_isomer * 100000 + mod ( heavyion_atomic_number, 100 ) * 1000 + heavyion_mass_number
               IJHION = IJHION * 100 + KXHEAV
               ionid  = IJHION
               CALL DCDION ( ionid )
               CALL SETION ( ionid )
               ILOFLK (NPFLKA) = IJHION
*           Flag this is prompt radiation
               LRADDC (NPFLKA) = .false.
*           Group number for "low" energy neutrons, set to 0 anyway
               IGROUP (NPFLKA) = 0
*           Parent radioactive isotope:
               IRDAZM (NPFLKA) = 0
               LFRPHN (NPFLKA) = .false.
            end if

*     Normal hadron:
         else
            ionid = particle_code
            ILOFLK (NPFLKA) = particle_code
*        Flag this is prompt radiation
            LRADDC (NPFLKA) = .false.
*        Group number for "low" energy neutrons, set to 0 anyway
            IGROUP (NPFLKA) = 0
*        Parent radioactive isotope:
            IRDAZM (NPFLKA) = 0
*        Check if it is a neutrino, if so force the interaction
*        (unless the relevant flag has been disabled)
            if ( lisnut (particle_code) .and. LNUFIN ) then
               LFRPHN (NPFLKA) = .true.
*        Not a neutrino
            else
               LFRPHN (NPFLKA) = .false.
            end if
         end if

         return

      end

*======================================================================*

      subroutine set_particle_age()

         use source_library
         use source_variables

         implicit none

         AGESTK (NPFLKA) = particle_age
         AKNSHR (NPFLKA) = kshort_component
         RADDLY (NPFLKA) = delayed_radioactive_decay

         return

      end

*======================================================================*

      subroutine set_particle_polarization()

         use source_library
         use source_variables

         implicit none

         TXPOL  (NPFLKA) = polarization_cosx
         TYPOL  (NPFLKA) = polarization_cosy
         TZPOL  (NPFLKA) = polarization_cosz

         return

      end

*======================================================================*

      subroutine set_particle_coordinates()

         use source_library
         use source_variables

         implicit none

         XFLK (NPFLKA) = coordinate_x
         YFLK (NPFLKA) = coordinate_y
         ZFLK (NPFLKA) = coordinate_z

         return

      end

*======================================================================*

      subroutine set_particle_direction()

         use source_library
         use source_variables

         implicit none

         double precision squared
         double precision normalization

         double precision, dimension (3,3) :: dummyaxes

         if ( direction_flag .eq. 0 ) then
            normalization = sqrt ( direction_cosx**2 + direction_cosy**2 + direction_cosz**2 )
            direction_cosx = direction_cosx / normalization
            direction_cosy = direction_cosy / normalization
            direction_cosz = direction_cosz / normalization
         else
            squared = ONEONE - direction_cosx**2 - direction_cosy**2

            if ( squared .lt. ZERZER ) then
               call FLABRT('set_particle_direction','Invalid direction values')
            end if

            if ( direction_flag .eq. 1 ) then
               direction_cosz = sqrt ( squared )
            else if ( direction_flag .eq. 2 ) then
               direction_cosz = -ONEONE * sqrt ( squared )
            else
               call FLABRT('set_particle_direction','Invalid direction flag')
            end if
         end if

         if ( ( divergence_x .lt. ANGLGB ) .and. ( divergence_y .lt. ANGLGB ) ) then
            divergence_x = ZERZER
            divergence_y = ZERZER
         else if ( divergence_x .lt. ANGLGB ) then
            divergence_x = ANGLGB
         else if ( divergence_y .lt. ANGLGB ) then
            divergence_y = ANGLGB
         end if

         call STBMDV ( .false., direction_cosx, direction_cosy, direction_cosz, divergence_x, divergence_y, dummyaxes,
     &      gaussian_divergence_logical_flag, .false. )

         return
      end

*======================================================================*

      subroutine set_particle_momentum_energy_weight( ionid )

         use source_library
         use source_variables

         implicit none

         integer ionid

*     Not (Radioactive) isotope:
         if ( .not. ( particle_code .eq. -2 .and. LRDBEA ) ) then
            if ( energy_logical_flag ) then
               if ( particle_code .eq. -2 ) then
                  TKEFLK (NPFLKA) = momentum_energy * AM (ionid) / AMUC12
               else
                  TKEFLK (NPFLKA) = momentum_energy
               end if
               PMOFLK (NPFLKA) = sqrt ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA) + TWOTWO * AM (ionid) ) )
            else
               if ( particle_code .eq. -2 ) then
                  PMOFLK (NPFLKA) = momentum_energy * AM (ionid) / AMUC12
               else
                  PMOFLK (NPFLKA) = momentum_energy
               end if
               TKEFLK (NPFLKA) = sqrt ( PMOFLK (NPFLKA)**2 + AM (ionid)**2 ) - AM (ionid)
            end if
         end if

         WTFLK  (NPFLKA) = particle_weight
         WEIPRI = WEIPRI + WTFLK (NPFLKA)

*     Calculate the total kinetic energy of the primaries:
*     (Radioactive) isotope:
         if ( particle_code .eq. -2 .and. LRDBEA ) then
*     Heavy ion:
         else if ( ILOFLK (NPFLKA) .eq. -2 .or. ILOFLK (NPFLKA) .gt. 100000 ) then
            TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
*     Standard particle:
         else if ( ILOFLK (NPFLKA) .ne. 0 ) then
            TKESUM = TKESUM + ( TKEFLK (NPFLKA) + AMDISC (ILOFLK(NPFLKA)) ) * WTFLK (NPFLKA)
         else
            TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
         end if

         return
      end

*======================================================================*

      subroutine search_starting_region()

         use source_library

         implicit none

         integer idisc

*     Here we ask for the region number of the hitting point.
*        NREG (NPFLKA) = ...
*     The following line makes the starting region search much more
*     robust if particles are starting very close to a boundary:
         call GEOCRS ( TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA) )
         call GEOREG ( XFLK  (NPFLKA), YFLK  (NPFLKA), ZFLK  (NPFLKA), NRGFLK(NPFLKA), idisc )
*     Do not change these cards:
         call GEOHSM ( NHSPNT (NPFLKA), 1, -11, MLATTC )
         NLATTC (NPFLKA) = MLATTC
         CMPATH (NPFLKA) = ZERZER
         call SOEVSV

         return
      end

*======================================================================*

*  Extra sampling functions and subroutines

*======================================================================*

      recursive function sample_flat_momentum_energy( min, max ) result( sample )

         use source_library

         implicit none

         double precision min, max, sample
         double precision FLRNDM, xdummy

         sample = min + (max - min) * FLRNDM( xdummy )

         if ( sample .le. ANGLGB ) then
            sample = sample_flat_momentum_energy ( min, max )
         end if

         return
      end

*======================================================================*

      recursive function sample_gaussian_momentum_energy( mean, fwhm ) result( sample )

         use source_library

         implicit none

         double precision mean, fwhm, sample, random

         call FLNRRN( random )

         do while ( random .gt. 12.D+00 )
            call FLNRRN ( random )
         end do

         sample = mean + fwhm * random / S2FWHM

         if ( sample .le. ANGLGB ) then
            sample = sample_gaussian_momentum_energy ( mean, fwhm )
         end if

         return
      end

*======================================================================*

      double precision function sample_maxwell_boltzmann_energy( temperature )

         use source_library

         implicit none

         double precision temperature
         double precision random1, random2, random3
         double precision FLRNDM, xdummy

         random1 = FLRNDM( xdummy )
         random2 = FLRNDM( xdummy )
         random3 = FLRNDM( xdummy )

         sample_maxwell_boltzmann_energy = -temperature * ( log( random1 ) + log( random2 ) * ( cos( PIPIPI * random3 / 2 ) )**2 )

         return
      end

*======================================================================*

      recursive function sample_histogram_momentum_energy( filename, energy_unit ) result ( sample )

         use source_library

         implicit none

         logical, save :: first_run = .true.
         character(*) :: filename, energy_unit
         integer ierr, io, line, bins, random_bin
         double precision sample
         double precision t1, t2, t3, FLRNDM, xdummy, energy_unit_to_scale_factor, scale_factor
         double precision, allocatable, dimension(:) :: e_min, e_max, intensity, probability
         double precision random_probability

         save e_min, e_max, intensity, probability, bins

         if ( first_run ) then
            write ( LUNOUT, * ) '# ------------------------------------------------'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Sampling from histogram requested'
            write ( LUNOUT, * ) '#'

            scale_factor = energy_unit_to_scale_factor(energy_unit)

            bins = 0

            call OAUXFI( filename, LUNRDB, 'OLD', ierr )

            if ( ierr .eq. 1 ) then
               call FLABRT('sample_histogram_momentum_energy', 'Unable to open histogram file')
            end if

            do
               read ( LUNRDB, *, iostat = io ) t1, t2, t3
               if ( io .eq. 0 ) then
                  bins = bins + 1
*           Failed to read line corretly, possible comment - skipping
               else if ( io .gt. 0 ) then
                  cycle
*           End of file
               else if ( io .lt. 0 ) then
                  exit
               end if
            end do

            allocate( e_min( bins ) )
            allocate( e_max( bins ) )
            allocate( intensity( bins ) )
            allocate( probability( bins ) )

            write ( LUNOUT, "(A, I6)" ) ' # Number of bins: ', bins
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Processed histogram:'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Emin [GeV],   Emax [GeV],    dN/dE'

            line = 0

            rewind( LUNRDB )

            do
               read ( LUNRDB, *, iostat = io ) t1, t2, t3
               if ( io .eq. 0 ) then
                  line = line + 1
                  e_min(line) = t1 * scale_factor
                  e_max(line) = t2 * scale_factor
                  intensity(line) = t3
                  write ( LUNOUT, "(3ES15.7)" ) e_min(line), e_max(line), intensity(line)
*           Failed to read line corretly, possible comment - skipping
               else if ( io .gt. 0 ) then
                  cycle
*           End of file
               else if ( io .lt. 0) then
                  exit
               end if
            end do

            do line = 1, bins
               intensity(line) = intensity(line) * (e_max(line) - e_min(line))
            end do

            probability(1) = intensity(1)

            do line = 2, bins
               probability(line) = probability(line - 1) + intensity(line)
            end do

            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# End of histogram'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# ------------------------------------------------'

            probability = probability / probability(bins)

            first_run = .false.
         end if

         random_probability = FLRNDM( xdummy )

         do line = 1, bins
            if (probability(line) .ge. random_probability) then
               random_bin = line
               exit
            end if
         end do

         sample = e_min(random_bin) + ( e_max(random_bin) - e_min(random_bin) ) * FLRNDM( xdummy )

         if ( sample .lt. ANGLGB ) then
            sample = sample_histogram_momentum_energy ( filename, energy_unit )
         end if

         return
      end

*======================================================================*

      recursive function sample_spectrum_momentum_energy( filename, energy_unit ) result ( sample )

         use source_library

         implicit none

         logical, save :: first_run = .true.
         character(*) :: filename, energy_unit
         integer ierr, io, line, bins, random_bin
         double precision sample
         double precision t1, t2, FLRNDM, xdummy, energy_unit_to_scale_factor, scale_factor
         double precision, allocatable, dimension(:) :: energy, intensity, probability
         double precision random_probability

         save energy, intensity, probability, bins

         if ( first_run ) then
            write ( LUNOUT, * ) '# ------------------------------------------------'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Sampling from spectrum requested'
            write ( LUNOUT, * ) '#'

            scale_factor = energy_unit_to_scale_factor(energy_unit)

            bins = -1

            call OAUXFI( filename, LUNRDB, 'OLD', ierr )

            if ( ierr .eq. 1 ) then
               call FLABRT('sample_spectrum_momentum_energy', 'Unable to open spectrum file')
            end if

            do
               read ( LUNRDB, *, iostat = io ) t1, t2
               if ( io .eq. 0 ) then
                  bins = bins + 1
*           Failed to read line corretly, possible comment - skipping
               else if ( io .gt. 0 ) then
                  cycle
*           End of file
               else if ( io .lt. 0 ) then
                  exit
               end if
            end do

            allocate( energy( 0:bins ) )
            allocate( intensity( 0:bins ) )
            allocate( probability( 0:bins ) )

            write ( LUNOUT, "(A, I6)" ) ' # Number of points: ', bins
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Processed spectrum:'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# E [GeV],       Intensity'

            line = 0

            rewind( LUNRDB )

            do
               read ( LUNRDB, *, iostat = io ) t1, t2
               if ( io .eq. 0 ) then

                  energy(line) = t1 * scale_factor
                  intensity(line) = t2

                  write ( LUNOUT, "(2ES15.7)" ) energy(line), intensity(line)

                  line = line + 1

*           Failed to read line corretly, possible comment - skipping
               else if ( io .gt. 0 ) then
                  cycle
*           End of file
               else if ( io .lt. 0) then
                  exit
               end if
            end do

            probability(0) = ZERZER

            do line = 1, bins
               probability(line) = (energy(line) - energy(line-1)) * (intensity(line) + intensity(line - 1)) + probability(line - 1)
            end do

            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# End of spectrum'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# ------------------------------------------------'

            probability = probability / probability(bins)

            first_run = .false.

         end if

         random_probability = FLRNDM( xdummy )

         do line = 1, bins
            if (probability(line) .ge. random_probability) then
               random_bin = line
               exit
            end if
         end do

         sample = (random_probability - probability(line - 1)) / (probability(line) - probability(line - 1))

         if (intensity(line) == intensity(line - 1)) then
            continue
         else
            sample = sample * (intensity(line) - intensity(line - 1)) * (intensity(line) + intensity(line - 1))
            sample = (sqrt(sample + intensity(line - 1)**2) - intensity(line - 1)) / (intensity(line) - intensity(line - 1))
         end if

         sample = sample * (energy(line) - energy(line - 1)) + energy(line - 1)

         if ( sample .lt. ANGLGB ) then
            sample = sample_spectrum_momentum_energy ( filename, energy_unit )
         end if

         return
      end

*======================================================================*

      subroutine read_phase_space_file( filename, energy_unit, length_unit, return_entry, sequential, nomore )

         use source_library
         use source_variables

         implicit none

         logical first_run, sequential
         character(*) :: filename, energy_unit, length_unit
         integer ierr, io, line, entries, nomore
         double precision FLRNDM, xdummy
         double precision energy_unit_to_scale_factor, energy_scale_factor
         double precision length_unit_to_scale_factor, length_scale_factor

         type(phase_space) :: return_entry, ent
         type(phase_space), allocatable, dimension(:) :: phase_space_data

         data first_run / .true. /

         save first_run, entries, line

         if ( first_run ) then
            write ( LUNOUT, * ) '# ------------------------------------------------'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Reading of phase space file requested'
            write ( LUNOUT, * ) '#'

            energy_scale_factor = energy_unit_to_scale_factor( energy_unit )
            length_scale_factor = length_unit_to_scale_factor( length_unit )

            entries = 0

            call OAUXFI( filename, LUNRDB, 'OLD', ierr )

            if ( ierr .eq. 1 ) then
               call FLABRT('read_phase_space_file','Unable to open phase space file')
            end if

            do
               read ( LUNRDB, *, iostat = io ) ent%pc, ent%m_e, ent%x, ent%y, ent%z, ent%u, ent%v, ent%w, ent%wei
               if ( io .eq. 0 ) then
                  entries = entries + 1
*           Failed to read line corretly, possible comment - skipping
               else if ( io .gt. 0 ) then
                  cycle
*           End of file
               else if ( io .lt. 0 ) then
                  exit
               end if
            end do

            allocate ( phase_space_data( entries ) )

            write ( LUNOUT, "(A, I10)" ) ' # Number of entries in phase space file: ', entries
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Processed phase space data:'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# Particle    E [GeV]        X [cm]         Y [cm]         Z [cm]         ',
     &         'Cosx           Cosy           Cosz           Weight'

            line = 0

            rewind( LUNRDB )

            do
               read ( LUNRDB, *, iostat = io ) ent%pc, ent%m_e, ent%x, ent%y, ent%z, ent%u, ent%v, ent%w, ent%wei
               if ( io .eq. 0 ) then
                  line = line + 1
                  phase_space_data(line)%pc  = ent%pc
                  phase_space_data(line)%m_e = ent%m_e * energy_scale_factor
                  phase_space_data(line)%x   = ent%x * length_scale_factor
                  phase_space_data(line)%y   = ent%y * length_scale_factor
                  phase_space_data(line)%z   = ent%z * length_scale_factor
                  phase_space_data(line)%u   = ent%u
                  phase_space_data(line)%v   = ent%v
                  phase_space_data(line)%w   = ent%w
                  phase_space_data(line)%wei = ent%wei

                  if ( line .le. 100 ) then
                     write ( LUNOUT, "(I11, 8ES15.7)" )
     &                  phase_space_data(line)%pc, phase_space_data(line)%m_e,
     &                  phase_space_data(line)%x, phase_space_data(line)%y, phase_space_data(line)%z,
     &                  phase_space_data(line)%u, phase_space_data(line)%v, phase_space_data(line)%w,
     &                  phase_space_data(line)%wei
                  end if
*           Failed to read line corretly, possible comment - skipping
               else if ( io .gt. 0 ) then
                  cycle
*           End of file
               else if ( io .lt. 0 ) then
                  exit
               end if
            end do

            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# End of phase space data (Only the first 100 entries shown)'
            write ( LUNOUT, * ) '#'
            write ( LUNOUT, * ) '# ------------------------------------------------'

            if ( sequential ) line = 0

            first_run = .false.
         end if

         if ( sequential ) then
            line = line + 1

            if ( line .gt. entries ) then
               nomore = 1

               write ( LUNOUT, * ) '# ------------------------------------------------'
               write ( LUNOUT, * ) '#'
               write ( LUNOUT, * ) '# All entries of the phase space data has been replayed'
               write ( LUNOUT, * ) '#'
               write ( LUNOUT, * ) '# ------------------------------------------------'

               return_entry = phase_space_data( line - 1 )
            else
               return_entry = phase_space_data( line )
            end if
         else
            return_entry = phase_space_data( floor( entries * FLRNDM( xdummy ) + 1 ) )
         end if

         return
      end

*======================================================================*

      double precision function energy_unit_to_scale_factor( energy_unit )

         use source_library

         implicit none

         character(*) :: energy_unit

         select case (energy_unit)
         case ('TeV')
            energy_unit_to_scale_factor = 1.0D+3
         case ('TeV/c')
            energy_unit_to_scale_factor = 1.0D+3
         case ('GeV')
            energy_unit_to_scale_factor = ONEONE
         case ('GeV/c')
            energy_unit_to_scale_factor = ONEONE
         case ('MeV')
            energy_unit_to_scale_factor = 1.0D-3
         case ('MeV/c')
            energy_unit_to_scale_factor = 1.0D-3
         case ('keV')
            energy_unit_to_scale_factor = 1.0D-6
         case ('keV/c')
            energy_unit_to_scale_factor = 1.0D-6
         case ('eV')
            energy_unit_to_scale_factor = 1.0D-9
         case ('eV/c')
            energy_unit_to_scale_factor = 1.0D-9
         case ('J')
            energy_unit_to_scale_factor = 6.241506D+9
         case default
            energy_unit_to_scale_factor = ZERZER
            call FLABRT ('energy_unit_to_scale_factor', 'Unknown energy unit')
         end select

         return
      end

*======================================================================*

      double precision function length_unit_to_scale_factor( length_unit )

         use source_library

         implicit none

         character(*) :: length_unit

         select case (length_unit)
         case ('km')
            length_unit_to_scale_factor = 1.0D+5
         case ('m')
            length_unit_to_scale_factor = 1.0D+2
         case ('cm')
            length_unit_to_scale_factor = ONEONE
         case ('mm')
            length_unit_to_scale_factor = 1.0D-1
         case default
            length_unit_to_scale_factor = ZERZER
            call FLABRT ('length_unit_to_scale_factor', 'Unknown length unit')
         end select

         return
      end

*======================================================================*

      subroutine sample_exponential_energy_weight( e_min, e_max, intensity_ratio, kinetic_energy, weight )

         use source_library

         implicit none

         double precision FLRNDM, xdummy

*     Input variables
         double precision e_min, e_max, intensity_ratio
*     Output variables
         double precision kinetic_energy, weight

         kinetic_energy = e_min + ( e_max - e_min ) * FLRNDM( xdummy )
         weight = exp( kinetic_energy * log( intensity_ratio ) / ( e_max - e_min ) )

         return
      end

*======================================================================*

      function sample_flat_distribution( min, max ) result( sample )

         use source_library

         implicit none

         double precision min, max, sample
         double precision FLRNDM, xdummy

         sample = min + (max - min) * FLRNDM( xdummy )

         return
      end

*======================================================================*

      function sample_gaussian_distribution( mean, fwhm ) result( sample )

         use source_library

         implicit none

         double precision mean, fwhm, sample, random

         call FLNRRN( random )

         sample = mean + fwhm * random / S2FWHM

         return
      end

*======================================================================*

      subroutine sample_annular_distribution( rmin, rmax, axis_1, axis_2 )

         use source_library

         implicit none

         double precision rmin, rmax, axis_1, axis_2
         double precision sin_theta, cos_theta
         double precision radius, rmin_squared, rmax_squared
         double precision constant
         double precision FLRNDM, xdummy

         call SFECFE( sin_theta, cos_theta )

         rmin_squared = rmin * rmin
         rmax_squared = rmax * rmax

         constant = TWOTWO / ( rmax_squared - rmin_squared )
         radius = sqrt( TWOTWO * FLRNDM( xdummy ) / constant + rmin_squared )

         axis_1 = axis_1 + radius * cos_theta
         axis_2 = axis_2 + radius * sin_theta

         return
      end

*======================================================================*

      subroutine sample_isotropic_direction( dirx, diry, dirz )

         use source_library

         implicit none

         double precision dirx, diry, dirz

         call RACO( dirx, diry, dirz )

         return
      end

*======================================================================*

      subroutine print_primary( max_lines )

         use source_library
         use source_variables

         integer, save :: counter = 0

         logical, save :: lfirst = .true.

         double precision m_e

         if ( lfirst ) then
            write ( *, * ) 'source_newgen.f - Debug output'
            write ( *, * ) '-------------------------------'

            if (energy_logical_flag) then
               write ( *, * ) 'Particle -E [GeV]        X [cm]         Y [cm]         Z [cm]         ',
     &            'Cosx           Cosy           Cosz           Weight'
            else
               write ( *, * ) 'Particle  Mom. [GeV/c]   X [cm]         Y [cm]         Z [cm]         ',
     &            'Cosx           Cosy           Cosz           Weight'
            endif

            lfirst = .false.
         end if

         if ( counter .ge. max_lines ) return

         m_e = momentum_energy
         if ( energy_logical_flag ) then
            m_e = -ONEONE * m_e
         end if

         write ( *, "(I9, 8ES15.7)" ) particle_code, m_e, coordinate_x, coordinate_y, coordinate_z,
     &      TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA), particle_weight

         counter = counter + 1

         return
      end

*======================================================================*
