diff --git a/src/tides/tides_getacch_pl.f90 b/src/tides/tides_getacch_pl.f90 deleted file mode 100644 index 4feb76221..000000000 --- a/src/tides/tides_getacch_pl.f90 +++ /dev/null @@ -1,65 +0,0 @@ -submodule(swiftest_classes) s_tides_kick_getacch - use swiftest -contains - - module subroutine tides_kick_getacch_pl(self, system) - !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton - !! - !! Calculated tidal torques from central body to any planet and from any planet to central body - !! planet - planet interactions are considered negligable. - !! This is a constant time lag model. - !! - !! Adapted from Mercury-T code from Bolmont et al. (2015) - !! - !! Reference: - !! Bolmont, E., Raymond, S.N., Leconte, J., Hersant, F., Correia, A.C.M., 2015. - !! Mercury-T : A new code to study tidally evolving multi-planet systems. - !! Applications to Kepler-62. A&A 583, A116. https://doi.org/10.1051/0004-6361/201525909 - implicit none - ! Arguments - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - ! Internals - integer(I4B) :: i - real(DP) :: rmag, vmag - real(DP), dimension(NDIM) :: r_unit, v_unit, h_unit, theta_unit, theta_dot, F_T - real(DP) :: Ftr, Ptopl, Ptocb, r5cbterm, r5plterm - - associate(pl => self, npl => self%nbody, cb => system%cb) - pl%atide(:,:) = 0.0_DP - cb%atide(:) = 0.0_DP - do i = 1, npl - rmag = norm2(pl%xh(:,i)) - vmag = norm2(pl%vh(:,i)) - r_unit(:) = pl%xh(:,i) / rmag - v_unit(:) = pl%vh(:,i) / vmag - h_unit(:) = r_unit(:) .cross. v_unit(:) - theta_unit(:) = h_unit(:) .cross. r_unit(:) - theta_dot = dot_product(pl%vh(:,i), theta_unit(:)) - - ! First calculate the tangential component of the force vector (eq. 5 & 6 of Bolmont et al. 2015) - ! The radial component is already computed in the obl_acc methods - r5cbterm = pl%Gmass(i)**2 * cb%k2 * cb%radius**5 - r5plterm = cb%Gmass**2 * pl%k2(i) * pl%radius(i)**5 - - Ptopl = 3 * r5plterm * pl%tlag(i) / rmag**7 - Ptocb = 3 * r5cbterm * cb%tlag / rmag**7 - - Ftr = -3 / rmag**7 * (r5cbterm + r5plterm) - 3 * vmag / rmag * (Ptocb + Ptopl) - - F_T(:) = (Ftr + (Ptocb + Ptopl) * dot_product(v_unit, r_unit) / rmag) * r_unit(:) & - + Ptopl * ((pl%rot(:,i) - theta_dot(:)) .cross. r_unit(:)) & - + Ptocb * ((cb%rot(:) - theta_dot(:)) .cross. r_unit(:)) - cb%atide(:) = cb%atide(:) + F_T(:) / cb%Gmass - pl%atide(:,i) = F_T(:) / pl%Gmass(i) - end do - - do i = 1, npl - pl%ah(:,i) = pl%ah(:,i) + pl%atide(:,i) + cb%atide(:) - end do - end associate - - return - end subroutine tides_kick_getacch_pl - -end submodule s_tides_kick_getacch \ No newline at end of file diff --git a/src/tides/tides_spin_step.f90 b/src/tides/tides_spin_step.f90 deleted file mode 100644 index 576aff8d7..000000000 --- a/src/tides/tides_spin_step.f90 +++ /dev/null @@ -1,132 +0,0 @@ -submodule(swiftest_classes) s_tides_step_spin - use swiftest - - type, extends(lambda_obj_tvar) :: tides_derivs_func - !! Base class for an lambda function object. This object takes no additional arguments other than the dependent variable x, an array of real numbers - procedure(tidederiv), pointer, nopass :: lambdaptr_tides_deriv - real(DP), dimension(:,:), allocatable :: xbeg - real(DP), dimension(:,:), allocatable :: xend - real(DP) :: dt - contains - generic :: init => tides_derivs_init - procedure :: evalt => tides_derivs_eval - procedure, nopass :: tides_derivs_init - end type - interface lambda_obj - module procedure tides_derivs_init - end interface - abstract interface - function tidederiv(x, t, dt, xbeg, xend) result(y) - ! Template for a 0 argument function - import DP, swiftest_nbody_system - real(DP), dimension(:), intent(in) :: x - real(DP), intent(in) :: t - real(DP), intent(in) :: dt - real(DP), dimension(:,:), intent(in) :: xbeg - real(DP), dimension(:,:), intent(in) :: xend - real(DP), dimension(:), allocatable :: y - end function - end interface - -contains - - module subroutine tides_step_spin_system(self, param, t, dt) - !! author: Jennifer L.L. Pouplin and David A. Minton - !! - !! Integrates the spin equations for central and massive bodies of the system subjected to tides. - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Simulation time - real(DP), intent(in) :: dt !! Current stepsize - ! Internals - real(DP), dimension(:), allocatable :: rot0, rot1 - real(DP) :: subt - real(DP), parameter :: tol=1e-6_DP !! Just a guess at the moment - real(DP) :: subdt - - associate(pl => self%pl, npl => self%pl%nbody, cb => self%cb) - allocate(rot0(NDIM*(npl+1))) - rot0 = [pack(pl%rot(:,1:npl),.true.), pack(cb%rot(:),.true.)] - ! Use this space call the ode_solver, passing tides_spin_derivs as the function: - subdt = dt / 20._DP - !rot1(:) = util_solve_rkf45(lambda_obj(tides_spin_derivs, subdt, pl%xbeg, pl%xend), rot0, dt, subdt tol) - ! Recover with unpack - !pl%rot(:,1:npl) = unpack(rot1... - !cb%rot(:) = unpack(rot1... - end associate - - return - end subroutine tides_step_spin_system - - - function tides_spin_derivs(rot_pl_cb, t, dt, xbeg, xend) result(drot) !! Need to add more arguments so we can pull in mass, radius, Ip, J2, etc... - !! author: Jennifer L.L. Pouplin and David A. Minton - !! - !! function used to calculate the derivatives that are fed to the ODE solver - implicit none - ! Arguments - real(DP), dimension(:,:), intent(in) :: rot_pl_cb !! Array of rotations. The last element is the central body, and all others are massive bodies - real(DP), intent(in) :: t !! Current time, which is used to interpolate the massive body positions - real(DP), intent(in) :: dt !! Total step size - real(DP), dimension(:,:), intent(in) :: xbeg - real(DP), dimension(:,:), intent(in) :: xend - ! Internals - real(DP), dimension(:,:), allocatable :: drot - real(DP), dimension(:), allocatable :: flatrot - real(DP), dimension(NDIM) :: N_Tcb, N_Rcb, N_Tpl, N_Rpl, xinterp - real(DP) :: C_cb, C_pl, r_dot_rot_cb, r_dot_rot_pl, rmag - integer(I4B) :: i, n - - - n = size(rot_pl_cb,2) - if (allocated(drot)) deallocate(drot) - allocate(drot, mold=rot_pl_cb) - drot(:,:) = 0.0_DP - do i = 1,n-1 - xinterp(:) = xbeg(:,i) + t / dt * (xend(:,i) - xbeg(:,i)) - ! Calculate Ncb and Npl as a function of xinterp - !drot(:,i) = -Mcb / (Mcb + Mpl(i)) * (N_Tpl + N_Rpl) - !drot(:,n) = drot(:,n) - Mcb / (Mcb + Mpl(i) * (N_Tcb + N_Rcb) - ! - end do - - return - end function tides_spin_derivs - - function tides_derivs_eval(self, x, t) result(y) - implicit none - ! Arguments - class(tides_derivs_func), intent(inout) :: self - real(DP), dimension(:), intent(in) :: x - real(DP), intent(in) :: t - ! Result - real(DP), dimension(:), allocatable :: y - if (associated(self%lambdaptr_tides_deriv)) then - y = self%lambdaptr_tides_deriv(x, t, self%dt, self%xbeg, self%xend) - else - error stop "Lambda function was not initialized" - end if - - return - end function tides_derivs_eval - - function tides_derivs_init(lambda, dt, xbeg, xend) result(f) - implicit none - ! Arguments - procedure(tidederiv) :: lambda - real(DP), intent(in) :: dt - real(DP), dimension(:,:), intent(in) :: xbeg - real(DP), dimension(:,:), intent(in) :: xend - ! Result - type(tides_derivs_func) :: f - f%lambdaptr_tides_deriv => lambda - f%dt = dt - allocate(f%xbeg, source = xbeg) - allocate(f%xend, source = xend) - - return - end function tides_derivs_init - -end submodule s_tides_step_spin \ No newline at end of file