Bug fix: particles emitted from a moving emitter exploded with verlet integration
[blender.git] / source / blender / blenkernel / intern / particle_system.c
1 /* particle_system.c
2  *
3  *
4  * $Id$
5  *
6  * ***** BEGIN GPL LICENSE BLOCK *****
7  *
8  * This program is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * as published by the Free Software Foundation; either version 2
11  * of the License, or (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software Foundation,
20  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
21  *
22  * The Original Code is Copyright (C) 2007 by Janne Karhu.
23  * All rights reserved.
24  *
25  * The Original Code is: all of this file.
26  *
27  * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn.
28  *
29  * ***** END GPL LICENSE BLOCK *****
30  */
31
32 /** \file blender/blenkernel/intern/particle_system.c
33  *  \ingroup bke
34  */
35
36
37 #include <stddef.h>
38 #include "BLI_storage.h" /* _LARGEFILE_SOURCE */
39
40 #include <stdlib.h>
41 #include <math.h>
42 #include <string.h>
43
44 #include "MEM_guardedalloc.h"
45
46 #include "DNA_anim_types.h"
47 #include "DNA_boid_types.h"
48 #include "DNA_particle_types.h"
49 #include "DNA_mesh_types.h"
50 #include "DNA_meshdata_types.h"
51 #include "DNA_modifier_types.h"
52 #include "DNA_object_force.h"
53 #include "DNA_object_types.h"
54 #include "DNA_material_types.h"
55 #include "DNA_curve_types.h"
56 #include "DNA_group_types.h"
57 #include "DNA_scene_types.h"
58 #include "DNA_texture_types.h"
59 #include "DNA_ipo_types.h" // XXX old animation system stuff... to be removed!
60 #include "DNA_listBase.h"
61
62 #include "BLI_edgehash.h"
63 #include "BLI_rand.h"
64 #include "BLI_jitter.h"
65 #include "BLI_math.h"
66 #include "BLI_blenlib.h"
67 #include "BLI_kdtree.h"
68 #include "BLI_kdopbvh.h"
69 #include "BLI_listbase.h"
70 #include "BLI_threads.h"
71 #include "BLI_storage.h" /* For _LARGEFILE64_SOURCE;  zlib needs this on some systems */
72 #include "BLI_utildefines.h"
73
74 #include "BKE_main.h"
75 #include "BKE_animsys.h"
76 #include "BKE_boids.h"
77 #include "BKE_cdderivedmesh.h"
78 #include "BKE_collision.h"
79 #include "BKE_displist.h"
80 #include "BKE_effect.h"
81 #include "BKE_particle.h"
82 #include "BKE_global.h"
83
84 #include "BKE_DerivedMesh.h"
85 #include "BKE_object.h"
86 #include "BKE_material.h"
87 #include "BKE_cloth.h"
88 #include "BKE_depsgraph.h"
89 #include "BKE_lattice.h"
90 #include "BKE_pointcache.h"
91 #include "BKE_mesh.h"
92 #include "BKE_modifier.h"
93 #include "BKE_scene.h"
94 #include "BKE_bvhutils.h"
95
96 #include "PIL_time.h"
97
98 #include "RE_shader_ext.h"
99
100 /* fluid sim particle import */
101 #ifndef DISABLE_ELBEEM
102 #include "DNA_object_fluidsim.h"
103 #include "LBM_fluidsim.h"
104 #include <zlib.h>
105 #include <string.h>
106
107 #ifdef WIN32
108 #ifndef snprintf
109 #define snprintf _snprintf
110 #endif
111 #endif
112
113 #endif // DISABLE_ELBEEM
114
115 /************************************************/
116 /*                      Reacting to system events                       */
117 /************************************************/
118
119 static int particles_are_dynamic(ParticleSystem *psys) {
120         if(psys->pointcache->flag & PTCACHE_BAKED)
121                 return 0;
122
123         if(psys->part->type == PART_HAIR)
124                 return psys->flag & PSYS_HAIR_DYNAMICS;
125         else
126                 return ELEM3(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
127 }
128
129 static int psys_get_current_display_percentage(ParticleSystem *psys)
130 {
131         ParticleSettings *part=psys->part;
132
133         if((psys->renderdata && !particles_are_dynamic(psys)) /* non-dynamic particles can be rendered fully */
134                 || (part->child_nbr && part->childtype) /* display percentage applies to children */
135                 || (psys->pointcache->flag & PTCACHE_BAKING)) /* baking is always done with full amount */
136                 return 100;
137
138         return psys->part->disp;
139 }
140
141 static int tot_particles(ParticleSystem *psys, PTCacheID *pid)
142 {
143         if(pid && psys->pointcache->flag & PTCACHE_EXTERNAL)
144                 return pid->cache->totpoint;
145         else if(psys->part->distr == PART_DISTR_GRID && psys->part->from != PART_FROM_VERT)
146                 return psys->part->grid_res * psys->part->grid_res * psys->part->grid_res - psys->totunexist;
147         else
148                 return psys->part->totpart - psys->totunexist;
149 }
150
151 void psys_reset(ParticleSystem *psys, int mode)
152 {
153         PARTICLE_P;
154
155         if(ELEM(mode, PSYS_RESET_ALL, PSYS_RESET_DEPSGRAPH)) {
156                 if(mode == PSYS_RESET_ALL || !(psys->flag & PSYS_EDITED)) {
157                         /* don't free if not absolutely necessary */
158                         if(psys->totpart != tot_particles(psys, NULL)) {
159                                 psys_free_particles(psys);
160                                 psys->totpart= 0;
161                         }
162
163                         psys->totkeyed= 0;
164                         psys->flag &= ~(PSYS_HAIR_DONE|PSYS_KEYED);
165
166                         if(psys->edit && psys->free_edit) {
167                                 psys->free_edit(psys->edit);
168                                 psys->edit = NULL;
169                                 psys->free_edit = NULL;
170                         }
171                 }
172         }
173         else if(mode == PSYS_RESET_CACHE_MISS) {
174                 /* set all particles to be skipped */
175                 LOOP_PARTICLES
176                         pa->flag |= PARS_NO_DISP;
177         }
178
179         /* reset children */
180         if(psys->child) {
181                 MEM_freeN(psys->child);
182                 psys->child= NULL;
183         }
184
185         psys->totchild= 0;
186
187         /* reset path cache */
188         psys_free_path_cache(psys, psys->edit);
189
190         /* reset point cache */
191         BKE_ptcache_invalidate(psys->pointcache);
192
193         if(psys->fluid_springs) {
194                 MEM_freeN(psys->fluid_springs);
195                 psys->fluid_springs = NULL;
196         }
197
198         psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
199 }
200
201 static void realloc_particles(ParticleSimulationData *sim, int new_totpart)
202 {
203         ParticleSystem *psys = sim->psys;
204         ParticleSettings *part = psys->part;
205         ParticleData *newpars = NULL;
206         BoidParticle *newboids = NULL;
207         PARTICLE_P;
208         int totpart, totsaved = 0;
209
210         if(new_totpart<0) {
211                 if(part->distr==PART_DISTR_GRID  && part->from != PART_FROM_VERT) {
212                         totpart= part->grid_res;
213                         totpart*=totpart*totpart;
214                 }
215                 else
216                         totpart=part->totpart;
217         }
218         else
219                 totpart=new_totpart;
220
221         if(totpart != psys->totpart) {
222                 if(psys->edit && psys->free_edit) {
223                         psys->free_edit(psys->edit);
224                         psys->edit = NULL;
225                         psys->free_edit = NULL;
226                 }
227
228                 if(totpart) {
229                         newpars= MEM_callocN(totpart*sizeof(ParticleData), "particles");
230                         if(newpars == NULL)
231                                 return;
232
233                         if(psys->part->phystype == PART_PHYS_BOIDS) {
234                                 newboids= MEM_callocN(totpart*sizeof(BoidParticle), "boid particles");
235
236                                 if(newboids == NULL) {
237                                          /* allocation error! */
238                                         if(newpars)
239                                                 MEM_freeN(newpars);
240                                         return;
241                                 }
242                         }
243                 }
244         
245                 if(psys->particles) {
246                         totsaved=MIN2(psys->totpart,totpart);
247                         /*save old pars*/
248                         if(totsaved) {
249                                 memcpy(newpars,psys->particles,totsaved*sizeof(ParticleData));
250
251                                 if(psys->particles->boid)
252                                         memcpy(newboids, psys->particles->boid, totsaved*sizeof(BoidParticle));
253                         }
254
255                         if(psys->particles->keys)
256                                 MEM_freeN(psys->particles->keys);
257
258                         if(psys->particles->boid)
259                                 MEM_freeN(psys->particles->boid);
260
261                         for(p=0, pa=newpars; p<totsaved; p++, pa++) {
262                                 if(pa->keys) {
263                                         pa->keys= NULL;
264                                         pa->totkey= 0;
265                                 }
266                         }
267
268                         for(p=totsaved, pa=psys->particles+totsaved; p<psys->totpart; p++, pa++)
269                                 if(pa->hair) MEM_freeN(pa->hair);
270
271                         MEM_freeN(psys->particles);
272                         psys_free_pdd(psys);
273                 }
274                 
275                 psys->particles=newpars;
276                 psys->totpart=totpart;
277
278                 if(newboids) {
279                         LOOP_PARTICLES
280                                 pa->boid = newboids++;
281                 }
282         }
283
284         if(psys->child) {
285                 MEM_freeN(psys->child);
286                 psys->child=NULL;
287                 psys->totchild=0;
288         }
289 }
290
291 static int get_psys_child_number(struct Scene *scene, ParticleSystem *psys)
292 {
293         int nbr;
294
295         if(!psys->part->childtype)
296                 return 0;
297
298         if(psys->renderdata)
299                 nbr= psys->part->ren_child_nbr;
300         else
301                 nbr= psys->part->child_nbr;
302
303         return get_render_child_particle_number(&scene->r, nbr);
304 }
305
306 static int get_psys_tot_child(struct Scene *scene, ParticleSystem *psys)
307 {
308         return psys->totpart*get_psys_child_number(scene, psys);
309 }
310
311 static void alloc_child_particles(ParticleSystem *psys, int tot)
312 {
313         if(psys->child){
314                 /* only re-allocate if we have to */
315                 if(psys->part->childtype && psys->totchild == tot) {
316                         memset(psys->child, 0, tot*sizeof(ChildParticle));
317                         return;
318                 }
319
320                 MEM_freeN(psys->child);
321                 psys->child=NULL;
322                 psys->totchild=0;
323         }
324
325         if(psys->part->childtype) {
326                 psys->totchild= tot;
327                 if(psys->totchild)
328                         psys->child= MEM_callocN(psys->totchild*sizeof(ChildParticle), "child_particles");
329         }
330 }
331
332 /************************************************/
333 /*                      Distribution                                            */
334 /************************************************/
335
336 void psys_calc_dmcache(Object *ob, DerivedMesh *dm, ParticleSystem *psys)
337 {
338         /* use for building derived mesh mapping info:
339
340            node: the allocated links - total derived mesh element count 
341            nodearray: the array of nodes aligned with the base mesh's elements, so
342                                   each original elements can reference its derived elements
343         */
344         Mesh *me= (Mesh*)ob->data;
345         PARTICLE_P;
346         
347         /* CACHE LOCATIONS */
348         if(!dm->deformedOnly) {
349                 /* Will use later to speed up subsurf/derivedmesh */
350                 LinkNode *node, *nodedmelem, **nodearray;
351                 int totdmelem, totelem, i, *origindex;
352
353                 if(psys->part->from == PART_FROM_VERT) {
354                         totdmelem= dm->getNumVerts(dm);
355                         totelem= me->totvert;
356                         origindex= dm->getVertDataArray(dm, CD_ORIGINDEX);
357                 }
358                 else { /* FROM_FACE/FROM_VOLUME */
359                         totdmelem= dm->getNumFaces(dm);
360                         totelem= me->totface;
361                         origindex= dm->getFaceDataArray(dm, CD_ORIGINDEX);
362                 }
363         
364                 nodedmelem= MEM_callocN(sizeof(LinkNode)*totdmelem, "psys node elems");
365                 nodearray= MEM_callocN(sizeof(LinkNode *)*totelem, "psys node array");
366                 
367                 for(i=0, node=nodedmelem; i<totdmelem; i++, origindex++, node++) {
368                         node->link= SET_INT_IN_POINTER(i);
369
370                         if(*origindex != -1) {
371                                 if(nodearray[*origindex]) {
372                                         /* prepend */
373                                         node->next = nodearray[*origindex];
374                                         nodearray[*origindex]= node;
375                                 }
376                                 else
377                                         nodearray[*origindex]= node;
378                         }
379                 }
380                 
381                 /* cache the verts/faces! */
382                 LOOP_PARTICLES {
383                         if(pa->num < 0) {
384                                 pa->num_dmcache = -1;
385                                 continue;
386                         }
387
388                         if(psys->part->from == PART_FROM_VERT) {
389                                 if(nodearray[pa->num])
390                                         pa->num_dmcache= GET_INT_FROM_POINTER(nodearray[pa->num]->link);
391                         }
392                         else { /* FROM_FACE/FROM_VOLUME */
393                                 /* Note that sometimes the pa->num is over the nodearray size, this is bad, maybe there is a better place to fix this,
394                                  * but for now passing NULL is OK. every face will be searched for the particle so its slower - Campbell */
395                                 pa->num_dmcache= psys_particle_dm_face_lookup(ob, dm, pa->num, pa->fuv, pa->num < totelem ? nodearray[pa->num] : NULL);
396                         }
397                 }
398
399                 MEM_freeN(nodearray);
400                 MEM_freeN(nodedmelem);
401         }
402         else {
403                 /* TODO PARTICLE, make the following line unnecessary, each function
404                  * should know to use the num or num_dmcache, set the num_dmcache to
405                  * an invalid value, just incase */
406                 
407                 LOOP_PARTICLES
408                         pa->num_dmcache = -1;
409         }
410 }
411
412 static void distribute_simple_children(Scene *scene, Object *ob, DerivedMesh *finaldm, ParticleSystem *psys)
413 {
414         ChildParticle *cpa = NULL;
415         int i, p;
416         int child_nbr= get_psys_child_number(scene, psys);
417         int totpart= get_psys_tot_child(scene, psys);
418
419         alloc_child_particles(psys, totpart);
420
421         cpa = psys->child;
422         for(i=0; i<child_nbr; i++){
423                 for(p=0; p<psys->totpart; p++,cpa++){
424                         float length=2.0;
425                         cpa->parent=p;
426                                         
427                         /* create even spherical distribution inside unit sphere */
428                         while(length>=1.0f){
429                                 cpa->fuv[0]=2.0f*BLI_frand()-1.0f;
430                                 cpa->fuv[1]=2.0f*BLI_frand()-1.0f;
431                                 cpa->fuv[2]=2.0f*BLI_frand()-1.0f;
432                                 length=len_v3(cpa->fuv);
433                         }
434
435                         cpa->num=-1;
436                 }
437         }
438         /* dmcache must be updated for parent particles if children from faces is used */
439         psys_calc_dmcache(ob, finaldm, psys);
440 }
441 static void distribute_grid(DerivedMesh *dm, ParticleSystem *psys)
442 {
443         ParticleData *pa=NULL;
444         float min[3], max[3], delta[3], d;
445         MVert *mv, *mvert = dm->getVertDataArray(dm,0);
446         int totvert=dm->getNumVerts(dm), from=psys->part->from;
447         int i, j, k, p, res=psys->part->grid_res, size[3], axis;
448
449         mv=mvert;
450
451         /* find bounding box of dm */
452         copy_v3_v3(min, mv->co);
453         copy_v3_v3(max, mv->co);
454         mv++;
455
456         for(i=1; i<totvert; i++, mv++){
457                 min[0]=MIN2(min[0],mv->co[0]);
458                 min[1]=MIN2(min[1],mv->co[1]);
459                 min[2]=MIN2(min[2],mv->co[2]);
460
461                 max[0]=MAX2(max[0],mv->co[0]);
462                 max[1]=MAX2(max[1],mv->co[1]);
463                 max[2]=MAX2(max[2],mv->co[2]);
464         }
465
466         VECSUB(delta,max,min);
467
468         /* determine major axis */
469         axis = (delta[0]>=delta[1]) ? 0 : ((delta[1]>=delta[2]) ? 1 : 2);
470          
471         d = delta[axis]/(float)res;
472
473         size[axis] = res;
474         size[(axis+1)%3] = (int)ceil(delta[(axis+1)%3]/d);
475         size[(axis+2)%3] = (int)ceil(delta[(axis+2)%3]/d);
476
477         /* float errors grrr.. */
478         size[(axis+1)%3] = MIN2(size[(axis+1)%3],res);
479         size[(axis+2)%3] = MIN2(size[(axis+2)%3],res);
480
481         size[0] = MAX2(size[0], 1);
482         size[1] = MAX2(size[1], 1);
483         size[2] = MAX2(size[2], 1);
484
485         /* no full offset for flat/thin objects */
486         min[0]+= d < delta[0] ? d/2.f : delta[0]/2.f;
487         min[1]+= d < delta[1] ? d/2.f : delta[1]/2.f;
488         min[2]+= d < delta[2] ? d/2.f : delta[2]/2.f;
489
490         for(i=0,p=0,pa=psys->particles; i<res; i++){
491                 for(j=0; j<res; j++){
492                         for(k=0; k<res; k++,p++,pa++){
493                                 pa->fuv[0] = min[0] + (float)i*d;
494                                 pa->fuv[1] = min[1] + (float)j*d;
495                                 pa->fuv[2] = min[2] + (float)k*d;
496                                 pa->flag |= PARS_UNEXIST;
497                                 pa->hair_index = 0; /* abused in volume calculation */
498                         }
499                 }
500         }
501
502         /* enable particles near verts/edges/faces/inside surface */
503         if(from==PART_FROM_VERT){
504                 float vec[3];
505
506                 pa=psys->particles;
507
508                 min[0] -= d/2.0f;
509                 min[1] -= d/2.0f;
510                 min[2] -= d/2.0f;
511
512                 for(i=0,mv=mvert; i<totvert; i++,mv++){
513                         sub_v3_v3v3(vec,mv->co,min);
514                         vec[0]/=delta[0];
515                         vec[1]/=delta[1];
516                         vec[2]/=delta[2];
517                         (pa     +((int)(vec[0]*(size[0]-1))*res
518                                 +(int)(vec[1]*(size[1]-1)))*res
519                                 +(int)(vec[2]*(size[2]-1)))->flag &= ~PARS_UNEXIST;
520                 }
521         }
522         else if(ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)){
523                 float co1[3], co2[3];
524
525                 MFace *mface= NULL, *mface_array;
526                 float v1[3], v2[3], v3[3], v4[4], lambda;
527                 int a, a1, a2, a0mul, a1mul, a2mul, totface;
528                 int amax= from==PART_FROM_FACE ? 3 : 1;
529
530                 totface=dm->getNumFaces(dm);
531                 mface_array= dm->getFaceDataArray(dm,CD_MFACE);
532                 
533                 for(a=0; a<amax; a++){
534                         if(a==0){ a0mul=res*res; a1mul=res; a2mul=1; }
535                         else if(a==1){ a0mul=res; a1mul=1; a2mul=res*res; }
536                         else{ a0mul=1; a1mul=res*res; a2mul=res; }
537
538                         for(a1=0; a1<size[(a+1)%3]; a1++){
539                                 for(a2=0; a2<size[(a+2)%3]; a2++){
540                                         mface= mface_array;
541
542                                         pa = psys->particles + a1*a1mul + a2*a2mul;
543                                         copy_v3_v3(co1, pa->fuv);
544                                         co1[a] -= d < delta[a] ? d/2.f : delta[a]/2.f;
545                                         copy_v3_v3(co2, co1);
546                                         co2[a] += delta[a] + 0.001f*d;
547                                         co1[a] -= 0.001f*d;
548                                         
549                                         /* lets intersect the faces */
550                                         for(i=0; i<totface; i++,mface++){
551                                                 copy_v3_v3(v1, mvert[mface->v1].co);
552                                                 copy_v3_v3(v2, mvert[mface->v2].co);
553                                                 copy_v3_v3(v3, mvert[mface->v3].co);
554
555                                                 if(isect_axial_line_tri_v3(a, co1, co2, v2, v3, v1, &lambda)){
556                                                         if(from==PART_FROM_FACE)
557                                                                 (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
558                                                         else /* store number of intersections */
559                                                                 (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
560                                                 }
561                                                 
562                                                 if(mface->v4){
563                                                         copy_v3_v3(v4, mvert[mface->v4].co);
564
565                                                         if(isect_axial_line_tri_v3(a, co1, co2, v4, v1, v3, &lambda)){
566                                                                 if(from==PART_FROM_FACE)
567                                                                         (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
568                                                                 else
569                                                                         (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
570                                                         }
571                                                 }
572                                         }
573
574                                         if(from==PART_FROM_VOLUME){
575                                                 int in=pa->hair_index%2;
576                                                 if(in) pa->hair_index++;
577                                                 for(i=0; i<size[0]; i++){
578                                                         if(in || (pa+i*a0mul)->hair_index%2)
579                                                                 (pa+i*a0mul)->flag &= ~PARS_UNEXIST;
580                                                         /* odd intersections == in->out / out->in */
581                                                         /* even intersections -> in stays same */
582                                                         in=(in + (pa+i*a0mul)->hair_index) % 2;
583                                                 }
584                                         }
585                                 }
586                         }
587                 }
588         }
589
590         if(psys->part->flag & PART_GRID_HEXAGONAL) {
591                 for(i=0,p=0,pa=psys->particles; i<res; i++){
592                         for(j=0; j<res; j++){
593                                 for(k=0; k<res; k++,p++,pa++){
594                                         if(j%2)
595                                                 pa->fuv[0] += d/2.f;
596
597                                         if(k%2) {
598                                                 pa->fuv[0] += d/2.f;
599                                                 pa->fuv[1] += d/2.f;
600                                         }
601                                 }
602                         }
603                 }
604         }
605
606         if(psys->part->flag & PART_GRID_INVERT){
607                 for(i=0; i<size[0]; i++){
608                         for(j=0; j<size[1]; j++){
609                                 pa=psys->particles + res*(i*res + j);
610                                 for(k=0; k<size[2]; k++, pa++){
611                                         pa->flag ^= PARS_UNEXIST;
612                                 }
613                         }
614                 }
615         }
616
617         if(psys->part->grid_rand > 0.f) {
618                 float rfac = d * psys->part->grid_rand;
619                 for(p=0,pa=psys->particles; p<psys->totpart; p++,pa++){
620                         if(pa->flag & PARS_UNEXIST)
621                                 continue;
622
623                         pa->fuv[0] += rfac * (PSYS_FRAND(p + 31) - 0.5f);
624                         pa->fuv[1] += rfac * (PSYS_FRAND(p + 32) - 0.5f);
625                         pa->fuv[2] += rfac * (PSYS_FRAND(p + 33) - 0.5f);
626                 }
627         }
628 }
629
630 /* modified copy from rayshade.c */
631 static void hammersley_create(float *out, int n, int seed, float amount)
632 {
633         RNG *rng;
634         double p, t, offs[2];
635         int k, kk;
636
637         rng = rng_new(31415926 + n + seed);
638         offs[0]= rng_getDouble(rng) + amount;
639         offs[1]= rng_getDouble(rng) + amount;
640         rng_free(rng);
641
642         for (k = 0; k < n; k++) {
643                 t = 0;
644                 for (p = 0.5, kk = k; kk; p *= 0.5, kk >>= 1)
645                         if (kk & 1) /* kk mod 2 = 1 */
646                                 t += p;
647     
648                 out[2*k + 0]= fmod((double)k/(double)n + offs[0], 1.0);
649                 out[2*k + 1]= fmod(t + offs[1], 1.0);
650         }
651 }
652
653 /* modified copy from effect.c */
654 static void init_mv_jit(float *jit, int num, int seed2, float amount)
655 {
656         RNG *rng;
657         float *jit2, x, rad1, rad2, rad3;
658         int i, num2;
659
660         if(num==0) return;
661
662         rad1= (float)(1.0/sqrt((float)num));
663         rad2= (float)(1.0/((float)num));
664         rad3= (float)sqrt((float)num)/((float)num);
665
666         rng = rng_new(31415926 + num + seed2);
667         x= 0;
668                 num2 = 2 * num;
669         for(i=0; i<num2; i+=2) {
670         
671                 jit[i]= x + amount*rad1*(0.5f - rng_getFloat(rng));
672                 jit[i+1]= i/(2.0f*num) + amount*rad1*(0.5f - rng_getFloat(rng));
673                 
674                 jit[i]-= (float)floor(jit[i]);
675                 jit[i+1]-= (float)floor(jit[i+1]);
676                 
677                 x+= rad3;
678                 x -= (float)floor(x);
679         }
680
681         jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
682
683         for (i=0 ; i<4 ; i++) {
684                 BLI_jitterate1(jit, jit2, num, rad1);
685                 BLI_jitterate1(jit, jit2, num, rad1);
686                 BLI_jitterate2(jit, jit2, num, rad2);
687         }
688         MEM_freeN(jit2);
689         rng_free(rng);
690 }
691
692 static void psys_uv_to_w(float u, float v, int quad, float *w)
693 {
694         float vert[4][3], co[3];
695
696         if(!quad) {
697                 if(u+v > 1.0f)
698                         v= 1.0f-v;
699                 else
700                         u= 1.0f-u;
701         }
702
703         vert[0][0]= 0.0f; vert[0][1]= 0.0f; vert[0][2]= 0.0f;
704         vert[1][0]= 1.0f; vert[1][1]= 0.0f; vert[1][2]= 0.0f;
705         vert[2][0]= 1.0f; vert[2][1]= 1.0f; vert[2][2]= 0.0f;
706
707         co[0]= u;
708         co[1]= v;
709         co[2]= 0.0f;
710
711         if(quad) {
712                 vert[3][0]= 0.0f; vert[3][1]= 1.0f; vert[3][2]= 0.0f;
713                 interp_weights_poly_v3( w,vert, 4, co);
714         }
715         else {
716                 interp_weights_poly_v3( w,vert, 3, co);
717                 w[3]= 0.0f;
718         }
719 }
720
721 /* Find the index in "sum" array before "value" is crossed. */
722 static int distribute_binary_search(float *sum, int n, float value)
723 {
724         int mid, low=0, high=n;
725
726         if(value == 0.f)
727                 return 0;
728
729         while(low <= high) {
730                 mid= (low + high)/2;
731                 
732                 if(sum[mid] < value && value <= sum[mid+1])
733                         return mid;
734                 
735                 if(sum[mid] >= value)
736                         high= mid - 1;
737                 else if(sum[mid] < value)
738                         low= mid + 1;
739                 else
740                         return mid;
741         }
742
743         return low;
744 }
745
746 /* the max number if calls to rng_* funcs within psys_thread_distribute_particle
747  * be sure to keep up to date if this changes */
748 #define PSYS_RND_DIST_SKIP 2
749
750 /* note: this function must be thread safe, for from == PART_FROM_CHILD */
751 #define ONLY_WORKING_WITH_PA_VERTS 0
752 static void distribute_threads_exec(ParticleThread *thread, ParticleData *pa, ChildParticle *cpa, int p)
753 {
754         ParticleThreadContext *ctx= thread->ctx;
755         Object *ob= ctx->sim.ob;
756         DerivedMesh *dm= ctx->dm;
757         float *v1, *v2, *v3, *v4, nor[3], orco1[3], co1[3], co2[3], nor1[3];
758         float cur_d, min_d, randu, randv;
759         int from= ctx->from;
760         int cfrom= ctx->cfrom;
761         int distr= ctx->distr;
762         int i, intersect, tot;
763         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
764
765         if(from == PART_FROM_VERT) {
766                 /* TODO_PARTICLE - use original index */
767                 pa->num= ctx->index[p];
768                 pa->fuv[0] = 1.0f;
769                 pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
770
771 #if ONLY_WORKING_WITH_PA_VERTS
772                 if(ctx->tree){
773                         KDTreeNearest ptn[3];
774                         int w, maxw;
775
776                         psys_particle_on_dm(ctx->dm,from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co1,0,0,0,orco1,0);
777                         transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
778                         maxw = BLI_kdtree_find_n_nearest(ctx->tree,3,orco1,NULL,ptn);
779
780                         for(w=0; w<maxw; w++){
781                                 pa->verts[w]=ptn->num;
782                         }
783                 }
784 #endif
785         }
786         else if(from == PART_FROM_FACE || from == PART_FROM_VOLUME) {
787                 MFace *mface;
788
789                 pa->num = i = ctx->index[p];
790                 mface = dm->getFaceData(dm,i,CD_MFACE);
791                 
792                 switch(distr){
793                 case PART_DISTR_JIT:
794                         if(ctx->jitlevel == 1) {
795                                 if(mface->v4)
796                                         psys_uv_to_w(0.5f, 0.5f, mface->v4, pa->fuv);
797                                 else
798                                         psys_uv_to_w(0.33333f, 0.33333f, mface->v4, pa->fuv);
799                         }
800                         else {
801                                 ctx->jitoff[i] = fmod(ctx->jitoff[i],(float)ctx->jitlevel);
802                                 psys_uv_to_w(ctx->jit[2*(int)ctx->jitoff[i]], ctx->jit[2*(int)ctx->jitoff[i]+1], mface->v4, pa->fuv);
803                                 ctx->jitoff[i]++;
804                         }
805                         break;
806                 case PART_DISTR_RAND:
807                         randu= rng_getFloat(thread->rng);
808                         randv= rng_getFloat(thread->rng);
809                         rng_skip_tot -= 2;
810
811                         psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
812                         break;
813                 }
814                 pa->foffset= 0.0f;
815                 
816                 /* experimental */
817                 if(from==PART_FROM_VOLUME){
818                         MVert *mvert=dm->getVertDataArray(dm,CD_MVERT);
819
820                         tot=dm->getNumFaces(dm);
821
822                         psys_interpolate_face(mvert,mface,0,0,pa->fuv,co1,nor,0,0,0,0);
823
824                         normalize_v3(nor);
825                         mul_v3_fl(nor,-100.0);
826
827                         VECADD(co2,co1,nor);
828
829                         min_d=2.0;
830                         intersect=0;
831
832                         for(i=0,mface=dm->getFaceDataArray(dm,CD_MFACE); i<tot; i++,mface++){
833                                 if(i==pa->num) continue;
834
835                                 v1=mvert[mface->v1].co;
836                                 v2=mvert[mface->v2].co;
837                                 v3=mvert[mface->v3].co;
838
839                                 if(isect_line_tri_v3(co1, co2, v2, v3, v1, &cur_d, 0)){
840                                         if(cur_d<min_d){
841                                                 min_d=cur_d;
842                                                 pa->foffset=cur_d*50.0f; /* to the middle of volume */
843                                                 intersect=1;
844                                         }
845                                 }
846                                 if(mface->v4){
847                                         v4=mvert[mface->v4].co;
848
849                                         if(isect_line_tri_v3(co1, co2, v4, v1, v3, &cur_d, 0)){
850                                                 if(cur_d<min_d){
851                                                         min_d=cur_d;
852                                                         pa->foffset=cur_d*50.0f; /* to the middle of volume */
853                                                         intersect=1;
854                                                 }
855                                         }
856                                 }
857                         }
858                         if(intersect==0)
859                                 pa->foffset=0.0;
860                         else switch(distr){
861                                 case PART_DISTR_JIT:
862                                         pa->foffset*= ctx->jit[p%(2*ctx->jitlevel)];
863                                         break;
864                                 case PART_DISTR_RAND:
865                                         pa->foffset*=BLI_frand();
866                                         break;
867                         }
868                 }
869         }
870         else if(from == PART_FROM_CHILD) {
871                 MFace *mf;
872
873                 if(ctx->index[p] < 0) {
874                         cpa->num=0;
875                         cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]=0.0f;
876                         cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
877                         return;
878                 }
879
880                 mf= dm->getFaceData(dm, ctx->index[p], CD_MFACE);
881
882                 randu= rng_getFloat(thread->rng);
883                 randv= rng_getFloat(thread->rng);
884                 rng_skip_tot -= 2;
885
886                 psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);
887
888                 cpa->num = ctx->index[p];
889
890                 if(ctx->tree){
891                         KDTreeNearest ptn[10];
892                         int w,maxw;//, do_seams;
893                         float maxd,mind,/*dd,*/totw=0.0;
894                         int parent[10];
895                         float pweight[10];
896
897                         psys_particle_on_dm(dm,cfrom,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co1,nor1,NULL,NULL,orco1,NULL);
898                         transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
899                         maxw = BLI_kdtree_find_n_nearest(ctx->tree,4,orco1,NULL,ptn);
900
901                         maxd=ptn[maxw-1].dist;
902                         mind=ptn[0].dist;
903                         
904                         /* the weights here could be done better */
905                         for(w=0; w<maxw; w++){
906                                 parent[w]=ptn[w].index;
907                                 pweight[w]=(float)pow(2.0,(double)(-6.0f*ptn[w].dist/maxd));
908                         }
909                         for(;w<10; w++){
910                                 parent[w]=-1;
911                                 pweight[w]=0.0f;
912                         }
913
914                         for(w=0,i=0; w<maxw && i<4; w++){
915                                 if(parent[w]>=0){
916                                         cpa->pa[i]=parent[w];
917                                         cpa->w[i]=pweight[w];
918                                         totw+=pweight[w];
919                                         i++;
920                                 }
921                         }
922                         for(;i<4; i++){
923                                 cpa->pa[i]=-1;
924                                 cpa->w[i]=0.0f;
925                         }
926
927                         if(totw>0.0f) for(w=0; w<4; w++)
928                                 cpa->w[w]/=totw;
929
930                         cpa->parent=cpa->pa[0];
931                 }
932         }
933
934         if(rng_skip_tot > 0) /* should never be below zero */
935                 rng_skip(thread->rng, rng_skip_tot);
936 }
937
938 static void *distribute_threads_exec_cb(void *data)
939 {
940         ParticleThread *thread= (ParticleThread*)data;
941         ParticleSystem *psys= thread->ctx->sim.psys;
942         ParticleData *pa;
943         ChildParticle *cpa;
944         int p, totpart;
945
946         if(thread->ctx->from == PART_FROM_CHILD) {
947                 totpart= psys->totchild;
948                 cpa= psys->child;
949
950                 for(p=0; p<totpart; p++, cpa++) {
951                         if(thread->ctx->skip) /* simplification skip */
952                                 rng_skip(thread->rng, PSYS_RND_DIST_SKIP * thread->ctx->skip[p]);
953
954                         if((p+thread->num) % thread->tot == 0)
955                                 distribute_threads_exec(thread, NULL, cpa, p);
956                         else /* thread skip */
957                                 rng_skip(thread->rng, PSYS_RND_DIST_SKIP);
958                 }
959         }
960         else {
961                 totpart= psys->totpart;
962                 pa= psys->particles + thread->num;
963                 for(p=thread->num; p<totpart; p+=thread->tot, pa+=thread->tot)
964                         distribute_threads_exec(thread, pa, NULL, p);
965         }
966
967         return 0;
968 }
969
970 /* not thread safe, but qsort doesn't take userdata argument */
971 static int *COMPARE_ORIG_INDEX = NULL;
972 static int distribute_compare_orig_index(const void *p1, const void *p2)
973 {
974         int index1 = COMPARE_ORIG_INDEX[*(const int*)p1];
975         int index2 = COMPARE_ORIG_INDEX[*(const int*)p2];
976
977         if(index1 < index2)
978                 return -1;
979         else if(index1 == index2) {
980                 /* this pointer comparison appears to make qsort stable for glibc,
981                  * and apparently on solaris too, makes the renders reproducable */
982                 if(p1 < p2)
983                         return -1;
984                 else if(p1 == p2)
985                         return 0;
986                 else
987                         return 1;
988         }
989         else
990                 return 1;
991 }
992
993 static void distribute_invalid(Scene *scene, ParticleSystem *psys, int from)
994 {
995         if(from == PART_FROM_CHILD) {
996                 ChildParticle *cpa;
997                 int p, totchild = get_psys_tot_child(scene, psys);
998
999                 if(psys->child && totchild) {
1000                         for(p=0,cpa=psys->child; p<totchild; p++,cpa++){
1001                                 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]= 0.0;
1002                                 cpa->foffset= 0.0f;
1003                                 cpa->parent=0;
1004                                 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
1005                                 cpa->num= -1;
1006                         }
1007                 }
1008         }
1009         else {
1010                 PARTICLE_P;
1011                 LOOP_PARTICLES {
1012                         pa->fuv[0]=pa->fuv[1]=pa->fuv[2]= pa->fuv[3]= 0.0;
1013                         pa->foffset= 0.0f;
1014                         pa->num= -1;
1015                 }
1016         }
1017 }
1018
1019 /* Creates a distribution of coordinates on a DerivedMesh       */
1020 /* This is to denote functionality that does not yet work with mesh - only derived mesh */
1021 static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, DerivedMesh *finaldm, int from)
1022 {
1023         ParticleThreadContext *ctx= threads[0].ctx;
1024         Object *ob= ctx->sim.ob;
1025         ParticleSystem *psys= ctx->sim.psys;
1026         ParticleData *pa=0, *tpars= 0;
1027         ParticleSettings *part;
1028         ParticleSeam *seams= 0;
1029         KDTree *tree=0;
1030         DerivedMesh *dm= NULL;
1031         float *jit= NULL;
1032         int i, seed, p=0, totthread= threads[0].tot;
1033         int cfrom=0;
1034         int totelem=0, totpart, *particle_element=0, children=0, totseam=0;
1035         int jitlevel= 1, distr;
1036         float *element_weight=NULL,*element_sum=NULL,*jitter_offset=NULL, *vweight=NULL;
1037         float cur, maxweight=0.0, tweight, totweight, inv_totweight, co[3], nor[3], orco[3], ornor[3];
1038         
1039         if(ELEM3(NULL, ob, psys, psys->part))
1040                 return 0;
1041
1042         part=psys->part;
1043         totpart=psys->totpart;
1044         if(totpart==0)
1045                 return 0;
1046
1047         if (!finaldm->deformedOnly && !finaldm->getFaceDataArray(finaldm, CD_ORIGINDEX)) {
1048                 printf("Can't create particles with the current modifier stack, disable destructive modifiers\n");
1049 // XXX          error("Can't paint with the current modifier stack, disable destructive modifiers");
1050                 return 0;
1051         }
1052
1053         /* First handle special cases */
1054         if(from == PART_FROM_CHILD) {
1055                 /* Simple children */
1056                 if(part->childtype != PART_CHILD_FACES) {
1057                         BLI_srandom(31415926 + psys->seed + psys->child_seed);
1058                         distribute_simple_children(scene, ob, finaldm, psys);
1059                         return 0;
1060                 }
1061         }
1062         else {
1063                 /* Grid distribution */
1064                 if(part->distr==PART_DISTR_GRID && from != PART_FROM_VERT){
1065                         BLI_srandom(31415926 + psys->seed);
1066                         dm= CDDM_from_mesh((Mesh*)ob->data, ob);
1067                         distribute_grid(dm,psys);
1068                         dm->release(dm);
1069                         return 0;
1070                 }
1071         }
1072         
1073         /* Create trees and original coordinates if needed */
1074         if(from == PART_FROM_CHILD) {
1075                 distr=PART_DISTR_RAND;
1076                 BLI_srandom(31415926 + psys->seed + psys->child_seed);
1077                 dm= finaldm;
1078                 children=1;
1079
1080                 tree=BLI_kdtree_new(totpart);
1081
1082                 for(p=0,pa=psys->particles; p<totpart; p++,pa++){
1083                         psys_particle_on_dm(dm,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co,nor,0,0,orco,ornor);
1084                         transform_mesh_orco_verts((Mesh*)ob->data, &orco, 1, 1);
1085                         BLI_kdtree_insert(tree, p, orco, ornor);
1086                 }
1087
1088                 BLI_kdtree_balance(tree);
1089
1090                 totpart = get_psys_tot_child(scene, psys);
1091                 cfrom = from = PART_FROM_FACE;
1092         }
1093         else {
1094                 distr = part->distr;
1095                 BLI_srandom(31415926 + psys->seed);
1096                 
1097                 dm= CDDM_from_mesh((Mesh*)ob->data, ob);
1098
1099                 /* we need orco for consistent distributions */
1100                 DM_add_vert_layer(dm, CD_ORCO, CD_ASSIGN, get_mesh_orco_verts(ob));
1101
1102                 if(from == PART_FROM_VERT) {
1103                         MVert *mv= dm->getVertDataArray(dm, CD_MVERT);
1104                         float (*orcodata)[3]= dm->getVertDataArray(dm, CD_ORCO);
1105                         int totvert = dm->getNumVerts(dm);
1106
1107                         tree=BLI_kdtree_new(totvert);
1108
1109                         for(p=0; p<totvert; p++) {
1110                                 if(orcodata) {
1111                                         VECCOPY(co,orcodata[p])
1112                                         transform_mesh_orco_verts((Mesh*)ob->data, &co, 1, 1);
1113                                 }
1114                                 else
1115                                         VECCOPY(co,mv[p].co)
1116                                 BLI_kdtree_insert(tree,p,co,NULL);
1117                         }
1118
1119                         BLI_kdtree_balance(tree);
1120                 }
1121         }
1122
1123         /* Get total number of emission elements and allocate needed arrays */
1124         totelem = (from == PART_FROM_VERT) ? dm->getNumVerts(dm) : dm->getNumFaces(dm);
1125
1126         if(totelem == 0){
1127                 distribute_invalid(scene, psys, children ? PART_FROM_CHILD : 0);
1128
1129                 if(G.f & G_DEBUG)
1130                         fprintf(stderr,"Particle distribution error: Nothing to emit from!\n");
1131
1132                 if(dm != finaldm) dm->release(dm);
1133                 return 0;
1134         }
1135
1136         element_weight  = MEM_callocN(sizeof(float)*totelem, "particle_distribution_weights");
1137         particle_element= MEM_callocN(sizeof(int)*totpart, "particle_distribution_indexes");
1138         element_sum             = MEM_callocN(sizeof(float)*(totelem+1), "particle_distribution_sum");
1139         jitter_offset   = MEM_callocN(sizeof(float)*totelem, "particle_distribution_jitoff");
1140
1141         /* Calculate weights from face areas */
1142         if((part->flag&PART_EDISTR || children) && from != PART_FROM_VERT){
1143                 MVert *v1, *v2, *v3, *v4;
1144                 float totarea=0.f, co1[3], co2[3], co3[3], co4[3];
1145                 float (*orcodata)[3];
1146                 
1147                 orcodata= dm->getVertDataArray(dm, CD_ORCO);
1148
1149                 for(i=0; i<totelem; i++){
1150                         MFace *mf=dm->getFaceData(dm,i,CD_MFACE);
1151
1152                         if(orcodata) {
1153                                 VECCOPY(co1, orcodata[mf->v1]);
1154                                 VECCOPY(co2, orcodata[mf->v2]);
1155                                 VECCOPY(co3, orcodata[mf->v3]);
1156                                 transform_mesh_orco_verts((Mesh*)ob->data, &co1, 1, 1);
1157                                 transform_mesh_orco_verts((Mesh*)ob->data, &co2, 1, 1);
1158                                 transform_mesh_orco_verts((Mesh*)ob->data, &co3, 1, 1);
1159                                 if(mf->v4) {
1160                                         VECCOPY(co4, orcodata[mf->v4]);
1161                                         transform_mesh_orco_verts((Mesh*)ob->data, &co4, 1, 1);
1162                                 }
1163                         }
1164                         else {
1165                                 v1= (MVert*)dm->getVertData(dm,mf->v1,CD_MVERT);
1166                                 v2= (MVert*)dm->getVertData(dm,mf->v2,CD_MVERT);
1167                                 v3= (MVert*)dm->getVertData(dm,mf->v3,CD_MVERT);
1168                                 VECCOPY(co1, v1->co);
1169                                 VECCOPY(co2, v2->co);
1170                                 VECCOPY(co3, v3->co);
1171                                 if(mf->v4) {
1172                                         v4= (MVert*)dm->getVertData(dm,mf->v4,CD_MVERT);
1173                                         VECCOPY(co4, v4->co);
1174                                 }
1175                         }
1176
1177                         cur = mf->v4 ? area_quad_v3(co1, co2, co3, co4) : area_tri_v3(co1, co2, co3);
1178                         
1179                         if(cur > maxweight)
1180                                 maxweight = cur;
1181
1182                         element_weight[i] = cur;
1183                         totarea += cur;
1184                 }
1185
1186                 for(i=0; i<totelem; i++)
1187                         element_weight[i] /= totarea;
1188
1189                 maxweight /= totarea;
1190         }
1191         else{
1192                 float min=1.0f/(float)(MIN2(totelem,totpart));
1193                 for(i=0; i<totelem; i++)
1194                         element_weight[i]=min;
1195                 maxweight=min;
1196         }
1197
1198         /* Calculate weights from vgroup */
1199         vweight = psys_cache_vgroup(dm,psys,PSYS_VG_DENSITY);
1200
1201         if(vweight){
1202                 if(from==PART_FROM_VERT) {
1203                         for(i=0;i<totelem; i++)
1204                                 element_weight[i]*=vweight[i];
1205                 }
1206                 else { /* PART_FROM_FACE / PART_FROM_VOLUME */
1207                         for(i=0;i<totelem; i++){
1208                                 MFace *mf=dm->getFaceData(dm,i,CD_MFACE);
1209                                 tweight = vweight[mf->v1] + vweight[mf->v2] + vweight[mf->v3];
1210                                 
1211                                 if(mf->v4) {
1212                                         tweight += vweight[mf->v4];
1213                                         tweight /= 4.0;
1214                                 }
1215                                 else {
1216                                         tweight /= 3.0;
1217                                 }
1218
1219                                 element_weight[i]*=tweight;
1220                         }
1221                 }
1222                 MEM_freeN(vweight);
1223         }
1224
1225         /* Calculate total weight of all elements */
1226         totweight= 0.0f;
1227         for(i=0;i<totelem; i++)
1228                 totweight += element_weight[i];
1229
1230         inv_totweight = (totweight > 0.f ? 1.f/totweight : 0.f);
1231
1232         /* Calculate cumulative weights */
1233         element_sum[0]= 0.0f;
1234         for(i=0; i<totelem; i++)
1235                 element_sum[i+1]= element_sum[i] + element_weight[i] * inv_totweight;
1236         
1237         /* Finally assign elements to particles */
1238         if((part->flag&PART_TRAND) || (part->simplify_flag&PART_SIMPLIFY_ENABLE)) {
1239                 float pos;
1240
1241                 for(p=0; p<totpart; p++) {
1242                         /* In theory element_sum[totelem] should be 1.0, but due to float errors this is not necessarily always true, so scale pos accordingly. */
1243                         pos= BLI_frand() * element_sum[totelem];
1244                         particle_element[p]= distribute_binary_search(element_sum, totelem, pos);
1245                         particle_element[p]= MIN2(totelem-1, particle_element[p]);
1246                         jitter_offset[particle_element[p]]= pos;
1247                 }
1248         }
1249         else {
1250                 double step, pos;
1251                 
1252                 step= (totpart < 2) ? 0.5 : 1.0/(double)totpart;
1253                 pos= 1e-16; /* tiny offset to avoid zero weight face */
1254                 i= 0;
1255
1256                 for(p=0; p<totpart; p++, pos+=step) {
1257                         while((i < totelem) && (pos > element_sum[i+1]))
1258                                 i++;
1259
1260                         particle_element[p]= MIN2(totelem-1, i);
1261
1262                         /* avoid zero weight face */
1263                         if(p == totpart-1 && element_weight[particle_element[p]] == 0.0f)
1264                                 particle_element[p]= particle_element[p-1];
1265
1266                         jitter_offset[particle_element[p]]= pos;
1267                 }
1268         }
1269
1270         MEM_freeN(element_sum);
1271
1272         /* For hair, sort by origindex (allows optimizations in rendering), */
1273         /* however with virtual parents the children need to be in random order. */
1274         if(part->type == PART_HAIR && !(part->childtype==PART_CHILD_FACES && part->parents!=0.0)) {
1275                 COMPARE_ORIG_INDEX = NULL;
1276
1277                 if(from == PART_FROM_VERT) {
1278                         if(dm->numVertData)
1279                                 COMPARE_ORIG_INDEX= dm->getVertDataArray(dm, CD_ORIGINDEX);
1280                 }
1281                 else {
1282                         if(dm->numFaceData)
1283                                 COMPARE_ORIG_INDEX= dm->getFaceDataArray(dm, CD_ORIGINDEX);
1284                 }
1285
1286                 if(COMPARE_ORIG_INDEX) {
1287                         qsort(particle_element, totpart, sizeof(int), distribute_compare_orig_index);
1288                         COMPARE_ORIG_INDEX = NULL;
1289                 }
1290         }
1291
1292         /* Create jittering if needed */
1293         if(distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
1294                 jitlevel= part->userjit;
1295                 
1296                 if(jitlevel == 0) {
1297                         jitlevel= totpart/totelem;
1298                         if(part->flag & PART_EDISTR) jitlevel*= 2;      /* looks better in general, not very scietific */
1299                         if(jitlevel<3) jitlevel= 3;
1300                 }
1301                 
1302                 jit= MEM_callocN((2+ jitlevel*2)*sizeof(float), "jit");
1303
1304                 /* for small amounts of particles we use regular jitter since it looks
1305                  * a bit better, for larger amounts we switch to hammersley sequence 
1306                  * because it is much faster */
1307                 if(jitlevel < 25)
1308                         init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
1309                 else
1310                         hammersley_create(jit, jitlevel+1, psys->seed, part->jitfac);
1311                 BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
1312         }
1313
1314         /* Setup things for threaded distribution */
1315         ctx->tree= tree;
1316         ctx->seams= seams;
1317         ctx->totseam= totseam;
1318         ctx->sim.psys= psys;
1319         ctx->index= particle_element;
1320         ctx->jit= jit;
1321         ctx->jitlevel= jitlevel;
1322         ctx->jitoff= jitter_offset;
1323         ctx->weight= element_weight;
1324         ctx->maxweight= maxweight;
1325         ctx->from= (children)? PART_FROM_CHILD: from;
1326         ctx->cfrom= cfrom;
1327         ctx->distr= distr;
1328         ctx->dm= dm;
1329         ctx->tpars= tpars;
1330
1331         if(children) {
1332                 totpart= psys_render_simplify_distribution(ctx, totpart);
1333                 alloc_child_particles(psys, totpart);
1334         }
1335
1336         if(!children || psys->totchild < 10000)
1337                 totthread= 1;
1338         
1339         seed= 31415926 + ctx->sim.psys->seed;
1340         for(i=0; i<totthread; i++) {
1341                 threads[i].rng= rng_new(seed);
1342                 threads[i].tot= totthread;
1343         }
1344
1345         return 1;
1346 }
1347
1348 static void distribute_particles_on_dm(ParticleSimulationData *sim, int from)
1349 {
1350         DerivedMesh *finaldm = sim->psmd->dm;
1351         ListBase threads;
1352         ParticleThread *pthreads;
1353         ParticleThreadContext *ctx;
1354         int i, totthread;
1355
1356         pthreads= psys_threads_create(sim);
1357
1358         if(!distribute_threads_init_data(pthreads, sim->scene, finaldm, from)) {
1359                 psys_threads_free(pthreads);
1360                 return;
1361         }
1362
1363         totthread= pthreads[0].tot;
1364         if(totthread > 1) {
1365                 BLI_init_threads(&threads, distribute_threads_exec_cb, totthread);
1366
1367                 for(i=0; i<totthread; i++)
1368                         BLI_insert_thread(&threads, &pthreads[i]);
1369
1370                 BLI_end_threads(&threads);
1371         }
1372         else
1373                 distribute_threads_exec_cb(&pthreads[0]);
1374
1375         psys_calc_dmcache(sim->ob, finaldm, sim->psys);
1376
1377         ctx= pthreads[0].ctx;
1378         if(ctx->dm != finaldm)
1379                 ctx->dm->release(ctx->dm);
1380
1381         psys_threads_free(pthreads);
1382 }
1383
1384 /* ready for future use, to emit particles without geometry */
1385 static void distribute_particles_on_shape(ParticleSimulationData *sim, int UNUSED(from))
1386 {
1387         distribute_invalid(sim->scene, sim->psys, 0);
1388
1389         fprintf(stderr,"Shape emission not yet possible!\n");
1390 }
1391
1392 static void distribute_particles(ParticleSimulationData *sim, int from)
1393 {
1394         PARTICLE_PSMD;
1395         int distr_error=0;
1396
1397         if(psmd){
1398                 if(psmd->dm)
1399                         distribute_particles_on_dm(sim, from);
1400                 else
1401                         distr_error=1;
1402         }
1403         else
1404                 distribute_particles_on_shape(sim, from);
1405
1406         if(distr_error){
1407                 distribute_invalid(sim->scene, sim->psys, from);
1408
1409                 fprintf(stderr,"Particle distribution error!\n");
1410         }
1411 }
1412
1413 /* threaded child particle distribution and path caching */
1414 ParticleThread *psys_threads_create(ParticleSimulationData *sim)
1415 {
1416         ParticleThread *threads;
1417         ParticleThreadContext *ctx;
1418         int i, totthread;
1419
1420         if(sim->scene->r.mode & R_FIXED_THREADS)
1421                 totthread= sim->scene->r.threads;
1422         else
1423                 totthread= BLI_system_thread_count();
1424         
1425         threads= MEM_callocN(sizeof(ParticleThread)*totthread, "ParticleThread");
1426         ctx= MEM_callocN(sizeof(ParticleThreadContext), "ParticleThreadContext");
1427
1428         ctx->sim = *sim;
1429         ctx->dm= ctx->sim.psmd->dm;
1430         ctx->ma= give_current_material(sim->ob, sim->psys->part->omat);
1431
1432         memset(threads, 0, sizeof(ParticleThread)*totthread);
1433
1434         for(i=0; i<totthread; i++) {
1435                 threads[i].ctx= ctx;
1436                 threads[i].num= i;
1437                 threads[i].tot= totthread;
1438         }
1439
1440         return threads;
1441 }
1442
1443 void psys_threads_free(ParticleThread *threads)
1444 {
1445         ParticleThreadContext *ctx= threads[0].ctx;
1446         int i, totthread= threads[0].tot;
1447
1448         /* path caching */
1449         if(ctx->vg_length)
1450                 MEM_freeN(ctx->vg_length);
1451         if(ctx->vg_clump)
1452                 MEM_freeN(ctx->vg_clump);
1453         if(ctx->vg_kink)
1454                 MEM_freeN(ctx->vg_kink);
1455         if(ctx->vg_rough1)
1456                 MEM_freeN(ctx->vg_rough1);
1457         if(ctx->vg_rough2)
1458                 MEM_freeN(ctx->vg_rough2);
1459         if(ctx->vg_roughe)
1460                 MEM_freeN(ctx->vg_roughe);
1461
1462         if(ctx->sim.psys->lattice){
1463                 end_latt_deform(ctx->sim.psys->lattice);
1464                 ctx->sim.psys->lattice= NULL;
1465         }
1466
1467         /* distribution */
1468         if(ctx->jit) MEM_freeN(ctx->jit);
1469         if(ctx->jitoff) MEM_freeN(ctx->jitoff);
1470         if(ctx->weight) MEM_freeN(ctx->weight);
1471         if(ctx->index) MEM_freeN(ctx->index);
1472         if(ctx->skip) MEM_freeN(ctx->skip);
1473         if(ctx->seams) MEM_freeN(ctx->seams);
1474         //if(ctx->vertpart) MEM_freeN(ctx->vertpart);
1475         BLI_kdtree_free(ctx->tree);
1476
1477         /* threads */
1478         for(i=0; i<totthread; i++) {
1479                 if(threads[i].rng)
1480                         rng_free(threads[i].rng);
1481                 if(threads[i].rng_path)
1482                         rng_free(threads[i].rng_path);
1483         }
1484
1485         MEM_freeN(ctx);
1486         MEM_freeN(threads);
1487 }
1488
1489 /* set particle parameters that don't change during particle's life */
1490 void initialize_particle(ParticleSimulationData *sim, ParticleData *pa, int p)
1491 {
1492         ParticleSystem *psys = sim->psys;
1493         ParticleSettings *part = psys->part;
1494         ParticleTexture ptex;
1495
1496         pa->flag &= ~PARS_UNEXIST;
1497
1498         if(part->type != PART_FLUID) {
1499                 psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);
1500                 
1501                 if(ptex.exist < PSYS_FRAND(p+125))
1502                         pa->flag |= PARS_UNEXIST;
1503
1504                 pa->time = (part->type == PART_HAIR) ? 0.f : part->sta + (part->end - part->sta)*ptex.time;
1505         }
1506
1507         pa->hair_index = 0;
1508         /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
1509         /* usage other than straight after distribute has to handle this index by itself - jahka*/
1510         //pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we dont have a derived mesh face */
1511 }
1512 static void initialize_all_particles(ParticleSimulationData *sim)
1513 {
1514         ParticleSystem *psys = sim->psys;
1515         PARTICLE_P;
1516
1517         psys->totunexist = 0;
1518
1519         LOOP_PARTICLES {
1520                 if((pa->flag & PARS_UNEXIST)==0)
1521                         initialize_particle(sim, pa, p);
1522
1523                 if(pa->flag & PARS_UNEXIST)
1524                         psys->totunexist++;
1525         }
1526
1527         /* Free unexisting particles. */
1528         if(psys->totpart && psys->totunexist == psys->totpart) {
1529                 if(psys->particles->boid)
1530                         MEM_freeN(psys->particles->boid);
1531
1532                 MEM_freeN(psys->particles);
1533                 psys->particles = NULL;
1534                 psys->totpart = psys->totunexist = 0;
1535         }
1536
1537         if(psys->totunexist) {
1538                 int newtotpart = psys->totpart - psys->totunexist;
1539                 ParticleData *npa, *newpars;
1540                 
1541                 npa = newpars = MEM_callocN(newtotpart * sizeof(ParticleData), "particles");
1542
1543                 for(p=0, pa=psys->particles; p<newtotpart; p++, pa++, npa++) {
1544                         while(pa->flag & PARS_UNEXIST)
1545                                 pa++;
1546
1547                         memcpy(npa, pa, sizeof(ParticleData));
1548                 }
1549
1550                 if(psys->particles->boid)
1551                         MEM_freeN(psys->particles->boid);
1552                 MEM_freeN(psys->particles);
1553                 psys->particles = newpars;
1554                 psys->totpart -= psys->totunexist;
1555
1556                 if(psys->particles->boid) {
1557                         BoidParticle *newboids = MEM_callocN(psys->totpart * sizeof(BoidParticle), "boid particles");
1558
1559                         LOOP_PARTICLES
1560                                 pa->boid = newboids++;
1561
1562                 }
1563         }
1564 }
1565 /* sets particle to the emitter surface with initial velocity & rotation */
1566 void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
1567 {
1568         Object *ob = sim->ob;
1569         ParticleSystem *psys = sim->psys;
1570         ParticleSettings *part;
1571         ParticleTexture ptex;
1572         float fac, phasefac, nor[3]={0,0,0},loc[3],vel[3]={0.0,0.0,0.0},rot[4],q2[4];
1573         float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3]={0.0,0.0,0.0};
1574         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};
1575         float q_phase[4];
1576         int p = pa - psys->particles;
1577         part=psys->part;
1578         
1579         /* get precise emitter matrix if particle is born */
1580         if(part->type!=PART_HAIR && dtime > 0.f && pa->time < cfra && pa->time >= sim->psys->cfra) {
1581                 /* we have to force RECALC_ANIM here since where_is_objec_time only does drivers */
1582                 while(ob) {
1583                         BKE_animsys_evaluate_animdata(&ob->id, ob->adt, pa->time, ADT_RECALC_ANIM);
1584                         ob = ob->parent;
1585                 }
1586                 ob = sim->ob;
1587                 where_is_object_time(sim->scene, ob, pa->time);
1588         }
1589
1590         /* get birth location from object               */
1591         if(part->tanfac != 0.f)
1592                 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0,0);
1593         else
1594                 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,0,0,0,0);
1595                 
1596         /* get possible textural influence */
1597         psys_get_texture(sim, pa, &ptex, PAMAP_IVEL|PAMAP_LIFE, cfra);
1598
1599         /* particles live in global space so    */
1600         /* let's convert:                                               */
1601         /* -location                                                    */
1602         mul_m4_v3(ob->obmat, loc);
1603                 
1604         /* -normal                                                              */
1605         mul_mat3_m4_v3(ob->obmat, nor);
1606         normalize_v3(nor);
1607
1608         /* -tangent                                                             */
1609         if(part->tanfac!=0.0){
1610                 //float phase=vg_rot?2.0f*(psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_rot)-0.5f):0.0f;
1611                 float phase=0.0f;
1612                 mul_v3_fl(vtan,-(float)cos(M_PI*(part->tanphase+phase)));
1613                 fac=-(float)sin(M_PI*(part->tanphase+phase));
1614                 VECADDFAC(vtan,vtan,utan,fac);
1615
1616                 mul_mat3_m4_v3(ob->obmat,vtan);
1617
1618                 VECCOPY(utan,nor);
1619                 mul_v3_fl(utan,dot_v3v3(vtan,nor));
1620                 VECSUB(vtan,vtan,utan);
1621                         
1622                 normalize_v3(vtan);
1623         }
1624                 
1625
1626         /* -velocity                                                    */
1627         if(part->randfac!=0.0){
1628                 r_vel[0] = 2.0f * (PSYS_FRAND(p + 10) - 0.5f);
1629                 r_vel[1] = 2.0f * (PSYS_FRAND(p + 11) - 0.5f);
1630                 r_vel[2] = 2.0f * (PSYS_FRAND(p + 12) - 0.5f);
1631
1632                 mul_mat3_m4_v3(ob->obmat, r_vel);
1633                 normalize_v3(r_vel);
1634         }
1635
1636         /* -angular velocity                                    */
1637         if(part->avemode==PART_AVE_RAND){
1638                 r_ave[0] = 2.0f * (PSYS_FRAND(p + 13) - 0.5f);
1639                 r_ave[1] = 2.0f * (PSYS_FRAND(p + 14) - 0.5f);
1640                 r_ave[2] = 2.0f * (PSYS_FRAND(p + 15) - 0.5f);
1641
1642                 mul_mat3_m4_v3(ob->obmat,r_ave);
1643                 normalize_v3(r_ave);
1644         }
1645                 
1646         /* -rotation                                                    */
1647         if(part->randrotfac != 0.0f){
1648                 r_rot[0] = 2.0f * (PSYS_FRAND(p + 16) - 0.5f);
1649                 r_rot[1] = 2.0f * (PSYS_FRAND(p + 17) - 0.5f);
1650                 r_rot[2] = 2.0f * (PSYS_FRAND(p + 18) - 0.5f);
1651                 r_rot[3] = 2.0f * (PSYS_FRAND(p + 19) - 0.5f);
1652                 normalize_qt(r_rot);
1653
1654                 mat4_to_quat(rot,ob->obmat);
1655                 mul_qt_qtqt(r_rot,r_rot,rot);
1656         }
1657 #if 0
1658         }
1659 #endif
1660
1661         if(part->phystype==PART_PHYS_BOIDS && pa->boid) {
1662                 BoidParticle *bpa = pa->boid;
1663                 float dvec[3], q[4], mat[3][3];
1664
1665                 copy_v3_v3(pa->state.co,loc);
1666
1667                 /* boids don't get any initial velocity  */
1668                 zero_v3(pa->state.vel);
1669
1670                 /* boids store direction in ave */
1671                 if(fabs(nor[2])==1.0f) {
1672                         sub_v3_v3v3(pa->state.ave, loc, ob->obmat[3]);
1673                         normalize_v3(pa->state.ave);
1674                 }
1675                 else {
1676                         VECCOPY(pa->state.ave, nor);
1677                 }
1678                 /* and gravity in r_ve */
1679                 bpa->gravity[0] = bpa->gravity[1] = 0.0f;
1680                 bpa->gravity[2] = -1.0f;
1681                 if((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY)
1682                         && sim->scene->physics_settings.gravity[2]!=0.0f)
1683                         bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
1684
1685                 /* calculate rotation matrix */
1686                 project_v3_v3v3(dvec, r_vel, pa->state.ave);
1687                 sub_v3_v3v3(mat[0], pa->state.ave, dvec);
1688                 normalize_v3(mat[0]);
1689                 negate_v3_v3(mat[2], r_vel);
1690                 normalize_v3(mat[2]);
1691                 cross_v3_v3v3(mat[1], mat[2], mat[0]);
1692                 
1693                 /* apply rotation */
1694                 mat3_to_quat_is_ok( q,mat);
1695                 copy_qt_qt(pa->state.rot, q);
1696
1697                 bpa->data.health = part->boids->health;
1698                 bpa->data.mode = eBoidMode_InAir;
1699                 bpa->data.state_id = ((BoidState*)part->boids->states.first)->id;
1700                 bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f;
1701         }
1702         else {
1703                 /* conversion done so now we apply new: */
1704                 /* -velocity from:                                              */
1705
1706                 /*              *reactions                                              */
1707                 if(dtime > 0.f){
1708                         sub_v3_v3v3(vel, pa->state.vel, pa->prev_state.vel);
1709                 }
1710
1711                 /*              *emitter velocity                               */
1712                 if(dtime != 0.f && part->obfac != 0.f){
1713                         sub_v3_v3v3(vel, loc, pa->state.co);
1714                         mul_v3_fl(vel, part->obfac/dtime);
1715                 }
1716                 
1717                 /*              *emitter normal                                 */
1718                 if(part->normfac != 0.f)
1719                         madd_v3_v3fl(vel, nor, part->normfac);
1720                 
1721                 /*              *emitter tangent                                */
1722                 if(sim->psmd && part->tanfac != 0.f)
1723                         madd_v3_v3fl(vel, vtan, part->tanfac);
1724
1725                 /*              *emitter object orientation             */
1726                 if(part->ob_vel[0] != 0.f) {
1727                         normalize_v3_v3(vec, ob->obmat[0]);
1728                         madd_v3_v3fl(vel, vec, part->ob_vel[0]);
1729                 }
1730                 if(part->ob_vel[1] != 0.f) {
1731                         normalize_v3_v3(vec, ob->obmat[1]);
1732                         madd_v3_v3fl(vel, vec, part->ob_vel[1]);
1733                 }
1734                 if(part->ob_vel[2] != 0.f) {
1735                         normalize_v3_v3(vec, ob->obmat[2]);
1736                         madd_v3_v3fl(vel, vec, part->ob_vel[2]);
1737                 }
1738
1739                 /*              *texture                                                */
1740                 /* TODO */
1741
1742                 /*              *random                                                 */
1743                 if(part->randfac != 0.f)
1744                         madd_v3_v3fl(vel, r_vel, part->randfac);
1745
1746                 /*              *particle                                               */
1747                 if(part->partfac != 0.f)
1748                         madd_v3_v3fl(vel, p_vel, part->partfac);
1749                 
1750                 mul_v3_v3fl(pa->state.vel, vel, ptex.ivel);
1751
1752                 /* -location from emitter                               */
1753                 copy_v3_v3(pa->state.co,loc);
1754
1755                 /* -rotation                                                    */
1756                 unit_qt(pa->state.rot);
1757
1758                 if(part->rotmode){
1759                         /* create vector into which rotation is aligned */
1760                         switch(part->rotmode){
1761                                 case PART_ROT_NOR:
1762                                         copy_v3_v3(rot_vec, nor);
1763                                         break;
1764                                 case PART_ROT_VEL:
1765                                         copy_v3_v3(rot_vec, vel);
1766                                         break;
1767                                 case PART_ROT_GLOB_X:
1768                                 case PART_ROT_GLOB_Y:
1769                                 case PART_ROT_GLOB_Z:
1770                                         rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
1771                                         break;
1772                                 case PART_ROT_OB_X:
1773                                 case PART_ROT_OB_Y:
1774                                 case PART_ROT_OB_Z:
1775                                         copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
1776                                         break;
1777                         }
1778                         
1779                         /* create rotation quat */
1780                         negate_v3(rot_vec);
1781                         vec_to_quat( q2,rot_vec, OB_POSX, OB_POSZ);
1782
1783                         /* randomize rotation quat */
1784                         if(part->randrotfac!=0.0f)
1785                                 interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
1786                         else
1787                                 copy_qt_qt(rot,q2);
1788
1789                         /* rotation phase */
1790                         phasefac = part->phasefac;
1791                         if(part->randphasefac != 0.0f)
1792                                 phasefac += part->randphasefac * PSYS_FRAND(p + 20);
1793                         axis_angle_to_quat( q_phase,x_vec, phasefac*(float)M_PI);
1794
1795                         /* combine base rotation & phase */
1796                         mul_qt_qtqt(pa->state.rot, rot, q_phase);
1797                 }
1798
1799                 /* -angular velocity                                    */
1800
1801                 zero_v3(pa->state.ave);
1802
1803                 if(part->avemode){
1804                         switch(part->avemode){
1805                                 case PART_AVE_SPIN:
1806                                         copy_v3_v3(pa->state.ave, vel);
1807                                         break;
1808                                 case PART_AVE_RAND:
1809                                         copy_v3_v3(pa->state.ave, r_ave);
1810                                         break;
1811                         }
1812                         normalize_v3(pa->state.ave);
1813                         mul_v3_fl(pa->state.ave,part->avefac);
1814                 }
1815         }
1816
1817
1818         if(part->type == PART_HAIR){
1819                 pa->lifetime = 100.0f;
1820         }
1821         else{
1822                 pa->lifetime = part->lifetime * ptex.life;
1823
1824                 if(part->randlife != 0.0)
1825                         pa->lifetime *= 1.0f - part->randlife * PSYS_FRAND(p + 21);
1826         }
1827
1828         pa->dietime = pa->time + pa->lifetime;
1829
1830         if(sim->psys->pointcache && sim->psys->pointcache->flag & PTCACHE_BAKED &&
1831                 sim->psys->pointcache->mem_cache.first) {
1832                 float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1833                 pa->dietime = MIN2(pa->dietime, dietime);
1834         }
1835
1836         if(pa->time > cfra)
1837                 pa->alive = PARS_UNBORN;
1838         else if(pa->dietime <= cfra)
1839                 pa->alive = PARS_DEAD;
1840         else
1841                 pa->alive = PARS_ALIVE;
1842
1843         pa->state.time = cfra;
1844 }
1845 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1846 {
1847         ParticleData *pa;
1848         int p, totpart=sim->psys->totpart;
1849         
1850         for(p=from, pa=sim->psys->particles+from; p<totpart; p++, pa++)
1851                 reset_particle(sim, pa, dtime, cfra);
1852 }
1853 /************************************************/
1854 /*                      Particle targets                                        */
1855 /************************************************/
1856 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1857 {
1858         ParticleSystem *psys = NULL;
1859
1860         if(pt->ob == NULL || pt->ob == ob)
1861                 psys = BLI_findlink(&ob->particlesystem, pt->psys-1);
1862         else
1863                 psys = BLI_findlink(&pt->ob->particlesystem, pt->psys-1);
1864
1865         if(psys)
1866                 pt->flag |= PTARGET_VALID;
1867         else
1868                 pt->flag &= ~PTARGET_VALID;
1869
1870         return psys;
1871 }
1872 /************************************************/
1873 /*                      Keyed particles                                         */
1874 /************************************************/
1875 /* Counts valid keyed targets */
1876 void psys_count_keyed_targets(ParticleSimulationData *sim)
1877 {
1878         ParticleSystem *psys = sim->psys, *kpsys;
1879         ParticleTarget *pt = psys->targets.first;
1880         int keys_valid = 1;
1881         psys->totkeyed = 0;
1882
1883         for(; pt; pt=pt->next) {
1884                 kpsys = psys_get_target_system(sim->ob, pt);
1885
1886                 if(kpsys && kpsys->totpart) {
1887                         psys->totkeyed += keys_valid;
1888                         if(psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
1889                                 psys->totkeyed += 1;
1890                 }
1891                 else {
1892                         keys_valid = 0;
1893                 }
1894         }
1895
1896         psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1897 }
1898
1899 static void set_keyed_keys(ParticleSimulationData *sim)
1900 {
1901         ParticleSystem *psys = sim->psys;
1902         ParticleSimulationData ksim= {0};
1903         ParticleTarget *pt;
1904         PARTICLE_P;
1905         ParticleKey *key;
1906         int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1907
1908         ksim.scene= sim->scene;
1909         
1910         /* no proper targets so let's clear and bail out */
1911         if(psys->totkeyed==0) {
1912                 free_keyed_keys(psys);
1913                 psys->flag &= ~PSYS_KEYED;
1914                 return;
1915         }
1916
1917         if(totpart && psys->particles->totkey != totkeys) {
1918                 free_keyed_keys(psys);
1919                 
1920                 key = MEM_callocN(totpart*totkeys*sizeof(ParticleKey), "Keyed keys");
1921                 
1922                 LOOP_PARTICLES {
1923                         pa->keys = key;
1924                         pa->totkey = totkeys;
1925                         key += totkeys;
1926                 }
1927         }
1928         
1929         psys->flag &= ~PSYS_KEYED;
1930
1931
1932         pt = psys->targets.first;
1933         for(k=0; k<totkeys; k++) {
1934                 ksim.ob = pt->ob ? pt->ob : sim->ob;
1935                 ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
1936
1937                 LOOP_PARTICLES {
1938                         key = pa->keys + k;
1939                         key->time = -1.0; /* use current time */
1940
1941                         psys_get_particle_state(&ksim, p%ksim.psys->totpart, key, 1);
1942
1943                         if(psys->flag & PSYS_KEYED_TIMING){
1944                                 key->time = pa->time + pt->time;
1945                                 if(pt->duration != 0.0f && k+1 < totkeys) {
1946                                         copy_particle_key(key+1, key, 1);
1947                                         (key+1)->time = pa->time + pt->time + pt->duration;
1948                                 }
1949                         }
1950                         else if(totkeys > 1)
1951                                 key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
1952                         else
1953                                 key->time = pa->time;
1954                 }
1955
1956                 if(psys->flag & PSYS_KEYED_TIMING && pt->duration!=0.0f)
1957                         k++;
1958
1959                 pt = (pt->next && pt->next->flag & PTARGET_VALID)? pt->next : psys->targets.first;
1960         }
1961
1962         psys->flag |= PSYS_KEYED;
1963 }
1964
1965 /************************************************/
1966 /*                      Point Cache                                                     */
1967 /************************************************/
1968 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
1969 {
1970         PointCache *cache = psys->pointcache;
1971
1972         if(cache->flag & PTCACHE_DISK_CACHE && cache->mem_cache.first == NULL) {
1973                 PTCacheID pid;
1974                 BKE_ptcache_id_from_particles(&pid, ob, psys);
1975                 cache->flag &= ~PTCACHE_DISK_CACHE;
1976                 BKE_ptcache_disk_to_mem(&pid);
1977                 cache->flag |= PTCACHE_DISK_CACHE;
1978         }
1979 }
1980 static void psys_clear_temp_pointcache(ParticleSystem *psys)
1981 {
1982         if(psys->pointcache->flag & PTCACHE_DISK_CACHE)
1983                 BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
1984 }
1985 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
1986 {
1987         ParticleSettings *part = psys->part;
1988
1989         *sfra = MAX2(1, (int)part->sta);
1990         *efra = MIN2((int)(part->end + part->lifetime + 1.0), scene->r.efra);
1991 }
1992
1993 /************************************************/
1994 /*                      Effectors                                                       */
1995 /************************************************/
1996 static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
1997 {
1998         if(psys) {
1999                 PARTICLE_P;
2000                 int totpart = 0;
2001
2002                 if(!psys->bvhtree || psys->bvhtree_frame != cfra) {
2003                         LOOP_SHOWN_PARTICLES {
2004                                 totpart++;
2005                         }
2006                         
2007                         BLI_bvhtree_free(psys->bvhtree);
2008                         psys->bvhtree = BLI_bvhtree_new(totpart, 0.0, 4, 6);
2009
2010                         LOOP_SHOWN_PARTICLES {
2011                                 if(pa->alive == PARS_ALIVE) {
2012                                         if(pa->state.time == cfra)
2013                                                 BLI_bvhtree_insert(psys->bvhtree, p, pa->prev_state.co, 1);
2014                                         else
2015                                                 BLI_bvhtree_insert(psys->bvhtree, p, pa->state.co, 1);
2016                                 }
2017                         }
2018                         BLI_bvhtree_balance(psys->bvhtree);
2019
2020                         psys->bvhtree_frame = cfra;
2021                 }
2022         }
2023 }
2024 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
2025 {
2026         if(psys) {
2027                 PARTICLE_P;
2028                 int totpart = 0;
2029
2030                 if(!psys->tree || psys->tree_frame != cfra) {
2031                         LOOP_SHOWN_PARTICLES {
2032                                 totpart++;
2033                         }
2034
2035                         BLI_kdtree_free(psys->tree);
2036                         psys->tree = BLI_kdtree_new(psys->totpart);
2037
2038                         LOOP_SHOWN_PARTICLES {
2039                                 if(pa->alive == PARS_ALIVE) {
2040                                         if(pa->state.time == cfra)
2041                                                 BLI_kdtree_insert(psys->tree, p, pa->prev_state.co, NULL);
2042                                         else
2043                                                 BLI_kdtree_insert(psys->tree, p, pa->state.co, NULL);
2044                                 }
2045                         }
2046                         BLI_kdtree_balance(psys->tree);
2047
2048                         psys->tree_frame = cfra;
2049                 }
2050         }
2051 }
2052
2053 static void psys_update_effectors(ParticleSimulationData *sim)
2054 {
2055         pdEndEffectors(&sim->psys->effectors);
2056         sim->psys->effectors = pdInitEffectors(sim->scene, sim->ob, sim->psys, sim->psys->part->effector_weights);
2057         precalc_guides(sim, sim->psys->effectors);
2058 }
2059
2060 static void integrate_particle(ParticleSettings *part, ParticleData *pa, float dtime, float *external_acceleration, void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse), void *forcedata)
2061 {
2062         ParticleKey states[5];
2063         float force[3],acceleration[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
2064         float pa_mass= (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
2065         int i, steps=1;
2066         int integrator = part->integrator;
2067
2068         copy_v3_v3(oldpos, pa->state.co);
2069         
2070         /* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
2071         if(pa->prev_state.time < 0.f && integrator == PART_INT_VERLET)
2072                 integrator = PART_INT_EULER;
2073
2074         switch(integrator){
2075                 case PART_INT_EULER:
2076                         steps=1;
2077                         break;
2078                 case PART_INT_MIDPOINT:
2079                         steps=2;
2080                         break;
2081                 case PART_INT_RK4:
2082                         steps=4;
2083                         break;
2084                 case PART_INT_VERLET:
2085                         steps=1;
2086                         break;
2087         }
2088
2089         copy_particle_key(states, &pa->state, 1);
2090
2091         states->time = 0.f;
2092
2093         for(i=0; i<steps; i++){
2094                 zero_v3(force);
2095                 zero_v3(impulse);
2096
2097                 force_func(forcedata, states+i, force, impulse);
2098
2099                 /* force to acceleration*/
2100                 mul_v3_v3fl(acceleration, force, 1.0f/pa_mass);
2101
2102                 if(external_acceleration)
2103                         add_v3_v3(acceleration, external_acceleration);
2104                 
2105                 /* calculate next state */
2106                 add_v3_v3(states[i].vel, impulse);
2107
2108                 switch(integrator){
2109                         case PART_INT_EULER:
2110                                 madd_v3_v3v3fl(pa->state.co, states->co, states->vel, dtime);
2111                                 madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2112                                 break;
2113                         case PART_INT_MIDPOINT:
2114                                 if(i==0){
2115                                         madd_v3_v3v3fl(states[1].co, states->co, states->vel, dtime*0.5f);
2116                                         madd_v3_v3v3fl(states[1].vel, states->vel, acceleration, dtime*0.5f);
2117                                         states[1].time = dtime*0.5f;
2118                                         /*fra=sim->psys->cfra+0.5f*dfra;*/
2119                                 }
2120                                 else{
2121                                         madd_v3_v3v3fl(pa->state.co, states->co, states[1].vel, dtime);
2122                                         madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2123                                 }
2124                                 break;
2125                         case PART_INT_RK4:
2126                                 switch(i){
2127                                         case 0:
2128                                                 copy_v3_v3(dx[0], states->vel);
2129                                                 mul_v3_fl(dx[0], dtime);
2130                                                 copy_v3_v3(dv[0], acceleration);
2131                                                 mul_v3_fl(dv[0], dtime);
2132
2133                                                 madd_v3_v3v3fl(states[1].co, states->co, dx[0], 0.5f);
2134                                                 madd_v3_v3v3fl(states[1].vel, states->vel, dv[0], 0.5f);
2135                                                 states[1].time = dtime*0.5f;
2136                                                 /*fra=sim->psys->cfra+0.5f*dfra;*/
2137                                                 break;
2138                                         case 1:
2139                                                 madd_v3_v3v3fl(dx[1], states->vel, dv[0], 0.5f);
2140                                                 mul_v3_fl(dx[1], dtime);
2141                                                 copy_v3_v3(dv[1], acceleration);
2142                                                 mul_v3_fl(dv[1], dtime);
2143
2144                                                 madd_v3_v3v3fl(states[2].co, states->co, dx[1], 0.5f);
2145                                                 madd_v3_v3v3fl(states[2].vel, states->vel, dv[1], 0.5f);
2146                                                 states[2].time = dtime*0.5f;
2147                                                 break;
2148                                         case 2:
2149                                                 madd_v3_v3v3fl(dx[2], states->vel, dv[1], 0.5f);
2150                                                 mul_v3_fl(dx[2], dtime);
2151                                                 copy_v3_v3(dv[2], acceleration);
2152                                                 mul_v3_fl(dv[2], dtime);
2153
2154                                                 add_v3_v3v3(states[3].co, states->co, dx[2]);
2155                                                 add_v3_v3v3(states[3].vel, states->vel, dv[2]);
2156                                                 states[3].time = dtime;
2157                                                 /*fra=cfra;*/
2158                                                 break;
2159                                         case 3:
2160                                                 add_v3_v3v3(dx[3], states->vel, dv[2]);
2161                                                 mul_v3_fl(dx[3], dtime);
2162                                                 copy_v3_v3(dv[3], acceleration);
2163                                                 mul_v3_fl(dv[3], dtime);
2164
2165                                                 madd_v3_v3v3fl(pa->state.co, states->co, dx[0], 1.0f/6.0f);
2166                                                 madd_v3_v3fl(pa->state.co, dx[1], 1.0f/3.0f);
2167                                                 madd_v3_v3fl(pa->state.co, dx[2], 1.0f/3.0f);
2168                                                 madd_v3_v3fl(pa->state.co, dx[3], 1.0f/6.0f);
2169
2170                                                 madd_v3_v3v3fl(pa->state.vel, states->vel, dv[0], 1.0f/6.0f);
2171                                                 madd_v3_v3fl(pa->state.vel, dv[1], 1.0f/3.0f);
2172                                                 madd_v3_v3fl(pa->state.vel, dv[2], 1.0f/3.0f);
2173                                                 madd_v3_v3fl(pa->state.vel, dv[3], 1.0f/6.0f);
2174                                 }
2175                                 break;
2176                         case PART_INT_VERLET:   /* Verlet integration */
2177                                 madd_v3_v3v3fl(pa->state.vel, pa->prev_state.vel, acceleration, dtime);
2178                                 madd_v3_v3v3fl(pa->state.co, pa->prev_state.co, pa->state.vel, dtime);
2179
2180                                 sub_v3_v3v3(pa->state.vel, pa->state.co, oldpos);
2181                                 mul_v3_fl(pa->state.vel, 1.0f/dtime);
2182                                 break;
2183                 }
2184         }
2185 }
2186
2187 /*********************************************************************************************************
2188                     SPH fluid physics 
2189
2190  In theory, there could be unlimited implementation of SPH simulators
2191
2192  This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
2193
2194  Titled: Particle-based Viscoelastic Fluid Simulation.
2195  Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
2196  Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
2197
2198  Presented at Siggraph, (2005)
2199
2200 ***********************************************************************************************************/
2201 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
2202 static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
2203 {
2204         /* Are more refs required? */
2205         if(psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
2206                 psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
2207                 psys->fluid_springs = (ParticleSpring*)MEM_callocN(psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
2208         }
2209         else if(psys->tot_fluidsprings == psys->alloc_fluidsprings) {
2210                 /* Double the number of refs allocated */
2211                 psys->alloc_fluidsprings *= 2;
2212                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
2213         }
2214
2215         memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
2216         psys->tot_fluidsprings++;
2217
2218         return psys->fluid_springs + psys->tot_fluidsprings - 1;
2219 }
2220 static void sph_spring_delete(ParticleSystem *psys, int j)
2221 {
2222         if (j != psys->tot_fluidsprings - 1)
2223                 psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
2224
2225         psys->tot_fluidsprings--;
2226
2227         if (psys->tot_fluidsprings < psys->alloc_fluidsprings/2 && psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE){
2228                 psys->alloc_fluidsprings /= 2;
2229                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs,  psys->alloc_fluidsprings * sizeof(ParticleSpring));
2230         }
2231 }
2232 static void sph_springs_modify(ParticleSystem *psys, float dtime){
2233         SPHFluidSettings *fluid = psys->part->fluid;
2234         ParticleData *pa1, *pa2;
2235         ParticleSpring *spring = psys->fluid_springs;
2236         
2237         float h, d, Rij[3], rij, Lij;
2238         int i;
2239
2240         float yield_ratio = fluid->yield_ratio;
2241         float plasticity = fluid->plasticity_constant;
2242         /* scale things according to dtime */
2243         float timefix = 25.f * dtime;
2244
2245         if((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
2246                 return;
2247
2248         /* Loop through the springs */
2249         for(i=0; i<psys->tot_fluidsprings; i++, spring++) {
2250                 pa1 = psys->particles + spring->particle_index[0];
2251                 pa2 = psys->particles + spring->particle_index[1];
2252
2253                 sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
2254                 rij = normalize_v3(Rij);
2255
2256                 /* adjust rest length */
2257                 Lij = spring->rest_length;
2258                 d = yield_ratio * timefix * Lij;
2259
2260                 if (rij > Lij + d) // Stretch
2261                         spring->rest_length += plasticity * (rij - Lij - d) * timefix;
2262                 else if(rij < Lij - d) // Compress
2263                         spring->rest_length -= plasticity * (Lij - d - rij) * timefix;
2264
2265                 h = 4.f*pa1->size;
2266
2267                 if(spring->rest_length > h)
2268                         spring->delete_flag = 1;
2269         }
2270
2271         /* Loop through springs backwaqrds - for efficient delete function */
2272         for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
2273                 if(psys->fluid_springs[i].delete_flag)
2274                         sph_spring_delete(psys, i);
2275         }
2276 }
2277 static EdgeHash *sph_springhash_build(ParticleSystem *psys)
2278 {
2279         EdgeHash *springhash = NULL;
2280         ParticleSpring *spring;
2281         int i = 0;
2282
2283         springhash = BLI_edgehash_new();
2284
2285         for(i=0, spring=psys->fluid_springs; i<psys->tot_fluidsprings; i++, spring++)
2286                 BLI_edgehash_insert(springhash, spring->particle_index[0], spring->particle_index[1], SET_INT_IN_POINTER(i+1));
2287
2288         return springhash;
2289 }
2290
2291 typedef struct SPHNeighbor
2292 {
2293         ParticleSystem *psys;
2294         int index;
2295 } SPHNeighbor;
2296 typedef struct SPHRangeData
2297 {
2298         SPHNeighbor neighbors[128];
2299         int tot_neighbors;
2300
2301         float density, near_density;
2302         float h;
2303
2304         ParticleSystem *npsys;
2305         ParticleData *pa;
2306
2307         float massfac;
2308         int use_size;
2309 } SPHRangeData;
2310 typedef struct SPHData {
2311         ParticleSystem *psys[10];
2312         ParticleData *pa;
2313         float mass;
2314         EdgeHash *eh;
2315         float *gravity;
2316 }SPHData;
2317 static void sph_density_accum_cb(void *userdata, int index, float squared_dist)
2318 {
2319         SPHRangeData *pfr = (SPHRangeData *)userdata;
2320         ParticleData *npa = pfr->npsys->particles + index;
2321         float q;
2322
2323         if(npa == pfr->pa || squared_dist < FLT_EPSILON)
2324                 return;
2325
2326         /* Ugh! One particle has over 128 neighbors! Really shouldn't happen,
2327          * but even if it does it shouldn't do any terrible harm if all are
2328          * not taken into account - jahka
2329          */
2330         if(pfr->tot_neighbors >= 128)
2331                 return;
2332         
2333         pfr->neighbors[pfr->tot_neighbors].index = index;
2334         pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
2335         pfr->tot_neighbors++;
2336
2337         q = (1.f - sqrt(squared_dist)/pfr->h) * pfr->massfac;
2338
2339         if(pfr->use_size)
2340                 q *= npa->size;
2341
2342         pfr->density += q*q;
2343         pfr->near_density += q*q*q;
2344 }
2345 static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
2346 {
2347         SPHData *sphdata = (SPHData *)sphdata_v;
2348         ParticleSystem **psys = sphdata->psys;
2349         ParticleData *pa = sphdata->pa;
2350         SPHFluidSettings *fluid = psys[0]->part->fluid;
2351         ParticleSpring *spring = NULL;
2352         SPHRangeData pfr;
2353         SPHNeighbor *pfn;
2354         float mass = sphdata->mass;
2355         float *gravity = sphdata->gravity;
2356         EdgeHash *springhash = sphdata->eh;
2357
2358         float q, u, rij, dv[3];
2359         float pressure, near_pressure;
2360
2361         float visc = fluid->viscosity_omega;
2362         float stiff_visc = fluid->viscosity_beta * (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
2363
2364         float inv_mass = 1.0f/mass;
2365         float spring_constant = fluid->spring_k;
2366         
2367         float h = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.f*pa->size : 1.f); /* 4.0 seems to be a pretty good value */
2368         float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f); /* 4.77 is an experimentally determined density factor */
2369         float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
2370
2371         float stiffness = fluid->stiffness_k;
2372         float stiffness_near_fac = fluid->stiffness_knear * (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);
2373
2374         ParticleData *npa;
2375         float vec[3];
2376         float vel[3];
2377         float co[3];
2378
2379         int i, spring_index, index = pa - psys[0]->particles;
2380
2381         pfr.tot_neighbors = 0;
2382         pfr.density = pfr.near_density = 0.f;
2383         pfr.h = h;
2384         pfr.pa = pa;
2385
2386         for(i=0; i<10 && psys[i]; i++) {
2387                 pfr.npsys = psys[i];
2388                 pfr.massfac = psys[i]->part->mass*inv_mass;
2389                 pfr.use_size = psys[i]->part->flag & PART_SIZEMASS;
2390
2391                 BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sph_density_accum_cb, &pfr);
2392         }
2393
2394         pressure =  stiffness * (pfr.density - rest_density);
2395         near_pressure = stiffness_near_fac * pfr.near_density;
2396
2397         pfn = pfr.neighbors;
2398         for(i=0; i<pfr.tot_neighbors; i++, pfn++) {
2399                 npa = pfn->psys->particles + pfn->index;
2400
2401                 madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
2402
2403                 sub_v3_v3v3(vec, co, state->co);
2404                 rij = normalize_v3(vec);
2405
2406                 q = (1.f - rij/h) * pfn->psys->part->mass * inv_mass;
2407
2408                 if(pfn->psys->part->flag & PART_SIZEMASS)
2409                         q *= npa->size;
2410
2411                 copy_v3_v3(vel, npa->prev_state.vel);
2412
2413                 /* Double Density Relaxation */
2414                 madd_v3_v3fl(force, vec, -(pressure + near_pressure*q)*q);
2415
2416                 /* Viscosity */
2417                 if(visc > 0.f   || stiff_visc > 0.f) {          
2418                         sub_v3_v3v3(dv, vel, state->vel);
2419                         u = dot_v3v3(vec, dv);
2420
2421                         if(u < 0.f && visc > 0.f)
2422                                 madd_v3_v3fl(force, vec, 0.5f * q * visc * u );
2423
2424                         if(u > 0.f && stiff_visc > 0.f)
2425                                 madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u );
2426                 }
2427
2428                 if(spring_constant > 0.f) {
2429                         /* Viscoelastic spring force */
2430                         if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
2431                                 spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, pfn->index));
2432
2433                                 if(spring_index) {
2434                                         spring = psys[0]->fluid_springs + spring_index - 1;
2435
2436                                         madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (spring->rest_length - rij));
2437                                 }
2438                                 else if(fluid->spring_frames == 0 || (pa->prev_state.time-pa->time) <= fluid->spring_frames){
2439                                         ParticleSpring temp_spring;
2440                                         temp_spring.particle_index[0] = index;
2441                                         temp_spring.particle_index[1] = pfn->index;
2442                                         temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length;
2443                                         temp_spring.delete_flag = 0;
2444                                                                 
2445                                         sph_spring_add(psys[0], &temp_spring);
2446                                 }
2447                         }
2448                         else {/* PART_SPRING_HOOKES - Hooke's spring force */
2449                                 madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (rest_length - rij));
2450                         }
2451                 }
2452         }
2453         
2454         /* Artificial buoyancy force in negative gravity direction  */
2455         if (fluid->buoyancy > 0.f && gravity)
2456                 madd_v3_v3fl(force, gravity, fluid->buoyancy * (pfr.density-rest_density));
2457 }
2458
2459 static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, float *gravity, EdgeHash *springhash){
2460         ParticleTarget *pt;
2461         int i;
2462
2463         ParticleSettings *part = sim->psys->part;
2464         // float timestep = psys_get_timestep(sim); // UNUSED
2465         float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
2466         float dtime = dfra*psys_get_timestep(sim);
2467         // int steps = 1; // UNUSED
2468         float effector_acceleration[3];
2469         SPHData sphdata;
2470
2471         sphdata.psys[0] = sim->psys;
2472         for(i=1, pt=sim->psys->targets.first; i<10; i++, pt=(pt?pt->next:NULL))
2473                 sphdata.psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL;
2474
2475         sphdata.pa = pa;
2476         sphdata.gravity = gravity;
2477         sphdata.mass = pa_mass;
2478         sphdata.eh = springhash;
2479
2480         /* restore previous state and treat gravity & effectors as external acceleration*/
2481         sub_v3_v3v3(effector_acceleration, pa->state.vel, pa->prev_state.vel);
2482         mul_v3_fl(effector_acceleration, 1.f/dtime);
2483
2484         copy_particle_key(&pa->state, &pa->prev_state, 0);
2485
2486         integrate_particle(part, pa, dtime, effector_acceleration, sph_force_cb, &sphdata);
2487 }
2488
2489 /************************************************/
2490 /*                      Basic physics                                           */
2491 /************************************************/
2492 typedef struct EfData
2493 {
2494         ParticleTexture ptex;
2495         ParticleSimulationData *sim;
2496         ParticleData *pa;
2497 } EfData;
2498 static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, float *impulse)
2499 {
2500         EfData *efdata = (EfData *)efdata_v;
2501         ParticleSimulationData *sim = efdata->sim;
2502         ParticleSettings *part = sim->psys->part;
2503         ParticleData *pa = efdata->pa;
2504         EffectedPoint epoint;
2505
2506         /* add effectors */
2507         pd_point_from_particle(efdata->sim, efdata->pa, state, &epoint);
2508         if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
2509                 pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
2510
2511         mul_v3_fl(force, efdata->ptex.field);
2512         mul_v3_fl(impulse, efdata->ptex.field);
2513
2514         /* calculate air-particle interaction */
2515         if(part->dragfac != 0.0f)
2516                 madd_v3_v3fl(force, state->vel, -part->dragfac * pa->size * pa->size * len_v3(state->vel));
2517
2518         /* brownian force */
2519         if(part->brownfac != 0.0){
2520                 force[0] += (BLI_frand()-0.5f) * part->brownfac;
2521                 force[1] += (BLI_frand()-0.5f) * part->brownfac;
2522                 force[2] += (BLI_frand()-0.5f) * part->brownfac;
2523         }
2524 }
2525 /* gathers all forces that effect particles and calculates a new state for the particle */
2526 static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
2527 {
2528         ParticleSettings *part = sim->psys->part;
2529         ParticleData *pa = sim->psys->particles + p;
2530         ParticleKey tkey;
2531         float dtime=dfra*psys_get_timestep(sim), time;
2532         float *gravity = NULL, gr[3];
2533         EfData efdata;
2534
2535         psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);
2536
2537         efdata.pa = pa;
2538         efdata.sim = sim;
2539
2540         /* add global acceleration (gravitation) */
2541         if(psys_uses_gravity(sim)
2542                 /* normal gravity is too strong for hair so it's disabled by default */
2543                 && (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
2544                 zero_v3(gr);
2545                 madd_v3_v3fl(gr, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * efdata.ptex.gravity);
2546                 gravity = gr;
2547         }
2548
2549         /* maintain angular velocity */
2550         copy_v3_v3(pa->state.ave, pa->prev_state.ave);
2551
2552         integrate_particle(part, pa, dtime, gravity, basic_force_cb, &efdata);
2553
2554         /* damp affects final velocity */
2555         if(part->dampfac != 0.f)
2556                 mul_v3_fl(pa->state.vel, 1.f - part->dampfac * efdata.ptex.damp);
2557
2558         //VECCOPY(pa->state.ave, states->ave);
2559
2560         /* finally we do guides */
2561         time=(cfra-pa->time)/pa->lifetime;
2562         CLAMP(time,0.0,1.0);
2563
2564         VECCOPY(tkey.co,pa->state.co);
2565         VECCOPY(tkey.vel,pa->state.vel);
2566         tkey.time=pa->state.time;
2567
2568         if(part->type != PART_HAIR) {
2569                 if(do_guides(sim->psys->effectors, &tkey, p, time)) {
2570                         VECCOPY(pa->state.co,tkey.co);
2571                         /* guides don't produce valid velocity */
2572                         VECSUB(pa->state.vel,tkey.co,pa->prev_state.co);
2573                         mul_v3_fl(pa->state.vel,1.0f/dtime);
2574                         pa->state.time=tkey.time;
2575                 }
2576         }
2577 }
2578 static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2579 {
2580         float rotfac, rot1[4], rot2[4]={1.0,0.0,0.0,0.0}, dtime=dfra*timestep;
2581
2582         if((part->flag & PART_ROT_DYN)==0){
2583                 if(part->avemode==PART_AVE_SPIN){
2584                         float angle;
2585                         float len1 = len_v3(pa->prev_state.vel);
2586                         float len2 = len_v3(pa->state.vel);
2587
2588                         if(len1==0.0f || len2==0.0f)
2589                                 pa->state.ave[0]=pa->state.ave[1]=pa->state.ave[2]=0.0f;
2590                         else{
2591                                 cross_v3_v3v3(pa->state.ave,pa->prev_state.vel,pa->state.vel);
2592                                 normalize_v3(pa->state.ave);
2593                                 angle=dot_v3v3(pa->prev_state.vel,pa->state.vel)/(len1*len2);
2594                                 mul_v3_fl(pa->state.ave,saacos(angle)/dtime);
2595                         }
2596
2597                         axis_angle_to_quat(rot2,pa->state.vel,dtime*part->avefac);
2598                 }
2599         }
2600
2601         rotfac=len_v3(pa->state.ave);
2602         if(rotfac==0.0){ /* unit_qt(in VecRotToQuat) doesn't give unit quat [1,0,0,0]?? */
2603                 rot1[0]=1.0;
2604                 rot1[1]=rot1[2]=rot1[3]=0;
2605         }
2606         else{
2607                 axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
2608         }
2609         mul_qt_qtqt(pa->state.rot,rot1,pa->prev_state.rot);
2610         mul_qt_qtqt(pa->state.rot,rot2,pa->state.rot);
2611
2612         /* keep rotation quat in good health */
2613         normalize_qt(pa->state.rot);
2614 }
2615
2616 /************************************************/
2617 /*                      Collisions                                                      */
2618 /************************************************/
2619 #define COLLISION_MAX_COLLISIONS        10
2620 #define COLLISION_MIN_RADIUS 0.001f
2621 #define COLLISION_MIN_DISTANCE 0.0001f
2622 #define COLLISION_ZERO 0.00001f
2623 typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
2624 static float nr_signed_distance_to_plane(float *p, float radius, ParticleCollisionElement *pce, float *nor)
2625 {
2626         float p0[3], e1[3], e2[3], d;
2627
2628         sub_v3_v3v3(e1, pce->x1, pce->x0);
2629         sub_v3_v3v3(e2, pce->x2, pce->x0);
2630         sub_v3_v3v3(p0, p, pce->x0);
2631
2632         cross_v3_v3v3(nor, e1, e2);
2633         normalize_v3(nor);
2634
2635         d = dot_v3v3(p0, nor);
2636
2637         if(pce->inv_nor == -1) {
2638                 if(d < 0.f)
2639                         pce->inv_nor = 1;
2640                 else
2641                         pce->inv_nor = 0;
2642         }
2643
2644         if(pce->inv_nor == 1) {
2645                 mul_v3_fl(nor, -1.f);
2646                 d = -d;
2647         }
2648
2649         return d - radius;
2650 }
2651 static float nr_distance_to_edge(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(nor))
2652 {
2653         float v0[3], v1[3], v2[3], c[3];
2654
2655         sub_v3_v3v3(v0, pce->x1, pce->x0);
2656         sub_v3_v3v3(v1, p, pce->x0);
2657         sub_v3_v3v3(v2, p, pce->x1);
2658
2659         cross_v3_v3v3(c, v1, v2);
2660
2661         return fabs(len_v3(c)/len_v3(v0)) - radius;
2662 }
2663 static float nr_distance_to_vert(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(nor))
2664 {
2665         return len_v3v3(p, pce->x0) - radius;
2666 }
2667 static void collision_interpolate_element(ParticleCollisionElement *pce, float t, float fac, ParticleCollision *col)
2668 {
2669         /* t is the current time for newton rhapson */
2670         /* fac is the starting factor for current collision iteration */
2671         /* the col->fac's are factors for the particle subframe step start and end during collision modifier step */
2672         float f = fac + t*(1.f-fac);
2673         float mul = col->fac1 + f * (col->fac2-col->fac1);
2674         if(pce->tot > 0) {
2675                 madd_v3_v3v3fl(pce->x0, pce->x[0], pce->v[0], mul);
2676
2677                 if(pce->tot > 1) {
2678                         madd_v3_v3v3fl(pce->x1, pce->x[1], pce->v[1], mul);
2679
2680                         if(pce->tot > 2)
2681                                 madd_v3_v3v3fl(pce->x2, pce->x[2], pce->v[2], mul);
2682                 }
2683         }
2684 }
2685 static void collision_point_velocity(ParticleCollisionElement *pce)
2686 {
2687         float v[3];
2688
2689         copy_v3_v3(pce->vel, pce->v[0]);
2690
2691         if(pce->tot > 1) {
2692                 sub_v3_v3v3(v, pce->v[1], pce->v[0]);
2693                 madd_v3_v3fl(pce->vel, v, pce->uv[0]);
2694
2695                 if(pce->tot > 2) {
2696                         sub_v3_v3v3(v, pce->v[2], pce->v[0]);
2697                         madd_v3_v3fl(pce->vel, v, pce->uv[1]);
2698                 }
2699         }
2700 }
2701 static float collision_point_distance_with_normal(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *nor)
2702 {
2703         if(fac >= 0.f)
2704                 collision_interpolate_element(pce, 0.f, fac, col);
2705
2706         switch(pce->tot) {
2707                 case 1:
2708                 {
2709                         sub_v3_v3v3(nor, p, pce->x0);
2710                         return normalize_v3(nor);
2711                 }
2712                 case 2:
2713                 {
2714                         float u, e[3], vec[3];
2715                         sub_v3_v3v3(e, pce->x1, pce->x0);
2716                         sub_v3_v3v3(vec, p, pce->x0);
2717                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2718
2719                         madd_v3_v3v3fl(nor, vec, e, -u);
2720                         return normalize_v3(nor);
2721                 }
2722                 case 3:
2723                         return nr_signed_distance_to_plane(p, 0.f, pce, nor);
2724         }
2725         return 0;
2726 }
2727 static void collision_point_on_surface(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *co)
2728 {
2729         collision_interpolate_element(pce, 0.f, fac, col);
2730
2731         switch(pce->tot) {
2732                 case 1:
2733                 {
2734                         sub_v3_v3v3(co, p, pce->x0);
2735                         normalize_v3(co);
2736                         madd_v3_v3v3fl(co, pce->x0, co, col->radius);
2737                         break;
2738                 }
2739                 case 2:
2740                 {
2741                         float u, e[3], vec[3], nor[3];
2742                         sub_v3_v3v3(e, pce->x1, pce->x0);
2743                         sub_v3_v3v3(vec, p, pce->x0);
2744                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2745
2746                         madd_v3_v3v3fl(nor, vec, e, -u);
2747                         normalize_v3(nor);
2748
2749                         madd_v3_v3v3fl(co, pce->x0, e, pce->uv[0]);
2750                         madd_v3_v3fl(co, nor, col->radius);
2751                         break;
2752                 }
2753                 case 3:
2754                 {
2755                                 float p0[3], e1[3], e2[3], nor[3];
2756
2757                                 sub_v3_v3v3(e1, pce->x1, pce->x0);
2758                                 sub_v3_v3v3(e2, pce->x2, pce->x0);
2759                                 sub_v3_v3v3(p0, p, pce->x0);
2760
2761                                 cross_v3_v3v3(nor, e1, e2);
2762                                 normalize_v3(nor);
2763
2764                                 if(pce->inv_nor == 1)
2765                                         mul_v3_fl(nor, -1.f);
2766
2767                                 madd_v3_v3v3fl(co, pce->x0, nor, col->radius);
2768                                 madd_v3_v3fl(co, e1, pce->uv[0]);
2769                                 madd_v3_v3fl(co, e2, pce->uv[1]);
2770                         break;
2771                 }
2772         }
2773 }
2774 /* find first root in range [0-1] starting from 0 */
2775 static float collision_newton_rhapson(ParticleCollision *col, float radius, ParticleCollisionElement *pce, NRDistanceFunc distance_func)
2776 {
2777         float t0, t1, d0, d1, dd, n[3];
2778         int iter;
2779
2780         pce->inv_nor = -1;
2781
2782         /* start from the beginning */
2783         t0 = 0.f;
2784         collision_interpolate_element(pce, t0, col->f, col);
2785         d0 = distance_func(col->co1, radius, pce, n);
2786         t1 = 0.001f;
2787         d1 = 0.f;
2788
2789         for(iter=0; iter<10; iter++) {//, itersum++) {
2790                 /* get current location */
2791                 collision_interpolate_element(pce, t1, col->f, col);
2792                 interp_v3_v3v3(pce->p, col->co1, col->co2, t1);
2793
2794                 d1 = distance_func(pce->p, radius, pce, n);
2795
2796                 /* no movement, so no collision */
2797                 if(d1 == d0) {
2798                         return -1.f;
2799                 }
2800
2801                 /* particle already inside face, so report collision */
2802                 if(iter == 0 && d0 < 0.f && d0 > -radius) {
2803                         copy_v3_v3(pce->p, col->co1);
2804                         copy_v3_v3(pce->nor, n);
2805                         pce->inside = 1;
2806                         return 0.f;
2807                 }
2808                 
2809                 dd = (t1-t0)/(d1-d0);
2810
2811                 t0 = t1;
2812                 d0 = d1;
2813
2814                 t1 -= d1*dd;
2815
2816                 /* particle movin away from plane could also mean a strangely rotating face, so check from end */
2817                 if(iter == 0 && t1 < 0.f) {
2818                         t0 = 1.f;
2819                         collision_interpolate_element(pce, t0, col->f, col);
2820                         d0 = distance_func(col->co2, radius, pce, n);
2821                         t1 = 0.999f;
2822                         d1 = 0.f;
2823
2824                         continue;
2825                 }
2826                 else if(iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f))
2827                         return -1.f;
2828
2829                 if(d1 <= COLLISION_ZERO && d1 >= -COLLISION_ZERO) {
2830                         if(t1 >= -COLLISION_ZERO && t1 <= 1.f) {
2831                                 if(distance_func == nr_signed_distance_to_plane)
2832                                         copy_v3_v3(pce->nor, n);
2833
2834                                 CLAMP(t1, 0.f, 1.f);
2835
2836                                 return t1;
2837                         }
2838                         else
2839                                 return -1.f;
2840                 }
2841         }
2842         return -1.0;
2843 }
2844 static int collision_sphere_to_tri(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2845 {
2846         ParticleCollisionElement *result = &col->pce;
2847         float ct, u, v;
2848
2849         pce->inv_nor = -1;
2850         pce->inside = 0;
2851
2852         ct = collision_newton_rhapson(col, radius, pce, nr_signed_distance_to_plane);
2853
2854         if(ct >= 0.f && ct < *t && (result->inside==0 || pce->inside==1) ) {
2855                 float e1[3], e2[3], p0[3];
2856                 float e1e1, e1e2, e1p0, e2e2, e2p0, inv;
2857
2858                 sub_v3_v3v3(e1, pce->x1, pce->x0);
2859                 sub_v3_v3v3(e2, pce->x2, pce->x0);
2860                 /* XXX: add radius correction here? */
2861                 sub_v3_v3v3(p0, pce->p, pce->x0);
2862
2863                 e1e1 = dot_v3v3(e1, e1);
2864                 e1e2 = dot_v3v3(e1, e2);
2865                 e1p0 = dot_v3v3(e1, p0);
2866                 e2e2 = dot_v3v3(e2, e2);
2867                 e2p0 = dot_v3v3(e2, p0);
2868
2869                 inv = 1.f/(e1e1 * e2e2 - e1e2 * e1e2);
2870                 u = (e2e2 * e1p0 - e1e2 * e2p0) * inv;
2871                 v = (e1e1 * e2p0 - e1e2 * e1p0) * inv;
2872
2873                 if(u>=0.f && u<=1.f && v>=0.f && u+v<=1.f) {
2874                         *result = *pce;
2875
2876                         /* normal already calculated in pce */
2877
2878                         result->uv[0] = u;
2879                         result->uv[1] = v;
2880
2881                         *t = ct;
2882                         return 1;
2883                 }
2884         }
2885         return 0;
2886 }
2887 static int collision_sphere_to_edges(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2888 {
2889         ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
2890         ParticleCollisionElement *result = &col->pce;
2891
2892         float ct;
2893         int i;
2894
2895         for(i=0; i<3; i++) {
2896                 /* in case of a quad, no need to check "edge" that goes through face twice */
2897                 if((pce->x[3] && i==2))
2898                         continue;
2899
2900                 cur = edge+i;
2901                 cur->x[0] = pce->x[i]; cur->x[1] = pce->x[(i+1)%3];
2902                 cur->v[0] = pce->v[i]; cur->v[1] = pce->v[(i+1)%3];
2903                 cur->tot = 2;
2904                 cur->inside = 0;
2905
2906                 ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_edge);
2907
2908                 if(ct >= 0.f && ct < *t) {
2909                         float u, e[3], vec[3];
2910
2911                         sub_v3_v3v3(e, cur->x1, cur->x0);
2912                         sub_v3_v3v3(vec, cur->p, cur->x0);
2913                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2914
2915                         if(u < 0.f || u > 1.f)
2916                                 break;
2917
2918                         *result = *cur;
2919
2920                         madd_v3_v3v3fl(result->nor, vec, e, -u);
2921                         normalize_v3(result->nor);
2922
2923                         result->uv[0] = u;
2924
2925                         
2926                         hit = cur;
2927                         *t = ct;
2928                 }
2929
2930         }
2931
2932         return hit != NULL;
2933 }
2934 static int collision_sphere_to_verts(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2935 {
2936         ParticleCollisionElement vert[3], *cur = NULL, *hit = NULL;
2937         ParticleCollisionElement *result = &col->pce;
2938
2939         float ct;
2940         int i;
2941
2942         for(i=0; i<3; i++) {
2943                 /* in case of quad, only check one vert the first time */
2944                 if(pce->x[3] && i != 1)
2945                         continue;
2946
2947                 cur = vert+i;
2948                 cur->x[0] = pce->x[i];
2949                 cur->v[0] = pce->v[i];
2950                 cur->tot = 1;
2951                 cur->inside = 0;
2952
2953                 ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_vert);
2954                 
2955                 if(ct >= 0.f && ct < *t) {
2956                         *result = *cur;
2957
2958                         sub_v3_v3v3(result->nor, cur->p, cur->x0);
2959                         normalize_v3(result->nor);
2960
2961                         hit = cur;
2962                         *t = ct;
2963                 }
2964
2965         }
2966
2967         return hit != NULL;
2968 }
2969 /* Callback for BVHTree near test */
2970 void BKE_psys_collision_neartest_cb(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
2971 {
2972         ParticleCollision *col = (ParticleCollision *) userdata;
2973         ParticleCollisionElement pce;
2974         MFace *face = col->md->mfaces + index;
2975         MVert *x = col->md->x;
2976         MVert *v = col->md->current_v;
2977         float t = hit->dist/col->original_ray_length;
2978         int collision = 0;
2979
2980         pce.x[0] = x[face->v1].co;
2981         pce.x[1] = x[face->v2].co;
2982         pce.x[2] = x[face->v3].co;
2983         pce.x[3] = face->v4 ? x[face->v4].co : NULL;
2984
2985         pce.v[0] = v[face->v1].co;
2986         pce.v[1] = v[face->v2].co;
2987         pce.v[2] = v[face->v3].co;
2988         pce.v[3] = face->v4 ? v[face->v4].co : NULL;
2989
2990         pce.tot = 3;
2991         pce.inside = 0;
2992
2993         do
2994         {       
2995                 collision = collision_sphere_to_tri(col, ray->radius, &pce, &t);
2996                 if(col->pce.inside == 0) {
2997                         collision += collision_sphere_to_edges(col, ray->radius, &pce, &t);
2998                         collision += collision_sphere_to_verts(col, ray->radius, &pce, &t);
2999                 }
3000
3001                 if(collision) {
3002                         hit->dist = col->original_ray_length * t;
3003                         hit->index = index;
3004                                 
3005                         collision_point_velocity(&col->pce);
3006
3007                         col->hit = col->current;
3008                 }
3009
3010                 pce.x[1] = pce.x[2];
3011                 pce.x[2] = pce.x[3];
3012                 pce.x[3] = NULL;
3013
3014                 pce.v[1] = pce.v[2];
3015                 pce.v[2] = pce.v[3];
3016                 pce.v[3] = NULL;
3017
3018         } while(pce.x[2]);
3019 }
3020 static int collision_detect(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, ListBase *colliders)
3021 {
3022         ColliderCache *coll;
3023         float ray_dir[3];
3024
3025         if(colliders->first == NULL)
3026                 return 0;
3027