module classes_dlc_visitor

use, intrinsic :: iso_fortran_env, only: DP => REAL64
use data_constants, only: num_dimensions, PI
use classes_periodic_box, only: Abstract_Periodic_Box
use classes_reciprocal_lattice, only: Abstract_Reciprocal_Lattice
use types_particle_wrapper, only: Concrete_Particle
use classes_dlc_weight, only: Abstract_DLC_Weight
use classes_dlc_structures, only: Abstract_DLC_Structures
use classes_structure_visitor, only: Abstract_Structure_Visitor
use procedures_dipolar_interactions_micro, only: set_fourier, set_exp_kz, reci_number_1_sym

implicit none

private

    type, extends(Abstract_Structure_Visitor), abstract, public :: Abstract_DLC_Visitor
    private
        class(Abstract_Periodic_Box), pointer :: periodic_box => null()
        integer :: reci_numbers(num_dimensions) = 0
        class(Abstract_DLC_Weight), pointer :: weight => null()
        class(Abstract_DLC_Structures), pointer :: structures => null()
    contains
        procedure :: construct => Abstract_construct
        procedure :: destroy => Abstract_destroy
        procedure :: target => Abstract_target
        procedure :: visit => Abstract_visit
        procedure :: visit_translation => Abstract_visit_translation
        procedure :: visit_transmutation => Abstract_visit_transmutation
        procedure :: visit_rotation => Abstract_visit_rotation
        procedure :: visit_add => Abstract_visit_add
        procedure :: visit_remove => Abstract_visit_remove
        procedure :: visit_switch => Abstract_visit_switch
        procedure, private :: visit_exchange => Abstract_visit_exchange
    end type Abstract_DLC_Visitor

    type, extends(Abstract_DLC_Visitor), public :: Concrete_DLC_Visitor

    end type Concrete_DLC_Visitor

    type, extends(Abstract_DLC_Visitor), public :: Null_DLC_Visitor
    contains
        procedure :: construct => Null_construct
        procedure :: destroy => Null_destroy
        procedure :: target => Null_target
        procedure :: visit => Null_visit
        procedure :: visit_translation => Null_visit_translation
        procedure :: visit_transmutation => Null_visit_transmutation
        procedure :: visit_switch => Null_visit_switch
        procedure, private :: visit_exchange => Null_visit_exchange
    end type Null_DLC_Visitor

contains

