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