Nugroho's blog.

Friday, December 19, 2014

3D Waterwave Simulation using Python

I used Numpy  Matplotlib with Animation and 3d Plot module on my OS X Yosemite.

The code is still messy and clearly not efficient (there's slow loop here and there) but it works, :)
Here The Result
The Code
`import numpy as npn = 8;g = 9.8;dt = 0.02;dx = 1.0;dy = 1.0;h = np.ones((n+2,n+2))u = np.zeros((n+2,n+2))v = np.zeros((n+2,n+2))hx = np.zeros((n+1,n+1))ux = np.zeros((n+1,n+1))vx = np.zeros((n+1,n+1))hy = np.zeros((n+1,n+1))uy = np.zeros((n+1,n+1))vy = np.zeros((n+1,n+1))nsteps = 0h[1,1] = .5;def reflective():    h[:,0] = h[:,1]    h[:,n+1] = h[:,n]    h[0,:] = h[1,:]    h[n+1,:] = h[n,:]    u[:,0] = u[:,1]    u[:,n+1] = u[:,n]    u[0,:] = -u[1,:]    u[n+1,:] = -u[n,:]    v[:,0] = -v[:,1]    v[:,n+1] = -v[:,n]    v[0,:] = v[1,:]    v[n+1,:] = v[n,:]def proses():    #hx = (h[1:,:]+h[:-1,:])/2-dt/(2*dx)*(u[1:,:]-u[:-1,:])    for i in range (n+1):        for j in range(n):            hx[i,j] = (h[i+1,j+1]+h[i,j+1])/2 - dt/(2*dx)*(u[i+1,j+1]-u[i,j+1])            ux[i,j] = (u[i+1,j+1]+u[i,j+1])/2- dt/(2*dx)*((pow(u[i+1,j+1],2)/h[i+1,j+1]+ g/2*pow(h[i+1,j+1],2))- (pow(u[i,j+1],2)/h[i,j+1]+ g/2*pow(h[i,j+1],2)))            vx[i,j] = (v[i+1,j+1]+v[i,j+1])/2 - dt/(2*dx)*((u[i+1,j+1]*v[i+1,j+1]/h[i+1,j+1]) - (u[i,j+1]*v[i,j+1]/h[i,j+1]))    for i in range (n):        for j in range(n+1):            hy[i,j] = (h[i+1,j+1]+h[i+1,j])/2 - dt/(2*dy)*(v[i+1,j+1]-v[i+1,j])            uy[i,j] = (u[i+1,j+1]+u[i+1,j])/2 - dt/(2*dy)*((v[i+1,j+1]*u[i+1,j+1]/h[i+1,j+1]) - (v[i+1,j]*u[i+1,j]/h[i+1,j]))            vy[i,j] = (v[i+1,j+1]+v[i+1,j])/2 - dt/(2*dy)*((pow(v[i+1,j+1],2)/h[i+1,j+1] + g/2*pow(h[i+1,j+1],2)) - (pow(v[i+1,j],2)/h[i+1,j] + g/2*pow(h[i+1,j],2)))        for i in range (1,n+1):        for j in range(1,n+1):            h[i,j] = h[i,j] - (dt/dx)*(ux[i,j-1]-ux[i-1,j-1]) - (dt/dy)*(vy[i-1,j]-vy[i-1,j-1])            u[i,j] = u[i,j] - (dt/dx)*((pow(ux[i,j-1],2)/hx[i,j-1] + g/2*pow(hx[i,j-1],2)) - (pow(ux[i-1,j-1],2)/hx[i-1,j-1] + g/2*pow(hx[i-1,j-1],2))) - (dt/dy)*((vy[i-1,j]*uy[i-1,j]/hy[i-1,j]) - (vy[i-1,j-1]*uy[i-1,j-1]/hy[i-1,j-1]))            v[i,j] = v[i,j] - (dt/dx)*((ux[i,j-1]*vx[i,j-1]/hx[i,j-1]) - (ux[i-1,j-1]*vx[i-1,j-1]/hx[i-1,j-1])) - (dt/dy)*((pow(vy[i-1,j],2)/hy[i-1,j] + g/2*pow(hy[i-1,j],2)) - (pow(vy[i-1,j-1],2)/hy[i-1,j-1] + g/2*pow(hy[i-1,j-1],2)))#dh = dt/dt*(ux[1:,:]-ux[:-1,:])+ dt/dy*(vy[:,1:]-vy[:,:-1])    reflective()    return h,u,v'''for i in range (17):    #print h    proses(1)'''import matplotlib.pyplot as pltfrom matplotlib import cmfrom matplotlib.ticker import LinearLocator, FormatStrFormatterfrom mpl_toolkits.mplot3d import Axes3Da = nx = np.arange(n+2)y = np.arange(n+2)x,y = np.meshgrid(x,y)fig = plt.figure()ax = fig.add_subplot(111, projection='3d')def plotset():    ax.set_xlim3d(0, a)    ax.set_ylim3d(0, a)    ax.set_zlim3d(0.5,1.5)    ax.set_autoscalez_on(False)    ax.zaxis.set_major_locator(LinearLocator(10))    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))    cset = ax.contour(x, y, h, zdir='x', offset=0 , cmap=cm.coolwarm)    cset = ax.contour(x, y, h, zdir='y', offset=n , cmap=cm.coolwarm)    cset = ax.contour(x, y, h, zdir='z', offset=.5, cmap=cm.coolwarm)plotset()surf = ax.plot_surface(x, y, h,rstride=1, cstride=1,cmap=cm.coolwarm,linewidth=0,antialiased=False, alpha=0.7)fig.colorbar(surf, shrink=0.5, aspect=5)from matplotlib import animationdef data(k,h,surf):    proses()    ax.clear()    plotset()    surf = ax.plot_surface(x, y, h,rstride=1, cstride=1,cmap=cm.coolwarm,linewidth=0,antialiased=False, alpha=0.7)    return surf,ani = animation.FuncAnimation(fig, data, fargs=(h,surf), interval=10, blit=False)#ani.save('air.mp4', bitrate=512)plt.show()`

and the snapshot

Thursday, December 18, 2014

3D Surface Plot Animation using Matplotlib in Python

And here's the animation
```import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

def data(i, z, line):
z = np.sin(x+y+i)
ax.clear()
line = ax.plot_surface(x, y, z,color= 'b')
return line,

n = 2.*np.pi
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

x = np.linspace(0,n,100)
y = np.linspace(0,n,100)
x,y = np.meshgrid(x,y)
z = np.sin(x+y)
line = ax.plot_surface(x, y, z,color= 'b')

ani = animation.FuncAnimation(fig, data, fargs=(z, line), interval=30, blit=False)

plt.show()```

The result

The snapshot

3D Surface Plot using Matplotlib in Python

It's slightly modified from before

`import numpy as npimport matplotlib.pyplot as pltimport matplotlib.animation as animationfrom mpl_toolkits.mplot3d import Axes3Dn = 2.*np.pifig = plt.figure()ax = fig.add_subplot(111, projection='3d')x = np.linspace(0,n,100)y = np.linspace(0,n,100)x,y = np.meshgrid(x,y)z = np.sin(x+y)line = ax.plot_surface(x, y, z,color= 'b')plt.show()`

the result

the snapshot

Wednesday, December 17, 2014

Matplotlib Animation in Python

Here is the update from before

`import numpy as npimport matplotlib.pyplot as pltimport matplotlib.animation as animationdef simData():    t_max = n    dt = 1./8    k = 0.0    t = np.linspace(0,t_max,100)    while k < t_max:        x = np.sin(np.pi*t+np.pi*k)        k = k + dt        yield x, tdef simPoints(simData):    x, t = simData[0], simData[1]    line.set_data(t, x)    return linen = 2.fig = plt.figure()ax = fig.add_subplot(111)line, = ax.plot([], [], 'b')ax.set_ylim(-1, 1)ax.set_xlim(0, n)ani = animation.FuncAnimation(fig, simPoints, simData, blit=False,\                              interval=100, repeat=True)plt.show()`

and the result

Tuesday, December 16, 2014

Playing with Matplotlib Animation in Python

Coding like this

`import numpy as npfrom matplotlib import pyplot as pltfrom matplotlib import animationfig = plt.figure()n = 10x = np.linspace(0,2*np.pi,100)def init():    passdef animate(k):    h = np.sin(x+np.pik)    plt.plot(x,h)ax = plt.axes(xlim=(0, 2*np.pi), ylim=(-1.1, 1.1))anim = animation.FuncAnimation(fig, animate,init_func=init,frames=360,interval=20,blit=False)plt.show()`

The result

Friday, December 5, 2014

Playing (again) with 'Home Made' Vector in Delphi

Here it is. I create a vector as new type, which is in itself is three dimension array.

Then I declared u as vector with three dimension;
u (h,i,j)

where h = 0, 1, 2  as physical component (eg: height, velocity, momentum)
i , j = 0, 1, 2, ..., n as row n column

So if we read u[0,1,1], it means height value at coordinate (1,1); u[1,1,1] is the velocity value; [2,1,1] is the momentum value at the same coordinate.

Trying some of properties of it. I found out that we can initialize all component of vector-u with this one line code

`u:=fu(h[i,j],i,j);`

so the component u(h,i,j) will filled. Notice that the function has vector (or in this case array) return value.

The code below show how I fill the value of component u(0, i, j)

`unit Unit1;interfaceuses  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,  Dialogs, StdCtrls;const n=3;type vector=array[0..2,0..n,0..n]of real;type  TForm1 = class(TForm)    Button1: TButton;    Memo1: TMemo;    function fu(a:real;i,j:integer):vector;    procedure FormCreate(Sender: TObject);  private    { Private declarations }  public    { Public declarations }  end;var  Form1: TForm1;implementation{\$R *.dfm}function tform1.fu(a:real;i,j:integer):vector;begin  fu[0,i,j]:=a;end;procedure TForm1.FormCreate(Sender: TObject);var i,j:integer;h:array[0..n,0..n]of real;u:vector;begin  for i:=0 to n do begin    for j:=0 to n do begin       h[i,j]:=1;       u:=fu(h[i,j],i,j);    end;  end;  memo1.Text:='';  memo1.Lines.Append('h[1,1]='+floattostr(h[1,1]));  memo1.Lines.Append('u[0,1,1]='+floattostr(u[0,1,1]));  memo1.Lines.Append('u[0,2,1]='+floattostr(u[0,2,1]));end;end.`
:)

