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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2006      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     New source for FLUKA9x-FLUKA200x:                                *
*                                                                      *
*     Created on 07 january 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on 03-mar-06     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.            *
*                                                                      *
*----------------------------------------------------------------------*

*----------------------------------------------------------------------*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      LOGICAL LFIRST
*
*  Version for isotropic Heavy ion source
*
      DIMENSION  CUMENE(0:58), ENELOW(0:58), ENEUP(0:58) 
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
*    Heavy ion energy group(Gev/n)
      DATA ENELOW /
     & 0.00000D-03, 1.00000D-04, 2.80808D-04, 4.61616D-04, 6.42424D-04,
     & 8.23232D-04, 1.00404D-03, 1.18485D-03, 1.36566D-03, 1.54646D-03,
     & 1.72727D-03, 1.90808D-03, 2.08889D-03, 2.26970D-03, 2.45051D-03,
     & 2.63131D-03, 2.81212D-03, 2.99293D-03, 3.17374D-03, 3.35455D-03,
     & 3.53535D-03, 3.71616D-03, 3.89697D-03, 4.07778D-03, 4.25859D-03,
     & 4.43939D-03, 4.62020D-03, 4.80101D-03, 4.98182D-03, 5.16263D-03,
     & 5.34343D-03, 5.52424D-03, 5.70505D-03, 5.88586D-03, 6.06667D-03,
     & 6.24747D-03, 6.42828D-03, 6.60909D-03, 6.78990D-03, 6.97071D-03,
     & 7.15152D-03, 7.33232D-03, 7.51313D-03, 7.69394D-03, 7.87475D-03,
     & 8.05556D-03, 8.23636D-03, 8.41717D-03, 8.59798D-03, 8.77879D-03,
     & 8.95960D-03, 9.14040D-03, 9.32121D-03, 9.50202D-03, 9.68283D-03,
     & 9.86364D-03, 1.00444D-02, 1.02253D-02, 1.04061D-02/

      DATA ENEUP /
     & 1.00000D-04, 2.80808D-04, 4.61616D-04, 6.42424D-04, 8.23232D-04,
     & 1.00404D-03, 1.18485D-03, 1.36566D-03, 1.54646D-03, 1.72727D-03,
     & 1.90808D-03, 2.08889D-03, 2.26970D-03, 2.45051D-03, 2.63131D-03,
     & 2.81212D-03, 2.99293D-03, 3.17374D-03, 3.35455D-03, 3.53535D-03,
     & 3.71616D-03, 3.89697D-03, 4.07778D-03, 4.25859D-03, 4.43939D-03,
     & 4.62020D-03, 4.80101D-03, 4.98182D-03, 5.16263D-03, 5.34343D-03,
     & 5.52424D-03, 5.70505D-03, 5.88586D-03, 6.06667D-03, 6.24747D-03,
     & 6.42828D-03, 6.60909D-03, 6.78990D-03, 6.97071D-03, 7.15152D-03,
     & 7.33232D-03, 7.51313D-03, 7.69394D-03, 7.87475D-03, 8.05556D-03,
     & 8.23636D-03, 8.41717D-03, 8.59798D-03, 8.77879D-03, 8.95960D-03,
     & 9.14040D-03, 9.32121D-03, 9.50202D-03, 9.68283D-03, 9.86364D-03,
     & 1.00444D-02, 1.02253D-02, 1.04061D-02, 1.05869D-02/

*    Normalized spectrum
      DATA CUMENE /
     & 3.22531E-02, 1.09893E+02, 2.34973E+02, 3.56235E+02, 4.64787E+02,
     & 5.57476E+02, 6.34223E+02, 6.96450E+02, 7.46176E+02, 7.85523E+02,
     & 8.16469E+02, 8.40738E+02, 8.59775E+02, 8.74754E+02, 8.86612E+02,
     & 8.96083E+02, 9.03736E+02, 9.10006E+02, 9.15227E+02, 9.19648E+02,
     & 9.23458E+02, 9.26798E+02, 9.29774E+02, 9.32463E+02, 9.34923E+02,
     & 9.37197E+02, 9.39317E+02, 9.41306E+02, 9.43181E+02, 9.44956E+02,
     & 9.46640E+02, 9.48239E+02, 9.49761E+02, 9.51208E+02, 9.52585E+02,
     & 9.53893E+02, 9.55134E+02, 9.56312E+02, 9.57426E+02, 9.58478E+02,
     & 9.59469E+02, 9.60400E+02, 9.61273E+02, 9.62087E+02, 9.62844E+02,
     & 9.63544E+02, 9.64188E+02, 9.64777E+02, 9.65312E+02, 9.65792E+02,
     & 9.66219E+02, 9.66593E+02, 9.66915E+02, 9.67185E+02, 9.67404E+02,
     & 9.67572E+02, 9.67691E+02, 9.67759E+02, 9.67779E+02/
    
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialization ***
	 XI = FLRNDM(DUMMY)
      DO K = 1 ,100
         IF (XI .LE. CUMENE(K)) THEN
            ENERGY = ((ENEUP(K) + ENELOW(K))/ TWOTWO) 
            GO TO 600
         END IF
      END DO
 600   CONTINUE 
      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 = 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
* 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)
      TKEFLK (NPFLKA) = ENERGY
*     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)
      COSTH=ONEONE-2.D0*RAN(ISEED)
      IF(COSTH.GT.ONEONE) COSTH=ONEONE
      IF(COSTH.LT.-ONEONE) COSTH=-ONEONE
      ANG=TWOPIP*RAN2(ISEED2)
      UBEAM=SQRT(ONEONE-COSTH*COSTH)*COS(ANG)
      VBEAM=SQRT(ONEONE-COSTH*COSTH)*SIN(ANG)
      WBEAM=COSTH
      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) = 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

