/* particle_system.c
 *
 *
 * $Id$
 *
 * ***** BEGIN GPL LICENSE BLOCK *****
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * 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.
 *
 * The Original Code is: all of this file.
 *
 * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn.
 *
 * ***** END GPL LICENSE BLOCK *****
 */

#include <stddef.h>
#include "BLI_storage.h" /* _LARGEFILE_SOURCE */

#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.h"
#include "DNA_object_types.h"
#include "DNA_material_types.h"
#include "DNA_curve_types.h"
#include "DNA_group_types.h"
#include "DNA_scene_types.h"
#include "DNA_texture_types.h"
#include "DNA_ipo_types.h" // XXX old animation system stuff... to be removed!
#include "DNA_listBase.h"

#include "BLI_edgehash.h"
#include "BLI_rand.h"
#include "BLI_jitter.h"
#include "BLI_math.h"
#include "BLI_blenlib.h"
#include "BLI_kdtree.h"
#include "BLI_kdopbvh.h"
#include "BLI_listbase.h"
#include "BLI_threads.h"
#include "BLI_storage.h" /* For _LARGEFILE64_SOURCE;  zlib needs this on some systems */
#include "BLI_utildefines.h"

#include "BKE_main.h"
#include "BKE_animsys.h"
#include "BKE_boids.h"
#include "BKE_cdderivedmesh.h"
#include "BKE_collision.h"
#include "BKE_displist.h"
#include "BKE_effect.h"
#include "BKE_particle.h"
#include "BKE_global.h"

#include "BKE_DerivedMesh.h"
#include "BKE_object.h"
#include "BKE_material.h"
#include "BKE_cloth.h"
#include "BKE_depsgraph.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 "PIL_time.h"

#include "RE_shader_ext.h"

/* fluid sim particle import */
#ifndef DISABLE_ELBEEM
#include "DNA_object_fluidsim.h"
#include "LBM_fluidsim.h"
#include <zlib.h>
#include <string.h>

#ifdef WIN32
#ifndef snprintf
#define snprintf _snprintf
#endif
#endif

#endif // DISABLE_ELBEEM

/************************************************/
/*			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 ELEM3(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
}

static int psys_get_current_display_percentage(ParticleSystem *psys)
{
	ParticleSettings *part=psys->part;

	if((psys->renderdata && !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 100;

	return psys->part->disp;
}

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

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

static int get_psys_child_number(struct Scene *scene, ParticleSystem *psys)
{
	int nbr;

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

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

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

static int get_psys_tot_child(struct Scene *scene, ParticleSystem *psys)
{
	return psys->totpart*get_psys_child_number(scene, psys);
}

static void alloc_child_particles(ParticleSystem *psys, int tot)
{
	if(psys->child){
		/* only re-allocate if we have to */
		if(psys->part->childtype && psys->totchild == tot) {
			memset(psys->child, 0, tot*sizeof(ChildParticle));
			return;
		}

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

	if(psys->part->childtype) {
		psys->totchild= tot;
		if(psys->totchild)
			psys->child= MEM_callocN(psys->totchild*sizeof(ChildParticle), "child_particles");
	}
}

void psys_calc_dmcache(Object *ob, DerivedMesh *dm, 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;
	PARTICLE_P;
	
	/* CACHE LOCATIONS */
	if(!dm->deformedOnly) {
		/* Will use later to speed up subsurf/derivedmesh */
		LinkNode *node, *nodedmelem, **nodearray;
		int totdmelem, totelem, i, *origindex;

		if(psys->part->from == PART_FROM_VERT) {
			totdmelem= dm->getNumVerts(dm);
			totelem= me->totvert;
			origindex= dm->getVertDataArray(dm, CD_ORIGINDEX);
		}
		else { /* FROM_FACE/FROM_VOLUME */
			totdmelem= dm->getNumFaces(dm);
			totelem= me->totface;
			origindex= dm->getFaceDataArray(dm, CD_ORIGINDEX);
		}
	
		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++, origindex++, node++) {
			node->link= SET_INT_IN_POINTER(i);

			if(*origindex != -1) {
				if(nodearray[*origindex]) {
					/* prepend */
					node->next = nodearray[*origindex];
					nodearray[*origindex]= node;
				}
				else
					nodearray[*origindex]= node;
			}
		}
		
		/* cache the verts/faces! */
		LOOP_PARTICLES {
			if(pa->num < 0) {
				pa->num_dmcache = -1;
				continue;
			}

			if(psys->part->from == PART_FROM_VERT) {
				if(nodearray[pa->num])
					pa->num_dmcache= GET_INT_FROM_POINTER(nodearray[pa->num]->link);
			}
			else { /* FROM_FACE/FROM_VOLUME */
				/* Note that sometimes the pa->num is over the nodearray size, this is bad, maybe there is a better place to fix this,
				 * but for now passing NULL is OK. every face will be searched for the particle so its slower - Campbell */
				pa->num_dmcache= psys_particle_dm_face_lookup(ob, dm, pa->num, pa->fuv, pa->num < totelem ? nodearray[pa->num] : NULL);
			}
		}

		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 incase */
		
		LOOP_PARTICLES
			pa->num_dmcache = -1;
	}
}

