propagator_mxll.F90

Go to the documentation of this file.

!! Copyright (C) 2019 R. Jestaedt, F. Bonafe, H. Appel
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
 
#include "global.h"
 
module propagator_mxll_oct_m
  use accel_oct_m
  use batch_ops_oct_m
  use batch_oct_m
  use comm_oct_m
  use cube_function_oct_m
  use cube_oct_m
  use derivatives_oct_m
  use density_oct_m
  use energy_calc_oct_m
  use energy_mxll_oct_m
  use exponential_oct_m
  use external_densities_oct_m
  use external_waves_oct_m
  use fft_oct_m
  use fourier_space_oct_m
  use grid_oct_m
  use global_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_mxll_oct_m
  use helmholtz_decomposition_m
  use index_oct_m
  use io_oct_m
  use io_function_oct_m
  use lalg_adv_oct_m
  use lalg_basic_oct_m
  use loct_math_oct_m
  use math_oct_m
  use maxwell_boundary_op_oct_m
  use maxwell_function_oct_m
  use mesh_oct_m
  use mesh_batch_oct_m
  use mesh_cube_parallel_map_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mpi_debug_oct_m
  use output_oct_m
  use parser_oct_m
  use par_vec_oct_m
  use profiling_oct_m
  use poisson_oct_m
  use quantity_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_mxll_oct_m
  use string_oct_m
  use tdfunction_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                                  &
    propagator_mxll_t,                       &
    propagator_mxll_init,                    &
    mxll_propagation_step,                   &
    mxll_propagate_leapfrog,                 &
    mxll_propagate_expgauss1,                &
    mxll_propagate_expgauss2,                &
    transform_rs_densities,                  &
    energy_mxll_calc,                        &
    energy_mxll_calc_batch,                  &
    spatial_constant_calculation,            &
    set_medium_rs_state,                     &
    plane_waves_in_box_calculation,          &
    plane_waves_boundaries_calculation,      &
    mirror_pmc_boundaries_calculation,       &
    mirror_pec_boundaries_calculation,       &
    mask_absorbing_boundaries,               &
    constant_boundaries_calculation,         &
    mxll_apply_boundaries
 
  ! The following routines are currently unused, but will be used in the near future.
  ! In order not to generate warnings about them, we declared them as public
  public ::                                  &
    calculate_matter_longitudinal_field,     &
    get_vector_pot_and_transverse_field,     &
    calculate_vector_potential
 
  type propagator_mxll_t
    logical             :: bc_add_ab_region  = .false.
    logical             :: bc_zero           = .false.
    logical             :: bc_constant       = .false.
    logical             :: bc_mirror_pec     = .false.
    logical             :: bc_mirror_pmc     = .false.
    logical             :: bc_plane_waves    = .false.
    logical             :: bc_medium         = .false.
    type(exponential_t) :: te
    logical             :: plane_waves_in_box
    integer             :: tr_etrs_approx
  end type propagator_mxll_t
 
  integer, public, parameter ::   &
    RS_TRANS_FORWARD  = 1,       &
    rs_trans_backward = 2
 
  integer, parameter ::    &
    MXWLL_ETRS_FULL  = 0,  &
    mxwll_etrs_const = 1
 