Wednesday, November 26, 2014

Returning Function as Array in Delphi

Do you wonder how to do vector operation in Delphi? No, of course, :).

We could go like this.

`function tform1.adv(a,b:real):real;begin  adv:=a+b;end;`

The problem is the return is real, which is single value only. We want a and b as vector. Wait...

How we define vector in Delphi? I don't know. I used to treat a vector in Delphi as array. So I coded it like this

`var a,b:array[0..1]of real;`

So far I had no problem. Lately, I am going crazy with overuse functions in Delphi, and trying operating vectors using function too.

But if I write the code like this
`function tform1.adv(a,b:array[0..1]of real):real;begin  adv:=a[0]+b[0];    {a[1]+b[1]?}end;`

It will only return one value. So I improvised by modify it

like this

`function tform1.adv(a,b:array[0..1]of real):array[0..1]ofreal;begin  adv[0]:=a[0]+b[0];    adv[1]:=a[1]+b[1];  end;`

But it won't compile. (it will give error message "identifier expected but ARRAY found"). So I try another approach
`type  vector=array[0..1] of real;function tform1.adv(a,b:vector):vector;begin  adv[0]:=a[0]+b[0];  adv[1]:=a[1]+b[1];end;`

It works, :).

Here my last night tinkering with "vector" in Delphi, :)