static void distribute_particles_in_grid(DerivedMesh *dm, ParticleSystem *psys)
{
	ParticleData *pa=NULL;
	float min[3], max[3], delta[3], d;
	MVert *mv, *mvert = dm->getVertDataArray(dm,0);
	int totvert=dm->getNumVerts(dm), from=psys->part->from;
	int i, j, k, p, res=psys->part->grid_res, size[3], axis;

	mv=mvert;

	/* find bounding box of dm */
	copy_v3_v3(min, mv->co);
	copy_v3_v3(max, mv->co);
	mv++;

	for(i=1; i<totvert; i++, mv++){
		min[0]=MIN2(min[0],mv->co[0]);
		min[1]=MIN2(min[1],mv->co[1]);
		min[2]=MIN2(min[2],mv->co[2]);

		max[0]=MAX2(max[0],mv->co[0]);
		max[1]=MAX2(max[1],mv->co[1]);
		max[2]=MAX2(max[2],mv->co[2]);
	}

	VECSUB(delta,max,min);

	/* determine major axis */
	axis = (delta[0]>=delta[1]) ? 0 : ((delta[1]>=delta[2]) ? 1 : 2);
	 
	d = delta[axis]/(float)res;

	size[axis] = res;
	size[(axis+1)%3] = (int)ceil(delta[(axis+1)%3]/d);
	size[(axis+2)%3] = (int)ceil(delta[(axis+2)%3]/d);

	/* float errors grrr.. */
	size[(axis+1)%3] = MIN2(size[(axis+1)%3],res);
	size[(axis+2)%3] = MIN2(size[(axis+2)%3],res);

	size[0] = MAX2(size[0], 1);
	size[1] = MAX2(size[1], 1);
	size[2] = MAX2(size[2], 1);

	/* no full offset for flat/thin objects */
	min[0]+= d < delta[0] ? d/2.f : delta[0]/2.f;
	min[1]+= d < delta[1] ? d/2.f : delta[1]/2.f;
	min[2]+= d < delta[2] ? d/2.f : delta[2]/2.f;

	for(i=0,p=0,pa=psys->particles; i<res; i++){
		for(j=0; j<res; j++){
			for(k=0; k<res; k++,p++,pa++){
				pa->fuv[0] = min[0] + (float)i*d;
				pa->fuv[1] = min[1] + (float)j*d;
				pa->fuv[2] = min[2] + (float)k*d;
				pa->flag |= PARS_UNEXIST;
				pa->hair_index = 0; /* abused in volume calculation */
			}
		}
	}

	/* enable particles near verts/edges/faces/inside surface */
	if(from==PART_FROM_VERT){
		float vec[3];

		pa=psys->particles;

		min[0] -= d/2.0f;
		min[1] -= d/2.0f;
		min[2] -= d/2.0f;

		for(i=0,mv=mvert; i<totvert; i++,mv++){
			sub_v3_v3v3(vec,mv->co,min);
			vec[0]/=delta[0];
			vec[1]/=delta[1];
			vec[2]/=delta[2];
			(pa	+((int)(vec[0]*(size[0]-1))*res
				+(int)(vec[1]*(size[1]-1)))*res
				+(int)(vec[2]*(size[2]-1)))->flag &= ~PARS_UNEXIST;
		}
	}
	else if(ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)){
		float co1[3], co2[3];

		MFace *mface= NULL, *mface_array;
		float v1[3], v2[3], v3[3], v4[4], lambda;
		int a, a1, a2, a0mul, a1mul, a2mul, totface;
		int amax= from==PART_FROM_FACE ? 3 : 1;

		totface=dm->getNumFaces(dm);
		mface_array= dm->getFaceDataArray(dm,CD_MFACE);
		
		for(a=0; a<amax; a++){
			if(a==0){ a0mul=res*res; a1mul=res; a2mul=1; }
			else if(a==1){ a0mul=res; a1mul=1; a2mul=res*res; }
			else{ a0mul=1; a1mul=res*res; a2mul=res; }

			for(a1=0; a1<size[(a+1)%3]; a1++){
				for(a2=0; a2<size[(a+2)%3]; a2++){
					mface= mface_array;

					pa = psys->particles + a1*a1mul + a2*a2mul;
					copy_v3_v3(co1, pa->fuv);
					co1[a] -= d < delta[a] ? d/2.f : delta[a]/2.f;
					copy_v3_v3(co2, co1);
					co2[a] += delta[a] + 0.001f*d;
					co1[a] -= 0.001f*d;
					
					/* lets intersect the faces */
					for(i=0; i<totface; i++,mface++){
						copy_v3_v3(v1, mvert[mface->v1].co);
						copy_v3_v3(v2, mvert[mface->v2].co);
						copy_v3_v3(v3, mvert[mface->v3].co);

						if(isect_axial_line_tri_v3(a, co1, co2, v2, v3, v1, &lambda)){
							if(from==PART_FROM_FACE)
								(pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
							else /* store number of intersections */
								(pa+(int)(lambda*size[a])*a0mul)->hair_index++;
						}
						
						if(mface->v4){
							copy_v3_v3(v4, mvert[mface->v4].co);

							if(isect_axial_line_tri_v3(a, co1, co2, v4, v1, v3, &lambda)){
								if(from==PART_FROM_FACE)
									(pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
								else
									(pa+(int)(lambda*size[a])*a0mul)->hair_index++;
							}
						}
					}

					if(from==PART_FROM_VOLUME){
						int in=pa->hair_index%2;
						if(in) pa->hair_index++;
						for(i=0; i<size[0]; i++){
							if(in || (pa+i*a0mul)->hair_index%2)
								(pa+i*a0mul)->flag &= ~PARS_UNEXIST;
							/* odd intersections == in->out / out->in */
							/* even intersections -> in stays same */
							in=(in + (pa+i*a0mul)->hair_index) % 2;
						}
					}
				}
			}
		}
	}

	if(psys->part->flag & PART_GRID_HEXAGONAL) {
		for(i=0,p=0,pa=psys->particles; i<res; i++){
			for(j=0; j<res; j++){
				for(k=0; k<res; k++,p++,pa++){
					if(j%2)
						pa->fuv[0] += d/2.f;

					if(k%2) {
						pa->fuv[0] += d/2.f;
						pa->fuv[1] += d/2.f;
					}
				}
			}
		}
	}

	if(psys->part->flag & PART_GRID_INVERT){
		for(i=0; i<size[0]; i++){
			for(j=0; j<size[1]; j++){
				pa=psys->particles + res*(i*res + j);
				for(k=0; k<size[2]; k++, pa++){
					pa->flag ^= PARS_UNEXIST;
				}
			}
		}
	}

	if(psys->part->grid_rand > 0.f) {
		float rfac = d * psys->part->grid_rand;
		for(p=0,pa=psys->particles; p<psys->totpart; p++,pa++){
			if(pa->flag & PARS_UNEXIST)
				continue;

			pa->fuv[0] += rfac * (PSYS_FRAND(p + 31) - 0.5f);
			pa->fuv[1] += rfac * (PSYS_FRAND(p + 32) - 0.5f);
			pa->fuv[2] += rfac * (PSYS_FRAND(p + 33) - 0.5f);
		}
	}
}

/* modified copy from rayshade.c */
static void hammersley_create(float *out, int n, int seed, float amount)
{
	RNG *rng;
	double p, t, offs[2];
	int k, kk;

	rng = rng_new(31415926 + n + seed);
	offs[0]= rng_getDouble(rng) + amount;
	offs[1]= rng_getDouble(rng) + amount;
	rng_free(rng);

	for (k = 0; k < n; k++) {
		t = 0;
		for (p = 0.5, kk = k; kk; p *= 0.5, kk >>= 1)
			if (kk & 1) /* kk mod 2 = 1 */
				t += p;
    
		out[2*k + 0]= fmod((double)k/(double)n + offs[0], 1.0);
		out[2*k + 1]= fmod(t + offs[1], 1.0);
	}
}

/* modified copy from effect.c */
static void init_mv_jit(float *jit, int num, int seed2, float amount)
{
	RNG *rng;
	float *jit2, x, rad1, rad2, rad3;
	int i, num2;

	if(num==0) return;

	rad1= (float)(1.0/sqrt((float)num));
	rad2= (float)(1.0/((float)num));
	rad3= (float)sqrt((float)num)/((float)num);

	rng = rng_new(31415926 + num + seed2);
	x= 0;
		num2 = 2 * num;
	for(i=0; i<num2; i+=2) {
	
		jit[i]= x + amount*rad1*(0.5f - rng_getFloat(rng));
		jit[i+1]= i/(2.0f*num) + amount*rad1*(0.5f - rng_getFloat(rng));
		
		jit[i]-= (float)floor(jit[i]);
		jit[i+1]-= (float)floor(jit[i+1]);
		
		x+= rad3;
		x -= (float)floor(x);
	}

	jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");

	for (i=0 ; i<4 ; i++) {
		BLI_jitterate1(jit, jit2, num, rad1);
		BLI_jitterate1(jit, jit2, num, rad1);
		BLI_jitterate2(jit, jit2, num, rad2);
	}
	MEM_freeN(jit2);
	rng_free(rng);
}

static void psys_uv_to_w(float u, float v, int quad, float *w)
{
	float vert[4][3], co[3];

	if(!quad) {
		if(u+v > 1.0f)
			v= 1.0f-v;
		else
			u= 1.0f-u;
	}

	vert[0][0]= 0.0f; vert[0][1]= 0.0f; vert[0][2]= 0.0f;
	vert[1][0]= 1.0f; vert[1][1]= 0.0f; vert[1][2]= 0.0f;
	vert[2][0]= 1.0f; vert[2][1]= 1.0f; vert[2][2]= 0.0f;

	co[0]= u;
	co[1]= v;
	co[2]= 0.0f;

	if(quad) {
		vert[3][0]= 0.0f; vert[3][1]= 1.0f; vert[3][2]= 0.0f;
		interp_weights_poly_v3( w,vert, 4, co);
	}
	else {
		interp_weights_poly_v3( w,vert, 3, co);
		w[3]= 0.0f;
	}
}

/* Find the index in "sum" array before "value" is crossed. */
static int binary_search_distribution(float *sum, int n, float value)
{
	int mid, low=0, high=n;

	if(value == 0.f)
		return 0;

	while(low <= high) {
		mid= (low + high)/2;
		
		if(sum[mid] < value && value <= sum[mid+1])
			return mid;
		
		if(sum[mid] >= value)
			high= mid - 1;
		else if(sum[mid] < value)
			low= mid + 1;
		else
			return mid;
	}

	return low;
}

/* the max number if calls to rng_* funcs within psys_thread_distribute_particle
 * be sure to keep up to date if this changes */
#define PSYS_RND_DIST_SKIP 2

/* note: this function must be thread safe, for from == PART_FROM_CHILD */
#define ONLY_WORKING_WITH_PA_VERTS 0
static void psys_thread_distribute_particle(ParticleThread *thread, ParticleData *pa, ChildParticle *cpa, int p)
{
	ParticleThreadContext *ctx= thread->ctx;
	Object *ob= ctx->sim.ob;
	DerivedMesh *dm= ctx->dm;
	ParticleData *tpa;
/*	ParticleSettings *part= ctx->sim.psys->part; */
	float *v1, *v2, *v3, *v4, nor[3], orco1[3], co1[3], co2[3], nor1[3];
	float cur_d, min_d, randu, randv;
	int from= ctx->from;
	int cfrom= ctx->cfrom;
	int distr= ctx->distr;
	int i, intersect, tot;
	int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */

	if(from == PART_FROM_VERT) {
		/* TODO_PARTICLE - use original index */
		pa->num= ctx->index[p];
		pa->fuv[0] = 1.0f;
		pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;

#if ONLY_WORKING_WITH_PA_VERTS
		if(ctx->tree){
			KDTreeNearest ptn[3];
			int w, maxw;

			psys_particle_on_dm(ctx->dm,from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co1,0,0,0,orco1,0);
			transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
			maxw = BLI_kdtree_find_n_nearest(ctx->tree,3,orco1,NULL,ptn);

			for(w=0; w<maxw; w++){
				pa->verts[w]=ptn->num;
			}
		}
#endif
	}
	else if(from == PART_FROM_FACE || from == PART_FROM_VOLUME) {
		MFace *mface;

		pa->num = i = ctx->index[p];
		mface = dm->getFaceData(dm,i,CD_MFACE);
		
		switch(distr){
		case PART_DISTR_JIT:
			ctx->jitoff[i] = fmod(ctx->jitoff[i],(float)ctx->jitlevel);
			psys_uv_to_w(ctx->jit[2*(int)ctx->jitoff[i]], ctx->jit[2*(int)ctx->jitoff[i]+1], mface->v4, pa->fuv);
			ctx->jitoff[i]++;
			break;
		case PART_DISTR_RAND:
			randu= rng_getFloat(thread->rng);
			randv= rng_getFloat(thread->rng);
			rng_skip_tot -= 2;

			psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
			break;
		}
		pa->foffset= 0.0f;
		
		/* experimental */
		if(from==PART_FROM_VOLUME){
			MVert *mvert=dm->getVertDataArray(dm,CD_MVERT);

			tot=dm->getNumFaces(dm);

			psys_interpolate_face(mvert,mface,0,0,pa->fuv,co1,nor,0,0,0,0);

			normalize_v3(nor);
			mul_v3_fl(nor,-100.0);

			VECADD(co2,co1,nor);

			min_d=2.0;
			intersect=0;

			for(i=0,mface=dm->getFaceDataArray(dm,CD_MFACE); i<tot; i++,mface++){
				if(i==pa->num) continue;

				v1=mvert[mface->v1].co;
				v2=mvert[mface->v2].co;
				v3=mvert[mface->v3].co;

				if(isect_line_tri_v3(co1, co2, v2, v3, v1, &cur_d, 0)){
					if(cur_d<min_d){
						min_d=cur_d;
						pa->foffset=cur_d*50.0f; /* to the middle of volume */
						intersect=1;
					}
				}
				if(mface->v4){
					v4=mvert[mface->v4].co;

					if(isect_line_tri_v3(co1, co2, v4, v1, v3, &cur_d, 0)){
						if(cur_d<min_d){
							min_d=cur_d;
							pa->foffset=cur_d*50.0f; /* to the middle of volume */
							intersect=1;
						}
					}
				}
			}
			if(intersect==0)
				pa->foffset=0.0;
			else switch(distr){
				case PART_DISTR_JIT:
					pa->foffset*= ctx->jit[p%(2*ctx->jitlevel)];
					break;
				case PART_DISTR_RAND:
					pa->foffset*=BLI_frand();
					break;
			}
		}
	}
	else if(from == PART_FROM_PARTICLE) {
		tpa=ctx->tpars+ctx->index[p];
		pa->num=ctx->index[p];
		pa->fuv[0]=tpa->fuv[0];
		pa->fuv[1]=tpa->fuv[1];
		/* abusing foffset a little for timing in near reaction */
		pa->foffset=ctx->weight[ctx->index[p]];
		ctx->weight[ctx->index[p]]+=ctx->maxweight;
	}
	else if(from == PART_FROM_CHILD) {
		MFace *mf;

		if(ctx->index[p] < 0) {
			cpa->num=0;
			cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]=0.0f;
			cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
			return;
		}

		mf= dm->getFaceData(dm, ctx->index[p], CD_MFACE);

		randu= rng_getFloat(thread->rng);
		randv= rng_getFloat(thread->rng);
		rng_skip_tot -= 2;

		psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);

		cpa->num = ctx->index[p];

		if(ctx->tree){
			KDTreeNearest ptn[10];
			int w,maxw;//, do_seams;
			float maxd,mind,/*dd,*/totw=0.0;
			int parent[10];
			float pweight[10];

			psys_particle_on_dm(dm,cfrom,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co1,nor1,NULL,NULL,orco1,NULL);
			transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
			maxw = BLI_kdtree_find_n_nearest(ctx->tree,4,orco1,NULL,ptn);

			maxd=ptn[maxw-1].dist;
			mind=ptn[0].dist;
			/*dd=maxd-mind;*/ /*UNUSED*/
			
			/* the weights here could be done better */
			for(w=0; w<maxw; w++){
				parent[w]=ptn[w].index;
				pweight[w]=(float)pow(2.0,(double)(-6.0f*ptn[w].dist/maxd));
			}
			for(;w<10; w++){
				parent[w]=-1;
				pweight[w]=0.0f;
			}

			for(w=0,i=0; w<maxw && i<4; w++){
				if(parent[w]>=0){
					cpa->pa[i]=parent[w];
					cpa->w[i]=pweight[w];
					totw+=pweight[w];
					i++;
				}
			}
			for(;i<4; i++){
				cpa->pa[i]=-1;
				cpa->w[i]=0.0f;
			}

			if(totw>0.0f) for(w=0; w<4; w++)
				cpa->w[w]/=totw;

			cpa->parent=cpa->pa[0];
		}
	}

	if(rng_skip_tot > 0) /* should never be below zero */
		rng_skip(thread->rng, rng_skip_tot);
}

static void *exec_distribution(void *data)
{
	ParticleThread *thread= (ParticleThread*)data;
	ParticleSystem *psys= thread->ctx->sim.psys;
	ParticleData *pa;
	ChildParticle *cpa;
	int p, totpart;

	if(thread->ctx->from == PART_FROM_CHILD) {
		totpart= psys->totchild;
		cpa= psys->child;

		for(p=0; p<totpart; p++, cpa++) {
			if(thread->ctx->skip) /* simplification skip */
				rng_skip(thread->rng, PSYS_RND_DIST_SKIP * thread->ctx->skip[p]);

			if((p+thread->num) % thread->tot == 0)
				psys_thread_distribute_particle(thread, NULL, cpa, p);
			else /* thread skip */
				rng_skip(thread->rng, PSYS_RND_DIST_SKIP);
		}
	}
	else {
		totpart= psys->totpart;
		pa= psys->particles + thread->num;
		for(p=thread->num; p<totpart; p+=thread->tot, pa+=thread->tot)
			psys_thread_distribute_particle(thread, pa, NULL, p);
	}

	return 0;
}

/* not thread safe, but qsort doesn't take userdata argument */
static int *COMPARE_ORIG_INDEX = NULL;
static int compare_orig_index(const void *p1, const void *p2)
{
	int index1 = COMPARE_ORIG_INDEX[*(const int*)p1];
	int index2 = COMPARE_ORIG_INDEX[*(const int*)p2];

	if(index1 < index2)
		return -1;
	else if(index1 == index2) {
		/* this pointer comparison appears to make qsort stable for glibc,
		 * and apparently on solaris too, makes the renders reproducable */
		if(p1 < p2)
			return -1;
		else if(p1 == p2)
			return 0;
		else
			return 1;
	}
	else
		return 1;
}

/* creates a distribution of coordinates on a DerivedMesh	*/
/*															*/
/* 1. lets check from what we are emitting					*/
/* 2. now we know that we have something to emit from so	*/
/*	  let's calculate some weights							*/
/* 2.1 from even distribution								*/
/* 2.2 and from vertex groups								*/
/* 3. next we determine the indexes of emitting thing that	*/
/*	  the particles will have								*/
/* 4. let's do jitter if we need it							*/
/* 5. now we're ready to set the indexes & distributions to	*/
/*	  the particles											*/
/* 6. and we're done!										*/

