# HG changeset patch
# User John Tsiombikas <nuclear@member.fsf.org>
# Date 1380432019 -10800
# Node ID 4deb0b12fe14de29b05220f321103c7d8db802d8
# Parent  1c601bf07b863f0d16dce45c607bde4db4172e68
wtf... corrupted heap?

diff -r 1c601bf07b86 -r 4deb0b12fe14 exporters/maxgoat/maxgoat.vcxproj
--- a/exporters/maxgoat/maxgoat.vcxproj	Sat Sep 28 20:36:55 2013 +0300
+++ b/exporters/maxgoat/maxgoat.vcxproj	Sun Sep 29 08:20:19 2013 +0300
@@ -33,7 +33,7 @@
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
     <ConfigurationType>DynamicLibrary</ConfigurationType>
     <UseDebugLibraries>true</UseDebugLibraries>
-    <PlatformToolset>v110</PlatformToolset>
+    <PlatformToolset>v100</PlatformToolset>
     <CharacterSet>Unicode</CharacterSet>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
@@ -46,7 +46,7 @@
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
     <ConfigurationType>DynamicLibrary</ConfigurationType>
     <UseDebugLibraries>false</UseDebugLibraries>
-    <PlatformToolset>v110</PlatformToolset>
+    <PlatformToolset>v100</PlatformToolset>
     <WholeProgramOptimization>true</WholeProgramOptimization>
     <CharacterSet>Unicode</CharacterSet>
   </PropertyGroup>
@@ -70,7 +70,7 @@
     <LinkIncremental>true</LinkIncremental>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
-    <LinkIncremental>true</LinkIncremental>
+    <LinkIncremental>false</LinkIncremental>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
     <LinkIncremental>false</LinkIncremental>
@@ -98,14 +98,14 @@
       </PrecompiledHeader>
       <WarningLevel>Level3</WarningLevel>
       <Optimization>Disabled</Optimization>
-      <PreprocessorDefinitions>WIN32;_DEBUG;_WINDOWS;_USRDLL;MAXGOAT_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
-      <DisableSpecificWarnings>4996</DisableSpecificWarnings>
+      <PreprocessorDefinitions>WIN32;_DEBUG;_WINDOWS;_USRDLL;MAXGOAT_EXPORTS;%(PreprocessorDefinitions);BUILD_MAXPLUGIN</PreprocessorDefinitions>
+      <DisableSpecificWarnings>4996;4244</DisableSpecificWarnings>
       <AdditionalIncludeDirectories>$(SolutionDir)\src</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
       <GenerateDebugInformation>true</GenerateDebugInformation>
-      <AdditionalDependencies>goat3d-x64.lib;libvmath-x64.lib;libanim-x64.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>goat3d-x64.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <AdditionalLibraryDirectories>$(TargetDir);%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
     </Link>
     <PostBuildEvent>
@@ -139,7 +139,7 @@
       <FunctionLevelLinking>true</FunctionLevelLinking>
       <IntrinsicFunctions>true</IntrinsicFunctions>
       <PreprocessorDefinitions>WIN32;NDEBUG;_WINDOWS;_USRDLL;MAXGOAT_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
-      <DisableSpecificWarnings>4996</DisableSpecificWarnings>
+      <DisableSpecificWarnings>4996;4244</DisableSpecificWarnings>
       <AdditionalIncludeDirectories>$(SolutionDir)\src</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
@@ -147,7 +147,7 @@
       <GenerateDebugInformation>true</GenerateDebugInformation>
       <EnableCOMDATFolding>true</EnableCOMDATFolding>
       <OptimizeReferences>true</OptimizeReferences>
-      <AdditionalDependencies>goat3d-x64.lib;libvmath-x64.lib;libanim-x64.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>goat3d-x64.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <AdditionalLibraryDirectories>$(TargetDir);%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
     </Link>
     <PostBuildEvent>
@@ -157,9 +157,6 @@
   <ItemGroup>
     <ClCompile Include="src\maxgoat.cc" />
   </ItemGroup>
-  <ItemGroup>
-    <ClInclude Include="src\config.h" />
-  </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
   </ImportGroup>
diff -r 1c601bf07b86 -r 4deb0b12fe14 exporters/maxgoat/maxgoat.vcxproj.filters
--- a/exporters/maxgoat/maxgoat.vcxproj.filters	Sat Sep 28 20:36:55 2013 +0300
+++ b/exporters/maxgoat/maxgoat.vcxproj.filters	Sun Sep 29 08:20:19 2013 +0300
@@ -11,9 +11,4 @@
       <Filter>src</Filter>
     </ClCompile>
   </ItemGroup>
-  <ItemGroup>
-    <ClInclude Include="src\config.h">
-      <Filter>src</Filter>
-    </ClInclude>
-  </ItemGroup>
 </Project>
\ No newline at end of file
diff -r 1c601bf07b86 -r 4deb0b12fe14 exporters/maxgoat/src/maxgoat.cc
--- a/exporters/maxgoat/src/maxgoat.cc	Sat Sep 28 20:36:55 2013 +0300
+++ b/exporters/maxgoat/src/maxgoat.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -1,6 +1,8 @@
 #include <stdio.h>
 #include <string.h>
+#include <stdlib.h>
 #include <errno.h>
+#include <map>
 #include <windows.h>
 #include <shlobj.h>
 #include "max.h"
@@ -33,6 +35,10 @@
 static HINSTANCE hinst;
 
 class GoatExporter : public SceneExport {
+private:
+	//std::map<IGameMaterial*, goat3d_material*> mtlmap;
+	//std::map<IGameNode*, goat3d_node*> nodemap;
+
 public:
 	IGameScene *igame;
 
@@ -49,9 +55,13 @@
 
 	int DoExport(const MCHAR *name, ExpInterface *eiface, Interface *iface, BOOL silent = FALSE, DWORD opt = 0);
 
-	void export_materials(goat3d *goat);
-	void export_meshes(goat3d *goat);
-	void process_node(goat3d *goat, IGameNode *maxnode);
+	void process_materials(goat3d *goat);
+
+	void process_node(goat3d *goat, goat3d_node *parent, IGameNode *maxnode);
+
+	void process_mesh(goat3d *goat, goat3d_mesh *mesh, IGameObject *maxobj);
+	void process_light(goat3d *goat, goat3d_light *light, IGameObject *maxobj);
+	void process_camera(goat3d *goat, goat3d_camera *cam, IGameObject *maxobj);
 };
 
 
@@ -62,7 +72,7 @@
 
 const TCHAR *GoatExporter::Ext(int n)
 {
-	return L"txt";
+	return L"xml";
 }
 
 const TCHAR *GoatExporter::LongDesc()
@@ -108,6 +118,9 @@
 int GoatExporter::DoExport(const MCHAR *name, ExpInterface *eiface, Interface *iface,
 		BOOL non_interactive, DWORD opt)
 {
+	//mtlmap.clear();
+	//nodemap.clear();
+
 	char fname[512];
 	wcstombs(fname, name, sizeof fname - 1);
 
@@ -122,20 +135,19 @@
 
 	goat3d *goat = goat3d_create();
 
-	export_materials(goat);
-	export_meshes(goat);
+	process_materials(goat);
+
+	// process all nodes
+	for(int i=0; i<igame->GetTopLevelNodeCount(); i++) {
+		IGameNode *node = igame->GetTopLevelNode(i);
+		process_node(goat, 0, node);
+	}
 
 	if(goat3d_save(goat, fname) == -1) {
 		goat3d_free(goat);
 		return IMPEXP_FAIL;
 	}
 
-	// process all nodes
-	for(int i=0; i<igame->GetTopLevelNodeCount(); i++) {
-		IGameNode *node = igame->GetTopLevelNode(i);
-		process_node(goat, node);
-	}
-
 	goat3d_free(goat);
 	return IMPEXP_SUCCESS;
 }
@@ -148,7 +160,7 @@
 	return str;
 }
 
-void GoatExporter::export_materials(goat3d *goat)
+void GoatExporter::process_materials(goat3d *goat)
 {
 	IGameProperty *prop;
 
@@ -231,19 +243,138 @@
 			}
 
 			goat3d_add_mtl(goat, mtl);
+			//mtlmap[maxmtl] = mtl;
 		}
 	}
 }
 
