Particles

=========

Merge of the famous particle patch by Janne Karhu, a full rewrite
of the Blender particle system. This includes:

- Emitter, Hair and Reactor particle types.
- Newtonian, Keyed and Boids physics.
- Various particle visualisation and rendering types.
- Vertex group and texture control for various properties.
- Interpolated child particles from parents.
- Hair editing with combing, growing, cutting, .. .
- Explode modifier.
- Harmonic, Magnetic fields, and multiple falloff types.

.. and lots of other things, some more info is here:

http://wiki.blender.org/index.php/BlenderDev/Particles_Rewrite
http://wiki.blender.org/index.php/BlenderDev/Particles_Rewrite_Doc

The new particle system cannot be backwards compatible. Old particle
systems are being converted to the new system, but will require
tweaking to get them looking the same as before.

Point Cache
===========

The new system to replace manual baking, based on automatic caching
on disk. This is currently used by softbodies and the particle system.

See the Cache API section on:
http://wiki.blender.org/index.php/BlenderDev/PhysicsSprint

Documentation
=============

These new features still need good docs for the release logs, help
for this is appreciated.

4 files changed, 958 insertions, 6 deletions

diff --git a/source/blender/blenlib/BLI_arithb.h b/source/blender/blenlib/BLI_arithb.h
index 32308f89a15..11a1fd277da 100644
--- a/source/blender/blenlib/BLI_arithb.h
+++ b/source/blender/blenlib/BLI_arithb.h
@@ -248,6 +248,7 @@ void VecMidf(float *v, float *v1, float *v2);
 void VecOrthoBasisf(float *v, float *v1, float *v2);
 
 float Vec2Lenf(float *v1, float *v2);
+float Vec2Length(float *v);
 void Vec2Mulf(float *v1, float f);
 void Vec2Addf(float *v, float *v1, float *v2);
 void Vec2Subf(float *v, float *v1, float *v2);
@@ -281,6 +282,11 @@ float AreaPoly3Dfl(int nr, float *verts, float *normal);
 extern short IsectLL2Df(float *v1, float *v2, float *v3, float *v4);
 extern short IsectLL2Ds(short *v1, short *v2, short *v3, short *v4);
 
+/*point in tri,  0 no intersection, 1 intersect */
+int IsectPT2Df(float pt[2], float v1[2], float v2[2], float v3[2]);
+/* point in quad,  0 no intersection, 1 intersect */
+int IsectPQ2Df(float pt[2], float v1[2], float v2[2], float v3[2], float v4[2]);
+
 /* interpolation weights of point in a triangle or quad, v4 may be NULL */
 void InterpWeightsQ3Dfl(float *v1, float *v2, float *v3, float *v4, float *co, float *w);
 /* interpolation weights of point in a polygon with >= 3 vertices */