/* This is to denote functionality that does not yet work with mesh - only derived mesh */
static int psys_threads_init_distribution(ParticleThread *threads, Scene *scene, DerivedMesh *finaldm, int from)
{
	ParticleThreadContext *ctx= threads[0].ctx;
	Object *ob= ctx->sim.ob;
	ParticleSystem *psys= ctx->sim.psys;
	Object *tob;
	ParticleData *pa=0, *tpars= 0;
	ParticleSettings *part;
	ParticleSystem *tpsys;
	ParticleSeam *seams= 0;
	ChildParticle *cpa=0;
	KDTree *tree=0;
	DerivedMesh *dm= NULL;
	float *jit= NULL;
	int i, seed, p=0, totthread= threads[0].tot;
	int /*no_distr=0,*/ cfrom=0;
	int tot=0, totpart, *index=0, children=0, totseam=0;
	//int *vertpart=0;
	int jitlevel= 1, distr;
	float *weight=0,*sum=0,*jitoff=0;
	float cur, maxweight=0.0, tweight, totweight, co[3], nor[3], orco[3], ornor[3];
	
	if(ob==0 || psys==0 || psys->part==0)
		return 0;

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

	if (!finaldm->deformedOnly && !finaldm->getFaceDataArray(finaldm, CD_ORIGINDEX)) {
		printf("Can't create particles with the current modifier stack, disable destructive modifiers\n");
// XXX		error("Can't paint with the current modifier stack, disable destructive modifiers");
		return 0;
	}

	BLI_srandom(31415926 + psys->seed);
	
	if(from==PART_FROM_CHILD){
		distr=PART_DISTR_RAND;
		BLI_srandom(31415926 + psys->seed + psys->child_seed);

		if(part->from!=PART_FROM_PARTICLE && part->childtype==PART_CHILD_FACES){
			dm= finaldm;
			children=1;

			tree=BLI_kdtree_new(totpart);

			for(p=0,pa=psys->particles; p<totpart; p++,pa++){
				psys_particle_on_dm(dm,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co,nor,0,0,orco,ornor);
				transform_mesh_orco_verts((Mesh*)ob->data, &orco, 1, 1);
				BLI_kdtree_insert(tree, p, orco, ornor);
			}

			BLI_kdtree_balance(tree);

			totpart=get_psys_tot_child(scene, psys);
			cfrom=from=PART_FROM_FACE;
		}
		else{
			/* no need to figure out distribution */
			int child_nbr= get_psys_child_number(scene, psys);

			totpart= get_psys_tot_child(scene, psys);
			alloc_child_particles(psys, totpart);
			cpa=psys->child;
			for(i=0; i<child_nbr; i++){
				for(p=0; p<psys->totpart; p++,cpa++){
					float length=2.0;
					cpa->parent=p;
					
					/* create even spherical distribution inside unit sphere */
					while(length>=1.0f){
						cpa->fuv[0]=2.0f*BLI_frand()-1.0f;
						cpa->fuv[1]=2.0f*BLI_frand()-1.0f;
						cpa->fuv[2]=2.0f*BLI_frand()-1.0f;
						length=len_v3(cpa->fuv);
					}

					cpa->num=-1;
				}
			}
			/* dmcache must be updated for parent particles if children from faces is used */
			psys_calc_dmcache(ob, finaldm, psys);

			return 0;
		}
	}
	else{
		dm= CDDM_from_mesh((Mesh*)ob->data, ob);

		/* special handling of grid distribution */
		if(part->distr==PART_DISTR_GRID && from != PART_FROM_VERT){
			distribute_particles_in_grid(dm,psys);
			dm->release(dm);
			return 0;
		}

		/* we need orco for consistent distributions */
		DM_add_vert_layer(dm, CD_ORCO, CD_ASSIGN, get_mesh_orco_verts(ob));

		distr=part->distr;

		if(from==PART_FROM_VERT){
			MVert *mv= dm->getVertDataArray(dm, CD_MVERT);
			float (*orcodata)[3]= dm->getVertDataArray(dm, CD_ORCO);
			int totvert = dm->getNumVerts(dm);

			tree=BLI_kdtree_new(totvert);

			for(p=0; p<totvert; p++){
				if(orcodata) {
					VECCOPY(co,orcodata[p])
					transform_mesh_orco_verts((Mesh*)ob->data, &co, 1, 1);
				}
				else
					VECCOPY(co,mv[p].co)
				BLI_kdtree_insert(tree,p,co,NULL);
			}

			BLI_kdtree_balance(tree);
		}
	}

	/* 1. */
	switch(from){
		case PART_FROM_VERT:
			tot = dm->getNumVerts(dm);
			break;
		case PART_FROM_VOLUME:
		case PART_FROM_FACE:
			tot = dm->getNumFaces(dm);
			break;
		case PART_FROM_PARTICLE:
			if(psys->target_ob)
				tob=psys->target_ob;
			else
				tob=ob;

			if((tpsys=BLI_findlink(&tob->particlesystem,psys->target_psys-1))){
				tpars=tpsys->particles;
				tot=tpsys->totpart;
			}
			break;
	}

	if(tot==0){
		/*no_distr=1;*/ /*UNUSED*/
		if(children){
			if(G.f & G_DEBUG)
				fprintf(stderr,"Particle child distribution error: Nothing to emit from!\n");
			if(psys->child) {
				for(p=0,cpa=psys->child; p<totpart; p++,cpa++){
					cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]= 0.0;
					cpa->foffset= 0.0f;
					cpa->parent=0;
					cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
					cpa->num= -1;
				}
			}
		}
		else {
			if(G.f & G_DEBUG)
				fprintf(stderr,"Particle distribution error: Nothing to emit from!\n");
			for(p=0,pa=psys->particles; p<totpart; p++,pa++){
				pa->fuv[0]=pa->fuv[1]=pa->fuv[2]= pa->fuv[3]= 0.0;
				pa->foffset= 0.0f;
				pa->num= -1;
			}
		}

		if(dm != finaldm) dm->release(dm);
		return 0;
	}

	/* 2. */

	weight=MEM_callocN(sizeof(float)*tot, "particle_distribution_weights");
	index=MEM_callocN(sizeof(int)*totpart, "particle_distribution_indexes");
	sum=MEM_callocN(sizeof(float)*(tot+1), "particle_distribution_sum");
	jitoff=MEM_callocN(sizeof(float)*tot, "particle_distribution_jitoff");

	/* 2.1 */
	if((part->flag&PART_EDISTR || children) && ELEM(from,PART_FROM_PARTICLE,PART_FROM_VERT)==0){
		MVert *v1, *v2, *v3, *v4;
		float totarea=0.0, co1[3], co2[3], co3[3], co4[3];
		float (*orcodata)[3];
		
		orcodata= dm->getVertDataArray(dm, CD_ORCO);

		for(i=0; i<tot; i++){
			MFace *mf=dm->getFaceData(dm,i,CD_MFACE);

			if(orcodata) {
				VECCOPY(co1, orcodata[mf->v1]);
				VECCOPY(co2, orcodata[mf->v2]);
				VECCOPY(co3, orcodata[mf->v3]);
				transform_mesh_orco_verts((Mesh*)ob->data, &co1, 1, 1);
				transform_mesh_orco_verts((Mesh*)ob->data, &co2, 1, 1);
				transform_mesh_orco_verts((Mesh*)ob->data, &co3, 1, 1);
			}
			else {
				v1= (MVert*)dm->getVertData(dm,mf->v1,CD_MVERT);
				v2= (MVert*)dm->getVertData(dm,mf->v2,CD_MVERT);
				v3= (MVert*)dm->getVertData(dm,mf->v3,CD_MVERT);
				VECCOPY(co1, v1->co);
				VECCOPY(co2, v2->co);
				VECCOPY(co3, v3->co);
			}

			if (mf->v4){
				if(orcodata) {
					VECCOPY(co4, orcodata[mf->v4]);
					transform_mesh_orco_verts((Mesh*)ob->data, &co4, 1, 1);
				}
				else {
					v4= (MVert*)dm->getVertData(dm,mf->v4,CD_MVERT);
					VECCOPY(co4, v4->co);
				}
				cur= area_quad_v3(co1, co2, co3, co4);
			}
			else
				cur= area_tri_v3(co1, co2, co3);
			
			if(cur>maxweight)
				maxweight=cur;

			weight[i]= cur;
			totarea+=cur;
		}

		for(i=0; i<tot; i++)
			weight[i] /= totarea;

		maxweight /= totarea;
	}
	else if(from==PART_FROM_PARTICLE){
		float val=(float)tot/(float)totpart;
		for(i=0; i<tot; i++)
			weight[i]=val;
		maxweight=val;
	}
	else{
		float min=1.0f/(float)(MIN2(tot,totpart));
		for(i=0; i<tot; i++)
			weight[i]=min;
		maxweight=min;
	}

	/* 2.2 */
	if(ELEM3(from,PART_FROM_VERT,PART_FROM_FACE,PART_FROM_VOLUME)){
		float *vweight= psys_cache_vgroup(dm,psys,PSYS_VG_DENSITY);

		if(vweight){
			if(from==PART_FROM_VERT) {
				for(i=0;i<tot; i++)
					weight[i]*=vweight[i];
			}
			else { /* PART_FROM_FACE / PART_FROM_VOLUME */
				for(i=0;i<tot; i++){
					MFace *mf=dm->getFaceData(dm,i,CD_MFACE);
					tweight = vweight[mf->v1] + vweight[mf->v2] + vweight[mf->v3];
				
					if(mf->v4) {
						tweight += vweight[mf->v4];
						tweight /= 4.0;
					}
					else {
						tweight /= 3.0;
					}

					weight[i]*=tweight;
				}
			}
			MEM_freeN(vweight);
		}
	}

	/* 3. */
	totweight= 0.0f;
	for(i=0;i<tot; i++)
		totweight += weight[i];

	if(totweight > 0.0f)
		totweight= 1.0f/totweight;

	sum[0]= 0.0f;
	for(i=0;i<tot; i++)
		sum[i+1]= sum[i]+weight[i]*totweight;
	
	if((part->flag&PART_TRAND) || (part->simplify_flag&PART_SIMPLIFY_ENABLE)) {
		float pos;

		for(p=0; p<totpart; p++) {
			/* In theory sys[tot] should be 1.0, but due to float errors this is not necessarily always true, so scale pos accordingly. */
			pos= BLI_frand() * sum[tot];
			index[p]= binary_search_distribution(sum, tot, pos);
			index[p]= MIN2(tot-1, index[p]);
			jitoff[index[p]]= pos;
		}
	}
	else {
		double step, pos;
		
		step= (totpart <= 1)? 0.5: 1.0/(totpart-1);
		pos= 1e-16f; /* tiny offset to avoid zero weight face */
		i= 0;

		for(p=0; p<totpart; p++, pos+=step) {
			while((i < tot) && (pos > sum[i+1]))
				i++;

			index[p]= MIN2(tot-1, i);

			/* avoid zero weight face */
			if(p == totpart-1 && weight[index[p]] == 0.0f)
				index[p]= index[p-1];

			jitoff[index[p]]= pos;
		}
	}

	MEM_freeN(sum);

	/* for hair, sort by origindex, allows optimizations in rendering */
	/* however with virtual parents the children need to be in random order */
	if(part->type == PART_HAIR && !(part->childtype==PART_CHILD_FACES && part->parents!=0.0)) {
		if(from != PART_FROM_PARTICLE) {
			COMPARE_ORIG_INDEX = NULL;

			if(from == PART_FROM_VERT) {
				if(dm->numVertData)
					COMPARE_ORIG_INDEX= dm->getVertDataArray(dm, CD_ORIGINDEX);
			}
			else {
				if(dm->numFaceData)
					COMPARE_ORIG_INDEX= dm->getFaceDataArray(dm, CD_ORIGINDEX);
			}

			if(COMPARE_ORIG_INDEX) {
				qsort(index, totpart, sizeof(int), compare_orig_index);
				COMPARE_ORIG_INDEX = NULL;
			}
		}
	}

	/* weights are no longer used except for FROM_PARTICLE, which needs them zeroed for indexing */
	if(from==PART_FROM_PARTICLE){
		for(i=0; i<tot; i++)
			weight[i]=0.0f;
	}

	/* 4. */
	if(distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
		jitlevel= part->userjit;
		
		if(jitlevel == 0) {
			jitlevel= totpart/tot;
			if(part->flag & PART_EDISTR) jitlevel*= 2;	/* looks better in general, not very scietific */
			if(jitlevel<3) jitlevel= 3;
		}
		
		jit= MEM_callocN((2+ jitlevel*2)*sizeof(float), "jit");

		/* for small amounts of particles we use regular jitter since it looks
		 * a bit better, for larger amounts we switch to hammersley sequence 
		 * because it is much faster */
		if(jitlevel < 25)
			init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
		else
			hammersley_create(jit, jitlevel+1, psys->seed, part->jitfac);
		BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
	}

	/* 5. */
	ctx->tree= tree;
	ctx->seams= seams;
	ctx->totseam= totseam;
	ctx->sim.psys= psys;
	ctx->index= index;
	ctx->jit= jit;
	ctx->jitlevel= jitlevel;
	ctx->jitoff= jitoff;
	ctx->weight= weight;
	ctx->maxweight= maxweight;
	ctx->from= (children)? PART_FROM_CHILD: from;
	ctx->cfrom= cfrom;
	ctx->distr= distr;
	ctx->dm= dm;
	ctx->tpars= tpars;

	if(children) {
		totpart= psys_render_simplify_distribution(ctx, totpart);
		alloc_child_particles(psys, totpart);
	}

	if(!children || psys->totchild < 10000)
		totthread= 1;
	
	seed= 31415926 + ctx->sim.psys->seed;
	for(i=0; i<totthread; i++) {
		threads[i].rng= rng_new(seed);
		threads[i].tot= totthread;
	}

	return 1;
}

