module mod_matrix_coeffs
    use mod_parameters
    use mod_variables
    implicit none
    type :: sph_coeffs
        real(8) :: mass_bal_rhojm, mass_bal_vjm1, mass_bal_vjp1
        real(8) :: const_p_vjm1(1:2), const_p_vjp1(1:2), const_p_sj0
        real(8) :: Mom__jm1_vjm1, Mom__jm1_sj0(1:2), Mom__jm1_sjm2(1:2), Mom__jm1_sj2(1:2), Mom__jm1_rho
        real(8) :: Mom__jp1_vjp1, Mom__jp1_sj0(1:2), Mom__jp1_sjp2(1:2), Mom__jp1_sj2(1:2), Mom__jp1_rho
        real(8) :: Const_jm2_vjm1(1:2)
        real(8) :: Const_j2__vjm1(1:2), Const_j2__vjp1(1:2)
        real(8) :: Const_jp2_vjp1(1:2)
    end type sph_coeffs

    type :: sph_boundary_coeffs
        real(8) :: Bound_jm1_s_0, Bound_jm1_s_j2, Bound_jm1_s_jm2
        real(8) :: Bound_jm1_vjm1, Bound_jm1_vjp1
        real(8) :: Bound_jp1_s_0, Bound_jp1_s_j2, Bound_jp1_s_jp2
        real(8) :: Bound_jp1_vjm1, Bound_jp1_vjp1
    end type sph_boundary_coeffs

    type :: tor_coeffs
        real(8) :: Mom__j_v, Mom__j_s_jm1(1:2), Mom__j_s_jp1(1:2)
        real(8) :: Const_jm1_vj(1:2), Const_jp1_vj(1:2)
    end type tor_coeffs

    type :: tor_boundary_coeffs
        real(8) :: Bound_j_s_jm1, Bound_j_s_jp1, Bound_j_vj
    end type tor_boundary_coeffs

    type(sph_coeffs)   :: sc
    type(sph_boundary_coeffs) :: sc_b
    type(tor_coeffs)     :: tc
    type(tor_boundary_coeffs) :: tc_b
contains
    

subroutine matrix_coeffs(j)
    implicit none
    real(8) :: j
    
    call spheriodal_coeffs(j)
    call spheriodal_boundary_coeffs(j)
    call toroidal_coeffs(j)
    call toroidal_boundary_coefss(j)
end subroutine matrix_coeffs

subroutine spheriodal_coeffs(j)
    implicit none
    real(8) :: j
    !! Balance of mass rho_t = -div(rho0 v)
    !! Below the coefficients of rho_t and vr
    sc%mass_bal_rhojm = 1.0d0 / time_step
    sc%mass_bal_vjm1  =  dsqrt(  j        / (2.0d0 * j + 1.0d0))
    sc%mass_bal_vjp1  = -dsqrt( (j+1.0d0) / (2.0d0 * j + 1.0d0))

    !! Constitutive - pressure if solid   p_t - vr rho0 g0 = -K div(v)
    !!                         if liquid  p   = -K * rho/rho0 + zeta * div(v)
    !!                            last term in the liquid case may be negligible
    !! below coefficients are just give for div(v)
    sc%const_p_vjm1(1) =  dsqrt(  j        / (2.0d0 * j + 1.0d0))
    sc%const_p_vjp1(1) = -dsqrt( (j+1.0d0) / (2.0d0 * j + 1.0d0))
    sc%const_p_vjm1(2) = -(j-1.0d0)*sc%const_p_vjm1(1)
    sc%const_p_vjp1(2) =  (j+2.0d0)*sc%const_p_vjp1(1)
    !! Balance of momentum rho0 v_t = div(sigma) - rho * g *er - rho0(nabla_V) - 2rho0 omega x v
    !!                         |         implicit            |  |          explicit/RHS         |
    !! below just the coefficients of implicit part
    ! !Y_{JM}^{J-1}
    sc%Mom__jm1_sj0(1)  = -dsqrt( j / (3.0d0 * (2.0d0 * j + 1.0d0)))
    sc%Mom__jm1_sjm2(1) =  dsqrt((j - 1.0d0) / (2.0d0 * j - 1.0d0))
    sc%Mom__jm1_sj2(1)  = -dsqrt( (j + 1.0d0)*(2.0d0 * j + 3.0d0) /(6.d0*(2.d0*j - 1.d0)*(2.d0*j + 1.d0) ))
    sc%Mom__jm1_sj0(2)  =  sc%Mom__jm1_sj0(1)  * (j + 1.0d0)
    sc%Mom__jm1_sjm2(2) = -sc%Mom__jm1_sjm2(1) * (j - 2.0d0)
    sc%Mom__jm1_sj2(2)  =  sc%Mom__jm1_sj2(1)  * (j + 1.0d0)
    sc%Mom__jm1_rho     = -dsqrt( j       /(2.0d0 * j + 1.0d0))

    ! !Y_{JM}^{J+1}
    sc%Mom__jp1_sj0(1)   =  dsqrt((j + 1.0d0)/((3.d0*(2.d0*j + 1.d0))))
    sc%Mom__jp1_sjp2(1)  = -dsqrt((j +2.d0) / (2.d0*j + 3.d0))
    sc%Mom__jp1_sj2(1)   =  dsqrt(( j*(2.d0*j - 1.d0) )/(6.d0*(2.d0*j + 3.d0)*(2.d0*j + 1.d0)))
    sc%Mom__jp1_sj0(2)   = -sc%Mom__jp1_sj0(1)  * j
    sc%Mom__jp1_sjp2(2)  =  sc%Mom__jp1_sjp2(1) * (j + 3.0d0)
    sc%Mom__jp1_sj2(2)   = -sc%Mom__jp1_sj2(1)  * j
    sc%Mom__jp1_rho      = -(-dsqrt((j+1.0d0)/(2.0d0 * j + 1.0d0)))

    !!Constitutive - Deviatoric (sigma^d) / eta + (sigma^d)_t / mu  = (nabla_v + nabla_v^T)
    ! ! Y_{JM}^{J-2,2}
    sc%Const_jm2_vjm1(1) = -dsqrt((j - 1.0d0) / (2.0d0 * j - 1.0d0))
    sc%Const_jm2_vjm1(2) =  sc%Const_jm2_vjm1(1) * j

    ! ! Y_{JM}^{J  ,2}
    sc%Const_j2__vjm1(1) =  dsqrt( ((j + 1.d0)*(2.d0*j + 3.d0))/(6.d0*(2.d0*j - 1.d0)*(2.d0*j + 1.d0)))
    sc%Const_j2__vjp1(1) = -dsqrt((j*(2.d0*j - 1.d0)) / (6.d0*(2.d0*j + 3.d0)*(2.d0*j + 1.d0)))
    sc%Const_j2__vjm1(2) = -sc%Const_j2__vjm1(1) * (j - 1.0d0)
    sc%Const_j2__vjp1(2) =  sc%Const_j2__vjp1(1) * (j + 2.0d0)

    ! ! Y_{JM}^{J+2,2}
    sc%Const_jp2_vjp1(1) =  dsqrt( (j + 2.0d0)  / (2.0d0 * j + 3.0d0) )
    sc%Const_jp2_vjp1(2) = -(j + 1.0d0) * sc%Const_jp2_vjp1(1)
