From 5045e39c327e43876330d94167383f61cd813d83 Mon Sep 17 00:00:00 2001 From: anand43 Date: Fri, 26 Jan 2024 12:21:48 -0500 Subject: [PATCH] Added an extra check for rotation matrix --- src/swiftest/swiftest_obl.f90 | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/swiftest/swiftest_obl.f90 b/src/swiftest/swiftest_obl.f90 index f8fa18105..f00a95d85 100644 --- a/src/swiftest/swiftest_obl.f90 +++ b/src/swiftest/swiftest_obl.f90 @@ -55,8 +55,8 @@ module subroutine swiftest_obl_rot_matrix(n, rot, rot_matrix, rot_matrix_inv) ! Internals real(DP) :: theta !! angle to rotate it through - real(DP), dimension(3) :: u, z_hat !! unit vector about which we rotate and z_hat - real(DP), dimension(3, 3) :: S_matrix !! rotation matrices + real(DP), dimension(3) :: u, z_hat, check !! unit vector about which we rotate, z_hat, and a check variable + real(DP), dimension(3, 3) :: S_matrix, temp !! rotation matrices, and a temporary variable integer :: i, j !! dummy variable ! Assumed that NDIM = 3 @@ -100,7 +100,18 @@ module subroutine swiftest_obl_rot_matrix(n, rot, rot_matrix, rot_matrix_inv) end do rot_matrix = matinv3(rot_matrix_inv) - + + ! Check that the correct rotation matrix is used + ! rot_matrix * rot should be in the z_hat direction + check = MATMUL(rot, rot_matrix) ! 1x3 matrix x 3x3 matrix + check = .unit. check(:) + + if(abs(check(1)) .gt. EPSILON(0.0_DP) .or. abs(check(2)) .gt. EPSILON(0.0_DP)) then + temp = rot_matrix + rot_matrix = rot_matrix_inv + rot_matrix_inv = temp + end if + return end subroutine swiftest_obl_rot_matrix