From: Sunil C (sunilc@barc.gov.in)
Date: Wed Mar 14 2007 - 13:02:14 CET
Hi All
I have the Am-Be neutron spectrum modeled in the source.f attached here,
based on what Nisy had done for Cf-252. I have few questions.
1. Is it necessary to bin the source spectrum strictly in accordance to
the group structure? What happens if I use the spectrum as given in IAEA
tecdoc 413 (2001), the lethargy converted to fluence and normalized to
one, but energy remaining exactly as in the book?
Such as
En (MeV) Fluence
7.94E-02 -
1.00E-01 -
1.25E-01 --
1.58E-01 -
1.99E-01 -
2.51E-01 -
3.16E-01 -
......
....
....
...
Should all the 72 groups be listed in the cumulative index? In this
case the numbers are zero below 7e-2 MeV and above 14.9 MeV. The
cumulative index will show lots of zero below 0.07 MeV and lots of 1.0
above 14.9MeV. In the source attached I have re-binned the spectrum from
the IAEA reference to match the group structure and have retained only
the first groups with 0.0 and 1.0 as cumulative index.
2. Does the neutron group 'NEUGRP' needs to be faithfully passed on from
the source subroutine to the main program? IS the cross section invoked
by this information or is there an internal check?
In the source.f here I am starting from group 3 by setting the 'NEUGRP'
index to 'NEUGRP+2'
And Thank you very much Nisy.
Best Regards
Sunil C
*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
SUBROUTINE SOURCE ( NOMORE )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*----------------------------------------------------------------------*
* Version for isotropic Cf source.
*
c DIMENSION CFCUM(0:72), ENEDGE(73), ANG(0:40), CUMANG(0:40)
DIMENSION CFCUM(0:40), ENEDGE(41), ANG(0:40), CUMANG(0:40)
*----------------------------------------------------------------------*
*
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. /
*----------------------------------------------------------------------*
* Neutron energy group boundaries
DATA ENEDGE /
c Am-Be neutron spectrum
& 1.4918E-02, 1.3499E-02, 1.2214E-02,
& 1.1052E-02, 1.0000E-02, 9.0484E-03, 8.1873E-03, 7.4082E-03,
& 6.7032E-03, 6.0653E-03, 5.4881E-03, 4.9659E-03, 4.4933E-03,
& 4.0657E-03, 3.6788E-03, 3.3287E-03, 3.0119E-03, 2.7253E-03,
& 2.4660E-03, 2.2313E-03, 2.0190E-03, 1.8268E-03, 1.6530E-03,
& 1.4957E-03, 1.3534E-03, 1.2246E-03, 1.1080E-03, 1.0026E-03,
& 9.0718E-04, 8.2085E-04, 7.4274E-04, 6.0810E-04, 4.9787E-04,
& 4.0762E-04, 3.3373E-04, 2.7324E-04, 2.2371E-04, 1.8316E-04,
& 1.4996E-04, 1.2277E-04, 8.6517E-05/
*
* Cumulative Am-Be spectrum in group strcuture
DATA CFCUM / 0.D0,
&1.871941772E-03,4.460515263E-03,
&7.681224425E-03,1.166466537E-02,1.654782161E-02,2.244066112E-02,
&2.947992860E-02,3.777423905E-02,4.725383086E-02,5.774037795E-02,
&6.890897217E-02,8.058700982E-02,9.281815275E-02,1.055467203E-01,
&1.186954005E-01,1.327763145E-01,1.492876327E-01,1.684029556E-01,
&1.901719153E-01,2.149087836E-01,2.455909710E-01,2.844440519E-01,
&3.313109456E-01,3.847880155E-01,4.455136331E-01,5.112235084E-01,
&5.824383965E-01,6.541031298E-01,7.240992794E-01,7.917544886E-01,
&8.561530386E-01,9.152071332E-01,9.564046266E-01,9.777466352E-01,
&9.901464874E-01,9.956112046E-01,9.972629018E-01,9.985685291E-01,
&9.994966257E-01,1.000000000E+00/
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
* Angle cosines
DATA ANG/ -1.0D0, -0.95D0, -0.9D0, -0.85D0, -0.8D0, -0.75D0,
& -0.7D0, -0.65D0, -0.6D0, -0.55D0, -0.5D0, -0.45D0,
& -0.4D0, -0.35D0, -0.3D0, -0.25D0, -0.2D0, -0.15D0,
& -0.1D0, -5.D-02, 0.0D0, 5.D-02, 0.1D0, 0.15D0,
& 0.2D0, 0.25D0, 0.3D0, 0.35D0, 0.4D0, 0.45D0,
& 0.5D0, 0.55D0, 0.6D0, 0.65D0, 0.7D0, 0.75D0,
& 0.8D0, 0.85D0, 0.9D0, 0.95D0, 1.0D0 /
* Cumulative angular distribution
DATA CUMANG / 0.0D0, 4.3787904871D-4, 9.309108266D-4,
& 1.5080318037D-3, 2.1835838416D-3, 2.9743548826D-3,
& 3.8909576476D-3, 4.8923191479D-3, 5.9862770681D-3,
& 7.1813938468D-3, 8.4870236975D-3, 9.9133858209D-3,
& 1.1471644393D-2, 1.3097707545D-2, 1.4776914040D-2,
& 1.6511000723D-2, 1.8301761200D-2, 2.0151047694D-2,
& 2.2060772967D-2, 2.4032912292D-2, 2.6069505497D-2,
& 2.9084203728D-2, 3.3546756876D-2, 4.0152519402D-2,
& 4.9930800378D-2, 6.4405250301D-2, 8.5831276094D-2,
& 0.1175474784D0, 0.1569148680D0, 0.1981130258D0,
& 0.2412270912D0, 0.2863461629D0, 0.3335634830D0,
& 0.3829766299D0, 0.4346877199D0, 0.4919529577D0,
& 0.5559660349D0, 0.6275220820D0, 0.7075099237D0,
& 0.8189694928D0, 1.0D0 /
*----------------------------------------------------------------------*
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
*----------------------------------------------------------------------*
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)') "Californium source. Version corrected"
write(lunout,*)
write(lunout,*)
write(lunout,'(a,132a)') ("*",i=1,132)
write(lunout,*)
END IF
* |
* +-------------------------------------------------------------------*
* Sample the energy group
XI = FLRNDM(DUMMY)
DO 500 K = 1, 40
IF(XI .LE. CFCUM(K)) THEN
NEUGRP = K+2
ENERGY = ENEDGE(K) -
& (XI-CFCUM(K-1))*(ENEDGE(K)-ENEDGE(K+1))/(CFCUM(K)-CFCUM(K-1))
GO TO 600
END IF
500 CONTINUE
STOP ' Failed to sample the energy group'
600 CONTINUE
kount=kount+1
*
* Sample the cosine of the polar angle
XI = FLRNDM(XI)
DO 300 K = 1, 40
IF(XI .LE. CUMANG(K)) THEN
COSTHE = ANG(K-1) +
& (XI - CUMANG(K-1))*(ANG(K)-ANG(K-1))/(CUMANG(K)-CUMANG(K-1))
GO TO 400
END IF
300 CONTINUE
STOP ' Failed to sample the polar angle cosine'
400 CONTINUE
SINTHE = SQRT(ONEONE - COSTHE**2)
IF(FLRNDM(COSTHE) .LE. HLFHLF) SINTHE = -SINTHE
* Samples the azimuthal angle
PHI = TWOPIP * FLRNDM(COSTHE)
COSPHI = COS(PHI)
* 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
* +-------------------------------------------------------------------*
* | Heavy ion:
IF ( IJBEAM .EQ. -2 ) THEN
IJHION = IPROZ * 1000 + IPROA
IJHION = IJHION * 100 + KXHEAV
IONID = IJHION
CALL DCDION ( IONID )
CALL SETION ( IONID )
ILOFLK (NPFLKA) = IJHION
* |
* +-------------------------------------------------------------------*
* | Normal hadron:
ELSE
IONID = IJBEAM
ILOFLK (NPFLKA) = IJBEAM
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
IGROUP (NPFLKA) = NEUGRP
* Kinetic energy of the particle (GeV)
TKEFLK (NPFLKA) = ENERGY
* Particle momentum
PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
* + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
* Cosines (tx,ty,tz) (make sure they are properly normalized)
TZFLK (NPFLKA) = COSTHE
TYFLK (NPFLKA) = SINTHE * COSPHI
TXFLK (NPFLKA) = SQRT(ONEONE-TZFLK(NPFLKA)**2-TYFLK(NPFLKA)**2)
IF(FLRNDM(COSPHI) .LE. HLFHLF) TXFLK (NPFLKA) = -TXFLK (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
* Flag this is prompt radiation
LRADDC (NPFLKA) = .FALSE.
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
This archive was generated by hypermail 2.1.6 : Thu Mar 15 2007 - 15:59:33 CET