blender-archive/source/blender/blenkernel/intern/particle_system.c

/*
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 * The Original Code is Copyright (C) 2007 by Janne Karhu.
 * All rights reserved.
 * Adaptive time step
 * Classical SPH
 * Copyright 2011-2012 AutoCRC
 */

/** \file
 * \ingroup bke
 */


#include <stddef.h>

#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "MEM_guardedalloc.h"

#include "DNA_anim_types.h"
#include "DNA_boid_types.h"
#include "DNA_particle_types.h"
#include "DNA_mesh_types.h"
#include "DNA_meshdata_types.h"
#include "DNA_modifier_types.h"
#include "DNA_object_force_types.h"
#include "DNA_object_types.h"
#include "DNA_curve_types.h"
#include "DNA_scene_types.h"
#include "DNA_texture_types.h"
#include "DNA_listBase.h"

#include "BLI_utildefines.h"
#include "BLI_edgehash.h"
#include "BLI_rand.h"
#include "BLI_math.h"
#include "BLI_blenlib.h"
#include "BLI_kdtree.h"
#include "BLI_kdopbvh.h"
#include "BLI_task.h"
#include "BLI_threads.h"
#include "BLI_linklist.h"
#include "BLI_string_utils.h"

#include "BKE_animsys.h"
#include "BKE_boids.h"
#include "BKE_collision.h"
#include "BKE_colortools.h"
#include "BKE_effect.h"
#include "BKE_library.h"
#include "BKE_library_query.h"
#include "BKE_particle.h"

#include "BKE_collection.h"
#include "BKE_object.h"
#include "BKE_material.h"
#include "BKE_cloth.h"
#include "BKE_lattice.h"
#include "BKE_pointcache.h"
#include "BKE_mesh.h"
#include "BKE_modifier.h"
#include "BKE_scene.h"
#include "BKE_bvhutils.h"

#include "DEG_depsgraph.h"
#include "DEG_depsgraph_physics.h"
#include "DEG_depsgraph_query.h"

#include "PIL_time.h"

#include "RE_shader_ext.h"

/* fluid sim particle import */
#ifdef WITH_MOD_FLUID
#include "DNA_object_fluidsim_types.h"
#include "LBM_fluidsim.h"
#include <zlib.h>
#include <string.h>

#endif // WITH_MOD_FLUID

static ThreadRWMutex psys_bvhtree_rwlock = BLI_RWLOCK_INITIALIZER;

/************************************************/
/*			Reacting to system events			*/
/************************************************/

static int particles_are_dynamic(ParticleSystem *psys)
{
	if (psys->pointcache->flag & PTCACHE_BAKED)
		return 0;

	if (psys->part->type == PART_HAIR)
		return psys->flag & PSYS_HAIR_DYNAMICS;
	else
		return ELEM(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
}

float psys_get_current_display_percentage(ParticleSystem *psys, const bool use_render_params)
{
	ParticleSettings *part=psys->part;

	if ((use_render_params && !particles_are_dynamic(psys)) ||  /* non-dynamic particles can be rendered fully */
	    (part->child_nbr && part->childtype)  ||    /* display percentage applies to children */
	    (psys->pointcache->flag & PTCACHE_BAKING))  /* baking is always done with full amount */
	{
		return 1.0f;
	}

	return psys->part->disp/100.0f;
}

static int tot_particles(ParticleSystem *psys, PTCacheID *pid)
{
	if (pid && psys->pointcache->flag & PTCACHE_EXTERNAL)
		return pid->cache->totpoint;
	else if (psys->part->distr == PART_DISTR_GRID && psys->part->from != PART_FROM_VERT)
		return psys->part->grid_res * psys->part->grid_res * psys->part->grid_res - psys->totunexist;
	else
		return psys->part->totpart - psys->totunexist;
}

void psys_reset(ParticleSystem *psys, int mode)
{
	PARTICLE_P;

	if (ELEM(mode, PSYS_RESET_ALL, PSYS_RESET_DEPSGRAPH)) {
		if (mode == PSYS_RESET_ALL || !(psys->flag & PSYS_EDITED)) {
			/* don't free if not absolutely necessary */
			if (psys->totpart != tot_particles(psys, NULL)) {
				psys_free_particles(psys);
				psys->totpart= 0;
			}

			psys->totkeyed= 0;
			psys->flag &= ~(PSYS_HAIR_DONE|PSYS_KEYED);

			if (psys->edit && psys->free_edit) {
				psys->free_edit(psys->edit);
				psys->edit = NULL;
				psys->free_edit = NULL;
			}
		}
	}
	else if (mode == PSYS_RESET_CACHE_MISS) {
		/* set all particles to be skipped */
		LOOP_PARTICLES {
			pa->flag |= PARS_NO_DISP;
		}
	}

	/* reset children */
	if (psys->child) {
		MEM_freeN(psys->child);
		psys->child= NULL;
	}

	psys->totchild= 0;

	/* reset path cache */
	psys_free_path_cache(psys, psys->edit);

	/* reset point cache */
	BKE_ptcache_invalidate(psys->pointcache);

	if (psys->fluid_springs) {
		MEM_freeN(psys->fluid_springs);
		psys->fluid_springs = NULL;
	}

	psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
}

void psys_unique_name(Object *object, ParticleSystem *psys, const char *defname)
{
	BLI_uniquename(&object->particlesystem, psys, defname, '.',
	    offsetof(ParticleSystem, name), sizeof(psys->name));
}

static void realloc_particles(ParticleSimulationData *sim, int new_totpart)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	ParticleData *newpars = NULL;
	BoidParticle *newboids = NULL;
	PARTICLE_P;
	int totpart, totsaved = 0;

	if (new_totpart<0) {
		if ((part->distr == PART_DISTR_GRID) && (part->from != PART_FROM_VERT)) {
			totpart= part->grid_res;
			totpart*=totpart*totpart;
		}
		else
			totpart=part->totpart;
	}
	else
		totpart=new_totpart;

	if (totpart != psys->totpart) {
		if (psys->edit && psys->free_edit) {
			psys->free_edit(psys->edit);
			psys->edit = NULL;
			psys->free_edit = NULL;
		}

		if (totpart) {
			newpars= MEM_callocN(totpart*sizeof(ParticleData), "particles");
			if (newpars == NULL)
				return;

			if (psys->part->phystype == PART_PHYS_BOIDS) {
				newboids= MEM_callocN(totpart*sizeof(BoidParticle), "boid particles");

				if (newboids == NULL) {
					/* allocation error! */
					if (newpars)
						MEM_freeN(newpars);
					return;
				}
			}
		}

		if (psys->particles) {
			totsaved=MIN2(psys->totpart,totpart);
			/*save old pars*/
			if (totsaved) {
				memcpy(newpars,psys->particles,totsaved*sizeof(ParticleData));

				if (psys->particles->boid)
					memcpy(newboids, psys->particles->boid, totsaved*sizeof(BoidParticle));
			}

			if (psys->particles->keys)
				MEM_freeN(psys->particles->keys);

			if (psys->particles->boid)
				MEM_freeN(psys->particles->boid);

			for (p=0, pa=newpars; p<totsaved; p++, pa++) {
				if (pa->keys) {
					pa->keys= NULL;
					pa->totkey= 0;
				}
			}

			for (p=totsaved, pa=psys->particles+totsaved; p<psys->totpart; p++, pa++)
				if (pa->hair) MEM_freeN(pa->hair);

			MEM_freeN(psys->particles);
			psys_free_pdd(psys);
		}

		psys->particles=newpars;
		psys->totpart=totpart;

		if (newboids) {
			LOOP_PARTICLES {
				pa->boid = newboids++;
			}
		}
	}

	if (psys->child) {
		MEM_freeN(psys->child);
		psys->child=NULL;
		psys->totchild=0;
	}
}

int psys_get_child_number(Scene *scene, ParticleSystem *psys, const bool use_render_params)
{
	int nbr;

	if (!psys->part->childtype)
		return 0;

	if (use_render_params)
		nbr= psys->part->ren_child_nbr;
	else
		nbr= psys->part->child_nbr;

	return get_render_child_particle_number(&scene->r, nbr, use_render_params);
}

int psys_get_tot_child(Scene *scene, ParticleSystem *psys, const bool use_render_params)
{
	return psys->totpart*psys_get_child_number(scene, psys, use_render_params);
}

/************************************************/
/*			Distribution						*/
/************************************************/

void psys_calc_dmcache(Object *ob, Mesh *mesh_final, Mesh *mesh_original, ParticleSystem *psys)
{
	/* use for building derived mesh mapping info:
	 *
	 * node: the allocated links - total derived mesh element count
	 * nodearray: the array of nodes aligned with the base mesh's elements, so
	 *            each original elements can reference its derived elements
	 */
	Mesh *me= (Mesh*)ob->data;
	bool use_modifier_stack= psys->part->use_modifier_stack;
	PARTICLE_P;

	/* CACHE LOCATIONS */
	if (!mesh_final->runtime.deformed_only) {
		/* Will use later to speed up subsurf/evaluated mesh. */
		LinkNode *node, *nodedmelem, **nodearray;
		int totdmelem, totelem, i, *origindex, *origindex_poly = NULL;

		if (psys->part->from == PART_FROM_VERT) {
			totdmelem = mesh_final->totvert;

			if (use_modifier_stack) {
				totelem= totdmelem;
				origindex= NULL;
			}
			else {
				totelem= me->totvert;
				origindex = CustomData_get_layer(&mesh_final->vdata, CD_ORIGINDEX);
			}
		}
		else { /* FROM_FACE/FROM_VOLUME */
			totdmelem= mesh_final->totface;

			if (use_modifier_stack) {
				totelem= totdmelem;
				origindex= NULL;
				origindex_poly= NULL;
			}
			else {
				totelem = mesh_original->totface;
				origindex = CustomData_get_layer(&mesh_final->fdata, CD_ORIGINDEX);

				/* for face lookups we need the poly origindex too */
				origindex_poly = CustomData_get_layer(&mesh_final->pdata, CD_ORIGINDEX);
				if (origindex_poly == NULL) {
					origindex= NULL;
				}
			}
		}

		nodedmelem= MEM_callocN(sizeof(LinkNode)*totdmelem, "psys node elems");
		nodearray= MEM_callocN(sizeof(LinkNode *)*totelem, "psys node array");

		for (i=0, node=nodedmelem; i<totdmelem; i++, node++) {
			int origindex_final;
			node->link = POINTER_FROM_INT(i);

			/* may be vertex or face origindex */
			if (use_modifier_stack) {
				origindex_final = i;
			}
			else {
				origindex_final = origindex ? origindex[i] : ORIGINDEX_NONE;

				/* if we have a poly source, do an index lookup */
				if (origindex_poly && origindex_final != ORIGINDEX_NONE) {
					origindex_final = origindex_poly[origindex_final];
				}
			}

			if (origindex_final != ORIGINDEX_NONE && origindex_final < totelem) {
				if (nodearray[origindex_final]) {
					/* prepend */
					node->next = nodearray[origindex_final];
					nodearray[origindex_final] = node;
				}
				else {
					nodearray[origindex_final] = node;
				}
			}
		}

		/* cache the verts/faces! */
		LOOP_PARTICLES {
			if (pa->num < 0) {
				pa->num_dmcache = DMCACHE_NOTFOUND;
				continue;
			}

			if (use_modifier_stack) {
				if (pa->num < totelem)
					pa->num_dmcache = DMCACHE_ISCHILD;
				else
					pa->num_dmcache = DMCACHE_NOTFOUND;
			}
			else {
				if (psys->part->from == PART_FROM_VERT) {
					if (pa->num < totelem && nodearray[pa->num])
						pa->num_dmcache= POINTER_AS_INT(nodearray[pa->num]->link);
					else
						pa->num_dmcache = DMCACHE_NOTFOUND;
				}
				else { /* FROM_FACE/FROM_VOLUME */
					pa->num_dmcache = psys_particle_dm_face_lookup(mesh_final, mesh_original, pa->num, pa->fuv, nodearray);
				}
			}
		}

		MEM_freeN(nodearray);
		MEM_freeN(nodedmelem);
	}
	else {
		/* TODO PARTICLE, make the following line unnecessary, each function
		 * should know to use the num or num_dmcache, set the num_dmcache to
		 * an invalid value, just in case */

		LOOP_PARTICLES {
			pa->num_dmcache = DMCACHE_NOTFOUND;
		}
	}
}

/* threaded child particle distribution and path caching */
void psys_thread_context_init(ParticleThreadContext *ctx, ParticleSimulationData *sim)
{
	memset(ctx, 0, sizeof(ParticleThreadContext));
	ctx->sim = *sim;
	ctx->mesh = ctx->sim.psmd->mesh_final;
	ctx->ma = give_current_material(sim->ob, sim->psys->part->omat);
}

#define MAX_PARTICLES_PER_TASK 256 /* XXX arbitrary - maybe use at least number of points instead for better balancing? */

BLI_INLINE int ceil_ii(int a, int b)
{
	return (a + b - 1) / b;
}

void psys_tasks_create(ParticleThreadContext *ctx, int startpart, int endpart, ParticleTask **r_tasks, int *r_numtasks)
{
	ParticleTask *tasks;
	int numtasks = ceil_ii((endpart - startpart), MAX_PARTICLES_PER_TASK);
	float particles_per_task = (float)(endpart - startpart) / (float)numtasks, p, pnext;
	int i;

	tasks = MEM_callocN(sizeof(ParticleTask) * numtasks, "ParticleThread");
	*r_numtasks = numtasks;
	*r_tasks = tasks;

	p = (float)startpart;
	for (i = 0; i < numtasks; i++, p = pnext) {
		pnext = p + particles_per_task;

		tasks[i].ctx = ctx;
		tasks[i].begin = (int)p;
		tasks[i].end = min_ii((int)pnext, endpart);
	}
}

void psys_tasks_free(ParticleTask *tasks, int numtasks)
{
	int i;

	/* threads */
	for (i = 0; i < numtasks; ++i) {
		if (tasks[i].rng)
			BLI_rng_free(tasks[i].rng);
		if (tasks[i].rng_path)
			BLI_rng_free(tasks[i].rng_path);
	}

	MEM_freeN(tasks);
}

void psys_thread_context_free(ParticleThreadContext *ctx)
{
	/* path caching */
	if (ctx->vg_length)
		MEM_freeN(ctx->vg_length);
	if (ctx->vg_clump)
		MEM_freeN(ctx->vg_clump);
	if (ctx->vg_kink)
		MEM_freeN(ctx->vg_kink);
	if (ctx->vg_rough1)
		MEM_freeN(ctx->vg_rough1);
	if (ctx->vg_rough2)
		MEM_freeN(ctx->vg_rough2);
	if (ctx->vg_roughe)
		MEM_freeN(ctx->vg_roughe);
	if (ctx->vg_twist)
		MEM_freeN(ctx->vg_twist);

	if (ctx->sim.psys->lattice_deform_data) {
		end_latt_deform(ctx->sim.psys->lattice_deform_data);
		ctx->sim.psys->lattice_deform_data = NULL;
	}

	/* distribution */
	if (ctx->jit) MEM_freeN(ctx->jit);
	if (ctx->jitoff) MEM_freeN(ctx->jitoff);
	if (ctx->weight) MEM_freeN(ctx->weight);
	if (ctx->index) MEM_freeN(ctx->index);
	if (ctx->seams) MEM_freeN(ctx->seams);
	//if (ctx->vertpart) MEM_freeN(ctx->vertpart);
	BLI_kdtree_3d_free(ctx->tree);

	if (ctx->clumpcurve != NULL) {
		curvemapping_free(ctx->clumpcurve);
	}
	if (ctx->roughcurve != NULL) {
		curvemapping_free(ctx->roughcurve);
	}
	if (ctx->twistcurve != NULL) {
		curvemapping_free(ctx->twistcurve);
	}
}

