Forced to get Force using Center of Mass

 I know. Using solely center of mass as main contributor in force calculation is bad idea in nearly homogenous n-body system.

 Well, I'm curious about that. So, just bite it.

 And, yeah, bad idea.

 The system is fine, generally.

 Each body move as usual, at a glimpse.

 But if we go detail about that. There's strange behavior here and there.


 On previous code with brute force calculation, we could saw, binary system consist of two body orbiting each other. But, with this center-of-mass-code, we see no one. A body that have very close neighbor don't make a twin system. It continues orbiting center of mass.

 As we have system with -body with same mass, it should be natural if two close-body will have greater gravity force than with any other body. It should have "priority"; orbiting each other.

 Ok, this "forced"-center-mass-code is weird, but it'll useful in other system. :)

from visual import *
from random import uniform,random

l       = 17.
dl      = .01
display(center=(0,0,0),background=(1,1,1),autoscale=False, width=600, height=600,forward=(-0.7,-0.7,-1))
distant_light(direction=(1,1,1), color=color.red)
#local_light(pos=(0,0,0), color=(0,0,1))
#sphere(pos=(0,0,0), color=(0,0,0), material=materials.emissive, opacity=0.9)

box(color=color.white, pos=(0,l/2,0),length=l,height=dl, width=l, opacity=0.3)
box(color=color.white, pos=(0,-l/2,0),length=l,height=dl, width=l, opacity=0.3)
box(color=color.white, pos=(l/2,0,0),length=dl,height=l, width=l, opacity=0.3)
box(color=color.white, pos=(-l/2,0,0),length=dl,height=l, width=l, opacity=0.3)
box(color=color.white, pos=(0,0,l/2),length=l,height=l, width=dl, opacity=0.3)

bola    = [] 
n       = 11
dv      = .1        #kecepatan setelah menabrak dinding
G       = 31.
for i in range(n):
    ball            = sphere (pos=(uniform(-7,7),uniform(-7,7),uniform(-7,7)), radius=.2, color=(random(),random(),random()))
    ball.v          = vector(uniform(-1,1),uniform(-1,1),uniform(-1,1))
    
    bola.append(ball)


dt = 1./32.
bola[0].pos=(1,1,1)


def hitungGaya(i,rcm):
    r           = bola[i].pos-rcm
    jarak       = mag(r)
    arah        = -norm(r)
    gaya        = G*1./(pow(jarak,2)+.1)*arah*n
    return gaya
def hitungPusatMassa():
    xcm = 0.
    ycm = 0.
    zcm = 0.
    for i in arange(n):
        xcm += bola[i].x
        ycm += bola[i].y
        zcm += bola[i].z
    xcm /= n
    ycm /= n
    zcm /= n
    return vector(xcm,ycm,zcm)
        
def proses():
    global bola
    rcm = hitungPusatMassa()
    for i in arange(n):
        gaya        = hitungGaya(i,rcm)
        r           = bola[i].pos
        #tambah gaya dari pusat
        '''
        jarak       = mag(bola[i].pos)
        arah        = -norm(bola[i].pos)
        gaya        += 10*G*1./(pow(jarak,2)+.05)*arah
        '''
        a           = gaya
        v           = bola[i].v
        v           += a*dt
        r           += v*dt
        #cek pantul ke tembok
        if bola[i].x > l/2.:
            bola[i].x   = l/2.
            v.x         *= -dv
        elif bola[i].x < -l/2.:
            bola[i].x   = -l/2.
            v.x         *= -dv
        elif bola[i].y > l/2.:
            bola[i].y   = l/2.
            v.y         *= -dv
        if bola[i].y < -l/2.:
            bola[i].y   = -l/2.
            v.y         *= -dv
        if bola[i].z > l/2.:
            bola[i].z   = l/2.
            v.z         *= -dv
        if bola[i].z < -l/2.:
            bola[i].z   = -l/2.
            v.z         *= -dv

        bola[i].pos = r
        
while 1:
    rate (10)
    proses()

.