*
*=== source ===========================================================*
*
      SUBROUTINE SOURCE ( NOMORE )

      INCLUDE 'dblprc.inc'
      INCLUDE 'dimpar.inc'
      INCLUDE 'iounit.inc'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 2003-2019:  CERN & INFN                            *
*     All Rights Reserved.                                             *
*                                                                      *
*     New source for FLUKA9x-FLUKA20xy:                                *
*                                                                      *
*     Created on 07 January 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*  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              *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE 'beamcm.inc'
      INCLUDE 'fheavy.inc'
      INCLUDE 'flkstk.inc'
      INCLUDE 'ioiocm.inc'
      INCLUDE 'ltclcm.inc'
      INCLUDE 'paprop.inc'
      INCLUDE 'sourcm.inc'
      INCLUDE 'sumcou.inc'
*
      LOGICAL LFIRST, LISNUT
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*  Statement function:
      LISNUT (IJ) = INDEX ( PRNAME (IJ), 'NEUTRI' ) .GT. 0
*======================================================================*
*                                                                      *
*                 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
*  |
*  +-------------------------------------------------------------------*
*  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 = IPROM  * 100000 + MOD ( IPROZ, 100 ) * 1000 + IPROA
         IJHION = IJHION * 100    + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         LFRPHN (NPFLKA) = .FALSE.
*  |
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      ELSE IF ( IJBEAM .EQ. -2 ) THEN
         IJHION = IPROM  * 100000 + MOD ( IPROZ, 100 ) * 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
*  |  Parent radioactive isotope:
         IRDAZM (NPFLKA) = 0
*  |  Particle age (s)
         AGESTK (NPFLKA) = +ZERZER
*  |  Kinetic energy of the particle (GeV)
         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) ) )
         LFRPHN (NPFLKA) = .FALSE.
*  |
*  +-------------------------------------------------------------------*
*  |  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
*  |  Parent radioactive isotope:
         IRDAZM (NPFLKA) = 0
*  |  Particle age (s)
         AGESTK (NPFLKA) = +ZERZER
*  |  Kinetic energy of the particle (GeV)
         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) ) )
*  |  +----------------------------------------------------------------*
*  |  |  Check if it is a neutrino, if so force the interaction
*  |  |  (unless the relevant flag has been disabled)
         IF ( LISNUT (IJBEAM) .AND. LNUFIN ) THEN
            LFRPHN (NPFLKA) = .TRUE.
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |  Not a neutrino
         ELSE
            LFRPHN (NPFLKA) = .FALSE.
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
      END IF
*  |
*  +-------------------------------------------------------------------*
*  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 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
      AKNSHR (NPFLKA) = -TWOTWO
*  Cosines (tx,ty,tz)
      CALL RACO(TXFLK(NPFLKA),TYFLK(NPFLKA),TZFLK(NPFLKA))
