diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 1e98ea6a6..3fac596b7 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -213,21 +213,23 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L ! Re-normalize position and velocity vectors by the fragment number so that for our initial guess we weight each ! fragment position by the mass and assume equipartition of energy for the velocity - r_col_norm = 1.5_DP * 2 * sum(rad_frag(:)) / theta / (nfrag - 2) + r_col_norm = 1.5_DP * 4 * sum(rad_frag(:)) / theta / (nfrag - 2) ! We will treat the first two fragments of the list as special cases. They get placed on the "poles" of the ! fragment distribution (aligned with the collisional system z-axis) - x_frag(:, 1) = -r_col_norm * z_col_unit(:) - x_frag(:, 2) = r_col_norm * z_col_unit(:) + x_frag(:, 1) = - z_col_unit(:) + x_frag(:, 2) = z_col_unit(:) + call random_number(x_frag(:,3:nfrag)) + x_frag(:,:) = x_frag(:,:) * r_col_norm ! The rest of the fragments will go in a ring on the x-y plane - do i = 3, nfrag - frag_vec(:) = [b2a * cos(theta * (i - 1)), sin(theta * (i - 1))] - frag_vec(:) = matmul(orientation(:,:), frag_vec(:)) - - r_frag_norm = r_col_norm - x_frag(:, i) = r_frag_norm * (frag_vec(1) * x_col_unit(:) + frag_vec(2) * y_col_unit(:)) - end do + !do i = 3, nfrag + ! frag_vec(:) = [b2a * cos(theta * (i - 1)), sin(theta * (i - 1))] + ! frag_vec(:) = matmul(orientation(:,:), frag_vec(:)) +! +! r_frag_norm = r_col_norm +! x_frag(:, i) = r_frag_norm * (frag_vec(1) * x_col_unit(:) + frag_vec(2) * y_col_unit(:)) +! end do call symba_frag_pos_com_adjust(xcom, m_frag, x_frag) call symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag)