Added baking for softbodies.
[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 #include "DNA_listBase.h"
40 #include "DNA_effect_types.h"
41 #include "DNA_object_types.h"
42 #include "DNA_object_force.h"
43 #include "DNA_mesh_types.h"
44 #include "DNA_meshdata_types.h"
45 #include "DNA_material_types.h"
46 #include "DNA_curve_types.h"
47 #include "DNA_key_types.h"
48 #include "DNA_texture_types.h"
49 #include "DNA_scene_types.h"
50 #include "DNA_lattice_types.h"
51 #include "DNA_ipo_types.h"
52
53 #include "BLI_blenlib.h"
54 #include "BLI_arithb.h"
55 #include "BLI_rand.h"
56
57 #include "BKE_action.h"
58 #include "BKE_bad_level_calls.h"
59 #include "BKE_blender.h"
60 #include "BKE_constraint.h"
61 #include "BKE_deform.h"
62 #include "BKE_displist.h"
63 #include "BKE_DerivedMesh.h"
64 #include "BKE_effect.h"
65 #include "BKE_global.h"
66 #include "BKE_ipo.h"
67 #include "BKE_key.h"
68 #include "BKE_lattice.h"
69 #include "BKE_mesh.h"
70 #include "BKE_material.h"
71 #include "BKE_main.h"
72 #include "BKE_object.h"
73 #include "BKE_screen.h"
74 #include "BKE_utildefines.h"
75
76 #include "render.h"  // externtex, bad level call (ton)
77
78
79 #ifdef HAVE_CONFIG_H
80 #include <config.h>
81 #endif
82
83 Effect *add_effect(int type)
84 {
85         Effect *eff=0;
86         BuildEff *bld;
87         PartEff *paf;
88         WaveEff *wav;
89         int a;
90         
91         switch(type) {
92         case EFF_BUILD:
93                 bld= MEM_callocN(sizeof(BuildEff), "neweff");
94                 eff= (Effect *)bld;
95                 
96                 bld->sfra= 1.0;
97                 bld->len= 100.0;
98                 break;
99                 
100         case EFF_PARTICLE:
101                 paf= MEM_callocN(sizeof(PartEff), "neweff");
102                 eff= (Effect *)paf;
103                 
104                 paf->sta= 1.0;
105                 paf->end= 100.0;
106                 paf->lifetime= 50.0;
107                 for(a=0; a<PAF_MAXMULT; a++) {
108                         paf->life[a]= 50.0;
109                         paf->child[a]= 4;
110                         paf->mat[a]= 1;
111                 }
112                 
113                 paf->totpart= 1000;
114                 paf->totkey= 8;
115                 paf->staticstep= 5;
116                 paf->defvec[2]= 1.0f;
117                 paf->nabla= 0.05f;
118                 
119                 break;
120                 
121         case EFF_WAVE:
122                 wav= MEM_callocN(sizeof(WaveEff), "neweff");
123                 eff= (Effect *)wav;
124                 
125                 wav->flag |= (WAV_X+WAV_Y+WAV_CYCL);
126                 
127                 wav->height= 0.5f;
128                 wav->width= 1.5f;
129                 wav->speed= 0.5f;
130                 wav->narrow= 1.5f;
131                 wav->lifetime= 0.0f;
132                 wav->damp= 10.0f;
133                 
134                 break;
135         }
136         
137         eff->type= eff->buttype= type;
138         eff->flag |= SELECT;
139         
140         return eff;
141 }
142
143 void free_effect(Effect *eff)
144 {
145         PartEff *paf;
146         
147         if(eff->type==EFF_PARTICLE) {
148                 paf= (PartEff *)eff;
149                 if(paf->keys) MEM_freeN(paf->keys);
150         }
151         MEM_freeN(eff); 
152 }
153
154
155 void free_effects(ListBase *lb)
156 {
157         Effect *eff;
158         
159         eff= lb->first;
160         while(eff) {
161                 BLI_remlink(lb, eff);
162                 free_effect(eff);
163                 eff= lb->first;
164         }
165 }
166
167 Effect *copy_effect(Effect *eff) 
168 {
169         Effect *effn;
170
171         effn= MEM_dupallocN(eff);
172         if(effn->type==EFF_PARTICLE) ((PartEff *)effn)->keys= 0;
173
174         return effn;    
175 }
176
177 void copy_act_effect(Object *ob)
178 {
179         /* return a copy of the active effect */
180         Effect *effn, *eff;
181         
182         eff= ob->effect.first;
183         while(eff) {
184                 if(eff->flag & SELECT) {
185                         
186                         effn= copy_effect(eff);
187                         BLI_addtail(&ob->effect, effn);
188                         
189                         eff->flag &= ~SELECT;
190                         return;
191                         
192                 }
193                 eff= eff->next;
194         }
195         
196         /* when it comes here: add new effect */
197         eff= add_effect(EFF_BUILD);
198         BLI_addtail(&ob->effect, eff);
199                         
200 }
201
202 void copy_effects(ListBase *lbn, ListBase *lb)
203 {
204         Effect *eff, *effn;
205
206         lbn->first= lbn->last= 0;
207
208         eff= lb->first;
209         while(eff) {
210                 effn= copy_effect(eff);
211                 BLI_addtail(lbn, effn);
212                 
213                 eff= eff->next;
214         }
215         
216 }
217
218 void deselectall_eff(Object *ob)
219 {
220         Effect *eff= ob->effect.first;
221         
222         while(eff) {
223                 eff->flag &= ~SELECT;
224                 eff= eff->next;
225         }
226 }
227
228 void set_buildvars(Object *ob, int *start, int *end)
229 {
230         BuildEff *bld;
231         float ctime;
232         
233         bld= ob->effect.first;
234         while(bld) {
235                 if(bld->type==EFF_BUILD) {
236                         ctime= bsystem_time(ob, 0, (float)G.scene->r.cfra, bld->sfra-1.0f);
237                         if(ctime < 0.0) {
238                                 *end= *start;
239                         }
240                         else if(ctime < bld->len) {
241                                 *end= *start+ (int)((*end - *start)*ctime/bld->len);
242                         }
243                         
244                         return;
245                 }
246                 bld= bld->next;
247         }
248 }
249
250 /* ***************** PARTICLES ***************** */
251
252 Particle *new_particle(PartEff *paf)
253 {
254         static Particle *pa;
255         static int cur;
256
257         /* we agree: when paf->keys==0: alloc */        
258         if(paf->keys==0) {
259                 pa= paf->keys= MEM_callocN( paf->totkey*paf->totpart*sizeof(Particle), "particlekeys" );
260                 cur= 0;
261         }
262         else {
263                 if(cur && cur<paf->totpart) pa+=paf->totkey;
264                 cur++;
265         }
266         return pa;
267 }
268
269 PartEff *give_parteff(Object *ob)
270 {
271         PartEff *paf;
272         
273         paf= ob->effect.first;
274         while(paf) {
275                 if(paf->type==EFF_PARTICLE) return paf;
276                 paf= paf->next;
277         }
278         return 0;
279 }
280
281 void where_is_particle(PartEff *paf, Particle *pa, float ctime, float *vec)
282 {
283         Particle *p[4];
284         float dt, t[4];
285         int a;
286         
287         if(paf->totkey==1) {
288                 VECCOPY(vec, pa->co);
289                 return;
290         }
291         
292         /* first find the first particlekey */
293         a= (int)((paf->totkey-1)*(ctime-pa->time)/pa->lifetime);
294         if(a>=paf->totkey) a= paf->totkey-1;
295         else if(a<0) a= 0;
296         
297         pa+= a;
298         
299         if(a>0) p[0]= pa-1; else p[0]= pa;
300         p[1]= pa;
301         
302         if(a+1<paf->totkey) p[2]= pa+1; else p[2]= pa;
303         if(a+2<paf->totkey) p[3]= pa+2; else p[3]= p[2];
304         
305         if(p[1]==p[2]) dt= 0.0;
306         else dt= (ctime-p[1]->time)/(p[2]->time - p[1]->time);
307
308         if(paf->flag & PAF_BSPLINE) set_four_ipo(dt, t, KEY_BSPLINE);
309         else set_four_ipo(dt, t, KEY_CARDINAL);
310
311         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];
312         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];
313         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];
314
315 }
316
317 void particle_tex(MTex *mtex, PartEff *paf, float *co, float *no)
318 {                               
319         float tin, tr, tg, tb, ta;
320         float old;
321         
322         externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
323
324         if(paf->texmap==PAF_TEXINT) {
325                 tin*= paf->texfac;
326                 no[0]+= tin*paf->defvec[0];
327                 no[1]+= tin*paf->defvec[1];
328                 no[2]+= tin*paf->defvec[2];
329         }
330         else if(paf->texmap==PAF_TEXRGB) {
331                 no[0]+= (tr-0.5f)*paf->texfac;
332                 no[1]+= (tg-0.5f)*paf->texfac;
333                 no[2]+= (tb-0.5f)*paf->texfac;
334         }
335         else {  /* PAF_TEXGRAD */
336                 
337                 old= tin;
338                 co[0]+= paf->nabla;
339                 externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
340                 no[0]+= (old-tin)*paf->texfac;
341                 
342                 co[0]-= paf->nabla;
343                 co[1]+= paf->nabla;
344                 externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
345                 no[1]+= (old-tin)*paf->texfac;
346                 
347                 co[1]-= paf->nabla;
348                 co[2]+= paf->nabla;
349                 externtex(mtex, co, &tin, &tr, &tg, &tb, &ta);
350                 no[2]+= (old-tin)*paf->texfac;
351                 
352         }
353 }
354
355 static int linetriangle(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *labda)
356 {
357         float p[3], s[3], d[3], e1[3], e2[3], q[3];
358         float a, f, u, v;
359         
360         VECSUB(e1, v1, v0);
361         VECSUB(e2, v2, v0);
362         VECSUB(d, p2, p1);
363         
364         Crossf(p, d, e2);
365         a = INPR(e1, p);
366         if ((a > -0.000001) && (a < 0.000001)) return 0;
367         f = 1.0f/a;
368         
369         VECSUB(s, p1, v0);
370         
371         Crossf(q, s, e1);
372         *labda = f * INPR(e2, q);
373         if ((*labda < 0.0)||(*labda > 1.0)) return 0;
374         
375         u = f * INPR(s, p);
376         if ((u < 0.0)||(u > 1.0)) return 0;
377         
378         v = f * INPR(d, q);
379         if ((v < 0.0)||((u + v) > 1.0)) return 0;
380         
381         return 1;
382 }
383
384 /*  -------- pdDoEffector() --------
385     generic force/speed system, now used for particles and softbodies
386         opco            = global coord, as input
387     force               = force accumulator
388     speed               = speed accumulator
389     cur_time    = in frames
390     par_layer   = layer the caller is in
391
392 */
393 void pdDoEffector(float *opco, float *force, float *speed, float cur_time, unsigned int par_layer,unsigned int flags)
394 {
395 /*
396         Modifies the force on a particle according to its
397         relation with the effector object
398         Different kind of effectors include:
399                 Forcefields: Gravity-like attractor
400                 (force power is related to the inverse of distance to the power of a falloff value)
401                 Vortex fields: swirling effectors
402                 (particles rotate around Z-axis of the object. otherwise, same relation as)
403                 (Forcefields, but this is not done through a force/acceleration)
404         
405 */
406         Object *ob;
407         Base *base;
408         PartDeflect *pd;
409         float vect_to_vert[3];
410         float force_vec[3];
411         float f_force, distance;
412         float *obloc;
413         float force_val, ffall_val;
414         short cur_frame;
415
416         /* Cycle through objects, get total of (1/(gravity_strength * dist^gravity_power)) */
417         /* Check for min distance here? (yes would be cool to add that, ton) */
418         
419         for(base = G.scene->base.first; base; base= base->next) {
420                 if( (base->lay & par_layer) && base->object->pd) {
421                         ob= base->object;
422                         pd= ob->pd;
423                         
424                         /* checking if to continue or not */
425                         if(pd->forcefield==0) continue;
426                         
427                         /* Get IPO force strength and fall off values here */
428                         if (has_ipo_code(ob->ipo, OB_PD_FSTR))
429                                 force_val = IPO_GetFloatValue(ob->ipo, OB_PD_FSTR, cur_time);
430                         else 
431                                 force_val = pd->f_strength;
432                         
433                         if (has_ipo_code(ob->ipo, OB_PD_FFALL)) 
434                                 ffall_val = IPO_GetFloatValue(ob->ipo, OB_PD_FFALL, cur_time);
435                         else 
436                                 ffall_val = pd->f_power;
437                         
438                         
439                         /* Need to set r.cfra for paths (investigate, ton) (uses ob->ctime now, ton) */
440                         if(ob->ctime!=cur_time) {
441                                 cur_frame = G.scene->r.cfra;
442                                 G.scene->r.cfra = (short)cur_time;
443                                 where_is_object_time(ob, cur_time);
444                                 G.scene->r.cfra = cur_frame;
445                         }
446                         
447                         /* use center of object for distance calculus */
448                         obloc= ob->obmat[3];
449                         VECSUB(vect_to_vert, obloc, opco);
450                         distance = VecLength(vect_to_vert);
451                         
452                         if((pd->flag & PFIELD_USEMAX) && distance>pd->maxdist)
453                                 ;       /* don't do anything */
454                         else if(pd->forcefield == PFIELD_WIND) {
455                                 VECCOPY(force_vec, ob->obmat[2]);
456                                 
457                                 /* wind works harder perpendicular to normal, would be nice for softbody later (ton) */
458                                 
459                                 /* Limit minimum distance to vertex so that */
460                                 /* the force is not too big */
461                                 if (distance < 0.001) distance = 0.001f;
462                                 f_force = (force_val)*(1/(1000 * (float)pow((double)distance, (double)ffall_val)));
463                                 if(flags &&PE_WIND_AS_SPEED){
464                                 speed[0] -= (force_vec[0] * f_force );
465                                 speed[1] -= (force_vec[1] * f_force );
466                                 speed[2] -= (force_vec[2] * f_force );
467                                 }
468                                 else{
469                                 force[0] += force_vec[0]*f_force;
470                                 force[1] += force_vec[1]*f_force;
471                                 force[2] += force_vec[2]*f_force;
472                                 }
473                         }
474                         else if(pd->forcefield == PFIELD_FORCE) {
475                                 
476                                 /* only use center of object */
477                                 obloc= ob->obmat[3];
478
479                                 /* Now calculate the gravitational force */
480                                 VECSUB(vect_to_vert, obloc, opco);
481                                 distance = VecLength(vect_to_vert);
482
483                                 /* Limit minimum distance to vertex so that */
484                                 /* the force is not too big */
485                                 if (distance < 0.001) distance = 0.001f;
486                                 f_force = (force_val)*(1/(1000 * (float)pow((double)distance, (double)ffall_val)));
487                                 force[0] += (vect_to_vert[0] * f_force );
488                                 force[1] += (vect_to_vert[1] * f_force );
489                                 force[2] += (vect_to_vert[2] * f_force );
490
491                         }
492                         else if(pd->forcefield == PFIELD_VORTEX) {
493
494                                 /* only use center of object */
495                                 obloc= ob->obmat[3];
496
497                                 /* Now calculate the vortex force */
498                                 VECSUB(vect_to_vert, obloc, opco);
499                                 distance = VecLength(vect_to_vert);
500
501                                 Crossf(force_vec, ob->obmat[2], vect_to_vert);
502                                 Normalise(force_vec);
503
504                                 /* Limit minimum distance to vertex so that */
505                                 /* the force is not too big */
506                                 if (distance < 0.001) distance = 0.001f;
507                                 f_force = (force_val)*(1/(100 * (float)pow((double)distance, (double)ffall_val)));
508                                 speed[0] -= (force_vec[0] * f_force );
509                                 speed[1] -= (force_vec[1] * f_force );
510                                 speed[2] -= (force_vec[2] * f_force );
511
512                         }
513                 }
514         }
515 }
516
517 static void cache_object_vertices(Object *ob)
518 {
519         Mesh *me;
520         MVert *mvert;
521         float *fp;
522         int a;
523         
524         me= ob->data;
525         if(me->totvert==0) return;
526
527         fp= ob->sumohandle= MEM_mallocN(3*sizeof(float)*me->totvert, "cache particles");
528         mvert= me->mvert;
529         a= me->totvert;
530         while(a--) {
531                 VECCOPY(fp, mvert->co);
532                 Mat4MulVecfl(ob->obmat, fp);
533                 mvert++;
534                 fp+= 3;
535         }
536 }
537
538 int pdDoDeflection(float opco[3], float npco[3], float opno[3],
539         float npno[3], float life, float force[3], int def_depth,
540         float cur_time, unsigned int par_layer, int *last_object,
541                 int *last_face, int *same_face)
542 {
543         /* Particle deflection code */
544         /* The code is in two sections: the first part checks whether a particle has            */
545         /* intersected a face of a deflector mesh, given its old and new co-ords, opco and npco */
546         /* and which face it hit first                                                          */
547         /* The second part calculates the new co-ordinates given that collision and updates     */
548         /* the new co-ordinates accordingly */
549         Base *base;
550         Object *ob, *deflection_object = NULL;
551         Mesh *def_mesh;
552         MFace *mface, *deflection_face = NULL;
553         float *v1, *v2, *v3, *v4, *vcache=NULL;
554         float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3];
555         float dv1[3], dv2[3], dv3[3];
556         float vect_to_int[3], refl_vel[3];
557         float d_intersect_co[3], d_intersect_vect[3], d_nvect[3], d_i_co_above[3];
558         float forcec[3];
559         float k_point3, dist_to_plane;
560         float first_dist, ref_plane_mag;
561         float dk_plane=0, dk_point1=0;
562         float icalctop, icalcbot, n_mag;
563         float mag_iv, x_m,y_m,z_m;
564         float damping, perm_thresh;
565         float perm_val, rdamp_val;
566         int a, deflected=0, deflected_now=0;
567         float t,t2, min_t;
568         float mat[3][3], obloc[3];
569         short cur_frame;
570         float time_before, time_after;
571         float force_mag_norm;
572         int d_object=0, d_face=0, ds_object=0, ds_face=0;
573
574         first_dist = 200000;
575         min_t = 200000;
576
577         /* The first part of the code, finding the first intersected face*/
578         base= G.scene->base.first;
579         while (base) {
580                 /*Only proceed for mesh object in same layer */
581                 if(base->object->type==OB_MESH && (base->lay & par_layer)) {
582                         ob= base->object;
583                         /* only with deflecting set */
584                         if(ob->pd && ob->pd->deflect) {
585                                 def_mesh= ob->data;
586                         
587                                 d_object = d_object + 1;
588
589                                 d_face = d_face + 1;
590                                 mface= def_mesh->mface;
591                                 a = def_mesh->totface;
592                                 
593                                 
594                                 if(ob->parent==NULL && ob->ipo==NULL) { // static
595                                         if(ob->sumohandle==NULL) cache_object_vertices(ob);
596                                         vcache= ob->sumohandle;
597                                 }
598                                 else {
599                                         /*Find out where the object is at this time*/
600                                         cur_frame = G.scene->r.cfra;
601                                         G.scene->r.cfra = (short)cur_time;
602                                         where_is_object_time(ob, cur_time);
603                                         G.scene->r.cfra = cur_frame;
604                                         
605                                         /*Pass the values from ob->obmat to mat*/
606                                         /*and the location values to obloc           */
607                                         Mat3CpyMat4(mat,ob->obmat);
608                                         obloc[0] = ob->obmat[3][0];
609                                         obloc[1] = ob->obmat[3][1];
610                                         obloc[2] = ob->obmat[3][2];
611                                         vcache= NULL;
612
613                                 }
614                                 
615                                 while (a--) {
616
617                                         if(vcache) {
618                                                 v1= vcache+ 3*(mface->v1);
619                                                 VECCOPY(nv1, v1);
620                                                 v1= vcache+ 3*(mface->v2);
621                                                 VECCOPY(nv2, v1);
622                                                 v1= vcache+ 3*(mface->v3);
623                                                 VECCOPY(nv3, v1);
624                                                 v1= vcache+ 3*(mface->v4);
625                                                 VECCOPY(nv4, v1);
626                                         }
627                                         else {
628                                                 /* Calculate the global co-ordinates of the vertices*/
629                                                 v1= (def_mesh->mvert+(mface->v1))->co;
630                                                 v2= (def_mesh->mvert+(mface->v2))->co;
631                                                 v3= (def_mesh->mvert+(mface->v3))->co;
632                                                 v4= (def_mesh->mvert+(mface->v4))->co;
633         
634                                                 VECCOPY(nv1, v1);
635                                                 VECCOPY(nv2, v2);
636                                                 VECCOPY(nv3, v3);
637                                                 VECCOPY(nv4, v4);
638         
639                                                 /*Apply the objects deformation matrix*/
640                                                 Mat3MulVecfl(mat, nv1);
641                                                 Mat3MulVecfl(mat, nv2);
642                                                 Mat3MulVecfl(mat, nv3);
643                                                 Mat3MulVecfl(mat, nv4);
644         
645                                                 VECADD(nv1, nv1, obloc);
646                                                 VECADD(nv2, nv2, obloc);
647                                                 VECADD(nv3, nv3, obloc);
648                                                 VECADD(nv4, nv4, obloc);
649                                         }
650                                         
651                                         deflected_now = 0;
652
653                                                 
654                                                 
655 //                                      t= 0.5; // this is labda of line, can use it optimize quad intersection
656 // sorry but no .. see below (BM)                                       
657                                         if( linetriangle(opco, npco, nv1, nv2, nv3, &t) ) {
658                                                 if (t < min_t) {
659                                                         deflected = 1;
660                                                         deflected_now = 1;
661                                                 }
662                                         }
663 //                                      else if (mface->v4 && (t>=0.0 && t<=1.0)) {
664 // no, you can't skip testing the other triangle
665 // it might give a smaller t on (close to) the edge .. this is numerics not esoteric maths :)
666 // note: the 2 triangles don't need to share a plane ! (BM)
667                                         if (mface->v4) {
668                                                 if( linetriangle(opco, npco, nv1, nv3, nv4, &t2) ) {
669                                                         if (t2 < min_t) {
670                                                                 deflected = 1;
671                                                                 deflected_now = 2;
672                                                         }
673                                                 }
674                                         }
675                                         
676                                         if ((deflected_now > 0) && ((t < min_t) ||(t2 < min_t))) {
677                         min_t = t;
678                         ds_object = d_object;
679                                                 ds_face = d_face;
680                                                 deflection_object = ob;
681                                                 deflection_face = mface;
682                                                 if (deflected_now==1) {
683                         min_t = t;
684                                                         VECCOPY(dv1, nv1);
685                                                         VECCOPY(dv2, nv2);
686                                                         VECCOPY(dv3, nv3);
687                                                 }
688                                                 else {
689                         min_t = t2;
690                                                         VECCOPY(dv1, nv1);
691                                                         VECCOPY(dv2, nv3);
692                                                         VECCOPY(dv3, nv4);
693                                                 }
694                                         }
695                                         mface++;
696                                 }
697                         }
698                 }
699                 base = base->next;
700         }
701
702
703         /* Here's the point to do the permeability calculation */
704         /* Set deflected to 0 if a random number is below the value */
705         /* Get the permeability IPO here*/
706         if (deflected) {
707                 
708                 if (has_ipo_code(deflection_object->ipo, OB_PD_PERM)) 
709                         perm_val = IPO_GetFloatValue(deflection_object->ipo, OB_PD_PERM, cur_time);
710                 else 
711                         perm_val = deflection_object->pd->pdef_perm;
712
713                 perm_thresh =  (float)BLI_drand() - perm_val;
714                 if (perm_thresh < 0 ) {
715                         deflected = 0;
716                 }
717         }
718
719         /* Now for the second part of the deflection code - work out the new speed */
720         /* and position of the particle if a collision occurred */
721         if (deflected) {
722         VECSUB(edge1, dv1, dv2);
723                 VECSUB(edge2, dv3, dv2);
724                 Crossf(d_nvect, edge2, edge1);
725                 n_mag = Normalise(d_nvect);
726                 dk_plane = INPR(d_nvect, nv1);
727                 dk_point1 = INPR(d_nvect,opco);
728
729                 VECSUB(d_intersect_vect, npco, opco);
730
731                 d_intersect_co[0] = opco[0] + (min_t * (npco[0] - opco[0]));
732                 d_intersect_co[1] = opco[1] + (min_t * (npco[1] - opco[1]));
733                 d_intersect_co[2] = opco[2] + (min_t * (npco[2] - opco[2]));
734                 
735                 d_i_co_above[0] = (d_intersect_co[0] + (0.001f * d_nvect[0]));
736                 d_i_co_above[1] = (d_intersect_co[1] + (0.001f * d_nvect[1]));
737                 d_i_co_above[2] = (d_intersect_co[2] + (0.001f * d_nvect[2]));
738                 mag_iv = Normalise(d_intersect_vect);
739                 VECCOPY(npco, d_intersect_co);
740                 
741                 VECSUB(vect_to_int, opco, d_intersect_co);
742                 first_dist = Normalise(vect_to_int);
743
744                 /* Work out the lengths of time before and after collision*/
745                 time_before = (life*(first_dist / (mag_iv)));
746                 time_after =  (life*((mag_iv - first_dist) / (mag_iv)));
747
748                 /* We have to recalculate what the speed would have been at the */
749                 /* point of collision, not the key frame time */
750                 npno[0]= opno[0] + time_before*force[0];
751                 npno[1]= opno[1] + time_before*force[1];
752                 npno[2]= opno[2] + time_before*force[2];
753
754
755                 /* Reflect the speed vector in the face */
756                 x_m = (2 * npno[0] * d_nvect[0]);
757                 y_m = (2 * npno[1] * d_nvect[1]);
758                 z_m = (2 * npno[2] * d_nvect[2]);
759                 refl_vel[0] = npno[0] - (d_nvect[0] * (x_m + y_m + z_m));
760                 refl_vel[1] = npno[1] - (d_nvect[1] * (x_m + y_m + z_m));
761                 refl_vel[2] = npno[2] - (d_nvect[2] * (x_m + y_m + z_m));
762
763                 /*A random variation in the damping factor........ */
764                 /*Get the IPO values for damping here*/
765                 
766                 if (has_ipo_code(deflection_object->ipo, OB_PD_SDAMP)) 
767                         damping = IPO_GetFloatValue(deflection_object->ipo, OB_PD_SDAMP, cur_time);
768                 else 
769                         damping = deflection_object->pd->pdef_damp;
770                 
771                 if (has_ipo_code(deflection_object->ipo, OB_PD_RDAMP)) 
772                         rdamp_val = IPO_GetFloatValue(deflection_object->ipo, OB_PD_RDAMP, cur_time);
773                 else 
774                         rdamp_val = deflection_object->pd->pdef_rdamp;
775
776                 damping = damping + ((1 - damping) * ((float)BLI_drand()*rdamp_val));
777                 damping = damping * damping;
778         ref_plane_mag = INPR(refl_vel,d_nvect);
779
780                 if (damping > 0.999) damping = 0.999f;
781
782                 /* Now add in the damping force - only damp in the direction of */
783                 /* the faces normal vector */
784                 npno[0] = (refl_vel[0] - (d_nvect[0] * ref_plane_mag * damping));
785                 npno[1] = (refl_vel[1] - (d_nvect[1] * ref_plane_mag * damping));
786                 npno[2] = (refl_vel[2] - (d_nvect[2] * ref_plane_mag * damping));
787
788                 /* Now reset opno */
789                 VECCOPY(opno,npno);
790                 VECCOPY(forcec, force);
791
792                 /* If the particle has bounced more than four times on the same */
793                 /* face within this cycle (depth > 4, same face > 4 )           */
794                 /* Then set the force to be only that component of the force    */
795                 /* in the same direction as the face normal                     */
796                 /* i.e. subtract the component of the force in the direction    */
797                 /* of the face normal from the actual force                     */
798                 if ((ds_object == *last_object) && (ds_face == *last_face)) {
799                         /* Increment same_face */
800                         *same_face = *same_face + 1;
801                         if ((*same_face > 3) && (def_depth > 3)) {
802                 force_mag_norm = INPR(forcec, d_nvect);
803                 forcec[0] = forcec[0] - (d_nvect[0] * force_mag_norm);
804                 forcec[1] = forcec[1] - (d_nvect[1] * force_mag_norm);
805                 forcec[2] = forcec[2] - (d_nvect[2] * force_mag_norm);
806                         }
807                 }
808                 else *same_face = 1;
809
810                 *last_object = ds_object;
811                 *last_face = ds_face;
812
813                 /* We have the particles speed at the point of collision    */
814                 /* Now we want the particles speed at the current key frame */
815
816                 npno[0]= npno[0] + time_after*forcec[0];
817                 npno[1]= npno[1] + time_after*forcec[1];
818                 npno[2]= npno[2] + time_after*forcec[2];
819
820                 /* Now we have to recalculate pa->co for the remainder*/
821                 /* of the time since the intersect*/
822                 npco[0]= npco[0] + time_after*npno[0];
823                 npco[1]= npco[1] + time_after*npno[1];
824                 npco[2]= npco[2] + time_after*npno[2];
825
826                 /* And set the old co-ordinates back to the point just above the intersection */
827                 VECCOPY(opco, d_i_co_above);
828
829                 /* Finally update the time */
830                 life = time_after;
831                 cur_time += time_before;
832
833                 /* The particle may have fallen through the face again by now!!*/
834                 /* So check if the particle has changed sides of the plane compared*/
835                 /* the co-ordinates at the last keyframe*/
836                 /* But only do this as a last resort, if we've got to the end of the */
837                 /* number of collisions allowed */
838                 if (def_depth==9) {
839                         k_point3 = INPR(d_nvect,npco);
840                         if (((dk_plane > k_point3) && (dk_plane < dk_point1))||((dk_plane < k_point3) && (dk_plane > dk_point1))) {
841
842                                 /* Yup, the pesky particle may have fallen through a hole!!! */
843                 /* So we'll cheat a bit and move the particle along the normal vector */
844                 /* until it's just the other side of the plane */
845                 icalctop = (dk_plane - d_nvect[0]*npco[0] - d_nvect[1]*npco[1] - d_nvect[2]*npco[2]);
846                 icalcbot = (d_nvect[0]*d_nvect[0] + d_nvect[1]*d_nvect[1] + d_nvect[2]*d_nvect[2]);
847                 dist_to_plane = icalctop / icalcbot;
848
849                 /*  Now just increase the distance a little to place */
850                 /* the point the other side of the plane */
851                 dist_to_plane *= 1.1f;
852                 npco[0]= npco[0] + (dist_to_plane * d_nvect[0]);
853                 npco[1]= npco[1] + (dist_to_plane * d_nvect[1]);
854                 npco[2]= npco[2] + (dist_to_plane * d_nvect[2]);
855
856                         }
857                 }
858         }
859         return deflected;
860 }
861
862 void make_particle_keys(int depth, int nr, PartEff *paf, Particle *part, float *force, int deform, MTex *mtex, unsigned int par_layer)
863 {
864         Particle *pa, *opa = NULL;
865         float damp, deltalife, life;
866         float cur_time;
867         float opco[3], opno[3], npco[3], npno[3], new_force[3], new_speed[3];
868         int b, rt1, rt2, deflected, deflection, finish_defs, def_count;
869         int last_ob, last_fc, same_fc;
870
871         damp= 1.0f-paf->damp;
872         pa= part;
873
874         /* start speed: random */
875         if(paf->randfac!=0.0) {
876                 pa->no[0]+= (float)(paf->randfac*( BLI_drand() -0.5));
877                 pa->no[1]+= (float)(paf->randfac*( BLI_drand() -0.5));
878                 pa->no[2]+= (float)(paf->randfac*( BLI_drand() -0.5));
879         }
880
881         /* start speed: texture */
882         if(mtex && paf->texfac!=0.0) {
883                 particle_tex(mtex, paf, pa->co, pa->no);
884         }
885
886         if(paf->totkey>1) deltalife= pa->lifetime/(paf->totkey-1);
887         else deltalife= pa->lifetime;
888
889         opa= pa;
890         pa++;
891
892         b= paf->totkey-1;
893         while(b--) {
894                 /* new time */
895                 pa->time= opa->time+deltalife;
896
897                 /* set initial variables                                */
898                 opco[0] = opa->co[0];
899                 opco[1] = opa->co[1];
900                 opco[2] = opa->co[2];
901
902                 new_force[0] = force[0];
903                 new_force[1] = force[1];
904                 new_force[2] = force[2];
905                 new_speed[0] = 0.0;
906                 new_speed[1] = 0.0;
907                 new_speed[2] = 0.0;
908
909                 /* Check force field */
910                 cur_time = pa->time;
911                 pdDoEffector(opco, new_force, new_speed, cur_time, par_layer,0);
912
913                 /* new location */
914                 pa->co[0]= opa->co[0] + deltalife * (opa->no[0] + new_speed[0] + 0.5f*new_force[0]);
915                 pa->co[1]= opa->co[1] + deltalife * (opa->no[1] + new_speed[1] + 0.5f*new_force[1]);
916                 pa->co[2]= opa->co[2] + deltalife * (opa->no[2] + new_speed[2] + 0.5f*new_force[2]);
917
918                 /* new speed */
919                 pa->no[0]= opa->no[0] + deltalife*new_force[0];
920                 pa->no[1]= opa->no[1] + deltalife*new_force[1];
921                 pa->no[2]= opa->no[2] + deltalife*new_force[2];
922
923                 /* Particle deflection code                             */
924                 deflection = 0;
925                 finish_defs = 1;
926                 def_count = 0;
927
928                 VECCOPY(opno, opa->no);
929                 VECCOPY(npco, pa->co);
930                 VECCOPY(npno, pa->no);
931
932                 life = deltalife;
933                 cur_time -= deltalife;
934
935                 last_ob = -1;
936                 last_fc = -1;
937                 same_fc = 0;
938
939                 /* First call the particle deflection check for the particle moving   */
940                 /* between the old co-ordinates and the new co-ordinates              */
941                 /* If a deflection occurs, call the code again, this time between the */
942                 /* intersection point and the updated new co-ordinates                */
943                 /* Bail out if we've done the calculation 10 times - this seems ok     */
944         /* for most scenes I've tested */
945                 while (finish_defs) {
946                         deflected =  pdDoDeflection(opco, npco, opno, npno, life, new_force,
947                                                         def_count, cur_time, par_layer,
948                                                         &last_ob, &last_fc, &same_fc);
949                         if (deflected) {
950                                 def_count = def_count + 1;
951                                 deflection = 1;
952                                 if (def_count==10) finish_defs = 0;
953                         }
954                         else {
955                                 finish_defs = 0;
956                         }
957                 }
958
959                 /* Only update the particle positions and speed if we had a deflection */
960                 if (deflection) {
961                         pa->co[0] = npco[0];
962                         pa->co[1] = npco[1];
963                         pa->co[2] = npco[2];
964                         pa->no[0] = npno[0];
965                         pa->no[1] = npno[1];
966                         pa->no[2] = npno[2];
967                 }
968
969
970                 /* speed: texture */
971                 if(mtex && paf->texfac!=0.0) {
972                         particle_tex(mtex, paf, pa->co, pa->no);
973                 }
974                 if(damp!=1.0) {
975                         pa->no[0]*= damp;
976                         pa->no[1]*= damp;
977                         pa->no[2]*= damp;
978                 }
979         
980
981
982                 opa= pa;
983                 pa++;
984                 /* opa is used later on too! */
985         }
986
987         if(deform) {
988                 /* deform all keys */
989                 pa= part;
990                 b= paf->totkey;
991                 while(b--) {
992                         calc_latt_deform(pa->co);
993                         pa++;
994                 }
995         }
996         
997         /* the big multiplication */
998         if(depth<PAF_MAXMULT && paf->mult[depth]!=0.0) {
999                 
1000                 /* new 'child' emerges from an average 'mult' part from 
1001                         the particles */
1002                 damp = (float)nr;
1003                 rt1= (int)(damp*paf->mult[depth]);
1004                 rt2= (int)((damp+1.0)*paf->mult[depth]);
1005                 if(rt1!=rt2) {
1006                         
1007                         for(b=0; b<paf->child[depth]; b++) {
1008                                 pa= new_particle(paf);
1009                                 *pa= *opa;
1010                                 pa->lifetime= paf->life[depth];
1011                                 if(paf->randlife!=0.0) {
1012                                         pa->lifetime*= 1.0f+ (float)(paf->randlife*( BLI_drand() - 0.5));
1013                                 }
1014                                 pa->mat_nr= paf->mat[depth];
1015
1016                                 make_particle_keys(depth+1, b, paf, pa, force, deform, mtex, par_layer);
1017                         }
1018                 }
1019         }
1020 }
1021
1022 void init_mv_jit(float *jit, int num,int seed2)
1023 {
1024         float *jit2, x, rad1, rad2, rad3;
1025         int i, num2;
1026
1027         if(num==0) return;
1028
1029         rad1= (float)(1.0/sqrt((float)num));
1030         rad2= (float)(1.0/((float)num));
1031         rad3= (float)sqrt((float)num)/((float)num);
1032
1033         BLI_srand(31415926 + num + seed2);
1034         x= 0;
1035         num2 = 2 * num;
1036         for(i=0; i<num2; i+=2) {
1037         
1038                 jit[i]= x+ (float)(rad1*(0.5-BLI_drand()));
1039                 jit[i+1]= ((float)i/2)/num +(float)(rad1*(0.5-BLI_drand()));
1040                 
1041                 jit[i]-= (float)floor(jit[i]);
1042                 jit[i+1]-= (float)floor(jit[i+1]);
1043                 
1044                 x+= rad3;
1045                 x -= (float)floor(x);
1046         }
1047
1048         jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
1049
1050         for (i=0 ; i<4 ; i++) {
1051                 RE_jitterate1(jit, jit2, num, rad1);
1052                 RE_jitterate1(jit, jit2, num, rad1);
1053                 RE_jitterate2(jit, jit2, num, rad2);
1054         }
1055         MEM_freeN(jit2);
1056 }
1057
1058
1059 static void give_mesh_mvert(Mesh *me, DispListMesh *dlm, int nr, float *co, short *no, int seed2)
1060 {
1061         static float *jit=0;
1062         static int jitlevel=1;
1063         MVert *mvert, *mvertbase=NULL;
1064         MFace *mface, *mfacebase=NULL;
1065         float u, v, *v1, *v2, *v3, *v4;
1066         int totface=0, totvert=0, curface, curjit;
1067         short *n1, *n2, *n3, *n4;
1068         
1069         /* signal */
1070         if(me==0) {
1071                 if(jit) MEM_freeN(jit);
1072                 jit= 0;
1073                 return;
1074         }
1075         
1076         if(dlm) {
1077                 mvertbase= dlm->mvert;
1078                 mfacebase= dlm->mface;
1079                 totface= dlm->totface;
1080                 totvert= dlm->totvert;
1081         }
1082         
1083         if(totvert==0) {
1084                 mvertbase= me->mvert;
1085                 mfacebase= me->mface;
1086                 totface= me->totface;
1087                 totvert= me->totvert;
1088         }
1089         
1090         if(totface==0 || nr<totvert) {
1091                 mvert= mvertbase + (nr % totvert);
1092                 VECCOPY(co, mvert->co);
1093                 VECCOPY(no, mvert->no);
1094         }
1095         else {
1096
1097                 nr-= totvert;
1098                 
1099                 if(jit==0) {
1100                         jitlevel= nr/totface;
1101                         if(jitlevel==0) jitlevel= 1;
1102                         if(jitlevel>100) jitlevel= 100;
1103
1104                         jit= MEM_callocN(2+ jitlevel*2*sizeof(float), "jit");
1105                         init_mv_jit(jit, jitlevel,seed2);
1106                         
1107                 }
1108
1109                 curjit= nr/totface;
1110                 curjit= curjit % jitlevel;
1111
1112                 curface= nr % totface;
1113                 
1114                 mface= mfacebase;
1115                 mface+= curface;
1116
1117                 v1= (mvertbase+(mface->v1))->co;
1118                 v2= (mvertbase+(mface->v2))->co;
1119                 n1= (mvertbase+(mface->v1))->no;
1120                 n2= (mvertbase+(mface->v2))->no;
1121                 if(mface->v3==0) {
1122                         v3= (mvertbase+(mface->v2))->co;
1123                         v4= (mvertbase+(mface->v1))->co;
1124                         n3= (mvertbase+(mface->v2))->no;
1125                         n4= (mvertbase+(mface->v1))->no;
1126                 }
1127                 else if(mface->v4==0) {
1128                         v3= (mvertbase+(mface->v3))->co;
1129                         v4= (mvertbase+(mface->v1))->co;
1130                         n3= (mvertbase+(mface->v3))->no;
1131                         n4= (mvertbase+(mface->v1))->no;
1132                 }
1133                 else {
1134                         v3= (mvertbase+(mface->v3))->co;
1135                         v4= (mvertbase+(mface->v4))->co;
1136                         n3= (mvertbase+(mface->v3))->no;
1137                         n4= (mvertbase+(mface->v4))->no;
1138                 }
1139
1140                 u= jit[2*curjit];
1141                 v= jit[2*curjit+1];
1142
1143                 co[0]= (float)((1.0-u)*(1.0-v)*v1[0] + (1.0-u)*(v)*v2[0] + (u)*(v)*v3[0] + (u)*(1.0-v)*v4[0]);
1144                 co[1]= (float)((1.0-u)*(1.0-v)*v1[1] + (1.0-u)*(v)*v2[1] + (u)*(v)*v3[1] + (u)*(1.0-v)*v4[1]);
1145                 co[2]= (float)((1.0-u)*(1.0-v)*v1[2] + (1.0-u)*(v)*v2[2] + (u)*(v)*v3[2] + (u)*(1.0-v)*v4[2]);
1146                 
1147                 no[0]= (short)((1.0-u)*(1.0-v)*n1[0] + (1.0-u)*(v)*n2[0] + (u)*(v)*n3[0] + (u)*(1.0-v)*n4[0]);
1148                 no[1]= (short)((1.0-u)*(1.0-v)*n1[1] + (1.0-u)*(v)*n2[1] + (u)*(v)*n3[1] + (u)*(1.0-v)*n4[1]);
1149                 no[2]= (short)((1.0-u)*(1.0-v)*n1[2] + (1.0-u)*(v)*n2[2] + (u)*(v)*n3[2] + (u)*(1.0-v)*n4[2]);
1150                 
1151         }
1152 }
1153
1154
1155 void build_particle_system(Object *ob)
1156 {
1157         Base *base;
1158         Object *par;
1159         PartEff *paf;
1160         Particle *pa;
1161         Mesh *me;
1162         MVert *mvert;
1163         MTex *mtexmove=0;
1164         Material *ma;
1165         DispListMesh *dlm;
1166         float framelenont, ftime, dtime, force[3], imat[3][3], vec[3];
1167         float fac, prevobmat[4][4], sfraont, co[3];
1168         int deform=0, a, cur, cfraont, cfralast, totpart;
1169         short no[3];
1170
1171         if(ob->type!=OB_MESH) return;
1172         me= ob->data;
1173         if(me->totvert==0) return;
1174
1175         ma= give_current_material(ob, 1);
1176         if(ma) {
1177                 mtexmove= ma->mtex[7];
1178         }
1179
1180         paf= give_parteff(ob);
1181         if(paf==NULL) return;
1182
1183         waitcursor(1);
1184
1185         disable_speed_curve(1);
1186         
1187         /* warning! we cannot call this when modifier is active! */
1188         mesh_modifier(ob, 's');
1189
1190         /* generate all particles */
1191         if(paf->keys) MEM_freeN(paf->keys);
1192         paf->keys= NULL;
1193         new_particle(paf);
1194
1195         /* reset deflector cache, sumohandle is free, but its still sorta abuse... (ton) */
1196         for(base= G.scene->base.first; base; base= base->next) {
1197                 base->object->sumohandle= NULL;
1198         }
1199
1200         cfraont= G.scene->r.cfra;
1201         cfralast= -1000;
1202         framelenont= G.scene->r.framelen;
1203         G.scene->r.framelen= 1.0;
1204         sfraont= ob->sf;
1205         ob->sf= 0.0;
1206
1207         /* mult generations? */
1208         totpart= paf->totpart;
1209         for(a=0; a<PAF_MAXMULT; a++) {
1210                 if(paf->mult[a]!=0.0) {
1211                         /* interessant formula! this way after 'x' generations the total is paf->totpart */
1212                         totpart= (int)(totpart / (1.0+paf->mult[a]*paf->child[a]));
1213                 }
1214                 else break;
1215         }
1216
1217         ftime= paf->sta;
1218         dtime= (paf->end - paf->sta)/totpart;
1219
1220         /* remember full hierarchy */
1221         par= ob;
1222         while(par) {
1223                 pushdata(par, sizeof(Object));
1224                 par= par->parent;
1225         }
1226         
1227         /* for static particles, calculate system on current frame */
1228         if(ma) do_mat_ipo(ma);
1229
1230         /* set it all at first frame */
1231         G.scene->r.cfra= cfralast= (int)floor(ftime);
1232         par= ob;
1233         while(par) {
1234                 /* do_ob_ipo(par); */
1235                 do_ob_key(par);
1236                 par= par->parent;
1237         }
1238         
1239         if((paf->flag & PAF_STATIC)==0) {
1240                 if(ma) do_mat_ipo(ma);  // nor for static
1241                 
1242                 where_is_object(ob);
1243                 Mat4CpyMat4(prevobmat, ob->obmat);
1244                 Mat4Invert(ob->imat, ob->obmat);
1245                 Mat3CpyMat4(imat, ob->imat);
1246         }
1247         else {
1248                 Mat4One(prevobmat);
1249                 Mat3One(imat);
1250         }
1251         
1252         BLI_srand(paf->seed);
1253         
1254         /* otherwise it goes way too fast */
1255         force[0]= paf->force[0]*0.05f;
1256         force[1]= paf->force[1]*0.05f;
1257         force[2]= paf->force[2]*0.05f;
1258         
1259         if( paf->flag & PAF_STATIC ) deform= 0;
1260         else {
1261                 deform= (ob->parent && ob->parent->type==OB_LATTICE);
1262                 if(deform) init_latt_deform(ob->parent, 0);
1263         }
1264         
1265         /* init */
1266         if (mesh_uses_displist(me)) {
1267                 DerivedMesh *dm = mesh_get_derived(ob);
1268
1269                 dlm = dm->convertToDispListMesh(dm);
1270         } else {
1271                 dlm = NULL;
1272         }
1273
1274         give_mesh_mvert(me, dlm, totpart, co, no, paf->seed);
1275         
1276         if(G.f & G_DEBUG) {
1277                 printf("\n");
1278                 printf("Calculating particles......... \n");
1279         }
1280         for(a=0; a<totpart; a++, ftime+=dtime) {
1281                 
1282                 pa= new_particle(paf);
1283                 pa->time= ftime;
1284                 
1285                 if(G.f & G_DEBUG) {
1286                         int b, c;
1287                         
1288                         c = totpart/100;
1289                         if (c==0){
1290                                 c = 1;
1291                         }
1292
1293                         b=(a%c);
1294                         if (b==0) {
1295                                 printf("\r Particle: %d / %d ", a, totpart);
1296                                 fflush(stdout);
1297                         }
1298                 }
1299                 /* set ob at correct time */
1300                 
1301                 if((paf->flag & PAF_STATIC)==0) {
1302
1303                         cur= (int)floor(ftime) + 1 ;            /* + 1 has a reason: (obmat/prevobmat) otherwise comet-tails start too late */
1304                         if(cfralast != cur) {
1305                                 G.scene->r.cfra= cfralast= cur;
1306         
1307                                 /* added later: blur? */
1308                                 bsystem_time(ob, ob->parent, (float)G.scene->r.cfra, 0.0);
1309
1310                                 par= ob;
1311                                 while(par) {
1312                                         /* do_ob_ipo(par); */
1313                                         par->ctime= -1234567.0;
1314                                         do_ob_key(par);
1315                                         if(par->type==OB_ARMATURE) {
1316                                                 do_all_actions(par);    // only does this object actions
1317                                                 clear_object_constraint_status(par);    // mysterious call, otherwise do_actions doesnt work???
1318                                         }
1319                                         par= par->parent;
1320                                 }
1321                                 if(ma) do_mat_ipo(ma);
1322                                 Mat4CpyMat4(prevobmat, ob->obmat);
1323                                 where_is_object(ob);
1324                                 Mat4Invert(ob->imat, ob->obmat);
1325                                 Mat3CpyMat4(imat, ob->imat);
1326                         }
1327                 }
1328                 /* get coordinates */
1329                 if(paf->flag & PAF_FACE) give_mesh_mvert(me, dlm, a, co, no, paf->seed);
1330                 else {
1331                         mvert= me->mvert + (a % me->totvert);
1332                         VECCOPY(co, mvert->co);
1333                         VECCOPY(no, mvert->no);
1334                 }
1335                 
1336                 VECCOPY(pa->co, co);
1337                 
1338                 if(paf->flag & PAF_STATIC);
1339                 else {
1340                         Mat4MulVecfl(ob->obmat, pa->co);
1341                 
1342                         VECCOPY(vec, co);
1343                         Mat4MulVecfl(prevobmat, vec);
1344                         
1345                         /* first start speed: object */
1346                         VECSUB(pa->no, pa->co, vec);
1347                         VecMulf(pa->no, paf->obfac);
1348                         
1349                         /* calculate the correct inter-frame */ 
1350                         fac= (ftime- (float)floor(ftime));
1351                         pa->co[0]= fac*pa->co[0] + (1.0f-fac)*vec[0];
1352                         pa->co[1]= fac*pa->co[1] + (1.0f-fac)*vec[1];
1353                         pa->co[2]= fac*pa->co[2] + (1.0f-fac)*vec[2];
1354                 }
1355
1356                 /* start speed: normal */
1357                 if(paf->normfac!=0.0) {
1358                         /* sp= mvert->no; */
1359                                 /* transpose ! */
1360                         vec[0]= imat[0][0]*no[0] + imat[0][1]*no[1] + imat[0][2]*no[2];
1361                         vec[1]= imat[1][0]*no[0] + imat[1][1]*no[1] + imat[1][2]*no[2];
1362                         vec[2]= imat[2][0]*no[0] + imat[2][1]*no[1] + imat[2][2]*no[2];         
1363                 
1364                         Normalise(vec);
1365                         VecMulf(vec, paf->normfac);
1366                         VECADD(pa->no, pa->no, vec);
1367                 }
1368                 pa->lifetime= paf->lifetime;
1369                 if(paf->randlife!=0.0) {
1370                         pa->lifetime*= 1.0f+ (float)(paf->randlife*( BLI_drand() - 0.5));
1371                 }
1372                 pa->mat_nr= 1;
1373                 
1374                 make_particle_keys(0, a, paf, pa, force, deform, mtexmove, ob->lay);
1375         }
1376         
1377         if(G.f & G_DEBUG) {
1378                 printf("\r Particle: %d / %d \n", totpart, totpart);
1379                 fflush(stdout);
1380         }       
1381         if(deform) end_latt_deform();
1382                 
1383         /* restore */
1384         G.scene->r.cfra= cfraont;
1385         G.scene->r.framelen= framelenont;
1386         give_mesh_mvert(0, 0, 0, 0, 0,paf->seed);
1387
1388         /* put hierarchy back */
1389         par= ob;
1390         while(par) {
1391                 popfirst(par);
1392                 /* do not do ob->ipo: keep insertkey */
1393                 do_ob_key(par);
1394                 
1395                 if(par->type==OB_ARMATURE) {
1396                         do_all_actions(par);    // only does this object actions
1397                         clear_object_constraint_status(par);    // mysterious call, otherwise do_actions doesnt work???
1398                 }
1399                 par= par->parent;
1400         }
1401
1402         /* reset deflector cache */
1403         for(base= G.scene->base.first; base; base= base->next) {
1404                 if(base->object->sumohandle) {
1405                         MEM_freeN(base->object->sumohandle);
1406                         base->object->sumohandle= NULL;
1407                 }
1408         }
1409
1410         /* restore: AFTER popfirst */
1411         ob->sf= sfraont;
1412
1413         if(ma) do_mat_ipo(ma);  // set back on current time
1414         disable_speed_curve(0);
1415         
1416         mesh_modifier(ob, 'e');
1417
1418         waitcursor(0);
1419
1420         if (dlm) {
1421                 displistmesh_free(dlm);
1422         }
1423 }
1424
1425 /* ************* WAVE **************** */
1426
1427 void calc_wave_deform(WaveEff *wav, float ctime, float *co)
1428 {
1429         /* co is in local coords */
1430         float lifefac, x, y, amplit;
1431         
1432         /* actually this should not happen */
1433         if((wav->flag & (WAV_X+WAV_Y))==0) return;      
1434
1435         lifefac= wav->height;
1436         
1437         if( wav->lifetime!=0.0) {
1438                 x= ctime - wav->timeoffs;
1439                 if(x>wav->lifetime) {
1440                         
1441                         lifefac= x-wav->lifetime;
1442                         
1443                         if(lifefac > wav->damp) lifefac= 0.0;
1444                         else lifefac= (float)(wav->height*(1.0 - sqrt(lifefac/wav->damp)));
1445                 }
1446         }
1447         if(lifefac==0.0) return;
1448
1449         x= co[0]-wav->startx;
1450         y= co[1]-wav->starty;
1451
1452         if(wav->flag & WAV_X) {
1453                 if(wav->flag & WAV_Y) amplit= (float)sqrt( (x*x + y*y));
1454                 else amplit= x;
1455         }
1456         else amplit= y;
1457         
1458         /* this way it makes nice circles */
1459         amplit-= (ctime-wav->timeoffs)*wav->speed;
1460
1461         if(wav->flag & WAV_CYCL) {
1462                 amplit = (float)fmod(amplit-wav->width, 2.0*wav->width) + wav->width;
1463         }
1464
1465         /* GAUSSIAN */
1466         
1467         if(amplit> -wav->width && amplit<wav->width) {
1468         
1469                 amplit = amplit*wav->narrow;
1470                 amplit= (float)(1.0/exp(amplit*amplit) - wav->minfac);
1471
1472                 co[2]+= lifefac*amplit;
1473         }
1474 }
1475
1476 /* return 1 if deformed
1477    Note: it works on mvert now, so assumes to be callied in modifier stack \
1478 */
1479 int object_wave(Object *ob)
1480 {
1481         WaveEff *wav;
1482         Mesh *me;
1483         MVert *mvert;
1484         float ctime;
1485         int a;
1486         
1487         /* is there a wave */
1488         wav= ob->effect.first;
1489         while(wav) {
1490                 if(wav->type==EFF_WAVE) break;
1491                 wav= wav->next;
1492         }
1493         if(wav==NULL) return 0;
1494         
1495         if(ob->type==OB_MESH) {
1496
1497                 ctime= bsystem_time(ob, 0, (float)G.scene->r.cfra, 0.0);
1498                 
1499                 me= ob->data;
1500
1501                 wav= ob->effect.first;
1502                 while(wav) {
1503                         if(wav->type==EFF_WAVE) {
1504                                 
1505                                 /* precalculate */
1506                                 wav->minfac= (float)(1.0/exp(wav->width*wav->narrow*wav->width*wav->narrow));
1507                                 if(wav->damp==0) wav->damp= 10.0f;
1508                                 
1509                                 mvert= me->mvert;
1510                                 
1511                                 for(a=0; a<me->totvert; a++, mvert++) {
1512                                         calc_wave_deform(wav, ctime, mvert->co);
1513                                 }
1514                         }
1515                         wav= wav->next;
1516                 }
1517         }
1518         return 1;
1519 }
1520
1521 int SoftBodyDetectCollision(float opco[3], float npco[3], float colco[3],
1522         float facenormal[3], float *damp, float force[3], int mode,
1523         float cur_time, unsigned int par_layer,struct Object *vertexowner)
1524 {
1525         Base *base;
1526         Object *ob, *deflection_object = NULL;
1527         Mesh *def_mesh;
1528         MFace *mface, *deflection_face = NULL;
1529         float *v1, *v2, *v3, *v4, *vcache=NULL;
1530         float mat[3][3];
1531         float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3],d_nvect[3], obloc[3];
1532         float dv1[3], dv2[3], dv3[3];
1533         float facedist,n_mag,t,t2, min_t,force_mag_norm;
1534         int a, deflected=0, deflected_now=0;
1535         short cur_frame;
1536         int d_object=0, d_face=0, ds_object=0, ds_face=0;
1537
1538 // i'm going to rearrange it to declaration rules when WIP is finished (BM)
1539         float u,v,len_u,len_v;
1540         float innerfacethickness = -0.5f;
1541         float outerfacethickness = 0.2f;
1542         float ee = 5.0f;
1543         float ff = 0.1f;
1544         float fa;
1545
1546         min_t = 200000;
1547
1548         /* The first part of the code, finding the first intersected face*/
1549         base= G.scene->base.first;
1550         while (base) {
1551                 /*Only proceed for mesh object in same layer */
1552                 if(base->object->type==OB_MESH && (base->lay & par_layer)) {
1553                         ob= base->object;
1554                         if((vertexowner) && (ob == vertexowner)){ 
1555                                 /* if vertexowner is given 
1556                                  * we don't want to check collision with owner object */ 
1557                                 base = base->next;
1558                                 continue;
1559                         }
1560                         /* only with deflecting set */
1561                         if(ob->pd && ob->pd->deflect) {
1562                                 def_mesh= ob->data;
1563                         
1564                                 d_object = d_object + 1;
1565
1566                                 d_face = d_face + 1;
1567                                 mface= def_mesh->mface;
1568                                 a = def_mesh->totface;
1569 /* need to have user control for that since it depends on model scale */
1570                                 
1571                                 innerfacethickness =-ob->pd->pdef_sbift;
1572                                 outerfacethickness =ob->pd->pdef_sboft;
1573                                                 fa = (ff*outerfacethickness-outerfacethickness);
1574                                                 fa *= fa;
1575                                                 fa = 1.0f/fa;
1576                                 
1577                                 if(ob->parent==NULL && ob->ipo==NULL) { // static
1578                                         if(ob->sumohandle==NULL) cache_object_vertices(ob);
1579                                         vcache= ob->sumohandle;
1580                                 }
1581                                 else {
1582                                         /*Find out where the object is at this time*/
1583                                         cur_frame = G.scene->r.cfra;
1584                                         G.scene->r.cfra = (short)cur_time;
1585                                         where_is_object_time(ob, cur_time);
1586                                         G.scene->r.cfra = cur_frame;
1587                                         
1588                                         /*Pass the values from ob->obmat to mat*/
1589                                         /*and the location values to obloc           */
1590                                         Mat3CpyMat4(mat,ob->obmat);
1591                                         obloc[0] = ob->obmat[3][0];
1592                                         obloc[1] = ob->obmat[3][1];
1593                                         obloc[2] = ob->obmat[3][2];
1594                                         /* not cachable */
1595                                         vcache= NULL;
1596
1597                                 }
1598                                 
1599                                 while (a--) {
1600
1601                                         if(vcache) {
1602                                                 v1= vcache+ 3*(mface->v1);
1603                                                 VECCOPY(nv1, v1);
1604                                                 v1= vcache+ 3*(mface->v2);
1605                                                 VECCOPY(nv2, v1);
1606                                                 v1= vcache+ 3*(mface->v3);
1607                                                 VECCOPY(nv3, v1);
1608                                                 v1= vcache+ 3*(mface->v4);
1609                                                 VECCOPY(nv4, v1);
1610                                         }
1611                                         else {
1612                                                 /* Calculate the global co-ordinates of the vertices*/
1613                                                 v1= (def_mesh->mvert+(mface->v1))->co;
1614                                                 v2= (def_mesh->mvert+(mface->v2))->co;
1615                                                 v3= (def_mesh->mvert+(mface->v3))->co;
1616                                                 v4= (def_mesh->mvert+(mface->v4))->co;
1617         
1618                                                 VECCOPY(nv1, v1);
1619                                                 VECCOPY(nv2, v2);
1620                                                 VECCOPY(nv3, v3);
1621                                                 VECCOPY(nv4, v4);
1622         
1623                                                 /*Apply the objects deformation matrix*/
1624                                                 Mat3MulVecfl(mat, nv1);
1625                                                 Mat3MulVecfl(mat, nv2);
1626                                                 Mat3MulVecfl(mat, nv3);
1627                                                 Mat3MulVecfl(mat, nv4);
1628         
1629                                                 VECADD(nv1, nv1, obloc);
1630                                                 VECADD(nv2, nv2, obloc);
1631                                                 VECADD(nv3, nv3, obloc);
1632                                                 VECADD(nv4, nv4, obloc);
1633                                         }
1634                                         
1635                                         deflected_now = 0;
1636
1637                                         if (mode == 1){ // face intrusion test
1638                                                 // switch origin to be nv2
1639                                                 VECSUB(edge1, nv1, nv2);
1640                                                 VECSUB(edge2, nv3, nv2);
1641                                                 VECSUB(dv1,opco,nv2); // abuse dv1 to have vertex in question at *origin* of triangle
1642                                                 
1643                                                 len_u=Normalise(edge1);
1644                                                 len_v=Normalise(edge2);
1645                                                 
1646                                                 u = Inpf(dv1,edge1)/len_u;
1647                                                 v = Inpf(dv1,edge2)/len_v;
1648                                                 if ( (u >= 0.0f) && (v >= 0.0f) && ((u+v) <= 1.0)){ // inside prims defined by triangle and normal
1649                                                         Crossf(d_nvect, edge2, edge1);
1650                                                         n_mag = Normalise(d_nvect);
1651                                                         // ok lets  add force 
1652                                                         facedist = Inpf(dv1,d_nvect);
1653                                                         if ((facedist > innerfacethickness) && (facedist < outerfacethickness)){
1654                                                                 //force_mag_norm =ee*(facedist - outerfacethickness)*(facedist - outerfacethickness);
1655                                                                 force_mag_norm =(float)exp(-ee*facedist);
1656                                                                 if (facedist > outerfacethickness*ff)
1657                                                                 force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness);
1658
1659                                                                 force[0] += force_mag_norm*d_nvect[0] ;
1660                                                                 force[1] += force_mag_norm*d_nvect[1] ;
1661                                                                 force[2] += force_mag_norm*d_nvect[2] ;
1662                                                                 *damp=ob->pd->pdef_sbdamp;
1663                                                                 deflected = 2;
1664
1665                                                                 colco[0] = nv2[0] + len_u*u*edge1[0] + len_v*v*edge2[0];
1666                                                                 colco[1] = nv2[1] + len_u*u*edge1[1] + len_v*v*edge2[1];
1667                                                                 colco[2] = nv2[2] + len_u*u*edge1[2] + len_v*v*edge2[2];
1668
1669                                                         }
1670                                                 }               
1671                                                 if (mface->v4){ // quad
1672                                                         // switch origin to be nv4
1673                                                         VECSUB(edge1, nv3, nv4);
1674                                                         VECSUB(edge2, nv1, nv4);
1675                                                         VECSUB(dv1,opco,nv4); // abuse dv1 to have vertex in question at *origin* of triangle
1676                                                         
1677                                                         len_u=Normalise(edge1);
1678                                                         len_v=Normalise(edge2);
1679                                                         
1680                                                         u = Inpf(dv1,edge1)/len_u;
1681                                                         v = Inpf(dv1,edge2)/len_v;
1682                                                         if ( (u >= 0.0f) && (v >= 0.0f) && ((u+v) <= 1.0)){ // inside prims defined by triangle and normal
1683                                                                 Crossf(d_nvect, edge2, edge1);
1684                                                                 n_mag = Normalise(d_nvect);
1685                                                                 // ok lets  add force 
1686                                                                 facedist = Inpf(dv1,d_nvect);
1687                                                                 if ((facedist > innerfacethickness) && (facedist < outerfacethickness)){
1688                                                                         //force_mag_norm =ee*(facedist - outerfacethickness)*(facedist - outerfacethickness);
1689                                                                 force_mag_norm =(float)exp(-ee*facedist);
1690                                                                 if (facedist > outerfacethickness*ff)
1691                                                                 force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness);
1692
1693                                                                         force[0] += force_mag_norm*d_nvect[0] ;
1694                                                                         force[1] += force_mag_norm*d_nvect[1] ;
1695                                                                         force[2] += force_mag_norm*d_nvect[2] ;
1696                                                                         *damp=ob->pd->pdef_sbdamp;
1697                                                                         deflected = 2;
1698                                                                 colco[0] = nv4[0] + len_u*u*edge1[0] + len_v*v*edge2[0];
1699                                                                 colco[1] = nv4[1] + len_u*u*edge1[1] + len_v*v*edge2[1];
1700                                                                 colco[2] = nv4[2] + len_u*u*edge1[2] + len_v*v*edge2[2];
1701
1702                                                                 }
1703                                                         }
1704                                                         
1705                                                 }
1706                                         }
1707                                         if (mode == 2){ // edge intrusion test
1708                                                 
1709                                                 
1710 //                                      t= 0.5; // this is labda of line, can use it optimize quad intersection
1711 // sorry but no .. see below (BM)                                       
1712                                         if( linetriangle(opco, npco, nv1, nv2, nv3, &t) ) {
1713                                                 if (t < min_t) {
1714                                                         deflected = 1;
1715                                                         deflected_now = 1;
1716                                                 }
1717                                         }
1718 //                                      else if (mface->v4 && (t>=0.0 && t<=1.0)) {
1719 // no, you can't skip testing the other triangle
1720 // it might give a smaller t on (close to) the edge .. this is numerics not esoteric maths :)
1721 // note: the 2 triangles don't need to share a plane ! (BM)
1722                                         if (mface->v4) {
1723                                                 if( linetriangle(opco, npco, nv1, nv3, nv4, &t2) ) {
1724                                                         if (t2 < min_t) {
1725                                                                 deflected = 1;
1726                                                                 deflected_now = 2;
1727                                                         }
1728                                                 }
1729                                         }
1730                                         
1731                                         if ((deflected_now > 0) && ((t < min_t) ||(t2 < min_t))) {
1732                         min_t = t;
1733                         ds_object = d_object;
1734                                                 ds_face = d_face;
1735                                                 deflection_object = ob;
1736                                                 deflection_face = mface;
1737                                                 if (deflected_now==1) {
1738                         min_t = t;
1739                                                         VECCOPY(dv1, nv1);
1740                                                         VECCOPY(dv2, nv2);
1741                                                         VECCOPY(dv3, nv3);
1742                                                 }
1743                                                 else {
1744                         min_t = t2;
1745                                                         VECCOPY(dv1, nv1);
1746                                                         VECCOPY(dv2, nv3);
1747                                                         VECCOPY(dv3, nv4);
1748                                                 }
1749                                         }
1750                                         } // not -100
1751                                         mface++;
1752                                 }
1753                         }
1754                 }
1755                 base = base->next;
1756         } // while (base)
1757
1758         if (mode == 1){ // face 
1759                 return deflected;
1760                                                 }
1761
1762         if (mode == 2){ // edge intrusion test
1763                 if (deflected) {
1764                         VECSUB(edge1, dv1, dv2);
1765                         VECSUB(edge2, dv3, dv2);
1766                         Crossf(d_nvect, edge2, edge1);
1767                         n_mag = Normalise(d_nvect);
1768                         // return point of intersection
1769                         colco[0] = opco[0] + (min_t * (npco[0] - opco[0]));
1770                         colco[1] = opco[1] + (min_t * (npco[1] - opco[1]));
1771                         colco[2] = opco[2] + (min_t * (npco[2] - opco[2]));
1772                         
1773                         VECCOPY(facenormal,d_nvect);
1774                         {
1775                                 facenormal[0] *= -1.0f;                 
1776                                 facenormal[1] *= -1.0f;                 
1777                                 facenormal[2] *= -1.0f;                 
1778                         }
1779                         
1780                         
1781                 }
1782         }
1783         return deflected;
1784 }
1785