module classes_des_self_component

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_permittivity, only: Abstract_Permittivity
use classes_component_dipole_moments, only: Abstract_Component_Dipole_Moments
use classes_des_convergence_parameter, only: Abstract_DES_Convergence_Parameter

implicit none

private

    type, abstract, public :: Abstract_DES_Self_Component
    private
        class(Abstract_Periodic_Box), pointer :: periodic_box => null()
        class(Abstract_Component_Dipole_Moments), pointer :: dipole_moments => null()
        real(DP) :: factor_denominator = 0._DP
        real(DP) :: alpha_x_box_edge = 0._DP
    contains
        procedure :: construct => Abstract_construct
        procedure :: destroy => Abstract_destroy
        procedure :: visit => Abstract_visit
        procedure :: meet => Abstract_meet
    end type Abstract_DES_Self_Component

    type, extends(Abstract_DES_Self_Component), public :: Concrete_DES_Self_Component

    end type Concrete_DES_Self_Component

    type, extends(Abstract_DES_Self_Component), public :: Null_DES_Self_Component
    contains
        procedure :: construct => Null_construct
        procedure :: destroy => Null_destroy
        procedure :: visit => Null_visit
        procedure :: meet => Null_meet
    end type Null_DES_Self_Component

    type, public :: DES_Self_Component_Wrapper
        class(Abstract_DES_Self_Component), allocatable :: component
    end type DES_Self_Component_Wrapper

contains

!implementation Abstract_DES_Self_Component

    subroutine Abstract_construct(this, periodic_box, permittivity, dipole_moments, alpha)
        class(Abstract_DES_Self_Component), intent(out) :: this
        class(Abstract_Periodic_Box), target, intent(in) :: periodic_box
        class(Abstract_Permittivity), intent(in) :: permittivity
        class(Abstract_Component_Dipole_Moments), target, intent(in) :: dipole_moments
        class(Abstract_DES_Convergence_Parameter), intent(in) :: alpha

        this%periodic_box => periodic_box
        this%dipole_moments => dipole_moments
        this%factor_denominator = 6._DP * permittivity%get() * PI**(3._DP/2._DP)
        this%alpha_x_box_edge = alpha%get_times_box_edge()
    end subroutine Abstract_construct

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

        this%dipole_moments => null()
        this%periodic_box => null()
    end subroutine Abstract_destroy

    !> \[ \frac{\alpha^3}{6\epsilon\pi^{3/2}} \sum_{\vec{\mu}} \vec{\mu}\cdot\vec{\mu} \]
    pure real(DP) function Abstract_visit(this) result(energy)
        class(Abstract_DES_Self_Component), intent(in) :: this

        real(DP) :: box_size(num_dimensions)
        real(DP) :: moment_squared_sum
        integer :: i_particle

        moment_squared_sum = 0._DP
        do i_particle = 1, this%dipole_moments%get_num()
            moment_squared_sum = moment_squared_sum + dot_product(this%dipole_moments%&
                get(i_particle), this%dipole_moments%get(i_particle))
        end do
        box_size = this%periodic_box%get_size()
        energy = (this%alpha_x_box_edge/box_size(1))**3 / this%factor_denominator * &
            moment_squared_sum
    end function Abstract_visit

    !> \[ u(\mu) = \frac{\alpha^3}{6\epsilon\pi^{3/2}} \vec{\mu}\cdot\vec{\mu} \]
    pure real(DP) function Abstract_meet(this, dipole_moment) result(energy)
        class(Abstract_DES_Self_Component), intent(in) :: this
        real(DP), intent(in) :: dipole_moment(:)

        real(DP) :: box_size(num_dimensions)

        box_size = this%periodic_box%get_size()
        energy = (this%alpha_x_box_edge/box_size(1))**3 / this%factor_denominator * &
            dot_product(dipole_moment, dipole_moment)
    end function Abstract_meet

!end implementation Abstract_DES_Self_Component

!implementation Null_DES_Self_Component

    subroutine Null_construct(this, periodic_box, permittivity, dipole_moments, alpha)
        class(Null_DES_Self_Component), intent(out) :: this
        class(Abstract_Periodic_Box), target, intent(in) :: periodic_box
        class(Abstract_Permittivity), intent(in) :: permittivity
        class(Abstract_Component_Dipole_Moments), target, intent(in) :: dipole_moments
        class(Abstract_DES_Convergence_Parameter), intent(in) :: alpha
    end subroutine Null_construct

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

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

    pure real(DP) function Null_meet(this, dipole_moment) result(energy)
        class(Null_DES_Self_Component), intent(in) :: this
        real(DP), intent(in) :: dipole_moment(:)
        energy = 0._DP
    end function Null_meet

!end implementation Null_DES_Self_Component

end module classes_des_self_component