!implementation Abstract_DLC_Visitor

    subroutine Abstract_construct(this, periodic_box, reciprocal_lattice, weight, structures)
        class(Abstract_DLC_Visitor), intent(out) :: this
        class(Abstract_Periodic_Box), target, intent(in) :: periodic_box
        class(Abstract_Reciprocal_Lattice), intent(in) :: reciprocal_lattice
        class(Abstract_DLC_Weight), intent(in) :: weight
        class(Abstract_DLC_Structures), intent(in) :: structures

        this%periodic_box => periodic_box
        this%reci_numbers = reciprocal_lattice%get_numbers()
        call this%target(weight, structures)
    end subroutine Abstract_construct

    subroutine Abstract_destroy(this)
        class(Abstract_DLC_Visitor), intent(inout) :: this

        this%structures => null()
        this%weight => null()
        this%periodic_box => null()
    end subroutine Abstract_destroy

    subroutine Abstract_target(this, weight, structures)
        class(Abstract_DLC_Visitor), intent(inout) :: this
        class(Abstract_DLC_Weight), target, intent(in) :: weight
        class(Abstract_DLC_Structures), target, intent(in) :: structures

        this%weight => weight
        this%structures => structures
    end subroutine Abstract_target

    !> \[
    !>      U = \sum_{\vec{k}_{1:2}} w(\vec{k}_{1:2})
    !>          \Re[S_+(\vec{k}_{1:2}) S_-(\vec{k}_{1:2})^\ast]
    !> \]
    !> where \( w(\vec{k}_{1:2}) \) is [[Abstract_DLC_Weight]] and
    !> \( S_\pm(\vec{k}_{1:2}) \) is [[Abstract_DLC_Structures]].
    pure real(DP) function Abstract_visit(this) result(energy)
        class(Abstract_DLC_Visitor), intent(in) :: this

        integer :: n_1, n_2

        energy = 0._DP
        do n_2 = 0, this%reci_numbers(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle
                energy = energy + this%weight%get(n_1, n_2) * &
                    real(this%structures%get_plus(n_1, n_2) * &
                        conjg(this%structures%get_minus(n_1, n_2)), DP)
            end do
        end do
        energy = 2._DP * energy ! symmetry: half wave vectors -> double energy
    end function Abstract_visit

    !> Energy delta when a particle is translated:
    !> \( (\vec{x}, \vec{\mu}) \to (\vec{x}^\prime, \vec{\mu}) \).
    !> \[
    !>      \Delta S_-^\ast(\vec{k}_{1:2}) =
    !>          (-k_{1:2} \mu_3 - i \vec{k}_{1:2} \cdot \vec{\mu}_{1:2})
    !>          \left(
    !>              e^{-k_{1:2} x^\prime_3} e^{-i \vec{k}_{1:2} \cdot \vec{x}^\prime_{1:2}} -
    !>              e^{-k_{1:2} x_3} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !>          \right)
    !> \]
    !> \[
    !>      \Delta S_+(\vec{k}_{1:2}) = (+k_{1:2} \mu_3 + i \vec{k}_{1:2} \cdot \vec{\mu}_{1:2})
    !>          \left(
    !>              e^{+k_{1:2} x^\prime_3} e^{+i \vec{k}_{1:2} \cdot \vec{x}^\prime_{1:2}} -
    !>              e^{+k_{1:2} x_3} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !>          \right)
    !> \]
    !> \[
    !>      \Delta S_+(\vec{k}_{1:2}) \Delta S_-^\ast(\vec{k}_{1:2}) =
    !>          \left[
    !>              -(k_{1:2} \mu_3)^2 + (\vec{k}_{1:2} \cdot \vec{\mu}_{1:2})^2 -
    !>              2 i k_{1:2} \mu_3 \vec{k}_{1:2} \cdot \vec{\mu}_{1:2}
    !>          \right]
    !>          \left(
    !>              2 - \\
    !>              e^{+k_{1:2} x^\prime_3} e^{+i \vec{k}_{1:2} \cdot \vec{x}^\prime_{1:2}}
    !>              e^{-k_{1:2} x_3} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{1:2}} - \\
    !>              e^{-k_{1:2} x^\prime_3} e^{-i \vec{k}_{1:2} \cdot \vec{x}^\prime_{1:2}}
    !>              e^{+k_{1:2} x_3} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !>          \right)
    !> \]
    pure real(DP) function Abstract_visit_translation(this, i_component, new_position, old) &
        result(delta_energy)
        class(Abstract_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: i_component
        real(DP), intent(in) :: new_position(:)
        type(Concrete_Particle), intent(in) :: old

        real(DP) :: real_part_1, real_part_2, real_part_3

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position_old, wave_1_x_position_new, wave_vector
        real(DP) :: wave_dot_moment_12, wave_x_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position_new
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_new_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_new_2
        real(DP) :: exp_kz_new
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_new_tab
        complex(DP) :: fourier_position_old
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_old_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_old_2
        real(DP) :: exp_kz_old
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_old_tab

        delta_energy = 0._DP
        if (.not.this%structures%is_dipolar(i_component)) return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position_old = 2._DP*PI * old%position(1:2) / surface_size
        call set_fourier(fourier_position_old_1, this%reci_numbers(1), wave_1_x_position_old(1))
        call set_fourier(fourier_position_old_2, this%reci_numbers(2), wave_1_x_position_old(2))
        call set_exp_kz(exp_kz_old_tab, surface_size, old%position(3))
        wave_1_x_position_new = 2._DP*PI * new_position(1:2) / surface_size
        call set_fourier(fourier_position_new_1, this%reci_numbers(1), wave_1_x_position_new(1))
        call set_fourier(fourier_position_new_2, this%reci_numbers(2), wave_1_x_position_new(2))
        call set_exp_kz(exp_kz_new_tab, surface_size, new_position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position_old = fourier_position_old_1(n_1) * fourier_position_old_2(n_2)
                exp_kz_old = exp_kz_old_tab(abs(n_1), abs(n_2))
                fourier_position_new = fourier_position_new_1(n_1) * fourier_position_new_2(n_2)
                exp_kz_new = exp_kz_new_tab(abs(n_1), abs(n_2))

                wave_dot_moment_12 = dot_product(wave_vector, old%dipole_moment(1:2))
                wave_x_moment_3 = norm2(wave_vector) * old%dipole_moment(3)

                real_part_1 = real(this%structures%get_plus(n_1, n_2) * &
                    cmplx(-wave_x_moment_3, -wave_dot_moment_12, DP) * &
                    conjg(fourier_position_new/exp_kz_new - fourier_position_old/exp_kz_old), DP)
                real_part_2 = real(conjg(this%structures%get_minus(n_1, n_2)) * &
                    cmplx(+wave_x_moment_3, +wave_dot_moment_12, DP) * &
                    (fourier_position_new * exp_kz_new - fourier_position_old * exp_kz_old), DP)
                real_part_3 = real(cmplx(-wave_x_moment_3**2 + wave_dot_moment_12**2, &
                    -2._DP*wave_x_moment_3 * wave_dot_moment_12, DP) * (2._DP - &
                    fourier_position_new*exp_kz_new * conjg(fourier_position_old)/exp_kz_old - &
                    fourier_position_old*exp_kz_old * conjg(fourier_position_new)/exp_kz_new), DP)

                delta_energy = delta_energy + this%weight%get(n_1, n_2) * &
                    (real_part_1 + real_part_2 + real_part_3)
            end do
        end do
        delta_energy = 2._DP * delta_energy ! symmetry: half wave vectors -> double energy
    end function Abstract_visit_translation

    !> Energy delta when a particle is transmuted:
    !> \( (\vec{x}, \vec{\mu}) \to (\vec{x}, \vec{\mu}^\prime) \).
    !> \[
    !>      \Delta S_-^\ast(\vec{k}_{1:2}) =
    !>          \left[
    !>              -k_{1:2} (\mu^\prime_3 - \mu_3) -
    !>              i \vec{k}_{1:2} \cdot (\vec{\mu}^\prime_{1:2} - \vec{\mu}_{1:2})
    !>          \right]
    !>          e^{-k_{1:2} x_3} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !> \]
    !> \[
    !>      \Delta S_+(\vec{k}_{1:2}) =
    !>          \left[
    !>              +k_{1:2} (\mu^\prime_3 - \mu_3) +
    !>              i \vec{k}_{1:2} \cdot (\vec{\mu}^\prime_{1:2} - \vec{\mu}_{1:2})
    !>          \right]
    !>          e^{+k_{1:2} x_3} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !> \]
    !> \[
    !>      \Re[\Delta S_+(\vec{k}_{1:2}) \Delta S_-^\ast(\vec{k}_{1:2})] =
    !>          -[k_{1:2} (\mu^\prime_3 - \mu_3)]^2 +
    !>          [\vec{k}_{1:2} \cdot (\vec{\mu}^\prime_{1:2} - \vec{\mu}_{1:2})]^2
    !> \]
    pure real(DP) function Abstract_visit_transmutation(this, ij_components, new_dipole_moment, &
        old) result(delta_energy)
        class(Abstract_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: ij_components(:)
        real(DP), intent(in) :: new_dipole_moment(:)
        type(Concrete_Particle), intent(in) :: old

        real(DP) :: real_part_1, real_part_2, real_part_3

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position, wave_vector
        real(DP) :: wave_dot_delta_moment_12, wave_delta_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_2
        real(DP) :: exp_kz
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_tab

        delta_energy = 0._DP
        if (.not.(this%structures%is_dipolar(ij_components(1)) .or. &
            this%structures%is_dipolar(ij_components(2)))) return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position = 2._DP*PI * old%position(1:2) / surface_size
        call set_fourier(fourier_position_1, this%reci_numbers(1), wave_1_x_position(1))
        call set_fourier(fourier_position_2, this%reci_numbers(2), wave_1_x_position(2))
        call set_exp_kz(exp_kz_tab, surface_size, old%position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position = fourier_position_1(n_1) * fourier_position_2(n_2)
                exp_kz = exp_kz_tab(abs(n_1), abs(n_2))

                wave_dot_delta_moment_12 = dot_product(wave_vector, new_dipole_moment(1:2) - old%&
                    dipole_moment(1:2))
                wave_delta_moment_3 = norm2(wave_vector) * (new_dipole_moment(3) - old%&
                    dipole_moment(3))

                real_part_1 = real(this%structures%get_plus(n_1, n_2) * &
                    cmplx(-wave_delta_moment_3, -wave_dot_delta_moment_12, DP) * &
                    conjg(fourier_position) / exp_kz, DP)
                real_part_2 = real(conjg(this%structures%get_minus(n_1, n_2)) * &
                    cmplx(+wave_delta_moment_3, +wave_dot_delta_moment_12, DP) * &
                    fourier_position * exp_kz, DP)
                real_part_3 = -wave_delta_moment_3**2 + wave_dot_delta_moment_12**2

                delta_energy = delta_energy + this%weight%get(n_1, n_2) * &
                    (real_part_1 + real_part_2 + real_part_3)
            end do
        end do
        delta_energy = 2._DP * delta_energy ! symmetry: half wave vectors -> double energy
    end function Abstract_visit_transmutation

    pure real(DP) function Abstract_visit_rotation(this, i_component, new_dipole_moment, old) &
        result(delta_energy)
        class(Abstract_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: i_component
        real(DP), intent(in) :: new_dipole_moment(:)
        type(Concrete_Particle), intent(in) :: old

        delta_energy = this%visit_transmutation([i_component, i_component], new_dipole_moment, old)
    end function Abstract_visit_rotation

    pure real(DP) function Abstract_visit_add(this, i_component, particle) result(delta_energy)
        class(Abstract_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle

        delta_energy = this%visit_exchange(i_component, particle, +1._DP)
    end function Abstract_visit_add

    pure real(DP) function Abstract_visit_remove(this, i_component, particle) result(delta_energy)
        class(Abstract_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle

        delta_energy = this%visit_exchange(i_component, particle, -1._DP)
    end function Abstract_visit_remove

    !> Energy delta when a particle of coordinates \( (\vec{x}, \vec{\mu}) \) is added
    !> (\( \pmb{+} \)) or removed (\( \pmb{-} \)).
    !> \[
    !>      \Delta S_-^\ast(\vec{k}_{1:2}) =
    !>          \pmb{\pm} (-k_{1:2} \mu_3 - i \vec{k}_{1:2} \cdot \vec{\mu}_{1:2})
    !>          e^{-k_{1:2} x_3} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !> \]
    !> \[
    !>      \Delta S_+(\vec{k}_{1:2}) =
    !>          \pmb{\pm} (+k_{1:2} \mu_3 + i \vec{k}_{1:2} \cdot \vec{\mu}_{1:2})
    !>          e^{+k_{1:2} x_3} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !> \]
    !> \[
    !>      \Re[\Delta S_+(\vec{k}_{1:2}) \Delta S_-^\ast(\vec{k}_{1:2})] =
    !>          -(k_{1:2} \mu_3)^2 + (\vec{k}_{1:2} \cdot \vec{\mu}_{1:2})^2
    !> \]
    pure real(DP) function Abstract_visit_exchange(this, i_component, particle, signed) &
        result(delta_energy)
        class(Abstract_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle
        real(DP), intent(in) :: signed

        real(DP) :: real_part_1, real_part_2, real_part_3

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position, wave_vector
        real(DP) :: wave_dot_moment_12, wave_x_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_2
        real(DP) :: exp_kz
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_tab

        delta_energy = 0._DP
        if (.not.this%structures%is_dipolar(i_component)) return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position = 2._DP*PI * particle%position(1:2) / surface_size
        call set_fourier(fourier_position_1, this%reci_numbers(1), wave_1_x_position(1))
        call set_fourier(fourier_position_2, this%reci_numbers(2), wave_1_x_position(2))
        call set_exp_kz(exp_kz_tab, surface_size, particle%position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position = fourier_position_1(n_1) * fourier_position_2(n_2)
                exp_kz = exp_kz_tab(abs(n_1), abs(n_2))
                wave_dot_moment_12 = dot_product(wave_vector, particle%dipole_moment(1:2))
                wave_x_moment_3 = norm2(wave_vector) * particle%dipole_moment(3)

                real_part_1 = signed * real(this%structures%get_plus(n_1, n_2) * &
                    cmplx(-wave_x_moment_3, -wave_dot_moment_12, DP) * &
                    conjg(fourier_position) / exp_kz, DP)
                real_part_2 = signed * real(conjg(this%structures%get_minus(n_1, n_2)) * &
                    cmplx(+wave_x_moment_3, +wave_dot_moment_12, DP) * &
                    fourier_position * exp_kz, DP)
                real_part_3 = -wave_x_moment_3**2 + wave_dot_moment_12**2

                delta_energy = delta_energy + this%weight%get(n_1, n_2) * &
                    (real_part_1 + real_part_2 + real_part_3)
            end do
        end do
        delta_energy = 2._DP * delta_energy ! symmetry: half wave vectors -> double energy
    end function Abstract_visit_exchange

    !> Energy delta when 2 particles of coordinates \( (\vec{x}_1, \vec{\mu}_1) \) and
    !> \( (\vec{x}_2, \vec{\mu}_2) \) switch their positions.
    !> \[
    !>      \Delta S_-^\ast(\vec{k}_{1:2}) =
    !>          \left[
    !>              -k_{1:2} (\mu_{1, 3} - \mu_{2, 3}) -
    !>              i \vec{k}_{1:2} \cdot (\vec{\mu}_{1, 1:2} - \vec{\mu}_{2, 1:2})
    !>          \right] \\
    !>          \left(
    !>              e^{-k_{1:2} x_{2, 3}} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{2, 1:2}} -
    !>              e^{-k_{1:2} x_{1, 3}} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{1, 1:2}}
    !>          \right)
    !> \]
    !> \[
    !>      \Delta S_+(\vec{k}_{1:2}) =
    !>          \left[
    !>              +k_{1:2} (\mu_{1, 3} - \mu_{2, 3}) +
    !>              i \vec{k}_{1:2} \cdot (\vec{\mu}_{1, 1:2} - \vec{\mu}_{2, 1:2})
    !>          \right] \\
    !>          \left(
    !>              e^{+k_{1:2} x_{2, 3}} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{2, 1:2}} -
    !>              e^{+k_{1:2} x_{1, 3}} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{1, 1:2}}
    !>          \right)
    !> \]
    !> \[
    !>      \Delta S_+(\vec{k}_{1:2}) \Delta S_-^\ast(\vec{k}_{1:2}) =
    !>          \left\{
    !>              -[k_{1:2} (\mu_{1, 3} - \mu_{2, 3})]^2 +
    !>              [\vec{k}_{1:2} \cdot (\vec{\mu}_{1, 1:2} - \vec{\mu}_{2, 1:2})]^2 - \\
    !>              2 i k_{1:2} (\mu_{1, 3} - \mu_{2, 3})
    !>              \vec{k}_{1:2} \cdot (\vec{\mu}_{1, 1:2} - \vec{\mu}_{2, 1:2})
    !>          \right\}
    !>          \left(
    !>              2 - \\
    !>              e^{+k_{1:2} x_{2, 3}} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{2, 1:2}}
    !>              e^{-k_{1:2} x_{1, 3}} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{1, 1:2}} - \\
    !>              e^{+k_{1:2} x_{1, 3}} e^{+i \vec{k}_{1:2} \cdot \vec{x}_{1, 1:2}}
    !>              e^{-k_{1:2} x_{2, 3}} e^{-i \vec{k}_{1:2} \cdot \vec{x}_{2, 1:2}}
    !>          \right)
    !> \]
    pure real(DP) function Abstract_visit_switch(this, ij_components, particles) &
        result(delta_energy)
        class(Abstract_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: ij_components(:)
        type(Concrete_Particle), intent(in) :: particles(:)

        real(DP) :: real_part_1, real_part_2, real_part_3

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position_1, wave_1_x_position_2, wave_vector
        real(DP) :: wave_dot_delta_moment_12, wave_delta_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position_1
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_1_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_1_2
        real(DP) :: exp_kz_1
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_1_tab
        complex(DP) :: fourier_position_2
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_2_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_2_2
        real(DP) :: exp_kz_2
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_2_tab

        delta_energy = 0._DP
        if (.not.(this%structures%is_dipolar(ij_components(1)) .or. &
            this%structures%is_dipolar(ij_components(2)))) return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position_1 = 2._DP*PI * particles(1)%position(1:2) / surface_size
        call set_fourier(fourier_position_1_1, this%reci_numbers(1), wave_1_x_position_1(1))
        call set_fourier(fourier_position_1_2, this%reci_numbers(2), wave_1_x_position_1(2))
        call set_exp_kz(exp_kz_1_tab, surface_size, particles(1)%position(3))
        wave_1_x_position_2 = 2._DP*PI * particles(2)%position(1:2) / surface_size
        call set_fourier(fourier_position_2_1, this%reci_numbers(1), wave_1_x_position_2(1))
        call set_fourier(fourier_position_2_2, this%reci_numbers(2), wave_1_x_position_2(2))
        call set_exp_kz(exp_kz_2_tab, surface_size, particles(2)%position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position_1 = fourier_position_1_1(n_1) * fourier_position_1_2(n_2)
                exp_kz_1 = exp_kz_1_tab(abs(n_1), abs(n_2))
                fourier_position_2 = fourier_position_2_1(n_1) * fourier_position_2_2(n_2)
                exp_kz_2 = exp_kz_2_tab(abs(n_1), abs(n_2))

                wave_dot_delta_moment_12 = dot_product(wave_vector, &
                    particles(1)%dipole_moment(1:2) - particles(2)%dipole_moment(1:2))
                wave_delta_moment_3 = norm2(wave_vector) * &
                    (particles(1)%dipole_moment(3) - particles(2)%dipole_moment(3))

                real_part_1 = real(this%structures%get_plus(n_1, n_2) * &
                    cmplx(-wave_delta_moment_3, -wave_dot_delta_moment_12, DP) * &
                    conjg(fourier_position_2 / exp_kz_2 - fourier_position_1 / exp_kz_1), DP)
                real_part_2 = real(conjg(this%structures%get_minus(n_1, n_2)) * &
                    cmplx(+wave_delta_moment_3, +wave_dot_delta_moment_12, DP) * &
                    (fourier_position_2 * exp_kz_2 - fourier_position_1 * exp_kz_1), DP)
                real_part_3 = real(cmplx(-wave_delta_moment_3**2 + wave_dot_delta_moment_12**2, &
                    -2._DP*wave_delta_moment_3 * wave_dot_delta_moment_12, DP) * (2._DP - &
                    fourier_position_2 * exp_kz_2 * conjg(fourier_position_1) / exp_kz_1 - &
                    fourier_position_1 * exp_kz_1 * conjg(fourier_position_2) / exp_kz_2), DP)

                delta_energy = delta_energy + this%weight%get(n_1, n_2) * &
                    (real_part_1 + real_part_2 + real_part_3)
            end do
        end do
        delta_energy = 2._DP * delta_energy ! symmetry: half wave vectors -> double energy
    end function Abstract_visit_switch

!end implementation Abstract_DLC_Visitor

!implementation Null_DLC_Visitor

    subroutine Null_construct(this, periodic_box, reciprocal_lattice, weight, structures)
        class(Null_DLC_Visitor), intent(out) :: this
        class(Abstract_Periodic_Box), target, intent(in) :: periodic_box
        class(Abstract_Reciprocal_Lattice), intent(in) :: reciprocal_lattice
        class(Abstract_DLC_Weight), intent(in) :: weight
        class(Abstract_DLC_Structures), intent(in) :: structures
    end subroutine Null_construct

    subroutine Null_destroy(this)
        class(Null_DLC_Visitor), intent(inout) :: this
    end subroutine Null_destroy

    subroutine Null_target(this, weight, structures)
        class(Null_DLC_Visitor), intent(inout) :: this
        class(Abstract_DLC_Weight), target, intent(in) :: weight
        class(Abstract_DLC_Structures), target, intent(in) :: structures
    end subroutine Null_target

    pure real(DP) function Null_visit(this) result(energy)
        class(Null_DLC_Visitor), intent(in) :: this
        energy = 0._DP
    end function Null_visit

    pure real(DP) function Null_visit_translation(this, i_component, new_position, old) &
        result(delta_energy)
        class(Null_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: i_component
        real(DP), intent(in) :: new_position(:)
        type(Concrete_Particle), intent(in) :: old
        delta_energy = 0._DP
    end function Null_visit_translation

    pure real(DP) function Null_visit_transmutation(this, ij_components, new_dipole_moment, old) &
        result(delta_energy)
        class(Null_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: ij_components(:)
        real(DP), intent(in) :: new_dipole_moment(:)
        type(Concrete_Particle), intent(in) :: old
        delta_energy = 0._DP
    end function Null_visit_transmutation

    pure real(DP) function Null_visit_exchange(this, i_component, particle, signed) &
        result(delta_energy)
        class(Null_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle
        real(DP), intent(in) :: signed
        delta_energy = 0._DP
    end function Null_visit_exchange

    pure real(DP) function Null_visit_switch(this, ij_components, particles) result(delta_energy)
        class(Null_DLC_Visitor), intent(in) :: this
        integer, intent(in) :: ij_components(:)
        type(Concrete_Particle), intent(in) :: particles(:)
        delta_energy = 0._DP
    end function Null_visit_switch

!end implementation Null_DLC_Visitor

end module classes_dlc_visitor