mover.f

Fortran的数学物理方程数值算法源程序。这是"Numerical Methods for Physics"第二版的Fortran源程序

      subroutine mover( x, v, npart, L, mpv, vwall, tau,
     &                  strikes, delv, seed )
      integer*4 MAXnpart, MAXncell
      parameter( MAXnpart = 10000, MAXncell = 500 )
      integer*4 npart, seed, strikes(2)
      real*8 x(MAXnpart), v(MAXnpart,3), L, mpv, vwall, tau, delv(2)

! mover - Function to move particles by free flight
!         Also handles collisions with walls
! Inputs
!    x        Positions of the particles
!    v        Velocities of the particles
!    npart    Number of particles in the system
!    L        System length
!    mpv      Most probable velocity off the wall
!    vwall    Wall velocities
!    tau      Time step
!    seed     Random number seed
! Outputs
!    x,v      Updated positions and velocities
!    strikes  Number of particles striking each wall
!    delv     Change of y-velocity at each wall
!    seed     Random number seed

      integer*4 i, flag
      real*8 x_old(MAXnpart), xwall(2), vw(2), direction(2), stdev
      real*8 vyInitial, dtr
      real*8 rand, randn

      !* Move all particles pretending walls are absent
      do i=1,npart
        x_old(i) = x(i)            ! Remember original position
        x(i) = x_old(i) + v(i,1)*tau
      enddo

      !* Check each particle to see if it strikes a wall
      strikes(1) = 0
      strikes(2) = 0
      delv(1) = 0.0
      delv(2) = 0.0
      xwall(1) = 0
      xwall(2) = L    ! Positions of walls
      vw(1) = -vwall
      vw(2) = vwall   ! Velocities of walls
      stdev = mpv/sqrt(2.0)

      direction(1) = 1
      direction(2) = -1  ! Direction of particle leaving wall
      do i=1,npart

        !* Test if particle strikes either wall
        flag = 0
        if( x(i) .le. 0 ) then
          flag=1        ! Particle strikes left wall
        else if( x(i) .ge. L ) then
          flag=2        ! Particle strikes right wall
        endif

        !* If particle strikes a wall, reset its position
        !  and velocity. Record velocity change.
        if( flag .ne. 0 )  then
          strikes(flag) = strikes(flag) + 1
          vyInitial = v(i,2)
          !* Reset velocity components as biased Maxwellian,
          !  Exponential dist. in x; Gaussian in y and z
          v(i,1) = direction(flag)*sqrt(-dlog(1.0-rand(seed))) * mpv
          v(i,2) = stdev*randn(seed) + vw(flag)  ! Add wall velocity
          v(i,3) = stdev*randn(seed)
          ! Time of flight after leaving wall
          dtr = tau*(x(i)-xwall(flag))/(x(i)-x_old(i))
          !* Reset position after leaving wall
          x(i) = xwall(flag) + v(i,1)*dtr
          !* Record velocity change for force measurement
          delv(flag) = delv(flag) + (v(i,2) - vyInitial)
        endif
      enddo

      return
      end