*=== 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)'

*
       LOGICAL LFIRST
*
       SAVE LFIRST
       DATA LFIRST / .TRUE. /
       DOUBLE PRECISION ZULIM, ZLLIM, RLIMIT, RANGLE, SPHI, STHE, SRV

*======================================================================*
*                                                                      *
*                 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
*  |
*  +-------------------------------------------------------------------*
       ZLLIM = 1982.125
       ZULIM = 0.01
       RLIMIT = 0.7
       RANGLE = 2*PIPIPI

*  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
* LSTOPP is the stack counter: of course any time source is called it
* must be =0
       IJBEAM = -1
       LSTOPP = LSTOPP + 1
*     Weight of optical photon
       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
*
*     Insert in POPTPH (LSTOPP) the proper energy for optical photon
*
       WVL2EN = 600E-9
       POPTPH (LSTOPP) = (((6.626D-34*299792458)/WVL2EN)/ELCMKS)/1D9
       DONEAR (LSTOPP) = ZERZER

* below should be a proper isotropic distribution
        STHE = ACOS(2*FLRNDM(TKESUM) - 1)
        SPHI = 2*PIPIPI*FLRNDM(TKESUM)
        SRV = 1
* convert random spherical coordinates to cartesian...

        TXOPPH (LSTOPP) = SRV*COS(SPHI)*SIN(STHE)

        TYOPPH (LSTOPP) = SRV*SIN(SPHI)*SIN(STHE)
        TZOPPH (LSTOPP) = SRV*COS(STHE)

* Cosine normlization to prevent a separate routine from normalizing
       TNORML=(SQRT(TXOPPH(LSTOPP)**2+TYOPPH(LSTOPP)**2
      &       +TZOPPH(LSTOPP)**2))

       TXOPPH (LSTOPP) = TXOPPH (LSTOPP) / TNORML
       TYOPPH (LSTOPP) = TYOPPH (LSTOPP) / TNORML
       TZOPPH (LSTOPP) = TZOPPH (LSTOPP) / TNORML


*     Injection coordinates of optical photon
*  Particle coordinates  FLRNDM(TKESUM)
       RLIMIT = 0.7*SQRT(FLRNDM(TKESUM))
       RANGLE = FLRNDM(TKESUM)*4*PIPIPI
       XOPTPH (LSTOPP)   = RLIMIT*COS(RANGLE)
       YOPTPH (LSTOPP) = RLIMIT*SIN(RANGLE)
       ZOPTPH (LSTOPP)  = ZLLIM + ZULIM*FLRNDM(TKESUM)

*  Set-up the polarization vector
       TXPOPP (LSTOPP) = -TWOTWO
       TYPOPP (LSTOPP) = ZERZER
       TZPOPP (LSTOPP) = ZERZER
*  age
       AGOPPH (LSTOPP) = ZERZER
* total path
       CMPOPP (LSTOPP) = ZERZER
* Particle generation
       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
       TKESUM = TKESUM + POPTPH (LSTOPP) * WTOPPH (LSTOPP)
*
       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
       CALL SOEVSV
       RETURN
*=== End of subroutine Source =========================================*
       END

