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

       INCLUDE '(DBLPRC)'
       INCLUDE '(DIMPAR)'
       INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2010      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  17-Oct-10    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)'
*
       LOGICAL LFIRST
*     final computed momentum (jba 10.15.12)
*     Is this necessary?
       DOUBLE PRECISION PMUFIN
       DOUBLE PRECISION PMUTST


*
       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
       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
*  ( This part here is uninteresting (jba))
*  +-------------------------------------------------------------------*
*  |  (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 .....

*  But we want to generate muons randomly, so instead...

* BTW, learn more about how source.f works here:
* http://fluka-course.web.cern.ch/fluka-course/
*        Lectures_pdf/c14_UserRoutines.pdf

* We don't care about absolute normalization

* Apologies to the coding gods for using a goto here.  Blame Fortran.

*  Here we begin making the angular distribution throw
*  vdepth = vertical depth (m.w.e.)
*  angtar = target for prob distribution (gets almost up to 1)
*  angcos = cosine of vertical angle, 1 = straight down
*  angpro = prob of getting that angle

       vdepth = 1585.0D+00

  97   continue
          write(*,*) "We've hit 97!"
          angtar = FLRNDM(XDUMMY)
          angcos = FLRNDM(XDUMMY)
          angpro = angcos**1.53D0 * EXP(-ONEONE * vdepth * 0.0008D0 *
      &                               (ONEONE/angcos - ONEONE))
          if (angpro .lt. angtar) then
             goto 97
          endif
*     Now we should have our vertical angle!

*     Now, compute the other angle direction with built-in function
*     We will have SINT become the Y-cosine, and COST become the X-cos
*     WBEAMU = z-cosine in our user routine
*     UBEAMU = x-cosine in our user routine
*     VBEAMU = y-cosine in our user routine
*
       CALL SFECFE(SINT,COST)
       WBEAMU = -ONEONE * angcos
       hozcos = SQRT( ONEONE - WBEAMU * WBEAMU )
*     Normalize the x and y cosines, carefully
       UBEAMU = COST * hozcos
       VBEAMU = SINT * hozcos
*     We should now have normalized angles!  For debugging....
*      WRITE(*,*) WBEAMU**2 + UBEAMU**2 + VBEAMU**2

*     Next, compute the slant depth for this particular muon...
*     For this formula, we will use km.w.e.
       sldpth = vdepth / angcos / 1000.0D0
*     Now, do the same MC thing, but now for energy
*     We will try this with energies from 0 to 2 TeV
*     We have some parameters from Mei + Hime
       embvar = 0.4D+00
       emgm   = 3.77D+00
       emepm  = 693.0D+00

98    Continue
          write(*,*) "We've hit 98!"
          emutry = FLRNDM(XDUMMY) * 2000.0D+00
          emutar = FLRNDM(XDUMMY)
          emnorm = EXP(-embvar * sldpth * (emgm - ONEONE))
      &             * (emepm *(ONEONE - EXP(-embvar*sldpth)))**(-emgm)
          emupro = EXP(-embvar * sldpth * (emgm - ONEONE)) / emnorm
      &    * (emutry + emepm *(ONEONE-EXP(-embvar*sldpth)))**(-emgm)
          if (emupro .lt. emutar) then
             goto 98
          endif
*     Now we should have the muon energy! Let's get a momentum.
       PMUFIN = SQRT(dble(emutry**TWOTWO) - dble(AM(IONID)**TWOTWO))
       PMUTST = PMUFIN
*      PMUFIN = 252.232229D+00
*      PMUFIN = 0.252232228711425506162414D+03

*     Next up, the muon position.  We will give it a random position
*     2.5 meters above the detector center, and we'll spread the
*     origin 3m each way in X and Y, so it will be in total a
*     6x6 m square, for 36 m^2 total.  This is important for
*     when we normalize to the actual muon flux
       XPOS = FLRNDM(XDUMMY) * 600.0D0 - 300.0D0
       YPOS = FLRNDM(XDUMMY) * 600.0D0 - 300.0D0
       ZPOS = 250.0D0

*  Particle generation (1 for primaries)
       LOFLK  (NPFLKA) = 1
*  User dependent flag:
       LOUSE  (NPFLKA) = 0
*  No channeling:
       LCHFLK (NPFLKA) = .FALSE.
       DCHFLK (NPFLKA) = ZERZER
*  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) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
       TKEFLK (NPFLKA) = SQRT ( PMUFIN**2 + AM (IONID)**2 ) - AM (IONID)
*  Particle momentum
*      PMOFLK (NPFLKA) = PBEAM
*      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
*     &                       + TWOTWO * AM (IONID) ) )
       PMOFLK (NPFLKA) = PMUFIN
*  Cosines (tx,ty,tz)
       TXFLK  (NPFLKA) = UBEAMU
       TYFLK  (NPFLKA) = VBEAMU
       TZFLK  (NPFLKA) = WBEAMU
*     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) = XBEAM
*      YFLK   (NPFLKA) = YBEAM
*      ZFLK   (NPFLKA) = ZBEAM
       XFLK   (NPFLKA) = XPOS
       YFLK   (NPFLKA) = YPOS
       ZFLK   (NPFLKA) = ZPOS
*  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

*  Unfortunate printout!
       write(*,*) "Just about done the source.f thing! Our params are..."
       write(*,*) XFLK(NPFLKA), YFLK(NPFLKA), ZFLK(NPFLKA), NPFLKA
       write(*,*) TXFLK(NPFLKA), TYFLK(NPFLKA), TZFLK(NPFLKA)
       write(*,40) PMUFIN, PMOFLK(NPFLKA), TKEFLK(NPFLKA)
*      write(*,40) PMUTST, PMOFLK(NPFLKA), TKEFLK(NPFLKA)
40    format (3d30.24)

*  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

