Nugroho's blog.: Forced to get Force using Center of Mass

Saturday, April 16, 2016

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()


.


No comments:

323f (5) amp (1) android (12) apple (7) arduino (18) art (1) assembler (21) astina (4) ATTiny (23) blackberry (4) camera (3) canon (2) cerita (2) computer (106) crazyness (11) debian (1) delphi (39) diary (286) flash (8) fortran (6) freebsd (6) google apps script (8) guitar (2) HTML5 (10) IFTTT (7) Instagram (7) internet (12) iOS (5) iPad (6) iPhone (5) java (1) javascript (1) keynote (2) LaTeX (6) lazarus (1) linux (29) lion (15) mac (28) macbook air (8) macbook pro (3) macOS (1) Math (3) mathematica (1) maverick (6) mazda (4) microcontroler (35) mountain lion (2) music (37) netbook (1) nugnux (6) os x (36) php (1) Physicist (29) Picture (3) programming (189) Python (109) S2 (13) software (7) Soliloquy (125) Ubuntu (5) unix (4) Video (8) wayang (3) yosemite (3)