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