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