*$ 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
*
c   defining and saving spectrum arrays
      DIMENSION ENEPOI(0:1000),ENEPRO(0:1000),ENECUM(0:1000)
      SAVE ENEPOI, ENEPRO, ENECUM
c   saving spectrum dimension
      SAVE IMAX
*
      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 ***
         CALL OAUXFI('spectrum.dat',LUNRDB,'OLD',IERR)
c   reading spectrum
         DO I=0,10000
           READ(LUNRDB,*,END=1972) ENEPOI(I),ENEPRO(I)
           IMAX=I
c    WRITE(*, 780) ENEPOI(I),' ',ENEPRO(I)
c    WRITE(*, 782) 'IMAX= ',IMAX
         ENDDO 
         STOP ' spectrum reading uncomplete!'
c  
 1972    CONTINUE
         ENECUM(0)=ZERZER
c   building cumulative spectrum
         DO I=1,IMAX
           ENECUM(I)=ENEPRO(I-1)+ENECUM(I-1)
c    WRITE(*, 781) 'ENECUM(I)= ', ENECUM(I)
         ENDDO
c   normalizing cumulative spectrum
         DO I=1,IMAX
           ENECUM(I)=ENECUM(I)/ENECUM(IMAX)
c    WRITE(*, 781) 'ENECUM(I)_norm= ', ENECUM(I)
         ENDDO
      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.
*  |
*  +-------------------------------------------------------------------*
*  |  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
c
c   sampling from the normalized cumulative spectrum
      XYZ=FLRNDM(DUMMY)
c      WRITE(*, 781) 'XYZ= ', XYZ
c   the sampled energy lies in the bin I (between ENEPOI(I-1) and ENEPOI(I))
c   now determining the energy inside the bin I according to a linear spectrum
      DO I=1,IMAX
        IF((XYZ.GE.ENECUM(I-1)).AND.(XYZ.LT.ENECUM(I))) THEN
           ESAMPLE=ENEPOI(I-1)+FLRNDM(DUMMY)*(ENEPOI(I)-ENEPOI(I-1))
c          WRITE(*, 781) 'ESAMPLE= ', ESAMPLE
c    WRITE(*, 781) 'ENEPOI(I-1)= ', ENEPOI(I-1)
c    WRITE(*, 781) 'ENEPOI(I)= ', ENEPOI(I)
        ENDIF
      ENDDO
* Kinetic energy of the particle (GeV)
      TKEFLK (NPFLKA) = ESAMPLE
* Particle momentum
      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (IONID) ) )
* Cosines (tx,ty,tz)
c   uniform angular distribution (isotropic source)
      WBEAM = -1.+2.*FLRNDM(DUMMY)
      VBEAM = SIN( XYZ*2.*PIPIPI ) * 
     &                       SQRT( 1. - WBEAM*WBEAM)
      UBEAM = COS( XYZ*2.*PIPIPI ) * 
     &                       SQRT( 1. - WBEAM*WBEAM)
      TXFLK  (NPFLKA) = UBEAM
      TYFLK  (NPFLKA) = VBEAM
      TZFLK  (NPFLKA) = WBEAM
c       WRITE(*, 998) ' UBEAM= ', UBEAM, ' VBEAM= ', VBEAM,  
c      &              ' WBEAM= ', WBEAM, ' END ', 
c      & UBEAM*UBEAM+VBEAM*VBEAM+WBEAM*WBEAM, ' ', XYZ
* Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
* Particle coordinates
c   uniform 3D cylindrical distribution where:
c   WHASOU(1) = Rmin; WHASOU(2) = Rmax;
c   WHASOU(3) = Hmin; WHASOU(4) = Hmax;
c   WHASOU(5) = THTmin; WHASOU(6) = THTmax;  !!! IN DEGREES
c   conversion in radians   
      DEGmin = WHASOU(5) * PIPIPI / 180.    
      DEGmax = WHASOU(6) * PIPIPI / 180.    
c
      Tht = DEGmin + (DEGmax - DEGmin) * FLRNDM(DUMMY)
      Hal = WHASOU(3) + (WHASOU(4) - WHASOU(3)) * FLRNDM(DUMMY)
      Rho = 0.5* (WHASOU(1)+SQRT(WHASOU(1)**2+
     &        4*FLRNDM(DUMMY)*WHASOU(2)*(WHASOU(2)-WHASOU(1))))
c       
      XBEAM = Rho * COS( Tht )      
      XFLK   (NPFLKA) = XBEAM
      YBEAM = Rho * SIN( Tht )   
      YFLK   (NPFLKA) = YBEAM
      ZBEAM = Hal    
      ZFLK   (NPFLKA) = ZBEAM
c      WRITE(*, 997) 
c     & ' WHASOU(1)= ', WHASOU(1), ' WHASOU(2)= ', WHASOU(2),  
c     & ' WHASOU(3)= ', WHASOU(3), ' WHASOU(4)= ', WHASOU(4),
c     & ' WHASOU(5)= ', WHASOU(5), ' WHASOU(6)= ', WHASOU(6),
c     & ' DEGmin= ', DEGmin, ' DEGmax= ', DEGmax,
c     & ' END'
      WRITE(*, 999) ' XBEAM= ', XBEAM, ' YBEAM= ', YBEAM,  
     &              ' ZBEAM= ', ZBEAM, ' END'
c 
  780 FORMAT (E12.6, A, E12.6)
  781 FORMAT (A, E12.6)
  782 FORMAT (A, I4)
  999 FORMAT (A, E12.6, A, E12.6, A, E12.6, A)
  998 FORMAT (A, E12.6, A, E12.6, A, E12.6, A, E12.6, A, E12.6)
  997 FORMAT (A, E12.6, A, E12.6, 
     &        A, E12.6, A, E12.6, 
     &        A, E12.6, A, E12.6,
     &        A, E12.6, A, E12.6,
     &        A)  
c
*  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
