CnlPepper
commented
4 years ago Hi @Mateasek,
yes I can confirm there is a bug there. The last sample is not being counted. Fortunately for most plasma or emission scenarios this will have little effect as the last sample is generally low values and the step size small which is probably why it was being missed.
We do need a test for this, but if I remember correctly I was planning on rewriting this code so never ended up implementing one. I'm going to fix the bug now and will build a test. The actual fix should be:
for interval in range(0, intervals):
# calculate location of sample point at the top of the interval
t = (interval + 1) * step
sample_point = new_point3d(
start.x + t * integration_direction.x,
start.y + t * integration_direction.y,
start.z + t * integration_direction.z
)Nice catch.
Mateasek
I have noticed an issue while working on unittests for the Thomson scattering module for Cherab. I constructed a uniform plasma and a laser with uniform power over its volume. This is setup is to easily test in a semi-analytical way outputs of spectrum traced by a single ray, to see if everything works as expected. The issue raised when I was trying to compare intensity of the traced spectrum to the analytical one. In this case the intensity can be calculated easily by scattering from a single point times the cross section length of the ray and the laser: I(x) * length. Although length used by the integrator and me was the same, the spectrum traced by the ray had always lower intensity than the testing one.
The problem I think is in the NumericalIntegrator and more concretely in the loop taking care of the trapezoid integration.
The trapezoidal integration integrates n intervals by evaluating the function at the interval ends. This gives n+1 functional values for n intervals
As far as what I understood from the code, the variable intervals in the integrator is the number of intervals n the spectrum integration should be separated into. In the method, the 0th functional value is evaluated and then the rest is evaluated and intervals are integrated in this loop defined by
for interval in range(1, intervals):, which leads to evaluation of n functional values and thus integration of n-1 intervals.When I changed the definition of the loop to
for interval in range(1, intervals + 1):, the results from my "semi-analytical" calculations and the tracing were a good match up to a small negligible error (numerical errors and differences in physical constants probably).I think nobody has noticed this before, because there was a large number of integration steps in every material and this made the error very small. Also when I go for smaller integration steps, the error vanishes.
I think adding a test for this would be a good idea. It could be done by making a new material for testing purposes, which would allow specification of the material with function3D and an emission spectrum. Then the same test I did for the TS model could be done in Raysect.