*$ 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
*
*  Version for isotropic Heavy ion source
*
       DIMENSION  CUMENE(0:100), ENELOW(0:99), ENEUP(0:99)
       SAVE LFIRST
       DATA LFIRST / .TRUE. /
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
*    Heavy ion energy group(Gev/n)
       DATA ENELOW /
      &  0.00000E+00, 1.01010E-03, 2.02020E-03, 3.03030E-03, 4.04040E-03,
      &  5.05051E-03, 6.06061E-03, 7.07071E-03, 8.08081E-03, 9.09091E-03,
      &  1.01010E-02, 1.11111E-02, 1.21212E-02, 1.31313E-02, 1.41414E-02,
      &  1.51515E-02, 1.61616E-02, 1.71717E-02, 1.81818E-02, 1.91919E-02,
      &  2.02020E-02, 2.12121E-02, 2.22222E-02, 2.32323E-02, 2.42424E-02,
      &  2.52525E-02, 2.62626E-02, 2.72727E-02, 2.82828E-02, 2.92929E-02,
      &  3.03030E-02, 3.13131E-02, 3.23232E-02, 3.33333E-02, 3.43434E-02,
      &  3.53535E-02, 3.63636E-02, 3.73737E-02, 3.83838E-02, 3.93939E-02,
      &  4.04040E-02, 4.14141E-02, 4.24242E-02, 4.34343E-02, 4.44444E-02,
      &  4.54545E-02, 4.64646E-02, 4.74747E-02, 4.84848E-02, 4.94949E-02,
      &  5.05051E-02, 5.15152E-02, 5.25253E-02, 5.35354E-02, 5.45455E-02,
      &  5.55556E-02, 5.65657E-02, 5.75758E-02, 5.85859E-02, 5.95960E-02,
      &  6.06061E-02, 6.16162E-02, 6.26263E-02, 6.36364E-02, 6.46465E-02,
      &  6.56566E-02, 6.66667E-02, 6.76768E-02, 6.86869E-02, 6.96970E-02,
      &  7.07071E-02, 7.17172E-02, 7.27273E-02, 7.37374E-02, 7.47475E-02,
      &  7.57576E-02, 7.67677E-02, 7.77778E-02, 7.87879E-02, 7.97980E-02,
      &  8.08081E-02, 8.18182E-02, 8.28283E-02, 8.38384E-02, 8.48485E-02,
      &  8.58586E-02, 8.68687E-02, 8.78788E-02, 8.88889E-02, 8.98990E-02,
      &  9.09091E-02, 9.19192E-02, 9.29293E-02, 9.39394E-02, 9.49495E-02,
      &  9.59596E-02, 9.69697E-02, 9.79798E-02, 9.89899E-02, 1.00000E-01/
       DATA ENEUP /
      &  1.01010E-03, 2.02020E-03, 3.03030E-03, 4.04040E-03, 5.05051E-03,
      &  6.06061E-03, 7.07071E-03, 8.08081E-03, 9.09091E-03, 1.01010E-02,
      &  1.11111E-02, 1.21212E-02, 1.31313E-02, 1.41414E-02, 1.51515E-02,
      &  1.61616E-02, 1.71717E-02, 1.81818E-02, 1.91919E-02, 2.02020E-02,
      &  2.12121E-02, 2.22222E-02, 2.32323E-02, 2.42424E-02, 2.52525E-02,
      &  2.62626E-02, 2.72727E-02, 2.82828E-02, 2.92929E-02, 3.03030E-02,
      &  3.13131E-02, 3.23232E-02, 3.33333E-02, 3.43434E-02, 3.53535E-02,
      &  3.63636E-02, 3.73737E-02, 3.83838E-02, 3.93939E-02, 4.04040E-02,
      &  4.14141E-02, 4.24242E-02, 4.34343E-02, 4.44444E-02, 4.54545E-02,
      &  4.64646E-02, 4.74747E-02, 4.84848E-02, 4.94949E-02, 5.05051E-02,
      &  5.15152E-02, 5.25253E-02, 5.35354E-02, 5.45455E-02, 5.55556E-02,
      &  5.65657E-02, 5.75758E-02, 5.85859E-02, 5.95960E-02, 6.06061E-02,
      &  6.16162E-02, 6.26263E-02, 6.36364E-02, 6.46465E-02, 6.56566E-02,
      &  6.66667E-02, 6.76768E-02, 6.86869E-02, 6.96970E-02, 7.07071E-02,
      &  7.17172E-02, 7.27273E-02, 7.37374E-02, 7.47475E-02, 7.57576E-02,
      &  7.67677E-02, 7.77778E-02, 7.87879E-02, 7.97980E-02, 8.08081E-02,
      &  8.18182E-02, 8.28283E-02, 8.38384E-02, 8.48485E-02, 8.58586E-02,
      &  8.68687E-02, 8.78788E-02, 8.88889E-02, 8.98990E-02, 9.09091E-02,
      &  9.19192E-02, 9.29293E-02, 9.39394E-02, 9.49495E-02, 9.59596E-02,
      &  9.69697E-02, 9.79798E-02, 9.89899E-02, 1.00000E-01, 1.01010E-01/
