diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index e91a19b7a..1e992b5f5 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -535,11 +535,10 @@ module subroutine collision_generate_disruption_vel_vec(collider) integer(I4B) :: i,j, loop logical :: lhitandrun, lnoncat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual - real(DP), dimension(NDIM,2) :: vbounce, vcloud - real(DP), dimension(2) :: vimp, mcloud - real(DP) :: vmag, vdisp, Lmag, Lscale, rotmag + real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof integer(I4B), dimension(collider%fragments%nbody) :: vsign - real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale + real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail + integer(I4B), parameter :: MAXLOOP = 100 associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) @@ -553,35 +552,11 @@ module subroutine collision_generate_disruption_vel_vec(collider) vsign(:) = 1 end where - ! The dominant velocities of the two clouds will be scaled by the CoM velocities of the two bodies - vimp(1:2) = .mag.impactors%vc(:,1:2) - - ! Now shift the CoM of each fragment cloud to what the origin body would have been in a bounce - vbounce(:,1) = -.mag.impactors%vc(:,1) * impactors%bounce_unit(:) - vbounce(:,2) = .mag.impactors%vc(:,2) * impactors%bounce_unit(:) - - ! Compute the current CoM of the fragment clouds - vcloud(:,:) = 0.0_DP - do concurrent(j = 1:2) - mcloud(j) = sum(fragments%mass(:), fragments%origin_body(:) == j) - do concurrent(i = 1:nfrag, fragments%origin_body(i) == j) - vcloud(:,j) = vcloud(:,j) + fragments%mass(i) * fragments%vc(:,i) - end do - vcloud(:,j) = vcloud(:,j) / mcloud(j) - end do - - ! Subtract off the difference between the cloud CoM velocity and the expected CoM velocity from bouncing - vcloud(:,:) = vcloud(:,:) - vbounce(:,:) - do concurrent(i = 1:nfrag) - j = fragments%origin_body(i) - fragments%vc(:,i) = fragments%vc(:,i) - vcloud(:,j) - end do - - ! Compute the velocity dispersion based on the escape velocity + ! The minimum fragment velocity will be set by the escape velocity if (lhitandrun) then - vdisp = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) + vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) else - vdisp = 2 * sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) + vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) end if ! Scale the magnitude of the velocity by the distance from the impact point @@ -590,44 +565,54 @@ module subroutine collision_generate_disruption_vel_vec(collider) rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) end do - vscale(:) = vscale(:)/maxval(vscale(:)) + vscale(:) = vscale(:)/minval(vscale(:)) ! Give the fragment velocities a random value that is scaled with fragment mass call random_number(mass_vscale) mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2 mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The 1/8 power is arbitrary. It just gives the velocity a small mass dependence - mass_vscale(:) = 2*mass_vscale(:) / maxval(mass_vscale(:)) + mass_vscale(:) = mass_vscale(:) / minval(mass_vscale(:)) ! Set the velocities of all fragments using all of the scale factors determined above - fragments%vc(:,1) = .mag.impactors%vc(:,1) * impactors%bounce_unit(:) - do concurrent(i = 2:nfrag) + do concurrent(i = 1:nfrag) j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) - vmag = .mag.fragments%vc(:,i) * vscale(i) * mass_vscale(i) + vmag = vesc * vscale(i) * mass_vscale(i) if (lhitandrun) then - fragments%vc(:,i) = vmag * 0.5_DP * impactors%bounce_unit(:) * vsign(i) + vrot(:) + fragments%vc(:,i) = vmag * impactors%bounce_unit(:) * vsign(i) + vrot(:) else ! Add more velocity dispersion to disruptions vs hit and runs. rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) vimp_unit(:) = .unit. rimp(:) - fragments%vc(:,i) = vmag * 0.5_DP * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:) + fragments%vc(:,i) = vmag * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:) end if end do - ! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body - call collider%set_coordinate_system() - call fragments%get_angular_momentum() + do loop = 1, MAXLOOP + call collider%set_coordinate_system() + call fragments%get_kinetic_energy() + ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape + ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%vmag(:) + ke_per_dof = -ke_residual/nfrag + do concurrent(i = 1:nfrag, ke_avail(i) > ke_per_dof) + fragments%vmag(i) = sqrt(2 * (fragments%ke_orbit_frag(i) - ke_per_dof)/fragments%mass(i)) + fragments%vc(:,i) = fragments%vmag(i) * fragments%v_unit(:,i) + end do + call fragments%get_kinetic_energy() + ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + + ! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body + call collider%set_coordinate_system() + call fragments%get_angular_momentum() + Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) + rotmag = .mag. fragments%rot(:,1) + fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) + rotmag = .mag. fragments%rot(:,1) - Lscale = fragments%mtot * sum(fragments%radius(:)) * sum(vimp(:)) - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) - Lmag = .mag.Lresidual(:)/Lscale - rotmag = .mag. fragments%rot(:,1) - fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) - rotmag = .mag. fragments%rot(:,1) + if (ke_residual >= 0.0_DP) exit - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) - Lmag = .mag.Lresidual(:)/Lscale + end do do concurrent(i = 1:nfrag) fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) @@ -639,12 +624,9 @@ module subroutine collision_generate_disruption_vel_vec(collider) end do impactors%vbcom(:) = impactors%vbcom(:) / fragments%mtot - ! Distribute any remaining angular momentum into fragments pin - call fragments%set_spins() - end associate return end subroutine collision_generate_disruption_vel_vec -end submodule s_collision_generate \ No newline at end of file +end submodule s_collision_generate diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 217e9a74e..c9d81322f 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -106,16 +106,19 @@ module collision real(DP), dimension(nbody) :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame real(DP), dimension(nbody) :: vmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame real(DP), dimension(nbody) :: rotmag !! Array of rotation magnitudes of individual fragments - real(DP), dimension(NDIM,nbody) :: v_r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(NDIM,nbody) :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(NDIM,nbody) :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: v_unit !! Array of velocity direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from - real(DP), dimension(NDIM) :: Lorbit !! Orbital angular momentum vector of all fragments - real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments - real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments - real(DP) :: ke_spin !! Spin kinetic energy of all fragments - real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories - real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories + real(DP), dimension(NDIM) :: Lorbit !! Orbital angular momentum vector of all fragments + real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments + real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments + real(DP) :: ke_spin !! Spin kinetic energy of all fragments + real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories + real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories + real(DP), dimension(nbody) :: ke_orbit_frag !! Orbital kinetic energy of each individual fragment + real(DP), dimension(nbody) :: ke_spin_frag !! Spin kinetic energy of each individual fragment contains procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index d7ae591f1..0d76af301 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -205,21 +205,23 @@ module subroutine collision_util_get_kinetic_energy(self) !! Calculates the current kinetic energy of the fragments implicit none ! Argument - class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + class(collision_fragments(*)), intent(inout) :: self !! Fragment system object ! Internals integer(I4B) :: i associate(fragments => self, nfrag => self%nbody) - fragments%ke_orbit = 0.0_DP - fragments%ke_spin = 0.0_DP + fragments%ke_orbit_frag(:) = 0.0_DP + fragments%ke_spin_frag(:) = 0.0_DP - do i = 1, nfrag - fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) - fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) + do concurrent(i = 1:nfrag) + fragments%ke_orbit_frag(i) = fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) + fragments%ke_spin_frag(i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) end do - fragments%ke_orbit = fragments%ke_orbit / 2 - fragments%ke_spin = fragments%ke_spin / 2 + fragments%ke_orbit_frag(:) = fragments%ke_orbit_frag(:) / 2 + fragments%ke_spin_frag(:) = fragments%ke_spin_frag(:) / 2 + fragments%ke_orbit = sum(fragments%ke_orbit_frag(:)) + fragments%ke_spin = sum(fragments%ke_spin_frag(:)) end associate @@ -374,9 +376,9 @@ module subroutine collision_util_reset_fragments(self) self%density(:) = 0.0_DP self%rc(:,:) = 0.0_DP self%vc(:,:) = 0.0_DP - self%v_r_unit(:,:) = 0.0_DP - self%v_t_unit(:,:) = 0.0_DP - self%v_n_unit(:,:) = 0.0_DP + self%r_unit(:,:) = 0.0_DP + self%t_unit(:,:) = 0.0_DP + self%n_unit(:,:) = 0.0_DP return end subroutine collision_util_reset_fragments @@ -453,9 +455,10 @@ module subroutine collision_util_set_coordinate_collider(self) fragments%vmag(:) = .mag. fragments%vc(:,:) ! Define the radial, normal, and tangential unit vectors for each individual fragment - fragments%v_r_unit(:,:) = .unit. fragments%rc(:,:) - fragments%v_t_unit(:,:) = .unit. fragments%vc(:,:) - fragments%v_n_unit(:,:) = .unit. (fragments%v_r_unit(:,:) .cross. fragments%v_t_unit(:,:)) + fragments%r_unit(:,:) = .unit. fragments%rc(:,:) + fragments%v_unit(:,:) = .unit. fragments%vc(:,:) + fragments%n_unit(:,:) = .unit. (fragments%rc(:,:) .cross. fragments%vc(:,:)) + fragments%t_unit(:,:) = .unit. (fragments%r_unit(:,:) .cross. fragments%r_unit(:,:)) end associate @@ -770,9 +773,9 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag) self%fragments%vc(:,:) = 0.0_DP self%fragments%rot(:,:) = 0.0_DP self%fragments%Ip(:,:) = 0.0_DP - self%fragments%v_r_unit(:,:) = 0.0_DP - self%fragments%v_t_unit(:,:) = 0.0_DP - self%fragments%v_n_unit(:,:) = 0.0_DP + self%fragments%r_unit(:,:) = 0.0_DP + self%fragments%t_unit(:,:) = 0.0_DP + self%fragments%n_unit(:,:) = 0.0_DP self%fragments%mass(:) = 0.0_DP self%fragments%radius(:) = 0.0_DP self%fragments%density(:) = 0.0_DP diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index cff097b4b..f8361be19 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -181,7 +181,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) Efunc = lambda_obj(E_objective_function) tol = TOL_INIT - fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:) + fragments%r_unit(:,:) = .unit. fragments%vc(:,:) fragments%vmag(:) = .mag. fragments%vc(:,1:nfrag) input_v(:) = fragments%vmag(1:nfrag) lfirst_Efunc = .true. @@ -191,7 +191,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) fval = E_objective_function(output_v) fragments%vmag(1:nfrag) = output_v(1:nfrag) do concurrent(i=1:nfrag) - fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i) + fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%r_unit(:,i) end do ! Set spins in order to conserve angular momentum @@ -226,7 +226,7 @@ function E_objective_function(val_input) result(fval) allocate(tmp_frag, source=fragments) tmp_frag%vmag(1:nfrag) = val_input(1:nfrag) do concurrent(i=1:nfrag) - tmp_frag%vc(:,i) = abs(tmp_frag%vmag(i)) * tmp_frag%v_r_unit(:,i) + tmp_frag%vc(:,i) = abs(tmp_frag%vmag(i)) * tmp_frag%r_unit(:,i) end do call collision_util_shift_vector_to_origin(tmp_frag%mass, tmp_frag%vc) @@ -303,8 +303,8 @@ subroutine fraggle_generate_tan_vel(collider, lfailure) tol = TOL_INIT lfirst_func = .true. do i = 1, nfrag - fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%v_t_unit(:,i)) - fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%v_r_unit(:,i)) + fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%t_unit(:,i)) + fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%r_unit(:,i)) end do allocate(v_t_initial, source=fragments%v_t_mag) do while(tol < TOL_MIN) @@ -322,8 +322,8 @@ subroutine fraggle_generate_tan_vel(collider, lfailure) fragments%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) - fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & - fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) + fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & + fragments%t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) do concurrent (i = 1:nfrag) fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:) end do @@ -393,10 +393,10 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) L_orb_others(:) = 0.0_DP do i = 1, nfrag if (i <= 2 * NDIM) then ! The tangential velocities of the first set of bodies will be the unknowns we will solve for to satisfy the constraints - A(1:3, i) = fragments%mass(i) * fragments%v_t_unit(:, i) - A(4:6, i) = fragments%mass(i) * fragments%rmag(i) * (fragments%v_r_unit(:, i) .cross. fragments%v_t_unit(:, i)) + A(1:3, i) = fragments%mass(i) * fragments%t_unit(:, i) + A(4:6, i) = fragments%mass(i) * fragments%rmag(i) * (fragments%r_unit(:, i) .cross. fragments%t_unit(:, i)) else if (present(v_t_mag_input)) then - vtmp(:) = v_t_mag_input(i - 6) * fragments%v_t_unit(:, i) + vtmp(:) = v_t_mag_input(i - 6) * fragments%t_unit(:, i) L_lin_others(:) = L_lin_others(:) + fragments%mass(i) * vtmp(:) L(:) = fragments%mass(i) * (fragments%rc(:, i) .cross. vtmp(:)) L_orb_others(:) = L_orb_others(:) + L(:) @@ -436,8 +436,8 @@ function tangential_objective_function(v_t_mag_input, lfailure) result(fval) lfailure = .false. v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lfailure=lfailure) - vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), v_t_new(1:nfrag), & - fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) + vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%r_unit(:,1:nfrag), v_t_new(1:nfrag), & + fragments%t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) do concurrent (i = 1:nfrag) vc(:,i) = vb(:,i) - impactors%vbcom(:) end do @@ -491,8 +491,8 @@ subroutine fraggle_generate_rad_vel(collider, lfailure) ke_radial = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit do i = 1, nfrag - fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%v_t_unit(:,i)) - fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%v_r_unit(:,i)) + fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%t_unit(:,i)) + fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%r_unit(:,i)) end do allocate(v_r_initial, source=fragments%v_r_mag) @@ -515,8 +515,8 @@ subroutine fraggle_generate_rad_vel(collider, lfailure) ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) fragments%ke_orbit = 0.0_DP - fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), & - fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) + fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%r_unit(:,1:nfrag), & + fragments%v_t_mag(1:nfrag), fragments%t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) do i = 1, nfrag fragments%vc(:, i) = fragments%vb(:, i) - impactors%vbcom(:) fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:, i), fragments%vc(:, i)) @@ -561,7 +561,7 @@ function radial_objective_function(v_r_mag_input) result(fval) select type(fragments => collider%fragments) class is (fraggle_fragments(*)) allocate(v_shift, mold=fragments%vb) - v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%v_r_unit, fragments%v_t_mag, fragments%v_t_unit, fragments%mass, impactors%vbcom) + v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%r_unit, fragments%v_t_mag, fragments%t_unit, fragments%mass, impactors%vbcom) do i = 1,fragments%nbody v_shift(:,i) = v_shift(:,i) - impactors%vbcom(:) rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2 diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 8d6964aa4..2c5ee9672 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -91,11 +91,11 @@ module subroutine fraggle_util_set_original_scale_factors(self) class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object end subroutine fraggle_util_set_original_scale_factors - module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) + module function fraggle_util_vmag_to_vb(v_r_mag, r_unit, v_t_mag, t_unit, m_frag, vcom) result(vb) implicit none real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint - real(DP), dimension(:,:), intent(in) :: v_r_unit, v_t_unit !! Radial and tangential unit vectors for each fragment + real(DP), dimension(:,:), intent(in) :: r_unit, t_unit !! Radial and tangential unit vectors for each fragment real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass real(DP), dimension(:,:), allocatable :: vb diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 94e245c0b..cff4fd832 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -52,9 +52,9 @@ module subroutine fraggle_util_reset_fragments(self) self%rb(:,:) = 0.0_DP self%vb(:,:) = 0.0_DP self%rot(:,:) = 0.0_DP - self%v_r_unit(:,:) = 0.0_DP - self%v_t_unit(:,:) = 0.0_DP - self%v_n_unit(:,:) = 0.0_DP + self%r_unit(:,:) = 0.0_DP + self%t_unit(:,:) = 0.0_DP + self%n_unit(:,:) = 0.0_DP self%rmag(:) = 0.0_DP self%rotmag(:) = 0.0_DP @@ -226,7 +226,7 @@ module subroutine fraggle_util_setup_fragments_system(self, nfrag) end subroutine fraggle_util_setup_fragments_system - module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) + module function fraggle_util_vmag_to_vb(v_r_mag, r_unit, v_t_mag, t_unit, m_frag, vcom) result(vb) !! Author: David A. Minton !! !! Converts radial and tangential velocity magnitudes into barycentric velocity @@ -234,7 +234,7 @@ module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_ ! Arguments real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint - real(DP), dimension(:,:), intent(in) :: v_r_unit, v_t_unit !! Radial and tangential unit vectors for each fragment + real(DP), dimension(:,:), intent(in) :: r_unit, t_unit !! Radial and tangential unit vectors for each fragment real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass ! Result @@ -242,17 +242,17 @@ module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_ ! Internals integer(I4B) :: i, nfrag - allocate(vb, mold=v_r_unit) + allocate(vb, mold=r_unit) ! Make sure the velocity magnitude stays positive nfrag = size(m_frag) do i = 1, nfrag - vb(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) + vb(:,i) = abs(v_r_mag(i)) * r_unit(:, i) end do ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass call collision_util_shift_vector_to_origin(m_frag, vb) do i = 1, nfrag - vb(:, i) = vb(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) + vb(:, i) = vb(:, i) + v_t_mag(i) * t_unit(:, i) + vcom(:) end do return