From c11cab0804c209ffee3fa0e82c6275464b698ee5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 8 Jan 2012 12:40:27 -0500 Subject: [PATCH 1/2] Fixed error in S(a,b) scattering with elastic Bragg peaks. Closes #65. --- src/cross_section.F90 | 34 +++++++++++++++++----------------- src/physics.F90 | 6 +++++- 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 18706fe1e8..9ec74e6b73 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -191,8 +191,8 @@ contains integer, intent(in) :: index_nuclide ! index into nuclides array integer, intent(in) :: index_sab ! index into sab_tables array - integer :: IE_sab ! index on S(a,b) energy grid - real(8) :: f_sab ! interp factor on S(a,b) energy grid + integer :: IE ! index on S(a,b) energy grid + real(8) :: f ! interp factor on S(a,b) energy grid real(8) :: inelastic ! S(a,b) inelastic cross section real(8) :: elastic ! S(a,b) elastic cross section type(SAB_Table), pointer :: sab => null() @@ -205,17 +205,17 @@ contains ! Get index and interpolation factor for inelastic grid if (p%E < sab % inelastic_e_in(1)) then - IE_sab = 1 - f_sab = ZERO + IE = 1 + f = ZERO else - IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) - f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & - (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) + IE = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) + f = (p%E - sab%inelastic_e_in(IE)) / & + (sab%inelastic_e_in(IE+1) - sab%inelastic_e_in(IE)) end if ! Calculate S(a,b) inelastic scattering cross section - inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & - sab % inelastic_sigma(IE_sab + 1) + inelastic = (ONE - f) * sab % inelastic_sigma(IE) + f * & + sab % inelastic_sigma(IE + 1) ! Check for elastic data if (p % E < sab % threshold_elastic) then @@ -229,24 +229,24 @@ contains ! cross section will be zero elastic = ZERO else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - elastic = sab % elastic_P(IE_sab) / p % E + IE = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + elastic = sab % elastic_P(IE) / p % E end if else ! Determine index on elastic energy grid if (p % E < sab % elastic_e_in(1)) then - IE_sab = 1 + IE = 1 else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + IE = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) end if ! Get interpolation factor for elastic grid - f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & - sab%elastic_e_in(IE_sab)) + f = (p%E - sab%elastic_e_in(IE))/(sab%elastic_e_in(IE+1) - & + sab%elastic_e_in(IE)) ! Calculate S(a,b) elastic scattering cross section - elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & - sab % elastic_P(IE_sab + 1) + elastic = (ONE - f) * sab % elastic_P(IE) + f * & + sab % elastic_P(IE + 1) end if else ! No elastic data diff --git a/src/physics.F90 b/src/physics.F90 index 767dcc2aa5..4d9a20d336 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -608,7 +608,11 @@ contains ! Sample a Bragg edge between 1 and i prob = prn() * sab % elastic_P(i+1) - k = binary_search(sab % elastic_P(1:i+1), i+1, prob) + if (prob < sab % elastic_P(1)) then + k = 1 + else + k = binary_search(sab % elastic_P(1:i+1), i+1, prob) + end if ! Characteristic scattering cosine for this Bragg egg mu = ONE - 2.0*sab % elastic_e_in(k) / p % E From 83a803b83f08029a45e9d636d5b4ad75045cfa12 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 9 Jan 2012 13:14:49 -0500 Subject: [PATCH 2/2] Fixed bug with probability tables (default target velocity). --- src/physics.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/physics.F90 b/src/physics.F90 index 4d9a20d336..b0f1e2e2cd 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -498,6 +498,8 @@ contains ! Sample velocity of target nucleus if (.not. micro_xs(index_nuclide) % use_ptable) then call sample_target_velocity(p, nuc, v_t) + else + v_t = ZERO end if ! Velocity of center-of-mass