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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2005      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 17-jun-05     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.            *
*                                                                      *
*                                                                      *
* Modified to load particles from a particle DB                        *
*                                                                      *
* Version 1.0	28.11.2005	       C. Theis                        *
* 
*                          
* Version 1.1   19.04.2011             L. Jaegerhofer                  *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
      INCLUDE '(CASLIM)'
*
      LOGICAL LFIRST
      CHARACTER*(16) FileName      
      LOGICAL bDebug
      LOGICAL LDIREC
      integer COUNTER, PARTNUM,NCASESOLD,HELITOT, nrpart, i
      common /grslbl/lblout
      common /helicom/WEITOTAL,HELITOT
      integer nngrs
      save COUNTER,PARTNUM,NCASESOLD,nngrs
    
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /

      double precision x, y , z, weight,cox, coy, coz, ekin
      double precision WEITOTAL, TX, TY, TZ
      integer ij
      SAVE TX, TY, TZ
      
       
*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialization ***
         write(90,*) 'ILOFLK (NPFLKA), XFLK(NPFLKA), YFLK(NPFLKA), 
     &   ZFLK(NPFLKA), TXFLK(NPFLKA), TYFLK(NPFLKA), TZFLK(NPFLKA),
     &   TKEFLK(NPFLKA)'
         write(90,*) 'ParticleID, PosX, PosY, PosZ, 
     &   DirX, DirX, DirX, EKin'
 
         bDebug = .FALSE.
* The first parameter indicates the number of iterations
* over the particle database
* 
* open debug log if WHAT(6) is set!	 
*
   	 IF( WHASOU(6) .NE. 0.0 ) THEN 
	   OPEN(98,FILE="debug.log", STATUS='UNKNOWN')          
           bDebug = .TRUE.
         END IF

	 IF( SDUSOU.EQ.'RUN2' ) THEN
  
           ctr=0
	   nngrs=0
           counter=0	 


* assemble filename & path        

           FileName = 'particles.dat'	      
* Open file & read number of primaries (== WEIPRI of 1st run)
	   WRITE(LUNOUT,*) 'Trying to open particle database...'
           OPEN(99,FILE=FileName,STATUS='OLD')
	   WRITE(LUNOUT,*) 'Opened file ', FileName
           READ(99,*,end=3000) NCASESOLD
             
           WRITE(LUNOUT,*) 
     &     '-------------------------------------------------------'
           WRITE(LUNOUT,*) 
     &     'Precalculated particle source implementation by Chris Theis'
           WRITE(LUNOUT,*) 
     &     'modified by Lukas Jaegerhofer'
           WRITE(LUNOUT,*) 
     &     'Version 2.1 12.04.2011'
           WRITE(LUNOUT,*) 
     &     '-------------------------------------------------------'
           WRITE(LUNOUT,*) 
     &     '| Reading from file ', FileName, ' Nr of primaries: ', 
     &     NCASESOLD
           WRITE(LUNOUT,*) 
     &     '-------------------------------------------------------'


	 END IF

* print some debug statements
         write (LUNOUT,*) 'w(1) w(2) w(3) w(4) w(5) w(6)'
         write (LUNOUT,*) whasou(1), whasou(2), whasou(3), whasou(4),
     &               whasou(5), whasou(6) 
         write (LUNOUT,*) 'w(7) w(8) w(9) w(10) w(11) w(12)'
         write (LUNOUT,*) whasou(7), whasou(8), whasou(9), whasou(10),
     &               whasou(11), whasou(12) 
         WRITE(LUNOUT,*) 
     &     '-------------------------------------------------------'

      END IF


*  |
*  +-------------------------------------------------------------------*
*  +-- READ PARTICLE DATA *
      IF( SDUSOU.EQ.'RUN2' ) THEN

1000      READ( 99, 1210, END = 3000, ERR = 1100 )
     &   ij , weight, x, y, z, cox, coy, coz, ekin

       


*
* Check if we are in the second run
* and if a special particle type was
* requested with WHAT(2)
*
          IF( SDUSOU.EQ.'RUN2' .AND. WHASOU(2) .NE. 0) THEN
            IF( WHASOU(2) .NE. IJ ) THEN
	      GOTO 1000
  	    END IF
          ENDIF

*
* debug = 1 ... write statement only during the 1st iteration
*
          IF( WHASOU(6) .EQ. 1.0 .AND. COUNTER .EQ. 0 ) THEN
            WRITE(98,1210) ij , weight, x, y, z, cox, coy, coz, ekin
	  END IF
*
* debug = 2 ... write statement for every reading
*  
          IF( WHASOU(6) .EQ. 2.0 ) THEN
            WRITE(98,1210) ij , weight, x, y, z, cox, coy, coz, ekin
	  END IF
	  
* Number of lines per particle database
         nngrs=nngrs+1

