subroutine particles(k0)
use varible
real bsigma
integer ::i,j,k,l,m,n,pnum
!-------------------求解顆粒處流體的速度-----------
do i=1,num
   fvx(i)=0.0
   fvy(i)=0.0
   do k=1,k0+1
      do L=1,knp
	     if(ndp(k,L)==0)then 
		   continue
		 else
           dis=sqrt((px(i)-Gx(k,L))**2+(py(i)-Gy(k,L))**2)
	       if(dis<sigma0(k,L))then
		     dis=sigma0(k,L)
		   else 
		     continue
		   end if
		   fvx(i)=fvx(i)+tao(k,L)*(-(py(i)-Gy(k,L)))/2/pi/dis/dis
		   fvy(i)=fvy(i)+tao(k,L)*(px(i)-Gx(k,L))/2/pi/dis/dis
		end if
	  end do
   end do
   do j=1,nall
      dis=sqrt((px(i)-bv(j,1))**2+(py(i)-bv(j,2))**2)
	  bsigma=max(abs(0.5*(bv(j+1,1)-bv(j,1))),abs(0.5*(bv(j+1,2)-bv(j,2))))
	  if(dis<bsigma)then
	    dis=bsigma
	  else
	    continue
	  end if
	  fvx(i)=fvx(i)+btao(j)*(-(py(i)-bv(j,2)))/2/pi/dis/dis
	  fvy(i)=fvy(i)+btao(j)*(px(i)-bv(j,1))/2/pi/dis/dis
   end do
   fvx(i)=U*cos(alpha)-fvx(i)
   fvy(i)=U*sin(alpha)-fvy(i)
end do
!-----------------------顆粒運動計算--------------------------

do i=1,num
   do j=1,6
      k1(j)=dfunc(ti,pvx(i),pvy(i),pvz(i),px(i),py(i),pz(i),j)
   end do
   do j=1,6
      k2(j)=dfunc(ti+dt/2,pvx(i)+k1(1)*dt/2,pvy(i)+k1(2)*dt/2,pvz(i)+k1(3)*dt/2,px(i)+k1(4)*dt/2,py(i)+k1(5)*dt/2,pz(i)+k1(6)*dt/2,j)
   end do
   do j=1,6
      k3(j)=dfunc(ti+dt/2,pvx(i)+k2(1)*dt/2,pvy(i)+k2(2)*dt/2,pvz(i)+k2(3)*dt/2,px(i)+k2(4)*dt/2,py(i)+k2(5)*dt/2,pz(i)+k2(6)*dt/2,j)
   end do
   do j=1,6
      k4(j)=dfunc(ti+dt/2,pvx(i)+k3(1)*dt/2,pvy(i)+k3(2)*dt/2,pvz(i)+k3(3)*dt/2,px(i)+k3(4)*dt/2,py(i)+k3(5)*dt/2,pz(i)+k3(6)*dt/2,j)
   end do
   pvx0=pvx(i)
   pvy0=pvy(i)
   pvz0=pvz(i)
   px0=px(i)
   py0=py(i)
   pz0=pz(i)
   pvx(i)=pvx(i)+dt*(k1(1)+2*k2(1)+2*k3(1)+k4(1))/6
   pvy(i)=pvy(i)+dt*(k1(2)+2*k2(2)+2*k3(2)+k4(2))/6
   pvz(i)=pvz(i)+dt*(k1(3)+2*k2(3)+2*k3(3)+k4(3))/6
   px(i)=px(i)+dt*(k1(4)+2*k2(4)+2*k3(4)+k4(4))/6
   py(i)=py(i)+dt*(k1(5)+2*k2(5)+2*k3(5)+k4(5))/6
   pz(i)=pz(i)+dt*(k1(6)+2*k2(6)+2*k3(6)+k4(6))/6
   ti=ti+dt
end do