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

        INCLUDE '(DBLPRC)'
        INCLUDE '(DIMPAR)'
        INCLUDE '(IOUNIT)'
*
        INCLUDE '(BEAMCM)'
        INCLUDE '(FHEAVY)'
        INCLUDE '(FLKSTK)'
        INCLUDE '(IOIOCM)'
        INCLUDE '(LTCLCM)'
        INCLUDE '(PAPROP)'
        INCLUDE '(SOURCM)'
        INCLUDE '(SUMCOU)'
*
        LOGICAL LFIRST
*
        PARAMETER (NCF = 33, NAMBE = 22)
        DIMENSION CF(0:NCF), AMBE(0:NAMBE), ENERCF(0:NCF), ENERAM(0:NAMBE)
        DIMENSION CUMUCF(0:NCF), CUMUAM(0:NAMBE)
        SAVE LFIRST
        DATA LFIRST / .TRUE. /
        DATA ENERCF / 4.64d-6,
       & 1.00d-5, 1.25d-5, 1.58d-5, 1.99d-5, 2.51d-5, 3.16d-5,
       & 3.98d-5, 5.01d-5, 6.30d-5, 7.94d-5, 1.00d-4, 1.25d-4, 1.58d-4,
       & 1.99d-4, 2.51d-4, 3.16d-4, 3.98d-4, 5.01d-4, 6.30d-4, 7.94d-4,
       & 1.00d-3, 1.25d-3, 1.58d-3, 1.99d-3, 2.51d-3, 3.16d-3, 3.98d-3,
       & 5.01d-3, 6.30d-3, 7.94d-3, 1.00d-2, 1.58d-2, 2.00d-2 /
        DATA ENERAM / 7.94d-5, 1.00d-4, 1.25d-4, 1.58d-4,
       & 1.99d-4, 2.51d-4, 3.16d-4, 3.98d-4, 5.01d-4, 6.30d-4, 7.94d-4,
       & 1.00d-3, 1.25d-3, 1.58d-3, 1.99d-3, 2.51d-3, 3.16d-3, 3.98d-3,
       & 5.01d-3, 6.30d-3, 7.94d-3, 1.00d-2, 1.58d-2 /
* spectra are in n/cm2 per unit lethargy!
        DATA CF / 0.d0,
       & 9.22d-4, 1.29d-3, 1.79d-3, 2.48d-3, 3.39d-3, 4.64d-3,
       & 6.38d-3, 8.81d-3, 1.21d-2, 1.69d-2, 2.37d-2, 3.26d-2, 4.45d-2,
       & 6.05d-2, 8.22d-2, 1.09d-1, 1.45d-1, 1.89d-1, 2.41d-1, 3.00d-1,
       & 3.59d-1, 4.14d-1, 4.52d-1, 4.62d-1, 4.34d-1, 3.66d-1, 2.71d-1,
       & 1.69d-1, 8.61d-2, 3.35d-2, 4.27d-3, 4.31d-4, 0.d0 /
        DATA AMBE / 0.d0,
       &   1.66d-2, 2.21d-2, 2.87d-2, 3.67d-2, 4.65d-2, 5.77d-2,
       &   7.06d-2, 8.48d-2, 9.61d-2, 1.06d-1, 1.18d-1, 1.27d-1, 1.81d-1,
       &   2.43d-1, 4.21d-1, 5.71d-1, 6.86d-1, 6.50d-1, 5.78d-1, 1.66d-1,
       &   1.72d-2, 0.d0 /
        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,*)
           IF(NINT(WHASOU(1)) .EQ. 1) THEN
              WRITE(LUNOUT,'(A)') "@ Source: Cf-252"
           ELSE IF(NINT(WHASOU(1)) .EQ. 2) THEN
              WRITE(LUNOUT,'(A)') "@ Source: Am-Be"
           ELSE
              STOP "Incorrect WHASOU(1)"
           END IF
           WRITE(LUNOUT,*)
           WRITE(LUNOUT,*)
           WRITE(LUNOUT,'(A,132A)') ("*",I=1,132)
           WRITE(LUNOUT,*)
           if(nint(whasou(1)) .eq. 1) then
*           Transform energies to lethargies
              do i = 0, ncf
                 enercf(i) = log(enercf(i))
              end do
              cumucf(0) = 0.d0
*           prepare the cumulative function
              do i = 1, ncf
                 cumucf(i) = cumucf(i-1) +
       &                     (cf(i-1)+cf(i))*(enercf(i)-enercf(i-1))*0.5d0
              end do
*               WRITE(88,*) "Lethargies, Cumulative, non normalized"
*               DO I = 0, NCF
*                  WRITE(88,'(I3,1X,1P,2G15.9)') I, ENERCF(I), CUMUCF(I)
*               END DO
*           normalize the cumulative function
              cunorm = cumucf(ncf)
              do i = 1, ncf
                 cumucf(i) = cumucf(i)/cunorm
              end do
*               WRITE(88,*) "Lethargies, Cumulative, normalized"
*               DO I = 0, NCF
*                  WRITE(88,'(I3,1X,1P,2G15.9)') I, ENERCF(I), CUMUCF(I)
*               END DO
           else if(nint(whasou(1)) .eq. 2) then
*           Transform energies to lethargies
              do i = 0, nambe
                 eneram(i) = log(eneram(i))
              end do
              cumuam(0) = 0.d0
*           prepare the cumulative function
              do i = 1, nambe
               cumuam(i) = cumuam(i-1) +
       &                 (ambe(i-1)+ambe(i))*(eneram(i)-eneram(i-1))*0.5d0
              end do
*           normalize the cumulative function
              cunorm = cumuam(nambe)
              do i = 1, nambe
                 cumuam(i) = cumuam(i)/cunorm
              end do
           end if
*         kount = 0
*      WRITE(88,*)
*     & "Sampled: xi, cum(i-1), cum(i), leth(i-1), leth(i), leth, energy"
        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
*      kount = kount + 1
*      if(kount .gt. 100) nomore = 1
        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). Sample from the cumulative
        xi = flrndm(dummy)
        if(nint(whasou(1)) .eq. 1) then
           do i = 1, ncf
              if(xi .gt. cumucf(i-1) .and. xi .le. cumucf(i)) then
                 enkin = enercf(i-1)+(enercf(i)-enercf(i-1))*
       &         (sqrt(xi) - sqrt(cumucf(i-1)))/
       &                               (sqrt(cumucf(i))-sqrt(cumucf(i-1)))
*               WRITE(88,'(3G12.4,1P,4G12.4)')
*     & XI,CUMUCF(I-1),CUMUCF(I),ENERCF(I-1),ENERCF(I),ENKIN,EXP(ENKIN)
                 go to 1
              end if
           end do
        else if(nint(whasou(1)) .eq. 2) then
           do i = 1, nambe
              if(xi .gt. cumuam(i-1) .and. xi .le. cumuam(i)) then
                 enkin = eneram(i-1)+(eneram(i)-eneram(i-1))*
       &        (sqrt(xi) - sqrt(cumuam(i-1)))/
       &                               (sqrt(cumuam(i))-sqrt(cumuam(i-1)))
                 go to 1
              end if
           end do
        end if
   1    continue
*     back from lethargy to energy
        tkeflk (npflka) = exp(enkin)
*     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
* 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