-void GoatExporter::export_meshes(goat3d *goat)
+void GoatExporter::process_node(goat3d *goat, goat3d_node *parent, IGameNode *maxnode)
 {
-	Tab<IGameNode*> meshes = igame->GetIGameNodeByType(IGameObject::IGAME_MESH);
+	goat3d_node *node = goat3d_create_node();
+	goat3d_add_node(goat, node);
 
-	for(int i=0; i<meshes.Count(); i++) {
-		const char *name = max_string(meshes[i]->GetName());
+	const char *name = max_string(maxnode->GetName());
+	if(name) {
+		goat3d_set_node_name(node, name);
+	}
+
+	// no animation yet, just get the static PRS
+	GMatrix maxmatrix = maxnode->GetObjectTM();
+	Point3 trans = maxmatrix.Translation();
+	Quat rot = maxmatrix.Rotation();
+	Point3 scale = maxmatrix.Scaling();
+
+	goat3d_set_node_position(node, trans.x, trans.y, trans.z, 0);
+	goat3d_set_node_rotation(node, rot.x, rot.y, rot.z, rot.w, 0);
+	goat3d_set_node_scaling(node, scale.x, scale.y, scale.z, 0);
+
+	IGameObject *maxobj = maxnode->GetIGameObject();
+	IGameObject::ObjectTypes type = maxobj->GetIGameType();
+
+	switch(type) {
+	case IGameObject::IGAME_MESH:
+		{
+			goat3d_mesh *mesh = goat3d_create_mesh();
+			if(name) goat3d_set_mesh_name(mesh, name);
+			goat3d_set_node_object(node, GOAT3D_NODE_MESH, mesh);
+
+			// get the node material and assign it to the mesh
+			IGameMaterial *maxmtl = maxnode->GetNodeMaterial();
+			goat3d_material *mtl = 0;//mtlmap[maxmtl];
+			if(mtl) {
+				goat3d_set_mesh_mtl(mesh, mtl);
+			}
+
+			process_mesh(goat, mesh, maxobj);
+		}
+		break;
+
+	case IGameObject::IGAME_LIGHT:
+		{
+			goat3d_light *light = goat3d_create_light();
+			//if(name) goat3d_set_light_name(light, name);
+			goat3d_set_node_object(node, GOAT3D_NODE_LIGHT, light);
+
+			process_light(goat, light, maxobj);
+		}
+		break;
+
+	case IGameObject::IGAME_CAMERA:
+		{
+			goat3d_camera *cam = goat3d_create_camera();
+			//if(name) goat3d_set_camera_name(camera, name);
+			goat3d_set_node_object(node, GOAT3D_NODE_CAMERA, cam);
+
+			process_camera(goat, cam, maxobj);
+		}
+		break;
+
+	default:
+		// otherwise don't assign an object, essentially treating it as a null node
+		break;
+	}
+		
+
+	for(int i=0; i<maxnode->GetChildCount(); i++) {
+		process_node(goat, node, maxnode->GetNodeChild(i));
 	}
 }
 
+void GoatExporter::process_mesh(goat3d *goat, goat3d_mesh *mesh, IGameObject *maxobj)
+{
+	IGameMesh *maxmesh = (IGameMesh*)maxobj;
+
+	maxmesh->SetCreateOptimizedNormalList();	// not needed any more according to docs
+	maxobj->InitializeData();
+
+	int num_verts = maxmesh->GetNumberOfVerts();
+	assert(maxmesh->GetNumberOfTexVerts() == num_verts);
+
+	float *vertices = new float[num_verts * 3];
+	float *normals = new float[num_verts * 3];
+	float *texcoords = new float[num_verts * 2];
+
+	for(int i=0; i<num_verts; i++) {
+		Point3 v = maxmesh->GetVertex(i, true);
+		vertices[i * 3] = v.x;
+		vertices[i * 3 + 1] = v.y;
+		vertices[i * 3 + 2] = v.z;
+	}
+
+	for(int i=0; i<maxmesh->GetNumberOfNormals(); i++) {
+		Point3 norm = maxmesh->GetNormal(i);
+
+		int vidx = maxmesh->GetNormalVertexIndex(i);
+		normals[vidx * 3] = norm.x;
+		normals[vidx * 3 + 1] = norm.y;
+		normals[vidx * 3 + 2] = norm.z;
+	}
+
+	for(int i=0; i<maxmesh->GetNumberOfTexVerts(); i++) {
+		Point3 tex = maxmesh->GetTexVertex(i);
+
+		texcoords[i * 2] = tex.x;
+		texcoords[i * 2 + 1] = tex.y;
+	}
+
+	goat3d_set_mesh_attribs(mesh, GOAT3D_MESH_ATTR_VERTEX, vertices, num_verts);
+	goat3d_set_mesh_attribs(mesh, GOAT3D_MESH_ATTR_NORMAL, normals, num_verts);
+	goat3d_set_mesh_attribs(mesh, GOAT3D_MESH_ATTR_TEXCOORD, texcoords, num_verts);
+
+	delete [] vertices;
+	delete [] normals;
+	delete [] texcoords;
+}
+
+void GoatExporter::process_light(goat3d *goat, goat3d_light *light, IGameObject *maxobj)
+{
+}
+
+void GoatExporter::process_camera(goat3d *goat, goat3d_camera *cam, IGameObject *maxobj)
+{
+}
+
+
 // ------------------------------------------
 
 class GoatClassDesc : public ClassDesc2 {
diff -r 1c601bf07b86 -r 4deb0b12fe14 goat3d.vcxproj
--- a/goat3d.vcxproj	Sat Sep 28 20:36:55 2013 +0300
+++ b/goat3d.vcxproj	Sun Sep 29 08:20:19 2013 +0300
@@ -19,6 +19,9 @@
     </ProjectConfiguration>
   </ItemGroup>
   <ItemGroup>
+    <ClInclude Include="libs\anim\anim.h" />
+    <ClInclude Include="libs\anim\dynarr.h" />
+    <ClInclude Include="libs\anim\track.h" />
     <ClInclude Include="libs\openctm\internal.h" />
     <ClInclude Include="libs\openctm\liblzma\Alloc.h" />
     <ClInclude Include="libs\openctm\liblzma\LzFind.h" />
@@ -31,6 +34,16 @@
     <ClInclude Include="libs\openctm\openctm.h" />
     <ClInclude Include="libs\openctm\openctmpp.h" />
     <ClInclude Include="libs\tinyxml2\tinyxml2.h" />
+    <ClInclude Include="libs\vmath\basis.h" />
+    <ClInclude Include="libs\vmath\geom.h" />
+    <ClInclude Include="libs\vmath\matrix.h" />
+    <ClInclude Include="libs\vmath\quat.h" />
+    <ClInclude Include="libs\vmath\ray.h" />
+    <ClInclude Include="libs\vmath\sphvec.h" />
+    <ClInclude Include="libs\vmath\vector.h" />
+    <ClInclude Include="libs\vmath\vmath.h" />
+    <ClInclude Include="libs\vmath\vmath_config.h" />
+    <ClInclude Include="libs\vmath\vmath_types.h" />
     <ClInclude Include="src\camera.h" />
     <ClInclude Include="src\chunk.h" />
     <ClInclude Include="src\goat3d.h" />
@@ -44,6 +57,9 @@
     <ClInclude Include="src\xform_node.h" />
   </ItemGroup>
   <ItemGroup>
+    <ClCompile Include="libs\anim\anim.c" />
+    <ClCompile Include="libs\anim\dynarr.c" />
+    <ClCompile Include="libs\anim\track.c" />
     <ClCompile Include="libs\openctm\compressMG1.c" />
     <ClCompile Include="libs\openctm\compressMG2.c" />
     <ClCompile Include="libs\openctm\compressRAW.c" />
@@ -55,6 +71,18 @@
     <ClCompile Include="libs\openctm\openctm.c" />
     <ClCompile Include="libs\openctm\stream.c" />
     <ClCompile Include="libs\tinyxml2\tinyxml2.cpp" />
+    <ClCompile Include="libs\vmath\basis.cc" />
+    <ClCompile Include="libs\vmath\basis_c.c" />
+    <ClCompile Include="libs\vmath\geom.c" />
+    <ClCompile Include="libs\vmath\matrix.cc" />
+    <ClCompile Include="libs\vmath\matrix_c.c" />
+    <ClCompile Include="libs\vmath\quat.cc" />
+    <ClCompile Include="libs\vmath\quat_c.c" />
+    <ClCompile Include="libs\vmath\ray.cc" />
+    <ClCompile Include="libs\vmath\ray_c.c" />
+    <ClCompile Include="libs\vmath\sphvec.cc" />
+    <ClCompile Include="libs\vmath\vector.cc" />
+    <ClCompile Include="libs\vmath\vmath.c" />
     <ClCompile Include="src\camera.cc" />
     <ClCompile Include="src\chunk.cc" />
     <ClCompile Include="src\goat3d.cc" />
@@ -69,6 +97,13 @@
     <ClCompile Include="src\scene.cc" />
     <ClCompile Include="src\xform_node.cc" />
   </ItemGroup>
+  <ItemGroup>
+    <None Include="libs\vmath\matrix.inl" />
+    <None Include="libs\vmath\quat.inl" />
+    <None Include="libs\vmath\ray.inl" />
+    <None Include="libs\vmath\vector.inl" />
+    <None Include="libs\vmath\vmath.inl" />
+  </ItemGroup>
   <PropertyGroup Label="Globals">
     <ProjectGuid>{86BF319B-9222-4805-918D-DC1B9F77BCEF}</ProjectGuid>
     <Keyword>Win32Proj</Keyword>
@@ -84,7 +119,7 @@
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
     <ConfigurationType>StaticLibrary</ConfigurationType>
     <UseDebugLibraries>true</UseDebugLibraries>
-    <PlatformToolset>v110</PlatformToolset>
+    <PlatformToolset>v100</PlatformToolset>
     <CharacterSet>Unicode</CharacterSet>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
@@ -97,7 +132,7 @@
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
     <ConfigurationType>StaticLibrary</ConfigurationType>
     <UseDebugLibraries>false</UseDebugLibraries>
-    <PlatformToolset>v110</PlatformToolset>
+    <PlatformToolset>v100</PlatformToolset>
     <WholeProgramOptimization>true</WholeProgramOptimization>
     <CharacterSet>Unicode</CharacterSet>
   </PropertyGroup>
@@ -143,14 +178,18 @@
       </PrecompiledHeader>
       <WarningLevel>Level3</WarningLevel>
       <Optimization>Disabled</Optimization>
-      <PreprocessorDefinitions>WIN32;_DEBUG;_LIB;%(PreprocessorDefinitions);OPENCTM_STATIC</PreprocessorDefinitions>
+      <PreprocessorDefinitions>WIN32;_DEBUG;_LIB;%(PreprocessorDefinitions);OPENCTM_STATIC;BUILD_MAXPLUGIN</PreprocessorDefinitions>
       <DisableSpecificWarnings>4244;4305;4996</DisableSpecificWarnings>
-      <AdditionalIncludeDirectories>$(SolutionDir)\libs\openctm;$(SolutionDir)\libs\openctm\liblzma;$(SolutionDir)\libs\tinyxml2</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>$(SolutionDir)\src;$(SolutionDir)\libs\openctm;$(SolutionDir)\libs\openctm\liblzma;$(SolutionDir)\libs\tinyxml2</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
       <GenerateDebugInformation>true</GenerateDebugInformation>
     </Link>
+    <Lib>
+      <AdditionalDependencies>
+      </AdditionalDependencies>
+    </Lib>
   </ItemDefinitionGroup>
   <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
     <ClCompile>
@@ -180,7 +219,7 @@
       <IntrinsicFunctions>true</IntrinsicFunctions>
       <PreprocessorDefinitions>WIN32;NDEBUG;_LIB;%(PreprocessorDefinitions);OPENCTM_STATIC</PreprocessorDefinitions>
       <DisableSpecificWarnings>4244;4305;4996</DisableSpecificWarnings>
-      <AdditionalIncludeDirectories>$(SolutionDir)\libs\openctm;$(SolutionDir)\libs\openctm\liblzma;$(SolutionDir)\libs\tinyxml2</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>$(SolutionDir)\src;$(SolutionDir)\libs\openctm;$(SolutionDir)\libs\openctm\liblzma;$(SolutionDir)\libs\tinyxml2</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
@@ -188,6 +227,10 @@
       <EnableCOMDATFolding>true</EnableCOMDATFolding>
       <OptimizeReferences>true</OptimizeReferences>
     </Link>
+    <Lib>
+      <AdditionalDependencies>
+      </AdditionalDependencies>
+    </Lib>
   </ItemDefinitionGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
diff -r 1c601bf07b86 -r 4deb0b12fe14 goat3d.vcxproj.filters
--- a/goat3d.vcxproj.filters	Sat Sep 28 20:36:55 2013 +0300
+++ b/goat3d.vcxproj.filters	Sun Sep 29 08:20:19 2013 +0300
@@ -17,6 +17,12 @@
     <Filter Include="libs\openctm\liblzma">
       <UniqueIdentifier>{18a85bea-ddb2-49a1-a7b1-bd11b23d5549}</UniqueIdentifier>
     </Filter>
+    <Filter Include="libs\vmath">
+      <UniqueIdentifier>{dc4ab435-af1e-4786-8011-a7a793e229eb}</UniqueIdentifier>
+    </Filter>
+    <Filter Include="libs\anim">
+      <UniqueIdentifier>{05d8010d-b4d9-4829-80ea-8c8d018e1e04}</UniqueIdentifier>
+    </Filter>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="src\camera.h">
@@ -88,6 +94,45 @@
     <ClInclude Include="libs\openctm\liblzma\Types.h">
       <Filter>libs\openctm\liblzma</Filter>
     </ClInclude>
+    <ClInclude Include="libs\anim\anim.h">
+      <Filter>libs\anim</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\anim\dynarr.h">
+      <Filter>libs\anim</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\anim\track.h">
+      <Filter>libs\anim</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\basis.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\geom.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\matrix.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\quat.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\ray.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\sphvec.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\vector.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\vmath.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\vmath_config.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
+    <ClInclude Include="libs\vmath\vmath_types.h">
+      <Filter>libs\vmath</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="src\goat3d.cc">
@@ -162,5 +207,67 @@
     <ClCompile Include="libs\openctm\liblzma\LzmaLib.c">
       <Filter>libs\openctm\liblzma</Filter>
     </ClCompile>
+    <ClCompile Include="libs\anim\anim.c">
+      <Filter>libs\anim</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\anim\dynarr.c">
+      <Filter>libs\anim</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\anim\track.c">
+      <Filter>libs\anim</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\basis.cc">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\basis_c.c">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\geom.c">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\matrix.cc">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\matrix_c.c">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\quat.cc">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\quat_c.c">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\ray.cc">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\ray_c.c">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\sphvec.cc">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\vector.cc">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+    <ClCompile Include="libs\vmath\vmath.c">
+      <Filter>libs\vmath</Filter>
+    </ClCompile>
+  </ItemGroup>
+  <ItemGroup>
+    <None Include="libs\vmath\matrix.inl">
+      <Filter>libs\vmath</Filter>
+    </None>
+    <None Include="libs\vmath\quat.inl">
+      <Filter>libs\vmath</Filter>
+    </None>
+    <None Include="libs\vmath\ray.inl">
+      <Filter>libs\vmath</Filter>
+    </None>
+    <None Include="libs\vmath\vector.inl">
+      <Filter>libs\vmath</Filter>
+    </None>
+    <None Include="libs\vmath\vmath.inl">
+      <Filter>libs\vmath</Filter>
+    </None>
   </ItemGroup>
 </Project>
\ No newline at end of file
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/anim/anim.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/anim/anim.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,455 @@
+#include <stdlib.h>
+#include <limits.h>
+#include <assert.h>
+#include "anim.h"
+#include "dynarr.h"
+
+#define ROT_USE_SLERP
+
+static void invalidate_cache(struct anm_node *node);
+
+int anm_init_node(struct anm_node *node)
+{
+	int i, j;
+	static const float defaults[] = {
+		0.0f, 0.0f, 0.0f,		/* default position */
+		0.0f, 0.0f, 0.0f, 1.0f,	/* default rotation quat */
+		1.0f, 1.0f, 1.0f		/* default scale factor */
+	};
+
+	memset(node, 0, sizeof *node);
+
+	for(i=0; i<ANM_NUM_TRACKS; i++) {
+		if(anm_init_track(node->tracks + i) == -1) {
+			for(j=0; j<i; j++) {
+				anm_destroy_track(node->tracks + i);
+			}
+		}
+		anm_set_track_default(node->tracks + i, defaults[i]);
+	}
+
+	node->cache.time = ANM_TIME_INVAL;
+	node->cache.inv_time = ANM_TIME_INVAL;
+	return 0;
+}
+
+void anm_destroy_node(struct anm_node *node)
+{
+	int i;
+	free(node->name);
+
+	for(i=0; i<ANM_NUM_TRACKS; i++) {
+		anm_destroy_track(node->tracks + i);
+	}
+}
+
+void anm_destroy_node_tree(struct anm_node *tree)
+{
+	struct anm_node *c, *tmp;
+
+	if(!tree) return;
+
+	c = tree->child;
+	while(c) {
+		tmp = c;
+		c = c->next;
+
+		anm_destroy_node_tree(tmp);
+	}
+	anm_destroy_node(tree);
+}
+
+struct anm_node *anm_create_node(void)
+{
+	struct anm_node *n;
+
+	if((n = malloc(sizeof *n))) {
+		if(anm_init_node(n) == -1) {
+			free(n);
+			return 0;
+		}
+	}
+	return n;
+}
+
+void anm_free_node(struct anm_node *node)
+{
+	anm_destroy_node(node);
+	free(node);
+}
+
+void anm_free_node_tree(struct anm_node *tree)
+{
+	struct anm_node *c, *tmp;
+
+	if(!tree) return;
+
+	c = tree->child;
+	while(c) {
+		tmp = c;
+		c = c->next;
+
+		anm_free_node_tree(tmp);
+	}
+
+	anm_free_node(tree);
+}
+
+int anm_set_node_name(struct anm_node *node, const char *name)
+{
+	char *str;
+
+	if(!(str = malloc(strlen(name) + 1))) {
+		return -1;
+	}
+	strcpy(str, name);
+	free(node->name);
+	node->name = str;
+	return 0;
+}
+
+const char *anm_get_node_name(struct anm_node *node)
+{
+	return node->name ? node->name : "";
+}
+
+void anm_set_interpolator(struct anm_node *node, enum anm_interpolator in)
+{
+	int i;
+
+	for(i=0; i<ANM_NUM_TRACKS; i++) {
+		anm_set_track_interpolator(node->tracks + i, in);
+	}
+	invalidate_cache(node);
+}
+
+void anm_set_extrapolator(struct anm_node *node, enum anm_extrapolator ex)
+{
+	int i;
+
+	for(i=0; i<ANM_NUM_TRACKS; i++) {
+		anm_set_track_extrapolator(node->tracks + i, ex);
+	}
+	invalidate_cache(node);
+}
+
+void anm_link_node(struct anm_node *p, struct anm_node *c)
+{
+	c->next = p->child;
+	p->child = c;
+
+	c->parent = p;
+	invalidate_cache(c);
+}
+
+int anm_unlink_node(struct anm_node *p, struct anm_node *c)
+{
+	struct anm_node *iter;
+
+	if(p->child == c) {
+		p->child = c->next;
+		c->next = 0;
+		invalidate_cache(c);
+		return 0;
+	}
+
+	iter = p->child;
+	while(iter->next) {
+		if(iter->next == c) {
+			iter->next = c->next;
+			c->next = 0;
+			invalidate_cache(c);
+			return 0;
+		}
+	}
+	return -1;
+}
+
+void anm_set_position(struct anm_node *node, vec3_t pos, anm_time_t tm)
+{
+	anm_set_value(node->tracks + ANM_TRACK_POS_X, tm, pos.x);
+	anm_set_value(node->tracks + ANM_TRACK_POS_Y, tm, pos.y);
+	anm_set_value(node->tracks + ANM_TRACK_POS_Z, tm, pos.z);
+	invalidate_cache(node);
+}
+
+vec3_t anm_get_node_position(struct anm_node *node, anm_time_t tm)
+{
+	vec3_t v;
+	v.x = anm_get_value(node->tracks + ANM_TRACK_POS_X, tm);
+	v.y = anm_get_value(node->tracks + ANM_TRACK_POS_Y, tm);
+	v.z = anm_get_value(node->tracks + ANM_TRACK_POS_Z, tm);
+	return v;
+}
+
+void anm_set_rotation(struct anm_node *node, quat_t rot, anm_time_t tm)
+{
+	anm_set_value(node->tracks + ANM_TRACK_ROT_X, tm, rot.x);
+	anm_set_value(node->tracks + ANM_TRACK_ROT_Y, tm, rot.y);
+	anm_set_value(node->tracks + ANM_TRACK_ROT_Z, tm, rot.z);
+	anm_set_value(node->tracks + ANM_TRACK_ROT_W, tm, rot.w);
+	invalidate_cache(node);
+}
+
+quat_t anm_get_node_rotation(struct anm_node *node, anm_time_t tm)
+{
+#ifndef ROT_USE_SLERP
+	quat_t q;
+	q.x = anm_get_value(node->tracks + ANM_TRACK_ROT_X, tm);
+	q.y = anm_get_value(node->tracks + ANM_TRACK_ROT_Y, tm);
+	q.z = anm_get_value(node->tracks + ANM_TRACK_ROT_Z, tm);
+	q.w = anm_get_value(node->tracks + ANM_TRACK_ROT_W, tm);
+	return q;
+#else
+	int idx0, idx1, last_idx;
+	anm_time_t tstart, tend;
+	float t, dt;
+	struct anm_track *track_x, *track_y, *track_z, *track_w;
+	quat_t q, q1, q2;
+
+	track_x = node->tracks + ANM_TRACK_ROT_X;
+	track_y = node->tracks + ANM_TRACK_ROT_Y;
+	track_z = node->tracks + ANM_TRACK_ROT_Z;
+	track_w = node->tracks + ANM_TRACK_ROT_W;
+
+	if(!track_x->count) {
+		q.x = track_x->def_val;
+		q.y = track_y->def_val;
+		q.z = track_z->def_val;
+		q.w = track_w->def_val;
+		return q;
+	}
+
+	last_idx = track_x->count - 1;
+
+	tstart = track_x->keys[0].time;
+	tend = track_x->keys[last_idx].time;
+
+	if(tstart == tend) {
+		q.x = track_x->keys[0].val;
+		q.y = track_y->keys[0].val;
+		q.z = track_z->keys[0].val;
+		q.w = track_w->keys[0].val;
+		return q;
+	}
+
+	tm = anm_remap_time(track_x, tm, tstart, tend);
+
+	idx0 = anm_get_key_interval(track_x, tm);
+	assert(idx0 >= 0 && idx0 < track_x->count);
+	idx1 = idx0 + 1;
+
+	if(idx0 == last_idx) {
+		q.x = track_x->keys[idx0].val;
+		q.y = track_y->keys[idx0].val;
+		q.z = track_z->keys[idx0].val;
+		q.w = track_w->keys[idx0].val;
+		return q;
+	}
+
+	dt = (float)(track_x->keys[idx1].time - track_x->keys[idx0].time);
+	t = (float)(tm - track_x->keys[idx0].time) / dt;
+
+	q1.x = track_x->keys[idx0].val;
+	q1.y = track_y->keys[idx0].val;
+	q1.z = track_z->keys[idx0].val;
+	q1.w = track_w->keys[idx0].val;
+
+	q2.x = track_x->keys[idx1].val;
+	q2.y = track_y->keys[idx1].val;
+	q2.z = track_z->keys[idx1].val;
+	q2.w = track_w->keys[idx1].val;
+
+	/*q1 = quat_normalize(q1);
+	q2 = quat_normalize(q2);*/
+
+	return quat_slerp(q1, q2, t);
+#endif
+}
+
+void anm_set_scaling(struct anm_node *node, vec3_t scl, anm_time_t tm)
+{
+	anm_set_value(node->tracks + ANM_TRACK_SCL_X, tm, scl.x);
+	anm_set_value(node->tracks + ANM_TRACK_SCL_Y, tm, scl.y);
+	anm_set_value(node->tracks + ANM_TRACK_SCL_Z, tm, scl.z);
+	invalidate_cache(node);
+}
+
+vec3_t anm_get_node_scaling(struct anm_node *node, anm_time_t tm)
+{
+	vec3_t v;
+	v.x = anm_get_value(node->tracks + ANM_TRACK_SCL_X, tm);
+	v.y = anm_get_value(node->tracks + ANM_TRACK_SCL_Y, tm);
+	v.z = anm_get_value(node->tracks + ANM_TRACK_SCL_Z, tm);
+	return v;
+}
+
+
+vec3_t anm_get_position(struct anm_node *node, anm_time_t tm)
+{
+	mat4_t xform;
+	vec3_t pos = {0.0, 0.0, 0.0};
+
+	if(!node->parent) {
+		return anm_get_node_position(node, tm);
+	}
+
+	anm_get_matrix(node, xform, tm);
+	return v3_transform(pos, xform);
+}
+
+quat_t anm_get_rotation(struct anm_node *node, anm_time_t tm)
+{
+	quat_t rot, prot;
+	rot = anm_get_node_rotation(node, tm);
+
+	if(!node->parent) {
+		return rot;
+	}
+
+	prot = anm_get_rotation(node->parent, tm);
+	return quat_mul(prot, rot);
+}
+
+vec3_t anm_get_scaling(struct anm_node *node, anm_time_t tm)
+{
+	vec3_t s, ps;
+	s = anm_get_node_scaling(node, tm);
+
+	if(!node->parent) {
+		return s;
+	}
+
+	ps = anm_get_scaling(node->parent, tm);
+	return v3_mul(s, ps);
+}
+
+void anm_set_pivot(struct anm_node *node, vec3_t piv)
+{
+	node->pivot = piv;
+}
+
+vec3_t anm_get_pivot(struct anm_node *node)
+{
+	return node->pivot;
+}
+
+void anm_get_node_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm)
+{
+	int i;
+	mat4_t rmat;
+	vec3_t pos, scale;
+	quat_t rot;
+
+	pos = anm_get_node_position(node, tm);
+	rot = anm_get_node_rotation(node, tm);
+	scale = anm_get_node_scaling(node, tm);
+
+	m4_set_translation(mat, node->pivot.x, node->pivot.y, node->pivot.z);
+
+	quat_to_mat4(rmat, rot);
+	for(i=0; i<3; i++) {
+		mat[i][0] = rmat[i][0];
+		mat[i][1] = rmat[i][1];
+		mat[i][2] = rmat[i][2];
+	}
+	/* this loop is equivalent to: m4_mult(mat, mat, rmat); */
+
+	mat[0][0] *= scale.x; mat[0][1] *= scale.y; mat[0][2] *= scale.z; mat[0][3] += pos.x;
+	mat[1][0] *= scale.x; mat[1][1] *= scale.y; mat[1][2] *= scale.z; mat[1][3] += pos.y;
+	mat[2][0] *= scale.x; mat[2][1] *= scale.y; mat[2][2] *= scale.z; mat[2][3] += pos.z;
+
+	m4_translate(mat, -node->pivot.x, -node->pivot.y, -node->pivot.z);
+
+	/* that's basically: pivot * rotation * translation * scaling * -pivot */
+}
+
+void anm_get_node_inv_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm)
+{
+	mat4_t tmp;
+	anm_get_node_matrix(node, tmp, tm);
+	m4_inverse(mat, tmp);
+}
+
+void anm_get_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm)
+{
+	if(node->cache.time != tm) {
+		anm_get_node_matrix(node, node->cache.matrix, tm);
+
+		if(node->parent) {
+			mat4_t parent_mat;
+
+			anm_get_matrix(node->parent, parent_mat, tm);
+			m4_mult(node->cache.matrix, parent_mat, node->cache.matrix);
+		}
+		node->cache.time = tm;
+	}
+	m4_copy(mat, node->cache.matrix);
+}
+
+void anm_get_inv_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm)
+{
+	if(node->cache.inv_time != tm) {
+		anm_get_matrix(node, mat, tm);
+		m4_inverse(node->cache.inv_matrix, mat);
+		node->cache.inv_time = tm;
+	}
+	m4_copy(mat, node->cache.inv_matrix);
+}
+
+anm_time_t anm_get_start_time(struct anm_node *node)
+{
+	int i;
+	struct anm_node *c;
+	anm_time_t res = LONG_MAX;
+
+	for(i=0; i<ANM_NUM_TRACKS; i++) {
+		if(node->tracks[i].count) {
+			anm_time_t tm = node->tracks[i].keys[0].time;
+			if(tm < res) {
+				res = tm;
+			}
+		}
+	}
+
+	c = node->child;
+	while(c) {
+		anm_time_t tm = anm_get_start_time(c);
+		if(tm < res) {
+			res = tm;
+		}
+		c = c->next;
+	}
+	return res;
+}
+
+anm_time_t anm_get_end_time(struct anm_node *node)
+{
+	int i;
+	struct anm_node *c;
+	anm_time_t res = LONG_MIN;
+
+	for(i=0; i<ANM_NUM_TRACKS; i++) {
+		if(node->tracks[i].count) {
+			anm_time_t tm = node->tracks[i].keys[node->tracks[i].count - 1].time;
+			if(tm > res) {
+				res = tm;
+			}
+		}
+	}
+
+	c = node->child;
+	while(c) {
+		anm_time_t tm = anm_get_end_time(c);
+		if(tm > res) {
+			res = tm;
+		}
+		c = c->next;
+	}
+	return res;
+}
+
+static void invalidate_cache(struct anm_node *node)
+{
+	node->cache.time = node->cache.inv_time = ANM_TIME_INVAL;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/anim/anim.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/anim/anim.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,112 @@
+#ifndef LIBANIM_H_
+#define LIBANIM_H_
+
+#include "config.h"
+
+#include <vmath/vector.h>
+#include <vmath/quat.h>
+#include <vmath/matrix.h>
+#include "track.h"
+
+enum {
+	ANM_TRACK_POS_X,
+	ANM_TRACK_POS_Y,
+	ANM_TRACK_POS_Z,
+
+	ANM_TRACK_ROT_X,
+	ANM_TRACK_ROT_Y,
+	ANM_TRACK_ROT_Z,
+	ANM_TRACK_ROT_W,
+
+	ANM_TRACK_SCL_X,
+	ANM_TRACK_SCL_Y,
+	ANM_TRACK_SCL_Z,
+
+	ANM_NUM_TRACKS
+};
+
+struct anm_node {
+	char *name;
+
+	struct anm_track tracks[ANM_NUM_TRACKS];
+	vec3_t pivot;
+
+	/* matrix cache */
+	struct mat_cache {
+		mat4_t matrix, inv_matrix;
+		anm_time_t time, inv_time;
+	} cache;
+
+	struct anm_node *parent;
+	struct anm_node *child;
+	struct anm_node *next;
+};
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* node constructor and destructor */
+int anm_init_node(struct anm_node *node);
+void anm_destroy_node(struct anm_node *node);
+
+/* recursively destroy an animation node tree */
+void anm_destroy_node_tree(struct anm_node *tree);
+
+/* helper functions to allocate/construct and destroy/free with
+ * a single call. They call anm_init_node and anm_destroy_node
+ * internally.
+ */
+struct anm_node *anm_create_node(void);
+void anm_free_node(struct anm_node *node);
+
+/* recursively destroy and free the nodes of a node tree */
+void anm_free_node_tree(struct anm_node *tree);
+
+int anm_set_node_name(struct anm_node *node, const char *name);
+const char *anm_get_node_name(struct anm_node *node);
+
+void anm_set_interpolator(struct anm_node *node, enum anm_interpolator in);
+void anm_set_extrapolator(struct anm_node *node, enum anm_extrapolator ex);
+
+/* link and unlink nodes with parent/child relations */
+void anm_link_node(struct anm_node *parent, struct anm_node *child);
+int anm_unlink_node(struct anm_node *parent, struct anm_node *child);
+
+void anm_set_position(struct anm_node *node, vec3_t pos, anm_time_t tm);
+vec3_t anm_get_node_position(struct anm_node *node, anm_time_t tm);
+
+void anm_set_rotation(struct anm_node *node, quat_t rot, anm_time_t tm);
+quat_t anm_get_node_rotation(struct anm_node *node, anm_time_t tm);
+
+void anm_set_scaling(struct anm_node *node, vec3_t scl, anm_time_t tm);
+vec3_t anm_get_node_scaling(struct anm_node *node, anm_time_t tm);
+
+/* these three return the full p/r/s taking hierarchy into account */
+vec3_t anm_get_position(struct anm_node *node, anm_time_t tm);
+quat_t anm_get_rotation(struct anm_node *node, anm_time_t tm);
+vec3_t anm_get_scaling(struct anm_node *node, anm_time_t tm);
+
+void anm_set_pivot(struct anm_node *node, vec3_t pivot);
+vec3_t anm_get_pivot(struct anm_node *node);
+
+/* these calculate the matrix and inverse matrix of this node alone */
+void anm_get_node_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm);
+void anm_get_node_inv_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm);
+
+/* These calculate the matrix and inverse matrix of this node taking hierarchy
+ * into account. The results are cached in thread-specific storage and returned
+ * if there's no change in time or tracks from the last query...
+ */
+void anm_get_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm);
+void anm_get_inv_matrix(struct anm_node *node, mat4_t mat, anm_time_t tm);
+
+/* those return the start and end times of the whole tree */
+anm_time_t anm_get_start_time(struct anm_node *node);
+anm_time_t anm_get_end_time(struct anm_node *node);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif	/* LIBANIM_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/anim/config.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/anim/config.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,6 @@
+#ifndef ANIM_CONFIG_H_
+#define ANIM_CONFIG_H_
+
+#undef ANIM_THREAD_SAFE
+
+#endif	/* ANIM_CONFIG_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/anim/dynarr.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/anim/dynarr.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,122 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "dynarr.h"
+
+/* The array descriptor keeps auxilliary information needed to manipulate
+ * the dynamic array. It's allocated adjacent to the array buffer.
+ */
+struct arrdesc {
+	int nelem, szelem;
+	int max_elem;
+	int bufsz;	/* not including the descriptor */
+};
+
+#define DESC(x)		((struct arrdesc*)((char*)(x) - sizeof(struct arrdesc)))
+
+void *dynarr_alloc(int elem, int szelem)
+{
+	struct arrdesc *desc;
+
+	if(!(desc = malloc(elem * szelem + sizeof *desc))) {
+		return 0;
+	}
+	desc->nelem = desc->max_elem = elem;
+	desc->szelem = szelem;
+	desc->bufsz = elem * szelem;
+	return (char*)desc + sizeof *desc;
+}
+
+void dynarr_free(void *da)
+{
+	if(da) {
+		free(DESC(da));
+	}
+}
+
+void *dynarr_resize(void *da, int elem)
+{
+	int newsz;
+	void *tmp;
+	struct arrdesc *desc;
+
+	if(!da) return 0;
+	desc = DESC(da);
+
+	newsz = desc->szelem * elem;
+
+	if(!(tmp = realloc(desc, newsz + sizeof *desc))) {
+		return 0;
+	}
+	desc = tmp;
+
+	desc->nelem = desc->max_elem = elem;
+	desc->bufsz = newsz;
+	return (char*)desc + sizeof *desc;
+}
+
+int dynarr_empty(void *da)
+{
+	return DESC(da)->nelem ? 0 : 1;
+}
+
+int dynarr_size(void *da)
+{
+	return DESC(da)->nelem;
+}
+
+
+/* stack semantics */
+void *dynarr_push(void *da, void *item)
+{
+	struct arrdesc *desc;
+	int nelem;
+
+	desc = DESC(da);
+	nelem = desc->nelem;
+
+	if(nelem >= desc->max_elem) {
+		/* need to resize */
+		struct arrdesc *tmp;
+		int newsz = desc->max_elem ? desc->max_elem * 2 : 1;
+
+		if(!(tmp = dynarr_resize(da, newsz))) {
+			fprintf(stderr, "failed to resize\n");
+			return da;
+		}
+		da = tmp;
+		desc = DESC(da);
+		desc->nelem = nelem;
+	}
+
+	memcpy((char*)da + desc->nelem++ * desc->szelem, item, desc->szelem);
+	return da;
+}
+
+void *dynarr_pop(void *da)
+{
+	struct arrdesc *desc;
+	int nelem;
+
+	desc = DESC(da);
+	nelem = desc->nelem;
+
+	if(!nelem) return da;
+
+	if(nelem <= desc->max_elem / 3) {
+		/* reclaim space */
+		struct arrdesc *tmp;
+		int newsz = desc->max_elem / 2;
+
+		if(!(tmp = dynarr_resize(da, newsz))) {
+			fprintf(stderr, "failed to resize\n");
+			return da;
+		}
+		da = tmp;
+		desc = DESC(da);
+		desc->nelem = nelem;
+	}
+	desc->nelem--;
+
+	return da;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/anim/dynarr.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/anim/dynarr.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,18 @@
+#ifndef DYNARR_H_
+#define DYNARR_H_
+
+#include "config.h"
+
+void *dynarr_alloc(int elem, int szelem);
+void dynarr_free(void *da);
+void *dynarr_resize(void *da, int elem);
+
+int dynarr_empty(void *da);
+int dynarr_size(void *da);
+
+/* stack semantics */
+void *dynarr_push(void *da, void *item);
+void *dynarr_pop(void *da);
+
+
+#endif	/* DYNARR_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/anim/track.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/anim/track.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,316 @@
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include "track.h"
+#include "dynarr.h"
+
+static int keycmp(const void *a, const void *b);
+static int find_prev_key(struct anm_keyframe *arr, int start, int end, anm_time_t tm);
+
+static float interp_step(float v0, float v1, float v2, float v3, float t);
+static float interp_linear(float v0, float v1, float v2, float v3, float t);
+static float interp_cubic(float v0, float v1, float v2, float v3, float t);
+
+static anm_time_t remap_extend(anm_time_t tm, anm_time_t start, anm_time_t end);
+static anm_time_t remap_clamp(anm_time_t tm, anm_time_t start, anm_time_t end);
+static anm_time_t remap_repeat(anm_time_t tm, anm_time_t start, anm_time_t end);
+static anm_time_t remap_pingpong(anm_time_t tm, anm_time_t start, anm_time_t end);
+
+/* XXX keep this in sync with enum anm_interpolator at track.h */
+static float (*interp[])(float, float, float, float, float) = {
+	interp_step,
+	interp_linear,
+	interp_cubic,
+	0
+};
+
+/* XXX keep this in sync with enum anm_extrapolator at track.h */
+static anm_time_t (*remap_time[])(anm_time_t, anm_time_t, anm_time_t) = {
+	remap_extend,
+	remap_clamp,
+	remap_repeat,
+	remap_pingpong,
+	0
+};
+
+int anm_init_track(struct anm_track *track)
+{
+	memset(track, 0, sizeof *track);
+
+	if(!(track->keys = dynarr_alloc(0, sizeof *track->keys))) {
+		return -1;
+	}
+	track->interp = ANM_INTERP_LINEAR;
+	track->extrap = ANM_EXTRAP_CLAMP;
+	return 0;
+}
+
+void anm_destroy_track(struct anm_track *track)
+{
+	dynarr_free(track->keys);
+}
+
+struct anm_track *anm_create_track(void)
+{
+	struct anm_track *track;
+
+	if((track = malloc(sizeof *track))) {
+		if(anm_init_track(track) == -1) {
+			free(track);
+			return 0;
+		}
+	}
+	return track;
+}
+
+void anm_free_track(struct anm_track *track)
+{
+	anm_destroy_track(track);
+	free(track);
+}
+
+void anm_copy_track(struct anm_track *dest, struct anm_track *src)
+{
+	free(dest->name);
+	if(dest->keys) {
+		dynarr_free(dest->keys);
+	}
+
+	if(src->name) {
+		dest->name = malloc(strlen(src->name) + 1);
+		strcpy(dest->name, src->name);
+	}
+
+	dest->count = src->count;
+	dest->keys = dynarr_alloc(src->count, sizeof *dest->keys);
+	memcpy(dest->keys, src->keys, src->count * sizeof *dest->keys);
+
+	dest->def_val = src->def_val;
+	dest->interp = src->interp;
+	dest->extrap = src->extrap;
+}
+
+int anm_set_track_name(struct anm_track *track, const char *name)
+{
+	char *tmp;
+
+	if(!(tmp = malloc(strlen(name) + 1))) {
+		return -1;
+	}
+	free(track->name);
+	track->name = tmp;
+	return 0;
+}
+
+const char *anm_get_track_name(struct anm_track *track)
+{
+	return track->name;
+}
+
+void anm_set_track_interpolator(struct anm_track *track, enum anm_interpolator in)
+{
+	track->interp = in;
+}
+
+void anm_set_track_extrapolator(struct anm_track *track, enum anm_extrapolator ex)
+{
+	track->extrap = ex;
+}
+
+anm_time_t anm_remap_time(struct anm_track *track, anm_time_t tm, anm_time_t start, anm_time_t end)
+{
+	return remap_time[track->extrap](tm, start, end);
+}
+
+void anm_set_track_default(struct anm_track *track, float def)
+{
+	track->def_val = def;
+}
+
+int anm_set_keyframe(struct anm_track *track, struct anm_keyframe *key)
+{
+	int idx = anm_get_key_interval(track, key->time);
+
+	/* if we got a valid keyframe index, compare them... */
+	if(idx >= 0 && idx < track->count && keycmp(key, track->keys + idx) == 0) {
+		/* ... it's the same key, just update the value */
+		track->keys[idx].val = key->val;
+	} else {
+		/* ... it's a new key, add it and re-sort them */
+		void *tmp;
+		if(!(tmp = dynarr_push(track->keys, key))) {
+			return -1;
+		}
+		track->keys = tmp;
+		/* TODO lazy qsort */
+		qsort(track->keys, ++track->count, sizeof *track->keys, keycmp);
+	}
+	return 0;
+}
+
+static int keycmp(const void *a, const void *b)
+{
+	return ((struct anm_keyframe*)a)->time - ((struct anm_keyframe*)b)->time;
+}
+
+struct anm_keyframe *anm_get_keyframe(struct anm_track *track, int idx)
+{
+	if(idx < 0 || idx >= track->count) {
+		return 0;
+	}
+	return track->keys + idx;
+}
+
+int anm_get_key_interval(struct anm_track *track, anm_time_t tm)
+{
+	int last;
+
+	if(!track->count || tm < track->keys[0].time) {
+		return -1;
+	}
+
+	last = track->count - 1;
+	if(tm > track->keys[last].time) {
+		return last;
+	}
+
+	return find_prev_key(track->keys, 0, last, tm);
+}
+
+static int find_prev_key(struct anm_keyframe *arr, int start, int end, anm_time_t tm)
+{
+	int mid;
+
+	if(end - start <= 1) {
+		return start;
+	}
+
+	mid = (start + end) / 2;
+	if(tm < arr[mid].time) {
+		return find_prev_key(arr, start, mid, tm);
+	}
+	if(tm > arr[mid].time) {
+		return find_prev_key(arr, mid, end, tm);
+	}
+	return mid;
+}
+
+int anm_set_value(struct anm_track *track, anm_time_t tm, float val)
+{
+	struct anm_keyframe key;
+	key.time = tm;
+	key.val = val;
+
+	return anm_set_keyframe(track, &key);
+}
+
+float anm_get_value(struct anm_track *track, anm_time_t tm)
+{
+	int idx0, idx1, last_idx;
+	anm_time_t tstart, tend;
+	float t, dt;
+	float v0, v1, v2, v3;
+
+	if(!track->count) {
+		return track->def_val;
+	}
+
+	last_idx = track->count - 1;
+
+	tstart = track->keys[0].time;
+	tend = track->keys[last_idx].time;
+
+	if(tstart == tend) {
+		return track->keys[0].val;
+	}
+
+	tm = remap_time[track->extrap](tm, tstart, tend);
+
+	idx0 = anm_get_key_interval(track, tm);
+	assert(idx0 >= 0 && idx0 < track->count);
+	idx1 = idx0 + 1;
+
+	if(idx0 == last_idx) {
+		return track->keys[idx0].val;
+	}
+
+	dt = (float)(track->keys[idx1].time - track->keys[idx0].time);
+	t = (float)(tm - track->keys[idx0].time) / dt;
+
+	v1 = track->keys[idx0].val;
+	v2 = track->keys[idx1].val;
+
+	/* get the neigboring values to allow for cubic interpolation */
+	v0 = idx0 > 0 ? track->keys[idx0 - 1].val : v1;
+	v3 = idx1 < last_idx ? track->keys[idx1 + 1].val : v2;
+
+	return interp[track->interp](v0, v1, v2, v3, t);
+}
+
+
+static float interp_step(float v0, float v1, float v2, float v3, float t)
+{
+	return v1;
+}
+
+static float interp_linear(float v0, float v1, float v2, float v3, float t)
+{
+	return v1 + (v2 - v1) * t;
+}
+
+static float interp_cubic(float a, float b, float c, float d, float t)
+{
+	float x, y, z, w;
+	float tsq = t * t;
+
+	x = -a + 3.0 * b - 3.0 * c + d;
+	y = 2.0 * a - 5.0 * b + 4.0 * c - d;
+	z = c - a;
+	w = 2.0 * b;
+
+	return 0.5 * (x * tsq * t + y * tsq + z * t + w);
+}
+
+static anm_time_t remap_extend(anm_time_t tm, anm_time_t start, anm_time_t end)
+{
+	return remap_repeat(tm, start, end);
+}
+
+static anm_time_t remap_clamp(anm_time_t tm, anm_time_t start, anm_time_t end)
+{
+	if(start == end) {
+		return start;
+	}
+	return tm < start ? start : (tm >= end ? end - 1 : tm);
+}
+
+static anm_time_t remap_repeat(anm_time_t tm, anm_time_t start, anm_time_t end)
+{
+	anm_time_t x, interv = end - start;
+
+	if(interv == 0) {
+		return start;
+	}
+
+	x = (tm - start) % interv;
+	if(x < 0) {
+		x += interv;
+	}
+	return x + start;
+
+	/*if(tm < start) {
+		while(tm < start) {
+			tm += interv;
+		}
+		return tm;
+	}
+	return (tm - start) % interv + start;*/
+}
+
+static anm_time_t remap_pingpong(anm_time_t tm, anm_time_t start, anm_time_t end)
+{
+	anm_time_t interv = end - start;
+	anm_time_t x = remap_repeat(tm, start, end + interv);
+
+	return x > end ? end + interv - x : x;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/anim/track.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/anim/track.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,101 @@
+/* An animation track defines the values of a single scalar over time
+ * and supports various interpolation and extrapolation modes.
+ */
+#ifndef LIBANIM_TRACK_H_
+#define LIBANIM_TRACK_H_
+
+#include <limits.h>
+#include "config.h"
+
+enum anm_interpolator {
+	ANM_INTERP_STEP,
+	ANM_INTERP_LINEAR,
+	ANM_INTERP_CUBIC
+};
+
+enum anm_extrapolator {
+	ANM_EXTRAP_EXTEND,	/* extend to infinity */
+	ANM_EXTRAP_CLAMP,	/* clamp to last value */
+	ANM_EXTRAP_REPEAT,	/* repeat motion */
+	ANM_EXTRAP_PINGPONG	/* repeat with mirroring */
+};
+
+typedef long anm_time_t;
+#define ANM_TIME_INVAL	LONG_MIN
+
+#define ANM_SEC2TM(x)	((anm_time_t)((x) * 1000))
+#define ANM_MSEC2TM(x)	((anm_time_t)(x))
+#define ANM_TM2SEC(x)	((x) / 1000.0)
+#define ANM_TM2MSEC(x)	(x)
+
+struct anm_keyframe {
+	anm_time_t time;
+	float val;
+};
+
+struct anm_track {
+	char *name;
+	int count;
+	struct anm_keyframe *keys;
+
+	float def_val;
+
+	enum anm_interpolator interp;
+	enum anm_extrapolator extrap;
+};
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* track constructor and destructor */
+int anm_init_track(struct anm_track *track);
+void anm_destroy_track(struct anm_track *track);
+
+/* helper functions that use anm_init_track and anm_destroy_track internally */
+struct anm_track *anm_create_track(void);
+void anm_free_track(struct anm_track *track);
+
+/* copies track src to dest
+ * XXX: dest must have been initialized first
+ */
+void anm_copy_track(struct anm_track *dest, struct anm_track *src);
+
+int anm_set_track_name(struct anm_track *track, const char *name);
+const char *anm_get_track_name(struct anm_track *track);
+
+void anm_set_track_interpolator(struct anm_track *track, enum anm_interpolator in);
+void anm_set_track_extrapolator(struct anm_track *track, enum anm_extrapolator ex);
+
+anm_time_t anm_remap_time(struct anm_track *track, anm_time_t tm, anm_time_t start, anm_time_t end);
+
+void anm_set_track_default(struct anm_track *track, float def);
+
+/* set or update a keyframe */
+int anm_set_keyframe(struct anm_track *track, struct anm_keyframe *key);
+
+/* get the idx-th keyframe, returns null if it doesn't exist */
+struct anm_keyframe *anm_get_keyframe(struct anm_track *track, int idx);
+
+/* Finds the 0-based index of the intra-keyframe interval which corresponds
+ * to the specified time. If the time falls exactly onto the N-th keyframe
+ * the function returns N.
+ *
+ * Special cases:
+ * - if the time is before the first keyframe -1 is returned.
+ * - if the time is after the last keyframe, the index of the last keyframe
+ *   is returned.
+ */
+int anm_get_key_interval(struct anm_track *track, anm_time_t tm);
+
+int anm_set_value(struct anm_track *track, anm_time_t tm, float val);
+
+/* evaluates and returns the value of the track for a particular time */
+float anm_get_value(struct anm_track *track, anm_time_t tm);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif	/* LIBANIM_TRACK_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/basis.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/basis.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,63 @@
+#include "basis.h"
+#include "vmath.h"
+
+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 -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/basis.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/basis.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,57 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_BASIS_H_
+#define VMATH_BASIS_H_
+
+#include "vector.h"
+#include "matrix.h"
+
+enum {
+	LEFT_HANDED,
+	RIGHT_HANDED
+};
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+void basis_matrix(mat4_t res, vec3_t i, vec3_t j, vec3_t k);
+void basis_matrix_dir(mat4_t res, vec3_t dir);
+
+#ifdef __cplusplus
+}	/* extern "C" */
+
+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;
+};
+#endif	/* __cplusplus */
+
+#endif	/* VMATH_BASIS_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/basis_c.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/basis_c.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,37 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "basis.h"
+#include "matrix.h"
+
+void basis_matrix(mat4_t res, vec3_t i, vec3_t j, vec3_t k)
+{
+	m4_identity(res);
+	m4_set_column(res, v4_cons(i.x, i.y, i.z, 1.0), 0);
+	m4_set_column(res, v4_cons(j.x, j.y, j.z, 1.0), 1);
+	m4_set_column(res, v4_cons(k.x, k.y, k.z, 1.0), 2);
+}
+
+void basis_matrix_dir(mat4_t res, vec3_t dir)
+{
+	vec3_t k = v3_normalize(dir);
+	vec3_t j = {0, 1, 0};
+	vec3_t i = v3_cross(j, k);
+	j = v3_cross(k, i);
+	basis_matrix(res, i, j, k);
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/geom.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/geom.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,150 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <math.h>
+#include "geom.h"
+#include "vector.h"
+
+plane_t plane_cons(scalar_t nx, scalar_t ny, scalar_t nz, scalar_t d)
+{
+	plane_t p;
+	p.norm.x = nx;
+	p.norm.y = ny;
+	p.norm.z = nz;
+	p.d = d;
+	return p;
+}
+
+plane_t plane_poly(vec3_t v0, vec3_t v1, vec3_t v2)
+{
+	vec3_t a, b, norm;
+
+	a = v3_sub(v1, v0);
+	b = v3_sub(v2, v0);
+	norm = v3_cross(a, b);
+	norm = v3_normalize(norm);
+
+	return plane_ptnorm(v0, norm);
+}
+
+plane_t plane_ptnorm(vec3_t pt, vec3_t normal)
+{
+	plane_t plane;
+
+	plane.norm = normal;
+	plane.d = v3_dot(pt, normal);
+
+	return plane;
+}
+
+plane_t plane_invert(plane_t p)
+{
+	p.norm = v3_neg(p.norm);
+	p.d = -p.d;
+	return p;
+}
+
+scalar_t plane_signed_dist(plane_t plane, vec3_t pt)
+{
+	vec3_t pp = plane_point(plane);
+	vec3_t pptopt = v3_sub(pt, pp);
+	return v3_dot(pptopt, plane.norm);
+}
+
+scalar_t plane_dist(plane_t plane, vec3_t pt)
+{
+	return fabs(plane_signed_dist(plane, pt));
+}
+
+vec3_t plane_point(plane_t plane)
+{
+	return v3_scale(plane.norm, plane.d);
+}
+
+int plane_ray_intersect(ray_t ray, plane_t plane, scalar_t *pos)
+{
+	vec3_t pt, orig_to_pt;
+	scalar_t ndotdir;
+
+	pt = plane_point(plane);
+	ndotdir = v3_dot(plane.norm, ray.dir);
+
+	if(fabs(ndotdir) < 1e-7) {
+		return 0;
+	}
+
+	if(pos) {
+		orig_to_pt = v3_sub(pt, ray.origin);
+		*pos = v3_dot(plane.norm, orig_to_pt) / ndotdir;
+	}
+	return 1;
+}
+
+sphere_t sphere_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t rad)
+{
+	sphere_t sph;
+	sph.pos.x = x;
+	sph.pos.y = y;
+	sph.pos.z = z;
+	sph.rad = rad;
+	return sph;
+}
+
+int sphere_ray_intersect(ray_t ray, sphere_t sph, scalar_t *pos)
+{
+	scalar_t a, b, c, d, sqrt_d, t1, t2, t;
+
+	a = v3_dot(ray.dir, ray.dir);
+	b = 2.0 * ray.dir.x * (ray.origin.x - sph.pos.x) +
+		2.0 * ray.dir.y * (ray.origin.y - sph.pos.y) +
+		2.0 * ray.dir.z * (ray.origin.z - sph.pos.z);
+	c = v3_dot(sph.pos, sph.pos) + v3_dot(ray.origin, ray.origin) +
+		2.0 * v3_dot(v3_neg(sph.pos), ray.origin) - sph.rad * sph.rad;
+
+	d = b * b - 4.0 * a * c;
+	if(d < 0.0) {
+		return 0;
+	}
+
+	sqrt_d = sqrt(d);
+	t1 = (-b + sqrt_d) / (2.0 * a);
+	t2 = (-b - sqrt_d) / (2.0 * a);
+
+	if(t1 < 1e-7 || t1 > 1.0) {
+		t1 = t2;
+	}
+	if(t2 < 1e-7 || t2 > 1.0) {
+		t2 = t1;
+	}
+	t = t1 < t2 ? t1 : t2;
+
+	if(t < 1e-7 || t > 1.0) {
+		return 0;
+	}
+
+	if(pos) {
+		*pos = t;
+	}
+	return 1;
+}
+
+int sphere_sphere_intersect(sphere_t sph1, sphere_t sph2, scalar_t *pos, scalar_t *rad)
+{
+	return -1;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/geom.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/geom.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,72 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2012 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#ifndef LIBVMATH_GEOM_H_
+#define LIBVMATH_GEOM_H_
+
+#include "vector.h"
+#include "ray.h"
+
+typedef struct {
+	vec3_t norm;
+	scalar_t d;
+} plane_t;
+
+typedef struct {
+	vec3_t pos;
+	scalar_t rad;
+} sphere_t;
+
+typedef struct {
+	vec3_t min, max;
+} aabox_t;
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* planes are good... you need planes, yes you do */
+plane_t plane_cons(scalar_t nx, scalar_t ny, scalar_t nz, scalar_t d);
+plane_t plane_poly(vec3_t v0, vec3_t v1, vec3_t v2);
+plane_t plane_ptnorm(vec3_t pt, vec3_t normal);
+
+plane_t plane_invert(plane_t p);
+
+scalar_t plane_signed_dist(plane_t plane, vec3_t pt);
+scalar_t plane_dist(plane_t plane, vec3_t pt);
+vec3_t plane_point(plane_t plane);
+
+int plane_ray_intersect(ray_t ray, plane_t plane, scalar_t *pos);
+
+/* spheres always come in handy */
+sphere_t sphere_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t rad);
+
+int sphere_ray_intersect(ray_t ray, sphere_t sph, scalar_t *pos);
+int sphere_sphere_intersect(sphere_t sph1, sphere_t sph2, scalar_t *pos, scalar_t *rad);
+
+#ifdef __cplusplus
+}
+
+/* TODO
+class Plane : public plane_t {
+public:
+};
+*/
+
+#endif
+
+#endif	/* LIBVMATH_GEOM_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/matrix.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/matrix.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,821 @@
+#include <cstdio>
+#include <cmath>
+#include "matrix.h"
+#include "vector.h"
+#include "quat.h"
+
+using namespace std;
+
+// ----------- Matrix3x3 --------------
+
+Matrix3x3 Matrix3x3::identity = Matrix3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
+
+Matrix3x3::Matrix3x3()
+{
+	*this = Matrix3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
+}
+
+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 Vector3 &ivec, const Vector3 &jvec, const Vector3 &kvec)
+{
+	set_row_vector(ivec, 0);
+	set_row_vector(jvec, 1);
+	set_row_vector(kvec, 2);
+}
+
+Matrix3x3::Matrix3x3(const mat3_t cmat)
+{
+	memcpy(m, cmat, sizeof(mat3_t));
+}
+
+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;
+}
+
+// Algorithm in Ken Shoemake's article in 1987 SIGGRAPH course notes
+// article "Quaternion Calculus and Fast Animation".
+// adapted from: http://www.geometrictools.com/LibMathematics/Algebra/Wm5Quaternion.inl
+Quaternion Matrix3x3::get_rotation_quat() const
+{
+	static const int next[3] = {1, 2, 0};
+
+	float quat[4];
+
+	scalar_t trace = m[0][0] + m[1][1] + m[2][2];
+	scalar_t root;
+
+	if(trace > 0.0f) {
+		// |w| > 1/2
+		root = sqrt(trace + 1.0f);	// 2w
+		quat[0] = 0.5f * root;
+		root = 0.5f / root;	// 1 / 4w
+		quat[1] = (m[2][1] - m[1][2]) * root;
+		quat[2] = (m[0][2] - m[2][0]) * root;
+		quat[3] = (m[1][0] - m[0][1]) * root;
+	} else {
+		// |w| <= 1/2
+		int i = 0;
+		if(m[1][1] > m[0][0]) {
+			i = 1;
+		}
+		if(m[2][2] > m[i][i]) {
+			i = 2;
+		}
+		int j = next[i];
+		int k = next[j];
+
+		root = sqrt(m[i][i] - m[j][j] - m[k][k] + 1.0f);
+		quat[i + 1] = 0.5f * root;
+		root = 0.5f / root;
+		quat[0] = (m[k][j] - m[j][k]) * root;
+		quat[j + 1] = (m[j][i] - m[i][j]) * root;
+		quat[k + 1] = (m[k][i] - m[i][k]) * root;
+	}
+	return Quaternion(quat[0], quat[1], quat[2], quat[3]);
+}
+
+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 = Matrix4x4(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
+
+Matrix4x4::Matrix4x4()
+{
+	*this = Matrix4x4(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
+}
+
+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;
+}
+
+Matrix4x4::Matrix4x4(const mat4_t cmat)
+{
+	memcpy(m, cmat, sizeof(mat4_t));
+}
+
+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];
+		}
+	}
+}
+
+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;
+}
+
+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++;
+	}
+}
+
+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);
+}
+
+Vector3 Matrix4x4::get_translation() const
+{
+	return Vector3(m[0][3], m[1][3], m[2][3]);
+}
+
+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;
+
+	Matrix4x4 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 *= 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::rotate(const Quaternion &quat)
+{
+	*this *= quat.get_rotation_matrix();
+}
+
+void Matrix4x4::set_rotation(const Quaternion &quat)
+{
+	*this = quat.get_rotation_matrix();
+}
+
+Quaternion Matrix4x4::get_rotation_quat() const
+{
+	Matrix3x3 mat3 = *this;
+	return mat3.get_rotation_quat();
+}
+
+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);
+}
+
+Vector3 Matrix4x4::get_scaling() const
+{
+	Vector3 vi = get_row_vector(0);
+	Vector3 vj = get_row_vector(1);
+	Vector3 vk = get_row_vector(2);
+
+	return Vector3(vi.length(), vj.length(), vk.length());
+}
+
+void Matrix4x4::set_perspective(float vfov, float aspect, float znear, float zfar)
+{
+	float f = 1.0f / tan(vfov * 0.5f);
+    float dz = znear - zfar;
+
+	reset_identity();
+
+	m[0][0] = f / aspect;
+    m[1][1] = f;
+    m[2][2] = (zfar + znear) / dz;
+    m[3][2] = -1.0f;
+    m[2][3] = 2.0f * zfar * znear / dz;
+    m[3][3] = 0.0f;
+}
+
+void Matrix4x4::set_orthographic(float left, float right, float bottom, float top, float znear, float zfar)
+{
+	float dx = right - left;
+	float dy = top - bottom;
+	float dz = zfar - znear;
+
+	reset_identity();
+
+	m[0][0] = 2.0 / dx;
+	m[1][1] = 2.0 / dy;
+	m[2][2] = -2.0 / dz;
+	m[0][3] = -(right + left) / dx;
+	m[1][3] = -(top + bottom) / dy;
+	m[2][3] = -(zfar + znear) / dz;
+}
+
+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 adj = adjoint();
+
+	return adj * (1.0f / determinant());
+}
+
+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 -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/matrix.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/matrix.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,260 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_MATRIX_H_
+#define VMATH_MATRIX_H_
+
+#include <stdio.h>
+#include "vmath_types.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+/* C matrix 3x3 functions */
+static inline void m3_identity(mat3_t m);
+static inline void m3_cons(mat3_t m,
+		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);
+static inline void m3_copy(mat3_t dest, mat3_t src);
+void m3_to_m4(mat4_t dest, mat3_t src);
+
+void m3_print(FILE *fp, mat3_t m);
+
+/* C matrix 4x4 functions */
+static inline void m4_identity(mat4_t m);
+static inline void m4_cons(mat4_t m,
+		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);
+static inline void m4_copy(mat4_t dest, mat4_t src);
+void m4_to_m3(mat3_t dest, mat4_t src);
+
+static inline void m4_mult(mat4_t res, mat4_t m1, mat4_t m2);
+
+void m4_set_translation(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
+void m4_translate(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
+
+void m4_rotate(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
+
+void m4_set_rotation_x(mat4_t m, scalar_t angle);
+void m4_rotate_x(mat4_t m, scalar_t angle);
+void m4_set_rotation_y(mat4_t m, scalar_t angle);
+void m4_rotate_y(mat4_t m, scalar_t angle);
+void m4_set_rotation_z(mat4_t m, scalar_t angle);
+void m4_rotate_z(mat4_t m, scalar_t angle);
+/* axis-angle rotation */
+void m4_set_rotation_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z);
+void m4_rotate_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z);
+/* concatentate a rotation quaternion */
+void m4_rotate_quat(mat4_t m, quat_t q);
+
+void m4_set_scaling(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
+void m4_scale(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
+
+static inline void m4_set_column(mat4_t m, vec4_t v, int idx);
+static inline void m4_set_row(mat4_t m, vec4_t v, int idx);
+
+void m4_transpose(mat4_t res, mat4_t m);
+scalar_t m4_determinant(mat4_t m);
+void m4_adjoint(mat4_t res, mat4_t m);
+void m4_inverse(mat4_t res, mat4_t m);
+
+void m4_print(FILE *fp, mat4_t m);
+
+#ifdef __cplusplus
+}
+
+/* when included from C++ source files, also define the matrix classes */
+#include <iostream>
+
+/** 3x3 matrix */
+class Matrix3x3 {
+public:
+	scalar_t m[3][3];
+
+	static Matrix3x3 identity;
+
+	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 Vector3 &ivec, const Vector3 &jvec, const Vector3 &kvec);
+	Matrix3x3(const mat3_t cmat);
+
+	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);
+	Quaternion get_rotation_quat() const;
+
+	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);
+};
+
+/* binary operations matrix (op) matrix */
+Matrix3x3 operator +(const Matrix3x3 &m1, const Matrix3x3 &m2);
+Matrix3x3 operator -(const Matrix3x3 &m1, const Matrix3x3 &m2);
+Matrix3x3 operator *(const Matrix3x3 &m1, const Matrix3x3 &m2);
+
+void operator +=(Matrix3x3 &m1, const Matrix3x3 &m2);
+void operator -=(Matrix3x3 &m1, const Matrix3x3 &m2);
+void operator *=(Matrix3x3 &m1, const Matrix3x3 &m2);
+
+/* binary operations matrix (op) scalar and scalar (op) matrix */
+Matrix3x3 operator *(const Matrix3x3 &mat, scalar_t scalar);
+Matrix3x3 operator *(scalar_t scalar, const Matrix3x3 &mat);
+
+void operator *=(Matrix3x3 &mat, scalar_t scalar);
+
+std::ostream &operator <<(std::ostream &out, const Matrix3x3 &mat);
+
+
+
+/** 4x4 matrix */
+class Matrix4x4 {
+public:
+	scalar_t m[4][4];
+
+	static Matrix4x4 identity;
+
+	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 mat4_t cmat);
+
+	Matrix4x4(const Matrix3x3 &mat3x3);
+
+	/* 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 inline 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);
+	Vector3 get_translation() const;	/* extract translation */
+
+	void rotate(const Vector3 &euler_angles);			/* 3d rotation with euler angles */
+	void rotate(const Vector3 &axis, scalar_t angle);	/* 3d axis/angle rotation */
+	void rotate(const Quaternion &quat);
+	void set_rotation(const Vector3 &euler_angles);
+	void set_rotation(const Vector3 &axis, scalar_t angle);
+	void set_rotation(const Quaternion &quat);
+	Quaternion get_rotation_quat() const;		/* extract rotation */
+
+	void scale(const Vector4 &scale_vec);
+	void set_scaling(const Vector4 &scale_vec);
+	Vector3 get_scaling() const;		/* extract scaling */
+
+	void set_perspective(float vfov, float aspect, float znear, float zfar);
+	void set_orthographic(float left, float right, float bottom, float top, float znear = -1.0, float zfar = 1.0);
+
+	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;
+
+	friend std::ostream &operator <<(std::ostream &out, const Matrix4x4 &mat);
+};
+
+/* binary operations matrix (op) matrix */
+Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2);
+Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2);
+inline Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2);
+
+void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2);
+void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2);
+inline void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2);
+
+/* binary operations matrix (op) scalar and scalar (op) matrix */
+Matrix4x4 operator *(const Matrix4x4 &mat, scalar_t scalar);
+Matrix4x4 operator *(scalar_t scalar, const Matrix4x4 &mat);
+
+void operator *=(Matrix4x4 &mat, scalar_t scalar);
+
+std::ostream &operator <<(std::ostream &out, const Matrix4x4 &mat);
+
+#endif	/* __cplusplus */
+
+#include "matrix.inl"
+
+#endif	/* VMATH_MATRIX_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/matrix.inl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/matrix.inl	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,200 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <string.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+/* C matrix 3x3 functions */
+static inline void m3_identity(mat3_t m)
+{
+	static const mat3_t id = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+	memcpy(m, id, sizeof id);
+}
+
+static inline void m3_cons(mat3_t m,
+		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;
+}
+
+static inline void m3_copy(mat3_t dest, mat3_t src)
+{
+	memcpy(dest, src, sizeof(mat3_t));
+}
+
+
+/* C matrix 4x4 functions */
+static inline void m4_identity(mat4_t m)
+{
+	static const mat4_t id = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
+	memcpy(m, id, sizeof id);
+}
+
+static inline void m4_cons(mat4_t m,
+		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;
+}
+
+static inline void m4_copy(mat4_t dest, mat4_t src)
+{
+	memcpy(dest, src, sizeof(mat4_t));
+}
+
+static inline void m4_mult(mat4_t res, mat4_t m1, mat4_t m2)
+{
+	mat4_t tmp;
+
+	/*
+	int i, j;
+	for(i=0; i<4; i++) {
+		for(j=0; j<4; j++) {
+			tmp[i][j] = m1[i][0] * m2[0][j] + m1[i][1] * m2[1][j] + m1[i][2] * m2[2][j] + m1[i][3] * m2[3][j];
+		}
+	}
+	*/
+
+	tmp[0][0] = m1[0][0] * m2[0][0] + m1[0][1] * m2[1][0] + m1[0][2] * m2[2][0] + m1[0][3] * m2[3][0];
+	tmp[0][1] = m1[0][0] * m2[0][1] + m1[0][1] * m2[1][1] + m1[0][2] * m2[2][1] + m1[0][3] * m2[3][1];
+	tmp[0][2] = m1[0][0] * m2[0][2] + m1[0][1] * m2[1][2] + m1[0][2] * m2[2][2] + m1[0][3] * m2[3][2];
+	tmp[0][3] = m1[0][0] * m2[0][3] + m1[0][1] * m2[1][3] + m1[0][2] * m2[2][3] + m1[0][3] * m2[3][3];
+
+	tmp[1][0] = m1[1][0] * m2[0][0] + m1[1][1] * m2[1][0] + m1[1][2] * m2[2][0] + m1[1][3] * m2[3][0];
+	tmp[1][1] = m1[1][0] * m2[0][1] + m1[1][1] * m2[1][1] + m1[1][2] * m2[2][1] + m1[1][3] * m2[3][1];
+	tmp[1][2] = m1[1][0] * m2[0][2] + m1[1][1] * m2[1][2] + m1[1][2] * m2[2][2] + m1[1][3] * m2[3][2];
+	tmp[1][3] = m1[1][0] * m2[0][3] + m1[1][1] * m2[1][3] + m1[1][2] * m2[2][3] + m1[1][3] * m2[3][3];
+
+	tmp[2][0] = m1[2][0] * m2[0][0] + m1[2][1] * m2[1][0] + m1[2][2] * m2[2][0] + m1[2][3] * m2[3][0];
+	tmp[2][1] = m1[2][0] * m2[0][1] + m1[2][1] * m2[1][1] + m1[2][2] * m2[2][1] + m1[2][3] * m2[3][1];
+	tmp[2][2] = m1[2][0] * m2[0][2] + m1[2][1] * m2[1][2] + m1[2][2] * m2[2][2] + m1[2][3] * m2[3][2];
+	tmp[2][3] = m1[2][0] * m2[0][3] + m1[2][1] * m2[1][3] + m1[2][2] * m2[2][3] + m1[2][3] * m2[3][3];
+
+	tmp[3][0] = m1[3][0] * m2[0][0] + m1[3][1] * m2[1][0] + m1[3][2] * m2[2][0] + m1[3][3] * m2[3][0];
+	tmp[3][1] = m1[3][0] * m2[0][1] + m1[3][1] * m2[1][1] + m1[3][2] * m2[2][1] + m1[3][3] * m2[3][1];
+	tmp[3][2] = m1[3][0] * m2[0][2] + m1[3][1] * m2[1][2] + m1[3][2] * m2[2][2] + m1[3][3] * m2[3][2];
+	tmp[3][3] = m1[3][0] * m2[0][3] + m1[3][1] * m2[1][3] + m1[3][2] * m2[2][3] + m1[3][3] * m2[3][3];
+
+	m4_copy(res, tmp);
+}
+
+static inline void m4_set_column(mat4_t m, vec4_t v, int idx)
+{
+	m[0][idx] = v.x;
+	m[1][idx] = v.y;
+	m[2][idx] = v.z;
+	m[3][idx] = v.w;
+}
+
+static inline void m4_set_row(mat4_t m, vec4_t v, int idx)
+{
+	m[idx][0] = v.x;
+	m[idx][1] = v.y;
+	m[idx][2] = v.z;
+	m[idx][3] = v.w;
+}
+
+#ifdef __cplusplus
+}	/* extern "C" */
+
+
+/* unrolled to hell and inline */
+inline Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2)
+{
+	Matrix4x4 res;
+
+	/*
+	for(i=0; i<4; i++) {
+		for(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];
+		}
+	}
+	*/
+
+	res.m[0][0] = m1.m[0][0] * m2.m[0][0] + m1.m[0][1] * m2.m[1][0] + m1.m[0][2] * m2.m[2][0] + m1.m[0][3] * m2.m[3][0];
+	res.m[0][1] = m1.m[0][0] * m2.m[0][1] + m1.m[0][1] * m2.m[1][1] + m1.m[0][2] * m2.m[2][1] + m1.m[0][3] * m2.m[3][1];
+	res.m[0][2] = m1.m[0][0] * m2.m[0][2] + m1.m[0][1] * m2.m[1][2] + m1.m[0][2] * m2.m[2][2] + m1.m[0][3] * m2.m[3][2];
+	res.m[0][3] = m1.m[0][0] * m2.m[0][3] + m1.m[0][1] * m2.m[1][3] + m1.m[0][2] * m2.m[2][3] + m1.m[0][3] * m2.m[3][3];
+
+	res.m[1][0] = m1.m[1][0] * m2.m[0][0] + m1.m[1][1] * m2.m[1][0] + m1.m[1][2] * m2.m[2][0] + m1.m[1][3] * m2.m[3][0];
+	res.m[1][1] = m1.m[1][0] * m2.m[0][1] + m1.m[1][1] * m2.m[1][1] + m1.m[1][2] * m2.m[2][1] + m1.m[1][3] * m2.m[3][1];
+	res.m[1][2] = m1.m[1][0] * m2.m[0][2] + m1.m[1][1] * m2.m[1][2] + m1.m[1][2] * m2.m[2][2] + m1.m[1][3] * m2.m[3][2];
+	res.m[1][3] = m1.m[1][0] * m2.m[0][3] + m1.m[1][1] * m2.m[1][3] + m1.m[1][2] * m2.m[2][3] + m1.m[1][3] * m2.m[3][3];
+
+	res.m[2][0] = m1.m[2][0] * m2.m[0][0] + m1.m[2][1] * m2.m[1][0] + m1.m[2][2] * m2.m[2][0] + m1.m[2][3] * m2.m[3][0];
+	res.m[2][1] = m1.m[2][0] * m2.m[0][1] + m1.m[2][1] * m2.m[1][1] + m1.m[2][2] * m2.m[2][1] + m1.m[2][3] * m2.m[3][1];
+	res.m[2][2] = m1.m[2][0] * m2.m[0][2] + m1.m[2][1] * m2.m[1][2] + m1.m[2][2] * m2.m[2][2] + m1.m[2][3] * m2.m[3][2];
+	res.m[2][3] = m1.m[2][0] * m2.m[0][3] + m1.m[2][1] * m2.m[1][3] + m1.m[2][2] * m2.m[2][3] + m1.m[2][3] * m2.m[3][3];
+
+	res.m[3][0] = m1.m[3][0] * m2.m[0][0] + m1.m[3][1] * m2.m[1][0] + m1.m[3][2] * m2.m[2][0] + m1.m[3][3] * m2.m[3][0];
+	res.m[3][1] = m1.m[3][0] * m2.m[0][1] + m1.m[3][1] * m2.m[1][1] + m1.m[3][2] * m2.m[2][1] + m1.m[3][3] * m2.m[3][1];
+	res.m[3][2] = m1.m[3][0] * m2.m[0][2] + m1.m[3][1] * m2.m[1][2] + m1.m[3][2] * m2.m[2][2] + m1.m[3][3] * m2.m[3][2];
+	res.m[3][3] = m1.m[3][0] * m2.m[0][3] + m1.m[3][1] * m2.m[1][3] + m1.m[3][2] * m2.m[2][3] + m1.m[3][3] * m2.m[3][3];
+
+	return res;
+}
+
+inline void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2)
+{
+	Matrix4x4 res = m1 * m2;
+	m1 = res;
+}
+
+
+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()
+{
+	*this = identity;
+}
+
+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()
+{
+	*this = identity;
+}
+#endif	/* __cplusplus */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/matrix_c.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/matrix_c.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,292 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <stdio.h>
+#include "matrix.h"
+#include "vector.h"
+#include "quat.h"
+
+void m3_to_m4(mat4_t dest, mat3_t src)
+{
+	int i, j;
+
+	memset(dest, 0, sizeof(mat4_t));
+	for(i=0; i<3; i++) {
+		for(j=0; j<3; j++) {
+			dest[i][j] = src[i][j];
+		}
+	}
+	dest[3][3] = 1.0;
+}
+
+void m3_print(FILE *fp, mat3_t m)
+{
+	int i;
+	for(i=0; i<3; i++) {
+		fprintf(fp, "[ %12.5f %12.5f %12.5f ]\n", (float)m[i][0], (float)m[i][1], (float)m[i][2]);
+	}
+}
+
+/* C matrix 4x4 functions */
+void m4_to_m3(mat3_t dest, mat4_t src)
+{
+	int i, j;
+	for(i=0; i<3; i++) {
+		for(j=0; j<3; j++) {
+			dest[i][j] = src[i][j];
+		}
+	}
+}
+
+void m4_set_translation(mat4_t m, scalar_t x, scalar_t y, scalar_t z)
+{
+	m4_identity(m);
+	m[0][3] = x;
+	m[1][3] = y;
+	m[2][3] = z;
+}
+
+void m4_translate(mat4_t m, scalar_t x, scalar_t y, scalar_t z)
+{
+	mat4_t tm;
+	m4_set_translation(tm, x, y, z);
+	m4_mult(m, m, tm);
+}
+
+void m4_rotate(mat4_t m, scalar_t x, scalar_t y, scalar_t z)
+{
+	m4_rotate_x(m, x);
+	m4_rotate_y(m, y);
+	m4_rotate_z(m, z);
+}
+
+void m4_set_rotation_x(mat4_t m, scalar_t angle)
+{
+	m4_identity(m);
+	m[1][1] = cos(angle); m[1][2] = -sin(angle);
+	m[2][1] = sin(angle); m[2][2] = cos(angle);
+}
+
+void m4_rotate_x(mat4_t m, scalar_t angle)
+{
+	mat4_t rm;
+	m4_set_rotation_x(m, angle);
+	m4_mult(m, m, rm);
+}
+
+void m4_set_rotation_y(mat4_t m, scalar_t angle)
+{
+	m4_identity(m);
+	m[0][0] = cos(angle); m[0][2] = sin(angle);
+	m[2][0] = -sin(angle); m[2][2] = cos(angle);
+}
+
+void m4_rotate_y(mat4_t m, scalar_t angle)
+{
+	mat4_t rm;
+	m4_set_rotation_y(rm, angle);
+	m4_mult(m, m, rm);
+}
+
+void m4_set_rotation_z(mat4_t m, scalar_t angle)
+{
+	m4_identity(m);
+	m[0][0] = cos(angle); m[0][1] = -sin(angle);
+	m[1][0] = sin(angle); m[1][1] = cos(angle);
+}
+
+void m4_rotate_z(mat4_t m, scalar_t angle)
+{
+	mat4_t rm;
+	m4_set_rotation_z(rm, angle);
+	m4_mult(m, m, rm);
+}
+
+void m4_set_rotation_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z)
+{
+	scalar_t sina = sin(angle);
+	scalar_t cosa = cos(angle);
+	scalar_t one_minus_cosa = 1.0 - cosa;
+	scalar_t nxsq = x * x;
+	scalar_t nysq = y * y;
+	scalar_t nzsq = z * z;
+
+	m[0][0] = nxsq + (1.0 - nxsq) * cosa;
+	m[0][1] = x * y * one_minus_cosa - z * sina;
+	m[0][2] = x * z * one_minus_cosa + y * sina;
+	m[1][0] = x * y * one_minus_cosa + z * sina;
+	m[1][1] = nysq + (1.0 - nysq) * cosa;
+	m[1][2] = y * z * one_minus_cosa - x * sina;
+	m[2][0] = x * z * one_minus_cosa - y * sina;
+	m[2][1] = y * z * one_minus_cosa + x * sina;
+	m[2][2] = nzsq + (1.0 - nzsq) * cosa;
+
+	/* the rest are identity */
+	m[3][0] = m[3][1] = m[3][2] = m[0][3] = m[1][3] = m[2][3] = 0.0;
+	m[3][3] = 1.0;
+}
+
+void m4_rotate_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z)
+{
+	mat4_t xform;
+	m4_set_rotation_axis(xform, angle, x, y, z);
+	m4_mult(m, m, xform);
+}
+
+void m4_rotate_quat(mat4_t m, quat_t q)
+{
+	mat4_t rm;
+	quat_to_mat4(rm, q);
+	m4_mult(m, m, rm);
+}
+
+void m4_scale(mat4_t m, scalar_t x, scalar_t y, scalar_t z)
+{
+	mat4_t sm;
+	m4_identity(sm);
+	sm[0][0] = x;
+	sm[1][1] = y;
+	sm[2][2] = z;
+	m4_mult(m, m, sm);
+}
+
+void m4_transpose(mat4_t res, mat4_t m)
+{
+	int i, j;
+	mat4_t tmp;
+	m4_copy(tmp, m);
+
+	for(i=0; i<4; i++) {
+		for(j=0; j<4; j++) {
+			res[i][j] = tmp[j][i];
+		}
+	}
+}
+
+scalar_t m4_determinant(mat4_t m)
+{
+	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;
+}
+
+void m4_adjoint(mat4_t res, mat4_t m)
+{
+	int i, j;
+	mat4_t coef;
+
+	coef[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[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[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[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[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[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[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[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[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[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[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[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[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[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[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[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]));
+
+	m4_transpose(res, coef);
+
+	for(i=0; i<4; i++) {
+		for(j=0; j<4; j++) {
+			res[i][j] = j % 2 ? -res[i][j] : res[i][j];
+			if(i % 2) res[i][j] = -res[i][j];
+		}
+	}
+}
+
+void m4_inverse(mat4_t res, mat4_t m)
+{
+	int i, j;
+	mat4_t adj;
+	scalar_t det;
+
+	m4_adjoint(adj, m);
+	det = m4_determinant(m);
+	
+	for(i=0; i<4; i++) {
+		for(j=0; j<4; j++) {
+			res[i][j] = adj[i][j] / det;
+		}
+	}
+}
+
+void m4_print(FILE *fp, mat4_t m)
+{
+	int i;
+	for(i=0; i<4; i++) {
+		fprintf(fp, "[ %12.5f %12.5f %12.5f %12.5f ]\n", (float)m[i][0], (float)m[i][1], (float)m[i][2], (float)m[i][3]);
+	}
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/quat.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/quat.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,212 @@
+#include "quat.h"
+#include "vmath.h"
+
+Quaternion::Quaternion()
+{
+	s = 1.0;
+	v.x = v.y = v.z = 0.0;
+}
+
+Quaternion::Quaternion(scalar_t s, const Vector3 &v)
+{
+	this->s = s;
+	this->v = v;
+}
+
+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(const Vector3 &axis, scalar_t angle)
+{
+	set_rotation(axis, angle);
+}
+
+Quaternion::Quaternion(const quat_t &quat)
+{
+	v.x = quat.x;
+	v.y = quat.y;
+	v.z = quat.z;
+	s = quat.w;
+}
+
+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 half_angle = angle / 2.0;
+	s = cos(half_angle);
+	v = axis * sin(half_angle);
+}
+
+void Quaternion::rotate(const Vector3 &axis, scalar_t angle)
+{
+	Quaternion q;
+	scalar_t half_angle = angle / 2.0;
+	q.s = cos(half_angle);
+	q.v = axis * sin(half_angle);
+
+	*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 &quat1, const Quaternion &q2, scalar_t t)
+{
+	Quaternion q1;
+	scalar_t dot = q1.s * q2.s + q1.v.x * q2.v.x + q1.v.y * q2.v.y + q1.v.z * q2.v.z;
+
+	if(dot < 0.0) {
+		/* make sure we interpolate across the shortest arc */
+		q1 = -quat1;
+		dot = -dot;
+	} else {
+		q1 = quat1;
+	}
+
+	/* clamp dot to [-1, 1] in order to avoid domain errors in acos due to
+	 * floating point imprecisions
+	 */
+	if(dot < -1.0) dot = -1.0;
+	if(dot > 1.0) dot = 1.0;
+
+	scalar_t angle = acos(dot);
+	scalar_t a, b;
+
+	scalar_t sin_angle = sin(angle);
+	if(fabs(sin_angle) < SMALL_NUMBER) {
+		/* for very small angles or completely opposite orientations
+		 * use linear interpolation to avoid div/zero (in the first case it makes sense,
+		 * the second case is pretty much undefined anyway I guess ...
+		 */
+		a = 1.0f - t;
+		b = t;
+	} else {
+		a = sin((1.0f - t) * angle) / sin_angle;
+		b = sin(t * angle) / sin_angle;
+	}
+
+	scalar_t x = q1.v.x * a + q2.v.x * b;
+	scalar_t y = q1.v.y * a + q2.v.y * b;
+	scalar_t z = q1.v.z * a + q2.v.z * b;
+	scalar_t s = q1.s * a + q2.s * b;
+
+	return Quaternion(s, Vector3(x, y, z));
+}
+
+
+std::ostream &operator <<(std::ostream &out, const Quaternion &q)
+{
+	out << "(" << q.s << ", " << q.v << ")";
+	return out;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/quat.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/quat.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,118 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_QUATERNION_H_
+#define VMATH_QUATERNION_H_
+
+#include <stdio.h>
+#include "vmath_types.h"
+#include "vector.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+#define quat_cons(s, x, y, z)	v4_cons(x, y, z, s)
+#define quat_vec(q)				v3_cons((q).x, (q).y, (q).z)
+#define quat_s(q)				((q).w)
+#define quat_identity()			quat_cons(1.0, 0.0, 0.0, 0.0)
+void quat_print(FILE *fp, quat_t q);
+
+#define quat_add		v4_add
+#define quat_sub		v4_sub
+#define quat_neg		v4_neg
+
+static inline quat_t quat_mul(quat_t q1, quat_t q2);
+
+static inline quat_t quat_conjugate(quat_t q);
+
+#define quat_length		v4_length
+#define quat_length_sq	v4_length_sq
+
+#define quat_normalize	v4_normalize
+static inline quat_t quat_inverse(quat_t q);
+
+quat_t quat_rotate(quat_t q, scalar_t angle, scalar_t x, scalar_t y, scalar_t z);
+quat_t quat_rotate_quat(quat_t q, quat_t rotq);
+
+static inline void quat_to_mat3(mat3_t res, quat_t q);
+static inline void quat_to_mat4(mat4_t res, quat_t q);
+
+#define quat_lerp quat_slerp
+quat_t quat_slerp(quat_t q1, quat_t q2, scalar_t t);
+
+
+#ifdef __cplusplus
+}	/* extern "C" */
+
+#include <iostream>
+
+/* 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(const quat_t &quat);
+
+	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;
+};
+
+Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t);
+inline Quaternion lerp(const Quaternion &q1, const Quaternion &q2, scalar_t t);
+
+std::ostream &operator <<(std::ostream &out, const Quaternion &q);
+
+#endif	/* __cplusplus */
+
+#include "quat.inl"
+
+#endif	/* VMATH_QUATERNION_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/quat.inl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/quat.inl	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,81 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "vector.h"
+#include "matrix.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+static inline quat_t quat_mul(quat_t q1, quat_t q2)
+{
+	quat_t res;
+	vec3_t v1 = quat_vec(q1);
+	vec3_t v2 = quat_vec(q2);
+
+	res.w = q1.w * q2.w - v3_dot(v1, v2);
+	/* resvec = v2 * q1 + v1 * q2 + cross(v1, v2) */
+	res.x = v2.x * q1.w + v1.x * q2.w + (v1.y * v2.z - v1.z * v2.y);
+	res.y = v2.y * q1.w + v1.y * q2.w + (v1.z * v2.x - v1.x * v2.z);
+	res.z = v2.z * q1.w + v1.z * q2.w + (v1.x * v2.y - v1.y * v2.x);
+	return res;
+}
+
+static inline quat_t quat_conjugate(quat_t q)
+{
+	q.x = -q.x;
+	q.y = -q.y;
+	q.z = -q.z;
+	return q;
+}
+
+static inline quat_t quat_inverse(quat_t q)
+{
+	scalar_t lensq = quat_length_sq(q);
+	q = quat_conjugate(q);
+	q.x /= lensq;
+	q.y /= lensq;
+	q.z /= lensq;
+	q.w /= lensq;
+	return q;
+}
+
+static inline void quat_to_mat3(mat3_t res, quat_t q)
+{
+	m3_cons(res, 1.0 - 2.0 * q.y*q.y - 2.0 * q.z*q.z, 2.0 * q.x * q.y - 2.0 * q.w * q.z,   2.0 * q.z * q.x + 2.0 * q.w * q.y,
+				 2.0 * q.x * q.y + 2.0 * q.w * q.z,   1.0 - 2.0 * q.x*q.x - 2.0 * q.z*q.z, 2.0 * q.y * q.z - 2.0 * q.w * q.x,
+				 2.0 * q.z * q.x - 2.0 * q.w * q.y,   2.0 * q.y * q.z + 2.0 * q.w * q.x,   1.0 - 2.0 * q.x*q.x - 2.0 * q.y*q.y);
+}
+
+static inline void quat_to_mat4(mat4_t res, quat_t q)
+{
+	m4_cons(res, 1.0 - 2.0 * q.y*q.y - 2.0 * q.z*q.z, 2.0 * q.x * q.y - 2.0 * q.w * q.z,   2.0 * q.z * q.x + 2.0 * q.w * q.y,   0,
+				 2.0 * q.x * q.y + 2.0 * q.w * q.z,   1.0 - 2.0 * q.x*q.x - 2.0 * q.z*q.z, 2.0 * q.y * q.z - 2.0 * q.w * q.x,   0,
+				 2.0 * q.z * q.x - 2.0 * q.w * q.y,   2.0 * q.y * q.z + 2.0 * q.w * q.x,   1.0 - 2.0 * q.x*q.x - 2.0 * q.y*q.y, 0,
+				 0, 0, 0, 1);
+}
+
+#ifdef __cplusplus
+}	/* extern "C" */
+
+inline Quaternion lerp(const Quaternion &a, const Quaternion &b, scalar_t t)
+{
+	return slerp(a, b, t);
+}
+#endif	/* __cplusplus */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/quat_c.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/quat_c.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,89 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <stdio.h>
+#include <math.h>
+#include "quat.h"
+
+void quat_print(FILE *fp, quat_t q)
+{
+	fprintf(fp, "([ %.4f %.4f %.4f ] %.4f)", q.x, q.y, q.z, q.w);
+}
+
+quat_t quat_rotate(quat_t q, scalar_t angle, scalar_t x, scalar_t y, scalar_t z)
+{
+	quat_t rq;
+	scalar_t half_angle = angle * 0.5;
+	scalar_t sin_half = sin(half_angle);
+
+	rq.w = cos(half_angle);
+	rq.x = x * sin_half;
+	rq.y = y * sin_half;
+	rq.z = z * sin_half;
+
+	return quat_mul(q, rq);
+}
+
+quat_t quat_rotate_quat(quat_t q, quat_t rotq)
+{
+	return quat_mul(quat_mul(rotq, q), quat_conjugate(rotq));
+}
+
+quat_t quat_slerp(quat_t q1, quat_t q2, scalar_t t)
+{
+	quat_t res;
+	scalar_t a, b, angle, sin_angle, dot;
+
+	dot = q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z;
+	if(dot < 0.0) {
+		/* make sure we interpolate across the shortest arc */
+		q1.x = -q1.x;
+		q1.y = -q1.y;
+		q1.z = -q1.z;
+		q1.w = -q1.w;
+		dot = -dot;
+	}
+
+	/* clamp dot to [-1, 1] in order to avoid domain errors in acos due to
+	 * floating point imprecisions
+	 */
+	if(dot < -1.0) dot = -1.0;
+	if(dot > 1.0) dot = 1.0;
+
+	angle = acos(dot);
+	sin_angle = sin(angle);
+
+	if(fabs(sin_angle) < SMALL_NUMBER) {
+		/* for very small angles or completely opposite orientations
+		 * use linear interpolation to avoid div/zero (in the first case it makes sense,
+		 * the second case is pretty much undefined anyway I guess ...
+		 */
+		a = 1.0f - t;
+		b = t;
+	} else {
+		a = sin((1.0f - t) * angle) / sin_angle;
+		b = sin(t * angle) / sin_angle;
+	}
+
+	res.x = q1.x * a + q2.x * b;
+	res.y = q1.y * a + q2.y * b;
+	res.z = q1.z * a + q2.z * b;
+	res.w = q1.w * a + q2.w * b;
+	return res;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/ray.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/ray.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,39 @@
+#include "ray.h"
+#include "vector.h"
+
+scalar_t Ray::env_ior = 1.0;
+
+Ray::Ray()
+{
+	ior = env_ior;
+	energy = 1.0;
+	time = 0;
+	iter = 0;
+}
+
+Ray::Ray(const Vector3 &origin, const Vector3 &dir)
+{
+	this->origin = origin;
+	this->dir = dir;
+	ior = env_ior;
+	energy = 1.0;
+	time = 0;
+	iter = 0;
+}
+
+void Ray::transform(const Matrix4x4 &xform)
+{
+	Matrix4x4 upper = xform;
+	upper[0][3] = upper[1][3] = upper[2][3] = upper[3][0] = upper[3][1] = upper[3][2] = 0.0;
+	upper[3][3] = 1.0;
+
+	dir.transform(upper);
+	origin.transform(xform);
+}
+
+Ray Ray::transformed(const Matrix4x4 &xform) const
+{
+	Ray foo = *this;
+	foo.transform(xform);
+	return foo;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/ray.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/ray.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,64 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_RAY_H_
+#define VMATH_RAY_H_
+
+#include "matrix.h"
+#include "vector.h"
+
+typedef struct {
+	vec3_t origin, dir;
+} ray_t;
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+static inline ray_t ray_cons(vec3_t origin, vec3_t dir);
+ray_t ray_transform(ray_t r, mat4_t m);
+
+#ifdef __cplusplus
+}	/* __cplusplus */
+
+#include <stack>
+
+class Ray {
+public:
+	static scalar_t env_ior;
+
+	Vector3 origin, dir;
+	scalar_t energy;
+	int iter;
+	scalar_t ior;
+	long time;
+
+	Ray();
+	Ray(const Vector3 &origin, const Vector3 &dir);
+
+	void transform(const Matrix4x4 &xform);
+	Ray transformed(const Matrix4x4 &xform) const;
+};
+
+inline Ray reflect_ray(const Ray &inray, const Vector3 &norm);
+inline Ray refract_ray(const Ray &inray, const Vector3 &norm, scalar_t from_ior, scalar_t to_ior);
+#endif	/* __cplusplus */
+
+#include "ray.inl"
+
+#endif	/* VMATH_RAY_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/ray.inl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/ray.inl	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,52 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+static inline ray_t ray_cons(vec3_t origin, vec3_t dir)
+{
+	ray_t r;
+	r.origin = origin;
+	r.dir = dir;
+	return r;
+}
+
+#ifdef __cplusplus
+}
+
+inline Ray reflect_ray(const Ray &inray, const Vector3 &norm)
+{
+	Ray ray = inray;
+	ray.dir = ray.dir.reflection(norm);
+	return ray;
+}
+
+inline Ray refract_ray(const Ray &inray, const Vector3 &norm, scalar_t from_ior, scalar_t to_ior)
+{
+	Ray ray = inray;
+	ray.dir = ray.dir.refraction(norm, from_ior, to_ior);
+
+	/* check TIR */
+	if(dot_product(ray.dir, norm) > 0.0) {
+		return reflect_ray(inray, norm);
+	}
+	return ray;
+}
+#endif	/* __cplusplus */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/ray_c.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/ray_c.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,36 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "ray.h"
+#include "vector.h"
+
+ray_t ray_transform(ray_t r, mat4_t xform)
+{
+	mat4_t upper;
+	vec3_t dir;
+
+	m4_copy(upper, xform);
+	upper[0][3] = upper[1][3] = upper[2][3] = upper[3][0] = upper[3][1] = upper[3][2] = 0.0;
+	upper[3][3] = 1.0;
+
+	dir = v3_sub(r.dir, r.origin);
+	dir = v3_transform(dir, upper);
+	r.origin = v3_transform(r.origin, xform);
+	r.dir = v3_add(dir, r.origin);
+	return r;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/sphvec.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/sphvec.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,27 @@
+#include "sphvec.h"
+#include "vector.h"
+
+/* theta: 0 <= theta <= 2pi, the angle around Y axis.
+ * phi: 0 <= phi <= pi, the angle from Y axis.
+ * 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 -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/sphvec.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/sphvec.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,36 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_SPHVEC_H_
+#define VMATH_SPHVEC_H_
+
+#include "vmath_types.h"
+
+#ifdef __cplusplus
+/* 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	/* __cplusplus */
+
+#endif	/* VMATH_SPHVEC_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vector.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vector.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,326 @@
+#include "vector.h"
+#include "vmath.h"
+
+// ---------- Vector2 -----------
+
+Vector2::Vector2(scalar_t x, scalar_t y)
+{
+	this->x = x;
+	this->y = y;
+}
+
+Vector2::Vector2(const vec2_t &vec)
+{
+	x = vec.x;
+	y = vec.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 vec3_t &vec)
+{
+	x = vec.x;
+	y = vec.y;
+	z = vec.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
+{
+	return refraction(normal, src_ior / dst_ior);
+}
+
+Vector3 Vector3::refraction(const Vector3 &normal, scalar_t ior) const
+{
+	scalar_t cos_inc = dot_product(*this, -normal);
+
+	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 vec4_t &vec)
+{
+	x = vec.x;
+	y = vec.y;
+	z = vec.z;
+	w = vec.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 -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vector.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vector.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,292 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_VECTOR_H_
+#define VMATH_VECTOR_H_
+
+#include <stdio.h>
+#include "vmath_types.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+/* C 2D vector functions */
+static inline vec2_t v2_cons(scalar_t x, scalar_t y);
+static inline void v2_print(FILE *fp, vec2_t v);
+
+static inline vec2_t v2_add(vec2_t v1, vec2_t v2);
+static inline vec2_t v2_sub(vec2_t v1, vec2_t v2);
+static inline vec2_t v2_scale(vec2_t v, scalar_t s);
+static inline scalar_t v2_dot(vec2_t v1, vec2_t v2);
+static inline scalar_t v2_length(vec2_t v);
+static inline scalar_t v2_length_sq(vec2_t v);
+static inline vec2_t v2_normalize(vec2_t v);
+
+static inline vec2_t v2_lerp(vec2_t v1, vec2_t v2, scalar_t t);
+
+/* C 3D vector functions */
+static inline vec3_t v3_cons(scalar_t x, scalar_t y, scalar_t z);
+static inline void v3_print(FILE *fp, vec3_t v);
+
+static inline vec3_t v3_add(vec3_t v1, vec3_t v2);
+static inline vec3_t v3_sub(vec3_t v1, vec3_t v2);
+static inline vec3_t v3_neg(vec3_t v);
+static inline vec3_t v3_mul(vec3_t v1, vec3_t v2);
+static inline vec3_t v3_scale(vec3_t v1, scalar_t s);
+static inline scalar_t v3_dot(vec3_t v1, vec3_t v2);
+static inline vec3_t v3_cross(vec3_t v1, vec3_t v2);
+static inline scalar_t v3_length(vec3_t v);
+static inline scalar_t v3_length_sq(vec3_t v);
+static inline vec3_t v3_normalize(vec3_t v);
+static inline vec3_t v3_transform(vec3_t v, mat4_t m);
+
+static inline vec3_t v3_rotate(vec3_t v, scalar_t x, scalar_t y, scalar_t z);
+static inline vec3_t v3_rotate_axis(vec3_t v, scalar_t angle, scalar_t x, scalar_t y, scalar_t z);
+static inline vec3_t v3_rotate_quat(vec3_t v, quat_t q);
+
+static inline vec3_t v3_reflect(vec3_t v, vec3_t n);
+
+static inline vec3_t v3_lerp(vec3_t v1, vec3_t v2, scalar_t t);
+
+/* C 4D vector functions */
+static inline vec4_t v4_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t w);
+static inline void v4_print(FILE *fp, vec4_t v);
+
+static inline vec4_t v4_add(vec4_t v1, vec4_t v2);
+static inline vec4_t v4_sub(vec4_t v1, vec4_t v2);
+static inline vec4_t v4_neg(vec4_t v);
+static inline vec4_t v4_mul(vec4_t v1, vec4_t v2);
+static inline vec4_t v4_scale(vec4_t v, scalar_t s);
+static inline scalar_t v4_dot(vec4_t v1, vec4_t v2);
+static inline scalar_t v4_length(vec4_t v);
+static inline scalar_t v4_length_sq(vec4_t v);
+static inline vec4_t v4_normalize(vec4_t v);
+static inline vec4_t v4_transform(vec4_t v, mat4_t m);
+
+#ifdef __cplusplus
+}	/* extern "C" */
+
+/* when included from C++ source files, also define the vector classes */
+#include <iostream>
+
+/** 2D Vector */
+class Vector2 {
+public:
+	scalar_t x, y;
+
+	Vector2(scalar_t x = 0.0, scalar_t y = 0.0);
+	Vector2(const vec2_t &vec);
+	Vector2(const Vector3 &vec);
+	Vector2(const Vector4 &vec);
+
+	inline scalar_t &operator [](int elem);
+	inline const scalar_t &operator [](int elem) const;
+
+	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;
+};
+
+/* unary operations */
+inline Vector2 operator -(const Vector2 &vec);
+
+/* binary vector (op) vector operations */
+inline scalar_t dot_product(const Vector2 &v1, const Vector2 &v2);
+
+inline Vector2 operator +(const Vector2 &v1, const Vector2 &v2);
+inline Vector2 operator -(const Vector2 &v1, const Vector2 &v2);
+inline Vector2 operator *(const Vector2 &v1, const Vector2 &v2);
+inline Vector2 operator /(const Vector2 &v1, const Vector2 &v2);
+inline bool operator ==(const Vector2 &v1, const Vector2 &v2);
+
+inline void operator +=(Vector2 &v1, const Vector2 &v2);
+inline void operator -=(Vector2 &v1, const Vector2 &v2);
+inline void operator *=(Vector2 &v1, const Vector2 &v2);
+inline void operator /=(Vector2 &v1, const Vector2 &v2);
+
+/* binary vector (op) scalar and scalar (op) vector operations */
+inline Vector2 operator +(const Vector2 &vec, scalar_t scalar);
+inline Vector2 operator +(scalar_t scalar, const Vector2 &vec);
+inline Vector2 operator -(const Vector2 &vec, scalar_t scalar);
+inline Vector2 operator *(const Vector2 &vec, scalar_t scalar);
+inline Vector2 operator *(scalar_t scalar, const Vector2 &vec);
+inline Vector2 operator /(const Vector2 &vec, scalar_t scalar);
+
+inline void operator +=(Vector2 &vec, scalar_t scalar);
+inline void operator -=(Vector2 &vec, scalar_t scalar);
+inline void operator *=(Vector2 &vec, scalar_t scalar);
+inline void operator /=(Vector2 &vec, scalar_t scalar);
+
+std::ostream &operator <<(std::ostream &out, const Vector2 &vec);
+
+inline Vector2 lerp(const Vector2 &a, const Vector2 &b, scalar_t t);
+inline Vector2 catmull_rom_spline(const Vector2 &v0, const Vector2 &v1,
+		const Vector2 &v2, const Vector2 &v3, scalar_t t);
+
+/* 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 vec3_t &vec);
+	Vector3(const Vector2 &vec);
+	Vector3(const Vector4 &vec);
+	Vector3(const SphVector &sph);
+
+	Vector3 &operator =(const SphVector &sph);
+
+	inline scalar_t &operator [](int elem);
+	inline const scalar_t &operator [](int elem) const;
+
+	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;
+	Vector3 refraction(const Vector3 &normal, scalar_t ior) const;
+};
+
+/* unary operations */
+inline Vector3 operator -(const Vector3 &vec);
+
+/* binary vector (op) vector operations */
+inline scalar_t dot_product(const Vector3 &v1, const Vector3 &v2);
+inline Vector3 cross_product(const Vector3 &v1, const Vector3 &v2);
+
+inline Vector3 operator +(const Vector3 &v1, const Vector3 &v2);
+inline Vector3 operator -(const Vector3 &v1, const Vector3 &v2);
+inline Vector3 operator *(const Vector3 &v1, const Vector3 &v2);
+inline Vector3 operator /(const Vector3 &v1, const Vector3 &v2);
+inline bool operator ==(const Vector3 &v1, const Vector3 &v2);
+
+inline void operator +=(Vector3 &v1, const Vector3 &v2);
+inline void operator -=(Vector3 &v1, const Vector3 &v2);
+inline void operator *=(Vector3 &v1, const Vector3 &v2);
+inline void operator /=(Vector3 &v1, const Vector3 &v2);
+
+/* binary vector (op) scalar and scalar (op) vector operations */
+inline Vector3 operator +(const Vector3 &vec, scalar_t scalar);
+inline Vector3 operator +(scalar_t scalar, const Vector3 &vec);
+inline Vector3 operator -(const Vector3 &vec, scalar_t scalar);
+inline Vector3 operator *(const Vector3 &vec, scalar_t scalar);
+inline Vector3 operator *(scalar_t scalar, const Vector3 &vec);
+inline Vector3 operator /(const Vector3 &vec, scalar_t scalar);
+
+inline void operator +=(Vector3 &vec, scalar_t scalar);
+inline void operator -=(Vector3 &vec, scalar_t scalar);
+inline void operator *=(Vector3 &vec, scalar_t scalar);
+inline void operator /=(Vector3 &vec, scalar_t scalar);
+
+std::ostream &operator <<(std::ostream &out, const Vector3 &vec);
+
+inline Vector3 lerp(const Vector3 &a, const Vector3 &b, scalar_t t);
+inline Vector3 catmull_rom_spline(const Vector3 &v0, const Vector3 &v1,
+		const Vector3 &v2, const Vector3 &v3, scalar_t t);
+
+/* 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 vec4_t &vec);
+	Vector4(const Vector2 &vec);
+	Vector4(const Vector3 &vec);
+
+	inline scalar_t &operator [](int elem);
+	inline const scalar_t &operator [](int elem) const;
+
+	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;
+};
+
+
+/* unary operations */
+inline Vector4 operator -(const Vector4 &vec);
+
+/* binary vector (op) vector operations */
+inline scalar_t dot_product(const Vector4 &v1, const Vector4 &v2);
+inline Vector4 cross_product(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3);
+
+inline Vector4 operator +(const Vector4 &v1, const Vector4 &v2);
+inline Vector4 operator -(const Vector4 &v1, const Vector4 &v2);
+inline Vector4 operator *(const Vector4 &v1, const Vector4 &v2);
+inline Vector4 operator /(const Vector4 &v1, const Vector4 &v2);
+inline bool operator ==(const Vector4 &v1, const Vector4 &v2);
+
+inline void operator +=(Vector4 &v1, const Vector4 &v2);
+inline void operator -=(Vector4 &v1, const Vector4 &v2);
+inline void operator *=(Vector4 &v1, const Vector4 &v2);
+inline void operator /=(Vector4 &v1, const Vector4 &v2);
+
+/* binary vector (op) scalar and scalar (op) vector operations */
+inline Vector4 operator +(const Vector4 &vec, scalar_t scalar);
+inline Vector4 operator +(scalar_t scalar, const Vector4 &vec);
+inline Vector4 operator -(const Vector4 &vec, scalar_t scalar);
+inline Vector4 operator *(const Vector4 &vec, scalar_t scalar);
+inline Vector4 operator *(scalar_t scalar, const Vector4 &vec);
+inline Vector4 operator /(const Vector4 &vec, scalar_t scalar);
+
+inline void operator +=(Vector4 &vec, scalar_t scalar);
+inline void operator -=(Vector4 &vec, scalar_t scalar);
+inline void operator *=(Vector4 &vec, scalar_t scalar);
+inline void operator /=(Vector4 &vec, scalar_t scalar);
+
+std::ostream &operator <<(std::ostream &out, const Vector4 &vec);
+
+inline Vector4 lerp(const Vector4 &v0, const Vector4 &v1, scalar_t t);
+inline Vector4 catmull_rom_spline(const Vector4 &v0, const Vector4 &v1,
+		const Vector4 &v2, const Vector4 &v3, scalar_t t);
+
+#endif	/* __cplusplus */
+
+#include "vector.inl"
+
+#endif	/* VMATH_VECTOR_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vector.inl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vector.inl	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,761 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <math.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+/* C 2D vector functions */
+static inline vec2_t v2_cons(scalar_t x, scalar_t y)
+{
+	vec2_t v;
+	v.x = x;
+	v.y = y;
+	return v;
+}
+
+static inline void v2_print(FILE *fp, vec2_t v)
+{
+	fprintf(fp, "[ %.4f %.4f ]", v.x, v.y);
+}
+
+static inline vec2_t v2_add(vec2_t v1, vec2_t v2)
+{
+	vec2_t res;
+	res.x = v1.x + v2.x;
+	res.y = v1.y + v2.y;
+	return res;
+}
+
+static inline vec2_t v2_sub(vec2_t v1, vec2_t v2)
+{
+	vec2_t res;
+	res.x = v1.x - v2.x;
+	res.y = v1.y - v2.y;
+	return res;
+}
+
+static inline vec2_t v2_scale(vec2_t v, scalar_t s)
+{
+	vec2_t res;
+	res.x = v.x * s;
+	res.y = v.y * s;
+	return res;
+}
+
+static inline scalar_t v2_dot(vec2_t v1, vec2_t v2)
+{
+	return v1.x * v2.x + v1.y * v2.y;
+}
+
+static inline scalar_t v2_length(vec2_t v)
+{
+	return sqrt(v.x * v.x + v.y * v.y);
+}
+
+static inline scalar_t v2_length_sq(vec2_t v)
+{
+	return v.x * v.x + v.y * v.y;
+}
+
+static inline vec2_t v2_normalize(vec2_t v)
+{
+	scalar_t len = (scalar_t)sqrt(v.x * v.x + v.y * v.y);
+	v.x /= len;
+	v.y /= len;
+	return v;
+}
+
+static inline vec2_t v2_lerp(vec2_t v1, vec2_t v2, scalar_t t)
+{
+	vec2_t res;
+	res.x = v1.x + (v2.x - v1.x) * t;
+	res.y = v1.y + (v2.y - v1.y) * t;
+	return res;
+}
+
+
+/* C 3D vector functions */
+static inline vec3_t v3_cons(scalar_t x, scalar_t y, scalar_t z)
+{
+	vec3_t v;
+	v.x = x;
+	v.y = y;
+	v.z = z;
+	return v;
+}
+
+static inline void v3_print(FILE *fp, vec3_t v)
+{
+	fprintf(fp, "[ %.4f %.4f %.4f ]", v.x, v.y, v.z);
+}
+
+static inline vec3_t v3_add(vec3_t v1, vec3_t v2)
+{
+	v1.x += v2.x;
+	v1.y += v2.y;
+	v1.z += v2.z;
+	return v1;
+}
+
+static inline vec3_t v3_sub(vec3_t v1, vec3_t v2)
+{
+	v1.x -= v2.x;
+	v1.y -= v2.y;
+	v1.z -= v2.z;
+	return v1;
+}
+
+static inline vec3_t v3_neg(vec3_t v)
+{
+	v.x = -v.x;
+	v.y = -v.y;
+	v.z = -v.z;
+	return v;
+}
+
+static inline vec3_t v3_mul(vec3_t v1, vec3_t v2)
+{
+	v1.x *= v2.x;
+	v1.y *= v2.y;
+	v1.z *= v2.z;
+	return v1;
+}
+
+static inline vec3_t v3_scale(vec3_t v1, scalar_t s)
+{
+	v1.x *= s;
+	v1.y *= s;
+	v1.z *= s;
+	return v1;
+}
+
+static inline scalar_t v3_dot(vec3_t v1, vec3_t v2)
+{
+	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
+}
+
+static inline vec3_t v3_cross(vec3_t v1, vec3_t v2)
+{
+	vec3_t v;
+	v.x = v1.y * v2.z - v1.z * v2.y;
+	v.y = v1.z * v2.x - v1.x * v2.z;
+	v.z = v1.x * v2.y - v1.y * v2.x;
+	return v;
+}
+
+static inline scalar_t v3_length(vec3_t v)
+{
+	return sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
+}
+
+static inline scalar_t v3_length_sq(vec3_t v)
+{
+	return v.x * v.x + v.y * v.y + v.z * v.z;
+}
+
+static inline vec3_t v3_normalize(vec3_t v)
+{
+	scalar_t len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
+	v.x /= len;
+	v.y /= len;
+	v.z /= len;
+	return v;
+}
+
+static inline vec3_t v3_transform(vec3_t v, mat4_t m)
+{
+	vec3_t res;
+	res.x = m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z + m[0][3];
+	res.y = m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z + m[1][3];
+	res.z = m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z + m[2][3];
+	return res;
+}
+
+static inline vec3_t v3_rotate(vec3_t v, scalar_t x, scalar_t y, scalar_t z)
+{
+	void m4_rotate(mat4_t, scalar_t, scalar_t, scalar_t);
+
+	mat4_t m = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
+	m4_rotate(m, x, y, z);
+	return v3_transform(v, m);
+}
+
+static inline vec3_t v3_rotate_axis(vec3_t v, scalar_t angle, scalar_t x, scalar_t y, scalar_t z)
+{
+	void m4_rotate_axis(mat4_t, scalar_t, scalar_t, scalar_t, scalar_t);
+
+	mat4_t m = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
+	m4_rotate_axis(m, angle, x, y, z);
+	return v3_transform(v, m);
+}
+
+static inline vec3_t v3_rotate_quat(vec3_t v, quat_t q)
+{
+	quat_t quat_rotate_quat(quat_t, quat_t);
+
+	quat_t vq = v4_cons(v.x, v.y, v.z, 0.0);
+	quat_t res = quat_rotate_quat(vq, q);
+	return v3_cons(res.x, res.y, res.z);
+}
+
+static inline vec3_t v3_reflect(vec3_t v, vec3_t n)
+{
+	scalar_t dot = v3_dot(v, n);
+	return v3_sub(v3_scale(n, dot * 2.0), v);
+}
+
+static inline vec3_t v3_lerp(vec3_t v1, vec3_t v2, scalar_t t)
+{
+	v1.x += (v2.x - v1.x) * t;
+	v1.y += (v2.y - v1.y) * t;
+	v1.z += (v2.z - v1.z) * t;
+	return v1;
+}
+
+/* C 4D vector functions */
+static inline vec4_t v4_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t w)
+{
+	vec4_t v;
+	v.x = x;
+	v.y = y;
+	v.z = z;
+	v.w = w;
+	return v;
+}
+
+static inline void v4_print(FILE *fp, vec4_t v)
+{
+	fprintf(fp, "[ %.4f %.4f %.4f %.4f ]", v.x, v.y, v.z, v.w);
+}
+
+static inline vec4_t v4_add(vec4_t v1, vec4_t v2)
+{
+	v1.x += v2.x;
+	v1.y += v2.y;
+	v1.z += v2.z;
+	v1.w += v2.w;
+	return v1;
+}
+
+static inline vec4_t v4_sub(vec4_t v1, vec4_t v2)
+{
+	v1.x -= v2.x;
+	v1.y -= v2.y;
+	v1.z -= v2.z;
+	v1.w -= v2.w;
+	return v1;
+}
+
+static inline vec4_t v4_neg(vec4_t v)
+{
+	v.x = -v.x;
+	v.y = -v.y;
+	v.z = -v.z;
+	v.w = -v.w;
+	return v;
+}
+
+static inline vec4_t v4_mul(vec4_t v1, vec4_t v2)
+{
+	v1.x *= v2.x;
+	v1.y *= v2.y;
+	v1.z *= v2.z;
+	v1.w *= v2.w;
+	return v1;
+}
+
+static inline vec4_t v4_scale(vec4_t v, scalar_t s)
+{
+	v.x *= s;
+	v.y *= s;
+	v.z *= s;
+	v.w *= s;
+	return v;
+}
+
+static inline scalar_t v4_dot(vec4_t v1, vec4_t v2)
+{
+	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z + v1.w * v2.w;
+}
+
+static inline scalar_t v4_length(vec4_t v)
+{
+	return sqrt(v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w);
+}
+
+static inline scalar_t v4_length_sq(vec4_t v)
+{
+	return v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w;
+}
+
+static inline vec4_t v4_normalize(vec4_t v)
+{
+	scalar_t len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w);
+	v.x /= len;
+	v.y /= len;
+	v.z /= len;
+	v.w /= len;
+	return v;
+}
+
+static inline vec4_t v4_transform(vec4_t v, mat4_t m)
+{
+	vec4_t res;
+	res.x = m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z + m[0][3] * v.w;
+	res.y = m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z + m[1][3] * v.w;
+	res.z = m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z + m[2][3] * v.w;
+	res.w = m[3][0] * v.x + m[3][1] * v.y + m[3][2] * v.z + m[3][3] * v.w;
+	return res;
+}
+
+#ifdef __cplusplus
+}	/* extern "C" */
+
+
+/* --------------- C++ part -------------- */
+
+inline scalar_t &Vector2::operator [](int elem) {
+	return elem ? y : x;
+}
+
+inline const scalar_t &Vector2::operator [](int elem) const {
+	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 -(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 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;
+}
+
+inline Vector2 lerp(const Vector2 &a, const Vector2 &b, scalar_t t)
+{
+	return a + (b - a) * t;
+}
+
+inline Vector2 catmull_rom_spline(const Vector2 &v0, const Vector2 &v1,
+		const Vector2 &v2, const Vector2 &v3, scalar_t t)
+{
+	scalar_t spline(scalar_t, scalar_t, scalar_t, scalar_t, scalar_t);
+	scalar_t x = spline(v0.x, v1.x, v2.x, v3.x, t);
+	scalar_t y = spline(v0.y, v1.y, v2.y, v3.y, t);
+	return Vector2(x, y);
+}
+
+
+/* ------------- Vector3 -------------- */
+
+inline scalar_t &Vector3::operator [](int elem) {
+	return elem ? (elem == 1 ? y : z) : x;
+}
+
+inline const scalar_t &Vector3::operator [](int elem) const {
+	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 *(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 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;
+}
+
+inline Vector3 lerp(const Vector3 &a, const Vector3 &b, scalar_t t) {
+	return a + (b - a) * t;
+}
+
+inline Vector3 catmull_rom_spline(const Vector3 &v0, const Vector3 &v1,
+		const Vector3 &v2, const Vector3 &v3, scalar_t t)
+{
+	scalar_t spline(scalar_t, scalar_t, scalar_t, scalar_t, scalar_t);
+	scalar_t x = spline(v0.x, v1.x, v2.x, v3.x, t);
+	scalar_t y = spline(v0.y, v1.y, v2.y, v3.y, t);
+	scalar_t z = spline(v0.z, v1.z, v2.z, v3.z, t);
+	return Vector3(x, y, z);
+}
+
+/* ----------- Vector4 ----------------- */
+
+inline scalar_t &Vector4::operator [](int elem) {
+	return elem ? (elem == 1 ? y : (elem == 2 ? z : w)) : x;
+}
+
+inline const scalar_t &Vector4::operator [](int elem) const {
+	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) {
+	scalar_t 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 *(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 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;
+}
+
+inline Vector4 lerp(const Vector4 &v0, const Vector4 &v1, scalar_t t)
+{
+	return v0 + (v1 - v0) * t;
+}
+
+inline Vector4 catmull_rom_spline(const Vector4 &v0, const Vector4 &v1,
+		const Vector4 &v2, const Vector4 &v3, scalar_t t)
+{
+	scalar_t spline(scalar_t, scalar_t, scalar_t, scalar_t, scalar_t);
+	scalar_t x = spline(v0.x, v1.x, v2.x, v3.x, t);
+	scalar_t y = spline(v0.y, v1.y, v2.y, v3.y, t);
+	scalar_t z = spline(v0.z, v1.z, v2.z, v3.z, t);
+	scalar_t w = spline(v0.w, v1.w, v2.w, v3.w, t);
+	return Vector4(x, y, z, w);
+}
+
+#endif	/* __cplusplus */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vmath.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vmath.c	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,336 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdlib.h>
+#include <math.h>
+#include "vmath.h"
+
+/** Numerical calculation of integrals using simpson's rule */
+scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples)
+{
+	int i;
+	scalar_t h = (high - low) / (scalar_t)samples;
+	scalar_t sum = 0.0;
+
+	for(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(M_E, exponent) / (sdev * sqrt(TWO_PI));
+}
+
+
+/** b-spline approximation */
+scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
+{
+	vec4_t tmp;
+	scalar_t tsq = t * t;
+
+	static mat4_t bspline_mat = {
+		{-1,  3, -3,  1},
+		{3, -6,  3,  0},
+		{-3,  0,  3,  0},
+		{1,  4,  1,  0}
+	};
+
+	tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), bspline_mat), 1.0 / 6.0);
+	return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);
+}
+
+/** Catmull-rom spline interpolation */
+scalar_t spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
+{
+	vec4_t tmp;
+	scalar_t tsq = t * t;
+
+	static mat4_t crspline_mat = {
+		{-1,  3, -3,  1},
+		{2, -5,  4, -1},
+		{-1,  0,  1,  0},
+		{0,  2,  0,  0}
+	};
+
+	tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), crspline_mat), 0.5);
+	return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);
+}
+
+/** Bezier interpolation */
+scalar_t bezier(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
+{
+	scalar_t omt, omt3, t3, f;
+	t3 = t * t * t;
+	omt = 1.0f - t;
+	omt3 = omt * omt * omt;
+	f = 3 * t * omt;
+
+	return (a * omt3) + (b * f * omt) + (c * f * t) + (d * t3);
+}
+
+/* ---- Ken Perlin's implementation of noise ---- */
+
+#define B	0x100
+#define BM	0xff
+#define N	0x1000
+#define NP	12   /* 2^N */
+#define NM	0xfff
+
+#define s_curve(t) (t * t * (3.0f - 2.0f * t))
+
+#define setup(elem, b0, b1, r0, r1) \
+	do {							\
+		scalar_t t = elem + N;		\
+		b0 = ((int)t) & BM;			\
+		b1 = (b0 + 1) & BM;			\
+		r0 = t - (int)t;			\
+		r1 = r0 - 1.0f;				\
+	} while(0)
+
+
+static int perm[B + B + 2];			/* permuted index from g_n onto themselves */
+static vec3_t grad3[B + B + 2];		/* 3D random gradients */
+static vec2_t grad2[B + B + 2];		/* 2D random gradients */
+static scalar_t grad1[B + B + 2];	/* 1D random ... slopes */
+static int tables_valid;
+
+static void init_noise()
+{
+	int i;
+
+	/* calculate random gradients */
+	for(i=0; i<B; i++) {
+		perm[i] = i;	/* .. and initialize permutation mapping to identity */
+
+		grad1[i] = (scalar_t)((rand() % (B + B)) - B) / B;
+
+		grad2[i].x = (scalar_t)((rand() % (B + B)) - B) / B;
+		grad2[i].y = (scalar_t)((rand() % (B + B)) - B) / B;
+		grad2[i] = v2_normalize(grad2[i]);
+
+		grad3[i].x = (scalar_t)((rand() % (B + B)) - B) / B;
+		grad3[i].y = (scalar_t)((rand() % (B + B)) - B) / B;
+		grad3[i].z = (scalar_t)((rand() % (B + B)) - B) / B;
+		grad3[i] = v3_normalize(grad3[i]);
+	}
+
+	/* permute indices by swapping them randomly */
+	for(i=0; i<B; i++) {
+		int rand_idx = rand() % B;
+
+		int tmp = perm[i];
+		perm[i] = perm[rand_idx];
+		perm[rand_idx] = tmp;
+	}
+
+	/* fill up the rest of the arrays by duplicating the existing gradients */
+	/* and permutations */
+	for(i=0; i<B+2; i++) {
+		perm[B + i] = perm[i];
+		grad1[B + i] = grad1[i];
+		grad2[B + i] = grad2[i];
+		grad3[B + i] = grad3[i];
+	}
+}
+
+scalar_t noise1(scalar_t x)
+{
+	int bx0, bx1;
+	scalar_t rx0, rx1, sx, u, v;
+
+	if(!tables_valid) {
+		init_noise();
+		tables_valid = 1;
+	}
+
+	setup(x, bx0, bx1, rx0, rx1);
+	sx = s_curve(rx0);
+	u = rx0 * grad1[perm[bx0]];
+	v = rx1 * grad1[perm[bx1]];
+
+	return lerp(u, v, sx);
+}
+
+scalar_t noise2(scalar_t x, scalar_t y)
+{
+	int i, j, b00, b10, b01, b11;
+	int bx0, bx1, by0, by1;
+	scalar_t rx0, rx1, ry0, ry1;
+	scalar_t sx, sy, u, v, a, b;
+
+	if(!tables_valid) {
+		init_noise();
+		tables_valid = 1;
+	}
+
+	setup(x, bx0, bx1, rx0, rx1);
+	setup(y, by0, by1, ry0, ry1);
+
+	i = perm[bx0];
+	j = perm[bx1];
+
+	b00 = perm[i + by0];
+	b10 = perm[j + by0];
+	b01 = perm[i + by1];
+	b11 = perm[j + by1];
+
+	/* calculate hermite inteprolating factors */
+	sx = s_curve(rx0);
+	sy = s_curve(ry0);
+
+	/* interpolate along the left edge */
+	u = v2_dot(grad2[b00], v2_cons(rx0, ry0));
+	v = v2_dot(grad2[b10], v2_cons(rx1, ry0));
+	a = lerp(u, v, sx);
+
+	/* interpolate along the right edge */
+	u = v2_dot(grad2[b01], v2_cons(rx0, ry1));
+	v = v2_dot(grad2[b11], v2_cons(rx1, ry1));
+	b = lerp(u, v, sx);
+
+	/* interpolate between them */
+	return lerp(a, b, sy);
+}
+
+scalar_t noise3(scalar_t x, scalar_t y, scalar_t z)
+{
+	int i, j;
+	int bx0, bx1, by0, by1, bz0, bz1;
+	int b00, b10, b01, b11;
+	scalar_t rx0, rx1, ry0, ry1, rz0, rz1;
+	scalar_t sx, sy, sz;
+	scalar_t u, v, a, b, c, d;
+
+	if(!tables_valid) {
+		init_noise();
+		tables_valid = 1;
+	}
+
+	setup(x, bx0, bx1, rx0, rx1);
+	setup(y, by0, by1, ry0, ry1);
+	setup(z, bz0, bz1, rz0, rz1);
+
+	i = perm[bx0];
+	j = perm[bx1];
+
+	b00 = perm[i + by0];
+	b10 = perm[j + by0];
+	b01 = perm[i + by1];
+	b11 = perm[j + by1];
+
+	/* calculate hermite interpolating factors */
+	sx = s_curve(rx0);
+	sy = s_curve(ry0);
+	sz = s_curve(rz0);
+
+	/* interpolate along the top slice of the cell */
+	u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz0));
+	v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz0));
+	a = lerp(u, v, sx);
+
+	u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz0));
+	v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz0));
+	b = lerp(u, v, sx);
+
+	c = lerp(a, b, sy);
+
+	/* interpolate along the bottom slice of the cell */
+	u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz1));
+	v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz1));
+	a = lerp(u, v, sx);
+
+	u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz1));
+	v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz1));
+	b = lerp(u, v, sx);
+
+	d = lerp(a, b, sy);
+
+	/* interpolate between slices */
+	return lerp(c, d, sz);
+}
+
+scalar_t fbm1(scalar_t x, int octaves)
+{
+	int i;
+	scalar_t res = 0.0f, freq = 1.0f;
+	for(i=0; i<octaves; i++) {
+		res += noise1(x * freq) / freq;
+		freq *= 2.0f;
+	}
+	return res;
+}
+
+scalar_t fbm2(scalar_t x, scalar_t y, int octaves)
+{
+	int i;
+	scalar_t res = 0.0f, freq = 1.0f;
+	for(i=0; i<octaves; i++) {
+		res += noise2(x * freq, y * freq) / freq;
+		freq *= 2.0f;
+	}
+	return res;
+}
+
+scalar_t fbm3(scalar_t x, scalar_t y, scalar_t z, int octaves)
+{
+	int i;
+	scalar_t res = 0.0f, freq = 1.0f;
+	for(i=0; i<octaves; i++) {
+		res += noise3(x * freq, y * freq, z * freq) / freq;
+		freq *= 2.0f;
+	}
+	return res;
+}
+
+scalar_t turbulence1(scalar_t x, int octaves)
+{
+	int i;
+	scalar_t res = 0.0f, freq = 1.0f;
+	for(i=0; i<octaves; i++) {
+		res += fabs(noise1(x * freq) / freq);
+		freq *= 2.0f;
+	}
+	return res;
+}
+
+scalar_t turbulence2(scalar_t x, scalar_t y, int octaves)
+{
+	int i;
+	scalar_t res = 0.0f, freq = 1.0f;
+	for(i=0; i<octaves; i++) {
+		res += fabs(noise2(x * freq, y * freq) / freq);
+		freq *= 2.0f;
+	}
+	return res;
+}
+
+scalar_t turbulence3(scalar_t x, scalar_t y, scalar_t z, int octaves)
+{
+	int i;
+	scalar_t res = 0.0f, freq = 1.0f;
+	for(i=0; i<octaves; i++) {
+		res += fabs(noise3(x * freq, y * freq, z * freq) / freq);
+		freq *= 2.0f;
+	}
+	return res;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vmath.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vmath.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,94 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2013 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_H_
+#define VMATH_H_
+
+#include <math.h>
+#include "vmath_types.h"
+
+#ifndef M_PI
+#define M_PI	PI
+#endif
+
+#ifndef M_E
+#define M_E				2.718281828459045
+#endif
+
+#define PI				3.141592653589793
+#define HALF_PI			1.570796326794897
+#define QUARTER_PI		0.785398163397448
+#define TWO_PI			6.283185307179586
+
+
+#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))
+
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+#define MAX(a, b)	((a) > (b) ? (a) : (b))
+
+#ifndef __GNUC__
+#define round(x)	((x) >= 0 ? (x) + 0.5 : (x) - 0.5)
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+static inline scalar_t smoothstep(float a, float b, float x);
+
+static inline scalar_t frand(scalar_t range);
+static inline vec3_t sphrand(scalar_t rad);
+
+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);
+
+static inline scalar_t lerp(scalar_t a, scalar_t b, scalar_t t);
+
+scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
+scalar_t spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
+scalar_t bezier(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
+
+scalar_t noise1(scalar_t x);
+scalar_t noise2(scalar_t x, scalar_t y);
+scalar_t noise3(scalar_t x, scalar_t y, scalar_t z);
+
+scalar_t fbm1(scalar_t x, int octaves);
+scalar_t fbm2(scalar_t x, scalar_t y, int octaves);
+scalar_t fbm3(scalar_t x, scalar_t y, scalar_t z, int octaves);
+
+scalar_t turbulence1(scalar_t x, int octaves);
+scalar_t turbulence2(scalar_t x, scalar_t y, int octaves);
+scalar_t turbulence3(scalar_t x, scalar_t y, scalar_t z, int octaves);
+
+#ifdef __cplusplus
+}
+#endif	/* __cplusplus */
+
+#include "vmath.inl"
+
+#include "vector.h"
+#include "matrix.h"
+#include "quat.h"
+#include "sphvec.h"
+#include "ray.h"
+#include "geom.h"
+
+#endif	/* VMATH_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vmath.inl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vmath.inl	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,56 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdlib.h>
+
+static inline scalar_t smoothstep(float a, float b, float x)
+{
+	if(x < a) return 0.0;
+	if(x >= b) return 1.0;
+
+	x = (x - a) / (b - a);
+	return x * x * (3.0 - 2.0 * x);
+}
+
+/** Generates a random number in [0, range) */
+static inline scalar_t frand(scalar_t range)
+{
+	return range * (scalar_t)rand() / (scalar_t)RAND_MAX;
+}
+
+/** Generates a random vector on the surface of a sphere */
+static inline vec3_t sphrand(scalar_t rad)
+{
+	scalar_t u = (scalar_t)rand() / RAND_MAX;
+	scalar_t v = (scalar_t)rand() / RAND_MAX;
+
+	scalar_t theta = 2.0 * M_PI * u;
+	scalar_t phi = acos(2.0 * v - 1.0);
+
+	vec3_t res;
+	res.x = rad * cos(theta) * sin(phi);
+	res.y = rad * sin(theta) * sin(phi);
+	res.z = rad * cos(phi);
+	return res;
+}
+
+/** linear interpolation */
+static inline scalar_t lerp(scalar_t a, scalar_t b, scalar_t t)
+{
+	return a + (b - a) * t;
+}
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vmath_config.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vmath_config.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,19 @@
+#ifndef VMATH_CONFIG_H_
+#define VMATH_CONFIG_H_
+
+#if (__STDC_VERSION__ < 199999)
+#if defined(__GNUC__) || defined(_MSC_VER)
+#define inline __inline
+#else
+#define inline
+
+#ifdef VECTOR_H_
+#warning "compiling vector operations without inline, performance might suffer"
+#endif	/* VECTOR_H_ */
+
+#endif	/* gcc/msvc */
+#endif	/* not C99 */
+
+#define SINGLE_PRECISION_MATH
+
+#endif	/* VMATH_CONFIG_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 libs/vmath/vmath_types.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libs/vmath/vmath_types.h	Sun Sep 29 08:20:19 2013 +0300
@@ -0,0 +1,58 @@
+/*
+libvmath - a vector math library
+Copyright (C) 2004-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published
+by the Free Software Foundation, either version 3 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VMATH_TYPES_H_
+#define VMATH_TYPES_H_
+
+#include "vmath_config.h"
+
+#define SMALL_NUMBER	1.e-4
+#define XSMALL_NUMBER	1.e-8
+#define ERROR_MARGIN	1.e-6
+
+
+#ifdef SINGLE_PRECISION_MATH
+typedef float scalar_t;
+#else
+typedef double scalar_t;
+#endif	/* floating point precision */
+
+/* vectors */
+typedef struct { scalar_t x, y; } vec2_t;
+typedef struct { scalar_t x, y, z; } vec3_t;
+typedef struct { scalar_t x, y, z, w; } vec4_t;
+
+/* quaternions */
+typedef vec4_t quat_t;
+
+/* matrices */
+typedef scalar_t mat3_t[3][3];
+typedef scalar_t mat4_t[4][4];
+
+
+#ifdef __cplusplus
+class Vector2;
+class Vector3;
+class Vector4;
+class Quaternion;
+class Matrix3x3;
+class Matrix4x4;
+class SphVector;
+#endif	/* __cplusplus */
+
+#endif	/* VMATH_TYPES_H_ */
diff -r 1c601bf07b86 -r 4deb0b12fe14 src/goat3d.cc
--- a/src/goat3d.cc	Sat Sep 28 20:36:55 2013 +0300
+++ b/src/goat3d.cc	Sun Sep 29 08:20:19 2013 +0300
@@ -25,6 +25,7 @@
 struct goat3d *goat3d_create(void)
 {
 	goat3d *goat = new goat3d;
+	goat->flags = 0;
 	goat->scn = new Scene;
 	return goat;
 }
