@@ -18,7 +18,7 @@ module m_helper
1818 public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, &
1919 & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, &
2020 & f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, &
21- & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product
21+ & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul
2222
2323contains
2424
@@ -297,6 +297,48 @@ contains
297297
298298 end function f_cross
299299
300+ !> Generate a unit vector uniformly distributed on the sphere from two random parameters.
301+ function f_unit_vector(theta, eta) result(vec)
302+
303+ $:GPU_ROUTINE(parallelism=' [seq]' )
304+ real(wp), intent(in) :: theta, eta
305+ real(wp) :: zeta, xi
306+ real(wp), dimension(3) :: vec
307+
308+ xi = 2._wp*pi*theta
309+ zeta = acos(2._wp*eta - 1._wp)
310+ vec(1) = sin(zeta)*cos(xi)
311+ vec(2) = sin(zeta)*sin(xi)
312+ vec(3) = cos(zeta)
313+
314+ end function f_unit_vector
315+
316+ !> Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
317+ subroutine s_prng(var, seed)
318+
319+ $:GPU_ROUTINE(parallelism=' [seq]' )
320+ integer, intent(inout) :: seed
321+ real(wp), intent(out) :: var
322+
323+ seed = mod(modmul(seed), modulus)
324+ var = seed/real(modulus, wp)
325+
326+ end subroutine s_prng
327+
328+ !> Compute a modular multiplication step for the linear congruential pseudo-random number generator.
329+ function modmul(a) result(val)
330+
331+ $:GPU_ROUTINE(parallelism=' [seq]' )
332+ integer, intent(in) :: a
333+ integer :: val
334+ real(wp) :: x, y
335+
336+ x = (multiplier/real(modulus, wp))*a + (increment/real(modulus, wp))
337+ y = nint((x - floor(x))*decimal_trim)/decimal_trim
338+ val = nint(y*modulus)
339+
340+ end function modmul
341+
300342 !> @brief Computes the cross product c = a x b of two 3D vectors.
301343 subroutine s_cross_product(a, b, c)
302344
0 commit comments