static void distribute_particles_on_dm(ParticleSimulationData *sim, int from)
{
	DerivedMesh *finaldm = sim->psmd->dm;
	ListBase threads;
	ParticleThread *pthreads;
	ParticleThreadContext *ctx;
	int i, totthread;

	pthreads= psys_threads_create(sim);

	if(!psys_threads_init_distribution(pthreads, sim->scene, finaldm, from)) {
		psys_threads_free(pthreads);
		return;
	}

	totthread= pthreads[0].tot;
	if(totthread > 1) {
		BLI_init_threads(&threads, exec_distribution, totthread);

		for(i=0; i<totthread; i++)
			BLI_insert_thread(&threads, &pthreads[i]);

		BLI_end_threads(&threads);
	}
	else
		exec_distribution(&pthreads[0]);

	psys_calc_dmcache(sim->ob, finaldm, sim->psys);

	ctx= pthreads[0].ctx;
	if(ctx->dm != finaldm)
		ctx->dm->release(ctx->dm);

	psys_threads_free(pthreads);
}

/* ready for future use, to emit particles without geometry */
static void distribute_particles_on_shape(ParticleSimulationData *sim, int UNUSED(from))
{
	ParticleSystem *psys = sim->psys;
	PARTICLE_P;

	fprintf(stderr,"Shape emission not yet possible!\n");

	LOOP_PARTICLES {
		pa->fuv[0]=pa->fuv[1]=pa->fuv[2]=pa->fuv[3]= 0.0;
		pa->foffset= 0.0f;
		pa->num= -1;
	}
}
static void distribute_particles(ParticleSimulationData *sim, int from)
{
	PARTICLE_PSMD;
	int distr_error=0;

	if(psmd){
		if(psmd->dm)
			distribute_particles_on_dm(sim, from);
		else
			distr_error=1;
	}
	else
		distribute_particles_on_shape(sim, from);

	if(distr_error){
		ParticleSystem *psys = sim->psys;
		PARTICLE_P;

		fprintf(stderr,"Particle distribution error!\n");

		LOOP_PARTICLES {
			pa->fuv[0]=pa->fuv[1]=pa->fuv[2]=pa->fuv[3]= 0.0;
			pa->foffset= 0.0f;
			pa->num= -1;
		}
	}
}

/* threaded child particle distribution and path caching */
ParticleThread *psys_threads_create(ParticleSimulationData *sim)
{
	ParticleThread *threads;
	ParticleThreadContext *ctx;
	int i, totthread;

	if(sim->scene->r.mode & R_FIXED_THREADS)
		totthread= sim->scene->r.threads;
	else
		totthread= BLI_system_thread_count();
	
	threads= MEM_callocN(sizeof(ParticleThread)*totthread, "ParticleThread");
	ctx= MEM_callocN(sizeof(ParticleThreadContext), "ParticleThreadContext");

	ctx->sim = *sim;
	ctx->dm= ctx->sim.psmd->dm;
	ctx->ma= give_current_material(sim->ob, sim->psys->part->omat);

	memset(threads, 0, sizeof(ParticleThread)*totthread);

	for(i=0; i<totthread; i++) {
		threads[i].ctx= ctx;
		threads[i].num= i;
		threads[i].tot= totthread;
	}

	return threads;
}

void psys_threads_free(ParticleThread *threads)
{
	ParticleThreadContext *ctx= threads[0].ctx;
	int i, totthread= threads[0].tot;

	/* 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->sim.psys->lattice){
		end_latt_deform(ctx->sim.psys->lattice);
		ctx->sim.psys->lattice= 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->skip) MEM_freeN(ctx->skip);
	if(ctx->seams) MEM_freeN(ctx->seams);
	//if(ctx->vertpart) MEM_freeN(ctx->vertpart);
	BLI_kdtree_free(ctx->tree);

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

	MEM_freeN(ctx);
	MEM_freeN(threads);
}

/* set particle parameters that don't change during particle's life */
void initialize_particle(ParticleSimulationData *sim, ParticleData *pa, int p)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	ParticleTexture ptex;

	pa->flag &= ~PARS_UNEXIST;

	if(part->from != PART_FROM_PARTICLE && part->type != PART_FLUID) {
		psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);
		
		if(ptex.exist < PSYS_FRAND(p+125))
			pa->flag |= PARS_UNEXIST;

		pa->time = (part->type == PART_HAIR) ? 0.f : part->sta + (part->end - part->sta)*ptex.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 dont have a derived mesh face */
}
static void initialize_all_particles(ParticleSimulationData *sim)
{
	ParticleSystem *psys = sim->psys;
	PARTICLE_P;

	psys->totunexist = 0;

	LOOP_PARTICLES {
		if((pa->flag & PARS_UNEXIST)==0)
			initialize_particle(sim, pa, p);

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

	/* Free unexisting particles. */
	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++;

		}
	}
}
/* sets particle to the emitter surface with initial velocity & rotation */
void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
{
	Object *ob = sim->ob;
	ParticleSystem *psys = sim->psys;
	ParticleSettings *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];
	int p = pa - psys->particles;
	part=psys->part;
	