`unit Unit1;interfaceuses  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,  Dialogs, Math;type  vector=array[0..1] of real;type  TForm1 = class(TForm)    procedure proses;    function mux(a:real;b:vector):vector;    function dot(a,b:vector):real;    function norm(a,b:vector):vector;    function adv(a,b:vector):vector;    function suv(a,b:vector):vector;    procedure FormCreate(Sender: TObject);  private    { Private declarations }  public    { Public declarations }  end;const  n=3;var  Form1: TForm1;  r,v:array[1..n] of vector;implementation{\$R *.dfm}function tform1.adv(a,b:vector):vector;begin  adv[0]:=a[0]+b[0];  adv[1]:=a[1]+b[1];end;function tform1.suv(a,b:vector):vector;begin  suv[0]:=a[0]-b[0];  suv[1]:=a[1]-b[1];end;function tform1.mux(a:real;b:vector):vector;begin  mux[0]:=a*b[0];  mux[1]:=a*b[1];end;function tform1.dot(a,b:vector):real;begin  dot:=a[0]*b[0]+a[1]*b[1];end;function tform1.norm(a,b:vector):vector;var mag,i,j:real;begin  i:=b[0]-a[0];  j:=b[1]-a[1];  mag:=sqrt(power(i,2)+power(j,2));  if mag<>0 then begin    norm[0]:=i/mag;    norm[1]:=j/mag;  end;end;procedure tform1.proses;var direction:vector;vi,vj,swap:real;i,j:integer;begin  j:=2;i:=1;  direction:=norm(r[j],r[i]);//call function  vi:=dot(v[i],direction);  vj:=dot(v[j],direction);  swap:=vj-vi;  v[i]:=adv(v[i],mux(swap,direction));  v[j]:=suv(v[j],mux(swap,direction));end;procedure TForm1.FormCreate(Sender: TObject);begin  proses;end;end.`
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)