*$ CREATE SOURCE.FOR *COPY SOURCE * *=== source ===========================================================* * SUBROUTINE SOURCE ( NOMORE ) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' DIMENSION ENER(47), CUMER(0:46) * DATA ENER / 0.0003d0, 0.0006d0, 0.0010d0, 0.0015d0, 0.0020d0, * & 0.0025d0, 0.0030d0, 0.0035d0,0.0040d0, 0.0050d0, 0.0055d0, * & 0.0060d0, 0.0066d0, 0.0070d0, 0.0075d0, 0.0080d0, 0.0085d0, * & 0.0090d0, 0.0095d0, 0.0100d0, 0.0105d0, 0.011d0 / * DATA CUMER / 0.0d0, 0.0598186049d0, 0.1316009308d0, 0.1964019333, * & 0.2512281958d0, 0.3145275176d0, 0.3825606142d0, 0.4732780397, * & 0.5455576108d0, 0.6190802952d0, 0.6952980488d0, 0.7388567337, * & 0.7808242337d0, 0.8167153967d0, 0.8566044116d0, 0.9014758239, * & 0.9401018358d0, 0.9550589733d0, 0.9675199395d0, 0.9824671321, * & 0.9949380432d0, 1.0D0 / DATA ENER / 0.00d0, 0.00105d0, 0.00116d0, & 0.00126d0, & 0.00135d0, 0.00157d0, 0.00171d0, 0.00187d0, 0.00201d0, & 0.00219d0, 0.00231d0, 0.00248d0, 0.00259d0, & 0.00273d0, 0.00281d0, 0.00289d0, 0.00300d0, & 0.00304d0, 0.00314d0, 0.00325d0, 0.00338d0, & 0.00355d0, 0.00388d0, 0.00421d0, 0.00455d0, & 0.00485d0, 0.00518d0, 0.00544d0, 0.00570d0, & 0.00598d0, 0.00626d0, 0.00654d0, 0.00682d0, & 0.00706d0, 0.00735d0, 0.00760d0, 0.00787d0, & 0.00810d0, 0.00833d0, 0.00863d0, 0.00890d0, & 0.00924d0, 0.00964d0, 0.00996d0, 0.01027d0, & 0.01065d0, 0.01084 / DATA CUMER / 0.0d0, 0.026992d0, 0.048099d0, & 0.070133d0, 0.091081d0, 0.110400d0, 0.127472d0, & 0.147095d0, 0.169269d0, 0.188885d0, 0.209738d0, & 0.231207d0, 0.255072d0, 0.282702d0, 0.312214d0, & 0.346600d0, 0.383538d0, 0.423649d0, 0.463371d0, & 0.502394d0, 0.537777d0, 0.571532d0, 0.600403d0, & 0.630662d0, 0.662927d0, 0.697043d0, 0.730379d0, & 0.756589d0, 0.778691d0, 0.793590d0, 0.808174d0, & 0.827165d0, 0.847543d0, 0.863815d0, 0.882171d0, & 0.900134d0, 0.920727d0, 0.940694d0, 0.955161d0, & 0.966603d0, 0.972390d0, 0.977783d0, 0.985104d0, & 0.992651d0, 0.997560d0, 0.999520d0, 1.0d0 / *----------------------------------------------------------------------* * * * 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 * * * *----------------------------------------------------------------------* * INCLUDE '(BEAMCM)' INCLUDE '(FHEAVY)' INCLUDE '(FLKSTK)' INCLUDE '(IOIOCM)' INCLUDE '(LTCLCM)' INCLUDE '(PAPROP)' INCLUDE '(SOURCM)' INCLUDE '(SUMCOU)' * LOGICAL LFIRST * SAVE LFIRST DATA LFIRST / .TRUE. / *======================================================================* * * * 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(lunout,*) write(lunout,'(a,132a)') ("*",i=1,132) write(lunout,*) write(lunout,*) write(lunout,'(a)') "Pu-Be source" write(lunout,*) write(lunout,*) write(lunout,'(a,132a)') ("*",i=1,132) write(lunout,*) END IF *------------------------------------------------------------------------* ** the energy of Pu be source is taken from the litrecture ** XI=FLRNDM(DUMMY) DO 400 K=1,46 if (XI .gt. CUMER(K-1) .and. XI .LE. CUMER(K) ) THEN ENERGY= ENER(K)+ & (XI-CUMER(K-1))*(ENER(K+1)-ENER(K))/(CUMER(K)-CUMER(K-1)) GO TO 500 end if 400 CONTINUE STOP "FAILED TO SAMPLE THE ENERGIES" 500 CONTINUE kount= kount+1 ** sample the polar angles COSTHE = TWOTWO * FLRNDM(XXX) - ONEONE SINTHE=SQRT(ONEONE-COSTHE**2) ** SAMPLE THE AZIMUTHAL ANGLE PHI=TWOPIP * FLRNDM(CCC) COSPHI=COS(PHI) * | * +-------------------------------------------------------------------* * 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) TKEFLK (NPFLKA)= ENERGY * 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) TZFLK (NPFLKA) = COSTHE TXFLK (NPFLKA) = SINTHE * COSPHI TYFLK (NPFLKA) = SQRT( ONEONE - TXFLK (NPFLKA)**2 & - TZFLK (NPFLKA)**2 ) IF (FLRNDM(COSPHI) .LE. HLFHLF) TYFLK (NPFLKA)= -TYFLK (NPFLKA) * Polarization cosines: TXPOL (NPFLKA) = -TWOTWO TYPOL (NPFLKA) = +ZERZER TZPOL (NPFLKA) = +ZERZER * Particle coordinates 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