Description
Description
RustBCA's library provides a handful of helper functions for code coupling. One of these, rotate_given_surface_normal, rotates the direction of a particle that strikes the wall in a simulation from the simulation's coordinate system into the RustBCA coordinate system (where +x is into the surface); the particle direction is denoted u and the normal vector n.
This function uses an algorithm that constructs a rotation matrix that aligns the global into-the-surface vector (-n) with the RustBCA into-the-surface vector, (1, 0, 0). This algorithm includes the following calculation:
R = I + V V^2/(1 + c)
Where R is the rotation matrix, I is the identity, V is a matrix constructed from the elements of the cross product of the global and local (i.e., RustBCA) into-the-surface vectors, and c is the dot product of the local and global into-the-surface vectors.
There is an obvious singularity when c = -1, which is the case for n = (1, 0, 0) and u = (-1, 0, 0). I believe this is related to gimbal lock. Currently, if c == -1, the code correctly ignores the algorithm and instead calculates R as a 180-degree rotation around the y axis (chosen arbitrarily without loss of generality). However, when c is very close to -1, the last term in the above equation blows up and the algorithm produces incorrect rotation matrices that result in "initial direction out of surface" error messages.
To Reproduce
>>> import numpy as np
>>> from libRustBCA import rotate_given_surface_normal_py
>>> nx = 1 - 1e-9
>>> ny = 1 - nx*nx
>>> ux1, uy1, uz1 = rotate_given_surface_normal_py(nx, ny, 0.0, -np.sqrt(2)/2, np.sqrt(2)/2, 0)
thread '<unnamed>' panicked at src/lib.rs:1398:5:
Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. n = (0.999999999, 0.000000001999999943436137, 0) u = (-0.707106775943907, -0.7071067811865476, 0)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
pyo3_runtime.PanicException: Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. n = (0.999999999, 0.000000001999999943436137, 0) u = (-0.707106775943907, -0.7071067811865476, 0)
Expected behavior
There should be a threshold for c being too close to -1 that should revert to a 180-degree rotation around y. There is probably a way to calculate an ideal threshold, but from testing the code functions normally when nx = 1 - 1e-6 and breaks when nx = 1 - 1e-9, so the threshold should be somewhere in between those two values.