Merge pull request #275 from njoy/fix/pn-law2

Photonuclear heating number for LAW=2

ReleaseNotes.md

@@ -4,6 +4,7 @@ Given here are some release notes for NJOY2016. Each release is made through a f
 ## [NJOY2016.69](https://github.com/njoy/NJOY2016/pull/xxx)
 This update fixes a number of minor issues:
   - PURR now writes Bondarenko data obtained from the probability tables to MF2 MT152 instead of the Bondarenko data obtained from the direct sampled cross sections (for very low dilutions, the Bondarenko data obtained using these two methods does not align, with the direct sampled data leading to extremely low P1 values). When comparing with the Bondarenko data at low dilutions obtained with UNRESR, the Bondarenko data obtained from the probability table directly seems to be the best.
+  - MF6 LAW=2 represents discrete two body scattering in which only angular distribution data is given (knowing that the outgoing energy of the secondary particle can be determined through kinematics when the angle is known). When calculating heating numbers based on LAW=2, ACER assumes that the yield of the secondary particle is 1, which is correct in all cases except when MT5 is used as a lumped reaction. Heating numbers in ACER for photonuclear files using LAW=2 in an MT5 entry are now correctly multiplied by the yield. A warning message is printed out whenever this situation is detected. Test 78 was added as part of this correction.
 
 ## [NJOY2016.68](https://github.com/njoy/NJOY2016/pull/264)
 This update fixes a number of minor issues:

src/acepn.f90

@@ -52,14 +52,14 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
    integer::iint,nn,kk,m,intt,last,lf,jnt,ja,jb,ipp,irr
    integer::lee,lle,nd,na,ncyc,ng,ig,nnr,nnp,mf,mt
    integer::ipt,ntrp,pxs,phn,mtrp,tyrp,lsigp,sigp,landp,andp,ldlwp,dlwp
-   integer::izarec,nl,iil,nexn,nexd,ki
+   integer::izarec,nl,iil,nexn,nexd,ki,ilaw2mt5
    integer::nle
    integer::imu,intmu,nmu
    real(kr)::emc2,e,enext,s,y,ynext,heat,en,ep,g,h,epl
    real(kr)::tneut,tphot,tprot,tdeut,ttrit,the3,the4,thresh
    real(kr)::ss,tt,ubar,sum,renorm,ebar,hh,u,theta,x,anorm
    real(kr)::ee,amass,avadd,avlab,avll,test,rkal,akal
-   real(kr)::eavi,avl,avcm,sign,dele,avav,zaid,gl,awp,awr,q
+   real(kr)::eavi,avl,avcm,sign,dele,avav,zaid,gl,awp,awr,q,yld
    real(kr)::av,del
    real(kr)::awprec,awpp
    integer,parameter::mmax=1000
@@ -78,6 +78,8 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
    character(66)::text
    emc2=amassn*amu*clight*clight/ev/emev
    tvn=1
+   ilaw2mt5=0
+   yld=1.
 
    nxsd=0
    jxsd=0
@@ -1095,8 +1097,8 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
                              sigfig(renorm*xss(nex+1+2*n+ii),9,0)  ! CDF(ii)
                         enddo
                         nex=nex+2+3*n       ! index for the next distribution
-                        e=xss(ie+iie)
-                        scr(llht+6+2*iie)=e
+                        e=xss(ie+iie)       ! e from xss, thus MeV
+                        scr(llht+6+2*iie)=e ! store it in scr, so energy in MeV
                         scr(llht+7+2*iie)=(awr-awpp)*(e+q)/awr
                      enddo
                      ! add in contribution to heating
@@ -1105,13 +1107,32 @@ subroutine acephn(nendf,npend,nace,ndir,matd,tempd,iprint,mcnpx,&
                      do ie=it,nes
                         e=xss(esz+ie-1)/emev
                         call terpa(h,e,en,idis,scr(llht),npp,nrr)
+                        !--get the yield and modify heating
+                        !--beware: e is already in MeV, scr(llht) has been
+                        !--modified to be in MeV already (see above)
+                        call terpa(yld,e*emev,en,idis,scr,npp,nrr)
+                        h=yld*h ! just in case, multiply by yield
+                        if (ilaw2mt5.eq.0.and.mt.eq.5.and.yld.ne.one) then
+                           ilaw2mt5=1
+                           write(text,'(''yield different from 1 for outgoing particle '',i5,''.'')')ip
+                           if (izarec.eq.ip) then
+                              call mess('acephn',&
+                                        'law=2 (discrete two-body scattering) used in mt=5 recoil',&
+                                        text)
+                           else
+                              call mess('acephn',&
+                                        'law=2 (discrete two-body scattering) used in mt=5',&
+                                        text)
+                           endif
+                        endif
                         ss=0
                         if (ie.ge.iaa) ss=xss(2+k+ie-iaa)
                         xss(phn+2+ie-it)=xss(phn+2+ie-it)+h*ss
                         if (xss(tot+ie-1).ne.zero)&
                           xss(thn+ie-1)=xss(thn+ie-1)&
                           +h*ss/xss(tot+ie-1)
                      enddo
+                     ilaw2mt5=0
                   else
                      call skip6(nin,0,0,scr,law)
                   endif

tests/78/CMakeLists.txt

@@ -0,0 +1,15 @@
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/input"
+               "${CMAKE_CURRENT_BINARY_DIR}/input" COPYONLY )
+
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/referenceTape31"
+               "${CMAKE_CURRENT_BINARY_DIR}/referenceTape31" COPYONLY )
+
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/referenceTape34"
+               "${CMAKE_CURRENT_BINARY_DIR}/referenceTape34" COPYONLY )
+
+configure_file("${RESOURCES}/g-002-He-3-IAEA.endf"
+               "${CMAKE_CURRENT_BINARY_DIR}/tape20" COPYONLY )
+
+add_test( NAME "Test78"
+          COMMAND ${Python3_EXECUTABLE} "../execute.py"
+          WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}" )

tests/78/input

@@ -0,0 +1,27 @@
+-- moder
+moder
+20 -30
+-- reconr
+reconr
+-30 -21/
+'He3 IAEA photonuclear data' /
+225 /
+0.001 /
+0/
+moder
+-21 90
+-- acer
+acer
+-30 -21 0 31 32                                                 / Card 1
+5 0 1 .31                                                       / Card 2
+'He3 IAEA photonuclear data'                                    / Card 3
+225                                                             / Card 12
+-- acer check/plot/correct
+acer
+0 31 33 34 35
+7 1 1 -1/
+'He3 IAEA photonuclear data - check tape'
+-- viewr plots from check tape
+viewr
+33 43
+stop