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