contains
 
  ! ---------------------------------------------------------
  subroutine propagator_mxll_init(gr, namespace, st, hm, tr)
    type(grid_t),                 intent(in)    :: gr
    type(namespace_t),            intent(in)    :: namespace
    type(states_mxll_t),          intent(inout) :: st
    type(hamiltonian_mxll_t),     intent(inout) :: hm
    type(propagator_mxll_t),      intent(inout) :: tr
 
    integer :: idim
 
    push_sub(propagator_mxll_init)
 
    call profiling_in("PROPAGATOR_MXLL_INIT")
 
    do idim = 1, 3
      select case (hm%bc%bc_type(idim))
      case (mxll_bc_zero)
        hm%bc_zero = .true.
        tr%bc_zero = .true.
      case (mxll_bc_constant)
        tr%bc_constant = .true.
        tr%bc_add_ab_region = .true.
        hm%bc_constant = .true.
        hm%bc_add_ab_region = .true.
      case (mxll_bc_mirror_pec)
        tr%bc_mirror_pec = .true.
        hm%bc_mirror_pec = .true.
      case (mxll_bc_mirror_pmc)
        tr%bc_mirror_pmc = .true.
        hm%bc_mirror_pmc = .true.
      case (mxll_bc_plane_waves)
        tr%bc_plane_waves = .true.
        tr%bc_add_ab_region = .true.
        hm%plane_waves = .true.
        hm%bc_plane_waves = .true.
        hm%bc_add_ab_region = .true.
      end select
    end do
 
    if (any(hm%bc%bc_type(1:3) == mxll_bc_constant)) then
      call td_function_mxll_init(st, namespace, hm)
      safe_allocate(st%rs_state_const(1:st%dim))
      st%rs_state_const = m_z0
    end if
 
    !%Variable MaxwellTDETRSApprox
    !%Type integer
    !%Default no
    !%Section Maxwell::TD Propagation
    !%Description
    !% Whether to perform  approximations to the ETRS propagator.
    !%Option no 0
    !% No approximations.
    !%Option const_steps 1
    !% Use constant current density.
    !%End
    call parse_variable(namespace, 'MaxwellTDETRSApprox', mxwll_etrs_full, tr%tr_etrs_approx)
 
    !%Variable MaxwellPlaneWavesInBox
    !%Type logical
    !%Default no
    !%Section Maxwell
    !%Description
    !% Analytic evaluation of the incoming waves inside the box,
    !% not doing any numerical propagation of Maxwells equations.
    !%End
    call parse_variable(namespace, 'MaxwellPlaneWavesInBox', .false., tr%plane_waves_in_box)
    if (tr%plane_waves_in_box .and. .not. hm%bc%do_plane_waves) then
      call external_waves_init(hm%bc%plane_wave, namespace)
    end if
 
    call messages_print_with_emphasis(namespace=namespace)
 
    !tr%te%exp = .true.
    call exponential_init(tr%te, namespace, full_batch=.true.) ! initialize Maxwell propagator
 
    call profiling_out("PROPAGATOR_MXLL_INIT")
 
    pop_sub(propagator_mxll_init)
  end subroutine propagator_mxll_init
 
  ! ---------------------------------------------------------
  subroutine mxll_propagation_step(hm, namespace, gr, space, st, tr, rs_stateb, ff_rs_inhom_t1, ff_rs_inhom_t2, time, dt)
    type(hamiltonian_mxll_t),   intent(inout) :: hm
    type(namespace_t),          intent(in)    :: namespace
    type(grid_t),               intent(inout) :: gr
    class(space_t),             intent(in)    :: space
    type(states_mxll_t),        intent(inout) :: st
    type(propagator_mxll_t),    intent(inout) :: tr
    type(batch_t),              intent(inout) :: rs_stateb
    complex(real64), contiguous,          intent(in)    :: ff_rs_inhom_t1(:,:)
    complex(real64), contiguous,          intent(in)    :: ff_rs_inhom_t2(:,:)
    real(real64),               intent(in)    :: time
    real(real64),               intent(in)    :: dt
 
    integer            :: ii, ff_dim, idim, istate, inter_steps
    real(real64)       :: inter_dt, inter_time
 
 
    logical            :: pml_check
    type(batch_t) :: ff_rs_stateb, ff_rs_state_pmlb
    type(batch_t) :: ff_rs_inhom_1b, ff_rs_inhom_2b, ff_rs_inhom_meanb
    complex(real64), allocatable :: rs_state(:, :)
 
    push_sub(mxll_propagation_step)
 
    call profiling_in('MXLL_PROPAGATOR_STEP')
    pml_check = .false.
 
    if (hm%ma_mx_coupling_apply) then
      message(1) = "Maxwell-matter coupling not implemented yet"
      call messages_fatal(1, namespace=namespace)
    end if
    safe_allocate(rs_state(gr%np, st%dim))
 
    if (tr%plane_waves_in_box) then
      rs_state = m_z0
      call plane_waves_in_box_calculation(hm%bc, time+dt, gr, gr%der, st, rs_state)
      call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
      safe_deallocate_a(rs_state)
      pop_sub(mxll_propagation_step)
      return
    end if
 
    do idim = 1, 3
      if (hm%bc%bc_ab_type(idim) == option__maxwellabsorbingboundaries__cpml) then
        pml_check = .true.
      end if
    end do
 
    ! this must be called only once, but can not be placed in init routines because PML parameters need to know about dt
    if (pml_check .and. .not. hm%bc%pml%parameters_initialized) &
      call bc_mxll_generate_pml_parameters(hm%bc, space, gr, hm%c_factor, dt)
 
    ff_dim = hm%dim
 
    ! intermediate step variables
    inter_steps   = 1
    inter_dt      = m_one / inter_steps * dt
 
    call zbatch_init(ff_rs_stateb, 1, 1, hm%dim, gr%np_part)
    if (st%pack_states) call ff_rs_stateb%do_pack(copy=.false.)
 
    if (pml_check) then
      call ff_rs_stateb%copy_to(ff_rs_state_pmlb)
    end if
 
    ! first step of Maxwell inhomogeneity propagation with constant current density
    if ((hm%ma_mx_coupling_apply .or. hm%current_density_ext_flag .or. hm%current_density_from_medium) .and. &
      tr%tr_etrs_approx == mxwll_etrs_const) then
      call ff_rs_stateb%copy_to(ff_rs_inhom_1b)
      call ff_rs_stateb%copy_to(ff_rs_inhom_2b)
      call ff_rs_stateb%copy_to(ff_rs_inhom_meanb)
 
      do istate = 1, hm%dim
        call batch_set_state(ff_rs_inhom_meanb, istate, gr%np, ff_rs_inhom_t1(:, istate))
        call batch_set_state(ff_rs_inhom_2b, istate, gr%np, ff_rs_inhom_t2(:, istate))
      end do
      call batch_axpy(gr%np, m_one, ff_rs_inhom_2b, ff_rs_inhom_meanb)
      call batch_scal(gr%np, m_half, ff_rs_inhom_meanb)
 
      ! inhomogeneity propagation
      call ff_rs_inhom_meanb%copy_data_to(gr%np, ff_rs_inhom_1b)
      call ff_rs_inhom_meanb%copy_data_to(gr%np, ff_rs_inhom_2b)
 
      call hamiltonian_mxll_update(hm, time=time)
      hm%cpml_hamiltonian = .false.
      call tr%te%apply_batch(namespace, gr, hm, ff_rs_inhom_2b, inter_dt)
 
      ! add term U(time+dt,time)J(time)
      call batch_axpy(gr%np, m_one, ff_rs_inhom_2b, ff_rs_inhom_1b)
      call ff_rs_inhom_meanb%copy_data_to(gr%np, ff_rs_inhom_2b)
      call hamiltonian_mxll_update(hm, time=time)
      call tr%te%apply_batch(namespace, gr, hm, ff_rs_inhom_2b, inter_dt*m_half)
      ! add term U(time+dt/2,time)J(time)
      call batch_axpy(gr%np, m_one, ff_rs_inhom_2b, ff_rs_inhom_1b)
      call ff_rs_inhom_meanb%copy_data_to(gr%np, ff_rs_inhom_2b)
      call hamiltonian_mxll_update(hm, time=time)
      call tr%te%apply_batch(namespace, gr, hm, ff_rs_inhom_2b, -inter_dt*m_half)
      ! add term U(time,time+dt/2)J(time)
      call batch_axpy(gr%np, m_one, ff_rs_inhom_2b, ff_rs_inhom_1b)
      call ff_rs_inhom_2b%end()
      call ff_rs_inhom_meanb%end()
    end if
 
    do ii = 1, inter_steps
 
      ! intermediate time
      inter_time = time + inter_dt * (ii-1)
 
      ! transformation of RS state into 3x3 or 4x4 representation
      call transform_rs_state_batch(hm, gr, st, rs_stateb, ff_rs_stateb, rs_trans_forward)
 
      ! RS state propagation
      call hamiltonian_mxll_update(hm, time=inter_time)
      if (pml_check) then
        call pml_propagation_stage_1_batch(hm, gr, st, tr, ff_rs_stateb, ff_rs_state_pmlb)
      end if
 
      hm%cpml_hamiltonian = pml_check
      call tr%te%apply_batch(namespace, gr, hm, ff_rs_stateb, dt)
      hm%cpml_hamiltonian = .false.
 
      if (pml_check) then
        call pml_propagation_stage_2_batch(hm, namespace, gr, st, tr, inter_time, inter_dt, m_zero, ff_rs_state_pmlb, ff_rs_stateb)
      end if
 
      !Below we add the contribution from the inhomogeneous terms
      if ((hm%ma_mx_coupling_apply) .or. hm%current_density_ext_flag .or. hm%current_density_from_medium) then
        if (tr%tr_etrs_approx == mxwll_etrs_full) then
          call ff_rs_stateb%copy_to(ff_rs_inhom_1b)
          call ff_rs_stateb%copy_to(ff_rs_inhom_2b)
          call ff_rs_stateb%copy_to(ff_rs_inhom_meanb)
 
          ! Interpolation of the external current
          do istate = 1, hm%dim
            call batch_set_state(ff_rs_inhom_meanb, istate, gr%np, ff_rs_inhom_t2(:, istate))
            call batch_set_state(ff_rs_inhom_1b, istate, gr%np, ff_rs_inhom_t1(:, istate))
          end do
          ! store t1 - t2 for the interpolation in mean
          call batch_axpy(gr%np, -m_one, ff_rs_inhom_1b, ff_rs_inhom_meanb)
          call ff_rs_inhom_1b%copy_data_to(gr%np, ff_rs_inhom_2b)
          call batch_axpy(gr%np, ii / real(inter_steps, real64) , ff_rs_inhom_meanb, ff_rs_inhom_2b)
          call batch_axpy(gr%np, (ii-1) / real(inter_steps, real64) , ff_rs_inhom_meanb, ff_rs_inhom_1b)
 
          hm%cpml_hamiltonian = .false.
          call tr%te%apply_batch(namespace, gr, hm, ff_rs_inhom_1b, inter_dt)
          ! add terms U(time+dt,time)J(time) and J(time+dt)
          call batch_axpy(gr%np, -m_fourth * inter_dt, ff_rs_inhom_1b, ff_rs_stateb)
          call batch_axpy(gr%np, -m_fourth * inter_dt, ff_rs_inhom_2b, ff_rs_stateb)
 
          do istate = 1, hm%dim
            call batch_set_state(ff_rs_inhom_1b, istate, gr%np, ff_rs_inhom_t1(:, istate))
            call batch_set_state(ff_rs_inhom_2b, istate, gr%np, ff_rs_inhom_t2(:, istate))
          end do
          call batch_axpy(gr%np, m_one, ff_rs_inhom_2b, ff_rs_inhom_1b)
          call batch_scal(gr%np, m_half, ff_rs_inhom_1b)
          call ff_rs_inhom_1b%copy_data_to(gr%np, ff_rs_inhom_2b)
 
          call tr%te%apply_batch(namespace, gr, hm, ff_rs_inhom_1b, inter_dt/m_two)
          call tr%te%apply_batch(namespace, gr, hm, ff_rs_inhom_2b, -inter_dt/m_two)
 
          ! add terms U(time+dt/2,time)J(time) and U(time,time+dt/2)J(time+dt)
          call batch_axpy(gr%np, -m_fourth * inter_dt, ff_rs_inhom_1b, ff_rs_stateb)
          call batch_axpy(gr%np, -m_fourth * inter_dt, ff_rs_inhom_2b, ff_rs_stateb)
 
          call ff_rs_inhom_1b%end()
          call ff_rs_inhom_2b%end()
          call ff_rs_inhom_meanb%end()
        else if (tr%tr_etrs_approx == mxwll_etrs_const) then
          call batch_axpy(gr%np, -m_fourth * inter_dt, ff_rs_inhom_1b, ff_rs_stateb)
        end if
      end if
 
      ! PML convolution function update
      if (pml_check) then
        call cpml_conv_function_update(hm, gr, ff_rs_state_pmlb)
      end if
 
      ! back transformation of RS state representation
      call transform_rs_state_batch(hm, gr, st, rs_stateb, ff_rs_stateb, rs_trans_backward)
 
      if (tr%bc_constant) then
        call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
        ! Propagation dt with H(inter_time+inter_dt) for constant boundaries
        if (st%rs_state_const_external) then
          call spatial_constant_calculation(tr%bc_constant, st, gr, hm, inter_time, inter_dt, m_zero, rs_state)
        end if
        call constant_boundaries_calculation(tr%bc_constant, hm%bc, hm, st, rs_state)
        call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
      end if
 
      ! Propagation dt with H(inter_time+inter_dt) for PEC mirror boundaries
      if (any(hm%bc%bc_type == mxll_bc_mirror_pec)) then
        call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
        call mirror_pec_boundaries_calculation(hm%bc, st, rs_state)
        call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
      end if
 
      ! Propagation dt with H(inter_time+inter_dt) for PMC mirror boundaries
      if (any(hm%bc%bc_type == mxll_bc_mirror_pmc)) then
        call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
        call mirror_pmc_boundaries_calculation(hm%bc, st, rs_state)
        call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
      end if
 
      ! Apply mask absorbing boundaries
      if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__mask)) then
        call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
        call mask_absorbing_boundaries(namespace, gr, hm, st, tr, inter_time, inter_dt, m_zero, rs_state)
        call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
      end if
 
      if (tr%bc_plane_waves) then
        ! Propagation dt with H(inter_time+inter_dt) for plane waves boundaries
        call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
        call plane_waves_boundaries_calculation(hm, st, gr, inter_time+inter_dt, m_zero, rs_state)
        call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
      end if
 
    end do
 
    if (tr%tr_etrs_approx == option__maxwelltdetrsapprox__const_steps) then
      call ff_rs_inhom_1b%end()
    end if
 
    call ff_rs_stateb%end()
 
    if (pml_check) then
      call ff_rs_state_pmlb%end()
    end if
 
    safe_deallocate_a(rs_state)
 
    call profiling_out('MXLL_PROPAGATOR_STEP')
 
    pop_sub(mxll_propagation_step)
  end subroutine mxll_propagation_step
 
  ! ---------------------------------------------------------
  subroutine mxll_propagate_leapfrog(hm, namespace, gr, st, tr, time, dt, counter)
    type(hamiltonian_mxll_t),   intent(inout) :: hm
    type(namespace_t),          intent(in)    :: namespace
    type(grid_t),               intent(inout) :: gr
    type(states_mxll_t),        intent(inout) :: st
    type(propagator_mxll_t),    intent(inout) :: tr
    real(real64),               intent(in)    :: time
    real(real64),               intent(in)    :: dt
    integer,                    intent(in)    :: counter
 
    type(batch_t) :: rs_state_tmpb
 
    push_sub_with_profile(mxll_propagate_leapfrog)
 
    call st%rs_stateb%copy_to(rs_state_tmpb)
 
    call hamiltonian_mxll_update(hm, time)
 
    ! do boundaries at the beginning
    call mxll_apply_boundaries(tr, st, hm, gr, namespace, time, dt, st%rs_stateb)
 
    ! update PML convolution values
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      call mxll_update_pml_simple(hm, st%rs_stateb)
    end if
 
    ! apply hamiltonian
    call hamiltonian_mxll_apply_simple(hm, namespace, gr, st%rs_stateb, rs_state_tmpb)
    call batch_scal(gr%np, -m_zi, rs_state_tmpb)
 
    ! add inhomogeneous terms
    call batch_xpay(gr%np, st%inhomogeneousb, m_one, rs_state_tmpb)
 
    if (counter == 0) then
      ! for the first step, we do one forward Euler step
      call batch_xpay(gr%np, st%rs_stateb, dt, rs_state_tmpb)
    else
      ! the leapfrog step depends on the previous state
      call batch_xpay(gr%np, st%rs_state_prevb, m_two*dt, rs_state_tmpb)
    end if
 
    ! save the current rs state
    call st%rs_stateb%copy_data_to(gr%np, st%rs_state_prevb)
    ! update to new timestep
    call rs_state_tmpb%copy_data_to(gr%np, st%rs_stateb)
 
    ! update PML convolution values
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      call mxll_copy_pml_simple(hm, st%rs_stateb)
    end if
 
    call rs_state_tmpb%end()
 
    pop_sub_with_profile(mxll_propagate_leapfrog)
  end subroutine mxll_propagate_leapfrog
 
  ! ---------------------------------------------------------
  subroutine mxll_propagate_expgauss1(hm, namespace, gr, st, tr, time, dt)
    type(hamiltonian_mxll_t),   intent(inout) :: hm
    type(namespace_t),          intent(in)    :: namespace
    type(grid_t),               intent(inout) :: gr
    type(states_mxll_t),        intent(inout) :: st
    type(propagator_mxll_t),    intent(inout) :: tr
    real(real64),               intent(in)    :: time
    real(real64),               intent(in)    :: dt
 
    type(batch_t) :: rs_state_tmpb
 
    push_sub_with_profile(mxll_propagate_expgauss1)
 
    call st%rs_stateb%copy_to(rs_state_tmpb)
 
    call hamiltonian_mxll_update(hm, time)
 
    ! do boundaries at the beginning (should be included in Hamiltonian?)
    call mxll_apply_boundaries(tr, st, hm, gr, namespace, time, &
      dt, st%rs_stateb)
 
    ! update PML convolution values
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      call mxll_update_pml_simple(hm, st%rs_stateb)
    end if
 
    ! accumulate -i H F_n - J_1 in rs_state_tmpb
    ! compute H F_n
    call hm%zapply(namespace, gr, st%rs_stateb, rs_state_tmpb)
    ! compute -i H F_n
    call batch_scal(gr%np, -m_zi, rs_state_tmpb)
    if (hm%current_density_ext_flag .or. hm%current_density_from_medium) then
      ! set J_1
      call mxll_set_batch(st%inhomogeneousb, st%rs_current_density_t1, gr%np, st%dim)
      ! accumulate -J_1 to rs_state_tmpb
      call batch_axpy(gr%np, -m_one, st%inhomogeneousb, rs_state_tmpb)
    end if
    ! compute phi_1
    call tr%te%apply_phi_batch(namespace, gr, hm, rs_state_tmpb, dt, 1)
    ! F_{n+1} = F_n + dt * phi_1 (...)
    call batch_axpy(gr%np, dt, rs_state_tmpb, st%rs_stateb)
 
    call rs_state_tmpb%end()
 
    ! update PML convolution values
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      call mxll_copy_pml_simple(hm, st%rs_stateb)
    end if
 
    pop_sub_with_profile(mxll_propagate_expgauss1)
  end subroutine mxll_propagate_expgauss1
 
  ! ---------------------------------------------------------
  subroutine mxll_propagate_expgauss2(hm, namespace, gr, st, tr, time, dt)
    type(hamiltonian_mxll_t),   intent(inout) :: hm
    type(namespace_t),          intent(in)    :: namespace
    type(grid_t),               intent(inout) :: gr
    type(states_mxll_t),        intent(inout) :: st
    type(propagator_mxll_t),    intent(inout) :: tr
    real(real64),               intent(in)    :: time
    real(real64),               intent(in)    :: dt
 
    type(batch_t) :: rs_state_tmpb
 
    push_sub_with_profile(mxll_propagate_expgauss2)
 
    call st%rs_stateb%copy_to(rs_state_tmpb)
 
    call hamiltonian_mxll_update(hm, time)
 
    ! do boundaries at the beginning (should be included in Hamiltonian?)
    call mxll_apply_boundaries(tr, st, hm, gr, namespace, time, &
      dt, st%rs_stateb)
 
    ! update PML convolution values
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      call mxll_update_pml_simple(hm, st%rs_stateb)
    end if
 
    ! accumulate -i H F_n - a_1 J_1 - a_2 J_2 in rs_state_tmpb
    ! compute H F_n
    call hm%zapply(namespace, gr, st%rs_stateb, rs_state_tmpb)
    ! compute -i H F_n
    call batch_scal(gr%np, -m_zi, rs_state_tmpb)
    if (hm%current_density_ext_flag .or. hm%current_density_from_medium) then
      ! set J_1
      call mxll_set_batch(st%inhomogeneousb, st%rs_current_density_t1, gr%np, st%dim)
      ! accumulate -a_1 J_1 to rs_state_tmpb
      call batch_axpy(gr%np, -m_half*(m_one+sqrt(m_three)), st%inhomogeneousb, rs_state_tmpb)
      ! set J_2
      call mxll_set_batch(st%inhomogeneousb, st%rs_current_density_t2, gr%np, st%dim)
      ! accumulate -a_2 J_2 to rs_state_tmpb
      call batch_axpy(gr%np, -m_half*(m_one-sqrt(m_three)), st%inhomogeneousb, rs_state_tmpb)
    end if
    ! compute phi_1
    call tr%te%apply_phi_batch(namespace, gr, hm, rs_state_tmpb, dt, 1)
    ! accumulate phi_1 term: F_{n+1} = F_n + dt * phi_1 (...)
    call batch_axpy(gr%np, dt, rs_state_tmpb, st%rs_stateb)
    if (hm%current_density_ext_flag .or. hm%current_density_from_medium) then
      call batch_set_zero(rs_state_tmpb)
      ! set J_1
      call mxll_set_batch(st%inhomogeneousb, st%rs_current_density_t1, gr%np, st%dim)
      ! accumulate -b_1 J_1 to rs_state_tmpb
      call batch_axpy(gr%np, sqrt(m_three), st%inhomogeneousb, rs_state_tmpb)
      ! set J_2
      call mxll_set_batch(st%inhomogeneousb, st%rs_current_density_t2, gr%np, st%dim)
      ! accumulate -b_2 J_2 to rs_state_tmpb
      call batch_axpy(gr%np, -sqrt(m_three), st%inhomogeneousb, rs_state_tmpb)
 
      ! compute phi_2
      call tr%te%apply_phi_batch(namespace, gr, hm, rs_state_tmpb, dt, 2)
      ! accumulate phi_2 term: F_{n+1} = F_n + dt * phi_2 (...)
      call batch_axpy(gr%np, dt, rs_state_tmpb, st%rs_stateb)
    end if
 
    ! update PML convolution values
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
      call mxll_copy_pml_simple(hm, st%rs_stateb)
    end if
 
    call rs_state_tmpb%end()
 
    pop_sub_with_profile(mxll_propagate_expgauss2)
  end subroutine mxll_propagate_expgauss2
 
  ! ---------------------------------------------------------
  subroutine set_medium_rs_state(st, gr, hm)
    type(states_mxll_t),      intent(inout) :: st
    type(grid_t),             intent(in)    :: gr
    type(hamiltonian_mxll_t), intent(in)    :: hm
 
    integer :: ip, ip_in, il, idim
 
    push_sub(set_medium_rs_state)
 
    assert(allocated(st%ep) .and. allocated(st%mu))
 
    call profiling_in('SET_MEDIUM_RS_STATE')
 
    if (hm%calc_medium_box) then
      do il = 1, size(hm%medium_boxes)
        assert(.not. hm%medium_boxes(il)%has_mapping)
        do ip = 1, hm%medium_boxes(il)%points_number
          if (abs(hm%medium_boxes(il)%c(ip)) <= m_epsilon) cycle
          st%ep(ip) = hm%medium_boxes(il)%ep(ip)
          st%mu(ip) = hm%medium_boxes(il)%mu(ip)
        end do
      end do
    end if
 
    do idim = 1, st%dim
      if (hm%bc%bc_type(idim) == mxll_bc_medium) then
        do ip_in = 1, hm%bc%medium(idim)%points_number
          ip = hm%bc%medium(idim)%points_map(ip_in)
          st%ep(ip) = hm%bc%medium(idim)%ep(ip_in)
          st%mu(ip) = hm%bc%medium(idim)%mu(ip_in)
        end do
      end if
    end do
 
    call profiling_out('SET_MEDIUM_RS_STATE')
 
    pop_sub(set_medium_rs_state)
  end subroutine set_medium_rs_state
 
  ! ---------------------------------------------------------
  subroutine transform_rs_state_batch(hm, gr, st, rs_stateb, ff_rs_stateb, sign)
    type(hamiltonian_mxll_t), intent(in)         :: hm
    type(grid_t),             intent(in)         :: gr
    type(states_mxll_t),      intent(in)         :: st
    type(batch_t),            intent(inout)      :: rs_stateb
    type(batch_t),            intent(inout)      :: ff_rs_stateb
    integer,                  intent(in)         :: sign
 
    complex(real64), allocatable :: rs_state(:,:)
    complex(real64), allocatable :: rs_state_tmp(:,:)
    integer :: ii, np
 
    push_sub(transform_rs_state_batch)
 
    call profiling_in('TRANSFORM_RS_STATE')
 
    assert(sign == rs_trans_forward .or. sign == rs_trans_backward)
 
    np = gr%np
    safe_allocate(rs_state(1:gr%np, 1:st%dim))
 
    if (hm%operator == faraday_ampere_medium) then
      if (sign == rs_trans_forward) then
        call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
        ! 3 to 6
        do ii = 1, 3
          call batch_set_state(ff_rs_stateb, ii, np, rs_state(:, ii))
          call batch_set_state(ff_rs_stateb, ii+3, np, conjg(rs_state(:, ii)))
        end do
      else
        ! 6 to 3
        safe_allocate(rs_state_tmp(1:gr%np, 1:st%dim))
        do ii = 1, 3
          call batch_get_state(ff_rs_stateb, ii, np, rs_state(:, ii))
          call batch_get_state(ff_rs_stateb, ii+3, np, rs_state_tmp(:, ii))
          rs_state(1:np, ii) = m_half * (rs_state(1:np, ii) + conjg(rs_state_tmp(1:np, ii)))
        end do
        call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
        safe_deallocate_a(rs_state_tmp)
      end if
    else
      if (sign == rs_trans_forward) then
        call rs_stateb%copy_data_to(gr%np, ff_rs_stateb)
      else
        call ff_rs_stateb%copy_data_to(gr%np, rs_stateb)
      end if
    end if
    safe_deallocate_a(rs_state)
 
    call profiling_out('TRANSFORM_RS_STATE')
 
    pop_sub(transform_rs_state_batch)
 
  end subroutine transform_rs_state_batch
 
  ! ---------------------------------------------------------
  subroutine transform_rs_densities(hm, mesh, rs_current_density, ff_density, sign)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    class(mesh_t),            intent(in)    :: mesh
    complex(real64),          intent(inout) :: rs_current_density(:,:)
    complex(real64),          intent(inout) :: ff_density(:,:)
    integer,                  intent(in)    :: sign
 
 
    assert(sign == rs_trans_forward .or. sign == rs_trans_backward)
    assert(size(rs_current_density, dim=2) == 3)
 
    push_sub(transform_rs_densities)
 
    call profiling_in('TRANSFORM_RS_DENSITIES')
 
    if (hm%operator == faraday_ampere_medium) then
      if (sign == rs_trans_forward) then
        call transform_rs_densities_to_6x6_rs_densities_forward(mesh, rs_current_density, ff_density)
      else
        call transform_rs_densities_to_6x6_rs_densities_backward(mesh, ff_density, rs_current_density)
      end if
    else
      if (sign == rs_trans_forward) then
        ff_density(1:mesh%np, 1:3) = rs_current_density(1:mesh%np, 1:3)
      else
        rs_current_density(1:mesh%np, 1:3) = ff_density(1:mesh%np, 1:3)
      end if
    end if
 
    call profiling_out('TRANSFORM_RS_DENSITIES')
 
    pop_sub(transform_rs_densities)
 
  end subroutine transform_rs_densities
 
  !----------------------------------------------------------
  subroutine transform_rs_densities_to_6x6_rs_densities_forward(mesh, rs_current_density, rs_density_6x6)
    class(mesh_t),            intent(in)    :: mesh
    complex(real64),          intent(in)    :: rs_current_density(:,:)
    complex(real64),          intent(inout) :: rs_density_6x6(:,:)
 
    integer :: ii
 
    assert(size(rs_current_density, dim=2) == 3)
    assert(size(rs_density_6x6, dim=2) == 6)
 
    ! no push_sub, called to frequently
    do ii = 1, 3
      rs_density_6x6(1:mesh%np, ii) = rs_current_density(1:mesh%np, ii)
      rs_density_6x6(1:mesh%np, ii+3) = rs_current_density(1:mesh%np, ii)
    end do
 
  end subroutine transform_rs_densities_to_6x6_rs_densities_forward
 
  !----------------------------------------------------------
  subroutine transform_rs_densities_to_6x6_rs_densities_backward(mesh, rs_density_6x6, rs_current_density)
    class(mesh_t),            intent(in)    :: mesh
    complex(real64),          intent(in)    :: rs_density_6x6(:,:)
    complex(real64),          intent(inout) :: rs_current_density(:,:)
 
    integer :: ii
 
    assert(size(rs_current_density, dim=2) == 3)
    assert(size(rs_density_6x6, dim=2) == 6)
 
    ! no push_sub, called to frequently
    do ii = 1, 3
      rs_current_density(1:mesh%np, ii) = m_half * &
        real(rs_density_6x6(1:mesh%np, ii) + rs_density_6x6(1:mesh%np, ii+3), real64)
    end do
 
  end subroutine transform_rs_densities_to_6x6_rs_densities_backward
 
  !----------------------------------------------------------
  subroutine calculate_matter_longitudinal_field(gr_mxll, st_mxll, hm_mxll, gr_elec, hm_elec, rs_state_matter)
    type(grid_t),                  intent(in)    :: gr_mxll
    type(states_mxll_t),           intent(in)    :: st_mxll
    type(hamiltonian_mxll_t),      intent(in)    :: hm_mxll
    type(grid_t),                  intent(in)    :: gr_elec
    type(hamiltonian_elec_t),      intent(in)    :: hm_elec
    complex(real64),               intent(inout) :: rs_state_matter(:,:)
 
    complex(real64), allocatable :: tmp_pot_mx_gr(:,:), tmp_grad_mx_gr(:,:)
 
    safe_allocate(tmp_pot_mx_gr(1:gr_mxll%np_part,1))
    safe_allocate(tmp_grad_mx_gr(1:gr_mxll%np,1:gr_mxll%box%dim))
    ! this subroutine needs the matter part
 
    push_sub(calculate_matter_longitudinal_field)
 
    tmp_pot_mx_gr(:,:) = m_zero
    tmp_grad_mx_gr(:,:) = m_zero
    call zderivatives_grad(gr_mxll%der, tmp_pot_mx_gr(:,1), tmp_grad_mx_gr(:,:), set_bc = .false.)
    tmp_grad_mx_gr = - tmp_grad_mx_gr
 
    rs_state_matter = m_z0
    call build_rs_state(real(tmp_grad_mx_gr(1:gr_mxll%np,:)), aimag(tmp_grad_mx_gr(1:gr_mxll%np,:)), st_mxll%rs_sign, &
      rs_state_matter(1:gr_mxll%np,:), gr_mxll, st_mxll%ep(1:gr_mxll%np), st_mxll%mu(1:gr_mxll%np), &
      gr_mxll%np)
 
    safe_deallocate_a(tmp_pot_mx_gr)
    safe_deallocate_a(tmp_grad_mx_gr)
 
    pop_sub(calculate_matter_longitudinal_field)
  end subroutine calculate_matter_longitudinal_field
 
  !----------------------------------------------------------
  subroutine get_vector_pot_and_transverse_field(namespace, gr_mxll, hm_mxll, st_mxll, tr_mxll, hm, &
    poisson_solver, helmholtz, field, transverse_field, vector_potential)
    type(namespace_t),               intent(in)    :: namespace
    type(grid_t),                    intent(in)    :: gr_mxll
    type(hamiltonian_mxll_t),        intent(in)    :: hm_mxll
    type(states_mxll_t),             intent(in)    :: st_mxll
    type(propagator_mxll_t),         intent(in)    :: tr_mxll
    type(hamiltonian_elec_t),        intent(in)    :: hm
    type(poisson_t),                 intent(in)    :: poisson_solver
    type(helmholtz_decomposition_t), intent(inout)    :: helmholtz
    complex(real64),                 intent(inout) :: field(:,:)
    complex(real64),                 intent(inout) :: transverse_field(:,:)
    real(real64),                    intent(inout) :: vector_potential(:,:)
 
    integer            :: np
    complex(real64), allocatable :: rs_state_plane_waves(:, :)
 
    push_sub(get_vector_pot_and_transverse_field)
 
    transverse_field = m_z0
    vector_potential = m_zero
 
    np = gr_mxll%np
 
    if (hm_mxll%ma_mx_coupling) then
 
      ! check what other transverse field methods are needed
 
      ! trans_calc_method == OPTION__MAXWELLTRANSFIELDCALCULATIONMETHOD__TRANS_FIELD_POISSON
      if (tr_mxll%bc_plane_waves .and. hm_mxll%plane_waves_apply) then
        safe_allocate(rs_state_plane_waves(1:gr_mxll%np, 1:st_mxll%dim))
        call mxll_get_batch(st_mxll%rs_state_plane_wavesb, rs_state_plane_waves, gr_mxll%np, st_mxll%dim)
      end if
 
      ! plane waves subtraction
      if (tr_mxll%bc_plane_waves .and. hm_mxll%plane_waves_apply) then
        transverse_field(1:np,:) = field(1:np,:) - rs_state_plane_waves(1:np,:)
      else
        transverse_field(1:np,:) = field(1:np,:)
      end if
      ! apply helmholtz decomposition for transverse field
      call helmholtz%get_trans_field(namespace, transverse_field, total_field=field)
      ! plane waves addition
      if (tr_mxll%bc_plane_waves .and. hm_mxll%plane_waves_apply) then
        transverse_field(1:np,:) = transverse_field(1:np,:) + rs_state_plane_waves(1:np,:)
        safe_deallocate_a(rs_state_plane_waves)
      end if
 
    else
 
      transverse_field(1:np,:) = field
 
    end if
 
 
    pop_sub(get_vector_pot_and_transverse_field)
 
  end subroutine  get_vector_pot_and_transverse_field
 
  ! ---------------------------------------------------------
  subroutine calculate_vector_potential(namespace, poisson_solver, gr, st, field, vector_potential)
    type(namespace_t),          intent(in)    :: namespace
    type(poisson_t),            intent(in)    :: poisson_solver
    type(grid_t),               intent(in)    :: gr
    type(states_mxll_t),        intent(in)    :: st
    complex(real64),            intent(in)    :: field(:,:)
    real(real64), contiguous,   intent(inout) :: vector_potential(:,:)
 
    integer :: idim
    real(real64), allocatable :: dtmp(:,:)
 
    safe_allocate(dtmp(1:gr%np_part,1:3))
 
    dtmp = m_zero
 
    call get_magnetic_field_state(field, gr, st%rs_sign, vector_potential, st%mu, gr%np_part)
    dtmp = vector_potential
    call dderivatives_curl(gr%der, dtmp, vector_potential, set_bc = .false.)
    do idim=1, st%dim
      call dpoisson_solve(poisson_solver, namespace, dtmp(:,idim), vector_potential(:,idim), .true.)
    end do
    vector_potential = m_one / (m_four * m_pi) * vector_potential
 
    safe_deallocate_a(dtmp)
 
  end subroutine calculate_vector_potential
 
  !----------------------------------------------------------
  subroutine energy_mxll_calc(gr, st, hm, energy_mxll, rs_field, rs_field_plane_waves)
    type(grid_t),             intent(in)  :: gr
    type(states_mxll_t),      intent(in)  :: st
    type(hamiltonian_mxll_t), intent(in)  :: hm
    type(energy_mxll_t),      intent(inout) :: energy_mxll
    complex(real64),          intent(in)  :: rs_field(:,:)
    complex(real64), optional,          intent(in)  :: rs_field_plane_waves(:,:)
 
    real(real64), allocatable :: energy_density(:), e_energy_density(:), b_energy_density(:), energy_density_plane_waves(:)
 
    push_sub(energy_mxll_calc)
 
    call profiling_in('ENERGY_MXLL_CALC')
 
    safe_allocate(energy_density(1:gr%np))
    safe_allocate(e_energy_density(1:gr%np))
    safe_allocate(b_energy_density(1:gr%np))
    if (present(rs_field_plane_waves) .and. hm%plane_waves) then
      safe_allocate(energy_density_plane_waves(1:gr%np))
    end if
 
    call energy_density_calc(gr, st, rs_field, energy_density, e_energy_density, &
      b_energy_density, hm%plane_waves, rs_field_plane_waves, energy_density_plane_waves)
    energy_mxll%energy    = dmf_integrate(gr, energy_density, mask=st%inner_points_mask)
    energy_mxll%e_energy  = dmf_integrate(gr, e_energy_density, mask=st%inner_points_mask)
    energy_mxll%b_energy  = dmf_integrate(gr, b_energy_density, mask=st%inner_points_mask)
    if (present(rs_field_plane_waves) .and. hm%plane_waves) then
      energy_mxll%energy_plane_waves = dmf_integrate(gr, energy_density_plane_waves, mask=st%inner_points_mask)
    else
      energy_mxll%energy_plane_waves = m_zero
    end if
 
    energy_mxll%boundaries = dmf_integrate(gr, energy_density, mask=st%boundary_points_mask)
 
    safe_deallocate_a(energy_density)
    safe_deallocate_a(e_energy_density)
    safe_deallocate_a(b_energy_density)
    if (present(rs_field_plane_waves) .and. hm%plane_waves) then
      safe_deallocate_a(energy_density_plane_waves)
    end if
 
    call profiling_out('ENERGY_MXLL_CALC')
 
    pop_sub(energy_mxll_calc)
  end subroutine energy_mxll_calc
 
  !----------------------------------------------------------
  subroutine energy_mxll_calc_batch(gr, st, hm, energy_mxll, rs_fieldb, rs_field_plane_wavesb)
    type(grid_t),             intent(in)    :: gr
    type(states_mxll_t),      intent(in)    :: st
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(energy_mxll_t),      intent(inout) :: energy_mxll
    type(batch_t),            intent(in)    :: rs_fieldb
    type(batch_t),            intent(in)    :: rs_field_plane_wavesb
 
    type(batch_t) :: e_fieldb, b_fieldb, e_field_innerb, b_field_innerb, rs_field_plane_waves_innerb
    real(real64) :: tmp(1:st%dim)
    complex(real64) :: ztmp(1:st%dim)
 
    push_sub(energy_mxll_calc_batch)
 
    call profiling_in('ENERGY_MXLL_CALC_BATCH')
 
    call dbatch_init(e_fieldb, 1, 1, st%dim, gr%np)
    if (st%pack_states) then
      call e_fieldb%do_pack(copy=.false.)
    end if
    call e_fieldb%copy_to(b_fieldb)
    call e_fieldb%copy_to(e_field_innerb)
    call e_fieldb%copy_to(b_field_innerb)
 
    call batch_split_complex(gr%np, rs_fieldb, e_fieldb, b_fieldb)
 
    ! subtract energy of inner points
    call batch_set_zero(e_field_innerb)
    call batch_set_zero(b_field_innerb)
    if (accel_is_enabled()) then
      call batch_copy_with_map(st%inner_points_number, st%buff_inner_points_map, e_fieldb, e_field_innerb)
      call batch_copy_with_map(st%inner_points_number, st%buff_inner_points_map, b_fieldb, b_field_innerb)
    else
      call batch_copy_with_map(st%inner_points_number, st%inner_points_map, e_fieldb, e_field_innerb)
      call batch_copy_with_map(st%inner_points_number, st%inner_points_map, b_fieldb, b_field_innerb)
    end if
    call dmesh_batch_dotp_vector(gr, e_field_innerb, e_field_innerb, tmp)
    energy_mxll%e_energy = sum(tmp)
    call dmesh_batch_dotp_vector(gr, b_field_innerb, b_field_innerb, tmp)
    energy_mxll%b_energy = sum(tmp)
    energy_mxll%energy = energy_mxll%e_energy + energy_mxll%b_energy
 
    call dmesh_batch_dotp_vector(gr, e_fieldb, e_fieldb, tmp)
    energy_mxll%boundaries = sum(tmp)
    call dmesh_batch_dotp_vector(gr, b_fieldb, b_fieldb, tmp)
    energy_mxll%boundaries = energy_mxll%boundaries + sum(tmp)
    energy_mxll%boundaries = energy_mxll%boundaries - energy_mxll%energy
 
    if (hm%plane_waves) then
      call rs_field_plane_wavesb%copy_to(rs_field_plane_waves_innerb)
      call batch_set_zero(rs_field_plane_waves_innerb)
      if (accel_is_enabled()) then
        call batch_copy_with_map(st%inner_points_number, st%buff_inner_points_map, &
          rs_field_plane_wavesb, rs_field_plane_waves_innerb)
      else
        call batch_copy_with_map(st%inner_points_number, st%inner_points_map, &
          rs_field_plane_wavesb, rs_field_plane_waves_innerb)
      end if
      call zmesh_batch_dotp_vector(gr, rs_field_plane_waves_innerb, rs_field_plane_waves_innerb, ztmp)
      energy_mxll%energy_plane_waves = sum(real(tmp, real64) )
      call rs_field_plane_waves_innerb%end()
    else
      energy_mxll%energy_plane_waves = m_zero
    end if
 
    call e_fieldb%end()
    call b_fieldb%end()
    call e_field_innerb%end()
    call b_field_innerb%end()
 
    call profiling_out('ENERGY_MXLL_CALC_BATCH')
 
    pop_sub(energy_mxll_calc_batch)
  end subroutine energy_mxll_calc_batch
 
  ! ---------------------------------------------------------
  subroutine mask_absorbing_boundaries(namespace, gr, hm, st, tr, time, dt, time_delay, rs_state)
    type(namespace_t),          intent(in)    :: namespace
    type(grid_t),               intent(in)    :: gr
    type(hamiltonian_mxll_t),   intent(inout) :: hm
    type(states_mxll_t),        intent(inout) :: st
    type(propagator_mxll_t),    intent(inout) :: tr
    real(real64),               intent(in)    :: time
    real(real64),               intent(in)    :: dt
    real(real64),               intent(in)    :: time_delay
    complex(real64),            intent(inout) :: rs_state(:,:)
 
    integer            :: ip, ip_in, idim
    logical            :: mask_check
 
    push_sub(mask_absorbing_boundaries)
 
    call profiling_in('MASK_ABSORBING_BOUNDARIES')
    mask_check = .false.
 
    do idim = 1, 3
      if (hm%bc%bc_ab_type(idim) == option__maxwellabsorbingboundaries__mask) then
        mask_check = .true.
      end if
    end do
 
    if (mask_check) then
      if (tr%bc_plane_waves .and. hm%plane_waves_apply) then
        call plane_waves_propagation(hm, tr, namespace, st, gr, time, dt, time_delay)
        call mxll_get_batch(st%rs_state_plane_wavesb, st%rs_state_plane_waves, gr%np, st%dim)
        rs_state = rs_state - st%rs_state_plane_waves
        call maxwell_mask(hm, rs_state)
        rs_state = rs_state + st%rs_state_plane_waves
      else if (tr%bc_constant .and. hm%spatial_constant_apply) then
        !call constant_at_absorbing_boundaries_calculation(st, hm%bc)
        call constant_boundaries_calculation(tr%bc_constant, hm%bc, hm, st, rs_state)
        do ip_in=1, hm%bc%constant_points_number
          ip = hm%bc%constant_points_map(ip_in)
          rs_state(ip,:) = rs_state(ip,:) - st%rs_state_const(:)
        end do
        call maxwell_mask(hm, rs_state)
        do ip_in=1, hm%bc%constant_points_number
          ip = hm%bc%constant_points_map(ip_in)
          rs_state(ip,:) = rs_state(ip,:) + st%rs_state_const(:)
        end do
      else
        call maxwell_mask(hm, rs_state)
      end if
    end if
 
    call profiling_out('MASK_ABSORBING_BOUNDARIES')
 
    pop_sub(mask_absorbing_boundaries)
  end subroutine mask_absorbing_boundaries
 
  ! ---------------------------------------------------------
  subroutine maxwell_mask(hm, rs_state)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    complex(real64),          intent(inout) :: rs_state(:,:)
 
    integer :: ip, ip_in, idim
 
    push_sub(maxwell_mask)
 
    call profiling_in('MAXWELL_MASK')
 
    do idim = 1, 3
      if (hm%bc%bc_ab_type(idim) == option__maxwellabsorbingboundaries__mask) then
        do ip_in = 1, hm%bc%mask_points_number(idim)
          ip = hm%bc%mask_points_map(ip_in,idim)
          rs_state(ip,:) = rs_state(ip,:) * hm%bc%mask(ip_in,idim)
        end do
      end if
    end do
 
    call profiling_out('MAXWELL_MASK')
 
    pop_sub(maxwell_mask)
  end subroutine maxwell_mask
 
  ! ---------------------------------------------------------
  subroutine pml_propagation_stage_1_batch(hm, gr, st, tr, ff_rs_stateb, ff_rs_state_pmlb)
    type(hamiltonian_mxll_t),   intent(inout) :: hm
    type(grid_t),               intent(in)    :: gr
    type(states_mxll_t),        intent(inout) :: st
    type(propagator_mxll_t),    intent(inout) :: tr
    type(batch_t),              intent(in)    :: ff_rs_stateb
    type(batch_t),              intent(inout) :: ff_rs_state_pmlb
 
    integer            :: ii
    complex(real64), allocatable :: rs_state_constant(:,:)
    type(batch_t) :: rs_state_constantb
 
    push_sub(pml_propagation_stage_1_batch)
 
    call profiling_in('PML_PROP_STAGE_1_BATCH')
 
    if (tr%bc_plane_waves .and. hm%plane_waves_apply) then
      call transform_rs_state_batch(hm, gr, st, st%rs_state_plane_wavesb, &
        ff_rs_state_pmlb, rs_trans_forward)
      call batch_xpay(gr%np, ff_rs_stateb, -m_one, ff_rs_state_pmlb)
    else if (tr%bc_constant .and. hm%spatial_constant_apply) then
      ! this could be optimized: right now we broadcast the constant value
      !  to the full mesh to be able to use the batch functions easily.
      ! in principle, we would need to do the transform only for one point
      !  and then subtract that value from all points of the state
      safe_allocate(rs_state_constant(1:gr%np,1:3))
      do ii = 1, 3
        rs_state_constant(1:gr%np, ii) = st%rs_state_const(ii)
      end do
      call ff_rs_stateb%copy_to(rs_state_constantb)
      call mxll_set_batch(rs_state_constantb, rs_state_constant, gr%np, 3)
 
      call transform_rs_state_batch(hm, gr, st, rs_state_constantb, &
        ff_rs_state_pmlb, rs_trans_forward)
      call batch_xpay(gr%np, ff_rs_stateb, -m_one, ff_rs_state_pmlb)
 
      call rs_state_constantb%end()
 
      safe_deallocate_a(rs_state_constant)
    else
      ! this copy should not be needed
      call ff_rs_stateb%copy_data_to(gr%np, ff_rs_state_pmlb)
    end if
 
    call profiling_out('PML_PROP_STAGE_1_BATCH')
 
    pop_sub(pml_propagation_stage_1_batch)
  end subroutine pml_propagation_stage_1_batch
 
  ! ---------------------------------------------------------
  subroutine pml_propagation_stage_2_batch(hm, namespace, gr, st, tr, time, dt, time_delay, ff_rs_state_pmlb, ff_rs_stateb)
    type(hamiltonian_mxll_t),   intent(inout) :: hm
    type(namespace_t),          intent(in)    :: namespace
    type(grid_t),               intent(in)    :: gr
    type(states_mxll_t),        intent(inout) :: st
    type(propagator_mxll_t),    intent(inout) :: tr
    real(real64),               intent(in)    :: time
    real(real64),               intent(in)    :: dt
    real(real64),               intent(in)    :: time_delay
    type(batch_t),              intent(inout) :: ff_rs_state_pmlb
    type(batch_t),              intent(inout) :: ff_rs_stateb
 
    integer            :: ii, ff_dim
    complex(real64), allocatable :: rs_state_constant(:,:), ff_rs_state_constant(:,:)
    type(batch_t) :: ff_rs_state_plane_wavesb, ff_rs_constantb, rs_state_constantb
 
    push_sub(pml_propagation_stage_2_batch)
 
    call profiling_in('PML_PROP_STAGE_2_BATCH')
 
    if (tr%bc_plane_waves .and. hm%plane_waves_apply) then
      hm%cpml_hamiltonian = .true.
      call tr%te%apply_batch(namespace, gr, hm, ff_rs_state_pmlb, dt)
      hm%cpml_hamiltonian = .false.
      call plane_waves_propagation(hm, tr, namespace, st, gr, time, dt, time_delay)
 
      call ff_rs_stateb%copy_to(ff_rs_state_plane_wavesb)
      call transform_rs_state_batch(hm, gr, st, st%rs_state_plane_wavesb, ff_rs_state_plane_wavesb, rs_trans_forward)
 
      if (ff_rs_stateb%status() == batch_device_packed) then
        ! use the map of points stored on the GPU in this case
        call batch_add_with_map(hm%bc%plane_wave%points_number, hm%bc%plane_wave%buff_map, &
          ff_rs_state_pmlb, ff_rs_state_plane_wavesb, ff_rs_stateb)
      else
        call batch_add_with_map(hm%bc%plane_wave%points_number, hm%bc%plane_wave%points_map, &
          ff_rs_state_pmlb, ff_rs_state_plane_wavesb, ff_rs_stateb)
      end if
 
      call ff_rs_state_plane_wavesb%end()
 
    else if (tr%bc_constant .and. hm%spatial_constant_apply) then
      hm%cpml_hamiltonian = .true.
      call tr%te%apply_batch(namespace, gr, hm, ff_rs_state_pmlb, dt)
      hm%cpml_hamiltonian = .false.
 
      call ff_rs_stateb%copy_to(ff_rs_constantb)
      ff_dim = ff_rs_stateb%nst_linear
      safe_allocate(rs_state_constant(1:gr%np, 1:st%dim))
      ! copy the value to the full mesh to be able to use batches
      ! this is in principle unneeded, but otherwise we could not use batches...
      do ii = 1, st%dim
        rs_state_constant(1:gr%np, ii) = st%rs_state_const(ii)
      end do
      call ff_rs_stateb%copy_to(rs_state_constantb)
      call mxll_set_batch(rs_state_constantb, rs_state_constant, gr%np, st%dim)
 
      call transform_rs_state_batch(hm, gr, st, rs_state_constantb, ff_rs_constantb, rs_trans_forward)
      if (ff_rs_stateb%status() == batch_device_packed) then
        ! use the map of points stored on the GPU in this case
        call batch_add_with_map(hm%bc%constant_points_number, hm%bc%buff_constant_points_map, &
          ff_rs_state_pmlb, ff_rs_constantb, ff_rs_stateb)
      else
        call batch_add_with_map(hm%bc%constant_points_number, hm%bc%constant_points_map, &
          ff_rs_state_pmlb, ff_rs_constantb, ff_rs_stateb)
      end if
 
      call ff_rs_constantb%end()
      call rs_state_constantb%end()
 
      safe_deallocate_a(rs_state_constant)
      safe_deallocate_a(ff_rs_state_constant)
    end if
 
    call profiling_out('PML_PROP_STAGE_2_BATCH')
 
    pop_sub(pml_propagation_stage_2_batch)
  end subroutine pml_propagation_stage_2_batch
 
  ! ---------------------------------------------------------
  subroutine cpml_conv_function_update(hm, gr, ff_rs_state_pmlb)
    type(hamiltonian_mxll_t), intent(inout) :: hm
    type(grid_t),             intent(in)    :: gr
    type(batch_t),            intent(inout) :: ff_rs_state_pmlb
 
 
    push_sub(cpml_conv_function_update)
 
    call profiling_in('CPML_CONV_FUNCTION_UPDATE')
 
    call cpml_conv_function_update_via_riemann_silberstein(hm, gr, ff_rs_state_pmlb)
 
    call profiling_out('CPML_CONV_FUNCTION_UPDATE')
 
    pop_sub(cpml_conv_function_update)
  end subroutine cpml_conv_function_update
 
  ! ---------------------------------------------------------
  subroutine cpml_conv_function_update_via_riemann_silberstein(hm, gr, ff_rs_state_pmlb)
    type(hamiltonian_mxll_t), intent(inout) :: hm
    type(grid_t),             intent(in)    :: gr
    type(batch_t),            intent(inout) :: ff_rs_state_pmlb
 
    integer :: ip, ip_in, np_part, rs_sign
    complex(real64) :: pml_a, pml_b, pml_g, grad
    integer :: pml_dir, field_dir, ifield, idir
    integer, parameter :: field_dirs(3, 2) = reshape([2, 3, 1, 3, 1, 2], [3, 2])
    logical :: with_medium
    type(batch_t) :: gradb(gr%der%dim)
    type(accel_kernel_t), save :: ker_pml
    integer :: wgsize
 
    push_sub(cpml_conv_function_update_via_riemann_silberstein)
 
    call profiling_in('CPML_CONV_UPDATE_VIA_RS')
 
    assert(hm%dim == 3 .or. hm%dim == 6)
 
    np_part = gr%np_part
    rs_sign = hm%rs_sign
 
    call zderivatives_batch_grad(gr%der, ff_rs_state_pmlb, gradb)
 
    with_medium = hm%dim == 6
 
    do pml_dir = 1, hm%st%dim
      select case (gradb(pml_dir)%status())
      case (batch_not_packed)
        do ip_in=1, hm%bc%pml%points_number
          ip       = hm%bc%pml%points_map(ip_in)
          pml_a = hm%bc%pml%a(ip_in,pml_dir)
          pml_b = hm%bc%pml%b(ip_in,pml_dir)
          do ifield = 1, 2
            field_dir = field_dirs(pml_dir, ifield)
            grad = gradb(pml_dir)%zff_linear(ip, field_dir)
            pml_g = hm%bc%pml%conv_plus(ip_in, pml_dir, field_dir)
            hm%bc%pml%conv_plus(ip_in, pml_dir, field_dir) = pml_a * grad + pml_b * pml_g
            if (with_medium) then
              grad = gradb(pml_dir)%zff_linear(ip, field_dir+3)
              pml_g = hm%bc%pml%conv_minus(ip_in, pml_dir, field_dir)
              hm%bc%pml%conv_minus(ip_in, pml_dir, field_dir) = pml_a * grad + pml_b * pml_g
            end if
          end do
        end do
      case (batch_packed)
        do ip_in=1, hm%bc%pml%points_number
          ip       = hm%bc%pml%points_map(ip_in)
          pml_a = hm%bc%pml%a(ip_in,pml_dir)
          pml_b = hm%bc%pml%b(ip_in,pml_dir)
          do ifield = 1, 2
            field_dir = field_dirs(pml_dir, ifield)
            grad = gradb(pml_dir)%zff_pack(field_dir, ip)
            pml_g = hm%bc%pml%conv_plus(ip_in, pml_dir, field_dir)
            hm%bc%pml%conv_plus(ip_in, pml_dir, field_dir) = pml_a * grad + pml_b * pml_g
            if (with_medium) then
              grad = gradb(pml_dir)%zff_pack(field_dir+3, ip)
              pml_g = hm%bc%pml%conv_minus(ip_in, pml_dir, field_dir)
              hm%bc%pml%conv_minus(ip_in, pml_dir, field_dir) = pml_a * grad + pml_b * pml_g
            end if
          end do
        end do
      case (batch_device_packed)
        call accel_kernel_start_call(ker_pml, 'pml.cl', 'pml_update_conv')
 
        if (with_medium) then
          call accel_set_kernel_arg(ker_pml, 0, 1_int32)
        else
          call accel_set_kernel_arg(ker_pml, 0, 0_int32)
        end if
        call accel_set_kernel_arg(ker_pml, 1, hm%bc%pml%points_number)
        call accel_set_kernel_arg(ker_pml, 2, pml_dir-1)
        call accel_set_kernel_arg(ker_pml, 3, hm%bc%pml%buff_map)
        call accel_set_kernel_arg(ker_pml, 4, gradb(pml_dir)%ff_device)
        call accel_set_kernel_arg(ker_pml, 5, log2(int(gradb(pml_dir)%pack_size(1), int32)))
        call accel_set_kernel_arg(ker_pml, 6, hm%bc%pml%buff_a)
        call accel_set_kernel_arg(ker_pml, 7, hm%bc%pml%buff_b)
        call accel_set_kernel_arg(ker_pml, 8, hm%bc%pml%buff_conv_plus)
        call accel_set_kernel_arg(ker_pml, 9, hm%bc%pml%buff_conv_minus)
 
        wgsize = accel_max_workgroup_size()
        call accel_kernel_run(ker_pml, (/ accel_padded_size(hm%bc%pml%points_number) /), (/ wgsize /))
      end select
    end do
 
    do idir = 1, gr%der%dim
      call gradb(idir)%end()
    end do
 
    if (accel_is_enabled()) then
      call accel_finish()
    end if
 
    call profiling_out('CPML_CONV_UPDATE_VIA_RS')
 
    pop_sub(cpml_conv_function_update_via_riemann_silberstein)
  end subroutine cpml_conv_function_update_via_riemann_silberstein
 
  ! ---------------------------------------------------------
  subroutine td_function_mxll_init(st, namespace, hm)
    type(states_mxll_t),      intent(inout) :: st
    type(namespace_t),        intent(in)    :: namespace
    type(hamiltonian_mxll_t), intent(inout)    :: hm
 
    type(block_t)        :: blk
    integer              :: il, nlines, idim, ncols, ierr
    real(real64)         :: e_field(st%dim), b_field(st%dim)
    character(len=1024)  :: mxf_expression
 
    push_sub(td_function_mxll_init)
 
    call profiling_in('TD_FUNCTION_MXLL_INIT')
 
    !%Variable UserDefinedConstantSpatialMaxwellField
    !%Type block
    !%Section Maxwell
    !%Description
    !% Define parameters of spatially constant field.
    !%
    !% Example:
    !%
    !% <tt>%UserDefinedConstantSpatialMaxwellFields
    !% <br>&nbsp;&nbsp;   plane_wave_parser      | E_x | E_y | E_z | B_x | B_y | B_z | "tdf_function"
    !% <br>%</tt>
    !%
    !% This block defines three components of E field, three components of B field, and reference to
    !% the TD function.
    !%
    !%End
 
    if (parse_block(namespace, 'UserDefinedConstantSpatialMaxwellField', blk) == 0) then
      st%rs_state_const_external = .true.
      nlines = parse_block_n(blk)
      safe_allocate(st%rs_state_const_td_function(1:nlines))
      safe_allocate(st%rs_state_const_amp(1:st%dim, 1:nlines))
      ! read all lines
      do il = 1, nlines
        e_field = m_zero
        b_field = m_zero
        ! Check that number of columns is five or six.
        ncols = parse_block_cols(blk, il - 1)
        if (ncols  /= 7) then
          message(1) = 'Each line in the UserDefinedConstantSpatialMaxwellField block must have'
          message(2) = 'seven columns.'
          call messages_fatal(2, namespace=namespace)
        end if
        do idim = 1, st%dim
          call parse_block_float( blk, il - 1, idim-1, e_field(idim))
        end do
        do idim = 1, st%dim
          call parse_block_float( blk, il - 1, idim+2, b_field(idim))
        end do
        call parse_block_string( blk, il - 1, 6, mxf_expression)
        call build_rs_vector(e_field, b_field, st%rs_sign, st%rs_state_const_amp(:,il))
        call tdf_read(st%rs_state_const_td_function(il), namespace, trim(mxf_expression), ierr)
      end do
    end if
    call parse_block_end(blk)
 
    !%Variable PropagateSpatialMaxwellField
    !%Type logical
    !%Default yes
    !%Section Maxwell::TD Propagation
    !%Description
    !% Allow for numerical propagation of Maxwells equations of spatially constant field.
    !% If set to no, do only analytic evaluation of the field inside the box.
    !%End
 
    call parse_variable(namespace, 'PropagateSpatialMaxwellField', .true., hm%spatial_constant_propagate)
 
    call profiling_out('TD_FUNCTION_MXLL_INIT')
 
    pop_sub(td_function_mxll_init)
  end subroutine td_function_mxll_init
 
  ! ---------------------------------------------------------
  subroutine spatial_constant_calculation(constant_calc, st, gr, hm, time, dt, delay, rs_state, set_initial_state)
    logical,                  intent(in)    :: constant_calc
    type(states_mxll_t),      intent(inout) :: st
    type(grid_t),             intent(in)    :: gr
    type(hamiltonian_mxll_t), intent(in)    :: hm
    real(real64),             intent(in)    :: time
    real(real64),             intent(in)    :: dt
    real(real64),             intent(in)    :: delay
    complex(real64),          intent(inout) :: rs_state(:,:)
    logical,        optional, intent(in)    :: set_initial_state
 
    integer :: ip, ic, icn
    real(real64)   :: tf_old, tf_new
    logical :: set_initial_state_
 
    push_sub(spatial_constant_calculation)
 
    call profiling_in('SPATIAL_CONSTANT_CALCULATION')
 
    set_initial_state_ = .false.
    if (present(set_initial_state)) set_initial_state_ = set_initial_state
 
    if (hm%spatial_constant_apply) then
      if (constant_calc) then
        icn = size(st%rs_state_const_td_function(:))
        st%rs_state_const(:) = m_z0
        do ic = 1, icn
          tf_old = tdf(st%rs_state_const_td_function(ic), time-delay-dt)
          tf_new = tdf(st%rs_state_const_td_function(ic), time-delay)
          do ip = 1, gr%np
            if (set_initial_state_ .or. (.not. hm%spatial_constant_propagate)) then
              rs_state(ip,:) = st%rs_state_const_amp(:,ic) * tf_new
            else
              rs_state(ip,:) = rs_state(ip,:) + st%rs_state_const_amp(:,ic) * (tf_new - tf_old)
            end if
          end do
          st%rs_state_const(:) = st%rs_state_const(:) + st%rs_state_const_amp(:, ic)
        end do
        st%rs_state_const(:) = st%rs_state_const(:) * tf_new
      end if
    end if
 
    call profiling_out('SPATIAL_CONSTANT_CALCULATION')
 
    pop_sub(spatial_constant_calculation)
  end subroutine spatial_constant_calculation
 
  ! ---------------------------------------------------------
  subroutine constant_boundaries_calculation(constant_calc, bc, hm, st, rs_state)
    logical,                   intent(in)    :: constant_calc
    type(bc_mxll_t),           intent(inout) :: bc
    type(states_mxll_t),       intent(in)    :: st
    type(hamiltonian_mxll_t),  intent(in)    :: hm
    complex(real64),           intent(inout) :: rs_state(:,:)
 
    integer :: ip_in, ip
 
    push_sub(constant_boundaries_calculation)
    call profiling_in('CONSTANT_BOUNDARIES_CALC')
 
    if (hm%spatial_constant_apply) then
      if (constant_calc) then
        do ip_in = 1, bc%constant_points_number
          ip = bc%constant_points_map(ip_in)
          rs_state(ip,:) = st%rs_state_const(:)
          bc%constant_rs_state(ip_in,:) = st%rs_state_const(:)
        end do
      end if
    end if
 
    call profiling_out('CONSTANT_BOUNDARIES_CALC')
 
    pop_sub(constant_boundaries_calculation)
  end subroutine constant_boundaries_calculation
 
  ! ---------------------------------------------------------
  subroutine mirror_pec_boundaries_calculation(bc, st, rs_state)
    type(bc_mxll_t),     intent(in)    :: bc
    type(states_mxll_t), intent(in)    :: st
    complex(real64),     intent(inout) :: rs_state(:,:)
 
    integer                    :: ip, ip_in, idim
    real(real64)               :: e_field(st%dim), b_field(st%dim)
 
    push_sub(mirror_pec_boundaries_calculation)
 
    do idim = 1, 3
      if (bc%bc_type(idim) == mxll_bc_mirror_pec) then
        do ip_in = 1, bc%mirror_points_number(idim)
          ip = bc%mirror_points_map(ip_in, idim)
          e_field(:) = m_zero
          call get_magnetic_field_vector(rs_state(ip,:), st%rs_sign, b_field(:), st%mu(ip))
          call build_rs_vector(e_field(:), b_field(:), st%rs_sign, rs_state(ip,:), st%ep(ip), st%mu(ip))
        end do
      end if
    end do
 
    pop_sub(mirror_pec_boundaries_calculation)
  end subroutine mirror_pec_boundaries_calculation
 
  ! ---------------------------------------------------------
  subroutine mirror_pmc_boundaries_calculation(bc, st, rs_state)
    type(bc_mxll_t),     intent(in)    :: bc
    type(states_mxll_t), intent(in)    :: st
    complex(real64),     intent(inout) :: rs_state(:,:)
 
    integer                    :: ip, ip_in, idim
    real(real64)               :: e_field(st%dim), b_field(st%dim)
 
    push_sub(mirror_pmc_boundaries_calculation)
 
    do idim = 1, 3
      if (bc%bc_type(idim) == mxll_bc_mirror_pmc) then
        do ip_in = 1, bc%mirror_points_number(idim)
          ip = bc%mirror_points_map(ip_in,idim)
          b_field(:) = m_zero
          call get_electric_field_vector(rs_state(ip,:), e_field(:), st%ep(ip))
          call build_rs_vector(e_field(:), b_field(:), st%rs_sign, rs_state(ip,:), st%ep(ip), st%mu(ip))
        end do
      end if
    end do
 
    pop_sub(mirror_pmc_boundaries_calculation)
  end subroutine mirror_pmc_boundaries_calculation
 
  ! ---------------------------------------------------------
  subroutine plane_waves_boundaries_calculation(hm, st, mesh, time, time_delay, rs_state)
    type(hamiltonian_mxll_t), intent(in)    :: hm
    type(states_mxll_t),      intent(in)    :: st
    class(mesh_t),            intent(in)    :: mesh
    real(real64),             intent(in)    :: time
    real(real64),             intent(in)    :: time_delay
    complex(real64),          intent(inout) :: rs_state(:,:)
 
    integer                    :: ip, ip_in, wn
    real(real64)               :: x_prop(mesh%box%dim), rr, vv(mesh%box%dim), k_vector(mesh%box%dim)
    real(real64)               :: k_vector_abs, nn
    complex(real64)            :: e0(mesh%box%dim)
    real(real64)               :: e_field(mesh%box%dim), b_field(mesh%box%dim)
    complex(real64)            :: rs_state_add(mesh%box%dim)
    complex(real64)            :: mx_func
 
    push_sub(plane_waves_boundaries_calculation)
 
    call profiling_in('PLANE_WAVES_BOUNDARIES_C')
 
    if (hm%plane_waves_apply) then
      do wn = 1, hm%bc%plane_wave%number
        k_vector(:) = hm%bc%plane_wave%k_vector(1:mesh%box%dim, wn)
        k_vector_abs = norm2(k_vector(1:mesh%box%dim))
        vv(:) = hm%bc%plane_wave%v_vector(1:mesh%box%dim, wn)
        e0(:) = hm%bc%plane_wave%e_field(1:mesh%box%dim, wn)
        do ip_in = 1, hm%bc%plane_wave%points_number
          ip = hm%bc%plane_wave%points_map(ip_in)
          if (wn == 1) rs_state(ip,:) = m_z0
          nn = sqrt(st%ep(ip)/p_ep*st%mu(ip)/p_mu)
          x_prop(1:mesh%box%dim) = mesh%x(ip,1:mesh%box%dim) - vv(1:mesh%box%dim) * (time - time_delay)
          rr = norm2(x_prop(1:mesh%box%dim))
          if (hm%bc%plane_wave%modus(wn) == option__maxwellincidentwaves__plane_wave_mx_function) then
            ! Temporary variable assigned due to macro line length
            mx_func = mxf(hm%bc%plane_wave%mx_function(wn), x_prop(1:mesh%box%dim))
            e_field(1:mesh%box%dim) = real(e0(1:mesh%box%dim) * mx_func, real64)
          end if
          b_field(1:3) = dcross_product(k_vector, e_field) / p_c / k_vector_abs
          call build_rs_vector(e_field, b_field, st%rs_sign, rs_state_add, st%ep(ip), st%mu(ip))
          rs_state(ip, :) =  rs_state(ip, :) + rs_state_add(:)
        end do
      end do
    else
      do ip_in = 1, hm%bc%plane_wave%points_number
        ip             = hm%bc%plane_wave%points_map(ip_in)
        rs_state(ip,:) = m_z0
      end do
    end if
 
    call profiling_out('PLANE_WAVES_BOUNDARIES_C')
 
    pop_sub(plane_waves_boundaries_calculation)
  end subroutine plane_waves_boundaries_calculation
 
  ! ---------------------------------------------------------
  subroutine plane_waves_propagation(hm, tr, namespace, st, gr, time, dt, time_delay)
    type(hamiltonian_mxll_t), intent(inout) :: hm
    type(propagator_mxll_t),  intent(inout) :: tr
    type(namespace_t),        intent(in)    :: namespace
    type(states_mxll_t),      intent(inout) :: st
    type(grid_t),             intent(in)    :: gr
    real(real64),             intent(in)    :: time
    real(real64),             intent(in)    :: dt
    real(real64),             intent(in)    :: time_delay
 
    type(batch_t) :: ff_rs_stateb
    integer :: ff_dim
 
    push_sub(plane_waves_propagation)
 
    call profiling_in('PLANE_WAVES_PROPAGATION')
 
    ff_dim = hm%dim
    call zbatch_init(ff_rs_stateb, 1, 1, hm%dim, gr%np_part)
    if (st%pack_states) call ff_rs_stateb%do_pack(copy=.false.)
 
    call transform_rs_state_batch(hm, gr, st, st%rs_state_plane_wavesb, ff_rs_stateb, rs_trans_forward)
 
    ! Time evolution of RS plane waves state without any coupling with H(inter_time)
    call hamiltonian_mxll_update(hm, time=time)
    hm%cpml_hamiltonian = .false.
    call tr%te%apply_batch(namespace, gr, hm, ff_rs_stateb, dt)
 
    call transform_rs_state_batch(hm, gr, st, st%rs_state_plane_wavesb, ff_rs_stateb, rs_trans_backward)
    call mxll_get_batch(st%rs_state_plane_wavesb, st%rs_state_plane_waves, gr%np, st%dim)
    call plane_waves_boundaries_calculation(hm, st, gr, time+dt, time_delay, st%rs_state_plane_waves)
    call mxll_set_batch(st%rs_state_plane_wavesb, st%rs_state_plane_waves, gr%np, st%dim)
    call ff_rs_stateb%end()
 
    call profiling_out('PLANE_WAVES_PROPAGATION')
    pop_sub(plane_waves_propagation)
  end subroutine plane_waves_propagation
 
  ! ---------------------------------------------------------
  subroutine plane_waves_in_box_calculation(bc, time, mesh, der, st, rs_state)
    type(bc_mxll_t),           intent(inout) :: bc
    real(real64),              intent(in)    :: time
    class(mesh_t),             intent(in)    :: mesh
    type(derivatives_t),       intent(in)    :: der
    type(states_mxll_t),       intent(in)    :: st
    complex(real64),           intent(inout) :: rs_state(:,:)
 
    real(real64)         :: e_field_total(mesh%np,st%dim), b_field_total(mesh%np,st%dim)
    complex(real64)      :: rs_state_add(mesh%np,st%dim)
 
    push_sub(plane_waves_in_box_calculation)
 
    call profiling_in('PLANE_WAVES_IN_BOX_CALCULATION')
 
    call external_waves_eval(bc%plane_wave, time, mesh, "E field", e_field_total)
    call external_waves_eval(bc%plane_wave, time, mesh, "B field", b_field_total, der=der)
 
    call build_rs_state(e_field_total,  b_field_total, st%rs_sign, &
      rs_state_add(1:mesh%np,:), mesh, st%ep, st%mu)
    rs_state(1:mesh%np,:) = rs_state(1:mesh%np,:) + rs_state_add(1:mesh%np,:)
 
    call profiling_out('PLANE_WAVES_IN_BOX_CALCULATION')
 
    pop_sub(plane_waves_in_box_calculation)
  end subroutine plane_waves_in_box_calculation
 
  ! ---------------------------------------------------------
  subroutine mxll_apply_boundaries(tr, st, hm, gr, namespace, time, dt, rs_stateb)
    type(propagator_mxll_t), intent(inout) :: tr
    type(states_mxll_t),     intent(inout) :: st
    type(hamiltonian_mxll_t),intent(inout) :: hm
    type(grid_t),            intent(in)    :: gr
    type(namespace_t),       intent(in)    :: namespace
    real(real64),            intent(in)    :: time
    real(real64),            intent(in)    :: dt
    type(batch_t),           intent(inout) :: rs_stateb
 
    complex(real64), allocatable :: rs_state(:, :)
 
    push_sub(mxll_apply_boundaries)
 
    safe_allocate(rs_state(gr%np, st%dim))
 
    if (tr%bc_constant) then
      call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
      ! Propagation dt with H(inter_time+inter_dt) for constant boundaries
      if (st%rs_state_const_external) then
        call spatial_constant_calculation(tr%bc_constant, st, gr, hm, time, dt, m_zero, rs_state)
      end if
      call constant_boundaries_calculation(tr%bc_constant, hm%bc, hm, st, rs_state)
      call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
    end if
 
    ! PEC mirror boundaries
    if (any(hm%bc%bc_type == mxll_bc_mirror_pec)) then
      call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
      call mirror_pec_boundaries_calculation(hm%bc, st, rs_state)
      call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
    end if
 
    ! PMC mirror boundaries
    if (any(hm%bc%bc_type == mxll_bc_mirror_pmc)) then
      call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
      call mirror_pmc_boundaries_calculation(hm%bc, st, rs_state)
      call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
    end if
 
    if (any(hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__mask)) then
      call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
      ! Apply mask absorbing boundaries
      call mask_absorbing_boundaries(namespace, gr, hm, st, tr, time, dt, m_zero, rs_state)
      call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
    end if
 
    if (tr%bc_plane_waves) then
      call mxll_get_batch(rs_stateb, rs_state, gr%np, st%dim)
      ! calculate plane waves boundaries at t
      call plane_waves_boundaries_calculation(hm, st, gr, time, m_zero, rs_state)
      call mxll_set_batch(rs_stateb, rs_state, gr%np, st%dim)
    end if
 
    safe_deallocate_a(rs_state)
 
    pop_sub(mxll_apply_boundaries)
  end subroutine mxll_apply_boundaries
 
end module propagator_mxll_oct_m