! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> This module gathers the definition of the analytical solution by
!! scattering of an electromagnetic wave at a dielectric sphere/cylinder. The
!! solution is given by means of Mie-series.
!! The solution is given in:
!! Cai, W., & Deng, S. (2003). An upwinding embedded boundary method
!! for Maxwell’s equations in media with material interfaces: 2D case.
!! Journal of Computational Physics, 190(1), 159–183.
!! doi:10.1016/S0021-9991(03)00269-9
!!
module tem_miescatter_module

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

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

  implicit none

  private

  !> Parameters of the solution for Mie-Scatter at
  !! dielectric cylinder (infinite height in z direction)
  type tem_miescatter_type
    !> The center of the cylinder (in the x-y plane)
    real(kind=rk) :: center(2) = [ 0.0_rk, 0.0_rk ]

    !> The radius of the cylinder
    real(kind=rk) :: radius = 1.0_rk

    !> Permeability (mu) of the background.
    real(kind=rk) :: permeability_background = 1.0_rk

    !> Permitivity (epsilon) of the background.
    real(kind=rk) :: permitivity_background = 1.0_rk

    !> Permeability (mu) of the cylinder.
    real(kind=rk) :: permeability_cylinder = 1.0_rk

    !> Permitivity (epsilon) of the cylinder.
    real(kind=rk) :: permitivity_cylinder = 1.0_rk

    !> Wave number for the background
    real(kind=rk) :: wavenumber_background

    !> Wave number for the cylinder
    real(kind=rk) :: wavenumber_cylinder

    !> Angular frequency ( = 2 * pi * f , where
    !! f is the original frequency in Hertz)
    real(kind=rk) :: omega = 2.0_rk * PI * 1.0_rk

  end type tem_miescatter_type


  !> Expansion data for the Mier scatter solution.
  type tem_mieexpansion_type

    !> The number of expansion coefficients
    integer :: nCoeffs

    !> Expansion coefficients inside the cylinder
    complex(kind=rk), allocatable :: c_tot(:)

    !> Expansion coefficients outside the cylinder (for the scattered field)
    complex(kind=rk), allocatable :: c_scat(:)

  end type tem_mieexpansion_type


  !> Parameters of the solution for Mie-Scatter at
  !! dielectric cylinder (infinite height in z direction).
  type tem_miescatter_field_type

    !> Parameter of the geometrical and material setup
    type(tem_miescatter_type) :: miescatter

    !> The expansion parameters for the scattered solution
    type(tem_mieexpansion_type) :: mieexpansion

  end type tem_miescatter_field_type

  interface tem_load_miescatter_displacementfieldz
    module procedure tem_load_miescatter
  end interface tem_load_miescatter_displacementfieldz

  interface tem_load_miescatter_magneticfieldx
    module procedure tem_load_miescatter
  end interface tem_load_miescatter_magneticfieldx

  interface tem_load_miescatter_magneticfieldy
    module procedure tem_load_miescatter
  end interface tem_load_miescatter_magneticfieldy

  public :: tem_miescatter_field_type,              &
    &       tem_load_miescatter_displacementfieldz, &
    &       tem_load_miescatter_magneticfieldx,     &
    &       tem_load_miescatter_magneticfieldy,     &
    &       tem_eval_miescatter_displz,             &
    &       tem_eval_miescatter_magnx,              &
    &       tem_eval_miescatter_magny

contains