@@ -142,6 +143,11 @@
 }
 
 // ---- materials ----
+void goat3d_add_mtl(struct goat3d *g, struct goat3d_material *mtl)
+{
+	g->scn->add_material(mtl);
+}
+
 struct goat3d_material *goat3d_create_mtl(void)
 {
 	return new goat3d_material;
@@ -197,12 +203,27 @@
 	return (*mtl)[attrib].map.c_str();
 }
 
-void goat3d_add_mtl(struct goat3d *g, struct goat3d_material *mtl)
+// ---- meshes ----
+void goat3d_add_mesh(struct goat3d *g, struct goat3d_mesh *mesh)
 {
-	g->scn->add_material(mtl);
+	g->scn->add_mesh(mesh);
 }
 
-// ---- meshes ----
+int goat3d_get_mesh_count(struct goat3d *g)
+{
+	return g->scn->get_mesh_count();
+}
+
+struct goat3d_mesh *goat3d_get_mesh(struct goat3d *g, int idx)
+{
+	return (goat3d_mesh*)g->scn->get_mesh(idx);
+}
+
+struct goat3d_mesh *goat3d_get_mesh_by_name(struct goat3d *g, const char *name)
+{
+	return (goat3d_mesh*)g->scn->get_mesh(name);
+}
+
 struct goat3d_mesh *goat3d_create_mesh(void)
 {
 	return new goat3d_mesh;
@@ -459,23 +480,92 @@
 	im_use[GOAT3D_MESH_ATTR_COLOR] = true;
 }
 
