*$ 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.            *
*                                                                      *
*----------------------------------------------------------------------*
*  Version for isotropic Heavy ion source
*
      DIMENSION  CUMENE(0:98), ENELOW(99), ENEUP(99) 
*----------------------------------------------------------------------*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      LOGICAL LFIRST
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
*    Heavy Ion energy group(Gev/n)
      DATA ENELOW /
     & 1.00000D-03, 1.14980D-03, 1.32190D-03, 1.51990D-03, 1.74750D-03,
     & 2.00920D-03, 2.31010D-03, 2.65610D-03, 3.05390D-03, 3.51120D-03,
     & 4.03700D-03, 4.64160D-03, 5.33670D-03, 6.13590D-03, 7.05480D-03,
     & 8.11130D-03, 9.32600D-03, 1.07230D-02, 1.23280D-02, 1.41750D-02,
     & 1.62980D-02, 1.87380D-02, 2.15440D-02, 2.47710D-02, 2.84800D-02,
     & 3.27450D-02, 3.76490D-02, 4.32880D-02, 4.97700D-02, 5.72240D-02,
     & 6.57930D-02, 7.56460D-02, 8.69750D-02, 1.00000D-01, 1.14980D-01,
     & 1.32190D-01, 1.51990D-01, 1.74750D-01, 2.00920D-01, 2.31010D-01,
     & 2.65610D-01, 3.05390D-01, 3.51120D-01, 4.03700D-01, 4.64160D-01,
     & 5.33670D-01, 6.13590D-01, 7.05480D-01, 8.11130D-01, 9.32600D-01,
     & 1.07230D+00, 1.23280D+00, 1.41750D+00, 1.62980D+00, 1.87380D+00,
     & 2.15440D+00, 2.47710D+00, 2.84800D+00, 3.27450D+00, 3.76490D+00,
     & 4.32880D+00, 4.97700D+00, 5.72240D+00, 6.57930D+00, 7.56460D+00,
     & 8.69750D+00, 1.00000D+01, 1.14980D+01, 1.32190D+01, 1.51990D+01,
     & 1.74750D+01, 2.00920D+01, 2.31010D+01, 2.65610D+01, 3.05390D+01,
     & 3.51120D+01, 4.03700D+01, 4.64160D+01, 5.33670D+01, 6.13590D+01,
     & 7.05480D+01, 8.11130D+01, 9.32600D+01, 1.07230D+02, 1.23280D+02,
     & 1.41750D+02, 1.62980D+02, 1.87380D+02, 2.15440D+02, 2.47710D+02,
     & 2.84800D+02, 3.27450D+02, 3.76490D+02, 4.32880D+02, 4.97700D+02,
     & 5.72240D+02, 6.57930D+02, 7.56460D+02, 8.69750D+02/
      DATA ENEUP /
     & 1.14980D-03, 1.32190D-03, 1.51990D-03, 1.74750D-03, 2.00920D-03,
     & 2.31010D-03, 2.65610D-03, 3.05390D-03, 3.51120D-03, 4.03700D-03,
     & 4.64160D-03, 5.33670D-03, 6.13590D-03, 7.05480D-03, 8.11130D-03,
     & 9.32600D-03, 1.07230D-02, 1.23280D-02, 1.41750D-02, 1.62980D-02,
     & 1.87380D-02, 2.15440D-02, 2.47710D-02, 2.84800D-02, 3.27450D-02,
     & 3.76490D-02, 4.32880D-02, 4.97700D-02, 5.72240D-02, 6.57930D-02,
     & 7.56460D-02, 8.69750D-02, 1.00000D-01, 1.14980D-01, 1.32190D-01,
     & 1.51990D-01, 1.74750D-01, 2.00920D-01, 2.31010D-01, 2.65610D-01,
     & 3.05390D-01, 3.51120D-01, 4.03700D-01, 4.64160D-01, 5.33670D-01,
     & 6.13590D-01, 7.05480D-01, 8.11130D-01, 9.32600D-01, 1.07230D+00,
     & 1.23280D+00, 1.41750D+00, 1.62980D+00, 1.87380D+00, 2.15440D+00,
     & 2.47710D+00, 2.84800D+00, 3.27450D+00, 3.76490D+00, 4.32880D+00,
     & 4.97700D+00, 5.72240D+00, 6.57930D+00, 7.56460D+00, 8.69750D+00,
     & 1.00000D+01, 1.14980D+01, 1.32190D+01, 1.51990D+01, 1.74750D+01,
     & 2.00920D+01, 2.31010D+01, 2.65610D+01, 3.05390D+01, 3.51120D+01,
     & 4.03700D+01, 4.64160D+01, 5.33670D+01, 6.13590D+01, 7.05480D+01,
     & 8.11130D+01, 9.32600D+01, 1.07230D+02, 1.23280D+02, 1.41750D+02,
     & 1.62980D+02, 1.87380D+02, 2.15440D+02, 2.47710D+02, 2.84800D+02,
     & 3.27450D+02, 3.76490D+02, 4.32880D+02, 4.97700D+02, 5.72240D+02,
     & 6.57930D+02, 7.56460D+02, 8.69750D+02, 1.00000D+03/