! ****************************************************************************** !
  !> load gauss pulse variables to set initial condition
  !!
  subroutine tem_load_miescatter(conf, thandle, me)
    ! ---------------------------------------------------------------------------
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> Global gauss pulse data type
    type(tem_miescatter_field_type), intent(out) :: me
    ! ---------------------------------------------------------------------------
    integer :: iError, iPos, cent_handle, back_handle, circ_handle
    real(kind=rk) :: frequency
    ! ---------------------------------------------------------------------------

    ! Read center cooridnates
    call aot_table_open( L       = conf,                                       &
      &                  parent  = thandle,                                    &
      &                  thandle = cent_handle,                                &
      &                  key     = 'center' )
    if ( cent_handle.eq.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*)'Not able to read center-table for Mie-Scatter '//    &
        &                'solution from lua, stopping ...'
      call tem_abort()
    end if
    do iPos = 1, 2
      call aot_get_val( L       = conf,                                        &
        &               thandle = cent_handle,                                 &
        &               pos     = iPos,                                        &
        &               val     = me%miescatter%center(iPos),                  &
        &               ErrCode = iError )
      if ( iError.ne.0 ) then
        write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
        write(logUnit(1),*)'Not able to read center for Mie-Scatter '//        &
          &                'solution from lua, stopping ...'
        call tem_abort()
      end if
    end do
    call aot_table_close(conf, cent_handle)


    ! Read radius of the cylinder
    call aot_get_val( L       = conf,                                          &
      &               thandle = thandle,                                       &
      &               key     = 'radius',                                      &
      &               val     = me%miescatter%radius,                          &
      &               ErrCode = iError )
    if ( iError.ne.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*)'Not able to read radius for Mie-Scatter '//          &
         &                'solution from lua, stopping ...'
      call tem_abort()
    end if


    ! Read background permeability and permitivity
    call aot_table_open( L       = conf,                                       &
      &                  parent  = thandle,                                    &
      &                  thandle = back_handle,                                &
      &                  key     = 'permeaPermit_background' )
    if ( back_handle.eq.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*) 'Not able to open permeaPermit_background-table for '&
                    & // 'Mie-Scatter solution from lua, stopping ...'
      call tem_abort()
    end if
    call aot_get_val( L       = conf,                                          &
      &               thandle = back_handle,                                   &
      &               pos     = 1,                                             &
      &               val     = me%miescatter%permeability_background,         &
      &               ErrCode = iError )
    if ( iError.ne.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*) 'Not able to read background permeability for ' //   &
                    & 'Mie-Scatter solution from lua, stopping ...'
      call tem_abort()
    end if
    call aot_get_val( L       = conf,                                          &
      &               thandle = back_handle,                                   &
      &               pos     = 2,                                             &
      &               val     = me%miescatter%permitivity_background,          &
      &               ErrCode = iError )
    if ( iError.ne.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*) 'Not able to read background permitivity for ' //    &
                    & 'Mie-Scatter solution from lua, stopping ...'
      call tem_abort()
    end if
    call aot_table_close(conf, back_handle)


    ! Read background permeability and permitivity
    call aot_table_open( L       = conf,                                       &
      &                  parent  = thandle,                                    &
      &                  thandle = circ_handle,                                &
      &                  key     = 'permeaPermit_cylinder' )
    if ( back_handle.eq.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*) 'Not able to open permeaPermit_cylinder-table for '//&
                    & 'Mie-Scatter solution from lua, stopping ...'
      call tem_abort()
    end if
    call aot_get_val( L       = conf,                                          &
      &               thandle = circ_handle,                                   &
      &               pos     = 1,                                             &
      &               val     = me%miescatter%permeability_cylinder,           &
      &               ErrCode = iError )
    if ( iError.ne.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*) 'Not able to read cylinder permeability for ' // &
                    & 'Mie-Scatter solution from lua, stopping ...'
      call tem_abort()
    end if
    call aot_get_val( L       = conf,                                          &
      &               thandle = circ_handle,                                   &
      &               pos     = 2,                                             &
      &               val     = me%miescatter%permitivity_cylinder,            &
      &               ErrCode = iError )
    if ( iError.ne.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*) 'Not able to read cylinder permitivity for ' // &
                    & 'Mie-Scatter solution from lua, stopping ...'
      call tem_abort()
    end if
    call aot_table_close(conf, circ_handle)


    ! Read the frequency
    call aot_get_val( L       = conf,                                          &
      &               thandle = thandle,                                       &
      &               key     = 'frequency',                                   &
      &               val     = frequency,                                     &
      &               ErrCode = iError )
    if ( iError.ne.0 ) then
      write(logUnit(1),*)'ERROR in tem_load_miescatter_displacementfieldz: '
      write(logUnit(1),*)'Not able to read frequency for Mie-Scatter '//       &
        &                'solution from lua, stopping ...'
      call tem_abort()
    end if
    me%miescatter%omega = 2.0_rk * PI * frequency

    ! Calculate the wave numbers inside and outside of the cylinder
    me%miescatter%wavenumber_cylinder = me%miescatter%omega &
    & * sqrt(me%miescatter%permeability_cylinder *          &
    &        me%miescatter%permitivity_cylinder             )
    me%miescatter%wavenumber_background = me%miescatter%omega &
    & * sqrt(me%miescatter%permeability_background *          &
    &        me%miescatter%permitivity_background             )

    ! Setup the parameters for the solution inside and outside the
    ! dielectric cylinder...
    call tem_init_data( conf, thandle, me%miescatter, me%mieexpansion )

  end subroutine tem_load_miescatter

  !> Init the expansion coefficients for the Mie-Scattering.
  subroutine tem_init_data( conf, thandle, miescatter, mieexpansion )
    ! --------------------------------------------------------------------------
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> The description of the Mie-Scattering setup.
    type(tem_miescatter_type), intent(in) :: miescatter
    !> The expansion coefficients to be filled
    type(tem_mieexpansion_type), intent(out) :: mieexpansion
    ! --------------------------------------------------------------------------
    integer :: nCoeffs, n
    ! --------------------------------------------------------------------------

    ! Get the number of expansion coefficients to be used ....
    nCoeffs = tem_get_ncoeffs_miescat(conf,thandle)
    mieexpansion%nCoeffs = nCoeffs

    ! Calculate the expansion coefficients inside the cylinder
    allocate(mieexpansion%c_tot(-nCoeffs:nCoeffs))
    innerLoop: do n = -nCoeffs, nCoeffs

      ! C_n^{tot}
      mieexpansion%c_tot(n)                                                    &
        & = (0.0_rk,1.0_rk)**(-n) *                                            &
        & ( (miescatter%wavenumber_background /                                &
        &               miescatter%permeability_background)                    &
        & * bessel_jn_derivative(n,                                            &
        &               miescatter%wavenumber_background * miescatter%radius)  &
        & * hankel2_n(n,miescatter%wavenumber_background * miescatter%radius)  &
        & - (miescatter%wavenumber_background /                                &
        &               miescatter%permeability_background)                    &
        & * hankel2_n_derivative(n,                                            &
        &               miescatter%wavenumber_background * miescatter%radius)  &
        & * bessel_jn(n,miescatter%wavenumber_background * miescatter%radius)  &
        & )                                                                    &
        & /                                                                    &
        & (                                                                    &
        &   (miescatter%wavenumber_cylinder /                                  &
        &               miescatter%permeability_cylinder)                      &
        & * bessel_jn_derivative(n,                                            &
        &               miescatter%wavenumber_cylinder * miescatter%radius)    &
        & * hankel2_n(n,miescatter%wavenumber_background * miescatter%radius)  &
        & - (miescatter%wavenumber_background /                                &
        &               miescatter%permeability_background)&
        & * hankel2_n_derivative(n,                                            &
        &               miescatter%wavenumber_background * miescatter%radius)  &
        & * bessel_jn(n,miescatter%wavenumber_cylinder * miescatter%radius)    &
        & )

    end do innerLoop

    ! Calculate the expansion coefficients outside of the cylinder, for the
    ! scattered field
    allocate(mieexpansion%c_scat(-nCoeffs:nCoeffs))
    outerLoop: do n = -nCoeffs, nCoeffs

      ! C_n^{scat}
      mieexpansion%c_scat(n) &
      & = (0.0_rk,1.0_rk)**(-n) *                                              &
      & ( (miescatter%wavenumber_background /                                  &
      &               miescatter%permeability_background)                      &
      & * bessel_jn_derivative(n,                                              &
      &               miescatter%wavenumber_background*miescatter%radius)      &
      & * bessel_jn(n,miescatter%wavenumber_cylinder*miescatter%radius)        &
      & - (miescatter%wavenumber_cylinder/miescatter%permeability_cylinder)    &
      & * bessel_jn_derivative(n,                                              &
      &               miescatter%wavenumber_cylinder*miescatter%radius)        &
      & * bessel_jn(n,miescatter%wavenumber_background*miescatter%radius)      &
      & )                                                                      &
      & /                                                                      &
      & (                                                                      &
      &   (miescatter%wavenumber_cylinder/miescatter%permeability_cylinder)    &
      & * bessel_jn_derivative(n,                                              &
      &               miescatter%wavenumber_cylinder*miescatter%radius)        &
      & * hankel2_n(n,miescatter%wavenumber_background*miescatter%radius)      &
      & - (miescatter%wavenumber_background /                                  &
      &               miescatter%permeability_background)                      &
      & * hankel2_n_derivative(n,                                              &
      &               miescatter%wavenumber_background*miescatter%radius)      &
      & * bessel_jn(n,miescatter%wavenumber_cylinder*miescatter%radius)        &
      & )

    end do outerLoop

  end subroutine tem_init_data


  !> Evaluate displacement field (z component) for Mie-Scattering of
  !! electromagnetic wave at dielectric cylinder.
  function tem_eval_miescatter_displz(me, coord, time, n) result(res)
    ! ---------------------------------------------------------------------------
    !> The function to evaluate
    type(tem_miescatter_field_type), intent(in) :: me
    !> Number of points to evaluate the function for.
    integer, intent(in) :: n
    !> 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)
    !> The time to evaluate the function at.
    real(kind=rk), intent( in )  :: time
    !> return value of the function
    real(kind=rk) :: res(n)
    ! ---------------------------------------------------------------------------
    integer :: iPoint, iCoeff
    ! Polar coordinate vector in x-y plane.
    ! First entry: radius \n
    ! Second entry: angle
    real(kind=rk) :: polar(2)
    complex(kind=rk) :: tmp
    ! ---------------------------------------------------------------------------

    do iPoint = 1, n

      ! Convert to polar coordinates (relative to the center of the
      ! cylinder.
      polar = cart2polar( coord(iPoint,1)-me%miescatter%center(1), &
                        & coord(iPoint,2)-me%miescatter%center(2)  )

      ! Inside the cylinder
      if(polar(1).le.me%miescatter%radius) then

        tmp = (0.0_rk, 0.0_rk)
        do iCoeff = -me%mieexpansion%nCoeffs , me%mieexpansion%nCoeffs
          tmp = tmp                                                           &
            & + me%mieexpansion%c_tot(iCoeff)                                 &
            & * bessel_jn(iCoeff, me%miescatter%wavenumber_cylinder*polar(1)) &
            & * exp( (0.0_rk,1.0_rk) * iCoeff * polar(2)  )
        end do
        ! Twiddle by phase and get the real part
        res(iPoint) = real( tmp * exp( (0.0_rk,1.0_rk) * me%miescatter%omega  &
          &                * time ) )

        ! Multiply with permitivity to get the displacement field
        res(iPoint) = me%miescatter%permitivity_cylinder * res(iPoint)

      ! Outside the cylinder
      else

        ! Add up incoming and scattered field
        tmp = (0.0_rk, 0.0_rk)
        do iCoeff = -me%mieexpansion%nCoeffs , me%mieexpansion%nCoeffs
          tmp = tmp +                                                     &
            & ( (0.0_rk,1.0_rk)**(-iCoeff)                                &
            &   * bessel_jn(iCoeff,                                       &
            &               me%miescatter%wavenumber_background*polar(1)) &
            &   + me%mieexpansion%c_scat(iCoeff)                          &
            &   * hankel2_n(iCoeff,                                       &
            &               me%miescatter%wavenumber_background*polar(1)) &
            & ) * exp( (0.0_rk,1.0_rk) * iCoeff * polar(2)  )
        end do
        ! Twiddle by phase and get the real part
        res(iPoint) = real(tmp * exp( (0.0_rk,1.0_rk) * &
          &                me%miescatter%omega * time ) )

        ! Multiply with permitivity to get the displacement field
        res(iPoint) = me%miescatter%permitivity_background * res(iPoint)

      end if

    end do

  end function tem_eval_miescatter_displz

  !> Evaluate magnetic field (x-component) for Mie-Scattering of
  !! electromagnetic wave at dielectric cylinder.
  function tem_eval_miescatter_magnx(me, coord, time, n) result(res)
    ! ---------------------------------------------------------------------------
    !> The function to evaluate
    type(tem_miescatter_field_type), intent(in) :: me
    !> Number of points to evaluate the function for.
    integer, intent(in) :: n
    !> 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)
    !> The time to evaluate the function at.
    real(kind=rk), intent( in )  :: time
    !> return value of the function
    real(kind=rk) :: res(n)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: H_r(n), H_theta(n), H_x, H_y, polar(2)
    integer :: iPoint
    ! ---------------------------------------------------------------------------

    ! Calculate magnetizing field in polar coordinate system
    ! ... radial component
    H_r = tem_eval_miescatter_magnradial(me, coord, time, n)
    !... angular component
    H_theta = tem_eval_miescatter_magnangular(me, coord, time, n)

    ! Convert from polar vector field to cartesian vector field
    ! and convert from magnetizing field (i.e. H) to magnetic field (i.e. B)
    do iPoint = 1, n

      ! Convert to polar coordinates (relative to the center of the
      ! cylinder.
      polar = cart2polar( coord(iPoint,1)-me%miescatter%center(1), &
                        & coord(iPoint,2)-me%miescatter%center(2)  )

      ! Vector field in cartesian coordinates
      H_x = H_r(iPoint) * cos(polar(2)) - H_theta(iPoint) * sin(polar(2))
      H_y = H_r(iPoint) * sin(polar(2)) + H_theta(iPoint) * cos(polar(2))

      ! Inside the cylinder
      if(polar(1) .le. me%miescatter%radius) then
        res(iPoint) = me%miescatter%permeability_cylinder * H_x
      ! Outside the cylinder
      else
        res(iPoint) = me%miescatter%permeability_background * H_x
      end if

    end do

  end function tem_eval_miescatter_magnx


  !> Evaluate magnetic field (y-component) for Mie-Scattering of electromagnetic
  !! wave at dielectric cylinder.
  function tem_eval_miescatter_magny(me, coord, time, n) result(res)
    ! ---------------------------------------------------------------------------
    !> The function to evaluate
    type(tem_miescatter_field_type), intent(in) :: me
    !> Number of points to evaluate the function for.
    integer, intent(in) :: n
    !> 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)
    !> The time to evaluate the function at.
    real(kind=rk), intent( in )  :: time
    !> return value of the function
    real(kind=rk) :: res(n)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: H_r(n), H_theta(n), H_x, H_y, polar(2)
    integer :: iPoint
    ! ---------------------------------------------------------------------------

    ! Calculate magnetizing field in polar coordinate system
    ! ... radial component
    H_r = tem_eval_miescatter_magnradial(me, coord, time, n)
    !... angular component
    H_theta = tem_eval_miescatter_magnangular(me, coord, time, n)

    ! Convert from polar vector field to cartesian vector field
    ! and convert from magnetizing field (i.e. H) to magnetic field (i.e. B)
    do iPoint = 1, n

      ! Convert to polar coordinates (relative to the center of the
      ! cylinder.
      polar = cart2polar( coord(iPoint,1)-me%miescatter%center(1), &
                        & coord(iPoint,2)-me%miescatter%center(2)  )

      ! Vector field in cartesian coordinates
      H_x = H_r(iPoint) * cos(polar(2)) - H_theta(iPoint) * sin(polar(2))
      H_y = H_r(iPoint) * sin(polar(2)) + H_theta(iPoint) * cos(polar(2))

      ! Inside the cylinder
      if(polar(1) .le. me%miescatter%radius) then
        res(iPoint) = me%miescatter%permeability_cylinder * H_y
      ! Outside the cylinder
      else
        res(iPoint) = me%miescatter%permeability_background * H_y
      end if

    end do

  end function tem_eval_miescatter_magny


  !> Evaluate magnetizing field (angular-component) for Mie-Scattering of
  !! electromagnetic wave at dielectric cylinder.
  function tem_eval_miescatter_magnangular(me, coord, time, n) result(res)
    ! ---------------------------------------------------------------------------
    !> The function to evaluate
    type(tem_miescatter_field_type), intent(in) :: me
    !> Number of points to evaluate the function for.
    integer, intent(in) :: n
    !> 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)
    !> The time to evaluate the function at.
    real(kind=rk), intent( in )  :: time
    !> return value of the function
    real(kind=rk) :: res(n)
    ! ---------------------------------------------------------------------------
    integer :: iPoint, iCoeff
    ! Polar coordinate vector in x-y plane.
    ! First entry: radius \n
    ! Second entry: angle
    real(kind=rk) :: polar(2)
    complex(kind=rk) :: tmp
    ! ---------------------------------------------------------------------------

    do iPoint = 1, n

      ! Convert to polar coordinates (relative to the center of the
      ! cylinder.
      polar = cart2polar( coord(iPoint,1)-me%miescatter%center(1), &
        &                 coord(iPoint,2)-me%miescatter%center(2)  )

      ! Inside the cylinder
      if(polar(1).le.me%miescatter%radius) then

        tmp = (0.0_rk, 0.0_rk)
        do iCoeff = -me%mieexpansion%nCoeffs , me%mieexpansion%nCoeffs
          tmp = tmp                                                   &
            & + me%mieexpansion%c_tot(iCoeff)                         &
            & * bessel_jn_derivative(iCoeff,                          &
            &             me%miescatter%wavenumber_cylinder*polar(1)) &
            & * exp( (0.0_rk,1.0_rk) * iCoeff * polar(2) )
        end do
        ! Twiddle by phase and get the real part
        res(iPoint) = real( tmp                                         &
          ! JZ: the paper has an additional factor of -1 here. However, I think
          ! that this factor is wrong at this place. So, I commented out
          ! the factor.
          ! & * (-1.0_rk) &
          & * (exp( (0.0_rk,1.0_rk)*me%miescatter%omega*time))          &
          & * (0.0_rk,-1.0_rk)*me%miescatter%wavenumber_cylinder        &
          & / (me%miescatter%omega*me%miescatter%permeability_cylinder) )

      ! Outside the cylinder
      else

        ! Add up incoming and scattered field
        tmp = (0.0_rk, 0.0_rk)
        do iCoeff = -me%mieexpansion%nCoeffs , me%mieexpansion%nCoeffs
          tmp = tmp +                                                     &
            & ( (0.0_rk,1.0_rk)**(-iCoeff)                                &
            &   * bessel_jn_derivative(iCoeff,                            &
            &               me%miescatter%wavenumber_background*polar(1)) &
            &   + me%mieexpansion%c_scat(iCoeff)                          &
            &   * hankel2_n_derivative(iCoeff,                            &
            &               me%miescatter%wavenumber_background*polar(1)) &
            & ) * exp( (0.0_rk,1.0_rk) * iCoeff * polar(2)  )
        end do
        ! Twiddle by phase and get the real part
        res(iPoint) = real( tmp                                            &
           ! JZ: the paper has an additional factor of -1 here. However, I think
           ! that this factor is wrong at this place. So, I commented out
           ! the factor.
           !      & * (-1.0_rk) &
           & * (exp( (0.0_rk,1.0_rk) * me%miescatter%omega * time ))       &
           & * (0.0_rk,-1.0_rk)*me%miescatter%wavenumber_background        &
           & / (me%miescatter%omega*me%miescatter%permeability_background) )

      end if

    end do

  end function tem_eval_miescatter_magnangular


  !> Evaluate magnetizing field (radial-component) for Mie-Scattering of
  !! electromagnetic wave at dielectric cylinder.
  function tem_eval_miescatter_magnradial(me, coord, time, n) result(res)
    ! ---------------------------------------------------------------------------
    !> The function to evaluate
    type(tem_miescatter_field_type), intent(in) :: me
    !> Number of points to evaluate the function for.
    integer, intent(in) :: n
    !> 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)
    !> The time to evaluate the function at.
    real(kind=rk), intent( in )  :: time
    !> return value of the function
    real(kind=rk) :: res(n)
    ! ---------------------------------------------------------------------------
    integer :: iPoint, iCoeff
    ! Polar coordinate vector in x-y plane.
    ! First entry: radius \n
    ! Second entry: angle
    real(kind=rk) :: polar(2)
    complex(kind=rk) :: tmp
    ! ---------------------------------------------------------------------------

    do iPoint = 1, n

      ! Convert to polar coordinates (relative to the center of the
      ! cylinder.
      polar = cart2polar( coord(iPoint,1)-me%miescatter%center(1), &
        &                 coord(iPoint,2)-me%miescatter%center(2)  )

      ! Inside the cylinder
      if(polar(1).le.me%miescatter%radius) then

        tmp = (0.0_rk, 0.0_rk)
        do iCoeff = -me%mieexpansion%nCoeffs , me%mieexpansion%nCoeffs
          tmp = tmp                                                           &
            & + me%mieexpansion%c_tot(iCoeff)                                 &
            & * bessel_jn(iCoeff, me%miescatter%wavenumber_cylinder*polar(1)) &
            & * (0.0_rk,1.0_rk) * iCoeff                                      &
            & * exp( (0.0_rk,1.0_rk) * iCoeff * polar(2) )
        end do
        ! Twiddle by phase and get the real part, first we check if distance from
        ! center is above a certain threshold (avoid division by zero).
        if(abs(polar(1)) > epsilon(1.0_rk) ) then
          res(iPoint) = real( tmp                                        &
            ! JZ: the paper has an additional factor of -1 here. However, I think
            ! that this factor is wrong at this place. So, I commented out
            ! the factor.
            !     & * (-1.0_rk) &
            & * (exp( (0.0_rk,1.0_rk)*me%miescatter%omega*time))         &
            & * (0.0_rk,1.0_rk)                                          &
            & / (me%miescatter%omega*me%miescatter%permeability_cylinder &
            &    *polar(1))                                              &
            & )
        else
          res(iPoint) = 0.0_rk
        end if

      ! Outside the cylinder
      else

        ! Add up incoming and scattered field
        tmp = (0.0_rk, 0.0_rk)
        do iCoeff = -me%mieexpansion%nCoeffs , me%mieexpansion%nCoeffs
          tmp = tmp +                                                     &
            & ( (0.0_rk,1.0_rk)**(-iCoeff)                                &
            &   * bessel_jn(iCoeff,                                       &
            &               me%miescatter%wavenumber_background*polar(1)) &
            &   + me%mieexpansion%c_scat(iCoeff)                          &
            &   * hankel2_n(iCoeff,                                       &
            &               me%miescatter%wavenumber_background*polar(1)) &
            & )                                                           &
            & * (0.0_rk,1.0_rk) * iCoeff                                  &
            & * exp( (0.0_rk,1.0_rk) * iCoeff * polar(2) )
        end do
        ! Twiddle by phase and get the real part, first we check if distance from
        ! center is above a certain threshold (avoid division by zero).
        if(abs(polar(1)) > epsilon(1.0_rk) ) then
          res(iPoint) = real( tmp                                          &
            ! JZ: the paper has an additional factor of -1 here. However, I think
            ! that this factor is wrong at this place. So, I commented out
            ! the factor.
            !     & * (-1.0_rk) &
            & * (exp( (0.0_rk,1.0_rk) * me%miescatter%omega * time ))      &
            & * (0.0_rk,1.0_rk)                                            &
            & / (me%miescatter%omega*me%miescatter%permeability_background &
            &    *polar(1))                                                &
            & )
        else
          res(iPoint) = 0.0_rk
        end if

      end if

    end do

  end function tem_eval_miescatter_magnradial

  !> Convert from cartesian coordinates (in the x-y plane) to
  !! polar coordinates (radius,angle)
  function cart2polar(x,y) result(polar)
    ! --------------------------------------------------------------------------
    !> X coordinate
    real(kind=rk) :: x
    !> Y coordinate
    real(kind=rk) :: y
    !> Polar coordinates, radius (first entry) and angle (second entry)
    real(kind=rk) :: polar(2)
    ! --------------------------------------------------------------------------
    ! --------------------------------------------------------------------------

    polar(1) = sqrt(x*x + y*y)

    ! Atan2 is not defined when both coordinates are zero. To cover this
    ! situation correctly, we define the angle to be 0.
    if(polar(1) > epsilon(1.0_rk) ) then
      polar(2) = atan2(y,x)
    else
      polar(2) = 0.0_rk
    end if

  end function cart2polar


  !> Compute derivative for Bessel function of first kind of order n.
  function bessel_jn_derivative(n,x) result(der)
    ! --------------------------------------------------------------------------
    !> The order of the function.
    integer :: n
    !> The evaluation point.
    real(kind=rk), intent(in) :: x
    !> The function value.
    real(kind=rk) :: der
    ! --------------------------------------------------------------------------
    ! --------------------------------------------------------------------------

    ! Calculate the derivative by recursive relation
    der = 0.5_rk*(bessel_jn(n-1,x)-bessel_jn(n+1,x))

  end function bessel_jn_derivative



  !> Compute derivative for Hankel function of second kind of order n.
  function hankel2_n_derivative(n,x) result(der)
    ! --------------------------------------------------------------------------
    !> The order of the function.
    integer :: n
    !> The evaluation point.
    real(kind=rk), intent(in) :: x
    !> The function value.
    complex(kind=rk) :: der
    ! --------------------------------------------------------------------------
    ! --------------------------------------------------------------------------

    ! Calculate the derivative by recursive relation
    der = 0.5_rk*( hankel2_n(n-1,x)-hankel2_n(n+1,x) )

  end function hankel2_n_derivative



  !> Compute Hankel function of second kind of order n.
  function hankel2_n(n,x) result(val)
    ! --------------------------------------------------------------------------
    !> The order of the function.
    integer :: n
    !> The evaluation point.
    real(kind=rk), intent(in) :: x
    !> The function value.
    complex(kind=rk) :: val
    ! --------------------------------------------------------------------------
    ! --------------------------------------------------------------------------

    val = bessel_jn(n, x) - (0.0_rk, 1.0_rk)*bessel_yn(n, x)

  end function hankel2_n


  function tem_get_ncoeffs_miescat(conf, thandle ) result(nCoeffs)
    ! --------------------------------------------------------------------------
    !> lua state type
    type(flu_State) :: conf
    !> aotus parent handle
    integer, intent(in) :: thandle
    !> The number of coefficients necessary for the Fourier series
    integer :: nCoeffs
    ! --------------------------------------------------------------------------
    integer :: iError
    ! --------------------------------------------------------------------------

    !> @todo JZ: Currently, we read the number of necessary coefficients
    !! from the lua script. However, there are also some automatic ways
    !! to estimate the number of coefficients necessary to obtain a certain
    !! numerical accuaracy, e.g. in [Wiscombe1980] .

    ! Read the frequency
    call aot_get_val( L       = conf,                                          &
      &               thandle = thandle,                                       &
      &               key     = 'nCoeffs',                                     &
      &               val     = nCoeffs,                                       &
      &               ErrCode = iError )
    if ( iError.ne.0 ) then
      call tem_abort( 'ERROR in tem_get_ncoeffs_miescat: Not able to read' &
        & // ' nCoeffs for Mie-Scatter solution from lua'                  )
    end if

  end function tem_get_ncoeffs_miescat

end module tem_miescatter_module
! ****************************************************************************** !

