*$ 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. * * * *----------------------------------------------------------------------* *----------------------------------------------------------------------* INCLUDE '(BEAMCM)' INCLUDE '(FHEAVY)' INCLUDE '(FLKSTK)' INCLUDE '(IOIOCM)' INCLUDE '(LTCLCM)' INCLUDE '(PAPROP)' INCLUDE '(SOURCM)' INCLUDE '(SUMCOU)' * LOGICAL LFIRST * * Version for isotropic Heavy ion source * DIMENSION CUMENE(0:100), ENELOW(0:99), ENEUP(0:99) SAVE LFIRST DATA LFIRST / .TRUE. / *======================================================================* * * * BASIC VERSION * * * *======================================================================* * Heavy ion energy group(Gev/n) DATA ENELOW / & 0.00000E+00, 1.00000E-03, 1.14980E-03, 1.32190E-03, 1.51990E-03, & 1.74750E-03, 2.00920E-03, 2.31010E-03, 2.65610E-03, 3.05390E-03, & 3.51120E-03, 4.03700E-03, 4.64160E-03, 5.33670E-03, 6.13590E-03, & 7.05480E-03, 8.11130E-03, 9.32600E-03, 1.07230E-02, 1.23280E-02, & 1.41750E-02, 1.62980E-02, 1.87380E-02, 2.15440E-02, 2.47710E-02, & 2.84800E-02, 3.27450E-02, 3.76490E-02, 4.32880E-02, 4.97700E-02, & 5.72240E-02, 6.57930E-02, 7.56460E-02, 8.69750E-02, 1.00000E-01, & 1.14980E-01, 1.32190E-01, 1.51990E-01, 1.74750E-01, 2.00920E-01, & 2.31010E-01, 2.65610E-01, 3.05390E-01, 3.51120E-01, 4.03700E-01, & 4.64160E-01, 5.33670E-01, 6.13590E-01, 7.05480E-01, 8.11130E-01, & 9.32600E-01, 1.07230E+00, 1.23280E+00, 1.41750E+00, 1.62980E+00, & 1.87380E+00, 2.15440E+00, 2.47710E+00, 2.84800E+00, 3.27450E+00, & 3.76490E+00, 4.32880E+00, 4.97700E+00, 5.72240E+00, 6.57930E+00, & 7.56460E+00, 8.69750E+00, 1.00000E+01, 1.14980E+01, 1.32190E+01, & 1.51990E+01, 1.74750E+01, 2.00920E+01, 2.31010E+01, 2.65610E+01, & 3.05390E+01, 3.51120E+01, 4.03700E+01, 4.64160E+01, 5.33670E+01, & 6.13590E+01, 7.05480E+01, 8.11130E+01, 9.32600E+01, 1.07230E+02, & 1.23280E+02, 1.41750E+02, 1.62980E+02, 1.87380E+02, 2.15440E+02, & 2.47710E+02, 2.84800E+02, 3.27450E+02, 3.76490E+02, 4.32880E+02, & 4.97700E+02, 5.72240E+02, 6.57930E+02, 7.56460E+02, 8.69750E+02/ DATA ENEUP / & 1.00000E-03, 1.14980E-03, 1.32190E-03, 1.51990E-03, 1.74750E-03, & 2.00920E-03, 2.31010E-03, 2.65610E-03, 3.05390E-03, 3.51120E-03, & 4.03700E-03, 4.64160E-03, 5.33670E-03, 6.13590E-03, 7.05480E-03, & 8.11130E-03, 9.32600E-03, 1.07230E-02, 1.23280E-02, 1.41750E-02, & 1.62980E-02, 1.87380E-02, 2.15440E-02, 2.47710E-02, 2.84800E-02, & 3.27450E-02, 3.76490E-02, 4.32880E-02, 4.97700E-02, 5.72240E-02, & 6.57930E-02, 7.56460E-02, 8.69750E-02, 1.00000E-01, 1.14980E-01, & 1.32190E-01, 1.51990E-01, 1.74750E-01, 2.00920E-01, 2.31010E-01, & 2.65610E-01, 3.05390E-01, 3.51120E-01, 4.03700E-01, 4.64160E-01, & 5.33670E-01, 6.13590E-01, 7.05480E-01, 8.11130E-01, 9.32600E-01, & 1.07230E+00, 1.23280E+00, 1.41750E+00, 1.62980E+00, 1.87380E+00, & 2.15440E+00, 2.47710E+00, 2.84800E+00, 3.27450E+00, 3.76490E+00, & 4.32880E+00, 4.97700E+00, 5.72240E+00, 6.57930E+00, 7.56460E+00, & 8.69750E+00, 1.00000E+01, 1.14980E+01, 1.32190E+01, 1.51990E+01, & 1.74750E+01, 2.00920E+01, 2.31010E+01, 2.65610E+01, 3.05390E+01, & 3.51120E+01, 4.03700E+01, 4.64160E+01, 5.33670E+01, 6.13590E+01, & 7.05480E+01, 8.11130E+01, 9.32600E+01, 1.07230E+02, 1.23280E+02, & 1.41750E+02, 1.62980E+02, 1.87380E+02, 2.15440E+02, 2.47710E+02, & 2.84800E+02, 3.27450E+02, 3.76490E+02, 4.32880E+02, 4.97700E+02, & 5.72240E+02, 6.57930E+02, 7.56460E+02, 8.69750E+02, 1.00000E+03/ * Normalized spectrum DATA CUMENE / & 0.00000E+00, & 9.18166E-06, 1.07735E-05, 1.28897E-05, 1.57070E-05, 1.94544E-05, & 2.44404E-05, 3.10745E-05, 3.99033E-05, 5.16524E-05, 6.72883E-05, & 8.81044E-05, 1.15825E-04, 1.52739E-04, 2.01908E-04, 2.67410E-04, & 3.54670E-04, 4.70903E-04, 6.25735E-04, 8.31649E-04, 1.10575E-03, & 1.46971E-03, 1.95215E-03, 2.59057E-03, 3.43300E-03, 4.54006E-03, & 5.98934E-03, 7.87688E-03, 1.03211E-02, 1.34642E-02, 1.74782E-02, & 2.25612E-02, 2.89433E-02, 3.68810E-02, 4.66546E-02, 5.85661E-02, & 7.29078E-02, 9.00005E-02, 1.10111E-01, 1.33484E-01, 1.60298E-01, & 1.90654E-01, 2.24543E-01, 2.61837E-01, 3.02291E-01, 3.45521E-01, & 3.91010E-01, 4.38135E-01, 4.86189E-01, 5.34412E-01, 5.82036E-01, & 6.28334E-01, 6.72592E-01, 7.14310E-01, 7.53000E-01, 7.88351E-01, & 8.20216E-01, 8.48552E-01, 8.73416E-01, 8.94982E-01, 9.13477E-01, & 9.29173E-01, 9.42361E-01, 9.53350E-01, 9.62428E-01, 9.69876E-01, & 9.75946E-01, 9.80862E-01, 9.84824E-01, 9.87999E-01, 9.90536E-01, & 9.92553E-01, 9.94153E-01, 9.95417E-01, 9.96414E-01, 9.97199E-01, & 9.97814E-01, 9.98297E-01, 9.98674E-01, 9.98969E-01, 9.99199E-01, & 9.99378E-01, 9.99518E-01, 9.99627E-01, 9.99712E-01, 9.99777E-01, & 9.99828E-01, 9.99868E-01, 9.99899E-01, 9.99923E-01, 9.99941E-01, & 9.99956E-01, 9.99967E-01, 9.99975E-01, 9.99982E-01, 9.99987E-01, & 9.99991E-01, 9.99994E-01, 9.99997E-01, 9.99999E-01, 1.00000E+00/ * +-------------------------------------------------------------------* NOMORE = 0 * +-------------------------------------------------------------------* IF BLOCK * | 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 Proton 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 * +-------------------------------------------------------------------* * +--------------------------------------------------------------------* * | (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) * Sample the energy group XI = FLRNDM(DUMMY) DO K = 0 ,99 IF (XI .GT. CUMENE(K) .AND. XI .LE. CUMENE(K+1)) THEN ENERGY = ((ENEUP(K) + ENELOW(K))/ TWOTWO) GO TO 600 END IF END DO 600 CONTINUE 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) ) ) * +-----------------------------------------------------------------* Call SFLOOD ( XXX, YYY, ZZZ, UXXX, VXXX, WZZZ ) UBEAM = UXXX VBEAM = VYYY WBEAM = WZZZ * 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 ( radius of sphere =10 cm) XBEAM = XXX * 10.D0 YBEAM = YYY * 10.D0 ZBEAM = ZZZ * 10.D0 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