7. Rotations and mirrors

probe turns its tetrahedron about an axis and reflects it across a plane, then checks what a rotation keeps and what a mirror changes: the volume of chapter 4, computed by a function of the program.

program probe
!< Tutorial, chapter 7: the tetrahedron rotated and mirrored.
use penf, only : R8P
use vecfor, only : ex, ey, ez, face_normal3, mirror_matrix, rotation_matrix, vector
implicit none
type(vector) :: o           ! the origin
type(vector) :: p(4)        ! the vertices
type(vector) :: s(4)        ! the vertices of a copy
integer      :: face(3,4)   ! the vertices of each face, counter-clockwise seen from outside
real(R8P)    :: pi          ! pi
real(R8P)    :: r(3,3)      ! a rotation matrix
real(R8P)    :: m(3,3)      ! a mirror matrix
integer      :: i           ! counter

p(1) = o ; p(2) = ex ; p(3) = ey ; p(4) = ez
face = reshape([1,3,2, 1,2,4, 1,4,3, 2,3,4], [3,4])
pi = acos(-1._R8P)
s = p
do i = 1, 4
   call s(i)%rotate(axis=ez, angle=pi / 2)     ! a quarter turn about z, counter-clockwise
enddo
call show('rotated about z', s)
r = rotation_matrix(axis=ex + ey + ez, angle=2 * pi / 3)
do i = 1, 4
   s(i) = r .matrix. p(i)                     ! a third of a turn about the diagonal
enddo
call show('rotated about (1,1,1)', s)
print '(A,F8.4)', 'volume ', volume(s)
s = p
do i = 1, 4
   call s(i)%mirror(normal=ex - ey)          ! mirrored across the plane x = y
enddo
call show('mirrored across x = y', s)
print '(A,F8.4)', 'volume ', volume(s)
m = mirror_matrix(normal=ex - ey)
s = p
do i = 1, 4
   call s(i)%mirror(matrix=m)
   call s(i)%mirror(matrix=m)                 ! twice: back where it was
enddo
call show('mirrored twice', s)

contains
   subroutine show(title, v)
   !< Print a title and the vertices, rounded to 4 decimals.
   character(*), intent(in) :: title ! What the vertices are.
   type(vector), intent(in) :: v(:)  ! The vertices.
   integer                  :: j     ! Counter.

   print '(A)', title
   do j = 1, size(v)
      print '(A,I0,3F8.4)', '  p', j, v(j)%x, v(j)%y, v(j)%z
   enddo
   endsubroutine show

   function volume(v) result(vol)
   !< The volume of the tetrahedron of vertices v, by the divergence theorem (chapter 4).
   type(vector), intent(in) :: v(4) ! The vertices.
   real(R8P)                :: vol  ! The volume.
   type(vector)             :: n    ! The normal of a face.
   integer                  :: f    ! Counter.

   vol = 0._R8P
   do f = 1, 4
      n = face_normal3(pt1=v(face(1,f)), pt2=v(face(2,f)), pt3=v(face(3,f)))
      vol = vol + (((v(face(1,f)) + v(face(2,f)) + v(face(3,f))) / 3) .dot. n) / 3
   enddo
   endfunction volume
endprogram probe

$ probe
rotated about z
  p1  0.0000  0.0000  0.0000
  p2  0.0000  1.0000  0.0000
  p3 -1.0000  0.0000  0.0000
  p4  0.0000  0.0000  1.0000
rotated about (1,1,1)
  p1  0.0000  0.0000  0.0000
  p2  0.0000  1.0000  0.0000
  p3  0.0000  0.0000  1.0000
  p4  1.0000  0.0000  0.0000
volume   0.1667
mirrored across x = y
  p1  0.0000  0.0000  0.0000
  p2  0.0000  1.0000  0.0000
  p3  1.0000  0.0000  0.0000
  p4  0.0000  0.0000  1.0000
volume  -0.1667
mirrored twice
  p1  0.0000  0.0000  0.0000
  p2  1.0000  0.0000  0.0000
  p3  0.0000  1.0000  0.0000
  p4  0.0000  0.0000  1.0000

A rotation about an axis

call v%rotate(axis=, angle=) turns v in place about the axis through the origin along axis, by angle radians, counter-clockwise when seen from the tip of the axis (the right-hand rule). The axis need not be a unit vector. A quarter turn about z takes p2 = ex to ey and p3 = ey to -ex; p4, on the axis, stays.

rotate is not elemental: it changes its vector in place, so a loop turns the vertices one by one.

A rotation matrix

rotation_matrix(axis=, angle=) returns the 3 x 3 matrix of the same rotation (Rodrigues' formula), and r .matrix. v applies a matrix to a vector, returning a new vector; call v%rotate(matrix=r) applies it in place. Compute the matrix once when many vectors turn the same way. A third of a turn about the diagonal (1, 1, 1) takes ex to ey, ey to ez and ez to ex: the tetrahedron onto itself, its vertices renamed, its volume unchanged.

A mirror

call v%mirror(normal=) reflects v across the plane through the origin orthogonal to normal, of any length. The plane x = y, of normal ex - ey, swaps x and y: p2 and p3 change places. The mirrored tetrahedron occupies the same space, but its faces, still listed in the same order, now turn clockwise seen from outside, and the volume of chapter 4 comes out negative: a mirror reverses orientation, a rotation does not.

mirror_matrix(normal=) returns the matrix of the reflection, the Householder matrix $I-2\hat{n}{\hat{n}}^T$, and call v%mirror(matrix=m) applies it in place. A reflection is its own inverse: twice, and every vertex is back.

A rotation about an axis that misses the origin, or a mirror across a plane that misses it, is a move to the origin, the transformation, and the move back (see The mirror image across a plane).

What you learned

rotate by axis and angle or by matrix, rotation_matrix, .matrix., mirror by normal or by matrix, mirror_matrix; rotations keep the orientation, mirrors reverse it. Reference: Rotations and mirrors.

Next: 8. Many vectors at once.