#if 0 /* deprecated code */
	if(part->from==PART_FROM_PARTICLE){
		float speed;
		ParticleSimulationData tsim= {0};
		tsim.scene= sim->scene;
		tsim.ob= psys->target_ob ? psys->target_ob : ob;
		tsim.psys = BLI_findlink(&tsim.ob->particlesystem, sim->psys->target_psys-1);

		state.time = pa->time;
		if(pa->num == -1)
			memset(&state, 0, sizeof(state));
		else
			psys_get_particle_state(&tsim, pa->num, &state, 1);
		psys_get_from_key(&state, loc, nor, rot, 0);

		mul_qt_v3(rot, vtan);
		mul_qt_v3(rot, utan);

		speed= normalize_v3_v3(p_vel, state.vel);
		mul_v3_fl(p_vel, dot_v3v3(r_vel, p_vel));
		VECSUB(p_vel, r_vel, p_vel);
		normalize_v3(p_vel);
		mul_v3_fl(p_vel, speed);

		VECCOPY(pa->fuv, loc); /* abusing pa->fuv (not used for "from particle") for storing emit location */
	}
	else{
#endif
	/* get precise emitter matrix if particle is born */
	if(part->type!=PART_HAIR && pa->time < cfra && pa->time >= sim->psys->cfra) {
		/* we have to force RECALC_ANIM here since where_is_objec_time only does drivers */
		BKE_animsys_evaluate_animdata(&sim->ob->id, sim->ob->adt, pa->time, ADT_RECALC_ANIM);
		where_is_object_time(sim->scene, sim->ob, pa->time);
	}

	/* get birth location from object		*/
	if(part->tanfac != 0.f)
		psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0,0);
	else
		psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,0,0,0,0);
		
	/* get possible textural influence */
	psys_get_texture(sim, pa, &ptex, PAMAP_IVEL|PAMAP_LIFE, 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(part->tanfac!=0.0){
		//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,-(float)cos(M_PI*(part->tanphase+phase)));
		fac=-(float)sin(M_PI*(part->tanphase+phase));
		VECADDFAC(vtan,vtan,utan,fac);

		mul_mat3_m4_v3(ob->obmat,vtan);

		VECCOPY(utan,nor);
		mul_v3_fl(utan,dot_v3v3(vtan,nor));
		VECSUB(vtan,vtan,utan);
			
		normalize_v3(vtan);
	}
		

	/* -velocity							*/
	if(part->randfac!=0.0){
		r_vel[0] = 2.0f * (PSYS_FRAND(p + 10) - 0.5f);
		r_vel[1] = 2.0f * (PSYS_FRAND(p + 11) - 0.5f);
		r_vel[2] = 2.0f * (PSYS_FRAND(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(p + 13) - 0.5f);
		r_ave[1] = 2.0f * (PSYS_FRAND(p + 14) - 0.5f);
		r_ave[2] = 2.0f * (PSYS_FRAND(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(p + 16) - 0.5f);
		r_rot[1] = 2.0f * (PSYS_FRAND(p + 17) - 0.5f);
		r_rot[2] = 2.0f * (PSYS_FRAND(p + 18) - 0.5f);
		r_rot[3] = 2.0f * (PSYS_FRAND(p + 19) - 0.5f);
		normalize_qt(r_rot);

		mat4_to_quat(rot,ob->obmat);
		mul_qt_qtqt(r_rot,r_rot,rot);
	}
#if 0
	}
#endif

	if(part->phystype==PART_PHYS_BOIDS && pa->boid) {
		BoidParticle *bpa = pa->boid;
		float dvec[3], q[4], mat[3][3];

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

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

		/* boids store direction in ave */
		if(fabs(nor[2])==1.0f) {
			sub_v3_v3v3(pa->state.ave, loc, ob->obmat[3]);
			normalize_v3(pa->state.ave);
		}
		else {
			VECCOPY(pa->state.ave, nor);
		}
		/* 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];

		/* calculate rotation matrix */
		project_v3_v3v3(dvec, r_vel, pa->state.ave);
		sub_v3_v3v3(mat[0], pa->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(pa->state.rot, q);

		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;
	}
	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, pa->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(pa->state.vel, vel, ptex.ivel);

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

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

		if(part->rotmode){
			/* create vector into which rotation is aligned */
			switch(part->rotmode){
				case PART_ROT_NOR:
					copy_v3_v3(rot_vec, nor);
					break;
				case PART_ROT_VEL:
					copy_v3_v3(rot_vec, vel);
					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;
					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]);
					break;
			}
			
			/* create rotation quat */
			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);

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

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

		/* -angular velocity					*/

		zero_v3(pa->state.ave);

		if(part->avemode){
			switch(part->avemode){
				case PART_AVE_SPIN:
					copy_v3_v3(pa->state.ave, vel);
					break;
				case PART_AVE_RAND:
					copy_v3_v3(pa->state.ave, r_ave);
					break;
			}
			normalize_v3(pa->state.ave);
			mul_v3_fl(pa->state.ave,part->avefac);
		}
	}


	if(part->type == PART_HAIR){
		pa->lifetime = 100.0f;
	}
	else{
		pa->lifetime = part->lifetime * ptex.life;

		if(part->randlife != 0.0)
			pa->lifetime *= 1.0f - part->randlife * PSYS_FRAND(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;

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

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

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

	psys->flag |= PSYS_KEYED;
}
/************************************************/
/*			Reactors							*/
/************************************************/
//static void push_reaction(ParticleSimulationData *sim, int pa_num, int event, ParticleKey *state)
//{
//	Object *rob;
//	ParticleSystem *rpsys;
//	ParticleSettings *rpart;
//	ParticleData *pa;
//	ListBase *lb=&sim->psys->effectors;
//	ParticleEffectorCache *ec;
//	ParticleReactEvent *re;
//
//	if(lb->first) for(ec = lb->first; ec; ec= ec->next){
//		if(ec->type & PSYS_EC_REACTOR){
//			/* all validity checks already done in add_to_effectors */
//			rob=ec->ob;
//			rpsys=BLI_findlink(&rob->particlesystem,ec->psys_nbr);
//			rpart=rpsys->part;
//			if(rpsys->part->reactevent==event){
//				pa=sim->psys->particles+pa_num;
//				re= MEM_callocN(sizeof(ParticleReactEvent), "react event");
//				re->event=event;
//				re->pa_num = pa_num;
//				re->ob = sim->ob;
//				re->psys = sim->psys;
//				re->size = pa->size;
//				copy_particle_key(&re->state,state,1);
//
//				switch(event){
//					case PART_EVENT_DEATH:
//						re->time=pa->dietime;
//						break;
//					case PART_EVENT_COLLIDE:
//						re->time=state->time;
//						break;
//					case PART_EVENT_NEAR:
//						re->time=state->time;
//						break;
//				}
//
//				BLI_addtail(&rpsys->reactevents, re);
//			}
//		}
//	}
//}
//static void react_to_events(ParticleSystem *psys, int pa_num)
//{
//	ParticleSettings *part=psys->part;
//	ParticleData *pa=psys->particles+pa_num;
//	ParticleReactEvent *re=psys->reactevents.first;
//	int birth=0;
//	float dist=0.0f;
//
//	for(re=psys->reactevents.first; re; re=re->next){
//		birth=0;
//		if(part->from==PART_FROM_PARTICLE){
//			if(pa->num==re->pa_num && pa->alive==PARS_UNBORN){
//				if(re->event==PART_EVENT_NEAR){
//					ParticleData *tpa = re->psys->particles+re->pa_num;
//					float pa_time=tpa->time + pa->foffset*tpa->lifetime;
//					if(re->time >= pa_time){
//						pa->time=pa_time;
//						pa->dietime=pa->time+pa->lifetime;
//					}
//				}
//				else{
//					pa->time=re->time;
//					pa->dietime=pa->time+pa->lifetime;
//				}
//			}
//		}
//		else{
//			dist=len_v3v3(pa->state.co, re->state.co);
//			if(dist <= re->size){
//				if(pa->alive==PARS_UNBORN){
//					pa->time=re->time;
//					pa->dietime=pa->time+pa->lifetime;
//					birth=1;
//				}
//				if(birth || part->flag&PART_REACT_MULTIPLE){
//					float vec[3];
//					VECSUB(vec,pa->state.co, re->state.co);
//					if(birth==0)
//						mul_v3_fl(vec,(float)pow(1.0f-dist/re->size,part->reactshape));
//					VECADDFAC(pa->state.vel,pa->state.vel,vec,part->reactfac);
//					VECADDFAC(pa->state.vel,pa->state.vel,re->state.vel,part->partfac);
//				}
//				if(birth)
//					mul_v3_fl(pa->state.vel,(float)pow(1.0f-dist/re->size,part->reactshape));
//			}
//		}
//	}
//}
//void psys_get_reactor_target(ParticleSimulationData *sim, Object **target_ob, ParticleSystem **target_psys)
//{
//	Object *tob;
//
//	tob = sim->psys->target_ob ? sim->psys->target_ob : sim->ob;
//	
//	*target_psys = BLI_findlink(&tob->particlesystem, sim->psys->target_psys-1);
//	if(*target_psys)
//		*target_ob=tob;
//	else
//		*target_ob=0;
//}
/************************************************/
/*			Point Cache							*/
/************************************************/
void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
{
	PointCache *cache = psys->pointcache;

	if(cache->flag & PTCACHE_DISK_CACHE && cache->mem_cache.first == NULL) {
		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 = MAX2(1, (int)part->sta);
	*efra = MIN2((int)(part->end + part->lifetime + 1.0), scene->r.efra);
}

/************************************************/
/*			Effectors							*/
/************************************************/
void psys_update_particle_tree(ParticleSystem *psys, float cfra)
{
	if(psys) {
		PARTICLE_P;

		if(!psys->tree || psys->tree_frame != cfra) {
			
			BLI_kdtree_free(psys->tree);

			psys->tree = BLI_kdtree_new(psys->totpart);
			
			LOOP_SHOWN_PARTICLES {
				if(pa->alive == PARS_ALIVE) {
					if(pa->state.time == cfra)
						BLI_kdtree_insert(psys->tree, p, pa->prev_state.co, NULL);
					else
						BLI_kdtree_insert(psys->tree, p, pa->state.co, NULL);
				}
			}
			BLI_kdtree_balance(psys->tree);

			psys->tree_frame = psys->cfra;
		}
	}
}

static void psys_update_effectors(ParticleSimulationData *sim)
{
	pdEndEffectors(&sim->psys->effectors);
	sim->psys->effectors = pdInitEffectors(sim->scene, sim->ob, sim->psys, sim->psys->part->effector_weights);
	precalc_guides(sim, sim->psys->effectors);
}

/*********************************************************************************************************
                    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 *add_fluid_spring(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  delete_fluid_spring(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 EdgeHash *build_fluid_springhash(ParticleSystem *psys)
{
	EdgeHash *springhash = NULL;
	ParticleSpring *spring;
	int i = 0;

	springhash = BLI_edgehash_new();

	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], SET_INT_IN_POINTER(i+1));

	return springhash;
}
static void particle_fluidsim(ParticleSystem *psys, int own_psys, ParticleData *pa, float dtime, float mass, float *gravity, EdgeHash *springhash)
{
	SPHFluidSettings *fluid = psys->part->fluid;
	KDTreeNearest *ptn = NULL;
	ParticleData *npa;
	ParticleSpring *spring = NULL;

	float temp[3];
	float q, q1, u, I, D, rij, d, Lij;
	float pressure_near, pressure;
	float p=0, pnear=0;

	float omega = fluid->viscosity_omega;
	float beta = fluid->viscosity_beta;
	float massfactor = 1.0f/mass;
	float spring_k = fluid->spring_k;
	float h = fluid->radius;
	float L = fluid->rest_length * fluid->radius;

	int n, neighbours = BLI_kdtree_range_search(psys->tree, h, pa->prev_state.co, NULL, &ptn);
	int spring_index = 0, index = own_psys ? pa - psys->particles : -1;

	/* pressure and near pressure */
	for(n=own_psys?1:0; n<neighbours; n++) {
		/* disregard particles at the exact same location */
		if(ptn[n].dist < FLT_EPSILON)
			continue;

		sub_v3_v3(ptn[n].co, pa->prev_state.co);
		mul_v3_fl(ptn[n].co, 1.f/ptn[n].dist);
		q = ptn[n].dist/h;

		if(q < 1.f) {
			q1 = 1.f - q;

			p += q1*q1;
			pnear += q1*q1*q1;
		}
	}

	p *= mass;
	pnear *= mass;
	pressure =  fluid->stiffness_k * (p - fluid->rest_density);
	pressure_near = fluid->stiffness_knear * pnear;

	/* main calculations */
	for(n=own_psys?1:0; n<neighbours; n++) {
		/* disregard particles at the exact same location */
		if(ptn[n].dist < FLT_EPSILON)
			continue;

		npa = psys->particles + ptn[n].index;

		rij = ptn[n].dist;
		q = rij/h;
		q1 = 1.f-q;

		/* Double Density Relaxation - Algorithm 2 (can't be thread safe!)*/
		D =  dtime * dtime * (pressure + pressure_near*q1)*q1 * 0.5f;
		madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
		if(own_psys)
			madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);

		if(index < ptn[n].index) {
			/* Viscosity - Algorithm 5 */
			if(omega > 0.f	|| beta > 0.f) {		
				sub_v3_v3v3(temp, pa->state.vel, npa->state.vel);
				u = dot_v3v3(ptn[n].co, temp);

				if (u > 0){
					I = dtime * (q1 * (omega * u + beta * u*u)) * 0.5f;
					madd_v3_v3fl(pa->state.vel, ptn[n].co, -I * massfactor);

					if(own_psys)
						madd_v3_v3fl(npa->state.vel, ptn[n].co, I * massfactor);
				}
			}

			if(spring_k > 0.f) {
				/* Viscoelastic spring force - Algorithm 4*/
				if (fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash){
					spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, ptn[n].index));

					if(spring_index) {
						spring = psys->fluid_springs + spring_index - 1;
					}
					else {
						ParticleSpring temp_spring;
						temp_spring.particle_index[0] = index;
						temp_spring.particle_index[1] = ptn[n].index;
						temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : L;
						temp_spring.delete_flag = 0;
						
						spring = add_fluid_spring(psys, &temp_spring);
					}

					Lij = spring->rest_length;
					d = fluid->yield_ratio * Lij;

					if (rij > Lij + d) // Stretch, 25 is just a multiplier for plasticity_constant value to counter default dtime of 1/25
						spring->rest_length += dtime * 25.f * fluid->plasticity_constant * (rij - Lij - d);
					else if(rij < Lij - d) // Compress
						spring->rest_length -= dtime * 25.f * fluid->plasticity_constant * (Lij - d - rij);
				}
				else { /* PART_SPRING_HOOKES - Hooke's spring force */
					/* L is a factor of radius */
					D = 0.5 * dtime * dtime * 10.f * fluid->spring_k * (1.f - L/h) * (L - rij);
 
					madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
					if(own_psys)
						madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
				}
			}
		}
	} 

	/* Artificial buoyancy force in negative gravity direction  */
	if (fluid->buoyancy >= 0.f && gravity) {
		float B = -dtime * dtime * fluid->buoyancy * (p - fluid->rest_density) * 0.5f;
		madd_v3_v3fl(pa->state.co, gravity, -B * massfactor);
	}

	if(ptn)
		MEM_freeN(ptn);
}

static void apply_particle_fluidsim(Object *ob, ParticleSystem *psys, ParticleData *pa, float dtime, float *gravity, EdgeHash *springhash){
	ParticleTarget *pt;

	particle_fluidsim(psys, 1, pa, dtime, psys->part->mass, gravity, springhash);
	
	/*----check other SPH systems (Multifluids) , each fluid has its own parameters---*/
	for(pt=psys->targets.first; pt; pt=pt->next) {
		ParticleSystem *epsys = psys_get_target_system(ob, pt);

		if(epsys)
			particle_fluidsim(epsys, 0, pa, dtime, psys->part->mass, gravity, NULL);
	}
	/*----------------------------------------------------------------*/	 	 
}

static void apply_fluid_springs(ParticleSystem *psys, float timestep){
	SPHFluidSettings *fluid = psys->part->fluid;
	ParticleData *pa1, *pa2;
	ParticleSpring *spring = psys->fluid_springs;
	
	float h = fluid->radius;
	float massfactor = 1.0f/psys->part->mass;
	float D, Rij[3], rij, Lij;
	int i;

	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++) {
		Lij = spring->rest_length;

		if (Lij > h) {
			spring->delete_flag = 1;
		}
		else {
			pa1 = psys->particles + spring->particle_index[0];
			pa2 = psys->particles + spring->particle_index[1];

			sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
			rij = normalize_v3(Rij);

			/* Calculate displacement and apply value */
			D =  0.5f * timestep * timestep * 10.f * fluid->spring_k * (1.f - Lij/h) * (Lij - rij);

			madd_v3_v3fl(pa1->state.co, Rij, -D * pa1->state.time * pa1->state.time * massfactor);
			madd_v3_v3fl(pa2->state.co, Rij, D * pa2->state.time * pa2->state.time * massfactor);
		}
	}

	/* Loop through springs backwaqrds - for efficient delete function */
	for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
		if(psys->fluid_springs[i].delete_flag)
			delete_fluid_spring(psys, i);
	}
}

