!> \~japanese
!! \brief ベクトル演算や行列演算を行う基礎的な線形代数の手続きを提供するモジュール.
!! \~english
!! \brief This module provides elementary linear algebra subroutines that perform the vector and matrix operations.
!! \~
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module igmcore_linear_algebra
  
  ! モジュール引用 ; Use statements
  !

  ! 種類型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision
  
  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: vec_length, vec_normalize
  public :: dot, cross
  public :: rotateX, rotateY, rotateZ

  interface vec_normalize
!    module procedure vec_normalize1
    module procedure vec_normalize2
  end interface vec_normalize

contains

  !
  !> \~japanese
  !! \brief ベクトルを(デカルト座標系の) X 軸回りに角度のパラメータで指定された量だけ回転させる.
  !!
  !! @param[in] vec   回転行列 R を作用させるベクトル \f$\Dvect{A}=(A_x, A_y, A_z)\f$.
  !! @param[in] angle ラジアンで指定した回転角度.
  !! @return          回転後のベクトル(\f$ R \Dvect{A} \f$).
  !!
  !! \~english
  !! \brief Rotates a vector around X-axis (in Cartesian coordinate system) the amount specified by the angle parameter.
  !!
  !! @param[in] vec    A vector on acted by rotation matrix R \f$\Dvect{A}=(A_x, A_y, A_z)\f$.
  !! @param[in] angle Rotation angle specified in radians.
  !! @return          The vector after rotation (\f$ R \Dvect{A} \f$).
  !!
  function rotateX(     &
    & vec, angle        & ! (in) 
    & ) result(rel_vec)

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in) :: vec(3)   ! 回転行列を作用させるベクトル(A_x, A_y, A_z)
    real(DP), intent(in) :: angle    ! 回転させる角度
    real(DP) rel_vec(3)              ! 回転行列を作用させた後のベクトル

    ! 実行文 ; Executable statements
    !
    
    ! 3次元ベクトル vec を X 軸回りに角度 angle だけ回転させ,
    ! その結果を rel_vec に格納する.
    rel_vec(1) = vec(1)
    rel_vec(2) = cos(angle) * vec(2) - sin(angle) * vec(3)
    rel_vec(3) = sin(angle) * vec(2) + cos(angle) * vec(3)
  end function rotateX

  !
  !> \~japanese
  !! \brief ベクトルを(デカルト座標系の) Y 軸回りに角度のパラメータで指定された量だけ回転させる.
  !!
  !! @param[in] vec   回転行列 R を作用させるベクトル \f$\Dvect{A}=(A_x, A_y, A_z)\f$.
  !! @param[in] angle ラジアンで指定した回転角度.
  !! @return          回転後のベクトル(\f$ R \Dvect{A} \f$).
  !!
  !! \~english
  !! \brief Rotates a vector around Y-axis (in Cartesian coordinate system) the amount specified by the angle parameter.
  !!
  !! @param[in] vec    A vector on acted by rotation matrix R \f$\Dvect{A}=(A_x, A_y, A_z)\f$.
  !! @param[in] angle Rotation angle specified in radians.
  !! @return          The vector after rotation (\f$ R \Dvect{A} \f$).
  !!
  function rotateY(      &
    & vec, angle         & ! (in)
    & ) result(rel_vec)

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in) :: vec(3)      ! 回転行列を作用させるベクトル\f$(A_x, A_y, A_z)\f$.
    real(DP), intent(in) :: angle       ! 回転させる角度
    real(DP) rel_vec(3)                 ! 回転行列を作用させた後のベクトル

    ! 実行文 ; Executable statements
    !
    
    ! 3次元ベクトル vec を Y 軸回りに角度 angle だけ回転させ,
    ! その結果を rel_vec に格納する.
    rel_vec(1) = cos(angle) * vec(1) + sin(angle) * vec(3)
    rel_vec(2) = vec(2)
    rel_vec(3) = - sin(angle) * vec(1) + cos(angle) * vec(3)

  end function rotateY

  !
  !> \~japanese
  !! \brief ベクトルを(デカルト座標系の) Z 軸回りに角度のパラメータで指定された量だけ回転させる.
  !!
  !! @param[in] vec   回転行列 R を作用させるベクトル \f$\Dvect{A}=(A_x, A_y, A_z)\f$.
  !! @param[in] angle ラジアンで指定した回転角度.
  !! @return          回転後のベクトル(\f$ R \Dvect{A} \f$).
  !!
  !! \~english
  !! \brief Rotates a vector around Z-axis (in Cartesian coordinate system) the amount specified by the angle parameter.
  !!
  !! @param[in] vec    A vector on acted by rotation matrix R \f$\Dvect{A}=(A_x, A_y, A_z)\f$.
  !! @param[in] angle Rotation angle specified in radians.
  !! @return          The vector after rotation (\f$ R \Dvect{A} \f$).
  !!
  function rotateZ( &
    & vec, angle    & ! (in) 
    & ) result(rel_vec)

    ! 宣言文 ; Declaration statements 
    !

    real(DP), intent(in) :: vec(3)      ! 回転行列を作用させるベクトル\f$(A_x, A_y, A_z)\f$
    real(DP), intent(in) :: angle       ! 回転させる角度
    real(DP) rel_vec(3)                 ! 回転行列を作用させた後のベクトル

    ! 実行文 ; Executable statements 
    !

    ! 3次元ベクトル vec をZ 軸回りに角度 angle だけ回転させ,
    ! その結果を rel_vec に格納する.
    rel_vec(1) = cos(angle) * vec(1) - sin(angle) * vec(2)
    rel_vec(2) = sin(angle) * vec(1) + cos(angle) * vec(2)
    rel_vec(3) = vec(3)

  end function rotateZ

  !
  !> \~japanese
  !! \brief ベクトルの長さを取得する.
  !!
  !! @param[in] vec ベクトル.
  !! @return        指定されたベクトルの長さ.
  !!
  !! \~english
  !! \brief Gets the length of a vector.
  !!
  !! @param[in] vec A vector.
  !! @return        The length of the specified vector.
  !!
  function vec_length( &
    & vec              &   ! (in)
    & ) result(length)

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in) :: vec(3)  ! 長さを取得したいベクトル
    real(DP) length                 ! ベクトルの長さ

    ! 実行文 ; Executable statements
    !

    length = sqrt(vec(1)**2 + vec(2)**2 + vec(3)**2)
    !length = sqrt(dot_product(vec(:), vec(:)))

  end function vec_length

  !
  !> \~japanese
  !! \brief 指定されたベクトルを規格化する.
  !!
  !! @param[in,out]  vec The vector that will normarize.
  !!
  !! \~english
  !! \brief Normalizes the specified vector.
  !!
  !! @param[in,out]  vec The vector that will normarize.
  !!
  subroutine vec_normalize1( &
    & vec                   &   ! (inout)
    & )

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(inout) :: vec(3)    ! 規格化したいベクトル

    ! 実行文 ; Executable statements
    !

    vec(1:3) = vec(1:3) / vec_length(vec)

  end subroutine vec_normalize1

  function vec_normalize2( &
    & vec                   &   ! (inout)
    & ) result(normalVec)

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in) :: vec(3)    ! 規格化したいベクトル
    real(DP) :: normalvec(3)

    ! 実行文 ; Executable statements
    !

    normalvec(1:3) = vec(1:3) / vec_length(vec)

  end function vec_normalize2

  !
  !> \~japanese
  !! \brief ベクトル\f$\Dvect{v}_1, \Dvect{v}_2\f$の外積 \f$\Dvect{v}_1 \times \Dvect{v}_2\f$ を計算する.
  !!
  !! @param[in]  v1 ベクトル\f$\Dvect{v}_1\f$.
  !! @param[in]  v2 ベクトル\f$\Dvect{v}_2\f$.
  !! @return   The vector in which is stored the result of this cross product.
  !!
  !! \~english
  !! \brief Calculates the cross product of vectors \f$\Dvect{v}_1\f$ and \f$ \Dvect{v}_2\f$.
  !!
  !! @param[in]  v1 Vector \f$\Dvect{v}_1\f$.
  !! @param[in]  v2 Vector \f$\Dvect{v}_2\f$.
  !! @return     The vector in which is stored the result of this cross product.
  !!
  function cross( &
    & v1, v2    &     ! (in)
    & ) result(v)

    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in) :: v1(3)   ! ベクトル\f$\Dvect{v}_1\f$. 
    real(DP), intent(in) :: v2(3)   ! ベクトル\f$\Dvect{v}_2\f$. 
    real(DP) v(3)                   ! \f$\Dvect{v}_1 \times \Dvect{v}_2\f$の計算結果を格納したベクトル

    ! 実行文 ; Executable statements
    !

    v(1) = v1(2)*v2(3) - v1(3)*v2(2)
    v(2) = v1(3)*v2(1) - v1(1)*v2(3)
    v(3) = v1(1)*v2(2) - v1(2)*v2(1)
  end function cross

  !
  !> \~japanese
  !! \brief ベクトル\f$\Dvect{v}_1, \Dvect{v}_2\f$の内積 \f$\Dvect{v}_1 \cdot \Dvect{v}_2\f$ を計算する.
  !!
  !! @param[in]  v1 ベクトル\f$\Dvect{v}_1\f$.
  !! @param[in]  v2 ベクトル\f$\Dvect{v}_2\f$.
  !! @return        内積の計算結果を格納したベクトル.
  !!
  !! \~english
  !! \brief Calculates the dot product of vectors \f$\Dvect{v}_1\f$ and \f$ \Dvect{v}_2\f$.
  !!
  !! @param[in]  v1 Vector \f$\Dvect{v}_1\f$.
  !! @param[in]  v2 Vector \f$\Dvect{v}_2\f$.
  !! @return        The vector in which is stored the result of this dot product.
  !!
  function dot( &
    & v1, v2    &    ! (in)
    &) result(val)


    ! 宣言文 ; Declaration statements
    !
    real(DP), intent(in) :: v1(3)  ! ベクトル\f$\Dvect{v}_1\f$
    real(DP), intent(in) :: v2(3)  ! ベクトル\f$\Dvect{v}_2\f$
    real(DP) val                   ! \f$\Dvect{v}_1 \cdot \Dvect{v}_2\f$の計算結果

    ! 実行文 ; Executable statements
    !

    val = v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3)

  end function dot

end module igmcore_linear_algebra
