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