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 np

n = 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 = 0
h[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 plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from mpl_toolkits.mplot3d import Axes3D
a = n
x = 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 animation


def 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 np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

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

plt.show()


the result


the snapshot



Wednesday, December 17, 2014

Matplotlib Animation in Python

Here is the update from before

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def 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, t

def simPoints(simData):
x, t = simData[0], simData[1]
line.set_data(t, x)
return line
n = 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 np
from matplotlib import pyplot as plt
from matplotlib import animation

fig = plt.figure()
n = 10
x = np.linspace(0,2*np.pi,100)



def init():
pass
def 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;

interface

uses
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;

interface

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