*$ 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
INTEGER NNN, IJ(NMAX)
DIMENSION XXX(NMAX), YYY(NMAX), ZZZ(NMAX)
DIMENSION UUU(NMAX), VVV(NMAX), WWW(NMAX)
DIMENSION ERG(NMAX), WGT(NMAX)
SAVE XXX, YYY, ZZZ
SAVE UUU, VVV, WWW
SAVE IJ, ERG, WGT
*======================================================================*
* *
* 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))
POSMIN = WHASOU(2)
POSMAX = WHASOU(3)
WTOT = ZERZER
NNN = 0
10 CONTINUE
READ( LUNRD, '(A)', ERR=9999, END=20 ) LINE
READ (LINE,*,ERR=10) I, X, Y, Z, U, V, W, E, WG
* | rejection criteria:
* | store only points of interest
IF ( POSMIN.LT.X .AND. X.LT.POSMAX ) THEN
NNN = NNN + 1
WTOT = WTOT + WG
IF (NNN.GT.NMAX) CALL FLABRT('SOURCE','Increase NMAX')
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
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
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:
IF ( IJ(N) .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 ( IJ(N) .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 = 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
* |
* +-------------------------------------------------------------------*
* 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:
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) = 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
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 =========================================*
9999 CONTINUE
CALL FLABRT('SOURCE','Error reading source file')
END