Nugroho's blog.: programming
Showing posts with label programming. Show all posts
Showing posts with label programming. Show all posts

Friday, April 14, 2017

Numeric Integration using Python and Pylab

Pylab is Python module that contain NumPy and MatPlotlib.

Let f(x) = x^2 and we want to integrate it from 0 to 1.

Numerically we could use square method or trapezoidal (which is almost better).


Delphi.


Still using Delphi?

Yup, it teach me about programming discipline, :)


Thursday, April 13, 2017

Flappy Bird Like in Delphi

Remember the infamous Flappy Bird? Yup, I will create the program based on that algorithm.



Wednesday, April 12, 2017

Newton Polynomial in Python

I wrote code for this in Delphi. This time I want to rewrite it in Python based on this wiki.

I use this set of data point
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)

and I use xc=3 for the test data.

It's obvious that these sets of data points have quadratic form and f(xc) must have value of 9.




The heart of code lay on this




def n(j,xc,x):
n = 1
for i in arange(j):
n *= (xc-x[i])

return n
def a(j,l,x,y):
if j==0:
return y[0]
elif j-l==1 :
return (y[j]-y[l])/(x[j]-x[l])
else:
return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])


def N(xc,x,y):
N = 0
for j in range(len(x)):
N += a(j,0,x,y)*n(j,xc,x)

return N


Look at the function a(j,l,x,y), that's recursive function to obtain Newton divided difference value.

the whole code is this
from pylab import *

def n(j,xc,x):
n = 1
for i in arange(j):
n *= (xc-x[i])

return n
def a(j,l,x,y):
if j==0:
return y[0]
elif j-l==1 :
return (y[j]-y[l])/(x[j]-x[l])
else:
return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])


def N(xc,x,y):
N = 0
for j in range(len(x)):
N += a(j,0,x,y)*n(j,xc,x)

return N

x = []
y = []
#initial value
x.append(0)
x.append(1)
x.append(2)
x.append(4)
x.append(5)

y.append(0)
y.append(1)
y.append(4)
y.append(16)
y.append(25)

#for testing
xc = 3

yc = N(xc,x,y)

print ''
print xc, yc
#plot
t = linspace(-7,7,100)
u = N(t,x,y)
plot(t,u)
grid(True)
show()



.


here's the graphics


with another sets of data points, I have another result
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)




Tuesday, April 11, 2017

Newton Polinomial.


Here's code for Newton's divided differences interpolation polynomial (quite mouthful huh, :) ).

The purpose of this method is to create a function (polynomial) that passes through given set of data points.




I read data point from several edit box.

procedure TForm1.Button3Click(Sender: TObject);
var i:integer;
begin
for i:=0 to n do begin
x[i]:=strToFloat(kx[i].Text);
y[i]:=strToFloat(ky[i].Text);
end;
xc:=strToFloat(kxc.Text);
yc:=fn(xc);
kyc.Text:=floatToStr(yc);
gambarNewton;
end;

.
kx and ky is tEdit created when button1 is clicked

procedure TForm1.Button1Click(Sender: TObject);
var i:integer;
begin
button2.Enabled:=true;
button3.Enabled:=true;
button4.Enabled:=true;
button5.Enabled:=true;

n:=strToInt(edit1.Text);
kxc:=tEdit.Create(form1); kyc:=tEdit.Create(form1);
kxc.Parent:=form1; kyc.Parent:=form1;
kxc.Left:=36; kyc.Left:=72;
kxc.Width:=36; kyc.Width:=36;
kxc.Text:='0,5';

for i:=0 to n do begin
kx[i]:=tEdit.Create(form1); ky[i]:=tEdit.Create(form1);
kx[i].Parent:=form1; ky[i].Parent:=form1;
kx[i].Top:=36+36*i; ky[i].Top:=36+36*i;
kx[i].Left:=36; ky[i].Left:=72;
kx[i].Width:=36; ky[i].Width:=36;
kx[i].Text:=intToStr(i); ky[i].Text:=intToStr(i);
end;
end;


xc is x coordinate where the corresponding y (yc) is obtained using Newton method by calling it

