LAW=33 is up to it again

[whaeck]

LAW=33 is an undocumented exit law in the ACE file, mainly used for outgoing charged particles.

Normally, the format of any exit law involves the energy values and probabilities of the law occurring. In the case of LAW=33, it appears that these values were printed in the wrong order:

• correct order: all energy values followed by all probabilities (equal to 1 since there are no partial laws)
• the actual order was E, probability, E, probability

Another issue here seems to be that when the reaction is MT102, ie a projectile absorption, the LDAT(1) value was assumed to be zero in some cases.

For example: tail of the ACE files for n+H2 and p+H2. The following is the tail of p+H2:

                   1                   0                  33                  10
                   0                   2   1.00000000000E-01   1.00000000000E+00
   1.50000000000E+02   1.00000000000E+00   0.00000000000E+00   3.33382296973E-01
                   1                 102

Input files that reproduce this problem: input.txt acer.inp.txt

[whaeck]

The first issue is fixed easily. The code on lines 9982 to 9985 (and a few similar spots) is as follows:

               xss(next+2)=sigfig(xss(esz+it-1),7,0)
               xss(next+3)=1
               xss(next+4)=sigfig(xss(esz+nes-1),7,0)
               xss(next+5)=1

Changing this to:

               xss(next+2)=sigfig(xss(esz+it-1),7,0)
               xss(next+3)=sigfig(xss(esz+nes-1),7,0)
               xss(next+4)=1
               xss(next+5)=1

fixes it.

This fix is available in the fix/cp33 branch.

[whaeck]

For the second issue, NJOY seems to assume that the Q value is zero for incident neutron and charged particle files when MT=102 (see lines 10116-10117):

                     xss(next)=0
                     xss(next+1)=awi/(awr+awi)

while for other reactions it uses the following (lines 10393-10397 for LAW=3 and LAW=4):

                     amass=awr/awi
                     aprime=awp/awi
                     xss(next)=sigfig((1+amass)*(-q)/amass,7,0)
                     xss(next+1)=&
                       sigfig(amass*(amass+1-aprime)/(1+amass)**2,7,0)

[whaeck]

Released in NJOY2016.66.