-void goat3d_add_mesh(struct goat3d *g, struct goat3d_mesh *mesh)
+/* lights */
+void goat3d_add_light(struct goat3d *g, struct goat3d_light *lt)
 {
-	g->scn->add_mesh(mesh);
+	g->scn->add_light(lt);
 }
 
-int goat3d_get_mesh_count(struct goat3d *g)
+int goat3d_get_light_count(struct goat3d *g)
 {
-	return g->scn->get_mesh_count();
+	return g->scn->get_light_count();
 }
 
-struct goat3d_mesh *goat3d_get_mesh(struct goat3d *g, int idx)
+struct goat3d_light *goat3d_get_light(struct goat3d *g, int idx)
 {
-	return (goat3d_mesh*)g->scn->get_mesh(idx);
+	return (goat3d_light*)g->scn->get_light(idx);
 }
 
+struct goat3d_light *goat3d_get_light_by_name(struct goat3d *g, const char *name)
+{
+	return (goat3d_light*)g->scn->get_light(name);
+}
+
+
+struct goat3d_light *goat3d_create_light(void)
+{
+	return new goat3d_light;
+}
+
+void goat3d_destroy_light(struct goat3d_light *lt)
+{
+	delete lt;
+}
+
+
+/* cameras */
+void goat3d_add_camera(struct goat3d *g, struct goat3d_camera *cam)
+{
+	g->scn->add_camera(cam);
+}
+
+int goat3d_get_camera_count(struct goat3d *g)
+{
+	return g->scn->get_camera_count();
+}
+
+struct goat3d_camera *goat3d_get_camera(struct goat3d *g, int idx)
+{
+	return (goat3d_camera*)g->scn->get_camera(idx);
+}
+
+struct goat3d_camera *goat3d_get_camera_by_name(struct goat3d *g, const char *name)
+{
+	return (goat3d_camera*)g->scn->get_camera(name);
+}
+
+struct goat3d_camera *goat3d_create_camera(void)
+{
+	return new goat3d_camera;
+}
+
+void goat3d_destroy_camera(struct goat3d_camera *cam)
+{
+	delete cam;
+}
+
+
 
 // node
