*$ 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. * * * *----------------------------------------------------------------------* * Version for isotropic Heavy ion source * DIMENSION CUMENE(0:98), ENELOW(99), ENEUP(99) *----------------------------------------------------------------------* 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 * * * *======================================================================* * Heavy Ion energy group(Gev/n) DATA ENELOW / & 1.00000D-03, 1.14980D-03, 1.32190D-03, 1.51990D-03, 1.74750D-03, & 2.00920D-03, 2.31010D-03, 2.65610D-03, 3.05390D-03, 3.51120D-03, & 4.03700D-03, 4.64160D-03, 5.33670D-03, 6.13590D-03, 7.05480D-03, & 8.11130D-03, 9.32600D-03, 1.07230D-02, 1.23280D-02, 1.41750D-02, & 1.62980D-02, 1.87380D-02, 2.15440D-02, 2.47710D-02, 2.84800D-02, & 3.27450D-02, 3.76490D-02, 4.32880D-02, 4.97700D-02, 5.72240D-02, & 6.57930D-02, 7.56460D-02, 8.69750D-02, 1.00000D-01, 1.14980D-01, & 1.32190D-01, 1.51990D-01, 1.74750D-01, 2.00920D-01, 2.31010D-01, & 2.65610D-01, 3.05390D-01, 3.51120D-01, 4.03700D-01, 4.64160D-01, & 5.33670D-01, 6.13590D-01, 7.05480D-01, 8.11130D-01, 9.32600D-01, & 1.07230D+00, 1.23280D+00, 1.41750D+00, 1.62980D+00, 1.87380D+00, & 2.15440D+00, 2.47710D+00, 2.84800D+00, 3.27450D+00, 3.76490D+00, & 4.32880D+00, 4.97700D+00, 5.72240D+00, 6.57930D+00, 7.56460D+00, & 8.69750D+00, 1.00000D+01, 1.14980D+01, 1.32190D+01, 1.51990D+01, & 1.74750D+01, 2.00920D+01, 2.31010D+01, 2.65610D+01, 3.05390D+01, & 3.51120D+01, 4.03700D+01, 4.64160D+01, 5.33670D+01, 6.13590D+01, & 7.05480D+01, 8.11130D+01, 9.32600D+01, 1.07230D+02, 1.23280D+02, & 1.41750D+02, 1.62980D+02, 1.87380D+02, 2.15440D+02, 2.47710D+02, & 2.84800D+02, 3.27450D+02, 3.76490D+02, 4.32880D+02, 4.97700D+02, & 5.72240D+02, 6.57930D+02, 7.56460D+02, 8.69750D+02/ DATA ENEUP / & 1.14980D-03, 1.32190D-03, 1.51990D-03, 1.74750D-03, 2.00920D-03, & 2.31010D-03, 2.65610D-03, 3.05390D-03, 3.51120D-03, 4.03700D-03, & 4.64160D-03, 5.33670D-03, 6.13590D-03, 7.05480D-03, 8.11130D-03, & 9.32600D-03, 1.07230D-02, 1.23280D-02, 1.41750D-02, 1.62980D-02, & 1.87380D-02, 2.15440D-02, 2.47710D-02, 2.84800D-02, 3.27450D-02, & 3.76490D-02, 4.32880D-02, 4.97700D-02, 5.72240D-02, 6.57930D-02, & 7.56460D-02, 8.69750D-02, 1.00000D-01, 1.14980D-01, 1.32190D-01, & 1.51990D-01, 1.74750D-01, 2.00920D-01, 2.31010D-01, 2.65610D-01, & 3.05390D-01, 3.51120D-01, 4.03700D-01, 4.64160D-01, 5.33670D-01, & 6.13590D-01, 7.05480D-01, 8.11130D-01, 9.32600D-01, 1.07230D+00, & 1.23280D+00, 1.41750D+00, 1.62980D+00, 1.87380D+00, 2.15440D+00, & 2.47710D+00, 2.84800D+00, 3.27450D+00, 3.76490D+00, 4.32880D+00, & 4.97700D+00, 5.72240D+00, 6.57930D+00, 7.56460D+00, 8.69750D+00, & 1.00000D+01, 1.14980D+01, 1.32190D+01, 1.51990D+01, 1.74750D+01, & 2.00920D+01, 2.31010D+01, 2.65610D+01, 3.05390D+01, 3.51120D+01, & 4.03700D+01, 4.64160D+01, 5.33670D+01, 6.13590D+01, 7.05480D+01, & 8.11130D+01, 9.32600D+01, 1.07230D+02, 1.23280D+02, 1.41750D+02, & 1.62980D+02, 1.87380D+02, 2.15440D+02, 2.47710D+02, 2.84800D+02, & 3.27450D+02, 3.76490D+02, 4.32880D+02, 4.97700D+02, 5.72240D+02, & 6.57930D+02, 7.56460D+02, 8.69750D+02, 1.00000D+03/ * Normalized spectrum DATA CUMENE / 0.D0, & 3.11051E-06, 5.47369E-06, 8.61672E-06, 1.27980E-05, 1.83609E-05, & 2.57626E-05, 3.56111E-05, 4.87150E-05, 6.61565E-05, 8.93771E-05, & 1.20292E-04, 1.61460E-04, 2.16289E-04, 2.89316E-04, 3.86564E-04, & 5.16088E-04, 6.88345E-04, 9.17632E-04, 1.22215E-03, 1.62593E-03, & 2.16056E-03, 2.86662E-03, 3.79548E-03, 5.01318E-03, 6.60198E-03, & 8.66372E-03, 1.13218E-02, 1.47264E-02, 1.90525E-02, 2.45053E-02, & 3.13167E-02, 3.97438E-02, 5.00691E-02, 6.25734E-02, 7.75723E-02, & 9.53420E-02, 1.16151E-01, 1.40217E-01, 1.67702E-01, 1.98673E-01, & 2.33097E-01, 2.70837E-01, 3.11620E-01, 3.55035E-01, 4.00566E-01, & 4.47584E-01, 4.95385E-01, 5.43224E-01, 5.90360E-01, 6.36028E-01, & 6.79661E-01, 7.20665E-01, 7.58623E-01, 7.93271E-01, 8.24463E-01, & 8.52155E-01, 8.76442E-01, 8.97490E-01, 9.15532E-01, 9.30832E-01, & 9.43692E-01, 9.54402E-01, 9.63254E-01, 9.70518E-01, 9.76439E-01, & 9.81239E-01, 9.85107E-01, 9.88213E-01, 9.90694E-01, 9.92670E-01, & 9.94238E-01, 9.95479E-01, 9.96458E-01, 9.97229E-01, 9.97836E-01, & 9.98311E-01, 9.98684E-01, 9.98976E-01, 9.99203E-01, 9.99381E-01, & 9.99520E-01, 9.99628E-01, 9.99712E-01, 9.99778E-01, 9.99829E-01, & 9.99869E-01, 9.99899E-01, 9.99923E-01, 9.99942E-01, 9.99956E-01, & 9.99968E-01, 9.99976E-01, 9.99983E-01, 9.99988E-01, 9.99992E-01, & 9.99996E-01, 9.99998E-01, 1.00000E+00/ * +-------------------------------------------------------------------* 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,*) '--->modified source.f for Silicon & in solar min' 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 * +-------------------------------------------------------------------* * Sample the energy group XI = FLRNDM(DUMMY) YI = FLRNDM(DUMMY1) DO 500 K = 1 ,99 IF (XI .LE. CUMENE(K)) THEN ENERGY = (YI * (ENEUP(K) - ENELOW(K)) + ENELOW(K)) * 28.D0 GO TO 600 END IF 500 CONTINUE STOP ' Failed to sample the energy group ' 600 CONTINUE count=count+1 NPFLKA = NPFLKA + 1 WTFLK (NPFLKA) = ONEONE WEIPRI = WEIPRI + WTFLK (NPFLKA) * +--------------------------------------------------------------------* * | (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. * | * +-------------------------------------------------------------------* * | Normal hadron: ELSE IONID = IJBEAM ILOFLK (NPFLKA) = IJBEAM * | Flag this is prompt radiation * LRADDC (NPFLKA) = .FALSE. 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) 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) ) ) * +-----------------------------------------------------------------* RANPHI = FLRNDM(U) PHID = TWOPIP * RANPHI RANCOS = FLRNDM(V) COSD = Oneone - Twotwo * RANCOS SIND = SQRT(Oneone - COSD * COSD) UBEAM = -SIND * COS(PHID) VBEAM = -SIND * SIN(PHID) WBEAM = -COSD * 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 * +-----------------------------------------------------------------* RAD = 285 ZS = 500 RANPHI = FLRNDM(P) PHIS = TWOPIP * RANPHI RANR = FLRNDM(R) RS = RAD * SQRT(RANR) XS = RS * COS(PHIS) YS = RS * SIN(PHIS) XBEAM = XS * COSD * COS(PHID) - YS * SIN(PHID) + ZS * SIND & * COS(PHID) YBEAM = XS * COSD * SIN(PHID) + YS * COS(PHID) + ZS * SIND & * SIN(PHID) ZBEAM = -XS * SIND + ZS * COSD * 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