! This module is free for anyone to use and modify. If you find any bugs
! or make additions please send me a patch: evan.mcclain@gatech.edu
module vectormath
  implicit none
  type vector
     double precision :: x, y, z
  end type vector
  interface operator (*)
     module procedure scalar_vec_mult
  end interface
  interface operator (*)
     module procedure vec_scalar_mult
  end interface
  interface operator (+)
     module procedure vec_add
  end interface
  interface operator (-)
     module procedure vec_sub
  end interface
contains
  function make_vec(a, b, c)
    double precision :: a, b, c
    type (vector) :: make_vec
    make_vec%x = a
    make_vec%y = b
    make_vec%z = c
  end function make_vec
  function scalar_vec_mult(a, b)
    double precision, intent (in) :: a
    type (vector), intent (in) :: b
    type (vector) :: scalar_vec_mult
    scalar_vec_mult%x = a * b%x
    scalar_vec_mult%y = a * b%y
    scalar_vec_mult%z = a * b%z
  end function scalar_vec_mult
  function vec_scalar_mult(a, b)
    double precision, intent (in) :: b
    type (vector), intent (in) :: a
    type (vector) :: vec_scalar_mult
    vec_scalar_mult%x = b * a%x
    vec_scalar_mult%y = b * a%y
    vec_scalar_mult%z = b * a%z
  end function vec_scalar_mult
  function vec_add(a, b)
    type (vector), intent (in) :: a, b
    type (vector) vec_add
    vec_add%x = a%x + b%x
    vec_add%y = a%y + b%y
    vec_add%z = a%z + b%z
  end function vec_add
  function vec_sub(a, b)
    type (vector), intent (in) :: a, b
    type (vector) vec_sub
    vec_sub%x = a%x - b%x
    vec_sub%y = a%y - b%y
    vec_sub%z = a%z - b%z
  end function vec_sub
  function cross_prod(a, b)
    type (vector), intent (in) :: a, b
    type (vector) :: cross_prod
    cross_prod%x = a%y * b%z - a%z * b%y
    cross_prod%y = a%z * b%x - a%x * b%z
    cross_prod%z = a%x * b%y - a%y * b%x
  end function cross_prod
  function dot_prod(a, b)
    type (vector), intent (in) :: a, b
    double precision :: dot_prod
    dot_prod = a%x * b%x + a%y * b%y + a%z * b%z
  end function dot_prod
  function vec_norm(a)
    type (vector), intent (in) :: a
    double precision :: vec_norm
    vec_norm = sqrt( a%x**2 + a%y**2 + a%z**2 )
  end function vec_norm
  function unit_vec(a)
    type (vector), intent (in) :: a
    type (vector) :: unit_vec
    double precision :: inv_norm
    if (vec_norm(a) == 0.0d0) then
       unit_vec = make_vec(0.0d0, 0.0d0, 0.0d0)
       print *, "Sorry, can't make unit vec out of zero vec :("
    else
       inv_norm = 1 / vec_norm(a)
       unit_vec%x = a%x * inv_norm
       unit_vec%y = a%y * inv_norm
       unit_vec%z = a%z * inv_norm
    end if
  end function unit_vec
end module vectormath