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

       INCLUDE '(DBLPRC)'
       INCLUDE '(DIMPAR)'
       INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2009      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 08-feb-09     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)'
       INCLUDE '(OPPHST)'
       INCLUDE '(TRACKR)'
*
       INCLUDE 'gbat.inc'
*
       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 ***
       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
       IJBEAM = -1
       IF (IJBEAM .EQ. -1) THEN
          LSTOPP = LSTOPP + 1
          WTOPPH (LSTOPP) = ONEONE
          WEIPRI = WEIPRI + WTOPPH (LSTOPP)
          NUMOPH = NUMOPH + 1
          IF ( NUMOPH .GT. 1000000000 ) THEN
             MUMOPH = MUMOPH + 1
             NUMOPH = NUMOPH - 1000000000
          END IF
          WOPTPH = WOPTPH + ONEONE
          POPTPH (LSTOPP) = 4.D-09
          DONEAR (LSTOPP) = ZERZER
* sample initial position
          XGBOPT = -0.5D+00 + ONEONE * FLRNDM(XDUMMY)
          YGBOPT = -5.0D+00 + 10.0D+00 * FLRNDM(XDUMMY)
          ZGBOPT = 20.000D+00
          XOPTPH (LSTOPP) = XGBOPT
          YOPTPH (LSTOPP) = YGBOPT
          ZOPTPH (LSTOPP) = ZGBOPT
* sample direction cosines
          CALL SFECFE ( SINPHI, COSPHI )
          COSTHE = FLRNDM(XDUMMY)
          SINTHE = SQRT ( ( ONEONE - COSTHE ) * ( ONEONE + COSTHE ) )
          UGBOPT = COSPHI * SINTHE
          VGBOPT = SINPHI * SINTHE
          WGBOPT = - COSTHE
          TXOPPH (LSTOPP) = UGBOPT
          TYOPPH (LSTOPP) = VGBOPT
          TZOPPH (LSTOPP) = WGBOPT
*  |  Set-up the polarization vector: the following for linear pola-
*  |  rization lying in the plane defined by the photon and the
*  |  original particle direction:
          TXPOPP (LSTOPP) = -TWOTWO
          TYPOPP (LSTOPP) = ZERZER
          TZPOPP (LSTOPP) = ZERZER
          AGOPPH (LSTOPP) = ZERZER
          CMPOPP (LSTOPP) = ZERZER
          LOOPPH (LSTOPP) = 1
          LOUOPP (LSTOPP) = LLOUSE
          DO 2100 ISPR = 1, MKBMX1
             SPAROK (ISPR,LSTOPP) = ZERZER
  2100    CONTINUE
          DO 2200 ISPR = 1, MKBMX2
             ISPORK (ISPR,LSTOPP) = 0
  2200    CONTINUE
          CALL GEOCRS ( TXOPPH (LSTOPP), TYOPPH (LSTOPP), TZOPPH (LSTOPP)
      $        )
          CALL GEOREG ( XOPTPH (LSTOPP), YOPTPH (LSTOPP), ZOPTPH (LSTOPP)
      $        ,NREGOP (LSTOPP), IDISC )
*  Do not change these cards:
          CALL GEOHSM ( IHSPNT, 1, -11, MLATTC )
          NLATOP (LSTOPP) = MLATTC
       ELSE
          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
* 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) = 6.626E-20
          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) ) )
*      PMOFLK (NPFLKA) = TKEFLK (NPFLKA)
* Cosines (tx,ty,tz)
          TXFLK  (NPFLKA) = -ONEONE
          TYFLK  (NPFLKA) = ZERZER
          TZFLK  (NPFLKA) = ZERZER
*     TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
*&                        - TYFLK (NPFLKA)**2 )
* Polarization cosines:
          TXPOL  (NPFLKA) = -TWOTWO
          TYPOL  (NPFLKA) = +ZERZER
          TZPOL  (NPFLKA) = +ZERZER
* Particle coordinates
          XFLK   (NPFLKA) = 2.0000
          YFLK   (NPFLKA) = -5.D+00 + 10.D+00*FLRNDM(XDUMMY)
          ZFLK   (NPFLKA) = 19.9999
*  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
       ENDIF
       CALL SOEVSV
       RETURN
*=== End of subroutine Source =========================================*
       END


