Completely refactored sph fluid particles. Only the very core of the algorithm remains

the same, but big changes have happened both on the outside and on the inside.

New UI:
* The old parameters were quite true to the underlying algorithm, but were quite obscure
  from a users point of view. Now there are only a few intuitive basic parameters that
  define the basic fluid behavior.
** By default particle size is now used to determine the interaction radius, rest
   density and spring rest lengths so that it's easy to get stable simulations by simply
   emitting particles for a few frames and adjusting the particle size (easy when the
   particle size is drawn) so that the fluid appears continuous (particles are touching
   eachother).
** Stiffness - in reality most fluids are very incompressible, but this is a very hard
   problem to solve with particle based fluid simulation so some compromises have to be
   made. So the bigger the stiffness parameter is the less the fluid will compress under
   stress, but the more substeps are needed for stable simulation.
** Viscosity - how much internal friction there is in the fluid. Large viscosities also
   smooth out instabilities, so less viscous fluids again need more substeps to remain
   stable.
** Buoancy - with high buoancy low pressure areas inside the fluid start to rise against
   gravity, and high pressure areas start to come down.

* In addition to these basic parameters there are separate advanced parameters that can
  either be tweaked relative to the basic parameters (or particle size) or defined
  independently.
** Repulsion - the stiffness parameter tries to keep the fluid density constant, but this
   can lead to small clumps of particles, so the repulsion keeps the particles better
   separated.
** Stiff viscosity - the normal viscosity only applies when particles are moving closer to 
   eachother to allow free flowing fluids. Stiff viscosity also applies smoothing to
   particles that are moving away from eachother.
** Interaction radius - by default this is 4 * particle size.
** Rest density - by default this is a density that the particles have when they're packed
   densely next to eachother.
** Spring rest length - by default this is 2 * particle size.

* There are also new options for 3d view particle coloring in the display panel to show
  particle velocity and acceleration. These make it easier to see what's happening in the
  fluid simulations, but can of course be used with other particles as well.

* Viscoelastic springs have some new options too. The plasticity can now be set to much
  higher values for instant deletion of springs as the elastic limit is exeeded. In addition
  to that there is an option to only create springs for a certain number of frames when a
  particle is born. These options give new possibilities for breaking viscoelastic fluids.

New in the code:
* Most of the fluids code is now thread safe, so when particle dynamics go threaded there
  will be a nice speed boost to fluids as well.
* Fluids now use a bvh-tree instead of a kd-tree for the neighbor lookups. The bvh-tree 
  implementation makes the code quite a bit cleaner and should also give a slight speed
  boost to the simulation too.
* Previously only force fields were calculated with the different integration methods, but
  now the fluid calculations are also done using the selected integration method, so there
  are again more choices in effecting simulation accuracy and stability. This change also
  included a nice cleanup of the whole particle integration code.

As the internals are pretty stirred up old particle fluid simulations will probably not
work correctly straight away, but with some tweaking the same level of control is still
available by not using the "relative versions" of the advanced parameters (by default these
are not used when loading old files).

6 files changed, 573 insertions, 529 deletions

diff --git a/source/blender/blenkernel/intern/particle.c b/source/blender/blenkernel/intern/particle.c
index e2a23f5bc6c..fcea0ffb6cf 100644
--- a/source/blender/blenkernel/intern/particle.c
+++ b/source/blender/blenkernel/intern/particle.c
@@ -571,6 +571,7 @@ void psys_free(Object *ob, ParticleSystem * psys)
 		
 		BLI_freelistN(&psys->targets);
 
+		BLI_bvhtree_free(psys->bvhtree);
 		BLI_kdtree_free(psys->tree);
  
 		if(psys->fluid_springs)
@@ -2703,7 +2704,7 @@ static void psys_thread_create_path(ParticleThread *thread, struct ChildParticle
 			sub_v3_v3v3((child-1)->vel, child->co, (child-2)->co);
 			mul_v3_fl((child-1)->vel, 0.5);
 
-			if(ctx->ma && (part->draw & PART_DRAW_MAT_COL))
+			if(ctx->ma && (part->draw_col == PART_DRAW_COL_MAT))
 				get_strand_normal(ctx->ma, ornor, cur_length, (child-1)->vel);
 		}
 
@@ -2722,7 +2723,7 @@ static void psys_thread_create_path(ParticleThread *thread, struct ChildParticle
 			cur_length = 0.0f;
 		}
 
