quaternion.f90

module quart

  implicit none
  type quaternion
     ! q(4) = (w, x, y, z) where w is the scalar part
     real :: q(4)
     ! Flags to tell us if we want to keep the quaternion unitised, and
     ! also if the quaternion represents an improper rotation (rotation + inversion)
     logical :: unit
     logical :: improper
  end type quaternion
contains

  elemental subroutine pnormalise_quaternion(q)

    ! Normalise the quaternion by dividing the quaternion by it's
    ! 'reduced norm'

    type (quaternion), intent(inout) :: q

    q%q = q%q / (pnormq(q))

  end subroutine pnormalise_quaternion

  elemental subroutine enormalise_quaternion(q)

    ! Normalise the quaternion by dividing the quaternion by it's
    ! 'reduced norm'

    type (quaternion), intent(inout) :: q

    q%q = q%q / (enormq(q))

  end subroutine enormalise_quaternion

  pure real function pnormq(q)

    ! The 'reduced norm' of a quaternion is just the square root of the
    ! arithmetic sum of the squares of it's elements

    real :: normq
    type (quaternion), intent(in) :: q

    pnormq = sqrt(sum(q%q**2))

  end function pnormq

  elemental real function enormq(q)

    ! The 'reduced norm' of a quaternion is just the square root of the
    ! arithmetic sum of the squares of it's elements

    real :: normq
    type (quaternion), intent(in) :: q

    enormq = sqrt(sum(q%q**2))

  end function enormq
end module quart

program test
    use quart

    type(quaternion) :: q(2)

    q(1)%q = [1,2,2,4]
    q(2)%q = [1,0,0,0]

    call pnormalise_quaternion(q)

    write(*,*) q(1)%q
    write(*,*) q(2)%q

    q(1)%q = [1,2,2,4]
    q(2)%q = [1,0,0,0]

    call enormalise_quaternion(q)

    write(*,*) q(1)%q
    write(*,*) q(2)%q


end program