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