*     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
*  Marine irradiator - 25 active pencils
*       open (unit=10, file = 'test.txt')
	RNDNUM = FLRNDM(DUMMY)
	IF (RNDNUM .LE. 0.025) THEN 
		XORIG = 43
		YORIG = 0.0
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.025 .AND. RNDNUM .LE. 0.05) THEN
		XORIG = 42.47
		YORIG = 6.73
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.05 .AND. RNDNUM .LE. 0.075) THEN
		XORIG = 40.9
		YORIG = 13.29
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.075 .AND. RNDNUM .LE. 0.1) THEN
		XORIG = 38.32
		YORIG = 19.52
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.1 .AND. RNDNUM .LE. 0.125) THEN
		XORIG = 34.79
		YORIG = 25.27
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.125 .AND. RNDNUM .LE. 0.15) THEN
		XORIG = 30.41
		YORIG = 30.4
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.15 .AND. RNDNUM .LE. 0.175) THEN
		XORIG = 25.28
		YORIG = 34.78
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.175 .AND. RNDNUM .LE. 0.2) THEN
		XORIG = 19.53
		YORIG = 38.31
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.2 .AND. RNDNUM .LE. 0.225) THEN
		XORIG = 13.3
		YORIG = 40.89
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.225 .AND. RNDNUM .LE. 0.25) THEN
		XORIG = 6.74
		YORIG = 42.47
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.25 .AND. RNDNUM .LE. 0.275) THEN
		XORIG = 0
		YORIG = 43
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.275 .AND. RNDNUM .LE. 0.3) THEN
		XORIG = -6.73
		YORIG = 42.47
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.3 .AND. RNDNUM .LE. 0.325) THEN
		XORIG = -13.29
		YORIG = 40.9
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.325 .AND. RNDNUM .LE. 0.35) THEN
		XORIG = -19.52
		YORIG = 38.32
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.35 .AND. RNDNUM .LE. 0.375) THEN
		XORIG = -25.27
		YORIG = 34.79
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.375 .AND. RNDNUM .LE. 0.4) THEN
		XORIG = -30.40
		YORIG = 30.41
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.4 .AND. RNDNUM .LE. 0.425) THEN
		XORIG = -34.78
		YORIG = 25.28
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.425 .AND. RNDNUM .LE. 0.45) THEN
		XORIG = -38.31
		YORIG = 19.53
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.45 .AND. RNDNUM .LE. 0.475) THEN
		XORIG = -40.89
		YORIG = 13.3
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.475 .AND. RNDNUM .LE. 0.5) THEN
		XORIG = -42.47
		YORIG = 6.73
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.5 .AND. RNDNUM .LE. 0.525) THEN
		XORIG = -43
		YORIG = 0
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.525 .AND. RNDNUM .LE. 0.55) THEN
		XORIG = -42.47
		YORIG = -6.73
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.55 .AND. RNDNUM .LE. 0.575) THEN
		XORIG = -40.9
		YORIG = -13.29
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.575 .AND. RNDNUM .LE. 0.6) THEN
		XORIG = -38.32
		YORIG = -19.5
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.6 .AND. RNDNUM .LE. 0.625) THEN
		XORIG = -34.79
		YORIG = -25.27		
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.625 .AND. RNDNUM .LE. 0.65) THEN
		XORIG = -30.41
		YORIG = -30.4
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.65 .AND. RNDNUM .LE. 0.675) THEN
		XORIG = -25.28
		YORIG = -34.78
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.675 .AND. RNDNUM .LE. 0.7) THEN
		XORIG = -19.53
		YORIG = -38.31
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.7 .AND. RNDNUM .LE. 0.725) THEN
		XORIG = -13.3
		YORIG = -40.89
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.725 .AND. RNDNUM .LE. 0.75) THEN
		XORIG = -6.73
		YORIG = -42.47
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.75 .AND. RNDNUM .LE. 0.775) THEN
		XORIG = 0
		YORIG = -43
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.775 .AND. RNDNUM .LE. 0.8) THEN
		XORIG = 6.73
		YORIG = -42.47
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.8 .AND. RNDNUM .LE. 0.825) THEN
		XORIG = 13.29
		YORIG = -40.9
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.825 .AND. RNDNUM .LE. 0.85) THEN
		XORIG = 19.52
		YORIG = -38.32
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.85 .AND. RNDNUM .LE. 0.875) THEN
		XORIG = 25.27
		YORIG = -34.79
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.875 .AND. RNDNUM .LE. 0.9) THEN
		XORIG = 30.4
		YORIG = -30.41
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.9 .AND. RNDNUM .LE. 0.925) THEN
		XORIG = 34.78
		YORIG = -25.28
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.925 .AND. RNDNUM .LE. 0.95) THEN
		XORIG = 38.31
		YORIG = -19.53
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20
	ELSE IF (RNDNUM .GT. 0.95 .AND. RNDNUM .LE. 0.975) THEN
		XORIG = 40.89
		YORIG = -13.3
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20

	ELSE 
		XORIG = 42.47		
                YORIG = -6.73
		RMIN = 0.0
		RMAX = 0.3
		ZMIN = -20
		ZMAX = 20	
	ENDIF
		R=SQRT(FLRNDM(DUMMY)) * (RMAX-RMIN) + RMIN
		THETA = FLRNDM(DUMMY) * TWOPIP
    		x = R * cos(THETA) + XORIG
		y = R * sin(THETA) + YORIG
		z = ZMIN + (ZMAX - ZMIN) * FLRNDM(DUMMY)
	      	XFLK (NPFLKA) =x
      		YFLK (NPFLKA) =y
      		ZFLK (NPFLKA) =z
*	WRITE (10, *) x,y,z
*  Calculate the total kinetic energy of the primaries: don't change
*  +-------------------------------------------------------------------*
*  |  (Radioactive) isotope:
      IF ( IJBEAM .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
*  |
*  +-------------------------------------------------------------------*
      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

