! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
! summary: module that holds datatypes and routines for coordinate transformations.
!!
module tem_coordinate_module

  ! include treelm modules
!  use atl_equation_module, only: scalar_variable, vector3_variable
  use tem_varSys_module,  only: tem_varSys_type
  use tem_logging_module, only: logUnit

  implicit none

  private

  public :: coordRotation_type
  public :: initcoordinateRotation
  public :: xToXAxes, yToXAxes, zToXAxes
  public :: dirToString

  !> constant representing the x direction
  integer, parameter :: xDir = 1
  !> constant representing the y direction
  integer, parameter :: yDir = 2
  !> constant representing the z direction
  integer, parameter :: zDir = 3

  !> Identity transformation
  integer, parameter :: xToXAxes = 1
  !> Transformation to transform a y axis into an x axis
  integer, parameter :: yToXAxes = 2
  !> Transformation to transform a z axis into an x axis
  integer, parameter :: zToXAxes = 3

  !>  summary: datatype to transform varibales given in one coordinate
  !! system to another coordinate system.
  !! This can be usefule e.g. for calculating the fluxes which is always
  !! done across the surface in positive x direction.
  type coordRotation_type
    !> The type of the rotation, see parameters above.
    integer :: rotationType
    !> Array of integers defining how to transform the variables of your
    !! state vector. Therefore the size of this array is the number of scalar
    !! variables of your equation system.
    integer, allocatable :: varTransformIndices(:)
    !> Array of integers defining how to transform the derivatives of your
    !! state vector. Therefore the size of this array is the total number
    !! of derivatives you store in your state vector.
    integer, allocatable :: derTransformIndices(:)
  end type coordRotation_type

  contains

! ****************************************************************************** !
  !>  summary: routine to specify a coordinate transformation for the state
  !! variables of your equation system.
  !!
  subroutine initCoordinateRotation( varSys, coordTrans, derivatives,          &
    &                                rotation, dimen )
    ! ---------------------------------------------------------------------------
    !> The variables to build the permutations for.
    type(tem_varSys_type), intent(in) :: varSys
    !> The rotation you want to obtain. Please have a look at the parameters
    !! of this module to find a valid input argument.
    integer, intent(in) :: coordTrans
    !> The number of derivatives (already multidimensional) you need for your
    !! equation. Zero means that we calculate cell values only, one means all
    !! first derivatives and so on.
    integer, intent(in) :: derivatives
    !> The coordinate rotation you want to initialize.
    type(coordRotation_type), intent(out) :: rotation
    !> The spatial dimension of the system
    integer, intent(in) :: dimen
    ! ---------------------------------------------------------------------------
    integer :: prevScalarVars, iVar
    ! ---------------------------------------------------------------------------

    rotation%rotationType = coordTrans

    ! get the number of scalar variables

    allocate( rotation%varTransformIndices(varSys%nScalars) )

    ! lets iterate over the variables and specify the rotation indices first.
    prevScalarVars = 0
    do iVar = 1, varSys%nStateVars
      call appendRotatedVariable( varSys%method%val(iVar)%nComponents, &
        &                         prevScalarVars, coordTrans, rotation )
      prevScalarVars = prevScalarVars + varSys%method%val(iVar)%nComponents
    end do

    ! iterate over all the derivatives and specifiy rotations for the
    ! derivatives.
    allocate( rotation%derTransformIndices(3*derivatives+1) )
    call appendDerivative( derivatives, coordTrans, rotation, dimen )

  end subroutine
! ****************************************************************************** !


