Nugroho's blog.: programming

## Pages

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 ndef 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 ndef 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 Nx = []y = []#initial valuex.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 testingxc = 3yc = N(xc,x,y)print ''print xc, yc#plott = 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;interfaceuses  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 turtlefrom random import uniformimport numpy as npturtle.shape("turtle")#turtle.speed(1)x = 0y = 0rmax=40for 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 turtlefrom random import uniformimport numpy as npturtle.shape("turtle")#turtle.speed(1)x = 0y = 0for 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 turtlefrom random import uniformturtle.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.

`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
`#codefrom visual import *from random import uniformdisplay(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,                width=600, height=600,  forward=(-.4,-.3,-1)) #arah kameradistant_light(direction=(1,1,1), color=color.red)Ball    = sphere(radius=2, length=4, opacity=.3)Bola    = []n       = 5for 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./16def 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   = vjdef 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, :)

`#codefrom visual import *from random import uniformdisplay(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,                width=600, height=600,  forward=(-.4,-.3,-1)) #arah kameradistant_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      = 1bola.x      = -1bola.z      = 1v           = vector(1,-1,0)dt          = 1./16r           = bola.posrc          = Cone.radiush           = 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 vdef 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 uniformdisplay(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,                width=600, height=600,  forward=(-.4,-.3,-1)) #arah kameradistant_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      = 1bola.x      = -1bola.z      = 1v           = vector(2,1,0)dt          = 1./16r           = bola.posrc          = Ball.radiusdef 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 - vpdef 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 uniformdisplay(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,                width=600, height=600,  forward=(-.4,-.3,-1)) #arah kameradistant_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      = 0bola.x      = 0bola.z      = 1v           = vector(2,0,0)dt          = 1./16r           = bola.posrc          = silinder.radiusdef 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 - vpdef 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)
silinder.rotate(angle=pi/2, axis=(0,0,1),origin=(0,0,0))

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, :)

`#codefrom visual import *n       = 13display(center=(n/2,0,0),background=(1,1,1), autoscale=False, range=(7),                width=600, height=600,  forward=(-.4,-.3,-1)) #arah kameradistant_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)#usikanl[0]        = 1.2 l[n-1]      = .9#posisi xposition    = 0for 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 kameradistant_light(direction=(1,1,1), color=color.red)n       = 1dt      = 1./8.dl      = ones(1)pegas   = []print dlfor 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 kameradistant_light(direction=(1,1,1), color=color.red)n       = 1dt      = 1./8.#dl      = ones(1)dl      = []dl.append(.1)pegas   = []print dlfor 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

`#codefrom visual import *from random import uniform,randomfrom 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 kameradistant_light(direction=(1,1,1), color=color.red)n       = 2dt      = 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.1x0      = l0x1      = l0+l1v0      = 0.v1      = 0.y = 1.print x1def 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,randomfrom 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 kameradistant_light(direction=(1,1,1), color=color.red)n       = 2dt      = 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.1l1      = 1.x0      = l0+w[0]/2v0      = 0.v1      = 0.y = 1.print pegas[0].lengthdef 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) 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) 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)