/************************************************/
/*			Newtonian physics					*/
/************************************************/
/* gathers all forces that effect particles and calculates a new state for the particle */
static void apply_particle_forces(ParticleSimulationData *sim, int p, float dfra, float cfra)
{
	ParticleSettings *part = sim->psys->part;
	ParticleData *pa = sim->psys->particles + p;
	EffectedPoint epoint;
	ParticleKey states[5], tkey;
	float timestep = psys_get_timestep(sim);
	float force[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
	float dtime=dfra*timestep, time, pa_mass=part->mass, fac /*, fra=sim->psys->cfra*/;
	int i, steps=1;
	ParticleTexture ptex;

	psys_get_texture(sim, pa, &ptex, PAMAP_PHYSICS, cfra);
	
	/* maintain angular velocity */
	VECCOPY(pa->state.ave,pa->prev_state.ave);
	VECCOPY(oldpos,pa->state.co);

	if(part->flag & PART_SIZEMASS)
		pa_mass*=pa->size;

	switch(part->integrator){
		case PART_INT_EULER:
			steps=1;
			break;
		case PART_INT_MIDPOINT:
			steps=2;
			break;
		case PART_INT_RK4:
			steps=4;
			break;
		case PART_INT_VERLET:
			steps=1;
			break;
	}

	copy_particle_key(states,&pa->state,1);

	for(i=0; i<steps; i++){
		force[0]=force[1]=force[2]=0.0;
		impulse[0]=impulse[1]=impulse[2]=0.0;
		/* add effectors */
		pd_point_from_particle(sim, pa, states+i, &epoint);
		if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
			pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);

		mul_v3_fl(force, ptex.field);
		mul_v3_fl(impulse, ptex.field);

		/* calculate air-particle interaction */
		if(part->dragfac!=0.0f){
			fac=-part->dragfac*pa->size*pa->size*len_v3(states[i].vel);
			VECADDFAC(force,force,states[i].vel,fac);
		}

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

		/* force to acceleration*/
		mul_v3_fl(force,1.0f/pa_mass);

		/* add global acceleration (gravitation) */
		if(psys_uses_gravity(sim)
			/* normal gravity is too strong for hair so it's disabled by default */
			&& (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
			madd_v3_v3fl(force, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * ptex.gravity);
		}
		
		/* calculate next state */
		VECADD(states[i].vel,states[i].vel,impulse);

		switch(part->integrator){
			case PART_INT_EULER:
				VECADDFAC(pa->state.co,states->co,states->vel,dtime);
				VECADDFAC(pa->state.vel,states->vel,force,dtime);
				break;
			case PART_INT_MIDPOINT:
				if(i==0){
					VECADDFAC(states[1].co,states->co,states->vel,dtime*0.5f);
					VECADDFAC(states[1].vel,states->vel,force,dtime*0.5f);
					/*fra=sim->psys->cfra+0.5f*dfra;*/
				}
				else{
					VECADDFAC(pa->state.co,states->co,states[1].vel,dtime);
					VECADDFAC(pa->state.vel,states->vel,force,dtime);
				}
				break;
			case PART_INT_RK4:
				switch(i){
					case 0:
						VECCOPY(dx[0],states->vel);
						mul_v3_fl(dx[0],dtime);
						VECCOPY(dv[0],force);
						mul_v3_fl(dv[0],dtime);

						VECADDFAC(states[1].co,states->co,dx[0],0.5f);
						VECADDFAC(states[1].vel,states->vel,dv[0],0.5f);
						/*fra=sim->psys->cfra+0.5f*dfra;*/
						break;
					case 1:
						VECADDFAC(dx[1],states->vel,dv[0],0.5f);
						mul_v3_fl(dx[1],dtime);
						VECCOPY(dv[1],force);
						mul_v3_fl(dv[1],dtime);

						VECADDFAC(states[2].co,states->co,dx[1],0.5f);
						VECADDFAC(states[2].vel,states->vel,dv[1],0.5f);
						break;
					case 2:
						VECADDFAC(dx[2],states->vel,dv[1],0.5f);
						mul_v3_fl(dx[2],dtime);
						VECCOPY(dv[2],force);
						mul_v3_fl(dv[2],dtime);

						VECADD(states[3].co,states->co,dx[2]);
						VECADD(states[3].vel,states->vel,dv[2]);
						/*fra=cfra;*/
						break;
					case 3:
						VECADD(dx[3],states->vel,dv[2]);
						mul_v3_fl(dx[3],dtime);
						VECCOPY(dv[3],force);
						mul_v3_fl(dv[3],dtime);

						VECADDFAC(pa->state.co,states->co,dx[0],1.0f/6.0f);
						VECADDFAC(pa->state.co,pa->state.co,dx[1],1.0f/3.0f);
						VECADDFAC(pa->state.co,pa->state.co,dx[2],1.0f/3.0f);
						VECADDFAC(pa->state.co,pa->state.co,dx[3],1.0f/6.0f);

						VECADDFAC(pa->state.vel,states->vel,dv[0],1.0f/6.0f);
						VECADDFAC(pa->state.vel,pa->state.vel,dv[1],1.0f/3.0f);
						VECADDFAC(pa->state.vel,pa->state.vel,dv[2],1.0f/3.0f);
						VECADDFAC(pa->state.vel,pa->state.vel,dv[3],1.0f/6.0f);
				}
				break;
			case PART_INT_VERLET:   /* Verlet integration */
				VECADDFAC(pa->state.vel,pa->state.vel,force,dtime);
				VECADDFAC(pa->state.co,pa->state.co,pa->state.vel,dtime);

				VECSUB(pa->state.vel,pa->state.co,oldpos);
				mul_v3_fl(pa->state.vel,1.0f/dtime);
				break;
		}
	}

	/* damp affects final velocity */
	if(part->dampfac != 0.f)
		mul_v3_fl(pa->state.vel, 1.f - part->dampfac * ptex.damp);

	VECCOPY(pa->state.ave, states->ave);

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

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

	if(part->type != PART_HAIR) {
		if(do_guides(sim->psys->effectors, &tkey, p, time)) {
			VECCOPY(pa->state.co,tkey.co);
			/* guides don't produce valid velocity */
			VECSUB(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 rotate_particle(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;

	if((part->flag & PART_ROT_DYN)==0){
		if(part->avemode==PART_AVE_SPIN){
			float angle;
			float len1 = len_v3(pa->prev_state.vel);
			float len2 = len_v3(pa->state.vel);

			if(len1==0.0f || len2==0.0f)
				pa->state.ave[0]=pa->state.ave[1]=pa->state.ave[2]=0.0f;
			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);
			}

			axis_angle_to_quat(rot2,pa->state.vel,dtime*part->avefac);
		}
	}

	rotfac=len_v3(pa->state.ave);
	if(rotfac==0.0){ /* unit_qt(in VecRotToQuat) doesn't give unit quat [1,0,0,0]?? */
		rot1[0]=1.0;
		rot1[1]=rot1[2]=rot1[3]=0;
	}
	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);
}

/* convert from triangle barycentric weights to quad mean value weights */
static void intersect_dm_quad_weights(float *v1, float *v2, float *v3, float *v4, float *w)
{
	float co[3], vert[4][3];

	VECCOPY(vert[0], v1);
	VECCOPY(vert[1], v2);
	VECCOPY(vert[2], v3);
	VECCOPY(vert[3], v4);

	co[0]= v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2] + v4[0]*w[3];
	co[1]= v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2] + v4[1]*w[3];
	co[2]= v1[2]*w[0] + v2[2]*w[1] + v3[2]*w[2] + v4[2]*w[3];

	interp_weights_poly_v3( w,vert, 4, co);
}

/* check intersection with a derivedmesh */
int psys_intersect_dm(Scene *scene, Object *ob, DerivedMesh *dm, float *vert_cos, float *co1, float* co2, float *min_d, int *min_face, float *min_w,
						  float *face_minmax, float *pa_minmax, float radius, float *ipoint)
{
	MFace *mface=0;
	MVert *mvert=0;
	int i, totface, intersect=0;
	float cur_d, cur_uv[2], v1[3], v2[3], v3[3], v4[3], min[3], max[3], p_min[3],p_max[3];
	float cur_ipoint[3];
	
	if(dm==0){
		psys_disable_all(ob);

		dm=mesh_get_derived_final(scene, ob, 0);
		if(dm==0)
			dm=mesh_get_derived_deform(scene, ob, 0);

		psys_enable_all(ob);

		if(dm==0)
			return 0;
	}

	

	if(pa_minmax==0){
		INIT_MINMAX(p_min,p_max);
		DO_MINMAX(co1,p_min,p_max);
		DO_MINMAX(co2,p_min,p_max);
	}
	else{
		VECCOPY(p_min,pa_minmax);
		VECCOPY(p_max,pa_minmax+3);
	}

	totface=dm->getNumFaces(dm);
	mface=dm->getFaceDataArray(dm,CD_MFACE);
	mvert=dm->getVertDataArray(dm,CD_MVERT);
	
	/* lets intersect the faces */
	for(i=0; i<totface; i++,mface++){
		if(vert_cos){
			VECCOPY(v1,vert_cos+3*mface->v1);
			VECCOPY(v2,vert_cos+3*mface->v2);
			VECCOPY(v3,vert_cos+3*mface->v3);
			if(mface->v4)
				VECCOPY(v4,vert_cos+3*mface->v4)
		}
		else{
			VECCOPY(v1,mvert[mface->v1].co);
			VECCOPY(v2,mvert[mface->v2].co);
			VECCOPY(v3,mvert[mface->v3].co);
			if(mface->v4)
				VECCOPY(v4,mvert[mface->v4].co)
		}

		if(face_minmax==0){
			INIT_MINMAX(min,max);
			DO_MINMAX(v1,min,max);
			DO_MINMAX(v2,min,max);
			DO_MINMAX(v3,min,max);
			if(mface->v4)
				DO_MINMAX(v4,min,max)
			if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
				continue;
		}
		else{
			VECCOPY(min, face_minmax+6*i);
			VECCOPY(max, face_minmax+6*i+3);
			if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
				continue;
		}

		if(radius>0.0f){
			if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v2, v3, v1, &cur_d, cur_ipoint)){
				if(cur_d<*min_d){
					*min_d=cur_d;
					VECCOPY(ipoint,cur_ipoint);
					*min_face=i;
					intersect=1;
				}
			}
			if(mface->v4){
				if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v4, v1, v3, &cur_d, cur_ipoint)){
					if(cur_d<*min_d){
						*min_d=cur_d;
						VECCOPY(ipoint,cur_ipoint);
						*min_face=i;
						intersect=1;
					}
				}
			}
		}
		else{
			if(isect_line_tri_v3(co1, co2, v1, v2, v3, &cur_d, cur_uv)){
				if(cur_d<*min_d){
					*min_d=cur_d;
					min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
					min_w[1]= cur_uv[0];
					min_w[2]= cur_uv[1];
					min_w[3]= 0.0f;
					if(mface->v4)
						intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
					*min_face=i;
					intersect=1;
				}
			}
			if(mface->v4){
				if(isect_line_tri_v3(co1, co2, v1, v3, v4, &cur_d, cur_uv)){
					if(cur_d<*min_d){
						*min_d=cur_d;
						min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
						min_w[1]= 0.0f;
						min_w[2]= cur_uv[0];
						min_w[3]= cur_uv[1];
						intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
						*min_face=i;
						intersect=1;
					}
				}
			}
		}
	}
	return intersect;
}

