external_waves.F90

Go to the documentation of this file.

!! Copyright (C) 2023 E.I. Albar, F. Bonafe and Heiko 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 external_waves_oct_m
  use accel_oct_m
  use box_sphere_oct_m
  use box_parallelepiped_oct_m
  use clock_oct_m
  use cube_function_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use ghost_interaction_oct_m
  use global_oct_m
  use grid_oct_m
  use index_oct_m
  use interaction_surrogate_oct_m
  use interaction_enum_oct_m
  use interaction_partner_oct_m
  use interaction_oct_m
  use, intrinsic :: iso_fortran_env
  use lorentz_force_oct_m
  use io_oct_m
  use io_function_oct_m
  use lasers_oct_m
  use loct_math_oct_m
  use maxwell_function_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mxll_b_field_to_matter_oct_m
  use mxll_e_field_to_matter_oct_m
  use mxll_vec_pot_to_matter_oct_m
  use par_vec_oct_m
  use parser_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use states_elec_oct_m
  use string_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use states_mxll_oct_m
 
  implicit none
 
  private
  public ::                    &
    external_waves_t,          &
    external_waves_init,       &
    external_waves_eval,       &
    bessel_source_eval,        &
    external_waves_end,        &
    load_external_waves
 
  type bessel_beam_t
    integer,  allocatable :: helicity(:)
    integer,  allocatable :: m_order(:)
    real(real64),    allocatable :: amp(:)
    real(real64),    allocatable :: theta_k(:)
    real(real64),    allocatable :: omega(:)
    real(real64),    allocatable :: shift(:,:)
    logical,  allocatable :: envelope(:)
    integer,  allocatable :: lin_dir(:)
  contains
    procedure :: init  => bessel_beam_init
    procedure :: function => bessel_beam_function
    final     :: bessel_beam_finalize
  end type bessel_beam_t
 
  type, extends(interaction_partner_t) :: external_waves_t
    integer                          :: points_number
    integer,             allocatable :: points_map(:)
    integer                          :: number
    integer,             allocatable :: modus(:)
    integer,             allocatable :: field_type(:)
    character(len=1024), allocatable :: e_field_string(:,:)
    real(real64),        allocatable :: k_vector(:,:)
    real(real64),        allocatable :: v_vector(:,:)
    complex(real64),     allocatable :: e_field(:,:)
    real(real64),        allocatable :: pw_phase(:)
    type(mxf_t),         allocatable :: mx_function(:)
    integer                          :: out_file
    logical                          :: output_from_point = .false. 
    real(real64), allocatable               :: selected_point_coordinate(:)
    real(real64), allocatable               :: selected_point_field(:)
    real(real64)                     :: c_factor
    type(accel_mem_t)                :: buff_map
    type(bessel_beam_t)                :: bessel
  contains
    procedure :: init_interaction_as_partner => external_waves_init_interaction_as_partner
    procedure :: update_quantity => external_waves_update_quantity
    procedure :: copy_quantities_to_interaction => external_waves_copy_quantities_to_interaction
    final :: external_waves_end
  end type external_waves_t
 
  interface external_waves_t
    module procedure external_waves_constructor
  end interface external_waves_t
 
