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