! 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_step use swiftest contains module subroutine ringmoons_step_restructure_seed(self,cb,ring,param) use, intrinsic :: ieee_exceptions implicit none class(ringmoons_seed), intent(inout) :: self class(ringmoons_cb), intent(inout) :: cb class(ringmoons_ring), intent(inout) :: ring class(swiftest_parameters), intent(in) :: param integer(I4B) :: i,j,inner_outer_sign,nbin,Nactive,rbin real(DP) :: a, delta_mass, mass_left logical :: open_space,destructo,spawnbin real(DP), parameter :: NOSPAWN_THRESHOLD_SIGMA = 0.1_DP !! Threshold surface mass density below which seeds cannot spawn in units of kg/m^2 real(DP), parameter :: dzone_width = 0.025_DP !! Width of the destruction zone as a fraction of the RRL distance integer(I4B) :: dzone_inner,dzone_outer ! inner and outer destruction zone bins real(DP) :: ndz,Lseed_orig,c,b,dr,Lring,deltaL real(DP) :: mass_min,R_min,dt,dtleft logical :: stepfail class(ringmoons_ring), allocatable :: seedring logical, dimension(:), allocatable :: lactive 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.) ! First convert any recently destroyed satellites into ring material destructo = .false. associate(seed => self) if (seed%nbody > 0) then do i = 1, seed%nbody if (seed%status(i) == ACTIVE) then if (seed%a(i) <= ring%RRL) then ! Destroy the satellite! write(*,*) 'We are on our way to destruction!' ! DESTRUCTION_EVENT = .true. ! DESTRUCTION_COUNTER = 0 seed%status(i) = INACTIVE end if end if end do do i = 1,seed%nbody if ((seed%status(i) == INACTIVE).and.(seed%mass(i) > 0.0_DP)) then destructo = .true. Lseed_orig = seed%mass(i) * sqrt(seed%mu(i) * seed%a(i)) allocate(seedring, source=ring) seedring%mass(:) = 0.0_DP mass_left = seed%mass(i) c = dzone_width * seed%a(i) ! Create an approximately Gaussian distribution of mass rbin = ring%find_bin(seed%a(i)) a = mass_left / (sqrt(2 * PI) * c) do j = 0,(ring%nbins - rbin) do inner_outer_sign = -1,1,2 nbin = rbin + inner_outer_sign * j if ((nbin > seedring%inside).and.(nbin < seedring%nbins).and.(mass_left > 0.0_DP)) then dr = 0.5_DP * seedring%X(nbin) * seedring%deltaX delta_mass = min(mass_left,a * dr * exp(-(seedring%r(nbin) - seed%a(i))**2 / (2 * c**2))) seedring%mass(nbin) = seedring%mass(nbin) + delta_mass mass_left = mass_left - delta_mass end if if (j == 0) exit end do if (mass_left == 0.0_DP) exit end do ! Offset in angular momentum Lring = sum(seedring%mass(:) * seedring%Iz(:) * seedring%nkep(:)) deltaL = Lseed_orig - Lring ! Apply a torque to the temporary ring to bring it back to the seed's original angular momentum dt = 1._DP seedring%nu(1:ring%nbins) = 1._DP / (16 * 12 * dt / (seedring%deltaX)**2) / seedring%X2(1:ring%nbins) seedring%sigma(:) = seedring%mass(:) / seedring%deltaA(:) where (seedring%mass(:) > 0.0_DP) seedring%Torque(:) = deltaL * seedring%mass(:) / sum(seedring%mass(:)) elsewhere seedring%Torque(:) = 0.0_DP end where dtleft = dt dt = seedring%get_dt(dt) do call seedring%step(cb,dt,param,stepfail) dtleft = dtleft - dt if (dtleft <= 0.0_DP) exit dt = min(dtleft,dt) end do ring%mass(:) = ring%mass(:) + seedring%mass(:) ring%sigma(:) = ring%mass(:) / ring%deltaA(:) seed%mass(i) = 0.0_DP call seedring%dealloc() end if end do allocate(lactive(size(seed%status))) where(seed%status(:) == ACTIVE) lactive(:) = .true. elsewhere lactive(:) = .false. end where Nactive = count(lactive(:)) seed%a(1:Nactive) = pack(seed%a(:),lactive(:)) seed%mass(1:Nactive) = pack(seed%mass(:),lactive(:)) seed%status(1:Nactive) = ACTIVE if (size(seed%status) > Nactive) then seed%status(Nactive+1:size(seed%status)) = INACTIVE seed%ringbin(Nactive+1:size(seed%status)) = -1 seed%a(Nactive+1:size(seed%status)) = 0.0_DP seed%mass(Nactive+1:size(seed%status)) = 0.0_DP end if seed%nbody = Nactive seed%ringbin(1:seed%nbody) = ring%find_bin(seed%a(1:seed%nbody)) end if ! Check to see if we can spawn seeds beyond the FRL do j = ring%iFRL,ring%nbins ! skip bins that don't have enough mass if (ring%sigma(j) * param%MU2KG/param%DU2M**2 < NOSPAWN_THRESHOLD_SIGMA) cycle ! skip bins that already have a seed if (seed%nbody > 0) then if (any(seed%ringbin(1:seed%nbody) == j)) cycle end if call seed%spawn(cb,ring,j,param) end do if (seed%nbody > 0) then where (seed%status(:) == ACTIVE) seed%ringbin(:) = ring%find_bin(seed%a(:)) elsewhere seed%ringbin(:) = -1 end where end if end associate call ieee_set_halting_mode(IEEE_ALL, fpe_halting_modes) return end subroutine ringmoons_step_restructure_seed module subroutine ringmoons_step_ring(self,cb,dt,param,stepfail) use, intrinsic :: ieee_exceptions implicit none ! Arguments class(ringmoons_ring), intent(inout) :: self class(ringmoons_cb), intent(in) :: cb real(DP), intent(in) :: dt class(swiftest_parameters), intent(in) :: param logical, intent(out) :: stepfail ! Internals real(DP),dimension(0:self%nbins+1) :: S,Snew,Sn1,Sn2,fac,artnu,L,dM1,dM2 integer(I4B) :: i,N,j,loop logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes integer(I4B), parameter :: MAX_LOOP_SOLVER = 10000 real(DP), parameter :: ARTNU_FAC = 1.0_DP / 16.0_DP real(DP), save :: mass_ring_orig real(DP) :: mass_ring logical, save:: first_step = .true. !! Artificial viscosity factor to prevent negative mass bins. if (first_step) then mass_ring_orig = sum(self%mass(:)) first_step = .false. end if ! 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.) call self%update(cb,param) stepfail = .false. associate(ring => self, N => self%nbins) ! calculate yarkovsky-schach torque if flag is set if (param%lyarkovsky_schach) then call ring%yark_schach_torque(cb, param, ring%Torque(:)) end if S(0) = 0.0_DP S(1:N) = ring%sigma(1:N) * ring%X(1:N) S(N+1) = 0.0_DP ring%Torque(0) = 0.0_DP ring%Torque(N+1) = 0.0_DP ring%mass(0) = 0.0_DP ring%mass(N+1) = 0.0_DP ring%sigma(0) = 0.0_DP ring%sigma(N+1) = 0.0_DP fac(:) = 12 * dt / (ring%deltaX)**2 / ring%X2(:) Sn1(1:N) = ring%nu(2:N+1) * S(2:N+1) - 2 * ring%nu(1:N) * S(1:N) + ring%nu(0:N-1) * S(0:N-1) Sn2(1:N) = ring%Torque(2:N+1) - ring%Torque(0:N-1) Sn2(1:N) = Sn2(1:N) * (1._DP / (3 * PI * sqrt(cb%Gmass))) Snew(1:N) = S(1:N) + fac(1:N) * (Sn1(1:N) - Sn2(1:N)) ! Prevent any bins from having negative mass by diffusing mass with an artificial viscosity loop = 1 do while (any(Snew(1:N) < -epsilon(1._DP) * maxval(S(1:N)))) Snew(0) = 0.0_DP Snew(N+1) = 0.0_DP artnu(:) = 0.0_DP where (Snew(1:N) < 0._DP) artnu(1:N) = ARTNU_FAC / fac(1:N) artnu(0:N-1) = ARTNU_FAC / fac(0:N-1) artnu(2:N+1) = ARTNU_FAC / fac(2:N+1) end where S(1:N) = Snew(1:N) Sn1(1:N) = artnu(2:N+1) * S(2:N+1) - 2 * artnu(1:N) * S(1:N) + artnu(0:N-1) * S(0:N-1) Snew(1:N) = S(1:N) + fac(1:N) * Sn1(1:N) loop = loop + 1 if (loop > MAX_LOOP_SOLVER) then ! Give up and set all negative mass bins to 0 where(Snew(1:N) < epsilon(1.0_DP) * maxval(fac(1:N) * Sn1(1:N))) Snew(1:N) = 0.0_DP end where end if end do ring%sigma(1:N) = Snew(1:N) / ring%X(1:N) ring%Gsigma(1:N) = param%GU * ring%sigma(1:N) ring%mass(1:N) = ring%sigma(1:N) * ring%deltaA(1:N) ring%Torque(:) = 0.0_DP mass_ring = sum(ring%mass(:)) end associate call ieee_set_halting_mode(IEEE_ALL, fpe_halting_modes) return end subroutine ringmoons_step_ring module subroutine ringmoons_step_seed(self, cb, ring, dt, param, stepfail) !! author: David A. Minton !! !! Step the seed forward in time and handle any accretion events that occur during the step. use, intrinsic :: ieee_exceptions implicit none ! Arguments class(ringmoons_seed), intent(inout) :: self class(ringmoons_cb), intent(inout) :: cb class(ringmoons_ring), intent(inout) :: ring real(DP), intent(in) :: dt class(swiftest_parameters), intent(in) :: param logical, intent(out) :: stepfail ! Internals !Runge-Kutta-Fehlberg parameters integer(I4B),parameter :: rkfo = 6 real(DP),dimension(6,5),parameter :: rkf45_btab = reshape( & ! Butcher tableau for Runge-Kutta-Fehlberg method [ 1./4., 1./4., 0., 0., 0., 0.,& 3./8., 3./32., 9./32., 0., 0., 0.,& 12./13., 1932./2197., -7200./2197., 7296./2197., 0., 0.,& 1., 439./216., -8., 3680./513., -845./4104., 0.,& 1./2., -8./27., 2., -3544./2565., 1859./4104., -11./40.], shape(rkf45_btab)) real(DP),dimension(6),parameter :: rkf5_coeff = [ 16./135., 0., 6656./12825., 28561./56430., -9./50., 2./55. ] real(DP),dimension(6),parameter :: rkf4_coeff = [ 25./216., 0., 1408./2565. , 2197./4104. , -1./5. , 0. ] integer(I4B) :: i, j, iRRL, nfz, ilo,ihi, rkn,rbin, loop, nloops, Ns real(DP) :: dadt, sigavg, sigsum, Li, Lj, Ls,dti,dtleft,dtmin real(DP) :: impact_b class(ringmoons_ring), allocatable :: iring class(ringmoons_seed), allocatable :: iseed real(DP) :: da,mleft,dm,mdisk real(DP),dimension(0:ring%nbins+1) :: dTorque_ring,mringi,mringf,Torquei,Torquef real(DP),dimension(0:ring%nbins+1) :: Tlind,Tring real(DP),dimension(0:ring%nbins+1,rkfo) :: kr,kL real(DP),dimension(self%nbody,rkfo) :: ka,km,kT real(DP),dimension(0:ring%nbins+1) :: Er,rscale,rmdot real(DP),dimension(self%nbody) :: Ea, Em,ascale,mscale real(DP),dimension(self%nbody) :: ai,af,mi,mf, dTtide,Ttidef integer(I4B) :: Nactive real(DP),dimension(0:ring%nbins+1) :: Lring_orig,Lring_now real(DP),dimension(self%nbody) :: Lseeds_orig,Lseeds_now,Lres, mdot,adot logical, dimension(:), allocatable :: lactive real(DP) :: Lr0,Ls0,Lp0,Lr1,Ls1,Lp1,Lorig,sarr,maxE logical :: chomped,goodstep real(DP),parameter :: DTMIN_FAC = 1e-16_DP integer(I4B) :: Nnegative_seed,Nnegative_ring,Nbig_error 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(seed => self, Ns => self%nbody) ! The first time through we will initialize the laplace coefficients in the torque calculation stepfail = .false. if (seed%nbody == 0) return dti = dt dtleft = dt dtmin = DTMIN_FAC * dt ! Save initial state of the seed allocate(iring, source=ring) allocate(iseed, source=self) ai(1:Ns) = seed%a(1:Ns) mi(1:Ns) = seed%mass(1:Ns) mringi(:) = ring%mass(:) Torquei(:) = ring%Torque(:) dTorque_ring(:) = 0.0_DP dTtide(:) = 0.0_DP Nnegative_seed = 0 Nnegative_ring = 0 Nbig_error = 0 steploop: do loop = 1, LOOPMAX nloops = loop if (loop == LOOPMAX) then stepfail = .true. write(*,*) 'max loop reached in seed_evolve' return end if stepfail = .false. ka(:,:) = 0._DP km(:,:) = 0._DP kr(:,:) = 0._DP kT(:,:) = 0.0_DP kL(:,:) = 0.0_DP goodstep = .true. do rkn = 1,rkfo ! Runge-Kutta steps iseed%Torque(:) = 0.0_DP iseed%Ttide(:) = 0.0_DP ! Allow seeds to have negative mass/semimajor axis during intermediate steps, but we'll check at the end to make ! sure that the final result isn't negative iseed%a(1:Ns) = ai(1:Ns) + matmul(ka(1:Ns,1:rkn-1), rkf45_btab(2:rkn,rkn-1)) iseed%mass(1:Ns) = mi(1:Ns) + matmul(km(1:Ns,1:rkn-1), rkf45_btab(2:rkn,rkn-1)) iring%mass(:) = mringi(:) + matmul(kr(:,1:rkn-1),rkf45_btab(2:rkn,rkn-1)) iring%sigma(:) = iring%mass(:) / iring%deltaA(:) iring%Gsigma(:) = param%GU * iring%sigma(:) iseed%ringbin(1:Ns) = iring%find_bin(seed%a(1:Ns)) call iseed%get_tidal_torque(cb,param) mdot(1:Ns) = ringmoons_dMdt_seed(iseed,iring,cb,param) Tlind(:) = iring%get_lindblad_torque(cb,iseed,param) ! Add the correction factor that comes from conservation of angular momentum between the seed and the ring mass ! that it consumes where(mdot(1:Ns) > VSMALL) iseed%Torque(1:Ns) = iseed%Torque(1:Ns) + mdot(1:Ns) * iring%Iz(iseed%ringbin(1:Ns)) * iring%nkep(iseed%ringbin(1:Ns)) end where adot(:) = ringmoons_dadt_seed(iseed,cb,mdot) ! Set the RKF45 coefficients for this step kr(iseed%ringbin(1:Ns),rkn) = kr(iseed%ringbin(1:Ns),rkn) - dti * mdot(1:Ns) ! Eat the ring km(:,rkn) = dti * mdot(:) ! Grow the seed ka(:,rkn) = dti * adot(:) ! Migrate the seed kT(:,rkn) = dti * iseed%Ttide(:) ! Tidal torque on the seed kL(:,rkn) = kL(:,rkn) + dti * Tlind(:) ! Lindblad torque on the ring end do ! Allow ring mass to go negative, as it will get filled in by sigma_solver mringf(:) = mringi(:) + matmul(kr(:,1:rkfo), rkf5_coeff(1:rkfo)) af(1:Ns) = ai(1:Ns) + matmul(ka(1:Ns,1:rkfo), rkf5_coeff(1:rkfo)) !Don't let seed semimajor axes go negative if (any(af(1:Ns) < 0.0_DP)) then Nnegative_seed = Nnegative_seed + 1 dti = 0.5_DP * dti cycle steploop end if mf(1:Ns) = mi(1:Ns) + matmul(km(1:Ns,1:rkfo), rkf5_coeff(1:rkfo)) ! Don't let final masses go netative. if (any(mf(1:Ns) < 0.0_DP)) then Nnegative_seed = Nnegative_seed + 1 dti = 0.5_DP * dti cycle steploop end if ! use the initial value and derivative for error scaling ascale(1:Ns) = abs(ai(1:Ns)) + abs(ka(1:Ns,1)) Ea(1:Ns) = abs(matmul(ka(1:Ns,:), (rkf5_coeff(:) - rkf4_coeff(:)))) maxE = maxval(Ea(1:Ns) / ascale(1:Ns)) / iseed%rkf_tol if ((maxE > 1.0_DP).and.(dti > dtmin)) then ! seed a error too high dti = 0.9_DP * dti / maxE**(0.25_DP) goodstep =.false. Nbig_error = Nbig_error + 1 cycle steploop end if mscale(1:Ns) = abs(mi(1:Ns)) + abs(km(1:Ns,1)) Em(1:Ns) = abs(matmul(km(1:Ns,:), (rkf5_coeff(:) - rkf4_coeff(:)))) maxE = max(maxE, maxval(Em(1:Ns) / mscale(1:Ns)) / iseed%rkf_tol) if (maxE > 1.0_DP) then if (dti > dtmin) then ! seed m error too high dti = 0.9_DP * dti / maxE**(0.25_DP) goodstep =.false. Nbig_error = Nbig_error + 1 cycle steploop else ! already at minimum step size sarr = 1.0_DP end if else if (maxE < 2e-4_DP) then ! error very low sarr = 5._DP else ! adjust step size based on error estimate sarr = max(0.90_DP / maxE**(0.25_DP),1._DP) end if ! save final state of seed and ring and average of torques Torquef(:) = matmul(kL(:,1:rkfo), rkf5_coeff(1:rkfo)) Ttidef(1:Ns) = matmul(kT(1:Ns,1:rkfo), rkf5_coeff(1:rkfo)) ai(1:Ns) = af(1:Ns) mi(1:Ns) = mf(1:Ns) mringi(:) = mringf(:) dTtide(1:Ns) = dTtide(1:Ns) + Ttidef(1:Ns) dTorque_ring(:) = dTorque_ring(:) + Torquef(:) dtleft = dtleft - dti if (dtleft <= 0.0_DP) exit steploop dti = min(sarr * dti,dtleft) end do steploop seed%a(1:Ns) = af(1:Ns) seed%mass(1:Ns) = mf(1:Ns) seed%Gmass(1:Ns) = param%GU * seed%mass(1:Ns) seed%mu(1:Ns) = cb%Gmass + seed%Gmass(1:Ns) ring%mass(:) = mringf(:) ring%sigma(:) = ring%mass(:) / ring%deltaA(:) call ring%reset(seed,cb,param) call ring%update(cb,param) cb%dL(3) = cb%dL(3) - sum(dTtide(1:Ns)) ring%Torque(:) = Torquei(:) + dTorque_ring(:) / dt cb%rot(3) = (cb%L0(3) + cb%dL(3)) / (cb%Ip(3) * cb%mass * cb%radius**2) * RAD2DEG seed%Torque(:) = 0.0_DP seed%Ttide(:) = 0.0_DP seed%ringbin(1:Ns) = ring%find_bin(seed%a(1:Ns)) !I'm hungry! What's there to eat?! Look for neighboring seed chomped = .false. do i = 1, Ns if (seed%status(i) == ACTIVE) then do j = i + 1, Ns if (seed%status(j) == ACTIVE) then impact_b = 0.5_DP*(seed%a(i)+seed%a(j)) * ((seed%mass(i)+seed%mass(j)) / (3 * cb%mass))**THIRD impact_b = impact_b * seed%feeding_zone_factor if (abs(seed%a(i) - seed%a(j)) < impact_b) then ! conserve both mass and angular momentum Li = seed%mass(i) * sqrt(seed%mu(i) * seed%a(i)) Lj = seed%mass(j) * sqrt(seed%mu(j) * seed%a(j)) seed%mass(i) = seed%mass(i) + seed%mass(j) seed%Gmass(i) = param%GU * seed%mass(i) seed%mu(i) = cb%Gmass + seed%Gmass(i) seed%a(i) = ((Li + Lj) / seed%mass(i))**2 / seed%mu(i) seed%density(i) = 0.5_DP * (seed%density(i) + seed%density(j)) seed%radius(i) = (3*seed%mass(i) / (4 * PI * seed%density(i)))**THIRD seed%rhill(i) = seed%a(i) * (seed%mass(i) / (3 * cb%mass))**THIRD ! deactivate particle seed%status(j) = INACTIVE seed%a(j) = 0.0_DP seed%mass(j) = 0.0_DP seed%mu(j) = 0.0_DP seed%Gmass(j) = 0.0_DP seed%radius(j) = 0.0_DP seed%density(j) = 0.0_DP seed%rhill(j) = 0.0_DP chomped = .true. end if end if end do end if end do if (chomped) then allocate(lactive(size(seed%status))) lactive(:) = seed%status(:) == ACTIVE Nactive = count(lactive(:)) seed%id(1:Nactive) = pack(seed%id(:),lactive(:)) seed%a(1:Nactive) = pack(seed%a(:),lactive(:)) seed%mass(1:Nactive) = pack(seed%mass(:),lactive(:)) seed%Gmass(1:Nactive) = pack(seed%Gmass(:),lactive(:)) seed%mu(1:Nactive) = pack(seed%mu(:),lactive(:)) seed%info(1:Nactive) = pack(seed%info(:), lactive(:)) seed%status(1:Nactive) = pack(seed%status(:),lactive(:)) seed%density(1:Nactive) = pack(seed%density(:),lactive(:)) seed%radius(1:Nactive) = pack(seed%radius(:),lactive(:)) seed%rhill(1:Nactive) = pack(seed%rhill(:),lactive(:)) seed%ringbin(1:Nactive) = ring%find_bin(seed%a(1:Nactive)) seed%nbody = Nactive if (size(seed%id) > Nactive) then seed%status(Nactive+1:size(seed%id)) = INACTIVE seed%id(Nactive+1:size(seed%id)) = 0 seed%a(Nactive+1:size(seed%id)) = 0.0_DP seed%mass(Nactive+1:size(seed%id)) = 0.0_DP seed%Gmass(Nactive+1:size(seed%id)) = 0.0_DP seed%mu(Nactive+1:size(seed%id)) = 0.0_DP seed%density(Nactive+1:size(seed%id)) = 0.0_DP seed%rhill(Nactive+1:size(seed%id)) = 0.0_DP seed%radius(Nactive+1:size(seed%id)) = 0.0_DP seed%ringbin(Nactive+1:size(seed%id)) = 0 end if end if stepfail = .false. end associate call ieee_set_halting_mode(IEEE_ALL, fpe_halting_modes) return end subroutine ringmoons_step_seed function ringmoons_dadt_seed(seed,cb,mdot) result(adot) use, intrinsic :: ieee_exceptions implicit none ! Arguments class(ringmoons_seed), intent(in) :: seed class(ringmoons_cb), intent(in) :: cb real(DP), dimension(:), intent(in) :: mdot real(DP), dimension(1:seed%nbody) :: adot 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(Ns => seed%nbody) adot(1:Ns) = seed%Torque(1:Ns) / (sqrt(seed%mu(1:Ns) * seed%a(1:Ns))) & - mdot(1:Ns) * (1._DP + seed%Gmass(1:Ns) / (2 *seed%mu(1:Ns))) adot = 2 * adot * seed%a(1:Ns) / seed%mass(1:Ns) end associate call ieee_set_halting_mode(IEEE_ALL, fpe_halting_modes) return end function ringmoons_dadt_seed function ringmoons_dMdt_seed(seed,ring,cb,param) result(mdot) use, intrinsic :: ieee_exceptions implicit none class(ringmoons_seed),intent(in) :: seed class(ringmoons_ring), intent(in) :: ring class(ringmoons_cb), intent(in) :: cb class(swiftest_parameters), intent(in) :: param real(DP), dimension(1:seed%nbody) :: mdot ! Internals integer(I4B) :: i real(DP), dimension(seed%nbody) :: C real(DP),parameter :: eff2 = 1e-7_DP ! This term gets the growth rate to match up closely to Andy's real(DP),parameter :: growth_exp = 4._DP / 3._DP real(DP), parameter :: SIGLIMIT = 1e-100_DP 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(Ns => seed%nbody) ! Executable code where((seed%a(1:Ns) > VSMALL).and.ring%mass(seed%ringbin(1:Ns))/seed%mass(1:Ns) > epsilon(1.0_DP)) C(1:Ns) = 12*PI**(2._DP/3._DP)*(3._DP/(4*seed%density(1:Ns)))**(1._DP/3._DP)*sqrt(param%GU)/sqrt(cb%mass) mdot(1:Ns) = C(1:Ns)*ring%sigma(seed%ringbin(1:Ns))/(eff2*sqrt(seed%a(1:Ns)))*abs(seed%mass(1:Ns))**(growth_exp) elsewhere mdot(1:Ns) = 0.0_DP end where end associate return call ieee_set_halting_mode(IEEE_ALL, fpe_halting_modes) end function ringmoons_dMdt_seed module subroutine ringmoons_step_system(self, param, t, dt) !! author: David A. Minton !! !! Step the full nbody system forward in time and call the SyMBA step !! !! Adapted from Andrew Hesselbrock's RING-MOONS Python scripts implicit none ! Arguments class(ringmoons_nbody_system), intent(inout) :: self !! Ringmoons nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Simulation time real(DP), intent(in) :: dt ! Internals real(DP), parameter :: DTMIN_FAC = 1e-16_DP integer(I4B), parameter :: SUBMAX = 2 logical :: stepfail real(DP) :: dtring,dtleft integer(I4B) :: i,loop,subcount,loopcount class(ringmoons_nbody_system), allocatable :: old_system allocate(old_system, mold=self) subcount = 0 dtleft = dt dtring = dtleft allocate(old_system%ring, source=self%ring) allocate(old_system%seed, source=self%seed) allocate(old_system%cb, source=self%cb) do loop = 1, LOOPMAX select type(cb => self%cb) class is (ringmoons_cb) associate(ring => self%ring, seed => self%seed) if (loop == LOOPMAX) then write(*,*) 'LOOPMAX reached in seed evolution. Ringmoons_step failed' call base_util_exit(FAILURE) end if ring%Torque(:) = 0.0_DP if (seed%nbody > 0) seed%Torque(:) = 0.0_DP call ring%update(cb,param) dtring = ring%get_dt(dtring) call cb%accrete(ring,seed,param,dtring) call seed%step(cb,ring,dtring,param,stepfail) if (stepfail) then dtring = dtring / SUBMAX subcount = 0 call ringmoons_step_restart(self,old_system) cycle end if call ring%step(cb,dtring,param,stepfail) if (stepfail) then if (dtring > dt * DTMIN_FAC) then dtring = dtring / submax subcount = 0 call ringmoons_step_restart(self,old_system) cycle end if end if call seed%restructure(cb,ring,param) ! Spawn new seed in any available bins outside the FRL subcount = subcount + 1 dtleft = dtleft - dtring ! Scale the change in the ring torques by the step size reduction in order to get the time-averaged Torque if (dtleft <= 0.0_DP) exit loopcount = loop if (subcount == 2 * submax) then dtring = min(dtleft,submax * dtring) subcount = 0 end if dtring = min(dtleft,dtring) deallocate(old_system%cb, old_system%ring, old_system%seed) allocate(old_system%ring, source=self%ring) allocate(old_system%seed, source=self%seed) allocate(old_system%cb, source=self%cb) end associate end select end do self%ring%t = t + dt ! Step the nbody system like normal call symba_step_system(self, param, t, dt) return end subroutine ringmoons_step_system subroutine ringmoons_step_restart(system, old_system) implicit none class(ringmoons_nbody_system), intent(inout) :: system class(ringmoons_nbody_system), intent(inout) :: old_system call move_alloc(old_system%ring, system%ring) call move_alloc(old_system%seed, system%seed) call move_alloc(old_system%cb, system%cb) allocate(old_system%ring, source=system%ring) allocate(old_system%seed, source=system%seed) allocate(old_system%cb, source=system%cb) end subroutine ringmoons_step_restart end submodule s_ringmoons_step