!===============================================
! crossprod.f90
!
! prompt for two real vectors of length 3,
! calculate cross product with function
!===============================================
function dotp(x,y)
  implicit none
  real :: dotp
  real, dimension(3) :: x, y
  dotp = sum(x*y)
end function dotp

function xprod(x,y)
  implicit none
  real, dimension(3) :: xprod, x, y
  xprod(1) = x(2)*y(3) - x(3)*y(2)
  xprod(2) =-x(3)*y(1) + x(1)*y(3)
  xprod(3) = x(1)*y(2) - x(2)*y(1)
end function xprod

program crossprod
  implicit none
  integer :: i
  real :: dotp
  real, dimension(3) :: a, b, c

  interface 
     function xprod(x,y)
     implicit none
     real, dimension(3) :: xprod, x, y
     end function xprod
  end interface

  !----------------------------
  ! input vectors
  !----------------------------
  print*,'Enter first vector'
  read*, a
  print*,'Enter second vector'
  read*, b

  !----------------------------
  ! calculate dot product
  !----------------------------
  c = xprod(a,b)

  !----------------------------
  ! check magnitude
  !----------------------------
  if(abs(dotp(c,c)) < 1.0e-6) then
     print*,'---------------------------------------------------'
     print*,'Warning: magnitude of cross product is less than 1e-6'
     print*,'---------------------------------------------------'
  endif

  !----------------------------
  ! print result
  !----------------------------
  print*,'cross product = ', c

end program crossprod