*$ 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 (NTR = 53, NAMBE = 22)
      DIMENSION TR(0:NTR), AMBE(0:NAMBE), ENERTR(0:NTR), ENERAM(0:NAMBE) 
      DIMENSION CUMUTR(0:NTR), CUMUAM(0:NAMBE)
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
      DATA ENERTR / 1.0d-15,
     & 1.0d-14, 5.00d-14, 1.00d-13, 5.00d-13, 1.00d-12, 2.00d-12,
     & 3.00d-12, 4.00d-12, 5.00d-12, 6.00d-12, 7.00d-12, 8.00d-12, 
     & 9.00d-12, 1.00d-11, 2.00d-11, 3.00d-11, 4.00d-11, 5.00d-11, 
     & 6.00d-11, 7.00d-11, 8.00d-11, 9.00d-11, 1.00d-10, 1.25d-10, 
     & 2.00d-10, 3.06d-10, 4.00d-10, 5.00d-10, 6.00d-10, 7.00d-10, 
     & 8.00d-10, 9.00d-10, 1.00d-9, 2.00d-9, 3.00d-9, 4.00d-9, 5.00d-9, 
     & 6.00d-9, 7.00d-9, 8.00d-9, 9.00d-9, 1.00d-8, 5.00d-8, 1.00d-7, 
     & 1.00d-6, 5.00d-6, 1.00d-5, 1.50d-4, 2.00d-4, 5.00d-4, 1.00d-3,
     & 2.00d-3, 3.00d-3 /
      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 per unit energy (GeV)!
      DATA TR / 0.d0,
     & 9.71d+6, 4.85d+7, 9.68d+7, 4.76d+8, 9.34d+8, 1.80d+9,
     & 2.59d+9, 3.32d+9, 3.99d+9, 4.60d+9, 5.16d+9, 5.66d+9, 6.13d+9,
     & 6.54d+9, 8.81d+9, 8.90d+9, 7.99d+9, 6.73d+9, 5.44d+9, 4.27d+9,
     & 3.29d+9, 2.49d+9, 1.87d+9, 9.83d+8, 7.16d+7, 4.77d+7, 3.58d+7, 
     & 2.86d+7, 2.39d+7, 2.05d+7, 1.79d+7, 1.59d+7, 1.43d+7, 7.16d+6,
     & 4.77d+6, 3.58d+6, 2.86d+6, 2.39d+6, 2.05d+6, 1.79d+6, 1.59d+6,
     & 1.42d+6, 2.86d+5, 1.43d+5, 1.43d+4, 2.86d+3, 1.43d+3, 9.54d+1,
     & 7.16d+1, 2.86d+1, 1.43d+1, 0.d09, 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: Triga"
         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, ntr
*               enertr(i) = log(enertr(i))
*            end do
            cumutr(0) = 0.d0
*           prepare the cumulative function
            do i = 1, ntr 
               cumutr(i) = cumutr(i-1) +
     &                     (tr(i-1)+tr(i))*(enertr(i)-enertr(i-1))*0.5d0
            end do
*               WRITE(88,*) "Lethargies, Cumulative, non normalized"
*               DO I = 0, NTR
*                  WRITE(88,'(I3,1X,1P,2G15.9)') I, ENERTR(I), CUMUTR(I)
*               END DO
*           normalize the cumulative function
            cunorm = cumutr(ntr)
            do i = 1, ntr 
               cumutr(i) = cumutr(i)/cunorm
            end do
*               WRITE(88,*) "Lethargies, Cumulative, normalized"
*               DO I = 0, NTR
*                  WRITE(88,'(I3,1X,1P,2G15.9)') I, ENERTR(I), CUMUTR(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: dont 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, ntr
            if(xi .gt. cumutr(i-1) .and. xi .le. cumutr(i)) then
               enkin = enertr(i-1)+(enertr(i)-enertr(i-1))*
     &         (sqrt(xi) - sqrt(cumutr(i-1)))/
     &                              (sqrt(cumutr(i))-sqrt(cumutr(i-1)))
*               WRITE(88,'(3G12.4,1P,4G12.4)') 
*     & XI,CUMUTR(I-1),CUMUTR(I),ENERTR(I-1),ENERTR(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) = 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: dont 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