yc=fn(xc)

function tform1.fn(xs:real):real;
var i:integer;fs:real;
begin
fs:=0;
for i:=0 to n do begin
fs:=fs+b(i,0)*c(xs,i);
end;
fn:=fs;
end;



the fn function call the two other function. The b function, a recursive contain divided difference like this

function tform1.b(i,j:integer):real;
begin
if i=0 then b:=y[0]
else if (i-j)=1 then
b:=(y[i]-y[j])/(x[i]-x[j])
else
b:=(b(i,j+1)-b(i-1,j))/(x[i]-x[j]);
end;

and c function, a recursive function (or you could rewrite it using simple for command)

function tform1.c(xs:real;i:integer):real;
begin
if i=0 then c:=1
else c:=(xs-x[i-1])*c(xs,i-1);
end;

and finally, draw the data and the function on image1

fprocedure tform1.gambarNewton;
var i,x0,y0:integer;px,py:real;
begin
x0:=image1.Width div 2; y0:=image1.Height div 2;
image1.Canvas.Brush.Color:=clLime;
image1.Canvas.Rectangle(0,0,image1.Width,image1.Height);
image1.Canvas.Brush.Color:=clWhite;
image1.Canvas.Pen.Color:=clBlack;
image1.Canvas.MoveTo(0,y0); image1.Canvas.LineTo(image1.Width,y0);
image1.Canvas.MoveTo(x0,0); image1.Canvas.LineTo(x0,image1.Height);
for i:=-300 to 300 do begin
px:=i/skala; py:=skala*fn(px);
image1.Canvas.Pixels[x0+i,y0-round(py)]:=clGreen;
end;
for i:=0 to n do begin
px:=x0+skala*x[i]; py:=y0-skala*y[i];
image1.Canvas.Ellipse(round(px)-7,round(py)-7,round(px)+7,round(py)+7);
end;
px:=x0+skala*xc; py:=y0-skala*yc;
image1.Canvas.Brush.Color:=clred;
image1.Canvas.Ellipse(round(px)-7,round(py)-7,round(px)+7,round(py)+7);
image1.Canvas.Brush.Color:=clwhite;
end;






Monday, April 10, 2017

Short Function.


Here's my implementation of function according to The Power of 10;

"Restrict functions to a single printed page."

As bonus, I didn't use global variable if possible. So if a function or procedure need a variable from others, it have to be passed using parameter on that function.

If we look at the code below, we know that it can be rewritten using a long single procedure or function. But according The Power of Ten, a function should be as short as possible so it could be printed in a single page.

So, instead one long multiple page function, I write/break it as several short-single-printed-page functions. :)





unit Unit1;

interface

uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls;

type
TForm1 = class(TForm)
Edit1: TEdit;
Edit2: TEdit;
Button1: TButton;
Edit3: TEdit;
Edit4: TEdit;
procedure proses;
function konversi(a:real;c,d:char):string;
function konversiC(a:real;d:char):string;
function konversiF(a:real;d:char):string;
function konversiR(a:real;d:char):string;
function konversiK(a:real;d:char):string;
procedure Button1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Button1Click(Sender: TObject);
begin
proses;
end;
procedure tform1.proses;
var
a:real;
s,b:string;
c,d:char;
begin
a:=strToFloat(edit1.Text);
s:=edit2.Text;
c:=s[1];
s:=edit4.Text;
d:=s[1];
b:=konversi(a,c,d);
edit3.Text:=b;
end;
function tform1.konversi(a:real;c,d:char):string;
begin
case c of
'C':konversi:=konversiC(a,d);
'F':konversi:=konversiF(a,d)
'R':konversi:=konversiR(a,d)
'K':konversi:=konversiK(a,d)
else konversi:='error';
end;
end;
function tform1.konversiC(a:real;d:char):string;
begin
case d of
'C':konversiC:=floatToStr(a);
'F':konversiC:=floatToStr(a*9/5+32);
'R':konversiC:=floatToStr(a*4/5);
'K':konversiC:=floatToStr(a+273);
else konversiC:='Error';
end;
end;
function tform1.konversiF(a:real;d:char):string;
begin
case d of
'C':konversiF:=floatToStr((a-32)*5/9);
'F':konversiF:=floatToStr(a);
'R':konversiF:=floatToStr((a-32)*4/9);
'K':konversiF:=floatToStr((a-32)*5/9+273);
else konversiF:='Error';
end;
end;
function tform1.konversiR(a:real;d:char):string;
begin
case d of
'C':konversiR:=floatToStr(a*5/4);
'F':konversiR:=floatToStr(a*9/4+32);
'R':konversiR:=floatToStr(a);
'K':konversiR:=floatToStr(a*5/4+273);
else konversiR:='Error';
end;
end;
function tform1.konversiK(a:real;d:char):string;
begin
case d of
'C':konversiK:=floatToStr(a-273);
'F':konversiK:=floatToStr((a-273)*9/5+32);
'R':konversiK:=floatToStr((a-273)*4/5);
'K':konversiK:=floatToStr(a);
else konversiK:='Error';
end;
end;