*    Normalized spectrum
       DATA CUMENE /
      &  0.00000E+00,
      &  2.33568E-01, 3.45591E-01, 4.28197E-01, 4.93572E-01, 5.47237E-01,
      &  5.92330E-01, 6.30854E-01, 6.64180E-01, 6.93297E-01, 7.18945E-01,
      &  7.41689E-01, 7.61976E-01, 7.80163E-01, 7.96537E-01, 8.11337E-01,
      &  8.24761E-01, 8.36975E-01, 8.48120E-01, 8.58315E-01, 8.67664E-01,
      &  8.76255E-01, 8.84166E-01, 8.91465E-01, 8.98211E-01, 9.04455E-01,
      &  9.10244E-01, 9.15620E-01, 9.20617E-01, 9.25269E-01, 9.29605E-01,
      &  9.33651E-01, 9.37431E-01, 9.40965E-01, 9.44273E-01, 9.47373E-01,
      &  9.50279E-01, 9.53007E-01, 9.55570E-01, 9.57979E-01, 9.60245E-01,
      &  9.62379E-01, 9.64389E-01, 9.66285E-01, 9.68073E-01, 9.69761E-01,
      &  9.71355E-01, 9.72862E-01, 9.74287E-01, 9.75636E-01, 9.76912E-01,
      &  9.78121E-01, 9.79267E-01, 9.80354E-01, 9.81384E-01, 9.82362E-01,
      &  9.83290E-01, 9.84172E-01, 9.85010E-01, 9.85807E-01, 9.86564E-01,
      &  9.87285E-01, 9.87971E-01, 9.88624E-01, 9.89245E-01, 9.89838E-01,
      &  9.90403E-01, 9.90941E-01, 9.91454E-01, 9.91944E-01, 9.92411E-01,
      &  9.92858E-01, 9.93284E-01, 9.93691E-01, 9.94079E-01, 9.94451E-01,
      &  9.94806E-01, 9.95146E-01, 9.95470E-01, 9.95781E-01, 9.96078E-01,
      &  9.96363E-01, 9.96635E-01, 9.96896E-01, 9.97146E-01, 9.97386E-01,
      &  9.97615E-01, 9.97835E-01, 9.98046E-01, 9.98248E-01, 9.98442E-01,
      &  9.98628E-01, 9.98806E-01, 9.98978E-01, 9.99142E-01, 9.99300E-01,
      &  9.99451E-01, 9.99597E-01, 9.99737E-01, 9.99871E-01, 1.00000E+00/
*  +-------------------------------------------------------------------*
       NOMORE = 0
*  +-------------------------------------------------------------------* IF BLOCK
*  |  First call initializations:
       IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
          TKESUM = ZERZER
          LFIRST = .FALSE.
          LUSSRC = .TRUE.
*  |  *** User initialization ***
          write(LUNOUT,*)
          write(LUNOUT,*) '--->modified source.f for Proton in solar min'
       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
* Kinetic energy of the particle (GeV)
* Sample the energy group
       XI = FLRNDM(DUMMY)
       DO K = 0 ,99
          IF (XI .GT. CUMENE(K) .AND. XI .LE. CUMENE(K+1)) THEN
             ENERGY = ((ENEUP(K) + ENELOW(K))/ TWOTWO)
             GO TO 600
          END IF
       END DO
  600   CONTINUE
       TKEFLK (NPFLKA) = ENERGY
*      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) ) )

*    +-----------------------------------------------------------------*
       Call SFLOOD ( XXX, YYY, ZZZ, UXXX, VXXX, WZZZ )
       UBEAM = UXXX
       VBEAM = VYYY
       WBEAM = WZZZ
* Cosines (tx,ty,tz)
       TXFLK  (NPFLKA) = UBEAM
       TYFLK  (NPFLKA) = VBEAM
       TZFLK  (NPFLKA) = WBEAM
*      TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
*&                        - TYFLK (NPFLKA)**2 )
* Polarization cosines:
       TXPOL  (NPFLKA) = -TWOTWO
       TYPOL  (NPFLKA) = +ZERZER
       TZPOL  (NPFLKA) = +ZERZER
* Particle coordinates ( radius of sphere =10 cm)
        XBEAM = XXX * 10.D0
        YBEAM = YYY * 10.D0
        ZBEAM = ZZZ * 10.D0
       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



