6
5. 绘制两个单位异号、同号点电荷在其所在平面的电势分布，两个点电荷分别在x轴线上处于x1??2和x2?2的位置。同号电荷要求画出V=3，2，1，0.75，0.5等势线，异号电荷画出V=3，2，1，0.75，0.5，-0.5,－0.75，－1，－2，－3的等势线，并在图中用圆圈标出正电荷，十字为负电荷（令同种电荷：
14??0?1）。
异种电荷：
7
6. 利用蒙特卡罗方法计算定积分
，请给出其求解原理与计算
步骤（具体模拟方法、计算程序和结果）。 程序如下：
program main implicit none
integer k,m,n,idum real x,y,s
real,external::ran1 n=100000 idum=-1 m=0
do k=1,n
x=ran1(idum) y=ran1(idum)
if(x*x-y>0) then m=m+1 endif end do
s=4.*float(m)/float(n) write(*,*) 's=',s end
8
FUNCTION ran1(idum)
INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL ran1,AM,EPS,RNMX PARAMETER
(IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy
DATA iv /NTAB*0/, iy /0/
if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum
ran1=min(AM*iy,RNMX) return END
结果为： s = 1.342240
7．模拟氢原子基态（n=1,l=0,m=0）的电子云。氢原子基态电子分布
几率密
度
函数是
D?4ra123e?2r/a1，其中a1?5.29?10?2nm,D
的最大值Dmax?1.1，r0?0.25nm是D的收敛点（给出模拟方法和模拟结果）。
模拟方法：1.模拟时先产生一个随机的电子轨道半径r = r0rand(1), 显
9
然有0 ≤r ≤ r0, 由r计算出D(r)。
2. 再产生一个随机的概率判据D0 = Dmaxrand(1), 显然 有0≤D0≤Dmax,
3. 然后进行判断，如果D(r) < D0,则舍弃它，反 之就计算一个随机的角度值，θ = 2лrand(1), 4. 最后得到的点的坐标是x = r cos θ; y = r sin θ。 程序如下：
program main implicit none integer n,idum,k
real x,y,r,r0,dr,d0,theta,a real,external :: ran1
open(2,file='F:\\fortran\\rand.txt')
n = 60000 idum = -1 r0 = 0.25 a = 0.0529
do k = 1, n
r = r0 * ran1(idum)
dr = 4./a**3*r**2*exp(-2/a*r) d0 = 1.1*ran1(idum) if(dr .gt. d0)then
theta = 2*acos(-1.)*ran1(idum) x = r*cos(theta) y = r*sin(theta) write(2,*) x,y endif end do
end program
FUNCTION ran1(idum)
INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL ran1,AM,EPS,RNMX PARAMETER
(IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
10
INTEGER j,k,iv(NTAB),iy SAVE iv,iy
DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ
11 idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM continue if (j.le.NTAB) iv(j)=idum iy=iv(1) endif
k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum
ran1=min(AM*iy,RNMX)
return END
11