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

.

Menjinakkan Gravitasi.

 Cari cara agar mereka jadi enak untuk dilihat.

 Di kode sebelumnya, bola sudah dicegah agar tidak melarikan diri dengan menambahkan kotak, mereka akan memantul jika menabrak kotak.

 Beres, sekarang tambahkan interaksi antar bola, interaksi sederhana, gravitasi.

 Ok, pake metode brutal, untuk tiap-tiap bola, hitung satu-satu interaksi gravitasi dengan seluruh bola, kemudian jumlahkan.





(ide, cari pusat masa seluruh bola, perlakukan interaksi gaya seperti kode sebelumnya)

 (saat nulis tadi, sambil mikir, ternyata itu ide buruk, tapi boleh juga, lihat saja nanti, bagus atau tidak)

 (tapi tidak sekarang, tidak masuk kode yang sekarang)

 (kode di posting ini tetap pake metode brutal)

 Hasilnya adalah, ..., ok, mereka tetap memantul saat menabrak tembok, tetapi di dalam kotak mereka bergerak dengan sangat cepat, bola liar.

 (yep, metode brutal menghasilkan hasil brutal, tetapi tidak selalu...)

 Ide licik, tambahkan koefisien restitusi di dinding, tumbukan dengan dinding tidak lenting sempurna, kecepatannya tinggal 10% setelah memantul

 Dengan demikian mereka lebih jinak di dalam kotak, :)

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=.3, 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):
    gaya    = vector(0,0,0)
    for j in arange(n):
        if i!=j:  
            r           = bola[i].pos-bola[j].pos
            jarak       = mag(r)
            arah        = -norm(r)
            gaya        += G*1./(pow(jarak,2)+.1)*arah
    return gaya
def proses():
    global bola
    for i in arange(n):
        gaya        = hitungGaya(i)
        #tambah gaya dari pusat
        r           = bola[i].pos
        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()

.

Bermain dengan Gaya Sentral di Python.

 Merangkak dari kode sebelumnya, dengan tambahan kotak yang mengungkung bola-bola yang mengitari titik asal dengan berbagai kecepatan. Bola-bola tersebut tak bisa menembus kotak, jika menabrak salah satu dinding kotak, maka dia akan memantul

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)

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

G       = 31.
for i in range(n):
    ball            = sphere (pos=(uniform(1,7),0,uniform(-7,7)), radius=.3, color=(random(),random(),random()))
    ball.v          = vector(0,uniform(1,3),0)
    
    bola.append(ball)


dt = 1./32.
bola[0].pos=(1,1,1)
while 1:
    rate (100)
    for i in arange(n):
        r           = bola[i].pos
        v           = bola[i].v
        jarak       = mag(bola[i].pos)
        arah        = -norm(bola[i].pos)
        gaya        = G*1./(pow(jarak,2)+.1)*arah
        a           = gaya
        v           += a*dt
        r           += v*dt
        #cek pantul ke tembok
        if bola[i].x > l/2.:
            bola[i].x   = l/2.
            v.x         *= -1
        elif bola[i].x < -l/2.:
            bola[i].x   = -l/2.
            v.x         *= -1
        elif bola[i].y > l/2.:
            bola[i].y   = l/2.
            v.y         *= -1
        if bola[i].y < -l/2.:
            bola[i].y   = -l/2.
            v.y         *= -1
        if bola[i].z > l/2.:
            bola[i].z   = l/2.
            v.z         *= -1
        if bola[i].z < -l/2.:
            bola[i].z   = -l/2.
            v.z         *= -1

        bola[i].pos = r
        

.

Playing with Central Force in Pyton with Visual Module

 I don't use n body calculation, every sphere only attracted by force from origin, dependent upon its position.

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.4,-0.3,-1))
distant_light(direction=(1,1,1), color=color.red)
box(color=color.white, pos=(0,0,0),length=l,height=dl, width=l, opacity=0.3)
box(color=color.white, pos=(0,0,0),length=dl,height=l, width=l, opacity=0.3)
box(color=color.white, pos=(0,0,0),length=l,height=l, width=dl, opacity=0.3)

bola    = [] 
n       =   11

