The vector returned by Beam.direction(x, y, z) does not depend on coordinates

[vsnever]

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 and BeamCXLine models, is rather strange:

core/cherab/core/beam/node.pyx

Lines 233 to 253 in 547e7d3

	cpdef Vector3D direction(self, double x, double y, double z):
	"""
	Calculates the beam direction vector at a point in space.
	Note the values of the beam outside of the beam envelope should be
	treated with caution.
	:param x: x coordinate in meters.
	:param y: y coordinate in meters.
	:param z: z coordinate in meters.
	:return: Direction vector.
	"""
	# 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 dx = tan(DEGREES_TO_RADIANS * self._divergence_x)
	cdef double dy = tan(DEGREES_TO_RADIANS * self._divergence_y)
	return new_vector3d(dx, dy, 1.0).normalise()

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:

        # calculate direction from divergence
        cdef double phi = atan2(y, x)
        cdef double dx = cos(phi) * tan(DEGREES_TO_RADIANS * self._divergence_x)
        cdef double dy = sin(phi) * tan(DEGREES_TO_RADIANS * self._divergence_y)
        return new_vector3d(dx, dy, 1.0).normalise()

However, I'm not an NBI expert, so I'm asking for advice here. @Mateasek, what do you think?

[Mateasek]

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.

[vsnever]

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

core/cherab/core/model/attenuator/singleray.pyx

Lines 109 to 125 in b10723a

	# 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

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?

[vsnever]

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:

[vsnever]

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)

[Mateasek]

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?

[vsnever]

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.

[vsnever]

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}$$

[vsnever]

Fixed in #433.