end.


.

Another Turtle in Circle

There's always another way to solve something.

So, I have another code for "Turtle in Circle" code, :)

In the script below, I use turtle position to determine if it's still inside circle or not. If it outside circle, instead of send it to zero position, I send it to random position inside circle.

import turtle
from random import uniform
import numpy as np

turtle.shape("turtle")
#turtle.speed(1)
x = 0
y = 0
rmax=40

for i in range (1,1000):
a = uniform (-90,90) #angle
turtle.left(a)
d = uniform (-75,75) #distance
x = turtle.xcor()+d*np.cos(a*np.pi/180)
y = turtle.ycor()+d*np.sin(a*np.pi/180)
r = np.sqrt(x*x+y*y)
if r>rmax:
turtle.setx(uniform(-rmax,rmax))
turtle.sety(uniform(-rmax,rmax))
x = 0
y = 0
else:
turtle.forward(d)

turtle.exitonclick()





.

Sunday, April 9, 2017

Turtle in Circle


I use previous code and improve it so the turtle could only move at certain circle area.

import turtle
from random import uniform
import numpy as np

turtle.shape("turtle")
#turtle.speed(1)
x = 0
y = 0

for i in range (1,1000):
a = uniform (-90,90) #angle
turtle.left(a)
d = uniform (-75,75) #distance
x += d*np.cos(np.pi*a/180)
y += d*np.sin(np.pi*a/180)
r = np.sqrt(x*x+y*y)
if r>40:
turtle.setx(0)
turtle.sety(0)
x = 0
y = 0
turtle.forward(d)

turtle.exitonclick()





.



Saturday, April 8, 2017

Random Turtle Movement.


I use turtle module, the standard module, in Python.

The turtle movement has random direction (angle), and random distance (forward/backward).

import turtle
from random import uniform

turtle.shape("turtle")
turtle.speed(1)
for i in range (1,100):
#random angle
a = uniform (-90,90)
turtle.left(a)
#random move
d = uniform (-100,100)
turtle.forward(d)

turtle.exitonclick()






.



Friday, April 7, 2017

Prime Number on Python

We've done it using Delphi, how about Python? Easy, :)


for i in range (2,200):
prime = True
for j in range (2,i):
if i%j==0:
prime=False
if prime==True:
print i


.

Prime Number on Delphi

Using Delphi to generate prime number? Okay, there's some code out there with

if (i mod 2 <>0) and (i mod 3 <>0) and (i mod 5 <>0) and (i mod 7 <>0) then i is prime.

Don't use that.

That's just for prime number below 100.

Use this instead.

var i,j:integer;
prime:boolean;
begin
for i:=2 to 200 do begin
prime:=true;
for j:=2 to i-1 do begin
if i mod j=0 then prime:=false;
end;
if prime=true then memo1.Lines.Append(intToStr(i));
end;
end;
.

Thursday, June 2, 2016

Collision.





 Here's the Code
#code
from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

Ball = sphere(radius=2, length=4, opacity=.3)

