From 2644514115434e6a9815e9884cead919ffc94599 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 11 Jun 2025 10:32:35 -0700 Subject: [PATCH 1/2] Draft of fixing numerical issues in rotate functions. --- src/lib.rs | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 3a9f53e..547661b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1502,10 +1502,10 @@ pub fn simple_compound_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1 #[cfg(feature = "parry3d")] #[no_mangle] pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut f64, uz: &mut f64) { - const DELTA: f64 = 1e-9; - let RUSTBCA_DIRECTION: Vector3:: = Vector3::::new(1.0, 0.0, 0.0); + //const DELTA: f64 = 1e-9; + //let RUSTBCA_DIRECTION: Vector3:: = Vector3::::new(1.0, 0.0, 0.0); - let into_surface = Vector3::new(-nx, -ny, -nz); + //let into_surface = Vector3::new(-nx, -ny, -nz); let direction = Vector3::new(*ux, *uy, *uz); //Rotation to local RustBCA coordinates from global @@ -1514,12 +1514,14 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu //That rotation is then applied to the particle direction, and can be undone later. //Algorithm is from here: //https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/180436#180436 - let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); - let c = into_surface.dot(&RUSTBCA_DIRECTION); - let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); + //let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); + //let c = into_surface.dot(&RUSTBCA_DIRECTION); + //let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - let rotation_matrix = if (c + 1.0).abs() > 1e-6 { - Matrix3::identity() + vx + vx*vx/(1. + c) + let s = ny*ny + nz*nz; + + let rotation_matrix = if nx > 0.0 { + Matrix3::::new(0.0, -ny, -nz, ny, nz*nz/s, -ny*nz/s, nz, -ny*nz/s, ny*ny/s) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation //around a non-x axis; y is chosen arbitrarily @@ -1531,12 +1533,12 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu // ux must not be exactly 1.0 to avoid gimbal lock in RustBCA // simple_bca does not normalize direction before proceeding, must be done manually assert!( - incident.x + DELTA > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. c={} n = ({}, {}, {}) u = ({}, {}, {}), unew = ({}, {}, {})", - (c + 1.0).abs(), nx, ny, nz, ux, uy, uz, incident.x, incident.y, incident.z + incident.x > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. n = ({}, {}, {}) u = ({}, {}, {}), unew = ({}, {}, {})", + nx, ny, nz, ux, uy, uz, incident.x, incident.y, incident.z ); - *ux = incident.x + DELTA; - *uy = incident.y - DELTA; + *ux = incident.x; + *uy = incident.y; *uz = incident.z; let mag = (ux.powf(2.) + uy.powf(2.) + uz.powf(2.)).sqrt(); From 34e31a5899199d0dc47afe1fee052491177ac99b Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 12 Jun 2025 11:21:52 -0700 Subject: [PATCH 2/2] Removed panic in rotate_given_surface_normal_py and replaced linear algebra with precomputed matrix from SymPy --- src/lib.rs | 23 ++++++----------------- 1 file changed, 6 insertions(+), 17 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 547661b..f0d7769 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1518,10 +1518,10 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu //let c = into_surface.dot(&RUSTBCA_DIRECTION); //let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - let s = ny*ny + nz*nz; + //let s = ny*ny + nz*nz; - let rotation_matrix = if nx > 0.0 { - Matrix3::::new(0.0, -ny, -nz, ny, nz*nz/s, -ny*nz/s, nz, -ny*nz/s, ny*ny/s) + let rotation_matrix = if (1.0 - nx).abs() > 0.0 { + Matrix3::::new(1. + (-ny*ny - nz*nz)/(1. - nx), -ny, -nz, ny, -ny*ny/(1. - nx) + 1., -ny*nz/(1. - nx), nz, -ny*nz/(1. - nx), -nz*nz/(1. - nx) + 1.) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation //around a non-x axis; y is chosen arbitrarily @@ -1530,13 +1530,6 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu let incident = rotation_matrix*direction; - // ux must not be exactly 1.0 to avoid gimbal lock in RustBCA - // simple_bca does not normalize direction before proceeding, must be done manually - assert!( - incident.x > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. n = ({}, {}, {}) u = ({}, {}, {}), unew = ({}, {}, {})", - nx, ny, nz, ux, uy, uz, incident.x, incident.y, incident.z - ); - *ux = incident.x; *uy = incident.y; *uz = incident.z; @@ -1617,7 +1610,7 @@ pub fn rotate_given_surface_normal_vec_py(nx: Vec, ny: Vec, nz: Vec = Vector3::::new(1.0, 0.0, 0.0); - let into_surface = Vector3::new(-nx, -ny, -nz); + //let into_surface = Vector3::new(-nx, -ny, -nz); let direction = Vector3::new(*ux, *uy, *uz); //Rotation to local RustBCA coordinates from global @@ -1625,12 +1618,8 @@ pub extern "C" fn rotate_back(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut //into-the-surface vector (1.0, 0.0, 0.0) onto the local into-the-surface vector (negative normal w.r.t. ray origin). //That rotation is then applied to the particle direction, and can be undone later. //Algorithm is from here: - //https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/180436#180436 - let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); - let c = into_surface.dot(&RUSTBCA_DIRECTION); - let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - let rotation_matrix = if c != -1.0 { - Matrix3::identity() + vx + vx*vx/(1. + c) + let rotation_matrix = if (1.0 - nx).abs() > 0.0 { + Matrix3::::new(1. + (-ny*ny - nz*nz)/(1. - nx), -ny, -nz, ny, -ny*ny/(1. - nx) + 1., -ny*nz/(1. - nx), nz, -ny*nz/(1. - nx), -nz*nz/(1. - nx) + 1.) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation //around a non-x axis; y is chosen arbitrarily