*    Normalized spectrum
      DATA CUMENE / 0.D0,
     & 3.11051E-06, 5.47369E-06, 8.61672E-06, 1.27980E-05, 1.83609E-05,
     & 2.57626E-05, 3.56111E-05, 4.87150E-05, 6.61565E-05, 8.93771E-05,
     & 1.20292E-04, 1.61460E-04, 2.16289E-04, 2.89316E-04, 3.86564E-04,
     & 5.16088E-04, 6.88345E-04, 9.17632E-04, 1.22215E-03, 1.62593E-03,
     & 2.16056E-03, 2.86662E-03, 3.79548E-03, 5.01318E-03, 6.60198E-03,
     & 8.66372E-03, 1.13218E-02, 1.47264E-02, 1.90525E-02, 2.45053E-02,
     & 3.13167E-02, 3.97438E-02, 5.00691E-02, 6.25734E-02, 7.75723E-02,
     & 9.53420E-02, 1.16151E-01, 1.40217E-01, 1.67702E-01, 1.98673E-01,
     & 2.33097E-01, 2.70837E-01, 3.11620E-01, 3.55035E-01, 4.00566E-01,
     & 4.47584E-01, 4.95385E-01, 5.43224E-01, 5.90360E-01, 6.36028E-01,
     & 6.79661E-01, 7.20665E-01, 7.58623E-01, 7.93271E-01, 8.24463E-01,
     & 8.52155E-01, 8.76442E-01, 8.97490E-01, 9.15532E-01, 9.30832E-01,
     & 9.43692E-01, 9.54402E-01, 9.63254E-01, 9.70518E-01, 9.76439E-01,
     & 9.81239E-01, 9.85107E-01, 9.88213E-01, 9.90694E-01, 9.92670E-01,
     & 9.94238E-01, 9.95479E-01, 9.96458E-01, 9.97229E-01, 9.97836E-01,
     & 9.98311E-01, 9.98684E-01, 9.98976E-01, 9.99203E-01, 9.99381E-01,
     & 9.99520E-01, 9.99628E-01, 9.99712E-01, 9.99778E-01, 9.99829E-01,
     & 9.99869E-01, 9.99899E-01, 9.99923E-01, 9.99942E-01, 9.99956E-01,
     & 9.99968E-01, 9.99976E-01, 9.99983E-01, 9.99988E-01, 9.99992E-01,
     & 9.99996E-01, 9.99998E-01, 1.00000E+00/
*  +-------------------------------------------------------------------*   
      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,*) '--->modified source.f for Silicon
     &                         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
*  +-------------------------------------------------------------------*
* Sample the energy group
      XI = FLRNDM(DUMMY)
      YI = FLRNDM(DUMMY1)
      DO 500 K = 1 ,99
         IF (XI .LE. CUMENE(K)) THEN
            ENERGY = (YI * (ENEUP(K) - ENELOW(K)) + ENELOW(K)) * 28.D0
            GO TO 600
         END IF
 500   CONTINUE
       STOP ' Failed to sample the energy group '
 600   CONTINUE
       count=count+1
      NPFLKA = NPFLKA + 1
      WTFLK (NPFLKA) = ONEONE
      WEIPRI = WEIPRI + WTFLK (NPFLKA)
*  +--------------------------------------------------------------------*
*  |  (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)
      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) ) )

*    +-----------------------------------------------------------------*
      RANPHI = FLRNDM(U)
      PHID = TWOPIP * RANPHI
      RANCOS = FLRNDM(V)
      COSD = Oneone - Twotwo * RANCOS
      SIND = SQRT(Oneone - COSD * COSD)
      UBEAM = -SIND * COS(PHID)
      VBEAM = -SIND * SIN(PHID)
      WBEAM = -COSD
* 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
*    +-----------------------------------------------------------------*
      RAD = 285
      ZS = 500
      RANPHI = FLRNDM(P)
      PHIS = TWOPIP * RANPHI
      RANR = FLRNDM(R)
      RS = RAD * SQRT(RANR)
      XS = RS * COS(PHIS)
      YS = RS * SIN(PHIS)
      XBEAM = XS * COSD * COS(PHID) - YS * SIN(PHID) + ZS * SIND
     &                                   * COS(PHID)
      YBEAM = XS * COSD * SIN(PHID) + YS * COS(PHID) + ZS * SIND
     &                                   * SIN(PHID)
      ZBEAM = -XS * SIND + ZS * COSD
* Particle coordinates
      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


