From: John Tsiombikas <nuclear@member.fsf.org>
Date: Fri, 20 Nov 2009 01:20:46 +0000 (+0000)
Subject: added summerhack
X-Git-Url: http://git.mutantstargoat.com/user/nuclear/?a=commitdiff_plain;h=a80992973feb47adfc1214ad1a3f7fc34eb345c0;p=summerhack

added summerhack
---

a80992973feb47adfc1214ad1a3f7fc34eb345c0
diff --git a/3dengfx.conf b/3dengfx.conf
new file mode 100644
index 0000000..0af8e71
--- /dev/null
+++ b/3dengfx.conf
@@ -0,0 +1,11 @@
+; The Lab Demos :: 3D engine configuration file
+; ---------------------------------------
+; this file is generated by config_tool (our GUI demo configuration
+; tool) so your best bet would be using that tool to configure the demo.
+; If that does not work for some reason, modify this file with care...
+
+fullscreen = true
+resolution = 1024x768
+bpp = dontcare
+zbuffer = dontcare
+stencil = dontcare
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..4292695
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,41 @@
+src := src/sumhack.cpp src/events.cpp
+bin := sumhack-thelab
+
+include src/parts/Makefile-part
+include src/sdlvf/Makefile-part
+
+obj := $(src:.cpp=.o) $(csrc:.c=.o)
+
+opt := -O1
+opt := -g
+CXXFLAGS := -ansi -pedantic -Wall $(opt) -Isrc/3dengfx/src `src/3dengfx/3dengfx-config --cflags`
+CFLAGS := -std=c89 -pedantic -Wall $(opt) `src/3dengfx/3dengfx-config --cflags`
+libs := src/3dengfx/lib3dengfx.a `src/3dengfx/3dengfx-config --libs-no-3dengfx` -lGL -lvorbisfile
+
+$(bin): $(obj) src/3dengfx/lib3dengfx.a data/tex_list
+	$(CXX) -o $@ $(obj) $(libs)
+
+src/3dengfx/lib3dengfx.a:
+	cd src/3dengfx;\
+	./configure --with-gfxlib=sdl --disable-ft --enable-opt --disable-debug && make
+
+data/tex_list: $(src)
+	tools/find_textures >$@
+
+include $(obj:.o=.d)
+
+%.d: %.cpp
+	@set -e; rm -f $@; $(CXX) -MM $(CXXFLAGS) $< > $@.$$$$; \
+	sed 's,\($*\)\.o[ :]*,\1.o $@ : ,g' < $@.$$$$ > $@; rm -f $@.$$$$
+
+%.d: %.c
+	@set -e; rm -f $@; $(CC) -MM $(CFLAGS) $< > $@.$$$$; \
+	sed 's,\($*\)\.o[ :]*,\1.o $@ : ,g' < $@.$$$$ > $@; rm -f $@.$$$$
+
+.PHONY: clean
+clean:
+	$(RM) $(obj) $(bin)
+
+.PHONY: cleandep
+cleandep:
+	find src \( -name '*.d' -o -name '*.d.*' \) -exec $(RM) '{}' ';'
diff --git a/README b/README
new file mode 100644
index 0000000..e37c14c
--- /dev/null
+++ b/README
@@ -0,0 +1,16 @@
+Summer Hack / The Lab Demos
+---------------------------
+Demo code: Nuclear
+3d engine: Nuclear, Samurai
+Sound sys: Savvy
+    Music: Tourmpadakis
+
+Greets to everyone who watches this small demoish thingy.
+
+About this thing...
+Needs a video card that doesn't suck for fillrate. No shaders or other
+high-end gimmicks are used this time.
+Released under the GPL, do whatever you wish with it.
+
+Contact info...
+http://thelab.demoscene.gr - nuclear@siggraph.org
diff --git a/data/curve00 b/data/curve00
new file mode 100644
index 0000000..b8a1f68
--- /dev/null
+++ b/data/curve00
@@ -0,0 +1,23 @@
+curve_3dengfx
+
+bspline
+19
+-44.089134 38.695248 1.000000
+-51.828182 25.327801 1.000000
+-54.833332 12.833332 1.000000
+-56.833332 -2.833332 1.000000
+-57.222065 -20.168427 1.000000
+-53.000763 -19.464882 1.000000
+-41.333332 -4.166668 1.000000
+-28.500000 0.166668 1.000000
+-25.000000 -1.833332 1.000000
+-32.166668 -16.833336 1.000000
+-31.166664 -22.000000 1.000000
+-8.833332 -4.333336 1.000000
+-10.833336 -0.666664 1.000000
+-17.999996 2.000000 1.000000
+-24.333332 -8.666664 1.000000
+-21.333332 -17.500000 1.000000
+-11.333336 -17.833336 1.000000
+-5.333332 -7.833332 1.000000
+10.500000 4.166664 1.000000
diff --git a/data/curve01 b/data/curve01
new file mode 100644
index 0000000..3d9e50a
--- /dev/null
+++ b/data/curve01
@@ -0,0 +1,11 @@
+curve_3dengfx
+
+bspline
+7
+-40.249992 -36.575001 1.000000
+-41.124996 -28.699993 1.000000
+-32.374996 24.850000 1.000000
+-30.274998 22.749998 1.000000
+-34.124996 2.974998 1.000000
+-38.499996 -11.200001 1.000000
+-45.499996 -19.774998 1.000000
diff --git a/data/curve02 b/data/curve02
new file mode 100644
index 0000000..de7f9fc
--- /dev/null
+++ b/data/curve02
@@ -0,0 +1,13 @@
+curve_3dengfx
+
+bspline
+9
+-17.500000 25.025000 1.000000
+-19.464882 28.141998 1.000000
+-53.938831 28.141998 1.000000
+-75.514366 5.862915 1.000000
+-74.810814 -5.159370 1.000000
+-60.036263 8.677116 1.000000
+-61.912399 -3.986778 1.000000
+-71.996613 -16.650681 1.000000
+-77.155983 -42.212997 1.000000
diff --git a/data/curve03 b/data/curve03
new file mode 100644
index 0000000..95a936b
--- /dev/null
+++ b/data/curve03
@@ -0,0 +1,44 @@
+curve_3dengfx
+
+bspline
+40
+16.666290 24.372211 1.000000
+13.978188 31.719719 1.000000
+11.290070 33.691002 1.000000
+8.960373 32.615753 1.000000
+3.763359 24.730627 1.000000
+1.254456 6.809883 1.000000
+5.734642 -1.075245 1.000000
+5.913841 -11.469280 1.000000
+1.971283 -17.562336 1.000000
+0.896042 -20.967274 1.000000
+4.480186 -18.816784 1.000000
+17.194740 -12.341385 1.000000
+26.485451 -17.749409 1.000000
+31.477478 -18.304085 1.000000
+40.510368 -14.308704 1.000000
+45.341881 -5.017342 1.000000
+41.439507 0.557480 1.000000
+34.365734 0.575157 1.000000
+29.764469 -3.163376 1.000000
+28.757938 -15.960651 1.000000
+38.104267 -16.823391 1.000000
+46.156490 -6.470535 1.000000
+44.862377 5.464008 1.000000
+44.431011 0.718948 1.000000
+42.995193 -11.710075 1.000000
+46.490738 -17.652496 1.000000
+48.937622 -17.128166 1.000000
+54.355721 -8.738857 1.000000
+58.026039 25.342695 1.000000
+56.802597 27.789574 1.000000
+55.579155 26.216579 1.000000
+53.831383 -3.495548 1.000000
+54.597290 -14.425369 1.000000
+59.892685 -15.520966 1.000000
+64.640274 -12.051575 1.000000
+66.831467 -2.921589 1.000000
+61.353485 1.278194 1.000000
+55.510284 -2.191193 1.000000
+51.675697 -18.990356 1.000000
+52.953896 -23.920540 1.000000
diff --git a/data/exibit.psys b/data/exibit.psys
new file mode 100644
index 0000000..dc1c7c6
--- /dev/null
+++ b/data/exibit.psys
@@ -0,0 +1,42 @@
+# --- particle system options ---
+# psize: <scalar>  (particle size)
+# psize-r: <scalar> (particle size range, unused currently)
+# life: <scalar> (particle life in seconds)
+# life-r: <scalar> (particle life range)
+# birth-rate: <scalar> (birth rate, particles/sec)
+# birth-rate-r: <scalar> (birth rate range)
+# grav: <vecror3> (gravity vector)
+# shoot: <vector3> (shoot vector, direction AND force)
+# shoot-r: <vector3> (shoot vector range)
+# friction: <scalar> (environment friction multiplier)
+# spawn_off: <vector3> (spawn particles in specified offset from the emitter position)
+# spawn_off-r: <vector3> (range of spawn-offset random vector)
+# tex: <string> (texture filename)
+# color_start: <vector4> (color at time of birth, alpha included)
+# color_end: <vector4> (color at time of death, alpha included)
+# color: <vector4> (specify the same color for start/end)
+# rot: <scalar> (particle rotation in rad/sec)
+# halo: <string> (halo texture filename)
+# halo_size: <scalar> (halo size)
+# halo_rot: <scalar> (halo rotation in rad/sec)
+# big_particles: <boolean> (requirement for big screen-space particles)
+
+psize: 700.0
+life: 3.5
+birth-rate: 50.0
+grav: 0, -0.05, 0
+shoot: 0, 0.0, 0
+shoot-r: 0.5, 0.5, 0.5;
+friction: 0.995
+spawn_off: 0, 150.0, 0
+spawn_off-r: 5, 5, 5
+tex: data/img/psys03.png
+color_start: 1.0, 0.61, 0.18, 1.0
+color_end: 1.0, 0.91, 0.18, 1.0
+glob_rot: 0.0
+rot: 0.0
+
+halo: data/img/psys03.png
+halo_size: 3000.0
+
+big_particles: true
diff --git a/data/galaxy.psys b/data/galaxy.psys
new file mode 100644
index 0000000..b7b02b6
--- /dev/null
+++ b/data/galaxy.psys
@@ -0,0 +1,42 @@
+# --- particle system options ---
+# psize: <scalar>  (particle size)
+# psize-r: <scalar> (particle size range, unused currently)
+# life: <scalar> (particle life in seconds)
+# life-r: <scalar> (particle life range)
+# birth-rate: <scalar> (birth rate, particles/sec)
+# birth-rate-r: <scalar> (birth rate range)
+# grav: <vecror3> (gravity vector)
+# shoot: <vector3> (shoot vector, direction AND force)
+# shoot-r: <vector3> (shoot vector range)
+# friction: <scalar> (environment friction multiplier)
+# spawn_off: <vector3> (spawn particles in specified offset from the emitter position)
+# spawn_off-r: <vector3> (range of spawn-offset random vector)
+# tex: <string> (texture filename)
+# color_start: <vector4> (color at time of birth, alpha included)
+# color_end: <vector4> (color at time of death, alpha included)
+# color: <vector4> (specify the same color for start/end)
+# rot: <scalar> (particle rotation in rad/sec)
+# halo: <string> (halo texture filename)
+# halo_size: <scalar> (halo size)
+# halo_rot: <scalar> (halo rotation in rad/sec)
+# big_particles: <boolean> (requirement for big screen-space particles)
+
+psize: 800.0
+life: 5.0
+birth-rate: 130.0
+grav: 0, 0, 0
+shoot: 0, 0, 0
+shoot-r: 3.0, 0.2, 0.1;
+friction: 0.995
+spawn_off: 0, 0, 0
+spawn_off-r: 10, 1.0, 3
+tex: data/img/psys03.png
+color_start: 0.25, 0.4, 1.0, 1.0
+color_end: 0.5, 0.8, 1.0, 0.1
+glob_rot: 0.0
+rot: 0.0
+
+halo: data/img/psys03.png
+halo_size: 2800.0
+
+big_particles: true
diff --git a/data/pillar.curve b/data/pillar.curve
new file mode 100644
index 0000000..5d0e3e9
--- /dev/null
+++ b/data/pillar.curve
@@ -0,0 +1,15 @@
+curve_3dengfx
+
+catmullrom
+11
+4.192291 0.684456 0.000000
+5.047863 3.935621 0.000000
+5.304535 11.122410 0.000000
+4.192291 20.704794 0.000000
+4.962307 22.073706 0.000000
+5.561207 22.672604 0.000000
+5.047863 23.271503 0.000000
+5.390091 23.442617 0.000000
+5.903435 23.870401 0.000000
+5.218979 24.298187 0.000000
+0.085560 24.300000 0.000000
diff --git a/data/stars.psys b/data/stars.psys
new file mode 100644
index 0000000..7f6cabd
--- /dev/null
+++ b/data/stars.psys
@@ -0,0 +1,38 @@
+# --- particle system options ---
+# psize: <scalar>  (particle size)
+# psize-r: <scalar> (particle size range, unused currently)
+# life: <scalar> (particle life in seconds)
+# life-r: <scalar> (particle life range)
+# birth-rate: <scalar> (birth rate, particles/sec)
+# birth-rate-r: <scalar> (birth rate range)
+# grav: <vecror3> (gravity vector)
+# shoot: <vector3> (shoot vector, direction AND force)
+# shoot-r: <vector3> (shoot vector range)
+# friction: <scalar> (environment friction multiplier)
+# spawn_off: <vector3> (spawn particles in specified offset from the emitter position)
+# spawn_off-r: <vector3> (range of spawn-offset random vector)
+# tex: <string> (texture filename)
+# color_start: <vector4> (color at time of birth, alpha included)
+# color_end: <vector4> (color at time of death, alpha included)
+# color: <vector4> (specify the same color for start/end)
+# rot: <scalar> (particle rotation in rad/sec)
+# halo: <string> (halo texture filename)
+# halo_size: <scalar> (halo size)
+# halo_rot: <scalar> (halo rotation in rad/sec)
+# big_particles: <boolean> (requirement for big screen-space particles)
+
+psize: 1000.0
+psize-r: 300.0
+life: 6.0
+birth-rate: 10.0
+grav: 0, 0, 0
+shoot: 0, 0, 0
+shoot-r: 0.0, 0.0, 0.0;
+friction: 1.0
+spawn_off: 0, 0, 0
+spawn_off-r: 500, 500, 500
+tex: data/img/psys03.png
+color_start: 0.6, 0.75, 1.0, 1.0
+color_end: 0.8, 0.75, 1.0, 0.0
+glob_rot: 0.0
+rot: 0.0
diff --git a/data/test2.psys b/data/test2.psys
new file mode 100644
index 0000000..dace8bb
--- /dev/null
+++ b/data/test2.psys
@@ -0,0 +1,36 @@
+# --- particle system options ---
+# psize: <scalar>  (particle size)
+# psize-r: <scalar> (particle size range, unused currently)
+# life: <scalar> (particle life in seconds)
+# life-r: <scalar> (particle life range)
+# birth-rate: <scalar> (birth rate, particles/sec)
+# birth-rate-r: <scalar> (birth rate range)
+# grav: <vecror3> (gravity vector)
+# shoot: <vector3> (shoot vector, direction AND force)
+# shoot-r: <vector3> (shoot vector range)
+# friction: <scalar> (environment friction multiplier)
+# spawn_off: <vector3> (spawn particles in specified offset from the emitter position)
+# spawn_off-r: <vector3> (range of spawn-offset random vector)
+# tex: <string> (texture filename)
+# color_start: <vector4> (color at time of birth, alpha included)
+# color_end: <vector4> (color at time of death, alpha included)
+# color: <vector4> (specify the same color for start/end)
+# rot: <scalar> (particle rotation in rad/sec)
+# halo: <string> (halo texture filename)
+# halo_size: <scalar> (halo size)
+# halo_rot: <scalar> (halo rotation in rad/sec)
+# big_particles: <boolean> (requirement for big screen-space particles)
+
+psize: 1000.0
+psize_end: 0.1
+life: 3.5
+birth-rate: 100.0
+grav: 0, 0, 0
+friction: 0
+tex: data/img/scl_blobs0000.png
+color_start: 0.25, 0.4, 1.0, 0.9
+color_end: 0.5, 0.8, 1.0, 0.2
+rot: 4.0
+glob_rot: 2.0
+
+big_particles: true
diff --git a/data/tex_list b/data/tex_list
new file mode 100644
index 0000000..2e4eb86
--- /dev/null
+++ b/data/tex_list
@@ -0,0 +1,25 @@
+data/img/benedeti.jpg
+data/img/brick1_base1.jpg
+data/img/cracks.jpg
+data/img/greets.png
+data/img/label_dolphin.png
+data/img/label_duck.png
+data/img/label_face.png
+data/img/label_teapot.png
+data/img/loading.jpg
+data/img/mud.jpg
+data/img/nebula_nx.jpg
+data/img/nebula_ny.jpg
+data/img/nebula_nz.jpg
+data/img/nebula_px.jpg
+data/img/nebula_py.jpg
+data/img/nebula_pz.jpg
+data/img/nx.jpg
+data/img/ny.jpg
+data/img/nz.jpg
+data/img/presents.png
+data/img/px.jpg
+data/img/py.jpg
+data/img/pz.jpg
+data/img/title.png
+data/img/volcano.jpg
diff --git a/data/thelab.psys b/data/thelab.psys
new file mode 100644
index 0000000..b6bb008
--- /dev/null
+++ b/data/thelab.psys
@@ -0,0 +1,47 @@
+# --- particle system options ---
+# psize: <scalar>  (particle size)
+# psize-r: <scalar> (particle size range, unused currently)
+# life: <scalar> (particle life in seconds)
+# life-r: <scalar> (particle life range)
+# birth-rate: <scalar> (birth rate, particles/sec)
+# birth-rate-r: <scalar> (birth rate range)
+# grav: <vecror3> (gravity vector)
+# shoot: <vector3> (shoot vector, direction AND force)
+# shoot-r: <vector3> (shoot vector range)
+# friction: <scalar> (environment friction multiplier)
+# spawn_off: <vector3> (spawn particles in specified offset from the emitter position)
+# spawn_off-r: <vector3> (range of spawn-offset random vector)
+# spawn_offset_curve: <string> (filename of a curve file)
+# tex: <string> (texture filename)
+# color_start: <vector4> (color at time of birth, alpha included)
+# color_end: <vector4> (color at time of death, alpha included)
+# color: <vector4> (specify the same color for start/end)
+# rot: <scalar> (particle rotation in rad/sec)
+# halo: <string> (halo texture filename)
+# halo_size: <scalar> (halo size)
+# halo_rot: <scalar> (halo rotation in rad/sec)
+# big_particles: <boolean> (requirement for big screen-space particles)
+
+psize: 450.0
+psize-r: 200.0
+life: 3.5
+birth-rate: 250.0
+#grav: 0, 0.0025, 0
+grav: 0, 0.0025, 0
+shoot: 0, 0, 0
+shoot-r: 0.1, 0.1, 0.1;
+friction: 0.995
+spawn_off: 0, 0, 0
+spawn_off-r: 0.7, 0.7, 0.7
+tex: data/img/psys03.png
+#tex: data/waterfall010002.png
+#tex: /home/nuclear/3ddata/textures/particles/sq_star.png
+
+#color_start: 0.25, 0.4, 1.0, 1.0
+#color_end: 0.5, 0.8, 1.0, 0.1
+#color_start: 1.0, 0.2, 0.3, 1.0
+#color_end: 0.3, 0.2, 1.0, 0.05
+color_start:	1, 0.78, 0.2, 1
+color_end:		1, 0.27, 0.2, 0.15
+
+big_particles: true
diff --git a/demoscript b/demoscript
new file mode 100644
index 0000000..776b728
--- /dev/null
+++ b/demoscript
@@ -0,0 +1,83 @@
+# 3dengfx demosystem script v1.2
+# ------------------------------
+# format: <time> <command> <part_name> [args]
+# If time has an 's' suffix it's seconds, otherwise milliseconds.
+# For a list of available commands, see the bottom of this file.
+
+# end of titles-start of galaxy flash
+0s		fx flash 36200 400
+
+# end of galaxy fade-out
+0s fx fade 56s 2s <0,0,0,0> 0 <0,0,0,1> 0
+
+# start of temple fade-in
+0s fx fade 58s 2s <0,0,0,1> 0 <0,0,0,0> 0
+
+#end of temple fade-out
+0s fx fade 87s 2s <0,0,0,0> 0 <0,0,0,1> 0
+
+# -- demo sequence --
+0s		rename_part distort zdistort
+0s		start_part	pipes
+36200	end_part	pipes
+36200	start_part	space
+58s		end_part	space
+58s		start_part	temple
+
+65s		start_part	zdistort
+65s		set_rtarget	temple	t0
+71283	set_rtarget	temple	fb
+71283	end_part	zdistort
+
+89s		end_part	temple
+89200	end
+#--------------------
+
+# -----------------
+# command reference
+# -----------------
+#
+# start_part <part>
+#		starts the specified demo-part.
+#
+# end_part <part>
+#		ends the specified demo-part.
+#
+# rename_part <part> <new name>
+#		renames a part.
+#		Note: overlapping parts are executed in lexicographical order.
+#
+# set_clear <part> true/false
+#		sets whether the specified part should clear the framebuffer
+#		before starting its graphics update.
+#
+# set_rtarget <part> <target>
+#		changes the render target of the specified part.
+#		<target> must be either fb (for framebuffer) or the name of a
+#		texture	register (t0, t1, t2, t3).
+#
+# fx <effect> <time> <duration> <effect specific arguments> ...
+#		schedules the specified post-processing effect for execution at
+#		the specified time, for the specified duration.
+#		Note: effects may overlapped, in which case are executed in the order
+#		that they are encountered in the demoscript.
+#
+# end
+#		ends the demo...
+#
+# ---------------------------------
+# Post-Processing Effects Reference
+# ---------------------------------
+#
+# neg <start> <duration>
+#		inverse video for the specified time interval.
+#
+# flash <peak-time> <duration> [color]
+#		screen flash, color is optional, default is white.
+#		Note: the flash is sinusoidal, with the peak occuring at the specified time.
+#
+# fade <start> <duration> [color1] [tex1] [color2] [tex2]
+#       crossfade from color/tex 1 to color/tex 2
+#
+# overlay <start> <duration> <texture> [shader]
+#		--- SUBJECT TO CHANGE ---
diff --git a/src/3dengfx b/src/3dengfx
new file mode 120000
index 0000000..86b5175
--- /dev/null
+++ b/src/3dengfx
@@ -0,0 +1 @@
+/home/nuclear/code/graphics/3dengfx_dev/3dengfx
\ No newline at end of file
diff --git a/src/events.cpp b/src/events.cpp
new file mode 100644
index 0000000..e67a736
--- /dev/null
+++ b/src/events.cpp
@@ -0,0 +1,76 @@
+#include "3dengfx/3dengfx.hpp"
+
+using namespace engfx_state;
+
+static float zoom_factor = 1.0;
+
+void key_handler(int key) {
+	TargetCamera *cam = dynamic_cast<TargetCamera*>(const_cast<Camera*>(view_mat_camera));
+	
+	switch(key) {
+	case 'a':
+		cam->zoom(0.8);
+		zoom_factor *= 0.8;
+		break;
+
+	case 'z':
+		cam->zoom(1.2);
+		zoom_factor *= 1.2;
+		break;
+
+	case '`':
+		screen_capture();
+		break;
+
+	case 27:
+		exit(0);
+
+	default:
+		break;
+	}
+}
+
+static int prev_x = -1;
+static int prev_y = -1;
+
+void motion_handler(int x, int y) {
+	if(prev_x != -1) {
+		TargetCamera *cam = dynamic_cast<TargetCamera*>(const_cast<Camera*>(view_mat_camera));
+		float dx = (float)(x - prev_x) * 0.01;
+		float dy = (float)(y - prev_y) * 1.0;
+		
+		cam->rotate(Vector3(0, -dx, 0));
+		cam->translate(Vector3(0, -dy * zoom_factor, 0));
+		
+		prev_x = x;
+		prev_y = y;
+	}
+}
+
+void bn_handler(int bn, int pressed, int x, int y) {
+	TargetCamera *cam = dynamic_cast<TargetCamera*>(const_cast<Camera*>(view_mat_camera));
+
+	switch(bn) {
+	case fxwt::BN_RIGHT:
+		if(pressed) {
+			prev_x = x;
+			prev_y = y;
+		} else {
+			prev_x = prev_y = -1;
+		}
+		break;
+
+	case fxwt::BN_WHEELUP:
+		cam->zoom(0.9);
+		zoom_factor *= 0.9;
+		break;
+
+	case fxwt::BN_WHEELDOWN:
+		cam->zoom(1.1);
+		zoom_factor *= 1.1;
+		break;
+
+	default:
+		break;
+	}
+}
diff --git a/src/events.hpp b/src/events.hpp
new file mode 100644
index 0000000..d21a09f
--- /dev/null
+++ b/src/events.hpp
@@ -0,0 +1,8 @@
+#ifndef _EVENTS_HPP_
+#define _EVENTS_HPP_
+
+void key_handler(int key);
+void motion_handler(int x, int y);
+void bn_handler(int bn, int pressed, int x, int y);
+
+#endif	// _EVENTS_HPP_
diff --git a/src/parts/.dist.hpp.swp b/src/parts/.dist.hpp.swp
new file mode 100644
index 0000000..8b934a2
Binary files /dev/null and b/src/parts/.dist.hpp.swp differ
diff --git a/src/parts/.space.cpp.swp b/src/parts/.space.cpp.swp
new file mode 100644
index 0000000..5429ac9
Binary files /dev/null and b/src/parts/.space.cpp.swp differ
diff --git a/src/parts/Makefile-part b/src/parts/Makefile-part
new file mode 100644
index 0000000..31e90db
--- /dev/null
+++ b/src/parts/Makefile-part
@@ -0,0 +1,4 @@
+src += src/parts/pipes.cpp \
+		src/parts/space.cpp \
+		src/parts/temple.cpp \
+		src/parts/dist.cpp
diff --git a/src/parts/dist.cpp b/src/parts/dist.cpp
new file mode 100644
index 0000000..3156130
--- /dev/null
+++ b/src/parts/dist.cpp
@@ -0,0 +1,136 @@
+#include <iostream>
+#include "dist.hpp"
+
+using namespace std;
+
+static Vertex *varray;
+static Index *iarray;
+static TexCoord *tc;
+
+static const int subd = 50;
+static const int vsz = subd + 2;
+static int vcount, icount;
+
+DistFx::DistFx() : Part("distort") {
+	set_clear(false);
+
+	TriMesh mesh;
+	create_plane(&mesh, Vector3(0, 0, -1), Vector2(1.0, 1.0), subd);
+
+	vcount = mesh.get_vertex_array()->get_count();
+	icount = mesh.get_index_array()->get_count();
+
+	varray = new Vertex[vcount];
+	Vertex *vptr = varray;
+	for(int i=0; i<vcount; i++) {
+		*vptr = mesh.get_vertex_array()->get_data()[i];
+		vptr->pos += Vector3(0.5, 0.5);
+		vptr++;
+	}
+
+	iarray = new Index[icount];
+	Index *iptr = iarray;
+	for(int i=0; i<icount; i++) {
+		*iptr++ = mesh.get_index_array()->get_data()[i];
+	}
+
+	tc = new TexCoord[vcount];
+	TexCoord *tptr = tc;
+	for(int y=0; y<vsz; y++) {
+		for(int x=0; x<vsz; x++) {
+			*tptr++ = TexCoord((float)x / (float)(vsz - 1), 1.0 - (float)y / (float)(vsz - 1));
+		}
+	}
+	
+}
+
+DistFx::~DistFx() {
+	delete [] varray;
+	delete [] iarray;
+}
+
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+#define MAX(a, b)	((a) > (b) ? (a) : (b))
+#define CLAMP(x, a, b)	MIN(MAX(x, a), b)
+
+void DistFx::draw_part() {
+	// update
+	float t = (float)time / 1000.0;
+	float cost = cos(t * 2.0);
+	float sint = sin(t * 2.0);
+	
+	Vertex *vptr = varray;
+	for(int y=0; y<vsz; y++) {
+		for(int x=0; x<vsz; x++) {
+			float u = (float)x / (float)(vsz - 1);
+			float v = (1.0 - (float)y / (float)(vsz - 1));
+
+			float pu = u * 4.0;
+			float pv = v * 4.0;
+
+			float du = sin(pu + cost) * 10.0 +
+						cos((pu + sint) * 2.0) * 5.0 +
+						sin((pu + cost) * 3.0) * 3.33 + 
+						cos((pv + sint) * 2.0) * 5.0 +
+						sin(pv + cost) * 10.0;
+
+			float dv = cos(pv + sint) * 10.0 +
+						sin((pv + cost) * 2.0) * 5.0 +
+						cos((pu + sint) * 3.0) * 3.33 +
+						sin((pv + cost) * 2.0) * 6.0;
+
+			du *= sin(t / 2.0) * 0.003 + 0.003;
+			dv *= sin(t / 2.0) * 0.003 + 0.003;
+
+			vptr->tex[0].u = CLAMP(u + du, 0.0, 1.0);
+			vptr->tex[0].v = CLAMP(v + dv, 0.0, 1.0);
+			vptr++;
+		}
+	}
+	
+	// draw
+	set_matrix(XFORM_TEXTURE, Matrix4x4::identity_matrix);
+	load_xform_matrices();
+
+	glMatrixMode(GL_MODELVIEW);
+	glLoadIdentity();
+	glMatrixMode(GL_PROJECTION);
+	glLoadIdentity();
+	glOrtho(0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
+
+	set_lighting(false);
+	set_zbuffering(false);
+	set_backface_culling(false);
+
+	set_texture(0, dsys::tex[0]);
+	enable_texture_unit(0);
+
+	glBegin(GL_TRIANGLES);
+	glColor4f(1.0, 1.0, 1.0, 0.5);
+	for(int i=0; i<icount; i++) {
+		Vertex *vptr = &varray[iarray[i]];
+		glTexCoord2f(tc[iarray[i]].u, tc[iarray[i]].v);
+		glVertex3f(vptr->pos.x, vptr->pos.y, vptr->pos.z);
+	}
+	glEnd();
+
+	set_alpha_blending(true);
+	set_blend_func(BLEND_SRC_ALPHA, BLEND_ONE_MINUS_SRC_ALPHA);
+
+	glBegin(GL_TRIANGLES);
+	glColor4f(1.0, 1.0, 1.0, 0.7);
+	for(int i=0; i<icount; i++) {
+		Vertex *vptr = &varray[iarray[i]];
+		glTexCoord2f(vptr->tex[0].u, vptr->tex[0].v);
+		glVertex3f(vptr->pos.x, vptr->pos.y, vptr->pos.z);
+	}
+	glEnd();
+
+	set_alpha_blending(false);
+
+	disable_texture_unit(0);
+	
+	set_lighting(true);
+	set_zbuffering(true);
+	set_backface_culling(true);
+}
diff --git a/src/parts/dist.hpp b/src/parts/dist.hpp
new file mode 100644
index 0000000..fda880e
--- /dev/null
+++ b/src/parts/dist.hpp
@@ -0,0 +1,16 @@
+#ifndef _DISTORTION_FX_HPP_
+#define _DISTORTION_FX_HPP_
+
+#include "3dengfx/3dengfx.hpp"
+#include "dsys/demosys.hpp"
+
+class DistFx : public dsys::Part {
+private:
+	virtual void draw_part();
+
+public:
+	DistFx();
+	virtual ~DistFx();
+};
+
+#endif	// _DISTORTION_FX_HPP_
diff --git a/src/parts/parts.hpp b/src/parts/parts.hpp
new file mode 100644
index 0000000..9cdb762
--- /dev/null
+++ b/src/parts/parts.hpp
@@ -0,0 +1,16 @@
+#ifndef _PARTS_HPP_
+#define _PARTS_HPP_
+
+#include "pipes.hpp"
+#include "space.hpp"
+#include "temple.hpp"
+#include "dist.hpp"
+
+#define INIT_PARTS()	{\
+	parts.push_back(new PipesPart());\
+	parts.push_back(new SpacePart());\
+	parts.push_back(new TemplePart());\
+	parts.push_back(new DistFx());\
+}
+
+#endif	// _PARTS_HPP_
diff --git a/src/parts/pipes.cpp b/src/parts/pipes.cpp
new file mode 100644
index 0000000..ec45a28
--- /dev/null
+++ b/src/parts/pipes.cpp
@@ -0,0 +1,302 @@
+#include "pipes.hpp"
+
+static void make_skycube(Scene *scene);
+static void make_title(Scene *scene);
+static void make_vol(Scene *scene);
+static void render_vol(unsigned long time);
+static void render_presents(float t);
+
+#define PSYS_COUNT	4
+static ParticleSystem *ps[PSYS_COUNT];
+static ParticleSysParams *psys_params[PSYS_COUNT];
+static bool psys_removed = false;
+
+static const float inner_sz = 30.0;
+static const float outer_sz = 70.0;
+static const int sph_count = 50;
+static Object *inner_sph, *vol_sph[sph_count];
+
+static Texture *presents, *summer_hack;
+
+#define TITLE
+#define VOL
+#define ANIM
+
+PipesPart::PipesPart() : ScenePart("pipes", new Scene) {
+#ifdef ANIM
+	TargetCamera *cam = new TargetCamera(Vector3(0, 20, -150));
+#else
+	TargetCamera *cam = new TargetCamera(Vector3(0, 10, -150));
+#endif
+	cam->set_fov(DEG_TO_RAD(50));
+	scene->add_camera(cam);
+
+	PointLight *lt = new PointLight(Vector3(10, 100, -100));
+	scene->add_light(lt);
+
+	scene->set_ambient_light(0.1);
+
+	make_skycube(scene);
+#ifdef TITLE
+	make_title(scene);
+#endif
+#ifdef VOL
+	make_vol(scene);
+#endif
+
+	presents = get_texture("data/img/presents.png");
+	summer_hack = get_texture("data/img/title.png");
+
+	MotionController cam_x, cam_y, cam_z;
+
+	cam_x = MotionController(CTRL_COS, TIME_FREE);
+	cam_x.set_sin_func(0.5, 40);
+	cam_x.set_control_axis(CTRL_X);
+
+	cam_z = MotionController(CTRL_SIN, TIME_FREE);
+	cam_z.set_sin_func(1.5, 100);
+	cam_z.set_control_axis(CTRL_Z);
+
+	cam_y = MotionController(CTRL_SIN, TIME_FREE);
+	cam_y.set_sin_func(0.35, 20);
+	cam_y.set_control_axis(CTRL_Y);
+
+#ifdef ANIM
+	cam->add_controller(cam_x, CTRL_TRANSLATION);
+	//cam->add_controller(cam_z, CTRL_TRANSLATION);
+	cam->add_controller(cam_y, CTRL_TRANSLATION);
+#endif
+
+	scene->set_background(0.5);
+
+	// make sure the volume and title textures are resident from the beginning
+	dsys::overlay(get_texture("data/img/cracks.jpg"), Vector2(0, 0), Vector2(0.1, 0.1), 0.0);
+	dsys::overlay(get_texture("data/img/volcano.jpg"), Vector2(0, 0), Vector2(0.1, 0.1), 0.0);
+	dsys::overlay(get_texture("data/img/title.png"), Vector2(0, 0), Vector2(0.1, 0.1), 0.0);
+	dsys::overlay(get_texture("data/img/presents.png"), Vector2(0, 0), Vector2(0.1, 0.1), 0.0);
+
+	//scene->render_sequence(0, 15000, 25);
+	//exit(0);
+}
+
+PipesPart::~PipesPart() {
+	//if(psys_removed) {
+		for(int i=0; i<PSYS_COUNT; i++) {
+			delete ps[i];
+		}
+	//}
+}
+
+static const float pexpl_time = half_pi * 8.0;
+static const float pexpl_full_time = pexpl_time + half_pi * 3.4;
+static const float pstop_time = pexpl_time + half_pi * 6.0;
+
+void PipesPart::draw_part() {
+	float t = (float)time / 1000.0f;
+	
+	ScenePart::draw_part();
+#ifdef VOL
+	render_vol(time);
+#endif
+
+#ifdef TITLE
+	if(t < pstop_time) {
+		float p = t / 4.0;
+		if(t > pexpl_time) {
+			static unsigned long prev_drop;
+			bool drop = (time - prev_drop > 33.3333) ? true : false;
+			
+			//Fuzzy val(0.0, 1.75 * (1.0 - sin(fmod(p, half_pi) + half_pi)));
+			Fuzzy val(0.0, 1.75 * sin(fmod(p, half_pi)));
+			for(int i=0; i<PSYS_COUNT; i++) {
+				psys_params[i]->shoot_dir = FuzzyVec3(val, val, val);
+
+				if(t > pexpl_full_time && drop) {
+					psys_params[i]->birth_rate.num *= 0.8;
+				}
+			}
+		}
+
+		p = t / 15.0;
+		for(int i=0; i<PSYS_COUNT; i++) {
+			float val = p > 0.5 ? 0.5 : p;
+			float range = val * 2.0;
+			psys_params[i]->spawn_offset_curve_area = Fuzzy(val, range);
+
+			ps[i]->update();
+			ps[i]->draw();
+		}
+	}/* else {
+		if(!psys_removed) {
+			//for(int i=0; i<PSYS_COUNT; i++) {
+			//	scene->remove_particle_sys(ps[i]);
+			//}
+			psys_removed = true;
+		}
+	}*/
+#endif	// TITLE
+
+	render_presents(t);
+}
+
+static void make_skycube(Scene *scene) {
+	const float size = 1000;
+	Object *face[6];
+	Texture *tex[6];
+
+	face[0] = new ObjPlane(Vector3(0, -1, 0), Vector2(size, size), 0);
+	face[0]->translate(Vector3(0, size / 2, 0));
+	tex[0] = get_texture("data/img/nebula_py.jpg");
+	
+	face[1] = new ObjPlane(Vector3(0, 1, 0), Vector2(size, size), 0);
+	face[1]->translate(Vector3(0, -size / 2, 0));
+	tex[1] = get_texture("data/img/nebula_ny.jpg");
+
+	face[2] = new ObjPlane(Vector3(0, 0, -1), Vector2(size, size), 0);
+	face[2]->translate(Vector3(0, 0, size / 2));
+	tex[2] = get_texture("data/img/nebula_pz.jpg");
+	
+	face[3] = new ObjPlane(Vector3(0, 0, 1), Vector2(size, size), 0);
+	face[3]->translate(Vector3(0, 0, -size / 2));
+	tex[3] = get_texture("data/img/nebula_nz.jpg");
+
+	face[4] = new ObjPlane(Vector3(-1, 0, 0), Vector2(size, size), 0);
+	face[4]->translate(Vector3(size / 2, 0, 0));
+	tex[4] = get_texture("data/img/nebula_px.jpg");
+	
+	face[5] = new ObjPlane(Vector3(1, 0, 0), Vector2(size, size), 0);
+	face[5]->translate(Vector3(-size / 2, 0, 0));
+	tex[5] = get_texture("data/img/nebula_nx.jpg");
+
+	for(int i=0; i<6; i++) {
+		face[i]->apply_xform();
+		face[i]->set_rotation(Vector3(0, quarter_pi, 0));
+	
+		Material *mat = face[i]->get_material_ptr();
+		mat->emissive_color = 1.0;
+		add_texture(tex[i]);
+		mat->set_texture(tex[i], TEXTYPE_DIFFUSE);
+		face[i]->set_texture_addressing(TEXADDR_CLAMP);
+		scene->add_object(face[i]);
+	}
+}
+
+
+void make_title(Scene *scene) {
+
+	for(int i=0; i<PSYS_COUNT; i++) {
+		char curve_fname[64];
+		sprintf(curve_fname, "data/curve%02d", i);
+		
+		Curve *curve = load_curve(curve_fname);
+		curve->set_arc_parametrization(true);
+		ps[i] = new ParticleSystem("data/thelab.psys");
+		psys_params[i] = ps[i]->get_params();
+		psys_params[i]->spawn_offset_curve = curve;
+		if(i == 2) {
+			psys_params[i]->birth_rate.num /= 1.5;
+		}
+		if(i == 1) {
+			psys_params[i]->birth_rate.num /= 6.0;
+		}
+		if(i == 3) {
+			psys_params[i]->birth_rate.num *= 1.3;
+		}
+		//scene->add_particle_sys(ps[i]);
+	}
+}
+
+static float vol_alpha[sph_count];
+
+static void make_vol(Scene *scene) {
+	Material *mat;
+
+	MotionController rot1(CTRL_LIN, TIME_FREE);
+	MotionController rot2(CTRL_SIN, TIME_FREE);
+	
+	rot1.set_control_axis(CTRL_XY);
+	rot1.set_slope(1.08);
+	rot2.set_control_axis(CTRL_YZ);
+	rot2.set_slope(1.439999);
+	rot2.set_sin_func(1.5, 1.5);
+
+	inner_sph = new ObjSphere(inner_sz, 4);
+	mat = inner_sph->get_material_ptr();
+	mat->set_texture(get_texture("data/img/volcano.jpg"), TEXTYPE_DIFFUSE);
+	inner_sph->add_controller(rot1, CTRL_ROTATION);
+	inner_sph->add_controller(rot2, CTRL_ROTATION);
+	//scene->add_object(obj);
+
+	for(int i=0; i<sph_count; i++) {
+		float t = (float)i / (float)sph_count;
+		vol_sph[i] = new ObjSphere(inner_sz + (outer_sz - inner_sz) * t, 3);
+		mat = vol_sph[i]->get_material_ptr();
+		mat->set_texture(get_texture("data/img/cracks.jpg"), TEXTYPE_DIFFUSE);
+		mat->alpha = vol_alpha[i] = 1.0 - t;
+		mat->emissive_color = 1.0;
+		vol_sph[i]->set_handle_blending(false);
+		vol_sph[i]->set_blending_mode(BLEND_SRC_ALPHA, BLEND_ONE);
+		vol_sph[i]->set_blending(true);
+		vol_sph[i]->add_controller(rot1, CTRL_ROTATION);
+		vol_sph[i]->add_controller(rot2, CTRL_ROTATION);
+		vol_sph[i]->set_zwrite(false);
+		//scene->add_object(obj);
+	}
+}
+
+static void render_vol(unsigned long time) {
+	static const float leave_time = 34.0;
+	
+	float t = (float)time / 1000.0;
+	if(t < pexpl_full_time) return;
+
+	float sp = (t - pexpl_full_time) / (pstop_time - pexpl_full_time);
+	if(t > pstop_time) sp = 1.0;
+
+	// update
+	inner_sph->get_material_ptr()->alpha = sp;
+	for(int i=0; i<sph_count + 1; i++) {
+		Object *obj = i < sph_count ? vol_sph[i] : inner_sph;
+
+		obj->set_scaling(Vector3(sp, sp, sp));
+		if(i < sph_count) obj->get_material_ptr()->alpha = vol_alpha[i] * sp;
+
+		if(t > leave_time) {
+			float lt = t - leave_time;
+			obj->set_position(Vector3(0, pow(lt * 5.0, 2.0), 0));
+		}
+	}
+	
+	// render
+	set_front_face(ORDER_CCW);
+	set_zbuffering(false);
+	for(int i=0; i<sph_count; i++) {
+		vol_sph[i]->render(time);
+	}
+	set_zbuffering(true);
+	set_front_face(ORDER_CW);
+
+	inner_sph->render(time);
+
+	set_zbuffering(false);
+	for(int i=0; i<sph_count; i++) {
+		vol_sph[i]->render(time);
+	}
+	set_zbuffering(true);
+}
+
+static void render_presents(float t) {
+	static const float pres_start = pexpl_full_time - 0.5;
+	static const Vector2 pres_sz(0.35, 0.12);
+	static const Vector2 hack_sz(0.55, 0.12);
+	
+	t -= pres_start;
+
+	if(t < 0.0 || t > pres_start + 5.0) return;
+
+	float x = pow((t * 0.2) - 1.0, 3.0);
+	Vector2 pres_pos(-x + 0.5, 0.3);
+	Vector2 hack_pos(x + 0.5, 0.7);
+	dsys::overlay(presents, pres_pos - pres_sz / 2.0, pres_pos + pres_sz / 2.0, 1.0);
+	dsys::overlay(summer_hack, hack_pos - hack_sz / 2.0, hack_pos + hack_sz / 2.0, 1.0);
+}
diff --git a/src/parts/pipes.hpp b/src/parts/pipes.hpp
new file mode 100644
index 0000000..7829013
--- /dev/null
+++ b/src/parts/pipes.hpp
@@ -0,0 +1,16 @@
+#ifndef _PIPES_PART_HPP_
+#define _PIPES_PART_HPP_
+
+#include "3dengfx/3dengfx.hpp"
+#include "dsys/demosys.hpp"
+
+class PipesPart : public dsys::ScenePart {
+private:
+	virtual void draw_part();
+
+public:
+	PipesPart();
+	virtual ~PipesPart();
+};
+
+#endif	// _PIPES_PART_HPP_
diff --git a/src/parts/space.cpp b/src/parts/space.cpp
new file mode 100644
index 0000000..4bf9ba9
--- /dev/null
+++ b/src/parts/space.cpp
@@ -0,0 +1,182 @@
+#include "space.hpp"
+
+static void make_skycube(Scene *scene);
+static void make_particles(Scene *scene);
+static void draw_design(float t);
+
+static const float des_sz = 0.1;
+static const int des_count = (int)(1.0 / des_sz);
+static const float des_offs = des_sz;
+
+#define ANIM
+
+SpacePart::SpacePart() : ScenePart("space", new Scene) {
+#ifdef ANIM
+	TargetCamera *cam = new TargetCamera(Vector3(0, 15, -5));
+#else
+	TargetCamera *cam = new TargetCamera(Vector3(0, 30, -250));
+#endif
+	cam->set_fov(DEG_TO_RAD(50));
+	scene->add_camera(cam);
+
+	PointLight *lt = new PointLight(Vector3(10, 100, -100));
+	scene->add_light(lt);
+
+	scene->set_ambient_light(0.1);
+
+	make_skycube(scene);
+	make_particles(scene);
+
+	MotionController cam_x, cam_y, cam_z;
+
+	cam_x = MotionController(CTRL_COS, TIME_FREE);
+	cam_x.set_sin_func(0.7, 250);
+	cam_x.set_control_axis(CTRL_X);
+
+	cam_z = MotionController(CTRL_SIN, TIME_FREE);
+	cam_z.set_sin_func(0.9, 120);
+	cam_z.set_control_axis(CTRL_Z);
+
+	cam_y = MotionController(CTRL_SIN, TIME_FREE);
+	cam_y.set_sin_func(0.5, 40);
+	cam_y.set_control_axis(CTRL_Y);
+
+#ifdef ANIM
+	cam->add_controller(cam_x, CTRL_TRANSLATION);
+	cam->add_controller(cam_z, CTRL_TRANSLATION);
+	cam->add_controller(cam_y, CTRL_TRANSLATION);
+#endif
+
+	scene->set_background(0.5);
+
+	//scene->render_sequence(0, 15000, 25);
+	//exit(0);
+}
+
+SpacePart::~SpacePart() {
+}
+
+void SpacePart::draw_part() {
+	float t = (float)time / 1000.0;
+
+	ScenePart::draw_part();
+	draw_design(t);
+}
+
+static void make_skycube(Scene *scene) {
+	const float size = 1000;
+	Object *face[6];
+	Texture *tex[6];
+
+	face[0] = new ObjPlane(Vector3(0, -1, 0), Vector2(size, size), 0);
+	face[0]->translate(Vector3(0, size / 2, 0));
+	tex[0] = get_texture("data/img/nebula_py.jpg");
+	
+	face[1] = new ObjPlane(Vector3(0, 1, 0), Vector2(size, size), 0);
+	face[1]->translate(Vector3(0, -size / 2, 0));
+	tex[1] = get_texture("data/img/nebula_ny.jpg");
+
+	face[2] = new ObjPlane(Vector3(0, 0, -1), Vector2(size, size), 0);
+	face[2]->translate(Vector3(0, 0, size / 2));
+	tex[2] = get_texture("data/img/nebula_pz.jpg");
+	
+	face[3] = new ObjPlane(Vector3(0, 0, 1), Vector2(size, size), 0);
+	face[3]->translate(Vector3(0, 0, -size / 2));
+	tex[3] = get_texture("data/img/nebula_nz.jpg");
+
+	face[4] = new ObjPlane(Vector3(-1, 0, 0), Vector2(size, size), 0);
+	face[4]->translate(Vector3(size / 2, 0, 0));
+	tex[4] = get_texture("data/img/nebula_px.jpg");
+	
+	face[5] = new ObjPlane(Vector3(1, 0, 0), Vector2(size, size), 0);
+	face[5]->translate(Vector3(-size / 2, 0, 0));
+	tex[5] = get_texture("data/img/nebula_nx.jpg");
+
+	for(int i=0; i<6; i++) {
+		Material *mat = face[i]->get_material_ptr();
+		mat->emissive_color = 1.0;
+		add_texture(tex[i]);
+		mat->set_texture(tex[i], TEXTYPE_DIFFUSE);
+		face[i]->set_texture_addressing(TEXADDR_CLAMP);
+		scene->add_object(face[i]);
+	}
+}
+
+static void make_particles(Scene *scene) {
+	MotionController mc(CTRL_LIN, TIME_FREE);
+	mc.set_slope(1);
+	mc.set_control_axis(CTRL_Y);
+	
+	ParticleSystem *galaxy = new ParticleSystem("data/galaxy.psys");
+	galaxy->add_controller(mc, CTRL_ROTATION);
+	scene->add_particle_sys(galaxy);
+
+	galaxy = new ParticleSystem("data/galaxy.psys");
+	mc.set_origin(half_pi);
+	galaxy->add_controller(mc, CTRL_ROTATION);
+	scene->add_particle_sys(galaxy);
+
+	ParticleSystem *stars = new ParticleSystem("data/stars.psys");
+	scene->add_particle_sys(stars);
+}
+
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+
+static void draw_design(float t) {
+	static const float des_fill = 3.0;
+	static const float yoffs = des_sz * 0.8;// + des_sz / 8.0;
+	static const Vector2 szvec(des_sz / 2, des_sz / 2);
+	static const float phase_offs = half_pi;
+	static const float x = des_sz;
+	static const float bar_len = 0.95;
+
+	static const float leave_start = 19.0;
+	static const float leave_dur = 4.0;
+	
+	float fill_t;
+	if(t < leave_start) {
+		fill_t = MIN(1.0, t / des_fill);
+	} else {
+		fill_t = 1.0 - (t - leave_start) / leave_dur;
+	}
+
+	set_zbuffering(false);
+	set_lighting(false);
+	set_alpha_blending(true);
+	set_blend_func(BLEND_ONE_MINUS_DST_COLOR, BLEND_ZERO);
+
+	for(int i=0; i<des_count; i++) {
+		if((float)i / (float)des_count > fill_t) break;
+		Vector2 pos(x + des_sz * (float)i, yoffs);
+	
+		glMatrixMode(GL_MODELVIEW);
+		glLoadIdentity();
+		glTranslatef(pos.x, pos.y, 0);
+		glRotatef(100.0 * cos(t * 3.0 + (float)i * phase_offs), 0, 0, 1);
+		glTranslatef(-pos.x, -pos.y, 0);
+		
+		float sz = sin(t * 4.0 + (float)i * phase_offs) * 0.25 + 0.5;
+		draw_scr_quad(pos - szvec * sz, pos + szvec * sz, 1.0, false);
+		draw_scr_quad(pos - szvec * sz * 0.7, pos + szvec * sz * 0.7, 1.0, false);
+
+		pos = Vector2(1.0, 1.0) - pos;
+		glMatrixMode(GL_MODELVIEW);
+		glLoadIdentity();
+		glTranslatef(pos.x, pos.y, 0);
+		glRotatef(100.0 * cos(t * 3.0 + (float)i * phase_offs), 0, 0, 1);
+		glTranslatef(-pos.x, -pos.y, 0);
+		
+		draw_scr_quad(pos - szvec * sz, pos + szvec * sz, 1.0, false);
+		draw_scr_quad(pos - szvec * sz * 0.7, pos + szvec * sz * 0.7, 1.0, false);
+	}
+
+	float bar_rest = 1.0 - bar_len;
+	float bar_start = bar_rest / 2.0;
+	float bar_end = bar_start + bar_len * fill_t;
+	draw_scr_quad(Vector2(bar_start, yoffs - 0.002), Vector2(bar_end, yoffs + 0.002));
+	draw_scr_quad(Vector2(1.0 - bar_start, 1.0 - yoffs + 0.002), Vector2(1.0 - bar_end, 1.0 - yoffs - 0.002));
+	
+	set_alpha_blending(false);
+	set_lighting(true);
+	set_zbuffering(true);
+}
diff --git a/src/parts/space.hpp b/src/parts/space.hpp
new file mode 100644
index 0000000..5220625
--- /dev/null
+++ b/src/parts/space.hpp
@@ -0,0 +1,16 @@
+#ifndef _SPACE_PART_HPP_
+#define _SPACE_PART_HPP_
+
+#include "3dengfx/3dengfx.hpp"
+#include "dsys/demosys.hpp"
+
+class SpacePart : public dsys::ScenePart {
+private:
+	virtual void draw_part();
+
+public:
+	SpacePart();
+	virtual ~SpacePart();
+};
+
+#endif	// _SPACE_PART_HPP_
diff --git a/src/parts/temple.cpp b/src/parts/temple.cpp
new file mode 100644
index 0000000..e2a58e1
--- /dev/null
+++ b/src/parts/temple.cpp
@@ -0,0 +1,575 @@
+#include <cassert>
+#include "temple.hpp"
+#include "common/err_msg.h"
+
+static void make_skycube(Scene *scene);
+static void make_temple(Scene *scene);
+static bool make_pillars(Scene *scene);
+static void make_blobs(Scene *scene);
+
+static void move_blobs(scalar_t t);
+static scalar_t eval_func(const Vector3 &vec, scalar_t t);
+static Vector3 eval_normals(const Vector3 &vec, scalar_t t);
+
+static const int pillar_udiv = 16;
+
+ParticleSysParams *pp[4];
+
+Texture *greets;
+
+#define BLOB_COUNT	4
+static scalar_t bint[BLOB_COUNT] = {30, 28, 24, 15};
+static Vector3 pos[BLOB_COUNT];
+static ScalarField *sfield;
+static Object *blob;
+
+#define BLOBS
+#define TEMPLE
+#define EXIBITS
+
+#define ANIM
+
+TemplePart::TemplePart() : ScenePart("temple", new Scene) {
+#ifdef ANIM
+	TargetCamera *cam = new TargetCamera(Vector3(0, 52, 0), Vector3(0, 50, 0));
+#else
+	TargetCamera *cam = new TargetCamera(Vector3(0, 60, -160), Vector3(0, 50, 0));
+#endif
+	cam->set_fov(DEG_TO_RAD(50));
+	scene->add_camera(cam);
+
+	PointLight *lt;
+	
+	lt = new PointLight(Vector3(-200, 250, -400));
+	lt->set_intensity(0.9);
+	scene->add_light(lt);
+
+	lt = new PointLight(Vector3(200, 20, 0));
+	lt->set_intensity(0.6);
+	scene->add_light(lt);
+
+	lt = new PointLight(Vector3(20, 500, -400));
+	lt->set_intensity(0.5);
+	scene->add_light(lt);
+
+	scene->set_ambient_light(0.1);
+
+	make_skycube(scene);
+#ifdef TEMPLE
+	make_temple(scene);
+#endif
+#ifdef EXIBITS
+	make_pillars(scene);
+#endif
+#ifdef BLOBS
+	make_blobs(scene);
+#endif
+	
+	MotionController ctrl;
+
+#ifdef ANIM
+	ctrl = MotionController(CTRL_COS, TIME_FREE);
+	ctrl.set_sin_func(0.58, 250);
+	ctrl.set_control_axis(CTRL_X);
+	cam->add_controller(ctrl, CTRL_TRANSLATION);
+
+	ctrl = MotionController(CTRL_SIN, TIME_FREE);
+	ctrl.set_sin_func(1.17, 125);
+	ctrl.set_control_axis(CTRL_X);
+	cam->add_controller(ctrl, CTRL_TRANSLATION);
+
+	ctrl = MotionController(CTRL_SIN, TIME_FREE);
+	ctrl.set_sin_func(0.49, 220);
+	ctrl.set_control_axis(CTRL_Z);
+	cam->add_controller(ctrl, CTRL_TRANSLATION);
+
+	ctrl = MotionController(CTRL_COS, TIME_FREE);
+	ctrl.set_sin_func(0.99, 105);
+	ctrl.set_control_axis(CTRL_Z);
+	cam->add_controller(ctrl, CTRL_TRANSLATION);
+
+	ctrl = MotionController(CTRL_SIN, TIME_FREE);
+	ctrl.set_sin_func(0.45, 50);
+	ctrl.set_control_axis(CTRL_Y);
+	cam->add_controller(ctrl, CTRL_TRANSLATION);
+#endif
+
+	greets = get_texture("data/img/greets.png");
+
+	scene->set_background(0.5);
+
+	// in order to precalculate the cubemaps...
+	scene->render(0);
+
+	//scene->render_sequence(0, 15000, 25);
+	//exit(0);
+}
+
+TemplePart::~TemplePart() {}
+
+static void make_skycube(Scene *scene) {
+	const float size = 5000;
+	Object *face[6];
+	Texture *tex[6];
+
+	face[0] = new ObjPlane(Vector3(0, -1, 0), Vector2(size, size), 0);
+	face[0]->translate(Vector3(0, size / 2, 0));
+	tex[0] = get_texture("data/img/py.jpg");
+	
+	face[1] = new ObjPlane(Vector3(0, 1, 0), Vector2(size, size), 0);
+	face[1]->translate(Vector3(0, -size / 2, 0));
+	tex[1] = get_texture("data/img/ny.jpg");
+
+	face[2] = new ObjPlane(Vector3(0, 0, -1), Vector2(size, size), 0);
+	face[2]->translate(Vector3(0, 0, size / 2));
+	tex[2] = get_texture("data/img/pz.jpg");
+	
+	face[3] = new ObjPlane(Vector3(0, 0, 1), Vector2(size, size), 0);
+	face[3]->translate(Vector3(0, 0, -size / 2));
+	tex[3] = get_texture("data/img/nz.jpg");
+
+	face[4] = new ObjPlane(Vector3(-1, 0, 0), Vector2(size, size), 0);
+	face[4]->translate(Vector3(size / 2, 0, 0));
+	tex[4] = get_texture("data/img/px.jpg");
+	
+	face[5] = new ObjPlane(Vector3(1, 0, 0), Vector2(size, size), 0);
+	face[5]->translate(Vector3(-size / 2, 0, 0));
+	tex[5] = get_texture("data/img/nx.jpg");
+
+	for(int i=0; i<6; i++) {
+		Material *mat = face[i]->get_material_ptr();
+		mat->emissive_color = 1.0;
+		add_texture(tex[i]);
+		mat->set_texture(tex[i], TEXTYPE_DIFFUSE);
+		face[i]->set_texture_addressing(TEXADDR_CLAMP);
+		scene->add_object(face[i]);
+	}
+}
+
+void make_temple(Scene *scene) {
+	const Vector3 pillar_body[] = {
+		Vector3(8, 4, 0),
+		Vector3(12, 40, 0),
+		Vector3(11, 70, 0),
+		Vector3(4, 106, 0)
+	};
+
+	const Vector3 pillar_pos[] = {
+		Vector3(-50, 0, -50),
+		Vector3(-50, 0, 50),
+		Vector3(50, 0, 50),
+		Vector3(50, 0, -50)
+	};
+
+
+	TriMesh pillar_mesh;
+	create_revolution(&pillar_mesh, pillar_body, 4, pillar_udiv, 10);
+
+	Material pillar_mat;
+	pillar_mat.specular_color = 1.0;
+	pillar_mat.specular_power = 70.0;
+	pillar_mat.set_texture(get_texture("data/img/benedeti.jpg"), TEXTYPE_DIFFUSE);
+
+	for(int i=0; i<4; i++) {
+		Object *obj = new Object;
+		obj->get_mesh() = pillar_mesh;
+		obj->set_position(pillar_pos[i]);
+		*obj->get_material_ptr() = pillar_mat;
+		scene->add_object(obj);
+
+		obj = new ObjCylinder(13.0, 5.0, true, pillar_udiv);
+		obj->set_position(pillar_pos[i] + Vector3(0, 6, 0));
+		*obj->get_material_ptr() = pillar_mat;
+		scene->add_object(obj);
+		
+		obj = new ObjCylinder(16.0, 5.0, true, pillar_udiv);
+		obj->set_position(pillar_pos[i] + Vector3(0, 2, 0));
+		*obj->get_material_ptr() = pillar_mat;
+		scene->add_object(obj);
+
+		obj = new ObjTorus(2.0, 5.0);
+		obj->set_position(pillar_pos[i] + Vector3(0, 102, 0));
+		*obj->get_material_ptr() = pillar_mat;
+		scene->add_object(obj);
+
+		obj = new ObjTorus(2.0, 7);
+		obj->set_position(pillar_pos[i] + Vector3(0, 105, 0));
+		*obj->get_material_ptr() = pillar_mat;
+		scene->add_object(obj);
+	}
+
+
+	// floor
+	Object *obj = new ObjCylinder(95.0, 5, true, 25);
+	obj->set_position(Vector3(0, -2.5, 0));
+	obj->get_material_ptr()->set_texture(get_texture("data/img/brick1_base1.jpg"), TEXTYPE_DIFFUSE);
+	scene->add_object(obj);
+
+	obj = new ObjCylinder(103.0, 5, true, 25);
+	obj->set_position(Vector3(0, -7.5, 0));
+	obj->get_material_ptr()->set_texture(get_texture("data/img/brick1_base1.jpg"), TEXTYPE_DIFFUSE);
+	scene->add_object(obj);
+
+	// ceiling
+	const Vector3 ceil[] = {
+		Vector3(90, 0, 0),
+		Vector3(92, 4, 0),
+		Vector3(98, 14, 0),
+		Vector3(95, 20, 0),
+		Vector3(60, 40, 0),
+		Vector3(20, 60, 0),
+		Vector3(18, 64, 0),
+		Vector3(23, 66, 0),
+		Vector3(26, 67, 0),
+		Vector3(25, 70, 0),
+		Vector3(10, 78, 0),
+		Vector3(4, 84, 0),
+		Vector3(2, 90, 0),
+		Vector3(0.1, 100, 0)
+	};
+
+	obj = new Object;
+	create_revolution(obj->get_mesh_ptr(), ceil, 14, 24, 28);
+	obj->set_pivot(Vector3(0, 40, 0));
+	obj->set_position(Vector3(0, 110, 0));
+	obj->calculate_normals();
+
+	Texture *cubemap = new Texture(64, 64, TEX_CUBE);
+	add_texture(cubemap);
+
+	Material *mat = obj->get_material_ptr();
+	mat->diffuse_color = mat->ambient_color = int_color(239, 190, 37) * 0.8;
+	mat->specular_color = int_color(255, 198, 43);
+	mat->specular_power = 40.0;
+	mat->set_texture(get_texture("data/img/brick1_base1.jpg"), TEXTYPE_DIFFUSE);
+	mat->set_texture(cubemap, TEXTYPE_ENVMAP);
+	mat->auto_refl = false;
+	mat->env_intensity = 0.3;
+	scene->add_object(obj);
+
+	obj = new ObjCylinder(94, 5, true, 25);
+	obj->set_position(Vector3(0, 107.5, 0));
+	*obj->get_material_ptr() = *mat;
+	scene->add_object(obj);
+
+	// sting
+	const Vector3 sting[] = {
+		Vector3(12, 0, 0),
+		Vector3(9, 4, 0),
+		Vector3(5, 6, 0),
+		Vector3(2, 10, 0),
+		Vector3(0.01, 34, 0)
+	};
+
+	cubemap = new Texture(16, 16, TEX_CUBE);
+	add_texture(cubemap);
+
+	obj = new Object;
+	create_revolution(obj->get_mesh_ptr(), sting, 5, 10, 10);
+	obj->set_pivot(Vector3(0, 5, 0));
+	mat = obj->get_material_ptr();
+	mat->diffuse_color = mat->ambient_color = 0.1;
+	mat->specular_color = 1.0;
+	mat->specular_power = 80.0;
+	mat->set_texture(cubemap, TEXTYPE_ENVMAP);
+	mat->auto_refl = false;
+	scene->add_object(obj);
+
+
+	// floor-floor
+	obj = new ObjPlane(Vector3(0, 1, 0), Vector2(1200, 1200), 6);
+	obj->set_position(Vector3(0, -10, 0));
+	mat = obj->get_material_ptr();
+	mat->set_texture(get_texture("data/img/mud.jpg"), TEXTYPE_DIFFUSE);
+	mat->tmat[TEXTYPE_DIFFUSE].set_scaling(Vector3(4.0, 4.0, 4.0));
+	scene->add_object(obj);
+}
+
+static bool make_pillars(Scene *scene) {
+	static const Vector3 ppos[] = {
+		Vector3(-150, -7.5, 0),
+		Vector3(0, -7.5, 150),
+		Vector3(150, -7.5, 0),
+		Vector3(0, -7.5, -150)
+	};
+
+	static const char *label_tex[] = {
+		"data/img/label_dolphin.png",
+		"data/img/label_duck.png",
+		"data/img/label_face.png",
+		"data/img/label_teapot.png",
+		0
+	};
+	
+	Curve *curve = load_curve("data/pillar.curve");
+
+	for(int i=0; i<4; i++) {
+		char name[128];
+		sprintf(name, "pillar%02d", i);
+		
+		Object *obj = new Object;
+		obj->name = name;
+		Material *mat = obj->get_material_ptr();
+		create_revolution(obj->get_mesh_ptr(), *curve, pillar_udiv, 10);
+		obj->set_scaling(Vector3(2.3, 1.8, 2.3));
+		obj->apply_xform();
+		obj->normalize_normals();
+		obj->set_position(ppos[i]);
+
+		mat->specular_color = 1.0;
+		mat->specular_power = 70.0;
+		mat->set_texture(get_texture("data/img/benedeti.jpg"), TEXTYPE_DIFFUSE);
+		scene->add_object(obj);
+
+		obj = new ObjCylinder(14.0, 3.0);
+		obj->set_position(ppos[i]);
+		mat = obj->get_material_ptr();
+		mat->specular_color = 1.0;
+		mat->specular_power = 70.0;
+		mat->set_texture(get_texture("data/img/benedeti.jpg"), TEXTYPE_DIFFUSE);
+		scene->add_object(obj);
+
+
+		Texture *cubemap = new Texture(64, 64, TEX_CUBE);
+		add_texture(cubemap);
+		
+		// put exibit
+		switch(i) {
+		case 1:
+			{
+				Scene *scene = load_scene("data/geom/duck.3ds");
+				if(!scene) {
+					error("failed to load duck.3ds");
+					return false;
+				}
+				obj = scene->get_object("Object03");
+				assert(obj);
+				scene->remove_object(obj);
+				delete scene;
+
+				obj->set_scaling(Vector3(0.009, 0.009, 0.009));
+				obj->set_auto_normalize(true);
+			
+				obj->set_position(ppos[i] + Vector3(0, 50, 0));
+				obj->set_rotation(Vector3(0, quarter_pi, 0));
+			}
+			break;
+
+		case 2:
+			{
+				Scene *scene = load_scene("data/geom/chrmface.3ds");
+				if(!scene) {
+					error("failed to load chrmface.3ds");
+					return false;
+				}
+				obj = scene->get_object("0main01");
+				assert(obj);
+				scene->remove_object(obj);
+				delete scene;
+
+				obj->set_scaling(Vector3(0.2, 0.2, 0.2));
+				obj->set_auto_normalize(true);
+
+				obj->set_position(ppos[i] + Vector3(4, 58, 0));				
+				obj->set_rotation(Vector3(0, half_pi, 0));
+
+				obj->get_material_ptr()->two_sided = true;
+			}
+			break;
+
+		case 0:
+			{
+				Scene *scene = load_scene("data/geom/dolphin.3ds");
+				if(!scene) {
+					error("failed to load dolphin.3ds");
+					return false;
+				}
+				obj = scene->get_object("Loft01");
+				assert(obj);
+				scene->remove_object(obj);
+				delete scene;
+
+				obj->set_scaling(Vector3(0.14, 0.14, 0.14));
+				obj->set_auto_normalize(true);
+
+				obj->set_position(ppos[i] + Vector3(0, 48, 0));
+				obj->set_rotation(Vector3(half_pi, half_pi, 0));
+			}
+			break;
+		
+		case 3:
+		default:
+			obj = new ObjTeapot(6.0);
+			obj->set_position(ppos[i] + Vector3(0, 43.7, 0));
+			break;
+		}
+
+		mat = obj->get_material_ptr();
+		mat->diffuse_color = mat->ambient_color = 0.1;
+		mat->specular_color = 1.0;
+		mat->specular_power = 80.0;
+		mat->set_texture(cubemap, TEXTYPE_ENVMAP);
+		mat->auto_refl = false;
+		
+		MotionController ctrl(CTRL_SIN, TIME_FREE);
+		ctrl.set_control_axis(CTRL_Y);
+		ctrl.set_sin_func(frand(0.5) + 2.3, 3.0, frand(pi));
+		obj->add_controller(ctrl, CTRL_TRANSLATION);
+		obj->translate(Vector3(0, 3.0, 0));
+		scene->add_object(obj);
+
+		ParticleSystem *label = new ParticleSystem("data/exibit.psys");
+		pp[i] = label->get_params();
+		pp[i]->halo = get_texture(label_tex[i]);
+		label->set_position(ppos[i] + Vector3(0, 80, 0));
+		scene->add_particle_sys(label);
+
+		ParticleSystem *trail = new ParticleSystem("data/test2.psys");
+		trail->set_position(ppos[i] + Vector3(0, 50, 0));
+
+		ctrl = MotionController(CTRL_SIN, TIME_FREE);
+		ctrl.set_control_axis(CTRL_X);
+		ctrl.set_sin_func(5.0, 30.0, (float)i * half_pi);
+		trail->add_controller(ctrl, CTRL_TRANSLATION);
+
+		ctrl = MotionController(CTRL_COS, TIME_FREE);
+		ctrl.set_control_axis(CTRL_Z);
+		ctrl.set_sin_func(5.0, 30.0, (float)i * half_pi);
+		trail->add_controller(ctrl, CTRL_TRANSLATION);
+
+		ctrl = MotionController(CTRL_SIN, TIME_FREE);
+		ctrl.set_control_axis(CTRL_Y);
+		ctrl.set_sin_func(1.8, 25.0, (float)i * half_pi);
+		trail->add_controller(ctrl, CTRL_TRANSLATION);
+		
+		scene->add_particle_sys(trail);
+	}
+
+	delete curve;
+	return true;
+}
+
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+#define MAX(a, b)	((a) > (b) ? (a) : (b))
+
+void TemplePart::draw_part() {
+	static const float greets_time = 18.0;
+	float t = (float)time / 1000.0f;
+
+#ifdef EXIBITS
+	float yoffs = MIN(150.0, t * 25.0);
+	for(int i=0; i<4; i++) {
+		pp[i]->spawn_offset = FuzzyVec3(0.0, yoffs, 0.0);
+	}
+#endif
+
+#ifdef BLOBS
+	move_blobs(t);
+	sfield->triangulate(blob->get_mesh_ptr(), 1.0, t, true);
+#endif
+
+	scene->render(time);
+
+	if(t > greets_time) {
+		t *= 1.8;
+
+		float a2 = sin(t);
+    	a2 = a2 > 0 ? -a2 : a2;
+    	a2 += 1.0f;
+    	a2 *= frand(0.8) + 0.2;
+
+    	float a1 = a2;//(cos(t)+1)/2;
+    	a1 *= (sin(t / 2.0) > 0) ? 1.0 : 0.0;     // modulate with step func
+    	a1 *= sin(t / 2.0);
+
+		Color col(1.0, 1.0, 1.0, MAX(a2, a1) + 0.4);
+	
+		set_matrix(XFORM_TEXTURE, Matrix4x4::identity_matrix);
+		load_xform_matrices();
+
+		for(int i=0; i<10; i++) {
+			set_alpha_blending(true);
+			set_blend_func(BLEND_SRC_ALPHA, BLEND_ONE);
+			float off = (float)i * 0.0025;
+			col.a *= 0.8;
+			dsys::overlay(greets, Vector2(-off, -off), Vector2(1 + off, 1 + off), col, 0, false);
+			set_alpha_blending(false);
+		}
+	}
+}
+
+static void move_blobs(scalar_t t) {
+	Vector3 offs = Vector3(0, 2.75 * cos(t*1.5) + 2.25, 0);
+	pos[0] = offs + Vector3(sin(t), cos(t), 0) * 8.0;
+	pos[1] = offs + Vector3(sin(t*2) + cos(t*4)/2, sin(t*1.5), cos(t)) * 8.0;
+	pos[2] = offs + Vector3(cos(t*1.8)/2 + sin(t)/4, cos(t*2) + sin(t), sin(t*1.3)) * 8.0;
+	pos[3] = offs + Vector3(sin(t*1.2)/3 + cos(t*1.8)/3.5, sin(cos(t*2) + t)*1.5, cos(t*2.0)/2) * 7.0;
+
+	//cyl = Vector2(cos(t), sin(t*2.0)) * 8.0;
+}
+
+static scalar_t eval_func(const Vector3 &vec, scalar_t t) {
+	scalar_t val = 0.0;
+	for(int i=0; i<BLOB_COUNT; i++) {
+		val += (bint[i] * 0.75) / (vec - pos[i]).length_sq();
+	}
+
+	return val;
+}
+
+//#define NORMAL_GRAD
+static Vector3 eval_normals(const Vector3 &vec, scalar_t t) {
+#ifdef NORMAL_GRAD
+	const scalar_t diff = 0.1;
+
+	Vector3 grad;
+	grad.x = eval_func(vec + Vector3(diff, 0, 0), t) - eval_func(vec + Vector3(-diff, 0, 0), t);
+	grad.y = eval_func(vec + Vector3(0, diff, 0), t) - eval_func(vec + Vector3(0, -diff, 0), t);
+	grad.z = eval_func(vec + Vector3(0, 0, diff), t) - eval_func(vec + Vector3(0, 0, -diff), t);
+
+	return -grad.normalized();
+#else
+	Vector3 normal(0, 0, 0);
+	for(int i=0; i<BLOB_COUNT; i++) {
+		scalar_t len_sq = (vec - pos[i]).length_sq();
+		scalar_t len_quad = len_sq * len_sq;
+		normal += (vec - pos[i]) * 2.0 * bint[i] / len_quad;
+	}
+
+	return normal.normalized();
+#endif
+}
+
+
+static void make_blobs(Scene *scene) {
+	// scalar field
+	const int grid = 23;
+	const scalar_t field_size = 40.0;
+	const scalar_t hsz = field_size / 2.0;
+
+	sfield = new ScalarField(grid, Vector3(-hsz, -hsz, -hsz), Vector3(hsz, hsz, hsz));
+	sfield->set_evaluator(eval_func);
+	sfield->set_normal_evaluator(eval_normals);
+	
+	// rendering details
+	Texture *cubemap = new Texture(64, 64, TEX_CUBE);
+	add_texture(cubemap);
+
+	blob = new Object;
+	blob->name = "blob";
+	blob->set_dynamic(true);
+	blob->set_position(Vector3(0, 50, 0));
+	blob->set_scaling(Vector3(2.8, 2.6, 2.8));
+	blob->set_auto_normalize(true);
+
+	Material *mat = blob->get_material_ptr();
+	mat->ambient_color = mat->diffuse_color = Color(0.7, 0.4, 0.2) * 0.4;
+	mat->specular_color = Color(0.9, 0.7, 0.6);
+	mat->specular_power = 60.0;
+	mat->set_texture(cubemap, TEXTYPE_ENVMAP);
+	mat->env_intensity = 0.6;
+	mat->auto_refl = false;
+	mat->two_sided = true;
+
+	scene->add_object(blob);
+}
diff --git a/src/parts/temple.hpp b/src/parts/temple.hpp
new file mode 100644
index 0000000..01702d9
--- /dev/null
+++ b/src/parts/temple.hpp
@@ -0,0 +1,16 @@
+#ifndef _TEMPLE_PART_HPP_
+#define _TEMPLE_PART_HPP_
+
+#include "3dengfx/3dengfx.hpp"
+#include "dsys/demosys.hpp"
+
+class TemplePart : public dsys::ScenePart {
+private:
+	virtual void draw_part();
+
+public:
+	TemplePart();
+	virtual ~TemplePart();
+};
+
+#endif	// _TEMPLE_PART_HPP_
diff --git a/src/sdlvf/Makefile-part b/src/sdlvf/Makefile-part
new file mode 100644
index 0000000..61fd073
--- /dev/null
+++ b/src/sdlvf/Makefile-part
@@ -0,0 +1 @@
+csrc += src/sdlvf/sdlvf.c
diff --git a/src/sdlvf/sdlvf.c b/src/sdlvf/sdlvf.c
new file mode 100644
index 0000000..c96bc17
--- /dev/null
+++ b/src/sdlvf/sdlvf.c
@@ -0,0 +1,202 @@
+/*
+ * sdlvf - the SDL/vorbisfile sound system
+ * Copyright (C) 2004  Vasilis Vasaitis <vvas@hal.csd.auth.gr>
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include <vorbis/codec.h>
+#include <vorbis/vorbisfile.h>
+#include <SDL.h>
+#include <string.h>
+#include "sdlvf.h"
+
+#define SDL_SAMPLES 2048
+#define VORBISFILE_BUFFER 4096
+
+#define MIN(x, y) ((x) < (y) ? (x) : (y))
+
+/*
+ * The error strings array, which must be kept in sync with the error
+ * codes in sdlvf.h.
+ */
+static char *audio_errors[] = {
+	"playing audio",
+	"audio playback has been stopped",
+	"could not open input file for reading",
+	"input file does not appear to be an valid bitstream",
+	"unable to get information about the current logical bitstream",
+	"unable to open the audio device"
+};
+
+/*
+ * The static data shared by all SDL/vorbisfile functions. Note that
+ * there is never a need for more than one instance of the sound
+ * system, so this a reasonable thing to do.
+ */
+static OggVorbis_File audio_vf;
+static volatile int audio_stopped, audio_reopen;
+
+/*
+ * This function is called by SDL when more audio data is needed.
+ */
+static void audio_callback(void *userdata, Uint8 *stream, int len)
+{
+	/* state information kept across calls */
+	static char buf[VORBISFILE_BUFFER];
+	static int buflen = 0;
+	static int curr = -1;
+	static int stopped = 0, reopen = 0;
+
+	/* local variables */
+	int filled = 0;
+
+	/* check for exceptional conditions */
+	if (stopped || reopen) {
+		audio_stopped = stopped;
+		audio_reopen = reopen;
+		stopped = reopen = 0;
+	}
+	if (audio_stopped || audio_reopen) return;
+
+	/* check for leftovers in our buffer */
+	if (buflen != 0) {
+		int copy = MIN(buflen, len);
+		memcpy(stream, buf, copy);
+		filled += copy;
+		buflen -= copy;
+		memmove(buf, buf + copy, buflen);
+	}
+
+	/* main loop */
+	while (filled < len) {
+		int needed = MIN(VORBISFILE_BUFFER, len - filled), read;
+		int prev = curr;
+		do read = ov_read(&audio_vf, stream + filled, needed, 0, 2, 1, &curr);
+		while (read < 0);
+		if (read == 0) {
+			stopped = 1;
+			break;
+		}
+		if (curr != prev && prev != -1) {
+			memcpy(buf, stream + filled, buflen = read);
+			memset(stream + filled, 0, read);
+			reopen = 1;
+			break;
+		}
+		filled += read;
+	}
+}
+
+/*
+ * Opens the audio device based on the parameters of the current
+ * logical bitstream of the Ogg Vorbis file.
+ */
+static int audio_open(void)
+{
+	vorbis_info *vi;
+	SDL_AudioSpec as;
+
+	audio_stopped = 0;
+	audio_reopen = 0;
+
+	if ((vi = ov_info(&audio_vf, -1)) == NULL)
+		return SDLVF_BADSTREAM;
+	as.freq = vi->rate;
+	as.format = AUDIO_S16;
+	as.channels = vi->channels;
+	as.samples = SDL_SAMPLES;
+	as.callback = audio_callback;
+	as.userdata = NULL;
+	if (SDL_OpenAudio(&as, NULL) == -1)
+		return SDLVF_NOAUDIO;
+	SDL_PauseAudio(0);
+
+	return SDLVF_PLAYING;
+}
+
+/*
+ * Closes the audio device; provided for naming consistency.
+ */
+static void audio_close(void)
+{
+	SDL_CloseAudio();
+}
+
+/*
+ * Initialises SDL/vorbisfile and starts playing the Ogg Vorbis file
+ * <fname>. Returns zero (SDLVF_PLAYING) on success; any other value
+ * indicates an error.
+ */
+int sdlvf_init(const char *fname)
+{
+	FILE *f;
+	int result;
+	if ((f = fopen(fname, "rb")) == NULL)
+		return SDLVF_BADFILE;
+	if (ov_open(f, &audio_vf, NULL, 0) != 0) {
+		fclose(f);
+		return SDLVF_BADOGG;
+	}
+	if ((result = audio_open()) != SDLVF_PLAYING)
+		ov_clear(&audio_vf);
+	return result;
+}
+
+/*
+ * Checks the status of SDL/vorbisfile, reopening the audio device if
+ * necessary. It also reports whether audio playback has stopped.
+ */
+int sdlvf_check(void)
+{
+	if (audio_stopped)
+		return SDLVF_STOPPED;
+	if (audio_reopen) {
+		audio_close();
+		return audio_open();
+	}
+	return SDLVF_PLAYING;
+}
+
+/*
+ * Seeks to the specified <time> in seconds, taking care of locking
+ * and avoiding clicks and pops.
+ */
+int sdlvf_seek(double time)
+{
+	int result;
+	SDL_LockAudio();
+	result = ov_time_seek_lap(&audio_vf, time);
+	SDL_UnlockAudio();
+	return result;
+}
+
+/*
+ * Shuts down SDL/vorbisfile, closing the audio device and freeing
+ * data structures.
+ */
+void sdlvf_done(void)
+{
+	audio_close();
+	ov_clear(&audio_vf);
+}
+
+/*
+ * Returns the equivalent error string for the supplied error code.
+ */
+char *sdlvf_strerror(int error)
+{
+	return audio_errors[error];
+}
diff --git a/src/sdlvf/sdlvf.h b/src/sdlvf/sdlvf.h
new file mode 100644
index 0000000..f308c42
--- /dev/null
+++ b/src/sdlvf/sdlvf.h
@@ -0,0 +1,81 @@
+/*
+ * sdlvf - the SDL/vorbisfile sound system
+ * Copyright (C) 2004  Vasilis Vasaitis <vvas@hal.csd.auth.gr>
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _SDLVF_H
+#define _SDLVF_H
+
+/*
+ * All the possible return values of sdlvf_init() and sdlvf_check().
+ */
+enum {
+	SDLVF_PLAYING,              /* no error, playing file */
+	SDLVF_STOPPED,              /* no error, playback has stopped */
+	SDLVF_BADFILE,              /* could not open file specified */
+	SDLVF_BADOGG,               /* file specified is not valid */
+	SDLVF_BADSTREAM,            /* could not query current stream */
+	SDLVF_NOAUDIO               /* could not open audio device */
+};
+
+/*
+ * Make sure that everything works under C++ as well.
+ */
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ * Initialises the SDL/vorbisfile sound system. The audio subsystem of
+ * SDL must have been initialised before calling this function. On
+ * success, the Ogg Vorbis file specified is played, and zero
+ * (SDLVF_PLAYING) is returned. On error, the appropriate error code
+ * is returned.
+ */
+int sdlvf_init(const char *);
+
+/*
+ * This function should be called regularly, i.e. multiple times per
+ * second. It makes sure that everything is going as it should. It
+ * also checks whether playback has ended, returning the appropriate
+ * value if this is the case.
+ */
+int sdlvf_check(void);
+
+/*
+ * Seeks to the specified position in seconds.
+ */
+int sdlvf_seek(double);
+
+/*
+ * Stops playback and shuts down the sound system.
+ */
+void sdlvf_done(void);
+
+/*
+ * Returns the equivalent error string for the supplied error code.
+ */
+char *sdlvf_strerror(int);
+
+/*
+ * Close the C++-supporting block if needed.
+ */
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* !defined(_SDLVF_H) */
diff --git a/src/sumhack.cpp b/src/sumhack.cpp
new file mode 100644
index 0000000..75c6fa3
--- /dev/null
+++ b/src/sumhack.cpp
@@ -0,0 +1,197 @@
+#include <iostream>
+#include <vector>
+#include <cstdio>
+#include <cstdlib>
+#include "3dengfx/3dengfx.hpp"
+#include "dsys/demosys.hpp"
+#include "common/err_msg.h"
+#include "parts/parts.hpp"
+#include "events.hpp"
+#include "sdlvf/sdlvf.h"
+
+using namespace std;
+
+bool init();
+void update_gfx();
+void clean_up();
+
+dsys::Part *part;
+ntimer timer;
+
+std::vector<dsys::Part*> parts;
+
+bool music = true;
+
+bool render = false;
+int render_fps = 25;
+const char *render_path = "frames";
+
+unsigned long avg_frame_time, frame_count;
+
+const char *help_str = "command line options\n------------------\n"
+	"-r\trender the demo as a sequence of targa files\n"
+	"-f <fps>\tset the rendering framerate for -r above\n"
+	"-p <path>\tset the directory in which to render the sequnce for -r.\n"
+	"-m\ttoggle music.\n"
+	"-h\tthis help screen\n\n";
+
+int main(int argc, char **argv) {
+	for(int i=1; i<argc; i++) {
+		if(argv[i][0] == '-' && argv[i][2] == 0) {
+			switch(argv[i][1]) {
+			case 'r':
+				render = true;
+				break;
+
+			case 'f':
+				if(isdigit(argv[++i][0])) {
+					render_fps = atoi(argv[i]);
+				} else {
+					cerr << "invalid argument -f " << argv[i] << endl;
+					return EXIT_FAILURE;
+				}
+				break;
+
+			case 'p':
+				render_path = argv[++i];
+				break;
+
+			case 'm':
+				music = !music;
+				break;
+
+			default:
+				cout << help_str;
+				return argv[i][1] == 'h' ? 0 : EXIT_FAILURE;
+			}
+		} else {
+			cout << help_str;
+			return EXIT_FAILURE;
+		}
+	}
+
+	if(!init()) {
+		return EXIT_FAILURE;
+	}
+
+	return fxwt::main_loop();
+}
+
+bool init() {
+	GraphicsInitParameters *gip;
+	if(!(gip = load_graphics_context_config("3dengfx.conf"))) {
+		return false;
+	}
+
+	if(!create_graphics_context(*gip)) {
+		return false;
+	}
+
+	if(gip->fullscreen) {
+		SDL_ShowCursor(0);
+	}
+
+	fxwt::set_window_title("The Lab Demos :: Summer Hack");
+	fxwt::set_display_handler(update_gfx);
+	fxwt::set_idle_handler(update_gfx);
+	fxwt::set_keyboard_handler(key_handler);
+	fxwt::set_motion_handler(motion_handler);
+	fxwt::set_button_handler(bn_handler);
+	atexit(clean_up);
+
+	//dsys::init();
+
+	// show loading screen
+	dsys::overlay(get_texture("data/img/loading.jpg"), Vector2(0, 0), Vector2(1, 1), 1.0);
+	flip();
+
+	// check if we have all the textures we need
+	FILE *fp = fopen("data/tex_list", "r");
+	if(!fp) {
+		error("couldn't find texture list, run this only from its directory");
+		return false;
+	}
+
+	char tex_name[512];
+	while(fgets(tex_name, 512, fp)) {
+		char *last_char = tex_name + strlen(tex_name) - 1;
+		if(*last_char == '\n') {
+			*last_char = 0;
+		}
+		if(!get_texture(tex_name)) {
+			error("failed to load texture: \"%s\"", tex_name);
+			return false;
+		}
+	}
+	fclose(fp);
+	
+	set_scene_data_path("data/img");
+
+	INIT_PARTS();
+
+	for(size_t i=0; i<parts.size(); i++) {
+		add_part(parts[i]);
+	}
+
+	if(!render) {
+		dsys::start_demo();
+	} else {
+		dsys::render_demo(render_fps, render_path);
+	}
+
+	if(music) {
+		sdlvf_init("data/music/red_herring.ogg");
+	}
+	
+	timer_reset(&timer);
+	timer_start(&timer);
+	
+	return true;
+}
+
+void clean_up() {
+	unsigned long time = timer_getmsec(&timer);
+	
+	if(music) {
+		sdlvf_done();
+	}
+
+	SDL_ShowCursor(1);
+
+	for(size_t i=0; i<parts.size(); i++) {
+		delete parts[i];
+	}
+	destroy_graphics_context();
+
+	if(frame_count) {
+		avg_frame_time /= frame_count;
+		printf("average frame time: %lu msec, average fps: %.2f\n", avg_frame_time, 1000.0 / avg_frame_time);
+	}
+
+	time /= 1000;
+	cout << "ran for ";
+	if(time < 60) {
+		cout << time << " seconds\n";
+	} else {
+		cout << time / 60 << ":" << time - (time / 60) * 60 << "\n";
+	}
+}
+
+void update_gfx() {
+	unsigned long frame_start = timer_getmsec(&timer);
+	
+	if(music) {
+		int sdlvf_err;
+		if((sdlvf_err = sdlvf_check()) != SDLVF_PLAYING) {
+			fprintf(stderr, "quit: %s", sdlvf_strerror(sdlvf_err));
+			exit(0);
+		}
+	}
+	
+	if(dsys::update_graphics() == -1) {
+		exit(0);
+	}
+
+	avg_frame_time += timer_getmsec(&timer) - frame_start;
+	frame_count++;
+}
diff --git a/tools/curve_draw/Makefile b/tools/curve_draw/Makefile
new file mode 100644
index 0000000..171e2c8
--- /dev/null
+++ b/tools/curve_draw/Makefile
@@ -0,0 +1,13 @@
+obj = curve_draw.o curves.o cursors.o gfx.o
+bin = curve_draw
+
+CXXFLAGS = -g -ansi -pedantic -Wall `sdl-config --cflags`
+
+include vmath/Makefile-part
+
+$(bin): $(obj)
+	$(CXX) -o $@ $(obj) `sdl-config --libs`
+
+.PHONY: clean
+clean:
+	$(RM) $(obj) $(bin)
diff --git a/tools/curve_draw/cursors.cc b/tools/curve_draw/cursors.cc
new file mode 100644
index 0000000..31b58f4
--- /dev/null
+++ b/tools/curve_draw/cursors.cc
@@ -0,0 +1,176 @@
+#include <SDL.h>
+#include "cursors.h"
+
+static SDL_Cursor *create_cursor(const char **image);
+
+SDL_Cursor *cursor_std, *cursor_cross, *cursor_x;
+
+static const char *cur_std[] = {
+  /* width height num_colors chars_per_pixel */
+  "    32    32        3            1",
+  /* colors */
+  "X c #000000",
+  ". c #ffffff",
+  "  c None",
+  /* pixels */
+  ".                               ",
+  "..                              ",
+  ".X.                             ",
+  ".XX.                            ",
+  ".XXX.                           ",
+  ".XXXX.                          ",
+  ".XXXXX.                         ",
+  ".XXXXXX.                        ",
+  ".XXXXXXX.                       ",
+  ".XXXXXXXX.                      ",
+  ".XXXXX.....                     ",
+  ".XX.XX.                         ",
+  ".X. .XX.                        ",
+  "..  .XX.                        ",
+  ".    .XX.                       ",
+  "     .XX.                       ",
+  "      .XX.                      ",
+  "      .XX.                      ",
+  "       ..                       ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "0,0"
+};
+
+static const char *cur_cross[] = {
+  /* width height num_colors chars_per_pixel */
+  "    32    32        3            1",
+  /* colors */
+  "X c #000000",
+  ". c #ffffff",
+  "  c None",
+  /* pixels */
+  "               .                ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "               .                ",
+  " ............     ............  ",
+  ".XXXXXXXXXXXX.   .XXXXXXXXXXXX. ",
+  " ............     ............  ",
+  "               .                ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "              .X.               ",
+  "               .                ",
+  "                                ",
+  "15,15"
+};
+
+static const char *cur_x[] = {
+  /* width height num_colors chars_per_pixel */
+  "    32    32        3            1",
+  /* colors */
+  "X c #000000",
+  ". c #ffffff",
+  "  c None",
+  /* pixels */
+  "                                ",
+  "                                ",
+  "  ..                       ..   ",
+  "  .X.                     .X.   ",
+  "   .X.                   .X.    ",
+  "    .X.                 .X.     ",
+  "     .X.               .X.      ",
+  "      .X.             .X.       ",
+  "       .X.           .X.        ",
+  "        .X.   ...   .X.         ",
+  "         .X. .XXX. .X.          ",
+  "          .X.X...X.X.           ",
+  "           .X.   .X.            ",
+  "          .X.     .X.           ",
+  "         .X.       .X.          ",
+  "         .X.       .X.          ",
+  "         .X.       .X.          ",
+  "          .X.     .X.           ",
+  "           .X.   .X.            ",
+  "          .X.X...X.X.           ",
+  "         .X. .XXX. .X.          ",
+  "        .X.   ...   .X.         ",
+  "       .X.           .X.        ",
+  "      .X.             .X.       ",
+  "     .X.               .X.      ",
+  "    .X.                 .X.     ",
+  "   .X.                   .X.    ",
+  "  .X.                     .X.   ",
+  "  ..                       ..   ",
+  "                                ",
+  "                                ",
+  "                                ",
+  "15,15"
+};
+
+
+
+void init_cursors() {
+	cursor_std = create_cursor(cur_std);
+	cursor_cross = create_cursor(cur_cross);
+	cursor_x = create_cursor(cur_x);
+}
+
+static SDL_Cursor *create_cursor(const char **image) {
+	int i, row, col;
+	Uint8 data[4*32];
+	Uint8 mask[4*32];
+	int hot_x, hot_y;
+
+	i = -1;
+	for(row=0; row<32; ++row) {
+		for (col=0; col<32; ++col) {
+			if(col % 8) {
+				data[i] <<= 1;
+				mask[i] <<= 1;
+			} else {
+				++i;
+				data[i] = mask[i] = 0;
+			}
+			switch(image[4+row][col]) {
+			case 'X':
+				data[i] |= 0x01;
+				mask[i] |= 0x01;
+				break;
+			case '.':
+				mask[i] |= 0x01;
+				break;
+			case ' ':
+				break;
+			}
+		}
+	}
+	sscanf(image[4+row], "%d,%d", &hot_x, &hot_y);
+	return SDL_CreateCursor(data, mask, 32, 32, hot_x, hot_y);
+}
diff --git a/tools/curve_draw/cursors.h b/tools/curve_draw/cursors.h
new file mode 100644
index 0000000..43e163d
--- /dev/null
+++ b/tools/curve_draw/cursors.h
@@ -0,0 +1,8 @@
+#ifndef _CURSORS_H_
+#define _CURSORS_H_
+
+extern SDL_Cursor *cursor_std, *cursor_cross, *cursor_x;
+
+void init_cursors();
+
+#endif	// _CURSORS_H_
diff --git a/tools/curve_draw/curve_draw.cc b/tools/curve_draw/curve_draw.cc
new file mode 100644
index 0000000..80ff0cd
--- /dev/null
+++ b/tools/curve_draw/curve_draw.cc
@@ -0,0 +1,461 @@
+#include <iostream>
+#include <vector>
+#include <SDL.h>
+#include "vmath/vmath.h"
+#include "curves.h"
+#include "cursors.h"
+#include "gfx.h"
+
+using namespace std;
+
+enum {
+	TOOL_SELECT,
+	TOOL_INSERT,
+	TOOL_REMOVE
+};
+
+#define TOOL_COUNT	3
+
+SDL_Cursor *tool_cursors[TOOL_COUNT];
+
+struct EditState {
+	Vector2 view_size;
+	unsigned int tool;
+	std::vector<Curve*> curves;
+	int sel_curve;
+	int sel_pt;
+	int sel_x, sel_y;
+	bool draw_wire;
+	bool draw_pt;
+};
+
+void reset_state(EditState *s = 0);
+void redraw();
+Vector2 world_to_dev(const Vector2 &pt);
+Vector2 dev_to_world(const Vector2 &pt);
+void key_handler(SDLKey key, bool pressed);
+void motion_handler(int x, int y);
+void bn_handler(int x, int y, int bn, bool pressed);
+
+
+int xsz = 800, ysz = 600;
+SDL_Surface *sdl_fb;
+bool done = false;
+EditState state;
+bool view_valid = false;
+
+int main(int argc, char **argv) {
+	const char *usage = "curve_draw [-s(--size) WxH] [curve(s) to load]\n";
+
+	reset_state();
+
+	for(int i=1; i<argc; i++) {
+		if(argv[i][0] == '-') {
+			if(!strcmp(argv[i], "-s") || !strcmp(argv[i], "--size")) {
+				if(!isdigit(argv[++i][0])) {
+					cerr << "invalid size argument: " << argv[i] << endl;
+					return EXIT_FAILURE;
+				}
+				xsz = atoi(argv[i]);
+
+				char *ptr = strchr(argv[i], 'x');
+				if(!ptr || !isdigit(*(ptr + 1))) {
+					cerr << "invalid size argument: " << argv[i] << endl;
+					return EXIT_FAILURE;
+				}
+				ysz = atoi(ptr + 1);
+				
+			} else if(!strcmp(argv[i], "-h") || !strcmp(argv[i], "--help")) {
+				fputs(usage, stdout);
+				return 0;
+			} else {
+				fputs(usage, stderr);
+				return EXIT_FAILURE;
+			}
+		} else {
+			Curve *curve = load_curve(argv[i]);
+			if(!curve) {
+				fprintf(stderr, "failed to load curve: %s\n", argv[i]);
+				return EXIT_FAILURE;
+			}
+			state.curves.push_back(curve);
+			view_valid = false;
+		}
+	}
+	
+	SDL_Init(SDL_INIT_VIDEO);
+	if(!(sdl_fb = SDL_SetVideoMode(xsz, ysz, 32, SDL_SWSURFACE))) {
+		cerr << "failed to set initialize graphics\n";
+		return EXIT_FAILURE;
+	}
+
+	SDL_WM_SetCaption("curvedraw - press h for usage help", 0);
+
+	init_cursors();
+
+	tool_cursors[TOOL_SELECT] = cursor_std;
+	tool_cursors[TOOL_INSERT] = cursor_cross;
+	tool_cursors[TOOL_REMOVE] = cursor_x;
+
+
+	SDL_Event event;
+	while(!done && SDL_WaitEvent(&event)) {
+		do {
+			switch(event.type) {
+			case SDL_QUIT:
+				done = true;
+				break;
+
+			case SDL_KEYDOWN:
+			case SDL_KEYUP:
+				key_handler(event.key.keysym.sym, event.type == SDL_KEYDOWN);
+				break;
+
+			case SDL_MOUSEMOTION:
+				motion_handler(event.motion.x, event.motion.y);
+				break;
+
+			case SDL_MOUSEBUTTONDOWN:
+			case SDL_MOUSEBUTTONUP:
+				bn_handler(event.button.x, event.button.y, event.button.button, event.button.state == SDL_PRESSED);
+				break;
+
+			default:
+				break;
+			}
+		} while(SDL_PollEvent(&event));
+
+		if(!view_valid) redraw();
+	}
+	
+	SDL_Quit();
+	return 0;
+}
+
+
+void reset_state(EditState *s) {
+	if(!s) s = &state;
+
+	s->view_size.y = 100;
+	s->view_size.x = s->view_size.y * ((float)xsz / (float)ysz);
+
+	s->tool = TOOL_SELECT;
+	s->sel_curve = -1;
+	s->sel_pt = -1;
+	s->sel_x = s->sel_y = -1;
+
+	s->draw_wire = false;
+	s->draw_pt = true;
+
+	for(size_t i=0; i<s->curves.size(); i++) {
+		delete s->curves[i];
+	}
+	s->curves.clear();
+}
+
+void redraw() {
+
+	SDL_FillRect(sdl_fb, 0, 0);
+	
+	if(SDL_MUSTLOCK(sdl_fb)) {
+		SDL_LockSurface(sdl_fb);
+	}
+	unsigned int *fb = (unsigned int*)sdl_fb->pixels;
+
+	memset(fb + xsz * ysz / 2, 0x18, xsz * 4);
+
+	unsigned int *ptr = fb + xsz / 2;
+	for(int i=0; i<ysz; i++) {
+		*ptr = 0x181818;
+		ptr += xsz;
+	}
+	
+
+	for(size_t i=0; i<state.curves.size(); i++) {
+		Vector2 prev;
+
+		// draw a line between control points
+		if(state.draw_wire) {
+			for(int j=0; j<state.curves[i]->get_point_count(); j++) {
+				Vector3 *ptptr = state.curves[i]->get_control_point(j);
+				Vector2 pt = world_to_dev(*ptptr);
+
+				if(j > 0) {
+					draw_line((int)prev.x, (int)prev.y, (int)pt.x, (int)pt.y, 0x404040, fb, xsz);
+				}
+				prev = pt;
+			}
+		}
+
+
+		// draw the spline as a polyline
+		const int draw_seg = 12 * state.curves[i]->get_segment_count();
+		for(int j=0; j<=draw_seg; j++) {
+			float t = (float)j / (float)draw_seg;
+			Vector2 pt = (*state.curves[i])(t);
+
+			pt = world_to_dev(pt);
+
+			if(j > 0) {
+				unsigned int col = state.sel_curve == (int)i ? 0x00ff00 : 0x0000ff;
+				draw_line((int)prev.x, (int)prev.y, (int)pt.x, (int)pt.y, col, fb, xsz);
+			}
+			
+			prev = pt;
+		}
+
+		// draw the control points
+		if(state.draw_pt) {
+			for(int j=0; j<state.curves[i]->get_point_count(); j++) {
+				Vector3 *ptptr = state.curves[i]->get_control_point(j);
+				Vector2 pt = world_to_dev(*ptptr);
+
+				unsigned int col = 0x808080;
+				if(state.sel_curve == (int)i) {
+					col = state.sel_pt == j ? 0xffff00 : 0xff0000;
+				}
+			
+				draw_point((int)pt.x, (int)pt.y, col, fb, xsz);
+			}
+		}
+	}
+
+	if(SDL_MUSTLOCK(sdl_fb)) {
+		SDL_UnlockSurface(sdl_fb);
+	}
+
+	SDL_Flip(sdl_fb);
+	view_valid = true;
+}
+
+Vector2 world_to_dev(const Vector2 &pt) {
+	Vector2 norm = (Vector2(1, -1) * pt + state.view_size / 2.0) / state.view_size;
+	return norm * Vector2(xsz, ysz);
+}
+
+Vector2 dev_to_world(const Vector2 &pt) {
+	Vector2 norm = pt / Vector2(xsz, ysz);
+	return Vector2(1, -1) * ((norm * state.view_size) - state.view_size / 2.0);
+}
+
+
+void key_handler(SDLKey key, bool pressed) {
+	switch(key) {
+	case 'q':
+		done = true;
+		break;
+
+	case SDLK_ESCAPE:
+		state.sel_curve = -1;
+		state.sel_pt = -1;
+		view_valid = false;
+		break;
+
+	case '1':
+		redraw();
+		break;
+
+	case 'w':
+		if(!pressed) {
+			state.draw_wire = !state.draw_wire;
+			view_valid = false;
+		}
+		break;
+
+	case 'p':
+		if(!pressed) {
+			state.draw_pt = !state.draw_pt;
+			view_valid = false;
+		}
+		break;
+
+	case SDLK_LEFT:
+		if(!pressed) {
+			if(state.sel_curve <= 0) {
+				state.sel_curve = state.curves.size() - 1;
+			} else {
+				state.sel_curve--;
+			}
+			state.sel_pt = -1;
+			view_valid = false;
+		}
+		break;
+
+	case SDLK_SPACE:
+	case SDLK_RIGHT:
+		if(!pressed) {	
+			if(state.sel_curve == -1 || state.sel_curve >= (int)state.curves.size() - 1) {
+				state.sel_curve = 0;
+			} else {
+				state.sel_curve++;
+			}
+			state.sel_pt = -1;
+			view_valid = false;
+		}
+		break;
+
+	case SDLK_RETURN:
+		if(!pressed && state.sel_curve != -1) {
+			Curve *current = state.curves[state.sel_curve];
+			Curve *curve = dynamic_cast<BSplineCurve*>(current) ? (Curve*)new CatmullRomSplineCurve : (Curve*)new BSplineCurve;
+
+			for(int i=0; i<current->get_point_count(); i++) {
+				Vector3 pt = *current->get_control_point(i);
+				pt.z = 0;
+				curve->add_control_point(pt);
+			}
+
+			delete current;
+			state.curves[state.sel_curve] = curve;
+			view_valid = false;
+		}
+		break;
+
+	case 'h':
+		if(!pressed) {
+			cout << "----- controls -----\n";
+			cout << "[keyboard]\n";
+			cout << "space - cycle through curves\n";
+			cout << "right - cycle through curves forward\n";
+			cout << "left  - cycle through curves backwards\n";
+			cout << "esc   - deselect selected curve\n";
+			cout << "enter - toggle curve type (bspline/catmull-rom)\n";
+			cout << "q     - quit\n";
+			cout << "s     - save all curves to current directory\n";
+			cout << "[mouse]\n";
+			cout << "left-button  - perform action (select tool also drags points around)\n";
+			cout << "right-button - cycle through tools (select/insert/remove)\n";
+			cout << "mousewheel   - zoom in/out\n";
+			cout << endl;
+		}
+		break;
+
+	case 's':
+		if(!pressed) {
+			char fname[128];
+			for(size_t i=0; i<state.curves.size(); i++) {
+				sprintf(fname, "curve%02d", (int)i);
+				save_curve(fname, state.curves[i]);
+			}
+			cout << state.curves.size() << " curves saved.\n";
+		}
+		break;
+
+	default:
+		break;
+	}
+}
+
+void motion_handler(int x, int y) {
+	Vector2 mouse = dev_to_world(Vector2(x, y));
+	
+	if(state.sel_x != -1) {
+		if(x != state.sel_x || y != state.sel_y) {
+			Curve *curve = state.curves[state.sel_curve];
+			Vector3 *ptptr = curve->get_control_point(state.sel_pt);
+			*ptptr = mouse;
+
+			state.sel_x = x;
+			state.sel_y = y;
+			view_valid = false;
+		}
+	}
+}
+
+void bn_handler(int x, int y, int bn, bool pressed) {
+	Vector2 mouse = dev_to_world(Vector2(x, y));
+
+	if(bn == 3 && pressed) {
+		state.tool = (state.tool + 1) % TOOL_COUNT;
+		SDL_SetCursor(tool_cursors[state.tool]);
+	}
+	
+	if(bn == 1 && state.tool == TOOL_SELECT) {
+		if(pressed) {
+			
+			if(state.sel_curve >= 0) {
+				Curve *curve = state.curves[state.sel_curve];
+				for(int i=0; i<curve->get_point_count(); i++) {
+					Vector3 *ptptr = curve->get_control_point(i);
+					Vector2 pt = *ptptr;
+				
+					if((pt - mouse).length() < state.view_size.y * 0.008) {
+						state.sel_pt = i;
+						state.sel_x = x;
+						state.sel_y = y;
+						view_valid = false;
+					}
+				}
+			}
+
+		} else {	// depressed
+			state.sel_x = -1;
+			state.sel_y = -1;
+		}
+	}
+
+	if(bn == 1 && state.tool == TOOL_REMOVE) {
+		static int press_x, press_y;
+		
+		if(pressed) {
+			press_x = x;
+			press_y = y;
+		} else {
+			if(press_x == x && press_y == y && state.sel_curve >= 0) {
+				Curve *curve = state.curves[state.sel_curve];
+				for(int i=0; i<curve->get_point_count(); i++) {
+					Vector3 *ptptr = curve->get_control_point(i);
+					Vector2 pt = *ptptr;
+				
+					if((pt - mouse).length() < state.view_size.y * 0.008) {
+						curve->remove_control_point(i);
+						if(state.sel_pt >= i) state.sel_pt--;
+						if(!curve->get_point_count()) {
+							delete curve;
+							state.curves.erase(state.curves.begin() + state.sel_curve);
+						}
+						view_valid = false;
+						break;
+					}
+				}
+			}
+		}
+	}
+
+	if(bn == 1 && state.tool == TOOL_INSERT) {
+		static int press_x, press_y;
+		
+		if(pressed) {
+			press_x = x;
+			press_y = y;
+		} else {
+			if(press_x == x && press_y == y) {
+				if(state.sel_curve == -1) {
+					// create a new curve...
+					Curve *curve = new CatmullRomSplineCurve;
+					curve->add_control_point(Vector3(mouse.x, mouse.y, 0.0));
+
+					state.sel_curve = state.curves.size();
+					state.sel_pt = 0;
+					state.curves.push_back(curve);
+				} else {
+					// add a control point to the selected curve
+					Curve *curve = state.curves[state.sel_curve];
+					state.sel_pt = curve->get_point_count();
+					curve->add_control_point(Vector3(mouse.x, mouse.y, 0.0));
+				}
+				view_valid = false;
+			}
+		}
+	}
+
+	if(bn == 4 && pressed) {
+		state.view_size *= 0.95;
+		view_valid = false;
+	}
+
+	if(bn == 5 && pressed) {
+		state.view_size *= 1.05;
+		view_valid = false;
+	}
+}
diff --git a/tools/curve_draw/curves.cc b/tools/curve_draw/curves.cc
new file mode 100644
index 0000000..b8dc659
--- /dev/null
+++ b/tools/curve_draw/curves.cc
@@ -0,0 +1,488 @@
+/*
+Copyright 2004 John Tsiombikas <nuclear@siggraph.org>
+
+This file is part of the graphics core library.
+
+the graphics core library is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+the graphics core library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the graphics core library; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*/
+
+/* 3D Curves
+ * 
+ * author: John Tsiombikas 2003
+ * modified:
+ * 		John Tsiombikas 2004
+ * 		Mihalis Georgoulopoulos 2004
+ */
+
+#include <cstdio>
+#include <cstring>
+#include <cmath>
+#include <cctype>
+#include <cassert>
+#include "curves.h"
+//#include "common/err_msg.h"
+
+Curve::Curve() {
+	arc_parametrize = false;
+	ease_curve = 0;
+	samples = 0;
+
+	set_ease_sample_count(100);
+}
+
+Curve::~Curve() {
+	delete [] samples;
+	
+}
+
+void Curve::set_arc_parametrization(bool state) {
+	arc_parametrize = state;
+}
+
+#define Param	0
+#define ArcLen	1
+
+void Curve::sample_arc_lengths() {
+	const int samplesPerSegment = 30;
+	sample_count = get_segment_count() * samplesPerSegment;
+
+	arc_parametrize = false;	// to be able to interpolate with the original values
+
+	samples = new Vector2[sample_count];
+	Vector3 prevpos;
+	scalar_t step = 1.0f / (scalar_t)(sample_count-1);
+	for(int i=0; i<sample_count; i++) {
+		scalar_t t = step * (scalar_t)i;
+		Vector3 pos = interpolate(t);
+		samples[i][Param] = t;
+		if(!i) {
+			samples[i][ArcLen] = 0.0f;
+		} else {
+			samples[i][ArcLen] = (pos - prevpos).length() + samples[i-1][ArcLen];
+		}
+		prevpos = pos;
+	}
+
+	// normalize arc lenghts
+	scalar_t maxlen = samples[sample_count-1][ArcLen];
+	for(int i=0; i<sample_count; i++) {
+		samples[i][ArcLen] /= maxlen;
+	}
+
+	arc_parametrize = true;
+}
+
+static int binary_search(Vector2 *array, scalar_t key, int begin, int end) {
+	int middle = begin + ((end - begin)>>1);
+
+	if(array[middle][ArcLen] == key) return middle;
+	if(end == begin) return middle;
+
+	if(key < array[middle][ArcLen]) return binary_search(array, key, begin, middle);
+	if(key > array[middle][ArcLen]) return binary_search(array, key, middle+1, end);
+	return -1;	// just to make the compiler shut the fuck up
+}
+
+scalar_t Curve::parametrize(scalar_t t) const {
+	if(!samples) const_cast<Curve*>(this)->sample_arc_lengths();
+
+	int samplepos = binary_search(samples, t, 0, sample_count);
+	scalar_t par = samples[samplepos][Param];
+	scalar_t len = samples[samplepos][ArcLen];
+
+	// XXX: I can't remember the significance of this condition, I had xsmall_number here
+	// previously and if t was 0.9999 it broke. I just changed the number blindly which fixed
+	// the breakage but I should investigate further at some point.
+	if((len - t) < 0.0005) return par;
+
+	if(len < t) {
+		if(!samplepos) return par;
+		scalar_t prevlen = samples[samplepos-1][ArcLen];
+		scalar_t prevpar = samples[samplepos-1][Param];
+		scalar_t p = (t - prevlen) / (len - prevlen);
+		return prevpar + (par - prevpar) * p;
+	} else {
+		if(samplepos >= sample_count) return par;
+		scalar_t nextlen = samples[samplepos+1][ArcLen];
+		scalar_t nextpar = samples[samplepos+1][Param];
+		scalar_t p = (t - len) / (nextlen - len);
+		return par + (nextpar - par) * p;
+	}
+
+	return par;	// not gonna happen
+}
+
+
+#define MIN(a, b) ((a) < (b) ? (a) : (b))
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+
+scalar_t Curve::ease(scalar_t t) const {
+	if(!ease_curve) return t;
+
+	const_cast<Curve*>(this)->ease_curve->set_arc_parametrization(true);
+	scalar_t et = ease_curve->interpolate(t).y;
+
+	return MIN(MAX(et, 0.0f), 1.0f);
+}
+
+
+void Curve::add_control_point(const Vector3 &cp) {
+	control_points.push_back(cp);
+	delete [] samples;
+	samples = 0;
+}
+
+void Curve::remove_control_point(int index) {
+	if(index < 0 || index >= control_points.size()) return;
+
+	ListNode<Vector3> *node = control_points.begin();
+	for(int i=0; i<index; i++) {
+		node = node->next;
+	}
+
+	control_points.erase(node);
+}
+
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+#define MAX(a, b)	((a) > (b) ? (a) : (b))
+
+Vector3 *Curve::get_control_point(int index) {
+	index = MAX(0, MIN(index, control_points.size() - 1));
+
+	ListNode<Vector3> *node = control_points.begin();
+	for(int i=0; i<index; i++) {
+		assert(node);
+		node = node->next;
+	}
+	
+	assert(node);
+	return &node->data;
+}
+
+int Curve::get_point_count() const {
+	return control_points.size();
+}
+
+void Curve::set_ease_curve(Curve *curve) {
+	ease_curve = curve;
+}
+
+void Curve::set_ease_sample_count(int count) {
+	ease_sample_count = count;
+	ease_step = (int)(1.0f / (scalar_t)ease_sample_count);
+}
+
+Vector3 Curve::operator ()(scalar_t t) const {
+	return interpolate(t);
+}
+
+///////////////// B-Spline implementation ////////////////////
+
+int BSplineCurve::get_segment_count() const {
+	return control_points.size() - 3;
+}
+
+Vector3 BSplineCurve::interpolate(scalar_t t) const {
+	if(t > 1.0) t = 1.0;
+	if(t < 0.0) t = 0.0;
+
+	if(control_points.size() < 4) return Vector3(0, 0, 0);
+
+	if(arc_parametrize) {
+		t = ease(parametrize(t));
+	}
+
+	// find the appropriate segment of the spline that t lies and calculate the piecewise parameter
+	t = (scalar_t)(control_points.size() - 3) * t;
+	int seg = (int)t;
+	t -= (scalar_t)floor(t);
+	if(seg >= get_segment_count()) {
+		seg = get_segment_count() - 1;
+		t = 1.0f;
+	}
+	
+	ListNode<Vector3> *iter = const_cast<BSplineCurve*>(this)->control_points.begin();
+	for(int i=0; i<seg; i++) iter = iter->next;
+
+	Vector3 Cp[4];
+	for(int i=0; i<4; i++) {
+        Cp[i] = iter->data;
+		iter = iter->next;
+	}
+
+	Vector3 res;
+	res.x = bspline(Cp[0].x, Cp[1].x, Cp[2].x, Cp[3].x, t);
+	res.y = bspline(Cp[0].y, Cp[1].y, Cp[2].y, Cp[3].y, t);
+	res.z = bspline(Cp[0].z, Cp[1].z, Cp[2].z, Cp[3].z, t);
+
+	return res;
+}
+
+//////////////// Catmull-Rom Spline implementation //////////////////
+
+int CatmullRomSplineCurve::get_segment_count() const {
+	return control_points.size() - 1;
+}
+
+Vector3 CatmullRomSplineCurve::interpolate(scalar_t t) const {
+	if(t > 1.0) t = 1.0;
+	if(t < 0.0) t = 0.0;
+
+	if(control_points.size() < 2) return Vector3(0, 0, 0);
+
+	if(arc_parametrize) {
+		t = ease(parametrize(t));
+	}
+
+	// find the appropriate segment of the spline that t lies and calculate the piecewise parameter
+	t = (scalar_t)(control_points.size() - 1) * t;
+	int seg = (int)t;
+	t -= (scalar_t)floor(t);
+	if(seg >= get_segment_count()) {
+		seg = get_segment_count() - 1;
+		t = 1.0f;
+	}
+
+	Vector3 Cp[4];
+	ListNode<Vector3> *iter = const_cast<CatmullRomSplineCurve*>(this)->control_points.begin();
+	for(int i=0; i<seg; i++) iter = iter->next;
+
+	Cp[1] = iter->data;
+	Cp[2] = iter->next->data;
+	
+	if(!seg) {
+		Cp[0] = Cp[1];
+	} else {
+		Cp[0] = iter->prev->data;
+	}
+	
+	if(seg == control_points.size() - 2) {
+		Cp[3] = Cp[2];
+	} else {
+		Cp[3] = iter->next->next->data;
+	}
+
+	Vector3 res;
+	res.x = catmull_rom_spline(Cp[0].x, Cp[1].x, Cp[2].x, Cp[3].x, t);
+	res.y = catmull_rom_spline(Cp[0].y, Cp[1].y, Cp[2].y, Cp[3].y, t);
+	res.z = catmull_rom_spline(Cp[0].z, Cp[1].z, Cp[2].z, Cp[3].z, t);
+
+	return res;
+}
+
+/* BezierSpline implementation - (MG)
+ */
+int BezierSpline::get_segment_count() const
+{
+	if (control_points.size() < 0) return 0;
+	return control_points.size() / 4;
+}
+
+Vector3 BezierSpline::interpolate(scalar_t t) const
+{
+	if (!get_segment_count()) return Vector3(0, 0, 0);
+
+	if (arc_parametrize)
+	{
+		t = ease(parametrize(t));
+	}
+ 
+	t = (t * get_segment_count());
+	int seg = (int) t;
+	t -= (scalar_t) floor(t);
+	if (seg >= get_segment_count())
+	{
+		seg = get_segment_count() - 1;
+		t = 1.0f;
+	}
+
+	seg *= 4;
+	ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
+	for (int i = 0; i < seg; i++) iter = iter->next;
+	
+	Vector3 Cp[4];
+	for (int i = 0; i < 4; i++)
+	{
+		Cp[i] = iter->data;
+		iter = iter->next;
+	}
+	
+	// interpolate
+	return bezier(Cp[seg], Cp[seg + 1], Cp[seg + 2], Cp[seg + 3], t);
+}
+
+Vector3 BezierSpline::get_tangent(scalar_t t)
+{	
+	if (!get_segment_count()) return Vector3(0, 0, 0);
+
+	if (arc_parametrize)
+	{
+		t = ease(parametrize(t));
+	}
+ 
+	t = (t * get_segment_count());
+	int seg = (int) t;
+	t -= (scalar_t) floor(t);
+	if (seg >= get_segment_count())
+	{
+		seg = get_segment_count() - 1;
+		t = 1.0f;
+	}
+
+	seg *= 4;
+	ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
+	for (int i = 0; i < seg; i++) iter = iter->next;
+	
+	Vector3 Cp[4];
+	for (int i = 0; i < 4; i++)
+	{
+		Cp[i] = iter->data;
+		iter = iter->next;
+	}
+	
+	// interpolate
+	return bezier_tangent(Cp[0], Cp[1], Cp[2], Cp[3], t);
+}
+
+Vector3 BezierSpline::get_control_point(int i) const
+{	
+	
+	ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
+	for (int j = 0; j < i; j++) 
+	{
+		if (!iter->next)
+		{
+			return Vector3(0, 0, 0);
+		}
+		iter = iter->next;
+	}
+	
+	return iter->data;
+}
+
+
+
+bool save_curve(const char *fname, const Curve *curve) {
+	FILE *fp = fopen(fname, "w");
+	if(!fp) {
+		fprintf(stderr, "failed to save the curve %s", curve->name.c_str());
+		return false;
+	}
+
+	fputs("curve_3dengfx\n", fp);
+	fputs(curve->name.c_str(), fp);
+	fputs("\n", fp);
+
+	if(dynamic_cast<const BSplineCurve*>(curve)) {
+		fputs("bspline\n", fp);
+	} else if(dynamic_cast<const CatmullRomSplineCurve*>(curve)) {
+		fputs("catmullrom\n", fp);
+	} else if(dynamic_cast<const BezierSpline*>(curve)) {
+		fputs("bezier\n", fp);
+	} else {
+		fprintf(stderr, "unknown spline type, save failed %s", curve->name.c_str());
+		fclose(fp);
+		remove(fname);
+		return false;
+	}
+
+	fprintf(fp, "%d\n", curve->control_points.size());
+
+	const ListNode<Vector3> *node = curve->control_points.begin();
+	while(node) {
+		fprintf(fp, "%f %f %f\n", node->data.x, node->data.y, node->data.z);
+		node = node->next;
+	}
+
+	fclose(fp);
+	return true;
+}
+
+Curve *load_curve(const char *fname) {
+	FILE *fp = fopen(fname, "r");
+	if(!fp) {
+		fprintf(stderr, "failed to open file %s", fname);
+		return 0;
+	}
+
+	char buffer[256];
+
+	fgets(buffer, 256, fp);
+	if(strcmp(buffer, "curve_3dengfx\n") != 0) {
+		fprintf(stderr, "load_curve failed, %s is not a curve file", fname);
+		fclose(fp);
+		return 0;
+	}
+
+	Curve *curve;
+
+	fgets(buffer, 256, fp);
+	std::string name = buffer;
+
+	fgets(buffer, 256, fp);
+	if(!strcmp(buffer, "bspline\n")) {
+		curve = new BSplineCurve;
+	} else if(!strcmp(buffer, "catmullrom\n")) {
+		curve = new CatmullRomSplineCurve;
+	} else /*if(!strcmp(buffer, "bezier"))*/ {
+		fprintf(stderr, "unsupported curve type (%s) or not a curve file", buffer);
+		fclose(fp);
+		return 0;
+	}
+
+	curve->name = name;
+
+	fgets(buffer, 256, fp);
+	if(!isdigit(buffer[0])) {
+		fprintf(stderr, "load_curve failed, %s is not a valid curve file (count: %s)", fname, buffer);
+		delete curve;
+		fclose(fp);
+		return 0;
+	}
+	int count = atoi(buffer);
+
+	int failed = count;
+	for(int i=0; i<count; i++, failed--) {
+		fgets(buffer, 256, fp);
+		if(!isdigit(buffer[0]) && buffer[0] != '.' && buffer[0] != '-') {
+			break;
+		}
+		float x = atof(buffer);
+
+		char *ptr = strchr(buffer, ' ');
+		if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
+			break;
+		}
+		float y = atof(++ptr);
+		
+		ptr = strchr(ptr, ' ');
+		if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
+			break;
+		}
+		float z = atof(++ptr);
+
+		curve->add_control_point(Vector3(x, y, z));
+	}
+
+	fclose(fp);
+
+	if(failed) {
+		fprintf(stderr, "load_curve failed to read the data, %s is not a valid curve file", fname);
+		delete curve;
+		return 0;
+	}
+
+	return curve;
+}
diff --git a/tools/curve_draw/curves.h b/tools/curve_draw/curves.h
new file mode 100644
index 0000000..08769ec
--- /dev/null
+++ b/tools/curve_draw/curves.h
@@ -0,0 +1,96 @@
+/*
+This file is part of the graphics core library.
+
+Copyright (c) 2004, 2005 John Tsiombikas <nuclear@siggraph.org>
+
+the graphics core library is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+the graphics core library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the graphics core library; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*/
+
+/* 3D Curves
+ * 
+ * author: John Tsiombikas 2003
+ * modified:
+ * 		John Tsiombikas 2004, 2005
+ * 		Mihalis Georgoulopoulos 2004
+ */
+
+#ifndef _CURVES_HPP_
+#define _CURVES_HPP_
+
+#include <string>
+#include "vmath/vmath.h"
+#include "linkedlist.h"
+
+
+class Curve {
+protected:
+	LinkedList<Vector3> control_points;
+	Vector2 *samples;	// used for parametrizing by arc length
+	int sample_count;
+	bool arc_parametrize;
+
+	Curve *ease_curve;	// ease in/out curve (1D, x&z discarded)
+	int ease_sample_count, ease_step;
+	
+	void sample_arc_lengths();
+	scalar_t parametrize(scalar_t t) const;
+	scalar_t ease(scalar_t t) const;
+
+public:
+	std::string name;
+
+	Curve();
+	virtual ~Curve();
+	virtual void add_control_point(const Vector3 &cp);
+	virtual void remove_control_point(int index);
+	virtual Vector3 *get_control_point(int index);
+
+	virtual int get_point_count() const;
+	virtual int get_segment_count() const = 0;
+	virtual void set_arc_parametrization(bool state);
+	virtual void set_ease_curve(Curve *curve);
+	virtual void set_ease_sample_count(int count);
+
+	virtual Vector3 interpolate(scalar_t t) const = 0;
+	virtual Vector3 operator ()(scalar_t t) const;
+
+	friend bool save_curve(const char *fname, const Curve *curve);
+};
+
+class BSplineCurve : public Curve {
+public:
+	virtual int get_segment_count() const;
+	virtual Vector3 interpolate(scalar_t t) const;
+};
+
+class CatmullRomSplineCurve : public Curve {
+public:
+	virtual int get_segment_count() const;
+	virtual Vector3 interpolate(scalar_t t) const;
+};
+
+class BezierSpline : public Curve {
+public:
+	virtual int get_segment_count() const;
+	virtual Vector3 interpolate(scalar_t t) const;
+
+	Vector3 get_control_point(int i) const;
+	Vector3 get_tangent(scalar_t t);
+};
+
+bool save_curve(const char *fname, const Curve *curve);
+Curve *load_curve(const char *fname);
+
+#endif	// _CURVES_HPP_
diff --git a/tools/curve_draw/gfx.cc b/tools/curve_draw/gfx.cc
new file mode 100644
index 0000000..2332f55
--- /dev/null
+++ b/tools/curve_draw/gfx.cc
@@ -0,0 +1,299 @@
+#include "gfx.h"
+
+#define RGB(r, g, b)	(((r) << 16) & 0xff0000) | (((g) << 8) & 0x00ff00) | ((b) & 0xff)
+
+static int clip_line(int *x1, int *y1, int *x2, int *y2);
+
+
+extern int xsz, ysz;
+	
+
+void draw_point(int x, int y, unsigned int col, unsigned int *ptr, int xsz) {
+	const int sz = 9;
+	const int pt_img[sz][sz] = {
+		{0,0,0,0,1,0,0,0,0},
+		{0,0,1,1,1,1,1,0,0},
+		{0,1,2,2,1,1,1,3,0},
+		{0,1,2,1,1,1,1,3,0},
+		{1,1,1,1,1,1,1,3,3},
+		{0,1,1,1,1,1,3,3,0},
+		{0,1,1,1,1,3,3,3,0},
+		{0,0,3,3,3,3,3,0,0},
+		{0,0,0,0,3,0,0,0,0}
+	};
+
+	if(x < sz/2 || x >= xsz - sz/2 || y < sz/2 || y >= ysz - sz/2) {
+		return;
+	}
+
+	int base_r = (col & 0xff0000) >> 16;
+	int base_g = (col & 0x00ff00) >> 8;
+	int base_b = (col & 0x0000ff);
+
+	unsigned int dark = RGB(base_r - base_r/4, base_g - base_g/4, base_b - base_b/4);;
+	
+	unsigned int ctable[] = {0, col, 0xffffff, dark};
+	ptr += (y - sz / 2) * xsz + (x - sz / 2);
+
+	for(int i=0; i<sz; i++) {
+		for(int j=0; j<sz; j++) {
+			if(pt_img[i][j]) {
+				*(ptr + j) = ctable[pt_img[i][j]];
+			}
+		}
+		ptr += xsz;
+	}
+}
+
+void draw_line(int x1, int y1, int x2, int y2, unsigned int col, unsigned int *ptr, int xsz) {
+	int dx, dy, dx2, dy2;
+	int x_inc, y_inc;
+	int error;
+	int i;
+
+	if(!clip_line(&x1, &y1, &x2, &y2)) return;
+
+	ptr += y1 * xsz + x1;
+	dx = x2 - x1;
+	dy = y2 - y1;
+	
+	if(dx >= 0) {
+		x_inc = 1;
+	} else {
+		x_inc = -1;
+		dx = -dx;
+	}
+	
+	if(dy >= 0) {
+		y_inc = xsz;
+	} else {
+		y_inc = -xsz;
+		dy = -dy;
+	}
+	
+	dx2 = dx << 1;
+	dy2 = dy << 1;
+
+	if(dx > dy) {
+		error = dy2 - dx;
+		for(i=0; i<=dx; i++) {
+			*ptr = col;
+			if(error >= 0) {
+				error -= dx2;
+				ptr += y_inc;
+			}
+			error += dy2;
+			ptr += x_inc;
+		}
+	} else {
+		error = dx2 - dy;
+		for(i=0;i<=dy;i++) {
+			*ptr = col;
+			if(error >= 0) {
+				error -= dy2;
+				ptr += x_inc;
+			}
+			error += dx2;
+			ptr += y_inc;
+		}
+	}
+}
+
+enum {
+	CLIP_C		= 0x0000,
+	CLIP_N		= 0x0008,
+	CLIP_S		= 0x0004,
+	CLIP_E		= 0x0002,
+	CLIP_W		= 0x0001,
+	CLIP_NE		= 0x000a,
+	CLIP_SE		= 0x0006,
+	CLIP_NW		= 0x0009,
+	CLIP_SW		= 0x0005
+};
+
+#define MIN_CLIP_X	0
+#define MIN_CLIP_Y	0
+#define MAX_CLIP_X	(xsz - 1)
+#define MAX_CLIP_Y	(ysz - 1)
+
+static int clip_line(int *x1, int *y1, int *x2, int *y2) {
+	int xc1 = *x1;
+	int yc1 = *y1;
+	int xc2 = *x2;
+	int yc2 = *y2;
+	int p1_code = 0;
+	int p2_code = 0;
+
+	// determine codes for p1 and p2
+	if(*y1 < MIN_CLIP_Y) {
+		p1_code |= CLIP_N;
+	} else {
+		if(*y1 > MAX_CLIP_Y) {
+			p1_code |= CLIP_S;
+		} else {
+			if(*x1 < MIN_CLIP_X) {
+				p1_code |= CLIP_W;
+			} else {
+				if(*x1 > MAX_CLIP_X) p1_code |= CLIP_E;
+			}
+		}
+	}
+
+	if(*y2 < MIN_CLIP_Y) {
+		p2_code |= CLIP_N;
+	} else {
+		if(*y2 > MAX_CLIP_Y) {
+			p2_code |= CLIP_S;
+		} else {
+			if(*x2 < MIN_CLIP_X) {
+				p2_code |= CLIP_W;
+			} else {
+				if(*x2 > MAX_CLIP_X) p2_code |= CLIP_E;
+			}
+		}
+	}
+
+	if((p1_code & p2_code)) return 0; // trivial rejection
+	if(p1_code == 0 && p2_code == 0) return 1; // test if totally visible
+
+	// find clip point for p1
+	switch(p1_code) {
+		case CLIP_C: break;
+		
+		case CLIP_N:
+			yc1 = MIN_CLIP_Y;
+			xc1 = *x1 + 0.5 + (MIN_CLIP_Y - *y1) * (*x2 - *x1) / (*y2 - *y1);
+			break;
+
+		case CLIP_S:
+			yc1 = MAX_CLIP_Y;
+			xc1 = *x1 + 0.5 + (MAX_CLIP_Y - *y1) * (*x2 - *x1) / (*y2 - *y1);
+			break;
+
+		case CLIP_W:
+			xc1 = MIN_CLIP_X;
+			yc1 = *y1 + 0.5 + (MIN_CLIP_X - *x1) * (*y2 - *y1) / (*x2 - *x1);
+			break;
+
+		case CLIP_E:
+			xc1 = MAX_CLIP_X;
+			yc1 = *y1 + 0.5 + (MAX_CLIP_X - *x1) * (*y2 - *y1) / (*x2 - *x1);
+			break;
+
+		case CLIP_NE:
+			yc1 = MIN_CLIP_Y;
+			xc1 = *x1 + 0.5 + (MIN_CLIP_Y - *y1) * (*x2 - *x1) / (*y2 - *y1);
+			if(xc1 < MIN_CLIP_X || xc1 > MAX_CLIP_X) {
+				xc1 = MAX_CLIP_X;
+				yc1 = *y1 + 0.5 + (MAX_CLIP_X - *x1) * (*y2 - *y1) / (*x2 - *x1);
+			}
+			break;
+
+		case CLIP_SE:
+			yc1 = MAX_CLIP_Y;
+			xc1 = *x1 + 0.5 + (MAX_CLIP_Y - *y1) * (*x2 - *x1) / (*y2 - *y1);
+			if(xc1 < MIN_CLIP_X || xc1 > MAX_CLIP_X) {
+				xc1 = MAX_CLIP_X;
+				yc1 = *y1 + 0.5 + (MAX_CLIP_X - *x1) * (*y2 - *y1) / (*x2 - *x1);
+			}
+			break;
+
+		case CLIP_NW:
+			yc1 = MIN_CLIP_Y;
+			xc1 = *x1 + 0.5 + (MIN_CLIP_Y - *y1) * (*x2 - *x1) / (*y2 - *y1);
+			if(xc1 < MIN_CLIP_X || xc1 > MAX_CLIP_X) {
+				xc1 = MIN_CLIP_X;
+				yc1 = *y1 + 0.5 + (MIN_CLIP_X - *x1) * (*y2 - *y1) / (*x2 - *x1);
+			}
+			break;
+
+		case CLIP_SW:
+			yc1 = MAX_CLIP_Y;
+			xc1 = *x1 + 0.5 + (MAX_CLIP_Y - *y1) * (*x2 - *x1) / (*y2 - *y1);
+			if(xc1 < MIN_CLIP_X || xc1 > MAX_CLIP_X) {
+				xc1 = MIN_CLIP_X;
+				yc1 = *y1 + 0.5 + (MIN_CLIP_X - *x1) * (*y2 - *y1) / (*x2 - *x1);
+			}
+			break;
+
+		default: break;
+	}
+
+	// find clip point for p2
+	switch(p2_code) {
+		case CLIP_C: break;
+		
+		case CLIP_N:
+			yc2 = MIN_CLIP_Y;
+			xc2 = *x2 + 0.5 + (MIN_CLIP_Y - *y2) * (*x1 - *x2) / (*y1 - *y2);
+			break;
+
+		case CLIP_S:
+			yc2 = MAX_CLIP_Y;
+			xc2 = *x2 + 0.5 + (MAX_CLIP_Y - *y2) * (*x1 - *x2) / (*y1 - *y2);
+			break;
+
+		case CLIP_W:
+			xc2 = MIN_CLIP_X;
+			yc2 = *y2 + 0.5 + (MIN_CLIP_X - *x2) * (*y1 - *y2) / (*x1 - *x2);
+			break;
+
+		case CLIP_E:
+			xc2 = MAX_CLIP_X;
+			yc2 = *y1 + 0.5 + (MAX_CLIP_X - *x2) * (*y1 - *y2) / (*x1 - *x2);
+			break;
+
+		case CLIP_NE:
+			yc2 = MIN_CLIP_Y;
+			xc2 = *x2 + 0.5 + (MIN_CLIP_Y - *y2) * (*x1 - *x2) / (*y1 - *y2);
+			if(xc2 < MIN_CLIP_X || xc2 > MAX_CLIP_X) {
+				xc2 = MAX_CLIP_X;
+				yc2 = *y2 + 0.5 + (MAX_CLIP_X - *x2) * (*y1 - *y2) / (*x1 - *x2);
+			}
+			break;
+
+		case CLIP_SE:
+			yc2 = MAX_CLIP_Y;
+			xc2 = *x2 + 0.5 + (MAX_CLIP_Y - *y2) * (*x1 - *x2) / (*y1 - *y2);
+			if(xc2 < MIN_CLIP_X || xc2 > MAX_CLIP_X) {
+				xc2 = MAX_CLIP_X;
+				yc2 = *y2 + 0.5 + (MAX_CLIP_X - *x2) * (*y1 - *y2) / (*x1 - *x2);
+			}
+			break;
+
+		case CLIP_NW:
+			yc2 = MIN_CLIP_Y;
+			xc2 = *x2 + 0.5 + (MIN_CLIP_Y - *y2) * (*x1 - *x2) / (*y1 - *y2);
+			if(xc2 < MIN_CLIP_X || xc2 > MAX_CLIP_X) {
+				xc2 = MIN_CLIP_X;
+				yc2 = *y2 + 0.5 + (MIN_CLIP_X - *x2) * (*y1 - *y2) / (*x1 - *x2);
+			}
+			break;
+
+		case CLIP_SW:
+			yc2 = MAX_CLIP_Y;
+			xc2 = *x2 + 0.5 + (MAX_CLIP_Y - *y2) * (*x1 - *x2) / (*y1 - *y2);
+			if(xc2 < MIN_CLIP_X || xc2 > MAX_CLIP_X) {
+				xc2 = MIN_CLIP_X;
+				yc2 = *y2 + 0.5 + (MIN_CLIP_X - *x2) * (*y1 - *y2) / (*x1 - *x2);
+			}
+			break;
+
+		default: break;
+	}
+
+	// bounds check
+	if(	(xc1 < MIN_CLIP_X) || (xc1 > MAX_CLIP_X) ||
+		(yc1 < MIN_CLIP_Y) || (yc1 > MAX_CLIP_Y) ||
+		(xc2 < MIN_CLIP_X) || (xc2 > MAX_CLIP_X) ||
+		(yc2 < MIN_CLIP_Y) || (yc2 > MAX_CLIP_X)  ) return 0;
+
+	*x1 = xc1;
+	*y1 = yc1;
+	*x2 = xc2;
+	*y2 = yc2;
+
+	return 1;
+}
+
+
diff --git a/tools/curve_draw/gfx.h b/tools/curve_draw/gfx.h
new file mode 100644
index 0000000..3c7668f
--- /dev/null
+++ b/tools/curve_draw/gfx.h
@@ -0,0 +1,7 @@
+#ifndef _GFX_H_
+#define _GFX_H_
+
+void draw_point(int x, int y, unsigned int col, unsigned int *ptr, int xsz);
+void draw_line(int x1, int y1, int x2, int y2, unsigned int col, unsigned int *ptr, int xsz);
+
+#endif	// _GFX_H_
diff --git a/tools/curve_draw/linkedlist.h b/tools/curve_draw/linkedlist.h
new file mode 100644
index 0000000..4cb380b
--- /dev/null
+++ b/tools/curve_draw/linkedlist.h
@@ -0,0 +1,260 @@
+/*
+Copyright 2004 John Tsiombikas <nuclear@siggraph.org>
+
+This file is part of the eternal demo.
+
+The eternal library is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+The eternal demo is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the eternal demo; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*/
+#ifndef _LINKEDLIST_H_
+#define _LINKEDLIST_H_
+
+template <class T>
+struct ListNode {
+	T data;
+	ListNode<T> *next, *prev;
+};
+
+
+template <class T>
+class LinkedList {
+private:
+	ListNode<T> *head, *tail;
+	int sz;
+
+public:
+
+	LinkedList();
+	~LinkedList();
+
+	inline ListNode<T> *begin();
+	inline const ListNode<T> *begin() const;
+	inline ListNode<T> *end();
+	inline const ListNode<T> *end() const;
+
+	void push_back(ListNode<T> *node);
+	void push_back(T data);
+
+	void insert(ListNode<T> *pos, ListNode<T> *node);
+	void insert(ListNode<T> *pos, T data);
+
+	void remove(ListNode<T> *node);
+	ListNode<T> *erase(ListNode<T> *node);
+
+	ListNode<T> *find(T key);
+
+	int count_nodes();
+	inline int size() const;
+
+	void operator =(LinkedList &rhs);
+};
+
+
+/////////// implementation //////////
+template <class T>
+LinkedList<T>::LinkedList() {
+	head = tail = 0;
+	sz = 0;
+}
+
+
+template <class T>
+LinkedList<T>::~LinkedList() {
+
+	while(head) {
+		erase(head);
+	}
+}
+
+template <class T>
+ListNode<T> *LinkedList<T>::begin() {
+	return head;
+}
+
+template <class T>
+const ListNode<T> *LinkedList<T>::begin() const {
+	return head;
+}
+
+template <class T>
+ListNode<T> *LinkedList<T>::end() {
+	return tail;
+}
+
+template <class T>
+const ListNode<T> *LinkedList<T>::end() const {
+	return tail;
+}
+
+template <class T>
+void LinkedList<T>::push_back(ListNode<T> *node) {
+
+	if(!head) {		// empty list
+		head = tail = node;
+		node->next = node->prev = 0;
+	} else {
+		tail->next = node;
+		node->prev = tail;
+		tail = node;
+		node->next = 0;
+	}
+
+	sz++;
+}
+
+template <class T>
+void LinkedList<T>::push_back(T data) {
+
+	ListNode<T> *node = new ListNode<T>;
+	node->data = data;
+
+	if(!head) {		// empty list
+		head = tail = node;
+		node->next = node->prev = 0;
+	} else {
+		tail->next = node;
+		node->prev = tail;
+		tail = node;
+		node->next = 0;
+	}
+
+	sz++;
+}
+
+template <class T>
+void LinkedList<T>::insert(ListNode<T> *pos, ListNode<T> *node) {
+
+	if(!head) {
+		head = tail = node;
+		node->next = node->prev = 0;
+	} else {
+		node->prev = pos->prev;
+		pos->prev = node;
+		node->next = pos;
+		if(head == pos) head = node; else node->prev->next = node;
+	}
+
+	sz++;
+}
+
+template <class T>
+void LinkedList<T>::insert(ListNode<T> *pos, T data) {
+
+	ListNode<T> *node = new ListNode<T>;
+	node->data = data;
+
+	if(!head) {
+		head = tail = node;
+		node->next = node->prev = 0;
+	} else {
+		node->prev = pos->prev;
+		pos->prev = node;
+		node->next = pos;
+		if(head == pos) head = node; else node->prev->next = node;
+	}
+
+	sz++;
+}
+
+template <class T>
+void LinkedList<T>::remove(ListNode<T> *node) {
+
+	if(!node) return;	// e.g. remove(head) on an empty list
+
+	if(node->prev) {
+		node->prev->next = node->next;
+	} else {
+		head = node->next;
+	}
+
+	if(node->next) {
+		node->next->prev = node->prev;
+	} else {
+		tail = node->prev;
+	}
+
+	sz--;
+}
+
+template <class T>
+ListNode<T> *LinkedList<T>::erase(ListNode<T> *node) {
+
+	if(!node) return 0;	// e.g. erase(head) on an empty list
+
+	if(node->prev) {
+		node->prev->next = node->next;
+	} else {
+		head = node->next;
+	}
+
+	if(node->next) {
+		node->next->prev = node->prev;
+	} else {
+		tail = node->prev;
+	}
+
+	ListNode<T> *destr = node;
+	node = node->next;
+	
+	delete destr;
+
+	sz--;
+
+	return node;
+}
+
+template <class T>
+ListNode<T> *LinkedList<T>::find(T key) {
+	
+	ListNode<T> *iter = head;
+	while(iter) {
+		if(iter->data == key) return iter;
+		iter = iter->next;
+	}
+
+	return 0;	// null pointer if not found
+}
+
+template <class T>
+int LinkedList<T>::count_nodes() {
+
+	sz = 0;
+
+	ListNode<T> *iter = head;
+	while(iter) {
+		sz++;
+		iter = iter->next;
+	}
+
+	return sz;
+}
+
+template <class T>
+int LinkedList<T>::size() const {
+	return sz;
+}
+
+
+template <class T>
+void LinkedList<T>::operator =(LinkedList<T> &rhs) {
+	
+	ListNode<T> *src = rhs.begin();
+	while(src) {
+		push_back(src->data);
+		src = src->next;
+	}
+}
+
+
+#endif	// _LINKEDLIST_H_
diff --git a/tools/curve_draw/vmath/Makefile-part b/tools/curve_draw/vmath/Makefile-part
new file mode 100644
index 0000000..444bfe4
--- /dev/null
+++ b/tools/curve_draw/vmath/Makefile-part
@@ -0,0 +1,5 @@
+obj +=	vmath/vmath.o\
+		vmath/matrix.o\
+		vmath/vector.o\
+		vmath/quaternion.o\
+		vmath/sphvec.o
diff --git a/tools/curve_draw/vmath/matrix.cc b/tools/curve_draw/vmath/matrix.cc
new file mode 100644
index 0000000..44d0b35
--- /dev/null
+++ b/tools/curve_draw/vmath/matrix.cc
@@ -0,0 +1,704 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file matrix.cc
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Matrix math.
+ */
+
+#include <cstdio>
+#include <cmath>
+#include "matrix.h"
+#include "vector.h"
+#include "quaternion.h"
+
+using namespace std;
+
+// ----------- Matrix3x3 --------------
+
+Matrix3x3 Matrix3x3::identity_matrix = Matrix3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
+
+Matrix3x3::Matrix3x3() {
+	*this = identity_matrix;
+}
+
+Matrix3x3::Matrix3x3(	scalar_t m11, scalar_t m12, scalar_t m13,
+						scalar_t m21, scalar_t m22, scalar_t m23,
+						scalar_t m31, scalar_t m32, scalar_t m33) {
+	m[0][0] = m11; m[0][1] = m12; m[0][2] = m13;
+	m[1][0] = m21; m[1][1] = m22; m[1][2] = m23;
+	m[2][0] = m31; m[2][1] = m32; m[2][2] = m33;
+}
+
+Matrix3x3::Matrix3x3(const Matrix4x4 &mat4x4) {
+	for(int i=0; i<3; i++) {
+		for(int j=0; j<3; j++) {
+			m[i][j] = mat4x4[i][j];
+		}
+	}
+}
+
+Matrix3x3 operator +(const Matrix3x3 &m1, const Matrix3x3 &m2) {
+	Matrix3x3 res;
+	const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
+	scalar_t *dest = res.m[0];
+	
+	for(int i=0; i<9; i++) {
+		*dest++ = *op1++ + *op2++;
+	}
+	return res;
+}
+
+Matrix3x3 operator -(const Matrix3x3 &m1, const Matrix3x3 &m2) {
+	Matrix3x3 res;
+	const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
+	scalar_t *dest = res.m[0];
+	
+	for(int i=0; i<9; i++) {
+		*dest++ = *op1++ - *op2++;
+	}
+	return res;
+}
+
+Matrix3x3 operator *(const Matrix3x3 &m1, const Matrix3x3 &m2) {
+	Matrix3x3 res;
+	for(int i=0; i<3; i++) {
+		for(int j=0; j<3; j++) {
+			res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j];
+		}
+	}
+	return res;
+}
+
+void operator +=(Matrix3x3 &m1, const Matrix3x3 &m2) {
+	scalar_t *op1 = m1.m[0];
+	const scalar_t *op2 = m2.m[0];
+	
+	for(int i=0; i<9; i++) {
+		*op1++ += *op2++;
+	}
+}
+
+void operator -=(Matrix3x3 &m1, const Matrix3x3 &m2) {
+	scalar_t *op1 = m1.m[0];
+	const scalar_t *op2 = m2.m[0];
+	
+	for(int i=0; i<9; i++) {
+		*op1++ -= *op2++;
+	}
+}
+
+void operator *=(Matrix3x3 &m1, const Matrix3x3 &m2) {
+	Matrix3x3 res;
+	for(int i=0; i<3; i++) {
+		for(int j=0; j<3; j++) {
+			res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j];
+		}
+	}
+	memcpy(m1.m, res.m, 9 * sizeof(scalar_t));
+}
+
+Matrix3x3 operator *(const Matrix3x3 &mat, scalar_t scalar) {
+	Matrix3x3 res;
+	const scalar_t *mptr = mat.m[0];
+	scalar_t *dptr = res.m[0];
+	
+	for(int i=0; i<9; i++) {
+		*dptr++ = *mptr++ * scalar;
+	}
+	return res;
+}
+
+Matrix3x3 operator *(scalar_t scalar, const Matrix3x3 &mat) {
+	Matrix3x3 res;
+	const scalar_t *mptr = mat.m[0];
+	scalar_t *dptr = res.m[0];
+	
+	for(int i=0; i<9; i++) {
+		*dptr++ = *mptr++ * scalar;
+	}
+	return res;
+}
+
+void operator *=(Matrix3x3 &mat, scalar_t scalar) {
+	scalar_t *mptr = mat.m[0];
+	
+	for(int i=0; i<9; i++) {
+		*mptr++ *= scalar;
+	}
+}
+
+void Matrix3x3::translate(const Vector2 &trans) {
+	Matrix3x3 tmat(1, 0, trans.x, 0, 1, trans.y, 0, 0, 1);
+	*this *= tmat;
+}
+
+void Matrix3x3::set_translation(const Vector2 &trans) {
+	*this = Matrix3x3(1, 0, trans.x, 0, 1, trans.y, 0, 0, 1);
+}
+
+void Matrix3x3::rotate(scalar_t angle) {
+	scalar_t cos_a = cos(angle);
+	scalar_t sin_a = sin(angle);
+	Matrix3x3 rmat(	cos_a, 	-sin_a, 	0,
+					sin_a, 	cos_a, 		0,
+					0,		0,			1);
+	*this *= rmat;
+}
+
+void Matrix3x3::set_rotation(scalar_t angle) {
+	scalar_t cos_a = cos(angle);
+	scalar_t sin_a = sin(angle);
+	*this = Matrix3x3(cos_a, -sin_a, 0, sin_a, cos_a, 0, 0, 0, 1);
+}
+
+void Matrix3x3::rotate(const Vector3 &euler_angles) {
+	Matrix3x3 xrot, yrot, zrot;
+	
+	xrot = Matrix3x3(	1,			0,					0,
+						0,	cos(euler_angles.x),	-sin(euler_angles.x),
+						0,	sin(euler_angles.x),	cos(euler_angles.x));
+	
+	yrot = Matrix3x3(	cos(euler_angles.y),	0,	sin(euler_angles.y),
+								0,				1,				0,
+						-sin(euler_angles.y),	0,	cos(euler_angles.y));
+	
+	zrot = Matrix3x3(	cos(euler_angles.z),	-sin(euler_angles.z),	0,
+						sin(euler_angles.z),	cos(euler_angles.z),	0,
+								0,						0,				1);
+	
+	*this *= xrot * yrot * zrot;
+}
+
+void Matrix3x3::set_rotation(const Vector3 &euler_angles) {
+	Matrix3x3 xrot, yrot, zrot;
+	
+	xrot = Matrix3x3(	1,			0,					0,
+						0,	cos(euler_angles.x),	-sin(euler_angles.x),
+						0,	sin(euler_angles.x),	cos(euler_angles.x));
+	
+	yrot = Matrix3x3(	cos(euler_angles.y),	0,	sin(euler_angles.y),
+								0,				1,				0,
+						-sin(euler_angles.y),	0,	cos(euler_angles.y));
+	
+	zrot = Matrix3x3(	cos(euler_angles.z),	-sin(euler_angles.z),	0,
+						sin(euler_angles.z),	cos(euler_angles.z),	0,
+								0,						0,				1);
+	
+	*this = xrot * yrot * zrot;
+}
+
+void Matrix3x3::rotate(const Vector3 &axis, scalar_t angle) {
+	scalar_t sina = (scalar_t)sin(angle);
+	scalar_t cosa = (scalar_t)cos(angle);
+	scalar_t invcosa = 1-cosa;
+	scalar_t nxsq = axis.x * axis.x;
+	scalar_t nysq = axis.y * axis.y;
+	scalar_t nzsq = axis.z * axis.z;
+
+	Matrix3x3 xform;
+	xform.m[0][0] = nxsq + (1-nxsq) * cosa;
+	xform.m[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
+	xform.m[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
+	xform.m[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
+	xform.m[1][1] = nysq + (1-nysq) * cosa;
+	xform.m[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
+	xform.m[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
+	xform.m[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
+	xform.m[2][2] = nzsq + (1-nzsq) * cosa;
+
+	*this *= xform;
+}
+
+void Matrix3x3::set_rotation(const Vector3 &axis, scalar_t angle) {
+	scalar_t sina = (scalar_t)sin(angle);
+	scalar_t cosa = (scalar_t)cos(angle);
+	scalar_t invcosa = 1-cosa;
+	scalar_t nxsq = axis.x * axis.x;
+	scalar_t nysq = axis.y * axis.y;
+	scalar_t nzsq = axis.z * axis.z;
+
+	reset_identity();
+	m[0][0] = nxsq + (1-nxsq) * cosa;
+	m[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
+	m[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
+	m[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
+	m[1][1] = nysq + (1-nysq) * cosa;
+	m[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
+	m[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
+	m[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
+	m[2][2] = nzsq + (1-nzsq) * cosa;
+}
+
+void Matrix3x3::scale(const Vector3 &scale_vec) {
+	Matrix3x3 smat(	scale_vec.x, 0, 0,
+					0, scale_vec.y, 0,
+					0, 0, scale_vec.z);
+	*this *= smat;
+}
+
+void Matrix3x3::set_scaling(const Vector3 &scale_vec) {
+	*this = Matrix3x3(	scale_vec.x, 0, 0,
+						0, scale_vec.y, 0,
+						0, 0, scale_vec.z);
+}
+
+void Matrix3x3::set_column_vector(const Vector3 &vec, unsigned int col_index) {
+	m[0][col_index] = vec.x;
+	m[1][col_index] = vec.y;
+	m[2][col_index] = vec.z;
+}
+
+void Matrix3x3::set_row_vector(const Vector3 &vec, unsigned int row_index) {
+	m[row_index][0] = vec.x;
+	m[row_index][1] = vec.y;
+	m[row_index][2] = vec.z;
+}
+
+Vector3 Matrix3x3::get_column_vector(unsigned int col_index) const {
+	return Vector3(m[0][col_index], m[1][col_index], m[2][col_index]);
+}
+
+Vector3 Matrix3x3::get_row_vector(unsigned int row_index) const {
+	return Vector3(m[row_index][0], m[row_index][1], m[row_index][2]);
+}
+
+void Matrix3x3::transpose() {
+	Matrix3x3 tmp = *this;
+	for(int i=0; i<3; i++) {
+		for(int j=0; j<3; j++) {
+			m[i][j] = tmp[j][i];
+		}
+	}
+}
+
+Matrix3x3 Matrix3x3::transposed() const {
+	Matrix3x3 res;
+	for(int i=0; i<3; i++) {
+		for(int j=0; j<3; j++) {
+			res[i][j] = m[j][i];
+		}
+	}
+	return res;
+}
+
+scalar_t Matrix3x3::determinant() const {
+	return 	m[0][0] * (m[1][1]*m[2][2] - m[1][2]*m[2][1]) -
+			m[0][1] * (m[1][0]*m[2][2] - m[1][2]*m[2][0]) +
+			m[0][2] * (m[1][0]*m[2][1] - m[1][1]*m[2][0]);
+}
+
+Matrix3x3 Matrix3x3::inverse() const {
+	// TODO: implement 3x3 inverse
+	return *this;
+}
+
+ostream &operator <<(ostream &out, const Matrix3x3 &mat) {
+	for(int i=0; i<3; i++) {
+		char str[100];
+		sprintf(str, "[ %12.5f %12.5f %12.5f ]\n", (float)mat.m[i][0], (float)mat.m[i][1], (float)mat.m[i][2]);
+		out << str;
+	}
+	return out;
+}
+
+
+
+////////////////////// Matrix4x4 implementation ///////////////////////
+
+Matrix4x4 Matrix4x4::identity_matrix = Matrix4x4(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
+
+Matrix4x4::Matrix4x4() {
+	*this = identity_matrix;
+}
+
+Matrix4x4::Matrix4x4(	scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14,
+						scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24,
+						scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34,
+						scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44) {
+	m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[0][3] = m14;
+	m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[1][3] = m24;
+	m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; m[2][3] = m34;
+	m[3][0] = m41; m[3][1] = m42; m[3][2] = m43; m[3][3] = m44;
+	
+	glmatrix = 0;
+}
+
+Matrix4x4::Matrix4x4(const Matrix3x3 &mat3x3) {
+	reset_identity();
+	for(int i=0; i<3; i++) {
+		for(int j=0; j<3; j++) {
+			m[i][j] = mat3x3[i][j];
+		}
+	}
+
+	glmatrix = 0;
+}
+
+Matrix4x4::~Matrix4x4() {
+	if(glmatrix) delete [] glmatrix;
+}
+
+Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2) {
+	Matrix4x4 res;
+	const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
+	scalar_t *dest = res.m[0];
+	
+	for(int i=0; i<16; i++) {
+		*dest++ = *op1++ + *op2++;
+	}
+	return res;
+}
+
+Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2) {
+	Matrix4x4 res;
+	const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
+	scalar_t *dest = res.m[0];
+	
+	for(int i=0; i<16; i++) {
+		*dest++ = *op1++ - *op2++;
+	}
+	return res;
+}
+
+Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2) {
+	Matrix4x4 res;
+	
+	for(int i=0; i<4; i++) {
+		for(int j=0; j<4; j++) {
+			res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j];
+		}
+	}
+
+	return res;
+}
+
+void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2) {
+	scalar_t *op1 = m1.m[0];
+	const scalar_t *op2 = m2.m[0];
+	
+	for(int i=0; i<16; i++) {
+		*op1++ += *op2++;
+	}
+}
+
+void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2) {
+	scalar_t *op1 = m1.m[0];
+	const scalar_t *op2 = m2.m[0];
+	
+	for(int i=0; i<16; i++) {
+		*op1++ -= *op2++;
+	}
+}
+
+void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2) {
+	Matrix4x4 res;
+	for(int i=0; i<4; i++) {
+		for(int j=0; j<4; j++) {
+			res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j];
+		}
+	}
+	memcpy(m1.m, res.m, 16 * sizeof(scalar_t));
+}
+
+Matrix4x4 operator *(const Matrix4x4 &mat, scalar_t scalar) {
+	Matrix4x4 res;
+	const scalar_t *mptr = mat.m[0];
+	scalar_t *dptr = res.m[0];
+	
+	for(int i=0; i<16; i++) {
+		*dptr++ = *mptr++ * scalar;
+	}
+	return res;
+}
+
+Matrix4x4 operator *(scalar_t scalar, const Matrix4x4 &mat) {
+	Matrix4x4 res;
+	const scalar_t *mptr = mat.m[0];
+	scalar_t *dptr = res.m[0];
+	
+	for(int i=0; i<16; i++) {
+		*dptr++ = *mptr++ * scalar;
+	}
+	return res;
+}
+
+void operator *=(Matrix4x4 &mat, scalar_t scalar) {
+	scalar_t *mptr = mat.m[0];
+	
+	for(int i=0; i<16; i++) {
+		*mptr++ *= scalar;
+	}
+}
+
+void Matrix4x4::translate(const Vector3 &trans) {
+	Matrix4x4 tmat(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1);
+	*this *= tmat;
+}
+
+void Matrix4x4::set_translation(const Vector3 &trans) {
+	*this = Matrix4x4(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1);
+}
+
+void Matrix4x4::rotate(const Vector3 &euler_angles) {
+	Matrix3x3 xrot, yrot, zrot;
+	
+	xrot = Matrix3x3(	1,			0,					0,
+						0,	cos(euler_angles.x),	-sin(euler_angles.x),
+						0,	sin(euler_angles.x),	cos(euler_angles.x));
+	
+	yrot = Matrix3x3(	cos(euler_angles.y),	0,	sin(euler_angles.y),
+								0,				1,				0,
+						-sin(euler_angles.y),	0,	cos(euler_angles.y));
+	
+	zrot = Matrix3x3(	cos(euler_angles.z),	-sin(euler_angles.z),	0,
+						sin(euler_angles.z),	cos(euler_angles.z),	0,
+								0,						0,				1);
+	
+	*this *= Matrix4x4(xrot * yrot * zrot);
+}
+
+void Matrix4x4::set_rotation(const Vector3 &euler_angles) {
+	Matrix3x3 xrot, yrot, zrot;
+	
+	xrot = Matrix3x3(	1,			0,					0,
+						0,	cos(euler_angles.x),	-sin(euler_angles.x),
+						0,	sin(euler_angles.x),	cos(euler_angles.x));
+	
+	yrot = Matrix3x3(	cos(euler_angles.y),	0,	sin(euler_angles.y),
+								0,				1,				0,
+						-sin(euler_angles.y),	0,	cos(euler_angles.y));
+	
+	zrot = Matrix3x3(	cos(euler_angles.z),	-sin(euler_angles.z),	0,
+						sin(euler_angles.z),	cos(euler_angles.z),	0,
+								0,						0,				1);
+	
+	*this = Matrix4x4(xrot * yrot * zrot);
+}
+
+void Matrix4x4::rotate(const Vector3 &axis, scalar_t angle) {
+	scalar_t sina = (scalar_t)sin(angle);
+	scalar_t cosa = (scalar_t)cos(angle);
+	scalar_t invcosa = 1-cosa;
+	scalar_t nxsq = axis.x * axis.x;
+	scalar_t nysq = axis.y * axis.y;
+	scalar_t nzsq = axis.z * axis.z;
+
+	Matrix3x3 xform;
+	xform[0][0] = nxsq + (1-nxsq) * cosa;
+	xform[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
+	xform[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
+	xform[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
+	xform[1][1] = nysq + (1-nysq) * cosa;
+	xform[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
+	xform[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
+	xform[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
+	xform[2][2] = nzsq + (1-nzsq) * cosa;
+
+	*this *= Matrix4x4(xform);
+}
+
+void Matrix4x4::set_rotation(const Vector3 &axis, scalar_t angle) {
+	scalar_t sina = (scalar_t)sin(angle);
+	scalar_t cosa = (scalar_t)cos(angle);
+	scalar_t invcosa = 1-cosa;
+	scalar_t nxsq = axis.x * axis.x;
+	scalar_t nysq = axis.y * axis.y;
+	scalar_t nzsq = axis.z * axis.z;
+
+	reset_identity();
+	m[0][0] = nxsq + (1-nxsq) * cosa;
+	m[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
+	m[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
+	m[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
+	m[1][1] = nysq + (1-nysq) * cosa;
+	m[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
+	m[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
+	m[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
+	m[2][2] = nzsq + (1-nzsq) * cosa;
+}
+
+void Matrix4x4::scale(const Vector4 &scale_vec) {
+	Matrix4x4 smat(	scale_vec.x, 0, 0, 0,
+					0, scale_vec.y, 0, 0,
+					0, 0, scale_vec.z, 0,
+					0, 0, 0, scale_vec.w);
+	*this *= smat;
+}
+
+void Matrix4x4::set_scaling(const Vector4 &scale_vec) {
+	*this = Matrix4x4(	scale_vec.x, 0, 0, 0,
+						0, scale_vec.y, 0, 0,
+						0, 0, scale_vec.z, 0,
+						0, 0, 0, scale_vec.w);
+}
+
+void Matrix4x4::set_column_vector(const Vector4 &vec, unsigned int col_index) {
+	m[0][col_index] = vec.x;
+	m[1][col_index] = vec.y;
+	m[2][col_index] = vec.z;
+	m[3][col_index] = vec.w;
+}
+
+void Matrix4x4::set_row_vector(const Vector4 &vec, unsigned int row_index) {
+	m[row_index][0] = vec.x;
+	m[row_index][1] = vec.y;
+	m[row_index][2] = vec.z;
+	m[row_index][3] = vec.w;
+}
+
+Vector4 Matrix4x4::get_column_vector(unsigned int col_index) const {
+	return Vector4(m[0][col_index], m[1][col_index], m[2][col_index], m[3][col_index]);
+}
+
+Vector4 Matrix4x4::get_row_vector(unsigned int row_index) const {
+	return Vector4(m[row_index][0], m[row_index][1], m[row_index][2], m[row_index][3]);
+}
+
+void Matrix4x4::transpose() {
+	Matrix4x4 tmp = *this;
+	for(int i=0; i<4; i++) {
+		for(int j=0; j<4; j++) {
+			m[i][j] = tmp[j][i];
+		}
+	}
+}
+
+Matrix4x4 Matrix4x4::transposed() const {
+	Matrix4x4 res;
+	for(int i=0; i<4; i++) {
+		for(int j=0; j<4; j++) {
+			res[i][j] = m[j][i];
+		}
+	}
+	return res;
+}
+
+scalar_t Matrix4x4::determinant() const {
+
+	scalar_t det11 =	(m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
+						(m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
+						(m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
+
+	scalar_t det12 =	(m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
+						(m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
+						(m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
+
+	scalar_t det13 =	(m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
+						(m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
+						(m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
+
+	scalar_t det14 =	(m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
+						(m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
+						(m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
+
+	return m[0][0] * det11 - m[0][1] * det12 + m[0][2] * det13 - m[0][3] * det14;
+}
+
+
+Matrix4x4 Matrix4x4::adjoint() const {
+
+	Matrix4x4 coef;
+
+	coef.m[0][0] =	(m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
+					(m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
+					(m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
+	coef.m[0][1] =	(m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
+					(m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
+					(m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
+	coef.m[0][2] =	(m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
+					(m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
+					(m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
+	coef.m[0][3] =	(m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
+					(m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
+					(m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
+
+	coef.m[1][0] =	(m[0][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
+					(m[0][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
+					(m[0][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
+	coef.m[1][1] =	(m[0][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
+					(m[0][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
+					(m[0][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
+	coef.m[1][2] =	(m[0][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
+					(m[0][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
+					(m[0][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
+	coef.m[1][3] =	(m[0][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
+					(m[0][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
+					(m[0][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
+
+	coef.m[2][0] =	(m[0][1] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) -
+					(m[0][2] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) +
+					(m[0][3] * (m[1][1] * m[3][2] - m[3][1] * m[1][2]));
+	coef.m[2][1] =	(m[0][0] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) -
+					(m[0][2] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) +
+					(m[0][3] * (m[1][0] * m[3][2] - m[3][0] * m[1][2]));
+	coef.m[2][2] =	(m[0][0] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) -
+					(m[0][1] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) +
+					(m[0][3] * (m[1][0] * m[3][1] - m[3][0] * m[1][1]));
+	coef.m[2][3] =	(m[0][0] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])) -
+					(m[0][1] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])) +
+					(m[0][2] * (m[1][0] * m[3][1] - m[3][0] * m[1][1]));
+
+	coef.m[3][0] =	(m[0][1] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) -
+					(m[0][2] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) +
+					(m[0][3] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]));
+	coef.m[3][1] =	(m[0][0] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) -
+					(m[0][2] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) +
+					(m[0][3] * (m[1][0] * m[2][2] - m[2][0] * m[1][2]));
+	coef.m[3][2] =	(m[0][0] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) -
+					(m[0][1] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) +
+					(m[0][3] * (m[1][0] * m[2][1] - m[2][0] * m[1][1]));
+	coef.m[3][3] =	(m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])) -
+					(m[0][1] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])) +
+					(m[0][2] * (m[1][0] * m[2][1] - m[2][0] * m[1][1]));
+
+	coef.transpose();
+
+	for(int i=0; i<4; i++) {
+		for(int j=0; j<4; j++) {
+			coef.m[i][j] = j%2 ? -coef.m[i][j] : coef.m[i][j];
+			if(i%2) coef.m[i][j] = -coef.m[i][j];
+		}
+	}
+
+	return coef;
+}
+
+Matrix4x4 Matrix4x4::inverse() const {
+
+	Matrix4x4 AdjMat = adjoint();
+
+	return AdjMat * (1.0f / determinant());
+}
+
+const scalar_t *Matrix4x4::opengl_matrix() const {
+	return (const scalar_t*)m;
+}
+
+ostream &operator <<(ostream &out, const Matrix4x4 &mat) {
+	for(int i=0; i<4; i++) {
+		char str[100];
+		sprintf(str, "[ %12.5f %12.5f %12.5f %12.5f ]\n", (float)mat.m[i][0], (float)mat.m[i][1], (float)mat.m[i][2], (float)mat.m[i][3]);
+		out << str;
+	}
+	return out;
+}
diff --git a/tools/curve_draw/vmath/matrix.h b/tools/curve_draw/vmath/matrix.h
new file mode 100644
index 0000000..6c7c465
--- /dev/null
+++ b/tools/curve_draw/vmath/matrix.h
@@ -0,0 +1,162 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file matrix.h
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Matrix math.
+ */
+
+#ifndef MATRIX_H_
+#define MATRIX_H_
+
+#include <iostream>
+#include "vmath_types.h"
+
+/** 3x3 matrix */
+class Matrix3x3 {
+private:
+	scalar_t m[3][3];
+public:
+	
+	static Matrix3x3 identity_matrix;
+
+	Matrix3x3();
+	Matrix3x3(	scalar_t m11, scalar_t m12, scalar_t m13,
+				scalar_t m21, scalar_t m22, scalar_t m23,
+				scalar_t m31, scalar_t m32, scalar_t m33);
+	
+	Matrix3x3(const Matrix4x4 &mat4x4);
+	
+	// binary operations matrix (op) matrix
+	friend Matrix3x3 operator +(const Matrix3x3 &m1, const Matrix3x3 &m2);
+	friend Matrix3x3 operator -(const Matrix3x3 &m1, const Matrix3x3 &m2);
+	friend Matrix3x3 operator *(const Matrix3x3 &m1, const Matrix3x3 &m2);
+	
+	friend void operator +=(Matrix3x3 &m1, const Matrix3x3 &m2);
+	friend void operator -=(Matrix3x3 &m1, const Matrix3x3 &m2);
+	friend void operator *=(Matrix3x3 &m1, const Matrix3x3 &m2);
+	
+	// binary operations matrix (op) scalar and scalar (op) matrix
+	friend Matrix3x3 operator *(const Matrix3x3 &mat, scalar_t scalar);
+	friend Matrix3x3 operator *(scalar_t scalar, const Matrix3x3 &mat);
+	
+	friend void operator *=(Matrix3x3 &mat, scalar_t scalar);
+	
+	inline scalar_t *operator [](int index);
+	inline const scalar_t *operator [](int index) const;
+	
+	inline void reset_identity();
+	
+	void translate(const Vector2 &trans);
+	void set_translation(const Vector2 &trans);
+	
+	void rotate(scalar_t angle);						// 2d rotation
+	void rotate(const Vector3 &euler_angles);			// 3d rotation with euler angles
+	void rotate(const Vector3 &axis, scalar_t angle);	// 3d axis/angle rotation
+	void set_rotation(scalar_t angle);
+	void set_rotation(const Vector3 &euler_angles);
+	void set_rotation(const Vector3 &axis, scalar_t angle);
+	
+	void scale(const Vector3 &scale_vec);
+	void set_scaling(const Vector3 &scale_vec);
+	
+	void set_column_vector(const Vector3 &vec, unsigned int col_index);
+	void set_row_vector(const Vector3 &vec, unsigned int row_index);
+	Vector3 get_column_vector(unsigned int col_index) const;
+	Vector3 get_row_vector(unsigned int row_index) const;
+
+	void transpose();
+	Matrix3x3 transposed() const;	
+	scalar_t determinant() const;
+	Matrix3x3 inverse() const;
+	
+	friend std::ostream &operator <<(std::ostream &out, const Matrix3x3 &mat);
+};
+
+/** 4x4 matrix */
+class Matrix4x4 {
+private:
+	scalar_t m[4][4];
+	mutable float *glmatrix;
+public:
+	
+	static Matrix4x4 identity_matrix;
+
+	Matrix4x4();
+	Matrix4x4(	scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14,
+				scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24,
+				scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34,
+				scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44);
+	
+	Matrix4x4(const Matrix3x3 &mat3x3);
+	~Matrix4x4();
+	
+	// binary operations matrix (op) matrix
+	friend Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2);
+	friend Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2);
+	friend Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2);
+	
+	friend void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2);
+	friend void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2);
+	friend void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2);
+	
+	// binary operations matrix (op) scalar and scalar (op) matrix
+	friend Matrix4x4 operator *(const Matrix4x4 &mat, scalar_t scalar);
+	friend Matrix4x4 operator *(scalar_t scalar, const Matrix4x4 &mat);
+	
+	friend void operator *=(Matrix4x4 &mat, scalar_t scalar);
+	
+	inline scalar_t *operator [](int index);
+	inline const scalar_t *operator [](int index) const;
+	
+	inline void reset_identity();
+	
+	void translate(const Vector3 &trans);
+	void set_translation(const Vector3 &trans);
+	
+	void rotate(const Vector3 &euler_angles);			// 3d rotation with euler angles
+	void rotate(const Vector3 &axis, scalar_t angle);	// 3d axis/angle rotation
+	void set_rotation(const Vector3 &euler_angles);
+	void set_rotation(const Vector3 &axis, scalar_t angle);
+	
+	void scale(const Vector4 &scale_vec);
+	void set_scaling(const Vector4 &scale_vec);
+	
+	void set_column_vector(const Vector4 &vec, unsigned int col_index);
+	void set_row_vector(const Vector4 &vec, unsigned int row_index);
+	Vector4 get_column_vector(unsigned int col_index) const;
+	Vector4 get_row_vector(unsigned int row_index) const;
+	
+	void transpose();
+	Matrix4x4 transposed() const;
+	scalar_t determinant() const;
+	Matrix4x4 adjoint() const;
+	Matrix4x4 inverse() const;
+	
+	const scalar_t *opengl_matrix() const;
+		
+	friend std::ostream &operator <<(std::ostream &out, const Matrix4x4 &mat);
+};
+
+#include "matrix.inl"
+
+#endif	// MATRIX_HPP_
diff --git a/tools/curve_draw/vmath/matrix.inl b/tools/curve_draw/vmath/matrix.inl
new file mode 100644
index 0000000..58e93b9
--- /dev/null
+++ b/tools/curve_draw/vmath/matrix.inl
@@ -0,0 +1,50 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file matrix.inl
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Matrix inline functions.
+ */
+
+inline scalar_t *Matrix3x3::operator [](int index) {
+	return m[index];
+}
+
+inline const scalar_t *Matrix3x3::operator [](int index) const {
+	return m[index];
+}
+
+inline void Matrix3x3::reset_identity() {
+	memcpy(this->m, identity_matrix.m, 9 * sizeof(scalar_t));
+}
+
+inline scalar_t *Matrix4x4::operator [](int index) {
+	return m[index];
+}
+
+inline const scalar_t *Matrix4x4::operator [](int index) const {
+	return m[index];
+}
+
+inline void Matrix4x4::reset_identity() {
+	memcpy(this->m, identity_matrix.m, 16 * sizeof(scalar_t));
+}
diff --git a/tools/curve_draw/vmath/quaternion.cc b/tools/curve_draw/vmath/quaternion.cc
new file mode 100644
index 0000000..830e503
--- /dev/null
+++ b/tools/curve_draw/vmath/quaternion.cc
@@ -0,0 +1,182 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file quaternion.cxx
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Quaternion math.
+ */
+
+#include "quaternion.h"
+
+Quaternion::Quaternion() {
+	s = 1.0;
+	v.x = v.y = v.z = 0.0;
+}
+
+Quaternion::Quaternion(scalar_t s, scalar_t x, scalar_t y, scalar_t z) {
+	v.x = x;
+	v.y = y;
+	v.z = z;
+	this->s = s;
+}
+
+Quaternion::Quaternion(scalar_t s, const Vector3 &v) {
+	this->s = s;
+	this->v = v;
+}
+
+Quaternion::Quaternion(const Vector3 &axis, scalar_t angle) {
+	set_rotation(axis, angle);
+}
+
+Quaternion Quaternion::operator +(const Quaternion &quat) const {
+	return Quaternion(s + quat.s, v + quat.v);
+}
+
+Quaternion Quaternion::operator -(const Quaternion &quat) const {
+	return Quaternion(s - quat.s, v - quat.v);
+}
+
+Quaternion Quaternion::operator -() const {
+	return Quaternion(-s, -v);
+}
+
+/** Quaternion Multiplication:
+ * Q1*Q2 = [s1*s2 - v1.v2,  s1*v2 + s2*v1 + v1(x)v2]
+ */
+Quaternion Quaternion::operator *(const Quaternion &quat) const {
+	Quaternion newq;	
+	newq.s = s * quat.s - dot_product(v, quat.v);
+	newq.v = quat.v * s + v * quat.s + cross_product(v, quat.v);	
+	return newq;
+}
+
+void Quaternion::operator +=(const Quaternion &quat) {
+	*this = Quaternion(s + quat.s, v + quat.v);
+}
+
+void Quaternion::operator -=(const Quaternion &quat) {
+	*this = Quaternion(s - quat.s, v - quat.v);
+}
+
+void Quaternion::operator *=(const Quaternion &quat) {
+	*this = *this * quat;
+}
+
+void Quaternion::reset_identity() {
+	s = 1.0;
+	v.x = v.y = v.z = 0.0;
+}
+
+Quaternion Quaternion::conjugate() const {
+	return Quaternion(s, -v);
+}
+
+scalar_t Quaternion::length() const {
+	return (scalar_t)sqrt(v.x*v.x + v.y*v.y + v.z*v.z + s*s);
+}
+
+/** Q * ~Q = ||Q||^2 */
+scalar_t Quaternion::length_sq() const {
+	return v.x*v.x + v.y*v.y + v.z*v.z + s*s;
+}
+
+void Quaternion::normalize() {
+	scalar_t len = (scalar_t)sqrt(v.x*v.x + v.y*v.y + v.z*v.z + s*s);
+	v.x /= len;
+	v.y /= len;
+	v.z /= len;
+	s /= len;
+}
+
+Quaternion Quaternion::normalized() const {
+	Quaternion nq = *this;
+	scalar_t len = (scalar_t)sqrt(v.x*v.x + v.y*v.y + v.z*v.z + s*s);
+	nq.v.x /= len;
+	nq.v.y /= len;
+	nq.v.z /= len;
+	nq.s /= len;
+	return nq;
+}
+
+/** Quaternion Inversion: Q^-1 = ~Q / ||Q||^2 */
+Quaternion Quaternion::inverse() const {
+	Quaternion inv = conjugate();
+	scalar_t lensq = length_sq();
+	inv.v /= lensq;
+	inv.s /= lensq;
+
+	return inv;
+}
+
+
+void Quaternion::set_rotation(const Vector3 &axis, scalar_t angle) {
+	scalar_t HalfAngle = angle / 2.0;
+	s = cos(HalfAngle);
+	v = axis * sin(HalfAngle);
+}
+
+void Quaternion::rotate(const Vector3 &axis, scalar_t angle) {
+	Quaternion q;
+	scalar_t HalfAngle = angle / 2.0;
+	q.s = cos(HalfAngle);
+	q.v = axis * sin(HalfAngle);
+
+	*this *= q;
+}
+
+void Quaternion::rotate(const Quaternion &q) {
+	*this = q * *this * q.conjugate();
+}
+
+Matrix3x3 Quaternion::get_rotation_matrix() const {
+	return Matrix3x3(	1.0 - 2.0 * v.y*v.y - 2.0 * v.z*v.z,	2.0 * v.x * v.y + 2.0 * s * v.z,		2.0 * v.z * v.x - 2.0 * s * v.y,
+						2.0 * v.x * v.y - 2.0 * s * v.z,		1.0 - 2.0 * v.x*v.x - 2.0 * v.z*v.z,	2.0 * v.y * v.z + 2.0 * s * v.x,
+						2.0 * v.z * v.x + 2.0 * s * v.y,		2.0 * v.y * v.z - 2.0 * s * v.x,		1.0 - 2.0 * v.x*v.x - 2.0 * v.y*v.y);
+}
+
+
+/** Spherical linear interpolation (slerp) */
+Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t) {
+	scalar_t dot = dot_product(Vector4(q1.s, q1.v.x, q1.v.y, q1.v.z), Vector4(q2.s, q2.v.x, q2.v.y, q2.v.z));
+	
+	if(fabs(1.0 - dot) < xsmall_number) return q1;	// avoids divisions by zero later on if q1 == q2
+	
+	if(dot >= 0.0f) {
+		scalar_t angle = acos(dot);
+		scalar_t coef1 = (angle * sin(1.0 - t)) / sin(angle);
+		scalar_t coef2 = sin(angle * t) / sin(angle);
+		return Quaternion(q1.s * coef1 + q2.s * coef2, q1.v * coef1 + q2.v * coef2).normalized();
+	} else {
+		scalar_t angle = acos(-dot);
+		scalar_t coef1 = (angle * sin(1.0 - t)) / sin(angle);
+		scalar_t coef2 = sin(angle * t) / sin(angle);
+		return Quaternion(q2.s * coef1 + q1.s * coef2, q2.v * coef1 + q1.v * coef2).normalized();
+	}
+}
+	
+
+
+std::ostream &operator <<(std::ostream &out, const Quaternion &q) {
+	out << "(" << q.s << ", " << q.v << ")";
+	return out;
+}
diff --git a/tools/curve_draw/vmath/quaternion.h b/tools/curve_draw/vmath/quaternion.h
new file mode 100644
index 0000000..f8fc3ad
--- /dev/null
+++ b/tools/curve_draw/vmath/quaternion.h
@@ -0,0 +1,83 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file quaternion.h
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Quaternion math.
+ */
+
+#ifndef QUATERNION_H_
+#define QUATERNION_H_
+
+#include <iostream>
+#include "vmath_types.h"
+#include "vector.h"
+
+/** Quaternion */
+class Quaternion {
+public:
+	scalar_t s;
+	Vector3 v;
+
+	Quaternion();
+	Quaternion(scalar_t s, const Vector3 &v);
+	Quaternion(scalar_t s, scalar_t x, scalar_t y, scalar_t z);
+	Quaternion(const Vector3 &axis, scalar_t angle);
+
+	Quaternion operator +(const Quaternion &quat) const;
+	Quaternion operator -(const Quaternion &quat) const;
+	Quaternion operator -() const;
+	Quaternion operator *(const Quaternion &quat) const;
+	
+	void operator +=(const Quaternion &quat);
+	void operator -=(const Quaternion &quat);
+	void operator *=(const Quaternion &quat);
+
+	void reset_identity();
+
+	Quaternion conjugate() const;
+
+	scalar_t length() const;
+	scalar_t length_sq() const;
+	
+	void normalize();
+	Quaternion normalized() const;
+
+	Quaternion inverse() const;
+
+	void set_rotation(const Vector3 &axis, scalar_t angle);
+	void rotate(const Vector3 &axis, scalar_t angle);
+	/* note: this is a totally different operation from the above
+	 * this treats the quaternion as signifying direction and rotates
+	 * it by a rotation quaternion by rot * q * rot'
+	 */
+	void rotate(const Quaternion &q);
+
+	Matrix3x3 get_rotation_matrix() const;
+	
+	friend Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t);
+	
+	friend std::ostream &operator <<(std::ostream &out, const Quaternion &q);
+};
+
+
+#endif	// QUATERNION_H_
diff --git a/tools/curve_draw/vmath/sphvec.cc b/tools/curve_draw/vmath/sphvec.cc
new file mode 100644
index 0000000..8521bd1
--- /dev/null
+++ b/tools/curve_draw/vmath/sphvec.cc
@@ -0,0 +1,55 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file sphvec.cc
+ * @author John Tsiombikas
+ * @date 7 July 2005
+ *
+ * Spherical coordinates.
+ */
+
+#include "sphvec.h"
+#include "vector.h"
+
+/**
+ * @param theta 0 <= theta <= 2pi, the angle around Y axis.
+ * @param phi 0 <= phi <= pi, the angle from Y axis.
+ * @param r Radius.
+ */
+SphVector::SphVector(scalar_t theta, scalar_t phi, scalar_t r) {
+	this->theta = theta;
+	this->phi = phi;
+	this->r = r;
+}
+
+/** Constructs a spherical coordinate vector from a cartesian vector */
+SphVector::SphVector(const Vector3 &cvec) {
+	*this = cvec;
+}
+
+/** Assignment operator that converts cartesian to spherical coords */
+SphVector &SphVector::operator =(const Vector3 &cvec) {
+	r = cvec.length();
+	//theta = atan2(cvec.y, cvec.x);
+	theta = atan2(cvec.z, cvec.x);
+	//phi = acos(cvec.z / r);
+	phi = acos(cvec.y / r);
+	return *this;
+}
diff --git a/tools/curve_draw/vmath/sphvec.h b/tools/curve_draw/vmath/sphvec.h
new file mode 100644
index 0000000..dc83c76
--- /dev/null
+++ b/tools/curve_draw/vmath/sphvec.h
@@ -0,0 +1,43 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file sphvec.h
+ * @author John Tsiombikas
+ * @date 7 July 2005
+ *
+ * Spherical coordinates.
+ */
+
+#ifndef SPHVEC_H_
+#define SPHVEC_H_
+
+#include "vmath_types.h"
+
+/** Vector in spherical coordinates */
+class SphVector {
+public:
+	scalar_t theta, phi, r;
+
+	SphVector(scalar_t theta = 0.0, scalar_t phi = 0.0, scalar_t r = 1.0);
+	SphVector(const Vector3 &cvec);
+	SphVector &operator =(const Vector3 &cvec);
+};
+
+#endif	// SPHVEC_H_
diff --git a/tools/curve_draw/vmath/vector.cc b/tools/curve_draw/vmath/vector.cc
new file mode 100644
index 0000000..7ce5697
--- /dev/null
+++ b/tools/curve_draw/vmath/vector.cc
@@ -0,0 +1,290 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file vector.cxx
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Vector math.
+ */
+
+#include "vector.h"
+#include "sphvec.h"
+
+// ---------- Vector2 -----------
+
+Vector2::Vector2(scalar_t x, scalar_t y) {
+	this->x = x;
+	this->y = y;
+}
+
+Vector2::Vector2(const Vector3 &vec) {
+	x = vec.x;
+	y = vec.y;
+}
+
+Vector2::Vector2(const Vector4 &vec) {
+	x = vec.x;
+	y = vec.y;
+}
+
+void Vector2::normalize() {
+	scalar_t len = length();
+	x /= len;
+	y /= len;
+}
+
+Vector2 Vector2::normalized() const {	
+	scalar_t len = length();
+	return Vector2(x / len, y / len);
+}
+
+void Vector2::transform(const Matrix3x3 &mat) {
+	scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2];
+	y = mat[1][0] * x + mat[1][1] * y + mat[1][2];
+	x = nx;
+}
+
+Vector2 Vector2::transformed(const Matrix3x3 &mat) const {
+	Vector2 vec;
+	vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2];
+	vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2];
+	return vec;
+}
+
+void Vector2::rotate(scalar_t angle) {
+	*this = Vector2(cos(angle) * x - sin(angle) * y, sin(angle) * x + cos(angle) * y);
+}
+
+Vector2 Vector2::rotated(scalar_t angle) const {
+	return Vector2(cos(angle) * x - sin(angle) * y, sin(angle) * x + cos(angle) * y);
+}
+	
+Vector2 Vector2::reflection(const Vector2 &normal) const {
+	return 2.0 * dot_product(*this, normal) * normal - *this;
+}
+
+Vector2 Vector2::refraction(const Vector2 &normal, scalar_t src_ior, scalar_t dst_ior) const {
+	// quick and dirty implementation :)
+	Vector3 v3refr = Vector3(this->x, this->y, 1.0).refraction(Vector3(this->x, this->y, 1), src_ior, dst_ior);
+	return Vector2(v3refr.x, v3refr.y);
+}
+
+std::ostream &operator <<(std::ostream &out, const Vector2 &vec) {
+	out << "[" << vec.x << " " << vec.y << "]";
+	return out;
+}
+
+
+
+// --------- Vector3 ----------
+
+Vector3::Vector3(scalar_t x, scalar_t y, scalar_t z) {
+	this->x = x;
+	this->y = y;
+	this->z = z;
+}
+
+Vector3::Vector3(const Vector2 &vec) {
+	x = vec.x;
+	y = vec.y;
+	z = 1;
+}
+
+Vector3::Vector3(const Vector4 &vec) {
+	x = vec.x;
+	y = vec.y;
+	z = vec.z;
+}
+
+Vector3::Vector3(const SphVector &sph) {
+	*this = sph;
+}
+
+Vector3 &Vector3::operator =(const SphVector &sph) {
+	x = sph.r * cos(sph.theta) * sin(sph.phi);
+	z = sph.r * sin(sph.theta) * sin(sph.phi);
+	y = sph.r * cos(sph.phi);
+	return *this;
+}
+
+void Vector3::normalize() {
+	scalar_t len = length();
+	x /= len;
+	y /= len;
+	z /= len;	
+}
+
+Vector3 Vector3::normalized() const {
+	scalar_t len = length();
+	return Vector3(x / len, y / len, z / len);
+}
+
+Vector3 Vector3::reflection(const Vector3 &normal) const {
+	return -(2.0 * dot_product(*this, normal) * normal - *this);
+}
+
+Vector3 Vector3::refraction(const Vector3 &normal, scalar_t src_ior, scalar_t dst_ior) const {
+	scalar_t cos_inc = dot_product(*this, -normal);
+	scalar_t ior = src_ior / dst_ior;
+
+	scalar_t radical = 1.0 + SQ(ior) * (SQ(cos_inc) - 1.0);
+
+	if(radical < 0.0) {		// total internal reflection
+		return reflection(normal);
+	}
+
+	scalar_t beta = ior * cos_inc - sqrt(radical);
+
+	return *this * ior + normal * beta;
+}
+
+void Vector3::transform(const Matrix3x3 &mat) {
+	scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z;
+	scalar_t ny = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z;
+	z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z;
+	x = nx;
+	y = ny;
+}
+
+Vector3 Vector3::transformed(const Matrix3x3 &mat) const {
+	Vector3 vec;
+	vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z;
+	vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z;
+	vec.z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z;
+	return vec;
+}
+
+void Vector3::transform(const Matrix4x4 &mat) {
+	scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3];
+	scalar_t ny = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3];
+	z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3];
+	x = nx;
+	y = ny;
+}
+
+Vector3 Vector3::transformed(const Matrix4x4 &mat) const {
+	Vector3 vec;
+	vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3];
+	vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3];
+	vec.z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3];
+	return vec;
+}
+
+void Vector3::transform(const Quaternion &quat) {
+	Quaternion vq(0.0f, *this);
+	vq = quat * vq * quat.inverse();
+	*this = vq.v;
+}
+
+Vector3 Vector3::transformed(const Quaternion &quat) const {
+	Quaternion vq(0.0f, *this);
+	vq = quat * vq * quat.inverse();
+	return vq.v;
+}
+
+void Vector3::rotate(const Vector3 &euler) {
+	Matrix4x4 rot;
+	rot.set_rotation(euler);
+	transform(rot);
+}
+
+Vector3 Vector3::rotated(const Vector3 &euler) const {
+	Matrix4x4 rot;
+	rot.set_rotation(euler);
+	return transformed(rot);
+}
+
+
+std::ostream &operator <<(std::ostream &out, const Vector3 &vec) {
+	out << "[" << vec.x << " " << vec.y << " " << vec.z << "]";
+	return out;
+}
+
+
+
+// -------------- Vector4 --------------
+Vector4::Vector4(scalar_t x, scalar_t y, scalar_t z, scalar_t w) {
+	this->x = x;
+	this->y = y;
+	this->z = z;
+	this->w = w;
+}
+
+Vector4::Vector4(const Vector2 &vec) {
+	x = vec.x;
+	y = vec.y;
+	z = 1;
+	w = 1;
+}
+
+Vector4::Vector4(const Vector3 &vec) {
+	x = vec.x;
+	y = vec.y;
+	z = vec.z;
+	w = 1;
+}
+
+void Vector4::normalize() {
+	scalar_t len = (scalar_t)sqrt(x*x + y*y + z*z + w*w);
+	x /= len;
+	y /= len;
+	z /= len;
+	w /= len;
+}
+
+Vector4 Vector4::normalized() const {
+	scalar_t len = (scalar_t)sqrt(x*x + y*y + z*z + w*w);
+	return Vector4(x / len, y / len, z / len, w / len);
+}
+
+void Vector4::transform(const Matrix4x4 &mat) {
+	scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3] * w;
+	scalar_t ny = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3] * w;
+	scalar_t nz = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3] * w;
+	w = mat[3][0] * x + mat[3][1] * y + mat[3][2] * z + mat[3][3] * w;
+	x = nx;
+	y = ny;
+	z = nz;
+}
+
+Vector4 Vector4::transformed(const Matrix4x4 &mat) const {
+	Vector4 vec;
+	vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3] * w;
+	vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3] * w;
+	vec.z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3] * w;
+	vec.w = mat[3][0] * x + mat[3][1] * y + mat[3][2] * z + mat[3][3] * w;
+	return vec;
+}
+
+// TODO: implement 4D vector reflection
+Vector4 Vector4::reflection(const Vector4 &normal) const {
+	return *this;
+}
+
+// TODO: implement 4D vector refraction
+Vector4 Vector4::refraction(const Vector4 &normal, scalar_t src_ior, scalar_t dst_ior) const {
+	return *this;
+}
+
+std::ostream &operator <<(std::ostream &out, const Vector4 &vec) {
+	out << "[" << vec.x << " " << vec.y << " " << vec.z << " " << vec.w << "]";
+	return out;
+}
diff --git a/tools/curve_draw/vmath/vector.h b/tools/curve_draw/vmath/vector.h
new file mode 100644
index 0000000..68b9c40
--- /dev/null
+++ b/tools/curve_draw/vmath/vector.h
@@ -0,0 +1,222 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file vector.h
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Vector math.
+ */
+
+#ifndef VECTOR_H_
+#define VECTOR_H_
+
+#include <iostream>
+#include "vmath_types.h"
+
+/** 2D Vector */
+class Vector2 {
+public:
+	scalar_t x, y;
+
+	Vector2(scalar_t x = 0.0, scalar_t y = 0.0);
+	Vector2(const Vector3 &vec);
+	Vector2(const Vector4 &vec);
+	
+	inline scalar_t &operator [](int elem);	
+	
+	// unary operations
+	friend inline Vector2 operator -(const Vector2 &vec);
+	
+	// binary vector (op) vector operations
+	friend inline scalar_t dot_product(const Vector2 &v1, const Vector2 &v2);
+	
+	friend inline Vector2 operator +(const Vector2 &v1, const Vector2 &v2);
+	friend inline Vector2 operator -(const Vector2 &v1, const Vector2 &v2);	
+	friend inline Vector2 operator *(const Vector2 &v1, const Vector2 &v2);
+	friend inline Vector2 operator /(const Vector2 &v1, const Vector2 &v2);
+	friend inline bool operator ==(const Vector2 &v1, const Vector2 &v2);
+	
+	friend inline void operator +=(Vector2 &v1, const Vector2 &v2);
+	friend inline void operator -=(Vector2 &v1, const Vector2 &v2);
+	friend inline void operator *=(Vector2 &v1, const Vector2 &v2);
+	friend inline void operator /=(Vector2 &v1, const Vector2 &v2);
+	
+	// binary vector (op) scalar and scalar (op) vector operations
+	friend inline Vector2 operator +(const Vector2 &vec, scalar_t scalar);
+	friend inline Vector2 operator +(scalar_t scalar, const Vector2 &vec);
+	friend inline Vector2 operator -(const Vector2 &vec, scalar_t scalar);
+	friend inline Vector2 operator -(scalar_t scalar, const Vector2 &vec);
+	friend inline Vector2 operator *(const Vector2 &vec, scalar_t scalar);
+	friend inline Vector2 operator *(scalar_t scalar, const Vector2 &vec);
+	friend inline Vector2 operator /(const Vector2 &vec, scalar_t scalar);
+	friend inline Vector2 operator /(scalar_t scalar, const Vector2 &vec);
+	
+	friend inline void operator +=(Vector2 &vec, scalar_t scalar);
+	friend inline void operator -=(Vector2 &vec, scalar_t scalar);
+	friend inline void operator *=(Vector2 &vec, scalar_t scalar);
+	friend inline void operator /=(Vector2 &vec, scalar_t scalar);
+	
+	inline scalar_t length() const;
+	inline scalar_t length_sq() const;
+	void normalize();
+	Vector2 normalized() const;
+
+	void transform(const Matrix3x3 &mat);
+	Vector2 transformed(const Matrix3x3 &mat) const;
+
+	void rotate(scalar_t angle);
+	Vector2 rotated(scalar_t angle) const;
+			
+	Vector2 reflection(const Vector2 &normal) const;
+	Vector2 refraction(const Vector2 &normal, scalar_t src_ior, scalar_t dst_ior) const;
+	
+	friend std::ostream &operator <<(std::ostream &out, const Vector2 &vec);
+};
+
+/** 3D Vector */
+class Vector3 {
+public:
+	scalar_t x, y, z;
+
+	Vector3(scalar_t x = 0.0, scalar_t y = 0.0, scalar_t z = 0.0);
+	Vector3(const Vector2 &vec);
+	Vector3(const Vector4 &vec);
+	Vector3(const SphVector &sph);
+
+	Vector3 &operator =(const SphVector &sph);
+	
+	inline scalar_t &operator [](int elem);	
+	
+	// unary operations
+	friend inline Vector3 operator -(const Vector3 &vec);
+	
+	// binary vector (op) vector operations
+	friend inline scalar_t dot_product(const Vector3 &v1, const Vector3 &v2);
+	friend inline Vector3 cross_product(const Vector3 &v1, const Vector3 &v2);
+	
+	friend inline Vector3 operator +(const Vector3 &v1, const Vector3 &v2);
+	friend inline Vector3 operator -(const Vector3 &v1, const Vector3 &v2);	
+	friend inline Vector3 operator *(const Vector3 &v1, const Vector3 &v2);	
+	friend inline Vector3 operator /(const Vector3 &v1, const Vector3 &v2);	
+	friend inline bool operator ==(const Vector3 &v1, const Vector3 &v2);
+	
+	friend inline void operator +=(Vector3 &v1, const Vector3 &v2);
+	friend inline void operator -=(Vector3 &v1, const Vector3 &v2);
+	friend inline void operator *=(Vector3 &v1, const Vector3 &v2);
+	friend inline void operator /=(Vector3 &v1, const Vector3 &v2);
+	
+	// binary vector (op) scalar and scalar (op) vector operations
+	friend inline Vector3 operator +(const Vector3 &vec, scalar_t scalar);
+	friend inline Vector3 operator +(scalar_t scalar, const Vector3 &vec);
+	friend inline Vector3 operator -(const Vector3 &vec, scalar_t scalar);
+	friend inline Vector3 operator -(scalar_t scalar, const Vector3 &vec);
+	friend inline Vector3 operator *(const Vector3 &vec, scalar_t scalar);
+	friend inline Vector3 operator *(scalar_t scalar, const Vector3 &vec);
+	friend inline Vector3 operator /(const Vector3 &vec, scalar_t scalar);
+	friend inline Vector3 operator /(scalar_t scalar, const Vector3 &vec);
+	
+	friend inline void operator +=(Vector3 &vec, scalar_t scalar);
+	friend inline void operator -=(Vector3 &vec, scalar_t scalar);
+	friend inline void operator *=(Vector3 &vec, scalar_t scalar);
+	friend inline void operator /=(Vector3 &vec, scalar_t scalar);
+	
+	inline scalar_t length() const;
+	inline scalar_t length_sq() const;
+	void normalize();
+	Vector3 normalized() const;
+	
+	void transform(const Matrix3x3 &mat);
+	Vector3 transformed(const Matrix3x3 &mat) const;
+	void transform(const Matrix4x4 &mat);
+	Vector3 transformed(const Matrix4x4 &mat) const;
+	void transform(const Quaternion &quat);
+	Vector3 transformed(const Quaternion &quat) const;
+
+	void rotate(const Vector3 &euler);
+	Vector3 rotated(const Vector3 &euler) const;
+	
+	Vector3 reflection(const Vector3 &normal) const;
+	Vector3 refraction(const Vector3 &normal, scalar_t src_ior, scalar_t dst_ior) const;
+	
+	friend std::ostream &operator <<(std::ostream &out, const Vector3 &vec);
+};
+
+/** 4D Vector */
+class Vector4 {
+public:
+	scalar_t x, y, z, w;
+
+	Vector4(scalar_t x = 0.0, scalar_t y = 0.0, scalar_t z = 0.0, scalar_t w = 0.0);
+	Vector4(const Vector2 &vec);
+	Vector4(const Vector3 &vec);
+	
+	inline scalar_t &operator [](int elem);	
+	
+	// unary operations
+	friend inline Vector4 operator -(const Vector4 &vec);
+		
+	// binary vector (op) vector operations
+	friend inline scalar_t dot_product(const Vector4 &v1, const Vector4 &v2);
+	friend inline Vector4 cross_product(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3);
+	
+	friend inline Vector4 operator +(const Vector4 &v1, const Vector4 &v2);
+	friend inline Vector4 operator -(const Vector4 &v1, const Vector4 &v2);	
+	friend inline Vector4 operator *(const Vector4 &v1, const Vector4 &v2);	
+	friend inline Vector4 operator /(const Vector4 &v1, const Vector4 &v2);	
+	friend inline bool operator ==(const Vector4 &v1, const Vector4 &v2);
+	
+	friend inline void operator +=(Vector4 &v1, const Vector4 &v2);
+	friend inline void operator -=(Vector4 &v1, const Vector4 &v2);
+	friend inline void operator *=(Vector4 &v1, const Vector4 &v2);
+	friend inline void operator /=(Vector4 &v1, const Vector4 &v2);
+	
+	// binary vector (op) scalar and scalar (op) vector operations
+	friend inline Vector4 operator +(const Vector4 &vec, scalar_t scalar);
+	friend inline Vector4 operator +(scalar_t scalar, const Vector4 &vec);
+	friend inline Vector4 operator -(const Vector4 &vec, scalar_t scalar);
+	friend inline Vector4 operator -(scalar_t scalar, const Vector4 &vec);
+	friend inline Vector4 operator *(const Vector4 &vec, scalar_t scalar);
+	friend inline Vector4 operator *(scalar_t scalar, const Vector4 &vec);
+	friend inline Vector4 operator /(const Vector4 &vec, scalar_t scalar);
+	friend inline Vector4 operator /(scalar_t scalar, const Vector4 &vec);
+	
+	friend inline void operator +=(Vector4 &vec, scalar_t scalar);
+	friend inline void operator -=(Vector4 &vec, scalar_t scalar);
+	friend inline void operator *=(Vector4 &vec, scalar_t scalar);
+	friend inline void operator /=(Vector4 &vec, scalar_t scalar);
+	
+	inline scalar_t length() const;
+	inline scalar_t length_sq() const;
+	void normalize();
+	Vector4 normalized() const;
+
+	void transform(const Matrix4x4 &mat);
+	Vector4 transformed(const Matrix4x4 &mat) const;
+		
+	Vector4 reflection(const Vector4 &normal) const;
+	Vector4 refraction(const Vector4 &normal, scalar_t src_ior, scalar_t dst_ior) const;
+	
+	friend std::ostream &operator <<(std::ostream &out, const Vector4 &vec);
+};
+
+#include "vector.inl"
+
+#endif	// VECTOR_H_
diff --git a/tools/curve_draw/vmath/vector.inl b/tools/curve_draw/vmath/vector.inl
new file mode 100644
index 0000000..8852223
--- /dev/null
+++ b/tools/curve_draw/vmath/vector.inl
@@ -0,0 +1,427 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#include <cmath>
+#include "vmath.h"
+
+/**
+ * @file vector.inl
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * Inline functions of vector classes.
+ */
+
+
+inline scalar_t &Vector2::operator [](int elem) {
+	return elem ? y : x;
+}
+
+inline Vector2 operator -(const Vector2 &vec) {
+	return Vector2(-vec.x, -vec.y);
+}
+
+inline scalar_t dot_product(const Vector2 &v1, const Vector2 &v2) {
+	return v1.x * v2.x + v1.y * v2.y;
+}
+
+inline Vector2 operator +(const Vector2 &v1, const Vector2 &v2) {
+	return Vector2(v1.x + v2.x, v1.y + v2.y);
+}
+
+inline Vector2 operator -(const Vector2 &v1, const Vector2 &v2) {
+	return Vector2(v1.x - v2.x, v1.y - v2.y);
+}
+
+inline Vector2 operator *(const Vector2 &v1, const Vector2 &v2) {
+	return Vector2(v1.x * v2.x, v1.y * v2.y);
+}
+
+inline Vector2 operator /(const Vector2 &v1, const Vector2 &v2) {
+	return Vector2(v1.x / v2.x, v1.y / v2.y);
+}
+
+inline bool operator ==(const Vector2 &v1, const Vector2 &v2) {
+	return (fabs(v1.x - v2.x) < xsmall_number) && (fabs(v1.y - v2.x) < xsmall_number);
+}
+
+inline void operator +=(Vector2 &v1, const Vector2 &v2) {
+	v1.x += v2.x;
+	v1.y += v2.y;
+}
+
+inline void operator -=(Vector2 &v1, const Vector2 &v2) {
+	v1.x -= v2.x;
+	v1.y -= v2.y;
+}
+
+inline void operator *=(Vector2 &v1, const Vector2 &v2) {
+	v1.x *= v2.x;
+	v1.y *= v2.y;
+}
+
+inline void operator /=(Vector2 &v1, const Vector2 &v2) {
+	v1.x /= v2.x;
+	v1.y /= v2.y;
+}
+
+inline Vector2 operator +(const Vector2 &vec, scalar_t scalar) {
+	return Vector2(vec.x + scalar, vec.y + scalar);
+}
+
+inline Vector2 operator +(scalar_t scalar, const Vector2 &vec) {
+	return Vector2(vec.x + scalar, vec.y + scalar);
+}
+
+inline Vector2 operator -(const Vector2 &vec, scalar_t scalar) {
+	return Vector2(vec.x - scalar, vec.y - scalar);
+}
+
+inline Vector2 operator -(scalar_t scalar, const Vector2 &vec) {
+	return Vector2(vec.x - scalar, vec.y - scalar);
+}
+
+inline Vector2 operator *(const Vector2 &vec, scalar_t scalar) {
+	return Vector2(vec.x * scalar, vec.y * scalar);
+}
+
+inline Vector2 operator *(scalar_t scalar, const Vector2 &vec) {
+	return Vector2(vec.x * scalar, vec.y * scalar);
+}
+
+inline Vector2 operator /(const Vector2 &vec, scalar_t scalar) {
+	return Vector2(vec.x / scalar, vec.y / scalar);
+}
+
+inline Vector2 operator /(scalar_t scalar, const Vector2 &vec) {
+	return Vector2(vec.x / scalar, vec.y / scalar);
+}
+
+inline void operator +=(Vector2 &vec, scalar_t scalar) {
+	vec.x += scalar;
+	vec.y += scalar;
+}
+
+inline void operator -=(Vector2 &vec, scalar_t scalar) {
+	vec.x -= scalar;
+	vec.y -= scalar;
+}
+
+inline void operator *=(Vector2 &vec, scalar_t scalar) {
+	vec.x *= scalar;
+	vec.y *= scalar;
+}
+
+inline void operator /=(Vector2 &vec, scalar_t scalar) {
+	vec.x /= scalar;
+	vec.y /= scalar;
+}
+
+inline scalar_t Vector2::length() const {
+	return sqrt(x*x + y*y);
+}
+
+inline scalar_t Vector2::length_sq() const {
+	return x*x + y*y;
+}
+
+
+
+// ------------- Vector3 --------------
+
+inline scalar_t &Vector3::operator [](int elem) {
+	return elem ? (elem == 1 ? y : z) : x;
+}
+
+// unary operations
+inline Vector3 operator -(const Vector3 &vec) {
+	return Vector3(-vec.x, -vec.y, -vec.z);
+}
+
+// binary vector (op) vector operations
+inline scalar_t dot_product(const Vector3 &v1, const Vector3 &v2) {
+	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
+}
+
+inline Vector3 cross_product(const Vector3 &v1, const Vector3 &v2) {
+	return Vector3(v1.y * v2.z - v1.z * v2.y,  v1.z * v2.x - v1.x * v2.z,  v1.x * v2.y - v1.y * v2.x);
+}
+
+
+inline Vector3 operator +(const Vector3 &v1, const Vector3 &v2) {
+	return Vector3(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
+}
+
+inline Vector3 operator -(const Vector3 &v1, const Vector3 &v2) {
+	return Vector3(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
+}
+
+inline Vector3 operator *(const Vector3 &v1, const Vector3 &v2) {
+	return Vector3(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z);
+}
+
+inline Vector3 operator /(const Vector3 &v1, const Vector3 &v2) {
+	return Vector3(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z);
+}
+
+inline bool operator ==(const Vector3 &v1, const Vector3 &v2) {
+	return (fabs(v1.x - v2.x) < xsmall_number) && (fabs(v1.y - v2.y) < xsmall_number) && (fabs(v1.z - v2.z) < xsmall_number);
+}
+
+inline void operator +=(Vector3 &v1, const Vector3 &v2) {
+	v1.x += v2.x;
+	v1.y += v2.y;
+	v1.z += v2.z;
+}
+
+inline void operator -=(Vector3 &v1, const Vector3 &v2) {
+	v1.x -= v2.x;
+	v1.y -= v2.y;
+	v1.z -= v2.z;
+}
+
+inline void operator *=(Vector3 &v1, const Vector3 &v2) {
+	v1.x *= v2.x;
+	v1.y *= v2.y;
+	v1.z *= v2.z;
+}
+
+inline void operator /=(Vector3 &v1, const Vector3 &v2) {
+	v1.x /= v2.x;
+	v1.y /= v2.y;
+	v1.z /= v2.z;
+}
+// binary vector (op) scalar and scalar (op) vector operations
+inline Vector3 operator +(const Vector3 &vec, scalar_t scalar) {
+	return Vector3(vec.x + scalar, vec.y + scalar, vec.z + scalar);
+}
+
+inline Vector3 operator +(scalar_t scalar, const Vector3 &vec) {
+	return Vector3(vec.x + scalar, vec.y + scalar, vec.z + scalar);
+}
+
+inline Vector3 operator -(const Vector3 &vec, scalar_t scalar) {
+	return Vector3(vec.x - scalar, vec.y - scalar, vec.z - scalar);
+}
+
+inline Vector3 operator -(scalar_t scalar, const Vector3 &vec) {
+	return Vector3(vec.x - scalar, vec.y - scalar, vec.z - scalar);
+}
+
+inline Vector3 operator *(const Vector3 &vec, scalar_t scalar) {
+	return Vector3(vec.x * scalar, vec.y * scalar, vec.z * scalar);
+}
+
+inline Vector3 operator *(scalar_t scalar, const Vector3 &vec) {
+	return Vector3(vec.x * scalar, vec.y * scalar, vec.z * scalar);
+}
+
+inline Vector3 operator /(const Vector3 &vec, scalar_t scalar) {
+	return Vector3(vec.x / scalar, vec.y / scalar, vec.z / scalar);
+}
+
+inline Vector3 operator /(scalar_t scalar, const Vector3 &vec) {
+	return Vector3(vec.x / scalar, vec.y / scalar, vec.z / scalar);
+}
+
+inline void operator +=(Vector3 &vec, scalar_t scalar) {
+	vec.x += scalar;
+	vec.y += scalar;
+	vec.z += scalar;
+}
+
+inline void operator -=(Vector3 &vec, scalar_t scalar) {
+	vec.x -= scalar;
+	vec.y -= scalar;
+	vec.z -= scalar;
+}
+
+inline void operator *=(Vector3 &vec, scalar_t scalar) {
+	vec.x *= scalar;
+	vec.y *= scalar;
+	vec.z *= scalar;
+}
+
+inline void operator /=(Vector3 &vec, scalar_t scalar) {
+	vec.x /= scalar;
+	vec.y /= scalar;
+	vec.z /= scalar;
+}
+
+inline scalar_t Vector3::length() const {
+	return sqrt(x*x + y*y + z*z);
+}
+inline scalar_t Vector3::length_sq() const {
+	return x*x + y*y + z*z;
+}
+
+
+
+// ----------- Vector4 -----------------
+
+inline scalar_t &Vector4::operator [](int elem) {
+	return elem ? (elem == 1 ? y : (elem == 2 ? z : w)) : x;
+}
+
+inline Vector4 operator -(const Vector4 &vec) {
+	return Vector4(-vec.x, -vec.y, -vec.z, -vec.w);
+}
+
+inline scalar_t dot_product(const Vector4 &v1, const Vector4 &v2) {
+	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z + v1.w * v2.w;
+}
+
+inline Vector4 cross_product(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3) {
+	float  A, B, C, D, E, F;       // Intermediate Values
+    Vector4 result;
+
+    // Calculate intermediate values.
+    A = (v2.x * v3.y) - (v2.y * v3.x);
+    B = (v2.x * v3.z) - (v2.z * v3.x);
+    C = (v2.x * v3.w) - (v2.w * v3.x);
+    D = (v2.y * v3.z) - (v2.z * v3.y);
+    E = (v2.y * v3.w) - (v2.w * v3.y);
+    F = (v2.z * v3.w) - (v2.w * v3.z);
+
+    // Calculate the result-vector components.
+    result.x =   (v1.y * F) - (v1.z * E) + (v1.w * D);
+    result.y = - (v1.x * F) + (v1.z * C) - (v1.w * B);
+    result.z =   (v1.x * E) - (v1.y * C) + (v1.w * A);
+    result.w = - (v1.x * D) + (v1.y * B) - (v1.z * A);
+    return result;
+}
+
+inline Vector4 operator +(const Vector4 &v1, const Vector4 &v2) {
+	return Vector4(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z, v1.w + v2.w);
+}
+
+inline Vector4 operator -(const Vector4 &v1, const Vector4 &v2) {
+	return Vector4(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.w - v2.w);
+}
+
+inline Vector4 operator *(const Vector4 &v1, const Vector4 &v2) {
+	return Vector4(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z, v1.w * v2.w);
+}
+
+inline Vector4 operator /(const Vector4 &v1, const Vector4 &v2) {
+	return Vector4(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z, v1.w / v2.w);
+}
+
+inline bool operator ==(const Vector4 &v1, const Vector4 &v2) {
+	return 	(fabs(v1.x - v2.x) < xsmall_number) && 
+			(fabs(v1.y - v2.y) < xsmall_number) && 
+			(fabs(v1.z - v2.z) < xsmall_number) &&
+			(fabs(v1.w - v2.w) < xsmall_number);
+}
+
+inline void operator +=(Vector4 &v1, const Vector4 &v2) {
+	v1.x += v2.x;
+	v1.y += v2.y;
+	v1.z += v2.z;
+	v1.w += v2.w;
+}
+
+inline void operator -=(Vector4 &v1, const Vector4 &v2) {
+	v1.x -= v2.x;
+	v1.y -= v2.y;
+	v1.z -= v2.z;
+	v1.w -= v2.w;
+}
+
+inline void operator *=(Vector4 &v1, const Vector4 &v2) {
+	v1.x *= v2.x;
+	v1.y *= v2.y;
+	v1.z *= v2.z;
+	v1.w *= v2.w;
+}
+
+inline void operator /=(Vector4 &v1, const Vector4 &v2) {
+	v1.x /= v2.x;
+	v1.y /= v2.y;
+	v1.z /= v2.z;
+	v1.w /= v2.w;
+}
+
+// binary vector (op) scalar and scalar (op) vector operations
+inline Vector4 operator +(const Vector4 &vec, scalar_t scalar) {
+	return Vector4(vec.x + scalar, vec.y + scalar, vec.z + scalar, vec.w + scalar);
+}
+
+inline Vector4 operator +(scalar_t scalar, const Vector4 &vec) {
+	return Vector4(vec.x + scalar, vec.y + scalar, vec.z + scalar, vec.w + scalar);
+}
+
+inline Vector4 operator -(const Vector4 &vec, scalar_t scalar) {
+	return Vector4(vec.x - scalar, vec.y - scalar, vec.z - scalar, vec.w - scalar);
+}
+
+inline Vector4 operator -(scalar_t scalar, const Vector4 &vec) {
+	return Vector4(vec.x - scalar, vec.y - scalar, vec.z - scalar, vec.w - scalar);
+}
+
+inline Vector4 operator *(const Vector4 &vec, scalar_t scalar) {
+	return Vector4(vec.x * scalar, vec.y * scalar, vec.z * scalar, vec.w * scalar);
+}
+
+inline Vector4 operator *(scalar_t scalar, const Vector4 &vec) {
+	return Vector4(vec.x * scalar, vec.y * scalar, vec.z * scalar, vec.w * scalar);
+}
+
+inline Vector4 operator /(const Vector4 &vec, scalar_t scalar) {
+	return Vector4(vec.x / scalar, vec.y / scalar, vec.z / scalar, vec.w / scalar);
+}
+
+inline Vector4 operator /(scalar_t scalar, const Vector4 &vec) {
+	return Vector4(vec.x / scalar, vec.y / scalar, vec.z / scalar, vec.w / scalar);
+}
+
+inline void operator +=(Vector4 &vec, scalar_t scalar) {
+	vec.x += scalar;
+	vec.y += scalar;
+	vec.z += scalar;
+	vec.w += scalar;
+}
+
+inline void operator -=(Vector4 &vec, scalar_t scalar) {
+	vec.x -= scalar;
+	vec.y -= scalar;
+	vec.z -= scalar;
+	vec.w -= scalar;
+}
+
+inline void operator *=(Vector4 &vec, scalar_t scalar) {
+	vec.x *= scalar;
+	vec.y *= scalar;
+	vec.z *= scalar;
+	vec.w *= scalar;
+}
+
+inline void operator /=(Vector4 &vec, scalar_t scalar) {
+	vec.x /= scalar;
+	vec.y /= scalar;
+	vec.z /= scalar;
+	vec.w /= scalar;
+}
+
+inline scalar_t Vector4::length() const {
+	return sqrt(x*x + y*y + z*z + w*w);
+}
+inline scalar_t Vector4::length_sq() const {
+	return x*x + y*y + z*z + w*w;
+}
diff --git a/tools/curve_draw/vmath/vmath.cc b/tools/curve_draw/vmath/vmath.cc
new file mode 100644
index 0000000..c856538
--- /dev/null
+++ b/tools/curve_draw/vmath/vmath.cc
@@ -0,0 +1,211 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/**
+ * @file vmath.cxx
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * vmath provides the mathematical foundations for 3D graphics.
+ * This specific file contains various constants and basic operations
+ * as well as the Basis class, which represents an orthogonal basis in
+ * 3 dimensional space.
+ */
+
+#include "vmath.h"
+
+const scalar_t e = 2.7182818;
+
+const scalar_t pi = 3.1415926535897932;
+const scalar_t half_pi = pi / 2.0;
+const scalar_t quarter_pi = pi / 4.0;
+const scalar_t two_pi = pi * 2.0;
+const scalar_t three_half_pi = 3.0 * pi / 2.0;
+
+const scalar_t xsmall_number = 1.e-8;
+const scalar_t small_number = 1.e-4;
+
+const scalar_t error_margin = 1.e-6;
+
+
+/** Generates a random number in [0, 1) */
+scalar_t frand(scalar_t range) {
+	return range * (float)rand() / (float)RAND_MAX;
+}
+
+/** Numerical calculation of integrals using simpson's rule */
+scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples)  {
+	scalar_t h = (high - low) / (scalar_t)samples;
+	scalar_t sum = 0.0;
+	
+	for(int i=0; i<samples+1; i++) {
+		scalar_t y = f((scalar_t)i * h + low);
+		sum += ((!i || i == samples) ? y : ((i%2) ? 4.0 * y : 2.0 * y)) * (h / 3.0);
+	}
+	return sum;
+}
+
+/** Gaussuan function */
+scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev) {
+	scalar_t exponent = -SQ(x - mean) / (2.0 * SQ(sdev));
+	
+	return 1.0 - -pow(e, exponent) / (sdev * sqrt(two_pi));
+}
+
+
+/** b-spline approximation */
+scalar_t bspline(const Vector4 &cpvec, scalar_t t) {
+	Matrix4x4 bspline_mat(	-1,  3, -3,  1,
+							 3, -6,  3,  0,
+							-3,  0,  3,  0,
+							 1,  4,  1,  0);
+	
+	scalar_t t_square = t * t;
+	scalar_t t_cube = t_square * t;
+	Vector4 params(t_cube, t_square, t, 1.0);
+
+	return dot_product(params, cpvec.transformed(bspline_mat) / 6.0);
+}
+
+/** Catmull-rom spline interpolation */
+scalar_t catmull_rom_spline(const Vector4 &cpvec, scalar_t t) {
+	Matrix4x4 crspline_mat(	-1,  3, -3,  1,
+							 2, -5,  4, -1,
+							-1,  0,  1,  0,
+							 0,  2,  0,  0);
+
+	scalar_t t_square = t * t;
+	scalar_t t_cube = t_square * t;
+	Vector4 params(t_cube, t_square, t, 1.0);
+
+	return dot_product(params, cpvec.transformed(crspline_mat) / 2.0);
+}
+
+/**
+ * @author Mihalis Georgoulopoulos
+ * 
+ * Bezier interpolation.
+ * returns an interpolated scalar value given 4 control values.
+ */
+scalar_t bezier(const Vector4 &cp, scalar_t t)
+{
+	static scalar_t omt, omt3, t3, f;
+	t3 = t * t * t;
+	omt = 1.0f - t;
+	omt3 = omt * omt * omt;
+	f = 3 * t * omt;
+
+	return (cp.x * omt3) + (cp.y * f * omt) + (cp.z * f * t) + (cp.w * t3);
+}
+
+/**
+ * @author Mihalis Georgoulopoulos
+ * 
+ * Bezier interpolation.
+ * Vector3 version of the bezier function.
+ */
+Vector3 bezier(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2, const Vector3 &p3, scalar_t t)
+{
+	static scalar_t omt, omt3, t3, f;
+	t3 = t * t * t;
+	omt = 1.0f - t;
+	omt3 = omt * omt * omt;
+	f = 3 * t * omt;
+
+	return (p0 * omt3) + (p1 * f * omt) + (p2 * f * t) + (p3 * t3);
+}
+
+/**
+ * @author Mihalis Georgoulopoulos
+ *
+ * Bezier tangent calculation.
+ * returns a vector tangent to the bezier curve at the specified point.
+ */
+Vector3 bezier_tangent(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2, const Vector3 &p3, scalar_t t)
+{
+	static scalar_t omt;
+	omt = 1.0f - t;
+	static Vector3 p4, p5, p6, p7, p8;
+	p4 = p0 * omt + p1 * t;
+	p5 = p1 * omt + p2 * t;
+	p6 = p2 * omt + p3 * t;
+
+	p7 = p4 * omt + p5 * t;
+	p8 = p5 * omt + p6 * t;
+
+	return p8 - p7;
+}
+
+Basis::Basis() {
+	i = Vector3(1, 0, 0);
+	j = Vector3(0, 1, 0);
+	k = Vector3(0, 0, 1);
+}
+
+Basis::Basis(const Vector3 &i, const Vector3 &j, const Vector3 &k) {
+	this->i = i;
+	this->j = j;
+	this->k = k;
+}
+
+Basis::Basis(const Vector3 &dir, bool left_handed) {
+	k = dir;
+	j = Vector3(0, 1, 0);
+	i = cross_product(j, k);
+	j = cross_product(k, i);
+}
+
+/** Rotate with euler angles */
+void Basis::rotate(scalar_t x, scalar_t y, scalar_t z) {
+	Matrix4x4 RotMat;
+	RotMat.set_rotation(Vector3(x, y, z));
+	i.transform(RotMat);
+	j.transform(RotMat);
+	k.transform(RotMat);
+}
+
+/** Rotate by axis-angle specification */
+void Basis::rotate(const Vector3 &axis, scalar_t angle) {
+	Quaternion q;
+	q.set_rotation(axis, angle);
+	i.transform(q);
+	j.transform(q);
+	k.transform(q);
+}
+
+/** Rotate with a 4x4 matrix */
+void Basis::rotate(const Matrix4x4 &mat) {
+	i.transform(mat);
+	j.transform(mat);
+	k.transform(mat);
+}
+
+/** Rotate with a quaternion */
+void Basis::rotate(const Quaternion &quat) {
+	i.transform(quat);
+	j.transform(quat);
+	k.transform(quat);
+}
+
+/** Extract a rotation matrix from the orthogonal basis */
+Matrix3x3 Basis::create_rotation_matrix() const {
+	return Matrix3x3(	i.x, j.x, k.x,
+						i.y, j.y, k.y,
+						i.z, j.z, k.z);
+}
diff --git a/tools/curve_draw/vmath/vmath.h b/tools/curve_draw/vmath/vmath.h
new file mode 100644
index 0000000..3349821
--- /dev/null
+++ b/tools/curve_draw/vmath/vmath.h
@@ -0,0 +1,103 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/** 
+ * @file vmath.h
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * vmath provides the mathematical foundations for 3D graphics.
+ * This specific file contains various constants and basic operations
+ * as well as the Basis class, which represents an orthogonal basis in
+ * 3 dimensional space.
+ */
+
+#ifndef VMATH_H_
+#define VMATH_H_
+
+#include "vmath_types.h"
+
+extern const scalar_t e;
+
+extern const scalar_t pi;
+extern const scalar_t half_pi;
+extern const scalar_t quarter_pi;
+extern const scalar_t two_pi;
+extern const scalar_t three_half_pi;
+
+extern const scalar_t xsmall_number;
+extern const scalar_t small_number;
+
+extern const scalar_t error_margin;
+
+
+// angle conversion
+#define RAD_TO_DEG(a) ((((scalar_t)a) * 360.0) / two_pi)
+#define DEG_TO_RAD(a) (((scalar_t)a) * (pi / 180.0))
+
+#define SQ(x) ((x) * (x))
+
+#ifndef __GNUC__
+#define round(x)	((x) >= 0 ? (x) + 0.5 : (x) - 0.5)
+#endif	// __GNUC__
+
+// -- mathematical & helper functions --
+scalar_t frand(scalar_t range);
+scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples);
+scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev);
+
+// -- interpolation and approximation --
+inline scalar_t lerp(scalar_t a, scalar_t b, scalar_t t);
+
+scalar_t bspline(const Vector4 &cpvec, scalar_t t);
+inline scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
+
+scalar_t catmull_rom_spline(const Vector4 &cpvec, scalar_t t);
+inline scalar_t catmull_rom_spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
+
+scalar_t bezier(const Vector4 &cp, scalar_t t);
+Vector3 bezier(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2, const Vector3 &p3, scalar_t t);
+Vector3 bezier_tangent(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2, const Vector3 &p3, scalar_t t);
+
+// -- actual class definitions --
+#include "vector.h"
+#include "matrix.h"
+#include "quaternion.h"
+#include "sphvec.h"
+
+/** Orthogonal basis in 3-space */
+class Basis {
+public:
+	Vector3 i, j, k;
+
+	Basis();
+	Basis(const Vector3 &i, const Vector3 &j, const Vector3 &k);
+	Basis(const Vector3 &dir, bool left_handed = true);
+
+	void rotate(scalar_t x, scalar_t y, scalar_t z);
+	void rotate(const Vector3 &axis, scalar_t angle);
+	void rotate(const Matrix4x4 &mat);
+	void rotate(const Quaternion &quat);
+
+	Matrix3x3 create_rotation_matrix() const;
+};
+
+#include "vmath.inl"
+
+#endif	// VMATH_H_
diff --git a/tools/curve_draw/vmath/vmath.inl b/tools/curve_draw/vmath/vmath.inl
new file mode 100644
index 0000000..687ad14
--- /dev/null
+++ b/tools/curve_draw/vmath/vmath.inl
@@ -0,0 +1,41 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+/** 
+ * @file vmath.inl
+ * @author John Tsiombikas
+ * @date 28 June 2005
+ *
+ * vmath inline functions.
+ */
+
+/** linear interpolation */
+inline scalar_t lerp(scalar_t a, scalar_t b, scalar_t t) {
+	return a + (b - a) * t;
+}
+
+/** b-spline approximation */
+inline scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t) {
+	return bspline(Vector4(a, b, c, d), t);
+}
+
+/** catmull-rom spline interpolation */
+inline scalar_t catmull_rom_spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t) {
+	return catmull_rom_spline(Vector4(a, b, c, d), t);
+}
diff --git a/tools/curve_draw/vmath/vmath_types.h b/tools/curve_draw/vmath/vmath_types.h
new file mode 100644
index 0000000..74fb48d
--- /dev/null
+++ b/tools/curve_draw/vmath/vmath_types.h
@@ -0,0 +1,36 @@
+/*
+This file is part of XRay, a photorealistic 3D rendering library.
+Copyright (C) 2005 John Tsiombikas
+
+XRay is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+XRay is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with XRay; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef VMATH_TYPES_H_
+#define VMATH_TYPES_H_
+
+//#include "types.h"
+typedef float scalar_t;
+
+// forward declarations of all vmath classes
+class Vector2;
+class Vector3;
+class Vector4;
+class SphVector;
+class Matrix3x3;
+class Matrix4x4;
+class Quaternion;
+class Basis;
+
+#endif	// VMATH_TYPES_H_
diff --git a/tools/find_textures b/tools/find_textures
new file mode 100755
index 0000000..57024b7
--- /dev/null
+++ b/tools/find_textures
@@ -0,0 +1,12 @@
+#!/bin/sh
+
+tmpfile=/tmp/find_tex_tmp_file
+files=`find . -name '*.cpp'`
+
+cat /dev/null >$tmpfile
+for i in $files; do
+	grep '".*"' $i | egrep '(jpg|png)' | \
+	sed 's/.*[(,\, ,\t]"//g' | sed 's/"[),\,,].*//g' >>$tmpfile
+done
+
+cat $tmpfile | sort | uniq