void particle_intersect_face(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
{
	ParticleCollision *col = (ParticleCollision *) userdata;
	MFace *face = col->md->mfaces + index;
	MVert *x = col->md->x;
	MVert *v = col->md->current_v;
	float vel[3], co1[3], co2[3], uv[2], ipoint[3], temp[3], t;
	float x0[3], x1[3], x2[3], x3[3];
	float *t0=x0, *t1=x1, *t2=x2, *t3=(face->v4 ? x3 : NULL);

	/* move collision face to start of timestep */
	madd_v3_v3v3fl(t0, x[face->v1].co, v[face->v1].co, col->cfra);
	madd_v3_v3v3fl(t1, x[face->v2].co, v[face->v2].co, col->cfra);
	madd_v3_v3v3fl(t2, x[face->v3].co, v[face->v3].co, col->cfra);
	if(t3)
		madd_v3_v3v3fl(t3, x[face->v4].co, v[face->v4].co, col->cfra);

	/* calculate average velocity of face */
	copy_v3_v3(vel, v[ face->v1 ].co);
	add_v3_v3(vel, v[ face->v2 ].co);
	add_v3_v3(vel, v[ face->v3 ].co);
	mul_v3_fl(vel, 0.33334f*col->dfra);

	/* substract face velocity, in other words convert to 
	   a coordinate system where only the particle moves */
	madd_v3_v3v3fl(co1, col->co1, vel, -col->f);
	sub_v3_v3v3(co2, col->co2, vel);

	do
	{	
		if(ray->radius == 0.0f) {
			if(isect_line_tri_v3(co1, co2, t0, t1, t2, &t, uv)) {
				if(t >= 0.0f && t < hit->dist/col->ray_len) {
					hit->dist = col->ray_len * t;
					hit->index = index;

					/* calculate normal that's facing the particle */
					normal_tri_v3( col->nor,t0, t1, t2);
					VECSUB(temp, co2, co1);
					if(dot_v3v3(col->nor, temp) > 0.0f)
						negate_v3(col->nor);

					VECCOPY(col->vel,vel);

					col->hit_ob = col->ob;
					col->hit_md = col->md;
				}
			}
		}
		else {
			if(isect_sweeping_sphere_tri_v3(co1, co2, ray->radius, t0, t1, t2, &t, ipoint)) {
				if(t >=0.0f && t < hit->dist/col->ray_len) {
					hit->dist = col->ray_len * t;
					hit->index = index;

					interp_v3_v3v3(temp, co1, co2, t);
					
					VECSUB(col->nor, temp, ipoint);
					normalize_v3(col->nor);

					VECCOPY(col->vel,vel);

					col->hit_ob = col->ob;
					col->hit_md = col->md;
				}
			}
		}

		t1 = t2;
		t2 = t3;
		t3 = NULL;

	} while(t2);
}
/* Particle - Mesh collision code
 * Features:
 * - point and swept sphere to mesh surface collisions
 * - moving colliders (but not yet rotating or deforming colliders)
 * - friction & damping
 * - angular momentum <-> linear momentum
 * - high accuracy by re-applying particle acceleration after collision
 * - behaves relatively well even if limit of 10 collisions per simulation step is exceeded
 * Main parts:
 * 1. check for all possible deflectors for closest intersection on particle path
 * 2. if deflection was found calculate new coordinates or kill the particle
 */
static void deflect_particle(ParticleSimulationData *sim, int p, float dfra, float cfra){
	Object *ground_ob = NULL;
	ParticleSettings *part = sim->psys->part;
	ParticleData *pa = sim->psys->particles + p;
	ParticleCollision col;
	ColliderCache *coll;
	BVHTreeRayHit hit;
	float ray_dir[3], acc[3];
	float radius = ((part->flag & PART_SIZE_DEFL)?pa->size:0.0f), boid_z = 0.0f;
	float timestep = psys_get_timestep(sim) * dfra;
	float inv_timestep = 1.0f/timestep;
	int deflections=0, max_deflections=10;

	/* get acceleration (from gravity, forcefields etc. to be re-applied after collision) */
	sub_v3_v3v3(acc, pa->state.vel, pa->prev_state.vel);
	mul_v3_fl(acc, inv_timestep);

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

	/* override for boids */
	if(part->phystype == PART_PHYS_BOIDS) {
		BoidParticle *bpa = pa->boid;
		radius = pa->size;
		boid_z = pa->state.co[2];
		ground_ob = bpa->ground;
	}

	/* 10 iterations to catch multiple deflections */
	if(sim->colliders) while(deflections < max_deflections){
		/* 1. */

		sub_v3_v3v3(ray_dir, col.co2, col.co1);
		hit.index = -1;
		hit.dist = col.ray_len = len_v3(ray_dir);

		col.cfra = fmod(cfra-dfra, 1.0f);
		col.dfra = dfra;

		/* even if particle is stationary we want to check for moving colliders */
		/* if hit.dist is zero the bvhtree_ray_cast will just ignore everything */
		if(hit.dist == 0.0f)
			hit.dist = col.ray_len = 0.000001f;

		for(coll = sim->colliders->first; coll; coll=coll->next){
			/* for boids: don't check with current ground object */
			if(coll->ob == ground_ob)
				continue;

			/* particles should not collide with emitter at birth */
			if(coll->ob == sim->ob && pa->time < cfra && pa->time >= sim->psys->cfra)
				continue;

			col.ob = coll->ob;
			col.md = coll->collmd;

			if(col.md && col.md->bvhtree)
				BLI_bvhtree_ray_cast(col.md->bvhtree, col.co1, ray_dir, radius, &hit, particle_intersect_face, &col);
		}

		/* 2. */
		if(hit.index>=0) {
			PartDeflect *pd = col.hit_ob->pd;
			float co[3];							/* point of collision */
			float x = hit.dist/col.ray_len;			/* location factor of collision between this iteration */
			float f = col.f + x * (1.0f - col.f);	/* time factor of collision between timestep */
			float dt1 = (f - col.f) * timestep;		/* time since previous collision (in seconds) */
			float dt2 = (1.0f - f) * timestep;		/* time left after collision (in seconds) */
			int through = (BLI_frand() < pd->pdef_perm) ? 1 : 0; /* did particle pass through the collision surface? */

			deflections++;
			
			interp_v3_v3v3(co, col.co1, col.co2, x);
			
			/* make sure we don't hit the current face again */
			/* TODO: could/should this be proportional to pa->size? */
			madd_v3_v3fl(co, col.nor, (through ? -0.0001f : 0.0001f));

			/* particle dies in collision */
			if(through == 0 && (part->flag & PART_DIE_ON_COL || pd->flag & PDEFLE_KILL_PART)) {
				pa->alive = PARS_DYING;
				pa->dietime = sim->psys->cfra + (cfra - sim->psys->cfra) * f;

				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;
			}
			/* figure out velocity and other data after collision */
			else {
				float v0[3];	/* velocity directly before collision to be modified into velocity directly after collision */
				float v0_nor[3];/* normal component of v0 */
				float v0_tan[3];/* tangential component of v0 */
				float vc_tan[3];/* tangential component of collision surface velocity */
				float check[3];
				float v0_dot, vc_dot, check_dot;
				float damp, frict;

				/* get exact velocity right before collision */
				madd_v3_v3v3fl(v0, col.ve1, acc, dt1);
				
				/* convert collider velocity from 1/framestep to 1/s */
				mul_v3_fl(col.vel, inv_timestep);
				
				/* get damping & friction factors */
				damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_frand() - 0.5f);
				CLAMP(damp,0.0,1.0);

				frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_frand() - 0.5f);
				CLAMP(frict,0.0,1.0);

				/* treat normal & tangent components separately */
				v0_dot = dot_v3v3(col.nor, v0);
				madd_v3_v3v3fl(v0_tan, v0, col.nor, -v0_dot);

				vc_dot = dot_v3v3(col.nor, col.vel);
				madd_v3_v3v3fl(vc_tan, col.vel, col.nor, -vc_dot);

				/* handle friction effects (tangential and angular velocity) */
				if(frict > 0.0f) {
					/* angular <-> linear velocity */
					if(part->flag & PART_ROT_DYN) {
						float vr_tan[3], v1_tan[3], ave[3];
					
						/* linear velocity of particle surface */
						cross_v3_v3v3(vr_tan, col.nor, pa->state.ave);
						mul_v3_fl(vr_tan, pa->size);

						/* 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 */
						mul_v3_v3fl(vr_tan, v1_tan, -1.0f);

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

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

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

				/* stickness was possibly added before, so cancel that before calculating new normal velocity */
				/* otherwise particles go flying out of the surface because of high reversed sticky velocity */
				if(v0_dot < 0.0f) {
					v0_dot += pd->pdef_stickness;
					if(v0_dot > 0.0f)
						v0_dot = 0.0f;
				}

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

				/* special case for object hitting the particle from behind */
				if(through==0 && ((vc_dot>0.0f && v0_dot>0.0f && vc_dot>v0_dot) || (vc_dot<0.0f && v0_dot<0.0f && vc_dot<v0_dot)))
					mul_v3_v3fl(v0_nor, col.nor, vc_dot);
				else
					mul_v3_v3fl(v0_nor, col.nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);

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

				/* keep boids above ground */
				if(part->phystype == PART_PHYS_BOIDS && part->boids->options & BOID_ALLOW_LAND) {
					BoidParticle *bpa = pa->boid;
					if(bpa->data.mode == eBoidMode_OnLand || co[2] <= boid_z) {
						co[2] = boid_z;
						v0[2] = 0.0f;
					}
				}
				
				if(deflections < max_deflections) {
					/* re-apply acceleration to final velocity and location */
					madd_v3_v3v3fl(pa->state.vel, v0, acc, dt2);
					madd_v3_v3v3fl(pa->state.co, co, v0, dt2);
					madd_v3_v3fl(pa->state.co, acc, 0.5f*dt2*dt2);

					/* make sure particle stays on the right side of the surface */
					sub_v3_v3v3(check, pa->state.co, co);
					/* (collision surface has moved during the time too) */
					madd_v3_v3fl(check, col.vel, -dt2);

					check_dot = dot_v3v3(check, col.nor);
					if((!through && check_dot < 0.0f) || (through && check_dot > 0.0f))
						madd_v3_v3fl(pa->state.co, col.nor, (through ? -0.0001f : 0.0001f) - check_dot);

					/* Stickness to surface */
					madd_v3_v3fl(pa->state.vel, col.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;
				}
				else {
					/* final chance to prevent failure, so stick to the surface and hope for the best */
					madd_v3_v3v3fl(pa->state.co, co, col.vel, dt2);
					copy_v3_v3(pa->state.vel, v0);
				}
			}
		}
		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)
{
	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 != get_psys_tot_child(sim->scene, psys)) || psys->recalc&PSYS_RECALC_RESET)
		alloc=1;

	if(alloc || psys->recalc&PSYS_RECALC_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(get_psys_tot_child(sim->scene, psys)) {
			/* 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->from!=PART_FROM_PARTICLE && part->childtype==PART_CHILD_FACES && part->parents!=0.0)
					psys_find_parents(sim);
			}
		}
		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(!psys->renderdata) {
		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->scene, 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);

		/* 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);
		}
	}
	else if(psys->pathcache)
		psys_free_path_cache(psys, NULL);
}

static void do_hair_dynamics(ParticleSimulationData *sim)
{
	ParticleSystem *psys = sim->psys;
	DerivedMesh *dm = psys->hair_in_dm;
	MVert *mvert = NULL;
	MEdge *medge = NULL;
	MDeformVert *dvert = NULL;
	HairKey *key;
	PARTICLE_P;
	int totpoint = 0;
	int totedge;
	int k;
	float hairmat[4][4];

	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_GOAL|CLOTH_SIMSETTINGS_FLAG_NO_SPRING_COMPRESS;
		psys->clmd->coll_parms->flags &= ~CLOTH_COLLSETTINGS_FLAG_SELF;
	}

	/* create a dm from hair vertices */
	LOOP_PARTICLES
		totpoint += pa->totkey;

	totedge = totpoint;
	totpoint += psys->totpart;

	if(dm && (totpoint != dm->getNumVerts(dm) || totedge != dm->getNumEdges(dm))) {
		dm->release(dm);
		dm = psys->hair_in_dm = NULL;
	}

	if(!dm) {
		dm = psys->hair_in_dm = CDDM_new(totpoint, totedge, 0);
		DM_add_vert_layer(dm, CD_MDEFORMVERT, CD_CALLOC, NULL);
	}

	mvert = CDDM_get_verts(dm);
	medge = CDDM_get_edges(dm);
	dvert = DM_get_vert_data_layer(dm, CD_MDEFORMVERT);

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

	/* make vgroup for pin roots etc.. */
	psys->particles->hair_index = 1;
	LOOP_PARTICLES {
		if(p)
			pa->hair_index = (pa-1)->hair_index + (pa-1)->totkey + 1;

		psys_mat_hair_to_object(sim->ob, sim->psmd->dm, psys->part->from, pa, hairmat);

		for(k=0, key=pa->hair; k<pa->totkey; k++,key++) {
			
			/* create fake root before actual root to resist bending */
			if(k==0) {
				float temp[3];
				VECSUB(temp, key->co, (key+1)->co);
				VECCOPY(mvert->co, key->co);
				VECADD(mvert->co, mvert->co, temp);
				mul_m4_v3(hairmat, mvert->co);
				mvert++;

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

				if(dvert) {
					if(!dvert->totweight) {
						dvert->dw = MEM_callocN (sizeof(MDeformWeight), "deformWeight");
						dvert->totweight = 1;
					}

					dvert->dw->weight = 1.0f;
					dvert++;
				}
			}

			VECCOPY(mvert->co, key->co);
			mul_m4_v3(hairmat, mvert->co);
			mvert++;
			
			if(k) {
				medge->v1 = pa->hair_index + k - 1;
				medge->v2 = pa->hair_index + k;
				medge++;
			}

			if(dvert) {
				if(!dvert->totweight) {
					dvert->dw = MEM_callocN (sizeof(MDeformWeight), "deformWeight");
					dvert->totweight = 1;
				}
				/* roots should be 1.0, the rest can be anything from 0.0 to 1.0 */
				dvert->dw->weight = key->weight;
				dvert++;
			}
		}
	}

	if(psys->hair_out_dm)
		psys->hair_out_dm->release(psys->hair_out_dm);

	psys->clmd->point_cache = psys->pointcache;
	psys->clmd->sim_parms->effector_weights = psys->part->effector_weights;

	psys->hair_out_dm = clothModifier_do(psys->clmd, sim->scene, sim->ob, dm);

	psys->clmd->sim_parms->effector_weights = NULL;
}
static void hair_step(ParticleSimulationData *sim, float cfra)
{
	ParticleSystem *psys = sim->psys;
/*	ParticleSettings *part = psys->part; */
	PARTICLE_P;
	float disp = (float)psys_get_current_display_percentage(psys)/100.0f;

	BLI_srandom(psys->seed);

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

	if(psys->recalc & PSYS_RECALC_RESET) {
		/* need this for changing subsurf levels */
		psys_calc_dmcache(sim->ob, sim->psmd->dm, 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);

	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= psys_get_lattice(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) {
			VECSUB(key->co, key->co, root->co);
			psys_vec_rot_to_face(sim->psmd->dm, 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)
			root->co[0] = root->co[1] = root->co[2] = 0.0f;
	}
}
/************************************************/
/*			System Core							*/
/************************************************/
/* 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.0){
		LOOP_EXISTING_PARTICLES {
			psys_get_texture(sim, pa, &ptex, PAMAP_SIZE, cfra);
			pa->size = part->size*ptex.size;
			if(part->randsize > 0.0)
				pa->size *= 1.0f - part->randsize * PSYS_FRAND(p + 1);

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

	BLI_srandom(31415926 + (int)cfra + psys->seed);

	psys_update_effectors(sim);

	if(part->type != PART_HAIR)
		sim->colliders = get_collider_cache(sim->scene, sim->ob, NULL);

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

			psys_update_particle_tree(psys, cfra);

			boids_precalc_rules(part, cfra);

			for(; pt; pt=pt->next) {
				if(pt->ob)
					psys_update_particle_tree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
			}
			break;
		}
		case PART_PHYS_FLUID:
		{
			ParticleTarget *pt = psys->targets.first;
			psys_update_particle_tree(psys, cfra);
			
			for(; pt; pt=pt->next) {  /* Updating others systems particle tree for fluid-fluid interaction */
				if(pt->ob)
					psys_update_particle_tree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra);
			}
			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.0)
			pa->size *= 1.0f - part->randsize * PSYS_FRAND(p + 1);

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

		/* 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 */
				apply_particle_forces(sim, p, pa->state.time, cfra);
	
				/* deflection */
				if(sim->colliders)
					deflect_particle(sim, p, pa->state.time, cfra);

				/* rotations */
				rotate_particle(part, pa, pa->state.time, timestep);
			}
			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)
						deflect_particle(sim, p, pa->state.time, cfra);
				}
			}
			break;
		}
		case PART_PHYS_FLUID:
		{
			EdgeHash *springhash = build_fluid_springhash(psys);
			float *gravity = NULL;

			if(psys_uses_gravity(sim))
				gravity = sim->scene->physics_settings.gravity;

			/* do global forces & effectors */
			LOOP_DYNAMIC_PARTICLES {
				apply_particle_forces(sim, p, pa->state.time, cfra);
				/* in fluids forces only effect velocity */
				copy_v3_v3(pa->state.co, pa->prev_state.co);
			}

			/* actual fluids calculations (not threadsafe!) */
			LOOP_DYNAMIC_PARTICLES {
				apply_particle_fluidsim(sim->ob, psys, pa, pa->state.time*timestep, gravity, springhash);
			}

			/* Apply springs to particles */
			apply_fluid_springs(psys, timestep);

			/* apply velocity, collisions and rotation */
			LOOP_DYNAMIC_PARTICLES {
				/* velocity holds forces and viscosity, so apply them before collisions */
				madd_v3_v3fl(pa->state.co, pa->state.vel, pa->state.time*timestep);

				/* calculate new velocity based on new-old location */
				sub_v3_v3v3(pa->state.vel, pa->state.co, pa->prev_state.co);
				mul_v3_fl(pa->state.vel, 1.f/(pa->state.time*timestep));

				if(sim->colliders)
					deflect_particle(sim, p, pa->state.time, cfra);
				
				/* SPH particles are not physical particles, just interpolation particles,  thus rotation has not a direct sense for them */	
				rotate_particle(part, pa, pa->state.time, timestep);  
			}

			if(springhash) {
				BLI_edgehash_free(springhash, NULL);
				springhash = NULL;
			}
			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;
	}

	free_collider_cache(&sim->colliders);
}
static void update_children(ParticleSimulationData *sim)
{
	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 != get_psys_tot_child(sim->scene, sim->psys))
			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)
{
	ParticleSystem *psys = sim->psys;
	ParticleSettings *part = psys->part;
	ParticleTexture ptex;
	PARTICLE_P;
	float disp, dietime;

	BLI_srandom(psys->seed);

	psys_update_effectors(sim);
	
	disp= (float)psys_get_current_display_percentage(psys)/100.0f;

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

		psys->lattice= psys_get_lattice(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){
			end_latt_deform(psys->lattice);
			psys->lattice= NULL;
		}

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

static void particles_fluid_step(ParticleSimulationData *sim, int UNUSED(cfra))
{	
	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 */
	#ifndef DISABLE_ELBEEM
	{
		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_snprintf(filename, sizeof(filename), "%sfluidsurface_particles_####.gz", fss->surfdataPath);
			
			BLI_path_abs(filename, G.main->name);
			BLI_path_frame(filename, curFrame, 0); // fixed #frame-no 

			gzf = gzopen(filename, "rb");
			if (!gzf) {
				snprintf(debugStrBuffer,256,"readFsPartData::error - Unable to open file for reading '%s' \n", filename); 
				// XXX bad level call elbeemDebugOut(debugStrBuffer);
				return;
			}
	
			gzread(gzf, &totpart, sizeof(totpart));
			totpart = (G.rendering)?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;
					}
	
					pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0f;
					pa->state.rot[0] = 1.0;
					pa->state.rot[1] = pa->state.rot[2] = pa->state.rot[3] = 0.0;
	
					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;
			snprintf(debugStrBuffer,256,"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
	}
	#endif // DISABLE_ELBEEM
}

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)
{
	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 && BKE_ptcache_get_continue_physics() == 0) {
		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 everythin on start frame */
		if((int)cfra == startframe) {
			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 & PSYS_RECALC_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);

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

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

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

		/* 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= (float)psys_get_current_display_percentage(psys)/100.0f;

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

	if(psys->totpart) {
		int dframe, subframe = 0, totframesback = 0, totsubframe = part->subframes+1;
		float fraction;
		
		/* 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);
		
		for(dframe=-totframesback; dframe<=0; dframe++) {
			/* ok now we're all set so let's go */
			for (subframe = 1; subframe <= totsubframe; subframe++) {
				fraction = (float)subframe/(float)totsubframe;
				dynamics_step(sim, cfra+dframe+fraction - 1.f);
				psys->cfra = cfra+dframe+fraction - 1.f;
			}
		}
		
	}
	