-		if(ctx->ma && (part->draw & PART_DRAW_MAT_COL)) {
+		if(ctx->ma && (part->draw_col == PART_DRAW_COL_MAT)) {
 			VECCOPY(child->col, &ctx->ma->r)
 			get_strand_normal(ctx->ma, ornor, cur_length, child->vel);
 		}
@@ -2907,7 +2908,7 @@ void psys_cache_paths(ParticleSimulationData *sim, float cfra)
 
 	psys->lattice = psys_get_lattice(sim);
 	ma= give_current_material(sim->ob, psys->part->omat);
-	if(ma && (psys->part->draw & PART_DRAW_MAT_COL))
+	if(ma && (psys->part->draw_col == PART_DRAW_COL_MAT))
 		VECCOPY(col, &ma->r)
 
 	if((psys->flag & PSYS_GLOBAL_HAIR)==0) {
@@ -3535,16 +3536,15 @@ static void default_particle_settings(ParticleSettings *part)
 	part->clength=1.0f;
 	part->clength_thres=0.0f;
 
-	part->draw= PART_DRAW_EMITTER|PART_DRAW_MAT_COL;
+	part->draw= PART_DRAW_EMITTER;
 	part->draw_line[0]=0.5;
 	part->path_start = 0.0f;
 	part->path_end = 1.0f;
 
 	part->keyed_loops = 1;
 
-#if 0 // XXX old animation system
-	part->ipo = NULL;
-#endif // XXX old animation system
+	part->color_vec_max = 1.f;
+	part->draw_col = PART_DRAW_COL_MAT;
 
 	part->simplify_refsize= 1920;
 	part->simplify_rate= 1.0f;
diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c
index b8535d36368..55d34f19cdc 100644
--- a/source/blender/blenkernel/intern/particle_system.c
+++ b/source/blender/blenkernel/intern/particle_system.c
@@ -1957,112 +1957,7 @@ static void set_keyed_keys(ParticleSimulationData *sim)
 
 	psys->flag |= PSYS_KEYED;
 }
-/************************************************/
-/*			Reactors							*/
-/************************************************/
-//static void push_reaction(ParticleSimulationData *sim, int pa_num, int event, ParticleKey *state)
-//{
-//	Object *rob;
-//	ParticleSystem *rpsys;
-//	ParticleSettings *rpart;
-//	ParticleData *pa;
-//	ListBase *lb=&sim->psys->effectors;
-//	ParticleEffectorCache *ec;
-//	ParticleReactEvent *re;
-//
-//	if(lb->first) for(ec = lb->first; ec; ec= ec->next){
-//		if(ec->type & PSYS_EC_REACTOR){
-//			/* all validity checks already done in add_to_effectors */
-//			rob=ec->ob;
-//			rpsys=BLI_findlink(&rob->particlesystem,ec->psys_nbr);
-//			rpart=rpsys->part;
-//			if(rpsys->part->reactevent==event){
-//				pa=sim->psys->particles+pa_num;
-//				re= MEM_callocN(sizeof(ParticleReactEvent), "react event");
-//				re->event=event;
-//				re->pa_num = pa_num;
-//				re->ob = sim->ob;
-//				re->psys = sim->psys;
-//				re->size = pa->size;
-//				copy_particle_key(&re->state,state,1);
-//
-//				switch(event){
-//					case PART_EVENT_DEATH:
-//						re->time=pa->dietime;
-//						break;
-//					case PART_EVENT_COLLIDE:
-//						re->time=state->time;
-//						break;
-//					case PART_EVENT_NEAR:
-//						re->time=state->time;
-//						break;
-//				}
-//
-//				BLI_addtail(&rpsys->reactevents, re);
-//			}
-//		}
-//	}
-//}
-//static void react_to_events(ParticleSystem *psys, int pa_num)
-//{
-//	ParticleSettings *part=psys->part;
-//	ParticleData *pa=psys->particles+pa_num;
-//	ParticleReactEvent *re=psys->reactevents.first;
-//	int birth=0;
-//	float dist=0.0f;
-//
-//	for(re=psys->reactevents.first; re; re=re->next){
-//		birth=0;
-//		if(part->from==PART_FROM_PARTICLE){
-//			if(pa->num==re->pa_num && pa->alive==PARS_UNBORN){
-//				if(re->event==PART_EVENT_NEAR){
-//					ParticleData *tpa = re->psys->particles+re->pa_num;
-//					float pa_time=tpa->time + pa->foffset*tpa->lifetime;
-//					if(re->time >= pa_time){
-//						pa->time=pa_time;
-//						pa->dietime=pa->time+pa->lifetime;
-//					}
-//				}
-//				else{
-//					pa->time=re->time;
-//					pa->dietime=pa->time+pa->lifetime;
-//				}
-//			}
-//		}
-//		else{
-//			dist=len_v3v3(pa->state.co, re->state.co);
-//			if(dist <= re->size){
-//				if(pa->alive==PARS_UNBORN){
-//					pa->time=re->time;
-//					pa->dietime=pa->time+pa->lifetime;
-//					birth=1;
-//				}
-//				if(birth || part->flag&PART_REACT_MULTIPLE){
-//					float vec[3];
-//					VECSUB(vec,pa->state.co, re->state.co);
-//					if(birth==0)
-//						mul_v3_fl(vec,(float)pow(1.0f-dist/re->size,part->reactshape));
-//					VECADDFAC(pa->state.vel,pa->state.vel,vec,part->reactfac);
-//					VECADDFAC(pa->state.vel,pa->state.vel,re->state.vel,part->partfac);
-//				}
-//				if(birth)
-//					mul_v3_fl(pa->state.vel,(float)pow(1.0f-dist/re->size,part->reactshape));
-//			}
-//		}
-//	}
-//}
-//void psys_get_reactor_target(ParticleSimulationData *sim, Object **target_ob, ParticleSystem **target_psys)
-//{
-//	Object *tob;
-//
-//	tob = sim->psys->target_ob ? sim->psys->target_ob : sim->ob;
-//	
-//	*target_psys = BLI_findlink(&tob->particlesystem, sim->psys->target_psys-1);
-//	if(*target_psys)
-//		*target_ob=tob;
-//	else
-//		*target_ob=0;
-//}
+
 /************************************************/
 /*			Point Cache							*/
 /************************************************/
@@ -2094,17 +1989,48 @@ void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra
 /************************************************/
 /*			Effectors							*/
 /************************************************/
+static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
+{
+	if(psys) {
+		PARTICLE_P;
+		int totpart = 0;
+
+		if(!psys->bvhtree || psys->bvhtree_frame != cfra) {
+			LOOP_SHOWN_PARTICLES {
+				totpart++;
+			}
+			
+			BLI_bvhtree_free(psys->bvhtree);
+			psys->bvhtree = BLI_bvhtree_new(totpart, 0.0, 4, 6);
+
+			LOOP_SHOWN_PARTICLES {
+				if(pa->alive == PARS_ALIVE) {
+					if(pa->state.time == cfra)
+						BLI_bvhtree_insert(psys->bvhtree, p, pa->prev_state.co, 1);
+					else
+						BLI_bvhtree_insert(psys->bvhtree, p, pa->state.co, 1);
+				}
+			}
+			BLI_bvhtree_balance(psys->bvhtree);
+
+			psys->bvhtree_frame = cfra;
+		}
+	}
+}
 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
 {
 	if(psys) {
 		PARTICLE_P;
+		int totpart = 0;
 
 		if(!psys->tree || psys->tree_frame != cfra) {
-			
-			BLI_kdtree_free(psys->tree);
+			LOOP_SHOWN_PARTICLES {
+				totpart++;
+			}
 
+			BLI_kdtree_free(psys->tree);
 			psys->tree = BLI_kdtree_new(psys->totpart);
-			
+
 			LOOP_SHOWN_PARTICLES {
 				if(pa->alive == PARS_ALIVE) {
 					if(pa->state.time == cfra)
@@ -2115,7 +2041,7 @@ void psys_update_particle_tree(ParticleSystem *psys, float cfra)
 			}
 			BLI_kdtree_balance(psys->tree);
 
-			psys->tree_frame = psys->cfra;
+			psys->tree_frame = cfra;
 		}
 	}
 }
@@ -2127,6 +2053,128 @@ static void psys_update_effectors(ParticleSimulationData *sim)
 	precalc_guides(sim, sim->psys->effectors);
 }
 
+static void integrate_particle(ParticleSettings *part, ParticleData *pa, float dtime, float *external_acceleration, void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse), void *forcedata)
+{
+	ParticleKey states[5];
+	float force[3],acceleration[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
+	float pa_mass= (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
+	int i, steps=1;
+	
+	copy_v3_v3(oldpos, pa->state.co);
+
+	switch(part->integrator){
+		case PART_INT_EULER:
+			steps=1;
+			break;
+		case PART_INT_MIDPOINT:
+			steps=2;
+			break;
+		case PART_INT_RK4:
+			steps=4;
+			break;
+		case PART_INT_VERLET:
+			steps=1;
+			break;
+	}
+
+	copy_particle_key(states, &pa->state, 1);
+
+	states->time = 0.f;
+
+	for(i=0; i<steps; i++){
+		zero_v3(force);
+		zero_v3(impulse);
+
+		force_func(forcedata, states+i, force, impulse);
+
+		/* force to acceleration*/
+		mul_v3_v3fl(acceleration, force, 1.0f/pa_mass);
+
+		if(external_acceleration)
+			add_v3_v3(acceleration, external_acceleration);
+		
+		/* calculate next state */
+		add_v3_v3(states[i].vel, impulse);
+
+		switch(part->integrator){
+			case PART_INT_EULER:
+				madd_v3_v3v3fl(pa->state.co, states->co, states->vel, dtime);
+				madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
+				break;
+			case PART_INT_MIDPOINT:
+				if(i==0){
+					madd_v3_v3v3fl(states[1].co, states->co, states->vel, dtime*0.5f);
+					madd_v3_v3v3fl(states[1].vel, states->vel, acceleration, dtime*0.5f);
+					states[1].time = dtime*0.5f;
+					/*fra=sim->psys->cfra+0.5f*dfra;*/
+				}
+				else{
+					madd_v3_v3v3fl(pa->state.co, states->co, states[1].vel, dtime);
+					madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
+				}
+				break;
+			case PART_INT_RK4:
+				switch(i){
+					case 0:
+						copy_v3_v3(dx[0], states->vel);
+						mul_v3_fl(dx[0], dtime);
+						copy_v3_v3(dv[0], acceleration);
+						mul_v3_fl(dv[0], dtime);
+
+						madd_v3_v3v3fl(states[1].co, states->co, dx[0], 0.5f);
+						madd_v3_v3v3fl(states[1].vel, states->vel, dv[0], 0.5f);
+						states[1].time = dtime*0.5f;
+						/*fra=sim->psys->cfra+0.5f*dfra;*/
+						break;
+					case 1:
+						madd_v3_v3v3fl(dx[1], states->vel, dv[0], 0.5f);
+						mul_v3_fl(dx[1], dtime);
+						copy_v3_v3(dv[1], acceleration);
+						mul_v3_fl(dv[1], dtime);
+
+						madd_v3_v3v3fl(states[2].co, states->co, dx[1], 0.5f);
+						madd_v3_v3v3fl(states[2].vel, states->vel, dv[1], 0.5f);
+						states[2].time = dtime*0.5f;
+						break;
+					case 2:
+						madd_v3_v3v3fl(dx[2], states->vel, dv[1], 0.5f);
+						mul_v3_fl(dx[2], dtime);
+						copy_v3_v3(dv[2], acceleration);
+						mul_v3_fl(dv[2], dtime);
+
+						add_v3_v3v3(states[3].co, states->co, dx[2]);
+						add_v3_v3v3(states[3].vel, states->vel, dv[2]);
+						states[3].time = dtime;
+						/*fra=cfra;*/
+						break;
+					case 3:
+						add_v3_v3v3(dx[3], states->vel, dv[2]);
+						mul_v3_fl(dx[3], dtime);
+						copy_v3_v3(dv[3], acceleration);
+						mul_v3_fl(dv[3], dtime);
+
+						madd_v3_v3v3fl(pa->state.co, states->co, dx[0], 1.0f/6.0f);
+						madd_v3_v3fl(pa->state.co, dx[1], 1.0f/3.0f);
+						madd_v3_v3fl(pa->state.co, dx[2], 1.0f/3.0f);
+						madd_v3_v3fl(pa->state.co, dx[3], 1.0f/6.0f);
+
+						madd_v3_v3v3fl(pa->state.vel, states->vel, dv[0], 1.0f/6.0f);
+						madd_v3_v3fl(pa->state.vel, dv[1], 1.0f/3.0f);
+						madd_v3_v3fl(pa->state.vel, dv[2], 1.0f/3.0f);
+						madd_v3_v3fl(pa->state.vel, dv[3], 1.0f/6.0f);
+				}
+				break;
+			case PART_INT_VERLET:   /* Verlet integration */
+				madd_v3_v3v3fl(pa->state.vel, pa->prev_state.vel, acceleration, dtime);
+				madd_v3_v3v3fl(pa->state.co, pa->prev_state.co, pa->state.vel, dtime);
+
+				sub_v3_v3v3(pa->state.vel, pa->state.co, oldpos);
+				mul_v3_fl(pa->state.vel, 1.0f/dtime);
+				break;
+		}
+	}
+}
+
 /*********************************************************************************************************
                     SPH fluid physics 
 
@@ -2142,7 +2190,7 @@ static void psys_update_effectors(ParticleSimulationData *sim)
 
 ***********************************************************************************************************/
 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
-static ParticleSpring *add_fluid_spring(ParticleSystem *psys, ParticleSpring *spring)
+static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
 {
 	/* Are more refs required? */
 	if(psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
@@ -2160,8 +2208,7 @@ static ParticleSpring *add_fluid_spring(ParticleSystem *psys, ParticleSpring *sp
 
 	return psys->fluid_springs + psys->tot_fluidsprings - 1;
 }
-
-static void  delete_fluid_spring(ParticleSystem *psys, int j)
+static void sph_spring_delete(ParticleSystem *psys, int j)
 {
 	if (j != psys->tot_fluidsprings - 1)
 		psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
@@ -2173,8 +2220,52 @@ static void  delete_fluid_spring(ParticleSystem *psys, int j)
 		psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs,  psys->alloc_fluidsprings * sizeof(ParticleSpring));
 	}
 }
+static void sph_springs_modify(ParticleSystem *psys, float dtime){
+	SPHFluidSettings *fluid = psys->part->fluid;
+	ParticleData *pa1, *pa2;
+	ParticleSpring *spring = psys->fluid_springs;
+	
+	float h, d, Rij[3], rij, Lij;
+	int i;
+
+	float yield_ratio = fluid->yield_ratio;
+	float plasticity = fluid->plasticity_constant;
+	/* scale things according to dtime */
+	float timefix = 25.f * dtime;
+
+	if((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
+		return;
+
+	/* Loop through the springs */
+	for(i=0; i<psys->tot_fluidsprings; i++, spring++) {
+		pa1 = psys->particles + spring->particle_index[0];
+		pa2 = psys->particles + spring->particle_index[1];
+
+		sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
+		rij = normalize_v3(Rij);
+
+		/* adjust rest length */
+		Lij = spring->rest_length;
+		d = yield_ratio * timefix * Lij;
+
+		if (rij > Lij + d) // Stretch
+			spring->rest_length += plasticity * (rij - Lij - d) * timefix;
+		else if(rij < Lij - d) // Compress
+			spring->rest_length -= plasticity * (Lij - d - rij) * timefix;
+
+		h = 4.f*pa1->size;
+
+		if(spring->rest_length > h)
+			spring->delete_flag = 1;
+	}
 
-static EdgeHash *build_fluid_springhash(ParticleSystem *psys)
+	/* Loop through springs backwaqrds - for efficient delete function */
+	for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
+		if(psys->fluid_springs[i].delete_flag)
+			sph_spring_delete(psys, i);
+	}
+}
+static EdgeHash *sph_springhash_build(ParticleSystem *psys)
 {
 	EdgeHash *springhash = NULL;
 	ParticleSpring *spring;
@@ -2187,347 +2278,275 @@ static EdgeHash *build_fluid_springhash(ParticleSystem *psys)
 
 	return springhash;
 }
-static void particle_fluidsim(ParticleSystem *psys, int own_psys, ParticleData *pa, float dtime, float mass, float *gravity, EdgeHash *springhash)
+
+typedef struct SPHNeighbor
 {
-	SPHFluidSettings *fluid = psys->part->fluid;
-	KDTreeNearest *ptn = NULL;
-	ParticleData *npa;
-	ParticleSpring *spring = NULL;
+	ParticleSystem *psys;
+	int index;
+} SPHNeighbor;
+typedef struct SPHRangeData
+{
+	SPHNeighbor neighbors[128];
+	int tot_neighbors;
 
-	float temp[3];
-	float q, q1, u, I, D, rij, d, Lij;
-	float pressure_near, pressure;
-	float p=0, pnear=0;
+	float density, near_density;
+	float h;
 
-	float omega = fluid->viscosity_omega;
-	float beta = fluid->viscosity_beta;
-	float massfactor = 1.0f/mass;
-	float spring_k = fluid->spring_k;
-	float h = fluid->radius;
-	float L = fluid->rest_length * fluid->radius;
+	ParticleSystem *npsys;
+	ParticleData *pa;
 
-	int n, neighbours = BLI_kdtree_range_search(psys->tree, h, pa->prev_state.co, NULL, &ptn);
-	int spring_index = 0, index = own_psys ? pa - psys->particles : -1;
+	float massfac;
+	int use_size;
+} SPHRangeData;
+typedef struct SPHData {
+	ParticleSystem *psys[10];
+	ParticleData *pa;
+	float mass;
+	EdgeHash *eh;
+	float *gravity;
+}SPHData;
+static void sph_density_accum_cb(void *userdata, int index, float squared_dist)
+{
+	SPHRangeData *pfr = (SPHRangeData *)userdata;
+	ParticleData *npa = pfr->npsys->particles + index;
+	float q;
 
-	/* pressure and near pressure */
-	for(n=own_psys?1:0; n<neighbours; n++) {
-		/* disregard particles at the exact same location */
-		if(ptn[n].dist < FLT_EPSILON)
-			continue;
+	if(npa == pfr->pa || squared_dist < FLT_EPSILON)
+		return;
 
-		sub_v3_v3(ptn[n].co, pa->prev_state.co);
-		mul_v3_fl(ptn[n].co, 1.f/ptn[n].dist);
-		q = ptn[n].dist/h;
+	/* Ugh! One particle has over 128 neighbors! Really shouldn't happen,
+	 * but even if it does it shouldn't do any terrible harm if all are
+	 * not taken into account - jahka
+	 */
+	if(pfr->tot_neighbors >= 128)
+		return;
+	
+	pfr->neighbors[pfr->tot_neighbors].index = index;
+	pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
+	pfr->tot_neighbors++;
 
-		if(q < 1.f) {
-			q1 = 1.f - q;
+	q = (1.f - sqrt(squared_dist)/pfr->h) * pfr->massfac;
 
-			p += q1*q1;
-			pnear += q1*q1*q1;
-		}
-	}
+	if(pfr->use_size)
+		q *= npa->size;
 
-	p *= mass;
-	pnear *= mass;
-	pressure =  fluid->stiffness_k * (p - fluid->rest_density);
-	pressure_near = fluid->stiffness_knear * pnear;
+	pfr->density += q*q;
+	pfr->near_density += q*q*q;
+}
+static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *impulse)
+{
+	SPHData *sphdata = (SPHData *)sphdata_v;
+	ParticleSystem **psys = sphdata->psys;
+	ParticleData *pa = sphdata->pa;
+	SPHFluidSettings *fluid = psys[0]->part->fluid;
+	ParticleSpring *spring = NULL;
+	SPHRangeData pfr;
+	SPHNeighbor *pfn;
+	float mass = sphdata->mass;
+	float *gravity = sphdata->gravity;
+	EdgeHash *springhash = sphdata->eh;
 
-	/* main calculations */
-	for(n=own_psys?1:0; n<neighbours; n++) {
-		/* disregard particles at the exact same location */
-		if(ptn[n].dist < FLT_EPSILON)
-			continue;
+	float q, u, rij, dv[3];
+	float pressure, near_pressure;
 
-		npa = psys->particles + ptn[n].index;
+	float visc = fluid->viscosity_omega;
+	float stiff_visc = fluid->viscosity_beta * (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
 
-		rij = ptn[n].dist;
-		q = rij/h;
-		q1 = 1.f-q;
+	float inv_mass = 1.0f/mass;
+	float spring_constant = fluid->spring_k;
+	
+	float h = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.f*pa->size : 1.f); /* 4.0 seems to be a pretty good value */
+	float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f); /* 4.77 is an experimentally determined density factor */
+	float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
 
-		/* Double Density Relaxation - Algorithm 2 (can't be thread safe!)*/
-		D =  dtime * dtime * (pressure + pressure_near*q1)*q1 * 0.5f;
-		madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
-		if(own_psys)
-			madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
+	float stiffness = fluid->stiffness_k;
+	float stiffness_near_fac = fluid->stiffness_knear * (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);
 
-		if(index < ptn[n].index) {
-			/* Viscosity - Algorithm 5 */
-			if(omega > 0.f	|| beta > 0.f) {		
-				sub_v3_v3v3(temp, pa->state.vel, npa->state.vel);
-				u = dot_v3v3(ptn[n].co, temp);
+	ParticleData *npa;
+	float vec[3];
+	float vel[3];
+	float co[3];
 
-				if (u > 0){
-					I = dtime * (q1 * (omega * u + beta * u*u)) * 0.5f;
-					madd_v3_v3fl(pa->state.vel, ptn[n].co, -I * massfactor);
+	int i, spring_index, index = pa - psys[0]->particles;
 
-					if(own_psys)
-						madd_v3_v3fl(npa->state.vel, ptn[n].co, I * massfactor);
-				}
-			}
+	pfr.tot_neighbors = 0;
+	pfr.density = pfr.near_density = 0.f;
+	pfr.h = h;
+	pfr.pa = pa;
 
-			if(spring_k > 0.f) {
-				/* Viscoelastic spring force - Algorithm 4*/
-				if (fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash){
-					spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, ptn[n].index));
+	for(i=0; i<10 && psys[i]; i++) {
+		pfr.npsys = psys[i];
+		pfr.massfac = psys[i]->part->mass*inv_mass;
+		pfr.use_size = psys[i]->part->flag & PART_SIZEMASS;
 
-					if(spring_index) {
-						spring = psys->fluid_springs + spring_index - 1;
-					}
-					else {
-						ParticleSpring temp_spring;
-						temp_spring.particle_index[0] = index;
-						temp_spring.particle_index[1] = ptn[n].index;
-						temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : L;
-						temp_spring.delete_flag = 0;
-						
-						spring = add_fluid_spring(psys, &temp_spring);
-					}
+		BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sph_density_accum_cb, &pfr);
+	}
 
-					Lij = spring->rest_length;
-					d = fluid->yield_ratio * Lij;
+	pressure =  stiffness * (pfr.density - rest_density);
+	near_pressure = stiffness_near_fac * pfr.near_density;
 
-					if (rij > Lij + d) // Stretch, 25 is just a multiplier for plasticity_constant value to counter default dtime of 1/25
-						spring->rest_length += dtime * 25.f * fluid->plasticity_constant * (rij - Lij - d);
-					else if(rij < Lij - d) // Compress
-						spring->rest_length -= dtime * 25.f * fluid->plasticity_constant * (Lij - d - rij);
-				}
-				else { /* PART_SPRING_HOOKES - Hooke's spring force */
-					/* L is a factor of radius */
-					D = 0.5 * dtime * dtime * 10.f * fluid->spring_k * (1.f - L/h) * (L - rij);
- 
-					madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
-					if(own_psys)
-						madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
-				}
-			}
-		}
-	} 
+	pfn = pfr.neighbors;
+	for(i=0; i<pfr.tot_neighbors; i++, pfn++) {
+		npa = pfn->psys->particles + pfn->index;
 
-	/* Artificial buoyancy force in negative gravity direction  */
-	if (fluid->buoyancy >= 0.f && gravity) {
-		float B = -dtime * dtime * fluid->buoyancy * (p - fluid->rest_density) * 0.5f;
-		madd_v3_v3fl(pa->state.co, gravity, -B * massfactor);
-	}
+		madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
 
-	if(ptn)
-		MEM_freeN(ptn);
-}
+		sub_v3_v3v3(vec, co, state->co);
+		rij = normalize_v3(vec);
 
