New particle collisions code:

* The old collisions code detected particle collisions by calculating the
  collision times analytically from the collision mesh faces. This was
  pretty accurate, but didn't support rotating/deforming faces at all, as
  the equations for these quickly become quite nasty.
* The new code uses a simple "distance to plane/edge/vert" function and
  iterates this with the Newton-Rhapson method to find the closest particle
  distance during a simulation step.
* The advantage in this is that the collision object can now move, rotate,
  scale or even deform freely and collisions are still detected reliably.
* For some extreme movements the calculation errors could stack up so much
  that the detection fails, but this can be easily fixed by increasing the
  particle size or simulation substeps.
* As a side note the algorithm doesn't really do point particles anymore,
  but uses a very small radius as the particle size when "size deflect" isn't
  selected.
* I've also updated the collision response code a bit, so now the particles
  shouldn't leak even from tight corners.

All in all the collisions code is now much cleaner and more robust than before!

2 files changed, 626 insertions, 406 deletions

diff --git a/source/blender/blenkernel/intern/boids.c b/source/blender/blenkernel/intern/boids.c
index 9fb7f6408ac..286cd3b2499 100644
--- a/source/blender/blenkernel/intern/boids.c
+++ b/source/blender/blenkernel/intern/boids.c
@@ -213,10 +213,8 @@ static int rule_avoid_collision(BoidRule *rule, BoidBrainData *bbd, BoidValues *
 		sub_v3_v3v3(ray_dir, col.co2, col.co1);
 		mul_v3_fl(ray_dir, acbr->look_ahead);
 		col.f = 0.0f;
-		col.cfra = fmod(bbd->cfra-bbd->dfra, 1.0f);
-		col.dfra = bbd->dfra;
 		hit.index = -1;
-		hit.dist = col.ray_len = len_v3(ray_dir);
+		hit.dist = col.original_ray_length = len_v3(ray_dir);
 
 		/* find out closest deflector object */
 		for(coll = bbd->sim->colliders->first; coll; coll=coll->next) {
@@ -224,18 +222,18 @@ static int rule_avoid_collision(BoidRule *rule, BoidBrainData *bbd, BoidValues *
 			if(coll->ob == bpa->ground)
 				continue;
 
-			col.ob = coll->ob;
+			col.current = coll->ob;
 			col.md = coll->collmd;
 
 			if(col.md && col.md->bvhtree)
-				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, BKE_psys_collision_neartest_cb, &col);
 		}
 		/* then avoid that object */
 		if(hit.index>=0) {
-			t = hit.dist/col.ray_len;
+			t = hit.dist/col.original_ray_length;
 
 			/* avoid head-on collision */
-			if(dot_v3v3(col.nor, pa->prev_state.ave) < -0.99) {
+			if(dot_v3v3(col.pce.nor, pa->prev_state.ave) < -0.99) {
 				/* don't know why, but uneven range [0.0,1.0] */
 				/* works much better than even [-1.0,1.0] */
 				bbd->wanted_co[0] = BLI_frand();
@@ -243,7 +241,7 @@ static int rule_avoid_collision(BoidRule *rule, BoidBrainData *bbd, BoidValues *
 				bbd->wanted_co[2] = BLI_frand();
 			}
 			else {
-				VECCOPY(bbd->wanted_co, col.nor);
+				copy_v3_v3(bbd->wanted_co, col.pce.nor);
 			}
 
 			mul_v3_fl(bbd->wanted_co, (1.0f - t) * val->personal_space * pa->size);
@@ -779,27 +777,26 @@ static Object *boid_find_ground(BoidBrainData *bbd, ParticleData *pa, float *gro
 		/* first try to find below boid */
 		copy_v3_v3(col.co1, pa->state.co);
 		sub_v3_v3v3(col.co2, pa->state.co, zvec);
-		sub_v3_v3(col.co2, zvec);
 		sub_v3_v3v3(ray_dir, col.co2, col.co1);
 		col.f = 0.0f;
-		col.cfra = fmod(bbd->cfra-bbd->dfra, 1.0f);
-		col.dfra = bbd->dfra;
 		hit.index = -1;
-		hit.dist = col.ray_len = len_v3(ray_dir);
+		hit.dist = col.original_ray_length = len_v3(ray_dir);
+		col.pce.inside = 0;
 
 		for(coll = bbd->sim->colliders->first; coll; coll = coll->next){
-			col.ob = coll->ob;
+			col.current = coll->ob;
 			col.md = coll->collmd;
+			col.fac1 = col.fac2 = 0.f;
 
 			if(col.md && col.md->bvhtree)
-				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, BKE_psys_collision_neartest_cb, &col);
 		}
 		/* then use that object */
 		if(hit.index>=0) {
-			t = hit.dist/col.ray_len;
+			t = hit.dist/col.original_ray_length;
 			interp_v3_v3v3(ground_co, col.co1, col.co2, t);
-			normalize_v3_v3(ground_nor, col.nor);
-			return col.hit_ob;
+			normalize_v3_v3(ground_nor, col.pce.nor);
+			return col.hit;
 		}
 
 		/* couldn't find below, so find upmost deflector object */
@@ -808,24 +805,22 @@ static Object *boid_find_ground(BoidBrainData *bbd, ParticleData *pa, float *gro
 		sub_v3_v3(col.co2, zvec);
 		sub_v3_v3v3(ray_dir, col.co2, col.co1);
 		col.f = 0.0f;
-		col.cfra = fmod(bbd->cfra-bbd->dfra, 1.0f);
-		col.dfra = bbd->dfra;
 		hit.index = -1;
-		hit.dist = col.ray_len = len_v3(ray_dir);
+		hit.dist = col.original_ray_length = len_v3(ray_dir);
 
 		for(coll = bbd->sim->colliders->first; coll; coll = coll->next){
-			col.ob = coll->ob;
+			col.current = coll->ob;
 			col.md = coll->collmd;
 
 			if(col.md && col.md->bvhtree)
-				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, BKE_psys_collision_neartest_cb, &col);
 		}
 		/* then use that object */
 		if(hit.index>=0) {
-			t = hit.dist/col.ray_len;
+			t = hit.dist/col.original_ray_length;
 			interp_v3_v3v3(ground_co, col.co1, col.co2, t);
-			normalize_v3_v3(ground_nor, col.nor);
-			return col.hit_ob;
+			normalize_v3_v3(ground_nor, col.pce.nor);
+			return col.hit;
 		}
 
 		/* default to z=0 */
diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c
index 5f7c62042ec..b5eb9b6ab29 100644
--- a/source/blender/blenkernel/intern/particle_system.c
+++ b/source/blender/blenkernel/intern/particle_system.c
@@ -2607,466 +2607,691 @@ static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, f
 /************************************************/
 /*			Collisions							*/
 /************************************************/
-/* convert from triangle barycentric weights to quad mean value weights */
-static void intersect_dm_quad_weights(float *v1, float *v2, float *v3, float *v4, float *w)
+#define COLLISION_MAX_COLLISIONS	10
+#define COLLISION_MIN_RADIUS 0.001f
+#define COLLISION_MIN_DISTANCE 0.0001f
+#define COLLISION_ZERO 0.00001f
+typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
+static float nr_signed_distance_to_plane(float *p, float radius, ParticleCollisionElement *pce, float *nor)
 {
-	float co[3], vert[4][3];
+	float p0[3], e1[3], e2[3], d;
 
-	VECCOPY(vert[0], v1);
-	VECCOPY(vert[1], v2);
-	VECCOPY(vert[2], v3);
-	VECCOPY(vert[3], v4);
+	sub_v3_v3v3(e1, pce->x1, pce->x0);
+	sub_v3_v3v3(e2, pce->x2, pce->x0);
+	sub_v3_v3v3(p0, p, pce->x0);
 
-	co[0]= v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2] + v4[0]*w[3];
-	co[1]= v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2] + v4[1]*w[3];
-	co[2]= v1[2]*w[0] + v2[2]*w[1] + v3[2]*w[2] + v4[2]*w[3];
+	cross_v3_v3v3(nor, e1, e2);
+	normalize_v3(nor);
 
-	interp_weights_poly_v3( w,vert, 4, co);
-}
+	d = dot_v3v3(p0, nor);
+
+	if(pce->inv_nor == -1) {
+		if(d < 0.f)
+			pce->inv_nor = 1;
+		else
+			pce->inv_nor = 0;
+	}
+
+	if(pce->inv_nor == 1) {
+		mul_v3_fl(nor, -1.f);
+		d = -d;
+	}
 
-/* check intersection with a derivedmesh */
-int psys_intersect_dm(Scene *scene, Object *ob, DerivedMesh *dm, float *vert_cos, float *co1, float* co2, float *min_d, int *min_face, float *min_w,
-						  float *face_minmax, float *pa_minmax, float radius, float *ipoint)
+	return d - radius;
+}
+static float nr_distance_to_edge(float *p, float radius, ParticleCollisionElement *pce, float *nor)
 {
-	MFace *mface=0;
-	MVert *mvert=0;
-	int i, totface, intersect=0;
-	float cur_d, cur_uv[2], v1[3], v2[3], v3[3], v4[3], min[3], max[3], p_min[3],p_max[3];
-	float cur_ipoint[3];
-	
-	if(dm==0){
-		psys_disable_all(ob);
+	float v0[3], v1[3], v2[3], c[3];
 
-		dm=mesh_get_derived_final(scene, ob, 0);
-		if(dm==0)
-			dm=mesh_get_derived_deform(scene, ob, 0);
+	sub_v3_v3v3(v0, pce->x1, pce->x0);
+	sub_v3_v3v3(v1, p, pce->x0);
+	sub_v3_v3v3(v2, p, pce->x1);
 
-		psys_enable_all(ob);
+	cross_v3_v3v3(c, v1, v2);
 
-		if(dm==0)
-			return 0;
+	return fabs(len_v3(c)/len_v3(v0)) - radius;
+}
+static float nr_distance_to_vert(float *p, float radius, ParticleCollisionElement *pce, float *nor)
+{
+	return len_v3v3(p, pce->x0) - radius;
+}
+static void collision_interpolate_element(ParticleCollisionElement *pce, float t, float fac, ParticleCollision *col)
+{
+	/* t is the current time for newton rhapson */
+	/* fac is the starting factor for current collision iteration */
+	/* the col->fac's are factors for the particle subframe step start and end during collision modifier step */
+	float f = fac + t*(1.f-fac);
+	float mul = col->fac1 + f * (col->fac2-col->fac1);
+	if(pce->tot > 0) {
+		madd_v3_v3v3fl(pce->x0, pce->x[0], pce->v[0], mul);
+
+		if(pce->tot > 1) {
+			madd_v3_v3v3fl(pce->x1, pce->x[1], pce->v[1], mul);
+
+			if(pce->tot > 2)
+				madd_v3_v3v3fl(pce->x2, pce->x[2], pce->v[2], mul);
+		}
 	}
+}
+static void collision_point_velocity(ParticleCollisionElement *pce)
+{
+	float v[3];
 
-	
+	copy_v3_v3(pce->vel, pce->v[0]);
 
-	if(pa_minmax==0){
-		INIT_MINMAX(p_min,p_max);
-		DO_MINMAX(co1,p_min,p_max);
-		DO_MINMAX(co2,p_min,p_max);
+	if(pce->tot > 1) {
+		sub_v3_v3v3(v, pce->v[1], pce->v[0]);
+		madd_v3_v3fl(pce->vel, v, pce->uv[0]);
+
+		if(pce->tot > 2) {
+			sub_v3_v3v3(v, pce->v[2], pce->v[0]);
+			madd_v3_v3fl(pce->vel, v, pce->uv[1]);
+		}
 	}
-	else{
-		VECCOPY(p_min,pa_minmax);
-		VECCOPY(p_max,pa_minmax+3);
+}
+static float collision_point_distance_with_normal(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *nor)
+{
+	if(fac >= 0.f)
+		collision_interpolate_element(pce, 0.f, fac, col);
+
+	switch(pce->tot) {
+		case 1:
+		{
+			sub_v3_v3v3(nor, p, pce->x0);
+			return normalize_v3(nor);
+		}
+		case 2:
+		{
+			float u, e[3], vec[3];
+			sub_v3_v3v3(e, pce->x1, pce->x0);
+			sub_v3_v3v3(vec, p, pce->x0);
+			u = dot_v3v3(vec, e) / dot_v3v3(e, e);
+
+			madd_v3_v3v3fl(nor, vec, e, -u);
+			return normalize_v3(nor);
+		}
+		case 3:
+			return nr_signed_distance_to_plane(p, 0.f, pce, nor);
 	}
+	return 0;
+}
+static void collision_point_on_surface(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *co)
+{
+	collision_interpolate_element(pce, 0.f, fac, col);
 
-	totface=dm->getNumFaces(dm);
-	mface=dm->getFaceDataArray(dm,CD_MFACE);
-	mvert=dm->getVertDataArray(dm,CD_MVERT);
-	
-	/* lets intersect the faces */
-	for(i=0; i<totface; i++,mface++){
-		if(vert_cos){
-			VECCOPY(v1,vert_cos+3*mface->v1);
-			VECCOPY(v2,vert_cos+3*mface->v2);
-			VECCOPY(v3,vert_cos+3*mface->v3);
-			if(mface->v4)
-				VECCOPY(v4,vert_cos+3*mface->v4)
+	switch(pce->tot) {
+		case 1:
+		{
+			sub_v3_v3v3(co, p, pce->x0);
+			normalize_v3(co);
+			madd_v3_v3v3fl(co, pce->x0, co, col->radius);
+			break;
 		}
-		else{
-			VECCOPY(v1,mvert[mface->v1].co);
-			VECCOPY(v2,mvert[mface->v2].co);
-			VECCOPY(v3,mvert[mface->v3].co);
-			if(mface->v4)
-				VECCOPY(v4,mvert[mface->v4].co)
+		case 2:
+		{
+			float u, e[3], vec[3], nor[3];
+			sub_v3_v3v3(e, pce->x1, pce->x0);
+			sub_v3_v3v3(vec, p, pce->x0);
+			u = dot_v3v3(vec, e) / dot_v3v3(e, e);
+
+			madd_v3_v3v3fl(nor, vec, e, -u);
+			normalize_v3(nor);
+
+			madd_v3_v3v3fl(co, pce->x0, e, pce->uv[0]);
+			madd_v3_v3fl(co, nor, col->radius);
+			break;
 		}
+		case 3:
+		{
+				float p0[3], e1[3], e2[3], nor[3];
 
-		if(face_minmax==0){
-			INIT_MINMAX(min,max);
-			DO_MINMAX(v1,min,max);
-			DO_MINMAX(v2,min,max);
-			DO_MINMAX(v3,min,max);
-			if(mface->v4)
-				DO_MINMAX(v4,min,max)
-			if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
-				continue;
+				sub_v3_v3v3(e1, pce->x1, pce->x0);
+				sub_v3_v3v3(e2, pce->x2, pce->x0);
+				sub_v3_v3v3(p0, p, pce->x0);
+
+				cross_v3_v3v3(nor, e1, e2);
+				normalize_v3(nor);
+
+				if(pce->inv_nor == 1)
+					mul_v3_fl(nor, -1.f);
+
+				madd_v3_v3v3fl(co, pce->x0, nor, col->radius);
+				madd_v3_v3fl(co, e1, pce->uv[0]);
+				madd_v3_v3fl(co, e2, pce->uv[1]);
+			break;
 		}
-		else{
-			VECCOPY(min, face_minmax+6*i);
-			VECCOPY(max, face_minmax+6*i+3);
-			if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
-				continue;
+	}
+}
+/* find first root in range [0-1] starting from 0 */
+static float collision_newton_rhapson(ParticleCollision *col, float radius, ParticleCollisionElement *pce, NRDistanceFunc distance_func)
+{
+	float t0, t1, d0, d1, dd, n[3];
+	int iter;
+
+	pce->inv_nor = -1;
+
+	/* start from the beginning */
+	t0 = 0.f;
+	collision_interpolate_element(pce, t0, col->f, col);
+	d0 = distance_func(col->co1, radius, pce, n);
+	t1 = 0.001f;
+	d1 = 0.f;
+
+	for(iter=0; iter<10; iter++) {//, itersum++) {
+		/* get current location */
+		collision_interpolate_element(pce, t1, col->f, col);
+		interp_v3_v3v3(pce->p, col->co1, col->co2, t1);
+
+		d1 = distance_func(pce->p, radius, pce, n);
+
+		/* no movement, so no collision */
+		if(d1 == d0) {
+			return -1.f;
 		}
 
-		if(radius>0.0f){
-			if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v2, v3, v1, &cur_d, cur_ipoint)){
-				if(cur_d<*min_d){
-					*min_d=cur_d;
-					VECCOPY(ipoint,cur_ipoint);
-					*min_face=i;
-					intersect=1;
-				}
-			}
-			if(mface->v4){
-				if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v4, v1, v3, &cur_d, cur_ipoint)){
-					if(cur_d<*min_d){
-						*min_d=cur_d;
-						VECCOPY(ipoint,cur_ipoint);
-						*min_face=i;
-						intersect=1;
-					}
-				}
-			}
+		/* particle already inside face, so report collision */
+		if(iter == 0 && d0 < 0.f && d0 > -radius) {
+			copy_v3_v3(pce->p, col->co1);
+			copy_v3_v3(pce->nor, n);
+			pce->inside = 1;
+			return 0.f;
 		}
-		else{
-			if(isect_line_tri_v3(co1, co2, v1, v2, v3, &cur_d, cur_uv)){
-				if(cur_d<*min_d){
-					*min_d=cur_d;
-					min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
-					min_w[1]= cur_uv[0];
-					min_w[2]= cur_uv[1];
-					min_w[3]= 0.0f;
-					if(mface->v4)
-						intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
-					*min_face=i;
-					intersect=1;
-				}
-			}
-			if(mface->v4){
-				if(isect_line_tri_v3(co1, co2, v1, v3, v4, &cur_d, cur_uv)){
-					if(cur_d<*min_d){
-						*min_d=cur_d;
-						min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
-						min_w[1]= 0.0f;
-						min_w[2]= cur_uv[0];
-						min_w[3]= cur_uv[1];
-						intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
-						*min_face=i;
-						intersect=1;
-					}
-				}
+		
+		dd = (t1-t0)/(d1-d0);
+
+		t0 = t1;
+		d0 = d1;
+
+		t1 -= d1*dd;
+
+		/* particle movin away from plane could also mean a strangely rotating face, so check from end */
+		if(iter == 0 && t1 < 0.f) {
+			t0 = 1.f;
+			collision_interpolate_element(pce, t0, col->f, col);
+			d0 = distance_func(col->co2, radius, pce, n);
+			t1 = 0.999f;
+			d1 = 0.f;
+
+			continue;
+		}
+		else if(iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f))
+			return -1.f;
+
+		if(d1 <= COLLISION_ZERO && d1 >= -COLLISION_ZERO) {
+			if(t1 >= -COLLISION_ZERO && t1 <= 1.f) {
+				if(distance_func == nr_signed_distance_to_plane)
+					copy_v3_v3(pce->nor, n);
+
+				CLAMP(t1, 0.f, 1.f);
+
+				return t1;
 			}
+			else
+				return -1.f;
 		}
 	}
-	return intersect;
+	return -1.0;
+}
+static int collision_sphere_to_tri(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
+{
+	ParticleCollisionElement *result = &col->pce;
+	float ct, u, v;
+
+	pce->inv_nor = -1;
+	pce->inside = 0;
+
+	ct = collision_newton_rhapson(col, radius, pce, nr_signed_distance_to_plane);
+
+	if(ct >= 0.f && ct < *t && (result->inside==0 || pce->inside==1) ) {
+		float e1[3], e2[3], p0[3];
+		float e1e1, e1e2, e1p0, e2e2, e2p0, inv;
+
+		sub_v3_v3v3(e1, pce->x1, pce->x0);
+		sub_v3_v3v3(e2, pce->x2, pce->x0);
+		/* XXX: add radius correction here? */
+		sub_v3_v3v3(p0, pce->p, pce->x0);
+
+		e1e1 = dot_v3v3(e1, e1);
+		e1e2 = dot_v3v3(e1, e2);
+		e1p0 = dot_v3v3(e1, p0);
+		e2e2 = dot_v3v3(e2, e2);
+		e2p0 = dot_v3v3(e2, p0);
+
+		inv = 1.f/(e1e1 * e2e2 - e1e2 * e1e2);
+		u = (e2e2 * e1p0 - e1e2 * e2p0) * inv;
+		v = (e1e1 * e2p0 - e1e2 * e1p0) * inv;
+
+		if(u>=0.f && u<=1.f && v>=0.f && u+v<=1.f) {
+			*result = *pce;
+
+			/* normal already calculated in pce */
+
+			result->uv[0] = u;
+			result->uv[1] = v;
+
+			*t = ct;
+			return 1;
+		}
+	}
+	return 0;
+}
+static int collision_sphere_to_edges(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t, int quad)
+{
+	ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
+	ParticleCollisionElement *result = &col->pce;
+
+	float ct;
+	int i;
+
+	for(i=0; i<3; i++) {
+		/* in case of a quad, no need to check "edge" that goes through face */
+		if((pce->x[3] && i==2) || (quad && i==0))
+			continue;
+
+		cur = edge+i;
+		cur->x[0] = pce->x[i]; cur->x[1] = pce->x[(i+1)%3];
+		cur->v[0] = pce->v[i]; cur->v[1] = pce->v[(i+1)%3];
+		cur->tot = 2;
+		cur->inside = 0;
+
+		ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_edge);
+
+		if(ct >= 0.f && ct < *t) {
+			float u, e[3], vec[3];
+
+			sub_v3_v3v3(e, cur->x1, cur->x0);
+			sub_v3_v3v3(vec, cur->p, cur->x0);
+			u = dot_v3v3(vec, e) / dot_v3v3(e, e);
+
+			if(u < 0.f || u > 1.f)
+				break;
+
+			*result = *cur;
+
+			madd_v3_v3v3fl(result->nor, vec, e, -u);
+			normalize_v3(result->nor);
+
+			result->uv[0] = u;
+
+			
+			hit = cur;
+			*t = ct;
+		}
+
+	}
+
+	return hit != NULL;
 }
+static int collision_sphere_to_verts(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
+{
+	ParticleCollisionElement vert[3], *cur = NULL, *hit = NULL;
+	ParticleCollisionElement *result = &col->pce;
+
+	float ct;
+	int i;
+
+	for(i=0; i<3; i++) {
+		/* in case of quad, only check one vert the first time */
+		if(pce->x[3] && i != 1)
+			continue;
+
+		cur = vert+i;
+		cur->x[0] = pce->x[i];
+		cur->v[0] = pce->v[i];
+		cur->tot = 1;
+		cur->inside = 0;
+
+		ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_vert);
+		
+		if(ct >= 0.f && ct < *t) {
+			*result = *cur;
+
+			sub_v3_v3v3(result->nor, cur->p, cur->x0);
+			normalize_v3(result->nor);
+
+			hit = cur;
+			*t = ct;
+		}
 
-void particle_intersect_face(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
+	}
+
+	return hit != NULL;
+}
+/* Callback for BVHTree near test */
+void BKE_psys_collision_neartest_cb(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
 {
 	ParticleCollision *col = (ParticleCollision *) userdata;
+	ParticleCollisionElement pce;
 	MFace *face = col->md->mfaces + index;
 	MVert *x = col->md->x;
 	MVert *v = col->md->current_v;
-	float vel[3], co1[3], co2[3], uv[2], ipoint[3], temp[3], t;
-	float x0[3], x1[3], x2[3], x3[3];
-	float *t0=x0, *t1=x1, *t2=x2, *t3=(face->v4 ? x3 : NULL);
-
-	/* move collision face to start of timestep */
-	madd_v3_v3v3fl(t0, x[face->v1].co, v[face->v1].co, col->cfra);
-	madd_v3_v3v3fl(t1, x[face->v2].co, v[face->v2].co, col->cfra);
-	madd_v3_v3v3fl(t2, x[face->v3].co, v[face->v3].co, col->cfra);
-	if(t3)
-		madd_v3_v3v3fl(t3, x[face->v4].co, v[face->v4].co, col->cfra);
-
-	/* calculate average velocity of face */
-	copy_v3_v3(vel, v[ face->v1 ].co);
-	add_v3_v3(vel, v[ face->v2 ].co);
-	add_v3_v3(vel, v[ face->v3 ].co);
-	mul_v3_fl(vel, 0.33334f*col->dfra);
-
-	/* substract face velocity, in other words convert to 
-	   a coordinate system where only the particle moves */
-	madd_v3_v3v3fl(co1, col->co1, vel, -col->f);
-	sub_v3_v3v3(co2, col->co2, vel);
+	float t = hit->dist/col->original_ray_length;
+	float uv[2] = {0.f,0.f};
+	int collision = 0, quad = 0;
+
+	pce.x[0] = x[face->v1].co;
+	pce.x[1] = x[face->v2].co;
+	pce.x[2] = x[face->v3].co;
+	pce.x[3] = face->v4 ? x[face->v4].co : NULL;
+
+	pce.v[0] = v[face->v1].co;
+	pce.v[1] = v[face->v2].co;
+	pce.v[2] = v[face->v3].co;
+	pce.v[3] = face->v4 ? v[face->v4].co : NULL;
+
+	pce.tot = 3;
+	pce.inside = 0;
 
 	do
 	{	
-		if(ray->radius == 0.0f) {
-			if(isect_line_tri_v3(co1, co2, t0, t1, t2, &t, uv)) {
-				if(t >= 0.0f && t < hit->dist/col->ray_len) {
-					hit->dist = col->ray_len * t;
-					hit->index = index;
-
-					/* calculate normal that's facing the particle */
-					normal_tri_v3( col->nor,t0, t1, t2);
-					VECSUB(temp, co2, co1);
-					if(dot_v3v3(col->nor, temp) > 0.0f)
-						negate_v3(col->nor);
-
-					VECCOPY(col->vel,vel);
-
-					col->hit_ob = col->ob;
-					col->hit_md = col->md;
-				}
-			}
+		collision = collision_sphere_to_tri(col, ray->radius, &pce, &t);
+		if(col->pce.inside == 0) {
+			collision += collision_sphere_to_edges(col, ray->radius, &pce, &t, quad);
+			collision += collision_sphere_to_verts(col, ray->radius, &pce, &t);
 		}
-		else {
-			if(isect_sweeping_sphere_tri_v3(co1, co2, ray->radius, t0, t1, t2, &t, ipoint)) {
-				if(t >=0.0f && t < hit->dist/col->ray_len) {
-					hit->dist = col->ray_len * t;
-					hit->index = index;
-
-					interp_v3_v3v3(temp, co1, co2, t);
-					
-					VECSUB(col->nor, temp, ipoint);
-					normalize_v3(col->nor);
 
-					VECCOPY(col->vel,vel);
+		if(collision) {
+			hit->dist = col->original_ray_length * t;
+			hit->index = index;
+				
+			collision_point_velocity(&col->pce);
 
-					col->hit_ob = col->ob;
-					col->hit_md = col->md;
-				}
-			}
+			col->hit = col->current;
 		}
 
-		t1 = t2;
-		t2 = t3;
-		t3 = NULL;
+		pce.x[1] = pce.x[2];
+		pce.x[2] = pce.x[3];
+		pce.x[3] = NULL;
 
-	} while(t2);
+		pce.v[1] = pce.v[2];
+		pce.v[2] = pce.v[3];
+		pce.v[3] = NULL;
+		quad++;
+
+	} while(pce.x[2]);
 }
-/* Particle - Mesh collision code
- * Features:
- * - point and swept sphere to mesh surface collisions
- * - moving colliders (but not yet rotating or deforming colliders)
- * - friction & damping
- * - angular momentum <-> linear momentum
- * - high accuracy by re-applying particle acceleration after collision
- * - behaves relatively well even if limit of 10 collisions per simulation step is exceeded
- * Main parts:
- * 1. check for all possible deflectors for closest intersection on particle path
- * 2. if deflection was found calculate new coordinates or kill the particle
- */
-static void deflect_particle(ParticleSimulationData *sim, int p, float dfra, float cfra){
-	Object *ground_ob = NULL;
-	ParticleSettings *part = sim->psys->part;
-	ParticleData *pa = sim->psys->particles + p;
-	ParticleCollision col;
+static int collision_detect(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, ListBase *colliders)
+{
 	ColliderCache *coll;
-	BVHTreeRayHit hit;
-	float ray_dir[3], acc[3];
-	float radius = ((part->flag & PART_SIZE_DEFL)?pa->size:0.0f), boid_z = 0.0f;
-	float timestep = psys_get_timestep(sim) * dfra;
-	float inv_timestep = 1.0f/timestep;
-	int deflections=0, max_deflections=10;
+	float ray_dir[3];
 
-	/* get acceleration (from gravity, forcefields etc. to be re-applied after collision) */
-	sub_v3_v3v3(acc, pa->state.vel, pa->prev_state.vel);
-	mul_v3_fl(acc, inv_timestep);
+	if(colliders->first == NULL)
+		return 0;
 
-	/* set values for first iteration */
-	copy_v3_v3(col.co1, pa->prev_state.co);
-	copy_v3_v3(col.co2, pa->state.co);
-	copy_v3_v3(col.ve1, pa->prev_state.vel);
-	copy_v3_v3(col.ve2, pa->state.vel);
-	col.f = 0.0f;
+	sub_v3_v3v3(ray_dir, col->co2, col->co1);
+	hit->index = -1;
+	hit->dist = col->original_ray_length = len_v3(ray_dir);
+	col->pce.inside = 0;
 
-	/* override for boids */
-	if(part->phystype == PART_PHYS_BOIDS) {
-		BoidParticle *bpa = pa->boid;
-		radius = pa->size;
-		boid_z = pa->state.co[2];
-		ground_ob = bpa->ground;
+	/* even if particle is stationary we want to check for moving colliders */
+	/* if hit.dist is zero the bvhtree_ray_cast will just ignore everything */
+	if(hit->dist == 0.0f)
+		hit->dist = col->original_ray_length = 0.000001f;
+
+	for(coll = colliders->first; coll; coll=coll->next){
+		/* for boids: don't check with current ground object */
+		if(coll->ob == col->skip)
+			continue;
+
+		/* particles should not collide with emitter at birth */
+		if(coll->ob == col->emitter && pa->time < col->cfra && pa->time >= col->old_cfra)
+			continue;
+
+		col->current = coll->ob;
+		col->md = coll->collmd;
+		col->fac1 = (col->old_cfra - coll->collmd->time_x) / (coll->collmd->time_xnew - coll->collmd->time_x);
+		col->fac2 = (col->cfra - coll->collmd->time_x) / (coll->collmd->time_xnew - coll->collmd->time_x);
+
+		if(col->md && col->md->bvhtree)
+			BLI_bvhtree_ray_cast(col->md->bvhtree, col->co1, ray_dir, col->radius, hit, BKE_psys_collision_neartest_cb, col);
 	}
 
-	/* 10 iterations to catch multiple deflections */
-	if(sim->colliders) while(deflections < max_deflections){
-		/* 1. */
+	return hit->index >= 0;
+}
+static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, int kill, int dynamic_rotation)
+{
+	ParticleCollisionElement *pce = &col->pce;
+	PartDeflect *pd = col->hit->pd;
+	float co[3];										/* point of collision */
+	float x = hit->dist/col->original_ray_length;		/* location factor of collision between this iteration */
+	float f = col->f + x * (1.0f - col->f);				/* time factor of collision between timestep */
+	float dt1 = (f - col->f) * col->total_time;			/* time since previous collision (in seconds) */
+	float dt2 = (1.0f - f) * col->total_time;			/* time left after collision (in seconds) */
+	int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */
+
+	/* calculate exact collision location */
+	interp_v3_v3v3(co, col->co1, col->co2, x);
+
+	/* particle dies in collision */
+	if(through == 0 && (kill || pd->flag & PDEFLE_KILL_PART)) {
+		pa->alive = PARS_DYING;
+		pa->dietime = col->old_cfra + (col->cfra - col->old_cfra) * f;
+
+		copy_v3_v3(pa->state.co, co);
+		interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, f);
+		interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, f);
+		interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, f);
+
+		/* particle is dead so we don't need to calculate further */
+		return 0;
+	}
+	/* figure out velocity and other data after collision */
+	else {
+		float v0[3];	/* velocity directly before collision to be modified into velocity directly after collision */
+		float v0_nor[3];/* normal component of v0 */
+		float v0_tan[3];/* tangential component of v0 */
+		float vc_tan[3];/* tangential component of collision surface velocity */
+		float v0_dot, vc_dot;
+		float damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_frand() - 0.5f);
+		float frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_frand() - 0.5f);
+		float distance, nor[3], dot;
+
+		CLAMP(damp,0.0,1.0);
+		CLAMP(frict,0.0,1.0);
+
+		/* get exact velocity right before collision */
+		madd_v3_v3v3fl(v0, col->ve1, col->acc, dt1);
+				
+		/* convert collider velocity from 1/framestep to 1/s TODO: here we assume 1 frame step for collision modifier */
+		mul_v3_fl(pce->vel, col->inv_timestep);
+
+		/* calculate tangential particle velocity */
+		v0_dot = dot_v3v3(pce->nor, v0);
+		madd_v3_v3v3fl(v0_tan, v0, pce->nor, -v0_dot);
+
+		/* calculate tangential collider velocity */
+		vc_dot = dot_v3v3(pce->nor, pce->vel);
+		madd_v3_v3v3fl(vc_tan, pce->vel, pce->nor, -vc_dot);
+
+		/* handle friction effects (tangential and angular velocity) */
+		if(frict > 0.0f) {
+			/* angular <-> linear velocity */
+			if(dynamic_rotation) {
+				float vr_tan[3], v1_tan[3], ave[3];
+					
+				/* linear velocity of particle surface */
+				cross_v3_v3v3(vr_tan, pce->nor, pa->state.ave);
+				mul_v3_fl(vr_tan, pa->size);
 
-		sub_v3_v3v3(ray_dir, col.co2, col.co1);
-		hit.index = -1;
-		hit.dist = col.ray_len = len_v3(ray_dir);
+				/* change to coordinates that move with the collision plane */
+				sub_v3_v3v3(v1_tan, v0_tan, vc_tan);
+						
+				/* The resulting velocity is a weighted average of particle cm & surface
+					* velocity. This weight (related to particle's moment of inertia) could
+					* be made a parameter for angular <-> linear conversion.
+					*/
+				madd_v3_v3fl(v1_tan, vr_tan, -0.4);
+				mul_v3_fl(v1_tan, 1.0f/1.4f); /* 1/(1+0.4) */
 
-		col.cfra = fmod(cfra-dfra, 1.0f);
-		col.dfra = dfra;
+				/* rolling friction is around 0.01 of sliding friction (could be made a parameter) */
+				mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);
 
-		/* even if particle is stationary we want to check for moving colliders */
-		/* if hit.dist is zero the bvhtree_ray_cast will just ignore everything */
-		if(hit.dist == 0.0f)
-			hit.dist = col.ray_len = 0.000001f;
+				/* surface_velocity is opposite to cm velocity */
+				mul_v3_v3fl(vr_tan, v1_tan, -1.0f);
 
-		for(coll = sim->colliders->first; coll; coll=coll->next){
-			/* for boids: don't check with current ground object */
-			if(coll->ob == ground_ob)
-				continue;
+				/* get back to global coordinates */
+				add_v3_v3(v1_tan, vc_tan);
 
-			/* particles should not collide with emitter at birth */
-			if(coll->ob == sim->ob && pa->time < cfra && pa->time >= sim->psys->cfra)
-				continue;
+				/* convert to angular velocity*/
+				cross_v3_v3v3(ave, vr_tan, pce->nor);
+				mul_v3_fl(ave, 1.0f/MAX2(pa->size, 0.001));
 
-			col.ob = coll->ob;
-			col.md = coll->collmd;
+				/* only friction will cause change in linear & angular velocity */
+				interp_v3_v3v3(pa->state.ave, pa->state.ave, ave, frict);
+				interp_v3_v3v3(v0_tan, v0_tan, v1_tan, frict);
+			}
+			else {
+				/* just basic friction (unphysical due to the friction model used in Blender) */
+				interp_v3_v3v3(v0_tan, v0_tan, vc_tan, frict);
+			}
+		}
 
-			if(col.md && col.md->bvhtree)
-				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
+		/* stickness was possibly added before, so cancel that before calculating new normal velocity */
+		/* otherwise particles go flying out of the surface because of high reversed sticky velocity */
+		if(v0_dot < 0.0f) {
+			v0_dot += pd->pdef_stickness;
+			if(v0_dot > 0.0f)
+				v0_dot = 0.0f;
 		}
 
-		/* 2. */
-		if(hit.index>=0) {
-			PartDeflect *pd = col.hit_ob->pd;
-			float co[3];							/* point of collision */
-			float x = hit.dist/col.ray_len;			/* location factor of collision between this iteration */
-			float f = col.f + x * (1.0f - col.f);	/* time factor of collision between timestep */
-			float dt1 = (f - col.f) * timestep;		/* time since previous collision (in seconds) */
-			float dt2 = (1.0f - f) * timestep;		/* time left after collision (in seconds) */
-			int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */
-
-			deflections++;
-			
-			interp_v3_v3v3(co, col.co1, col.co2, x);
-			
-			/* make sure we don't hit the current face again */
-			/* TODO: could/should this be proportional to pa->size? */
-			madd_v3_v3fl(co, col.nor, (through ? -0.0001f : 0.0001f));
+		/* damping and flipping of velocity around normal */
+		v0_dot *= 1.0f - damp;
+		vc_dot *= through ? damp : 1.0f;
 
-			/* particle dies in collision */
-			if(through == 0 && (part->flag & PART_DIE_ON_COL || pd->flag & PDEFLE_KILL_PART)) {
-				pa->alive = PARS_DYING;
-				pa->dietime = sim->psys->cfra + (cfra - sim->psys->cfra) * f;
+		/* calculate normal particle velocity */
+		/* special case for object hitting the particle from behind */
+		if(through==0 && ((vc_dot>0.0f && v0_dot>0.0f && vc_dot>v0_dot) || (vc_dot<0.0f && v0_dot<0.0f && vc_dot<v0_dot)))
+			mul_v3_v3fl(v0_nor, pce->nor, vc_dot);
+		else if(v0_dot > 0.f)
+			mul_v3_v3fl(v0_nor, pce->nor, vc_dot + (through ? -1.0f : 1.0f) * v0_dot);
+		else
+			mul_v3_v3fl(v0_nor, pce->nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);
 
-				copy_v3_v3(pa->state.co, co);
-				interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, f);
-				interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, f);
-				interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, f);
+		/* combine components together again */
+		add_v3_v3v3(v0, v0_nor, v0_tan);
 
-				/* particle is dead so we don't need to calculate further */
-				return;
+		if(col->boid) {
+			/* keep boids above ground */
+			BoidParticle *bpa = pa->boid;
+			if(bpa->data.mode == eBoidMode_OnLand || co[2] <= col->boid_z) {
+				co[2] = col->boid_z;
+				v0[2] = 0.0f;
 			}
-			/* figure out velocity and other data after collision */
-			else {
-				float v0[3];	/* velocity directly before collision to be modified into velocity directly after collision */
-				float v0_nor[3];/* normal component of v0 */
-				float v0_tan[3];/* tangential component of v0 */
-				float vc_tan[3];/* tangential component of collision surface velocity */
-				float check[3];
-				float v0_dot, vc_dot, check_dot;
-				float damp, frict;
-
-				/* get exact velocity right before collision */
-				madd_v3_v3v3fl(v0, col.ve1, acc, dt1);
-				
-				/* convert collider velocity from 1/framestep to 1/s */
-				mul_v3_fl(col.vel, inv_timestep);
-				
-				/* get damping & friction factors */
-				damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_frand() - 0.5f);
-				CLAMP(damp,0.0,1.0);
+		}
+		
+		/* re-apply acceleration to final location and velocity */
+		madd_v3_v3v3fl(pa->state.co, co, v0, dt2);
+		madd_v3_v3fl(pa->state.co, col->acc, 0.5f*dt2*dt2);
+		madd_v3_v3v3fl(pa->state.vel, v0, col->acc, dt2);
+
+		/* make sure particle stays on the right side of the surface */
+		if(!through) {
+			distance = collision_point_distance_with_normal(co, pce, -1.f, col, nor);
+			
+			if(distance < col->radius + COLLISION_MIN_DISTANCE)
+				madd_v3_v3fl(co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
 
-				frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_frand() - 0.5f);
-				CLAMP(frict,0.0,1.0);
+			dot = dot_v3v3(nor, v0);
+			if(dot < 0.f)
+				madd_v3_v3fl(v0, nor, -dot);
 
-				/* treat normal & tangent components separately */
-				v0_dot = dot_v3v3(col.nor, v0);
-				madd_v3_v3v3fl(v0_tan, v0, col.nor, -v0_dot);
+			distance = collision_point_distance_with_normal(pa->state.co, pce, 1.f, col, nor);
 
-				vc_dot = dot_v3v3(col.nor, col.vel);
-				madd_v3_v3v3fl(vc_tan, col.vel, col.nor, -vc_dot);
+			if(distance < col->radius + COLLISION_MIN_DISTANCE)
+				madd_v3_v3fl(pa->state.co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
 
-				/* handle friction effects (tangential and angular velocity) */
-				if(frict > 0.0f) {
-					/* angular <-> linear velocity */
-					if(part->flag & PART_ROT_DYN) {
-						float vr_tan[3], v1_tan[3], ave[3];
-					
-						/* linear velocity of particle surface */
-						cross_v3_v3v3(vr_tan, col.nor, pa->state.ave);
-						mul_v3_fl(vr_tan, pa->size);
+			dot = dot_v3v3(nor, pa->state.vel);
+			if(dot < 0.f)
+				madd_v3_v3fl(pa->state.vel, nor, -dot);
+		}
 
-						/* change to coordinates that move with the collision plane */
-						sub_v3_v3v3(v1_tan, v0_tan, vc_tan);
-						
-						/* The resulting velocity is a weighted average of particle cm & surface
-						 * velocity. This weight (related to particle's moment of inertia) could
-						 * be made a parameter for angular <-> linear conversion.
-						 */
-						madd_v3_v3fl(v1_tan, vr_tan, -0.4);
-						mul_v3_fl(v1_tan, 1.0f/1.4f); /* 1/(1+0.4) */
+		/* add stickness to surface */
+		madd_v3_v3fl(pa->state.vel, pce->nor, -pd->pdef_stickness);
 
-						/* rolling friction is around 0.01 of sliding friction (could be made a parameter) */
-						mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);
+		/* set coordinates for next iteration */
+		copy_v3_v3(col->co1, co);
+		copy_v3_v3(col->co2, pa->state.co);
 
-						/* surface_velocity is opposite to cm velocity */
-						mul_v3_v3fl(vr_tan, v1_tan, -1.0f);
+		copy_v3_v3(col->ve1, v0);
+		copy_v3_v3(col->ve2, pa->state.vel);
 
-						/* get back to global coordinates */
-						add_v3_v3(v1_tan, vc_tan);
+		col->f = f;
+	}
 
-						/* convert to angular velocity*/
-						cross_v3_v3v3(ave, vr_tan, col.nor);
-						mul_v3_fl(ave, 1.0f/MAX2(pa->size, 0.001));
+	col->prev = col->hit;
+	col->prev_index = hit->index;
 
-						/* only friction will cause change in linear & angular velocity */
-						interp_v3_v3v3(pa->state.ave, pa->state.ave, ave, frict);
-						interp_v3_v3v3(v0_tan, v0_tan, v1_tan, frict);
-					}
-					else {
-						/* just basic friction (unphysical due to the friction model used in Blender) */
-						interp_v3_v3v3(v0_tan, v0_tan, vc_tan, frict);
-					}
-				}
+	return 1;
+}
+static void collision_fail(ParticleData *pa, ParticleCollision *col)
+{
+	/* final chance to prevent total failure, so stick to the surface and hope for the best */
+	collision_point_on_surface(col->co1, &col->pce, 1.f, col, pa->state.co);
 
-				/* stickness was possibly added before, so cancel that before calculating new normal velocity */
-				/* otherwise particles go flying out of the surface because of high reversed sticky velocity */
-				if(v0_dot < 0.0f) {
-					v0_dot += pd->pdef_stickness;
-					if(v0_dot > 0.0f)
-						v0_dot = 0.0f;
-				}
+	copy_v3_v3(pa->state.vel, col->pce.vel);
+	mul_v3_fl(pa->state.vel, col->inv_timestep);
 
-				/* damping and flipping of velocity around normal */
-				v0_dot *= 1.0f - damp;
-				vc_dot *= through ? damp : 1.0f;
 
-				/* special case for object hitting the particle from behind */
-				if(through==0 && ((vc_dot>0.0f && v0_dot>0.0f && vc_dot>v0_dot) || (vc_dot<0.0f && v0_dot<0.0f && vc_dot<v0_dot)))
-					mul_v3_v3fl(v0_nor, col.nor, vc_dot);
-				else
-					mul_v3_v3fl(v0_nor, col.nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);
+	/* printf("max iterations\n"); */
+}
 
-				/* combine components together again */
-				add_v3_v3v3(v0, v0_nor, v0_tan);
+/* Particle - Mesh collision detection and response
+ * Features:
+ * -friction and damping
+ * -angular momentum <-> linear momentum
+ * -high accuracy by re-applying particle acceleration after collision
+ * -handles moving, rotating and deforming meshes
+ * -uses Newton-Rhapson iteration to find the collisions
+ * -handles spherical particles and (nearly) point like particles
+ */
+static void collision_check(ParticleSimulationData *sim, int p, float dfra, float cfra){
+	Object *ground_ob = NULL;
+	ParticleSettings *part = sim->psys->part;
+	ParticleData *pa = sim->psys->particles + p;
+	ParticleCollision col;
+	BVHTreeRayHit hit;
+	int collision_count=0;
 
-				/* keep boids above ground */
-				if(part->phystype == PART_PHYS_BOIDS && part->boids->options & BOID_ALLOW_LAND) {
-					BoidParticle *bpa = pa->boid;
-					if(bpa->data.mode == eBoidMode_OnLand || co[2] <= boid_z) {
-						co[2] = boid_z;
-						v0[2] = 0.0f;
-					}
-				}
-				
-				if(deflections < max_deflections) {
-					/* re-apply acceleration to final velocity and location */
-					madd_v3_v3v3fl(pa->state.vel, v0, acc, dt2);
-					madd_v3_v3v3fl(pa->state.co, co, v0, dt2);
-					madd_v3_v3fl(pa->state.co, acc, 0.5f*dt2*dt2);
+	float timestep = psys_get_timestep(sim);
 
-					/* make sure particle stays on the right side of the surface */
-					sub_v3_v3v3(check, pa->state.co, co);
-					/* (collision surface has moved during the time too) */
-					madd_v3_v3fl(check, col.vel, -dt2);
+	memset(&col, 0, sizeof(ParticleCollision));
 
-					check_dot = dot_v3v3(check, col.nor);
-					if((!through && check_dot < 0.0f) || (through && check_dot > 0.0f))
-						madd_v3_v3fl(pa->state.co, col.nor, (through ? -0.0001f : 0.0001f) - check_dot);
+	col.total_time = timestep * dfra;
+	col.inv_timestep = 1.0f/timestep;
 
-					/* Stickness to surface */
-					madd_v3_v3fl(pa->state.vel, col.nor, -pd->pdef_stickness);
+	col.cfra = cfra;
+	col.old_cfra = sim->psys->cfra;
 
-					/* set coordinates for next iteration */
-					copy_v3_v3(col.co1, co);
-					copy_v3_v3(col.co2, pa->state.co);
+	/* get acceleration (from gravity, forcefields etc. to be re-applied in collision response) */
+	sub_v3_v3v3(col.acc, pa->state.vel, pa->prev_state.vel);
+	mul_v3_fl(col.acc, 1.f/col.total_time);
 
-					copy_v3_v3(col.ve1, v0);
-					copy_v3_v3(col.ve2, pa->state.vel);
+	/* set values for first iteration */
+	copy_v3_v3(col.co1, pa->prev_state.co);
+	copy_v3_v3(col.co2, pa->state.co);
+	copy_v3_v3(col.ve1, pa->prev_state.vel);
+	copy_v3_v3(col.ve2, pa->state.vel);
+	col.f = 0.0f;
 
-					col.f = f;
-				}
-				else {
-					/* final chance to prevent failure, so stick to the surface and hope for the best */
-					madd_v3_v3v3fl(pa->state.co, co, col.vel, dt2);
-					copy_v3_v3(pa->state.vel, v0);
-				}
-			}
+	col.radius = ((part->flag & PART_SIZE_DEFL) || (part->phystype == PART_PHYS_BOIDS)) ? pa->size : COLLISION_MIN_RADIUS;
+
+	/* override for boids */
+	if(part->phystype == PART_PHYS_BOIDS && part->boids->options & BOID_ALLOW_LAND) {
+		col.boid = 1;
+		col.boid_z = pa->state.co[2];
+		col.skip = pa->boid->ground;
+	}
+
+	/* 10 iterations to catch multiple collisions */
+	while(collision_count < COLLISION_MAX_COLLISIONS){
+		if(collision_detect(pa, &col, &hit, sim->colliders)) {
+			
+			collision_count++;
+
+			if(collision_count == COLLISION_MAX_COLLISIONS)
+				collision_fail(pa, &col);
+			else if(collision_response(pa, &col, &hit, part->flag & PART_DIE_ON_COL, part->flag & PART_ROT_DYN)==0)
+				return;
 		}
 		else
 			return;
@@ -3454,7 +3679,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 	
 				/* deflection */
 				if(sim->colliders)
-					deflect_particle(sim, p, pa->state.time, cfra);
+					collision_check(sim, p, pa->state.time, cfra);
 
 				/* rotations */
 				basic_rotate(part, pa, pa->state.time, timestep);
@@ -3473,7 +3698,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 
 					/* deflection */
 					if(sim->colliders)
-						deflect_particle(sim, p, pa->state.time, cfra);
+						collision_check(sim, p, pa->state.time, cfra);
 				}
 			}
 			break;
@@ -3494,7 +3719,7 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 				sph_integrate(sim, pa, pa->state.time, gravity, springhash);
 
 				if(sim->colliders)
-					deflect_particle(sim, p, pa->state.time, cfra);
+					collision_check(sim, p, pa->state.time, cfra);
 				
 				/* SPH particles are not physical particles, just interpolation particles,  thus rotation has not a direct sense for them */	
 				basic_rotate(part, pa, pa->state.time, timestep);