+void goat3d_add_node(struct goat3d *g, struct goat3d_node *node)
+{
+	g->scn->add_node(node);
+}
+
+int goat3d_get_node_count(struct goat3d *g)
+{
+	return g->scn->get_node_count();
+}
+
+struct goat3d_node *goat3d_get_node(struct goat3d *g, int idx)
+{
+	return (goat3d_node*)g->scn->get_node(idx);
+}
+
+struct goat3d_node *goat3d_get_node_by_name(struct goat3d *g, const char *name)
+{
+	return (goat3d_node*)g->scn->get_node(name);
+}
 
 struct goat3d_node *goat3d_create_node(void)
 {
diff -r 1c601bf07b86 -r 4deb0b12fe14 src/goat3d.h
--- a/src/goat3d.h	Sat Sep 28 20:36:55 2013 +0300
+++ b/src/goat3d.h	Sun Sep 29 08:20:19 2013 +0300
@@ -90,6 +90,8 @@
 const float *goat3d_get_ambient(const struct goat3d *g);
 
 /* materials */
+void goat3d_add_mtl(struct goat3d *g, struct goat3d_material *mtl);
+
 struct goat3d_material *goat3d_create_mtl(void);
 void goat3d_destroy_mtl(struct goat3d_material *mtl);
 
@@ -105,9 +107,12 @@
 void goat3d_set_mtl_attrib_map(struct goat3d_material *mtl, const char *attrib, const char *mapname);
 const char *goat3d_get_mtl_attrib_map(struct goat3d_material *mtl, const char *attrib);
 
-void goat3d_add_mtl(struct goat3d *g, struct goat3d_material *mtl);
+/* meshes */
+void goat3d_add_mesh(struct goat3d *g, struct goat3d_mesh *mesh);
+int goat3d_get_mesh_count(struct goat3d *g);
+struct goat3d_mesh *goat3d_get_mesh(struct goat3d *g, int idx);
+struct goat3d_mesh *goat3d_get_mesh_by_name(struct goat3d *g, const char *name);
 
-/* meshes */
 struct goat3d_mesh *goat3d_create_mesh(void);
 void goat3d_destroy_mesh(struct goat3d_mesh *mesh);
 
@@ -160,13 +165,32 @@
 void goat3d_color3f(float x, float y, float z);
 void goat3d_color4f(float x, float y, float z, float w);
 
-void goat3d_add_mesh(struct goat3d *g, struct goat3d_mesh *mesh);
+/* lights */
+void goat3d_add_light(struct goat3d *g, struct goat3d_light *lt);
+int goat3d_get_light_count(struct goat3d *g);
+struct goat3d_light *goat3d_get_light(struct goat3d *g, int idx);
+struct goat3d_light *goat3d_get_light_by_name(struct goat3d *g, const char *name);
 
-int goat3d_get_mesh_count(struct goat3d *g);
-struct goat3d_mesh *goat3d_get_mesh(struct goat3d *g, int idx);
+struct goat3d_light *goat3d_create_light(void);
+void goat3d_destroy_light(struct goat3d_light *lt);
+
+/* cameras */
+void goat3d_add_camera(struct goat3d *g, struct goat3d_camera *cam);
+int goat3d_get_camera_count(struct goat3d *g);
+struct goat3d_camera *goat3d_get_camera(struct goat3d *g, int idx);
+struct goat3d_camera *goat3d_get_camera_by_name(struct goat3d *g, const char *name);
+
+struct goat3d_camera *goat3d_create_camera(void);
+void goat3d_destroy_camera(struct goat3d_camera *cam);
 
 /* nodes */
+void goat3d_add_node(struct goat3d *g, struct goat3d_node *node);
+int goat3d_get_node_count(struct goat3d *g);
+struct goat3d_node *goat3d_get_node(struct goat3d *g, int idx);
+struct goat3d_node *goat3d_get_node_by_name(struct goat3d *g, const char *name);
+
 struct goat3d_node *goat3d_create_node(void);
+void goat3d_destroy_node(struct goat3d_node *node);
 
 void goat3d_set_node_name(struct goat3d_node *node, const char *name);
 const char *goat3d_get_node_name(const struct goat3d_node *node);