-static void apply_particle_fluidsim(Object *ob, ParticleSystem *psys, ParticleData *pa, float dtime, float *gravity, EdgeHash *springhash){
-	ParticleTarget *pt;
+		q = (1.f - rij/h) * pfn->psys->part->mass * inv_mass;
 
-	particle_fluidsim(psys, 1, pa, dtime, psys->part->mass, gravity, springhash);
-	
-	/*----check other SPH systems (Multifluids) , each fluid has its own parameters---*/
-	for(pt=psys->targets.first; pt; pt=pt->next) {
-		ParticleSystem *epsys = psys_get_target_system(ob, pt);
+		if(pfn->psys->part->flag & PART_SIZEMASS)
+			q *= npa->size;
 
-		if(epsys)
-			particle_fluidsim(epsys, 0, pa, dtime, psys->part->mass, gravity, NULL);
-	}
-	/*----------------------------------------------------------------*/	 	 
-}
+		copy_v3_v3(vel, npa->prev_state.vel);
 
-static void apply_fluid_springs(ParticleSystem *psys, float timestep){
-	SPHFluidSettings *fluid = psys->part->fluid;
-	ParticleData *pa1, *pa2;
-	ParticleSpring *spring = psys->fluid_springs;
-	
-	float h = fluid->radius;
-	float massfactor = 1.0f/psys->part->mass;
-	float D, Rij[3], rij, Lij;
-	int i;
+		/* Double Density Relaxation */
+		madd_v3_v3fl(force, vec, -(pressure + near_pressure*q)*q);
 
-	if((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
-		return;
+		/* Viscosity */
+		if(visc > 0.f	|| stiff_visc > 0.f) {		
+			sub_v3_v3v3(dv, vel, state->vel);
+			u = dot_v3v3(vec, dv);
 
-	/* Loop through the springs */
-	for(i=0; i<psys->tot_fluidsprings; i++, spring++) {
-		Lij = spring->rest_length;
+			if(u < 0.f && visc > 0.f)
+				madd_v3_v3fl(force, vec, 0.5f * q * visc * u );
 
-		if (Lij > h) {
-			spring->delete_flag = 1;
+			if(u > 0.f && stiff_visc > 0.f)
+				madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u );
 		}
-		else {
-			pa1 = psys->particles + spring->particle_index[0];
-			pa2 = psys->particles + spring->particle_index[1];
 
-			sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
-			rij = normalize_v3(Rij);
+		if(spring_constant > 0.f) {
+			/* Viscoelastic spring force */
+			if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
+				spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, pfn->index));
 
-			/* Calculate displacement and apply value */
-			D =  0.5f * timestep * timestep * 10.f * fluid->spring_k * (1.f - Lij/h) * (Lij - rij);
+				if(spring_index) {
+					spring = psys[0]->fluid_springs + spring_index - 1;
 
-			madd_v3_v3fl(pa1->state.co, Rij, -D * pa1->state.time * pa1->state.time * massfactor);
-			madd_v3_v3fl(pa2->state.co, Rij, D * pa2->state.time * pa2->state.time * massfactor);
+					madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (spring->rest_length - rij));
+				}
+				else if(fluid->spring_frames == 0 || (pa->prev_state.time-pa->time) <= fluid->spring_frames){
+					ParticleSpring temp_spring;
+					temp_spring.particle_index[0] = index;
+					temp_spring.particle_index[1] = pfn->index;
+					temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length;
+					temp_spring.delete_flag = 0;
+								
+					sph_spring_add(psys[0], &temp_spring);
+				}
+			}
+			else {/* PART_SPRING_HOOKES - Hooke's spring force */
+				madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (rest_length - rij));
+			}
 		}
 	}
