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