static void initialize_particle_texture(ParticleSimulationData *sim, ParticleData *pa, int p)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	ParticleTexture ptex;

	psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);

	switch (part->type) {
		case PART_EMITTER:
			if (ptex.exist < psys_frand(psys, p + 125)) {
				pa->flag |= PARS_UNEXIST;
			}
			pa->time = part->sta + (part->end - part->sta)*ptex.time;
			break;
		case PART_HAIR:
			if (ptex.exist < psys_frand(psys, p + 125)) {
				pa->flag |= PARS_UNEXIST;
			}
			pa->time = 0.f;
			break;
		case PART_FLUID:
			break;
	}
}

/* set particle parameters that don't change during particle's life */
void initialize_particle(ParticleSimulationData *sim, ParticleData *pa)
{
	ParticleSettings *part = sim->psys->part;
	float birth_time = (float)(pa - sim->psys->particles) / (float)sim->psys->totpart;

	pa->flag &= ~PARS_UNEXIST;
	pa->time = part->sta + (part->end - part->sta) * birth_time;

	pa->hair_index = 0;
	/* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
	/* usage other than straight after distribute has to handle this index by itself - jahka*/
	//pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we don't have a derived mesh face */
}

static void initialize_all_particles(ParticleSimulationData *sim)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	/* Grid distributionsets UNEXIST flag, need to take care of
	 * it here because later this flag is being reset.
	 *
	 * We can't do it for any distribution, because it'll then
	 * conflict with texture influence, which does not free
	 * unexisting particles and only sets flag.
	 *
	 * It's not so bad, because only grid distribution sets
	 * UNEXIST flag.
	 */
	const bool emit_from_volume_grid = (part->distr == PART_DISTR_GRID) &&
	                                   (!ELEM(part->from, PART_FROM_VERT, PART_FROM_CHILD));
	PARTICLE_P;
	LOOP_PARTICLES {
		if (!(emit_from_volume_grid && (pa->flag & PARS_UNEXIST) != 0)) {
			initialize_particle(sim, pa);
		}
	}
}

static void free_unexisting_particles(ParticleSimulationData *sim)
{
	ParticleSystem *psys = sim->psys;
	PARTICLE_P;

	psys->totunexist = 0;

	LOOP_PARTICLES {
		if (pa->flag & PARS_UNEXIST) {
			psys->totunexist++;
		}
	}

	if (psys->totpart && psys->totunexist == psys->totpart) {
		if (psys->particles->boid)
			MEM_freeN(psys->particles->boid);

		MEM_freeN(psys->particles);
		psys->particles = NULL;
		psys->totpart = psys->totunexist = 0;
	}

	if (psys->totunexist) {
		int newtotpart = psys->totpart - psys->totunexist;
		ParticleData *npa, *newpars;

		npa = newpars = MEM_callocN(newtotpart * sizeof(ParticleData), "particles");

		for (p=0, pa=psys->particles; p<newtotpart; p++, pa++, npa++) {
			while (pa->flag & PARS_UNEXIST)
				pa++;

			memcpy(npa, pa, sizeof(ParticleData));
		}

		if (psys->particles->boid)
			MEM_freeN(psys->particles->boid);
		MEM_freeN(psys->particles);
		psys->particles = newpars;
		psys->totpart -= psys->totunexist;

		if (psys->particles->boid) {
			BoidParticle *newboids = MEM_callocN(psys->totpart * sizeof(BoidParticle), "boid particles");

			LOOP_PARTICLES {
				pa->boid = newboids++;
			}

		}
	}
}

static void get_angular_velocity_vector(short avemode, ParticleKey *state, float vec[3])
{
	switch (avemode) {
		case PART_AVE_VELOCITY:
			copy_v3_v3(vec, state->vel);
			break;
		case PART_AVE_HORIZONTAL:
		{
			float zvec[3];
			zvec[0] = zvec[1] = 0;
			zvec[2] = 1.f;
			cross_v3_v3v3(vec, state->vel, zvec);
			break;
		}
		case PART_AVE_VERTICAL:
		{
			float zvec[3], temp[3];
			zvec[0] = zvec[1] = 0;
			zvec[2] = 1.f;
			cross_v3_v3v3(temp, state->vel, zvec);
			cross_v3_v3v3(vec, temp, state->vel);
			break;
		}
		case PART_AVE_GLOBAL_X:
			vec[0] = 1.f;
			vec[1] = vec[2] = 0;
			break;
		case PART_AVE_GLOBAL_Y:
			vec[1] = 1.f;
			vec[0] = vec[2] = 0;
			break;
		case PART_AVE_GLOBAL_Z:
			vec[2] = 1.f;
			vec[0] = vec[1] = 0;
			break;
	}
}

void psys_get_birth_coords(ParticleSimulationData *sim, ParticleData *pa, ParticleKey *state, float dtime, float cfra)
{
	Object *ob = sim->ob;
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	ParticleTexture ptex;
	float fac, phasefac, nor[3] = {0,0,0},loc[3],vel[3] = {0.0,0.0,0.0},rot[4],q2[4];
	float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3] = {0.0,0.0,0.0};
	float x_vec[3] = {1.0,0.0,0.0}, utan[3] = {0.0,1.0,0.0}, vtan[3] = {0.0,0.0,1.0}, rot_vec[3] = {0.0,0.0,0.0};
	float q_phase[4];

	const bool use_boids = ((part->phystype == PART_PHYS_BOIDS) &&
	                        (pa->boid != NULL));
	const bool use_tangents = ((use_boids == false) &&
	                           ((part->tanfac != 0.0f) || (part->rotmode == PART_ROT_NOR_TAN)));

	int p = pa - psys->particles;

	/* get birth location from object		*/
	if (use_tangents)
		psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0);
	else
		psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,0,0,0);

	/* get possible textural influence */
	psys_get_texture(sim, pa, &ptex, PAMAP_IVEL, cfra);

	/* particles live in global space so	*/
	/* let's convert:						*/
	/* -location							*/
	mul_m4_v3(ob->obmat, loc);

	/* -normal								*/
	mul_mat3_m4_v3(ob->obmat, nor);
	normalize_v3(nor);

	/* -tangent								*/
	if (use_tangents) {
		//float phase=vg_rot?2.0f*(psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_rot)-0.5f):0.0f;
		float phase=0.0f;
		mul_v3_fl(vtan,-cosf((float)M_PI*(part->tanphase+phase)));
		fac= -sinf((float)M_PI*(part->tanphase+phase));
		madd_v3_v3fl(vtan, utan, fac);

		mul_mat3_m4_v3(ob->obmat,vtan);

		copy_v3_v3(utan, nor);
		mul_v3_fl(utan,dot_v3v3(vtan,nor));
		sub_v3_v3(vtan, utan);

		normalize_v3(vtan);
	}


	/* -velocity (boids need this even if there's no random velocity) */
	if (part->randfac != 0.0f || (part->phystype==PART_PHYS_BOIDS && pa->boid)) {
		r_vel[0] = 2.0f * (psys_frand(psys, p + 10) - 0.5f);
		r_vel[1] = 2.0f * (psys_frand(psys, p + 11) - 0.5f);
		r_vel[2] = 2.0f * (psys_frand(psys, p + 12) - 0.5f);

		mul_mat3_m4_v3(ob->obmat, r_vel);
		normalize_v3(r_vel);
	}

	/* -angular velocity					*/
	if (part->avemode==PART_AVE_RAND) {
		r_ave[0] = 2.0f * (psys_frand(psys, p + 13) - 0.5f);
		r_ave[1] = 2.0f * (psys_frand(psys, p + 14) - 0.5f);
		r_ave[2] = 2.0f * (psys_frand(psys, p + 15) - 0.5f);

		mul_mat3_m4_v3(ob->obmat,r_ave);
		normalize_v3(r_ave);
	}

	/* -rotation							*/
	if (part->randrotfac != 0.0f) {
		r_rot[0] = 2.0f * (psys_frand(psys, p + 16) - 0.5f);
		r_rot[1] = 2.0f * (psys_frand(psys, p + 17) - 0.5f);
		r_rot[2] = 2.0f * (psys_frand(psys, p + 18) - 0.5f);
		r_rot[3] = 2.0f * (psys_frand(psys, p + 19) - 0.5f);
		normalize_qt(r_rot);

		mat4_to_quat(rot,ob->obmat);
		mul_qt_qtqt(r_rot,r_rot,rot);
	}

	if (use_boids) {
		float dvec[3], q[4], mat[3][3];

		copy_v3_v3(state->co,loc);

		/* boids don't get any initial velocity  */
		zero_v3(state->vel);

		/* boids store direction in ave */
		if (fabsf(nor[2])==1.0f) {
			sub_v3_v3v3(state->ave, loc, ob->obmat[3]);
			normalize_v3(state->ave);
		}
		else {
			copy_v3_v3(state->ave, nor);
		}

		/* calculate rotation matrix */
		project_v3_v3v3(dvec, r_vel, state->ave);
		sub_v3_v3v3(mat[0], state->ave, dvec);
		normalize_v3(mat[0]);
		negate_v3_v3(mat[2], r_vel);
		normalize_v3(mat[2]);
		cross_v3_v3v3(mat[1], mat[2], mat[0]);

		/* apply rotation */
		mat3_to_quat_is_ok( q,mat);
		copy_qt_qt(state->rot, q);
	}
	else {
		/* conversion done so now we apply new:	*/
		/* -velocity from:						*/

		/*		*reactions						*/
		if (dtime > 0.f) {
			sub_v3_v3v3(vel, pa->state.vel, pa->prev_state.vel);
		}

		/*		*emitter velocity				*/
		if (dtime != 0.f && part->obfac != 0.f) {
			sub_v3_v3v3(vel, loc, state->co);
			mul_v3_fl(vel, part->obfac/dtime);
		}

		/*		*emitter normal					*/
		if (part->normfac != 0.f)
			madd_v3_v3fl(vel, nor, part->normfac);

		/*		*emitter tangent				*/
		if (sim->psmd && part->tanfac != 0.f)
			madd_v3_v3fl(vel, vtan, part->tanfac);

		/*		*emitter object orientation		*/
		if (part->ob_vel[0] != 0.f) {
			normalize_v3_v3(vec, ob->obmat[0]);
			madd_v3_v3fl(vel, vec, part->ob_vel[0]);
		}
		if (part->ob_vel[1] != 0.f) {
			normalize_v3_v3(vec, ob->obmat[1]);
			madd_v3_v3fl(vel, vec, part->ob_vel[1]);
		}
		if (part->ob_vel[2] != 0.f) {
			normalize_v3_v3(vec, ob->obmat[2]);
			madd_v3_v3fl(vel, vec, part->ob_vel[2]);
		}

		/*		*texture						*/
		/* TODO	*/

		/*		*random							*/
		if (part->randfac != 0.f)
			madd_v3_v3fl(vel, r_vel, part->randfac);

		/*		*particle						*/
		if (part->partfac != 0.f)
			madd_v3_v3fl(vel, p_vel, part->partfac);

		mul_v3_v3fl(state->vel, vel, ptex.ivel);

		/* -location from emitter				*/
		copy_v3_v3(state->co,loc);

		/* -rotation							*/
		unit_qt(state->rot);

		if (part->rotmode) {
			bool use_global_space;

			/* create vector into which rotation is aligned */
			switch (part->rotmode) {
				case PART_ROT_NOR:
				case PART_ROT_NOR_TAN:
					copy_v3_v3(rot_vec, nor);
					use_global_space = false;
					break;
				case PART_ROT_VEL:
					copy_v3_v3(rot_vec, vel);
					use_global_space = true;
					break;
				case PART_ROT_GLOB_X:
				case PART_ROT_GLOB_Y:
				case PART_ROT_GLOB_Z:
					rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
					use_global_space = true;
					break;
				case PART_ROT_OB_X:
				case PART_ROT_OB_Y:
				case PART_ROT_OB_Z:
					copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
					use_global_space = false;
					break;
				default:
					use_global_space = true;
					break;
			}

			/* create rotation quat */


			if (use_global_space) {
				negate_v3(rot_vec);
				vec_to_quat(q2, rot_vec, OB_POSX, OB_POSZ);

				/* randomize rotation quat */
				if (part->randrotfac != 0.0f) {
					interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
				}
				else {
					copy_qt_qt(rot, q2);
				}
			}
			else {
				/* calculate rotation in local-space */
				float q_obmat[4];
				float q_imat[4];

				mat4_to_quat(q_obmat, ob->obmat);
				invert_qt_qt_normalized(q_imat, q_obmat);


				if (part->rotmode != PART_ROT_NOR_TAN) {
					float rot_vec_local[3];

					/* rot_vec */
					negate_v3(rot_vec);
					copy_v3_v3(rot_vec_local, rot_vec);
					mul_qt_v3(q_imat, rot_vec_local);
					normalize_v3(rot_vec_local);

					vec_to_quat(q2, rot_vec_local, OB_POSX, OB_POSZ);
				}
				else {
					/* (part->rotmode == PART_ROT_NOR_TAN) */
					float tmat[3][3];

					/* note: utan_local is not taken from 'utan', we calculate from rot_vec/vtan */
					/* note: it looks like rotation phase may be applied twice (once with vtan, again below)
					 * however this isn't the case - campbell */
					float *rot_vec_local = tmat[0];
					float *vtan_local    = tmat[1];
					float *utan_local    = tmat[2];

					/* use tangents */
					BLI_assert(use_tangents == true);

					/* rot_vec */
					copy_v3_v3(rot_vec_local, rot_vec);
					mul_qt_v3(q_imat, rot_vec_local);

					/* vtan_local */
					copy_v3_v3(vtan_local, vtan);  /* flips, cant use */
					mul_qt_v3(q_imat, vtan_local);

					/* ensure orthogonal matrix (rot_vec aligned) */
					cross_v3_v3v3(utan_local, vtan_local, rot_vec_local);
					cross_v3_v3v3(vtan_local, utan_local, rot_vec_local);

					/* note: no need to normalize */
					mat3_to_quat(q2, tmat);
				}

				/* randomize rotation quat */
				if (part->randrotfac != 0.0f) {
					mul_qt_qtqt(r_rot, r_rot, q_imat);
					interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
				}
				else {
					copy_qt_qt(rot, q2);
				}

				mul_qt_qtqt(rot, q_obmat, rot);
			}

			/* rotation phase */
			phasefac = part->phasefac;
			if (part->randphasefac != 0.0f)
				phasefac += part->randphasefac * psys_frand(psys, p + 20);
			axis_angle_to_quat( q_phase,x_vec, phasefac*(float)M_PI);

			/* combine base rotation & phase */
			mul_qt_qtqt(state->rot, rot, q_phase);
		}

		/* -angular velocity					*/

		zero_v3(state->ave);

		if (part->avemode) {
			if (part->avemode == PART_AVE_RAND)
				copy_v3_v3(state->ave, r_ave);
			else
				get_angular_velocity_vector(part->avemode, state, state->ave);

			normalize_v3(state->ave);
			mul_v3_fl(state->ave, part->avefac);
		}
	}
}

