An N-body College Reunion

原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/an-n-body-college-reunion

9007 Rate this article:

5.0

An N-body College Reunion

Jim Pendleton Thursday, May 21, 2015

A little over 30 years ago, the professor in my sophomore quantum mechanics class gave us a free range assignment. We could use a computer to solve any sort of physics problem. It need not be related to quantum in any way, which was just as well since I've never grokked the wave function in fullness.

This was back in the day when students training for careers in basic sciences were expected to know how to write their own analysis applications, and not mindlessly press buttons on spreadsheets or web pages. In no small part because the latter hadn't been invented yet.

I chose a problem in astrophysics. At the time there was academic research being performed on simplified galactic collisions, and the plots in the scientific journals were fascinating to me. N-body problems were well beyond my grasp, unless N was one or two.  

I had no access to supercomputers, and my classical mechanics training was only at the freshman level.  

But I liked pretty pictures.

It was challenging at the time to debug Fortran PLOT10 programs when one needed to plan for additional wait time in line at the campus computing center for the output and it took an hour for a Calcomp plotter to produce the (garbage) graph the (garbage) code was generating.

Unfortunately, my attempt burned down, fell over, and sank into the swamp. I vaguely recall I received an orange slice and a participant’s ribbon for the effort. Thanks for not failing me, Prof. Dutta.

But a decade later, and with access to IDL, I reattempted the feat during a competition to produce demos for IDL 5.0.

The results were more satisfying visually, if not more realistic physically.

I've left the code below more or less intact. It’s a main application.  Dig the DEVICE, DECOMPOSED=0 and parentheses instead of square brackets – 100% OG, and I don’t mean Object Graphics.

Essentially, two point masses represent galactic centers and 1000 mass-less test points (“stars”) are given initial circular paths around those centers. Both galactic centers then move relative to their center of mass with some initial conditions intended to provoke a visually appealing, if physically bogus, result.

I wrote a number of variants on this theme, color-coded the stars by speed relative to the center of mass, added other “galactic” masses, etc. If you find this interesting, try varying initial parameters such as the relative masses and velocities.

; initialization section

G = 4.4999e-12 ; pc^3 y^-2 Msol^-1

epsilon = 1.e5 ; yr

snapshot = 16.e6 ;yr

nsnapshots = 16

maxradius = 15d3 ; pc

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 ; Msol

x(0) = initialsep/2.*cos(!pi/4.)

y(0) = initialsep/2.*sin(!pi/4.)

m(1) = 9.d9 ;Msol

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.

; execution section

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

;

;   Erase the previous frame

;

closedcircles,sym=1

oplot,oldx(0:NGalaxies-1),oldy(0:NGalaxies-1), color = 0

!psym = 3

oplot, oldx, oldy, color = 0

;

;   Draw the galaxy centers

;

closedcircles,sym=1

for i = 0, NGalaxies - 1 Do Begin

oplot,[x(i),x(i)],[y(i),y(i)],color=galaxycolors(i)

endfor

;

;   Draw the stars

;

!psym = 3

for i = 0, NGalaxies - 1 Do Begin

indices = indgen(nsuns) + (i*nsuns) + Ngalaxies

oplot, x(indices), y(indices), color=galaxycolors(i)

endfor

;

;   Save the coordinates

;

oldx = x

oldy = y

empty

endelse

endfor

end

You will also need the following support routine.

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