-
-	/* Loop through springs backwaqrds - for efficient delete function */
-	for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
-		if(psys->fluid_springs[i].delete_flag)
-			delete_fluid_spring(psys, i);
-	}
+	
+	/* Artificial buoyancy force in negative gravity direction  */
+	if (fluid->buoyancy > 0.f && gravity)
+		madd_v3_v3fl(force, gravity, fluid->buoyancy * (pfr.density-rest_density));
 }
 
-/************************************************/
-/*			Newtonian physics					*/
-/************************************************/
-/* gathers all forces that effect particles and calculates a new state for the particle */
-static void apply_particle_forces(ParticleSimulationData *sim, int p, float dfra, float cfra)
-{
+static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, float *gravity, EdgeHash *springhash){
+	ParticleTarget *pt;
+	int i;
+
 	ParticleSettings *part = sim->psys->part;
-	ParticleData *pa = sim->psys->particles + p;
-	EffectedPoint epoint;
-	ParticleKey states[5], tkey;
 	float timestep = psys_get_timestep(sim);
-	float force[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
-	float dtime=dfra*timestep, time, pa_mass=part->mass, fac /*, fra=sim->psys->cfra*/;
-	int i, steps=1;
-	ParticleTexture ptex;
+	float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
+	float dtime = dfra*psys_get_timestep(sim);
+	int steps = 1;
+	float effector_acceleration[3];
+	SPHData sphdata;
 
-	psys_get_texture(sim, pa, &ptex, PAMAP_PHYSICS, cfra);
-	
-	/* maintain angular velocity */
-	VECCOPY(pa->state.ave,pa->prev_state.ave);
-	VECCOPY(oldpos,pa->state.co);
+	sphdata.psys[0] = sim->psys;
+	for(i=1, pt=sim->psys->targets.first; i<10; i++, pt=(pt?pt->next:NULL))
+		sphdata.psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL;
 
-	if(part->flag & PART_SIZEMASS)
-		pa_mass*=pa->size;
+	sphdata.pa = pa;
+	sphdata.gravity = gravity;
+	sphdata.mass = pa_mass;
+	sphdata.eh = springhash;
 
-	switch(part->integrator){
-		case PART_INT_EULER:
-			steps=1;
-			break;
-		case PART_INT_MIDPOINT:
-			steps=2;
-			break;
-		case PART_INT_RK4:
-			steps=4;
-			break;
-		case PART_INT_VERLET:
-			steps=1;
-			break;
-	}
+	/* restore previous state and treat gravity & effectors as external acceleration*/
+	sub_v3_v3v3(effector_acceleration, pa->state.vel, pa->prev_state.vel);
+	mul_v3_fl(effector_acceleration, 1.f/dtime);
 
-	copy_particle_key(states,&pa->state,1);
+	copy_particle_key(&pa->state, &pa->prev_state, 0);
 
-	for(i=0; i<steps; i++){
-		force[0]=force[1]=force[2]=0.0;
-		impulse[0]=impulse[1]=impulse[2]=0.0;
-		/* add effectors */
-		pd_point_from_particle(sim, pa, states+i, &epoint);
-		if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
-			pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
-
-		mul_v3_fl(force, ptex.field);
-		mul_v3_fl(impulse, ptex.field);
-
-		/* calculate air-particle interaction */
-		if(part->dragfac!=0.0f){
-			fac=-part->dragfac*pa->size*pa->size*len_v3(states[i].vel);
-			VECADDFAC(force,force,states[i].vel,fac);
-		}
+	integrate_particle(part, pa, dtime, effector_acceleration, sph_force_cb, &sphdata);
+}
 
-		/* brownian force */
-		if(part->brownfac!=0.0){
-			force[0]+=(BLI_frand()-0.5f)*part->brownfac;
-			force[1]+=(BLI_frand()-0.5f)*part->brownfac;
-			force[2]+=(BLI_frand()-0.5f)*part->brownfac;
-		}
+/************************************************/
+/*			Basic physics						*/
+/************************************************/
+typedef struct EfData
+{
+	ParticleTexture ptex;
+	ParticleSimulationData *sim;
+	ParticleData *pa;
+} EfData;
+static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, float *impulse)
+{
+	EfData *efdata = (EfData *)efdata_v;
+	ParticleSimulationData *sim = efdata->sim;
+	ParticleSettings *part = sim->psys->part;
+	ParticleData *pa = efdata->pa;
+	EffectedPoint epoint;
 
-		/* force to acceleration*/
-		mul_v3_fl(force,1.0f/pa_mass);
+	/* add effectors */
+	pd_point_from_particle(efdata->sim, efdata->pa, state, &epoint);
+	if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
+		pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
 
-		/* add global acceleration (gravitation) */
-		if(psys_uses_gravity(sim)
-			/* normal gravity is too strong for hair so it's disabled by default */
-			&& (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
-			madd_v3_v3fl(force, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * ptex.gravity);
-		}
-		
-		/* calculate next state */
-		VECADD(states[i].vel,states[i].vel,impulse);
+	mul_v3_fl(force, efdata->ptex.field);
+	mul_v3_fl(impulse, efdata->ptex.field);
 
-		switch(part->integrator){
-			case PART_INT_EULER:
-				VECADDFAC(pa->state.co,states->co,states->vel,dtime);
-				VECADDFAC(pa->state.vel,states->vel,force,dtime);
-				break;
-			case PART_INT_MIDPOINT:
-				if(i==0){
-					VECADDFAC(states[1].co,states->co,states->vel,dtime*0.5f);
-					VECADDFAC(states[1].vel,states->vel,force,dtime*0.5f);
-					/*fra=sim->psys->cfra+0.5f*dfra;*/
-				}
-				else{
-					VECADDFAC(pa->state.co,states->co,states[1].vel,dtime);
-					VECADDFAC(pa->state.vel,states->vel,force,dtime);
-				}
-				break;
-			case PART_INT_RK4:
-				switch(i){
-					case 0:
-						VECCOPY(dx[0],states->vel);
-						mul_v3_fl(dx[0],dtime);
-						VECCOPY(dv[0],force);
-						mul_v3_fl(dv[0],dtime);
+	/* calculate air-particle interaction */
+	if(part->dragfac != 0.0f)
+		madd_v3_v3fl(force, state->vel, -part->dragfac * pa->size * pa->size * len_v3(state->vel));
 
-						VECADDFAC(states[1].co,states->co,dx[0],0.5f);
-						VECADDFAC(states[1].vel,states->vel,dv[0],0.5f);
-						/*fra=sim->psys->cfra+0.5f*dfra;*/
-						break;
-					case 1:
-						VECADDFAC(dx[1],states->vel,dv[0],0.5f);
-						mul_v3_fl(dx[1],dtime);
-						VECCOPY(dv[1],force);
-						mul_v3_fl(dv[1],dtime);
+	/* brownian force */
+	if(part->brownfac != 0.0){
+		force[0] += (BLI_frand()-0.5f) * part->brownfac;
+		force[1] += (BLI_frand()-0.5f) * part->brownfac;
+		force[2] += (BLI_frand()-0.5f) * part->brownfac;
+	}
+}
+/* gathers all forces that effect particles and calculates a new state for the particle */
+static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
+{
+	ParticleSettings *part = sim->psys->part;
+	ParticleData *pa = sim->psys->particles + p;
+	ParticleKey tkey;
+	float dtime=dfra*psys_get_timestep(sim), time;
+	float *gravity = NULL, gr[3];
+	EfData efdata;
 
-						VECADDFAC(states[2].co,states->co,dx[1],0.5f);
-						VECADDFAC(states[2].vel,states->vel,dv[1],0.5f);
-						break;
-					case 2:
-						VECADDFAC(dx[2],states->vel,dv[1],0.5f);
-						mul_v3_fl(dx[2],dtime);
-						VECCOPY(dv[2],force);
-						mul_v3_fl(dv[2],dtime);
+	psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);
 
-						VECADD(states[3].co,states->co,dx[2]);
-						VECADD(states[3].vel,states->vel,dv[2]);
-						/*fra=cfra;*/
-						break;
-					case 3:
-						VECADD(dx[3],states->vel,dv[2]);
-						mul_v3_fl(dx[3],dtime);
-						VECCOPY(dv[3],force);
-						mul_v3_fl(dv[3],dtime);
-
-						VECADDFAC(pa->state.co,states->co,dx[0],1.0f/6.0f);
-						VECADDFAC(pa->state.co,pa->state.co,dx[1],1.0f/3.0f);
-						VECADDFAC(pa->state.co,pa->state.co,dx[2],1.0f/3.0f);
-						VECADDFAC(pa->state.co,pa->state.co,dx[3],1.0f/6.0f);
-
-						VECADDFAC(pa->state.vel,states->vel,dv[0],1.0f/6.0f);
-						VECADDFAC(pa->state.vel,pa->state.vel,dv[1],1.0f/3.0f);
-						VECADDFAC(pa->state.vel,pa->state.vel,dv[2],1.0f/3.0f);
-						VECADDFAC(pa->state.vel,pa->state.vel,dv[3],1.0f/6.0f);
-				}
-				break;
-			case PART_INT_VERLET:   /* Verlet integration */
-				VECADDFAC(pa->state.vel,pa->state.vel,force,dtime);
-				VECADDFAC(pa->state.co,pa->state.co,pa->state.vel,dtime);
+	efdata.pa = pa;
+	efdata.sim = sim;
 
-				VECSUB(pa->state.vel,pa->state.co,oldpos);
-				mul_v3_fl(pa->state.vel,1.0f/dtime);
-				break;
-		}
+	/* add global acceleration (gravitation) */
+	if(psys_uses_gravity(sim)
+		/* normal gravity is too strong for hair so it's disabled by default */
+		&& (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
+		zero_v3(gr);
+		madd_v3_v3fl(gr, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * efdata.ptex.gravity);
+		gravity = gr;
 	}
 
+	/* maintain angular velocity */
+	copy_v3_v3(pa->state.ave, pa->prev_state.ave);
+
+	integrate_particle(part, pa, dtime, gravity, basic_force_cb, &efdata);
+
 	/* damp affects final velocity */
 	if(part->dampfac != 0.f)
-		mul_v3_fl(pa->state.vel, 1.f - part->dampfac * ptex.damp);
+		mul_v3_fl(pa->state.vel, 1.f - part->dampfac * efdata.ptex.damp);
 
-	VECCOPY(pa->state.ave, states->ave);
+	//VECCOPY(pa->state.ave, states->ave);
 
 	/* finally we do guides */
 	time=(cfra-pa->time)/pa->lifetime;
@@ -2547,7 +2566,7 @@ static void apply_particle_forces(ParticleSimulationData *sim, int p, float dfra
 		}
 	}
 }
-static void rotate_particle(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
+static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
 {
 	float rotfac, rot1[4], rot2[4]={1.0,0.0,0.0,0.0}, dtime=dfra*timestep;
 
@@ -2585,6 +2604,9 @@ static void rotate_particle(ParticleSettings *part, ParticleData *pa, float dfra
 	normalize_qt(pa->state.rot);
 }
 
+/************************************************/
+/*			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)
 {
@@ -3372,11 +3394,11 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 		case PART_PHYS_FLUID:
 		{
 			ParticleTarget *pt = psys->targets.first;
-			psys_update_particle_tree(psys, cfra);
+			psys_update_particle_bvhtree(psys, psys->cfra);
 			
 			for(; pt; pt=pt->next) {  /* Updating others systems particle tree for fluid-fluid interaction */
 				if(pt->ob)
-					psys_update_particle_tree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
+					psys_update_particle_bvhtree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), psys->cfra);
 			}
 			break;
 		}
@@ -3428,14 +3450,14 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 		{
 			LOOP_DYNAMIC_PARTICLES {
 				/* do global forces & effectors */
-				apply_particle_forces(sim, p, pa->state.time, cfra);
+				basic_integrate(sim, p, pa->state.time, cfra);
 	
 				/* deflection */
 				if(sim->colliders)
 					deflect_particle(sim, p, pa->state.time, cfra);
 
 				/* rotations */
-				rotate_particle(part, pa, pa->state.time, timestep);
+				basic_rotate(part, pa, pa->state.time, timestep);
 			}
 			break;
 		}
@@ -3458,43 +3480,28 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
 		}
 		case PART_PHYS_FLUID:
 		{
-			EdgeHash *springhash = build_fluid_springhash(psys);
+			EdgeHash *springhash = sph_springhash_build(psys);
 			float *gravity = NULL;
 
 			if(psys_uses_gravity(sim))
 				gravity = sim->scene->physics_settings.gravity;
 
-			/* do global forces & effectors */
-			LOOP_DYNAMIC_PARTICLES {
-				apply_particle_forces(sim, p, pa->state.time, cfra);
-				/* in fluids forces only effect velocity */
-				copy_v3_v3(pa->state.co, pa->prev_state.co);
-			}
-
-			/* actual fluids calculations (not threadsafe!) */
-			LOOP_DYNAMIC_PARTICLES {
-				apply_particle_fluidsim(sim->ob, psys, pa, pa->state.time*timestep, gravity, springhash);
-			}
-
-			/* Apply springs to particles */
-			apply_fluid_springs(psys, timestep);
-
-			/* apply velocity, collisions and rotation */
 			LOOP_DYNAMIC_PARTICLES {
-				/* velocity holds forces and viscosity, so apply them before collisions */
-				madd_v3_v3fl(pa->state.co, pa->state.vel, pa->state.time*timestep);
+				/* do global forces & effectors */
+				basic_integrate(sim, p, pa->state.time, cfra);
 
-				/* calculate new velocity based on new-old location */
-				sub_v3_v3v3(pa->state.vel, pa->state.co, pa->prev_state.co);
-				mul_v3_fl(pa->state.vel, 1.f/(pa->state.time*timestep));
+				/* actual fluids calculations */
+				sph_integrate(sim, pa, pa->state.time, gravity, springhash);
 
 				if(sim->colliders)
 					deflect_particle(sim, p, pa->state.time, cfra);
 				
 				/* SPH particles are not physical particles, just interpolation particles,  thus rotation has not a direct sense for them */	
-				rotate_particle(part, pa, pa->state.time, timestep);  
+				basic_rotate(part, pa, pa->state.time, timestep);  
 			}
 
+			sph_springs_modify(psys, timestep);
+
 			if(springhash) {
 				BLI_edgehash_free(springhash, NULL);
 				springhash = NULL;
@@ -3898,17 +3905,18 @@ void psys_check_boid_data(ParticleSystem *psys)
 static void fluid_default_settings(ParticleSettings *part){
 	SPHFluidSettings *fluid = part->fluid;
 
-	fluid->radius = 0.5f;
 	fluid->spring_k = 0.f;
 	fluid->plasticity_constant = 0.1f;
 	fluid->yield_ratio = 0.1f;
-	fluid->rest_length = 0.5f;
+	fluid->rest_length = 1.f;
 	fluid->viscosity_omega = 2.f;
-	fluid->viscosity_beta = 0.f;
-	fluid->stiffness_k = 0.1f;
-	fluid->stiffness_knear = 0.05f;
-	fluid->rest_density = 10.f;
+	fluid->viscosity_beta = 0.1f;
+	fluid->stiffness_k = 1.f;
+	fluid->stiffness_knear = 1.f;
+	fluid->rest_density = 1.f;
 	fluid->buoyancy = 0.f;
+	fluid->radius = 1.f;
+	fluid->flag |= SPH_FAC_REPULSION|SPH_FAC_DENSITY|SPH_FAC_RADIUS|SPH_FAC_VISCOSITY|SPH_FAC_REST_LENGTH;
 }
 
 static void psys_prepare_physics(ParticleSimulationData *sim)
diff --git a/source/blender/blenloader/intern/readfile.c b/source/blender/blenloader/intern/readfile.c
index cd8ea50b0ba..d5da98499bc 100644
--- a/source/blender/blenloader/intern/readfile.c
+++ b/source/blender/blenloader/intern/readfile.c
@@ -3330,6 +3330,7 @@ static void direct_link_particlesystems(FileData *fd, ListBase *particles)
 		}
 
 		psys->tree = NULL;
+		psys->bvhtree = NULL;
 	}
 	return;
 }
@@ -11491,7 +11492,7 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
 		bScreen *sc;
 		Brush *brush;
 		Object *ob;
-
+		ParticleSettings *part;
 		Material *mat;
 		int tex_nr, transp_tex;
 		
@@ -11539,6 +11540,12 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
 				}
 			}
 		}
+
+		/* particle draw color from material */
+		for(part = main->particle.first; part; part = part->id.next) {
+			if(part->draw & PART_DRAW_MAT_COL)
+				part->draw_col = PART_DRAW_COL_MAT;
+		}
 	}
 
 	/* put compatibility code here until next subversion bump */
diff --git a/source/blender/editors/space_view3d/drawobject.c b/source/blender/editors/space_view3d/drawobject.c
index bd3c140c647..4979c3c2a2a 100644
--- a/source/blender/editors/space_view3d/drawobject.c
+++ b/source/blender/editors/space_view3d/drawobject.c
@@ -3465,7 +3465,7 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 	Material *ma;
 	float vel[3], imat[4][4];
 	float timestep, pixsize=1.0, pa_size, r_tilt, r_length;
-	float pa_time, pa_birthtime, pa_dietime, pa_health;
+	float pa_time, pa_birthtime, pa_dietime, pa_health, intensity;
 	float cfra;
 	float ma_r=0.0f, ma_g=0.0f, ma_b=0.0f;
 	int a, totpart, totpoint=0, totve=0, drawn, draw_as, totchild=0;
@@ -3529,7 +3529,7 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 
 	if(v3d->zbuf) glDepthMask(1);
 
