*
*----------------------------------------------------------------------*
*                                                                      *
*     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
      ! ====================================

      double precision random_variable
      double precision photon_electron_proportion


      ! ==============================
      ! 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
      ! ==============================


      ! 1. Accessing variables from the SOURCE card
      ! ===========================================
      ! Values set on the SOURCE card can be accessed with the following variables:
      !    Numerical values (double precision):
      !       WHASOU(1), WHATSOU(2), .., WHASOU(18)
      !    SDUM text (8 character):
      !       SDUSOU


      ! 2. Primary particle
      ! ===================

      ! 2.1. Particle code
      ! ------------------
      ! FLUKA particle code of the primary
      !    See section 5.1 of the FLUKA manual for the list of particle code
      ! Default:
      !    Particle code of the primary defined on the BEAM card if present, otherwise 1 (proton)

*      particle_code = ...

      random_variable = FLRNDM(XDUMMY)
      photon_electron_proportion = WHASOU(1)
      if ( random_variable .lt. photon_electron_proportion ) then
         ! 3 := Electron
         particle_code = 3
      else 
         ! 7 := Photon
         particle_code = 7
      end if


      ! 3. Particle momentum / energy and weight
      ! ========================================

      ! 3.4. Sampling functions and subroutines
      ! ---------------------------------------

      ! 3.4.2. Gaussian distribution
      ! ------

*      momentum_energy = sample_gaussian_momentum_energy( [mean], [fwhm] )

      energy_logical_flag = .true.
      if ( random_variable .lt. photon_electron_proportion ) then
         ! Electron, Maximum = 40 MeV, FWHM = 5 MeV
         momentum_energy = sample_gaussian_momentum_energy( 40.0D-3, 5.0D-3 )
      else 
         ! Photon, Maximum = 10 MeV, FWHM = 2 MeV
         momentum_energy = sample_gaussian_momentum_energy( 10.0D-3, 5.0D-3 )
      end if


      ! 9. Debugging
      ! ============

      ! 9.1. Debug logical flag
      ! -----------------------
      ! Enables or disables the printout of the beam parameters for debugging
      ! Possible values:
      !    .false. : Debug output disabled (Default)
      !    .true.  : Debug output enabled

      debug_logical_flag = .true.
      
      ! 9.2. Debug lines
      ! ----------------
      ! Sets the maximum number of lines printed in the debug output
      ! Default:
      !    100

*      debug_lines = ...

      debug_lines = 50


      ! ==============================================
      ! 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