contains
 
  function external_waves_constructor(namespace) result(this)
    class(external_waves_t), pointer :: this
    type(namespace_t), intent(in) :: namespace
 
    integer :: iq
    character(len=:), allocatable :: quantities(:)
 
    push_sub(external_waves_constructor)
 
    safe_allocate(this)
 
    this%namespace = namespace_t("ExternalSource", parent=namespace)
 
    message(1) = 'Plane-wave is currently always 3D and non-periodic.'
    call messages_warning(1)
    call external_waves_init(this, this%namespace)
 
    quantities = [character(16) :: "E field", "vector potential", "B field"]
    do iq = 1, size(quantities)
      call this%quantities%add(quantity_t(quantities(iq), always_available = .true., &
        updated_on_demand = .true., iteration = clock_t()))
    end do
 
    this%supported_interactions_as_partner = [mxll_e_field_to_matter, mxll_b_field_to_matter, &
      mxll_vec_pot_to_matter]
 
    pop_sub(external_waves_constructor)
  end function external_waves_constructor
 
  ! ---------------------------------------------------------
  subroutine external_waves_init_interaction_as_partner(partner, interaction)
    class(external_waves_t),    intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(external_waves_init_interaction_as_partner)
 
    select type (interaction)
    type is (lorentz_force_t)
      ! Nothing to be initialized
    type is (mxll_e_field_to_matter_t)
      ! Nothing to be initialized
    type is (mxll_vec_pot_to_matter_t)
      ! Nothing to be initialized
    type is (mxll_b_field_to_matter_t)
      ! Nothing to be initialized
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
    pop_sub(external_waves_init_interaction_as_partner)
  end subroutine external_waves_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine external_waves_update_quantity(this, label)
    class(external_waves_t),  intent(inout) :: this
    character(len=*),         intent(in)    :: label
 
    push_sub(external_waves_update_quantity)
 
    select case (label)
    case ("E field", "B field", "vector potential")
      ! We will not update the quantities here because they are computed
      ! on-the-fly when copying them to the corresponding interaction (see
      ! copy_quantities_to_interaction routine)
    case default
      message(1) = "Incompatible quantity."
      call messages_fatal(1, namespace=this%namespace)
    end select
 
    pop_sub(external_waves_update_quantity)
  end subroutine external_waves_update_quantity
 
  ! ---------------------------------------------------------
  subroutine external_waves_copy_quantities_to_interaction(partner, interaction)
    class(external_waves_t), intent(inout) :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    class(quantity_t), pointer :: quantity
 
    push_sub(external_waves_copy_quantities_to_interaction)
 
    select type(interaction)
    class is (mxll_e_field_to_matter_t)
      quantity => partner%quantities%get("E field")
      interaction%system_field = m_zero
      call external_waves_eval(partner, quantity%iteration%value(), interaction%system_gr, "E field", &
        interaction%system_field)
      call interaction%do_mapping()
 
    class is (mxll_vec_pot_to_matter_t)
      quantity => partner%quantities%get("vector potential")
      interaction%system_field = m_zero
      call external_waves_eval(partner, quantity%iteration%value(), interaction%system_gr, &
        "vector potential", interaction%system_field)
      call interaction%do_mapping()
 
    class is (mxll_b_field_to_matter_t)
      quantity => partner%quantities%get("B field")
      interaction%system_field = m_zero
      call external_waves_eval(partner, quantity%iteration%value(), interaction%system_gr, "B field", &
        interaction%system_field, der=interaction%system_gr%der)
      call interaction%do_mapping()
 
    class default
      message(1) = "Incompatible interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(external_waves_copy_quantities_to_interaction)
  end subroutine external_waves_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  ! Load the external source for the multisystem framework
  subroutine load_external_waves(partners, namespace)
    class(partner_list_t), intent(inout)  :: partners
    type(namespace_t),     intent(in)     :: namespace
 
    logical :: has_source
 
    push_sub(load_external_waves)
 
    !%Variable AnalyticalExternalSource
    !%Type logical
    !%Default no
    !%Section Maxwell
    !%Description
    !% This means the analytical evaluation of formula will be used, Maxwell propagation will not be used.
    !%End
    call parse_variable(namespace, 'AnalyticalExternalSource', .false., has_source)
 
    if (has_source) then
      call partners%add(external_waves_t(namespace))
    end if
 
    pop_sub(load_external_waves)
  end subroutine load_external_waves
 
 
  ! ---------------------------------------------------------
  subroutine external_waves_init(external_waves, namespace)
    type(external_waves_t),    intent(inout) :: external_waves
    type(namespace_t),         intent(in)    :: namespace
    type(block_t)        :: blk
    integer              :: il, nlines, ncols, iex_norm, idim
    integer, parameter   :: sys_dim = 3
    real(real64)         :: k_vector(sys_dim), velocity(sys_dim), x_pos(sys_dim)
    real(real64)         :: x_norm, dummy(sys_dim), k_dot_e , test_limit, k_norm, output_pos(3)
    complex(real64)      :: e_field(sys_dim)
    character(len=1024)  :: k_string(sys_dim)
    character(len=1), dimension(sys_dim), parameter :: dims = ["x", "y", "z"]
    character(len=1024)  :: mxf_expression
 
    push_sub(external_waves_init)
 
    call profiling_in('EXTERNAL_WAVES_INIT')
 
    test_limit = 10.0e-9_real64
 
    !%Variable ExternalSourceBesselOutput
    !%Type block
    !%Section Maxwell
    !%Description
    !%  The ExternalSourceBesselOutput block allows to output analytically calculated fields at a
    !%  particular point in space. The columns denote the x, y, and z coordinate of the point.
    !%  Please be aware that ExternalSource lives on the grid of the system that it is applied to.
    !%  Therefore, it might not be evaluated at every point in space. When comparing, please be sure
    !%  to check the log and compare if your required point in space matches the evaluated position.
    !%
    !% <tt>%ExternalSourceBesselOutput
    !% <br>&nbsp;&nbsp;    -1.0 | 2.0 |  4.0
    !% <br>%</tt>
    !%
    !%End
 
    if (parse_block(namespace, 'ExternalSourceBesselOutput', blk) == 0) then
      nlines = parse_block_n(blk)
      if (nlines > 1 ) then
        message(2) = 'ExternalSource output is limited to one point.'
        call messages_fatal(1, namespace=namespace)
      end if
      ncols = parse_block_cols(blk,0)
      if (ncols /= 3 ) then
        message(1) = 'ExternalSourceBesselOutput must have 3 columns.'
        call messages_fatal(1, namespace=namespace)
      end if
      external_waves%output_from_point= .true.
      safe_allocate(external_waves%selected_point_coordinate(1:3))
      safe_allocate(external_waves%selected_point_field(1:3))
 
      do idim = 1, 3
        call parse_block_float(blk, 0, idim-1, output_pos(idim), units_inp%length)
      end do
      external_waves%selected_point_coordinate(1:3) = output_pos(1:3)
      external_waves%selected_point_field(1:3)  = m_zero
 
      call parse_block_end(blk)
      call io_mkdir('ExternalSource')
      external_waves%out_file  = io_open('bessel_source_at_point', namespace=namespace, action='write')
      write(external_waves%out_file, '(12a) ')  '#    Time (a.u.)    ' , '  Field-x  ' , '  Field-y  ' , '   Field-z   '
 
    else
      external_waves%output_from_point= .false.
 
    end if
 
    ! This variable is documented in hamiltonian_mxll.F90
    call parse_variable(namespace, 'SpeedOfLightFactor', m_one, external_waves%c_factor)
 
    !%Variable MaxwellIncidentWaves
    !%Type block
    !%Section Maxwell
    !%Description
    !% The initial electromagnetic fields can be set by the user
    !% with the <tt>MaxwellIncidentWaves</tt> block variable.
    !% The electromagnetic fields have to fulfill the
    !% Maxwells equations in vacuum.
    !% For a Maxwell propagation, setting the electric field is sufficient,
    !% the magnetic field (for plane waves) will be calculated from it as 1/(c.|k|) . (k x E).
    !%
    !% Example:
    !%
    !% <tt>%MaxwellIncidentWaves
    !% <br>&nbsp;&nbsp;   plane_wave_parser      | "field_type"| "k1x" | "k1y" | "k1z" | "E1x" | "E1z" | "E1x"
    !% <br>&nbsp;&nbsp;   plane_wave_mx_function | "field_type"| "E4x" | "E4y" | "E4z" | mx_envelope_name | phase
    !% <br>&nbsp;&nbsp;   bessel_function | "field_type"| A_0 | m | omega | helicity | <math>\theta_{k}</math> | mx_envelope_name | lin_dir
    !% <br>%</tt>
    !%
    !% Field type can be "electric_field" or "vector_potential". Note that in order to couple to an
    !% electronic system, the <tt>MaxwellCouplingMode</tt> variable needs to be set to a coupling type
    !% compatible with the requested field type ("electric_field" is compatible with length gauge,
    !% while "vector_potential" is compatible with velocity gauge and full minimal coupling).
    !% Otherwise, the field will not be calculated or applied to the electronic Hamiltonian.
    !%
    !%Option plane_wave_parser 0
    !% Parser input modus
    !%Option plane_wave_mx_function 1
    !% The incident wave envelope is defined by an mx_function
    !%Option bessel_function 2
    !% The incident source is a generalized Bessel beam, parametrized by its amplitude, opening angle, helicity, order and frequency.
    !% This beam is a solution of Maxwell equations, and inherently circularly polarized and is parametrized by its amplitude,
    !% opening angle, helicity, order and frequency.
    !% Please keep in mind, if you set linear polarization lin_dir,
    !% you will obtain a linearly polarized Bessel beam.
    !%End
 
    if (parse_block(namespace, 'MaxwellIncidentWaves', blk) == 0) then
 
      call messages_print_with_emphasis(msg='Substitution of the electromagnetic incident waves', namespace=namespace)
 
      ! find out how many lines (i.e. states) the block has
      nlines = parse_block_n(blk)
 
      external_waves%number = nlines
      safe_allocate(external_waves%modus(1:nlines))
      safe_allocate(external_waves%e_field_string(1:3, 1:nlines))
      safe_allocate(external_waves%e_field(1:3, 1:nlines))
      safe_allocate(external_waves%k_vector(1:3, 1:nlines))
      safe_allocate(external_waves%v_vector(1:3, 1:nlines))
      safe_allocate(external_waves%mx_function(1:nlines))
      safe_allocate(external_waves%field_type(1:nlines))
      safe_allocate(external_waves%pw_phase(1:nlines))
      external_waves%pw_phase = m_zero
 
      call external_waves%bessel%init(nlines, 3)
 
      do il = 1, nlines
        ncols = parse_block_cols(blk, il - 1)
        if ((ncols < 5) .or. (ncols > 9)) then
          message(1) = 'Each line in the MaxwellIncidentWaves block must have five to nine columns.'
          call messages_fatal(1, namespace=namespace)
        end if
 
        ! check input modus e.g. parser of defined functions
        call parse_block_integer(blk, il - 1, 0, external_waves%modus(il))
        call parse_block_integer(blk, il - 1, 1, external_waves%field_type(il))
 
        ! parse formula string
        if (external_waves%modus(il) == option__maxwellincidentwaves__plane_wave_parser) then
          do idim = 1, 3
            call parse_block_string( blk, il - 1, idim + 1, k_string(idim))
            call parse_block_string( blk, il - 1, 3 + idim + 1, external_waves%e_field_string(idim, il))
          end do
          write(message(1), '(a,i2,a) ') 'Substituting electromagnetic incident wave ', il, ' with the expressions: '
          call messages_info(1, namespace=namespace)
          do idim = 1, 3
            write(message(idim), '(6a)')     '  Wave vector k('//dims(idim)//')   = ', trim(k_string(idim))
            write(message(idim+1), '(2a)')     '  E-field('//dims(idim)//') for t_0 = ', trim(external_waves%e_field_string(idim, il))
          end do
          call messages_info(6, namespace=namespace)
 
          do idim = 1, 3
            call conv_to_c_string(k_string(idim))
            call conv_to_c_string(external_waves%e_field_string(idim, il))
          end do
 
          x_pos(:) = m_zero
          x_norm    = m_zero
          do idim = 1, 3
            call parse_expression(k_vector(idim), dummy(idim), idim, x_pos, x_norm, m_zero, k_string(idim))
          end do
 
          k_vector = units_to_atomic(unit_one/units_inp%length, k_vector)
          k_norm = norm2(k_vector)
 
          velocity(:)    = k_vector(:) / k_norm * p_c * external_waves%c_factor
          external_waves%k_vector(:,il) = k_vector(:)
          external_waves%v_vector(:,il) = velocity(:)
 
        else if (external_waves%modus(il) == option__maxwellincidentwaves__plane_wave_mx_function) then
          do idim = 1, 3
            call parse_block_cmplx( blk, il - 1, idim + 1, e_field(idim))
          end do
          call parse_block_string( blk, il - 1, 3 + 2, mxf_expression)
 
          write(message(1), '(a,i2) ') 'Substituting electromagnetic incident wave ', il
          write(message(2), '(a)'    ) 'with the expression: '
          call messages_info(2, namespace=namespace)
 
          do idim = 1, 3
            write(message(idim), '(a,f9.4,sp,f9.4,"i")') '  E-field('//trim(dims(idim))//') complex amplitude  = ', &
              real(e_field(idim)), aimag(e_field(idim))
          end do
          write(message(4), '(2a)'    )      '  Maxwell wave function name = ', trim(mxf_expression)
          call messages_info(4, namespace=namespace)
 
          call mxf_read(external_waves%mx_function(il), namespace, trim(mxf_expression), iex_norm)
          if (iex_norm /= 0) then
            write(message(1),'(3A)') 'Ex_norm in the ""', trim(mxf_expression), &
              '"" field defined in the MaxwellIncidentWaves block'
            call messages_fatal(1, namespace=namespace)
          end if
          if (parse_block_cols(blk, il-1) == 7) then
            call parse_block_float( blk, il - 1, 3 + 3 , external_waves%pw_phase(il))
          end if
          e_field  = units_to_atomic(units_inp%energy/units_inp%length, e_field)
          k_vector(1:3) = external_waves%mx_function(il)%k_vector(1:3)
          k_norm = norm2(k_vector)
 
          k_dot_e = real(dot_product(k_vector, e_field), real64)
          if (abs(k_dot_e) > test_limit) then
            write(message(1), '(a) ') 'The wave vector k or its electric field are not perpendicular. '
            write(message(2), '(a,f8.3,a)' ) 'Their dot product yields the magnitude', abs(k_dot_e) , ' while'
            write(message(3), '(a,f8.3,a)' ) 'tolerance is ', test_limit ,'.'
            call messages_fatal(3, namespace=namespace)
          end if
          if (k_norm < 1e-10_real64) then
            message(1) = 'The k vector is not set correctly: |k|~0 .'
            call messages_fatal(1, namespace=namespace)
          end if
 
          external_waves%e_field(:,il)  = e_field(:)
          external_waves%k_vector(:,il) = k_vector(:)
          external_waves%v_vector(:,il) = k_vector(:) / k_norm * p_c * external_waves%c_factor
 
        else if (external_waves%modus(il) == option__maxwellincidentwaves__bessel_function) then
          call parse_block_float( blk, il - 1, 2 , external_waves%bessel%amp(il))
          call parse_block_integer( blk, il - 1, 3 , external_waves%bessel%m_order(il))
          call parse_block_float( blk, il - 1, 4 , external_waves%bessel%omega(il))
          call parse_block_integer( blk, il - 1, 5 , external_waves%bessel%helicity(il))
          call parse_block_float( blk, il - 1, 6 , external_waves%bessel%theta_k(il))
          if (parse_block_cols(blk, il-1) > 7) then
            call parse_block_string( blk, il - 1, 7 , mxf_expression)
            external_waves%bessel%envelope(il) = .true.
            call mxf_read(external_waves%mx_function(il), namespace, trim(mxf_expression), iex_norm)
          end if
          if (parse_block_cols(blk, il-1) == 9) then
            call parse_block_integer( blk, il - 1, 8 , external_waves%bessel%lin_dir(il))
          end if
 
          write(message(1), '(a,i2) ') 'Incident Bessel Beam', il
          call messages_info(1, namespace=namespace)
 
          if (abs(external_waves%bessel%helicity(il)) /= 1) then
            write(message(1),'(A)') 'Helicity has to be either +1 or -1 !'
            call messages_fatal(1, namespace=namespace)
          end if
 
          write(message(1), '(a,f5.3)'    )  ' Bessel Amplitude ', external_waves%bessel%amp(il)
          write(message(2), '(a,i2)'    )  ' Bessel Order m', external_waves%bessel%m_order(il)
          write(message(3), '(a,f5.3)'    )  ' Bessel Frequency ', external_waves%bessel%omega(il)
          write(message(4), '(a,i2)'    )  ' Bessel Helicity ', external_waves%bessel%helicity(il)
          write(message(5), '(a,f5.3)'    )  ' Bessel Opening Angle ', external_waves%bessel%theta_k(il)
          call messages_info(4, namespace=namespace)
 
          if (external_waves%bessel%lin_dir(il)/= 0) then
            write(message(5), '(a,i2)'    )  ' Bessel is Linearly Polarized in Direction : ', external_waves%bessel%lin_dir(il)
            call messages_info(4, namespace=namespace)
          end if
 
        end if
      end do
 
      call parse_block_end(blk)
 
      call messages_print_with_emphasis(namespace=namespace)
    else
      external_waves%number = 0
 
    end if
 
    !%Variable BesselBeamAxisShift
    !%Type block
    !%Section Maxwell
    !%Description
    !% The BesselBeamAxisShift block allows to shift the Bessel Beam, which is centered at (0,0,0) as default.
    !% Selected position point will be used as the new center of the Bessel Beam.
    !% When defining a BesselBeamAxisShift, please make sure to define a shift for each Bessel source you use,
    !% then it is possible to tell which source is shifted according to which BesselShift, respectively.
    !% <tt>%BesselBeamAxisShift
    !% <br>&nbsp;&nbsp;    0.0 | 2.0 |  5.0
    !% <br>%</tt>
    !%
    !%End
 
    if (parse_block(namespace, 'BesselBeamAxisShift', blk) == 0) then
      nlines = parse_block_n(blk)
      ncols = parse_block_cols(blk,0)
      if (ncols /= 3 ) then
        message(1) = 'BesselBeamAxisShift must have 3 columns.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      do il = 1, nlines
        do idim = 1, 3
          call parse_block_float(blk, il-1, idim-1, external_waves%bessel%shift(il, idim), units_inp%length)
        end do
      end do
 
      call parse_block_end(blk)
    end if
 
    call profiling_out('EXTERNAL_WAVES_INIT')
 
    pop_sub(external_waves_init)
  end subroutine external_waves_init
 
  ! ---------------------------------------------------------
  subroutine external_waves_end(external_waves)
    type(external_waves_t),   intent(inout) :: external_waves
 
    push_sub(external_waves_end)
 
    if (external_waves%output_from_point) then
      call io_close(external_waves%out_file)
      safe_deallocate_a(external_waves%selected_point_coordinate)
      safe_deallocate_a(external_waves%selected_point_field)
    end if
 
    safe_deallocate_a(external_waves%bessel%shift)
    safe_deallocate_a(external_waves%points_map)
    safe_deallocate_a(external_waves%modus)
    safe_deallocate_a(external_waves%e_field_string)
    safe_deallocate_a(external_waves%k_vector)
    safe_deallocate_a(external_waves%v_vector)
    safe_deallocate_a(external_waves%e_field)
    safe_deallocate_a(external_waves%mx_function)
    safe_deallocate_a(external_waves%pw_phase)
 
    if (accel_is_enabled()) then
      call accel_release_buffer(external_waves%buff_map)
    end if
 
    pop_sub(external_waves_end)
  end subroutine external_waves_end
 
  ! ---------------------------------------------------------
  subroutine external_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
    type(external_waves_t),    intent(inout) :: external_waves
    real(real64),              intent(in)    :: time
    class(mesh_t),             intent(in)    :: mesh
    character(len=*),          intent(in)    :: type_of_field
    real(real64),              intent(out)   :: out_field_total(:, :)
    type(derivatives_t), optional, intent(in):: der
 
 
    push_sub(external_waves_eval)
 
    call profiling_in('EXTERNAL_WAVES_EVAL')
 
    out_field_total = m_zero
 
    call plane_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der=der)
    call bessel_source_eval(external_waves, time, mesh, type_of_field, out_field_total, der=der)
 
    call profiling_out('EXTERNAL_WAVES_EVAL')
 
    pop_sub(external_waves_eval)
  end subroutine external_waves_eval
 
  ! ---------------------------------------------------------
  subroutine plane_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
    type(external_waves_t),    intent(inout) :: external_waves
    real(real64),              intent(in)    :: time
    class(mesh_t),             intent(in)    :: mesh
    character(len=*),          intent(in)    :: type_of_field
    real(real64),              intent(out)   :: out_field_total(:, :)
    type(derivatives_t), optional, intent(in):: der
 
    integer              :: wn
    real(real64), allocatable   :: pw_field(:,:), ztmp(:,:), b_field_aux(:,:)
    real(real64)         :: p_c_
    integer, allocatable :: indices_pw_parser(:)
    integer, allocatable :: indices_mx_ftc(:)
    integer :: n_plane_waves, n_points
 
    push_sub(plane_waves_eval)
 
    call profiling_in('PLANE_WAVES_EVAL')
 
    indices_pw_parser = pack([(wn, wn = 1,external_waves%number)], &
      external_waves%modus == option__maxwellincidentwaves__plane_wave_parser)
 
    indices_mx_ftc = pack([(wn, wn = 1,external_waves%number)], &
      external_waves%modus == option__maxwellincidentwaves__plane_wave_mx_function)
 
    n_plane_waves = size(indices_pw_parser) + size(indices_mx_ftc)
 
    p_c_ = p_c * external_waves%c_factor
 
    if (n_plane_waves == 0) then
      call profiling_out('PLANE_WAVES_EVAL')
      pop_sub(plane_waves_eval)
      return
    end if
 
    if (type_of_field == "B field" .and. any(external_waves%field_type == e_field_vector_potential)) then
      assert(present(der))
      safe_allocate(ztmp(mesh%np, size(out_field_total, dim=2)))
      n_points = mesh%np_part
    else
      n_points = mesh%np
    end if
    safe_allocate(pw_field(n_points, size(out_field_total, dim=2)))
    pw_field(:,:) = m_zero
 
    ! The E_field (or A_field, rescaled later) we calculate always
    do wn = 1, external_waves%number
 
      select case(external_waves%modus(wn))
      case (option__maxwellincidentwaves__plane_wave_parser)
        call pw_parsed_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
 
      case (option__maxwellincidentwaves__plane_wave_mx_function)
        call pw_mx_function_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
      end select
 
      select case (external_waves%field_type(wn))
 
      case(e_field_electric)
 
        select case (type_of_field)
        case ("E field")
          out_field_total(1:mesh%np,:) = out_field_total(1:mesh%np,:) + pw_field(1:mesh%np,:)
        case ("vector potential")
          call messages_not_implemented("Calculation of a vector potential from a plane wave specified as electric field")
        case ("B field")
          safe_allocate(b_field_aux(1:mesh%np, 1:mesh%box%dim))
          call get_pw_b_field(external_waves, mesh, wn, pw_field, b_field_aux)
          out_field_total(:,:) = out_field_total(:,:) + b_field_aux(:,:)
          safe_deallocate_a(b_field_aux)
        end select
 
      case(e_field_vector_potential)
 
        select case (type_of_field)
        case ("E field")
          call messages_not_implemented("Calculation of an electric field from a plane wave specified as vector potential")
        case ("vector potential")
          out_field_total(1:mesh%np,:) = out_field_total(1:mesh%np,:) - m_one/p_c_ * pw_field(1:mesh%np,1:3)
        case ("B field")
          call dderivatives_curl(der, pw_field(1:mesh%np_part,1:3), ztmp(1:mesh%np,1:3), set_bc = .false.)
          out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c_ * ztmp(1:mesh%np, 1:3)
        end select
 
      end select
    end do
 
    safe_deallocate_a(pw_field)
    safe_deallocate_a(ztmp)
    call profiling_out('PLANE_WAVES_EVAL')
 
    pop_sub(plane_waves_eval)
 
  end subroutine plane_waves_eval
 
  ! ---------------------------------------------------------
  subroutine pw_parsed_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
    type(external_waves_t),    intent(inout) :: external_waves
    integer,                   intent(in)    :: wn
    real(real64),              intent(in)    :: time
    class(mesh_t),             intent(in)    :: mesh
    integer,                   intent(in)    :: n_points
    real(real64),              intent(out)   :: pw_field(:,:)
 
    real(real64)         :: x_prop(3), x_norm
    real(real64)         :: velocity_time(3)
    real(real64)         :: parsed_field(3)
    real(real64)         :: dummy(3)
    integer              :: idim, ip
 
    velocity_time(:) = external_waves%v_vector(1:3, wn) * time
    do idim = 1, 3
      do ip = 1, n_points
        x_prop = mesh%x(ip, :) - velocity_time
        x_norm = norm2(x_prop(1:3))
        call parse_expression(parsed_field(idim), dummy(idim), 3, x_prop, x_norm, m_zero, &
          external_waves%e_field_string(idim, wn))
        pw_field(ip, idim) = units_to_atomic(units_inp%energy/units_inp%length, parsed_field(idim))
      end do
    end do
 
  end subroutine pw_parsed_evaluation
 
  ! ---------------------------------------------------------
  subroutine pw_mx_function_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
    type(external_waves_t),    intent(inout) :: external_waves
    integer,                   intent(in)    :: wn
    real(real64),              intent(in)    :: time
    class(mesh_t),             intent(in)    :: mesh
    integer,                   intent(in)    :: n_points
    real(real64),              intent(out)   :: pw_field(:,:)
 
    real(real64)         :: x_prop(3), x_norm
    real(real64)         :: velocity_time(3)
    complex(real64)      :: efield_ip(3)
    complex(real64)      :: e0(3)
    integer              :: ip
 
    velocity_time(:) = external_waves%v_vector(1:3, wn) * time
    e0(:) = external_waves%e_field(1:3, wn)
    do ip = 1, n_points
      x_prop = mesh%x(ip, :) - velocity_time
      x_norm = norm2(x_prop(1:3))
      efield_ip = mxf(external_waves%mx_function(wn), x_prop, external_waves%pw_phase(wn))
      pw_field(ip, :) = real(e0(1:3) * efield_ip, real64)
    end do
 
  end subroutine pw_mx_function_evaluation
 
  ! ---------------------------------------------------------
  subroutine get_pw_b_field(external_waves, mesh, pwidx, e_field, b_field)
    type(external_waves_t),  intent(in)  :: external_waves
    class(mesh_t),           intent(in)  :: mesh
    real(real64),            intent(in)  :: e_field(:,:)
    real(real64),            intent(out) :: b_field(:,:)
    integer,                 intent(in)  :: pwidx
 
    real(real64)   :: k_vector(3), k_vector_abs
    real(real64)   :: velocity(3)
    real(real64)   :: P_c_
    complex(real64)   :: e0(3)
    integer :: ip
 
    velocity = external_waves%v_vector(1:3, pwidx)
    k_vector = external_waves%k_vector(1:3, pwidx)
    k_vector_abs = norm2(k_vector(1:3))
    e0 = external_waves%e_field(1:3, pwidx)
    p_c_ = p_c * external_waves%c_factor
 
    b_field = m_zero
    do ip = 1, mesh%np
      b_field(ip, :) = m_one/(p_c_ * k_vector_abs) * dcross_product(k_vector, e_field(ip, :))
    end do
 
  end subroutine get_pw_b_field
 
  ! ---------------------------------------------------------
  subroutine bessel_source_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
    type(external_waves_t),    intent(inout) :: external_waves
    real(real64),              intent(in)    :: time
    class(mesh_t),             intent(in)    :: mesh
    character(len=*),          intent(in)    :: type_of_field
    real(real64),              intent(out)   :: out_field_total(:, :)
    type(derivatives_t), optional, intent(in):: der
 
    real(real64)         :: dmin, omega, k_vector(3), c_factor
    integer              :: iline, wn, pos_index, n_points, rankmin
    real(real64), allocatable   :: shift(:,:)
    complex(real64), allocatable   :: bessel_field_total(:,:), ztmp(:,:), vec_pot(:,:)
    integer, allocatable :: indices_bessel_ftc(:)
    type(mxf_t)          :: envelope_mxf
 
    push_sub(bessel_source_eval)
 
    call profiling_in('BESSEL_SOURCE_EVAL')
 
    indices_bessel_ftc = pack([(wn, wn = 1,external_waves%number)], &
      external_waves%modus == option__maxwellincidentwaves__bessel_function)
 
    if (size(indices_bessel_ftc) == 0) then
      call profiling_out('BESSEL_SOURCE_EVAL')
      pop_sub(bessel_source_eval)
      return
    end if
 
    ! Check if the BesselBeamAxisShift is defined for every incoming Bessel Beam.
    if (allocated(external_waves%bessel%shift) .and. &
      size(external_waves%bessel%shift(:,1)) /= size(indices_bessel_ftc)) then
      message(1) = 'Number of BesselBeamAxisShift defined in input file'
      message(2) = 'does not match the number of Bessel beams.'
      call messages_fatal(2)
    end if
 
    safe_allocate(shift(size(indices_bessel_ftc), 3))
    if (allocated(external_waves%bessel%shift)) then
      shift = external_waves%bessel%shift
    else
      shift = m_zero
    end if
 
    if (type_of_field == "B field") then
      assert(present(der))
      safe_allocate(vec_pot(mesh%np_part, size(out_field_total, dim=2)))
      safe_allocate(ztmp(size(out_field_total, dim=1), size(out_field_total, dim=2)))
      n_points = mesh%np_part ! needed for curl
    else
      n_points = mesh%np
    end if
 
 
    safe_allocate(bessel_field_total(1:n_points, 1:3))
    bessel_field_total = m_zero
 
    do iline = 1, size(indices_bessel_ftc)
      wn = indices_bessel_ftc(iline)
      omega = external_waves%bessel%omega(wn)
      k_vector = external_waves%mx_function(wn)%k_vector
      c_factor =  external_waves%c_factor
      envelope_mxf = external_waves%mx_function(wn)
 
      call external_waves%bessel%function(wn, shift, mesh, n_points, time, k_vector, c_factor, envelope_mxf, bessel_field_total)
 
      select case (external_waves%field_type(wn))
 
      case(e_field_vector_potential)
        ! interpreting bessel_field_total as a vector potential (as requested by the user)
        select case (type_of_field)
        case ("E field")
          out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) + real(m_zi*omega*bessel_field_total(1:mesh%np,1:3))
        case ("vector potential")
          ! For the vector potential, we multiply by -1/c becuase of the electronic Hamiltonian
          ! being in Gaussian units
          out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c * real(bessel_field_total(1:mesh%np,1:3))
        case ("B field")
          call zderivatives_curl(der, bessel_field_total(1:mesh%np_part,1:3), ztmp(1:mesh%np,1:3), set_bc = .false.)
          out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c * real(ztmp(1:mesh%np, 1:3))
        end select
 
      case(e_field_electric)
        ! interpreting bessel_field_total as an electric field (as requested by the user)
        select case (type_of_field)
        case ("E field")
          out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) + real(bessel_field_total(1:mesh%np,1:3))
        case ("vector potential")
          ! We calculate the vector potential as real(E/i*omega),
          ! and convert it to the proper units by multiplying by -1/c
          out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c * &
            real(bessel_field_total(1:mesh%np,1:3)/M_zI/omega)
        case ("B field")
          vec_pot(1:mesh%np_part,1:3) = - m_one/p_c * real(bessel_field_total(1:mesh%np_part,1:3)/m_zi/omega)
          call zderivatives_curl(der, vec_pot(1:mesh%np_part,1:3), ztmp(1:mesh%np,1:3), set_bc = .false.)
          out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - real(ztmp(1:mesh%np, 1:3))
        end select
 
      end select
    end do
 
    if (external_waves%output_from_point) then
      pos_index = mesh_nearest_point(mesh, external_waves%selected_point_coordinate(1:3), dmin, rankmin)
      if (mesh%mpi_grp%rank == rankmin) then
        external_waves%selected_point_field(:) = out_field_total(pos_index,:)
        write(external_waves%out_file, "(4F14.8, 4x)") time, external_waves%selected_point_field(:)
      end if
    end if
 
    safe_deallocate_a(shift)
    safe_deallocate_a(ztmp)
    safe_deallocate_a(vec_pot)
    safe_deallocate_a(bessel_field_total)
    call profiling_out('BESSEL_SOURCE_EVAL')
 
    pop_sub(bessel_source_eval)
 
  end subroutine bessel_source_eval
 
  ! ---------------------------------------------------------
  subroutine bessel_beam_function(this, iline, shift, mesh, n_points, time, k_vector, c_factor, envelope_mxf, bessel_field)
    class(bessel_beam_t), intent(in)   :: this
    integer,              intent(in)   :: iline
    real(real64),         intent(in)   :: shift(:,:), time, k_vector(3), c_factor
    class(mesh_t),        intent(in)   :: mesh
    integer,              intent(in)   :: n_points
    type(mxf_t),          intent(in)   :: envelope_mxf
    complex(real64),      intent(out)  :: bessel_field(:,:)
 
    real(real64)   :: pos(3), temp, temp2, temp3, rho, phi_rho, wigner(3)
    real(real64)   :: hel, theta, omega, amp, kappa, proj, k_norm, velocity_time(3), x_prop(3)
    complex(real64)   :: efield_ip(3)
    real(real64)   :: bessel_plus, bessel_minus
    integer :: ip, mm, pol
 
    assert(iline <= size(this%omega))
    hel = real(this%helicity(iline), real64)
    theta = this%theta_k(iline)
    mm = this%m_order(iline)
    amp = this%amp(iline) / sqrt(m_two)
    omega = this%omega(iline)
    proj = omega * cos(theta) / p_c             ! k_z
    kappa = sqrt(omega**2 - (proj* p_c)**2)     ! parse omega
    ! Set Wigner Coefficients from theta
    wigner(1) = hel * sin(theta) / sqrt(m_two)  ! mu = 0
    wigner(2) = 0.5 * (1 + hel * cos(theta))    ! mu = 1
    wigner(3) = 0.5 * (1 - hel * cos(theta))    ! mu = -1
    proj = omega * cos(theta) / p_c             ! k_z
    pol = this%lin_dir(iline)  ! Incoming polarization corresponding to beam in question
 
    do ip = 1, n_points
      pos(:) = mesh%x(ip, :) -  shift(iline,:)
      rho = norm2(pos(1:2))
      phi_rho = atan2(pos(2) , pos(1))
      temp = proj * pos(3) + phi_rho * (mm + 1) - omega*time ! temp, temp2 and temp3 should be unitless
      temp2 = proj * pos(3) + phi_rho * (mm - 1) - omega*time
      temp3 = proj * pos(3) + phi_rho * mm - omega*time
      bessel_plus = loct_bessel(mm+1, kappa * rho / p_c)
      bessel_minus = loct_bessel(mm-1, kappa * rho / p_c)
 
      ! Calculate complex Ax component, if generalized bessel OR x -polarized bessel
      if (pol /= 2) then
        bessel_field(ip, 1) = amp * (exp(m_zi*temp) * wigner(3) * bessel_plus + exp(m_zi*temp2) * wigner(2) * bessel_minus)
      end if
      ! Calculate complex Ay component if generalized bessel OR y -polarized bessel
      if (pol/=1) then
        bessel_field(ip, 2) = m_zi * amp * (-exp(m_zi*temp) * wigner(3) * bessel_plus + &
          exp(m_zi*temp2) * wigner(2) * bessel_minus)
      end if
      ! Calculate complex Az component, only iff generalized Bessel
      if (pol == 0) then
        bessel_field(ip, 3) = - m_zi * amp * sqrt(m_two) * wigner(1) * loct_bessel(mm, kappa * rho / p_c) * exp(m_zi*temp3)
      end if
 
      if (this%envelope(iline)) then
        k_norm = norm2(k_vector)
        velocity_time = k_vector * p_c * c_factor * time  / k_norm
        x_prop(:) = pos(:) - velocity_time(:)
        efield_ip = mxf_envelope_eval(envelope_mxf, x_prop)
        bessel_field(ip, :) = bessel_field(ip, :) * real(efield_ip, real64)
      end if
 
    end do
 
  end subroutine bessel_beam_function
 
  subroutine bessel_beam_init(this, nlines, dim)
    class(bessel_beam_t), intent(out) :: this
    integer, intent(in) :: nlines
    integer, intent(in) :: dim
 
    safe_allocate(this%amp(1: nlines))
    safe_allocate(this%omega(1:nlines))
    safe_allocate(this%theta_k(1:nlines))
    safe_allocate(this%m_order(1:nlines))
    safe_allocate(this%helicity(1:nlines))
    safe_allocate(this%shift(1:nlines, 1:dim))
    safe_allocate(this%envelope(1:nlines))
    safe_allocate(this%lin_dir(1:nlines))
    this%amp = m_zero
    this%omega = m_zero
    this%theta_k = m_zero
    this%m_order = m_zero
    this%helicity = m_zero
    this%shift = m_zero
    this%lin_dir = m_zero
    this%envelope = .false.
 
  end subroutine bessel_beam_init
 
  subroutine bessel_beam_finalize(this)
    type(bessel_beam_t), intent(inout) :: this
 
    safe_deallocate_a(this%amp)
    safe_deallocate_a(this%omega)
    safe_deallocate_a(this%theta_k)
    safe_deallocate_a(this%m_order)
    safe_deallocate_a(this%helicity)
    safe_deallocate_a(this%shift)
    safe_deallocate_a(this%lin_dir)
    safe_deallocate_a(this%envelope)
 
  end subroutine bessel_beam_finalize
 
end module external_waves_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: