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