*       IJBEAM = ij

	 
         if (ekin .LE. 0) goto 1100
* set the particle type

   
         GOTO 1200
1100     WRITE(*,*) '+----------------------------------------+'
         WRITE(*,*) '| Error while reading from file!         |'
         WRITE(*,*) '+----------------------------------------+'
         WRITE(*, fmt = 1210)
     &   ij , weight, x, y, z, cox, coy, coz, ekin
*  if the end of the particle database is reached or an error occurred
*  we stop
         STOP 'Error (Ekin <= 0). The system rolled over & died!'
*
1200     END IF
1210     FORMAT (1x, i3, 1x, E8.2, 7(1x, E12.5) )

*  +-------------------------------------------------------------------*
*load 2 particles at the stack - later it will be the array of particles 
*scored in particles.dat originating from the same primary particle.

      nrpart = 2
      
      do 10 i = 1, nrpart

      NPFLKA = NPFLKA + 1


* set the particle type
         IF ( i.EQ.1 ) THEN
	 IJBEAM = 1
        WRITE(LUNOUT,*) 'i: ', i, ' PartID: ', IJBEAM
* set random directionL
         END IF
         IF ( i.EQ.1 ) THEN
         call RACO (TXX, TYY, TZZ)
	 TX = TXX
	 TY = TYY
	 TZ = TZZ
	 END IF
*         ELSE 
	 if ( i.EQ.2 ) THEN
	 IJBEAM = -4
         WRITE(LUNOUT,*) 'i: ', i, ' PartID: ', IJBEAM
        ENDIF
*      IJBEAM = -4	 
         

* Wt is the weight of the particle
*
* If WHAT(1) != 0 then the particle weight is taken
* from what was read from the file. Furthermore 
      IF( SDUSOU.EQ.'RUN2' ) THEN
        WTFLK( NPFLKA) = weight
        WEIPRI=NCASESOLD*WHASOU(1)
      ELSE 
        WTFLK  (NPFLKA) = ONEONE
	WEIPRI = WEIPRI + WTFLK (NPFLKA)
      END IF
            
* Particle type (1=proton.....). Ijbeam is the type set by the BEAM
* card
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      IF ( IJBEAM .EQ. -2 ) THEN
         IJHION = IPROZ  * 1000 + IPROA
         IJHION = IJHION * 100 + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         ILOFLK (NPFLKA) = IJHION
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID  = IJBEAM
         ILOFLK (NPFLKA) = IJBEAM
         WRITE(LUNOUT,*) 'ionID: ', IONID, ' ILOFLK: ', ILOFLK (NPFLKA)
      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


      IF( SDUSOU.EQ.'RUN2' ) THEN
* Read the kinetic energy directly
        IF ( i.EQ.1 ) THEN
	TKEFLK(NPFLKA) = 0.000573
        PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
        ENDIF
        IF ( i.EQ.2 ) THEN
	TKEFLK(NPFLKA) = 0.000191
        PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
        ENDIF	

      ELSE
* Standard behavior - calculate kinetic energy from the beam card
        TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
* Particle momentum
        PMOFLK (NPFLKA) = PBEAM
      END IF


* Cosines (tx,ty,tz)
*
* ---> DON'T forget that the cosine MUST be normalized!!!
* The normalization via the division could lead to trouble on 64bit machines
* under some circumstances, but for the time being 
*       
      IF( SDUSOU.EQ.'RUN2' ) THEN

*        COSSUM = SQRT( COX * COX + COY * COY + COZ * COZ )
*        TXFLK  (NPFLKA) = COX / COSSUM
*        TYFLK  (NPFLKA) = COY / COSSUM
*        TZFLK  (NPFLKA) = COZ / COSSUM
        IF ( i.EQ.1 ) THEN
	 TXFLK  (NPFLKA) = TX
         TYFLK  (NPFLKA) = TY
         TZFLK  (NPFLKA) = TZ
	ENDIF
        IF ( i.EQ.2 ) THEN
         TXFLK  (NPFLKA) = -TX
         TYFLK  (NPFLKA) = -TY
         TZFLK  (NPFLKA) = -TZ
	ENDIF
      ELSE 
* Take particle direction from BEAM card although there is a SOURCE card!
       LDIREC = .FALSE.
       CALL STBMDV( LDIREC ) 
      END IF
      
*     TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
*    &                       - TYFLK (NPFLKA)**2 )
* Polarization cosines:
* Take particle polarization from BEAM card although there is a SOURCE card!
      CALL STBMPL
* Particle coordinates
      IF( SDUSOU.EQ.'RUN2' )  THEN
        IF ( i.EQ.1 ) THEN
        XFLK   (NPFLKA) = x
        YFLK   (NPFLKA) = y
        ZFLK   (NPFLKA) = z
	ENDIF
        IF ( i.EQ.2 ) THEN
        XFLK   (NPFLKA) = x
        YFLK   (NPFLKA) = y
        ZFLK   (NPFLKA) = z
        ENDIF
