*$ 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
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /

       PARAMETER (COEFC = 1.230E+12)
      PARAMETER (COEFK = 2.115E+00)
      PARAMETER (ALPHA = 0.2815E+00)
      PARAMETER (IBNUM = 999)      ! number of energy bins

      DIMENSION ENERGY(0:IBNUM)    ! Initial Energy Of Each Bin
      DIMENSION SPECTR(0:IBNUM)    ! differential fluence function
      DIMENSION CUMUL(0:IBNUM)     ! Cumulative function
*======================================================================*
*                                                                      *
*                 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(*,*)
         WRITE(*,*) "source for Al slab/ monodirectional beam"
         WRITE(*,*) "Version 1 called"
         WRITE(*,*)
 

         DO I = 0, IBNUM
!          WRITE(*,*) I
          ENERGY(I) = 1.0E-3 + 1.0E-3 * I
!          WRITE(*,*) I, ENERGY(I)
          SPECTR(I) = COEFC * COEFK * ALPHA *
     &                EXP(-COEFK * (ENERGY(I) ** ALPHA))* 
     &                ENERGY(I) ** (ALPHA - 1.D0)   
*          WRITE(*,*) I, ENERGY(I), SPECTR(I)   
         END DO

!         WRITE(*,*) "enegry, spectr initialized"
         CUMUL(0) = ZERZER

* calculate the cumulative function

         DO I = 1, IBNUM
          CUMUL(I) = CUMUL(I-1) + HLFHLF * (SPECTR(I) + SPECTR(I-1))
     &               * (ENERGY(I) - ENERGY(I-1))
*          WRITE(*,*) I, ENERGY(I), CUMUL(I)
         END DO
!         WRITE(*,*) "cumulative function initialized", CUMUL(IBNUM)
* normalize the cumulative function to 1

         SUM = CUMUL(IBNUM)
         DO I = 1, IBNUM
           CUMUL(I) = CUMUL(I)/SUM
*           WRITE(*,*) I, ENERGY(I), CUMUL(I)
         END DO 
    
!         WRITE(*,*) "cumulative function norm", CUMUL(IBNUM)

      END IF

*      IJBEAM = 1
      
*  |
*  +-------------------------------------------------------------------*
*  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.
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID = IJBEAM
         ILOFLK (NPFLKA) = IJBEAM
*  |  Flag this is prompt radiation
         LRADDC (NPFLKA) = .FALSE.
      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
* Group number for "low" energy neutrons, set to 0 anyway
      IGROUP (NPFLKA) = 0
* Kinetic energy of the particle (GeV)

!      TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
!
!  Sampling from the distribution 
!
      XI = FLRNDM(DUMMY)
*      WRITE(*,*) " XI", XI
       DO 111 I = 0, IBNUM - 1
          IF ( XI .LT. CUMUL(I+1) ) THEN 
* come here when cumul(i) =< xi < cumul(i+1). Interpolate
*            write(*,*) ENERGY(I), ENERGY(I+1)
            SAMPLE = ENERGY(I) + (ENERGY(I+1) - ENERGY(I))
     &               * (XI - CUMUL(I)) / (CUMUL(I+1) - CUMUL(I))
*            WRITE(*, *) "Sample", SAMPLE
           GO TO 222
          END IF
111    CONTINUE
222    CONTINUE    
       TKEFLK (NPFLKA) = SAMPLE
      
!
!
* Particle momentum
*      PMOFLK (NPFLKA) = PBEAM
!
      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (IONID) ) )
*      WRITE(*, *) " momentum", PMOFLK (NPFLKA)
!
* Cosines (tx,ty,tz)
!      TXFLK  (NPFLKA) = UBEAM
!      TYFLK  (NPFLKA) = VBEAM
!      TZFLK  (NPFLKA) = WBEAM

      TXFLK  (NPFLKA) = 0.D0
      TYFLK  (NPFLKA) = 0.D0
      TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
     &                       - TYFLK (NPFLKA)**2 )
*      WRITE(*, *) " cos z", TZFLK  (NPFLKA)
! 
* Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
* Particle coordinates 
!      XFLK   (NPFLKA) = XBEAM
!      YFLK   (NPFLKA) = YBEAM
!      ZFLK   (NPFLKA) = ZBEAM

      XFLK   (NPFLKA) = -40.D0 + 80.D0 * FLRNDM(DUMMY) 
      YFLK   (NPFLKA) = -40.D0 + 80.D0 * FLRNDM(DUMMY)
      ZFLK   (NPFLKA) = -10.D0
*      WRITE(*, *) " cordinates",XFLK(NPFLKA),YFLK(NPFLKA),ZFLK (NPFLKA) 
*      WRITE(*, *) "***************************************************"
*  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
*      WRITE(*, *) " total kinetic energy", TKESUM
      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
*      WRITE(*, *) " just before exiting tracker" 
      RETURN
*=== End of subroutine Source =========================================*
      END

