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