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