N体问题的同学会
原文链接：https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/an-n-body-college-reunion
9007 给本文评分：
5.0
N体问题的同学会
Jim Pendleton 2015年5月21日 星期四
大约30多年前，在我大二的量子力学课上，教授给我们布置了一项自由选题作业。我们可以用计算机解决任何物理问题，完全不必与量子力学相关，这对我来说正好，因为我从未完全理解过波函数。
在那个年代，为从事基础科学事业而受训的学生，都被期望懂得如何编写自己的分析应用，而不是盲目地点击电子表格或网页上的按钮。这在很大程度上是因为后者当时尚未被发明出来。
我选择了一个天体物理学问题。当时有关于简化星系碰撞的学术研究在进行，科学期刊上的那些图表让我着迷。N体问题远超出我的理解范围，除非N等于1或2。
我那时接触不到超级计算机，而我的经典力学训练也只停留在大一水平。
但我喜欢漂亮的图片。
那时候调试 Fortran PLOT10 程序颇具挑战性，因为你需要在校园计算中心排队等待额外的时间才能看到输出结果，而 Calcomp 绘图仪则需要一个小时才能生成（糟糕的）代码所产生的（糟糕的）图表。
不幸的是，我的尝试以失败、崩溃和沉入沼泽而告终。我隐约记得我因此获得了一片橙子和一张参与奖绶带。感谢您没有让我挂科，Dutta教授。
但十年后，借助IDL，在为IDL 5.0制作演示程序的竞赛中，我重新尝试了这一壮举。
结果在视觉上更加令人满意，尽管在物理上可能更不真实。
我大致原封不动地保留了下面的代码。这是一个主应用程序。注意里面的 DEVICE, DECOMPOSED=0 和圆括号而非方括号——这是100%的原版代码，我指的不是面向对象图形。
本质上，两个质点代表星系中心，1000个无质量测试点（“恒星”）被赋予围绕这些中心的初始圆形轨道。然后，两个星系中心在它们的质心系中运动，具有一些旨在产生视觉吸引人（尽管在物理上可能站不住脚）的结果的初始条件。
我基于这个主题写了很多变体，根据相对于质心的速度为恒星着色，添加了其他“星系”质量等等。如果你对此感兴趣，可以尝试改变相对质量和速度等初始参数。
; 初始化部分
G = 4.4999e-12 ; 秒差距^3 年^-2 太阳质量^-1
epsilon = 1.e5 ; 年
snapshot = 16.e6 ; 年
nsnapshots = 16
maxradius = 15d3 ; 秒差距
initialsep = 2.*maxradius
Ngalaxies = 2
GalaxyColors = [150, 200]
Nsuns = 1000
vx = dblarr(Ngalaxies*(Nsuns + 1)) & vy = vx & vz = vx & x = vx & y = vx
z = vx & m = vx & rdist= vx & thetadist = vx & ax = vx & ay = vx & az = vx
m(0) = 10.d9 ; 太阳质量
x(0) = initialsep/2.*cos(!pi/4.)
y(0) = initialsep/2.*sin(!pi/4.)
m(1) = 9.d9 ; 太阳质量
x(1) = -x(0)
y(1) = -y(0)
speed = 2*sqrt(G*m(1)/initialsep)
vx(0) = -speed*sin(!pi/4.)/2.
vy(0) = speed*cos(!pi/4.)/2.
vx(1) = -vx(0)*sqrt(m(0)/m(1))/3.
vy(1) = -vy(0)*sqrt(m(0)/m(1))/3.
; 执行部分
device, decomposed = 0
iseed = 138987
for i = Ngalaxies, Ngalaxies*(Nsuns + 1) - 1 do begin
rdist(i) = randomu(iseed) * maxradius
thetadist(i) = randomu(iseed)*2*!pi
m(i) = 1.
endfor
for j = 0, ngalaxies - 1 do begin
for i = ngalaxies + j*nsuns, ngalaxies + (j + 1)*nsuns - 1 do begin
x(i) = rdist(i)*cos(thetadist(i)) + x(j)
y(i) = rdist(i)*sin(thetadist(i)) + y(j)
z(i) = 0.
speed = sqrt(G*m(j)/max([Abs(rdist(i)),1.e-30]))
vx(i) = -speed*sin(thetadist(i)) + vx(j)
vy(i) = speed*cos(thetadist(i)) + vy(j)
vz(i) = 0.
endfor
endfor
loadct,15
!ignore = -1
!x.style = 2
!y.style = 2
window, 1, xsize=700, ysize=700, xpos = 0, ypos = 0
wset, 1
for t = -epsilon, snapshot*(nsnapshots - 1), epsilon do begin
ax(*) = 0 & ay = ax & az = ax
for j = 0, NGalaxies - 1 do begin
r = sqrt((x - x(j))^2 + (y - y(j))^2 + (z - z(j))^2)
r3 = r^3 > 1.e-20
Gprod = -G*m(j)/r3
ax = Gprod*(x-x(j)) + ax
ay = Gprod*(y-y(j)) + ay
az = Gprod*(z-z(j)) + az
endfor
vx = vx + epsilon*ax
vy = vy + epsilon*ay
vz = vz + epsilon*az
x = x + epsilon*vx
y = y + epsilon*vy
z = z + epsilon*vz
if (t eq -epsilon) then begin
rs = 1.5*initialsep
plot, [rs, -rs], [rs, -rs], color = 0
closedcircles,sym=1
for i = 0, NGalaxies - 1 Do Begin
oplot,[x(i),x(i)],[y(i),y(i)],color=galaxycolors(i)
endfor
!psym = 3
for i = 0, NGalaxies - 1 Do Begin
indices = indgen(nsuns) + (i*nsuns) + Ngalaxies
oplot, x(indices), y(indices), color=galaxycolors(i)
endfor
oldx = x
oldy = y
endif else begin
;
; 擦除前一帧
;
closedcircles,sym=1
oplot,oldx(0:NGalaxies-1),oldy(0:NGalaxies-1), color = 0
!psym = 3
oplot, oldx, oldy, color = 0
;
; 绘制星系中心
;
closedcircles,sym=1
for i = 0, NGalaxies - 1 Do Begin
oplot,[x(i),x(i)],[y(i),y(i)],color=galaxycolors(i)
endfor
;
; 绘制恒星
;
!psym = 3
for i = 0, NGalaxies - 1 Do Begin
indices = indgen(nsuns) + (i*nsuns) + Ngalaxies
oplot, x(indices), y(indices), color=galaxycolors(i)
endfor
;
; 保存坐标
;
oldx = x
oldy = y
empty
endelse
endfor
end
你还需要以下辅助程序。
Pro ClosedCircles, Symbol_Radius = Symbol_Radius, Color = Color
On_Error, 2
Ang = (360./47.*Findgen(48))/!RaDeg
X = Cos(Ang)
Y = Sin(Ang)
If (N_Elements(Symbol_Radius) ne 0) then Begin
X = X * Symbol_Radius
Y = Y * Symbol_Radius
EndIf
If (N_Elements(Color) eq 0) then Begin
UserSym, X, Y, /Fill
EndIf Else Begin
UserSym, X, Y, /Fill, Color = Color
EndElse
!P.PSym = 8
Return
End
User-Defined ENVITask That Masks Some Output Pixels Computer Based Training