*$ 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 * * * *----------------------------------------------------------------------* * * * DIMENSION CUMPR(43), ENEDGE(43) *----------------------------------------------------------------------* * INCLUDE '(BEAMCM)' INCLUDE '(CASLIM)' INCLUDE '(FHEAVY)' INCLUDE '(FLKSTK)' INCLUDE '(IOIOCM)' INCLUDE '(LTCLCM)' INCLUDE '(PAPROP)' INCLUDE '(SOURCM)' INCLUDE '(SUMCOU)' * LOGICAL LFIRST * SAVE LFIRST DATA LFIRST / .TRUE. / *----------------------------------------------------------------------* * * Proton energy group boundaries !43 energies DATA ENEDGE / & 0.11521, & 0.11655, & 0.11787, & 0.11918, & 0.12045, & 0.12171, & 0.12296, & 0.12413, & 0.12539, & 0.12666, & 0.12788, & 0.12909, & 0.13028, & 0.13139, & 0.13261, & 0.13383, & 0.13501, & 0.13618, & 0.13733, & 0.13845, & 0.13955, & 0.14073, & 0.14189, & 0.14300, & 0.14409, & 0.14524, & 0.14635, & 0.14745, & 0.14855, & 0.14962, & 0.15066, & 0.15177, & 0.15287, & 0.15386, & 0.15497, & 0.15710, & 0.15920, & 0.16127, & 0.16332, & 0.16534, & 0.16738, & 0.16938, & 0.17137/ * * Cumulative spectrum !43 weights DATA CUMPR / & 0.01054430535851470, & 0.01431470005844220, & 0.02235293146569520, & 0.02786211997428630, & 0.03479711441176390, & 0.04115148724087290, & 0.04810046058817790, & 0.05417059411222460, & 0.06234335477433980, & 0.06898244954716260, & 0.07635045114455820, & 0.08403102422993940, & 0.09107531571584640, & 0.09881893629872820, & 0.10689405545350300, & 0.11551085231704100, & 0.12307495747963500, & 0.13255567582233000, & 0.14031266983304600, & 0.14992378276338300, & 0.15826892365338800, & 0.16868329100665900, & 0.17827383468625500, & 0.18867374702026500, & 0.19804611723494700, & 0.21038811454095700, & 0.21948664581049600, & 0.23350115883350400, & 0.24206807468343800, & 0.25798986891144900, & 0.26550517262357100, & 0.28579903901021900, & 0.29170908517045800, & 0.31179289083864600, & 0.33021480199225000, & 0.36701408972104000, & 0.40310877546159700, & 0.44540179368478600, & 0.49381078238080200, & 0.55065625008642900, & 0.62417519316010700, & 0.71312514880381800, & 1/ * *======================================================================* * * * 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 * | * +-------------------------------------------------------------------* * * In order to cover the target volume in depth the energy must be * sampled from the ENEDGE, with cumulative weights from CUMPR * Sample the energy group XI = FLRNDM(DUMMY) !FLUKA random number generator, 1 is not included DO 500 K = 1, 43 !500 statement label IF(XI .LE. CUMPR(K)) THEN !.LE. meaning <= ENERGY = ENEDGE(K) GO TO 600 END IF 500 CONTINUE STOP 'Failed to sample the energy group' 600 CONTINUE * * 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 * +-------------------------------------------------------------------* * | (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 ..... * 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) * With Gaussian sigma spread 0.1% initial energy CALL FLNRRN(RGAUSS) ! which returns a normal gaussian number RGAUSS TKEFLK (NPFLKA) = ENERGY * (1 + 0.004246609001440 * RGAUSS) !WRITE(101,*)RGAUSS !TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID) * Particle momentum !PMOFLK (NPFLKA) = PBEAM PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA) & + TWOTWO * AM (IONID) ) ) * 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 XFLK (NPFLKA) = (FLRNDM(DUMMY)-0.5)*2.0*(10) YFLK (NPFLKA) = (FLRNDM(DUMMY)-0.5)*2.0*(10) ZFLK (NPFLKA) = -0.1 !XFLK (NPFLKA) = XBEAM !YFLK (NPFLKA) = YBEAM !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 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