书中的图有些误导人，因为那是一个中心排斥势的散射图，而Lennard-Jones势函数在散射中应该是起到了吸引的作用。假设Lennard-Jones势函数中的能量是量子化的。当入射粒子能量和瞄准距离都合适时，可能被俘获。也就是形成共振现象。在此处所谓的奇点便是指的这个。当粒子接近到
此时，径向能量为：
那么，
也就是书中提到的，被积分式接近零，使得积分出现奇点。
书中给出的公式为：
其中第一项的奇点是
点击(此处)折叠或打开
- !
- 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
当碰撞参数很小时，出现大角度散射现象；当碰撞参数变大时，散射角逐渐减小。因此计算结果也是符合物理直觉的。