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