Fixed all gcc 4 warnings in blenkernel. Found 2 potentially harmful
[blender.git] / source / blender / blenkernel / intern / effect.c
1 /*  effect.c
2  * 
3  * 
4  * $Id$
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) 2001-2002 by NaN Holding BV.
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 <math.h>
36 #include <stdlib.h>
37
38 #include "MEM_guardedalloc.h"
39
40 #include "DNA_curve_types.h"
41 #include "DNA_effect_types.h"
42 #include "DNA_group_types.h"
43 #include "DNA_ipo_types.h"
44 #include "DNA_key_types.h"
45 #include "DNA_lattice_types.h"
46 #include "DNA_listBase.h"
47 #include "DNA_mesh_types.h"
48 #include "DNA_meshdata_types.h"
49 #include "DNA_material_types.h"
50 #include "DNA_object_types.h"
51 #include "DNA_object_force.h"
52 #include "DNA_texture_types.h"
53 #include "DNA_scene_types.h"
54
55 #include "BLI_arithb.h"
56 #include "BLI_blenlib.h"
57 #include "BLI_jitter.h"
58 #include "BLI_rand.h"
59
60 #include "BKE_action.h"
61 #include "BKE_anim.h"           /* needed for where_on_path */
62 #include "BKE_armature.h"
63 #include "BKE_bad_level_calls.h"
64 #include "BKE_blender.h"
65 #include "BKE_constraint.h"
66 #include "BKE_deform.h"
67 #include "BKE_depsgraph.h"
68 #include "BKE_displist.h"
69 #include "BKE_DerivedMesh.h"
70 #include "BKE_effect.h"
71 #include "BKE_global.h"
72 #include "BKE_group.h"
73 #include "BKE_ipo.h"
74 #include "BKE_key.h"
75 #include "BKE_lattice.h"
76 #include "BKE_mesh.h"
77 #include "BKE_material.h"
78 #include "BKE_main.h"
79 #include "BKE_object.h"
80 #include "BKE_scene.h"
81 #include "BKE_screen.h"
82 #include "BKE_utildefines.h"
83
84 #include "PIL_time.h"
85 #include "RE_render_ext.h"
86
87 /* fluid sim particle import */
88 #ifndef DISABLE_ELBEEM
89 #include "DNA_object_fluidsim.h"
90 #include "LBM_fluidsim.h"
91 #include "elbeem.h"
92 #include <zlib.h>
93 #include <string.h>
94 #endif // DISABLE_ELBEEM
95
96 /* temporal struct, used for reading return of mesh_get_mapped_verts_nors() */
97 typedef struct VeNoCo {
98         float co[3], no[3];
99 } VeNoCo;
100
101 Effect *add_effect(int type)
102 {
103         Effect *eff=0;
104         PartEff *paf;
105         int a;
106         
107         switch(type) {
108         case EFF_PARTICLE:
109                 paf= MEM_callocN(sizeof(PartEff), "neweff");
110                 eff= (Effect *)paf;
111                 
112                 paf->sta= 1.0;
113                 paf->end= 100.0;
114                 paf->lifetime= 50.0;
115                 for(a=0; a<PAF_MAXMULT; a++) {
116                         paf->life[a]= 50.0;
117                         paf->child[a]= 4;
118                         paf->mat[a]= 1;
119                 }
120                 
121                 paf->totpart= 1000;
122                 paf->totkey= 8;
123                 paf->staticstep= 5;
124                 paf->defvec[2]= 1.0f;
125                 paf->nabla= 0.05f;
126                 paf->disp = 100;
127                 paf->speedtex = 8;
128                 paf->omat = 1;
129                 paf->flag= PAF_FACE;
130
131                 break;
132         }
133         
134         eff->type= eff->buttype= type;
135         eff->flag |= SELECT;
136         
137         return eff;
138 }
139
140 void free_effect(Effect *eff)
141 {
142         PartEff *paf;
143         
144         if(eff->type==EFF_PARTICLE) {
145                 paf= (PartEff *)eff;
146                 if(paf->keys) MEM_freeN(paf->keys);
147         }
148         MEM_freeN(eff); 
149 }
150
151
152 void free_effects(ListBase *lb)
153 {
154         Effect *eff;
155         
156         eff= lb->first;
157         while(eff) {
158                 BLI_remlink(lb, eff);
159                 free_effect(eff);
160                 eff= lb->first;
161         }
162 }
163
164 Effect *copy_effect(Effect *eff) 
165 {
166         Effect *effn;
167
168         effn= MEM_dupallocN(eff);
169         if(effn->type==EFF_PARTICLE) ((PartEff *)effn)->keys= 0;
170
171         return effn;    
172 }
173
174 void copy_act_effect(Object *ob)
175 {
176         /* return a copy of the active effect */
177         Effect *effn, *eff;
178         
179         eff= ob->effect.first;
180         while(eff) {
181                 if(eff->flag & SELECT) {
182                         
183                         effn= copy_effect(eff);
184                         BLI_addtail(&ob->effect, effn);
185                         
186                         eff->flag &= ~SELECT;
187                         return;
188                         
189                 }
190                 eff= eff->next;
191         }
192         
193         /* when it comes here: add new effect */
194         eff= add_effect(EFF_PARTICLE);
195         BLI_addtail(&ob->effect, eff);
196                         
197 }
198
199 void copy_effects(ListBase *lbn, ListBase *lb)
200 {
201         Effect *eff, *effn;
202
203         lbn->first= lbn->last= 0;
204
205         eff= lb->first;
206         while(eff) {
207                 effn= copy_effect(eff);
208                 BLI_addtail(lbn, effn);
209                 
210                 eff= eff->next;
211         }
212         
213 }
214
215 void deselectall_eff(Object *ob)
216 {
217         Effect *eff= ob->effect.first;
218         
219         while(eff) {
220                 eff->flag &= ~SELECT;
221                 eff= eff->next;
222         }
223 }
224
225 /* ***************** PARTICLES ***************** */
226
227 static Particle *new_particle(PartEff *paf)
228 {
229         static Particle *pa;
230         static int cur;
231
232         /* we agree: when paf->keys==0: alloc */        
233         if(paf->keys==NULL) {
234                 pa= paf->keys= MEM_callocN( paf->totkey*paf->totpart*sizeof(Particle), "particlekeys" );
235                 cur= 0;
236         }
237         else {
238                 if(cur && cur<paf->totpart) pa+=paf->totkey;
239                 cur++;
240         }
241         return pa;
242 }
243
244 PartEff *give_parteff(Object *ob)
245 {
246         PartEff *paf;
247         
248         paf= ob->effect.first;
249         while(paf) {
250                 if(paf->type==EFF_PARTICLE) return paf;
251                 paf= paf->next;
252         }
253         return 0;
254 }
255
256 void where_is_particle(PartEff *paf, Particle *pa, float ctime, float *vec)
257 {
258         Particle *p[4];
259         float dt, t[4];
260         int a;
261         
262         if(paf->totkey==1 || ctime < pa->time) {
263                 VECCOPY(vec, pa->co);
264                 return;
265         }
266         
267         /* first find the first particlekey */
268         a= (int)((paf->totkey-1)*(ctime-pa->time)/pa->lifetime);
269         if(a>=paf->totkey) a= paf->totkey-1;
270         else if(a<0) a= 0;
271         
272         pa+= a;
273         
274         if(a>0) p[0]= pa-1; else p[0]= pa;
275         p[1]= pa;
276         
277         if(a+1<paf->totkey) p[2]= pa+1; else p[2]= pa;
278         if(a+2<paf->totkey) p[3]= pa+2; else p[3]= p[2];
279         
280         if(p[1]==p[2] || p[2]->time == p[1]->time) dt= 0.0;
281         else dt= (ctime-p[1]->time)/(p[2]->time - p[1]->time);
282
283         if(paf->flag & PAF_BSPLINE) set_four_ipo(dt, t, KEY_BSPLINE);
284         else set_four_ipo(dt, t, KEY_CARDINAL);
285
286         vec[0]= t[0]*p[0]->co[0] + t[1]*p[1]->co[0] + t[2]*p[2]->co[0] + t[3]*p[3]->co[0];
287         vec[1]= t[0]*p[0]->co[1] + t[1]*p[1]->co[1] + t[2]*p[2]->co[1] + t[3]*p[3]->co[1];
288         vec[2]= t[0]*p[0]->co[2] + t[1]*p[1]->co[2] + t[2]*p[2]->co[2] + t[3]*p[3]->co[2];
289
290 }
291
292 static void particle_tex(MTex *mtex, PartEff *paf, float *co, float *no)
293 {                               
294         float tin, tr, tg, tb, ta;
295         float old;
296         
297         externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
298
299         if(paf->texmap==PAF_TEXINT) {
300                 tin*= paf->texfac;
301                 no[0]+= tin*paf->defvec[0];
302                 no[1]+= tin*paf->defvec[1];
303                 no[2]+= tin*paf->defvec[2];
304         }
305         else if(paf->texmap==PAF_TEXRGB) {
306                 no[0]+= (tr-0.5f)*paf->texfac;
307                 no[1]+= (tg-0.5f)*paf->texfac;
308                 no[2]+= (tb-0.5f)*paf->texfac;
309         }
310         else {  /* PAF_TEXGRAD */
311                 
312                 old= tin;
313                 co[0]+= paf->nabla;
314                 externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
315                 no[0]+= (old-tin)*paf->texfac;
316                 
317                 co[0]-= paf->nabla;
318                 co[1]+= paf->nabla;
319                 externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
320                 no[1]+= (old-tin)*paf->texfac;
321                 
322                 co[1]-= paf->nabla;
323                 co[2]+= paf->nabla;
324                 externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
325                 no[2]+= (old-tin)*paf->texfac;
326                 
327         }
328 }
329
330 /* -------------------------- Effectors ------------------ */
331
332 static void add_to_effectorcache(ListBase *lb, Object *ob, Object *obsrc)
333 {
334         pEffectorCache *ec;
335         PartDeflect *pd= ob->pd;
336                         
337         if(pd->forcefield == PFIELD_GUIDE) {
338                 if(ob->type==OB_CURVE && obsrc->type==OB_MESH) {        /* guides only do mesh particles */
339                         Curve *cu= ob->data;
340                         if(cu->flag & CU_PATH) {
341                                 if(cu->path==NULL || cu->path->data==NULL)
342                                         makeDispListCurveTypes(ob, 0);
343                                 if(cu->path && cu->path->data) {
344                                         ec= MEM_callocN(sizeof(pEffectorCache), "effector cache");
345                                         ec->ob= ob;
346                                         BLI_addtail(lb, ec);
347                                 }
348                         }
349                 }
350         }
351         else if(pd->forcefield) {
352                 ec= MEM_callocN(sizeof(pEffectorCache), "effector cache");
353                 ec->ob= ob;
354                 BLI_addtail(lb, ec);
355         }
356 }
357
358 /* returns ListBase handle with objects taking part in the effecting */
359 ListBase *pdInitEffectors(Object *obsrc, Group *group)
360 {
361         static ListBase listb={NULL, NULL};
362         pEffectorCache *ec;
363         Base *base;
364         unsigned int layer= obsrc->lay;
365         
366         if(group) {
367                 GroupObject *go;
368                 
369                 for(go= group->gobject.first; go; go= go->next) {
370                         if( (go->ob->lay & layer) && go->ob->pd && go->ob!=obsrc) {
371                                 add_to_effectorcache(&listb, go->ob, obsrc);
372                         }
373                 }
374         }
375         else {
376                 for(base = G.scene->base.first; base; base= base->next) {
377                         if( (base->lay & layer) && base->object->pd && base->object!=obsrc) {
378                                 add_to_effectorcache(&listb, base->object, obsrc);
379                         }
380                 }
381         }
382         
383         /* make a full copy */
384         for(ec= listb.first; ec; ec= ec->next) {
385                 ec->obcopy= *(ec->ob);
386         }
387
388         if(listb.first)
389                 return &listb;
390         
391         return NULL;
392 }
393
394 void pdEndEffectors(ListBase *lb)
395 {
396         if(lb) {
397                 pEffectorCache *ec;
398                 /* restore full copy */
399                 for(ec= lb->first; ec; ec= ec->next)
400                         *(ec->ob)= ec->obcopy;
401
402                 BLI_freelistN(lb);
403         }
404 }
405
406 /* local for this c file, only for guides now */
407 static void precalc_effectors(Object *ob, PartEff *paf, Particle *pa, ListBase *lb)
408 {
409         pEffectorCache *ec;
410         
411         for(ec= lb->first; ec; ec= ec->next) {
412                 PartDeflect *pd= ec->ob->pd;
413                 
414                 ec->oldspeed[0]= ec->oldspeed[1]= ec->oldspeed[2]= 0.0f;
415                 
416                 if(pd->forcefield==PFIELD_GUIDE && ec->ob->type==OB_CURVE) {
417                         float vec[4], dir[3];
418                         
419                         if(!(paf->flag & PAF_STATIC))
420                                 where_is_object_time(ec->ob, pa->time);
421                         
422                         /* scale corrects speed vector to curve size */
423                         if(paf->totkey>1) ec->scale= (paf->totkey-1)/pa->lifetime;
424                         else ec->scale= 1.0f;
425                         
426                         /* time_scale is for random life */
427                         if(pa->lifetime>paf->lifetime)
428                                 ec->time_scale= paf->lifetime/pa->lifetime;
429                         else
430                                 ec->time_scale= pa->lifetime/paf->lifetime;
431
432                         /* distance of first path point to particle origin */
433                         where_on_path(ec->ob, 0.0f, vec, dir);
434                         VECCOPY(ec->oldloc, vec);       /* store local coord for differences */
435                         Mat4MulVecfl(ec->ob->obmat, vec);
436                         
437                         /* for static we need to move to global space */
438                         if(paf->flag & PAF_STATIC) {
439                                 VECCOPY(dir, pa->co);
440                                 Mat4MulVecfl(ob->obmat, dir);
441                                 ec->guide_dist= VecLenf(vec, dir);
442                         }
443                         else 
444                                 ec->guide_dist= VecLenf(vec, pa->co);
445                 }
446         }
447 }
448
449
450 /*  -------- pdDoEffectors() --------
451     generic force/speed system, now used for particles and softbodies
452         lb                      = listbase with objects that take part in effecting
453         opco            = global coord, as input
454     force               = force accumulator
455     speed               = actual current speed which can be altered
456         cur_time        = "external" time in frames, is constant for static particles
457         loc_time        = "local" time in frames, range <0-1> for the lifetime of particle
458     par_layer   = layer the caller is in
459         flags           = only used for softbody wind now
460         guide           = old speed of particle
461
462 */
463 void pdDoEffectors(ListBase *lb, float *opco, float *force, float *speed, float cur_time, float loc_time, unsigned int flags)
464 {
465 /*
466         Modifies the force on a particle according to its
467         relation with the effector object
468         Different kind of effectors include:
469                 Forcefields: Gravity-like attractor
470                 (force power is related to the inverse of distance to the power of a falloff value)
471                 Vortex fields: swirling effectors
472                 (particles rotate around Z-axis of the object. otherwise, same relation as)
473                 (Forcefields, but this is not done through a force/acceleration)
474                 Guide: particles on a path
475                 (particles are guided along a curve bezier or old nurbs)
476                 (is independent of other effectors)
477 */
478         Object *ob;
479         pEffectorCache *ec;
480         PartDeflect *pd;
481         float vect_to_vert[3];
482         float f_force, force_vec[3];
483         float *obloc;
484         float distance, force_val, ffall_val;
485         float guidecollect[3], guidedist= 0.0f;
486         int cur_frame;
487         
488         guidecollect[0]= guidecollect[1]= guidecollect[2]=0.0f;
489
490         /* Cycle through collected objects, get total of (1/(gravity_strength * dist^gravity_power)) */
491         /* Check for min distance here? (yes would be cool to add that, ton) */
492         
493         for(ec = lb->first; ec; ec= ec->next) {
494                 /* object effectors were fully checked to be OK to evaluate! */
495                 ob= ec->ob;
496                 pd= ob->pd;
497                         
498                 /* Get IPO force strength and fall off values here */
499                 if (has_ipo_code(ob->ipo, OB_PD_FSTR))
500                         force_val = IPO_GetFloatValue(ob->ipo, OB_PD_FSTR, cur_time);
501                 else 
502                         force_val = pd->f_strength;
503                 
504                 if (has_ipo_code(ob->ipo, OB_PD_FFALL)) 
505                         ffall_val = IPO_GetFloatValue(ob->ipo, OB_PD_FFALL, cur_time);
506                 else 
507                         ffall_val = pd->f_power;
508                         
509                 /* Need to set r.cfra for paths (investigate, ton) (uses ob->ctime now, ton) */
510                 if(ob->ctime!=cur_time) {
511                         cur_frame = G.scene->r.cfra;
512                         G.scene->r.cfra = (int)cur_time;
513                         where_is_object_time(ob, cur_time);
514                         G.scene->r.cfra = cur_frame;
515                 }
516                         
517                 /* use center of object for distance calculus */
518                 obloc= ob->obmat[3];
519                 VECSUB(vect_to_vert, obloc, opco);
520                 distance = VecLength(vect_to_vert);
521                         
522                 if((pd->flag & PFIELD_USEMAX) && distance>pd->maxdist && pd->forcefield != PFIELD_GUIDE)
523                         ;       /* don't do anything */
524                 else if(pd->forcefield == PFIELD_WIND) {
525                         VECCOPY(force_vec, ob->obmat[2]);
526                         
527                         /* wind works harder perpendicular to normal, would be nice for softbody later (ton) */
528                         
529                         /* Limit minimum distance to vertex so that */
530                         /* the force is not too big */
531                         if (distance < 0.001) distance = 0.001f;
532                         f_force = (force_val)*(1/(1000 * (float)pow((double)distance, (double)ffall_val)));
533                         /* this option for softbody only */
534                         if(flags && PE_WIND_AS_SPEED){
535                                 speed[0] -= (force_vec[0] * f_force );
536                                 speed[1] -= (force_vec[1] * f_force );
537                                 speed[2] -= (force_vec[2] * f_force );
538                         }
539                         else{
540                                 force[0] += force_vec[0]*f_force;
541                                 force[1] += force_vec[1]*f_force;
542                                 force[2] += force_vec[2]*f_force;
543                         }
544                 }
545                 else if(pd->forcefield == PFIELD_FORCE) {
546                         
547                         /* only use center of object */
548                         obloc= ob->obmat[3];
549
550                         /* Now calculate the gravitational force */
551                         VECSUB(vect_to_vert, obloc, opco);
552                         distance = VecLength(vect_to_vert);
553
554                         /* Limit minimum distance to vertex so that */
555                         /* the force is not too big */
556                         if (distance < 0.001) distance = 0.001f;
557                         f_force = (force_val)*(1.0/(1000.0 * (float)pow((double)distance, (double)ffall_val)));
558                         force[0] += (vect_to_vert[0] * f_force );
559                         force[1] += (vect_to_vert[1] * f_force );
560                         force[2] += (vect_to_vert[2] * f_force );
561                 }
562                 else if(pd->forcefield == PFIELD_VORTEX) {
563                         float vortexvec[3];
564                         
565                         /* only use center of object */
566                         obloc= ob->obmat[3];
567
568                         /* Now calculate the vortex force */
569                         VECSUB(vect_to_vert, obloc, opco);
570                         distance = VecLength(vect_to_vert);
571
572                         Crossf(force_vec, ob->obmat[2], vect_to_vert);
573                         Normalise(force_vec);
574
575                         /* Limit minimum distance to vertex so that */
576                         /* the force is not too big */
577                         if (distance < 0.001) distance = 0.001f;
578                         f_force = (force_val)*(1.0/(100.0 * (float)pow((double)distance, (double)ffall_val)));
579                         vortexvec[0]= -(force_vec[0] * f_force );
580                         vortexvec[1]= -(force_vec[1] * f_force );
581                         vortexvec[2]= -(force_vec[2] * f_force );
582                         
583                         /* this option for softbody only */
584                         if(flags &&PE_WIND_AS_SPEED) {
585                                 speed[0]+= vortexvec[0];
586                                 speed[1]+= vortexvec[1];
587                                 speed[2]+= vortexvec[2];
588                         }
589                         else {
590                                 /* since vortex alters the speed, we have to correct for the previous vortex result */
591                                 speed[0]+= vortexvec[0] - ec->oldspeed[0];
592                                 speed[1]+= vortexvec[1] - ec->oldspeed[1];
593                                 speed[2]+= vortexvec[2] - ec->oldspeed[2];
594                                 
595                                 VECCOPY(ec->oldspeed, vortexvec);
596                         }
597                 }
598                 else if(pd->forcefield == PFIELD_GUIDE) {
599                         float guidevec[4], guidedir[3];
600                         float mindist= force_val; /* force_val is actually mindist in the UI */
601                         
602                         distance= ec->guide_dist;
603                         
604                         /* WARNING: bails out with continue here */
605                         if((pd->flag & PFIELD_USEMAX) && distance>pd->maxdist) continue;
606                         
607                         /* calculate contribution factor for this guide */
608                         if(distance<=mindist) f_force= 1.0f;
609                         else if(pd->flag & PFIELD_USEMAX) {
610                                 if(distance>pd->maxdist || mindist>=pd->maxdist) f_force= 0.0f;
611                                 else {
612                                         f_force= 1.0f - (distance-mindist)/(pd->maxdist - mindist);
613                                         if(ffall_val!=0.0f)
614                                                 f_force = (float)pow(f_force, ffall_val+1.0);
615                                 }
616                         }
617                         else {
618                                 f_force= 1.0f/(1.0f + distance-mindist);
619                                 if(ffall_val!=0.0f)
620                                         f_force = (float)pow(f_force, ffall_val+1.0);
621                         }
622                         
623                         /* now derive path point from loc_time */
624                         if(pd->flag & PFIELD_GUIDE_PATH_ADD)
625                                 where_on_path(ob, f_force*loc_time*ec->time_scale, guidevec, guidedir);
626                         else
627                                 where_on_path(ob, loc_time*ec->time_scale, guidevec, guidedir);
628                         
629                         VECSUB(guidedir, guidevec, ec->oldloc);
630                         VECCOPY(ec->oldloc, guidevec);
631                         
632                         Mat4Mul3Vecfl(ob->obmat, guidedir);
633                         VecMulf(guidedir, ec->scale);   /* correction for lifetime and speed */
634                         
635                         /* we subtract the speed we gave it previous step */
636                         VECCOPY(guidevec, guidedir);
637                         VECSUB(guidedir, guidedir, ec->oldspeed);
638                         VECCOPY(ec->oldspeed, guidevec);
639                         
640                         /* if it fully contributes, we stop */
641                         if(f_force==1.0) {
642                                 VECCOPY(guidecollect, guidedir);
643                                 guidedist= 1.0f;
644                                 break;
645                         }
646                         else if(guidedist<1.0f) {
647                                 VecMulf(guidedir, f_force);
648                                 VECADD(guidecollect, guidecollect, guidedir);
649                                 guidedist += f_force;
650                         }                                       
651                 }
652         }
653
654         /* all guides are accumulated here */
655         if(guidedist!=0.0f) {
656                 if(guidedist!=1.0f) VecMulf(guidecollect, 1.0f/guidedist);
657                 VECADD(speed, speed, guidecollect);
658         }
659 }
660
661 static void cache_object_vertices(Object *ob)
662 {
663         Mesh *me;
664         MVert *mvert;
665         float *fp;
666         int a;
667         
668         me= ob->data;
669         if(me->totvert==0) return;
670
671         fp= ob->sumohandle= MEM_mallocN(3*sizeof(float)*me->totvert, "cache particles");
672         mvert= me->mvert;
673         a= me->totvert;
674         while(a--) {
675                 VECCOPY(fp, mvert->co);
676                 Mat4MulVecfl(ob->obmat, fp);
677                 mvert++;
678                 fp+= 3;
679         }
680 }
681
682 static int pdDoDeflection(RNG *rng, float opco[3], float npco[3], float opno[3],
683         float npno[3], float life, float force[3], int def_depth,
684         float cur_time, unsigned int par_layer, int *last_object,
685                 int *last_face, int *same_face)
686 {
687         /* Particle deflection code */
688         /* The code is in two sections: the first part checks whether a particle has            */
689         /* intersected a face of a deflector mesh, given its old and new co-ords, opco and npco */
690         /* and which face it hit first                                                          */
691         /* The second part calculates the new co-ordinates given that collision and updates     */
692         /* the new co-ordinates accordingly */
693         Base *base;
694         Object *ob, *deflection_object = NULL;
695         Mesh *def_mesh;
696         MFace *mface, *deflection_face = NULL;
697         float *v1, *v2, *v3, *v4, *vcache=NULL;
698         float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3];
699         float dv1[3] = {0}, dv2[3] = {0}, dv3[3] = {0};
700         float vect_to_int[3], refl_vel[3];
701         float d_intersect_co[3], d_intersect_vect[3], d_nvect[3], d_i_co_above[3];
702         float forcec[3];
703         float k_point3, dist_to_plane;
704         float first_dist, ref_plane_mag;
705         float dk_plane=0, dk_point1=0;
706         float icalctop, icalcbot, n_mag;
707         float mag_iv, x_m,y_m,z_m;
708         float damping, perm_thresh;
709         float perm_val, rdamp_val;
710         int a, deflected=0, deflected_now=0;
711         float t,t2, min_t;
712         float mat[3][3], obloc[3] = {0};
713         int cur_frame;
714         float time_before, time_after;
715         float force_mag_norm;
716         int d_object=0, d_face=0, ds_object=0, ds_face=0;
717
718         first_dist = 200000;
719         min_t = 200000;
720
721         /* The first part of the code, finding the first intersected face*/
722         base= G.scene->base.first;
723         while (base) {
724                 /*Only proceed for mesh object in same layer */
725                 if(base->object->type==OB_MESH && (base->lay & par_layer)) {
726                         ob= base->object;
727                         /* only with deflecting set */
728                         if(ob->pd && ob->pd->deflect) {
729                                 def_mesh= ob->data;
730                         
731                                 d_object = d_object + 1;
732
733                                 d_face = d_face + 1;
734                                 mface= def_mesh->mface;
735                                 a = def_mesh->totface;
736                                 
737                                 
738                                 if(ob->parent==NULL && ob->ipo==NULL) { // static
739                                         if(ob->sumohandle==NULL) cache_object_vertices(ob);
740                                         vcache= ob->sumohandle;
741                                 }
742                                 else {
743                                         /*Find out where the object is at this time*/
744                                         cur_frame = G.scene->r.cfra;
745                                         G.scene->r.cfra = (int)cur_time;
746                                         where_is_object_time(ob, cur_time);
747                                         G.scene->r.cfra = cur_frame;
748                                         
749                                         /*Pass the values from ob->obmat to mat*/
750                                         /*and the location values to obloc           */
751                                         Mat3CpyMat4(mat,ob->obmat);
752                                         obloc[0] = ob->obmat[3][0];
753                                         obloc[1] = ob->obmat[3][1];
754                                         obloc[2] = ob->obmat[3][2];
755                                         vcache= NULL;
756
757                                 }
758                                 
759                                 while (a--) {
760
761                                         if(vcache) {
762                                                 v1= vcache+ 3*(mface->v1);
763                                                 VECCOPY(nv1, v1);
764                                                 v1= vcache+ 3*(mface->v2);
765                                                 VECCOPY(nv2, v1);
766                                                 v1= vcache+ 3*(mface->v3);
767                                                 VECCOPY(nv3, v1);
768                                                 v1= vcache+ 3*(mface->v4);
769                                                 VECCOPY(nv4, v1);
770                                         }
771                                         else {
772                                                 /* Calculate the global co-ordinates of the vertices*/
773                                                 v1= (def_mesh->mvert+(mface->v1))->co;
774                                                 v2= (def_mesh->mvert+(mface->v2))->co;
775                                                 v3= (def_mesh->mvert+(mface->v3))->co;
776                                                 v4= (def_mesh->mvert+(mface->v4))->co;
777         
778                                                 VECCOPY(nv1, v1);
779                                                 VECCOPY(nv2, v2);
780                                                 VECCOPY(nv3, v3);
781                                                 VECCOPY(nv4, v4);
782         
783                                                 /*Apply the objects deformation matrix*/
784                                                 Mat3MulVecfl(mat, nv1);
785                                                 Mat3MulVecfl(mat, nv2);
786                                                 Mat3MulVecfl(mat, nv3);
787                                                 Mat3MulVecfl(mat, nv4);
788         
789                                                 VECADD(nv1, nv1, obloc);
790                                                 VECADD(nv2, nv2, obloc);
791                                                 VECADD(nv3, nv3, obloc);
792                                                 VECADD(nv4, nv4, obloc);
793                                         }
794                                         
795                                         deflected_now = 0;
796
797                                                 
798                                                 
799 //                                      t= 0.5; // this is labda of line, can use it optimize quad intersection
800 // sorry but no .. see below (BM)                                       
801                                         if( LineIntersectsTriangle(opco, npco, nv1, nv2, nv3, &t) ) {
802                                                 if (t < min_t) {
803                                                         deflected = 1;
804                                                         deflected_now = 1;
805                                                 }
806                                         }
807 //                                      else if (mface->v4 && (t>=0.0 && t<=1.0)) {
808 // no, you can't skip testing the other triangle
809 // it might give a smaller t on (close to) the edge .. this is numerics not esoteric maths :)
810 // note: the 2 triangles don't need to share a plane ! (BM)
811                                         if (mface->v4) {
812                                                 if( LineIntersectsTriangle(opco, npco, nv1, nv3, nv4, &t2) ) {
813                                                         if (t2 < min_t) {
814                                                                 deflected = 1;
815                                                                 deflected_now = 2;
816                                                         }
817                                                 }
818                                         }
819                                         
820                                         if ((deflected_now > 0) && ((t < min_t) ||(t2 < min_t))) {
821                         min_t = t;
822                         ds_object = d_object;
823                                                 ds_face = d_face;
824                                                 deflection_object = ob;
825                                                 deflection_face = mface;
826                                                 if (deflected_now==1) {
827                         min_t = t;
828                                                         VECCOPY(dv1, nv1);
829                                                         VECCOPY(dv2, nv2);
830                                                         VECCOPY(dv3, nv3);
831                                                 }
832                                                 else {
833                         min_t = t2;
834                                                         VECCOPY(dv1, nv1);
835                                                         VECCOPY(dv2, nv3);
836                                                         VECCOPY(dv3, nv4);
837                                                 }
838                                         }
839                                         mface++;
840                                 }
841                         }
842                 }
843                 base = base->next;
844         }
845
846
847         /* Here's the point to do the permeability calculation */
848         /* Set deflected to 0 if a random number is below the value */
849         /* Get the permeability IPO here*/
850         if (deflected) {
851                 
852                 if (has_ipo_code(deflection_object->ipo, OB_PD_PERM)) 
853                         perm_val = IPO_GetFloatValue(deflection_object->ipo, OB_PD_PERM, cur_time);
854                 else 
855                         perm_val = deflection_object->pd->pdef_perm;
856
857                 perm_thresh = rng_getFloat(rng) - perm_val;
858                 if (perm_thresh < 0 ) {
859                         deflected = 0;
860                 }
861         }
862
863         /* Now for the second part of the deflection code - work out the new speed */
864         /* and position of the particle if a collision occurred */
865         if (deflected) {
866         VECSUB(edge1, dv1, dv2);
867                 VECSUB(edge2, dv3, dv2);
868                 Crossf(d_nvect, edge2, edge1);
869                 n_mag = Normalise(d_nvect);
870                 dk_plane = INPR(d_nvect, nv1);
871                 dk_point1 = INPR(d_nvect,opco);
872
873                 VECSUB(d_intersect_vect, npco, opco);
874
875                 d_intersect_co[0] = opco[0] + (min_t * (npco[0] - opco[0]));
876                 d_intersect_co[1] = opco[1] + (min_t * (npco[1] - opco[1]));
877                 d_intersect_co[2] = opco[2] + (min_t * (npco[2] - opco[2]));
878                 
879                 d_i_co_above[0] = (d_intersect_co[0] + (0.001f * d_nvect[0]));
880                 d_i_co_above[1] = (d_intersect_co[1] + (0.001f * d_nvect[1]));
881                 d_i_co_above[2] = (d_intersect_co[2] + (0.001f * d_nvect[2]));
882                 mag_iv = Normalise(d_intersect_vect);
883                 VECCOPY(npco, d_intersect_co);
884                 
885                 VECSUB(vect_to_int, opco, d_intersect_co);
886                 first_dist = Normalise(vect_to_int);
887
888                 /* Work out the lengths of time before and after collision*/
889                 time_before = (life*(first_dist / (mag_iv)));
890                 time_after =  (life*((mag_iv - first_dist) / (mag_iv)));
891
892                 /* We have to recalculate what the speed would have been at the */
893                 /* point of collision, not the key frame time */
894                 npno[0]= opno[0] + time_before*force[0];
895                 npno[1]= opno[1] + time_before*force[1];
896                 npno[2]= opno[2] + time_before*force[2];
897
898
899                 /* Reflect the speed vector in the face */
900                 x_m = (2 * npno[0] * d_nvect[0]);
901                 y_m = (2 * npno[1] * d_nvect[1]);
902                 z_m = (2 * npno[2] * d_nvect[2]);
903                 refl_vel[0] = npno[0] - (d_nvect[0] * (x_m + y_m + z_m));
904                 refl_vel[1] = npno[1] - (d_nvect[1] * (x_m + y_m + z_m));
905                 refl_vel[2] = npno[2] - (d_nvect[2] * (x_m + y_m + z_m));
906
907                 /*A random variation in the damping factor........ */
908                 /*Get the IPO values for damping here*/
909                 
910                 if (has_ipo_code(deflection_object->ipo, OB_PD_SDAMP)) 
911                         damping = IPO_GetFloatValue(deflection_object->ipo, OB_PD_SDAMP, cur_time);
912                 else 
913                         damping = deflection_object->pd->pdef_damp;
914                 
915                 if (has_ipo_code(deflection_object->ipo, OB_PD_RDAMP)) 
916                         rdamp_val = IPO_GetFloatValue(deflection_object->ipo, OB_PD_RDAMP, cur_time);
917                 else 
918                         rdamp_val = deflection_object->pd->pdef_rdamp;
919
920                 damping = damping + ((1.0f - damping) * rng_getFloat(rng) *rdamp_val);
921                 damping = damping * damping;
922         ref_plane_mag = INPR(refl_vel,d_nvect);
923
924                 if (damping > 0.999) damping = 0.999f;
925
926                 /* Now add in the damping force - only damp in the direction of */
927                 /* the faces normal vector */
928                 npno[0] = (refl_vel[0] - (d_nvect[0] * ref_plane_mag * damping));
929                 npno[1] = (refl_vel[1] - (d_nvect[1] * ref_plane_mag * damping));
930                 npno[2] = (refl_vel[2] - (d_nvect[2] * ref_plane_mag * damping));
931
932                 /* Now reset opno */
933                 VECCOPY(opno,npno);
934                 VECCOPY(forcec, force);
935
936                 /* If the particle has bounced more than four times on the same */
937                 /* face within this cycle (depth > 4, same face > 4 )           */
938                 /* Then set the force to be only that component of the force    */
939                 /* in the same direction as the face normal                     */
940                 /* i.e. subtract the component of the force in the direction    */
941                 /* of the face normal from the actual force                     */
942                 if ((ds_object == *last_object) && (ds_face == *last_face)) {
943                         /* Increment same_face */
944                         *same_face = *same_face + 1;
945                         if ((*same_face > 3) && (def_depth > 3)) {
946                 force_mag_norm = INPR(forcec, d_nvect);
947                 forcec[0] = forcec[0] - (d_nvect[0] * force_mag_norm);
948                 forcec[1] = forcec[1] - (d_nvect[1] * force_mag_norm);
949                 forcec[2] = forcec[2] - (d_nvect[2] * force_mag_norm);
950                         }
951                 }
952                 else *same_face = 1;
953
954                 *last_object = ds_object;
955                 *last_face = ds_face;
956
957                 /* We have the particles speed at the point of collision    */
958                 /* Now we want the particles speed at the current key frame */
959
960                 npno[0]= npno[0] + time_after*forcec[0];
961                 npno[1]= npno[1] + time_after*forcec[1];
962                 npno[2]= npno[2] + time_after*forcec[2];
963
964                 /* Now we have to recalculate pa->co for the remainder*/
965                 /* of the time since the intersect*/
966                 npco[0]= npco[0] + time_after*npno[0];
967                 npco[1]= npco[1] + time_after*npno[1];
968                 npco[2]= npco[2] + time_after*npno[2];
969
970                 /* And set the old co-ordinates back to the point just above the intersection */
971                 VECCOPY(opco, d_i_co_above);
972
973                 /* Finally update the time */
974                 life = time_after;
975                 cur_time += time_before;
976
977                 /* The particle may have fallen through the face again by now!!*/
978                 /* So check if the particle has changed sides of the plane compared*/
979                 /* the co-ordinates at the last keyframe*/
980                 /* But only do this as a last resort, if we've got to the end of the */
981                 /* number of collisions allowed */
982                 if (def_depth==9) {
983                         k_point3 = INPR(d_nvect,npco);
984                         if (((dk_plane > k_point3) && (dk_plane < dk_point1))||((dk_plane < k_point3) && (dk_plane > dk_point1))) {
985
986                                 /* Yup, the pesky particle may have fallen through a hole!!! */
987                 /* So we'll cheat a bit and move the particle along the normal vector */
988                 /* until it's just the other side of the plane */
989                 icalctop = (dk_plane - d_nvect[0]*npco[0] - d_nvect[1]*npco[1] - d_nvect[2]*npco[2]);
990                 icalcbot = (d_nvect[0]*d_nvect[0] + d_nvect[1]*d_nvect[1] + d_nvect[2]*d_nvect[2]);
991                 dist_to_plane = icalctop / icalcbot;
992
993                 /*  Now just increase the distance a little to place */
994                 /* the point the other side of the plane */
995                 dist_to_plane *= 1.1f;
996                 npco[0]= npco[0] + (dist_to_plane * d_nvect[0]);
997                 npco[1]= npco[1] + (dist_to_plane * d_nvect[1]);
998                 npco[2]= npco[2] + (dist_to_plane * d_nvect[2]);
999
1000                         }
1001                 }
1002         }
1003         return deflected;
1004 }
1005
1006 /*
1007         rng= random number generator
1008         ob = object that spawns the particles
1009         depth = for fireworks
1010         nr = index nr of current particle
1011         paf = the particle system
1012         part = current particle
1013         force = force vector
1014         deform = flag to indicate lattice deform
1015  */
1016 static void make_particle_keys(RNG *rng, Object *ob, int depth, int nr, PartEff *paf, Particle *part, float *force, int deform, MTex *mtex, ListBase *effectorbase)
1017 {
1018         Particle *pa, *opa = NULL;
1019         float damp, deltalife, life;
1020         float cur_time;
1021         float opco[3], opno[3], npco[3], npno[3], new_force[3], new_speed[3];
1022         int b, rt1, rt2, deflected, deflection, finish_defs, def_count;
1023         int last_ob, last_fc, same_fc;
1024
1025         damp= 1.0f-paf->damp;
1026         pa= part;
1027
1028         /* start speed: random */
1029         if(paf->randfac!=0.0) {
1030                 pa->no[0]+= paf->randfac*(rng_getFloat(rng) - 0.5f);
1031                 pa->no[1]+= paf->randfac*(rng_getFloat(rng) - 0.5f);
1032                 pa->no[2]+= paf->randfac*(rng_getFloat(rng) - 0.5f);
1033         }
1034
1035         /* start speed: texture */
1036         if(mtex && paf->texfac!=0.0) {
1037                 particle_tex(mtex, paf, pa->co, pa->no);
1038         }
1039
1040         /* effectors here? */
1041         if(effectorbase)
1042                 precalc_effectors(ob, paf, pa, effectorbase);
1043         
1044         if(paf->totkey>1) deltalife= pa->lifetime/(paf->totkey-1);
1045         else deltalife= pa->lifetime;
1046
1047         /* longer lifetime results in longer distance covered */
1048         VecMulf(pa->no, deltalife);
1049         
1050         opa= pa;
1051         pa++;
1052
1053         for(b=1; b<paf->totkey; b++) {
1054
1055                 /* new time */
1056                 pa->time= opa->time+deltalife;
1057                 cur_time = pa->time;
1058
1059                 /* set initial variables                                */
1060                 VECCOPY(opco, opa->co);
1061                 VECCOPY(new_force, force);
1062                 VECCOPY(new_speed, opa->no);
1063                 VecMulf(new_speed, 1.0f/deltalife);
1064                 //new_speed[0] = new_speed[1] = new_speed[2] = 0.0f;
1065
1066                 /* handle differences between static (local coords, fixed frame) and dynamic */
1067                 if(effectorbase) {
1068                         float loc_time= ((float)b)/(float)(paf->totkey-1);
1069                         
1070                         if(paf->flag & PAF_STATIC) {
1071                                 float opco1[3], new_force1[3];
1072                                 
1073                                 /* move co and force to global coords */
1074                                 VECCOPY(opco1, opco);
1075                                 Mat4MulVecfl(ob->obmat, opco1);
1076                                 VECCOPY(new_force1, new_force);
1077                                 Mat4Mul3Vecfl(ob->obmat, new_force1);
1078                                 Mat4Mul3Vecfl(ob->obmat, new_speed);
1079                                 
1080                                 cur_time = G.scene->r.cfra;
1081                                 
1082                                 /* force fields */
1083                                 pdDoEffectors(effectorbase, opco1, new_force1, new_speed, cur_time, loc_time, 0);
1084                                 
1085                                 /* move co, force and newspeed back to local */
1086                                 VECCOPY(opco, opco1);
1087                                 Mat4MulVecfl(ob->imat, opco);
1088                                 VECCOPY(new_force, new_force1);
1089                                 Mat4Mul3Vecfl(ob->imat, new_force);
1090                                 Mat4Mul3Vecfl(ob->imat, new_speed);
1091                         }
1092                         else {
1093                                  /* force fields */
1094                                 pdDoEffectors(effectorbase, opco, new_force, new_speed, cur_time, loc_time, 0);
1095                         }
1096                 }
1097                 
1098                 /* new speed */
1099                 pa->no[0]= deltalife * (new_speed[0] + new_force[0]);
1100                 pa->no[1]= deltalife * (new_speed[1] + new_force[1]);
1101                 pa->no[2]= deltalife * (new_speed[2] + new_force[2]);
1102                 
1103                 /* new location */
1104                 pa->co[0]= opa->co[0] + pa->no[0];
1105                 pa->co[1]= opa->co[1] + pa->no[1];
1106                 pa->co[2]= opa->co[2] + pa->no[2];
1107
1108                 /* Particle deflection code */
1109                 if((paf->flag & PAF_STATIC)==0) {
1110                         deflection = 0;
1111                         finish_defs = 1;
1112                         def_count = 0;
1113
1114                         VECCOPY(opno, opa->no);
1115                         VECCOPY(npco, pa->co);
1116                         VECCOPY(npno, pa->no);
1117
1118                         life = deltalife;
1119                         cur_time -= deltalife;
1120
1121                         last_ob = -1;
1122                         last_fc = -1;
1123                         same_fc = 0;
1124
1125                         /* First call the particle deflection check for the particle moving   */
1126                         /* between the old co-ordinates and the new co-ordinates              */
1127                         /* If a deflection occurs, call the code again, this time between the */
1128                         /* intersection point and the updated new co-ordinates                */
1129                         /* Bail out if we've done the calculation 10 times - this seems ok     */
1130                         /* for most scenes I've tested */
1131                         while (finish_defs) {
1132                                 deflected =  pdDoDeflection(rng, opco, npco, opno, npno, life, new_force,
1133                                                                 def_count, cur_time, ob->lay,
1134                                                                 &last_ob, &last_fc, &same_fc);
1135                                 if (deflected) {
1136                                         def_count = def_count + 1;
1137                                         deflection = 1;
1138                                         if (def_count==10) finish_defs = 0;
1139                                 }
1140                                 else {
1141                                         finish_defs = 0;
1142                                 }
1143                         }
1144
1145                         /* Only update the particle positions and speed if we had a deflection */
1146                         if (deflection) {
1147                                 pa->co[0] = npco[0];
1148                                 pa->co[1] = npco[1];
1149                                 pa->co[2] = npco[2];
1150                                 pa->no[0] = npno[0];
1151                                 pa->no[1] = npno[1];
1152                                 pa->no[2] = npno[2];
1153                         }
1154                 }
1155                 
1156                 /* speed: texture */
1157                 if(mtex && paf->texfac!=0.0) {
1158                         particle_tex(mtex, paf, pa->co, pa->no);
1159                 }
1160                 if(damp!=1.0) {
1161                         pa->no[0]*= damp;
1162                         pa->no[1]*= damp;
1163                         pa->no[2]*= damp;
1164                 }
1165                 
1166                 opa= pa;
1167                 pa++;
1168                 /* opa is used later on too! */
1169         }
1170
1171         if(deform) {
1172                 /* deform all keys */
1173                 pa= part;
1174                 b= paf->totkey;
1175                 while(b--) {
1176                         calc_latt_deform(pa->co, 1.0f);
1177                         pa++;
1178                 }
1179         }
1180         
1181         /* the big multiplication */
1182         if(depth<PAF_MAXMULT && paf->mult[depth]!=0.0) {
1183                 
1184                 /* new 'child' emerges from an average 'mult' part from 
1185                         the particles */
1186                 damp = (float)nr;
1187                 rt1= (int)(damp*paf->mult[depth]);
1188                 rt2= (int)((damp+1.0)*paf->mult[depth]);
1189                 if(rt1!=rt2) {
1190                         
1191                         for(b=0; b<paf->child[depth]; b++) {
1192                                 pa= new_particle(paf);
1193                                 *pa= *opa;
1194                                 pa->lifetime= paf->life[depth];
1195                                 if(paf->randlife!=0.0) {
1196                                         pa->lifetime*= 1.0f + paf->randlife*(rng_getFloat(rng) - 0.5f);
1197                                 }
1198                                 pa->mat_nr= paf->mat[depth];
1199
1200                                 make_particle_keys(rng, ob, depth+1, b, paf, pa, force, deform, mtex, effectorbase);
1201                         }
1202                 }
1203         }
1204 }
1205
1206 static void init_mv_jit(float *jit, int num, int seed2)
1207 {
1208         RNG *rng;
1209         float *jit2, x, rad1, rad2, rad3;
1210         int i, num2;
1211
1212         if(num==0) return;
1213
1214         rad1= (float)(1.0/sqrt((float)num));
1215         rad2= (float)(1.0/((float)num));
1216         rad3= (float)sqrt((float)num)/((float)num);
1217
1218         rng = rng_new(31415926 + num + seed2);
1219         x= 0;
1220         num2 = 2 * num;
1221         for(i=0; i<num2; i+=2) {
1222         
1223                 jit[i]= x + rad1*(0.5f - rng_getFloat(rng));
1224                 jit[i+1]= i/(2.0f*num) + rad1*(0.5f - rng_getFloat(rng));
1225                 
1226                 jit[i]-= (float)floor(jit[i]);
1227                 jit[i+1]-= (float)floor(jit[i+1]);
1228                 
1229                 x+= rad3;
1230                 x -= (float)floor(x);
1231         }
1232
1233         jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
1234
1235         for (i=0 ; i<4 ; i++) {
1236                 BLI_jitterate1(jit, jit2, num, rad1);
1237                 BLI_jitterate1(jit, jit2, num, rad1);
1238                 BLI_jitterate2(jit, jit2, num, rad2);
1239         }
1240         MEM_freeN(jit2);
1241         rng_free(rng);
1242 }
1243
1244 #define JIT_RAND        32
1245
1246 /* for a position within a face, tot is total amount of faces */
1247 static void give_mesh_particle_coord(PartEff *paf, VeNoCo *noco, MFace *mface, int partnr, int subnr, float *co, float *no)
1248 {
1249         static float *jit= NULL;
1250         static float *trands= NULL;
1251         static int jitlevel= 1;
1252         float *v1, *v2, *v3, *v4;
1253         float u, v;
1254         float *n1, *n2, *n3, *n4;
1255         
1256         /* free signal */
1257         if(paf==NULL) {
1258                 if(jit) MEM_freeN(jit);
1259                 jit= NULL;
1260                 if(trands) MEM_freeN(trands);
1261                 trands= NULL;
1262                 return;
1263         }
1264         
1265         /* first time initialize jitter or trand, partnr then is total amount of particles, subnr total amount of faces */
1266         if(trands==NULL && jit==NULL) {
1267                 RNG *rng = rng_new(31415926 + paf->seed);
1268                 int i, tot;
1269
1270                 if(paf->flag & PAF_TRAND)
1271                         tot= partnr;
1272                 else 
1273                         tot= JIT_RAND;  /* arbitrary... allows JIT_RAND times more particles in a face for jittered distro */
1274                         
1275                 trands= MEM_callocN(2+2*tot*sizeof(float), "trands");
1276                 for(i=0; i<tot; i++) {
1277                         trands[2*i]= rng_getFloat(rng);
1278                         trands[2*i+1]= rng_getFloat(rng);
1279                 }
1280                 rng_free(rng);
1281
1282                 if((paf->flag & PAF_TRAND)==0) {
1283                         jitlevel= paf->userjit;
1284                         
1285                         if(jitlevel == 0) {
1286                                 jitlevel= partnr/subnr;
1287                                 if(paf->flag & PAF_EDISTR) jitlevel*= 2;        /* looks better in general, not very scietific */
1288                                 if(jitlevel<3) jitlevel= 3;
1289                                 if(jitlevel>100) jitlevel= 100;
1290                         }
1291                         
1292                         jit= MEM_callocN(2+ jitlevel*2*sizeof(float), "jit");
1293                         init_mv_jit(jit, jitlevel, paf->seed);
1294                         BLI_array_randomize(jit, 2*sizeof(float), jitlevel, paf->seed); /* for custom jit or even distribution */
1295                 }
1296                 return;
1297         }
1298         
1299         if(paf->flag & PAF_TRAND) {
1300                 u= trands[2*partnr];
1301                 v= trands[2*partnr+1];
1302         }
1303         else {
1304                 /* jittered distribution gets fixed random offset */
1305                 if(subnr>=jitlevel) {
1306                         int jitrand= (subnr/jitlevel) % JIT_RAND;
1307                 
1308                         subnr %= jitlevel;
1309                         u= jit[2*subnr] + trands[2*jitrand];
1310                         v= jit[2*subnr+1] + trands[2*jitrand+1];
1311                         if(u > 1.0f) u-= 1.0f;
1312                         if(v > 1.0f) v-= 1.0f;
1313                 }
1314                 else {
1315                         u= jit[2*subnr];
1316                         v= jit[2*subnr+1];
1317                 }
1318         }
1319         
1320         v1= (noco+(mface->v1))->co;
1321         v2= (noco+(mface->v2))->co;
1322         v3= (noco+(mface->v3))->co;
1323         n1= (noco+(mface->v1))->no;
1324         n2= (noco+(mface->v2))->no;
1325         n3= (noco+(mface->v3))->no;
1326         
1327         if(mface->v4) {
1328                 float uv= u*v;
1329                 float muv= (1.0f-u)*(v);
1330                 float umv= (u)*(1.0f-v);
1331                 float mumv= (1.0f-u)*(1.0f-v);
1332                 
1333                 v4= (noco+(mface->v4))->co;
1334                 n4= (noco+(mface->v4))->no;
1335                 
1336                 co[0]= mumv*v1[0] + muv*v2[0] + uv*v3[0] + umv*v4[0];
1337                 co[1]= mumv*v1[1] + muv*v2[1] + uv*v3[1] + umv*v4[1];
1338                 co[2]= mumv*v1[2] + muv*v2[2] + uv*v3[2] + umv*v4[2];
1339
1340                 no[0]= mumv*n1[0] + muv*n2[0] + uv*n3[0] + umv*n4[0];
1341                 no[1]= mumv*n1[1] + muv*n2[1] + uv*n3[1] + umv*n4[1];
1342                 no[2]= mumv*n1[2] + muv*n2[2] + uv*n3[2] + umv*n4[2];
1343         }
1344         else {
1345                 /* mirror triangle uv coordinates when on other side */
1346                 if(u + v > 1.0f) {
1347                         u= 1.0f-u;
1348                         v= 1.0f-v;
1349                 }
1350                 co[0]= v1[0] + u*(v3[0]-v1[0]) + v*(v2[0]-v1[0]);
1351                 co[1]= v1[1] + u*(v3[1]-v1[1]) + v*(v2[1]-v1[1]);
1352                 co[2]= v1[2] + u*(v3[2]-v1[2]) + v*(v2[2]-v1[2]);
1353                 
1354                 no[0]= n1[0] + u*(n3[0]-n1[0]) + v*(n2[0]-n1[0]);
1355                 no[1]= n1[1] + u*(n3[1]-n1[1]) + v*(n2[1]-n1[1]);
1356                 no[2]= n1[2] + u*(n3[2]-n1[2]) + v*(n2[2]-n1[2]);
1357         }
1358 }
1359
1360
1361 /* Gets a MDeformVert's weight in group (0 if not in group) */
1362 /* note; this call could be in mesh.c or deform.c, but OK... it's in armature.c too! (ton) */
1363 static float vert_weight(MDeformVert *dvert, int group)
1364 {
1365         MDeformWeight *dw;
1366         int i;
1367         
1368         if(dvert) {
1369                 dw= dvert->dw;
1370                 for(i= dvert->totweight; i>0; i--, dw++) {
1371                         if(dw->def_nr == group) return dw->weight;
1372                         if(i==1) break; /*otherwise dw will point to somewhere it shouldn't*/
1373                 }
1374         }
1375         return 0.0;
1376 }
1377
1378 /* Gets a faces average weight in a group, helper for below, face and weights are always set */
1379 static float face_weight(MFace *face, float *weights)
1380 {
1381         float tweight;
1382         
1383         tweight = weights[face->v1] + weights[face->v2] + weights[face->v3];
1384         
1385         if(face->v4) {
1386                 tweight += weights[face->v4];
1387                 tweight /= 4.0;
1388         }
1389         else {
1390                 tweight /= 3.0;
1391         }
1392
1393         return tweight;
1394 }
1395
1396 /* helper function for build_particle_system() */
1397 static void make_weight_tables(PartEff *paf, Mesh *me, int totpart, VeNoCo *vertlist, int totvert, MFace *facelist, int totface, float **vweights, float **fweights)
1398 {
1399         MFace *mface;
1400         float *foweights=NULL, *voweights=NULL;
1401         float totvweight=0.0f, totfweight=0.0f;
1402         int a;
1403         
1404         if((paf->flag & PAF_FACE)==0) totface= 0;
1405         
1406         /* collect emitting vertices & faces if vert groups used */
1407         if(paf->vertgroup && me->dvert) {
1408                 
1409                 /* allocate weights array for all vertices, also for lookup of faces later on. note it's a malloc */
1410                 *vweights= voweights= MEM_mallocN( totvert*sizeof(float), "pafvoweights" );
1411                 totvweight= 0.0f;
1412                 for(a=0; a<totvert; a++) {
1413                         voweights[a]= vert_weight(me->dvert+a, paf->vertgroup-1);
1414                         totvweight+= voweights[a];
1415                 }
1416                 
1417                 if(totface) {
1418                         /* allocate weights array for faces, note it's a malloc */
1419                         *fweights= foweights= MEM_mallocN(totface*sizeof(float), "paffoweights" );
1420                         for(a=0, mface=facelist; a<totface; a++, mface++) {
1421                                 foweights[a] = face_weight(mface, voweights);
1422                         }
1423                 }
1424         }
1425         
1426         /* make weights for faces or for even area distribution */
1427         if(totface && (paf->flag & PAF_EDISTR)) {
1428                 float maxfarea= 0.0f, curfarea;
1429                 
1430                 /* two cases for area distro, second case we already have group weights */
1431                 if(foweights==NULL) {
1432                         /* allocate weights array for faces, note it's a malloc */
1433                         *fweights= foweights= MEM_mallocN(totface*sizeof(float), "paffoweights" );
1434                         
1435                         for(a=0, mface=facelist; a<totface; a++, mface++) {
1436                                 if (mface->v4)
1437                                         curfarea= AreaQ3Dfl(vertlist[mface->v1].co, vertlist[mface->v2].co, vertlist[mface->v3].co, vertlist[mface->v4].co);
1438                                 else
1439                                         curfarea= AreaT3Dfl(vertlist[mface->v1].co,  vertlist[mface->v2].co, vertlist[mface->v3].co);
1440                                 if(curfarea>maxfarea)
1441                                         maxfarea = curfarea;
1442                                 foweights[a]= curfarea;
1443                         }
1444                 }
1445                 else {
1446                         for(a=0, mface=facelist; a<totface; a++, mface++) {
1447                                 if(foweights[a]!=0.0f) {
1448                                         if (mface->v4)
1449                                                 curfarea= AreaQ3Dfl(vertlist[mface->v1].co, vertlist[mface->v2].co, vertlist[mface->v3].co, vertlist[mface->v4].co);
1450                                         else
1451                                                 curfarea= AreaT3Dfl(vertlist[mface->v1].co,  vertlist[mface->v2].co, vertlist[mface->v3].co);
1452                                         if(curfarea>maxfarea)
1453                                                 maxfarea = curfarea;
1454                                         foweights[a]*= curfarea;
1455                                 }
1456                         }
1457                 }
1458                 
1459                 /* normalize weights for max face area, calculate tot */
1460                 if(maxfarea!=0.0f) {
1461                         maxfarea= 1.0f/maxfarea;
1462                         for(a=0; a< totface; a++) {
1463                                 if(foweights[a]!=0.0) {
1464                                         foweights[a] *= maxfarea;
1465                                         totfweight+= foweights[a];
1466                                 }
1467                         }
1468                 }
1469         }
1470         else if(foweights) {
1471                 /* only add totfweight value */
1472                 for(a=0; a< totface; a++) {
1473                         if(foweights[a]!=0.0) {
1474                                 totfweight+= foweights[a];
1475                         }
1476                 }
1477         }
1478         
1479         /* if weight arrays, we turn these arrays into the amount of particles */
1480         if(totvert && voweights) {
1481                 float mult= (float)totpart/totvweight;
1482                 
1483                 for(a=0; a< totvert; a++) {
1484                         if(voweights[a]!=0.0)
1485                                 voweights[a] *= mult;
1486                 }
1487         }
1488         
1489         if(totface && foweights) {
1490                 float mult= (float)totpart/totfweight;
1491                 
1492                 for(a=0; a< totface; a++) {
1493                         if(foweights[a]!=0.0)
1494                                 foweights[a] *= mult;
1495                 }
1496         }
1497 }
1498
1499 /* helper function for build_particle_system() */
1500 static void make_length_tables(PartEff *paf, Mesh *me, int totvert, MFace *facelist, int totface, float **vlengths, float **flengths)
1501 {
1502         MFace *mface;
1503         float *folengths=NULL, *volengths=NULL;
1504         int a;
1505         
1506         if((paf->flag & PAF_FACE)==0) totface= 0;
1507         
1508         /* collect emitting vertices & faces if vert groups used */
1509         if(paf->vertgroup_v && me->dvert) {
1510                 
1511                 /* allocate lengths array for all vertices, also for lookup of faces later on. note it's a malloc */
1512                 *vlengths= volengths= MEM_mallocN( totvert*sizeof(float), "pafvolengths" );
1513                 for(a=0; a<totvert; a++) {
1514                         volengths[a]= vert_weight(me->dvert+a, paf->vertgroup_v-1);
1515                 }
1516                 
1517                 if(totface) {
1518                         /* allocate lengths array for faces, note it's a malloc */
1519                         *flengths= folengths= MEM_mallocN(totface*sizeof(float), "paffolengths" );
1520                         for(a=0, mface=facelist; a<totface; a++, mface++) {
1521                                 folengths[a] = face_weight(mface, volengths);
1522                         }
1523                 }
1524         }
1525 }
1526
1527
1528 /* for paf start to end, store all matrices for objects */
1529 typedef struct pMatrixCache {
1530         float obmat[4][4];
1531         float imat[3][3];
1532 } pMatrixCache;
1533
1534
1535 /* WARN: this function stores data in ob->id.idnew! */
1536 static pMatrixCache *cache_object_matrices(Object *ob, int start, int end)
1537 {
1538         pMatrixCache *mcache, *mc;
1539         Group *group= NULL;
1540         Object *obcopy;
1541         Base *base;
1542         float framelenold, cfrao, sfo;
1543         
1544         /* object can be linked in group... stupid exception */
1545         if(NULL==object_in_scene(ob, G.scene))
1546                 group= find_group(ob);
1547         
1548         mcache= mc= MEM_mallocN( (end-start+1)*sizeof(pMatrixCache), "ob matrix cache");
1549         
1550         framelenold= G.scene->r.framelen;
1551         G.scene->r.framelen= 1.0f;
1552         cfrao= G.scene->r.cfra;
1553         sfo= ob->sf;
1554         ob->sf= 0.0f;
1555
1556         /* clear storage */
1557         for(obcopy= G.main->object.first; obcopy; obcopy= obcopy->id.next) 
1558                 obcopy->id.newid= NULL;
1559         
1560         /* all objects get tagged recalc that influence this object (does group too) */
1561         /* another hack; while transform you cannot call this, it sets own recalc flags */
1562         if(G.moving==0) {
1563                 ob->recalc |= OB_RECALC_OB; /* make sure a recalc gets flushed */
1564                 DAG_object_update_flags(G.scene, ob, -1);
1565         }
1566         
1567         for(G.scene->r.cfra= start; G.scene->r.cfra<=end; G.scene->r.cfra++, mc++) {
1568                 
1569                 if(group) {
1570                         GroupObject *go;
1571
1572                         for(go= group->gobject.first; go; go= go->next) {
1573                                 if(go->ob->recalc) {
1574                                         where_is_object(go->ob);
1575                                         
1576                                         do_ob_key(go->ob);
1577                                         if(go->ob->type==OB_ARMATURE) {
1578                                                 do_all_pose_actions(go->ob);    // only does this object actions
1579                                                 where_is_pose(go->ob);
1580                                         }
1581                                 }
1582                         }
1583                 }
1584                 else {
1585                         for(base= G.scene->base.first; base; base= base->next) {
1586                                 if(base->object->recalc) {
1587                                         if(base->object->id.newid==NULL)
1588                                                 base->object->id.newid= MEM_dupallocN(base->object);
1589                                         
1590                                         where_is_object(base->object);
1591                                         
1592                                         do_ob_key(base->object);
1593                                         if(base->object->type==OB_ARMATURE) {
1594                                                 do_all_pose_actions(base->object);      // only does this object actions
1595                                                 where_is_pose(base->object);
1596                                         }
1597                                 }
1598                         }
1599                 }               
1600                 Mat4CpyMat4(mc->obmat, ob->obmat);
1601                 Mat4Invert(ob->imat, ob->obmat);
1602                 Mat3CpyMat4(mc->imat, ob->imat);
1603                 Mat3Transp(mc->imat);
1604         }
1605         
1606         /* restore */
1607         G.scene->r.cfra= cfrao;
1608         G.scene->r.framelen= framelenold;
1609         ob->sf= sfo;
1610         
1611         if(group) {
1612                 GroupObject *go;
1613                 
1614                 for(go= group->gobject.first; go; go= go->next) {
1615                         if(go->ob->recalc) {
1616                                 where_is_object(go->ob);
1617                                 
1618                                 do_ob_key(go->ob);
1619                                 if(go->ob->type==OB_ARMATURE) {
1620                                         do_all_pose_actions(go->ob);    // only does this object actions
1621                                         where_is_pose(go->ob);
1622                                 }
1623                         }
1624                 }
1625         }
1626         else {
1627                 for(base= G.scene->base.first; base; base= base->next) {
1628                         if(base->object->recalc) {
1629                                 
1630                                 if(base->object->id.newid) {
1631                                         obcopy= (Object *)base->object->id.newid;
1632                                         *(base->object) = *(obcopy); 
1633                                         MEM_freeN(obcopy);
1634                                         base->object->id.newid= NULL;
1635                                 }
1636                                 
1637                                 do_ob_key(base->object);
1638                                 if(base->object->type==OB_ARMATURE) {
1639                                         do_all_pose_actions(base->object);      // only does this object actions
1640                                         where_is_pose(base->object);
1641                                 }
1642                         }
1643                 }
1644         }       
1645         return mcache;
1646 }
1647
1648 /* for fluidsim win32 debug messages */
1649 #if defined(WIN32) && (!(defined snprintf))
1650 #define snprintf _snprintf
1651 #endif
1652
1653 /* main particle building function 
1654    one day particles should become dynamic (realtime) with the current method as a 'bake' (ton) */
1655 void build_particle_system(Object *ob)
1656 {
1657         RNG *rng;
1658         PartEff *paf;
1659         Particle *pa;
1660         Mesh *me;
1661         Base *base;
1662         MTex *mtexmove=0, *mtextime=0;
1663         Material *ma;
1664         MFace *facelist= NULL;
1665         pMatrixCache *mcache=NULL, *mcnow, *mcprev;
1666         ListBase *effectorbase;
1667         VeNoCo *vertexcosnos;
1668         double startseconds= PIL_check_seconds_timer();
1669         float ftime, dtime, force[3], vec[3], fac, co[3], no[3];
1670         float *voweights= NULL, *foweights= NULL, maxw=1.0f;
1671         float *volengths= NULL, *folengths= NULL;
1672         int deform=0, a, totpart, paf_sta, paf_end;
1673         int waitcursor_set= 0, totvert, totface, curface, curvert;
1674         int readMask, activeParts, fileParts;
1675         
1676         /* return conditions */
1677         if(ob->type!=OB_MESH) return;
1678         me= ob->data;
1679
1680         paf= give_parteff(ob);
1681         if(paf==NULL) return;
1682         
1683         if(G.rendering==0 && paf->disp==0) return;
1684         
1685         if(paf->keys) MEM_freeN(paf->keys);     /* free as early as possible, for returns */
1686         paf->keys= NULL;
1687         
1688         printf("build particles\n");
1689         
1690         /* fluid sim particle import handling, actual loading */
1691         #ifndef DISABLE_ELBEEM
1692         if( (1) && (ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) && 
1693             (ob->fluidsimSettings) && 
1694                   (ob->fluidsimSettings->type == OB_FLUIDSIM_PARTICLE)) {
1695                 char *suffix  = "fluidsurface_particles_#";
1696                 char *suffix2 = ".gz";
1697                 char filename[256];
1698                 char debugStrBuffer[256];
1699                 int  curFrame = G.scene->r.cfra -1; // warning - sync with derived mesh fsmesh loading
1700                 int  j, numFileParts;
1701                 gzFile gzf;
1702                 float vel[3];
1703
1704                 if(ob==G.obedit) { // off...
1705                         paf->totpart = 1;
1706                         return;
1707                 }
1708
1709                 // ok, start loading
1710                 strcpy(filename, ob->fluidsimSettings->surfdataPath);
1711                 strcat(filename, suffix);
1712                 BLI_convertstringcode(filename, G.sce, curFrame); // fixed #frame-no 
1713                 strcat(filename, suffix2);
1714
1715                 gzf = gzopen(filename, "rb");
1716                 if (!gzf) {
1717                         snprintf(debugStrBuffer,256,"readFsPartData::error - Unable to open file for reading '%s' \n", filename); 
1718                         //elbeemDebugOut(debugStrBuffer);
1719                         paf->totpart = 1;
1720                         return;
1721                 }
1722
1723                 gzread(gzf, &totpart, sizeof(totpart));
1724                 numFileParts = totpart;
1725                 totpart = (G.rendering)?totpart:(paf->disp*totpart)/100;
1726                 paf->totpart= totpart;
1727                 paf->totkey= 1;
1728                 /* initialize particles */
1729                 new_particle(paf); 
1730                 ftime = 0.0; // unused...
1731
1732                 // set up reading mask
1733                 readMask = ob->fluidsimSettings->typeFlags;
1734                 activeParts=0;
1735                 fileParts=0;
1736                 
1737                 for(a=0; a<totpart; a++) {
1738                         int ptype=0;
1739                         short shsize=0;
1740                         float convertSize=0.0;
1741                         gzread(gzf, &ptype, sizeof( ptype )); 
1742                         if(ptype&readMask) {
1743                                 activeParts++;
1744                                 pa= new_particle(paf);
1745                                 pa->time= ftime;
1746                                 pa->lifetime= ftime + 10000.; // add large number to make sure they are displayed, G.scene->r.efra +1.0;
1747                                 pa->co[0] = 0.0;
1748                                 pa->co[1] = 
1749                                 pa->co[2] = 1.0*(float)a / (float)totpart;
1750                                 pa->no[0] = pa->no[1] = pa->no[2] = 0.0;
1751                                 pa->mat_nr= paf->omat;
1752                                 gzread(gzf, &convertSize, sizeof( float )); 
1753                                 // convert range of  1.0-10.0 to shorts 1000-10000)
1754                                 shsize = (short)(convertSize*1000.0);
1755                                 pa->rt = shsize;
1756                                 //if(a<200) fprintf(stderr,"SREAD %f %d %d \n",convertSize,shsize,pa->rt);
1757
1758                                 for(j=0; j<3; j++) {
1759                                         float wrf;
1760                                         gzread(gzf, &wrf, sizeof( wrf )); 
1761                                         pa->co[j] = wrf;
1762                                         //fprintf(stderr,"Rj%d ",j);
1763                                 }
1764                                 for(j=0; j<3; j++) {
1765                                         float wrf;
1766                                         gzread(gzf, &wrf, sizeof( wrf )); 
1767                                         vel[j] = wrf;
1768                                 }
1769                                 //if(a<25) fprintf(stderr,"FSPARTICLE debug set %s , a%d = %f,%f,%f , life=%f \n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
1770                         } else {
1771                                 // skip...
1772                                 for(j=0; j<2*3+1; j++) {
1773                                         float wrf; gzread(gzf, &wrf, sizeof( wrf )); 
1774                                 }
1775                         }
1776                         fileParts++;
1777                 }
1778                 gzclose( gzf );
1779
1780                 totpart = paf->totpart = activeParts;
1781                 snprintf(debugStrBuffer,256,"readFsPartData::done - particles:%d, active:%d, file:%d, mask:%d  \n", paf->totpart,activeParts,fileParts,readMask);
1782                 elbeemDebugOut(debugStrBuffer);
1783                 return;
1784         } // fluid sim particles done
1785         #endif // DISABLE_ELBEEM
1786         
1787         if(paf->end < paf->sta) return;
1788         
1789         if( (paf->flag & PAF_OFACE) && (paf->flag & PAF_FACE)==0) return;
1790         
1791         if(me->totvert==0) return;
1792         
1793         if(ob==G.obedit) return;
1794         totpart= (G.rendering)?paf->totpart:(paf->disp*paf->totpart)/100;
1795         if(totpart==0) return;
1796         
1797         /* No returns after this line! */
1798         
1799         /* material */
1800         ma= give_current_material(ob, paf->omat);
1801         if(ma) {
1802                 if(paf->speedtex)
1803                         mtexmove= ma->mtex[paf->speedtex-1];
1804                 mtextime= ma->mtex[paf->timetex-1];
1805         }
1806
1807         disable_speed_curve(1); /* check this... */
1808
1809         /* initialize particles */
1810         new_particle(paf);
1811
1812         /* reset deflector cache, sumohandle is free, but its still sorta abuse... (ton) */
1813         for(base= G.scene->base.first; base; base= base->next)
1814                 base->object->sumohandle= NULL;
1815
1816         /* all object positions from start to end */
1817         paf_sta= (int)floor(paf->sta);
1818         paf_end= (int)ceil(paf->end);
1819         if((paf->flag & PAF_STATIC)==0)
1820                 mcache= cache_object_matrices(ob, paf_sta, paf_end);
1821         
1822         /* mult generations? */
1823         for(a=0; a<PAF_MAXMULT; a++) {
1824                 if(paf->mult[a]!=0.0) {
1825                         /* interesting formula! this way after 'x' generations the total is paf->totpart */
1826                         totpart= (int)(totpart / (1.0+paf->mult[a]*paf->child[a]));
1827                 }
1828                 else break;
1829         }
1830
1831         /* for static particles, calculate system on current frame (? ton) */
1832         if(ma) do_mat_ipo(ma);
1833         
1834         /* matrix invert for static too */
1835         Mat4Invert(ob->imat, ob->obmat);
1836         Mat4CpyMat4(paf->imat, ob->imat);       /* used for duplicators */
1837         
1838         /* new random generator */
1839         rng = rng_new(paf->seed);
1840         
1841         /* otherwise it goes way too fast */
1842         force[0]= paf->force[0]*0.05f;
1843         force[1]= paf->force[1]*0.05f;
1844         force[2]= paf->force[2]*0.05f;
1845         
1846         if( paf->flag & PAF_STATIC ) deform= 0;
1847         else {
1848                 deform= (ob->parent && ob->parent->type==OB_LATTICE && ob->partype==PARSKEL);
1849                 if(deform) init_latt_deform(ob->parent, 0);
1850         }
1851         
1852         /* get the effectors */
1853         effectorbase= pdInitEffectors(ob, paf->group);
1854         
1855         /* init geometry, return is 6 x float * me->totvert in size */
1856         vertexcosnos= (VeNoCo *)mesh_get_mapped_verts_nors(ob);
1857         facelist= me->mface;
1858         totvert= me->totvert;
1859         totface= me->totface;
1860         
1861         /* if vertexweights or even distribution, it makes weight tables, also checks where it emits from */
1862         make_weight_tables(paf, me, totpart, vertexcosnos, totvert, facelist, totface, &voweights, &foweights);
1863         
1864         /* vertexweights can define lengths too */
1865         make_length_tables(paf, me, totvert, facelist, totface, &volengths, &folengths);
1866         
1867         /* now define where to emit from, if there are face weights we skip vertices */
1868         if(paf->flag & PAF_OFACE) totvert= 0;
1869         if((paf->flag & PAF_FACE)==0) totface= 0;
1870         if(foweights) totvert= 0;
1871         
1872         /* initialize give_mesh_particle_coord */
1873         if(totface)
1874                 give_mesh_particle_coord(paf, vertexcosnos, facelist, totpart, totface, NULL, NULL);
1875         
1876         /* correction for face timing when using weighted average */
1877         if(totface && foweights) {
1878                 maxw= (paf->end-paf->sta)/foweights[0];
1879         }
1880         else if(totvert && voweights) {
1881                 maxw= (paf->end-paf->sta)/voweights[0];
1882         }
1883         
1884         /* for loop below */
1885         if (paf->flag & PAF_STATIC) {
1886                 ftime = G.scene->r.cfra;
1887                 dtime= 0.0f;
1888         } else {
1889                 ftime= paf->sta;
1890                 dtime= (paf->end - paf->sta)/(float)totpart;
1891         }
1892         
1893         curface= curvert= 0;
1894         for(a=0; a<totpart; a++, ftime+=dtime) {
1895                 
1896                 /* we set waitcursor only when a half second expired, particles now are realtime updated */
1897                 if(waitcursor_set==0 && (a % 256)==255) {
1898                         double seconds= PIL_check_seconds_timer();
1899                         if(seconds - startseconds > 0.5) {
1900                                 waitcursor(1);
1901                                 waitcursor_set= 1;
1902                         }
1903                 }
1904                 
1905                 pa= new_particle(paf);
1906                 pa->time= ftime;
1907
1908                 /* get coordinates from faces, only when vertices set to zero */
1909                 if(totvert==0 && totface) {
1910                         int curjit;
1911                         
1912                         /* use weight table, we have to do faces in order to be able to use jitter table... */
1913                         if(foweights) {
1914                                 
1915                                 if(foweights[curface] < 1.0f) {
1916                                         float remainder= 0.0f;
1917                                         
1918                                         while(remainder + foweights[curface] < 1.0f && curface<totface-1) {
1919                                                 remainder += foweights[curface];
1920                                                 curface++;
1921                                         }
1922                                         /* if this is the last face, the foweights[] can be zero, so we don't add a particle extra */
1923                                         if(curface!=totface-1)
1924                                                 foweights[curface] += remainder;
1925                                         
1926                                         maxw= (paf->end-paf->sta)/foweights[curface];
1927                                 }
1928
1929                                 if(foweights[curface]==0.0f)
1930                                         break;  /* WARN skips here out of particle generating */
1931                                 else {
1932                                         if(foweights[curface] >= 1.0f)          /* note the >= here, this because of the < 1.0f above, it otherwise will stick to 1 face forever */
1933                                                 foweights[curface] -= 1.0f;
1934                                         
1935                                         curjit= (int) foweights[curface];
1936                                         give_mesh_particle_coord(paf, vertexcosnos, facelist+curface, a, curjit, co, no);
1937                                         
1938                                         /* time correction to make particles appear evenly, maxw does interframe (0-1) */
1939                                         pa->time= paf->sta + maxw*foweights[curface];
1940                                 }
1941                         }
1942                         else {
1943                                 curface= a % totface;
1944                                 curjit= a/totface;
1945                                 give_mesh_particle_coord(paf, vertexcosnos, facelist+curface, a, curjit, co, no);
1946                         }
1947                 }
1948                 /* get coordinates from vertices */
1949                 if(totvert) {
1950                         /* use weight table */
1951                         if(voweights) {
1952                                 
1953                                 if(voweights[curvert] < 1.0f) {
1954                                         float remainder= 0.0f;
1955                                         
1956                                         while(remainder + voweights[curvert] < 1.0f && curvert<totvert-1) {
1957                                                 remainder += voweights[curvert];
1958                                                 curvert++;
1959                                         }
1960                                         voweights[curvert] += remainder;
1961                                         maxw= (paf->end-paf->sta)/voweights[curvert];
1962                                 }
1963                                 
1964                                 if(voweights[curvert]==0.0f)
1965                                         break;  /* WARN skips here out of particle generating */
1966                                 else {
1967                                         if(voweights[curvert] > 1.0f)
1968                                                 voweights[curvert] -= 1.0f;
1969                                         
1970                                         /* time correction to make particles appear evenly */
1971                                         pa->time= paf->sta + maxw*voweights[curvert];
1972                                 }
1973                         }
1974                         else {
1975                                 curvert= a % totvert;
1976                                 if(a >= totvert && totface)
1977                                         totvert= 0;
1978                         }
1979                         
1980                         VECCOPY(co, vertexcosnos[curvert].co);
1981                         VECCOPY(no, vertexcosnos[curvert].no);
1982                 }
1983                 
1984                 VECCOPY(pa->co, co);
1985                 
1986                 /* dynamic options */
1987                 if((paf->flag & PAF_STATIC)==0) {
1988                         int cur;
1989                         
1990                         /* particle retiming with texture */
1991                         if(mtextime && (paf->flag2 & PAF_TEXTIME)) {
1992                                 float tin, tr, tg, tb, ta, orco[3];
1993                                 
1994                                 /* calculate normalized orco */
1995                                 orco[0] = (co[0]-me->loc[0])/me->size[0];
1996                                 orco[1] = (co[1]-me->loc[1])/me->size[1];
1997                                 orco[2] = (co[2]-me->loc[2])/me->size[2];
1998                                 externtex(mtextime, orco, &tin, &tr, &tg, &tb, &ta);
1999                                 
2000                                 if(paf->flag2neg & PAF_TEXTIME)
2001                                         pa->time = paf->sta + (paf->end - paf->sta)*tin;
2002                                 else
2003                                         pa->time = paf->sta + (paf->end - paf->sta)*(1.0f-tin);
2004                         }
2005
2006                         /* set ob at correct time, we use cached matrices */
2007                         cur= (int)floor(pa->time) + 1 ;         /* + 1 has a reason: (obmat/prevobmat) otherwise comet-tails start too late */
2008                         
2009                         if(cur <= paf_end) mcnow= mcache + cur - paf_sta;
2010                         else mcnow= mcache + paf_end - paf_sta + 1;
2011                         
2012                         if(cur > paf_sta) mcprev= mcnow-1;
2013                         else mcprev= mcache;
2014                         
2015                         /* move to global space */
2016                         Mat4MulVecfl(mcnow->obmat, pa->co);
2017                 
2018                         VECCOPY(vec, co);
2019                         Mat4MulVecfl(mcprev->obmat, vec);
2020                         
2021                         /* first start speed: object */
2022                         VECSUB(pa->no, pa->co, vec);
2023
2024                         VecMulf(pa->no, paf->obfac);
2025                         
2026                         /* calculate the correct inter-frame */ 
2027                         fac= (pa->time- (float)floor(pa->time));
2028                         pa->co[0]= fac*pa->co[0] + (1.0f-fac)*vec[0];
2029                         pa->co[1]= fac*pa->co[1] + (1.0f-fac)*vec[1];
2030                         pa->co[2]= fac*pa->co[2] + (1.0f-fac)*vec[2];
2031
2032                         /* start speed: normal */
2033                         if(paf->normfac!=0.0) {
2034                                 /* imat is transpose ! */
2035                                 VECCOPY(vec, no);
2036                                 Mat3MulVecfl(mcnow->imat, vec);
2037                         
2038                                 Normalise(vec);
2039                                 VecMulf(vec, paf->normfac);
2040                                 VECADD(pa->no, pa->no, vec);
2041                         }
2042                 }
2043                 else {
2044                         if(paf->normfac!=0.0) {
2045                                 VECCOPY(pa->no, no);
2046                                 Normalise(pa->no);
2047                                 VecMulf(pa->no, paf->normfac);
2048                         }
2049                 }
2050                 
2051                 pa->lifetime= paf->lifetime;
2052                 if(paf->randlife!=0.0) {
2053                         pa->lifetime*= 1.0f + paf->randlife*(rng_getFloat(rng) - 0.5f);
2054                 }
2055                 pa->mat_nr= paf->omat;
2056                 
2057                 if(folengths)
2058                         pa->lifetime*= folengths[curface];
2059
2060                 make_particle_keys(rng, ob, 0, a, paf, pa, force, deform, mtexmove, effectorbase);
2061         }
2062         
2063         /* free stuff */
2064         give_mesh_particle_coord(NULL, NULL, NULL, 0, 0, NULL, NULL);
2065         MEM_freeN(vertexcosnos);
2066         if(voweights) MEM_freeN(voweights);
2067         if(foweights) MEM_freeN(foweights);
2068         if(volengths) MEM_freeN(volengths);
2069         if(folengths) MEM_freeN(folengths);
2070         if(mcache) MEM_freeN(mcache);
2071         rng_free(rng);
2072
2073         if(deform) end_latt_deform();
2074         
2075         if(effectorbase)
2076                 pdEndEffectors(effectorbase);   
2077         
2078         /* reset deflector cache */
2079         for(base= G.scene->base.first; base; base= base->next) {
2080                 if(base->object->sumohandle) {
2081                         MEM_freeN(base->object->sumohandle);
2082                         base->object->sumohandle= NULL;
2083                 }
2084         }
2085
2086         disable_speed_curve(0);
2087         
2088         if(waitcursor_set) waitcursor(0);
2089 }
2090