Computational Physics: Fortran Version——Project I Lennard-Jones中心势散射

　　这个项目给人的启发也是不小的，例如在此处截断和关于奇点的讨论。尤其是后一点，很有意思，在此详细论述一下。
　　书中的图有些误导人，因为那是一个中心排斥势的散射图，而Lennard-Jones势函数在散射中应该是起到了吸引的作用。假设Lennard-Jones势函数中的能量是量子化的。当入射粒子能量和瞄准距离都合适时，可能被俘获。也就是形成共振现象。在此处所谓的奇点便是指的这个。当粒子接近到时，被俘获。此时，粒子能量分为径向能量和切向能量。被俘获时，径向能量与在该距离下的Lennard-Jones势相抵消，而只剩余切向能量。从而粒子只能绕行很多圈。
　　此时，径向能量为：
　　
　　那么，,两边再同时除以E,即可得到如下关系式：
　　
　　也就是书中提到的，被积分式接近零，使得积分出现奇点。
　　书中给出的公式为：
　　
　　其中第一项的奇点是，即下域；第二项的奇点便是
　　

!
module param
  implicit none
  real,parameter :: pi=3.1415926
  real,parameter :: E=1.0 ! incident energy in unit of V0
  real,parameter :: Bmin=0.1 ! minimum of impact parameter
  real,parameter :: Bmax=2.4 ! maximum of impact parameter
  real,parameter :: Rmax=2.5 ! maximum of radius
  real,parameter :: Tolr=0.0001 ! tolerance for turning point search
  integer,parameter :: Nb=20 ! number of steps in impact parameter b
  integer,parameter :: Npts=40 ! number of quadrature points
  real,parameter :: db=(Bmax-Bmin)/(Nb-1)
  integer,parameter :: BnumMax=100

end module
!

!
Program main
! Project 1: Scattering by the Lennard-Jones potential
!
  use param
  call archon

stop

end program
!

!
subroutine archon
!
! Global variables
  use param
! Local variables
  implicit none
  real :: Impact(0:BnumMax)=0.0
  integer :: Ib=0
  real :: Angle(0:BnumMax)=0.0
  real :: Int1,Int2,Rmin

  do Ib=0,(Nb-1),1
    Impact(Ib)=Bmin+Ib*db
    call Theta(Impact(Ib),Int1,Int2,Rmin)
    write(*,*)"Int1:",Int1,"Int2:",Int2
    Angle(Ib)=2.0*Impact(Ib)*(Int1-Int2)*180.0/pi
    write(*,*)"b(",Ib,"):",Impact(Ib),"Angle(",Ib,"):",Angle(Ib)
  end do

return

end subroutine
!

!
subroutine Theta(b,Int1,Int2,Rmin)
!
! Global variables
  use param
! Input/Output variables
  real :: b ! impact parameter (Input)
  real :: Int1 ! the 1st integral (Output)
  real :: Int2 ! the 2nd integral (Output)
  real :: Rmin ! radius of the closes approach (Output)
! Local variables
  real :: r=1.0
  real :: dr=0.2
  real,external :: pot,fun,fun1,fun2
  real :: u
  real :: total=0.0
  integer :: Iu
  real :: umax
  real :: step

  Rmin=Rmax

  do while(dr>Tolr)
    Rmin=Rmin-dr
    if(fun(Rmin,b)<0)then
      Rmin=Rmin+dr
      dr=dr/2.0
    end if
  end do
!
  umax=sqrt(Rmax-b)
  step=umax/Npts
  total=0.0
  do Iu=1,Npts,1
    u=step*(Iu-0.5)
    r=u**2+b
    total=total+u*fun1(r,b)
  end do
  Int1=2.0*step*total
!
  total=0.0
  umax=sqrt(Rmax-Rmin)
  step=umax/Npts
  do Iu=1,Npts,1
    u=step*(Iu-0.5)
    r=u**2+Rmin
    total=total+u*fun2(r,b)
  end do
  Int2=2.0*step*total

return
end subroutine
!

!
real function pot(r)
  implicit none
  real :: r
  pot=4.0*(1/r**12-1/r**6)
  return
end function
!

!
real function fun(r,b)
  use param
  implicit none
  real :: r,b
  real,external :: pot
  fun=1-b**2/r**2-pot(r)/E
return
end function
!

!
real function fun1(r,b)
  implicit none
  real :: r,b
  fun1=1/r**2*(1/sqrt(1-b**2/r**2))
return
end function
!

!
real function fun2(r,b)
  use param
  implicit none
  real :: r,b
  real,external :: pot
  fun2=1/r**2*(1/sqrt(1-b**2/r**2-pot(r)/E))
return
end function
!

　　得到的计算结果为：
 b(0):  0.100000001     Angle(0):   168.910172    
 b(1):  0.221052647     Angle(1):   169.855896    
 b(2):  0.342105299     Angle(2):   164.271286    
 b(3):  0.463157952     Angle(3):   158.649170    
 b(4):  0.584210575     Angle(4):   152.974335    
 b(5):  0.705263257     Angle(5):   147.230621    
 b(6):  0.826315880     Angle(6):   141.401428    
 b(7):  0.947368503     Angle(7):   135.465668    
 b(8):   1.06842113     Angle(8):   129.401596    
 b(9):   1.18947387     Angle(9):   123.180473    
 b(10):   1.31052649     Angle(10):   116.771988    
 b(11):   1.43157911     Angle(11):   110.134575    
 b(12):   1.55263174     Angle(12):   103.215302    
 b(13):   1.67368436     Angle(13):   95.9486313    
 b(14):   1.79473698     Angle(14):   88.2398682    
 b(15):   1.91578972     Angle(15):   79.9521637    
 b(16):   2.03684235     Angle(16):   70.8765564    
 b(17):   2.15789509     Angle(17):   60.6551018    
 b(18):   2.27894759     Angle(18):   48.5505600    
 b(19):   2.40000033     Angle(19):   32.5175247 
　　当碰撞参数很小时，出现大角度散射现象；当碰撞参数变大时，散射角逐渐减小。因此计算结果也是符合物理直觉的。