*$ 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.            *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      LOGICAL LFIRST
      LOGICAL bDebug
*
      SAVE LFIRST
      SAVE bDebug
      DATA LFIRST / .TRUE. /
      DATA bDebug / .FALSE. /
      DOUBLE PRECISION EKINSPECT,WSPECT
*======================================================================*
*                                                                      *
*                 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(*,*) '##################################################'
	 WRITE(*,*) '##'	 
         WRITE(*,*) '## Reading source spectrum from spectrum.data'
	 WRITE(*,*) '## V 1.0 28.06.2006 Helmut Vincke, Chris Theis'
	 WRITE(*,*) '##'	 	 
         WRITE(*,*) '##################################################'
         IF( WHASOU(4) .NE. 0.0 ) THEN
	   open(51,file='debug.log',status='UNKNOWN',
     &     form = 'formatted')	   
	   bDebug = .TRUE.
	 ELSE
	   bDebug = .FALSE.
	 END IF
	 
      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
*  Lstack is the stack counter: of course any time source is called it
*  must be =0
      NPFLKA = NPFLKA + 1
* KIETIC ENERGY AND WEIGHT OF THE PARTICLE 
      CALL SPECT(EKINSPECT, WSPECT)
      IF (WHASOU(1) .EQ. 0.0) THEN
        WHASOU(1) = 1.0
      ENDIF
      TKEFLK (NPFLKA) = EKINSPECT*WHASOU(1)
      WTFLK  (NPFLKA) = WSPECT
      WEIPRI = WEIPRI + WTFLK (NPFLKA)
* 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
      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) ,, ENERGY HANDLED ABOVE
*      TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
* Particle momentum
*      PMOFLK (NPFLKA) = PBEAM
      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
* 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 
* sample X & Y as a rectangular area of WHASOU(2) * WHASOU(3) with XBEAM & YBEAM
* being the center.
      IF( WHASOU(2) .LT. 0.0 .OR. WHASOU(3) .LT. 0.0 ) THEN
*        WRITE(*,*) 'ERROR: Source params 2 & 3 are invalid!!!'
*        WRITE(*,*) 'They should define the length of the rectangular beam spot!'	
        STOP       
      END IF

      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
*  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
      
      
      IF( bDebug .EQV. .TRUE.) THEN
        WRITE( 51,1000) IJBEAM, WEIPRI, XFLK(NPFLKA),YFLK(NPFLKA),
     &   ZFLK(NPFLKA), TXFLK(NPFLKA), TYFLK(NPFLKA), TZFLK(NPFLKA),
     &   TKEFLK(NPFLKA)
1000      format (1x, i3, 1x, E8.2, 7(1x, E12.5) )     
      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
      RETURN
*=== End of subroutine Source =========================================*
      END

      SUBROUTINE SPECT(EKINSPECT,WSPECT)
*
* provides ENERGY AND WEIGHT ACCORDING TO A SPECTRUM GIVEN IN FILE
*
      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
      INTEGER TOTMAX,MAX
      PARAMETER (TOTMAX=1000)
      DOUBLE PRECISION ESPECTRUM(TOTMAX),WSPECTRUM(TOTMAX),EKINSPECT,WSPECT
      DOUBLE PRECISION RNCURR,FLRNDM
      LOGICAL LOADED
      CHARACTER STERN*1
      SAVE MAX
      DATA MAX /0/
      SAVE LOADED
      DATA LOADED /.FALSE./
      SAVE ESPECTRUM,WSPECTRUM
      DATA ESPECTRUM/TOTMAX*0.0/
      DATA WSPECTRUM/TOTMAX*0.0/
*         
      IF (.NOT.LOADED) THEN
       open(50,file='../ambe12.data',status='OLD')
*         
C
C    IF * ==> COMMENT LINE
C     
       K=0
  150  CONTINUE 
       READ(50,'(A1)',END=200,ERR=1111) STERN
       IF (STERN .EQ. '*') THEN
        GOTO 150
       ENDIF
C
C    OTHERWISE PROGRAM READS THE LINE AS EXPECTED (2 FLOAT NUMBERS)
C
       BACKSPACE 50
       K = K+1
       READ(50,*,END=200,ERR=1111) Espectrum(K),Wspectrum(K)
       IF(Espectrum(K).EQ.0.0) THEN
         K=K-1
         GOTO 150
       ENDIF
       MAX = K
C
C    CHECK IF ALLOCATED SPACE FOR SPECTRUM IS BIG ENOUGH
C
       IF (MAX .GT. TOTMAX) THEN 
         WRITE(*,*)"TOTAL NUMBER OF SPECTRUM ENTRIES IS BIGGER THAN"
	 WRITE(*,*)"ALLOCATED SPACE. CHANGE TOTMAX IN SSPECT.f"
	 STOP
       ENDIF
       GOTO 150
 200   CONTINUE  
       LOADED = .TRUE.
      ENDIF
 250  CONTINUE  
C
C  GAMBLE ENERGY AND CHOOSE THE ACCORDING WEIGHT (=INTENSITY)
C
      RNCURR=FLRNDM(RNCURR)
      EKINSPECT=RNCURR*(ESPECTRUM(MAX)-ESPECTRUM(1))+ESPECTRUM(1)
      DO i=1,MAX
       IF (EKINSPECT .LT. ESPECTRUM(i)) THEN
C
C  CHECK IF WEIGHT VALUES ARE NEGATIVE OR EQUAL 0
C
         IF(WSPECTRUM(i-1) .EQ. 0.0D0) THEN
	   GOTO 250
         ELSEIF(WSPECTRUM(i-1) .LE. 0.0D0) THEN
	   WRITE(*,*) "NEGATIVE WEIGHT VALUES IN SPECTRUM"
	   WRITE(*,*) " RUN STOPS. PLEASE CHECK THE DATA"
	   STOP
	 ELSE
	   WSPECT = WSPECTRUM(i-1)
           GOTO 300
	 ENDIF
       ENDIF
      ENDDO
  300 CONTINUE       
      GOTO 1500   
 1111 CONTINUE
      CALL READERR()
 1500 CONTINUE
      END
      
C===============================================================================      
      SUBROUTINE READERR()
C
C FALLS PROGRAMM EINEN FEHLER WAEHREND DES READINGS FINDET DANN SPRINGT ES 
C HIERHER UM DIE UMGEBUNG DES FEHLERS AUF DEM SCHIRM AUSZUDRUCKEN. 
C
      CHARACTER CDUM*100
      BACKSPACE (50)
      READ(50,'(A100)') CDUM
      WRITE(*,'(A40)') ' PROGRAM HAS PROBLEMS WITH LINE BELOW'
      WRITE(*,*) 
      WRITE(*,'(A100)') CDUM
      WRITE(*,*) 
      WRITE(*,*) 'PROGRAM STOPS!!!!'
      STOP
      END
      
