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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      DIMENSION ENER(47), CUMER(0:46)
*       DATA ENER / 0.0003d0, 0.0006d0, 0.0010d0, 0.0015d0, 0.0020d0,
*     & 0.0025d0, 0.0030d0, 0.0035d0,0.0040d0, 0.0050d0, 0.0055d0, 
*     & 0.0060d0, 0.0066d0, 0.0070d0, 0.0075d0, 0.0080d0, 0.0085d0,
*     & 0.0090d0, 0.0095d0, 0.0100d0, 0.0105d0, 0.011d0  /
*       DATA CUMER / 0.0d0, 0.0598186049d0, 0.1316009308d0, 0.1964019333,
*     & 0.2512281958d0, 0.3145275176d0, 0.3825606142d0, 0.4732780397,
*     & 0.5455576108d0, 0.6190802952d0, 0.6952980488d0, 0.7388567337,
*     & 0.7808242337d0, 0.8167153967d0, 0.8566044116d0, 0.9014758239,
*     & 0.9401018358d0, 0.9550589733d0, 0.9675199395d0, 0.9824671321,
*     & 0.9949380432d0, 1.0D0 /
       DATA ENER / 0.00d0, 0.00105d0, 0.00116d0,
     & 0.00126d0, 
     & 0.00135d0, 0.00157d0, 0.00171d0, 0.00187d0, 0.00201d0,
     & 0.00219d0, 0.00231d0, 0.00248d0, 0.00259d0,
     & 0.00273d0, 0.00281d0, 0.00289d0, 0.00300d0,
     & 0.00304d0, 0.00314d0, 0.00325d0, 0.00338d0,
     & 0.00355d0, 0.00388d0, 0.00421d0, 0.00455d0,
     & 0.00485d0, 0.00518d0, 0.00544d0, 0.00570d0,
     & 0.00598d0, 0.00626d0, 0.00654d0, 0.00682d0,
     & 0.00706d0, 0.00735d0, 0.00760d0, 0.00787d0,
     & 0.00810d0, 0.00833d0, 0.00863d0, 0.00890d0,
     & 0.00924d0, 0.00964d0, 0.00996d0, 0.01027d0,
     & 0.01065d0, 0.01084 /
       DATA CUMER / 0.0d0, 0.026992d0, 0.048099d0,
     & 0.070133d0, 0.091081d0, 0.110400d0, 0.127472d0,
     & 0.147095d0, 0.169269d0, 0.188885d0, 0.209738d0,
     & 0.231207d0, 0.255072d0, 0.282702d0, 0.312214d0,
     & 0.346600d0, 0.383538d0, 0.423649d0, 0.463371d0,
     & 0.502394d0, 0.537777d0, 0.571532d0, 0.600403d0,
     & 0.630662d0, 0.662927d0, 0.697043d0, 0.730379d0,
     & 0.756589d0, 0.778691d0, 0.793590d0, 0.808174d0,
     & 0.827165d0, 0.847543d0, 0.863815d0, 0.882171d0,
     & 0.900134d0, 0.920727d0, 0.940694d0, 0.955161d0,
     & 0.966603d0, 0.972390d0, 0.977783d0, 0.985104d0,
     & 0.992651d0, 0.997560d0, 0.999520d0, 1.0d0 /
*----------------------------------------------------------------------*
*                                                                      *
*     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              *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      LOGICAL LFIRST
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialization ***
        write(lunout,*)
        write(lunout,'(a,132a)') ("*",i=1,132)
        write(lunout,*)
        write(lunout,*)
        write(lunout,'(a)') "Pu-Be source"
        write(lunout,*)
        write(lunout,*)
        write(lunout,'(a,132a)') ("*",i=1,132)
        write(lunout,*)
       END IF
*------------------------------------------------------------------------*
** the energy of Pu be source is taken from the litrecture             **
          XI=FLRNDM(DUMMY)
          DO 400 K=1,46
         if (XI .gt. CUMER(K-1) .and. XI .LE. CUMER(K) ) THEN
         ENERGY= ENER(K)+
     & (XI-CUMER(K-1))*(ENER(K+1)-ENER(K))/(CUMER(K)-CUMER(K-1))
         GO TO 500
          end if
 400   CONTINUE
       STOP "FAILED TO SAMPLE THE ENERGIES"
 500   CONTINUE
       kount= kount+1
** sample the polar angles
         COSTHE = TWOTWO * FLRNDM(XXX) - ONEONE
         SINTHE=SQRT(ONEONE-COSTHE**2)
**         SAMPLE THE AZIMUTHAL ANGLE
         PHI=TWOPIP * FLRNDM(CCC)
         COSPHI=COS(PHI)      
*  |
*  +-------------------------------------------------------------------*
*  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)
       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)
      TZFLK  (NPFLKA) = COSTHE
      TXFLK  (NPFLKA) = SINTHE * COSPHI
      TYFLK  (NPFLKA) = SQRT( ONEONE - TXFLK (NPFLKA)**2
     &                       - TZFLK (NPFLKA)**2 )
      IF (FLRNDM(COSPHI) .LE. HLFHLF) TYFLK (NPFLKA)= -TYFLK (NPFLKA)
*  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

