=bmesh= merge from trunk at r36529
[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->getNumTessFaces(dm);
360                         totelem= me->totface;
361                         origindex= dm->getTessFaceDataArray(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->getNumTessFaces(dm);
531                 mface=mface_array=dm->getTessFaceDataArray(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) + (double)amount;
639         offs[1]= rng_getDouble(rng) + (double)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.0f/sqrtf((float)num));
663         rad2= (float)(1.0f/((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->getTessFaceData(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->getNumTessFaces(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->getTessFaceDataArray(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->getTessFaceData(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.0f;
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->getTessFaceDataArray(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->getNumTessFaces(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->getTessFaceData(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->getTessFaceData(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.0f;
1214                                 }
1215                                 else {
1216                                         tweight /= 3.0f;
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.0f)) {
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->getTessFaceDataArray(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.0f){
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((float)M_PI*(part->tanphase+phase)));
1613                 fac=-(float)sin((float)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.0f){
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(fabsf(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         if(part->type == PART_HAIR){
1818                 pa->lifetime = 100.0f;
1819         }
1820         else{
1821                 pa->lifetime = part->lifetime * ptex.life;
1822
1823                 if(part->randlife != 0.0f)
1824                         pa->lifetime *= 1.0f - part->randlife * PSYS_FRAND(p + 21);
1825         }
1826
1827         pa->dietime = pa->time + pa->lifetime;
1828
1829         if(sim->psys->pointcache && sim->psys->pointcache->flag & PTCACHE_BAKED &&
1830                 sim->psys->pointcache->mem_cache.first) {
1831                 float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1832                 pa->dietime = MIN2(pa->dietime, dietime);
1833         }
1834
1835         if(pa->time > cfra)
1836                 pa->alive = PARS_UNBORN;
1837         else if(pa->dietime <= cfra)
1838                 pa->alive = PARS_DEAD;
1839         else
1840                 pa->alive = PARS_ALIVE;
1841
1842         pa->state.time = cfra;
1843 }
1844 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1845 {
1846         ParticleData *pa;
1847         int p, totpart=sim->psys->totpart;
1848         
1849         for(p=from, pa=sim->psys->particles+from; p<totpart; p++, pa++)
1850                 reset_particle(sim, pa, dtime, cfra);
1851 }
1852 /************************************************/
1853 /*                      Particle targets                                        */
1854 /************************************************/
1855 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1856 {
1857         ParticleSystem *psys = NULL;
1858
1859         if(pt->ob == NULL || pt->ob == ob)
1860                 psys = BLI_findlink(&ob->particlesystem, pt->psys-1);
1861         else
1862                 psys = BLI_findlink(&pt->ob->particlesystem, pt->psys-1);
1863
1864         if(psys)
1865                 pt->flag |= PTARGET_VALID;
1866         else
1867                 pt->flag &= ~PTARGET_VALID;
1868
1869         return psys;
1870 }
1871 /************************************************/
1872 /*                      Keyed particles                                         */
1873 /************************************************/
1874 /* Counts valid keyed targets */
1875 void psys_count_keyed_targets(ParticleSimulationData *sim)
1876 {
1877         ParticleSystem *psys = sim->psys, *kpsys;
1878         ParticleTarget *pt = psys->targets.first;
1879         int keys_valid = 1;
1880         psys->totkeyed = 0;
1881
1882         for(; pt; pt=pt->next) {
1883                 kpsys = psys_get_target_system(sim->ob, pt);
1884
1885                 if(kpsys && kpsys->totpart) {
1886                         psys->totkeyed += keys_valid;
1887                         if(psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
1888                                 psys->totkeyed += 1;
1889                 }
1890                 else {
1891                         keys_valid = 0;
1892                 }
1893         }
1894
1895         psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1896 }
1897
1898 static void set_keyed_keys(ParticleSimulationData *sim)
1899 {
1900         ParticleSystem *psys = sim->psys;
1901         ParticleSimulationData ksim= {0};
1902         ParticleTarget *pt;
1903         PARTICLE_P;
1904         ParticleKey *key;
1905         int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1906
1907         ksim.scene= sim->scene;
1908         
1909         /* no proper targets so let's clear and bail out */
1910         if(psys->totkeyed==0) {
1911                 free_keyed_keys(psys);
1912                 psys->flag &= ~PSYS_KEYED;
1913                 return;
1914         }
1915
1916         if(totpart && psys->particles->totkey != totkeys) {
1917                 free_keyed_keys(psys);
1918                 
1919                 key = MEM_callocN(totpart*totkeys*sizeof(ParticleKey), "Keyed keys");
1920                 
1921                 LOOP_PARTICLES {
1922                         pa->keys = key;
1923                         pa->totkey = totkeys;
1924                         key += totkeys;
1925                 }
1926         }
1927         
1928         psys->flag &= ~PSYS_KEYED;
1929
1930
1931         pt = psys->targets.first;
1932         for(k=0; k<totkeys; k++) {
1933                 ksim.ob = pt->ob ? pt->ob : sim->ob;
1934                 ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
1935
1936                 LOOP_PARTICLES {
1937                         key = pa->keys + k;
1938                         key->time = -1.0; /* use current time */
1939
1940                         psys_get_particle_state(&ksim, p%ksim.psys->totpart, key, 1);
1941
1942                         if(psys->flag & PSYS_KEYED_TIMING){
1943                                 key->time = pa->time + pt->time;
1944                                 if(pt->duration != 0.0f && k+1 < totkeys) {
1945                                         copy_particle_key(key+1, key, 1);
1946                                         (key+1)->time = pa->time + pt->time + pt->duration;
1947                                 }
1948                         }
1949                         else if(totkeys > 1)
1950                                 key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
1951                         else
1952                                 key->time = pa->time;
1953                 }
1954
1955                 if(psys->flag & PSYS_KEYED_TIMING && pt->duration!=0.0f)
1956                         k++;
1957
1958                 pt = (pt->next && pt->next->flag & PTARGET_VALID)? pt->next : psys->targets.first;
1959         }
1960
1961         psys->flag |= PSYS_KEYED;
1962 }
1963
1964 /************************************************/
1965 /*                      Point Cache                                                     */
1966 /************************************************/
1967 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
1968 {
1969         PointCache *cache = psys->pointcache;
1970
1971         if(cache->flag & PTCACHE_DISK_CACHE && cache->mem_cache.first == NULL) {
1972                 PTCacheID pid;
1973                 BKE_ptcache_id_from_particles(&pid, ob, psys);
1974                 cache->flag &= ~PTCACHE_DISK_CACHE;
1975                 BKE_ptcache_disk_to_mem(&pid);
1976                 cache->flag |= PTCACHE_DISK_CACHE;
1977         }
1978 }
1979 static void psys_clear_temp_pointcache(ParticleSystem *psys)
1980 {
1981         if(psys->pointcache->flag & PTCACHE_DISK_CACHE)
1982                 BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
1983 }
1984 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
1985 {
1986         ParticleSettings *part = psys->part;
1987
1988         *sfra = MAX2(1, (int)part->sta);
1989         *efra = MIN2((int)(part->end + part->lifetime + 1.0f), scene->r.efra);
1990 }
1991
1992 /************************************************/
1993 /*                      Effectors                                                       */
1994 /************************************************/
1995 static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
1996 {
1997         if(psys) {
1998                 PARTICLE_P;
1999                 int totpart = 0;
2000
2001                 if(!psys->bvhtree || psys->bvhtree_frame != cfra) {
2002                         LOOP_SHOWN_PARTICLES {
2003                                 totpart++;
2004                         }
2005                         
2006                         BLI_bvhtree_free(psys->bvhtree);
2007                         psys->bvhtree = BLI_bvhtree_new(totpart, 0.0, 4, 6);
2008
2009                         LOOP_SHOWN_PARTICLES {
2010                                 if(pa->alive == PARS_ALIVE) {
2011                                         if(pa->state.time == cfra)
2012                                                 BLI_bvhtree_insert(psys->bvhtree, p, pa->prev_state.co, 1);
2013                                         else
2014                                                 BLI_bvhtree_insert(psys->bvhtree, p, pa->state.co, 1);
2015                                 }
2016                         }
2017                         BLI_bvhtree_balance(psys->bvhtree);
2018
2019                         psys->bvhtree_frame = cfra;
2020                 }
2021         }
2022 }
2023 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
2024 {
2025         if(psys) {
2026                 PARTICLE_P;
2027                 int totpart = 0;
2028
2029                 if(!psys->tree || psys->tree_frame != cfra) {
2030                         LOOP_SHOWN_PARTICLES {
2031                                 totpart++;
2032                         }
2033
2034                         BLI_kdtree_free(psys->tree);
2035                         psys->tree = BLI_kdtree_new(psys->totpart);
2036
2037                         LOOP_SHOWN_PARTICLES {
2038                                 if(pa->alive == PARS_ALIVE) {
2039                                         if(pa->state.time == cfra)
2040                                                 BLI_kdtree_insert(psys->tree, p, pa->prev_state.co, NULL);
2041                                         else
2042                                                 BLI_kdtree_insert(psys->tree, p, pa->state.co, NULL);
2043                                 }
2044                         }
2045                         BLI_kdtree_balance(psys->tree);
2046
2047                         psys->tree_frame = cfra;
2048                 }
2049         }
2050 }
2051
2052 static void psys_update_effectors(ParticleSimulationData *sim)
2053 {
2054         pdEndEffectors(&sim->psys->effectors);
2055         sim->psys->effectors = pdInitEffectors(sim->scene, sim->ob, sim->psys, sim->psys->part->effector_weights);
2056         precalc_guides(sim, sim->psys->effectors);
2057 }
2058
2059 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)
2060 {
2061         ParticleKey states[5];
2062         float force[3],acceleration[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
2063         float pa_mass= (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
2064         int i, steps=1;
2065         int integrator = part->integrator;
2066
2067         copy_v3_v3(oldpos, pa->state.co);
2068         
2069         /* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
2070         if(pa->prev_state.time < 0.f && integrator == PART_INT_VERLET)
2071                 integrator = PART_INT_EULER;
2072
2073         switch(integrator){
2074                 case PART_INT_EULER:
2075                         steps=1;
2076                         break;
2077                 case PART_INT_MIDPOINT:
2078                         steps=2;
2079                         break;
2080                 case PART_INT_RK4:
2081                         steps=4;
2082                         break;
2083                 case PART_INT_VERLET:
2084                         steps=1;
2085                         break;
2086         }
2087
2088         copy_particle_key(states, &pa->state, 1);
2089
2090         states->time = 0.f;
2091
2092         for(i=0; i<steps; i++){
2093                 zero_v3(force);
2094                 zero_v3(impulse);
2095
2096                 force_func(forcedata, states+i, force, impulse);
2097
2098                 /* force to acceleration*/
2099                 mul_v3_v3fl(acceleration, force, 1.0f/pa_mass);
2100
2101                 if(external_acceleration)
2102                         add_v3_v3(acceleration, external_acceleration);
2103                 
2104                 /* calculate next state */
2105                 add_v3_v3(states[i].vel, impulse);
2106
2107                 switch(integrator){
2108                         case PART_INT_EULER:
2109                                 madd_v3_v3v3fl(pa->state.co, states->co, states->vel, dtime);
2110                                 madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2111                                 break;
2112                         case PART_INT_MIDPOINT:
2113                                 if(i==0){
2114                                         madd_v3_v3v3fl(states[1].co, states->co, states->vel, dtime*0.5f);
2115                                         madd_v3_v3v3fl(states[1].vel, states->vel, acceleration, dtime*0.5f);
2116                                         states[1].time = dtime*0.5f;
2117                                         /*fra=sim->psys->cfra+0.5f*dfra;*/
2118                                 }
2119                                 else{
2120                                         madd_v3_v3v3fl(pa->state.co, states->co, states[1].vel, dtime);
2121                                         madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2122                                 }
2123                                 break;
2124                         case PART_INT_RK4:
2125                                 switch(i){
2126                                         case 0:
2127                                                 copy_v3_v3(dx[0], states->vel);
2128                                                 mul_v3_fl(dx[0], dtime);
2129                                                 copy_v3_v3(dv[0], acceleration);
2130                                                 mul_v3_fl(dv[0], dtime);
2131
2132                                                 madd_v3_v3v3fl(states[1].co, states->co, dx[0], 0.5f);
2133                                                 madd_v3_v3v3fl(states[1].vel, states->vel, dv[0], 0.5f);
2134                                                 states[1].time = dtime*0.5f;
2135                                                 /*fra=sim->psys->cfra+0.5f*dfra;*/
2136                                                 break;
2137                                         case 1:
2138                                                 madd_v3_v3v3fl(dx[1], states->vel, dv[0], 0.5f);
2139                                                 mul_v3_fl(dx[1], dtime);
2140                                                 copy_v3_v3(dv[1], acceleration);
2141                                                 mul_v3_fl(dv[1], dtime);
2142
2143                                                 madd_v3_v3v3fl(states[2].co, states->co, dx[1], 0.5f);
2144                                                 madd_v3_v3v3fl(states[2].vel, states->vel, dv[1], 0.5f);
2145                                                 states[2].time = dtime*0.5f;
2146                                                 break;
2147                                         case 2:
2148                                                 madd_v3_v3v3fl(dx[2], states->vel, dv[1], 0.5f);
2149                                                 mul_v3_fl(dx[2], dtime);
2150                                                 copy_v3_v3(dv[2], acceleration);
2151                                                 mul_v3_fl(dv[2], dtime);
2152
2153                                                 add_v3_v3v3(states[3].co, states->co, dx[2]);
2154                                                 add_v3_v3v3(states[3].vel, states->vel, dv[2]);
2155                                                 states[3].time = dtime;
2156                                                 /*fra=cfra;*/
2157                                                 break;
2158                                         case 3:
2159                                                 add_v3_v3v3(dx[3], states->vel, dv[2]);
2160                                                 mul_v3_fl(dx[3], dtime);
2161                                                 copy_v3_v3(dv[3], acceleration);
2162                                                 mul_v3_fl(dv[3], dtime);
2163
2164                                                 madd_v3_v3v3fl(pa->state.co, states->co, dx[0], 1.0f/6.0f);
2165                                                 madd_v3_v3fl(pa->state.co, dx[1], 1.0f/3.0f);
2166                                                 madd_v3_v3fl(pa->state.co, dx[2], 1.0f/3.0f);
2167                                                 madd_v3_v3fl(pa->state.co, dx[3], 1.0f/6.0f);
2168
2169                                                 madd_v3_v3v3fl(pa->state.vel, states->vel, dv[0], 1.0f/6.0f);
2170                                                 madd_v3_v3fl(pa->state.vel, dv[1], 1.0f/3.0f);
2171                                                 madd_v3_v3fl(pa->state.vel, dv[2], 1.0f/3.0f);
2172                                                 madd_v3_v3fl(pa->state.vel, dv[3], 1.0f/6.0f);
2173                                 }
2174                                 break;
2175                         case PART_INT_VERLET:   /* Verlet integration */
2176                                 madd_v3_v3v3fl(pa->state.vel, pa->prev_state.vel, acceleration, dtime);
2177                                 madd_v3_v3v3fl(pa->state.co, pa->prev_state.co, pa->state.vel, dtime);
2178
2179                                 sub_v3_v3v3(pa->state.vel, pa->state.co, oldpos);
2180                                 mul_v3_fl(pa->state.vel, 1.0f/dtime);
2181                                 break;
2182                 }
2183         }
2184 }
2185
2186 /*********************************************************************************************************
2187                     SPH fluid physics 
2188
2189  In theory, there could be unlimited implementation of SPH simulators
2190
2191  This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
2192
2193  Titled: Particle-based Viscoelastic Fluid Simulation.
2194  Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
2195  Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
2196
2197  Presented at Siggraph, (2005)
2198
2199 ***********************************************************************************************************/
2200 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
2201 static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
2202 {
2203         /* Are more refs required? */
2204         if(psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
2205                 psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
2206                 psys->fluid_springs = (ParticleSpring*)MEM_callocN(psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
2207         }
2208         else if(psys->tot_fluidsprings == psys->alloc_fluidsprings) {
2209                 /* Double the number of refs allocated */
2210                 psys->alloc_fluidsprings *= 2;
2211                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
2212         }
2213
2214         memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
2215         psys->tot_fluidsprings++;
2216
2217         return psys->fluid_springs + psys->tot_fluidsprings - 1;
2218 }
2219 static void sph_spring_delete(ParticleSystem *psys, int j)
2220 {
2221         if (j != psys->tot_fluidsprings - 1)
2222                 psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
2223
2224         psys->tot_fluidsprings--;
2225
2226         if (psys->tot_fluidsprings < psys->alloc_fluidsprings/2 && psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE){
2227                 psys->alloc_fluidsprings /= 2;
2228                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs,  psys->alloc_fluidsprings * sizeof(ParticleSpring));
2229         }
2230 }
2231 static void sph_springs_modify(ParticleSystem *psys, float dtime){
2232         SPHFluidSettings *fluid = psys->part->fluid;
2233         ParticleData *pa1, *pa2;
2234         ParticleSpring *spring = psys->fluid_springs;
2235         
2236         float h, d, Rij[3], rij, Lij;
2237         int i;
2238
2239         float yield_ratio = fluid->yield_ratio;
2240         float plasticity = fluid->plasticity_constant;
2241         /* scale things according to dtime */
2242         float timefix = 25.f * dtime;
2243
2244         if((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
2245                 return;
2246
2247         /* Loop through the springs */
2248         for(i=0; i<psys->tot_fluidsprings; i++, spring++) {
2249                 pa1 = psys->particles + spring->particle_index[0];
2250                 pa2 = psys->particles + spring->particle_index[1];
2251
2252                 sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
2253                 rij = normalize_v3(Rij);
2254
2255                 /* adjust rest length */
2256                 Lij = spring->rest_length;
2257                 d = yield_ratio * timefix * Lij;
2258
2259                 if (rij > Lij + d) // Stretch
2260                         spring->rest_length += plasticity * (rij - Lij - d) * timefix;
2261                 else if(rij < Lij - d) // Compress
2262                         spring->rest_length -= plasticity * (Lij - d - rij) * timefix;
2263
2264                 h = 4.f*pa1->size;
2265
2266                 if(spring->rest_length > h)
2267                         spring->delete_flag = 1;
2268         }
2269
2270         /* Loop through springs backwaqrds - for efficient delete function */
2271         for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
2272                 if(psys->fluid_springs[i].delete_flag)
2273                         sph_spring_delete(psys, i);
2274         }
2275 }
2276 static EdgeHash *sph_springhash_build(ParticleSystem *psys)
2277 {
2278         EdgeHash *springhash = NULL;
2279         ParticleSpring *spring;
2280         int i = 0;
2281
2282         springhash = BLI_edgehash_new();
2283
2284         for(i=0, spring=psys->fluid_springs; i<psys->tot_fluidsprings; i++, spring++)
2285                 BLI_edgehash_insert(springhash, spring->particle_index[0], spring->particle_index[1], SET_INT_IN_POINTER(i+1));
2286
2287         return springhash;
2288 }
2289
2290 typedef struct SPHNeighbor
2291 {
2292         ParticleSystem *psys;
2293         int index;
2294 } SPHNeighbor;
2295 typedef struct SPHRangeData
2296 {
2297         SPHNeighbor neighbors[128];
2298         int tot_neighbors;
2299
2300         float density, near_density;
2301         float h;
2302
2303         ParticleSystem *npsys;
2304         ParticleData *pa;
2305
2306         float massfac;
2307         int use_size;
2308 } SPHRangeData;
2309 typedef struct SPHData {
2310         ParticleSystem *psys[10];
2311         ParticleData *pa;
2312         float mass;
2313         EdgeHash *eh;
2314         float *gravity;
2315 }SPHData;
2316 static void sph_density_accum_cb(void *userdata, int index, float squared_dist)
2317 {
2318         SPHRangeData *pfr = (SPHRangeData *)userdata;
2319         ParticleData *npa = pfr->npsys->particles + index;
2320         float q;
2321
2322         if(npa == pfr->pa || squared_dist < FLT_EPSILON)
2323                 return;
2324
2325         /* Ugh! One particle has over 128 neighbors! Really shouldn't happen,
2326          * but even if it does it shouldn't do any terrible harm if all are
2327          * not taken into account - jahka
2328          */
2329         if(pfr->tot_neighbors >= 128)
2330                 return;
2331         
2332         pfr->neighbors[pfr->tot_neighbors].index = index;
2333         pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
2334         pfr->tot_neighbors++;
2335
2336         q = (1.f - sqrtf(squared_dist)/pfr->h) * pfr->massfac;
2337
2338         if(pfr->use_size)
2339                 q *= npa->size;
2340
2341         pfr->density += q*q;
2342         pfr->near_density += q*q*q;
2343 }
2344 static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
2345 {
2346         SPHData *sphdata = (SPHData *)sphdata_v;
2347         ParticleSystem **psys = sphdata->psys;
2348         ParticleData *pa = sphdata->pa;
2349         SPHFluidSettings *fluid = psys[0]->part->fluid;
2350         ParticleSpring *spring = NULL;
2351         SPHRangeData pfr;
2352         SPHNeighbor *pfn;
2353         float mass = sphdata->mass;
2354         float *gravity = sphdata->gravity;
2355         EdgeHash *springhash = sphdata->eh;
2356
2357         float q, u, rij, dv[3];
2358         float pressure, near_pressure;
2359
2360         float visc = fluid->viscosity_omega;
2361         float stiff_visc = fluid->viscosity_beta * (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
2362
2363         float inv_mass = 1.0f/mass;
2364         float spring_constant = fluid->spring_k;
2365         
2366         float h = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.f*pa->size : 1.f); /* 4.0 seems to be a pretty good value */
2367         float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f); /* 4.77 is an experimentally determined density factor */
2368         float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
2369
2370         float stiffness = fluid->stiffness_k;
2371         float stiffness_near_fac = fluid->stiffness_knear * (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);
2372
2373         ParticleData *npa;
2374         float vec[3];
2375         float vel[3];
2376         float co[3];
2377
2378         int i, spring_index, index = pa - psys[0]->particles;
2379
2380         pfr.tot_neighbors = 0;
2381         pfr.density = pfr.near_density = 0.f;
2382         pfr.h = h;
2383         pfr.pa = pa;
2384
2385         for(i=0; i<10 && psys[i]; i++) {
2386                 pfr.npsys = psys[i];
2387                 pfr.massfac = psys[i]->part->mass*inv_mass;
2388                 pfr.use_size = psys[i]->part->flag & PART_SIZEMASS;
2389
2390                 BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sph_density_accum_cb, &pfr);
2391         }
2392
2393         pressure =  stiffness * (pfr.density - rest_density);
2394         near_pressure = stiffness_near_fac * pfr.near_density;
2395
2396         pfn = pfr.neighbors;
2397         for(i=0; i<pfr.tot_neighbors; i++, pfn++) {
2398                 npa = pfn->psys->particles + pfn->index;
2399
2400                 madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
2401
2402                 sub_v3_v3v3(vec, co, state->co);
2403                 rij = normalize_v3(vec);
2404
2405                 q = (1.f - rij/h) * pfn->psys->part->mass * inv_mass;
2406
2407                 if(pfn->psys->part->flag & PART_SIZEMASS)
2408                         q *= npa->size;
2409
2410                 copy_v3_v3(vel, npa->prev_state.vel);
2411
2412                 /* Double Density Relaxation */
2413                 madd_v3_v3fl(force, vec, -(pressure + near_pressure*q)*q);
2414
2415                 /* Viscosity */
2416                 if(visc > 0.f   || stiff_visc > 0.f) {          
2417                         sub_v3_v3v3(dv, vel, state->vel);
2418                         u = dot_v3v3(vec, dv);
2419
2420                         if(u < 0.f && visc > 0.f)
2421                                 madd_v3_v3fl(force, vec, 0.5f * q * visc * u );
2422
2423                         if(u > 0.f && stiff_visc > 0.f)
2424                                 madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u );
2425                 }
2426
2427                 if(spring_constant > 0.f) {
2428                         /* Viscoelastic spring force */
2429                         if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
2430                                 spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, pfn->index));
2431
2432                                 if(spring_index) {
2433                                         spring = psys[0]->fluid_springs + spring_index - 1;
2434
2435                                         madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (spring->rest_length - rij));
2436                                 }
2437                                 else if(fluid->spring_frames == 0 || (pa->prev_state.time-pa->time) <= fluid->spring_frames){
2438                                         ParticleSpring temp_spring;
2439                                         temp_spring.particle_index[0] = index;
2440                                         temp_spring.particle_index[1] = pfn->index;
2441                                         temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length;
2442                                         temp_spring.delete_flag = 0;
2443                                                                 
2444                                         sph_spring_add(psys[0], &temp_spring);
2445                                 }
2446                         }
2447                         else {/* PART_SPRING_HOOKES - Hooke's spring force */
2448                                 madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (rest_length - rij));
2449                         }
2450                 }
2451         }
2452         
2453         /* Artificial buoyancy force in negative gravity direction  */
2454         if (fluid->buoyancy > 0.f && gravity)
2455                 madd_v3_v3fl(force, gravity, fluid->buoyancy * (pfr.density-rest_density));
2456 }
2457
2458 static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, float *gravity, EdgeHash *springhash){
2459         ParticleTarget *pt;
2460         int i;
2461
2462         ParticleSettings *part = sim->psys->part;
2463         // float timestep = psys_get_timestep(sim); // UNUSED
2464         float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
2465         float dtime = dfra*psys_get_timestep(sim);
2466         // int steps = 1; // UNUSED
2467         float effector_acceleration[3];
2468         SPHData sphdata;
2469
2470         sphdata.psys[0] = sim->psys;
2471         for(i=1, pt=sim->psys->targets.first; i<10; i++, pt=(pt?pt->next:NULL))
2472                 sphdata.psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL;
2473
2474         sphdata.pa = pa;
2475         sphdata.gravity = gravity;
2476         sphdata.mass = pa_mass;
2477         sphdata.eh = springhash;
2478
2479         /* restore previous state and treat gravity & effectors as external acceleration*/
2480         sub_v3_v3v3(effector_acceleration, pa->state.vel, pa->prev_state.vel);
2481         mul_v3_fl(effector_acceleration, 1.f/dtime);
2482
2483         copy_particle_key(&pa->state, &pa->prev_state, 0);
2484
2485         integrate_particle(part, pa, dtime, effector_acceleration, sph_force_cb, &sphdata);
2486 }
2487
2488 /************************************************/
2489 /*                      Basic physics                                           */
2490 /************************************************/
2491 typedef struct EfData
2492 {
2493         ParticleTexture ptex;
2494         ParticleSimulationData *sim;
2495         ParticleData *pa;
2496 } EfData;
2497 static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, float *impulse)
2498 {
2499         EfData *efdata = (EfData *)efdata_v;
2500         ParticleSimulationData *sim = efdata->sim;
2501         ParticleSettings *part = sim->psys->part;
2502         ParticleData *pa = efdata->pa;
2503         EffectedPoint epoint;
2504
2505         /* add effectors */
2506         pd_point_from_particle(efdata->sim, efdata->pa, state, &epoint);
2507         if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
2508                 pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
2509
2510         mul_v3_fl(force, efdata->ptex.field);
2511         mul_v3_fl(impulse, efdata->ptex.field);
2512
2513         /* calculate air-particle interaction */
2514         if(part->dragfac != 0.0f)
2515                 madd_v3_v3fl(force, state->vel, -part->dragfac * pa->size * pa->size * len_v3(state->vel));
2516
2517         /* brownian force */
2518         if(part->brownfac != 0.0f){
2519                 force[0] += (BLI_frand()-0.5f) * part->brownfac;
2520                 force[1] += (BLI_frand()-0.5f) * part->brownfac;
2521                 force[2] += (BLI_frand()-0.5f) * part->brownfac;
2522         }
2523 }
2524 /* gathers all forces that effect particles and calculates a new state for the particle */
2525 static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
2526 {
2527         ParticleSettings *part = sim->psys->part;
2528         ParticleData *pa = sim->psys->particles + p;
2529         ParticleKey tkey;
2530         float dtime=dfra*psys_get_timestep(sim), time;
2531         float *gravity = NULL, gr[3];
2532         EfData efdata;
2533
2534         psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);
2535
2536         efdata.pa = pa;
2537         efdata.sim = sim;
2538
2539         /* add global acceleration (gravitation) */
2540         if(psys_uses_gravity(sim)
2541                 /* normal gravity is too strong for hair so it's disabled by default */
2542                 && (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
2543                 zero_v3(gr);
2544                 madd_v3_v3fl(gr, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * efdata.ptex.gravity);
2545                 gravity = gr;
2546         }
2547
2548         /* maintain angular velocity */
2549         copy_v3_v3(pa->state.ave, pa->prev_state.ave);
2550
2551         integrate_particle(part, pa, dtime, gravity, basic_force_cb, &efdata);
2552
2553         /* damp affects final velocity */
2554         if(part->dampfac != 0.f)
2555                 mul_v3_fl(pa->state.vel, 1.f - part->dampfac * efdata.ptex.damp * 25.f * dtime);
2556
2557         //VECCOPY(pa->state.ave, states->ave);
2558
2559         /* finally we do guides */
2560         time=(cfra-pa->time)/pa->lifetime;
2561         CLAMP(time, 0.0f, 1.0f);
2562
2563         VECCOPY(tkey.co,pa->state.co);
2564         VECCOPY(tkey.vel,pa->state.vel);
2565         tkey.time=pa->state.time;
2566
2567         if(part->type != PART_HAIR) {
2568                 if(do_guides(sim->psys->effectors, &tkey, p, time)) {
2569                         VECCOPY(pa->state.co,tkey.co);
2570                         /* guides don't produce valid velocity */
2571                         VECSUB(pa->state.vel,tkey.co,pa->prev_state.co);
2572                         mul_v3_fl(pa->state.vel,1.0f/dtime);
2573                         pa->state.time=tkey.time;
2574                 }
2575         }
2576 }
2577 static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2578 {
2579         float rotfac, rot1[4], rot2[4]={1.0,0.0,0.0,0.0}, dtime=dfra*timestep;
2580
2581         if((part->flag & PART_ROT_DYN)==0){
2582                 if(part->avemode==PART_AVE_SPIN){
2583                         float angle;
2584                         float len1 = len_v3(pa->prev_state.vel);
2585                         float len2 = len_v3(pa->state.vel);
2586
2587                         if(len1==0.0f || len2==0.0f)
2588                                 pa->state.ave[0]=pa->state.ave[1]=pa->state.ave[2]=0.0f;
2589                         else{
2590                                 cross_v3_v3v3(pa->state.ave,pa->prev_state.vel,pa->state.vel);
2591                                 normalize_v3(pa->state.ave);
2592                                 angle=dot_v3v3(pa->prev_state.vel,pa->state.vel)/(len1*len2);
2593                                 mul_v3_fl(pa->state.ave,saacos(angle)/dtime);
2594                         }
2595
2596                         axis_angle_to_quat(rot2,pa->state.vel,dtime*part->avefac);
2597                 }
2598         }
2599
2600         rotfac=len_v3(pa->state.ave);
2601         if(rotfac == 0.0f){ /* unit_qt(in VecRotToQuat) doesn't give unit quat [1,0,0,0]?? */
2602                 rot1[0]=1.0f;
2603                 rot1[1]=rot1[2]=rot1[3]=0;
2604         }
2605         else{
2606                 axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
2607         }
2608         mul_qt_qtqt(pa->state.rot,rot1,pa->prev_state.rot);
2609         mul_qt_qtqt(pa->state.rot,rot2,pa->state.rot);
2610
2611         /* keep rotation quat in good health */
2612         normalize_qt(pa->state.rot);
2613 }
2614
2615 /************************************************/
2616 /*                      Collisions                                                      */
2617 /************************************************/
2618 #define COLLISION_MAX_COLLISIONS        10
2619 #define COLLISION_MIN_RADIUS 0.001f
2620 #define COLLISION_MIN_DISTANCE 0.0001f
2621 #define COLLISION_ZERO 0.00001f
2622 typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
2623 static float nr_signed_distance_to_plane(float *p, float radius, ParticleCollisionElement *pce, float *nor)
2624 {
2625         float p0[3], e1[3], e2[3], d;
2626
2627         sub_v3_v3v3(e1, pce->x1, pce->x0);
2628         sub_v3_v3v3(e2, pce->x2, pce->x0);
2629         sub_v3_v3v3(p0, p, pce->x0);
2630
2631         cross_v3_v3v3(nor, e1, e2);
2632         normalize_v3(nor);
2633
2634         d = dot_v3v3(p0, nor);
2635
2636         if(pce->inv_nor == -1) {
2637                 if(d < 0.f)
2638                         pce->inv_nor = 1;
2639                 else
2640                         pce->inv_nor = 0;
2641         }
2642
2643         if(pce->inv_nor == 1) {
2644                 mul_v3_fl(nor, -1.f);
2645                 d = -d;
2646         }
2647
2648         return d - radius;
2649 }
2650 static float nr_distance_to_edge(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(nor))
2651 {
2652         float v0[3], v1[3], v2[3], c[3];
2653
2654         sub_v3_v3v3(v0, pce->x1, pce->x0);
2655         sub_v3_v3v3(v1, p, pce->x0);
2656         sub_v3_v3v3(v2, p, pce->x1);
2657
2658         cross_v3_v3v3(c, v1, v2);
2659
2660         return fabsf(len_v3(c)/len_v3(v0)) - radius;
2661 }
2662 static float nr_distance_to_vert(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(nor))
2663 {
2664         return len_v3v3(p, pce->x0) - radius;
2665 }
2666 static void collision_interpolate_element(ParticleCollisionElement *pce, float t, float fac, ParticleCollision *col)
2667 {
2668         /* t is the current time for newton rhapson */
2669         /* fac is the starting factor for current collision iteration */
2670         /* the col->fac's are factors for the particle subframe step start and end during collision modifier step */
2671         float f = fac + t*(1.f-fac);
2672         float mul = col->fac1 + f * (col->fac2-col->fac1);
2673         if(pce->tot > 0) {
2674                 madd_v3_v3v3fl(pce->x0, pce->x[0], pce->v[0], mul);
2675
2676                 if(pce->tot > 1) {
2677                         madd_v3_v3v3fl(pce->x1, pce->x[1], pce->v[1], mul);
2678
2679                         if(pce->tot > 2)
2680                                 madd_v3_v3v3fl(pce->x2, pce->x[2], pce->v[2], mul);
2681                 }
2682         }
2683 }
2684 static void collision_point_velocity(ParticleCollisionElement *pce)
2685 {
2686         float v[3];
2687
2688         copy_v3_v3(pce->vel, pce->v[0]);
2689
2690         if(pce->tot > 1) {
2691                 sub_v3_v3v3(v, pce->v[1], pce->v[0]);
2692                 madd_v3_v3fl(pce->vel, v, pce->uv[0]);
2693
2694                 if(pce->tot > 2) {
2695                         sub_v3_v3v3(v, pce->v[2], pce->v[0]);
2696                         madd_v3_v3fl(pce->vel, v, pce->uv[1]);
2697                 }
2698         }
2699 }
2700 static float collision_point_distance_with_normal(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *nor)
2701 {
2702         if(fac >= 0.f)
2703                 collision_interpolate_element(pce, 0.f, fac, col);
2704
2705         switch(pce->tot) {
2706                 case 1:
2707                 {
2708                         sub_v3_v3v3(nor, p, pce->x0);
2709                         return normalize_v3(nor);
2710                 }
2711                 case 2:
2712                 {
2713                         float u, e[3], vec[3];
2714                         sub_v3_v3v3(e, pce->x1, pce->x0);
2715                         sub_v3_v3v3(vec, p, pce->x0);
2716                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2717
2718                         madd_v3_v3v3fl(nor, vec, e, -u);
2719                         return normalize_v3(nor);
2720                 }
2721                 case 3:
2722                         return nr_signed_distance_to_plane(p, 0.f, pce, nor);
2723         }
2724         return 0;
2725 }
2726 static void collision_point_on_surface(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *co)
2727 {
2728         collision_interpolate_element(pce, 0.f, fac, col);
2729
2730         switch(pce->tot) {
2731                 case 1:
2732                 {
2733                         sub_v3_v3v3(co, p, pce->x0);
2734                         normalize_v3(co);
2735                         madd_v3_v3v3fl(co, pce->x0, co, col->radius);
2736                         break;
2737                 }
2738                 case 2:
2739                 {
2740                         float u, e[3], vec[3], nor[3];
2741                         sub_v3_v3v3(e, pce->x1, pce->x0);
2742                         sub_v3_v3v3(vec, p, pce->x0);
2743                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2744
2745                         madd_v3_v3v3fl(nor, vec, e, -u);
2746                         normalize_v3(nor);
2747
2748                         madd_v3_v3v3fl(co, pce->x0, e, pce->uv[0]);
2749                         madd_v3_v3fl(co, nor, col->radius);
2750                         break;
2751                 }
2752                 case 3:
2753                 {
2754                                 float p0[3], e1[3], e2[3], nor[3];
2755
2756                                 sub_v3_v3v3(e1, pce->x1, pce->x0);
2757                                 sub_v3_v3v3(e2, pce->x2, pce->x0);
2758                                 sub_v3_v3v3(p0, p, pce->x0);
2759
2760                                 cross_v3_v3v3(nor, e1, e2);
2761                                 normalize_v3(nor);
2762
2763                                 if(pce->inv_nor == 1)
2764                                         mul_v3_fl(nor, -1.f);
2765
2766                                 madd_v3_v3v3fl(co, pce->x0, nor, col->radius);
2767                                 madd_v3_v3fl(co, e1, pce->uv[0]);
2768                                 madd_v3_v3fl(co, e2, pce->uv[1]);
2769                         break;
2770                 }
2771         }
2772 }
2773 /* find first root in range [0-1] starting from 0 */
2774 static float collision_newton_rhapson(ParticleCollision *col, float radius, ParticleCollisionElement *pce, NRDistanceFunc distance_func)
2775 {
2776         float t0, t1, d0, d1, dd, n[3];
2777         int iter;
2778
2779         pce->inv_nor = -1;
2780
2781         /* start from the beginning */
2782         t0 = 0.f;
2783         collision_interpolate_element(pce, t0, col->f, col);
2784         d0 = distance_func(col->co1, radius, pce, n);
2785         t1 = 0.001f;
2786         d1 = 0.f;
2787
2788         for(iter=0; iter<10; iter++) {//, itersum++) {
2789                 /* get current location */
2790                 collision_interpolate_element(pce, t1, col->f, col);
2791                 interp_v3_v3v3(pce->p, col->co1, col->co2, t1);
2792
2793                 d1 = distance_func(pce->p, radius, pce, n);
2794
2795                 /* no movement, so no collision */
2796                 if(d1 == d0) {
2797                         return -1.f;
2798                 }
2799
2800                 /* particle already inside face, so report collision */
2801                 if(iter == 0 && d0 < 0.f && d0 > -radius) {
2802                         copy_v3_v3(pce->p, col->co1);
2803                         copy_v3_v3(pce->nor, n);
2804                         pce->inside = 1;
2805                         return 0.f;
2806                 }
2807                 
2808                 dd = (t1-t0)/(d1-d0);
2809
2810                 t0 = t1;
2811                 d0 = d1;
2812
2813                 t1 -= d1*dd;
2814
2815                 /* particle movin away from plane could also mean a strangely rotating face, so check from end */
2816                 if(iter == 0 && t1 < 0.f) {
2817                         t0 = 1.f;
2818                         collision_interpolate_element(pce, t0, col->f, col);
2819                         d0 = distance_func(col->co2, radius, pce, n);
2820                         t1 = 0.999f;
2821                         d1 = 0.f;
2822
2823                         continue;
2824                 }
2825                 else if(iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f))
2826                         return -1.f;
2827
2828                 if(d1 <= COLLISION_ZERO && d1 >= -COLLISION_ZERO) {
2829                         if(t1 >= -COLLISION_ZERO && t1 <= 1.f) {
2830                                 if(distance_func == nr_signed_distance_to_plane)
2831                                         copy_v3_v3(pce->nor, n);
2832
2833                                 CLAMP(t1, 0.f, 1.f);
2834
2835                                 return t1;
2836                         }
2837                         else
2838                                 return -1.f;
2839                 }
2840         }
2841         return -1.0;
2842 }
2843 static int collision_sphere_to_tri(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2844 {
2845         ParticleCollisionElement *result = &col->pce;
2846         float ct, u, v;
2847
2848         pce->inv_nor = -1;
2849         pce->inside = 0;
2850
2851         ct = collision_newton_rhapson(col, radius, pce, nr_signed_distance_to_plane);
2852
2853         if(ct >= 0.f && ct < *t && (result->inside==0 || pce->inside==1) ) {
2854                 float e1[3], e2[3], p0[3];
2855                 float e1e1, e1e2, e1p0, e2e2, e2p0, inv;
2856
2857                 sub_v3_v3v3(e1, pce->x1, pce->x0);
2858                 sub_v3_v3v3(e2, pce->x2, pce->x0);
2859                 /* XXX: add radius correction here? */
2860                 sub_v3_v3v3(p0, pce->p, pce->x0);
2861
2862                 e1e1 = dot_v3v3(e1, e1);
2863                 e1e2 = dot_v3v3(e1, e2);
2864                 e1p0 = dot_v3v3(e1, p0);
2865                 e2e2 = dot_v3v3(e2, e2);
2866                 e2p0 = dot_v3v3(e2, p0);
2867
2868                 inv = 1.f/(e1e1 * e2e2 - e1e2 * e1e2);
2869                 u = (e2e2 * e1p0 - e1e2 * e2p0) * inv;
2870                 v = (e1e1 * e2p0 - e1e2 * e1p0) * inv;
2871
2872                 if(u>=0.f && u<=1.f && v>=0.f && u+v<=1.f) {
2873                         *result = *pce;
2874
2875                         /* normal already calculated in pce */
2876
2877                         result->uv[0] = u;
2878                         result->uv[1] = v;
2879
2880                         *t = ct;
2881                         return 1;
2882                 }
2883         }
2884         return 0;
2885 }
2886 static int collision_sphere_to_edges(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2887 {
2888         ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
2889         ParticleCollisionElement *result = &col->pce;
2890
2891         float ct;
2892         int i;
2893
2894         for(i=0; i<3; i++) {
2895                 /* in case of a quad, no need to check "edge" that goes through face twice */
2896                 if((pce->x[3] && i==2))
2897                         continue;
2898
2899                 cur = edge+i;
2900                 cur->x[0] = pce->x[i]; cur->x[1] = pce->x[(i+1)%3];
2901                 cur->v[0] = pce->v[i]; cur->v[1] = pce->v[(i+1)%3];
2902                 cur->tot = 2;
2903                 cur->inside = 0;
2904
2905                 ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_edge);
2906
2907                 if(ct >= 0.f && ct < *t) {
2908                         float u, e[3], vec[3];
2909
2910                         sub_v3_v3v3(e, cur->x1, cur->x0);
2911                         sub_v3_v3v3(vec, cur->p, cur->x0);
2912                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2913
2914                         if(u < 0.f || u > 1.f)
2915                                 break;
2916
2917                         *result = *cur;
2918
2919                         madd_v3_v3v3fl(result->nor, vec, e, -u);
2920                         normalize_v3(result->nor);
2921
2922                         result->uv[0] = u;
2923
2924                         
2925                         hit = cur;
2926                         *t = ct;
2927                 }
2928
2929         }
2930
2931         return hit != NULL;
2932 }
2933 static int collision_sphere_to_verts(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2934 {
2935         ParticleCollisionElement vert[3], *cur = NULL, *hit = NULL;
2936         ParticleCollisionElement *result = &col->pce;
2937
2938         float ct;
2939         int i;
2940
2941         for(i=0; i<3; i++) {
2942                 /* in case of quad, only check one vert the first time */
2943                 if(pce->x[3] && i != 1)
2944                         continue;
2945
2946                 cur = vert+i;
2947                 cur->x[0] = pce->x[i];
2948                 cur->v[0] = pce->v[i];
2949                 cur->tot = 1;
2950                 cur->inside = 0;
2951
2952                 ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_vert);
2953                 
2954                 if(ct >= 0.f && ct < *t) {
2955                         *result = *cur;
2956
2957                         sub_v3_v3v3(result->nor, cur->p, cur->x0);
2958                         normalize_v3(result->nor);
2959
2960                         hit = cur;
2961                         *t = ct;
2962                 }
2963
2964         }
2965
2966         return hit != NULL;
2967 }
2968 /* Callback for BVHTree near test */
2969 void BKE_psys_collision_neartest_cb(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
2970 {
2971         ParticleCollision *col = (ParticleCollision *) userdata;
2972         ParticleCollisionElement pce;
2973         MFace *face = col->md->mfaces + index;
2974         MVert *x = col->md->x;
2975         MVert *v = col->md->current_v;
2976         float t = hit->dist/col->original_ray_length;
2977         int collision = 0;
2978
2979         pce.x[0] = x[face->v1].co;
2980         pce.x[1] = x[face->v2].co;
2981         pce.x[2] = x[face->v3].co;
2982         pce.x[3] = face->v4 ? x[face->v4].co : NULL;
2983
2984         pce.v[0] = v[face->v1].co;
2985         pce.v[1] = v[face->v2].co;
2986         pce.v[2] = v[face->v3].co;
2987         pce.v[3] = face->v4 ? v[face->v4].co : NULL;
2988
2989         pce.tot = 3;
2990         pce.inside = 0;
2991         pce.index = index;
2992
2993         /* don't collide with same face again */
2994         if(col->hit == col->current && col->pce.index == index && col->pce.tot == 3)
2995                 return;
2996
2997         do
2998         {       
2999                 collision = collision_sphere_to_tri(col, ray->radius, &pce, &t);
3000                 if(col->pce.inside == 0) {
3001                         collision += collision_sphere_to_edges(col, ray->radius, &pce, &t);
3002                         collision += collision_sphere_to_verts(col, ray->radius, &pce, &t);
3003                 }
3004
3005                 if(collision) {
3006                         hit->dist = col->original_ray_length * t;
3007                         hit->index = index;
3008                                 
3009                         collision_point_velocity(&col->pce);
3010
3011                         col->hit = col->current;
3012                 }
3013
3014                 pce.x[1] = pce.x[2];
3015                 pce.x[2] = pce.x[3];
3016                 pce.x[3] = NULL;
3017
3018                 pce.v[1] = pce.v[2];
3019                 pce.v[2] = pce.v[3];
3020                 pce.v[3] = NULL;
3021
3022         } while(pce.x[2]);
3023 }
3024 static int collision_detect(ParticleData *pa, ParticleCollision *col, BVHTreeRayHit *hit, ListBase *colliders)
3025 {
3026         ColliderCache *coll;