N-Body merged velocity update

Help
1

Hi, I made a code to simulate N-Body collision and merging, the problem is after merging the body should continue moving but it stops instead

for example
if two bodies colliding where
Both Bodies : x velocity 5 and y velocity 0
and they are placed in a vertical line with some distance between them

they should collide & merge and one of them should sum all the data (mass, radius, position and velocity) while the other lose their mass and keep tracing the first one

(Im using gnuplot to show results) i will also link the portal to load it

program NBody
implicit none 

double precision :: pi,tmax, t, dt,rij,G
double precision , allocatable :: acc0(:,:)
double precision , allocatable :: posi(:,:), v(:,:),acc(:,:),M(:),r(:)
integer , allocatable :: link(:)
integer :: i, j, N ,lines, counter


N=2
allocate (posi(1:N,1:2), v(1:N,1:2),acc(1:N,1:2),M(1:N),r(1:N),link(1:N),acc0(1:N,1:2))



pi=4.0d0*atan(1.0d0) !pi
G=10d0!gravity
M(1:N)=1d0 !mass
r(1:N)=1d0 !radius
link(1:N)=0 !collision array

tmax=8d0 !max time
t=0d0 !time
dt=0.001d0 !timestep
lines=200d0 !(tmax/dt) !lines for printing


posi(1:N,1)=0.0d0 !(/-5d0,5d0/) !x position
posi(1:N,2)=(/-5d0,5d0/) !y position

v(1:N,1)=10d0 !velocity x
v(1:N,2)=0d0 !velocity y


acc(1:N,1)=0d0 !acceleration x
acc(1:N,2)=0d0 !acceleration y



open (unit=1 , file="results.txt")



do while (t<tmax) !time loop           
     counter=counter+1

     do i=1,N !position update loop
         if (link(i).NE.0) then !object position tracing
             posi(i,1:2)=posi(link(i),1:2)
         else !position update
             posi(i,1:2)=posi(i,1:2)+v(i,1:2)*dt+(0.5d0*acc(i,1:2)*dt**2d0)
         endif
     enddo


     do  i=1,N !position and acceleration update with collision case
	     !acc(i,1:2)=0d0 !acceleration
	     do j=1,N 
	   	     if (i==j) cycle 
			     rij=sqrt(  (posi(j,1)-posi(i,1))**2d0+ (posi(j,2)-posi(i,2))**2d0  )
				 !rij=sqrt(sum((posi(j,1:2)-posi(i,1:2))**2))
				 
				 acc0(i,1:2)=acc(i,1:2)
				 acc(i,1:2)=acc(i,1:2)+ G*M(j)*(posi(j,1:2)-posi(i,1:2))/rij**3d0
				 
			 if(rij<r(i)+r(j)) then !detect collision
			     link(j)=i !tracing
			     r(i)=r(i)+r(j)
			     v(i,1:2)=(M(i)*v(i,1:2)+M(j)*v(j,1:2))/(M(i)+M(j))
				 M(i)=M(i)+M(j)
				 M(j)=0d0
				 r(j)=0d0
				 posi(i,1:2)=(posi(i,1:2)+posi(j,1:2))/2.0d0
				 posi(j,1:2)=posi(i,1:2)
			 endif
			 
			
	     enddo
     enddo

     do i=1,N !velocity update loop
         if (link(i).NE.0) then !object velocity tracing
	         v(i,1:2)=v(link(i),1:2)
         else !velocity update
	         v(i,1:2)=v(i,1:2)+0.5d0*(acc0(i,1:2)+acc(i,1:2))*dt
         endif
     enddo


    ! do i=1,N !border collision loop
	!    if(posi(i,1)>=50d0)
	!	     v(i,1)=-1d0*v(i,1)
	!	 if (posi(i,2)>=50d0)
	!	     v(i,2)=-1d0*v(i,2)
	!	 endif
	! enddo
	 
	 

     if(mod(counter,(lines/100))==1) write (1,*) t,posi(1:N,1:2) 


     t=t+dt !time update
end do 

close (1)
deallocate (posi,v,acc,M,link,r)
   
end program NBody

and here’s the gnuplot portal

reset
set term qt
set xrange [-10:10]
set yrange [-10:10]
set size square
do for [i=1:400:1]{
plot  "results.txt" every ::i::i u 2:4 w p pt 7 ps 0.8 ,"results.txt" every ::1::i u 2:4 w l ,  "results.txt" every ::1::i u 3:5 w l, "results.txt" every ::i::i u 3:5 w p pt 7 ps 0.8
}
2

There are two things wrong in the following assignment.

The first error is that this violates Newton’s Second Law.

The second error is that rij can be zero after a collision has taken place, and that causes division by zero to occur.

A comment about the whole post: readers are placed in an uncomfortable position when you have a physical model that has not been described, and they have to infer features of the model from the code, and that code itself has errors.

When two bodies touch, you arbitrarily set the radius and mass of the merged object as the sums of the radii and masses prior to collision. This raises the question as to the relation between mass and radius. In 2D, mass ~ square(radius) and, in 3D, mass ~ cube(radius), but neither of these relations is being satisfied. So, what is the model being used?