@@ -358,7 +364,13 @@ void LocQuatSizeToMat4(float mat[][4], float loc[3], float quat[4], float size[3
 void tubemap(float x, float y, float z, float *u, float *v);
 void spheremap(float x, float y, float z, float *u, float *v);
 
-int LineIntersectsTriangle(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda);
+int LineIntersectsTriangle(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv);
+int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, float v0[3], float v1[3], float v2[3], float *lambda, float *ipoint);
+int AxialLineIntersectsTriangle(int axis, float co1[3], float co2[3], float v0[3], float v1[3], float v2[3], float *lambda);
+int AabbIntersectAabb(float min1[3], float max1[3], float min2[3], float max2[3]);
+void VecfCubicInterpol(float *x1, float *v1, float *x2, float *v2, float t, float *x, float *v);
+void PointInQuad2DUV(float v0[2], float v1[2], float v2[2], float v3[2], float pt[2], float *uv);
+void PointInFace2DUV(int isquad, float v0[2], float v1[2], float v2[2], float v3[2], float pt[2], float *uv);
 int point_in_tri_prism(float p[3], float v1[3], float v2[3], float v3[3]);
 
 float lambda_cp_line_ex(float p[3], float l1[3], float l2[3], float cp[3]);
diff --git a/source/blender/blenlib/BLI_kdtree.h b/source/blender/blenlib/BLI_kdtree.h
new file mode 100644
index 00000000000..d32e85c150c
--- /dev/null
+++ b/source/blender/blenlib/BLI_kdtree.h
@@ -0,0 +1,63 @@
+/**
+ * A kd-tree for nearest neighbour search.
+ * 
+ * $Id$
+ *
+ * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version. The Blender
+ * Foundation also sells licenses for use in proprietary software under
+ * the Blender License.  See http://www.blender.org/BL/ for information
+ * about this.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: none of this file.
+ *
+ * Contributor(s): Janne Karhu
+ *                 Brecht Van Lommel
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ */
+ 
+#ifndef BLI_KDTREE_H
+#define BLI_KDTREE_H
+
+struct KDTree;
+typedef struct KDTree KDTree;
+
+typedef struct KDTreeNearest {
+	int index;
+	float dist;
+	float co[3];
+} KDTreeNearest;
+
+/* Creates or free a kdtree */
+KDTree* BLI_kdtree_new(int maxsize);
+void BLI_kdtree_free(KDTree *tree);
+
+/* Construction: first insert points, then call balance. Normal is optional. */
+void BLI_kdtree_insert(KDTree *tree, int index, float *co, float *nor);
+void BLI_kdtree_balance(KDTree *tree);
+
+/* Find nearest returns index, and -1 if no node is found.
+ * Find n nearest returns number of points found, with results in nearest.
+ * Normal is optional. */
+int	BLI_kdtree_find_nearest(KDTree *tree, float *co, float *nor, KDTreeNearest *nearest);
+int	BLI_kdtree_find_n_nearest(KDTree *tree, int n, float *co, float *nor, KDTreeNearest *nearest);
+
+#endif
+
diff --git a/source/blender/blenlib/intern/BLI_kdtree.c b/source/blender/blenlib/intern/BLI_kdtree.c
new file mode 100644
index 00000000000..dbb9e95d379
--- /dev/null
+++ b/source/blender/blenlib/intern/BLI_kdtree.c
@@ -0,0 +1,327 @@
+/**
+ * $Id$
+ *
+ * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version. The Blender
+ * Foundation also sells licenses for use in proprietary software under
+ * the Blender License.  See http://www.blender.org/BL/ for information
+ * about this.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: none of this file.
+ *
+ * Contributor(s): Janne Karhu
+ *                 Brecht Van Lommel
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ */
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_arithb.h"
+#include "BLI_kdtree.h"
+
+#define SWAP(type, a, b) { type sw_ap; sw_ap=(a); (a)=(b); (b)=sw_ap; }
+
+typedef struct KDTreeNode {
+	struct KDTreeNode *left, *right;
+	float co[3], nor[3];
+	int index;
+	short d;
+} KDTreeNode;
+
+struct KDTree {
+	KDTreeNode *nodes;
+	int totnode;
+	KDTreeNode *root;
+};
+
+KDTree *BLI_kdtree_new(int maxsize)
+{
+	KDTree *tree;
+
+	tree= MEM_callocN(sizeof(KDTree), "KDTree");
+	tree->nodes= MEM_callocN(sizeof(KDTreeNode)*maxsize, "KDTreeNode");
+	tree->totnode= 0;
+
+	return tree;
+}
+
+void BLI_kdtree_free(KDTree *tree)
+{
+	if(tree) {
+		MEM_freeN(tree->nodes);
+		MEM_freeN(tree);
+	}
+}
+
+void BLI_kdtree_insert(KDTree *tree, int index, float *co, float *nor)
+{
+	KDTreeNode *node= &tree->nodes[tree->totnode++];
+
+	node->index= index;
+	VecCopyf(node->co, co);
+	if(nor) VecCopyf(node->nor, nor);
+}
+
+static KDTreeNode *kdtree_balance(KDTreeNode *nodes, int totnode, int axis)
+{
+	KDTreeNode *node;
+	float co;
+	int left, right, median, i, j;
+
+	if(totnode <= 0)
+		return NULL;
+	else if(totnode == 1)
+		return nodes;
+	
+	/* quicksort style sorting around median */
+	left= 0;
+	right= totnode-1;
+	median= totnode/2;
+
+	while(right > left) {
+		co= nodes[right].co[axis];
+		i= left-1;
+		j= right;
+
+		while(1) {
+			while(nodes[++i].co[axis] < co);
+			while(nodes[--j].co[axis] > co && j>left);
+
+			if(i >= j) break;
+			SWAP(KDTreeNode, nodes[i], nodes[j]);
+		}
+
+		SWAP(KDTreeNode, nodes[i], nodes[right]);
+		if(i >= median)
+			right= i-1;
+		if(i <= median)
+			left= i+1;
+	}
+
+	/* set node and sort subnodes */
+	node= &nodes[median];
+	node->d= axis;
+	node->left= kdtree_balance(nodes, median, (axis+1)%3);
+	node->right= kdtree_balance(nodes+median+1, (totnode-(median+1)), (axis+1)%3);
+
+	return node;
+}
+
+void BLI_kdtree_balance(KDTree *tree)
+{
+	tree->root= kdtree_balance(tree->nodes, tree->totnode, 0);
+}
+
+static float squared_distance(float *v1, float *v2, float *n1, float *n2)
+{
+	float d[3], dist;
+
+	d[0]= v2[0]-v1[0];
+	d[1]= v2[1]-v1[1];
+	d[2]= v2[2]-v1[2];
+
+	dist= d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+
+	if(n1 && n2 && n1[0]*n2[0] + n1[1]*n2[1] + n1[2]*n2[2] < 0.0f)
+		dist *= 10.0f;
+
+	return dist;
+}
+
+int	BLI_kdtree_find_nearest(KDTree *tree, float *co, float *nor, KDTreeNearest *nearest)
+{
+	KDTreeNode *root, *node, *min_node;
+	KDTreeNode **stack, *defaultstack[100];
+	float min_dist, cur_dist;
+	int totstack, cur=0;
+
+	if(!tree->root)
+		return -1;
+
+	stack= defaultstack;
+	totstack= 100;
+
+	root= tree->root;
+	min_node= root;
+	min_dist= squared_distance(root->co,co,root->nor,nor);
+
+	if(root->left)
+		stack[cur++]=root->left;
+
+	if(root->right)
+		stack[cur++]=root->right;
+	
+	while(cur--){
+		node=stack[cur];
+
+		cur_dist = node->co[node->d] - co[node->d];
+
+		if(cur_dist<0.0){
+			cur_dist= -cur_dist*cur_dist;
+
+			if(-cur_dist<min_dist){
+				cur_dist=squared_distance(node->co,co,node->nor,nor);
+				if(cur_dist<min_dist){
+					min_dist=cur_dist;
+					min_node=node;
+				}
+				if(node->left)
+					stack[cur++]=node->left;
+			}
+			if(node->right)
+				stack[cur++]=node->right;
+		}
+		else{
+			cur_dist= cur_dist*cur_dist;
+
+			if(cur_dist<min_dist){
+				cur_dist=squared_distance(node->co,co,node->nor,nor);
+				if(cur_dist<min_dist){
+					min_dist=cur_dist;
+					min_node=node;
+				}
+				if(node->right)
+					stack[cur++]=node->right;
+			}
+			if(node->left)
+				stack[cur++]=node->left;
+		}
+		if(cur+3 > totstack){
+			KDTreeNode **temp=MEM_callocN((totstack+100)*sizeof(KDTreeNode*), "psys_treestack");
+			memcpy(temp,stack,totstack*sizeof(KDTreeNode*));
+			if(stack != defaultstack)
+				MEM_freeN(stack);
+			stack=temp;
+			totstack+=100;
+		}
+	}
+
+	if(nearest) {
+		nearest->index= min_node->index;
+		nearest->dist= sqrt(min_dist);
+		VecCopyf(nearest->co, min_node->co);
+	}
+
+	if(stack != defaultstack)
+		MEM_freeN(stack);
+
+	return min_node->index;
+}
+
+static void add_nearest(KDTreeNearest *ptn, int *found, int n, int index, float dist, float *co)
+{
+	int i;
+
+	if(*found<n) (*found)++;
+
+	for(i=*found-1; i>0; i--) {
+		if(dist >= ptn[i-1].dist)
+			break;
+		else
+			ptn[i]= ptn[i-1];
+	}
+
+	ptn[i].index= index;
+	ptn[i].dist= dist;
+	VecCopyf(ptn[i].co, co);
+}
+
+/* finds the nearest n entries in tree to specified coordinates */
+int	BLI_kdtree_find_n_nearest(KDTree *tree, int n, float *co, float *nor, KDTreeNearest *nearest)
+{
+	KDTreeNode *root, *node=0;
+	KDTreeNode **stack, *defaultstack[100];
+	float cur_dist;
+	int i, totstack, cur=0, found=0;
+
+	if(!tree->root)
+		return 0;
+
+	stack= defaultstack;
+	totstack= 100;
+
+	root= tree->root;
+
+	cur_dist= squared_distance(root->co,co,root->nor,nor);
+	add_nearest(nearest,&found,n,root->index,cur_dist,root->co);
+	
+	if(root->left)
+		stack[cur++]=root->left;
+
+	if(root->right)
+		stack[cur++]=root->right;
+
+	while(cur--){
+		node=stack[cur];
+
+		cur_dist = node->co[node->d] - co[node->d];
+
+		if(cur_dist<0.0){
+			cur_dist= -cur_dist*cur_dist;
+
+			if(found<n || -cur_dist<nearest[found-1].dist){
+				cur_dist=squared_distance(node->co,co,node->nor,nor);
+
+				if(found<n || cur_dist<nearest[found-1].dist)
+					add_nearest(nearest,&found,n,node->index,cur_dist,node->co);
+
+				if(node->left)
+					stack[cur++]=node->left;
+			}
+			if(node->right)
+				stack[cur++]=node->right;
+		}
+		else{
+			cur_dist= cur_dist*cur_dist;
+
+			if(found<n || cur_dist<nearest[found-1].dist){
+				cur_dist=squared_distance(node->co,co,node->nor,nor);
+				if(found<n || cur_dist<nearest[found-1].dist)
+					add_nearest(nearest,&found,n,node->index,cur_dist,node->co);
+
+				if(node->right)
+					stack[cur++]=node->right;
+			}
+			if(node->left)
+				stack[cur++]=node->left;
+		}
+		if(cur+3 > totstack){
+			KDTreeNode **temp=MEM_callocN((totstack+100)*sizeof(KDTreeNode*), "psys_treestack");
+			memcpy(temp,stack,totstack*sizeof(KDTreeNode*));
+			if(stack != defaultstack)
+				MEM_freeN(stack);
+			stack=temp;
+			totstack+=100;
+		}
+	}
+
+	for(i=0; i<found; i++)
+		nearest[i].dist= sqrt(nearest[i].dist);
+
+	if(stack != defaultstack)
+		MEM_freeN(stack);
+
+	return found;
+}
+
diff --git a/source/blender/blenlib/intern/arithb.c b/source/blender/blenlib/intern/arithb.c
index 2888a684af1..15ea647361f 100644
--- a/source/blender/blenlib/intern/arithb.c
+++ b/source/blender/blenlib/intern/arithb.c
@@ -2421,6 +2421,98 @@ short IsectLL2Df(float *v1, float *v2, float *v3, float *v4)
 	return 0;
 }
 