/* recursively evaluate emitter parent anim at cfra */
static void evaluate_emitter_anim(struct Depsgraph *depsgraph, Scene *scene, Object *ob, float cfra)
{
	if (ob->parent)
		evaluate_emitter_anim(depsgraph, scene, ob->parent, cfra);

	BKE_object_where_is_calc_time(depsgraph, scene, ob, cfra);
}

/* sets particle to the emitter surface with initial velocity & rotation */
void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part;
	ParticleTexture ptex;
	int p = pa - psys->particles;
	part=psys->part;

	/* get precise emitter matrix if particle is born */
	if (part->type != PART_HAIR && dtime > 0.f && pa->time < cfra && pa->time >= sim->psys->cfra) {
		evaluate_emitter_anim(sim->depsgraph, sim->scene, sim->ob, pa->time);

		psys->flag |= PSYS_OB_ANIM_RESTORE;
	}

	psys_get_birth_coords(sim, pa, &pa->state, dtime, cfra);

	/* Initialize particle settings which depends on texture.
	 *
	 * We could only do it now because we'll need to know coordinate
	 * before sampling the texture.
	 */
	initialize_particle_texture(sim, pa, p);

	if (part->phystype==PART_PHYS_BOIDS && pa->boid) {
		BoidParticle *bpa = pa->boid;

		/* and gravity in r_ve */
		bpa->gravity[0] = bpa->gravity[1] = 0.0f;
		bpa->gravity[2] = -1.0f;
		if ((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) &&
		    (sim->scene->physics_settings.gravity[2] != 0.0f))
		{
			bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
		}

		bpa->data.health = part->boids->health;
		bpa->data.mode = eBoidMode_InAir;
		bpa->data.state_id = ((BoidState*)part->boids->states.first)->id;
		bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f;
	}

	if (part->type == PART_HAIR) {
		pa->lifetime = 100.0f;
	}
	else {
		/* initialize the lifetime, in case the texture coordinates
		 * are from Particles/Strands, which would cause undefined values
		 */
		pa->lifetime = part->lifetime * (1.0f - part->randlife * psys_frand(psys, p + 21));
		pa->dietime = pa->time + pa->lifetime;

		/* get possible textural influence */
		psys_get_texture(sim, pa, &ptex, PAMAP_LIFE, cfra);

		pa->lifetime = part->lifetime * ptex.life;

		if (part->randlife != 0.0f)
			pa->lifetime *= 1.0f - part->randlife * psys_frand(psys, p + 21);
	}

	pa->dietime = pa->time + pa->lifetime;

	if ((sim->psys->pointcache) &&
	    (sim->psys->pointcache->flag & PTCACHE_BAKED) &&
	    (sim->psys->pointcache->mem_cache.first))
	{
		float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
		pa->dietime = MIN2(pa->dietime, dietime);
	}

	if (pa->time > cfra)
		pa->alive = PARS_UNBORN;
	else if (pa->dietime <= cfra)
		pa->alive = PARS_DEAD;
	else
		pa->alive = PARS_ALIVE;

	pa->state.time = cfra;
}
static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
{
	ParticleData *pa;
	int p, totpart=sim->psys->totpart;

	for (p=from, pa=sim->psys->particles+from; p<totpart; p++, pa++)
		reset_particle(sim, pa, dtime, cfra);
}
/************************************************/
/*			Particle targets					*/
/************************************************/
ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
{
	ParticleSystem *psys = NULL;

	if (pt->ob == NULL || pt->ob == ob)
		psys = BLI_findlink(&ob->particlesystem, pt->psys-1);
	else
		psys = BLI_findlink(&pt->ob->particlesystem, pt->psys-1);

	if (psys)
		pt->flag |= PTARGET_VALID;
	else
		pt->flag &= ~PTARGET_VALID;

	return psys;
}
/************************************************/
/*			Keyed particles						*/
/************************************************/
/* Counts valid keyed targets */
void psys_count_keyed_targets(ParticleSimulationData *sim)
{
	ParticleSystem *psys = sim->psys, *kpsys;
	ParticleTarget *pt = psys->targets.first;
	int keys_valid = 1;
	psys->totkeyed = 0;

	for (; pt; pt=pt->next) {
		kpsys = psys_get_target_system(sim->ob, pt);

		if (kpsys && kpsys->totpart) {
			psys->totkeyed += keys_valid;
			if (psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
				psys->totkeyed += 1;
		}
		else {
			keys_valid = 0;
		}
	}

	psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
}

static void set_keyed_keys(ParticleSimulationData *sim)
{
	ParticleSystem *psys = sim->psys;
	ParticleSimulationData ksim= {0};
	ParticleTarget *pt;
	PARTICLE_P;
	ParticleKey *key;
	int totpart = psys->totpart, k, totkeys = psys->totkeyed;
	int keyed_flag = 0;

	ksim.depsgraph = sim->depsgraph;
	ksim.scene = sim->scene;

	/* no proper targets so let's clear and bail out */
	if (psys->totkeyed==0) {
		free_keyed_keys(psys);
		psys->flag &= ~PSYS_KEYED;
		return;
	}

	if (totpart && psys->particles->totkey != totkeys) {
		free_keyed_keys(psys);

		key = MEM_callocN(totpart*totkeys*sizeof(ParticleKey), "Keyed keys");

		LOOP_PARTICLES {
			pa->keys = key;
			pa->totkey = totkeys;
			key += totkeys;
		}
	}

	psys->flag &= ~PSYS_KEYED;


	pt = psys->targets.first;
	for (k=0; k<totkeys; k++) {
		ksim.ob = pt->ob ? pt->ob : sim->ob;
		ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
		keyed_flag = (ksim.psys->flag & PSYS_KEYED);
		ksim.psys->flag &= ~PSYS_KEYED;

		LOOP_PARTICLES {
			key = pa->keys + k;
			key->time = -1.0; /* use current time */

			psys_get_particle_state(&ksim, p%ksim.psys->totpart, key, 1);

			if (psys->flag & PSYS_KEYED_TIMING) {
				key->time = pa->time + pt->time;
				if (pt->duration != 0.0f && k+1 < totkeys) {
					copy_particle_key(key+1, key, 1);
					(key+1)->time = pa->time + pt->time + pt->duration;
				}
			}
			else if (totkeys > 1)
				key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
			else
				key->time = pa->time;
		}

		if (psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
			k++;

		ksim.psys->flag |= keyed_flag;

		pt = (pt->next && pt->next->flag & PTARGET_VALID) ? pt->next : psys->targets.first;
	}

	psys->flag |= PSYS_KEYED;
}

/************************************************/
/*			Point Cache							*/
/************************************************/
void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
{
	PointCache *cache = psys->pointcache;

	if (cache->flag & PTCACHE_DISK_CACHE && BLI_listbase_is_empty(&cache->mem_cache)) {
		PTCacheID pid;
		BKE_ptcache_id_from_particles(&pid, ob, psys);
		cache->flag &= ~PTCACHE_DISK_CACHE;
		BKE_ptcache_disk_to_mem(&pid);
		cache->flag |= PTCACHE_DISK_CACHE;
	}
}
static void psys_clear_temp_pointcache(ParticleSystem *psys)
{
	if (psys->pointcache->flag & PTCACHE_DISK_CACHE)
		BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
}
void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
{
	ParticleSettings *part = psys->part;

	*sfra = max_ii(1, (int)part->sta);
	*efra = min_ii((int)(part->end + part->lifetime + 1.0f), max_ii(scene->r.pefra, scene->r.efra));
}

/************************************************/
/*			Effectors							*/
/************************************************/
static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
{
	if (psys) {
		PARTICLE_P;
		int totpart = 0;
		bool need_rebuild;

		BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_READ);
		need_rebuild = !psys->bvhtree || psys->bvhtree_frame != cfra;
		BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);

		if (need_rebuild) {
			LOOP_SHOWN_PARTICLES {
				totpart++;
			}

			BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_WRITE);

			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;

			BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
		}
	}
}
void psys_update_particle_tree(ParticleSystem *psys, float cfra)
{
	if (psys) {
		PARTICLE_P;
		int totpart = 0;

		if (!psys->tree || psys->tree_frame != cfra) {
			LOOP_SHOWN_PARTICLES {
				totpart++;
			}

			BLI_kdtree_3d_free(psys->tree);
			psys->tree = BLI_kdtree_3d_new(psys->totpart);

			LOOP_SHOWN_PARTICLES {
				if (pa->alive == PARS_ALIVE) {
					if (pa->state.time == cfra)
						BLI_kdtree_3d_insert(psys->tree, p, pa->prev_state.co);
					else
						BLI_kdtree_3d_insert(psys->tree, p, pa->state.co);
				}
			}
			BLI_kdtree_3d_balance(psys->tree);

			psys->tree_frame = cfra;
		}
	}
}

static void psys_update_effectors(ParticleSimulationData *sim)
{
	BKE_effectors_free(sim->psys->effectors);
	sim->psys->effectors = BKE_effectors_create(sim->depsgraph,
	                                            sim->ob, sim->psys,
	                                            sim->psys->part->effector_weights);
	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)
{
#define ZERO_F43 {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}

	ParticleKey states[5];
	float force[3], acceleration[3], impulse[3], dx[4][3] = ZERO_F43, dv[4][3] = ZERO_F43, oldpos[3];
	float pa_mass= (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
	int i, steps=1;
	int integrator = part->integrator;

#undef ZERO_F43

	copy_v3_v3(oldpos, pa->state.co);

	/* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
	if (pa->prev_state.time < 0.f && integrator == PART_INT_VERLET)
		integrator = PART_INT_EULER;

	switch (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;
	}

	for (i=0; i<steps; i++) {
		copy_particle_key(states + i, &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 (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
 *
 * In theory, there could be unlimited implementation of SPH simulators
 *
 * This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
 *
 * Titled: Particle-based Viscoelastic Fluid Simulation.
 * Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
 * Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
 *
 * Presented at Siggraph, (2005)
 *
 * ********************************************************************************************************/
#define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
{
	/* Are more refs required? */
	if (psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
		psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
		psys->fluid_springs = (ParticleSpring*)MEM_callocN(psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
	}
	else if (psys->tot_fluidsprings == psys->alloc_fluidsprings) {
		/* Double the number of refs allocated */
		psys->alloc_fluidsprings *= 2;
		psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
	}

	memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
	psys->tot_fluidsprings++;

	return psys->fluid_springs + psys->tot_fluidsprings - 1;
}
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];

	psys->tot_fluidsprings--;

	if (psys->tot_fluidsprings < psys->alloc_fluidsprings/2 && psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE) {
		psys->alloc_fluidsprings /= 2;
		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;
	}

	/* 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;
	int i = 0;

	springhash = BLI_edgehash_new_ex(__func__, psys->tot_fluidsprings);

	for (i=0, spring=psys->fluid_springs; i<psys->tot_fluidsprings; i++, spring++)
		BLI_edgehash_insert(springhash, spring->particle_index[0], spring->particle_index[1], POINTER_FROM_INT(i+1));

	return springhash;
}

#define SPH_NEIGHBORS 512
typedef struct SPHNeighbor {
	ParticleSystem *psys;
	int index;
} SPHNeighbor;

typedef struct SPHRangeData {
	SPHNeighbor neighbors[SPH_NEIGHBORS];
	int tot_neighbors;

	float* data;

	ParticleSystem *npsys;
	ParticleData *pa;

	float h;
	float mass;
	float massfac;
	int use_size;
} SPHRangeData;

static void sph_evaluate_func(BVHTree *tree, ParticleSystem **psys, float co[3], SPHRangeData *pfr, float interaction_radius, BVHTree_RangeQuery callback)
{
	int i;

	pfr->tot_neighbors = 0;

	for (i=0; i < 10 && psys[i]; i++) {
		pfr->npsys    = psys[i];
		pfr->massfac  = psys[i]->part->mass / pfr->mass;
		pfr->use_size = psys[i]->part->flag & PART_SIZEMASS;

		if (tree) {
			BLI_bvhtree_range_query(tree, co, interaction_radius, callback, pfr);
			break;
		}
		else {
			BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_READ);

			BLI_bvhtree_range_query(psys[i]->bvhtree, co, interaction_radius, callback, pfr);

			BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
		}
	}
}
static void sph_density_accum_cb(void *userdata, int index, const float co[3], float squared_dist)
{
	SPHRangeData *pfr = (SPHRangeData *)userdata;
	ParticleData *npa = pfr->npsys->particles + index;
	float q;
	float dist;

	UNUSED_VARS(co);

	if (npa == pfr->pa || squared_dist < FLT_EPSILON)
		return;

	/* Ugh! One particle has too many neighbors! If some aren't taken into
	 * account, the forces will be biased by the tree search order. This
	 * effectively adds energy to the system, and results in a churning motion.
	 * But, we have to stop somewhere, and it's not the end of the world.
	 * - jahka and z0r
	 */
	if (pfr->tot_neighbors >= SPH_NEIGHBORS)
		return;

	pfr->neighbors[pfr->tot_neighbors].index = index;
	pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
	pfr->tot_neighbors++;

	dist = sqrtf(squared_dist);
	q = (1.f - dist/pfr->h) * pfr->massfac;

	if (pfr->use_size)
		q *= npa->size;

	pfr->data[0] += q*q;
	pfr->data[1] += q*q*q;
}

/*
 * Find the Courant number for an SPH particle (used for adaptive time step).
 */
static void sph_particle_courant(SPHData *sphdata, SPHRangeData *pfr)
{
	ParticleData *pa, *npa;
	int i;
	float flow[3], offset[3], dist;

	zero_v3(flow);

	dist = 0.0f;
	if (pfr->tot_neighbors > 0) {
		pa = pfr->pa;
		for (i=0; i < pfr->tot_neighbors; i++) {
			npa = pfr->neighbors[i].psys->particles + pfr->neighbors[i].index;
			sub_v3_v3v3(offset, pa->prev_state.co, npa->prev_state.co);
			dist += len_v3(offset);
			add_v3_v3(flow, npa->prev_state.vel);
		}
		dist += sphdata->psys[0]->part->fluid->radius; // TODO: remove this? - z0r
		sphdata->element_size = dist / pfr->tot_neighbors;
		mul_v3_v3fl(sphdata->flow, flow, 1.0f / pfr->tot_neighbors);
	}
	else {
		sphdata->element_size = FLT_MAX;
		copy_v3_v3(sphdata->flow, flow);
	}
}
static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(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 *gravity = sphdata->gravity;
	EdgeHash *springhash = sphdata->eh;

	float q, u, rij, dv[3];
	float pressure, near_pressure;

	float visc = fluid->viscosity_omega;
	float stiff_visc = fluid->viscosity_beta * (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);

	float inv_mass = 1.0f / sphdata->mass;
	float spring_constant = fluid->spring_k;

	/* 4.0 seems to be a pretty good value */
	float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
	float h = interaction_radius * sphdata->hfac;
	/* 4.77 is an experimentally determined density factor */
	float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f);
	float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);

	float stiffness = fluid->stiffness_k;
	float stiffness_near_fac = fluid->stiffness_knear * (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);

	ParticleData *npa;
	float vec[3];
	float vel[3];
	float co[3];
	float data[2];
	float density, near_density;

	int i, spring_index, index = pa - psys[0]->particles;

	data[0] = data[1] = 0;
	pfr.data = data;
	pfr.h = h;
	pfr.pa = pa;
	pfr.mass = sphdata->mass;

	sph_evaluate_func( NULL, psys, state->co, &pfr, interaction_radius, sph_density_accum_cb);

	density = data[0];
	near_density = data[1];

	pressure =  stiffness * (density - rest_density);
	near_pressure = stiffness_near_fac * near_density;

	pfn = pfr.neighbors;
	for (i=0; i<pfr.tot_neighbors; i++, pfn++) {
		npa = pfn->psys->particles + pfn->index;

		madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);

		sub_v3_v3v3(vec, co, state->co);
		rij = normalize_v3(vec);

		q = (1.f - rij/h) * pfn->psys->part->mass * inv_mass;

		if (pfn->psys->part->flag & PART_SIZEMASS)
			q *= npa->size;

		copy_v3_v3(vel, npa->prev_state.vel);

		/* Double Density Relaxation */
		madd_v3_v3fl(force, vec, -(pressure + near_pressure*q)*q);

		/* Viscosity */
		if (visc > 0.f  || stiff_visc > 0.f) {
			sub_v3_v3v3(dv, vel, state->vel);
			u = dot_v3v3(vec, dv);

			if (u < 0.f && visc > 0.f)
				madd_v3_v3fl(force, vec, 0.5f * q * visc * u );

			if (u > 0.f && stiff_visc > 0.f)
				madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u );
		}

		if (spring_constant > 0.f) {
			/* Viscoelastic spring force */
			if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
				/* BLI_edgehash_lookup appears to be thread-safe. - z0r */
				spring_index = POINTER_AS_INT(BLI_edgehash_lookup(springhash, index, pfn->index));

				if (spring_index) {
					spring = psys[0]->fluid_springs + spring_index - 1;

					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 is not thread-safe. - z0r */
					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));
			}
		}
	}

	/* Artificial buoyancy force in negative gravity direction  */
	if (fluid->buoyancy > 0.f && gravity)
		madd_v3_v3fl(force, gravity, fluid->buoyancy * (density-rest_density));

	if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF)
		sph_particle_courant(sphdata, &pfr);
	sphdata->pass++;
}

static void sphclassical_density_accum_cb(void *userdata, int index, const float co[3], float UNUSED(squared_dist))
{
	SPHRangeData *pfr = (SPHRangeData *)userdata;
	ParticleData *npa = pfr->npsys->particles + index;
	float q;
	float qfac = 21.0f / (256.f * (float)M_PI);
	float rij, rij_h;
	float vec[3];

	/* Exclude particles that are more than 2h away. Can't use squared_dist here
	 * because it is not accurate enough. Use current state, i.e. the output of
	 * basic_integrate() - z0r */
	sub_v3_v3v3(vec, npa->state.co, co);
	rij = len_v3(vec);
	rij_h = rij / pfr->h;
	if (rij_h > 2.0f)
		return;

	/* Smoothing factor. Utilise the Wendland kernel. gnuplot:
	 *     q1(x) = (2.0 - x)**4 * ( 1.0 + 2.0 * x)
	 *     plot [0:2] q1(x) */
	q  = qfac / pow3f(pfr->h) * pow4f(2.0f - rij_h) * ( 1.0f + 2.0f * rij_h);
	q *= pfr->npsys->part->mass;

	if (pfr->use_size)
		q *= pfr->pa->size;

	pfr->data[0] += q;
	pfr->data[1] += q / npa->sphdensity;
}

static void sphclassical_neighbour_accum_cb(void *userdata, int index, const float co[3], float UNUSED(squared_dist))
{
	SPHRangeData *pfr = (SPHRangeData *)userdata;
	ParticleData *npa = pfr->npsys->particles + index;
	float rij, rij_h;
	float vec[3];

	if (pfr->tot_neighbors >= SPH_NEIGHBORS)
		return;

	/* Exclude particles that are more than 2h away. Can't use squared_dist here
	 * because it is not accurate enough. Use current state, i.e. the output of
	 * basic_integrate() - z0r */
	sub_v3_v3v3(vec, npa->state.co, co);
	rij = len_v3(vec);
	rij_h = rij / pfr->h;
	if (rij_h > 2.0f)
		return;

	pfr->neighbors[pfr->tot_neighbors].index = index;
	pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
	pfr->tot_neighbors++;
}
static void sphclassical_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
{
	SPHData *sphdata = (SPHData *)sphdata_v;
	ParticleSystem **psys = sphdata->psys;
	ParticleData *pa = sphdata->pa;
	SPHFluidSettings *fluid = psys[0]->part->fluid;
	SPHRangeData pfr;
	SPHNeighbor *pfn;
	float *gravity = sphdata->gravity;

	float dq, u, rij, dv[3];
	float pressure, npressure;

	float visc = fluid->viscosity_omega;

	float interaction_radius;
	float h, hinv;
	/* 4.77 is an experimentally determined density factor */
	float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.0f);

	// Use speed of sound squared
	float stiffness = pow2f(fluid->stiffness_k);

	ParticleData *npa;
	float vec[3];
	float co[3];
	float pressureTerm;

	int i;

	float qfac2 = 42.0f / (256.0f * (float)M_PI);
	float rij_h;

	/* 4.0 here is to be consistent with previous formulation/interface */
	interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
	h = interaction_radius * sphdata->hfac;
	hinv = 1.0f / h;

	pfr.h = h;
	pfr.pa = pa;

	sph_evaluate_func(NULL, psys, state->co, &pfr, interaction_radius, sphclassical_neighbour_accum_cb);
	pressure =  stiffness * (pow7f(pa->sphdensity / rest_density) - 1.0f);

	/* multiply by mass so that we return a force, not accel */
	qfac2 *= sphdata->mass / pow3f(pfr.h);

	pfn = pfr.neighbors;
	for (i = 0; i < pfr.tot_neighbors; i++, pfn++) {
		npa = pfn->psys->particles + pfn->index;
		if (npa == pa) {
			/* we do not contribute to ourselves */
			continue;
		}

		/* Find vector to neighbor. Exclude particles that are more than 2h
		 * away. Can't use current state here because it may have changed on
		 * another thread - so do own mini integration. Unlike basic_integrate,
		 * SPH integration depends on neighboring particles. - z0r */
		madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
		sub_v3_v3v3(vec, co, state->co);
		rij = normalize_v3(vec);
		rij_h = rij / pfr.h;
		if (rij_h > 2.0f)
			continue;

		npressure = stiffness * (pow7f(npa->sphdensity / rest_density) - 1.0f);

		/* First derivative of smoothing factor. Utilise the Wendland kernel.
		 * gnuplot:
		 *     q2(x) = 2.0 * (2.0 - x)**4 - 4.0 * (2.0 - x)**3 * (1.0 + 2.0 * x)
		 *     plot [0:2] q2(x)
		 * Particles > 2h away are excluded above. */
		dq = qfac2 * (2.0f * pow4f(2.0f - rij_h) - 4.0f * pow3f(2.0f - rij_h) * (1.0f + 2.0f * rij_h)  );

		if (pfn->psys->part->flag & PART_SIZEMASS)
			dq *= npa->size;

		pressureTerm = pressure / pow2f(pa->sphdensity) + npressure / pow2f(npa->sphdensity);

		/* Note that 'minus' is removed, because vec = vecBA, not vecAB.
		 * This applies to the viscosity calculation below, too. */
		madd_v3_v3fl(force, vec, pressureTerm * dq);

		/* Viscosity */
		if (visc > 0.0f) {
			sub_v3_v3v3(dv, npa->prev_state.vel, pa->prev_state.vel);
			u = dot_v3v3(vec, dv);
			/* Apply parameters */
			u *= -dq * hinv * visc / (0.5f * npa->sphdensity + 0.5f * pa->sphdensity);
			madd_v3_v3fl(force, vec, u);
		}
	}

	/* Artificial buoyancy force in negative gravity direction  */
	if (fluid->buoyancy > 0.f && gravity)
		madd_v3_v3fl(force, gravity, fluid->buoyancy * (pa->sphdensity - rest_density));

	if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF)
		sph_particle_courant(sphdata, &pfr);
	sphdata->pass++;
}

static void sphclassical_calc_dens(ParticleData *pa, float UNUSED(dfra), SPHData *sphdata)
{
	ParticleSystem **psys = sphdata->psys;
	SPHFluidSettings *fluid = psys[0]->part->fluid;
	/* 4.0 seems to be a pretty good value */
	float interaction_radius  = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
	SPHRangeData pfr;
	float data[2];

	data[0] = 0;
	data[1] = 0;
	pfr.data = data;
	pfr.h = interaction_radius * sphdata->hfac;
	pfr.pa = pa;
	pfr.mass = sphdata->mass;

	sph_evaluate_func( NULL, psys, pa->state.co, &pfr, interaction_radius, sphclassical_density_accum_cb);
	pa->sphdensity = min_ff(max_ff(data[0], fluid->rest_density * 0.9f), fluid->rest_density * 1.1f);
}

void psys_sph_init(ParticleSimulationData *sim, SPHData *sphdata)
{
	ParticleTarget *pt;
	int i;

	// Add other coupled particle systems.
	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 (psys_uses_gravity(sim))
		sphdata->gravity = sim->scene->physics_settings.gravity;
	else
		sphdata->gravity = NULL;
	sphdata->eh = sph_springhash_build(sim->psys);

	// These per-particle values should be overridden later, but just for
	// completeness we give them default values now.
	sphdata->pa = NULL;
	sphdata->mass = 1.0f;

	if (sim->psys->part->fluid->solver == SPH_SOLVER_DDR) {
		sphdata->force_cb = sph_force_cb;
		sphdata->density_cb = sph_density_accum_cb;
		sphdata->hfac = 1.0f;
	}
	else {
		/* SPH_SOLVER_CLASSICAL */
		sphdata->force_cb = sphclassical_force_cb;
		sphdata->density_cb = sphclassical_density_accum_cb;
		sphdata->hfac = 0.5f;
	}

}

void psys_sph_finalise(SPHData *sphdata)
{
	if (sphdata->eh) {
		BLI_edgehash_free(sphdata->eh, NULL);
		sphdata->eh = NULL;
	}
}
/* Sample the density field at a point in space. */
void psys_sph_density(BVHTree *tree, SPHData *sphdata, float co[3], float vars[2])
{
	ParticleSystem **psys = sphdata->psys;
	SPHFluidSettings *fluid = psys[0]->part->fluid;
	/* 4.0 seems to be a pretty good value */
	float interaction_radius  = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
	SPHRangeData pfr;
	float density[2];

	density[0] = density[1] = 0.0f;
	pfr.data = density;
	pfr.h = interaction_radius * sphdata->hfac;
	pfr.mass = sphdata->mass;

	sph_evaluate_func(tree, psys, co, &pfr, interaction_radius, sphdata->density_cb);

	vars[0] = pfr.data[0];
	vars[1] = pfr.data[1];
}

static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, SPHData *sphdata)
{
	ParticleSettings *part = sim->psys->part;
	// float timestep = psys_get_timestep(sim); // UNUSED
	float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
	float dtime = dfra*psys_get_timestep(sim);
	// int steps = 1; // UNUSED
	float effector_acceleration[3];

	sphdata->pa = pa;
	sphdata->mass = pa_mass;
	sphdata->pass = 0;
	//sphdata.element_size and sphdata.flow are set in the callback.

	/* 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(&pa->state, &pa->prev_state, 0);

	integrate_particle(part, pa, dtime, effector_acceleration, sphdata->force_cb, sphdata);
}

/************************************************/
/*			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;
	RNG *rng = sim->rng;

	/* 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)
		BKE_effectors_apply(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);

	mul_v3_fl(force, efdata->ptex.field);
	mul_v3_fl(impulse, efdata->ptex.field);

	/* 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));

	/* brownian force */
	if (part->brownfac != 0.0f) {
		force[0] += (BLI_rng_get_float(rng)-0.5f) * part->brownfac;
		force[1] += (BLI_rng_get_float(rng)-0.5f) * part->brownfac;
		force[2] += (BLI_rng_get_float(rng)-0.5f) * part->brownfac;
	}

	if (part->flag & PART_ROT_DYN && epoint.ave)
		copy_v3_v3(pa->state.ave, epoint.ave);
}
/* 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;

	psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);

	efdata.pa = pa;
	efdata.sim = sim;

	/* 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 * efdata.ptex.damp * 25.f * dtime);

	//copy_v3_v3(pa->state.ave, states->ave);

	/* finally we do guides */
	time=(cfra-pa->time)/pa->lifetime;
	CLAMP(time, 0.0f, 1.0f);

	copy_v3_v3(tkey.co,pa->state.co);
	copy_v3_v3(tkey.vel,pa->state.vel);
	tkey.time=pa->state.time;

	if (part->type != PART_HAIR) {
		if (do_guides(sim->depsgraph, sim->psys->part, sim->psys->effectors, &tkey, p, time)) {
			copy_v3_v3(pa->state.co,tkey.co);
			/* guides don't produce valid velocity */
			sub_v3_v3v3(pa->state.vel, tkey.co, pa->prev_state.co);
			mul_v3_fl(pa->state.vel,1.0f/dtime);
			pa->state.time=tkey.time;
		}
	}
}
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, extrotfac;

	if ((part->flag & PART_ROTATIONS) == 0) {
		unit_qt(pa->state.rot);
		return;
	}

	if (part->flag & PART_ROT_DYN) {
		extrotfac = len_v3(pa->state.ave);
	}
	else {
		extrotfac = 0.0f;
	}

	if ((part->flag & PART_ROT_DYN) && ELEM(part->avemode, PART_AVE_VELOCITY, PART_AVE_HORIZONTAL, PART_AVE_VERTICAL)) {
		float angle;
		float len1 = len_v3(pa->prev_state.vel);
		float len2 = len_v3(pa->state.vel);
		float vec[3];

		if (len1 == 0.0f || len2 == 0.0f) {
			zero_v3(pa->state.ave);
		}
		else {
			cross_v3_v3v3(pa->state.ave, pa->prev_state.vel, pa->state.vel);
			normalize_v3(pa->state.ave);
			angle = dot_v3v3(pa->prev_state.vel, pa->state.vel) / (len1 * len2);
			mul_v3_fl(pa->state.ave, saacos(angle) / dtime);
		}

		get_angular_velocity_vector(part->avemode, &pa->state, vec);
		axis_angle_to_quat(rot2, vec, dtime*part->avefac);
	}

	rotfac = len_v3(pa->state.ave);
	if (rotfac == 0.0f || (part->flag & PART_ROT_DYN)==0 || extrotfac == 0.0f) {
		unit_qt(rot1);
	}
	else {
		axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
	}
	mul_qt_qtqt(pa->state.rot,rot1,pa->prev_state.rot);
	mul_qt_qtqt(pa->state.rot,rot2,pa->state.rot);

	/* keep rotation quat in good health */
	normalize_qt(pa->state.rot);
}

/************************************************
 * Collisions
 *
 * The algorithm is roughly:
 *  1. Use a BVH tree to search for faces that a particle may collide with.
 *  2. Use Newton's method to find the exact time at which the collision occurs.
 *     https://en.wikipedia.org/wiki/Newton's_method
 *
 ************************************************/
#define COLLISION_MIN_RADIUS 0.001f
#define COLLISION_MIN_DISTANCE 0.0001f
#define COLLISION_ZERO 0.00001f
#define COLLISION_INIT_STEP 0.00008f
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 p0[3], e1[3], e2[3], d;

	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);

	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) {
		negate_v3(nor);
		d = -d;
	}

	return d - radius;
}
static float nr_distance_to_edge(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(nor))
{
	float v0[3], v1[3], v2[3], c[3];

	sub_v3_v3v3(v0, pce->x1, pce->x0);
	sub_v3_v3v3(v1, p, pce->x0);
	sub_v3_v3v3(v2, p, pce->x1);

	cross_v3_v3v3(c, v1, v2);

	return fabsf(len_v3(c)/len_v3(v0)) - radius;
}
static float nr_distance_to_vert(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(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 (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]);
		}
	}
}
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);

	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;
		}
		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];

			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)
				negate_v3(nor);

			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;
		}
	}
}
/* 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, dt_init, d0, d1, dd, n[3];
	int iter;

	pce->inv_nor = -1;

	if (col->inv_total_time > 0.0f) {
		/* Initial step size should be small, but not too small or floating point
		 * precision errors will appear. - z0r */
		dt_init = COLLISION_INIT_STEP * col->inv_total_time;
	}
	else {
		dt_init = 0.001f;
	}

	/* start from the beginning */
	t0 = 0.f;
	collision_interpolate_element(pce, t0, col->f, col);
	d0 = distance_func(col->co1, radius, pce, n);
	t1 = dt_init;
	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);

		/* 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;
		}

		/* Zero gradient (no movement relative to element). Can't step from
		 * here. */
		if (d1 == d0) {
			/* If first iteration, try from other end where the gradient may be
			 * greater. Note: code duplicated below. */
			if (iter == 0) {
				t0 = 1.f;
				collision_interpolate_element(pce, t0, col->f, col);
				d0 = distance_func(col->co2, radius, pce, n);
				t1 = 1.0f - dt_init;
				d1 = 0.f;
				continue;
			}
			else
				return -1.f;
		}

		dd = (t1-t0)/(d1-d0);

		t0 = t1;
		d0 = d1;

		t1 -= d1*dd;

		/* Particle moving away from plane could also mean a strangely rotating
		 * face, so check from end. Note: code duplicated above. */
		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 = 1.0f - dt_init;
			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 -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)
{
	ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
	ParticleCollisionElement *result = &col->pce;

	float ct;
	int i;

	for (i=0; i<3; i++) {
		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++) {
		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;
		}

	}

	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;
	const MVertTri *vt = &col->md->tri[index];
	MVert *x = col->md->x;
	MVert *v = col->md->current_v;
	float t = hit->dist/col->original_ray_length;
	int collision = 0;

	pce.x[0] = x[vt->tri[0]].co;
	pce.x[1] = x[vt->tri[1]].co;
	pce.x[2] = x[vt->tri[2]].co;

	pce.v[0] = v[vt->tri[0]].co;
	pce.v[1] = v[vt->tri[1]].co;
	pce.v[2] = v[vt->tri[2]].co;

	pce.tot = 3;
	pce.inside = 0;
	pce.index = index;

	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);
		collision += collision_sphere_to_verts(col, ray->radius, &pce, &t);
	}

	if (collision) {
		hit->dist = col->original_ray_length * t;
		hit->index = index;

		collision_point_velocity(&col->pce);

		col->hit = col->current;
	}
}
static int collision_detect(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, ListBase *colliders)
{
	const int raycast_flag = BVH_RAYCAST_DEFAULT & ~(BVH_RAYCAST_WATERTIGHT);
	ColliderCache *coll;
	float ray_dir[3];

	if (BLI_listbase_is_empty(colliders))
		return 0;

	sub_v3_v3v3(ray_dir, col->co2, col->co1);
	hit->index = -1;
	hit->dist = col->original_ray_length = normalize_v3(ray_dir);
	col->pce.inside = 0;

	/* 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; also skip if permeated */
		bool skip = false;

		for (int i = 0; i < col->skip_count; i++) {
			if (coll->ob == col->skip[i]) {
				skip = true;
				break;
			}
		}

		if (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_ex(
			        col->md->bvhtree, col->co1, ray_dir, col->radius, hit,
			        BKE_psys_collision_neartest_cb, col, raycast_flag);
		}
	}

	return hit->index >= 0;
}
static int collision_response(ParticleSimulationData *sim, ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, int kill, int dynamic_rotation)
{
	ParticleCollisionElement *pce = &col->pce;
	PartDeflect *pd = col->hit->pd;
	RNG *rng = sim->rng;
	/* point of collision */
	float co[3];
	/* location factor of collision between this iteration */
	float x = hit->dist/col->original_ray_length;
	/* time factor of collision between timestep */
	float f = col->f + x * (1.0f - col->f);
	/* time since previous collision (in seconds) */
	float dt1 = (f - col->f) * col->total_time;
	/* time left after collision (in seconds) */
	float dt2 = (1.0f - f) * col->total_time;
	/* did particle pass through the collision surface? */
	int through = (BLI_rng_get_float(rng) < pd->pdef_perm) ? 1 : 0;

	/* 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 {
		/* velocity directly before collision to be modified into velocity directly after collision */
		float v0[3];
		/* normal component of v0 */
		float v0_nor[3];
		/* tangential component of v0 */
		float v0_tan[3];
		/* tangential component of collision surface velocity */
		float vc_tan[3];
		float v0_dot, vc_dot;
		float damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_rng_get_float(rng) - 0.5f);
		float frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_rng_get_float(rng) - 0.5f);
		float distance, nor[3], dot;

		CLAMP(damp,0.0f, 1.0f);
		CLAMP(frict,0.0f, 1.0f);

		/* 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);

				/* 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) */

				/* rolling friction is around 0.01 of sliding friction
				 * (could be made a parameter) */
				mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);

				/* surface_velocity is opposite to cm velocity */
				negate_v3_v3(vr_tan, v1_tan);

				/* get back to global coordinates */
				add_v3_v3(v1_tan, vc_tan);

				/* convert to angular velocity*/
				cross_v3_v3v3(ave, vr_tan, pce->nor);
				mul_v3_fl(ave, 1.0f/MAX2(pa->size, 0.001f));

				/* 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);
			}
		}

		/* stickiness 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;
		}

		/* damping and flipping of velocity around normal */
		v0_dot *= 1.0f - damp;
		vc_dot *= through ? damp : 1.0f;

		/* 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 + v0_dot);
		else
			mul_v3_v3fl(v0_nor, pce->nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);

		/* combine components together again */
		add_v3_v3v3(v0, v0_nor, v0_tan);

		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;
			}
		}

		/* 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);

			dot = dot_v3v3(nor, v0);
			if (dot < 0.f)
				madd_v3_v3fl(v0, nor, -dot);

			distance = collision_point_distance_with_normal(pa->state.co, pce, 1.f, col, nor);

			if (distance < col->radius + COLLISION_MIN_DISTANCE)
				madd_v3_v3fl(pa->state.co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);

			dot = dot_v3v3(nor, pa->state.vel);
			if (dot < 0.f)
				madd_v3_v3fl(pa->state.vel, nor, -dot);
		}

		/* add stickiness to surface */
		madd_v3_v3fl(pa->state.vel, pce->nor, -pd->pdef_stickness);

		/* set coordinates for next iteration */
		copy_v3_v3(col->co1, co);
		copy_v3_v3(col->co2, pa->state.co);

		copy_v3_v3(col->ve1, v0);
		copy_v3_v3(col->ve2, pa->state.vel);

		col->f = f;
	}

	/* if permeability random roll succeeded, disable collider for this sim step */
	if (through) {
		col->skip[col->skip_count++] = col->hit;
	}

	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);

	copy_v3_v3(pa->state.vel, col->pce.vel);
	mul_v3_fl(pa->state.vel, col->inv_timestep);


	/* printf("max iterations\n"); */
}

/* 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)
{
	ParticleSettings *part = sim->psys->part;
	ParticleData *pa = sim->psys->particles + p;
	ParticleCollision col;
	BVHTreeRayHit hit;
	int collision_count=0;

	float timestep = psys_get_timestep(sim);

	memset(&col, 0, sizeof(ParticleCollision));

	col.total_time = timestep * dfra;
	col.inv_total_time = 1.0f/col.total_time;
	col.inv_timestep = 1.0f/timestep;

	col.cfra = cfra;
	col.old_cfra = sim->psys->cfra;

	/* 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);

	/* 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.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[col.skip_count++] = pa->boid->ground;
	}

	/* 10 iterations to catch multiple collisions */
	while (collision_count < PARTICLE_COLLISION_MAX_COLLISIONS) {
		if (collision_detect(pa, &col, &hit, sim->colliders)) {

			collision_count++;

			if (collision_count == PARTICLE_COLLISION_MAX_COLLISIONS)
				collision_fail(pa, &col);
			else if (collision_response(sim, pa, &col, &hit, part->flag & PART_DIE_ON_COL, part->flag & PART_ROT_DYN)==0)
				return;
		}
		else
			return;
	}
}
/************************************************/
/*			Hair								*/
/************************************************/
/* check if path cache or children need updating and do it if needed */
static void psys_update_path_cache(ParticleSimulationData *sim, float cfra, const bool use_render_params)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	ParticleEditSettings *pset = &sim->scene->toolsettings->particle;
	int distr=0, alloc=0, skip=0;

	if ((psys->part->childtype && psys->totchild != psys_get_tot_child(sim->scene, psys, use_render_params)) || psys->recalc&ID_RECALC_PSYS_RESET)
		alloc=1;

	if (alloc || psys->recalc&ID_RECALC_PSYS_CHILD || (psys->vgroup[PSYS_VG_DENSITY] && (sim->ob && sim->ob->mode & OB_MODE_WEIGHT_PAINT)))
		distr=1;

	if (distr) {
		if (alloc)
			realloc_particles(sim, sim->psys->totpart);

		if (psys_get_tot_child(sim->scene, psys, use_render_params)) {
			/* don't generate children while computing the hair keys */
			if (!(psys->part->type == PART_HAIR) || (psys->flag & PSYS_HAIR_DONE)) {
				distribute_particles(sim, PART_FROM_CHILD);

				if (part->childtype==PART_CHILD_FACES && part->parents != 0.0f)
					psys_find_parents(sim, use_render_params);
			}
		}
		else
			psys_free_children(psys);
	}

	if ((part->type==PART_HAIR || psys->flag&PSYS_KEYED || psys->pointcache->flag & PTCACHE_BAKED)==0)
		skip = 1; /* only hair, keyed and baked stuff can have paths */
	else if (part->ren_as != PART_DRAW_PATH && !(part->type==PART_HAIR && ELEM(part->ren_as, PART_DRAW_OB, PART_DRAW_GR)))
		skip = 1; /* particle visualization must be set as path */
	else if (DEG_get_mode(sim->depsgraph) != DAG_EVAL_RENDER) {
		if (part->draw_as != PART_DRAW_REND)
			skip = 1; /* draw visualization */
		else if (psys->pointcache->flag & PTCACHE_BAKING)
			skip = 1; /* no need to cache paths while baking dynamics */

		else if (psys_in_edit_mode(sim->depsgraph, psys)) {
			if ((pset->flag & PE_DRAW_PART)==0)
				skip = 1;
			else if (part->childtype==0 && (psys->flag & PSYS_HAIR_DYNAMICS && psys->pointcache->flag & PTCACHE_BAKED)==0)
				skip = 1; /* in edit mode paths are needed for child particles and dynamic hair */
		}
	}

	if (!skip) {
		psys_cache_paths(sim, cfra, use_render_params);

		/* for render, child particle paths are computed on the fly */
		if (part->childtype) {
			if (!psys->totchild)
				skip = 1;
			else if (psys->part->type == PART_HAIR && (psys->flag & PSYS_HAIR_DONE)==0)
				skip = 1;

			if (!skip)
				psys_cache_child_paths(sim, cfra, 0, use_render_params);
		}
	}
	else if (psys->pathcache)
		psys_free_path_cache(psys, NULL);
}

static bool psys_hair_use_simulation(ParticleData *pa, float max_length)
{
	/* Minimum segment length relative to average length.
	 * Hairs with segments below this length will be excluded from the simulation,
	 * because otherwise the solver will become unstable.
	 * The hair system should always make sure the hair segments have reasonable length ratios,
	 * but this can happen in old files when e.g. cutting hair.
	 */
	const float min_length = 0.1f * max_length;

	HairKey *key;
	int k;

	if (pa->totkey < 2)
		return false;

	for (k=1, key=pa->hair+1; k<pa->totkey; k++,key++) {
		float length = len_v3v3(key->co, (key-1)->co);
		if (length < min_length)
			return false;
	}

	return true;
}

static MDeformVert *hair_set_pinning(MDeformVert *dvert, float weight)
{
	if (dvert) {
		if (!dvert->totweight) {
			dvert->dw = MEM_callocN(sizeof(MDeformWeight), "deformWeight");
			dvert->totweight = 1;
		}

		dvert->dw->weight = weight;
		dvert++;
	}
	return dvert;
}

static void hair_create_input_mesh(ParticleSimulationData *sim, int totpoint, int totedge, Mesh **r_mesh, ClothHairData **r_hairdata)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	Mesh *mesh;
	ClothHairData *hairdata;
	MVert *mvert;
	MEdge *medge;
	MDeformVert *dvert;
	HairKey *key;
	PARTICLE_P;
	int k, hair_index;
	float hairmat[4][4];
	float max_length;
	float hair_radius;

	mesh = *r_mesh;
	if (!mesh) {
		*r_mesh = mesh = BKE_mesh_new_nomain(totpoint, totedge, 0, 0, 0);
		CustomData_add_layer(&mesh->vdata, CD_MDEFORMVERT, CD_CALLOC, NULL, mesh->totvert);
		BKE_mesh_update_customdata_pointers(mesh, false);
	}
	mvert = mesh->mvert;
	medge = mesh->medge;
	dvert = mesh->dvert;

	hairdata = *r_hairdata;
	if (!hairdata) {
		*r_hairdata = hairdata = MEM_mallocN(sizeof(ClothHairData) * totpoint, "hair data");
	}

	/* calculate maximum segment length */
	max_length = 0.0f;
	LOOP_PARTICLES {
		if (!(pa->flag & PARS_UNEXIST)) {
			for (k=1, key=pa->hair+1; k<pa->totkey; k++,key++) {
				float length = len_v3v3(key->co, (key-1)->co);
				if (max_length < length)
					max_length = length;
			}
		}
	}

	psys->clmd->sim_parms->vgroup_mass = 1;

	/* XXX placeholder for more flexible future hair settings */
	hair_radius = part->size;

	/* make vgroup for pin roots etc.. */
	hair_index = 1;
	LOOP_PARTICLES {
		if (!(pa->flag & PARS_UNEXIST)) {
			float root_mat[4][4];
			float bending_stiffness;
			bool use_hair;

			pa->hair_index = hair_index;
			use_hair = psys_hair_use_simulation(pa, max_length);

			psys_mat_hair_to_object(sim->ob, sim->psmd->mesh_final, psys->part->from, pa, hairmat);
			mul_m4_m4m4(root_mat, sim->ob->obmat, hairmat);
			normalize_m4(root_mat);

			bending_stiffness = CLAMPIS(1.0f - part->bending_random * psys_frand(psys, p + 666), 0.0f, 1.0f);

			for (k=0, key=pa->hair; k<pa->totkey; k++,key++) {
				ClothHairData *hair;
				float *co, *co_next;

				co = key->co;
				co_next = (key+1)->co;

				/* create fake root before actual root to resist bending */
				if (k==0) {
					hair = &psys->clmd->hairdata[pa->hair_index - 1];
					copy_v3_v3(hair->loc, root_mat[3]);
					copy_m3_m4(hair->rot, root_mat);

					hair->radius = hair_radius;
					hair->bending_stiffness = bending_stiffness;

					add_v3_v3v3(mvert->co, co, co);
					sub_v3_v3(mvert->co, co_next);
					mul_m4_v3(hairmat, mvert->co);

					medge->v1 = pa->hair_index - 1;
					medge->v2 = pa->hair_index;

					dvert = hair_set_pinning(dvert, 1.0f);

					mvert++;
					medge++;
				}

				/* store root transform in cloth data */
				hair = &psys->clmd->hairdata[pa->hair_index + k];
				copy_v3_v3(hair->loc, root_mat[3]);
				copy_m3_m4(hair->rot, root_mat);

				hair->radius = hair_radius;
				hair->bending_stiffness = bending_stiffness;

				copy_v3_v3(mvert->co, co);
				mul_m4_v3(hairmat, mvert->co);

				if (k) {
					medge->v1 = pa->hair_index + k - 1;
					medge->v2 = pa->hair_index + k;
				}

				/* roots and disabled hairs should be 1.0, the rest can be anything from 0.0 to 1.0 */
				if (use_hair)
					dvert = hair_set_pinning(dvert, key->weight);
				else
					dvert = hair_set_pinning(dvert, 1.0f);

				mvert++;
				if (k)
					medge++;
			}

			hair_index += pa->totkey + 1;
		}
	}
}

static void do_hair_dynamics(ParticleSimulationData *sim)
{
	ParticleSystem *psys = sim->psys;
	PARTICLE_P;
	EffectorWeights *clmd_effweights;
	int totpoint;
	int totedge;
	float (*deformedVerts)[3];
	bool realloc_roots;

	if (!psys->clmd) {
		psys->clmd = (ClothModifierData*)modifier_new(eModifierType_Cloth);
		psys->clmd->sim_parms->goalspring = 0.0f;
		psys->clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_RESIST_SPRING_COMPRESS;
		psys->clmd->coll_parms->flags &= ~CLOTH_COLLSETTINGS_FLAG_SELF;
	}

	/* count simulated points */
	totpoint = 0;
	totedge = 0;
	LOOP_PARTICLES {
		if (!(pa->flag & PARS_UNEXIST)) {
			/* "out" dm contains all hairs */
			totedge += pa->totkey;
			totpoint += pa->totkey + 1; /* +1 for virtual root point */
		}
	}

	/* whether hair root info array has to be reallocated */
	realloc_roots = false;
	if (psys->hair_in_mesh) {
		Mesh *mesh = psys->hair_in_mesh;
		if (totpoint != mesh->totvert || totedge != mesh->totedge) {
			BKE_id_free(NULL, mesh);
			psys->hair_in_mesh = NULL;
			realloc_roots = true;
		}
	}

	if (!psys->hair_in_mesh || !psys->clmd->hairdata || realloc_roots) {
		if (psys->clmd->hairdata) {
			MEM_freeN(psys->clmd->hairdata);
			psys->clmd->hairdata = NULL;
		}
	}

	hair_create_input_mesh(sim, totpoint, totedge, &psys->hair_in_mesh, &psys->clmd->hairdata);

	if (psys->hair_out_mesh)
		BKE_id_free(NULL, psys->hair_out_mesh);

	psys->clmd->point_cache = psys->pointcache;
	/* for hair sim we replace the internal cloth effector weights temporarily
	 * to use the particle settings
	 */
	clmd_effweights = psys->clmd->sim_parms->effector_weights;
	psys->clmd->sim_parms->effector_weights = psys->part->effector_weights;

	BKE_id_copy_ex(NULL, &psys->hair_in_mesh->id, (ID **)&psys->hair_out_mesh, LIB_ID_COPY_LOCALIZE);
	deformedVerts = BKE_mesh_vertexCos_get(psys->hair_out_mesh, NULL);
	clothModifier_do(psys->clmd, sim->depsgraph, sim->scene, sim->ob, psys->hair_in_mesh, deformedVerts);
	BKE_mesh_apply_vert_coords(psys->hair_out_mesh, deformedVerts);

	MEM_freeN(deformedVerts);

	/* restore cloth effector weights */
	psys->clmd->sim_parms->effector_weights = clmd_effweights;
}
static void hair_step(ParticleSimulationData *sim, float cfra, const bool use_render_params)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	PARTICLE_P;
	float disp = psys_get_current_display_percentage(psys, use_render_params);

	LOOP_PARTICLES {
		pa->size = part->size;
		if (part->randsize > 0.0f)
			pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);

		if (psys_frand(psys, p) > disp)
			pa->flag |= PARS_NO_DISP;
		else
			pa->flag &= ~PARS_NO_DISP;
	}

	if (psys->recalc & ID_RECALC_PSYS_RESET) {
		/* need this for changing subsurf levels */
		psys_calc_dmcache(sim->ob, sim->psmd->mesh_final, sim->psmd->mesh_original, psys);

		if (psys->clmd)
			cloth_free_modifier(psys->clmd);
	}

	/* dynamics with cloth simulation, psys->particles can be NULL with 0 particles [#25519] */
	if (psys->part->type==PART_HAIR && psys->flag & PSYS_HAIR_DYNAMICS && psys->particles)
		do_hair_dynamics(sim);

	/* following lines were removed r29079 but cause bug [#22811], see report for details */
	psys_update_effectors(sim);
	psys_update_path_cache(sim, cfra, use_render_params);

	psys->flag |= PSYS_HAIR_UPDATED;
}

static void save_hair(ParticleSimulationData *sim, float UNUSED(cfra))
{
	Object *ob = sim->ob;
	ParticleSystem *psys = sim->psys;
	HairKey *key, *root;
	PARTICLE_P;

	invert_m4_m4(ob->imat, ob->obmat);

	psys->lattice_deform_data= psys_create_lattice_deform_data(sim);

	if (psys->totpart==0) return;

	/* save new keys for elements if needed */
	LOOP_PARTICLES {
		/* first time alloc */
		if (pa->totkey==0 || pa->hair==NULL) {
			pa->hair = MEM_callocN((psys->part->hair_step + 1) * sizeof(HairKey), "HairKeys");
			pa->totkey = 0;
		}

		key = root = pa->hair;
		key += pa->totkey;

		/* convert from global to geometry space */
		copy_v3_v3(key->co, pa->state.co);
		mul_m4_v3(ob->imat, key->co);

		if (pa->totkey) {
			sub_v3_v3(key->co, root->co);
			psys_vec_rot_to_face(sim->psmd->mesh_final, pa, key->co);
		}

		key->time = pa->state.time;

		key->weight = 1.0f - key->time / 100.0f;

		pa->totkey++;

		/* root is always in the origin of hair space so we set it to be so after the last key is saved*/
		if (pa->totkey == psys->part->hair_step + 1) {
			zero_v3(root->co);
		}

	}
}

/* Code for an adaptive time step based on the Courant-Friedrichs-Lewy
 * condition. */
static const float MIN_TIMESTEP = 1.0f / 101.0f;
/* Tolerance of 1.5 means the last subframe neither favors growing nor
 * shrinking (e.g if it were 1.3, the last subframe would tend to be too
 * small). */
static const float TIMESTEP_EXPANSION_FACTOR = 0.1f;
static const float TIMESTEP_EXPANSION_TOLERANCE = 1.5f;

/* Calculate the speed of the particle relative to the local scale of the
 * simulation. This should be called once per particle during a simulation
 * step, after the velocity has been updated. element_size defines the scale of
 * the simulation, and is typically the distance to neighboring particles. */
static void update_courant_num(ParticleSimulationData *sim, ParticleData *pa,
                               float dtime, SPHData *sphdata, SpinLock *spin)
{
	float relative_vel[3];

	sub_v3_v3v3(relative_vel, pa->prev_state.vel, sphdata->flow);

	const float courant_num = len_v3(relative_vel) * dtime / sphdata->element_size;
	if (sim->courant_num < courant_num) {
		BLI_spin_lock(spin);
		if (sim->courant_num < courant_num) {
			sim->courant_num = courant_num;
		}
		BLI_spin_unlock(spin);
	}
}
static float get_base_time_step(ParticleSettings *part)
{
	return 1.0f / (float) (part->subframes + 1);
}
/* Update time step size to suit current conditions. */
static void update_timestep(ParticleSystem *psys, ParticleSimulationData *sim)
{
	float dt_target;
	if (sim->courant_num == 0.0f)
		dt_target = 1.0f;
	else
		dt_target = psys->dt_frac * (psys->part->courant_target / sim->courant_num);

	/* Make sure the time step is reasonable. For some reason, the CLAMP macro
	 * doesn't work here. The time step becomes too large. - z0r */
	if (dt_target < MIN_TIMESTEP)
		dt_target = MIN_TIMESTEP;
	else if (dt_target > get_base_time_step(psys->part))
		dt_target = get_base_time_step(psys->part);

	/* Decrease time step instantly, but increase slowly. */
	if (dt_target > psys->dt_frac)
		psys->dt_frac = interpf(dt_target, psys->dt_frac, TIMESTEP_EXPANSION_FACTOR);
	else
		psys->dt_frac = dt_target;
}

static float sync_timestep(ParticleSystem *psys, float t_frac)
{
	/* Sync with frame end if it's close. */
	if (t_frac == 1.0f)
		return psys->dt_frac;
	else if (t_frac + (psys->dt_frac * TIMESTEP_EXPANSION_TOLERANCE) >= 1.0f)
		return 1.0f - t_frac;
	else
		return psys->dt_frac;
}

/************************************************/
/*			System Core							*/
/************************************************/

typedef struct DynamicStepSolverTaskData {
	ParticleSimulationData *sim;

	float cfra;
	float timestep;
	float dtime;

	SpinLock spin;
} DynamicStepSolverTaskData;

static void dynamics_step_sph_ddr_task_cb_ex(
        void *__restrict userdata,
        const int p,
        const ParallelRangeTLS *__restrict tls)
{
	DynamicStepSolverTaskData *data = userdata;
	ParticleSimulationData *sim = data->sim;
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;

	SPHData *sphdata = tls->userdata_chunk;

	ParticleData *pa;

	if ((pa = psys->particles + p)->state.time <= 0.0f) {
		return;
	}

	/* do global forces & effectors */
	basic_integrate(sim, p, pa->state.time, data->cfra);

	/* actual fluids calculations */
	sph_integrate(sim, pa, pa->state.time, sphdata);

	if (sim->colliders)
		collision_check(sim, p, pa->state.time, data->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, data->timestep);

	if (part->time_flag & PART_TIME_AUTOSF) {
		update_courant_num(sim, pa, data->dtime, sphdata, &data->spin);
	}
}

static void dynamics_step_sph_classical_basic_integrate_task_cb_ex(
        void *__restrict userdata,
        const int p,
        const ParallelRangeTLS *__restrict UNUSED(tls))
{
	DynamicStepSolverTaskData *data = userdata;
	ParticleSimulationData *sim = data->sim;
	ParticleSystem *psys = sim->psys;

	ParticleData *pa;

	if ((pa = psys->particles + p)->state.time <= 0.0f) {
		return;
	}

	basic_integrate(sim, p, pa->state.time, data->cfra);
}

static void dynamics_step_sph_classical_calc_density_task_cb_ex(
        void *__restrict userdata,
        const int p,
        const ParallelRangeTLS *__restrict tls)
{
	DynamicStepSolverTaskData *data = userdata;
	ParticleSimulationData *sim = data->sim;
	ParticleSystem *psys = sim->psys;

	SPHData *sphdata = tls->userdata_chunk;

	ParticleData *pa;

	if ((pa = psys->particles + p)->state.time <= 0.0f) {
		return;
	}

	sphclassical_calc_dens(pa, pa->state.time, sphdata);
}

static void dynamics_step_sph_classical_integrate_task_cb_ex(
        void *__restrict userdata,
        const int p,
        const ParallelRangeTLS *__restrict tls)
{
	DynamicStepSolverTaskData *data = userdata;
	ParticleSimulationData *sim = data->sim;
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;

	SPHData *sphdata = tls->userdata_chunk;

	ParticleData *pa;

	if ((pa = psys->particles + p)->state.time <= 0.0f) {
		return;
	}

	/* actual fluids calculations */
	sph_integrate(sim, pa, pa->state.time, sphdata);

	if (sim->colliders)
		collision_check(sim, p, pa->state.time, data->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, data->timestep);

	if (part->time_flag & PART_TIME_AUTOSF) {
		update_courant_num(sim, pa, data->dtime, sphdata, &data->spin);
	}
}

/* unbaked particles are calculated dynamically */
static void dynamics_step(ParticleSimulationData *sim, float cfra)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part=psys->part;
	BoidBrainData bbd;
	ParticleTexture ptex;
	PARTICLE_P;
	float timestep;
	/* frame & time changes */
	float dfra, dtime;
	float birthtime, dietime;

	/* where have we gone in time since last time */
	dfra= cfra - psys->cfra;

	timestep = psys_get_timestep(sim);
	dtime= dfra*timestep;

	if (dfra < 0.0f) {
		LOOP_EXISTING_PARTICLES {
			psys_get_texture(sim, pa, &ptex, PAMAP_SIZE, cfra);
			pa->size = part->size*ptex.size;
			if (part->randsize > 0.0f)
				pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);

			reset_particle(sim, pa, dtime, cfra);
		}
		return;
	}

	/* for now do both, boids us 'rng' */
	sim->rng = BLI_rng_new_srandom(31415926 + (int)cfra + psys->seed);

	psys_update_effectors(sim);

	if (part->type != PART_HAIR)
		sim->colliders = BKE_collider_cache_create(sim->depsgraph, sim->ob, part->collision_group);

	/* initialize physics type specific stuff */
	switch (part->phystype) {
		case PART_PHYS_BOIDS:
		{
			ParticleTarget *pt = psys->targets.first;
			bbd.sim = sim;
			bbd.part = part;
			bbd.cfra = cfra;
			bbd.dfra = dfra;
			bbd.timestep = timestep;
			bbd.rng = sim->rng;

			psys_update_particle_tree(psys, cfra);

			boids_precalc_rules(part, cfra);

			for (; pt; pt=pt->next) {
				ParticleSystem *psys_target = psys_get_target_system(sim->ob, pt);
				if (psys_target && psys_target != psys) {
					psys_update_particle_tree(psys_target, cfra);
				}
			}
			break;
		}
		case PART_PHYS_FLUID:
		{
			ParticleTarget *pt = psys->targets.first;
			psys_update_particle_bvhtree(psys, cfra);

			for (; pt; pt=pt->next) {  /* Updating others systems particle tree for fluid-fluid interaction */
				if (pt->ob)
					psys_update_particle_bvhtree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
			}
			break;
		}
	}
	/* initialize all particles for dynamics */
	LOOP_SHOWN_PARTICLES {
		copy_particle_key(&pa->prev_state,&pa->state,1);

		psys_get_texture(sim, pa, &ptex, PAMAP_SIZE, cfra);

		pa->size = part->size*ptex.size;
		if (part->randsize > 0.0f)
			pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);

		birthtime = pa->time;
		dietime = pa->dietime;

		/* store this, so we can do multiple loops over particles */
		pa->state.time = dfra;

		if (dietime <= cfra && psys->cfra < dietime) {
			/* particle dies some time between this and last step */
			pa->state.time = dietime - ((birthtime > psys->cfra) ? birthtime : psys->cfra);
			pa->alive = PARS_DYING;
		}
		else if (birthtime <= cfra && birthtime >= psys->cfra) {
			/* particle is born some time between this and last step*/
			reset_particle(sim, pa, dfra*timestep, cfra);
			pa->alive = PARS_ALIVE;
			pa->state.time = cfra - birthtime;
		}
		else if (dietime < cfra) {
			/* nothing to be done when particle is dead */
		}

		/* only reset unborn particles if they're shown or if the particle is born soon*/
		if (pa->alive==PARS_UNBORN && (part->flag & PART_UNBORN || (cfra + psys->pointcache->step > pa->time))) {
			reset_particle(sim, pa, dtime, cfra);
		}
		else if (part->phystype == PART_PHYS_NO) {
			reset_particle(sim, pa, dtime, cfra);
		}

		if (ELEM(pa->alive, PARS_ALIVE, PARS_DYING)==0 || (pa->flag & (PARS_UNEXIST|PARS_NO_DISP)))
			pa->state.time = -1.f;
	}

	switch (part->phystype) {
		case PART_PHYS_NEWTON:
		{
			LOOP_DYNAMIC_PARTICLES {
				/* do global forces & effectors */
				basic_integrate(sim, p, pa->state.time, cfra);

				/* deflection */
				if (sim->colliders)
					collision_check(sim, p, pa->state.time, cfra);

				/* rotations */
				basic_rotate(part, pa, pa->state.time, timestep);
			}
			break;
		}
		case PART_PHYS_BOIDS:
		{
			LOOP_DYNAMIC_PARTICLES {
				bbd.goal_ob = NULL;

				boid_brain(&bbd, p, pa);

				if (pa->alive != PARS_DYING) {
					boid_body(&bbd, pa);

					/* deflection */
					if (sim->colliders)
						collision_check(sim, p, pa->state.time, cfra);
				}
			}
			break;
		}
		case PART_PHYS_FLUID:
		{
			SPHData sphdata;
			psys_sph_init(sim, &sphdata);

			DynamicStepSolverTaskData task_data = {
			    .sim = sim, .cfra = cfra, .timestep = timestep, .dtime = dtime,
			};

			BLI_spin_init(&task_data.spin);

			if (part->fluid->solver == SPH_SOLVER_DDR) {
				/* Apply SPH forces using double-density relaxation algorithm
				 * (Clavat et. al.) */

				ParallelRangeSettings settings;
				BLI_parallel_range_settings_defaults(&settings);
				settings.use_threading = (psys->totpart > 100);
				settings.userdata_chunk = &sphdata;
				settings.userdata_chunk_size = sizeof(sphdata);
				BLI_task_parallel_range(
				        0, psys->totpart,
				        &task_data,
				        dynamics_step_sph_ddr_task_cb_ex,
				        &settings);

				sph_springs_modify(psys, timestep);
			}
			else {
				/* SPH_SOLVER_CLASSICAL */
				/* Apply SPH forces using classical algorithm (due to Gingold
				 * and Monaghan). Note that, unlike double-density relaxation,
				 * this algorithm is separated into distinct loops. */

				{
					ParallelRangeSettings settings;
					BLI_parallel_range_settings_defaults(&settings);
					settings.use_threading = (psys->totpart > 100);
					BLI_task_parallel_range(
					        0, psys->totpart,
					        &task_data,
					        dynamics_step_sph_classical_basic_integrate_task_cb_ex,
					        &settings);
				}

				/* calculate summation density */
				/* Note that we could avoid copying sphdata for each thread here (it's only read here),
				 * but doubt this would gain us anything except confusion... */
				{
					ParallelRangeSettings settings;
					BLI_parallel_range_settings_defaults(&settings);
					settings.use_threading = (psys->totpart > 100);
					settings.userdata_chunk = &sphdata;
					settings.userdata_chunk_size = sizeof(sphdata);
					BLI_task_parallel_range(
					        0, psys->totpart,
					        &task_data,
					        dynamics_step_sph_classical_calc_density_task_cb_ex,
					        &settings);
				}

				/* do global forces & effectors */
				{
					ParallelRangeSettings settings;
					BLI_parallel_range_settings_defaults(&settings);
					settings.use_threading = (psys->totpart > 100);
					settings.userdata_chunk = &sphdata;
					settings.userdata_chunk_size = sizeof(sphdata);
					BLI_task_parallel_range(
					        0, psys->totpart,
					        &task_data,
					        dynamics_step_sph_classical_integrate_task_cb_ex,
					        &settings);
				}
			}

			BLI_spin_end(&task_data.spin);

			psys_sph_finalise(&sphdata);
			break;
		}
	}

	/* finalize particle state and time after dynamics */
	LOOP_DYNAMIC_PARTICLES {
		if (pa->alive == PARS_DYING) {
			pa->alive=PARS_DEAD;
			pa->state.time=pa->dietime;
		}
		else
			pa->state.time=cfra;
	}

	BKE_collider_cache_free(&sim->colliders);
	BLI_rng_free(sim->rng);
	sim->rng = NULL;
}

static void update_children(ParticleSimulationData *sim, const bool use_render_params)
{
	if ((sim->psys->part->type == PART_HAIR) && (sim->psys->flag & PSYS_HAIR_DONE)==0)
	/* don't generate children while growing hair - waste of time */
		psys_free_children(sim->psys);
	else if (sim->psys->part->childtype) {
		if (sim->psys->totchild != psys_get_tot_child(sim->scene, sim->psys, use_render_params))
			distribute_particles(sim, PART_FROM_CHILD);
		else {
			/* Children are up to date, nothing to do. */
		}
	}
	else
		psys_free_children(sim->psys);
}
/* updates cached particles' alive & other flags etc..*/
static void cached_step(ParticleSimulationData *sim, float cfra, const bool use_render_params)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	ParticleTexture ptex;
	PARTICLE_P;
	float disp, dietime;

	psys_update_effectors(sim);

	disp= psys_get_current_display_percentage(psys, use_render_params);

	LOOP_PARTICLES {
		psys_get_texture(sim, pa, &ptex, PAMAP_SIZE, cfra);
		pa->size = part->size*ptex.size;
		if (part->randsize > 0.0f)
			pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);

		psys->lattice_deform_data = psys_create_lattice_deform_data(sim);

		dietime = pa->dietime;

		/* update alive status and push events */
		if (pa->time > cfra) {
			pa->alive = PARS_UNBORN;
			if (part->flag & PART_UNBORN && (psys->pointcache->flag & PTCACHE_EXTERNAL) == 0)
				reset_particle(sim, pa, 0.0f, cfra);
		}
		else if (dietime <= cfra)
			pa->alive = PARS_DEAD;
		else
			pa->alive = PARS_ALIVE;

		if (psys->lattice_deform_data) {
			end_latt_deform(psys->lattice_deform_data);
			psys->lattice_deform_data = NULL;
		}

		if (psys_frand(psys, p) > disp)
			pa->flag |= PARS_NO_DISP;
		else
			pa->flag &= ~PARS_NO_DISP;
	}
}

static void particles_fluid_step(
        ParticleSimulationData *sim, int UNUSED(cfra), const bool use_render_params)
{
	ParticleSystem *psys = sim->psys;
	if (psys->particles) {
		MEM_freeN(psys->particles);
		psys->particles = 0;
		psys->totpart = 0;
	}

	/* fluid sim particle import handling, actual loading of particles from file */
#ifdef WITH_MOD_FLUID
	{
		FluidsimModifierData *fluidmd = (FluidsimModifierData *)modifiers_findByType(sim->ob, eModifierType_Fluidsim);

		if ( fluidmd && fluidmd->fss) {
			FluidsimSettings *fss= fluidmd->fss;
			ParticleSettings *part = psys->part;
			ParticleData *pa=NULL;
			char filename[256];
			char debugStrBuffer[256];
			int  curFrame = sim->scene->r.cfra -1; // warning - sync with derived mesh fsmesh loading
			int  p, j, totpart;
			int readMask, activeParts = 0, fileParts = 0;
			gzFile gzf;

// XXX			if (ob==G.obedit) // off...
//				return;

			// ok, start loading
			BLI_join_dirfile(filename, sizeof(filename), fss->surfdataPath, OB_FLUIDSIM_SURF_PARTICLES_FNAME);

			BLI_path_abs(filename, modifier_path_relbase_from_global(sim->ob));

			BLI_path_frame(filename, curFrame, 0); // fixed #frame-no

			gzf = BLI_gzopen(filename, "rb");
			if (!gzf) {
				BLI_snprintf(debugStrBuffer, sizeof(debugStrBuffer),"readFsPartData::error - Unable to open file for reading '%s'\n", filename);
				// XXX bad level call elbeemDebugOut(debugStrBuffer);
				return;
			}

			gzread(gzf, &totpart, sizeof(totpart));
			totpart = (use_render_params) ? totpart:(part->disp*totpart) / 100;

			part->totpart= totpart;
			part->sta=part->end = 1.0f;
			part->lifetime = sim->scene->r.efra + 1;

			/* allocate particles */
			realloc_particles(sim, part->totpart);

			// set up reading mask
			readMask = fss->typeFlags;

			for (p=0, pa=psys->particles; p<totpart; p++, pa++) {
				int ptype=0;

				gzread(gzf, &ptype, sizeof( ptype ));
				if (ptype & readMask) {
					activeParts++;

					gzread(gzf, &(pa->size), sizeof(float));

					pa->size /= 10.0f;

					for (j=0; j<3; j++) {
						float wrf;
						gzread(gzf, &wrf, sizeof( wrf ));
						pa->state.co[j] = wrf;
						//fprintf(stderr,"Rj%d ",j);
					}
					for (j=0; j<3; j++) {
						float wrf;
						gzread(gzf, &wrf, sizeof( wrf ));
						pa->state.vel[j] = wrf;
					}

					zero_v3(pa->state.ave);
					unit_qt(pa->state.rot);

					pa->time = 1.f;
					pa->dietime = sim->scene->r.efra + 1;
					pa->lifetime = sim->scene->r.efra;
					pa->alive = PARS_ALIVE;
					//if (a < 25) fprintf(stderr,"FSPARTICLE debug set %s, a%d = %f,%f,%f, life=%f\n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
				}
				else {
					// skip...
					for (j=0; j<2*3+1; j++) {
						float wrf; gzread(gzf, &wrf, sizeof( wrf ));
					}
				}
				fileParts++;
			}
			gzclose(gzf);

			totpart = psys->totpart = activeParts;
			BLI_snprintf(debugStrBuffer,sizeof(debugStrBuffer),"readFsPartData::done - particles:%d, active:%d, file:%d, mask:%d\n", psys->totpart,activeParts,fileParts,readMask);
			// bad level call
			// XXX elbeemDebugOut(debugStrBuffer);

		} // fluid sim particles done
	}
#else
	UNUSED_VARS(use_render_params);
#endif // WITH_MOD_FLUID
}

static int emit_particles(ParticleSimulationData *sim, PTCacheID *pid, float UNUSED(cfra))
{
	ParticleSystem *psys = sim->psys;
	int oldtotpart = psys->totpart;
	int totpart = tot_particles(psys, pid);

	if (totpart != oldtotpart)
		realloc_particles(sim, totpart);

	return totpart - oldtotpart;
}

/* Calculates the next state for all particles of the system
 * In particles code most fra-ending are frames, time-ending are fra*timestep (seconds)
 * 1. Emit particles
 * 2. Check cache (if used) and return if frame is cached
 * 3. Do dynamics
 * 4. Save to cache */
static void system_step(ParticleSimulationData *sim, float cfra, const bool use_render_params)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	PointCache *cache = psys->pointcache;
	PTCacheID ptcacheid, *pid = NULL;
	PARTICLE_P;
	float disp, cache_cfra = cfra; /*, *vg_vel= 0, *vg_tan= 0, *vg_rot= 0, *vg_size= 0; */
	int startframe = 0, endframe = 100, oldtotpart = 0;

	/* cache shouldn't be used for hair or "continue physics" */
	if (part->type != PART_HAIR) {
		psys_clear_temp_pointcache(psys);

		/* set suitable cache range automatically */
		if ((cache->flag & (PTCACHE_BAKING|PTCACHE_BAKED))==0)
			psys_get_pointcache_start_end(sim->scene, psys, &cache->startframe, &cache->endframe);

		pid = &ptcacheid;
		BKE_ptcache_id_from_particles(pid, sim->ob, psys);

		BKE_ptcache_id_time(pid, sim->scene, 0.0f, &startframe, &endframe, NULL);

		/* clear everything on start frame, or when psys needs full reset! */
		if ((cfra == startframe) || (psys->recalc & ID_RECALC_PSYS_RESET)) {
			BKE_ptcache_id_reset(sim->scene, pid, PTCACHE_RESET_OUTDATED);
			BKE_ptcache_validate(cache, startframe);
			cache->flag &= ~PTCACHE_REDO_NEEDED;
		}

		CLAMP(cache_cfra, startframe, endframe);
	}

/* 1. emit particles and redo particles if needed */
	oldtotpart = psys->totpart;
	if (emit_particles(sim, pid, cfra) || psys->recalc & ID_RECALC_PSYS_RESET) {
		distribute_particles(sim, part->from);
		initialize_all_particles(sim);
		/* reset only just created particles (on startframe all particles are recreated) */
		reset_all_particles(sim, 0.0, cfra, oldtotpart);
		free_unexisting_particles(sim);

		if (psys->fluid_springs) {
			MEM_freeN(psys->fluid_springs);
			psys->fluid_springs = NULL;
		}

		psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;

		/* flag for possible explode modifiers after this system */
		sim->psmd->flag |= eParticleSystemFlag_Pars;

		BKE_ptcache_id_clear(pid, PTCACHE_CLEAR_AFTER, cfra);
	}

/* 2. try to read from the cache */
	if (pid) {
		int cache_result = BKE_ptcache_read(pid, cache_cfra, true);

		if (ELEM(cache_result, PTCACHE_READ_EXACT, PTCACHE_READ_INTERPOLATED)) {
			cached_step(sim, cfra, use_render_params);
			update_children(sim, use_render_params);
			psys_update_path_cache(sim, cfra, use_render_params);

			BKE_ptcache_validate(cache, (int)cache_cfra);

			if (cache_result == PTCACHE_READ_INTERPOLATED && cache->flag & PTCACHE_REDO_NEEDED)
				BKE_ptcache_write(pid, (int)cache_cfra);

			return;
		}
		/* Cache is supposed to be baked, but no data was found so bail out */
		else if (cache->flag & PTCACHE_BAKED) {
			psys_reset(psys, PSYS_RESET_CACHE_MISS);
			return;
		}
		else if (cache_result == PTCACHE_READ_OLD) {
			psys->cfra = (float)cache->simframe;
			cached_step(sim, psys->cfra, use_render_params);
		}

		/* if on second frame, write cache for first frame */
		if (psys->cfra == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact==0))
			BKE_ptcache_write(pid, startframe);
	}
	else
		BKE_ptcache_invalidate(cache);

