! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> This module gathers the various predefined initial conditions
!!
module tem_spongeLayer_module

  ! include treelm modules
  use env_module,         only: rk
  use tem_param_module,   only: PI
  use tem_aux_module,     only: tem_abort
  use tem_logging_module, only: logUnit

  ! include aotus modules
  use aotus_module,       only: flu_State, aot_get_val,           &
    &                           aoterr_NonExistent, aoterr_Fatal
  use aot_table_module,   only: aot_table_open, aot_table_close,  &
   &                            aot_table_length, aot_get_val

  implicit none
  private

  public :: load_spongeLayer, spongeLayer_type, tem_spongelayer_for

  !> This type contains datas to define spongeLayer
  type spongeLayer_type
    !> Sponge Plane origin
    real(kind=rk) :: plane_origin(3)
    !> Sponge Plane normal
    real(kind=rk) :: plane_normal(3)
    !> Damp factor for the sponge Layer
    real(kind=rk) :: dampFactor
    !> damping exponent for the sponge layer
    real(kind=rk) :: dampExponent
    !> target state
    real(kind=rk), allocatable :: targetState(:)
  end type spongeLayer_type


contains

  subroutine load_spongeLayer(conf, thandle, me, ndim)
    ! ---------------------------------------------------------------------------
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> Global spongeLayer data type
    type(spongeLayer_type), intent(out) :: me
    !> number of spatial dimensions
    integer, intent(in) :: ndim
    ! ---------------------------------------------------------------------------
    integer :: cent_handle
    integer :: i
    integer :: iError
    integer, allocatable :: vError(:)
    ! ---------------------------------------------------------------------------
    ! Plane_origin
    call aot_table_open( L       = conf,                                       &
      &                  parent  = thandle,                                    &
      &                  thandle = cent_handle,                                &
      &                  key     = 'plane_origin' )
    if (aot_table_length(L=conf, thandle=cent_handle) == 3) then
      do i=1,3
        call aot_get_val( L       = conf,                                      &
          &               thandle = cent_handle,                               &
          &               pos     = i,                                         &
          &               val     = me%plane_origin(i),                        &
          &               ErrCode = iError )
      end do
    else
      write(*,*) 'ERROR while reading Plane origin of sponge Layer,'
      write(*,*) 'should have 3 entries for each coordinate!'
      STOP
    end if
    call aot_table_close(conf, cent_handle)
    write(logUnit(1),*) ' * plane_origin =', me%plane_origin



    ! Plane_normal
    call aot_table_open( L       = conf,                                       &
      &                  parent  = thandle,                                    &
      &                  thandle = cent_handle,                                &
      &                  key     = 'plane_normal' )
    if (aot_table_length(L=conf, thandle=cent_handle) == 3) then
      do i=1,3
        call aot_get_val( L       = conf,                                      &
          &               thandle = cent_handle,                               &
          &               pos     = i,                                         &
          &               val     = me%plane_normal(i),                        &
          &               ErrCode = iError )
      end do
    else
      write(*,*) 'ERROR while reading Plane origin of sponge Layer,'
      write(*,*) 'should have 3 entries for each coordinate!'
      STOP
    end if
    call aot_table_close(conf, cent_handle)
    write(logUnit(1),*) ' * plane_normal =', me%plane_normal

    !damp_factor
    call aot_get_val( L       = conf,                                          &
      &               thandle = thandle,                                       &
      &               key     = 'damp_factor',                                 &
      &               val     = me%dampFactor,                                 &
      &               ErrCode = iError )
    write(logUnit(1),*) ' * Damping_factor =', me%dampFactor

    !damp_exponent
    call aot_get_val( L       = conf,                                          &
      &               thandle = thandle,                                       &
      &               key     = 'damp_exponent',                               &
      &               val     = me%dampExponent,                               &
      &               ErrCode = iError )
    if (iError .ne. 0) then
      me%dampExponent = 1.0
    endif

    write(logUnit(1),*) ' * Damping_exponent =', me%dampExponent

    ! target_state
    call aot_table_open( L       = conf,                                       &
      &                  parent  = thandle,                                    &
      &                  thandle = cent_handle,                                &
      &                  key     = 'target_state' )

    select case (ndim)
    case (3)
      allocate(me%targetState(5))
      allocate(vError(5))
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key     = 'density',          &
           &            val     = me%targetState(1),  &
           &            ErrCode = vError(1)           )
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key     = 'velocityX',        &
           &            val     = me%targetState(2),  &
           &            ErrCode = vError(2)           )
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key    = 'velocityY',         &
           &            val     = me%targetState(3),  &
           &            ErrCode = vError(3)           )
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key     = 'velocityZ',        &
           &            val     = me%targetState(4),  &
           &            ErrCode = vError(4)           )
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key     = 'pressure',         &
           &            val     = me%targetState(5),  &
           &            ErrCode = vError(5)           )
      if (any(btest(vError, aoterr_Fatal))) then
        write(*,*) 'FATAL Error occured, when loading target state '// &
          &        'make sure to have: density,velocityX, velocityY,'//&
          &        'velocityZ,pressure! Aborting'
        call tem_abort()
      end if

    case (2)
      allocate(me%targetState(4))
      allocate(vError(4))
      call aot_get_val( L       = conf,              &
           &            thandle = cent_handle,       &
           &            key    = 'density',          &
           &            val     = me%targetState(1), &
           &            ErrCode = vError(1)          )
      call aot_get_val( L       = conf,              &
           &            thandle = cent_handle,       &
           &            key    = 'velocityX',        &
           &            val     = me%targetState(2), &
           &            ErrCode = vError(2)          )
      call aot_get_val( L       = conf,              &
           &            thandle = cent_handle,       &
           &            key    = 'velocityY',        &
           &            val     = me%targetState(3), &
           &            ErrCode = vError(3)          )
      call aot_get_val( L       = conf,              &
           &            thandle = cent_handle,       &
           &            key    = 'pressure',         &
           &            val     = me%targetState(4), &
           &            ErrCode = vError(4)          )
      if (any(btest(vError, aoterr_Fatal))) then
        write(*,*) 'FATAL Error occured, when loading target state '// &
          &        'make sure to have: density,velocityX, velocityY,'//&
          &        'pressure! Aborting'
        call tem_abort()
      end if
    case (1)
      allocate(me%targetState(3))
      allocate(vError(3))
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key     = 'density',          &
           &            val     = me%targetState(1),  &
           &            ErrCode = vError(1)            )
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key     = 'velocityX',        &
           &            val     = me%targetState(2),  &
           &            ErrCode = vError(2)           )
      call aot_get_val( L       = conf,               &
           &            thandle = cent_handle,        &
           &            key     = 'pressure',         &
           &            val     = me%targetState(3),  &
           &            ErrCode = vError(3)           )
      if (any(btest(vError, aoterr_Fatal))) then
        write(*,*) 'FATAL Error occured, when loading target state '// &
          &        'make sure to have: density,velocityX, pressure!'// &
          &        'Aborting'
        call tem_abort()
      end if
    end select

    call aot_table_close(conf, cent_handle)
    write(logUnit(1),*) ' * target_state =', me%targetState

    deallocate(vError)

  end subroutine load_spongeLayer

  !> This function calculates the sigma for the spongelayer and fills up
  !  the res with the target state
  function tem_spongelayer_for(me, nComp, coord, n)  &
    &                           result(res)
    !> Spacetime function to evaluate
    type(spongeLayer_type) :: me
    !> Number of arrays to return
    integer, intent(in) :: n
    !> Number of entrys in each array
    integer, intent(in) :: ncomp
    !> barycentric Ids of an elements.
    !! 1st index goes over number of elements and
    !! 2nd index goes over x,y,z coordinates
    real(kind=rk), intent( in ) :: coord(n,3)
    !> return value
    real(kind=rk) :: res(n,ncomp)
    ! ---------------------------------------------------------------------------
    integer :: i
    real(kind=rk) :: sigma, origin(3), normal(3), vec(3), proj_len
    ! ---------------------------------------------------------------------------
    origin(:) = me%plane_origin
    normal(:) = me%plane_normal

    do i = 1,n
      vec (:) = coord(i,:) - origin(:)
      proj_len = (vec(1)*normal(1)+ vec(2)*normal(2)+vec(3)*normal(3))/   &
                 (normal(1)**2 + normal(2)**2 + normal(3)**2)
      sigma = me%dampFactor*((proj_len)**me%dampExponent)
      if (proj_len .gt. 0) then
        res(i,1) = sigma
      else
        res(i,1) = 0.0
      endif

    enddo

    if(ncomp > 1) then
      do i = 1,n
        res(i,2:) = me%targetState(:)
      enddo
    end if

  end function tem_spongelayer_for

end module tem_spongeLayer_module
