*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
      SUBROUTINE SOURCE ( NOMORE )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2010      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     New source for FLUKA9x-FLUKA20xy:                                *
*                                                                      *
*     Created on 07 January 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on  17-Oct-10    by    Alfredo Ferrari               *
*                                                                      *
*  This is just an example of a possible user written source routine.  *
*  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              *
*                                                                      *
*----------------------------------------------------------------------*
*
*
*
      DIMENSION CUMPR(43), ENEDGE(43)
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(CASLIM)' 
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      LOGICAL LFIRST
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*----------------------------------------------------------------------*
*
* Proton energy group boundaries !43 energies
      DATA ENEDGE / 
     & 0.11521,
     & 0.11655,
     & 0.11787,
     & 0.11918,
     & 0.12045,
     & 0.12171,
     & 0.12296,
     & 0.12413,
     & 0.12539,
     & 0.12666,
     & 0.12788,
     & 0.12909,
     & 0.13028,
     & 0.13139,
     & 0.13261,
     & 0.13383,
     & 0.13501,
     & 0.13618,
     & 0.13733,
     & 0.13845,
     & 0.13955,
     & 0.14073,
     & 0.14189,
     & 0.14300,
     & 0.14409,
     & 0.14524,
     & 0.14635,
     & 0.14745,
     & 0.14855,
     & 0.14962,
     & 0.15066,
     & 0.15177,
     & 0.15287,
     & 0.15386,
     & 0.15497,
     & 0.15710,
     & 0.15920,
     & 0.16127,
     & 0.16332,
     & 0.16534,
     & 0.16738,
     & 0.16938,
     & 0.17137/
* 
* Cumulative spectrum !43 weights
      DATA CUMPR / 
     & 0.01054430535851470, 
     & 0.01431470005844220,
     & 0.02235293146569520,
     & 0.02786211997428630,
     & 0.03479711441176390,
     & 0.04115148724087290,
     & 0.04810046058817790,
     & 0.05417059411222460,
     & 0.06234335477433980,
     & 0.06898244954716260,
     & 0.07635045114455820,
     & 0.08403102422993940,
     & 0.09107531571584640,
     & 0.09881893629872820,
     & 0.10689405545350300,
     & 0.11551085231704100,
     & 0.12307495747963500,
     & 0.13255567582233000,
     & 0.14031266983304600,
     & 0.14992378276338300,
     & 0.15826892365338800,
     & 0.16868329100665900,
     & 0.17827383468625500,
     & 0.18867374702026500,
     & 0.19804611723494700,
     & 0.21038811454095700,
     & 0.21948664581049600,
     & 0.23350115883350400,
     & 0.24206807468343800,
     & 0.25798986891144900,
     & 0.26550517262357100,
     & 0.28579903901021900,
     & 0.29170908517045800,
     & 0.31179289083864600,
     & 0.33021480199225000,
     & 0.36701408972104000,
     & 0.40310877546159700,
     & 0.44540179368478600,
     & 0.49381078238080200,
     & 0.55065625008642900,
     & 0.62417519316010700,
     & 0.71312514880381800,
     & 1/
*
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialization ***
      END IF
*  |
*  +-------------------------------------------------------------------*
*
* In order to cover the target volume in depth the energy must be
* sampled from the ENEDGE, with cumulative weights from CUMPR
* Sample the energy group
      XI = FLRNDM(DUMMY) !FLUKA random number generator, 1 is not included
      DO 500 K = 1, 43 !500 statement label
         IF(XI .LE. CUMPR(K)) THEN !.LE. meaning <=
           ENERGY = ENEDGE(K)
           GO TO 600
         END IF
 500  CONTINUE
      STOP 'Failed to sample the energy group'
 600  CONTINUE
*
*  Push one source particle to the stack. Note that you could as well
*  push many but this way we reserve a maximum amount of space in the
*  stack for the secondaries to be generated
*  Npflka is the stack counter: of course any time source is called it
*  must be =0
      NPFLKA = NPFLKA + 1
*  Wt is the weight of the particle
      WTFLK  (NPFLKA) = ONEONE
      WEIPRI = WEIPRI + WTFLK (NPFLKA)
*  Particle type (1=proton.....). Ijbeam is the type set by the BEAM
*  card
*  +-------------------------------------------------------------------*
*  |  (Radioactive) isotope:
      IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
         IARES  = IPROA
         IZRES  = IPROZ
         IISRES = IPROM
         CALL STISBM ( IARES, IZRES, IISRES )
         IJHION = IPROZ  * 1000 + IPROA
         IJHION = IJHION * 100 + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
*  |
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      ELSE IF ( IJBEAM .EQ. -2 ) THEN
         IJHION = IPROZ  * 1000 + IPROA
         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
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID = IJBEAM
         ILOFLK (NPFLKA) = IJBEAM
*  |  Flag this is prompt radiation
         LRADDC (NPFLKA) = .FALSE.
*  |  Group number for "low" energy neutrons, set to 0 anyway
         IGROUP (NPFLKA) = 0
      END IF
*  |
*  +-------------------------------------------------------------------*
*  From this point .....
*  Particle generation (1 for primaries)
      LOFLK  (NPFLKA) = 1
*  User dependent flag:
      LOUSE  (NPFLKA) = 0
*  No channeling:
      LCHFLK (NPFLKA) = .FALSE.
      DCHFLK (NPFLKA) = ZERZER
*  User dependent spare variables:
      DO 100 ISPR = 1, MKBMX1
         SPAREK (ISPR,NPFLKA) = ZERZER
 100  CONTINUE
*  User dependent spare flags:
      DO 200 ISPR = 1, MKBMX2
         ISPARK (ISPR,NPFLKA) = 0
 200  CONTINUE
*  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
*  Particle age (s)
      AGESTK (NPFLKA) = +ZERZER
      AKNSHR (NPFLKA) = -TWOTWO
*  Kinetic energy of the particle (GeV)
*  With Gaussian sigma spread 0.1% initial energy
      CALL FLNRRN(RGAUSS) ! which returns a normal gaussian number RGAUSS
      TKEFLK (NPFLKA) = ENERGY * (1 + 0.004246609001440 * RGAUSS)
      !WRITE(101,*)RGAUSS
      !TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
*  Particle momentum
      !PMOFLK (NPFLKA) = PBEAM
      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (IONID) ) )
*  Cosines (tx,ty,tz)
      TXFLK  (NPFLKA) = UBEAM
      TYFLK  (NPFLKA) = VBEAM
      TZFLK  (NPFLKA) = WBEAM
*     TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
*    &                       - TYFLK (NPFLKA)**2 )
*  Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
* Particle coordinates
      XFLK   (NPFLKA) = (FLRNDM(DUMMY)-0.5)*2.0*(10)
      YFLK   (NPFLKA) = (FLRNDM(DUMMY)-0.5)*2.0*(10)
      ZFLK   (NPFLKA) = -0.1
      !XFLK   (NPFLKA) = XBEAM
      !YFLK   (NPFLKA) = YBEAM
      !ZFLK   (NPFLKA) = ZBEAM
*  Calculate the total kinetic energy of the primaries: don't change
      IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. ILOFLK (NPFLKA) .GT. 100000 )
     &   THEN
         TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
      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
      RADDLY (NPFLKA) = ZERZER
*  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 of subroutine Source =========================================*
      END