/* 3. do dynamics */
	/* set particles to be not calculated TODO: can't work with pointcache */
	disp= psys_get_current_display_percentage(psys, use_render_params);

	LOOP_PARTICLES {
		if (psys_frand(psys, p) > disp)
			pa->flag |= PARS_NO_DISP;
		else
			pa->flag &= ~PARS_NO_DISP;
	}

	if (psys->totpart) {
		int dframe, totframesback = 0;
		float t_frac, dt_frac;

		/* handle negative frame start at the first frame by doing
		 * all the steps before the first frame */
		if ((int)cfra == startframe && part->sta < startframe)
			totframesback = (startframe - (int)part->sta);

		if (!(part->time_flag & PART_TIME_AUTOSF)) {
			/* Constant time step */
			psys->dt_frac = get_base_time_step(part);
		}
		else if ((int)cfra == startframe) {
			/* Variable time step; initialise to subframes */
			psys->dt_frac = get_base_time_step(part);
		}
		else if (psys->dt_frac < MIN_TIMESTEP) {
			/* Variable time step; subsequent frames */
			psys->dt_frac = MIN_TIMESTEP;
		}

		for (dframe=-totframesback; dframe<=0; dframe++) {
			/* simulate each subframe */
			dt_frac = psys->dt_frac;
			for (t_frac = dt_frac; t_frac <= 1.0f; t_frac += dt_frac) {
				sim->courant_num = 0.0f;
				dynamics_step(sim, cfra+dframe+t_frac - 1.f);
				psys->cfra = cfra+dframe+t_frac - 1.f;

				if (part->time_flag & PART_TIME_AUTOSF)
					update_timestep(psys, sim);
				/* Even without AUTOSF dt_frac may not add up to 1.0 due to float precision. */
				dt_frac = sync_timestep(psys, t_frac);
			}
		}
	}

/* 4. only write cache starting from second frame */
	if (pid) {
		BKE_ptcache_validate(cache, (int)cache_cfra);
		if ((int)cache_cfra != startframe)
			BKE_ptcache_write(pid, (int)cache_cfra);
	}

	update_children(sim, use_render_params);

/* cleanup */
	if (psys->lattice_deform_data) {
		end_latt_deform(psys->lattice_deform_data);
		psys->lattice_deform_data = NULL;
	}
}

/* system type has changed so set sensible defaults and clear non applicable flags */
void psys_changed_type(Object *ob, ParticleSystem *psys)
{
	ParticleSettings *part = psys->part;
	PTCacheID pid;

	BKE_ptcache_id_from_particles(&pid, ob, psys);

	if (part->phystype != PART_PHYS_KEYED)
		psys->flag &= ~PSYS_KEYED;

	if (part->type == PART_HAIR) {
		if (ELEM(part->ren_as, PART_DRAW_NOT, PART_DRAW_PATH, PART_DRAW_OB, PART_DRAW_GR)==0)
			part->ren_as = PART_DRAW_PATH;

		if (part->distr == PART_DISTR_GRID)
			part->distr = PART_DISTR_JIT;

		if (ELEM(part->draw_as, PART_DRAW_NOT, PART_DRAW_REND, PART_DRAW_PATH)==0)
			part->draw_as = PART_DRAW_REND;

		CLAMP(part->path_start, 0.0f, 100.0f);
		CLAMP(part->path_end, 0.0f, 100.0f);

		BKE_ptcache_id_clear(&pid, PTCACHE_CLEAR_ALL, 0);
	}
	else {
		free_hair(ob, psys, 1);

		CLAMP(part->path_start, 0.0f, MAX2(100.0f, part->end + part->lifetime));
		CLAMP(part->path_end, 0.0f, MAX2(100.0f, part->end + part->lifetime));
	}

	psys_reset(psys, PSYS_RESET_ALL);
}
void psys_check_boid_data(ParticleSystem *psys)
{
		BoidParticle *bpa;
		PARTICLE_P;

		pa = psys->particles;

		if (!pa)
			return;

		if (psys->part && psys->part->phystype==PART_PHYS_BOIDS) {
			if (!pa->boid) {
				bpa = MEM_callocN(psys->totpart * sizeof(BoidParticle), "Boid Data");

				LOOP_PARTICLES {
					pa->boid = bpa++;
				}
			}
		}
		else if (pa->boid) {
			MEM_freeN(pa->boid);
			LOOP_PARTICLES {
				pa->boid = NULL;
			}
		}
}

void BKE_particlesettings_fluid_default_settings(ParticleSettings *part)
{
	SPHFluidSettings *fluid = part->fluid;

	fluid->spring_k = 0.f;
	fluid->plasticity_constant = 0.1f;
	fluid->yield_ratio = 0.1f;
	fluid->rest_length = 1.f;
	fluid->viscosity_omega = 2.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)
{
	ParticleSettings *part = sim->psys->part;

	if (ELEM(part->phystype, PART_PHYS_NO, PART_PHYS_KEYED)) {
		PTCacheID pid;
		BKE_ptcache_id_from_particles(&pid, sim->ob, sim->psys);
		BKE_ptcache_id_clear(&pid, PTCACHE_CLEAR_ALL, 0);
	}
	else {
		free_keyed_keys(sim->psys);
		sim->psys->flag &= ~PSYS_KEYED;
	}

	/* RNA Update must ensure this is true. */
	if (part->phystype == PART_PHYS_BOIDS) {
		BLI_assert(part->boids != NULL);
	}
	else if (part->phystype == PART_PHYS_FLUID) {
		BLI_assert(part->fluid != NULL);
	}

	psys_check_boid_data(sim->psys);
}
static int hair_needs_recalc(ParticleSystem *psys)
{
	if (!(psys->flag & PSYS_EDITED) && (!psys->edit || !psys->edit->edited) &&
	    ((psys->flag & PSYS_HAIR_DONE)==0 || psys->recalc & ID_RECALC_PSYS_RESET || (psys->part->flag & PART_HAIR_REGROW && !psys->edit)))
	{
		return 1;
	}

	return 0;
}

static ParticleSettings *particle_settings_localize(ParticleSettings *particle_settings)
{
	ParticleSettings *particle_settings_local;
	BKE_id_copy_ex(NULL,
	               (ID *)&particle_settings->id,
	               (ID **)&particle_settings_local,
	               LIB_ID_COPY_LOCALIZE);
	return particle_settings_local;
}

static void particle_settings_free_local(ParticleSettings *particle_settings)
{
	BKE_libblock_free_datablock(&particle_settings->id, 0);
	BKE_libblock_free_data(&particle_settings->id, false);
	MEM_freeN(particle_settings);
}

/* main particle update call, checks that things are ok on the large scale and
 * then advances in to actual particle calculations depending on particle type */
void particle_system_update(struct Depsgraph *depsgraph, Scene *scene, Object *ob, ParticleSystem *psys, const bool use_render_params)
{
	ParticleSimulationData sim= {0};
	ParticleSettings *part = psys->part;
	ParticleSystem *psys_orig = psys_orig_get(psys);
	float cfra;
	ParticleSystemModifierData *psmd = psys_get_modifier(ob, psys);

	/* drawdata is outdated after ANY change */
	if (psys->pdd) psys->pdd->flag &= ~PARTICLE_DRAW_DATA_UPDATED;

	if (!psys_check_enabled(ob, psys, use_render_params))
		return;

	cfra = DEG_get_ctime(depsgraph);

	sim.depsgraph = depsgraph;
	sim.scene = scene;
	sim.ob = ob;
	sim.psys = psys;
	sim.psmd = psmd;

	/* system was already updated from modifier stack */
	if (sim.psmd->flag & eParticleSystemFlag_psys_updated) {
		sim.psmd->flag &= ~eParticleSystemFlag_psys_updated;
		/* make sure it really was updated to cfra */
		if (psys->cfra == cfra)
			return;
	}

	if (!sim.psmd->mesh_final)
		return;

	if (part->from != PART_FROM_VERT) {
		BKE_mesh_tessface_ensure(sim.psmd->mesh_final);
	}

	/* to verify if we need to restore object afterwards */
	psys->flag &= ~PSYS_OB_ANIM_RESTORE;

	if (psys->recalc & ID_RECALC_PSYS_RESET)
		psys->totunexist = 0;

	/* setup necessary physics type dependent additional data if it doesn't yet exist */
	psys_prepare_physics(&sim);

	switch (part->type) {
		case PART_HAIR:
		{
			/* nothing to do so bail out early */
			if (psys->totpart == 0 && part->totpart == 0) {
				psys_free_path_cache(psys, NULL);
				free_hair(ob, psys, 0);
				psys->flag |= PSYS_HAIR_DONE;
			}
			/* (re-)create hair */
			else if (hair_needs_recalc(psys)) {
				float hcfra=0.0f;
				int i, recalc = psys->recalc;

				free_hair(ob, psys, 0);

				if (psys_orig->edit && psys_orig->free_edit) {
					psys_orig->free_edit(psys_orig->edit);
					psys_orig->edit = NULL;
					psys_orig->free_edit = NULL;
				}

				/* first step is negative so particles get killed and reset */
				psys->cfra= 1.0f;

				ParticleSettings *part_local = part;
				if ((part->flag & PART_HAIR_REGROW) == 0) {
					part_local = particle_settings_localize(part);
					psys->part = part_local;
				}

				for (i=0; i<=part->hair_step; i++) {
					hcfra=100.0f*(float)i/(float)psys->part->hair_step;
					if ((part->flag & PART_HAIR_REGROW)==0)
						BKE_animsys_evaluate_animdata(depsgraph, scene, &part_local->id, part_local->adt, hcfra, ADT_RECALC_ANIM);
					system_step(&sim, hcfra, use_render_params);
					psys->cfra = hcfra;
					psys->recalc = 0;
					save_hair(&sim, hcfra);
				}

				if (part_local != part) {
					particle_settings_free_local(part_local);
					psys->part = part;
				}

				psys->flag |= PSYS_HAIR_DONE;
				psys->recalc = recalc;
			}
			else if (psys->flag & PSYS_EDITED)
				psys->flag |= PSYS_HAIR_DONE;

			if (psys->flag & PSYS_HAIR_DONE)
				hair_step(&sim, cfra, use_render_params);
			break;
		}
		case PART_FLUID:
		{
			particles_fluid_step(&sim, (int)cfra, use_render_params);
			break;
		}
		default:
		{
			switch (part->phystype) {
				case PART_PHYS_NO:
				case PART_PHYS_KEYED:
				{
					PARTICLE_P;
					float disp = psys_get_current_display_percentage(psys, use_render_params);
					bool free_unexisting = false;

					/* Particles without dynamics haven't been reset yet because they don't use pointcache */
					if (psys->recalc & ID_RECALC_PSYS_RESET)
						psys_reset(psys, PSYS_RESET_ALL);

					if (emit_particles(&sim, NULL, cfra) || (psys->recalc & ID_RECALC_PSYS_RESET)) {
						free_keyed_keys(psys);
						distribute_particles(&sim, part->from);
						initialize_all_particles(&sim);
						free_unexisting = true;

						/* flag for possible explode modifiers after this system */
						sim.psmd->flag |= eParticleSystemFlag_Pars;
					}

					LOOP_EXISTING_PARTICLES {
						pa->size = part->size;
						if (part->randsize > 0.0f)
							pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);

						reset_particle(&sim, pa, 0.0, cfra);

						if (psys_frand(psys, p) > disp)
							pa->flag |= PARS_NO_DISP;
						else
							pa->flag &= ~PARS_NO_DISP;
					}

					/* free unexisting after resetting particles */
					if (free_unexisting)
						free_unexisting_particles(&sim);

					if (part->phystype == PART_PHYS_KEYED) {
						psys_count_keyed_targets(&sim);
						set_keyed_keys(&sim);
						psys_update_path_cache(&sim, (int)cfra, use_render_params);
					}
					break;
				}
				default:
				{
					/* the main dynamic particle system step */
					system_step(&sim, cfra, use_render_params);
					break;
				}
			}
			break;
		}
	}

	/* make sure emitter is left at correct time (particle emission can change this) */
	if (psys->flag & PSYS_OB_ANIM_RESTORE) {
		evaluate_emitter_anim(depsgraph, scene, ob, cfra);
		psys->flag &= ~PSYS_OB_ANIM_RESTORE;
	}

	if (psys_orig->edit) {
		psys_orig->edit->flags |= PT_CACHE_EDIT_UPDATE_PARTICLE_FROM_EVAL;
	}

	psys->cfra = cfra;
	psys->recalc = 0;

	if (DEG_is_active(depsgraph)) {
		if (psys_orig != psys) {
			if (psys_orig->edit != NULL &&
			    psys_orig->edit->psys == psys_orig)
			{
				psys_orig->edit->psys_eval = psys;
				psys_orig->edit->psmd_eval = psmd;
			}
			psys_orig->flag = (psys->flag & ~PSYS_SHARED_CACHES);
			psys_orig->cfra = psys->cfra;
			psys_orig->recalc = psys->recalc;
		}
	}

	/* save matrix for duplicators, at rendertime the actual dupliobject's matrix is used so don't update! */
	invert_m4_m4(psys->imat, ob->obmat);

	BKE_particle_batch_cache_dirty_tag(psys, BKE_PARTICLE_BATCH_DIRTY_ALL);
}

/* ID looper */

void BKE_particlesystem_id_loop(ParticleSystem *psys, ParticleSystemIDFunc func, void *userdata)
{
	ParticleTarget *pt;

	func(psys, (ID **)&psys->part, userdata, IDWALK_CB_USER | IDWALK_CB_NEVER_NULL);
	func(psys, (ID **)&psys->target_ob, userdata, IDWALK_CB_NOP);
	func(psys, (ID **)&psys->parent, userdata, IDWALK_CB_NOP);

	for (pt = psys->targets.first; pt; pt = pt->next) {
		func(psys, (ID **)&pt->ob, userdata, IDWALK_CB_NOP);
	}

	/* Even though psys->part should never be NULL, this can happen as an exception during deletion.
	 * See ID_REMAP_SKIP/FORCE/FLAG_NEVER_NULL_USAGE in BKE_library_remap. */
	if (psys->part && psys->part->phystype == PART_PHYS_BOIDS) {
		ParticleData *pa;
		int p;

		for (p = 0, pa = psys->particles; p < psys->totpart; p++, pa++) {
			func(psys, (ID **)&pa->boid->ground, userdata, IDWALK_CB_NOP);
		}
	}
}

void BKE_particlesystem_reset_all(struct Object *object)
{
	for (ModifierData *md = object->modifiers.first; md != NULL; md = md->next) {
		if (md->type != eModifierType_ParticleSystem) {
			continue;
		}
		ParticleSystemModifierData *psmd = (ParticleSystemModifierData *)md;
		ParticleSystem *psys = psmd->psys;
		psys->recalc |= ID_RECALC_PSYS_RESET;
	}
}

/* **** Depsgraph evaluation **** */

void BKE_particle_settings_eval_reset(
        struct Depsgraph *depsgraph,
        ParticleSettings *particle_settings)
{
	DEG_debug_print_eval(depsgraph,
	                     __func__,
	                     particle_settings->id.name,
	                     particle_settings);
	particle_settings->id.recalc |= ID_RECALC_PSYS_RESET;
}

void BKE_particle_system_eval_init(struct Depsgraph *depsgraph,
                                   Object *object)
{
	DEG_debug_print_eval(depsgraph, __func__, object->id.name, object);
	for (ParticleSystem *psys = object->particlesystem.first;
	     psys != NULL;
	     psys = psys->next)
	{
		psys->recalc |= (psys->part->id.recalc & ID_RECALC_PSYS_ALL);
	}
}