+/*
+-1: colliniar
+ 1: intersection
+
+*/
+short IsectLLPt2Df(float x0,float y0,float x1,float y1,
+					 float x2,float y2,float x3,float y3, float *xi,float *yi)
+
+{
+	/*
+	this function computes the intersection of the sent lines
+	and returns the intersection point, note that the function assumes
+	the lines intersect. the function can handle vertical as well
+	as horizontal lines. note the function isn't very clever, it simply
+	applies the math, but we don't need speed since this is a
+	pre-processing step
+	*/
+	float c1,c2, // constants of linear equations
+	det_inv,  // the inverse of the determinant of the coefficient
+	m1,m2;    // the slopes of each line
+	/*
+	compute slopes, note the cludge for infinity, however, this will
+	be close enough
+	*/
+	if ( fabs( x1-x0 ) > 0.000001 )
+		m1 = ( y1-y0 ) / ( x1-x0 );
+	else
+		return -1; /*m1 = ( float ) 1e+10;*/   // close enough to infinity
+
+	if ( fabs( x3-x2 ) > 0.000001 )
+		m2 = ( y3-y2 ) / ( x3-x2 );
+	else
+		return -1; /*m2 = ( float ) 1e+10;*/   // close enough to infinity
+
+	if (fabs(m1-m2) < 0.000001)
+		return -1; /* paralelle lines */
+	
+// compute constants
+
+	c1 = ( y0-m1*x0 );
+	c2 = ( y2-m2*x2 );
+
+// compute the inverse of the determinate
+
+	det_inv = 1.0f / ( -m1 + m2 );
+
+// use Kramers rule to compute xi and yi
+
+	*xi= ( ( -c2 + c1 ) *det_inv );
+	*yi= ( ( m2*c1 - m1*c2 ) *det_inv );
+	
+	return 1; 
+} // end Intersect_Lines
+
+#define SIDE_OF_LINE(pa,pb,pp)	((pa[0]-pp[0])*(pb[1]-pp[1]))-((pb[0]-pp[0])*(pa[1]-pp[1]))
+#define ISECT_EPSILON 1e-6
+
+/* point in tri */
+int IsectPT2Df(float pt[2], float v1[2], float v2[2], float v3[2])
+{
+	if ((SIDE_OF_LINE(v1,v2,pt)>=-ISECT_EPSILON) &&
+		(SIDE_OF_LINE(v2,v3,pt)>=-ISECT_EPSILON) &&
+		(SIDE_OF_LINE(v3,v1,pt)>=-ISECT_EPSILON))
+		return 1;
+	else {
+		return 0;
+	}
+}
+/* point in quad - only convex quads */
+int IsectPQ2Df(float pt[2], float v1[2], float v2[2], float v3[2], float v4[2])
+{
+	if ((SIDE_OF_LINE(v1,v2,pt)>=-ISECT_EPSILON) &&
+		(SIDE_OF_LINE(v2,v3,pt)>=-ISECT_EPSILON) &&
+		(SIDE_OF_LINE(v3,v4,pt)>=-ISECT_EPSILON) &&
+		(SIDE_OF_LINE(v4,v1,pt)>=-ISECT_EPSILON))
+		return 1;
+	else
+		return 0;
+}
+
+
+
+			/* copied from Geometry.c - todo - move to arithb.c or some other generic place we can reuse */
+#define SIDE_OF_LINE(pa,pb,pp)	((pa[0]-pp[0])*(pb[1]-pp[1]))-((pb[0]-pp[0])*(pa[1]-pp[1]))
+#define POINT_IN_TRI(p0,p1,p2,p3)	((SIDE_OF_LINE(p1,p2,p0)>=0) && (SIDE_OF_LINE(p2,p3,p0)>=0) && (SIDE_OF_LINE(p3,p1,p0)>=0))
+
+/**
+ * 
+ * @param min 
+ * @param max 
+ * @param vec 
+ */
 void MinMax3(float *min, float *max, float *vec)
 {
 	if(min[0]>vec[0]) min[0]= vec[0];
@@ -3035,6 +3127,11 @@ float Vec2Lenf(float *v1, float *v2)
 	return (float)sqrt(x*x+y*y);
 }
 