/* 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);

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

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

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

	if(part->from == PART_FROM_PARTICLE) {
		//if(part->type != PART_REACTOR)
		part->from = PART_FROM_FACE;
		if(part->distr == PART_DISTR_GRID && part->from != PART_FROM_VERT)
			part->distr = PART_DISTR_JIT;
	}

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

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

		if(ELEM3(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(sim->ob, sim->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(sim->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;
		}
}

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

	fluid->radius = 0.5f;
	fluid->spring_k = 0.f;
	fluid->plasticity_constant = 0.1f;
	fluid->yield_ratio = 0.1f;
	fluid->rest_length = 0.5f;
	fluid->viscosity_omega = 2.f;
	fluid->viscosity_beta = 0.f;
	fluid->stiffness_k = 0.1f;
	fluid->stiffness_knear = 0.05f;
	fluid->rest_density = 10.f;
	fluid->buoyancy = 0.f;
}

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

	if(part->phystype == PART_PHYS_BOIDS && part->boids == NULL) {
		BoidState *state;

		part->boids = MEM_callocN(sizeof(BoidSettings), "Boid Settings");
		boid_default_settings(part->boids);

		state = boid_new_state(part->boids);
		BLI_addtail(&state->rules, boid_new_rule(eBoidRuleType_Separate));
		BLI_addtail(&state->rules, boid_new_rule(eBoidRuleType_Flock));

		((BoidRule*)state->rules.first)->flag |= BOIDRULE_CURRENT;

		state->flag |= BOIDSTATE_CURRENT;
		BLI_addtail(&part->boids->states, state);
	}
	else if(part->phystype == PART_PHYS_FLUID && part->fluid == NULL) {
		part->fluid = MEM_callocN(sizeof(SPHFluidSettings), "SPH Fluid Settings");
		fluid_default_settings(part);
	}

	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 & PSYS_RECALC_RESET || (psys->part->flag & PART_HAIR_REGROW && !psys->edit))) {
		return 1;
	}

	return 0;
}

/* 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(Scene *scene, Object *ob, ParticleSystem *psys)
{
	ParticleSimulationData sim= {0};
	ParticleSettings *part = psys->part;
	float cfra;

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

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

	cfra= BKE_curframe(scene);

	sim.scene= scene;
	sim.ob= ob;
	sim.psys= psys;
	sim.psmd= psys_get_modifier(ob, psys);

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

	/* execute drivers only, as animation has already been done */
	BKE_animsys_evaluate_animdata(&part->id, part->adt, cfra, ADT_RECALC_DRIVERS);

	if(psys->recalc & PSYS_RECALC_TYPE)
		psys_changed_type(&sim);

	if(psys->recalc & PSYS_RECALC_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);
			}
			/* (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->edit && psys->free_edit) {
					psys->free_edit(psys->edit);
					psys->edit = NULL;
					psys->free_edit = NULL;
				}

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

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

				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);
			break;
		}
		case PART_FLUID:
		{
			particles_fluid_step(&sim, (int)cfra);
			break;
		}
		default:
		{
			switch(part->phystype) {
				case PART_PHYS_NO:
				case PART_PHYS_KEYED:
				{
					PARTICLE_P;
					float disp = (float)psys_get_current_display_percentage(psys)/100.0f;

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

					if(emit_particles(&sim, NULL, cfra) || (psys->recalc & PSYS_RECALC_RESET)) {
						free_keyed_keys(psys);
						distribute_particles(&sim, part->from);
						initialize_all_particles(&sim);
					}

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

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

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

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

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

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