Bola = []
n = 5
for i in arange (n):
bola = sphere(color=color.green,radius=uniform(.2,.73))
bola.pos = vector(uniform(-1.5,1.5),uniform(-1.5,1.5),uniform(-1.5,1.5))
bola.v = vector(uniform(-1,1),uniform(-1,1),uniform(-1,1))
Bola.append(bola)

dt = 1./16

def pantul():
global Bola
for bola in Bola:
r = bola.pos
v = bola.v
if mag(r)>=Ball.radius:
r = 1.9*norm(r)
vp = (dot(v,norm(r)))*norm(r)
vr = v-vp
v = vr - vp
bola.r = r
bola.v = v
for i in arange (n-1):
for j in range(i+1, n):
ri = Bola[i].pos
rj = Bola[j].pos
vi = Bola[i].v
vj = Bola[j].v
rc = rj-ri
if Bola[i].radius+Bola[j].radius>mag(rc):
vpi = dot(vi,norm(rc))*norm(rc)
vri = vi-vpi
vpj = -dot(vj,norm(rc))*norm(rc)
vrj = vj-vpj

vi = vpj+vri
vj = vpi+vrj

Bola[i].v = vi
Bola[i].v = vj

def proses():
for bola in Bola:
r = bola.pos
v = bola.v
a = vector(0,0,0)
v += a*dt
r += v*dt

bola.pos = r

pantul()

while 1:
rate(37)
proses()

.

Sunday, May 29, 2016

Bouncing Ball inside a Cone


  I use vector projection and rejection to calculate velocity after bouncing the side of cone, :)





#code
from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

Cone = cone(pos = (0,0,0), axis=(0,5,0), radius = 3, opacity = .2)


bola = sphere(color=color.green,radius=.2)
bola.y = 1
bola.x = -1
bola.z = 1

v = vector(1,-1,0)
dt = 1./16
r = bola.pos
rc = Cone.radius
h = vector(Cone.axis)

def pantul():
global r,v
#tumbukan dengan lantai
if r.y<0:
r.y = 0
v.y *= -1

rp = vector(r.x,0,r.z)
hb = h.y - r.y
rmaks = hb/h.y*rc
c = h-rmaks*norm(rp) #vektor garis singgung
#selimut kerucut dengan bidang singgung
#tumbukan dengan selimut kerucut
if mag(rp)>rmaks:
rp = norm(rp)*rmaks
r = vector(rp.x,r.y,rp.z)
vp = dot(v,norm(c))*norm(c)
v = 2*vp-v
print v



def proses():
global r,v
a = vector(0,0,0)
v += a*dt
r += v*dt

bola.pos = r

pantul()

while 1:
rate(37)
proses()

.

Wednesday, May 25, 2016

Bouncing Ball inside Sphere


from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

Ball = sphere(radius=2, length=4, opacity=.3)


bola = sphere(color=color.green,radius=.2)
bola.y = 1
bola.x = -1
bola.z = 1

v = vector(2,1,0)
dt = 1./16
r = bola.pos
rc = Ball.radius

def pantul():
global r,v
if mag(r)>=rc:
r = 1.99*norm(r)
vp = (dot(v,norm(r)))*norm(r)
vr = v-vp
v = vr - vp

def proses():
global r,v
a = vector(0,0,0)
v += a*dt
r += v*dt

bola.pos = r

pantul()

while 1:
rate(37)
proses()


.


Tuesday, May 24, 2016

Bouncing inside Cylinder




from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

silinder = cylinder(radius=2, length=4, opacity=.3)
silinder.rotate(angle=pi/2, axis=(0,0,1),origin=(0,0,0))


bola = sphere(color=color.green,radius=.2)
bola.y = 0
bola.x = 0
bola.z = 1

v = vector(2,0,0)
dt = 1./16
r = bola.pos
rc = silinder.radius

def pantul():
global r,v
if mag(r)>=rc:
r = 1.9*norm(r)
vp = (dot(v,norm(r)))*norm(r)
vr = v-vp
v = vr - vp

def proses():
global r,v
a = vector(0,0,0)
v += a*dt
r += v*dt

bola.pos = r

pantul()

while 1:
rate(7)
proses()

.

Thursday, May 12, 2016