+float Vec2Length(float *v)
+{
+	return (float)sqrt(v[0]*v[0] + v[1]*v[1]);
+}
+
 void Vec2Mulf(float *v1, float f)
 {
 	v1[0]*= f;
@@ -3501,7 +3598,7 @@ void mul_v3_v3m4(float *v1, float *v2, float mat[][4])
    test if the line starting at p1 ending at p2 intersects the triangle v0..v2
    return non zero if it does 
 */
-int LineIntersectsTriangle(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda)
+int LineIntersectsTriangle(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv)
 {
 
 	float p[3], s[3], d[3], e1[3], e2[3], q[3];
@@ -3527,15 +3624,318 @@ int LineIntersectsTriangle(float p1[3], float p2[3], float v0[3], float v1[3], f
 	
 	v = f * Inpf(d, q);
 	if ((v < 0.0)||((u + v) > 1.0)) return 0;
+
+	if(uv) {
+		uv[0]= u;
+		uv[1]= v;
+	}
 	
 	return 1;
 }
 
+/* Adapted from the paper by Kasper Fauerby */
+/* "Improved Collision detection and Response" */
+int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, float v0[3], float v1[3], float v2[3], float *lambda, float *ipoint)
+{
+	float e1[3], e2[3], e3[3], point[3], vel[3], /*dist[3],*/ nor[3], temp[3], bv[3];
+	float a, b, c, d, e, x, y, z, t, t0, t1, radius2=radius*radius;
+	float elen2,edotv,edotbv,nordotv,vel2;
+	int embedded_in_plane=0, found_by_sweep=0;
 
-/*
-find closest point to p on line through l1,l2
-and return lambda, where (0 <= lambda <= 1) when cp is in the line segement l1,l2
-*/
+	VecSubf(e1,v1,v0);
+	VecSubf(e2,v2,v0);
+	VecSubf(vel,p2,p1);
+
+/*---test plane of tri---*/
+	Crossf(nor,e1,e2);
+	Normalize(nor);
+	/* flip normal */
+	if(Inpf(nor,vel)>0.0f) VecMulf(nor,-1.0f);
+	
+	a=Inpf(p1,nor)-Inpf(v0,nor);
+
+	nordotv=Inpf(nor,vel);
+
+	if ((nordotv > -0.000001) && (nordotv < 0.000001)) {
+		if(fabs(a)>=1.0f)
+			return 0;
+		else{
+			embedded_in_plane=1;
+			t0=0.0f;
+			t1=1.0f;
+		}
+	}
+	else{
+		t0=(radius-a)/nordotv;
+		t1=(-radius-a)/nordotv;
+		/* make t0<t1 */
+		if(t0>t1){b=t1; t1=t0; t0=b;}
+
+		if(t0>1.0f || t1<0.0f) return 0;
+
+		/* clamp to [0,1] */
+		t0=(t0<0.0f)?0.0f:((t0>1.0f)?1.0:t0);
+		t1=(t1<0.0f)?0.0f:((t1>1.0f)?1.0:t1);
+	}
+
+/*---test inside of tri---*/
+	if(embedded_in_plane==0){
+		/* plane intersection point */
+		VecCopyf(point,vel);
+		VecMulf(point,t0);
+		VecAddf(point,point,p1);
+		VecCopyf(temp,nor);
+		VecMulf(temp,radius);
+		VecSubf(point,point,temp);
+
+		/* is the point in the tri? */
+		a=Inpf(e1,e1);
+		b=Inpf(e1,e2);
+		c=Inpf(e2,e2);
+
+		VecSubf(temp,point,v0);
+		d=Inpf(temp,e1);
+		e=Inpf(temp,e2);
+		
+		x=d*c-e*b;
+		y=e*a-d*b;
+		z=x+y-(a*c-b*b);
+
+		if(( ((unsigned int)z)& ~(((unsigned int)x)|((unsigned int)y)) ) & 0x80000000){
+			*lambda=t0;
+			VecCopyf(ipoint,point);
+			return 1;
+		}
+	}
+
+	*lambda=1.0f;
+/*---test points---*/
+	a=vel2=Inpf(vel,vel);
+
+	/*v0*/
+	VecSubf(temp,p1,v0);
+	b=2.0f*Inpf(vel,temp);
+	c=Inpf(temp,temp)-radius2;
+	d=b*b-4*a*c;
+
+	if(d>=0.0f){
+		if(d==0.0f)
+			t=-b/2*a;
+		else{
+			z=sqrt(d);
+			x=(-b-z)*0.5/a;
+			y=(-b+z)*0.5/a;
+			t=x<y?x:y;
+		}
+
+		if(t>0.0 && t < *lambda){
+			*lambda=t;
+			VecCopyf(ipoint,v0);
+			found_by_sweep=1;
+		}
+	}
+
+	/*v1*/
+	VecSubf(temp,p1,v1);
+	b=2.0f*Inpf(vel,temp);
+	c=Inpf(temp,temp)-radius2;
+	d=b*b-4*a*c;
+
+	if(d>=0.0f){
+		if(d==0.0f)
+			t=-b/2*a;
+		else{
+			z=sqrt(d);
+			x=(-b-z)*0.5/a;
+			y=(-b+z)*0.5/a;
+			t=x<y?x:y;
+		}
+
+		if(t>0.0 && t < *lambda){
+			*lambda=t;
+			VecCopyf(ipoint,v1);
+			found_by_sweep=1;
+		}
+	}
+	/*v2*/
+	VecSubf(temp,p1,v2);
+	b=2.0f*Inpf(vel,temp);
+	c=Inpf(temp,temp)-radius2;
+	d=b*b-4*a*c;
+
+	if(d>=0.0f){
+		if(d==0.0f)
+			t=-b/2*a;
+		else{
+			z=sqrt(d);
+			x=(-b-z)*0.5/a;
+			y=(-b+z)*0.5/a;
+			t=x<y?x:y;
+		}
+
+		if(t>0.0 && t < *lambda){
+			*lambda=t;
+			VecCopyf(ipoint,v2);
+			found_by_sweep=1;
+		}
+	}
+
+/*---test edges---*/
+	/*e1*/
+	VecSubf(bv,v0,p1);
+	elen2 = Inpf(e1,e1);
+	edotv = Inpf(e1,vel);
+	edotbv = Inpf(e1,bv);
+
+	a=elen2*(-Inpf(vel,vel))+edotv*edotv;
+	b=2.0f*(elen2*Inpf(vel,bv)-edotv*edotbv);
+	c=elen2*(radius2-Inpf(bv,bv))+edotbv*edotbv;
+	d=b*b-4*a*c;
+	if(d>=0.0f){
+		if(d==0.0f)
+			t=-b/2*a;
+		else{
+			z=sqrt(d);
+			x=(-b-z)*0.5/a;
+			y=(-b+z)*0.5/a;
+			t=x<y?x:y;
+		}
+
+		e=(edotv*t-edotbv)/elen2;
+
+		if((e>=0.0f) && (e<=1.0f)){
+			if(t>0.0 && t < *lambda){
+				*lambda=t;
+				VecCopyf(ipoint,e1);
+				VecMulf(ipoint,e);
+				VecAddf(ipoint,ipoint,v0);
+				found_by_sweep=1;
+			}
+		}
+	}
+
+	/*e2*/
+	/*bv is same*/
+	elen2 = Inpf(e2,e2);
+	edotv = Inpf(e2,vel);
+	edotbv = Inpf(e2,bv);
+
+	a=elen2*(-Inpf(vel,vel))+edotv*edotv;
+	b=2.0f*(elen2*Inpf(vel,bv)-edotv*edotbv);
+	c=elen2*(radius2-Inpf(bv,bv))+edotbv*edotbv;
+	d=b*b-4*a*c;
+	if(d>=0.0f){
+		if(d==0.0f)
+			t=-b/2*a;
+		else{
+			z=sqrt(d);
+			x=(-b-z)*0.5/a;
+			y=(-b+z)*0.5/a;
+			t=x<y?x:y;
+		}
+
+		e=(edotv*t-edotbv)/elen2;
+
+		if((e>=0.0f) && (e<=1.0f)){
+			if(t>0.0 && t < *lambda){
+				*lambda=t;
+				VecCopyf(ipoint,e2);
+				VecMulf(ipoint,e);
+				VecAddf(ipoint,ipoint,v0);
+				found_by_sweep=1;
+			}
+		}
+	}
+
+	/*e3*/
+	VecSubf(e3,v2,v1);
+	VecSubf(bv,v1,p1);
+	elen2 = Inpf(e3,e3);
+	edotv = Inpf(e3,vel);
+	edotbv = Inpf(e3,bv);
+
+	a=elen2*(-Inpf(vel,vel))+edotv*edotv;
+	b=2.0f*(elen2*Inpf(vel,bv)-edotv*edotbv);
+	c=elen2*(radius2-Inpf(bv,bv))+edotbv*edotbv;
+	d=b*b-4*a*c;
+	if(d>=0.0f){
+		if(d==0.0f)
+			t=-b/2*a;
+		else{
+			z=sqrt(d);
+			x=(-b-z)*0.5/a;
+			y=(-b+z)*0.5/a;
+			t=x<y?x:y;
+		}
+
+		e=(edotv*t-edotbv)/elen2;
+
+		if((e>=0.0f) && (e<=1.0f)){
+			if(t>0.0 && t < *lambda){
+				*lambda=t;
+				VecCopyf(ipoint,e3);
+				VecMulf(ipoint,e);
+				VecAddf(ipoint,ipoint,v1);
+				found_by_sweep=1;
+			}
+		}
+	}
+
+	return found_by_sweep;
+}
+int AxialLineIntersectsTriangle(int axis, float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda)
+{
+	float p[3], e1[3], e2[3];
+	float u, v, f;
+	int a0=axis, a1=(axis+1)%3, a2=(axis+2)%3;
+
+	//return LineIntersectsTriangle(p1,p2,v0,v1,v2,lambda);
+
+	///* first a simple bounding box test */
+	//if(MIN3(v0[a1],v1[a1],v2[a1]) > p1[a1]) return 0;
+	//if(MIN3(v0[a2],v1[a2],v2[a2]) > p1[a2]) return 0;
+	//if(MAX3(v0[a1],v1[a1],v2[a1]) < p1[a1]) return 0;
+	//if(MAX3(v0[a2],v1[a2],v2[a2]) < p1[a2]) return 0;
+
+	///* then a full intersection test */
+	
+	VecSubf(e1,v1,v0);
+	VecSubf(e2,v2,v0);
+	VecSubf(p,v0,p1);
+
+	f= (e2[a1]*e1[a2]-e2[a2]*e1[a1]);
+	if ((f > -0.000001) && (f < 0.000001)) return 0;
+
+	v= (p[a2]*e1[a1]-p[a1]*e1[a2])/f;
+	if ((v < 0.0)||(v > 1.0)) return 0;
+	
+	f= e1[a1];
+	if((f > -0.000001) && (f < 0.000001)){
+		f= e1[a2];
+		if((f > -0.000001) && (f < 0.000001)) return 0;
+		u= (-p[a2]-v*e2[a2])/f;
+	}
+	else
+		u= (-p[a1]-v*e2[a1])/f;
+
+	if ((u < 0.0)||((u + v) > 1.0)) return 0;
+
+	*lambda = (p[a0]+u*e1[a0]+v*e2[a0])/(p2[a0]-p1[a0]);
+
+	if ((*lambda < 0.0)||(*lambda > 1.0)) return 0;
+
+	return 1;
+}
+
+int AabbIntersectAabb(float min1[3], float max1[3], float min2[3], float max2[3])
+{
+	return (min1[0]<max2[0] && min1[1]<max2[1] && min1[2]<max2[2] &&
+	        min2[0]<max1[0] && min2[1]<max1[1] && min2[2]<max1[2]);
+}
+
+/* find closest point to p on line through l1,l2 and return lambda,
+ * where (0 <= lambda <= 1) when cp is in the line segement l1,l2
+ */
 float lambda_cp_line_ex(float p[3], float l1[3], float l2[3], float cp[3])
 {
 	float h[3],u[3],lambda;
@@ -3547,6 +3947,7 @@ float lambda_cp_line_ex(float p[3], float l1[3], float l2[3], float cp[3])
 	cp[2] = l1[2] + u[2] * lambda;
 	return lambda;
 }
+
 /* little sister we only need to know lambda */
 float lambda_cp_line(float p[3], float l1[3], float l2[3])
 {
@@ -3556,7 +3957,156 @@ float lambda_cp_line(float p[3], float l1[3], float l2[3])
 	return(Inpf(u,h)/Inpf(u,u));
 }
 
+/* Similar to LineIntersectsTriangleUV, except it operates on a quad and in 2d, assumes point is in quad */
+void PointInQuad2DUV(float v0[2], float v1[2], float v2[2], float v3[2], float pt[2], float *uv)
+{
+	float x0,y0, x1,y1, wtot, v2d[2], w1, w2;
+	
+	/* used for paralelle lines */
+	float pt3d[3], l1[3], l2[3], pt_on_line[3];
+	
+	/* compute 2 edges  of the quad  intersection point */
+	if (IsectLLPt2Df(v0[0],v0[1],v1[0],v1[1],  v2[0],v2[1],v3[0],v3[1], &x0,&y0) == 1) {
+		/* the intersection point between the quad-edge intersection and the point in the quad we want the uv's for */
+		/* should never be paralle !! */
+		/*printf("\tnot paralelle 1\n");*/
+		IsectLLPt2Df(pt[0],pt[1],x0,y0,  v0[0],v0[1],v3[0],v3[1], &x1,&y1);
+		
+		/* Get the weights from the new intersection point, to each edge */
+		v2d[0] = x1-v0[0];
+		v2d[1] = y1-v0[1];
+		w1 = Vec2Length(v2d);
+		
+		v2d[0] = x1-v3[0]; /* some but for the other vert */
+		v2d[1] = y1-v3[1];
+		w2 = Vec2Length(v2d);
+		wtot = w1+w2;
+		/*w1 = w1/wtot;*/
+		/*w2 = w2/wtot;*/
+		uv[0] = w1/wtot;
+	} else {
+		/* lines are paralelle, lambda_cp_line_ex is 3d grrr */
+		/*printf("\tparalelle1\n");*/
+		pt3d[0] = pt[0];
+		pt3d[1] = pt[1];
+		pt3d[2] = l1[2] = l2[2] = 0.0f;
+		
+		l1[0] = v0[0]; l1[1] = v0[1];
+		l2[0] = v1[0]; l2[1] = v1[1];
+		lambda_cp_line_ex(pt3d, l1, l2, pt_on_line);
+		v2d[0] = pt[0]-pt_on_line[0]; /* same, for the other vert */
+		v2d[1] = pt[1]-pt_on_line[1];
+		w1 = Vec2Length(v2d);
+		
+		l1[0] = v2[0]; l1[1] = v2[1];
+		l2[0] = v3[0]; l2[1] = v3[1];
+		lambda_cp_line_ex(pt3d, l1, l2, pt_on_line);
+		v2d[0] = pt[0]-pt_on_line[0]; /* same, for the other vert */
+		v2d[1] = pt[1]-pt_on_line[1];
+		w2 = Vec2Length(v2d);
+		wtot = w1+w2;
+		uv[0] = w1/wtot;	
+	}
+	
+	/* Same as above to calc the uv[1] value, alternate calculation */
+	
+	if (IsectLLPt2Df(v0[0],v0[1],v3[0],v3[1],  v1[0],v1[1],v2[0],v2[1], &x0,&y0) == 1) { /* was v0,v1  v2,v3  now v0,v3  v1,v2*/
+		/* never paralle if above was not */
+		/*printf("\tnot paralelle2\n");*/
+		IsectLLPt2Df(pt[0],pt[1],x0,y0,  v0[0],v0[1],v1[0],v1[1], &x1,&y1);/* was v0,v3  now v0,v1*/
+		
+		v2d[0] = x1-v0[0];
+		v2d[1] = y1-v0[1];
+		w1 = Vec2Length(v2d);
+		
+		v2d[0] = x1-v1[0];
+		v2d[1] = y1-v1[1];
+		w2 = Vec2Length(v2d);
+		wtot = w1+w2;
+		uv[1] = w1/wtot;
+	} else {
+		/* lines are paralelle, lambda_cp_line_ex is 3d grrr */
+		/*printf("\tparalelle2\n");*/
+		pt3d[0] = pt[0];
+		pt3d[1] = pt[1];
+		pt3d[2] = l1[2] = l2[2] = 0.0f;
+		
+		
+		l1[0] = v0[0]; l1[1] = v0[1];
+		l2[0] = v3[0]; l2[1] = v3[1];
+		lambda_cp_line_ex(pt3d, l1, l2, pt_on_line);
+		v2d[0] = pt[0]-pt_on_line[0]; /* some but for the other vert */
+		v2d[1] = pt[1]-pt_on_line[1];
+		w1 = Vec2Length(v2d);
+		
+		l1[0] = v1[0]; l1[1] = v1[1];
+		l2[0] = v2[0]; l2[1] = v2[1];
+		lambda_cp_line_ex(pt3d, l1, l2, pt_on_line);
+		v2d[0] = pt[0]-pt_on_line[0]; /* some but for the other vert */
+		v2d[1] = pt[1]-pt_on_line[1];
+		w2 = Vec2Length(v2d);
+		wtot = w1+w2;
+		uv[1] = w1/wtot;
+	}
+	/* may need to flip UV's here */
+}
+
+/* same as above but does tri's and quads, tri's are a bit of a hack */
+void PointInFace2DUV(int isquad, float v0[2], float v1[2], float v2[2], float v3[2], float pt[2], float *uv)
+{
+	if (isquad) {
+		PointInQuad2DUV(v0, v1, v2, v3, pt, uv);
+	}
+	else {
+		/* not for quads, use for our abuse of LineIntersectsTriangleUV */
+		float p1_3d[3], p2_3d[3], v0_3d[3], v1_3d[3], v2_3d[3], lambda;
+			
+		p1_3d[0] = p2_3d[0] = uv[0];
+		p1_3d[1] = p2_3d[1] = uv[1];
+		p1_3d[2] = 1.0f;
+		p2_3d[2] = -1.0f;
+		v0_3d[2] = v1_3d[2] = v2_3d[2] = 0.0;
+		
+		/* generate a new fuv, (this is possibly a non optimal solution,
+		 * since we only need 2d calculation but use 3d func's)
+		 * 
+		 * this method makes an imaginary triangle in 2d space using the UV's from the derived mesh face
+		 * Then find new uv coords using the fuv and this face with LineIntersectsTriangleUV.
+		 * This means the new values will be correct in relation to the derived meshes face. 
+		 */
+		Vec2Copyf(v0_3d, v0);
+		Vec2Copyf(v1_3d, v1);
+		Vec2Copyf(v2_3d, v2);
+		
+		/* Doing this in 3D is not nice */
+		LineIntersectsTriangle(p1_3d, p2_3d, v0_3d, v1_3d, v2_3d, &lambda, uv);
+	}
+}
+
+/* (x1,v1)(t1=0)------(x2,v2)(t2=1), 0<t<1 --> (x,v)(t) */
+void VecfCubicInterpol(float *x1, float *v1, float *x2, float *v2, float t, float *x, float *v)
+{
+	float a[3],b[3];
+	float t2= t*t;
+	float t3= t2*t;
 
+	/* cubic interpolation */
+	a[0]= v1[0] + v2[0] + 2*(x1[0] - x2[0]);
+	a[1]= v1[1] + v2[1] + 2*(x1[1] - x2[1]);
+	a[2]= v1[2] + v2[2] + 2*(x1[2] - x2[2]);
+
+	b[0]= -2*v1[0] - v2[0] - 3*(x1[0] - x2[0]);
+	b[1]= -2*v1[1] - v2[1] - 3*(x1[1] - x2[1]);
+	b[2]= -2*v1[2] - v2[2] - 3*(x1[2] - x2[2]);
+
+	x[0]= a[0]*t3 + b[0]*t2 + v1[0]*t + x1[0];
+	x[1]= a[1]*t3 + b[1]*t2 + v1[1]*t + x1[1];
+	x[2]= a[2]*t3 + b[2]*t2 + v1[2]*t + x1[2];
+
+	v[0]= 3*a[0]*t2 + 2*b[0]*t + v1[0];
+	v[1]= 3*a[1]*t2 + 2*b[1]*t + v1[1];
+	v[2]= 3*a[2]*t2 + 2*b[2]*t + v1[2];
+}
 
 int point_in_slice(float p[3], float v1[3], float l1[3], float l2[3])
 {