end subroutine spheriodal_coeffs

subroutine spheriodal_boundary_coeffs(j)
    implicit none
    real(8) :: j
    ! !BOUNDARY CONDITIONS Y_{JM}^{J-1}
    sc_b%Bound_jm1_s_0   = -dsqrt(j / (3.d0*(2.d0*j +1.d0)) )
    sc_b%Bound_jm1_s_j2  = -dsqrt( ((j + 1.d0)*(2.d0*j + 3.d0)) / (6.d0*(2.d0*j + 1.d0)*(2.d0*j -1.0d0)) )
    sc_b%Bound_jm1_s_jm2 =  dsqrt( (j - 1.d0)/(2.d0*j -1.d0) )

    sc_b%Bound_jm1_vjm1 =                        (j / (2.0d0 * j + 1.0d0))
    sc_b%Bound_jm1_vjp1 = - dsqrt(j * (j + 1.0d0) ) / (2.0d0 * j + 1.0d0)
    ! !BOUNDARY CONDITIONS Y_{JM}^{J+1}
    sc_b%Bound_jp1_s_0   =  dsqrt( (j + 1.d0) / (3.d0*(2.d0*j + 1.d0)) )
    sc_b%Bound_jp1_s_j2  =  dsqrt( (  j*(2.d0*j - 1.d0)) / (6.d0*(2.d0*j + 1.d0)*(2.d0*j + 3.d0)) )
    sc_b%Bound_jp1_s_jp2 = -dsqrt( (  j + 2.d0) / (2.d0*j + 3.d0) )

    sc_b%Bound_jp1_vjm1 = - dsqrt(j * (j + 1.0d0) ) / (2.0d0 * j + 1.0d0)
    sc_b%Bound_jp1_vjp1 =             (j + 1.0d0)   / (2.0d0 * j + 1.0d0)
end subroutine spheriodal_boundary_coeffs

subroutine toroidal_coeffs(j)
    implicit none
    real(8) :: j

    tc%Mom__j_s_jm1(1) =  dsqrt((j - 1.0d0)/(2.0d0*(2.0d0 * j + 1.0d0)))
    tc%Mom__j_s_jm1(2) = -(j - 1.0d0) * tc%Mom__j_s_jm1(1)
    tc%Mom__j_s_jp1(1) = -dsqrt((j + 2.0d0)/(2.0d0*(2.0d0 * j + 1.0d0)))
    tc%Mom__j_s_jp1(2) =  (j + 2.0d0) * tc%Mom__j_s_jp1(1)


    tc%Const_jm1_vj(1) = -dsqrt((j - 1.0d0)/(2.0d0*(2.0d0 * j + 1.0d0)))
    tc%Const_jm1_vj(2) =  (j + 1.0d0) * tc%Const_jm1_vj(1)

    tc%Const_jp1_vj(1) =  dsqrt((j + 2.0d0)/(2.0d0*(2.0d0 * j + 1.0d0)))
    tc%Const_jp1_vj(2) = - j * tc%Const_jp1_vj(1)
end subroutine toroidal_coeffs

subroutine toroidal_boundary_coefss(j)
    implicit none
    real(8) :: j
    
    tc_b%Bound_j_s_jm1 =  dsqrt((j-1.d0)/(2.d0*(2.d0*j + 1.d0)))
    tc_b%Bound_j_s_jp1 = -dsqrt((j+2.d0)/(2.d0*(2.d0*j + 1.d0)))
end subroutine toroidal_boundary_coefss
end module mod_matrix_coeffs