ac06f75737ff714f014a3eab2968997272c22364
[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= 0.0f;
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                         if(p == totpart-1 && weight[index[p]] == 0.0f)
1150                                 index[p]= index[p-1];
1151
1152                         jitoff[index[p]]= pos;
1153                 }
1154         }
1155
1156         MEM_freeN(sum);
1157
1158         /* weights are no longer used except for FROM_PARTICLE, which needs them zeroed for indexing */
1159         if(from==PART_FROM_PARTICLE){
1160                 for(i=0; i<tot; i++)
1161                         weight[i]=0.0f;
1162         }
1163
1164         /* 4. */
1165         if(distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
1166                 jitlevel= part->userjit;
1167                 
1168                 if(jitlevel == 0) {
1169                         jitlevel= totpart/tot;
1170                         if(part->flag & PART_EDISTR) jitlevel*= 2;      /* looks better in general, not very scietific */
1171                         if(jitlevel<3) jitlevel= 3;
1172                         //if(jitlevel>100) jitlevel= 100;
1173                 }
1174                 
1175                 jit= MEM_callocN(2+ jitlevel*2*sizeof(float), "jit");
1176
1177                 init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
1178                 BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
1179         }
1180
1181         /* 5. */
1182         ctx->tree= tree;
1183         ctx->seams= seams;
1184         ctx->totseam= totseam;
1185         ctx->psys= psys;
1186         ctx->index= index;
1187         ctx->jit= jit;
1188         ctx->jitlevel= jitlevel;
1189         ctx->jitoff= jitoff;
1190         ctx->weight= weight;
1191         ctx->maxweight= maxweight;
1192         ctx->from= (children)? PART_FROM_CHILD: from;
1193         ctx->cfrom= cfrom;
1194         ctx->distr= distr;
1195         ctx->dm= dm;
1196         ctx->tpars= tpars;
1197
1198         if(children) {
1199                 totpart= psys_render_simplify_distribution(ctx, totpart);
1200                 alloc_child_particles(psys, totpart);
1201         }
1202
1203         if(!children || psys->totchild < 10000)
1204                 totthread= 1;
1205         
1206         seed= 31415926 + ctx->psys->seed;
1207         for(i=0; i<totthread; i++) {
1208                 threads[i].rng= rng_new(seed);
1209                 threads[i].tot= totthread;
1210         }
1211
1212         return 1;
1213 }
1214
1215 static void distribute_particles_on_dm(DerivedMesh *finaldm, Object *ob, ParticleSystem *psys, int from)
1216 {
1217         ListBase threads;
1218         ParticleThread *pthreads;
1219         ParticleThreadContext *ctx;
1220         int i, totthread;
1221
1222         pthreads= psys_threads_create(ob, psys, G.scene->r.threads);
1223
1224         if(!psys_threads_init_distribution(pthreads, finaldm, from)) {
1225                 psys_threads_free(pthreads);
1226                 return;
1227         }
1228
1229         totthread= pthreads[0].tot;
1230         if(totthread > 1) {
1231                 BLI_init_threads(&threads, exec_distribution, totthread);
1232
1233                 for(i=0; i<totthread; i++)
1234                         BLI_insert_thread(&threads, &pthreads[i]);
1235
1236                 BLI_end_threads(&threads);
1237         }
1238         else
1239                 exec_distribution(&pthreads[0]);
1240
1241         if (from == PART_FROM_FACE)
1242                 psys_calc_dmfaces(ob, finaldm, psys);
1243
1244         ctx= pthreads[0].ctx;
1245         if(ctx->dm != finaldm)
1246                 ctx->dm->release(ctx->dm);
1247
1248         psys_threads_free(pthreads);
1249 }
1250
1251 /* ready for future use, to emit particles without geometry */
1252 static void distribute_particles_on_shape(Object *ob, ParticleSystem *psys, int from)
1253 {
1254         ParticleData *pa;
1255         int totpart=psys->totpart, p;
1256
1257         fprintf(stderr,"Shape emission not yet possible!\n");
1258
1259         for(p=0,pa=psys->particles; p<totpart; p++,pa++){
1260                 pa->fuv[0]=pa->fuv[1]=pa->fuv[2]=pa->fuv[3]= 0.0;
1261                 pa->foffset= 0.0f;
1262                 pa->num= -1;
1263         }
1264 }
1265 static void distribute_particles(Object *ob, ParticleSystem *psys, int from)
1266 {
1267         ParticleSystemModifierData *psmd=0;
1268         int distr_error=0;
1269         psmd=psys_get_modifier(ob,psys);
1270
1271         if(psmd){
1272                 if(psmd->dm)
1273                         distribute_particles_on_dm(psmd->dm,ob,psys,from);
1274                 else
1275                         distr_error=1;
1276         }
1277         else
1278                 distribute_particles_on_shape(ob,psys,from);
1279
1280         if(distr_error){
1281                 ParticleData *pa;
1282                 int totpart=psys->totpart, p;
1283
1284                 fprintf(stderr,"Particle distribution error!\n");
1285
1286                 for(p=0,pa=psys->particles; p<totpart; p++,pa++){
1287                         pa->fuv[0]=pa->fuv[1]=pa->fuv[2]=pa->fuv[3]= 0.0;
1288                         pa->foffset= 0.0f;
1289                         pa->num= -1;
1290                 }
1291         }
1292 }
1293
1294 /* threaded child particle distribution and path caching */
1295 ParticleThread *psys_threads_create(struct Object *ob, struct ParticleSystem *psys, int totthread)
1296 {
1297         ParticleThread *threads;
1298         ParticleThreadContext *ctx;
1299         int i;
1300         
1301         threads= MEM_callocN(sizeof(ParticleThread)*totthread, "ParticleThread");
1302         ctx= MEM_callocN(sizeof(ParticleThreadContext), "ParticleThreadContext");
1303
1304         ctx->ob= ob;
1305         ctx->psys= psys;
1306         ctx->psmd= psys_get_modifier(ob, psys);
1307         ctx->dm= ctx->psmd->dm;
1308         ctx->ma= give_current_material(ob, psys->part->omat);
1309
1310         memset(threads, 0, sizeof(ParticleThread)*totthread);
1311
1312         for(i=0; i<totthread; i++) {
1313                 threads[i].ctx= ctx;
1314                 threads[i].num= i;
1315                 threads[i].tot= totthread;
1316         }
1317
1318         return threads;
1319 }
1320
1321 void psys_threads_free(ParticleThread *threads)
1322 {
1323         ParticleThreadContext *ctx= threads[0].ctx;
1324         int i, totthread= threads[0].tot;
1325
1326         /* path caching */
1327         if(ctx->vg_length)
1328                 MEM_freeN(ctx->vg_length);
1329         if(ctx->vg_clump)
1330                 MEM_freeN(ctx->vg_clump);
1331         if(ctx->vg_kink)
1332                 MEM_freeN(ctx->vg_kink);
1333         if(ctx->vg_rough1)
1334                 MEM_freeN(ctx->vg_rough1);
1335         if(ctx->vg_rough2)
1336                 MEM_freeN(ctx->vg_roughe);
1337         if(ctx->vg_roughe)
1338                 MEM_freeN(ctx->vg_roughe);
1339
1340         if(ctx->psys->lattice){
1341                 end_latt_deform();
1342                 ctx->psys->lattice=0;
1343         }
1344
1345         /* distribution */
1346         if(ctx->jit) MEM_freeN(ctx->jit);
1347         if(ctx->jitoff) MEM_freeN(ctx->jitoff);
1348         if(ctx->weight) MEM_freeN(ctx->weight);
1349         if(ctx->index) MEM_freeN(ctx->index);
1350         if(ctx->skip) MEM_freeN(ctx->skip);
1351         if(ctx->seams) MEM_freeN(ctx->seams);
1352         //if(ctx->vertpart) MEM_freeN(ctx->vertpart);
1353         BLI_kdtree_free(ctx->tree);
1354
1355         /* threads */
1356         for(i=0; i<totthread; i++) {
1357                 if(threads[i].rng)
1358                         rng_free(threads[i].rng);
1359                 if(threads[i].rng_path)
1360                         rng_free(threads[i].rng_path);
1361         }
1362
1363         MEM_freeN(ctx);
1364         MEM_freeN(threads);
1365 }
1366
1367 /* set particle parameters that don't change during particle's life */
1368 void initialize_particle(ParticleData *pa, int p, Object *ob, ParticleSystem *psys, ParticleSystemModifierData *psmd)
1369 {
1370         ParticleSettings *part;
1371         ParticleTexture ptex;
1372         Material *ma=0;
1373         IpoCurve *icu=0;
1374         int totpart;
1375         float rand,length;
1376
1377         part=psys->part;
1378
1379         totpart=psys->totpart;
1380
1381         ptex.life=ptex.size=ptex.exist=ptex.length=1.0;
1382         ptex.time=(float)p/(float)totpart;
1383
1384         BLI_srandom(psys->seed+p);
1385
1386         if(part->from!=PART_FROM_PARTICLE){
1387                 ma=give_current_material(ob,part->omat);
1388
1389                 /* TODO: needs some work to make most blendtypes generally usefull */
1390                 psys_get_texture(ob,ma,psmd,psys,pa,&ptex,MAP_PA_INIT);
1391         }
1392         
1393         pa->lifetime= part->lifetime*ptex.life;
1394
1395         if(part->type==PART_HAIR)
1396                 pa->time=0.0f;
1397         else if(part->type==PART_REACTOR && (part->flag&PART_REACT_STA_END)==0)
1398                 pa->time=MAXFRAMEF;
1399         else{
1400                 //icu=find_ipocurve(psys->part->ipo,PART_EMIT_TIME);
1401                 //if(icu){
1402                 //      calc_icu(icu,100*ptex.time);
1403                 //      ptex.time=icu->curval;
1404                 //}
1405
1406                 pa->time= part->sta + (part->end - part->sta)*ptex.time;
1407         }
1408
1409
1410         if(part->type==PART_HAIR){
1411                 pa->lifetime=100.0f;
1412         }
1413         else{
1414                 icu=find_ipocurve(psys->part->ipo,PART_EMIT_LIFE);
1415                 if(icu){
1416                         calc_icu(icu,100*ptex.time);
1417                         pa->lifetime*=icu->curval;
1418                 }
1419
1420         /* need to get every rand even if we don't use them so that randoms don't affect eachother */
1421                 rand= BLI_frand();
1422                 if(part->randlife!=0.0)
1423                         pa->lifetime*= 1.0f - part->randlife*rand;
1424         }
1425
1426         pa->dietime= pa->time+pa->lifetime;
1427
1428         pa->sizemul= BLI_frand();
1429
1430         rand= BLI_frand();
1431
1432         /* while loops are to have a spherical distribution (avoid cubic distribution) */
1433         length=2.0f;
1434         while(length>1.0){
1435                 pa->r_ve[0]=2.0f*(BLI_frand()-0.5f);
1436                 pa->r_ve[1]=2.0f*(BLI_frand()-0.5f);
1437                 pa->r_ve[2]=2.0f*(BLI_frand()-0.5f);
1438                 length=VecLength(pa->r_ve);
1439         }
1440
1441         length=2.0f;
1442         while(length>1.0){
1443                 pa->r_ave[0]=2.0f*(BLI_frand()-0.5f);
1444                 pa->r_ave[1]=2.0f*(BLI_frand()-0.5f);
1445                 pa->r_ave[2]=2.0f*(BLI_frand()-0.5f);
1446                 length=VecLength(pa->r_ave);
1447         }
1448
1449         pa->r_rot[0]=2.0f*(BLI_frand()-0.5f);
1450         pa->r_rot[1]=2.0f*(BLI_frand()-0.5f);
1451         pa->r_rot[2]=2.0f*(BLI_frand()-0.5f);
1452         pa->r_rot[3]=2.0f*(BLI_frand()-0.5f);
1453
1454         NormalQuat(pa->r_rot);
1455
1456         if(part->distr!=PART_DISTR_GRID){
1457                 /* any unique random number will do (r_ave[0]) */
1458                 if(ptex.exist < 0.5*(1.0+pa->r_ave[0]))
1459                         pa->flag |= PARS_UNEXIST;
1460                 else
1461                         pa->flag &= ~PARS_UNEXIST;
1462         }
1463
1464         pa->loop=0;
1465         /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
1466         /* usage other than straight after distribute has to handle this index by itself - jahka*/
1467         //pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we dont have a derived mesh face */
1468 }
1469 static void initialize_all_particles(Object *ob, ParticleSystem *psys, ParticleSystemModifierData *psmd)
1470 {
1471         IpoCurve *icu=0;
1472         ParticleData *pa;
1473         int p, totpart=psys->totpart;
1474
1475         for(p=0, pa=psys->particles; p<totpart; p++, pa++)
1476                 initialize_particle(pa,p,ob,psys,psmd);
1477         
1478         /* store the derived mesh face index for each particle */
1479         icu=find_ipocurve(psys->part->ipo,PART_EMIT_FREQ);
1480         if(icu){
1481                 float time=psys->part->sta, end=psys->part->end;
1482                 float v1, v2, a=0.0f, t1,t2, d;
1483
1484                 p=0;
1485                 pa=psys->particles;
1486
1487                 calc_icu(icu,time);
1488                 v1=icu->curval;
1489                 if(v1<0.0f) v1=0.0f;
1490
1491                 calc_icu(icu,time+1.0f);
1492                 v2=icu->curval;
1493                 if(v2<0.0f) v2=0.0f;
1494
1495                 for(p=0, pa=psys->particles; p<totpart && time<end; p++, pa++){
1496                         while(a+0.5f*(v1+v2) < (float)(p+1) && time<end){
1497                                 a+=0.5f*(v1+v2);
1498                                 v1=v2;
1499                                 time++;
1500                                 calc_icu(icu,time+1.0f);
1501                                 v2=icu->curval;
1502                         }
1503                         if(time<end){
1504                                 if(v1==v2){
1505                                         pa->time=time+((float)(p+1)-a)/v1;
1506                                 }
1507                                 else{
1508                                         d=(float)sqrt(v1*v1-2.0f*(v2-v1)*(a-(float)(p+1)));
1509                                         t1=(-v1+d)/(v2-v1);
1510                                         t2=(-v1-d)/(v2-v1);
1511
1512                                         /* the root between 0-1 is the correct one */
1513                                         if(t1>0.0f && t1<=1.0f)
1514                                                 pa->time=time+t1;
1515                                         else
1516                                                 pa->time=time+t2;
1517                                 }
1518                         }
1519
1520                         pa->dietime = pa->time+pa->lifetime;
1521                         pa->flag &= ~PARS_UNEXIST;
1522                 }
1523                 for(; p<totpart; p++, pa++){
1524                         pa->flag |= PARS_UNEXIST;
1525                 }
1526         }
1527 }
1528 /* sets particle to the emitter surface with initial velocity & rotation */
1529 void reset_particle(ParticleData *pa, ParticleSystem *psys, ParticleSystemModifierData *psmd, Object *ob,
1530                                         float dtime, float cfra, float *vg_vel, float *vg_tan, float *vg_rot)
1531 {
1532         ParticleSettings *part;
1533         ParticleTexture ptex;
1534         ParticleKey state;
1535         IpoCurve *icu=0;
1536         float fac, rotfac, phasefac, nor[3]={0,0,0},loc[3],tloc[3],vel[3]={0.0,0.0,0.0},rot[4],*q2=0;
1537         float r_vel[3],r_ave[3],r_rot[4],p_vel[3]={0.0,0.0,0.0};
1538         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};
1539         float q_one[4]={1.0,0.0,0.0,0.0}, q_phase[4];
1540         part=psys->part;
1541
1542         ptex.ivel=1.0;
1543         
1544         if(part->from==PART_FROM_PARTICLE){
1545                 Object *tob;
1546                 ParticleSystem *tpsys=0;
1547                 float speed;
1548
1549                 tob=psys->target_ob;
1550                 if(tob==0)
1551                         tob=ob;
1552
1553                 tpsys=BLI_findlink(&tob->particlesystem,psys->target_psys-1);
1554                 
1555                 /*TODO: get precise location of particle at birth*/
1556
1557                 state.time=cfra;
1558                 psys_get_particle_state(tob,tpsys,pa->num,&state,1);
1559                 psys_get_from_key(&state,loc,nor,rot,0);
1560
1561                 QuatMulVecf(rot,vtan);
1562                 QuatMulVecf(rot,utan);
1563                 VECCOPY(r_vel,pa->r_ve);
1564                 VECCOPY(r_rot,pa->r_rot);
1565                 VECCOPY(r_ave,pa->r_ave);
1566
1567                 VECCOPY(p_vel,state.vel);
1568                 speed=Normalize(p_vel);
1569                 VecMulf(p_vel,Inpf(pa->r_ve,p_vel));
1570                 VECSUB(p_vel,pa->r_ve,p_vel);
1571                 Normalize(p_vel);
1572                 VecMulf(p_vel,speed);
1573         }
1574         else{
1575                 /* get precise emitter matrix if particle is born */
1576                 if(part->type!=PART_HAIR && pa->time < cfra && pa->time >= psys->cfra)
1577                         where_is_object_time(ob,pa->time);
1578
1579                 /* get birth location from object               */
1580                 psys_particle_on_emitter(ob,psmd,part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0,0);
1581                 
1582                 /* save local coordinates for later             */
1583                 VECCOPY(tloc,loc);
1584                 
1585                 /* get possible textural influence */
1586                 psys_get_texture(ob,give_current_material(ob,part->omat),psmd,psys,pa,&ptex,MAP_PA_IVEL);
1587
1588                 if(vg_vel){
1589                         ptex.ivel*=psys_interpolate_value_from_verts(psmd->dm,part->from,pa->num,pa->fuv,vg_vel);
1590                 }
1591
1592                 /* particles live in global space so    */
1593                 /* let's convert:                                               */
1594                 /* -location                                                    */
1595                 Mat4MulVecfl(ob->obmat,loc);
1596                 
1597                 /* -normal                                                              */
1598                 VECADD(nor,tloc,nor);
1599                 Mat4MulVecfl(ob->obmat,nor);
1600                 VECSUB(nor,nor,loc);
1601                 Normalize(nor);
1602
1603                 /* -tangent                                                             */
1604                 if(part->tanfac!=0.0){
1605                         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;
1606                         VecMulf(vtan,-(float)cos(M_PI*(part->tanphase+phase)));
1607                         fac=-(float)sin(M_PI*(part->tanphase+phase));
1608                         VECADDFAC(vtan,vtan,utan,fac);
1609
1610                         VECADD(vtan,tloc,vtan);
1611                         Mat4MulVecfl(ob->obmat,vtan);
1612                         VECSUB(vtan,vtan,loc);
1613
1614                         VECCOPY(utan,nor);
1615                         VecMulf(utan,Inpf(vtan,nor));
1616                         VECSUB(vtan,vtan,utan);
1617                         
1618                         Normalize(vtan);
1619                 }
1620                 
1621
1622                 /* -velocity                                                    */
1623                 if(part->randfac!=0.0){
1624                         VECADD(r_vel,tloc,pa->r_ve);
1625                         Mat4MulVecfl(ob->obmat,r_vel);
1626                         VECSUB(r_vel,r_vel,loc);
1627                         Normalize(r_vel);
1628                 }
1629
1630                 /* -angular velocity                                    */
1631                 if(part->avemode==PART_AVE_RAND){
1632                         VECADD(r_ave,tloc,pa->r_ave);
1633                         Mat4MulVecfl(ob->obmat,r_ave);
1634                         VECSUB(r_ave,r_ave,loc);
1635                         Normalize(r_ave);
1636                 }
1637                 
1638                 /* -rotation                                                    */
1639                 if(part->randrotfac != 0.0f){
1640                         QUATCOPY(r_rot,pa->r_rot);
1641                         Mat4ToQuat(ob->obmat,rot);
1642                         QuatMul(r_rot,r_rot,rot);
1643                 }
1644         }
1645         /* conversion done so now we apply new: */
1646         /* -velocity from:                                              */
1647         /*              *emitter velocity                               */
1648         if(dtime!=0.0 && part->obfac!=0.0){
1649                 VECSUB(vel,loc,pa->state.co);
1650                 VecMulf(vel,part->obfac/dtime);
1651         }
1652         
1653         /*              *emitter normal                                 */
1654         if(part->normfac!=0.0)
1655                 VECADDFAC(vel,vel,nor,part->normfac);
1656         
1657         /*              *emitter tangent                                */
1658         if(part->tanfac!=0.0)
1659                 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));
1660
1661         /*              *texture                                                */
1662         /* TODO */
1663
1664         /*              *random                                                 */
1665         if(part->randfac!=0.0)
1666                 VECADDFAC(vel,vel,r_vel,part->randfac);
1667
1668         /*              *particle                                               */
1669         if(part->partfac!=0.0)
1670                 VECADDFAC(vel,vel,p_vel,part->partfac);
1671
1672         icu=find_ipocurve(psys->part->ipo,PART_EMIT_VEL);
1673         if(icu){
1674                 calc_icu(icu,100*((pa->time-part->sta)/(part->end-part->sta)));
1675                 ptex.ivel*=icu->curval;
1676         }
1677
1678         VecMulf(vel,ptex.ivel);
1679         
1680         VECCOPY(pa->state.vel,vel);
1681
1682         /* -location from emitter                               */
1683         VECCOPY(pa->state.co,loc);
1684
1685         /* -rotation                                                    */
1686         pa->state.rot[0]=1.0;
1687         pa->state.rot[1]=pa->state.rot[2]=pa->state.rot[3]=0.0;
1688
1689         if(part->rotmode){
1690                 /* create vector into which rotation is aligned */
1691                 switch(part->rotmode){
1692                         case PART_ROT_NOR:
1693                                 VecCopyf(rot_vec, nor);
1694                                 break;
1695                         case PART_ROT_VEL:
1696                                 VecCopyf(rot_vec, vel);
1697                                 break;
1698                         case PART_ROT_GLOB_X:
1699                         case PART_ROT_GLOB_Y:
1700                         case PART_ROT_GLOB_Z:
1701                                 rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
1702                                 break;
1703                         case PART_ROT_OB_X:
1704                         case PART_ROT_OB_Y:
1705                         case PART_ROT_OB_Z:
1706                                 VecCopyf(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
1707                                 break;
1708                 }
1709                 
1710                 /* create rotation quat */
1711                 VecMulf(rot_vec,-1.0);
1712                 q2= vectoquat(rot_vec, OB_POSX, OB_POSZ);
1713
1714                 /* randomize rotation quat */
1715                 if(part->randrotfac!=0.0f)
1716                         QuatInterpol(rot, q2, r_rot, part->randrotfac);
1717                 else
1718                         QuatCopy(rot,q2);
1719
1720                 /* rotation phase */
1721                 phasefac = part->phasefac;
1722                 if(part->randphasefac != 0.0f) /* abuse r_ave[0] as a random number */
1723                         phasefac += part->randphasefac * pa->r_ave[0];
1724                 VecRotToQuat(x_vec, phasefac*(float)M_PI, q_phase);
1725
1726                 /* combine base rotation & phase */
1727                 QuatMul(pa->state.rot, rot, q_phase);
1728         }
1729
1730         /* -angular velocity                                    */
1731
1732         pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0;
1733
1734         if(part->avemode){
1735                 switch(part->avemode){
1736                         case PART_AVE_SPIN:
1737                                 VECCOPY(pa->state.ave,vel);
1738                                 break;
1739                         case PART_AVE_RAND:
1740                                 VECCOPY(pa->state.ave,r_ave);
1741                                 break;
1742                 }
1743                 Normalize(pa->state.ave);
1744                 VecMulf(pa->state.ave,part->avefac);
1745
1746                 icu=find_ipocurve(psys->part->ipo,PART_EMIT_AVE);
1747                 if(icu){
1748                         calc_icu(icu,100*((pa->time-part->sta)/(part->end-part->sta)));
1749                         VecMulf(pa->state.ave,icu->curval);
1750                 }
1751         }
1752
1753         pa->dietime = pa->time + pa->lifetime;
1754
1755         if(pa->time >= cfra)
1756                 pa->alive = PARS_UNBORN;
1757
1758         pa->state.time = cfra;
1759
1760         pa->stick_ob = 0;
1761         pa->flag &= ~PARS_STICKY;
1762 }
1763 static void reset_all_particles(Object *ob, ParticleSystem *psys, ParticleSystemModifierData *psmd, float dtime, float cfra, int from)
1764 {
1765         ParticleData *pa;
1766         int p, totpart=psys->totpart;
1767         float *vg_vel=psys_cache_vgroup(psmd->dm,psys,PSYS_VG_VEL);
1768         float *vg_tan=psys_cache_vgroup(psmd->dm,psys,PSYS_VG_TAN);
1769         float *vg_rot=psys_cache_vgroup(psmd->dm,psys,PSYS_VG_ROT);
1770         
1771         for(p=from, pa=psys->particles+from; p<totpart; p++, pa++)
1772                 reset_particle(pa, psys, psmd, ob, dtime, cfra, vg_vel, vg_tan, vg_rot);
1773
1774         if(vg_vel)
1775                 MEM_freeN(vg_vel);
1776 }
1777 /************************************************/
1778 /*                      Keyed particles                                         */
1779 /************************************************/
1780 /* a bit of an unintuitive function :) counts objects in a keyed chain and returns 1 if some of them were selected (used in drawing) */
1781 int psys_count_keyed_targets(Object *ob, ParticleSystem *psys)
1782 {
1783         ParticleSystem *kpsys=psys,*tpsys;
1784         ParticleSettings *tpart;
1785         Object *kob=ob,*tob;
1786         int select=ob->flag&SELECT;
1787         short totkeyed=0;
1788         Base *base;
1789
1790         ListBase lb;
1791         lb.first=lb.last=0;
1792
1793         tob=psys->keyed_ob;
1794         while(tob){
1795                 if((tpsys=BLI_findlink(&tob->particlesystem,kpsys->keyed_psys-1))){
1796                         tpart=tpsys->part;
1797
1798                         if(tpart->phystype==PART_PHYS_KEYED){
1799                                 if(lb.first){
1800                                         for(base=lb.first;base;base=base->next){
1801                                                 if(tob==base->object){
1802                                                         fprintf(stderr,"Error: loop in keyed chain!\n");
1803                                                         BLI_freelistN(&lb);
1804                                                         return select;
1805                                                 }
1806                                         }
1807                                 }
1808                                 base=MEM_callocN(sizeof(Base), "keyed base");
1809                                 base->object=tob;
1810                                 BLI_addtail(&lb,base);
1811
1812                                 if(tob->flag&SELECT)
1813                                         select++;
1814                                 kob=tob;
1815                                 kpsys=tpsys;
1816                                 tob=tpsys->keyed_ob;
1817                                 totkeyed++;
1818                         }
1819                         else{
1820                                 tob=0;
1821                                 totkeyed++;
1822                         }
1823                 }
1824                 else
1825                         tob=0;
1826         }
1827         psys->totkeyed=totkeyed;
1828         BLI_freelistN(&lb);
1829         return select;
1830 }
1831 void set_keyed_keys(Object *ob, ParticleSystem *psys)
1832 {
1833         Object *kob = ob;
1834         ParticleSystem *kpsys = psys;
1835         ParticleData *pa;
1836         ParticleKey state;
1837         int totpart = psys->totpart, i, k, totkeys = psys->totkeyed + 1;
1838         float prevtime, nexttime, keyedtime;
1839
1840         /* no proper targets so let's clear and bail out */
1841         if(psys->totkeyed==0){
1842                 free_keyed_keys(psys);
1843                 psys->flag &= ~PSYS_KEYED;
1844                 return;
1845         }
1846
1847         if(totpart && psys->particles->totkey != totkeys){
1848                 free_keyed_keys(psys);
1849                 
1850                 psys->particles->keys = MEM_callocN(psys->totpart * totkeys * sizeof(ParticleKey),"Keyed keys");
1851
1852                 psys->particles->totkey = totkeys;
1853                 
1854                 for(i=1, pa=psys->particles+1; i<totpart; i++,pa++){
1855                         pa->keys = (pa-1)->keys + totkeys;
1856                         pa->totkey = totkeys;
1857                 }
1858         }
1859         
1860         psys->flag &= ~PSYS_KEYED;
1861         state.time=-1.0;
1862
1863         for(k=0; k<totkeys; k++){
1864                 for(i=0,pa=psys->particles; i<totpart; i++, pa++){
1865                         psys_get_particle_state(kob, kpsys, i%kpsys->totpart, pa->keys + k, 1);
1866
1867                         if(k==0)
1868                                 pa->keys->time = pa->time;
1869                         else if(k==totkeys-1)
1870                                 (pa->keys + k)->time = pa->time + pa->lifetime;
1871                         else{
1872                                 if(psys->flag & PSYS_KEYED_TIME){
1873                                         prevtime = (pa->keys + k - 1)->time;
1874                                         nexttime = pa->time + pa->lifetime;
1875                                         keyedtime = kpsys->part->keyed_time;
1876                                         (pa->keys + k)->time = (1.0f - keyedtime) * prevtime + keyedtime * nexttime;
1877                                 }
1878                                 else
1879                                         (pa->keys+k)->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
1880                         }
1881                 }
1882                 if(kpsys->keyed_ob){
1883                         kob = kpsys->keyed_ob;
1884                         kpsys = BLI_findlink(&kob->particlesystem, kpsys->keyed_psys - 1);
1885                 }
1886         }
1887
1888         psys->flag |= PSYS_KEYED;
1889 }
1890 /************************************************/
1891 /*                      Reactors                                                        */
1892 /************************************************/
1893 static void push_reaction(Object* ob, ParticleSystem *psys, int pa_num, int event, ParticleKey *state)
1894 {
1895         Object *rob;
1896         ParticleSystem *rpsys;
1897         ParticleSettings *rpart;
1898         ParticleData *pa;
1899         ListBase *lb=&psys->effectors;
1900         ParticleEffectorCache *ec;
1901         ParticleReactEvent *re;
1902
1903         if(lb->first) for(ec = lb->first; ec; ec= ec->next){
1904                 if(ec->type & PSYS_EC_REACTOR){
1905                         /* all validity checks already done in add_to_effectors */
1906                         rob=ec->ob;
1907                         rpsys=BLI_findlink(&rob->particlesystem,ec->psys_nbr);
1908                         rpart=rpsys->part;
1909                         if(rpsys->part->reactevent==event){
1910                                 pa=psys->particles+pa_num;
1911                                 re= MEM_callocN(sizeof(ParticleReactEvent), "react event");
1912                                 re->event=event;
1913                                 re->pa_num = pa_num;
1914                                 re->ob = ob;
1915                                 re->psys = psys;
1916                                 re->size = pa->size;
1917                                 copy_particle_key(&re->state,state,1);
1918
1919                                 switch(event){
1920                                         case PART_EVENT_DEATH:
1921                                                 re->time=pa->dietime;
1922                                                 break;
1923                                         case PART_EVENT_COLLIDE:
1924                                                 re->time=state->time;
1925                                                 break;
1926                                         case PART_EVENT_NEAR:
1927                                                 re->time=state->time;
1928                                                 break;
1929                                 }
1930
1931                                 BLI_addtail(&rpsys->reactevents, re);
1932                         }
1933                 }
1934         }
1935 }
1936 static void react_to_events(ParticleSystem *psys, int pa_num)
1937 {
1938         ParticleSettings *part=psys->part;
1939         ParticleData *pa=psys->particles+pa_num;
1940         ParticleReactEvent *re=psys->reactevents.first;
1941         int birth=0;
1942         float dist=0.0f;
1943
1944         for(re=psys->reactevents.first; re; re=re->next){
1945                 birth=0;
1946                 if(part->from==PART_FROM_PARTICLE){
1947                         if(pa->num==re->pa_num){
1948                                 if(re->event==PART_EVENT_NEAR){
1949                                         ParticleData *tpa = re->psys->particles+re->pa_num;
1950                                         float pa_time=tpa->time + pa->foffset*tpa->lifetime;
1951                                         if(re->time > pa_time){
1952                                                 pa->alive=PARS_ALIVE;
1953                                                 pa->time=pa_time;
1954                                                 pa->dietime=pa->time+pa->lifetime;
1955                                         }
1956                                 }
1957                                 else{
1958                                         if(pa->alive==PARS_UNBORN){
1959                                                 pa->alive=PARS_ALIVE;
1960                                                 pa->time=re->time;
1961                                                 pa->dietime=pa->time+pa->lifetime;
1962                                         }
1963                                 }
1964                         }
1965                 }
1966                 else{
1967                         dist=VecLenf(pa->state.co, re->state.co);
1968                         if(dist <= re->size){
1969                                 if(pa->alive==PARS_UNBORN){
1970                                         pa->alive=PARS_ALIVE;
1971                                         pa->time=re->time;
1972                                         pa->dietime=pa->time+pa->lifetime;
1973                                         birth=1;
1974                                 }
1975                                 if(birth || part->flag&PART_REACT_MULTIPLE){
1976                                         float vec[3];
1977                                         VECSUB(vec,pa->state.co, re->state.co);
1978                                         if(birth==0)
1979                                                 VecMulf(vec,(float)pow(1.0f-dist/re->size,part->reactshape));
1980                                         VECADDFAC(pa->state.vel,pa->state.vel,vec,part->reactfac);
1981                                         VECADDFAC(pa->state.vel,pa->state.vel,re->state.vel,part->partfac);
1982                                 }
1983                                 if(birth)
1984                                         VecMulf(pa->state.vel,(float)pow(1.0f-dist/re->size,part->reactshape));
1985                         }
1986                 }
1987         }
1988 }
1989 void psys_get_reactor_target(Object *ob, ParticleSystem *psys, Object **target_ob, ParticleSystem **target_psys)
1990 {
1991         Object *tob;
1992
1993         tob=psys->target_ob;
1994         if(tob==0)
1995                 tob=ob;
1996         
1997         *target_psys=BLI_findlink(&tob->particlesystem,psys->target_psys-1);
1998         if(*target_psys)
1999                 *target_ob=tob;
2000         else
2001                 *target_ob=0;
2002 }
2003 /************************************************/
2004 /*                      Point Cache                                                     */
2005 /************************************************/
2006 void clear_particles_from_cache(Object *ob, ParticleSystem *psys, int cfra)
2007 {
2008         ParticleSystemModifierData *psmd = psys_get_modifier(ob,psys);
2009         int stack_index = modifiers_indexInObject(ob,(ModifierData*)psmd);
2010
2011         BKE_ptcache_id_clear((ID *)ob, PTCACHE_CLEAR_ALL, cfra, stack_index);
2012 }
2013 static void write_particles_to_cache(Object *ob, ParticleSystem *psys, int cfra)
2014 {
2015         FILE *fp = NULL;
2016         ParticleSystemModifierData *psmd = psys_get_modifier(ob,psys);
2017         ParticleData *pa;
2018         int stack_index = modifiers_indexInObject(ob,(ModifierData*)psmd);
2019         int i, totpart = psys->totpart;
2020
2021         if(totpart == 0) return;
2022
2023         fp = BKE_ptcache_id_fopen((ID *)ob, 'w', cfra, stack_index);
2024         if(!fp) return;
2025
2026         for(i=0, pa=psys->particles; i<totpart; i++, pa++)
2027                 fwrite(&pa->state, sizeof(ParticleKey), 1, fp);
2028         
2029         fclose(fp);
2030 }
2031 static int get_particles_from_cache(Object *ob, ParticleSystem *psys, int cfra)
2032 {
2033         FILE *fp = NULL;
2034         ParticleSystemModifierData *psmd = psys_get_modifier(ob,psys);
2035         ParticleData *pa;
2036         int stack_index = modifiers_indexInObject(ob,(ModifierData*)psmd);
2037         int i, totpart = psys->totpart, ret = 1;
2038
2039         if(totpart == 0) return 0;
2040
2041         fp = BKE_ptcache_id_fopen((ID *)ob, 'r', cfra, stack_index);
2042         if(!fp)
2043                 ret = 0;
2044         else {
2045                 for(i=0, pa=psys->particles; i<totpart; i++, pa++)
2046                         if((fread(&pa->state, sizeof(ParticleKey), 1, fp)) != 1) {
2047                                 ret = 0;
2048                                 break;
2049                         }
2050                 
2051                 fclose(fp);
2052         }
2053
2054         return ret;
2055 }
2056 /************************************************/
2057 /*                      Effectors                                                       */
2058 /************************************************/
2059 static float effector_falloff(PartDeflect *pd, float *eff_velocity, float *vec_to_part)
2060 {
2061         float eff_dir[3], temp[3];
2062         float falloff=1.0, fac, r_fac;
2063         
2064         VecCopyf(eff_dir,eff_velocity);
2065         Normalize(eff_dir);
2066
2067         if(pd->flag & PFIELD_POSZ && Inpf(eff_dir,vec_to_part)<0.0f)
2068                 falloff=0.0f;
2069         else switch(pd->falloff){
2070                 case PFIELD_FALL_SPHERE:
2071                         fac=VecLength(vec_to_part);
2072                         if(pd->flag&PFIELD_USEMAX && fac>pd->maxdist){
2073                                 falloff=0.0f;
2074                                 break;
2075                         }
2076
2077                         if(pd->flag & PFIELD_USEMIN){
2078                                 if(fac>pd->mindist)
2079                                         fac+=1.0f-pd->mindist;
2080                                 else
2081                                         fac=1.0f;
2082                         }
2083                         else if(fac<0.001)
2084                                 fac=0.001f;
2085
2086                         falloff=1.0f/(float)pow((double)fac,(double)pd->f_power);
2087                         break;
2088
2089                 case PFIELD_FALL_TUBE:
2090                         fac=Inpf(vec_to_part,eff_dir);
2091                         if(pd->flag&PFIELD_USEMAX && ABS(fac)>pd->maxdist){
2092                                 falloff=0.0f;
2093                                 break;
2094                         }
2095
2096                         VECADDFAC(temp,vec_to_part,eff_dir,-fac);
2097                         r_fac=VecLength(temp);
2098                         if(pd->flag&PFIELD_USEMAXR && r_fac>pd->maxrad){
2099                                 falloff=0.0f;
2100                                 break;
2101                         }
2102
2103                         fac=ABS(fac);
2104
2105
2106                         if(pd->flag & PFIELD_USEMIN){
2107                                 if(fac>pd->mindist)
2108                                         fac+=1.0f-pd->mindist;
2109                                 else
2110                                         fac=1.0f;
2111                         }
2112                         else if(fac<0.001)
2113                                 fac=0.001f;
2114
2115                         if(pd->flag & PFIELD_USEMINR){
2116                                 if(r_fac>pd->minrad)
2117                                         r_fac+=1.0f-pd->minrad;
2118                                 else
2119                                         r_fac=1.0f;
2120                         }
2121                         else if(r_fac<0.001)
2122                                 r_fac=0.001f;
2123
2124                         falloff=1.0f/((float)pow((double)fac,(double)pd->f_power)*(float)pow((double)r_fac,(double)pd->f_power_r));
2125
2126                         break;
2127                 case PFIELD_FALL_CONE:
2128                         fac=Inpf(vec_to_part,eff_dir);
2129                         if(pd->flag&PFIELD_USEMAX && ABS(fac)>pd->maxdist){
2130                                 falloff=0.0f;
2131                                 break;
2132                         }
2133                         r_fac=saacos(fac/VecLength(vec_to_part))*180.0f/(float)M_PI;
2134                         if(pd->flag&PFIELD_USEMAXR && r_fac>pd->maxrad){
2135                                 falloff=0.0f;
2136                                 break;
2137                         }
2138
2139                         if(pd->flag & PFIELD_USEMIN){
2140                                 if(fac>pd->mindist)
2141                                         fac+=1.0f-pd->mindist;
2142                                 else
2143                                         fac=1.0f;
2144                         }
2145                         else if(fac<0.001)
2146                                 fac=0.001f;
2147
2148                         if(pd->flag & PFIELD_USEMINR){
2149                                 if(r_fac>pd->minrad)
2150                                         r_fac+=1.0f-pd->minrad;
2151                                 else
2152                                         r_fac=1.0f;
2153                         }
2154                         else if(r_fac<0.001)
2155                                 r_fac=0.001f;
2156
2157                         falloff=1.0f/((float)pow((double)fac,(double)pd->f_power)*(float)pow((double)r_fac,(double)pd->f_power_r));
2158
2159                         break;
2160 //              case PFIELD_FALL_INSIDE:
2161                                 //for(i=0; i<totface; i++,mface++){
2162                                 //      VECCOPY(v1,mvert[mface->v1].co);
2163                                 //      VECCOPY(v2,mvert[mface->v2].co);
2164                                 //      VECCOPY(v3,mvert[mface->v3].co);
2165
2166                                 //      if(AxialLineIntersectsTriangle(a,co1, co2, v2, v3, v1, &lambda)){
2167                                 //              if(from==PART_FROM_FACE)
2168                                 //                      (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
2169                                 //              else /* store number of intersections */
2170                                 //                      (pa+(int)(lambda*size[a])*a0mul)->loop++;
2171                                 //      }
2172                                 //      
2173                                 //      if(mface->v4){
2174                                 //              VECCOPY(v4,mvert[mface->v4].co);
2175
2176                                 //              if(AxialLineIntersectsTriangle(a,co1, co2, v4, v1, v3, &lambda)){
2177                                 //                      if(from==PART_FROM_FACE)
2178                                 //                              (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
2179                                 //                      else
2180                                 //                              (pa+(int)(lambda*size[a])*a0mul)->loop++;
2181                                 //              }
2182                                 //      }
2183                                 //}
2184
2185 //                      break;
2186         }
2187
2188         return falloff;
2189 }
2190 static void do_physical_effector(short type, float force_val, float distance, float falloff, float size, float damp,
2191                                                         float *eff_velocity, float *vec_to_part, float *velocity, float *field, int planar)
2192 {
2193         float mag_vec[3]={0,0,0};
2194         float temp[3], temp2[3];
2195         float eff_vel[3];
2196
2197         VecCopyf(eff_vel,eff_velocity);
2198         Normalize(eff_vel);
2199
2200         switch(type){
2201                 case PFIELD_WIND:
2202                         VECCOPY(mag_vec,eff_vel);
2203
2204                         VecMulf(mag_vec,force_val*falloff);
2205                         VecAddf(field,field,mag_vec);
2206                         break;
2207
2208                 case PFIELD_FORCE:
2209                         if(planar)
2210                                 Projf(mag_vec,vec_to_part,eff_vel);
2211                         else
2212                                 VecCopyf(mag_vec,vec_to_part);
2213
2214                         VecMulf(mag_vec,force_val*falloff);
2215                         VecAddf(field,field,mag_vec);
2216                         break;
2217
2218                 case PFIELD_VORTEX:
2219                         Crossf(mag_vec,eff_vel,vec_to_part);
2220                         Normalize(mag_vec);
2221
2222                         VecMulf(mag_vec,force_val*distance*falloff);
2223                         VecAddf(field,field,mag_vec);
2224
2225                         break;
2226                 case PFIELD_MAGNET:
2227                         if(planar)
2228                                 VecCopyf(temp,eff_vel);
2229                         else
2230                                 /* magnetic field of a moving charge */
2231                                 Crossf(temp,eff_vel,vec_to_part);
2232
2233                         Crossf(temp2,velocity,temp);
2234                         VecAddf(mag_vec,mag_vec,temp2);
2235
2236                         VecMulf(mag_vec,force_val*falloff);
2237                         VecAddf(field,field,mag_vec);
2238                         break;
2239                 case PFIELD_HARMONIC:
2240                         if(planar)
2241                                 Projf(mag_vec,vec_to_part,eff_vel);
2242                         else
2243                                 VecCopyf(mag_vec,vec_to_part);
2244
2245                         VecMulf(mag_vec,force_val*falloff);
2246                         VecSubf(field,field,mag_vec);
2247
2248                         VecCopyf(mag_vec,velocity);
2249                         /* 1.9 is an experimental value to get critical damping at damp=1.0 */
2250                         VecMulf(mag_vec,damp*1.9f*(float)sqrt(force_val));
2251                         VecSubf(field,field,mag_vec);
2252                         break;
2253                 case PFIELD_NUCLEAR:
2254                         /*pow here is root of cosine expression below*/
2255                         //rad=(float)pow(2.0,-1.0/power)*distance/size;
2256                         //VECCOPY(mag_vec,vec_to_part);
2257                         //Normalize(mag_vec);
2258                         //VecMulf(mag_vec,(float)cos(3.0*M_PI/2.0*(1.0-1.0/(pow(rad,power)+1.0)))/(rad+0.2f));
2259                         //VECADDFAC(field,field,mag_vec,force_val);
2260                         break;
2261         }
2262 }
2263 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)
2264 {
2265         TexResult result[4];
2266         float tex_co[3], strength, mag_vec[3];
2267         int i;
2268
2269         if(tex==0) return;
2270
2271         for(i=0; i<4; i++)
2272                 result[i].nor=0;
2273
2274         strength= force_val*falloff;///(float)pow((double)distance,(double)power);
2275
2276         VECCOPY(tex_co,pa_co);
2277
2278         if(is_2d){
2279                 float fac=-Inpf(tex_co,obmat[2]);
2280                 VECADDFAC(tex_co,tex_co,obmat[2],fac);
2281         }
2282
2283         if(object){
2284                 VecSubf(tex_co,tex_co,obmat[3]);
2285                 Mat4Mul3Vecfl(obmat,tex_co);
2286         }
2287
2288         multitex_ext(tex, tex_co, NULL,NULL, 1, result);
2289
2290         if(mode==PFIELD_TEX_RGB){
2291                 mag_vec[0]= (0.5f-result->tr)*strength;
2292                 mag_vec[1]= (0.5f-result->tg)*strength;
2293                 mag_vec[2]= (0.5f-result->tb)*strength;
2294         }
2295         else{
2296                 strength/=nabla;
2297
2298                 tex_co[0]+= nabla;
2299                 multitex_ext(tex, tex_co, NULL,NULL, 1, result+1);
2300
2301                 tex_co[0]-= nabla;
2302                 tex_co[1]+= nabla;
2303                 multitex_ext(tex, tex_co, NULL,NULL, 1, result+2);
2304
2305                 tex_co[1]-= nabla;
2306                 tex_co[2]+= nabla;
2307                 multitex_ext(tex, tex_co, NULL,NULL, 1, result+3);
2308
2309                 if(mode==PFIELD_TEX_GRAD){
2310                         mag_vec[0]= (result[0].tin-result[1].tin)*strength;
2311                         mag_vec[1]= (result[0].tin-result[2].tin)*strength;
2312                         mag_vec[2]= (result[0].tin-result[3].tin)*strength;
2313                 }
2314                 else{ /*PFIELD_TEX_CURL*/
2315                         float dbdy,dgdz,drdz,dbdx,dgdx,drdy;
2316
2317                         dbdy= result[2].tb-result[0].tb;
2318                         dgdz= result[3].tg-result[0].tg;
2319                         drdz= result[3].tr-result[0].tr;
2320                         dbdx= result[1].tb-result[0].tb;
2321                         dgdx= result[1].tg-result[0].tg;
2322                         drdy= result[2].tr-result[0].tr;
2323
2324                         mag_vec[0]=(dbdy-dgdz)*strength;
2325                         mag_vec[1]=(drdz-dbdx)*strength;
2326                         mag_vec[2]=(dgdx-drdy)*strength;
2327                 }
2328         }
2329
2330         if(is_2d){
2331                 float fac=-Inpf(mag_vec,obmat[2]);
2332                 VECADDFAC(mag_vec,mag_vec,obmat[2],fac);
2333         }
2334
2335         VecAddf(field,field,mag_vec);
2336 }
2337 static void add_to_effectors(ListBase *lb, Object *ob, Object *obsrc, ParticleSystem *psys)
2338 {
2339         ParticleEffectorCache *ec;
2340         PartDeflect *pd= ob->pd;
2341         short type=0,i;
2342
2343         if(pd && ob != obsrc){
2344                 if(pd->forcefield == PFIELD_GUIDE) {
2345                         if(ob->type==OB_CURVE) {
2346                                 Curve *cu= ob->data;
2347                                 if(cu->flag & CU_PATH) {
2348                                         if(cu->path==NULL || cu->path->data==NULL)
2349                                                 makeDispListCurveTypes(ob, 0);
2350                                         if(cu->path && cu->path->data) {
2351                                                 type |= PSYS_EC_EFFECTOR;
2352                                         }
2353                                 }
2354                         }
2355                 }
2356                 else if(pd->forcefield)
2357                         type |= PSYS_EC_EFFECTOR;
2358         }
2359         
2360         if(pd && pd->deflect)
2361                 type |= PSYS_EC_DEFLECT;
2362
2363         if(type){
2364                 ec= MEM_callocN(sizeof(ParticleEffectorCache), "effector cache");
2365                 ec->ob= ob;
2366                 ec->type=type;
2367                 ec->distances=0;
2368                 ec->locations=0;
2369                 BLI_addtail(lb, ec);
2370         }
2371
2372         type=0;
2373
2374         /* add particles as different effectors */
2375         if(ob->particlesystem.first){
2376                 ParticleSystem *epsys=ob->particlesystem.first;
2377                 ParticleSettings *epart=0;
2378                 Object *tob;
2379
2380                 for(i=0; epsys; epsys=epsys->next,i++){
2381                         type=0;
2382                         if(epsys!=psys){
2383                                 epart=epsys->part;
2384
2385                                 if(epsys->part->pd && epsys->part->pd->forcefield)
2386                                         type=PSYS_EC_PARTICLE;
2387
2388                                 if(epart->type==PART_REACTOR) {
2389                                         tob=epsys->target_ob;
2390                                         if(tob==0)
2391                                                 tob=ob;
2392                                         if(BLI_findlink(&tob->particlesystem,epsys->target_psys-1)==psys)
2393                                                 type|=PSYS_EC_REACTOR;
2394                                 }
2395
2396                                 if(type){
2397                                         ec= MEM_callocN(sizeof(ParticleEffectorCache), "effector cache");
2398                                         ec->ob= ob;
2399                                         ec->type=type;
2400                                         ec->psys_nbr=i;
2401                                         BLI_addtail(lb, ec);
2402                                 }
2403                         }
2404                 }
2405                                 
2406         }
2407 }
2408 void psys_init_effectors(Object *obsrc, Group *group, ParticleSystem *psys)
2409 {
2410         ListBase *listb=&psys->effectors;
2411         Base *base;
2412         unsigned int layer= obsrc->lay;
2413
2414         listb->first=listb->last=0;
2415         
2416         if(group) {
2417                 GroupObject *go;
2418                 
2419                 for(go= group->gobject.first; go; go= go->next) {
2420                         if( (go->ob->lay & layer) && (go->ob->pd || go->ob->particlesystem.first)) {
2421                                 add_to_effectors(listb, go->ob, obsrc, psys);
2422                         }
2423                 }
2424         }
2425         else {
2426                 for(base = G.scene->base.first; base; base= base->next) {
2427                         if( (base->lay & layer) && (base->object->pd || base->object->particlesystem.first)) {
2428                                 add_to_effectors(listb, base->object, obsrc, psys);
2429                         }
2430                 }
2431         }
2432 }
2433
2434 void psys_end_effectors(ParticleSystem *psys)
2435 {
2436         ListBase *lb=&psys->effectors;
2437         if(lb->first) {
2438                 ParticleEffectorCache *ec;
2439                 for(ec= lb->first; ec; ec= ec->next){
2440                         if(ec->distances)
2441                                 MEM_freeN(ec->distances);
2442
2443                         if(ec->locations)
2444                                 MEM_freeN(ec->locations);
2445
2446                         if(ec->face_minmax)
2447                                 MEM_freeN(ec->face_minmax);
2448
2449                         if(ec->vert_cos)
2450                                 MEM_freeN(ec->vert_cos);
2451
2452                         if(ec->tree)
2453                                 BLI_kdtree_free(ec->tree);
2454                 }
2455
2456                 BLI_freelistN(lb);
2457         }
2458 }
2459
2460 static void precalc_effectors(Object *ob, ParticleSystem *psys, ParticleSystemModifierData *psmd)
2461 {
2462         ListBase *lb=&psys->effectors;
2463         ParticleEffectorCache *ec;
2464         ParticleSettings *part=psys->part;
2465         ParticleData *pa;
2466         float vec2[3],loc[3],*co=0;
2467         int p,totpart,totvert;
2468         
2469         for(ec= lb->first; ec; ec= ec->next) {
2470                 PartDeflect *pd= ec->ob->pd;
2471                 
2472                 if(ec->type==PSYS_EC_EFFECTOR && pd->forcefield==PFIELD_GUIDE && ec->ob->type==OB_CURVE 
2473                         && part->phystype!=PART_PHYS_BOIDS) {
2474                         float vec[4];
2475
2476                         where_on_path(ec->ob, 0.0, vec, vec2);
2477
2478                         Mat4MulVecfl(ec->ob->obmat,vec);
2479                         Mat4Mul3Vecfl(ec->ob->obmat,vec2);
2480
2481                         QUATCOPY(ec->firstloc,vec);
2482                         VECCOPY(ec->firstdir,vec2);
2483
2484                         totpart=psys->totpart;
2485
2486                         if(totpart){
2487                                 ec->distances=MEM_callocN(totpart*sizeof(float),"particle distances");
2488                                 ec->locations=MEM_callocN(totpart*3*sizeof(float),"particle locations");
2489
2490                                 for(p=0,pa=psys->particles; p<totpart; p++, pa++){
2491                                         psys_particle_on_emitter(ob,psmd,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,loc,0,0,0,0,0);
2492                                         Mat4MulVecfl(ob->obmat,loc);
2493                                         ec->distances[p]=VecLenf(loc,vec);
2494                                         VECSUB(loc,loc,vec);
2495                                         VECCOPY(ec->locations+3*p,loc);
2496                                 }
2497                         }
2498                 }
2499                 else if(ec->type==PSYS_EC_DEFLECT){
2500                         DerivedMesh *dm;
2501                         MFace *mface=0;
2502                         MVert *mvert=0;
2503                         int i, totface;
2504                         float v1[3],v2[3],v3[3],v4[4], *min, *max;
2505
2506                         if(ob==ec->ob)
2507                                 dm=psmd->dm;
2508                         else{
2509                                 psys_disable_all(ec->ob);
2510
2511                                 dm=mesh_get_derived_final(ec->ob,0);
2512                                 
2513                                 psys_enable_all(ec->ob);
2514                         }
2515
2516                         if(dm){
2517                                 totvert=dm->getNumVerts(dm);
2518                                 totface=dm->getNumFaces(dm);
2519                                 mface=dm->getFaceDataArray(dm,CD_MFACE);
2520                                 mvert=dm->getVertDataArray(dm,CD_MVERT);
2521
2522                                 /* Decide which is faster to calculate by the amount of*/
2523                                 /* matrice multiplications needed to convert spaces. */
2524                                 /* With size deflect we have to convert allways because */
2525                                 /* the object can be scaled nonuniformly (sphere->ellipsoid). */
2526                                 if(totvert<2*psys->totpart || part->flag & PART_SIZE_DEFL){
2527                                         co=ec->vert_cos=MEM_callocN(sizeof(float)*3*totvert,"Particle deflection vert cos");
2528                                         /* convert vert coordinates to global (particle) coordinates */
2529                                         for(i=0; i<totvert; i++, co+=3){
2530                                                 VECCOPY(co,mvert[i].co);
2531                                                 Mat4MulVecfl(ec->ob->obmat,co);
2532                                         }
2533                                         co=ec->vert_cos;
2534                                 }
2535                                 else
2536                                         ec->vert_cos=0;
2537
2538                                 INIT_MINMAX(ec->ob_minmax,ec->ob_minmax+3);
2539
2540                                 min=ec->face_minmax=MEM_callocN(sizeof(float)*6*totface,"Particle deflection face minmax");
2541                                 max=min+3;
2542
2543                                 for(i=0; i<totface; i++,mface++,min+=6,max+=6){
2544                                         if(co){
2545                                                 VECCOPY(v1,co+3*mface->v1);
2546                                                 VECCOPY(v2,co+3*mface->v2);
2547                                                 VECCOPY(v3,co+3*mface->v3);
2548                                         }
2549                                         else{
2550                                                 VECCOPY(v1,mvert[mface->v1].co);
2551                                                 VECCOPY(v2,mvert[mface->v2].co);
2552                                                 VECCOPY(v3,mvert[mface->v3].co);
2553                                         }
2554                                         INIT_MINMAX(min,max);
2555                                         DO_MINMAX(v1,min,max);
2556                                         DO_MINMAX(v2,min,max);
2557                                         DO_MINMAX(v3,min,max);
2558
2559                                         if(mface->v4){
2560                                                 if(co){
2561                                                         VECCOPY(v4,co+3*mface->v4);
2562                                                 }
2563                                                 else{
2564                                                         VECCOPY(v4,mvert[mface->v4].co);
2565                                                 }
2566                                                 DO_MINMAX(v4,min,max);
2567                                         }
2568
2569                                         DO_MINMAX(min,ec->ob_minmax,ec->ob_minmax+3);
2570                                         DO_MINMAX(max,ec->ob_minmax,ec->ob_minmax+3);
2571                                 }
2572                         }
2573                         else
2574                                 ec->face_minmax=0;
2575                 }
2576                 else if(ec->type==PSYS_EC_PARTICLE){
2577                         if(psys->part->phystype==PART_PHYS_BOIDS){
2578                                 Object *eob = ec->ob;
2579                                 ParticleSystem *epsys;
2580                                 ParticleSettings *epart;
2581                                 ParticleData *epa;
2582                                 ParticleKey state;
2583                                 PartDeflect *pd;
2584                                 int totepart, p;
2585                                 epsys= BLI_findlink(&eob->particlesystem,ec->psys_nbr);
2586                                 epart= epsys->part;
2587                                 pd= epart->pd;
2588                                 totepart= epsys->totpart;
2589                                 if(pd->forcefield==PFIELD_FORCE && totepart){
2590                                         KDTree *tree;
2591
2592                                         tree=BLI_kdtree_new(totepart);
2593                                         ec->tree=tree;
2594
2595                                         for(p=0, epa=epsys->particles; p<totepart; p++,epa++)
2596                                                 if(epa->alive==PARS_ALIVE && psys_get_particle_state(eob,epsys,p,&state,0))
2597                                                         BLI_kdtree_insert(tree, p, state.co, NULL);
2598
2599                                         BLI_kdtree_balance(tree);
2600                                 }
2601                         }
2602                 }
2603         }
2604 }
2605
2606
2607 /* calculate forces that all effectors apply to a particle*/
2608 static void do_effectors(int pa_no, ParticleData *pa, ParticleKey *state, Object *ob, ParticleSystem *psys, float *force_field, float *vel,float framestep, float cfra)
2609 {
2610         Object *eob;
2611         ParticleSystem *epsys;
2612         ParticleSettings *epart;
2613         ParticleData *epa;
2614         ParticleKey estate;
2615         PartDeflect *pd;
2616         ListBase *lb=&psys->effectors;
2617         ParticleEffectorCache *ec;
2618         float distance, vec_to_part[3];
2619         float falloff;
2620         int p;
2621
2622         /* check all effector objects for interaction */
2623         if(lb->first){
2624                 for(ec = lb->first; ec; ec= ec->next){
2625                         eob= ec->ob;
2626                         if(ec->type & PSYS_EC_EFFECTOR){
2627                                 pd=eob->pd;
2628                                 if(psys->part->type!=PART_HAIR && psys->part->integrator)
2629                                         where_is_object_time(eob,cfra);
2630                                 /* Get IPO force strength and fall off values here */
2631                                 //if (has_ipo_code(eob->ipo, OB_PD_FSTR))
2632                                 //      force_val = IPO_GetFloatValue(eob->ipo, OB_PD_FSTR, cfra);
2633                                 //else 
2634                                 //      force_val = pd->f_strength;
2635                                 
2636                                 //if (has_ipo_code(eob->ipo, OB_PD_FFALL)) 
2637                                 //      ffall_val = IPO_GetFloatValue(eob->ipo, OB_PD_FFALL, cfra);
2638                                 //else 
2639                                 //      ffall_val = pd->f_power;
2640
2641                                 //if (has_ipo_code(eob->ipo, OB_PD_FMAXD)) 
2642                                 //      maxdist = IPO_GetFloatValue(eob->ipo, OB_PD_FMAXD, cfra);
2643                                 //else 
2644                                 //      maxdist = pd->maxdist;
2645
2646                                 /* use center of object for distance calculus */
2647                                 //obloc= eob->obmat[3];
2648                                 VecSubf(vec_to_part, state->co, eob->obmat[3]);
2649                                 distance = VecLength(vec_to_part);
2650
2651                                 falloff=effector_falloff(pd,eob->obmat[2],vec_to_part);
2652
2653                                 if(falloff<=0.0f)
2654                                         ;       /* don't do anything */
2655                                 else if(pd->forcefield==PFIELD_TEXTURE)
2656                                         do_texture_effector(pd->tex, pd->tex_mode, pd->flag&PFIELD_TEX_2D, pd->tex_nabla,
2657                                                                                 pd->flag & PFIELD_TEX_OBJECT, state->co, eob->obmat,
2658                                                                                 pd->f_strength, falloff, force_field);
2659                                 else
2660                                         do_physical_effector(pd->forcefield,pd->f_strength,distance,
2661                                                                                 falloff,pd->f_dist,pd->f_damp,eob->obmat[2],vec_to_part,
2662                                                                                 pa->state.vel,force_field,pd->flag&PFIELD_PLANAR);
2663                         }
2664                         if(ec->type & PSYS_EC_PARTICLE){
2665                                 int totepart;
2666                                 epsys= BLI_findlink(&eob->particlesystem,ec->psys_nbr);
2667                                 epart= epsys->part;
2668                                 pd= epart->pd;
2669                                 totepart= epsys->totpart;
2670
2671                                 if(pd->forcefield==PFIELD_HARMONIC){
2672                                         /* every particle is mapped to only one harmonic effector particle */
2673                                         p= pa_no%epsys->totpart;
2674                                         totepart= p+1;
2675                                 }
2676                                 else{
2677                                         p=0;
2678                                 }
2679
2680                                 epsys->lattice=psys_get_lattice(ob,psys);
2681
2682                                 for(; p<totepart; p++){
2683                                         epa = epsys->particles + p;
2684                                         estate.time=-1.0;
2685                                         if(psys_get_particle_state(eob,epsys,p,&estate,0)){
2686                                                 VECSUB(vec_to_part, state->co, estate.co);
2687                                                 distance = VecLength(vec_to_part);
2688
2689                                                 //if(pd->forcefield==PFIELD_HARMONIC){
2690                                                 //      //if(cfra < epa->time + radius){ /* radius is fade-in in ui */
2691                                                 //      //      eforce*=(cfra-epa->time)/radius;
2692                                                 //      //}
2693                                                 //}
2694                                                 //else{
2695                                                 //      /* Limit minimum distance to effector particle so that */
2696                                                 //      /* the force is not too big */
2697                                                 //      if (distance < 0.001) distance = 0.001f;
2698                                                 //}
2699                                                 
2700                                                 falloff=effector_falloff(pd,estate.vel,vec_to_part);
2701
2702                                                 if(falloff<=0.0f)
2703                                                         ;       /* don't do anything */
2704                                                 else
2705                                                         do_physical_effector(pd->forcefield,pd->f_strength,distance,
2706                                                         falloff,epart->size,pd->f_damp,estate.vel,vec_to_part,
2707                                                         state->vel,force_field,0);
2708                                         }
2709                                         else if(pd->forcefield==PFIELD_HARMONIC && cfra-framestep <= epa->dietime && cfra>epa->dietime){
2710                                                 /* first step after key release */
2711                                                 psys_get_particle_state(eob,epsys,p,&estate,1);
2712                                                 VECADD(vel,vel,estate.vel);
2713                                                 /* TODO: add rotation handling here too */
2714                                         }
2715                                 }
2716
2717                                 if(epsys->lattice){
2718                                         end_latt_deform();
2719                                         epsys->lattice=0;
2720                                 }
2721                         }
2722                 }
2723         }
2724 }
2725
2726 /************************************************/
2727 /*                      Newtonian physics                                       */
2728 /************************************************/
2729 /* gathers all forces that effect particles and calculates a new state for the particle */
2730 static void apply_particle_forces(int pa_no, ParticleData *pa, Object *ob, ParticleSystem *psys, ParticleSettings *part, float timestep, float dfra, float cfra, ParticleKey *state)
2731 {
2732         ParticleKey states[5], tkey;
2733         float force[3],tvel[3],dx[4][3],dv[4][3];
2734         float dtime=dfra*timestep, time, pa_mass=part->mass, fac, fra=psys->cfra;
2735         int i, steps=1;
2736         
2737         /* maintain angular velocity */
2738         VECCOPY(state->ave,pa->state.ave);
2739
2740         if(part->flag & PART_SIZEMASS)
2741                 pa_mass*=pa->size;
2742
2743         switch(part->integrator){
2744                 case PART_INT_EULER:
2745                         steps=1;
2746                         break;
2747                 case PART_INT_MIDPOINT:
2748                         steps=2;
2749                         break;
2750                 case PART_INT_RK4:
2751                         steps=4;
2752                         break;
2753         }
2754
2755         copy_particle_key(states,&pa->state,1);
2756
2757         for(i=0; i<steps; i++){
2758                 force[0]=force[1]=force[2]=0.0;
2759                 tvel[0]=tvel[1]=tvel[2]=0.0;
2760                 /* add effectors */
2761                 do_effectors(pa_no,pa,states+i,ob,psys,force,tvel,dfra,fra);
2762
2763                 /* calculate air-particle interaction */
2764                 if(part->dragfac!=0.0f){
2765                         fac=-part->dragfac*pa->size*pa->size*VecLength(states[i].vel);
2766                         VECADDFAC(force,force,states[i].vel,fac);
2767                 }
2768
2769                 /* brownian force */
2770                 if(part->brownfac!=0.0){
2771                         force[0]+=(BLI_frand()-0.5f)*part->brownfac;
2772                         force[1]+=(BLI_frand()-0.5f)*part->brownfac;
2773                         force[2]+=(BLI_frand()-0.5f)*part->brownfac;
2774                 }
2775
2776                 /* force to acceleration*/
2777                 VecMulf(force,1.0f/pa_mass);
2778
2779                 /* add global acceleration (gravitation) */
2780                 VECADD(force,force,part->acc);
2781
2782                 //VecMulf(force,dtime);
2783                 
2784                 /* calculate next state */
2785                 VECADD(states[i].vel,states[i].vel,tvel);
2786
2787                 //VecMulf(force,0.5f*dt);
2788                 switch(part->integrator){
2789                         case PART_INT_EULER:
2790                                 VECADDFAC(state->co,states->co,states->vel,dtime);
2791                                 VECADDFAC(state->vel,states->vel,force,dtime);
2792                                 break;
2793                         case PART_INT_MIDPOINT:
2794                                 if(i==0){
2795                                         VECADDFAC(states[1].co,states->co,states->vel,dtime*0.5f);
2796                                         VECADDFAC(states[1].vel,states->vel,force,dtime*0.5f);
2797                                         fra=psys->cfra+0.5f*dfra;
2798                                 }
2799                                 else{
2800                                         VECADDFAC(state->co,states->co,states[1].vel,dtime);
2801                                         VECADDFAC(state->vel,states->vel,force,dtime);
2802                                 }
2803                                 break;
2804                         case PART_INT_RK4:
2805                                 switch(i){
2806                                         case 0:
2807                                                 VECCOPY(dx[0],states->vel);
2808                                                 VecMulf(dx[0],dtime);
2809                                                 VECCOPY(dv[0],force);
2810                                                 VecMulf(dv[0],dtime);
2811
2812                                                 VECADDFAC(states[1].co,states->co,dx[0],0.5f);
2813                                                 VECADDFAC(states[1].vel,states->vel,dv[0],0.5f);
2814                                                 fra=psys->cfra+0.5f*dfra;
2815                                                 break;
2816                                         case 1:
2817                                                 VECADDFAC(dx[1],states->vel,dv[0],0.5f);
2818                                                 VecMulf(dx[1],dtime);
2819                                                 VECCOPY(dv[1],force);
2820                                                 VecMulf(dv[1],dtime);
2821
2822                                                 VECADDFAC(states[2].co,states->co,dx[1],0.5f);
2823                                                 VECADDFAC(states[2].vel,states->vel,dv[1],0.5f);
2824                                                 break;
2825                                         case 2:
2826                                                 VECADDFAC(dx[2],states->vel,dv[1],0.5f);
2827                                                 VecMulf(dx[2],dtime);
2828                                                 VECCOPY(dv[2],force);
2829                                                 VecMulf(dv[2],dtime);
2830
2831                                                 VECADD(states[3].co,states->co,dx[2]);
2832                                                 VECADD(states[3].vel,states->vel,dv[2]);
2833                                                 fra=cfra;
2834                                                 break;
2835                                         case 3:
2836                                                 VECADD(dx[3],states->vel,dv[2]);
2837                                                 VecMulf(dx[3],dtime);