Closed vsnever closed 3 months ago
Hi @vsnever,
yes, this looks better! The particles moving in the x-z plane should indeed have a 0 velocity component in the y direction for the thin beam approximation, since it assumes all the particles originate in a single point... Thanks for spotting this.
Hi Matej, thanks. I was trying to figure out if it was possible to calculate the bulk velocity more accurately, and I found another thing that seems strange to me. In SingleRayAttenuator.density(x, y, z)
we have:
# calculate beam width
sigma_x = self._beam.get_sigma() + z * self._tanxdiv
sigma_y = self._beam.get_sigma() + z * self._tanydiv
# normalised radius squared
norm_radius_sqr = ((x / sigma_x)**2 + (y / sigma_y)**2)
# clamp low densities to zero (beam models can skip their calculation if density is zero)
# comparison is done using the squared radius to avoid a costly square root
if self.clamp_to_zero:
if norm_radius_sqr > self._clamp_sigma_sqr:
return 0.0
# bi-variate Gaussian distribution (normalised)
gaussian_sample = exp(-0.5 * norm_radius_sqr) / (2 * M_PI * sigma_x * sigma_y)
return self._density.evaluate(z) * gaussian_sample
https://github.com/cherab/core/blob/b10723a2c03989a57f9872f0dfef87e9f33ccadc/cherab/core/model/attenuator/singleray.pyx#L109-L125
Suppose we have a point-like source of particles at (x', y', 0) emitting normally with divergence angles ɑx, ɑy, then at (x, y, z) we'll have the following distribution:
$$f_1(x, y, z) = \frac{exp\left(-\frac{1}{2}\left(\frac{(x-x')^2}{(ztg(a_x))^2}+\frac{(y-y')^2}{(ztg(a_y))^2}\right)\right)}{2\pi z^2tg(a_x)tg(a_y)}$$
But the point-like sources at z=0 plane are distributed normally as well:
$$\frac{exp\left(-\frac{x'^2 + y'^2}{2\sigma^2}\right)}{2\pi \sigma^2}$$
By averaging over the ensemble of sources at z=0, we'll get the convolution of two Gaussians:
$$f(x, y, z) = \int{-\infty}^{\infty}\int{-\infty}^{\infty}\frac{exp\left(-\frac{1}{2}\left(\frac{(x-x')^2}{(ztg(a_x))^2}+\frac{(y-y')^2}{(ztg(a_y))^2}\right)\right)}{2\pi z^2tg(a_x)tg(a_y)}\frac{exp\left(-\frac{x'^2 + y'^2}{2\sigma^2}\right)}{2\pi \sigma^2}dx'dy',$$
which is equal to:
$$f(x, y, z) = \frac{exp\left(-\frac{1}{2}\left(\frac{x^2}{\sigma_x^2}+\frac{y^2}{\sigma_y^2}\right)\right)}{2\pi \sigma_x\sigma_y},\hspace{0.5cm}\sigma_x = \sqrt{\sigma^2 + (ztg(a_x))^2}\hspace{0.5cm}\sigma_y = \sqrt{\sigma^2 + (ztg(a_y))^2}.$$
But in SingleRayAttenuator.density(x, y, z)
we have:
$$\sigma_x =\sigma + ztg(a_x)\hspace{0.5cm}\sigma_y = \sigma + ztg(a_y)$$
So is this a mistake, or am I misunderstanding the meaning of σ in the beam description?
If I am correct with the above expressions, then the direction can be calculated as follows:
# if behind the beam just return the beam axis (for want of a better value)
if z <= 0:
return self.BEAM_AXIS
# calculate direction from divergence
cdef double z_tanx_sqr = z * z * self._tanxdiv * self._tanxdiv
cdef double z_tany_sqr = z * z * self._tanydiv * self._tanydiv
cdef double sigma_sqr = self._sigma * self._sigma
cdef double sigma_x_sqr = sigma_sqr + z_tanx_sqr
cdef double sigma_y_sqr = sigma_sqr + z_tany_sqr
cdef double ex = x * z_tanx_sqr / sigma_x_sqr
cdef double ey = y * z_tany_sqr / sigma_y_sqr
return new_vector3d(ex, ey, z).normalise()
This is for σ=0.1 m and 5° divergence:
I was contacted by Péter Bálazs, who works on synthetic MSE diagnostics for ITER. He also came across this issue with the beam direction, which prevents natural broadening of the MSE line shape due to beam divergence and results in noticeable difference with the SoS code.
@Mateasek, @jacklovell, is it worth fixing this in the current beam model or we should wait until the new beam architecture is released? I think that if we fix beam direction calculation as described above, we also should change how beam's sigma is calculated for consistency. Namely, we also should change this:
sigma_x = self._beam.get_sigma() + z * self._tanxdiv
sigma_y = self._beam.get_sigma() + z * self._tanydiv
to this:
sigma_x = sqrt(self._beam.get_sigma()**2 + (z * self._tanxdiv)**2)
sigma_y = sqrt(self._beam.get_sigma()**2 + (z * self._tanxdiv)**2)
This looks good. I don't have anything against implementing it for the current beam. Also, we will have a thin beam model with the new architecture, so we will use this approximation also.
Should the beam warn the users that the beam direction definition has changed? Implementing this will change results. Or is in enough to have it in change log?
Should the beam warn the users that the beam direction definition has changed? Implementing this will change results. Or is in enough to have it in change log?
I think it's enough to add this to the changelog, because we are improving the model and the results will become more accurate. Previously, we have repeatedly improved models without generating warnings.
Should the beam warn the users that the beam direction definition has changed? Implementing this will change results. Or is in enough to have it in change log?
Also, not only will the beam direction change, which actually only affects the BES, but also the beam density, because in SingleRayAttenuator.density(x, y, z)
the
$$\sigma_x =\sigma + ztg(a_x)\hspace{0.5cm}\sigma_y = \sigma + ztg(a_y)$$
will be replaced with:
$$\sigma_x =\sqrt{\sigma^2 + (ztg(a_x))^2}\hspace{0.5cm}\sigma_y = \sqrt{\sigma^2 + (ztg(a_y))^2}$$
Fixed in #433.
I was working on updating the MSE model and found that the calculation of the direction of the beam at a point, which is used as a direction of beam bulk velocity in the
BeamEmissionLine
andBeamCXLine
models, is rather strange: https://github.com/cherab/core/blob/547e7d3b5ae979ae4bf7403c6ba1893c66b80bd3/cherab/core/beam/node.pyx#L233-L253 The vector returned by this function is independent of the given x, y coordinates and always points in the same direction, determined by the beam's divergence angles.I think that if the model takes into account beam divergence, then at least the polar angle should be taken into account:
However, I'm not an NBI expert, so I'm asking for advice here. @Mateasek, what do you think?