Particles todo item: particle textures
[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         ParticleSystem *psys = sim->psys;
1523         ParticleSettings *part = psys->part;
1524         ParticleTexture ptex;
1525
1526         pa->flag &= ~PARS_UNEXIST;
1527
1528         if(part->from != PART_FROM_PARTICLE && part->type != PART_FLUID) {
1529                 psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);
1530                 
1531                 if(ptex.exist < PSYS_FRAND(p+125))
1532                         pa->flag |= PARS_UNEXIST;
1533
1534                 pa->time = (part->type == PART_HAIR) ? 0.f : part->sta + (part->end - part->sta)*ptex.time;
1535         }
1536
1537         pa->hair_index = 0;
1538         /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
1539         /* usage other than straight after distribute has to handle this index by itself - jahka*/
1540         //pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we dont have a derived mesh face */
1541 }
1542 static void initialize_all_particles(ParticleSimulationData *sim)
1543 {
1544         ParticleSystem *psys = sim->psys;
1545         PARTICLE_P;
1546
1547         psys->totunexist = 0;
1548
1549         LOOP_PARTICLES {
1550                 if((pa->flag & PARS_UNEXIST)==0)
1551                         initialize_particle(sim, pa, p);
1552
1553                 if(pa->flag & PARS_UNEXIST)
1554                         psys->totunexist++;
1555         }
1556
1557         /* Free unexisting particles. */
1558         if(psys->totpart && psys->totunexist == psys->totpart) {
1559                 if(psys->particles->boid)
1560                         MEM_freeN(psys->particles->boid);
1561
1562                 MEM_freeN(psys->particles);
1563                 psys->particles = NULL;
1564                 psys->totpart = psys->totunexist = 0;
1565         }
1566
1567         if(psys->totunexist) {
1568                 int newtotpart = psys->totpart - psys->totunexist;
1569                 ParticleData *npa, *newpars;
1570                 
1571                 npa = newpars = MEM_callocN(newtotpart * sizeof(ParticleData), "particles");
1572
1573                 for(p=0, pa=psys->particles; p<newtotpart; p++, pa++, npa++) {
1574                         while(pa->flag & PARS_UNEXIST)
1575                                 pa++;
1576
1577                         memcpy(npa, pa, sizeof(ParticleData));
1578                 }
1579
1580                 if(psys->particles->boid)
1581                         MEM_freeN(psys->particles->boid);
1582                 MEM_freeN(psys->particles);
1583                 psys->particles = newpars;
1584                 psys->totpart -= psys->totunexist;
1585
1586                 if(psys->particles->boid) {
1587                         BoidParticle *newboids = MEM_callocN(psys->totpart * sizeof(BoidParticle), "boid particles");
1588
1589                         LOOP_PARTICLES
1590                                 pa->boid = newboids++;
1591
1592                 }
1593         }
1594 }
1595 /* sets particle to the emitter surface with initial velocity & rotation */
1596 void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
1597 {
1598         Object *ob = sim->ob;
1599         ParticleSystem *psys = sim->psys;
1600         ParticleSettings *part;
1601         ParticleTexture ptex;
1602         float fac, phasefac, nor[3]={0,0,0},loc[3],vel[3]={0.0,0.0,0.0},rot[4],q2[4];
1603         float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3]={0.0,0.0,0.0};
1604         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};
1605         float q_phase[4];
1606         int p = pa - psys->particles;
1607         part=psys->part;
1608         
1609 #if 0 /* deprecated code */
1610         if(part->from==PART_FROM_PARTICLE){
1611                 float speed;
1612                 ParticleSimulationData tsim= {0};
1613                 tsim.scene= sim->scene;
1614                 tsim.ob= psys->target_ob ? psys->target_ob : ob;
1615                 tsim.psys = BLI_findlink(&tsim.ob->particlesystem, sim->psys->target_psys-1);
1616
1617                 state.time = pa->time;
1618                 if(pa->num == -1)
1619                         memset(&state, 0, sizeof(state));
1620                 else
1621                         psys_get_particle_state(&tsim, pa->num, &state, 1);
1622                 psys_get_from_key(&state, loc, nor, rot, 0);
1623
1624                 mul_qt_v3(rot, vtan);
1625                 mul_qt_v3(rot, utan);
1626
1627                 speed= normalize_v3_v3(p_vel, state.vel);
1628                 mul_v3_fl(p_vel, dot_v3v3(r_vel, p_vel));
1629                 VECSUB(p_vel, r_vel, p_vel);
1630                 normalize_v3(p_vel);
1631                 mul_v3_fl(p_vel, speed);
1632
1633                 VECCOPY(pa->fuv, loc); /* abusing pa->fuv (not used for "from particle") for storing emit location */
1634         }
1635         else{
1636 #endif
1637         /* get precise emitter matrix if particle is born */
1638         if(part->type!=PART_HAIR && pa->time < cfra && pa->time >= sim->psys->cfra) {
1639                 /* we have to force RECALC_ANIM here since where_is_objec_time only does drivers */
1640                 BKE_animsys_evaluate_animdata(&sim->ob->id, sim->ob->adt, pa->time, ADT_RECALC_ANIM);
1641                 where_is_object_time(sim->scene, sim->ob, pa->time);
1642         }
1643
1644         /* get birth location from object               */
1645         if(part->tanfac != 0.f)
1646                 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0,0);
1647         else
1648                 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,0,0,0,0);
1649                 
1650         /* get possible textural influence */
1651         psys_get_texture(sim, pa, &ptex, PAMAP_IVEL|PAMAP_LIFE, cfra);
1652
1653         /* particles live in global space so    */
1654         /* let's convert:                                               */
1655         /* -location                                                    */
1656         mul_m4_v3(ob->obmat, loc);
1657                 
1658         /* -normal                                                              */
1659         mul_mat3_m4_v3(ob->obmat, nor);
1660         normalize_v3(nor);
1661
1662         /* -tangent                                                             */
1663         if(part->tanfac!=0.0){
1664                 //float phase=vg_rot?2.0f*(psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_rot)-0.5f):0.0f;
1665                 float phase=0.0f;
1666                 mul_v3_fl(vtan,-(float)cos(M_PI*(part->tanphase+phase)));
1667                 fac=-(float)sin(M_PI*(part->tanphase+phase));
1668                 VECADDFAC(vtan,vtan,utan,fac);
1669
1670                 mul_mat3_m4_v3(ob->obmat,vtan);
1671
1672                 VECCOPY(utan,nor);
1673                 mul_v3_fl(utan,dot_v3v3(vtan,nor));
1674                 VECSUB(vtan,vtan,utan);
1675                         
1676                 normalize_v3(vtan);
1677         }
1678                 
1679
1680         /* -velocity                                                    */
1681         if(part->randfac!=0.0){
1682                 r_vel[0] = 2.0f * (PSYS_FRAND(p + 10) - 0.5f);
1683                 r_vel[1] = 2.0f * (PSYS_FRAND(p + 11) - 0.5f);
1684                 r_vel[2] = 2.0f * (PSYS_FRAND(p + 12) - 0.5f);
1685
1686                 mul_mat3_m4_v3(ob->obmat, r_vel);
1687                 normalize_v3(r_vel);
1688         }
1689
1690         /* -angular velocity                                    */
1691         if(part->avemode==PART_AVE_RAND){
1692                 r_ave[0] = 2.0f * (PSYS_FRAND(p + 13) - 0.5f);
1693                 r_ave[1] = 2.0f * (PSYS_FRAND(p + 14) - 0.5f);
1694                 r_ave[2] = 2.0f * (PSYS_FRAND(p + 15) - 0.5f);
1695
1696                 mul_mat3_m4_v3(ob->obmat,r_ave);
1697                 normalize_v3(r_ave);
1698         }
1699                 
1700         /* -rotation                                                    */
1701         if(part->randrotfac != 0.0f){
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                 mat4_to_quat(rot,ob->obmat);
1709                 mul_qt_qtqt(r_rot,r_rot,rot);
1710         }
1711 #if 0
1712         }
1713 #endif
1714
1715         if(part->phystype==PART_PHYS_BOIDS && pa->boid) {
1716                 BoidParticle *bpa = pa->boid;
1717                 float dvec[3], q[4], mat[3][3];
1718
1719                 copy_v3_v3(pa->state.co,loc);
1720
1721                 /* boids don't get any initial velocity  */
1722                 zero_v3(pa->state.vel);
1723
1724                 /* boids store direction in ave */
1725                 if(fabs(nor[2])==1.0f) {
1726                         sub_v3_v3v3(pa->state.ave, loc, ob->obmat[3]);
1727                         normalize_v3(pa->state.ave);
1728                 }
1729                 else {
1730                         VECCOPY(pa->state.ave, nor);
1731                 }
1732                 /* and gravity in r_ve */
1733                 bpa->gravity[0] = bpa->gravity[1] = 0.0f;
1734                 bpa->gravity[2] = -1.0f;
1735                 if((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY)
1736                         && sim->scene->physics_settings.gravity[2]!=0.0f)
1737                         bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
1738
1739                 /* calculate rotation matrix */
1740                 project_v3_v3v3(dvec, r_vel, pa->state.ave);
1741                 sub_v3_v3v3(mat[0], pa->state.ave, dvec);
1742                 normalize_v3(mat[0]);
1743                 negate_v3_v3(mat[2], r_vel);
1744                 normalize_v3(mat[2]);
1745                 cross_v3_v3v3(mat[1], mat[2], mat[0]);
1746                 
1747                 /* apply rotation */
1748                 mat3_to_quat_is_ok( q,mat);
1749                 copy_qt_qt(pa->state.rot, q);
1750
1751                 bpa->data.health = part->boids->health;
1752                 bpa->data.mode = eBoidMode_InAir;
1753                 bpa->data.state_id = ((BoidState*)part->boids->states.first)->id;
1754                 bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f;
1755         }
1756         else {
1757                 /* conversion done so now we apply new: */
1758                 /* -velocity from:                                              */
1759
1760                 /*              *reactions                                              */
1761                 if(dtime > 0.f){
1762                         sub_v3_v3v3(vel, pa->state.vel, pa->prev_state.vel);
1763                 }
1764
1765                 /*              *emitter velocity                               */
1766                 if(dtime != 0.f && part->obfac != 0.f){
1767                         sub_v3_v3v3(vel, loc, pa->state.co);
1768                         mul_v3_fl(vel, part->obfac/dtime);
1769                 }
1770                 
1771                 /*              *emitter normal                                 */
1772                 if(part->normfac != 0.f)
1773                         madd_v3_v3fl(vel, nor, part->normfac);
1774                 
1775                 /*              *emitter tangent                                */
1776                 if(sim->psmd && part->tanfac != 0.f)
1777                         madd_v3_v3fl(vel, vtan, part->tanfac);
1778
1779                 /*              *emitter object orientation             */
1780                 if(part->ob_vel[0] != 0.f) {
1781                         normalize_v3_v3(vec, ob->obmat[0]);
1782                         madd_v3_v3fl(vel, vec, part->ob_vel[0]);
1783                 }
1784                 if(part->ob_vel[1] != 0.f) {
1785                         normalize_v3_v3(vec, ob->obmat[1]);
1786                         madd_v3_v3fl(vel, vec, part->ob_vel[1]);
1787                 }
1788                 if(part->ob_vel[2] != 0.f) {
1789                         normalize_v3_v3(vec, ob->obmat[2]);
1790                         madd_v3_v3fl(vel, vec, part->ob_vel[2]);
1791                 }
1792
1793                 /*              *texture                                                */
1794                 /* TODO */
1795
1796                 /*              *random                                                 */
1797                 if(part->randfac != 0.f)
1798                         madd_v3_v3fl(vel, r_vel, part->randfac);
1799
1800                 /*              *particle                                               */
1801                 if(part->partfac != 0.f)
1802                         madd_v3_v3fl(vel, p_vel, part->partfac);
1803                 
1804                 mul_v3_v3fl(pa->state.vel, vel, ptex.ivel);
1805
1806                 /* -location from emitter                               */
1807                 copy_v3_v3(pa->state.co,loc);
1808
1809                 /* -rotation                                                    */
1810                 unit_qt(pa->state.rot);
1811
1812                 if(part->rotmode){
1813                         /* create vector into which rotation is aligned */
1814                         switch(part->rotmode){
1815                                 case PART_ROT_NOR:
1816                                         copy_v3_v3(rot_vec, nor);
1817                                         break;
1818                                 case PART_ROT_VEL:
1819                                         copy_v3_v3(rot_vec, vel);
1820                                         break;
1821                                 case PART_ROT_GLOB_X:
1822                                 case PART_ROT_GLOB_Y:
1823                                 case PART_ROT_GLOB_Z:
1824                                         rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
1825                                         break;
1826                                 case PART_ROT_OB_X:
1827                                 case PART_ROT_OB_Y:
1828                                 case PART_ROT_OB_Z:
1829                                         copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
1830                                         break;
1831                         }
1832                         
1833                         /* create rotation quat */
1834                         negate_v3(rot_vec);
1835                         vec_to_quat( q2,rot_vec, OB_POSX, OB_POSZ);
1836
1837                         /* randomize rotation quat */
1838                         if(part->randrotfac!=0.0f)
1839                                 interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
1840                         else
1841                                 copy_qt_qt(rot,q2);
1842
1843                         /* rotation phase */
1844                         phasefac = part->phasefac;
1845                         if(part->randphasefac != 0.0f)
1846                                 phasefac += part->randphasefac * PSYS_FRAND(p + 20);
1847                         axis_angle_to_quat( q_phase,x_vec, phasefac*(float)M_PI);
1848
1849                         /* combine base rotation & phase */
1850                         mul_qt_qtqt(pa->state.rot, rot, q_phase);
1851                 }
1852
1853                 /* -angular velocity                                    */
1854
1855                 zero_v3(pa->state.ave);
1856
1857                 if(part->avemode){
1858                         switch(part->avemode){
1859                                 case PART_AVE_SPIN:
1860                                         copy_v3_v3(pa->state.ave, vel);
1861                                         break;
1862                                 case PART_AVE_RAND:
1863                                         copy_v3_v3(pa->state.ave, r_ave);
1864                                         break;
1865                         }
1866                         normalize_v3(pa->state.ave);
1867                         mul_v3_fl(pa->state.ave,part->avefac);
1868                 }
1869         }
1870
1871
1872         if(part->type == PART_HAIR){
1873                 pa->lifetime = 100.0f;
1874         }
1875         else{
1876                 pa->lifetime = part->lifetime * ptex.life;
1877
1878                 if(part->randlife != 0.0)
1879                         pa->lifetime *= 1.0f - part->randlife * PSYS_FRAND(p + 21);
1880         }
1881
1882         pa->dietime = pa->time + pa->lifetime;
1883
1884         if(sim->psys->pointcache && sim->psys->pointcache->flag & PTCACHE_BAKED &&
1885                 sim->psys->pointcache->mem_cache.first) {
1886                 float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1887                 pa->dietime = MIN2(pa->dietime, dietime);
1888         }
1889
1890         if(pa->time > cfra)
1891                 pa->alive = PARS_UNBORN;
1892         else if(pa->dietime <= cfra)
1893                 pa->alive = PARS_DEAD;
1894         else
1895                 pa->alive = PARS_ALIVE;
1896
1897         pa->state.time = cfra;
1898 }
1899 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1900 {
1901         ParticleData *pa;
1902         int p, totpart=sim->psys->totpart;
1903         
1904         for(p=from, pa=sim->psys->particles+from; p<totpart; p++, pa++)
1905                 reset_particle(sim, pa, dtime, cfra);
1906 }
1907 /************************************************/
1908 /*                      Particle targets                                        */
1909 /************************************************/
1910 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1911 {
1912         ParticleSystem *psys = NULL;
1913
1914         if(pt->ob == NULL || pt->ob == ob)
1915                 psys = BLI_findlink(&ob->particlesystem, pt->psys-1);
1916         else
1917                 psys = BLI_findlink(&pt->ob->particlesystem, pt->psys-1);
1918
1919         if(psys)
1920                 pt->flag |= PTARGET_VALID;
1921         else
1922                 pt->flag &= ~PTARGET_VALID;
1923
1924         return psys;
1925 }
1926 /************************************************/
1927 /*                      Keyed particles                                         */
1928 /************************************************/
1929 /* Counts valid keyed targets */
1930 void psys_count_keyed_targets(ParticleSimulationData *sim)
1931 {
1932         ParticleSystem *psys = sim->psys, *kpsys;
1933         ParticleTarget *pt = psys->targets.first;
1934         int keys_valid = 1;
1935         psys->totkeyed = 0;
1936
1937         for(; pt; pt=pt->next) {
1938                 kpsys = psys_get_target_system(sim->ob, pt);
1939
1940                 if(kpsys && kpsys->totpart) {
1941                         psys->totkeyed += keys_valid;
1942                         if(psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
1943                                 psys->totkeyed += 1;
1944                 }
1945                 else {
1946                         keys_valid = 0;
1947                 }
1948         }
1949
1950         psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1951 }
1952
1953 static void set_keyed_keys(ParticleSimulationData *sim)
1954 {
1955         ParticleSystem *psys = sim->psys;
1956         ParticleSimulationData ksim= {0};
1957         ParticleTarget *pt;
1958         PARTICLE_P;
1959         ParticleKey *key;
1960         int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1961
1962         ksim.scene= sim->scene;
1963         
1964         /* no proper targets so let's clear and bail out */
1965         if(psys->totkeyed==0) {
1966                 free_keyed_keys(psys);
1967                 psys->flag &= ~PSYS_KEYED;
1968                 return;
1969         }
1970
1971         if(totpart && psys->particles->totkey != totkeys) {
1972                 free_keyed_keys(psys);
1973                 
1974                 key = MEM_callocN(totpart*totkeys*sizeof(ParticleKey), "Keyed keys");
1975                 
1976                 LOOP_PARTICLES {
1977                         pa->keys = key;
1978                         pa->totkey = totkeys;
1979                         key += totkeys;
1980                 }
1981         }
1982         
1983         psys->flag &= ~PSYS_KEYED;
1984
1985
1986         pt = psys->targets.first;
1987         for(k=0; k<totkeys; k++) {
1988                 ksim.ob = pt->ob ? pt->ob : sim->ob;
1989                 ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
1990
1991                 LOOP_PARTICLES {
1992                         key = pa->keys + k;
1993                         key->time = -1.0; /* use current time */
1994
1995                         psys_get_particle_state(&ksim, p%ksim.psys->totpart, key, 1);
1996
1997                         if(psys->flag & PSYS_KEYED_TIMING){
1998                                 key->time = pa->time + pt->time;
1999                                 if(pt->duration != 0.0f && k+1 < totkeys) {
2000                                         copy_particle_key(key+1, key, 1);
2001                                         (key+1)->time = pa->time + pt->time + pt->duration;
2002                                 }
2003                         }
2004                         else if(totkeys > 1)
2005                                 key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
2006                         else
2007                                 key->time = pa->time;
2008                 }
2009
2010                 if(psys->flag & PSYS_KEYED_TIMING && pt->duration!=0.0f)
2011                         k++;
2012
2013                 pt = (pt->next && pt->next->flag & PTARGET_VALID)? pt->next : psys->targets.first;
2014         }
2015
2016         psys->flag |= PSYS_KEYED;
2017 }
2018 /************************************************/
2019 /*                      Reactors                                                        */
2020 /************************************************/
2021 //static void push_reaction(ParticleSimulationData *sim, int pa_num, int event, ParticleKey *state)
2022 //{
2023 //      Object *rob;
2024 //      ParticleSystem *rpsys;
2025 //      ParticleSettings *rpart;
2026 //      ParticleData *pa;
2027 //      ListBase *lb=&sim->psys->effectors;
2028 //      ParticleEffectorCache *ec;
2029 //      ParticleReactEvent *re;
2030 //
2031 //      if(lb->first) for(ec = lb->first; ec; ec= ec->next){
2032 //              if(ec->type & PSYS_EC_REACTOR){
2033 //                      /* all validity checks already done in add_to_effectors */
2034 //                      rob=ec->ob;
2035 //                      rpsys=BLI_findlink(&rob->particlesystem,ec->psys_nbr);
2036 //                      rpart=rpsys->part;
2037 //                      if(rpsys->part->reactevent==event){
2038 //                              pa=sim->psys->particles+pa_num;
2039 //                              re= MEM_callocN(sizeof(ParticleReactEvent), "react event");
2040 //                              re->event=event;
2041 //                              re->pa_num = pa_num;
2042 //                              re->ob = sim->ob;
2043 //                              re->psys = sim->psys;
2044 //                              re->size = pa->size;
2045 //                              copy_particle_key(&re->state,state,1);
2046 //
2047 //                              switch(event){
2048 //                                      case PART_EVENT_DEATH:
2049 //                                              re->time=pa->dietime;
2050 //                                              break;
2051 //                                      case PART_EVENT_COLLIDE:
2052 //                                              re->time=state->time;
2053 //                                              break;
2054 //                                      case PART_EVENT_NEAR:
2055 //                                              re->time=state->time;
2056 //                                              break;
2057 //                              }
2058 //
2059 //                              BLI_addtail(&rpsys->reactevents, re);
2060 //                      }
2061 //              }
2062 //      }
2063 //}
2064 //static void react_to_events(ParticleSystem *psys, int pa_num)
2065 //{
2066 //      ParticleSettings *part=psys->part;
2067 //      ParticleData *pa=psys->particles+pa_num;
2068 //      ParticleReactEvent *re=psys->reactevents.first;
2069 //      int birth=0;
2070 //      float dist=0.0f;
2071 //
2072 //      for(re=psys->reactevents.first; re; re=re->next){
2073 //              birth=0;
2074 //              if(part->from==PART_FROM_PARTICLE){
2075 //                      if(pa->num==re->pa_num && pa->alive==PARS_UNBORN){
2076 //                              if(re->event==PART_EVENT_NEAR){
2077 //                                      ParticleData *tpa = re->psys->particles+re->pa_num;
2078 //                                      float pa_time=tpa->time + pa->foffset*tpa->lifetime;
2079 //                                      if(re->time >= pa_time){
2080 //                                              pa->time=pa_time;
2081 //                                              pa->dietime=pa->time+pa->lifetime;
2082 //                                      }
2083 //                              }
2084 //                              else{
2085 //                                      pa->time=re->time;
2086 //                                      pa->dietime=pa->time+pa->lifetime;
2087 //                              }
2088 //                      }
2089 //              }
2090 //              else{
2091 //                      dist=len_v3v3(pa->state.co, re->state.co);
2092 //                      if(dist <= re->size){
2093 //                              if(pa->alive==PARS_UNBORN){
2094 //                                      pa->time=re->time;
2095 //                                      pa->dietime=pa->time+pa->lifetime;
2096 //                                      birth=1;
2097 //                              }
2098 //                              if(birth || part->flag&PART_REACT_MULTIPLE){
2099 //                                      float vec[3];
2100 //                                      VECSUB(vec,pa->state.co, re->state.co);
2101 //                                      if(birth==0)
2102 //                                              mul_v3_fl(vec,(float)pow(1.0f-dist/re->size,part->reactshape));
2103 //                                      VECADDFAC(pa->state.vel,pa->state.vel,vec,part->reactfac);
2104 //                                      VECADDFAC(pa->state.vel,pa->state.vel,re->state.vel,part->partfac);
2105 //                              }
2106 //                              if(birth)
2107 //                                      mul_v3_fl(pa->state.vel,(float)pow(1.0f-dist/re->size,part->reactshape));
2108 //                      }
2109 //              }
2110 //      }
2111 //}
2112 //void psys_get_reactor_target(ParticleSimulationData *sim, Object **target_ob, ParticleSystem **target_psys)
2113 //{
2114 //      Object *tob;
2115 //
2116 //      tob = sim->psys->target_ob ? sim->psys->target_ob : sim->ob;
2117 //      
2118 //      *target_psys = BLI_findlink(&tob->particlesystem, sim->psys->target_psys-1);
2119 //      if(*target_psys)
2120 //              *target_ob=tob;
2121 //      else
2122 //              *target_ob=0;
2123 //}
2124 /************************************************/
2125 /*                      Point Cache                                                     */
2126 /************************************************/
2127 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
2128 {
2129         PointCache *cache = psys->pointcache;
2130
2131         if(cache->flag & PTCACHE_DISK_CACHE && cache->mem_cache.first == NULL) {
2132                 PTCacheID pid;
2133                 BKE_ptcache_id_from_particles(&pid, ob, psys);
2134                 cache->flag &= ~PTCACHE_DISK_CACHE;
2135                 BKE_ptcache_disk_to_mem(&pid);
2136                 cache->flag |= PTCACHE_DISK_CACHE;
2137         }
2138 }
2139 static void psys_clear_temp_pointcache(ParticleSystem *psys)
2140 {
2141         if(psys->pointcache->flag & PTCACHE_DISK_CACHE)
2142                 BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
2143 }
2144 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
2145 {
2146         ParticleSettings *part = psys->part;
2147
2148         *sfra = MAX2(1, (int)part->sta);
2149         *efra = MIN2((int)(part->end + part->lifetime + 1.0), scene->r.efra);
2150 }
2151
2152 /************************************************/
2153 /*                      Effectors                                                       */
2154 /************************************************/
2155 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
2156 {
2157         if(psys) {
2158                 PARTICLE_P;
2159
2160                 if(!psys->tree || psys->tree_frame != cfra) {
2161                         
2162                         BLI_kdtree_free(psys->tree);
2163
2164                         psys->tree = BLI_kdtree_new(psys->totpart);
2165                         
2166                         LOOP_SHOWN_PARTICLES {
2167                                 if(pa->alive == PARS_ALIVE) {
2168                                         if(pa->state.time == cfra)
2169                                                 BLI_kdtree_insert(psys->tree, p, pa->prev_state.co, NULL);
2170                                         else
2171                                                 BLI_kdtree_insert(psys->tree, p, pa->state.co, NULL);
2172                                 }
2173                         }
2174                         BLI_kdtree_balance(psys->tree);
2175
2176                         psys->tree_frame = psys->cfra;
2177                 }
2178         }
2179 }
2180
2181 static void psys_update_effectors(ParticleSimulationData *sim)
2182 {
2183         pdEndEffectors(&sim->psys->effectors);
2184         sim->psys->effectors = pdInitEffectors(sim->scene, sim->ob, sim->psys, sim->psys->part->effector_weights);
2185         precalc_guides(sim, sim->psys->effectors);
2186 }
2187
2188 /*********************************************************************************************************
2189                     SPH fluid physics 
2190
2191  In theory, there could be unlimited implementation of SPH simulators
2192
2193  This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
2194
2195  Titled: Particle-based Viscoelastic Fluid Simulation.
2196  Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
2197  Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
2198
2199  Presented at Siggraph, (2005)
2200
2201 ***********************************************************************************************************/
2202 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
2203 ParticleSpring *add_fluid_spring(ParticleSystem *psys, ParticleSpring *spring)
2204 {
2205         /* Are more refs required? */
2206         if(psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
2207                 psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
2208                 psys->fluid_springs = (ParticleSpring*)MEM_callocN(psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
2209         }
2210         else if(psys->tot_fluidsprings == psys->alloc_fluidsprings) {
2211                 /* Double the number of refs allocated */
2212                 psys->alloc_fluidsprings *= 2;
2213                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
2214         }
2215
2216         memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
2217         psys->tot_fluidsprings++;
2218
2219         return psys->fluid_springs + psys->tot_fluidsprings - 1;
2220 }
2221
2222 void  delete_fluid_spring(ParticleSystem *psys, int j)
2223 {
2224         if (j != psys->tot_fluidsprings - 1)
2225                 psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
2226
2227         psys->tot_fluidsprings--;
2228
2229         if (psys->tot_fluidsprings < psys->alloc_fluidsprings/2 && psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE){
2230                 psys->alloc_fluidsprings /= 2;
2231                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs,  psys->alloc_fluidsprings * sizeof(ParticleSpring));
2232         }
2233 }
2234
2235 EdgeHash *build_fluid_springhash(ParticleSystem *psys)
2236 {
2237         EdgeHash *springhash = NULL;
2238         ParticleSpring *spring;
2239         int i = 0;
2240
2241         springhash = BLI_edgehash_new();
2242
2243         for(i=0, spring=psys->fluid_springs; i<psys->tot_fluidsprings; i++, spring++)
2244                 BLI_edgehash_insert(springhash, spring->particle_index[0], spring->particle_index[1], SET_INT_IN_POINTER(i+1));
2245
2246         return springhash;
2247 }
2248 static void particle_fluidsim(ParticleSystem *psys, int own_psys, ParticleData *pa, float dtime, float mass, float *gravity, EdgeHash *springhash)
2249 {
2250         SPHFluidSettings *fluid = psys->part->fluid;
2251         KDTreeNearest *ptn = NULL;
2252         ParticleData *npa;
2253         ParticleSpring *spring = NULL;
2254
2255         float temp[3];
2256         float q, q1, u, I, D, rij, d, Lij;
2257         float pressure_near, pressure;
2258         float p=0, pnear=0;
2259
2260         float omega = fluid->viscosity_omega;
2261         float beta = fluid->viscosity_beta;
2262         float massfactor = 1.0f/mass;
2263         float spring_k = fluid->spring_k;
2264         float h = fluid->radius;
2265         float L = fluid->rest_length * fluid->radius;
2266
2267         int n, neighbours = BLI_kdtree_range_search(psys->tree, h, pa->prev_state.co, NULL, &ptn);
2268         int spring_index = 0, index = own_psys ? pa - psys->particles : -1;
2269
2270         /* pressure and near pressure */
2271         for(n=own_psys?1:0; n<neighbours; n++) {
2272                 sub_v3_v3(ptn[n].co, pa->prev_state.co);
2273                 mul_v3_fl(ptn[n].co, 1.f/ptn[n].dist);
2274                 q = ptn[n].dist/h;
2275
2276                 if(q < 1.f) {
2277                         q1 = 1.f - q;
2278
2279                         p += q1*q1;
2280                         pnear += q1*q1*q1;
2281                 }
2282         }
2283
2284         p *= mass;
2285         pnear *= mass;
2286         pressure =  fluid->stiffness_k * (p - fluid->rest_density);
2287         pressure_near = fluid->stiffness_knear * pnear;
2288
2289         /* main calculations */
2290         for(n=own_psys?1:0; n<neighbours; n++) {
2291                 npa = psys->particles + ptn[n].index;
2292
2293                 rij = ptn[n].dist;
2294                 q = rij/h;
2295                 q1 = 1.f-q;
2296
2297                 /* Double Density Relaxation - Algorithm 2 (can't be thread safe!)*/
2298                 D =  dtime * dtime * (pressure + pressure_near*q1)*q1 * 0.5f;
2299                 madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
2300                 if(own_psys)
2301                         madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
2302
2303                 if(index < ptn[n].index) {
2304                         /* Viscosity - Algorithm 5 */
2305                         if(omega > 0.f  || beta > 0.f) {                
2306                                 sub_v3_v3v3(temp, pa->state.vel, npa->state.vel);
2307                                 u = dot_v3v3(ptn[n].co, temp);
2308
2309                                 if (u > 0){
2310                                         I = dtime * (q1 * (omega * u + beta * u*u)) * 0.5f;
2311                                         madd_v3_v3fl(pa->state.vel, ptn[n].co, -I * massfactor);
2312
2313                                         if(own_psys)
2314                                                 madd_v3_v3fl(npa->state.vel, ptn[n].co, I * massfactor);
2315                                 }
2316                         }
2317
2318                         if(spring_k > 0.f) {
2319                                 /* Viscoelastic spring force - Algorithm 4*/
2320                                 if (fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash){
2321                                         spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, ptn[n].index));
2322
2323                                         if(spring_index) {
2324                                                 spring = psys->fluid_springs + spring_index - 1;
2325                                         }
2326                                         else {
2327                                                 ParticleSpring temp_spring;
2328                                                 temp_spring.particle_index[0] = index;
2329                                                 temp_spring.particle_index[1] = ptn[n].index;
2330                                                 temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : L;
2331                                                 temp_spring.delete_flag = 0;
2332                                                 
2333                                                 spring = add_fluid_spring(psys, &temp_spring);
2334                                         }
2335
2336                                         Lij = spring->rest_length;
2337                                         d = fluid->yield_ratio * Lij;
2338
2339                                         if (rij > Lij + d) // Stretch, 25 is just a multiplier for plasticity_constant value to counter default dtime of 1/25
2340                                                 spring->rest_length += dtime * 25.f * fluid->plasticity_constant * (rij - Lij - d);
2341                                         else if(rij < Lij - d) // Compress
2342                                                 spring->rest_length -= dtime * 25.f * fluid->plasticity_constant * (Lij - d - rij);
2343                                 }
2344                                 else { /* PART_SPRING_HOOKES - Hooke's spring force */
2345                                         /* L is a factor of radius */
2346                                         D = 0.5 * dtime * dtime * 10.f * fluid->spring_k * (1.f - L/h) * (L - rij);
2347  
2348                                         madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
2349                                         if(own_psys)
2350                                                 madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
2351                                 }
2352                         }
2353                 }
2354         } 
2355
2356         /* Artificial buoyancy force in negative gravity direction  */
2357         if (fluid->buoyancy >= 0.f && gravity) {
2358                 float B = -dtime * dtime * fluid->buoyancy * (p - fluid->rest_density) * 0.5f;
2359                 madd_v3_v3fl(pa->state.co, gravity, -B * massfactor);
2360         }
2361
2362         if(ptn)
2363                 MEM_freeN(ptn);
2364 }
2365
2366 static void apply_particle_fluidsim(Object *ob, ParticleSystem *psys, ParticleData *pa, float dtime, float *gravity, EdgeHash *springhash){
2367         ParticleTarget *pt;
2368
2369         particle_fluidsim(psys, 1, pa, dtime, psys->part->mass, gravity, springhash);
2370         
2371         /*----check other SPH systems (Multifluids) , each fluid has its own parameters---*/
2372         for(pt=psys->targets.first; pt; pt=pt->next) {
2373                 ParticleSystem *epsys = psys_get_target_system(ob, pt);
2374
2375                 if(epsys)
2376                         particle_fluidsim(epsys, 0, pa, dtime, psys->part->mass, gravity, NULL);
2377         }
2378         /*----------------------------------------------------------------*/             
2379 }
2380
2381 static void apply_fluid_springs(ParticleSystem *psys, float timestep){
2382         SPHFluidSettings *fluid = psys->part->fluid;
2383         ParticleData *pa1, *pa2;
2384         ParticleSpring *spring = psys->fluid_springs;
2385         
2386         float h = fluid->radius;
2387         float massfactor = 1.0f/psys->part->mass;
2388         float D, Rij[3], rij, Lij;
2389         int i;
2390
2391         if((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
2392                 return;
2393
2394         /* Loop through the springs */
2395         for(i=0; i<psys->tot_fluidsprings; i++, spring++) {
2396                 Lij = spring->rest_length;
2397
2398                 if (Lij > h) {
2399                         spring->delete_flag = 1;
2400                 }
2401                 else {
2402                         pa1 = psys->particles + spring->particle_index[0];
2403                         pa2 = psys->particles + spring->particle_index[1];
2404
2405                         sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
2406                         rij = normalize_v3(Rij);
2407
2408                         /* Calculate displacement and apply value */
2409                         D =  0.5f * timestep * timestep * 10.f * fluid->spring_k * (1.f - Lij/h) * (Lij - rij);
2410
2411                         madd_v3_v3fl(pa1->state.co, Rij, -D * pa1->state.time * pa1->state.time * massfactor);
2412                         madd_v3_v3fl(pa2->state.co, Rij, D * pa2->state.time * pa2->state.time * massfactor);
2413                 }
2414         }
2415
2416         /* Loop through springs backwaqrds - for efficient delete function */
2417         for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
2418                 if(psys->fluid_springs[i].delete_flag)
2419                         delete_fluid_spring(psys, i);
2420         }
2421 }
2422
2423 /************************************************/
2424 /*                      Newtonian physics                                       */
2425 /************************************************/
2426 /* gathers all forces that effect particles and calculates a new state for the particle */
2427 static void apply_particle_forces(ParticleSimulationData *sim, int p, float dfra, float cfra)
2428 {
2429         ParticleSettings *part = sim->psys->part;
2430         ParticleData *pa = sim->psys->particles + p;
2431         EffectedPoint epoint;
2432         ParticleKey states[5], tkey;
2433         float timestep = psys_get_timestep(sim);
2434         float force[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
2435         float dtime=dfra*timestep, time, pa_mass=part->mass, fac /*, fra=sim->psys->cfra*/;
2436         int i, steps=1;
2437         ParticleTexture ptex;
2438
2439         psys_get_texture(sim, pa, &ptex, PAMAP_PHYSICS, cfra);
2440         
2441         /* maintain angular velocity */
2442         VECCOPY(pa->state.ave,pa->prev_state.ave);
2443         VECCOPY(oldpos,pa->state.co);
2444
2445         if(part->flag & PART_SIZEMASS)
2446                 pa_mass*=pa->size;
2447
2448         switch(part->integrator){
2449                 case PART_INT_EULER:
2450                         steps=1;
2451                         break;
2452                 case PART_INT_MIDPOINT:
2453                         steps=2;
2454                         break;
2455                 case PART_INT_RK4:
2456                         steps=4;
2457                         break;
2458                 case PART_INT_VERLET:
2459                         steps=1;
2460                         break;
2461         }
2462
2463         copy_particle_key(states,&pa->state,1);
2464
2465         for(i=0; i<steps; i++){
2466                 force[0]=force[1]=force[2]=0.0;
2467                 impulse[0]=impulse[1]=impulse[2]=0.0;
2468                 /* add effectors */
2469                 pd_point_from_particle(sim, pa, states+i, &epoint);
2470                 if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
2471                         pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
2472
2473                 mul_v3_fl(force, ptex.field);
2474                 mul_v3_fl(impulse, ptex.field);
2475
2476                 /* calculate air-particle interaction */
2477                 if(part->dragfac!=0.0f){
2478                         fac=-part->dragfac*pa->size*pa->size*len_v3(states[i].vel);
2479                         VECADDFAC(force,force,states[i].vel,fac);
2480                 }
2481
2482                 /* brownian force */
2483                 if(part->brownfac!=0.0){
2484                         force[0]+=(BLI_frand()-0.5f)*part->brownfac;
2485                         force[1]+=(BLI_frand()-0.5f)*part->brownfac;
2486                         force[2]+=(BLI_frand()-0.5f)*part->brownfac;
2487                 }
2488
2489                 /* force to acceleration*/
2490                 mul_v3_fl(force,1.0f/pa_mass);
2491
2492                 /* add global acceleration (gravitation) */
2493                 if(psys_uses_gravity(sim)
2494                         /* normal gravity is too strong for hair so it's disabled by default */
2495                         && (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
2496                         madd_v3_v3fl(force, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * ptex.gravity);
2497                 }
2498                 
2499                 /* calculate next state */
2500                 VECADD(states[i].vel,states[i].vel,impulse);
2501
2502                 switch(part->integrator){
2503                         case PART_INT_EULER:
2504                                 VECADDFAC(pa->state.co,states->co,states->vel,dtime);
2505                                 VECADDFAC(pa->state.vel,states->vel,force,dtime);
2506                                 break;
2507                         case PART_INT_MIDPOINT:
2508                                 if(i==0){
2509                                         VECADDFAC(states[1].co,states->co,states->vel,dtime*0.5f);
2510                                         VECADDFAC(states[1].vel,states->vel,force,dtime*0.5f);
2511                                         /*fra=sim->psys->cfra+0.5f*dfra;*/
2512                                 }
2513                                 else{
2514                                         VECADDFAC(pa->state.co,states->co,states[1].vel,dtime);
2515                                         VECADDFAC(pa->state.vel,states->vel,force,dtime);
2516                                 }
2517                                 break;
2518                         case PART_INT_RK4:
2519                                 switch(i){
2520                                         case 0:
2521                                                 VECCOPY(dx[0],states->vel);
2522                                                 mul_v3_fl(dx[0],dtime);
2523                                                 VECCOPY(dv[0],force);
2524                                                 mul_v3_fl(dv[0],dtime);
2525
2526                                                 VECADDFAC(states[1].co,states->co,dx[0],0.5f);
2527                                                 VECADDFAC(states[1].vel,states->vel,dv[0],0.5f);
2528                                                 /*fra=sim->psys->cfra+0.5f*dfra;*/
2529                                                 break;
2530                                         case 1:
2531                                                 VECADDFAC(dx[1],states->vel,dv[0],0.5f);
2532                                                 mul_v3_fl(dx[1],dtime);
2533                                                 VECCOPY(dv[1],force);
2534                                                 mul_v3_fl(dv[1],dtime);
2535
2536                                                 VECADDFAC(states[2].co,states->co,dx[1],0.5f);
2537                                                 VECADDFAC(states[2].vel,states->vel,dv[1],0.5f);
2538                                                 break;
2539                                         case 2:
2540                                                 VECADDFAC(dx[2],states->vel,dv[1],0.5f);
2541                                                 mul_v3_fl(dx[2],dtime);
2542                                                 VECCOPY(dv[2],force);
2543                                                 mul_v3_fl(dv[2],dtime);
2544
2545                                                 VECADD(states[3].co,states->co,dx[2]);
2546                                                 VECADD(states[3].vel,states->vel,dv[2]);
2547                                                 /*fra=cfra;*/
2548                                                 break;
2549                                         case 3:
2550                                                 VECADD(dx[3],states->vel,dv[2]);
2551                                                 mul_v3_fl(dx[3],dtime);
2552                                                 VECCOPY(dv[3],force);
2553                                                 mul_v3_fl(dv[3],dtime);
2554
2555                                                 VECADDFAC(pa->state.co,states->co,dx[0],1.0f/6.0f);
2556                                                 VECADDFAC(pa->state.co,pa->state.co,dx[1],1.0f/3.0f);
2557                                                 VECADDFAC(pa->state.co,pa->state.co,dx[2],1.0f/3.0f);
2558                                                 VECADDFAC(pa->state.co,pa->state.co,dx[3],1.0f/6.0f);
2559
2560                                                 VECADDFAC(pa->state.vel,states->vel,dv[0],1.0f/6.0f);
2561                                                 VECADDFAC(pa->state.vel,pa->state.vel,dv[1],1.0f/3.0f);
2562                                                 VECADDFAC(pa->state.vel,pa->state.vel,dv[2],1.0f/3.0f);
2563                                                 VECADDFAC(pa->state.vel,pa->state.vel,dv[3],1.0f/6.0f);
2564                                 }
2565                                 break;
2566                         case PART_INT_VERLET:   /* Verlet integration */
2567                                 VECADDFAC(pa->state.vel,pa->state.vel,force,dtime);
2568                                 VECADDFAC(pa->state.co,pa->state.co,pa->state.vel,dtime);
2569
2570                                 VECSUB(pa->state.vel,pa->state.co,oldpos);
2571                                 mul_v3_fl(pa->state.vel,1.0f/dtime);
2572                                 break;
2573                 }
2574         }
2575
2576         /* damp affects final velocity */
2577         if(part->dampfac != 0.f)
2578                 mul_v3_fl(pa->state.vel, 1.f - part->dampfac * ptex.damp);
2579
2580         VECCOPY(pa->state.ave, states->ave);
2581
2582         /* finally we do guides */
2583         time=(cfra-pa->time)/pa->lifetime;
2584         CLAMP(time,0.0,1.0);
2585
2586         VECCOPY(tkey.co,pa->state.co);
2587         VECCOPY(tkey.vel,pa->state.vel);
2588         tkey.time=pa->state.time;
2589
2590         if(part->type != PART_HAIR) {
2591                 if(do_guides(sim->psys->effectors, &tkey, p, time)) {
2592                         VECCOPY(pa->state.co,tkey.co);
2593                         /* guides don't produce valid velocity */
2594                         VECSUB(pa->state.vel,tkey.co,pa->prev_state.co);
2595                         mul_v3_fl(pa->state.vel,1.0f/dtime);
2596                         pa->state.time=tkey.time;
2597                 }
2598         }
2599 }
2600 static void rotate_particle(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2601 {
2602         float rotfac, rot1[4], rot2[4]={1.0,0.0,0.0,0.0}, dtime=dfra*timestep;
2603
2604         if((part->flag & PART_ROT_DYN)==0){
2605                 if(part->avemode==PART_AVE_SPIN){
2606                         float angle;
2607                         float len1 = len_v3(pa->prev_state.vel);
2608                         float len2 = len_v3(pa->state.vel);
2609
2610                         if(len1==0.0f || len2==0.0f)
2611                                 pa->state.ave[0]=pa->state.ave[1]=pa->state.ave[2]=0.0f;
2612                         else{
2613                                 cross_v3_v3v3(pa->state.ave,pa->prev_state.vel,pa->state.vel);
2614                                 normalize_v3(pa->state.ave);
2615                                 angle=dot_v3v3(pa->prev_state.vel,pa->state.vel)/(len1*len2);
2616                                 mul_v3_fl(pa->state.ave,saacos(angle)/dtime);
2617                         }
2618
2619                         axis_angle_to_quat(rot2,pa->state.vel,dtime*part->avefac);
2620                 }
2621         }
2622
2623         rotfac=len_v3(pa->state.ave);
2624         if(rotfac==0.0){ /* unit_qt(in VecRotToQuat) doesn't give unit quat [1,0,0,0]?? */
2625                 rot1[0]=1.0;
2626                 rot1[1]=rot1[2]=rot1[3]=0;
2627         }
2628         else{
2629                 axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
2630         }
2631         mul_qt_qtqt(pa->state.rot,rot1,pa->prev_state.rot);
2632         mul_qt_qtqt(pa->state.rot,rot2,pa->state.rot);
2633
2634         /* keep rotation quat in good health */
2635         normalize_qt(pa->state.rot);
2636 }
2637
2638 /* convert from triangle barycentric weights to quad mean value weights */
2639 static void intersect_dm_quad_weights(float *v1, float *v2, float *v3, float *v4, float *w)
2640 {
2641         float co[3], vert[4][3];
2642
2643         VECCOPY(vert[0], v1);
2644         VECCOPY(vert[1], v2);
2645         VECCOPY(vert[2], v3);
2646         VECCOPY(vert[3], v4);
2647
2648         co[0]= v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2] + v4[0]*w[3];
2649         co[1]= v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2] + v4[1]*w[3];
2650         co[2]= v1[2]*w[0] + v2[2]*w[1] + v3[2]*w[2] + v4[2]*w[3];
2651
2652         interp_weights_poly_v3( w,vert, 4, co);
2653 }
2654
2655 /* check intersection with a derivedmesh */
2656 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,
2657                                                   float *face_minmax, float *pa_minmax, float radius, float *ipoint)
2658 {
2659         MFace *mface=0;
2660         MVert *mvert=0;
2661         int i, totface, intersect=0;
2662         float cur_d, cur_uv[2], v1[3], v2[3], v3[3], v4[3], min[3], max[3], p_min[3],p_max[3];
2663         float cur_ipoint[3];
2664         
2665         if(dm==0){
2666                 psys_disable_all(ob);
2667
2668                 dm=mesh_get_derived_final(scene, ob, 0);
2669                 if(dm==0)
2670                         dm=mesh_get_derived_deform(scene, ob, 0);
2671
2672                 psys_enable_all(ob);
2673
2674                 if(dm==0)
2675                         return 0;
2676         }
2677
2678         
2679
2680         if(pa_minmax==0){
2681                 INIT_MINMAX(p_min,p_max);
2682                 DO_MINMAX(co1,p_min,p_max);
2683                 DO_MINMAX(co2,p_min,p_max);
2684         }
2685         else{
2686                 VECCOPY(p_min,pa_minmax);
2687                 VECCOPY(p_max,pa_minmax+3);
2688         }
2689
2690         totface=dm->getNumFaces(dm);
2691         mface=dm->getFaceDataArray(dm,CD_MFACE);
2692         mvert=dm->getVertDataArray(dm,CD_MVERT);
2693         
2694         /* lets intersect the faces */
2695         for(i=0; i<totface; i++,mface++){
2696                 if(vert_cos){
2697                         VECCOPY(v1,vert_cos+3*mface->v1);
2698                         VECCOPY(v2,vert_cos+3*mface->v2);
2699                         VECCOPY(v3,vert_cos+3*mface->v3);
2700                         if(mface->v4)
2701                                 VECCOPY(v4,vert_cos+3*mface->v4)
2702                 }
2703                 else{
2704                         VECCOPY(v1,mvert[mface->v1].co);
2705                         VECCOPY(v2,mvert[mface->v2].co);
2706                         VECCOPY(v3,mvert[mface->v3].co);
2707                         if(mface->v4)
2708                                 VECCOPY(v4,mvert[mface->v4].co)
2709                 }
2710
2711                 if(face_minmax==0){
2712                         INIT_MINMAX(min,max);
2713                         DO_MINMAX(v1,min,max);
2714                         DO_MINMAX(v2,min,max);
2715                         DO_MINMAX(v3,min,max);
2716                         if(mface->v4)
2717                                 DO_MINMAX(v4,min,max)
2718                         if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
2719                                 continue;
2720                 }
2721                 else{
2722                         VECCOPY(min, face_minmax+6*i);
2723                         VECCOPY(max, face_minmax+6*i+3);
2724                         if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
2725                                 continue;
2726                 }
2727
2728                 if(radius>0.0f){
2729                         if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v2, v3, v1, &cur_d, cur_ipoint)){
2730                                 if(cur_d<*min_d){
2731                                         *min_d=cur_d;
2732                                         VECCOPY(ipoint,cur_ipoint);
2733                                         *min_face=i;
2734                                         intersect=1;
2735                                 }
2736                         }
2737                         if(mface->v4){
2738                                 if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v4, v1, v3, &cur_d, cur_ipoint)){
2739                                         if(cur_d<*min_d){
2740                                                 *min_d=cur_d;
2741                                                 VECCOPY(ipoint,cur_ipoint);
2742                                                 *min_face=i;
2743                                                 intersect=1;
2744                                         }
2745                                 }
2746                         }
2747                 }
2748                 else{
2749                         if(isect_line_tri_v3(co1, co2, v1, v2, v3, &cur_d, cur_uv)){
2750                                 if(cur_d<*min_d){
2751                                         *min_d=cur_d;
2752                                         min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
2753                                         min_w[1]= cur_uv[0];
2754                                         min_w[2]= cur_uv[1];
2755                                         min_w[3]= 0.0f;
2756                                         if(mface->v4)
2757                                                 intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
2758                                         *min_face=i;
2759                                         intersect=1;
2760                                 }
2761                         }
2762                         if(mface->v4){
2763                                 if(isect_line_tri_v3(co1, co2, v1, v3, v4, &cur_d, cur_uv)){
2764                                         if(cur_d<*min_d){
2765                                                 *min_d=cur_d;
2766                                                 min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
2767                                                 min_w[1]= 0.0f;
2768                                                 min_w[2]= cur_uv[0];
2769                                                 min_w[3]= cur_uv[1];
2770                                                 intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
2771                                                 *min_face=i;
2772                                                 intersect=1;
2773                                         }
2774                                 }
2775                         }
2776                 }
2777         }
2778         return intersect;
2779 }
2780
2781 void particle_intersect_face(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
2782 {
2783         ParticleCollision *col = (ParticleCollision *) userdata;
2784         MFace *face = col->md->mfaces + index;
2785         MVert *x = col->md->x;
2786         MVert *v = col->md->current_v;
2787         float vel[3], co1[3], co2[3], uv[2], ipoint[3], temp[3], t;
2788         float x0[3], x1[3], x2[3], x3[3];
2789         float *t0=x0, *t1=x1, *t2=x2, *t3=(face->v4 ? x3 : NULL);
2790
2791         /* move collision face to start of timestep */
2792         madd_v3_v3v3fl(t0, x[face->v1].co, v[face->v1].co, col->cfra);
2793         madd_v3_v3v3fl(t1, x[face->v2].co, v[face->v2].co, col->cfra);
2794         madd_v3_v3v3fl(t2, x[face->v3].co, v[face->v3].co, col->cfra);
2795         if(t3)
2796                 madd_v3_v3v3fl(t3, x[face->v4].co, v[face->v4].co, col->cfra);
2797
2798         /* calculate average velocity of face */
2799         copy_v3_v3(vel, v[ face->v1 ].co);
2800         add_v3_v3(vel, v[ face->v2 ].co);
2801         add_v3_v3(vel, v[ face->v3 ].co);
2802         mul_v3_fl(vel, 0.33334f*col->dfra);
2803
2804         /* substract face velocity, in other words convert to 
2805            a coordinate system where only the particle moves */
2806         madd_v3_v3v3fl(co1, col->co1, vel, -col->f);
2807         sub_v3_v3v3(co2, col->co2, vel);
2808
2809         do
2810         {       
2811                 if(ray->radius == 0.0f) {
2812                         if(isect_line_tri_v3(co1, co2, t0, t1, t2, &t, uv)) {
2813                                 if(t >= 0.0f && t < hit->dist/col->ray_len) {
2814                                         hit->dist = col->ray_len * t;
2815                                         hit->index = index;
2816
2817                                         /* calculate normal that's facing the particle */
2818                                         normal_tri_v3( col->nor,t0, t1, t2);
2819                                         VECSUB(temp, co2, co1);
2820                                         if(dot_v3v3(col->nor, temp) > 0.0f)
2821                                                 negate_v3(col->nor);
2822
2823                                         VECCOPY(col->vel,vel);
2824
2825                                         col->hit_ob = col->ob;
2826                                         col->hit_md = col->md;
2827                                 }
2828                         }
2829                 }
2830                 else {
2831                         if(isect_sweeping_sphere_tri_v3(co1, co2, ray->radius, t0, t1, t2, &t, ipoint)) {
2832                                 if(t >=0.0f && t < hit->dist/col->ray_len) {
2833                                         hit->dist = col->ray_len * t;
2834                                         hit->index = index;
2835
2836                                         interp_v3_v3v3(temp, co1, co2, t);
2837                                         
2838                                         VECSUB(col->nor, temp, ipoint);
2839                                         normalize_v3(col->nor);
2840
2841                                         VECCOPY(col->vel,vel);
2842
2843                                         col->hit_ob = col->ob;
2844                                         col->hit_md = col->md;
2845                                 }
2846                         }
2847                 }
2848
2849                 t1 = t2;
2850                 t2 = t3;
2851                 t3 = NULL;
2852
2853         } while(t2);
2854 }
2855 /* Particle - Mesh collision code
2856  * Features:
2857  * - point and swept sphere to mesh surface collisions
2858  * - moving colliders (but not yet rotating or deforming colliders)
2859  * - friction & damping
2860  * - angular momentum <-> linear momentum
2861  * - high accuracy by re-applying particle acceleration after collision
2862  * - behaves relatively well even if limit of 10 collisions per simulation step is exceeded
2863  * Main parts:
2864  * 1. check for all possible deflectors for closest intersection on particle path
2865  * 2. if deflection was found calculate new coordinates or kill the particle
2866  */
2867 static void deflect_particle(ParticleSimulationData *sim, int p, float dfra, float cfra){
2868         Object *ground_ob = NULL;
2869         ParticleSettings *part = sim->psys->part;
2870         ParticleData *pa = sim->psys->particles + p;
2871         ParticleCollision col;
2872         ColliderCache *coll;