*for debugging writing in fort.90 file:
        write(90,*) ILOFLK (NPFLKA), XFLK(NPFLKA), YFLK(NPFLKA), 
     &  ZFLK(NPFLKA), TXFLK(NPFLKA), TYFLK(NPFLKA), TZFLK(NPFLKA),
     &  TKEFLK(NPFLKA)

* Sampling the beam in a rectangular field with lengths What(2) x What(3) 
* around BEAM coordinates
      ELSE 
       IF( WHASOU(2) .EQ. ZERZER .AND. WHASOU(3) .EQ. ZERZER ) THEN 	
       XFLK   (NPFLKA) = XBEAM
       YFLK   (NPFLKA) = YBEAM
       ZFLK   (NPFLKA) = ZBEAM
 

* Take particle polarization from BEAM card although there is a SOURCE card!
* Sampling the beam in a rectangular field with lengths What(2) x What(3) 
* around BEAM coordinates
 
	ELSE IF( WHASOU(2) .GT. ZERZER .AND. WHASOU(3) .GT. ZERZER) THEN
       RAN1 = FLRNDM(RAN1)
       RAN2 = FLRNDM(RAN2)      
       XFLK   (NPFLKA) = XBEAM + (WHASOU(2) * RAN1) - (WHASOU(2) / 2.0)
       YFLK   (NPFLKA) = YBEAM + (WHASOU(3) * RAN2) - (WHASOU(3) / 2.0)
       ZFLK   (NPFLKA) = ZBEAM
*        write(90,*) XFLK(NPFLKA), YFLK(NPFLKA), ZFLK(NPFLKA), 
*     &  DX, DY, DZ
*50      FORMAT("8 1 ", 6F12.6," 1")
            ELSE
             WRITE(*,*) 'ERROR: Source params 2 & 3 are invalid!!!'
             WRITE(*,*) 'They should define the length of the 
     &         rectangular beam spot, '	
             WRITE(*,*) 'or they must both equal 0 for a point source! '	
            STOP       
        END IF
       END IF

*  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
             

*  Flag this is prompt radiation
      LRADDC (NPFLKA) = .FALSE.
      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

*      
* debug = 2 - write current particle parameters
*
      if ( WHASOU(6) .GE. 2.0 ) THEN
        WRITE(98,2122) ILOFLK (NPFLKA),WTFLK (NPFLKA),
     &   XFLK (NPFLKA), YFLK (NPFLKA), ZFLK (NPFLKA),
     &   TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK(NPFLKA), 
     &   ekin, TKEFLK(NPFLKA), PMOFLK (NPFLKA), AM (IONID)
        WRITE(98,*) '----------------------------------------',NCASE

      endif   
2122  format (1x, i3, 1x, E8.2, 9(1x, E12.5) )

10    continue
2500  RETURN

* this part decides what to do in case we reach the end of the primary file.
* depending of the number given in what(1) in card source the repetition of 
* reading and using the primary particle file is performed. Also the total 
* number and the total weight to use for the normalisation is calculated here.
*       
3000  continue

      IF( SDUSOU.EQ.'RUN2' ) THEN 
* Sanity check. We have at least ONE iteration over the file
        IF(WHASOU(1).EQ.0.0) WHASOU(1)=1.0
        counter=counter+1

* After the first run set the internal counter of the total number of
* particles that are to be processed.
* This can only be done after the first iteration over the particle DB
* because we have counted the number of entries.

        IF (COUNTER.EQ.1)THEN
         partnum=nngrs
         HELITOT=nint(WHASOU(1)*partnum)
	 WRITE(LUNOUT,*) 'NCASES =', ncases
         ncases=HELITOT
	 WRITE(LUNOUT,*) 'NEW NCASES =', ncases	 	 

        END IF

*
* Perform the normalization to NrOfIterations over the particle file times the
* weight of the primaries of the original 1st run.
*
        IF( NCASES.EQ.NCASE) THEN
          WEIPRI=NCASESOLD*WHASOU(1)
* we have reached the max. nr of particles --> let's bail out	

          NOMORE = 1

        END IF
	
	IF( NCASE .LT. NCASES ) THEN
          REWIND(99)
          READ(99,*,END=3000) NCASESOLD

          GO TO 1000
	END IF

*
* If we only have one iteration over the particle DB then NCASES might report
* one particle too much, although the weight of the primaries is okay and also
* the number of particles that are processed are correct. In order to avoid
* confusion in the output files where the use would read e.g. 11 particles 
* processed with a weight of 10, we manually fix this.
*	
	IF( NCASE .GT. NCASES .AND. COUNTER .EQ. 1) THEN
	  NCASE = NCASES

	END IF	
	
      END IF    
*=== End of subroutine Source =========================================*
      END
