! Copyright 2026 - The Minton Group at Purdue University
! This file is part of Swiftest.
! Swiftest 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 3 of the License, or (at your option) any later version.
! Swiftest 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 Swiftest. 
! If not, see: https://www.gnu.org/licenses. 

submodule(ringmoons) s_ringmoons_util
    use swiftest
contains

    module subroutine ringmoons_util_accrete_cb(self,ring,seed,param,dt)
        use, intrinsic :: ieee_exceptions
        implicit none

        ! Arguments
        class(ringmoons_cb),        intent(inout) :: self
            !! Ringmoons central body object
        class(ringmoons_ring), intent(inout) :: ring
            !! Ringmoons ring object
        class(ringmoons_seed), intent(inout) :: seed
            !! Ringmoons seed obje3ct
        class(swiftest_parameters), intent(in) :: param
            !! Current run configuration parameters
        real(DP),intent(in)                :: dt
            !! Current time step size

        ! Internals
        integer(I4B) :: i,j,iin
        real(DP) :: Lplanet, Lring, Ltot,Rnew,Mnew, Lorig,Mring,dGMtot,Lnow,rfac
        real(DP),dimension(seed%nbody) :: afac
        real(DP),dimension(0:ring%nbins+1)        :: mtmp,Lring_orig,Lring_now,dL
        real(DP) :: Lp0,Ls0,Lp1,Ls1,Lr0,Lr1,drot0,drot1
        logical, dimension(size(IEEE_ALL))      :: fpe_halting_modes

        ! Guard against underflow errors when rings surface mass density gets too small
        call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes)
        call ieee_set_halting_mode(ieee_underflow, .false.)
        
        associate(cb => self)
            ring%inside = ring%find_bin(cb%radius)
            dGMtot = sum(param%GU*ring%mass(0:ring%inside))
            if (dGMtot <= epsilon(1.0_DP)*cb%dGM) return
            !Add ring mass and angular momentum to planet
            Lring_orig(:) = ring%mass(:) * ring%Iz(:) * ring%nkep(:) 
            Lring = sum(Lring_orig(0:ring%inside))
            cb%dGM = cb%dGM + dGMtot 
            cb%dL(3) = Lring + cb%dL(3) 
            cb%Gmass = cb%dGM  + cb%GM0
            cb%mass = cb%Gmass / param%GU
            cb%radius = (3*cb%mass/(4*PI*cb%density))**THIRD
            cb%dR= cb%radius - cb%R0
            drot0 = cb%L0(3) * RAD2DEG / (cb%Ip(3) * cb%mass * cb%radius**2)  
            drot1 = cb%dL(3) * RAD2DEG / (cb%Ip(3) * cb%mass * cb%radius**2)
            cb%rot(3) = drot1 + drot0

            ring%mass(0:ring%inside) = 0.0_DP
            ring%sigma(0:ring%inside) = 0.0_DP

            if (seed%nbody > 0) then
                afac(1:seed%nbody) = 1._DP - dGMtot / seed%mu(1:seed%nbody)
                seed%a(1:seed%nbody) = seed%a(1:seed%nbody) * afac(1:seed%nbody)
            end if

            ! update body-dependent parameters as needed
            rfac = 1._DP - dGMtot / cb%Gmass
            ring%r_outer = ring%r_outer * rfac
            ring%r_inner = ring%r_inner * rfac

            ! Save the mass so that we can correct for the change in geometry
            mtmp(:) = ring%mass(:)
            call ring%reset(seed,cb,param)
            ring%mass(:) = mtmp(:)
            ring%sigma(:) = ring%mass(:) / ring%deltaA(:)
            ring%Gsigma(:) = param%GU * ring%sigma(:)
            
            ! Any difference in angular momentum in each ring bin will result in a torque in that bin
            Lring_now(:) = ring%mass(:) * ring%Iz(:) * ring%nkep(:) 
            dL(0:ring%inside) = 0.0_DP
            dL(ring%nbins+1) = 0.0_DP
            dL(ring%inside+1:ring%nbins) = (Lring_now(ring%inside+1:ring%nbins) - Lring_orig(ring%inside+1:ring%nbins)) / dt
            ring%Torque(:) = ring%Torque(:) - dL(:) 
        end associate

        call ieee_set_halting_mode(IEEE_ALL, fpe_halting_modes)
        return
    end subroutine ringmoons_util_accrete_cb

    module subroutine ringmoons_util_dealloc_ring(self)
        !! author: David A. Minton
        !!
        !! Deallocates all allocatabale arrays
        implicit none
        ! Arguments
        class(ringmoons_ring),  intent(inout) :: self 
            !! Ringmoons ring object

        self%nbins = 0
        if (allocated(self%r))       deallocate(self%r)
        if (allocated(self%X))       deallocate(self%X)
        if (allocated(self%X2))      deallocate(self%X2)
        if (allocated(self%r_hstar)) deallocate(self%r_hstar)
        if (allocated(self%deltaA))  deallocate(self%deltaA)
        if (allocated(self%sigma))   deallocate(self%sigma)
        if (allocated(self%Gsigma))  deallocate(self%Gsigma)
        if (allocated(self%mass))    deallocate(self%mass)
        if (allocated(self%tau))     deallocate(self%tau)
        if (allocated(self%nu))      deallocate(self%nu)
        if (allocated(self%Q))       deallocate(self%Q)
        if (allocated(self%Iz))      deallocate(self%Iz)
        if (allocated(self%nkep))    deallocate(self%nkep)
        if (allocated(self%Torque))  deallocate(self%Torque)
        if (allocated(self%r_p))     deallocate(self%r_p)
        if (allocated(self%m_p))     deallocate(self%m_p)
        if (allocated(self%rho_p))   deallocate(self%rho_p)
        if (allocated(self%vrel_p))  deallocate(self%vrel_p)
        if (allocated(self%Y_21))    deallocate(self%Y_21)
        if (allocated(self%delta))   deallocate(self%delta)

        return
    end subroutine ringmoons_util_dealloc_ring

    module subroutine ringmoons_util_dealloc_seed(self)
        !! author: David A. Minton
        !!
        !! Deallocates all allocatabale arrays
        implicit none
        ! Arguments
        class(ringmoons_seed),  intent(inout) :: self 
            !! Ringmoons ring object

        self%nbody = 0
        if (allocated(self%status))  deallocate(self%status)
        if (allocated(self%id))      deallocate(self%id)
        if (allocated(self%info))    deallocate(self%info)
        if (allocated(self%a))       deallocate(self%a)
        if (allocated(self%mu))      deallocate(self%mu)
        if (allocated(self%mass))    deallocate(self%mass)
        if (allocated(self%Gmass))   deallocate(self%Gmass)
        if (allocated(self%rhill))   deallocate(self%rhill)
        if (allocated(self%radius))  deallocate(self%radius)
        if (allocated(self%density)) deallocate(self%density)
        if (allocated(self%ringbin)) deallocate(self%ringbin)
        if (allocated(self%Torque))  deallocate(self%Torque)
        if (allocated(self%Ttide))   deallocate(self%Ttide)

        return
    end subroutine ringmoons_util_dealloc_seed

    module subroutine ringmoons_util_dealloc_system(self)
        !! author: David A. Minton
        !!
        !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer
        implicit none
        class(ringmoons_nbody_system), intent(inout) :: self
            !! Ringmoons nbody system object to deallocate
        if (allocated(self%ring)) deallocate(self%ring)
        if (allocated(self%seed)) deallocate(self%seed)
        call self%symba_nbody_system%dealloc()

        return
    end subroutine ringmoons_util_dealloc_system

    pure elemental module function ringmoons_util_find_bin(self,r) result(bin)
        !! author: David A. Minton
        !!
        !! Returns the bin containing radius r from the input ring.
        implicit none
        class(ringmoons_ring), intent(in)      :: self
            !! Ringmoons ring object
        real(DP), intent(in)                   :: r
            !! Radial distance at which to search for the bin
        integer(I4B)                           :: bin
            !! The bin containing radial distance r

        if (r > self%r_outer) then
            bin = self%nbins+1
        else if (r < self%r_inner) then
            bin = 0
        else
            bin = ceiling(2 * (sqrt(r) - sqrt(self%r_inner)) / self%deltaX) 
        end if
        
        return
    end function ringmoons_util_find_bin

    module function ringmoons_util_get_dt_ring(self,dtin) result(dtout)
        !! author: David A. Minton
        !!
        !! Calculates the maximum stable timestep for the surface mass density evolution that is not larger than dtin.
        !!
        !! Adapted from Andrew Hesselbrock's  RING-MOONS Python scripts
        implicit none

        ! Arguments
        class(ringmoons_ring), intent(in)  :: self
        real(DP), intent(in)               :: dtin
        real(DP)                           :: dtout

        ! Internals
        integer(I4B)                           :: i
        real(DP)                               :: sig_max,nu_max
        real(DP)                               :: torque_term
       
        associate(ring => self)
            ! Start with viscous stability
            dtout = dtin
    
            nu_max = max(maxval(ring%nu(:) / ring%X2(:)),0.0_DP)
    
            if (nu_max > 0.0_DP) then
                sig_max = 16 * (12._DP / (ring%deltaX)**2) * nu_max
                dtout = min(dtin,(sig_max)**(-1))
            end if
        end associate
        
        return
    end function ringmoons_util_get_dt_ring

    module subroutine ringmoons_util_update_ring(self,cb,param)
        !! author: David A. Minton
        !!
        !! Updates the ring velocity dispersion, Toomre parameter, and viscosity values
        implicit none
        class(ringmoons_ring), intent(inout) :: self
        class(ringmoons_cb), intent(in) :: cb
        class(swiftest_parameters), intent(in) :: param
        ! Internals
        integer(I4B) :: i

        associate(ring => self)       
            call ring%set_velocity_dispersion(cb,param)
            where(ring%sigma(:) > 1000 * VSMALL) 
                ring%Q(:) = ring%nkep(:) * ring%vrel_p(:) / (3.36_DP * ring%Gsigma(:))
                ring%tau(:) = PI * ring%r_p(:)**2 * ring%sigma(:) / ring%m_p(:)
                ring%nu(:) = ringmoons_viscosity(ring%Gsigma(:), ring%m_p(:), (ring%vrel_p(:))**2, &
                                                ring%r_p(:), ring%r_hstar(:), ring%Q(:), ring%tau(:), ring%nkep(:))
            elsewhere
                ring%Q(:) = huge(1._DP) / 10._DP
                ring%tau(:) = 0.0_DP
                ring%nu(:) = 0.0_DP
            end where
        end associate

        return
    end subroutine ringmoons_util_update_ring

    module subroutine ringmoons_util_reset_ring(self,seed,cb,param)
        !! author: David A. Minton
        !!
        !!  Resets ring torques and recomputes all dimensional quantities, such as ring extent and limits based on the current
        !! surface mass density and central body properties.
        use, intrinsic :: ieee_exceptions
        implicit none
        ! Arguments
        class(ringmoons_ring), intent(inout) :: self
        class(ringmoons_seed), intent(inout) :: seed
        class(ringmoons_cb), intent(in) :: cb
        class(swiftest_parameters), intent(in) :: param
        ! Internals
        integer(I4B)                        :: i, iFRL, iRRL
        real(DP)                            :: Xlo, rho_p
        logical, dimension(size(IEEE_ALL))      :: fpe_halting_modes

        ! Guard against underflow errors when rings surface mass density gets too small
        call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes)
        call ieee_set_halting_mode(ieee_underflow, .false.)
        associate(ring => self)
            ring%nu = 0.0_DP
            ring%X_inner = 2 * sqrt(ring%r_inner)
            ring%X_outer = 2 * sqrt(ring%r_outer)
            ring%deltaX = (ring%X_outer - ring%X_inner) / ring%nbins
            where (ring%r_p(:) > tiny(1.0_DP))
                ring%rho_p(:) = ring%m_p(:) / ((4.0_DP / 3.0_DP) * PI * ring%r_p(:)**3)
            end where
            ring%iFRL = ring%nbins / 2
            ring%iRRL = ring%nbins / 2
            do i = 1, 10
                if (ring%rho_p(ring%iFRL) > tiny(1.0_DP)) then
                    rho_p = ring%rho_p(ring%iFRL)
                else
                    rho_p = ring%rho_p(ring%nbins / 2)
                end if
                ring%FRL = 2.456_DP * cb%radius * (cb%density / rho_p)**(1._DP / 3._DP)

                ! end if
                if (ring%rho_p(ring%iRRL) > tiny(1.0_DP)) then
                    rho_p = ring%rho_p(ring%iRRL)
                else
                    rho_p = ring%rho_p(ring%nbins / 2)
                end if
                ring%RRL = 1.44_DP  * cb%radius * (cb%density / rho_p)**(1._DP / 3._DP)
                iFRL = ring%find_bin(ring%FRL)
                iRRL = ring%find_bin(ring%RRL)
                if ((iFRL == ring%iFRL) .and. (iRRL == ring%iRRL)) exit
                ring%iFRL = iFRL
                ring%iRRL = iRRL
            end do

            do i = 0,ring%nbins + 1
                ! Set up X coordinate system (see Bath & Pringle 1981)
                Xlo = ring%X_inner + ring%deltaX * (i - 1)
        
                ring%X(i) = Xlo + 0.5_DP * ring%deltaX
            end do
            ring%X2(:) = ring%X(:)**2
            ring%deltaA(:) = 0.25_DP * PI * ring%X(:)**3 * ring%deltaX 
            ring%r(:) = 0.25_DP * (ring%X(:))**2

            ! Specific moment of inertia of the ring bin and ring angular velocity
            ring%Iz(:) = (ring%r(:))**2
            ring%nkep(:) = sqrt(cb%Gmass / ring%r(:)**3)
            ring%mass(:) = ring%sigma(:) * ring%deltaA(:)
            ring%Gsigma(:) = param%GU * ring%sigma(:)
            
            ring%Torque(:) = 0.0_DP

            if (seed%nbody > 0) then
                seed%rhill(:) = seed%a(:) * (seed%mass(:) / cb%mass / 3)**THIRD
                seed%mu(:) = cb%Gmass + seed%Gmass(:)
                where (seed%status(:) == ACTIVE)
                    seed%ringbin(:)   = ring%find_bin(seed%a(:))
                elsewhere
                    seed%ringbin(:)   = 0
                end where
            end if

        end associate
        call ieee_set_halting_mode(IEEE_ALL, fpe_halting_modes)
        return
    end subroutine ringmoons_util_reset_ring

    module subroutine ringmoons_util_setup_ring(self, n, param)
        !! author: David A. Minton
        !!
        !! Constructor for the ringmoons_ring class. 
        !! Allocates space for all bins and initializes all components with a value. 
        implicit none
        class(ringmoons_ring),      intent(inout) :: self  
            !! Ringmoons ring object
        integer(I4B),               intent(in)    :: n     
            !! Number of bins to allocate space for
        class(swiftest_parameters), intent(in)    :: param 
            !! Current run configuration parameters

        if (n < 0) return

        call self%dealloc()
        self%nbins = n

        if (n == 0) return

        allocate(self%r(0:n+1))
        allocate(self%X(0:n+1))
        allocate(self%X2(0:n+1))
        allocate(self%r_hstar(0:n+1))
        allocate(self%deltaA(0:n+1))
        allocate(self%sigma(0:n+1))
        allocate(self%Gsigma(0:n+1))
        allocate(self%mass(0:n+1))
        allocate(self%tau(0:n+1))
        allocate(self%nu(0:n+1))
        allocate(self%Q(0:n+1))
        allocate(self%Iz(0:n+1))
        allocate(self%nkep(0:n+1))
        allocate(self%Torque(0:n+1))
        allocate(self%r_p(0:n+1))
        allocate(self%m_p(0:n+1))
        allocate(self%rho_p(0:n+1))
        allocate(self%vrel_p(0:n+1))
        allocate(self%Y_21(0:n+1))
        allocate(self%delta(0:n+1))

        self%r(:) = 0.0_DP
        self%X(:) = 0.0_DP
        self%X2(:) = 0.0_DP
        self%r_hstar(:) = 0.0_DP
        self%deltaA(:) = 0.0_DP
        self%sigma(:) = 0.0_DP
        self%Gsigma(:) = 0.0_DP
        self%mass(:) = 0.0_DP
        self%tau(:) = 0.0_DP
        self%nu(:) = 0.0_DP
        self%Q(:) = 0.0_DP
        self%Iz(:) = 0.0_DP
        self%nkep(:) = 0.0_DP
        self%Torque(:) = 0.0_DP
        self%r_p(:) = 0.0_DP
        self%m_p(:) = 0.0_DP
        self%rho_p(:) = 0.0_DP
        self%vrel_p(:) = 0.0_DP
        self%Y_21(:) = 0.0_DP
        self%delta(:) = 0.0_DP

        return
    end subroutine ringmoons_util_setup_ring

    module subroutine ringmoons_util_setup_seed(self, n, param)
        !! author: David A. Minton
        !!
        !! Constructor for the ringmoons_seed class. 
        !! Allocates space for all bins and initializes all components with a value. 
        implicit none
        class(ringmoons_seed),     intent(inout) :: self  
            !! Ringmoons seed object
        integer(I4B),               intent(in)    :: n
            !! Number of seeds to allocate space for
        class(swiftest_parameters), intent(in)    :: param 
            !! Current run configuration parameters
        integer(I4B) :: i

        if (n < 0) return
        call self%dealloc()

        self%nbody = n
        if (n == 0) return

        allocate(self%id(n))
        allocate(swiftest_particle_info :: self%info(n))
        allocate(self%status(n))
        allocate(self%a(n))
        allocate(self%mass(n))
        allocate(self%Gmass(n))
        allocate(self%mu(n))
        allocate(self%rhill(n))
        allocate(self%radius(n))
        allocate(self%density(n))
        allocate(self%ringbin(n))
        allocate(self%Torque(n))
        allocate(self%Ttide(n))

        do i = 1, n
            call self%info(i)%set_value(&
                name = "UNNAMED", &
                particle_type = "UNKNOWN", &
                status = "INACTIVE", & 
                origin_type = "UNKNOWN", &
                collision_id = 0, &
                origin_time = -huge(1.0_DP), & 
                origin_rh = [0.0_DP, 0.0_DP, 0.0_DP], &
                origin_vh = [0.0_DP, 0.0_DP, 0.0_DP], &
                discard_time = huge(1.0_DP), & 
                discard_rh = [0.0_DP, 0.0_DP, 0.0_DP], &
                discard_vh = [0.0_DP, 0.0_DP, 0.0_DP], &
                discard_body_id = -1  &
            )
        end do
        self%id(:) = 0
        self%status(:) = INACTIVE
        self%a(:) = 0.0_DP
        self%mass(:) = 0.0_DP
        self%Gmass(:) = 0.0_DP
        self%mu(:) = 0.0_DP
        self%rhill(:) = 0.0_DP
        self%radius(:) = 0.0_DP
        self%density(:) = 0.0_DP
        self%ringbin(:) = 0
        self%Torque(:) = 0.0_DP
        self%Ttide(:) = 0.0_DP

        return
    end subroutine ringmoons_util_setup_seed

    module subroutine ringmoons_util_setup_initialize_system(self, system_history, param)
        !! author: David A. Minton
        !!
        !! Initialize an SyMBA nbody system from files and sets up the encounter and collision structures
        !! 
        implicit none
        ! Arguments
        class(ringmoons_nbody_system), intent(inout) :: self 
            !! SyMBA nbody_system object
        class(swiftest_storage),allocatable, intent(inout) :: system_history 
            !! Stores the system history between output dumps
        class(swiftest_parameters), intent(inout) :: param 
            !! Current run configuration parameters 

        ! Call parent method
        select type(cb => self%cb)
        class is (ringmoons_cb)
        associate(nbody_system => self, ring=>self%ring, seed=>self%seed)
            call symba_util_setup_initialize_system(nbody_system, system_history, param)
            ring%nc%file_name = param%ring_file
            call ring%read_frame(self%t,param)
            call seed%read_frame(self%t,system_history%nc,param)
            call ring%reset(seed,cb,param)
        end associate
        end select

        return
    end subroutine ringmoons_util_setup_initialize_system

    module subroutine ringmoons_util_spawn_seed(self, cb, ring, ring_bin, param)
        implicit none
        ! Arguments
        class(ringmoons_seed),          intent(inout) :: self
        class(ringmoons_ring),          intent(inout) :: ring
        class(ringmoons_cb),            intent(in)    :: cb
        integer(I4B),                   intent(in)    :: ring_bin
        class(swiftest_parameters),     intent(in)    :: param
        ! Internals
        integer(I4B)        :: i,j,seed_bin,nfz
        type(ringmoons_seed)  :: new_seed
        character(*), parameter :: SEEDFMT = '("Seed",I0.7)'
        character(NAMELEN) :: newname
        real(DP) :: Lr0,Lr1,Ls,ker0,ker1,per0,per1,kes,pes,bes

        associate(seed => self, maxid => self%maxid)
            seed_bin = seed%nbody + 1 
            call new_seed%setup(seed_bin, param)
            if (seed%nbody > 0) then
                ! Copy over the old seed properties to the new
                new_seed%status(1:seed%nbody) = seed%status(1:seed%nbody)
                new_seed%id(1:seed%nbody) = seed%id(1:seed%nbody)
                new_seed%info(1:seed%nbody) = seed%info(1:seed%nbody)
                new_seed%a(1:seed%nbody) = seed%a(1:seed%nbody)
                new_seed%mu(1:seed%nbody) = seed%mu(1:seed%nbody)
                new_seed%mass(1:seed%nbody) = seed%mass(1:seed%nbody)
                new_seed%Gmass(1:seed%nbody) = seed%Gmass(1:seed%nbody)
                new_seed%density(1:seed%nbody) = seed%density(1:seed%nbody)
                new_seed%radius(1:seed%nbody) = seed%radius(1:seed%nbody)
                new_seed%rhill(1:seed%nbody) = seed%rhill(1:seed%nbody)
                new_seed%ringbin(1:seed%nbody) = seed%ringbin(1:seed%nbody)
                new_seed%Torque(1:seed%nbody) = seed%Torque(1:seed%nbody)
                new_seed%Ttide(1:seed%nbody) = seed%Ttide(1:seed%nbody)
                new_seed%ringbin(1:seed%nbody) = seed%ringbin(1:seed%nbody)
            end if
            call seed%setup(seed_bin, param)
            seed%id = new_seed%id
            seed%status = new_seed%status
            seed%info = new_seed%info
            seed%a = new_seed%a
            seed%mass = new_seed%mass
            seed%Gmass = new_seed%Gmass
            seed%mu = new_seed%mu
            seed%density = new_seed%density
            seed%radius = new_seed%radius
            seed%rhill = new_seed%rhill
            seed%ringbin = new_seed%ringbin
            seed%Torque = new_seed%Torque
            seed%Ttide = new_seed%Ttide
            seed%ringbin = new_seed%ringbin

            i = seed_bin 
            j = ring_bin
            seed%ringbin(i) = ring_bin
            seed%id(i) = -1 ! Assign a temporary id to the seed and only update it it survives long enough to be recorded
            seed%status(i) = ACTIVE
            seed%mass(i) = min(ring%mass(j),ring%m_p(j))
            seed%Gmass(i) = param%GU * seed%mass(i)
            seed%mu(i) = cb%Gmass + seed%Gmass(i)
            seed%a(i) = (ring%Iz(j) * ring%nkep(j))**2 / seed%mu(i)
            seed%density(i) = ring%rho_p(j)
            seed%radius(i) = (3 * seed%mass(i) / (4 * PI * seed%density(i)))**(1._DP / 3._DP)
            seed%rhill(i) = seed%a(i) * (seed%mass(i) / cb%mass / 3)**THIRD 
            write(newname, SEEDFMT) seed%id(i)
            call seed%info(i)%set_value(name=newname, particle_type=SEED_TYPE_NAME, status="ACTIVE")

            ! Testing angular momentum conservation and energy loss
            ! Lr0 = ring%mass(j) * ring%Iz(j) * ring%nkep(j)
            ! ker0 = 0.5_DP * ring%mass(j) * ring%Iz(j) * ring%nkep(j)**2
            ! per0 = -cb%Gmass*ring%mass(j) / ring%r(j)

            ! Take away the mass from the ring
            ring%mass(j) = ring%mass(j) - seed%mass(i)
            ring%sigma(j) = ring%mass(j) / ring%deltaA(j)
            ring%Gsigma(j) = param%GU * ring%sigma(j)
            seed%ringbin(i) = ring%find_bin(seed%a(i))

            ! Testing angular momentum conservation and energy loss
            ! ker1 = 0.5_DP * ring%mass(j) * ring%Iz(j) * ring%nkep(j)**2
            ! per1 = -cb%Gmass*ring%mass(j) / ring%r(j)
            ! kes = 0.5*seed%mass(i)*seed%mu(i)/seed%a(i)
            ! pes = -cb%Gmass*seed%mass(i)/seed%a(i)
            ! bes = -3*seed%Gmass(i)*seed%mass(i)/seed%a(i)
            ! Lr1 = ring%mass(j) * ring%Iz(j) * ring%nkep(j)
            ! Ls = seed%mass(i) * sqrt(seed%mu(i) * seed%a(i))

        end associate

        return
    end subroutine ringmoons_util_spawn_seed

    module subroutine ringmoons_util_velocity_dispersion_ring(self,cb,param)
        implicit none
        class(ringmoons_ring), intent(inout) :: self
        class(ringmoons_cb),   intent(in)    :: cb
        class(swiftest_parameters), intent(in) :: param
        ! Internals
        integer(I4B)                         :: i
        real(DP),dimension(0:self%nbins+1)   :: kappa_rhstar,eta_rhstar

        associate(ring => self)
            where(ring%mass(:) > N_DISK_FACTOR * ring%m_p(:))
                ring%r_hstar(:) = ring%r(:) * (2 * ring%m_p(:) /(3._DP * cb%mass))**(1._DP/3._DP) / (2 * ring%r_p(:)) 
                ! See Salmon et al. 2010 for this
                kappa_rhstar(:) = ringmoons_transition_function(ring%r_hstar(:))
                eta_rhstar(:) = 1._DP - kappa_rhstar(:)
                ring%vrel_p(:) = kappa_rhstar(:) * sqrt(param%GU * ring%m_p(:) / ring%r_p(:)) + eta_rhstar(:) * &
                                                    (2 * ring%r_p(:) * ring%nkep(:))
            end where
        end associate

        return
    end subroutine ringmoons_util_velocity_dispersion_ring

    elemental pure function ringmoons_viscosity(Gsigma, m_p, v2_p, r_p, r_hstar, Q, tau, w) result(nu)
        ! Arguments
        real(DP),intent(in) :: Gsigma, m_p, v2_p, r_p, r_hstar, Q, tau, w
        real(DP) :: nu
        ! Internals
        real(DP)       :: kappa_Q,eta_Q,y
        real(DP)       :: nu_trans_stable,nu_grav_stable,nu_trans_unstable,nu_grav_unstable
        real(DP)       :: nu_trans,nu_grav,nu_coll

        nu_trans_unstable = 13 * r_hstar**5 * Gsigma**2 / w**3
        nu_trans_stable = (v2_p / (2 * w)) * (0.46_DP * tau / (1._DP + tau**2))

        nu_grav_stable = 0.0_DP
        nu_grav_unstable = nu_trans_unstable

        y = 0.25_DP * Q 
        kappa_Q = ringmoons_transition_function(y)
        eta_Q = 1._DP - kappa_Q

        nu_trans = kappa_Q * nu_trans_stable + eta_Q * nu_trans_unstable
        nu_grav  = kappa_Q * nu_grav_stable  + eta_Q * nu_grav_unstable
        nu_coll = r_p**2 * w * tau

        nu = nu_trans + nu_grav + nu_coll
        return
    end function ringmoons_viscosity


    elemental pure module function ringmoons_transition_function(yin) result(kappa)
        implicit none
        ! Arguments
        real(DP),intent(in) ::yin
        real(DP) :: kappa
        ! Internals
        real(DP) :: y
        y = min(max(yin,epsilon(1._DP)),1.0_DP-epsilon(1._DP))
        kappa =  0.5_DP * (1._DP + tanh((2 * y - 1._DP) / (y * (1._DP - y)))) 
        return

    end function ringmoons_transition_function

end submodule s_ringmoons_util