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