! ****************************************************************************** !
  !>  summary: subroutine to append a rotated derivative to coordRotation type.
  !!
  subroutine appendDerivative( derivatives , coordTrans, rotation, dimen )
    ! ---------------------------------------------------------------------------
    !> The number of derivatives of your equation (inclunding the zeroth order
    !! derivative).
    integer, intent(in) :: derivatives
    !> The coordinate transformation you apply.
    integer, intent(in) :: coordTrans
    !> The coordinate rotation you want to initialize.
    type(coordRotation_type), intent(inout) :: rotation
    !> The spatial dimension of the system
    integer, intent(in) :: dimen
    ! ---------------------------------------------------------------------------

    select case(dimen)

    ! 3D
    case(3)
      select case(derivatives)
      case(1)
        select case(coordTrans)
        case(xToXAxes)
          rotation%derTransformIndices(1:4) = [1,2,3,4]
        case(yToXAxes)
          rotation%derTransformIndices(1:4) = [1,3,4,2]
        case(zToXAxes)
          rotation%derTransformIndices(1:4) = [1,4,2,3]
        case default
          write(logUnit(1),*)'unknown coordinate rotation for derivatives, '//   &
            &            'stopping....'
          stop
        end select
      case(0)
        select case(coordTrans)
        case(xToXAxes)
          rotation%derTransformIndices(1) = 1
        case(yToXAxes)
          rotation%derTransformIndices(1) = 1
        case(zToXAxes)
          rotation%derTransformIndices(1) = 1
        case default
          write(logUnit(1),*)'unknown coordinate rotation for derivatives, '//   &
            &            'stopping....'
          stop
        end select
      case default
        write(logUnit(1),*)'not able to create coordinate rotations for this '// &
          &            'order of derivatives, stopping...'
        stop
      end select

   ! 2D
   case(2)
      select case(derivatives)
      case(1)
        select case(coordTrans)
        case(xToXAxes)
          rotation%derTransformIndices(1:3) = [1,2,3]
        case(yToXAxes)
          rotation%derTransformIndices(1:3) = [1,3,2]
        case default
          write(logUnit(1),*)'unknown coordinate rotation for derivatives, '//   &
            &            'stopping....'
          stop
        end select
      case(0)
        select case(coordTrans)
        case(xToXAxes)
          rotation%derTransformIndices(1) = 1
        case(yToXAxes)
          rotation%derTransformIndices(1) = 1
        case default
          write(logUnit(1),*)'unknown coordinate rotation for derivatives, '//   &
            &            'stopping....'
          stop
        end select
      case default
        write(logUnit(1),*)'not able to create coordinate rotations for this '// &
          &            'order of derivatives, stopping...'
        stop
      end select
   ! 1D
   case(1)
     rotation%derTransformIndices(:) = 1
    case default
      write(logUnit(1),*)'not able to create coordinate rotations for this '// &
          &            'spatial dimension, stopping...'
      stop
    end select

  end subroutine
! ****************************************************************************** !


! ****************************************************************************** !
  !>  summary: routine to append a rotated varibale to corrdRoation type.
  !!
  subroutine appendRotatedVariable( nComponents, prevScalarVars, coordTrans,   &
    &                               rotation )
    ! ---------------------------------------------------------------------------
    !> nComponents of variable you want to append.
    integer, intent(in) :: nComponents
    !> The number of scalar variable you append before you append this variable.
    integer, intent(in) :: prevScalarVars
    !> The coordinate transformation you apply.
    integer, intent(in) :: coordTrans
    !> The coordinate rotation you want to initialize.
    type(coordRotation_type), intent(inout) :: rotation
    ! ---------------------------------------------------------------------------
    integer, allocatable :: varIndices(:)
    ! ---------------------------------------------------------------------------

    allocate( varIndices(nComponents) )

    select case(nComponents)
    case(1)
      varIndices = rotateScalar()
      rotation%varTransformIndices(prevScalarVars+1:prevScalarVars+nComponents)&
        &             = varIndices + prevScalarVars
    case(2)
      varIndices = rotateVector2(coordTrans)
      rotation%varTransformIndices(prevScalarVars+1:prevScalarVars+nComponents)&
        &             = varIndices + prevScalarVars
    case(3)
      varIndices = rotateVector3(coordTrans)
      rotation%varTransformIndices(prevScalarVars+1:prevScalarVars+nComponents)&
        &             = varIndices + prevScalarVars
    case(9)
      varIndices = rotateTensor3(coordTrans)
      rotation%varTransformIndices(prevScalarVars+1:prevScalarVars+nComponents)&
        &             = varIndices + prevScalarVars

    case default
      write(logUnit(1),*)'not able to transform this variable to a new '//     &
        &            'coordinate system, stopping....'
      stop
    end select

  end subroutine
