*$ 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.            *
*                                                                      *
*       Output variables:                                              *
*                                                                      *
*              Nomore = if > 0 the run will be terminated              *
*                                                                      *
*  Read source term from a file                                        *
*     Author: Vasilis.Vlachoudis@cern.ch                               *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      LOGICAL LFIRST

      PARAMETER (NMAX=1000000)
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /

      CHARACTER*250 LINE
      CHARACTER*8 SRCDUM

*     Particle counter, particle ID
      INTEGER    NNN, IJ(NMAX), PTYPE
*     Coordinates
      DIMENSION  XXX(NMAX), YYY(NMAX), ZZZ(NMAX)
*     Directional cosines
      DIMENSION  UUU(NMAX), VVV(NMAX), WWW(NMAX)
*     Energy, weight
      DIMENSION  ERG(NMAX), WGT(NMAX)
      
      SAVE XXX, YYY, ZZZ
      SAVE UUU, VVV, WWW
      SAVE IJ, ERG, WGT
*      mypath="/home/cliu/build/fluka-2011/test/cosmic/second-step/"
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.

*  |  *** User initialization ***
*  |  WHASOU(1) = UNIT number (File name is linked with OPEN)
         LUNRD = NINT(WHASOU(1))
         PTYPE = NINT(WHASOU(2))

         SRCDUM = SDUSOU
*         IF (SRCDUM .EQ. 'PLANE') THEN
*             CALL OAUXFI('tempfile', LUNRD, 'OLD', IERR)
*         ELSE IF (SRCDUM .EQ. 'CIRCULAR') THEN
*             CALL OAUXFI('tempcirc', LUNRD, 'OLD', IERR)
*         ELSE
*             CALL OAUXFI('tempfile', LUNRD, 'OLD', IERR)
*         END IF
*         open(unit=81,file='tempfile', form='FORMATTED',status='OLD')
         WTOT = ZERZER
         NNN = 0
 10      CONTINUE
            READ( LUNRD, '(A)', ERR=9999, END=20 ) LINE
*            WRITE(*,*) LINE

* I: particle ID | X,Y,Z (x,y,z) | U,V,W (u,v,w) | E: energy | WG: weight
            READ (LINE,*,ERR=10) I, X, Y, Z, U, V, W, E, WG
            IF (PTYPE .NE. I .AND. PTYPE .NE. 201) THEN
                WRITE(*,*) 'Wrong type of particle:',I,'instead of',
     & PTYPE
                GOTO 10
            END IF
*            WRITE(*,*) PTYPE, I
            IF( SRCDUM .EQ. 'SFLOOD') THEN
                CALL SFLOOD(X,Y,Z,U,V,W)
                IF (Z .LT. 0) THEN
                    Z = Z*(-1)
                    Y = Y*(-1)
                    X = X*(-1)
                    U = MODULO(PI - U,2*PI)
                    V = MODULO(PI - V,2*PI)
                    W = MODULO(PI - W,2*PI)
                END IF
                X = X*10
                Y = Y*10
                Z = Z*10
            END IF

*  |        Rejection criteria:
*  |        store only points of interest
*            IF (PTYPE .EQ. I .OR. PTYPE .EQ. 201) THEN
*            WRITE(*,*) "PASSED"
            NNN = NNN + 1
            IF (NNN.GT.NMAX) CALL FLABRT('SOURCE','Increase NMAX')
            WTOT = WTOT + WG

               IJ(NNN)  = I

               XXX(NNN) = X
               YYY(NNN) = Y
               ZZZ(NNN) = Z

*  |  Normalize direction to 1.0
               UVW = SQRT(U**2 + V**2 + W**2)
               UUU(NNN) = U / UVW
               VVV(NNN) = V / UVW
               WWW(NNN) = W / UVW

               ERG(NNN) = E
               WGT(NNN) = WG
*            ENDIF
         GOTO 10
 20      CONTINUE
         WRITE(*,*) NNN
         IF (NNN.EQ.0) CALL FLABRT('SOURCE','Error reading file')
         WRITE (LUNOUT,*)
         WRITE (LUNOUT,*) '*** rdsource: ',NNN,' particles loaded',
     &' - total weight: ',WTOT
         WRITE (LUNOUT,*)
      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
*  +-------------------------------------------------------------------*
*  | Choose a random particle
*      RNDSIG = FLRNDM (RNDSIG)
*      N = INT(NNN*RNDSIG)+1
*****      NN=0
310   continue
      N=N+1
      WRITE(*,*) "TESTAA", PTYPE, IJ(N), N, NPFLKA
**      IF (PTYPE .NE. IJ(N) .AND. PTYPE .NE. 201) THEN 
**      GOTO 310
**      ENDIF      
*****      IF (N .GT. NNN) THEN
*****          write(*,'(A)') 'exceeded array size, leaving subroutine...'
*****          GOTO 777
*****      END IF
      write(*,'(A)') 'New particle loaded'
