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