*$ CREATE SOURCE.FOR *COPY SOURCE * *=== source ===========================================================* * SUBROUTINE SOURCE ( NOMORE ) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' DIMENSION ENRG(0:22) DIMENSION CUMPR(0:22) * *----------------------------------------------------------------------* * * * Copyright (C) 1990-2009 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 08-feb-09 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. / *======================================================================* C neutron spectrum of Am-Be from TRS 318 DATA ENRG/ & 7.94E-05, 1.00E-04, 1.26E-04, 1.58E-04, 2.00E-04, & 2.51E-04, 3.16E-04, 3.98E-04, 5.01E-04, 6.31E-04, & 7.94E-04, 1.00E-03, 1.26E-03, 1.58E-03, 2.00E-03, & 2.51E-03, 3.16E-03, 3.98E-03, 5.01E-03, 6.31E-03, & 7.94E-03, 1.00E-02, 1.26E-02 / C CUMULATIVE PROBABAILITY DATA CUMPR/ & 0, 3.34E-03, 7.55E-03, 1.29E-02, 1.96E-02, & 2.81E-02, 3.88E-02, 5.14E-02, 6.62E-02, 8.36E-02, & 1.04E-01, 1.28E-01, 1.56E-01, 1.85E-01, 2.27E-01, & 2.85E-01, 3.81E-01, 5.17E-01, 6.76E-01, 8.22E-01, & 9.56E-01, 9.97E-01, 1.00E+00 / * * * BASIC VERSION * * * *======================================================================* NOMORE = 0 * +-------------------------------------------------------------------* * | First call initializations: c open(99,file='sourceout.txt',status='unknown') IF ( LFIRST ) THEN * | *** The following 3 cards are mandatory *** TKESUM = ZERZER LFIRST = .FALSE. LUSSRC = .TRUE. * | *** User initialization *** C write(lunout,*) C write(lunout,'(a,132a)') ("*",i=1,132) C write(lunout,*) C write(lunout,*) write(lunout,'(a)') "Thermal neutron source" C write(lunout,*) C write(lunout,*) C write(lunout,'(a,132a)') ("*",i=1,132) C 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 * 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 * 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) xi = flrndm(dummy) do i = 1, 23 if(xi .gt. CUMPR(i-1) .and. xi .le. CUMPR(i)) then ENERGY= ENRG(i-1)+(ENRG(i)-ENRG(i-1))* & (xi - CUMPR(i-1))/(CUMPR(i)-CUMPR(i-1)) c write(99,*) 'Energy and group are ', ENERGY, i go to 1 end if end do 1 continue C TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID) * Particle momentum ENERGY = 2.5E-11 TKEFLK(NPFLKA) = ENERGY C PMOFLK (NPFLKA) = PBEAM PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA) & + TWOTWO * AM (IONID) ) ) * Cosines (tx,ty,tz) * Cosines (tx,ty,tz) COSTH = ONEONE-0.105572809*FLRNDM(TTT) IF(COSTH.GT.ONEONE) COSTH=ONEONE IF(COSTH.LT.-ONEONE) COSTH=-ONEONE PHI = TWOPIP*FLRNDM(SSS) WBEAM=SQRT(ONEONE-COSTH*COSTH)*COS(PHI) VBEAM=SQRT(ONEONE-COSTH*COSTH)*SIN(PHI) UBEAM=COSTH 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 X = -10.0 Y = 0.0 Z = 5 XBEAM = X YBEAM = Y ZBEAM = Z C XBEAM = 20.0+ (50-20)*FLRNDM(XXX) C YBEAM = 20.0+ (50-20)*FLRNDM(YYY) C ZBEAM = 20.0+ (50-20)*FLRNDM(ZZZ) 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