-	if((ma) && (part->draw&PART_DRAW_MAT_COL)) {
+	if((ma) && (part->draw_col == PART_DRAW_COL_MAT)) {
 		rgb_float_to_byte(&(ma->r), tcol);
 
 		ma_r = ma->r;
@@ -3630,6 +3630,10 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 		normalize_v3(imat[1]);
 	}
 
+	if(ELEM3(draw_as, PART_DRAW_DOT, PART_DRAW_CROSS, PART_DRAW_LINE)
+		&& part->draw_col > PART_DRAW_COL_MAT)
+		create_cdata = 1;
+
 	if(!create_cdata && pdd && pdd->cdata) {
 		MEM_freeN(pdd->cdata);
 		pdd->cdata = pdd->cd = NULL;
@@ -3672,7 +3676,7 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 		if(create_cdata && !pdd->cdata)
 			pdd->cdata = MEM_callocN(tot_vec_size, "particle_cdata");
 		if(create_ndata && !pdd->ndata)
-			pdd->ndata = MEM_callocN(tot_vec_size, "particle_vdata");
+			pdd->ndata = MEM_callocN(tot_vec_size, "particle_ndata");
 
 		if(part->draw & PART_DRAW_VEL && draw_as != PART_DRAW_LINE) {
 			if(!pdd->vedata)
@@ -3727,65 +3731,26 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 				else
 					pa_health = -1.0;
 
-#if 0 // XXX old animation system	
-				if((part->flag&PART_ABS_TIME)==0){			
-					if(ma && ma->ipo){
-						IpoCurve *icu;
-
-						/* correction for lifetime */
-						calc_ipo(ma->ipo, 100.0f*pa_time);
-
-						for(icu = ma->ipo->curve.first; icu; icu=icu->next) {
-							if(icu->adrcode == MA_COL_R)
-								ma_r = icu->curval;
-							else if(icu->adrcode == MA_COL_G)
-								ma_g = icu->curval;
-							else if(icu->adrcode == MA_COL_B)
-								ma_b = icu->curval;
-						}
-					}
-					if(part->ipo) {
-						IpoCurve *icu;
-
-						/* correction for lifetime */
-						calc_ipo(part->ipo, 100*pa_time);
+				r_tilt = 2.0f*(PSYS_FRAND(a + 21) - 0.5f);
+				r_length = PSYS_FRAND(a + 22);
 
-						for(icu = part->ipo->curve.first; icu; icu=icu->next) {
-							if(icu->adrcode == PART_SIZE)
-								pa_size = icu->curval;
-						}
+				if(part->draw_col > PART_DRAW_COL_MAT) {
+					switch(part->draw_col) {
+						case PART_DRAW_COL_VEL:
+							intensity = len_v3(pa->state.vel)/part->color_vec_max;
+							break;
+						case PART_DRAW_COL_ACC:
+							intensity = len_v3v3(pa->state.vel, pa->prev_state.vel)/((pa->state.time-pa->prev_state.time)*part->color_vec_max);
+							break;
 					}
+					CLAMP(intensity, 0.f, 1.f);
+					weight_to_rgb(intensity, &ma_r, &ma_g, &ma_b);
 				}
-#endif // XXX old animation system
-
-				r_tilt = 2.0f*(PSYS_FRAND(a + 21) - 0.5f);
-				r_length = PSYS_FRAND(a + 22);
 			}
 			else{
 				ChildParticle *cpa= &psys->child[a-totpart];
 
 				pa_time=psys_get_child_time(psys,cpa,cfra,&pa_birthtime,&pa_dietime);
-
-#if 0 // XXX old animation system
-				if((part->flag&PART_ABS_TIME)==0) {
-					if(ma && ma->ipo){
-						IpoCurve *icu;
-
-						/* correction for lifetime */
-						calc_ipo(ma->ipo, 100.0f*pa_time);
-
-						for(icu = ma->ipo->curve.first; icu; icu=icu->next) {
-							if(icu->adrcode == MA_COL_R)
-								ma_r = icu->curval;
-							else if(icu->adrcode == MA_COL_G)
-								ma_g = icu->curval;
-							else if(icu->adrcode == MA_COL_B)
-								ma_b = icu->curval;
-						}
-					}
-				}
-#endif // XXX old animation system
-
 				pa_size=psys_get_child_size(psys,cpa,cfra,NULL);
 
 				pa_health = -1.0;
@@ -3911,7 +3876,7 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 		if (1) { //ob_dt > OB_WIRE) {
 			glEnableClientState(GL_NORMAL_ARRAY);
 
-			if(part->draw&PART_DRAW_MAT_COL)
+			if(part->draw_col == PART_DRAW_COL_MAT)
 				glEnableClientState(GL_COLOR_ARRAY);
 
 			glEnable(GL_LIGHTING);
@@ -3940,7 +3905,7 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 
 				if(1) { //ob_dt > OB_WIRE) {
 					glNormalPointer(GL_FLOAT, sizeof(ParticleCacheKey), path->vel);
-					if(part->draw&PART_DRAW_MAT_COL)
+					if(part->draw_col == PART_DRAW_COL_MAT)
 						glColorPointer(3, GL_FLOAT, sizeof(ParticleCacheKey), path->col);
 				}
 
@@ -3956,7 +3921,7 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 
 			if(1) { //ob_dt > OB_WIRE) {
 				glNormalPointer(GL_FLOAT, sizeof(ParticleCacheKey), path->vel);
-				if(part->draw&PART_DRAW_MAT_COL)
+				if(part->draw_col == PART_DRAW_COL_MAT)
 					glColorPointer(3, GL_FLOAT, sizeof(ParticleCacheKey), path->col);
 			}
 
@@ -3966,7 +3931,7 @@ static void draw_new_particle_system(Scene *scene, View3D *v3d, RegionView3D *rv
 
 		/* restore & clean up */
 		if(1) { //ob_dt > OB_WIRE) {
-			if(part->draw&PART_DRAW_MAT_COL)
+			if(part->draw_col == PART_DRAW_COL_MAT)
 				glDisable(GL_COLOR_ARRAY);
 			glDisable(GL_COLOR_MATERIAL);
 		}
diff --git a/source/blender/makesdna/DNA_particle_types.h b/source/blender/makesdna/DNA_particle_types.h
index b146d49365c..09053848b28 100644
--- a/source/blender/makesdna/DNA_particle_types.h
+++ b/source/blender/makesdna/DNA_particle_types.h
@@ -123,16 +123,23 @@ typedef struct ParticleData {
 
 typedef struct SPHFluidSettings {
 	/*Particle Fluid*/
-	float spring_k, radius, rest_length, plasticity_constant, yield_ratio;
+	float radius, spring_k, rest_length;
+	float plasticity_constant, yield_ratio;
+	float plasticity_balance, yield_balance;
 	float viscosity_omega, viscosity_beta;
 	float stiffness_k, stiffness_knear, rest_density;
 	float buoyancy;
-	int flag, pad;
+	int flag, spring_frames;
 } SPHFluidSettings;
 
 /* fluid->flag */
 #define SPH_VISCOELASTIC_SPRINGS	1
 #define SPH_CURRENT_REST_LENGTH		2
+#define SPH_FAC_REPULSION			4
+#define SPH_FAC_DENSITY				8
+#define SPH_FAC_RADIUS				16
+#define SPH_FAC_VISCOSITY			32
+#define SPH_FAC_REST_LENGTH			64
 
 typedef struct ParticleSettings {
 	ID id;
@@ -143,12 +150,12 @@ typedef struct ParticleSettings {
 
 	struct EffectorWeights *effector_weights;
 
-	int flag;
+	int flag, rt;
 	short type, from, distr, texact;
 	/* physics modes */
 	short phystype, rotmode, avemode, reactevent;
 	short draw, draw_as, draw_size, childtype;
-	short ren_as, subframes;
+	short ren_as, subframes, draw_col;
 	/* number of path segments, power of 2 except */
 	short draw_step, ren_step;
 	short hair_step, keys_step;
@@ -157,12 +164,15 @@ typedef struct ParticleSettings {
 	short adapt_angle, adapt_pix;
 
 	short disp, omat, interpolation, rotfrom, integrator;
-	short kink, kink_axis, rt2;
+	short kink, kink_axis;
 
 	/* billboards */
 	short bb_align, bb_uv_split, bb_anim, bb_split_offset;
 	float bb_tilt, bb_rand_tilt, bb_offset[2];
 
+	/* draw color */
+	float color_vec_max;
+
 	/* simplification */
 	short simplify_flag, simplify_refsize;
 	float simplify_rate, simplify_transition;
@@ -249,9 +259,9 @@ typedef struct ParticleSystem{				/* note, make sure all (runtime) are NULL's in
 	char name[32];							/* particle system name */
 	
 	float imat[4][4];	/* used for duplicators */
-	float cfra, tree_frame;
+	float cfra, tree_frame, bvhtree_frame;
 	int seed, child_seed;
-	int flag, totpart, totunexist, totchild, totcached, totchildcache, rt;
+	int flag, totpart, totunexist, totchild, totcached, totchildcache;
 	short recalc, target_psys, totkeyed, bakespace;
 
 	char bb_uvname[3][32];					/* billboard uv name */
@@ -271,7 +281,8 @@ typedef struct ParticleSystem{				/* note, make sure all (runtime) are NULL's in
 	ParticleSpring *fluid_springs;
 	int tot_fluidsprings, alloc_fluidsprings;
 
-	struct KDTree *tree;					/* used for interactions with self and other systems */
+	struct KDTree *tree;								/* used for interactions with self and other systems */
+	struct BVHTree *bvhtree;								/* used for interactions with self and other systems */
 
 	struct ParticleDrawData *pdd;
 
@@ -375,10 +386,17 @@ typedef struct ParticleSystem{				/* note, make sure all (runtime) are NULL's in
 #define PART_DRAW_RAND_GR	1024
 #define PART_DRAW_REN_ADAPT	2048
 #define PART_DRAW_VEL_LENGTH	(1<<12)
-#define PART_DRAW_MAT_COL		(1<<13)
+#define PART_DRAW_MAT_COL		(1<<13) /* deprecated, but used in do_versions */
 #define PART_DRAW_WHOLE_GR		(1<<14)
 #define PART_DRAW_REN_STRAND	(1<<15)
 
+/* part->draw_col */
+#define PART_DRAW_COL_NONE		0
+#define PART_DRAW_COL_MAT		1
+#define PART_DRAW_COL_VEL		2
+#define PART_DRAW_COL_ACC		3
+
+
 /* part->simplify_flag */
 #define PART_SIMPLIFY_ENABLE	1
 #define PART_SIMPLIFY_VIEWPORT	2
diff --git a/source/blender/makesrna/intern/rna_particle.c b/source/blender/makesrna/intern/rna_particle.c
index 57d6e72a102..069f674bf4a 100644
--- a/source/blender/makesrna/intern/rna_particle.c
+++ b/source/blender/makesrna/intern/rna_particle.c
@@ -1073,10 +1073,9 @@ static void rna_def_fluid_settings(BlenderRNA *brna)
 	RNA_def_property_ui_text(prop, "Interaction Radius", "Fluid interaction radius");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
 
-	/* Hidden in ui to give a little easier user experience. */
  	prop= RNA_def_property(srna, "rest_length", PROP_FLOAT, PROP_NONE);
- 	RNA_def_property_range(prop, 0.0f, 1.0f);
-	RNA_def_property_ui_text(prop, "Rest Length", "Spring rest length (factor of interaction radius)");
+ 	RNA_def_property_range(prop, 0.0f, 2.0f);
+	RNA_def_property_ui_text(prop, "Rest Length", "Spring rest length (factor of particle radius)");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
 
 	prop= RNA_def_property(srna, "use_viscoelastic_springs", PROP_BOOLEAN, PROP_NONE);
@@ -1086,12 +1085,12 @@ static void rna_def_fluid_settings(BlenderRNA *brna)
 
 	prop= RNA_def_property(srna, "use_initial_rest_length", PROP_BOOLEAN, PROP_NONE);
 	RNA_def_property_boolean_sdna(prop, NULL, "flag", SPH_CURRENT_REST_LENGTH);
-	RNA_def_property_ui_text(prop, "Initial Rest Length", "Use the initial length as spring rest length instead of interaction radius/2");
+	RNA_def_property_ui_text(prop, "Initial Rest Length", "Use the initial length as spring rest length instead of 2 * particle size");
 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
 
  	prop= RNA_def_property(srna, "plasticity", PROP_FLOAT, PROP_NONE);
  	RNA_def_property_float_sdna(prop, NULL, "plasticity_constant");
- 	RNA_def_property_range(prop, 0.0f, 1.0f);
+ 	RNA_def_property_range(prop, 0.0f, 100.0f);
  	RNA_def_property_ui_text(prop, "Plasticity", "How much the spring rest length can change after the elastic limit is crossed");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
 
@@ -1100,6 +1099,11 @@ static void rna_def_fluid_settings(BlenderRNA *brna)
  	RNA_def_property_range(prop, 0.0f, 1.0f);
  	RNA_def_property_ui_text(prop, "Elastic Limit", "How much the spring has to be stretched/compressed in order to change it's rest length");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+ 	prop= RNA_def_property(srna, "spring_frames", PROP_INT, PROP_NONE);
+ 	RNA_def_property_range(prop, 0.0f, 100.0f);
+ 	RNA_def_property_ui_text(prop, "Spring Frames", "Create springs for this number of frames since particles birth (0 is always)");
+ 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
  
  	/* Viscosity */
  	prop= RNA_def_property(srna, "linear_viscosity", PROP_FLOAT, PROP_NONE);
@@ -1109,42 +1113,69 @@ static void rna_def_fluid_settings(BlenderRNA *brna)
  	RNA_def_property_ui_text(prop, "Viscosity", "Linear viscosity");
 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
   
- 	prop= RNA_def_property(srna, "square_viscosity", PROP_FLOAT, PROP_NONE);
+ 	prop= RNA_def_property(srna, "stiff_viscosity", PROP_FLOAT, PROP_NONE);
 	RNA_def_property_float_sdna(prop, NULL, "viscosity_beta");
 	RNA_def_property_range(prop, 0.0f, 100.0f);
-	RNA_def_property_ui_range(prop, 0.0f, 10.0f, 1, 3);
- 	RNA_def_property_ui_text(prop, "Square viscosity", "Square viscosity");
+	RNA_def_property_ui_range(prop, 0.0f, 2.0f, 1, 3);
+ 	RNA_def_property_ui_text(prop, "Stiff viscosity", "Creates viscosity for expanding fluid)");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
 
 	/* Double density relaxation */
- 	prop= RNA_def_property(srna, "density_force", PROP_FLOAT, PROP_NONE);
+ 	prop= RNA_def_property(srna, "stiffness", PROP_FLOAT, PROP_NONE);
  	RNA_def_property_float_sdna(prop, NULL, "stiffness_k");
  	RNA_def_property_range(prop, 0.0f, 100.0f);
 	RNA_def_property_ui_range(prop, 0.0f, 10.0f, 1, 3);
- 	RNA_def_property_ui_text(prop, "Density Force", "How strongly the fluid tends to rest density");
+ 	RNA_def_property_ui_text(prop, "Stiffness", "How incompressible the fluid is");
 	RNA_def_property_update(prop, 0, "rna_Particle_reset");
   
- 	prop= RNA_def_property(srna, "repulsion_force", PROP_FLOAT, PROP_NONE);
+ 	prop= RNA_def_property(srna, "repulsion", PROP_FLOAT, PROP_NONE);
  	RNA_def_property_float_sdna(prop, NULL, "stiffness_knear");
 	RNA_def_property_range(prop, 0.0f, 100.0f);
-	RNA_def_property_ui_range(prop, 0.0f, 10.0f, 1, 3);
- 	RNA_def_property_ui_text(prop, "Repulsion", "How strongly the fluid tries to keep from clustering");
+	RNA_def_property_ui_range(prop, 0.0f, 2.0f, 1, 3);
+ 	RNA_def_property_ui_text(prop, "Repulsion Factor", "How strongly the fluid tries to keep from clustering (factor of stiffness)");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
   
  	prop= RNA_def_property(srna, "rest_density", PROP_FLOAT, PROP_NONE);
 	RNA_def_property_float_sdna(prop, NULL, "rest_density");
- 	RNA_def_property_range(prop, 0.0f, 1000.0f);
-	RNA_def_property_ui_range(prop, 0.0f, 100.0f, 1, 3);
- 	RNA_def_property_ui_text(prop, "Rest Density", "Rest density of the fluid");
+ 	RNA_def_property_range(prop, 0.0f, 100.0f);
+	RNA_def_property_ui_range(prop, 0.0f, 2.0f, 1, 3);
+ 	RNA_def_property_ui_text(prop, "Rest Density", "Fluid rest density");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
 	
 	/* Buoyancy */
 	prop= RNA_def_property(srna, "buoyancy", PROP_FLOAT, PROP_NONE);
  	RNA_def_property_float_sdna(prop, NULL, "buoyancy");
- 	RNA_def_property_range(prop, 0.0f, 1.0f);
+ 	RNA_def_property_range(prop, 0.0f, 10.0f);
+	RNA_def_property_ui_range(prop, 0.0f, 1.0f, 1, 3);
  	RNA_def_property_ui_text(prop, "Buoyancy", "Artificial buoyancy force in negative gravity direction based on pressure differences inside the fluid");
  	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+	/* Factor flags */
 	
+	prop= RNA_def_property(srna, "factor_repulsion", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", SPH_FAC_REPULSION);
+	RNA_def_property_ui_text(prop, "Factor Repulsion", "Repulsion is a factor of stiffness");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+	prop= RNA_def_property(srna, "factor_density", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", SPH_FAC_DENSITY);
+	RNA_def_property_ui_text(prop, "Factor Density", "Density is calculated as a factor of default density (depends on particle size)");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+	prop= RNA_def_property(srna, "factor_radius", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", SPH_FAC_RADIUS);
+	RNA_def_property_ui_text(prop, "Factor Radius", "Interaction radius is a factor of 4 * particle size");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+	prop= RNA_def_property(srna, "factor_stiff_viscosity", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", SPH_FAC_VISCOSITY);
+	RNA_def_property_ui_text(prop, "Factor Stiff Viscosity", "Stiff viscosity is a factor of normal viscosity");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
+	prop= RNA_def_property(srna, "factor_rest_length", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "flag", SPH_FAC_REST_LENGTH);
+	RNA_def_property_ui_text(prop, "Factor Rest Length", "Spring rest length is a factor of 2 * particle size");
+	RNA_def_property_update(prop, 0, "rna_Particle_reset");
 }
 
 static void rna_def_particle_settings_mtex(BlenderRNA *brna)
@@ -1483,6 +1514,14 @@ static void rna_def_particle_settings(BlenderRNA *brna)
 		{0, NULL, 0, NULL, NULL}
 	};
 
+	static EnumPropertyItem draw_col_items[] = {
+		{PART_DRAW_COL_NONE, "NONE", 0, "None", ""},
+		{PART_DRAW_COL_MAT, "MATERIAL", 0, "Material", ""},
+		{PART_DRAW_COL_VEL, "VELOCITY", 0, "Velocity", ""},
+		{PART_DRAW_COL_ACC, "ACCELERATION", 0, "Acceleration", ""},
+		{0, NULL, 0, NULL, NULL}
+	};
+
 	srna= RNA_def_struct(brna, "ParticleSettings", "ID");
 	RNA_def_struct_ui_text(srna, "Particle Settings", "Particle settings, reusable by multiple particle systems");
 	RNA_def_struct_ui_icon(srna, ICON_PARTICLE_DATA);
@@ -1725,11 +1764,6 @@ static void rna_def_particle_settings(BlenderRNA *brna)
 	RNA_def_property_ui_text(prop, "Speed", "Multiply line length by particle speed");
 	RNA_def_property_update(prop, 0, "rna_Particle_redo");
 
-	prop= RNA_def_property(srna, "show_material_color", PROP_BOOLEAN, PROP_NONE);
-	RNA_def_property_boolean_sdna(prop, NULL, "draw", PART_DRAW_MAT_COL);
-	RNA_def_property_ui_text(prop, "Material Color", "Draw particles using material's diffuse color");
-	RNA_def_property_update(prop, 0, "rna_Particle_redo");
-
 	prop= RNA_def_property(srna, "use_whole_group", PROP_BOOLEAN, PROP_NONE);
 	RNA_def_property_boolean_sdna(prop, NULL, "draw", PART_DRAW_WHOLE_GR);
 	RNA_def_property_ui_text(prop, "Whole Group", "Use whole group at once");
@@ -1754,6 +1788,12 @@ static void rna_def_particle_settings(BlenderRNA *brna)
 	RNA_def_property_ui_text(prop, "Particle Rendering", "How particles are rendered");
 	RNA_def_property_update(prop, 0, "rna_Particle_redo");
 
+	prop= RNA_def_property(srna, "draw_color", PROP_ENUM, PROP_NONE);
+	RNA_def_property_enum_sdna(prop, NULL, "draw_col");
+	RNA_def_property_enum_items(prop, draw_col_items);
+	RNA_def_property_ui_text(prop, "Draw Color", "Draw additional particle data as a color");
+	RNA_def_property_update(prop, 0, "rna_Particle_redo");
+
 	prop= RNA_def_property(srna, "draw_size", PROP_INT, PROP_NONE);
 	RNA_def_property_range(prop, 0, 1000);
 	RNA_def_property_ui_range(prop, 0, 100, 1, 0);
@@ -1865,6 +1905,12 @@ static void rna_def_particle_settings(BlenderRNA *brna)
 	RNA_def_property_ui_text(prop, "Tilt", "Tilt of the billboards");
 	RNA_def_property_update(prop, 0, "rna_Particle_redo");
 
+	prop= RNA_def_property(srna, "color_maximum", PROP_FLOAT, PROP_NONE);
+	RNA_def_property_float_sdna(prop, NULL, "color_vec_max");
+	RNA_def_property_range(prop, 0.01f, 100.0f);
+	RNA_def_property_ui_text(prop, "Color Maximum", "Maximum length of the particle color vector");
+	RNA_def_property_update(prop, 0, "rna_Particle_redo");
+
 	prop= RNA_def_property(srna, "billboard_tilt_random", PROP_FLOAT, PROP_NONE);
 	RNA_def_property_float_sdna(prop, NULL, "bb_rand_tilt");
 	RNA_def_property_range(prop, 0.0f, 1.0f);