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