diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 88c973aee..45b5acfbc 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -283,7 +283,7 @@ end subroutine symba_io_dump_particle_info module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none - class(symba_parameters), intent(inout) :: self !! Collection of SyMBA parameters + class(symba_parameters), intent(inout) :: self !! Current run configuration parameters with SyMBA additionss integer, intent(in) :: unit !! File unit number character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. @@ -294,7 +294,7 @@ end subroutine symba_io_param_reader module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, iomsg) implicit none - class(symba_parameters),intent(in) :: self !! Collection of SyMBA parameters + class(symba_parameters),intent(in) :: self !! Current run configuration parameters with SyMBA additions integer, intent(in) :: unit !! File unit number character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 13c874b0f..edb02dee5 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -143,7 +143,7 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt end function symba_collision_check_one - function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag) + function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag) !! author: David A. Minton !! !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all family members, @@ -151,6 +151,7 @@ function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass, implicit none ! Arguments class(symba_pl), intent(inout) :: pl !! SyMBA massive body object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision integer(I4B), dimension(:), allocatable, intent(out) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(NDIM,2), intent(out) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision @@ -209,6 +210,8 @@ function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass, family = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] fam_size = count(pl%status(family(:)) == ACTIVE) family = pack(family(:), pl%status(family(:)) == ACTIVE) + L_spin(:,:) = 0.0_DP + Ip(:,:) = 0.0_DP ! Find the barycenter of each body along with its children, if it has any do j = 1, 2 @@ -216,10 +219,8 @@ function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass, v(:, j) = pl%vb(:, idx_parent(j)) Ip(:, j) = mass(j) * pl%Ip(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) - L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j)) + if (param%lrotation) L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j)) - ! Use conservation of energy and angular momentum to adjust velocities and distances to be those equivalent to where they would be when the mutual radii are touching - !call vector_adjust(mass(j), x(:,j), v(:,j)) if (nchild(j) > 0) then do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties idx_child = parent_child_index_array(j)%idx(i + 1) @@ -230,31 +231,33 @@ function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass, volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 volume(j) = volume(j) + volchild ! Get angular momentum of the child-parent pair and add that to the spin - xcom(:) = (mass(j) * x(:,j) + mchild * xchild(:)) / (mass(j) + mchild) - vcom(:) = (mass(j) * v(:,j) + mchild * vchild(:)) / (mass(j) + mchild) - xc(:) = x(:, j) - xcom(:) - vc(:) = v(:, j) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - L_spin(:, j) = L_spin(:, j) + mass(j) * xcrossv(:) - - xc(:) = xchild(:) - xcom(:) - vc(:) = vchild(:) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - L_spin(:, j) = L_spin(:, j) + mchild * xcrossv(:) - ! Add the child's spin - L_spin(:, j) = L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child) + if (param%lrotation) then + xcom(:) = (mass(j) * x(:,j) + mchild * xchild(:)) / (mass(j) + mchild) + vcom(:) = (mass(j) * v(:,j) + mchild * vchild(:)) / (mass(j) + mchild) + xc(:) = x(:, j) - xcom(:) + vc(:) = v(:, j) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + L_spin(:, j) = L_spin(:, j) + mass(j) * xcrossv(:) + + xc(:) = xchild(:) - xcom(:) + vc(:) = vchild(:) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + L_spin(:, j) = L_spin(:, j) + mchild * xcrossv(:) + + L_spin(:, j) = L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child) + Ip(:, j) = Ip(:, j) + mchild * pl%Ip(:, idx_child) + end if ! Merge the child and parent mass(j) = mass(j) + mchild x(:, j) = xcom(:) v(:, j) = vcom(:) - Ip(:, j) = Ip(:, j) + mchild * pl%Ip(:, idx_child) end do end if density(j) = mass(j) / volume(j) radius(j) = ((3 * mass(j)) / (density(j) * 4 * pi))**(1.0_DP / 3.0_DP) - Ip(:, j) = Ip(:, j) / mass(j) + if (param%lrotation) Ip(:, j) = Ip(:, j) / mass(j) end do lflag = .true. @@ -418,6 +421,26 @@ module subroutine symba_collision_resolve_mergers(self, system, param) class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions ! Internals + integer(I4B) :: i + logical :: lgoodcollision + integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + + associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) + select type(pl => system%pl) + class is (symba_pl) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) + if (.not. lgoodcollision) cycle + if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved + !call symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + end do + end select + end associate return end subroutine symba_collision_resolve_mergers