Particle Simulation

particles.f03

program particles
    ! ten particles in a field of 10000 x 10000
    ! base velocity is 10 x 10 

    integer, parameter      :: p = 10,                         &
                               fx = 10000, fy = 10000,         &
                               bx = 10, by = 10,               &
                               steps = 10000
    real, parameter         :: radius = 5
    real, dimension(p, 2)   :: position, velocity
    real                    :: num
     
    ! initialize the random number seed
    call reseed()    

    do i = 1, p
        call random_number(num)
        position(i, 1) = num * fx
		
        call random_number(num)
        position(i, 2) = num * fy

        call random_number(num)
        velocity(i, 1) = num * bx

        call random_number(num)
        velocity(i, 2) = num * by
        
        print *, "Particle ", i, " at (", position(i, 1), ", ", position(i, 2),") heading (",   &
                                          velocity(i, 1), ", ", velocity(i, 2), ")"
    end do
    
    call step(position, velocity)

contains

    subroutine reseed()
        integer                             :: i, n, clock
        integer, dimension(:), allocatable  :: seed

        call random_seed(size = n)
        allocate(seed(n))
        
        call system_clock(count = clock)
        
        seed = clock + 37 * (/ (i - 1, i = 1, n) /)
        call random_seed(put = seed)
        
        deallocate(seed)
    end subroutine reseed

    subroutine step(positions, velocities)
        real, dimension(p,2), intent(inout) :: positions, velocities
        
        ! verify lengths (? unnecessary?)
        if (size(positions, 1) .ne. size(velocities, 1)) then
            print *, "The sizes of two arrays are not the same"
        end if
        
        ! move the particles to the new position 
        do i = 1, p
            positions(i, :) = positions(i, :) + velocities(i, :)
!            positions(i, 1) = positions(i, 1) + velocities(i, 1)
!            positions(i, 2) = positions(i, 2) + velocities(i, 2)
        end do
               
        ! check for collisions :: only works for binary collisions...
        do 10 i = 1, p
            do 20 j = i, p
                if (collision(positions(i, :), positions(j, :))) then
                    
                        
20          continue
10      continue      
      
    end subroutine step
    
    function distance(p1, p2)
        real, dimension(2), intent(in) :: p1, p2, origin
        real                           :: distance
        
        origin = p2 - p1
        distance = sqrt(origin(1)**2 + origin(2)**2)
    end function distance 
    
    function collision(p1, p2)
        real, dimension(2), intent(in) :: p1, p2
        logical                        :: collision
        
        collision = (distance(p1, p2) .le. radius)
    end function collision
        

end program particles