Bounce Over Spherical Surface










#code
from visual import *
from random import uniform


display(center=(0,2,0),background=(1,1,1), autoscale=False, range=7.5,        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

ball        = sphere(radius=2, color=color.red, opacity = .5)
r2          = ball.radius
silinder    = cylinder(radius=2, length=4, opacity=.3)
silinder.rotate(angle=pi/2, axis=(0,0,1),origin=(0,0,0))


bola        = sphere(color=color.green,radius=.2)
bola.y      = 3
bola.x      = uniform(-1,1)
bola.z      = uniform(-1,1)


v           = vector(0,2,0)
dt          = 1./8.
r           = bola.pos

def pantul():
    global r,v
    print v
    if mag(r)<r2:
        print mag(r)
        arah    = norm(r)
        dv      = dot(v,arah)
        v       -= dv*arah
        r       = (r2+.2)*arah
    
    

def proses():
    global r,v
    a   = vector(0,-1,0)
    v   += a*dt
    r   += v*dt

    bola.pos  = r

    pantul()

    

while 1:
    rate(11)
    proses()

.




Tuesday, May 10, 2016

N-Spring System

 Using Visual Python

 I like the result, :)

#code
from visual import *

n = 13
display(center=(n/2,0,0),background=(1,1,1), autoscale=False, range=(7),
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.orange)
dt = 1./8.
dx = []
k = []
m = []
l0 = []
l = []
v = []
x = []
for i in arange(n):
dx.append(.1)
k.append(1.)
m.append(1.)
l0.append(1.)
l.append(1.)
v.append(0.)
x.append(0.)
pegas = []
kotak = []

for i in arange(n):
spring = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
pegas.append(spring)
ko = box(pos=(0,0,0), width=.5, height=.5, length= .5, color= color.green)
kotak.append(ko)

box(pos=(-1,.64,0), width=n, height=2, length= 2, color=color.black)
box(pos=(n/2.,-.36,0), width=n, height=.2, length= n, color=color.white,opacity=.9)

#usikan
l[0] = 1.2
l[n-1] = .9
#posisi x
position = 0
for i in arange(n):
position +=l[i]
x[i] = position
kotak[i].x = x[i]


def updatePegas():
global l
for i in arange(n):
if i!=0:
pegas[i].x = x[i-1]
l[i] = x[i]-x[i-1]
else:
l[i] = x[i]

kotak[i].x = x[i]
pegas[i].length = l[i]

def proses():
for i in arange(n):
dx[i] = l[i]-l0[i]
f0 = -k[i]*dx[i]
if i<n-1:
dx[i+1] = l[i+1]-l0[i+1]
f1 = -k[i+1]*dx[i+1]
a = (f0-f1)/m[i]
else:
a = f0/m[i]
v[i] += a*dt
x[i] += v[i]*dt
updatePegas()
while 1:
rate (19)
proses()

.

Monday, May 9, 2016

Here's the Culprit



 In Visual Python, helix object will generate error if helix.length = some array like the code below. I use dl, an array, for the length value
from visual import *

display(center=(1,0,0),background=(1,1,1), autoscale=False, range=(2,2,2),
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n = 1
dt = 1./8.
dl = ones(1)

pegas = []

print dl


for i in arange(n):
spring = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
pegas.append(spring)

def updatePegas(l):
pegas[0].length = l



def proses():
global dl
l = pegas[0].length
if l>2:
l = 2
dl[0] *= -1
elif l<.5:
l = .5
dl[0] *= -1
l += dl[0]
updatePegas(l)

while 1:
rate (19)
proses()








.



And the result is
Traceback (most recent call last):
File "springList.py", line 40, in
proses()
File "springList.py", line 36, in proses
updatePegas(l)
File "springList.py", line 22, in updatePegas
pegas[0].length = l
File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/VPython-6.11-py2.7-macosx-10.6-intel.egg/visual_common/primitives.py", line 850, in set_length
self.__frame.axis = self.__axis.norm()
AttributeError: 'numpy.ndarray' object has no attribute 'norm'


If we change dl from array to list, like the code below, everything is suddenly OK, :)
from visual import *

