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