*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2020:  CERN                                        *
*     All Rights Reserved.                                             *
*                                                                      *
*     Source routine or FLUKA 4:                                       *
*                                                                      *
*     Created on 24 September 2020 by David Horvath & Roberto Versaci  *
*                                             ELI Beamlines            *
*                                                                      *
*     Modified on 23 May 2022      by         David Horvath            *
*                                             ELI Beamlines            *
*                                                                      *
*  This is a simplified user written source routine utilizing a        *
*  separate source routine library.                                    *
*                                                                      *
*  It is intended as an alternative new-user-friendly version of the   *
*  source.f routine. Existing FLUKA 4 source routines remain           *
*  compatible.                                                         *
*                                                                      *
*  Note that the beam card still has some meaning - in the scoring the *
*  maximum momentum used in deciding the binning is taken from the     *
*  beam momentum.  Other beam card parameters are obsolete.            *
*                                                                      *
*       Output variables:                                              *
*                                                                      *
*              nomore = if > 0 the run will be terminated              *
*                                                                      *
*----------------------------------------------------------------------*
*
*  Quick start guide:
*  ------------------
*
*  This user source routine template aims to modernize the legacy routine
*  by implementing modern Fortran conventions and to provide built in
*  sampling functions.
*
*  The users only need to change / add code between the BEGINNING and END
*  marks, one section for declaration of user variables, and one for
*  assigning values to the beam parameters.
*
*  By default there is no user variable defined, and all code lines for
*  parameter assignment are commented out. These comments start with the
*  symbol: '*'. To enable one, the '*' needs to be deleted.
*
*  (Note: In Fortran each code line should start in column 7 or further in.)
*
*  Every beam parameter has a default value based on the FLUKA input file.
*  A parameter assignment should only be used if the default value has to be
*  modified.
*
*  There are three ways to assign a value to a parameter:
*
*     1. Direct assignment: A parameter is equal to a value. For example:
*
*           momentum_energy = 0.1D0
*
*        If the parameter defined as a double precision, then the assigned
*        value should be represented as double precision as well so as to
*        not loose numerical precision. To do this a 'D' exponential mark
*        must be used.
*
*     2. Using a sampling function: A parameter is assigned to a value,
*        which is calculated by a separate function. For example:
*
*           coordinate_x = sample_flat_distribution( [min], [max] )
*
*        The parameters with the '[' and ']' brackets need to be replaced
*        with numbers, or user variables containing the desired values, like:
*
*           coordinate_x = sample_flat_distribution( -1.0D0, 1.0D0 )
*
*     3. Using a sampling subroutine: They are similar to function, but they
*        are not returning values directly; instead they modify the variables
*        in their argument list. For example:
*
*           call sample_annular_distribution( [rmin], [rmax], coordinate_x, coordinate_y )
*
*        The example above has two input parameters with brackets, and two
*        output parameters (without bracket). The input parameters have to be
*        provided, as for functions. The output parameter names usually don't
*        need to be changed, but there are cases, where a subset of possible
*        output parameters has to be selected.
*
*  For further details see the FLUKA manual.
*
*----------------------------------------------------------------------*

      module source_variables

         implicit none

         integer, save :: particle_code
         integer, save :: heavyion_atomic_number, heavyion_mass_number, heavyion_isomer

         double precision, save :: momentum_energy, particle_weight
         logical, save :: energy_logical_flag

         double precision, save :: divergence_x, divergence_y
         logical, save :: gaussian_divergence_logical_flag

         double precision, save :: coordinate_x, coordinate_y, coordinate_z

         integer, save :: direction_flag
         double precision, save :: direction_cosx, direction_cosy, direction_cosz

         double precision, save :: polarization_cosx, polarization_cosy, polarization_cosz

         double precision, save :: particle_age
         double precision, save :: kshort_component
         double precision, save :: delayed_radioactive_decay

      end module source_variables

      include 'source_library.inc'

*=== Source ===========================================================*

      subroutine SOURCE ( nomore )

      use source_library
      use source_variables

      implicit none

      double precision xdummy
      integer :: nomore, debug_lines = 100
      logical, save :: first_run = .true., debug_logical_flag = .false.

      type(phase_space) phase_space_entry

      ! Function declarations

      double precision sample_flat_momentum_energy
      double precision sample_gaussian_momentum_energy
      double precision sample_maxwell_boltzmann_energy
      double precision sample_histogram_momentum_energy
      double precision sample_spectrum_momentum_energy

      double precision sample_gaussian_distribution
      double precision sample_flat_distribution

      double precision FLRNDM

      ! ====================================
      ! BEGINNING of user declared variables
      ! ====================================



      ! ==============================
      ! END of user declared variables
      ! ==============================

      nomore = 0

      call initialization()

      if ( first_run ) then

         ! ==================================
         ! BEGINNING of custom initialization
         ! ==================================



         ! ============================
         ! END of custom initialization
         ! ============================

         first_run = .false.
      end if

      ! ==============================
      ! BEGINNING of customizable code
      ! ==============================

      call read_phase_space_file( "phase_space_file.txt", "GeV/c", "cm", phase_space_entry, .true., nomore )

      particle_code   = phase_space_entry%pc
      momentum_energy = phase_space_entry%m_e

      coordinate_x = phase_space_entry%x
      coordinate_y = phase_space_entry%y
      coordinate_z = phase_space_entry%z

      direction_cosx = phase_space_entry%u
      direction_cosy = phase_space_entry%v
      direction_cosz = phase_space_entry%w

      particle_weight = phase_space_entry%wei

      ! ==============================================
      ! END of customizable code - Do not change below
      ! ==============================================

      if ( nomore .eq. 0 ) then

         call set_primary()

         if ( debug_logical_flag ) call print_primary( debug_lines )

      end if

      return
*=== End of subroutine Source =========================================*
      end
