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