display(center=(1,0,0),background=(1,1,1), autoscale=False, range=(2,2,2),
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n = 1
dt = 1./8.
#dl = ones(1)
dl = []
dl.append(.1)
pegas = []

print dl


for i in arange(n):
spring = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
pegas.append(spring)

def updatePegas(l):
pegas[0].length = l



def proses():
global dl
l = pegas[0].length
if l>2:
l = 2
dl[0] *= -1
elif l<.5:
l = .5
dl[0] *= -1
l += dl[0]
updatePegas(l)

while 1:
rate (19)
proses()








.

Gonna rewrite the code.


Sunday, May 8, 2016

Double Spring System

Akhirnya.

With Visual Python module

I couldn't use list for spring length since it'll trigger some error for helix object. It's very unfortunate because it would come handy as we expand the number of spring and mass.

Anyway, here's the code

#code
from visual import *
from random import uniform,random
from visual.controls import *

display(center=(1,0,0),background=(1,1,1), autoscale=False, range=(2,2,2),
width=600, height=600, forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n = 2
dt = 1./8.
pegas = []
kotak = []
box(pos=(-1,0,0), width=2, height=2, length= 2, color=color.black)
box(pos=(0,-.36,0), width=2, height=.2, length= 5, color=color.black,opacity=.3)
for i in arange(n):
spring = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
pegas.append(spring)
ko = box(pos=(0,0,0), width=.5, height=.5, length= .5, color= color.green)
kotak.append(ko)

k0 = 1.
k1 = 1.

m0 = 1.
m1 = 1.

l00 = 1.
l01 = 1.

l0 = 1.
l1 = 1.1

x0 = l0
x1 = l0+l1

v0 = 0.
v1 = 0.


y = 1.
print x1
def updatePegas():
global x0,x1
kotak[0].x = x0
pegas[0].length = l0

kotak[1].x = x1
pegas[1].x = l0
pegas[1].length = l1



def proses():
global l0,v0,x0,l1,v1,x1
#untuk m0
dx0 = l0-l00
f0 = -k0*dx0
dx1 = l1-l01
f1 = -k1*dx1
a0 = (f0-f1)/m0
v0 += a0*dt
x0 += v0*dt
l0 = x0
#untuk m1
a1 = f1/m1
v1 += a1*dt
x1 += v1*dt
l1 = x1-x0


updatePegas()

while 1:
rate (39)
y += .1
proses()








.


Saturday, May 7, 2016

Normal Mode


 Planned to write coupled oscillator, using double spring.

 For some reason, Vpython refused to set length of helix with some value from array or list or any value derived from them. It only accept plain number or  number from simple variable (like a = 3. ). It hugely messed the whole script as I used l = [] for helix length.

 So, rewriting the code, create l0 and l1 manually.

 Didn't have energy to code the rest. So, at the moment, just call it "normal mode" coupled oscillator, heheh...

#code

from visual import *
from random import uniform,random
from visual.controls import *

display(center=(0,0,0),background=(1,1,1), autoscale=False, range=(5,5,3),
width=600, height=600, forward=(-1.4,-1.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n = 2
dt = 1./8.
pegas = []
kotak = []
w = ones(n)
w /= 2.
for i in arange(n):
spring = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
pegas.append(spring)
ko = box(pos=(0,0,0), width=.5, height=.5, length= .5, color= color.green)
kotak.append(ko)

k0 = 1.
m0 = 1.

l00 = 1.
l01 = 1.

l0 = 1.1
l1 = 1.

x0 = l0+w[0]/2

v0 = 0.
v1 = 0.


y = 1.
print pegas[0].length
def updatePegas():
global x0
kotak[0].x = x0
pegas[0].length = l0

x1 = l0+w[0]+l1+w[1]/2
kotak[1].x = x1
pegas[1].x = l0+w[0]



def proses():
global l0,v0,x0
dx = l0-l00
f0 = -k0*dx
a0 = f0/m0
v0 += a0*dt
l0 += v0*dt
x0 = l0+w[0]/2


updatePegas()

while 1:
rate (19)
y += .1
proses()








.



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)