Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Changes to the positining that mainly affects hit and runs to increas…
Browse files Browse the repository at this point in the history
…e the success rate of fragmentations after hit and run.
  • Loading branch information
daminton committed Feb 21, 2023
1 parent 5cd7c93 commit abc674a
Showing 1 changed file with 28 additions and 25 deletions.
53 changes: 28 additions & 25 deletions src/fraggle/fraggle_generate.f90
Original file line number Diff line number Diff line change
Expand Up @@ -326,20 +326,17 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
rdistance = 0.5_DP * sum(impactors%radius(:))
istart = 3
else
rdistance = 2 * impactors%radius(2)
rdistance = impactors%radius(2)
istart = 3
end if

if (lhitandrun) then
mass_rscale(:) = 1.0_DP
else
! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point
! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be.
call random_number(mass_rscale)
mass_rscale(:) = (mass_rscale(:) + 1.0_DP) / 2
mass_rscale(:) = mass_rscale(:) * (fragments%mtot / fragments%mass(1:nfrag))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence
mass_rscale(:) = mass_rscale(:) / maxval(mass_rscale(:))
end if
mass_rscale(1:istart-1) = 1.0_DP
! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point
! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be.
call random_number(mass_rscale(istart:nfrag))
mass_rscale(istart:nfrag) = (mass_rscale(istart:nfrag) + 1.0_DP) / 2
mass_rscale(istart:nfrag) = mass_rscale(istart:nfrag) * (sum(fragments%mass(istart:nfrag)) / fragments%mass(istart:nfrag))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence
mass_rscale(istart:nfrag) = mass_rscale(istart:nfrag) / minval(mass_rscale(istart:nfrag))

loverlap(:) = .true.
do loop = 1, MAXLOOP
Expand All @@ -353,41 +350,47 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
else ! Keep the first and second bodies at approximately their original location, but so as not to be overlapping
fragment_cloud_center(:,1) = impactors%rcimp(:) - rbuffer * max(fragments%radius(1),impactors%radius(1)) * impactors%y_unit(:)
fragment_cloud_center(:,2) = impactors%rcimp(:) + rbuffer * max(fragments%radius(2),impactors%radius(2)) * impactors%y_unit(:)
fragment_cloud_radius(:) = cloud_size_scale_factor * rdistance
fragments%rc(:,1) = fragment_cloud_center(:,1)
fragments%rc(:,2) = fragment_cloud_center(:,2)
fragment_cloud_radius(:) = rdistance / pack_density
fragments%rc(:,1:2) = fragment_cloud_center(:,1:2)
end if

do i = 1, nfrag
if (loverlap(i)) then
call random_number(phi(i))
call random_number(theta(i))
call random_number(u(i))
phi(i) = TWOPI * phi(i)
theta(i) = asin(2 * theta(i) - 1.0_DP)
end if
end do

! Randomly place the n>2 fragments inside their cloud until none are overlapping
do concurrent(i = istart:nfrag, loverlap(i))
j = fragments%origin_body(i)

! Make a random cloud
phi(i) = TWOPI * phi(i)
theta(i) = acos(2 * theta(i) - 1.0_DP)

! Scale the cloud size
fragments%rmag(i) = fragment_cloud_radius(j) * mass_rscale(i) * u(i)**(THIRD)
if (lhitandrun) then
fragments%rmag(i) = fragment_cloud_radius(j) * u(i)**(THIRD)
else
fragments%rmag(i) = fragment_cloud_radius(j) * mass_rscale(i) * u(i)**(THIRD)
end if

! Position the fragment in a random point within the cloud
fragments%rc(1,i) = fragments%rmag(i) * sin(theta(i)) * cos(phi(i))
fragments%rc(2,i) = fragments%rmag(i) * sin(theta(i)) * sin(phi(i))
fragments%rc(3,i) = fragments%rmag(i) * cos(theta(i))

! Shift to the cloud center coordinates
fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j)

! Stretch out the hit and run cloud along the flight trajectory
if (lhitandrun) then
fragments%rc(:,i) = fragments%rc(:,i) + cloud_size_scale_factor * rdistance * u(i) * impactors%bounce_unit(:)
fragments%rc(:,i) = fragments%rc(:,i) * (1.0_DP + fragment_cloud_radius(j) * mass_rscale(i) * impactors%bounce_unit(:))
end if

fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j)

if (lhitandrun) then
fragments%rc(:,i) = fragments%rc(:,i) + (fragment_cloud_radius(j) + maxval(mass_rscale)) * impactors%bounce_unit(:) ! Shift the stretched cloud downrange
else
! Make sure that the fragments are positioned away from the impact point
direction = dot_product(fragments%rc(:,i) - impactors%rcimp(:), fragment_cloud_center(:,j) - impactors%rcimp(:))
Expand Down Expand Up @@ -549,7 +552,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
real(DP) :: vmax_guess
integer(I4B), parameter :: MAXLOOP = 25
integer(I4B), parameter :: MAXTRY = 10
integer(I4B), parameter :: MAXANGMTM = 1000
integer(I4B), parameter :: MAXANGMTM = 30000
class(collision_fraggle), allocatable :: collider_local
character(len=STRMAX) :: message

Expand Down Expand Up @@ -682,11 +685,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
if (all(dL_metric(:) <= 1.0_DP)) exit angmtm

do i = istart, fragments%nbody
dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody))
dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody))
call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i))
call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)
fragments%vmag(i) = .mag.fragments%vc(:,i)
end do
call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)
fragments%vmag(:) = .mag.fragments%vc(:,:)
end do angmtm

call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)
Expand Down

0 comments on commit abc674a

Please sign in to comment.