Merge pull request #139 from njoy/fix/cpt

Fix/cpt

ReleaseNotes.md

@@ -1,6 +1,16 @@
 # Release Notes—NJOY2016
 Given here are some release notes for NJOY2016. Each release is made through a formal [Pull Request](https://github.com/njoy/NJOY2016/pulls) made on GitHub. There are links in this document that point to each of those Pull Requests, where you can see in great details the changes that were made. Often the Pull Requests are made in response to an [issue](https://github.com/njoy/NJOY2016/issues). In such cases, links to those issues are also given.
 
+## [NJOY2016.52](https://github.com/njoy/NJOY2016/pull/139)
+An error was corrected in the calculation of the Coulomb elastic scattering cross section for incident charged particles. The issue was detected when producing an ACE file for alpha on alpha with ACER but the error also existed in GROUPR. All instances have been corrected. Tests were added for 5 cases to detect this issue in the future.
+
+This pull request addresses issue [\#138](https://github.com/njoy/NJOY2016/issues/138).
+
+## [NJOY2016.51](https://github.com/njoy/NJOY2016/pull/127)
+Changes were made to ACER to increase the number of ZA values that go into the ACE file for a thermal scattering file from 3 to 16 values. This was an issue when one tried to specify the ZA values for a thermal scattering file such as ZrH or SiO2 (the additional ZA values were added by hand as required).
+
+This pull request addresses issue [\#126](https://github.com/njoy/NJOY2016/issues/126).
+
 ## [NJOY2016.50](https://github.com/njoy/NJOY2016/pull/136)
 Changes were made to the ACER routines for photoatomic data to fix an infinite loop when using ENDF/B-VIII.0 data and to fix the scratch array size to accommodate the current size of the ENDF data. An additional test was added to detect this problem in the future.
 

docs/testDescription.md

@@ -258,3 +258,18 @@ title: NJOY2016 Test Descriptions
 
   This test is added following the change in how additional za values are handled in `ACER`. Previously, only 3 were allowed even though thermal scattering ACE file sometimes need more of them (e.g. Zr in ZrH). The user can now specify up to 16 values (the actual number of za values that the ACE file can store). The changes were made to be backwards compatible (test 25 ensures this).
 
+## Test Problem 50-54
+
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/50/input)]
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/51/input)]
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/52/input)]
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/53/input)]
+[[input](https://raw.githubusercontent.com/njoy/NJOY2016/master/tests/54/input)]
+
+  These tests were added following issue [\#138](https://github.com/njoy/NJOY2016/issues/91) in which a charged particle ACE file is produced with NaN values due to an error in the Coulomb elastic scattering cross section for identical particles. These 5 tests cover most of the different possibilities we may encounter:
+- test 50: LAW=5 LTP=12 for identical particles with a spin s = 0 (this file produces NaN values)
+- test 51: LAW=5 LTP=12 for different particles (as expected, this remains the same before and after the fix)
+- test 52: LAW=5 LTP=1 for identical particles with a spin s = 0.5 (as expected, this remains the same before and after the fix)
+- test 53: LAW=5 LTP=1 for identical particles with a spin s = 1 (this changes due to the fix but the original file does no have NaN values, as expected)
+- test 54: LAW=5 LTP=1 for different particles (as expected, this remains the same before and after the fix)
+

src/acefc.f90

@@ -6427,7 +6427,7 @@ subroutine acecpe(next,scr,nin,awr,awp,spi,ne,lidp,ie,il,nes)
                if (lidp.eq.0) sigc=(eta**2/wn**2)/(1-amuu)**2
                if (lidp.eq.1) sigc=((2*eta**2/wn**2)&
                  /(1-amuu**2))*((1+amuu**2)/(1-amuu**2)&
-                 +(-1**i2s)*cos(eta*log((1+amuu)/(1-amuu)))&
+                 +((-1)**i2s)*cos(eta*log((1+amuu)/(1-amuu)))&
                  /(2*spi+1))
                if (ltp.lt.12) pmu=pmu-sigc
                if (iterp.eq.1) then
@@ -7776,7 +7776,7 @@ subroutine ptlegc(c,awp,izai,awt,iza,spi)
    real(kr)::dy,dm,xm,yt,test,check,f,diff,ym,dco
    integer,parameter::imax=20
    real(kr)::x(imax),y(imax)
-   integer,parameter::maxang=2000
+   integer,parameter::maxang=4000
    real(kr)::aco(maxang),cprob(maxang)
    real(kr),parameter::tol1=.001e0_kr
    real(kr),parameter::tol2=.01e0_kr
@@ -7920,7 +7920,7 @@ subroutine coul(x,y,c,awp,izap,awt,izat,spi)
    if (lidp.eq.0) sigc=(eta**2/wn**2)/(1-x)**2
    if (lidp.eq.1) sigc=((2*eta**2/wn**2)&
      /(1-x**2))*((1+x**2)/(1-x**2)&
-     +(-1**i2s)*cos(eta*log((1+x)/(1-x)))/(2*spi+1))
+     +((-1)**i2s)*cos(eta*log((1+x)/(1-x)))/(2*spi+1))
    nt=nint(c(6))
    np=2*nt
    call legndr(x,p,np)

src/groupr.f90

@@ -7714,7 +7714,7 @@ subroutine getdis(e,enext,idisc,yld,ff,nl,ng,iglo,nq,nin,nlg)
          if (lidp.eq.0) sigc=(eta**2/wn**2)/(1-wqp)**2
          if (lidp.eq.1) sigc=((2*eta**2/wn**2)&
            /(1-wqp**2))*((1+wqp**2)/(1-wqp**2)&
-           +(-1**i2s)*cos(eta*log((1+wqp)/(1-wqp)))/(2*spi+1))
+           +((-1)**i2s)*cos(eta*log((1+wqp)/(1-wqp)))/(2*spi+1))
          if (lidp.eq.0) prob=prob/(1-wqp)
          if (lidp.eq.1) prob=prob/(1-wqp*wqp)
          prob=prob+sigc

src/vers.f90

@@ -3,7 +3,7 @@ module version
    ! These values are updated during the NJOY revision-control process.
    implicit none
    private
-   character(8),public::vers='2016.51'
-   character(8),public::vday='01Nov19'
+   character(8),public::vers='2016.52'
+   character(8),public::vday='02Nov19'
 end module version
 

tests/50/CMakeLists.txt

@@ -0,0 +1,16 @@
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/input"
+               "${CMAKE_CURRENT_BINARY_DIR}/input" COPYONLY )
+
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/referenceTape33"
+               "${CMAKE_CURRENT_BINARY_DIR}/referenceTape33" COPYONLY )
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/referenceTape34"
+               "${CMAKE_CURRENT_BINARY_DIR}/referenceTape34" COPYONLY )
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/referenceTape35"
+               "${CMAKE_CURRENT_BINARY_DIR}/referenceTape35" COPYONLY )
+
+configure_file("${RESOURCES}/a-002_He_004-ENDF8.0.endf"
+               "${CMAKE_CURRENT_BINARY_DIR}/tape20" COPYONLY )
+
+add_test( NAME "Test50"
+          COMMAND ${PYTHON_EXECUTABLE} "../execute.py"
+          WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}" )

tests/50/input

@@ -0,0 +1,20 @@
+-- MF6 MT2 LAW=5 LTP=12 LIDP=1 (identical particles) s=0 
+moder
+ 20 -21
+moder
+ -21 -22
+acer
+ -21 -22 0 31 32
+ 1 1 1 .10 /
+ 'ENDF/B-VIII.0 alpha + He4' /
+ 228 0.0 /
+ /
+ /
+acer
+ 0 31 33 34 35
+ 7 1 1 -1 /
+ 'ENDF/B-VIII.0 alpha + He4' /
+viewr
+ 33 43
+stop
+