*****      IF (PTYPE .EQ. IJ(N) .OR. PTYPE .EQ. 201) THEN
*****          NN = NN+1
          WRITE (LUNOUT,*) '*** N=',N, " ",IJ(N)," ",ERG(N)," ",WGT(N)
*  |
*  +-------------------------------------------------------------------*
*  Wt is the weight of the particle
          WTFLK  (NPFLKA) = WGT(N)
          WEIPRI = WEIPRI + WTFLK (NPFLKA)
*  Particle type (1=proton.....). Ijbeam is the type set by the BEAM
*  card
*  +-------------------------------------------------------------------*
*  |  (Radioactive) isotope:
          write(*,'(A)') 'Entering case by case'
          IF ( IJ(N) .EQ. -2 .AND. LRDBEA) THEN
          write(*,'(A)') 'FLAG1'
             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 ( IJ(N) .EQ. -2) THEN
          write(*,'(A)') 'FLAG2'
             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 
          write(*,'(A)') 'FLAG3'
             IONID = IJ(N)
             ILOFLK (NPFLKA) = IJ(N)
*  |  Flag this is prompt radiation
             LRADDC (NPFLKA) = .FALSE.
*  |  Group number for "low" energy neutrons, set to 0 anyway
             IGROUP (NPFLKA) = 0
          END IF
*  |
*  +-------------------------------------------------------------------*
          write(*,'(A)') 'FLAG4'
*  From this point .....
*  Particle generation (1 for primaries)
          LOFLK  (NPFLKA) = 1
*  User dependent flag:
          LOUSE  (NPFLKA) = 0
*  No channeling:
**      KCHFLK (NPFLKA) = 0
**      ECRFLK (NPFLKA) = ZERZER
*  Extra infos:
**      INFSTK (NPFLKA) = 0
**      LNFSTK (NPFLKA) = 0
**      ANFSTK (NPFLKA) = ZERZER
*  Parent variables:
**      IPRSTK (NPFLKA) = 0
**      EKPSTK (NPFLKA) = ZERZER
*  User dependent spare variables:
          write(*,'(A)') 'FLAG5'
          DO 100 ISPR = 1, MKBMX1
             SPAREK (ISPR,NPFLKA) = ZERZER
 100  CONTINUE
*  User dependent spare flags:
          write(*,'(A)') 'FLAG6'
          DO 200 ISPR = 1, MKBMX2
             ISPARK (ISPR,NPFLKA) = 0
 200  CONTINUE
*  Save the track number of the stack particle:
      write(*,'(A)') 'FLAG7'
          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)
          write(*,'(A)') 'FLAG8'
          TKEFLK (NPFLKA) = ERG(N)
*  Particle momentum
*      PMOFLK (NPFLKA) = PBEAM
          PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (IONID) ) )
*  Cosines (tx,ty,tz)
          TXFLK  (NPFLKA) = UUU(N)
          TYFLK  (NPFLKA) = VVV(N)
          TZFLK  (NPFLKA) = WWW(N)
*  Polarization cosines:
          TXPOL  (NPFLKA) = -TWOTWO
          TYPOL  (NPFLKA) = +ZERZER
          TZPOL  (NPFLKA) = +ZERZER
*  Particle coordinates
          XFLK   (NPFLKA) = XXX(N) !+ XBEAM
          YFLK   (NPFLKA) = YYY(N) !+ YBEAM
          ZFLK   (NPFLKA) = ZZZ(N) !+ ZBEAM

*  Calculate the total kinetic energy of the primaries: don't change
          write(*,'(A)') 'FLAG9'
          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:
          write(*,'(A)') 'FLAG10'
          CALL GEOCRS ( TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA) )
          CALL GEOREG ( XFLK  (NPFLKA), YFLK  (NPFLKA), ZFLK  (NPFLKA),
     &              NRGFLK(NPFLKA), IDISC )
*  Do not change these cards:
          write(*,'(A)') 'FLAG11'
          CALL GEOHSM ( NHSPNT (NPFLKA), 1, -11, MLATTC )
          NLATTC (NPFLKA) = MLATTC
          CMPATH (NPFLKA) = ZERZER
          CALL SOEVSV
*****      ELSE
*****          write(*,'(A)') 'wrong particle type'
*****          GOTO 310
*****      END IF
*****777   CONTINUE
      write(*,'(A)') 'end of subroutine'
      write(*,'(A)') '-------------------------------------------------'
      RETURN
*=== End of subroutine Source =========================================*
 9999 CONTINUE
      CALL FLABRT('SOURCE','Error reading source file')
      END
