! =============================================================================
!> This submodule defines a cylindrical equilibrium in which a Suydam surface
!! is present, such that this gives rises to Suydam cluster modes.
!! The geometry can be overriden using the parfile.
!!
!! This equilibrium is taken from
!! _Nijboer, R. J., Holst, B., Poedts, S., & Goedbloed, J. P. (1997).
!! Calculating magnetohydrodynamic flow spectra. Computer physics communications, 106(1-2), 39-52_.
!! @note Default values are given by
!!
!! - k2 = 1
!! - k3 = -1.2
!! - cte_rho0 = 1 : used to set the density value.
!! - cte_p0 = 0.05 : used to set the background pressure.
!! - cte_v02 = 0 : sets a constant flow \(v_\theta\).
!! - cte_v03 = 0.14 : prefactor for flow profile \(v_z\).
!! - p1 = 0.1 : constant that appears in magnetic field and pressure components.
!! - alpha = 2 : used in the Bessel functions.
!!
!! and can all be changed in the parfile. @endnote
submodule (mod_equilibrium) smod_equil_suydam_cluster
use mod_equilibrium_params, only: cte_rho0, cte_v02, cte_v03, cte_p0, p1, alpha
implicit none
contains
!> Sets the equilibrium
module procedure suydam_cluster_eq
if (settings%equilibrium%use_defaults) then ! LCOV_EXCL_START
call settings%grid%set_geometry("cylindrical")
call settings%grid%set_grid_boundaries(0.0_dp, 1.0_dp)
call settings%physics%enable_flow()
cte_rho0 = 1.0_dp
cte_v02 = 0.0_dp
cte_v03 = 0.14_dp
cte_p0 = 0.05_dp
p1 = 0.1_dp
alpha = 2.0_dp
k2 = 1.0_dp
k3 = -1.2_dp
end if ! LCOV_EXCL_STOP
call background%set_density_funcs(rho0_func=rho0)
call background%set_velocity_2_funcs(v02_func=v02)
call background%set_velocity_3_funcs(v03_func=v03, dv03_func=dv03)
call background%set_temperature_funcs(T0_func=T0, dT0_func=dT0)
call background%set_magnetic_2_funcs(B02_func=B02, dB02_func=dB02)
call background%set_magnetic_3_funcs(B03_func=B03, dB03_func=dB03)
end procedure suydam_cluster_eq
real(dp) function rho0()
rho0 = cte_rho0
end function rho0
real(dp) function T0(r)
real(dp), intent(in) :: r
T0 = (cte_p0 + 0.5_dp * p1 * J0(r)**2) / cte_rho0
end function T0
real(dp) function dT0(r)
real(dp), intent(in) :: r
dT0 = p1 * J0(r) * DJ0(r) / cte_rho0
end function dT0
real(dp) function v02()
v02 = cte_v02
end function v02
real(dp) function v03(r)
real(dp), intent(in) :: r
v03 = cte_v03 * (1.0_dp - r**2)
end function v03
real(dp) function dv03(r)
real(dp), intent(in) :: r
dv03 = -2.0_dp * cte_v03 * r
end function dv03
real(dp) function B02(r)
real(dp), intent(in) :: r
B02 = J1(r)
end function B02
real(dp) function dB02(r)
real(dp), intent(in) :: r
dB02 = DJ1(r)
end function dB02
real(dp) function B03(r)
real(dp), intent(in) :: r
B03 = sqrt(1.0_dp - p1) * J0(r)
end function B03
real(dp) function dB03(r)
real(dp), intent(in) :: r
dB03 = -alpha * sqrt(1.0_dp - p1) * J1(r)
end function dB03
real(dp) function J0(r)
real(dp), intent(in) :: r
J0 = bessel_jn(0, alpha * r)
end function J0
real(dp) function J1(r)
real(dp), intent(in) :: r
J1 = bessel_jn(1, alpha * r)
end function J1
real(dp) function DJ0(r)
real(dp), intent(in) :: r
DJ0 = -alpha * J1(r)
end function DJ0
real(dp) function DJ1(r)
real(dp), intent(in) :: r
DJ1 = alpha * (0.5_dp * J0(r) - 0.5_dp * bessel_jn(2, alpha * r))
end function DJ1
end submodule smod_equil_suydam_cluster