Hello Xilong Tong,
thank you very much for the extra files and your description. First of all, I want to re-assure you that I have been working on this problem since your message, and I am looking into any issues.
So far, I was able to check that your source routine is indeed correct. It does not interpolate between the given data points, which might provide you with smoother results, especially near cos(theta) = 0, but that would be only a minor correction for average quantities like neutron fluxes.
Please confirm that in MCNP you are also not interpolating in angle.
I found some errors in the source routine, which ultimately do not affect the results, but could be good to fix for clarity:
- when you set the names of the files to read, you are setting twice file number 9
- more serious, but inconsequential, you use the trapezoid rule for integration
probability(i,line)=(energy(i,line)-energy(i,line-1))*(intensity(i,line)+intensity(i,line-1))+probability(i,line-1)
which is missing a critical factor of 0.5 when finding the sum of the intensities. Ultimately this is not changing the results because your cdf starts from 0 and you are renormalizing at the end, but fixing it will make your code more stable to modifications.
Anyway, I have tested your sampling routine and everything appears to work (as long as what you wanted is a DT generator pointing in the y-direction)
My next step will be to see whether I can spot any reason for the enhanced fluxes at low energy in FLUKA. I will let you know if there is anything else I need,
Stefano