G       = 31.
for i in range(n):
    ball            = sphere (pos=(uniform(1,7),0,uniform(-7,7)), radius=.3, color=(random(),random(),random()))
    ball.v          = vector(0,uniform(1,3),0)
    
    bola.append(ball)


dt = 1./32.
bola[0].pos=(1,1,1)
while 1:
    rate (100)
    for i in arange(n):
        r           = bola[i].pos
        v           = bola[i].v
        jarak       = mag(bola[i].pos)
        arah        = -norm(bola[i].pos)
        gaya        = G*1./(pow(jarak,2)+.1)*arah
        a           = gaya
        v           += a*dt
        r           += v*dt

        bola[i].pos = r
        

.

Wave Simulation Using Visual Python.

 I used previous code (with matplotlib animation).

 Removed the matplotlib part, swap it with vpython, :)

"""
Gelombang
"""
import numpy as np                          #untuk operasi array
from visual import *

#variabel
n   = 39
x   = np.arange(0., 1.,1./n )   #array dari 0 s.d 1 berjarak 1/n
y   = np.zeros(n)               #array sejumlah n isinya 0
y1  = np.zeros(n)               # y1 : array [0..n] of real
y2  = np.zeros(n)
y1  = np.exp(-1*np.power(10*x-3,2))
y2  = np.exp(-1*np.power(10*x-3,2))

r2  = 1./512   

display(center=(.5,0,0),background=(1,1,1))
#kotak
dindingKiri   = box (pos=(0.,0,0), length=.01, height=1, width=1, color=color.green)
dindingKanan   = box (pos=(1.,0,0), length=.01, height=1, width=1, color=color.blue)
bola    = []
for i in range(n):
    ball            = sphere (pos=(x[i],y2[i],0), radius=.01, color=color.red)
    bola.append(ball)

#nilai awal fungsi gaussian

#print y2


def proses():
    #hitung nilai baru
    for i in np.arange(1,n-1):
        y[i]    =   2*(1-r2)*y1[i]-y2[i]+r2*(y1[i+1]+y1[i-1])
    #geser
    y2[:]   =   y1[:]
    y1[:]   =   y[:]
    #print y2

    return y2

while 1:
    rate(100)
    proses()
    for i in np.arange(n):
        bola[i].y   = y2[i]

.

                          

Hello (Again) Visual Python, :D .

Menyapa kembali mainan lama, :)

from visual import *

floor   = box (pos=(0,0,0), length=4, height=0.5, width=4, color=color.blue)
bola    = [] 
n   =   2
for i in range(n):
    ball            = sphere (pos=(0,4,0), radius=1, color=color.red)
    ball.v          = vector(0,-1,0)

    bola.append(ball)

dt = 0.01
bola[1].pos=(1,1,1)
while 1:
    rate (100)
    bola[0].pos = bola[0].pos + bola[0].v*dt
    if bola[0].y < bola[0].radius:
        bola[0].v.y = abs(bola[0].v.y)
    else:
        bola[0].v.y = bola[0].v.y - 9.8*dt

.

 

Lorenz Attractor 3D Scatter Plot using Python with Matplotlib

 After failed with regular plot() syntax, I have luck with the scatter() one, :)

"""
Cluster
"""
import numpy as np                          #untuk operasi array
import matplotlib.pyplot as plt             #untuk gambar grafik
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation    #untuk menggerakkan grafik
from itertools import cycle
randColor = cycle('bgrcmk').next
randMarker= cycle('.,ov^<>12348sp*hH+xDd|_').next

#fig, ax = plt.subplots()
fig = plt.figure()
ax  = fig.add_subplot(111,projection='3d')

x   = 1.
y   = 1.
z   = 1.

ax.scatter(x, y, z )

def animate(i):
    global x,y,z,n
    
    dt  = 1./64.
    s   = 10.
    b   = 8./3.
    r   = 28
    
    
    
    xdot    = s * (y-x)
    ydot    = x*r -x*z -y
    zdot    = x*y -b*z

    x       = x+xdot*dt
    y       = y+ydot*dt
    z       = z+zdot*dt


    ax.scatter(x,y,z, c=randColor(), marker=randMarker())

ani = animation.FuncAnimation(fig, animate, frames=2000, interval=10, blit=False)
#ani.save('LorenzAtrractor3D.mp4',bitrate=1024)
plt.show()

.