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