! ****************************************************************************** !


! ****************************************************************************** !
  !>  summary: rotate a vector in 3D by a given rotation.
  !!
  function rotateVector3(coordTrans) result(rotationIndices)
    ! ---------------------------------------------------------------------------
    !> The coordinate transformation you apply.
    integer, intent(in) :: coordTrans
    !> Rotation indices for the given transformation
    integer :: rotationIndices(3)
    ! ---------------------------------------------------------------------------

    select case(coordTrans)
    case(xToXAxes)
      ! this is just the identity transformation
      rotationIndices = [1,2,3]
    case(yToXAxes)
      rotationIndices = [2,3,1]
    case(zToXAxes)
      rotationIndices = [3,1,2]
    case default
      write(logUnit(1),*)'unknown coordinate rotation, stopping....'
      stop
    end select

  end function
! ****************************************************************************** !

! ****************************************************************************** !
  !>  summary: rotate a tensor in 3D by a given rotation.
  !!
  function rotateTensor3(coordTrans) result(rotationIndices)
    ! ---------------------------------------------------------------------------
    !> The coordinate transformation you apply.
    integer, intent(in) :: coordTrans
    !> Rotation indices for the given transformation
    integer :: rotationIndices(9)
    ! ---------------------------------------------------------------------------

    select case(coordTrans)
    case(xToXAxes)
      ! this is just the identity transformation
      rotationIndices = [1,2,3,4,5,6,7,8,9]
    case(yToXAxes)
      rotationIndices = [5,6,4,8,9,7,2,3,1]
    case(zToXAxes)
      rotationIndices = [9,7,8,3,1,2,6,4,5]
    case default
      write(logUnit(1),*)'unknown coordinate rotation, stopping....'
      stop
    end select

  end function
! ****************************************************************************** !

! ****************************************************************************** !
  !>  summary: rotate a vector in 2D by a given rotation.
  !!
  function rotateVector2(coordTrans) result(rotationIndices)
    ! ---------------------------------------------------------------------------
    !> The coordinate transformation you apply.
    integer, intent(in) :: coordTrans
    !> Rotation indices for the given transformation
    integer :: rotationIndices(2)
    ! ---------------------------------------------------------------------------

    select case(coordTrans)
    case(xToXAxes)
      ! this is just the identity transformation
      rotationIndices = [1,2]
    case(yToXAxes)
      rotationIndices = [2,1]
    case(zToXAxes)
      ! rotation to z axes does not make sense.
      ! Nevertheless we define it for convenience.
      rotationIndices = [1,2]
    case default
      write(logUnit(1),*)'unknown coordinate rotation, stopping....'
      stop
    end select

  end function
! ****************************************************************************** !


! ****************************************************************************** !
  !>  summary: rotate a scalar?
  !!
  function rotateScalar() result(rotationIndices)
    ! ---------------------------------------------------------------------------
    !> The coordinate transformation you apply.
    ! integer, intent(in) :: coordTrans
    !> Rotation indices for the given transformation
    integer :: rotationIndices(1)
    ! ---------------------------------------------------------------------------
    rotationIndices = [1]

  end function rotateScalar
! ****************************************************************************** !


! ****************************************************************************** !
  !>  summary: function to convert a direction to a string.
  !!
  function dirToString(direction) result(dirAsChar)
    ! ---------------------------------------------------------------------------
    !> direction to convert
    integer, intent(in) :: direction
    !> direction as string
    character(len=1) :: dirAsChar
    ! ---------------------------------------------------------------------------
    select case(direction)
    case(xDir)
      dirAsChar = 'X'
    case(yDir)
      dirAsChar = 'Y'
    case(zDir)
      dirAsChar = 'Z'
    case default
      write(logUnit(1),*)'unknown direction for conversion to string.'
      write(logUnit(1),*)'stopping...'
      stop
    end select

  end function
! ****************************************************************************** !


end module tem_coordinate_module
! ****************************************************************************** !
