Added baking for softbodies.
[blender.git] / source / blender / blenkernel / intern / softbody.c
1 /*  softbody.c      
2  * 
3  * $Id: 
4  *
5  * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
6  *
7  * This program is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License
9  * as published by the Free Software Foundation; either version 2
10  * of the License, or (at your option) any later version. The Blender
11  * Foundation also sells licenses for use in proprietary software under
12  * the Blender License.  See http://www.blender.org/BL/ for information
13  * about this.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with this program; if not, write to the Free Software Foundation,
22  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23  *
24  * The Original Code is Copyright (C) Blender Foundation
25  * All rights reserved.
26  *
27  * The Original Code is: all of this file.
28  *
29  * Contributor(s): none yet.
30  *
31  * ***** END GPL/BL DUAL LICENSE BLOCK *****
32  */
33
34 /*
35 ******
36 variables on the UI for now
37
38         float mediafrict;  friction to env 
39         float nodemass;   softbody mass of *vertex* 
40         float grav;        softbody amount of gravitaion to apply 
41         
42         float goalspring;  softbody goal springs 
43         float goalfrict;   softbody goal springs friction 
44         float mingoal;     quick limits for goal 
45         float maxgoal;
46
47         float inspring;   softbody inner springs 
48         float infrict;     softbody inner springs friction 
49
50 *****
51 */
52
53
54 #include <math.h>
55 #include <stdlib.h>
56 #include <string.h>
57
58 #include "MEM_guardedalloc.h"
59
60 /* types */
61 #include "DNA_curve_types.h"
62 #include "DNA_object_types.h"
63 #include "DNA_object_force.h"   /* here is the softbody struct */
64 #include "DNA_key_types.h"
65 #include "DNA_mesh_types.h"
66 #include "DNA_meshdata_types.h"
67 #include "DNA_lattice_types.h"
68 #include "DNA_scene_types.h"
69
70 #include "BLI_blenlib.h"
71 #include "BLI_arithb.h"
72
73 #include "BKE_displist.h"
74 #include "BKE_effect.h"
75 #include "BKE_global.h"
76 #include "BKE_key.h"
77 #include "BKE_object.h"
78 #include "BKE_softbody.h"
79 #include "BKE_utildefines.h"
80
81 #include  "BIF_editdeform.h"
82
83 extern bDeformGroup *get_named_vertexgroup(Object *ob, char *name);
84 extern int  get_defgroup_num (Object *ob, bDeformGroup        *dg);
85
86
87 /* ********** soft body engine ******* */
88 #define SOFTGOALSNAP  0.999f 
89 // if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp
90 // removes *unnecessary* stiffnes from ODE system
91 #define HEUNWARNLIMIT 1 // 50 would be fine i think for detecting severe *stiff* stuff
92
93
94 float SoftHeunTol = 1.0f; // humm .. this should be calculated from sb parameters and sizes
95
96 /* local prototypes */
97 static void free_softbody_intern(SoftBody *sb);
98
99
100 /*+++ frame based timing +++*/
101
102 //physical unit of force is [kg * m / sec^2]
103
104 static float sb_grav_force_scale(Object *ob)
105 // since unit of g is [m/sec^2] and F = mass * g we rescale unit mass of node to 1 gramm
106 // put it to a function here, so we can add user options later without touching simulation code
107 {
108         return (0.001f);
109 }
110
111 static float sb_fric_force_scale(Object *ob)
112 // rescaling unit of drag [1 / sec] to somehow reasonable
113 // put it to a function here, so we can add user options later without touching simulation code
114 {
115         return (0.01f);
116 }
117
118 static float sb_time_scale(Object *ob)
119 // defining the frames to *real* time relation
120 {
121         SoftBody *sb= ob->soft; // is supposed to be there
122         if (sb){
123                 return(sb->physics_speed); //hrms .. this could be IPO as well :) 
124                 // estimated range [0.001 sluggish slug - 100.0 very fast (i hope ODE solver can handle that)]
125                 // 1 approx = a unit 1 pendulum at g = 9.8 [earth conditions]  has period 65 frames
126         // theory would give a 50 frames period .. so there must be something inaccurate .. looking for that (BM) 
127         }
128         return (1.0f);
129         /* 
130         this would be frames/sec independant timing assuming 25 fps is default
131         but does not work very well with NLA
132                 return (25.0f/G.scene->r.frs_sec)
133         */
134 }
135 /*--- frame based timing ---*/
136
137
138 static int count_mesh_quads(Mesh *me)
139 {
140         int a,result = 0;
141         MFace *mface= me->mface;
142         
143         if(mface) {
144                 for(a=me->totface; a>0; a--, mface++) {
145                         if(mface->v4) result++;
146                 }
147         }       
148         return result;
149 }
150
151 static void add_mesh_quad_diag_springs(Object *ob)
152 {
153         Mesh *me= ob->data;
154         MFace *mface= me->mface;
155         BodyPoint *bp;
156         BodySpring *bs, *bs_new;
157         int a ;
158         
159         if (ob->soft){
160                 int nofquads;
161                 
162                 nofquads = count_mesh_quads(me);
163                 if (nofquads) {
164                         /* resize spring-array to hold additional quad springs */
165                         bs_new= MEM_callocN( (ob->soft->totspring + nofquads *2 )*sizeof(BodySpring), "bodyspring");
166                         memcpy(bs_new,ob->soft->bspring,(ob->soft->totspring )*sizeof(BodySpring));
167                         
168                         if(ob->soft->bspring)
169                                 MEM_freeN(ob->soft->bspring); /* do this before reassigning the pointer  or have a 1st class memory leak */
170                         ob->soft->bspring = bs_new; 
171                         
172                         /* fill the tail */
173                         a = 0;
174                         bs = bs_new+ob->soft->totspring;
175                         bp= ob->soft->bpoint;
176                         if(mface ) {
177                                 for(a=me->totface; a>0; a--, mface++) {
178                                         if(mface->v4) {
179                                                 bs->v1= mface->v1;
180                                                 bs->v2= mface->v3;
181                                                 bs->strength= 1.0;
182                                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
183                                                 bs++;
184                                                 bs->v1= mface->v2;
185                                                 bs->v2= mface->v4;
186                                                 bs->strength= 1.0;
187                                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
188                                                 bs++;
189                                                 
190                                         }
191                                 }       
192                         }
193                         
194             /* now we can announce new springs */
195                         ob->soft->totspring += nofquads *2;
196                 }
197         }
198 }
199
200
201 static void add_bp_springlist(BodyPoint *bp,int springID)
202 {
203         int *newlist;
204         
205         if (bp->springs == NULL) {
206                 bp->springs = MEM_callocN( sizeof(int), "bpsprings");
207                 bp->springs[0] = springID;
208                 bp->nofsprings = 1;
209         }
210         else {
211                 bp->nofsprings++;
212                 newlist = MEM_callocN(bp->nofsprings * sizeof(int), "bpsprings");
213                 memcpy(newlist,bp->springs,(bp->nofsprings-1)* sizeof(int));
214                 MEM_freeN(bp->springs);
215                 bp->springs = newlist;
216                 bp->springs[bp->nofsprings-1] = springID;
217         }
218 }
219
220 /* do this once when sb is build
221 it is O(N^2) so scanning for springs every iteration is too expensive
222 */
223 static void build_bps_springlist(Object *ob)
224 {
225         SoftBody *sb= ob->soft; // is supposed to be there
226         BodyPoint *bp;  
227         BodySpring *bs; 
228         int a,b;
229         
230         if (sb==NULL) return; // paranoya check
231         
232         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
233                 /* scan for attached inner springs */   
234                 for(b=sb->totspring, bs= sb->bspring; b>0; b--, bs++) {
235                         if (( (sb->totpoint-a) == bs->v1) ){ 
236                                 add_bp_springlist(bp,sb->totspring -b);
237                         }
238                         if (( (sb->totpoint-a) == bs->v2) ){ 
239                                 add_bp_springlist(bp,sb->totspring -b);
240                         }
241                 }//for springs
242                 // if (bp->nofsprings) printf(" node %d has %d spring links\n",a,bp->nofsprings);
243         }//for bp               
244 }
245
246
247 /* creates new softbody if didn't exist yet, makes new points and springs arrays */
248 /* called in mesh_to_softbody */
249 static void renew_softbody(Object *ob, int totpoint, int totspring)  
250 {
251         SoftBody *sb;
252         
253         if(ob->soft==NULL) ob->soft= sbNew();
254         else free_softbody_intern(ob->soft);
255         sb= ob->soft;
256            
257         if(totpoint) {
258                 sb->totpoint= totpoint;
259                 sb->totspring= totspring;
260                 
261                 sb->bpoint= MEM_mallocN( totpoint*sizeof(BodyPoint), "bodypoint");
262                 if(totspring) 
263                         sb->bspring= MEM_mallocN( totspring*sizeof(BodySpring), "bodyspring");
264         }
265 }
266
267 static void free_softbody_baked(SoftBody *sb)
268 {
269         SBVertex *key;
270         int k;
271         
272         for(k=0; k<sb->totkey; k++) {
273                 key= *(sb->keys + k);
274                 if(key) MEM_freeN(key);
275         }
276         if(sb->keys) MEM_freeN(sb->keys);
277         
278         sb->keys= NULL;
279         sb->totkey= 0;
280         
281 }
282
283 /* only frees internal data */
284 static void free_softbody_intern(SoftBody *sb)
285 {
286         if(sb) {
287                 int a;
288                 BodyPoint *bp;
289                 
290                 if(sb->bpoint){
291                         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
292                                 /* free spring list */ 
293                                 if (bp->springs != NULL) {
294                                         MEM_freeN(bp->springs);
295                                 }
296                         }
297                         MEM_freeN(sb->bpoint);
298                 }
299                 
300                 if(sb->bspring) MEM_freeN(sb->bspring);
301                 
302                 sb->totpoint= sb->totspring= 0;
303                 sb->bpoint= NULL;
304                 sb->bspring= NULL;
305                 
306                 free_softbody_baked(sb);
307         }
308 }
309
310
311 /* ************ dynamics ********** */
312
313 /* aye this belongs to arith.c */
314 static void Vec3PlusStVec(float *v, float s, float *v1)
315 {
316         v[0] += s*v1[0];
317         v[1] += s*v1[1];
318         v[2] += s*v1[2];
319 }
320
321 static int sb_deflect_face(Object *ob,float *actpos, float *futurepos,float *collisionpos, float *facenormal,float *force,float *cf ,float *bounce)
322 {
323         int deflected;
324         float s_actpos[3], s_futurepos[3];
325         VECCOPY(s_actpos,actpos);
326         if(futurepos)
327         VECCOPY(s_futurepos,futurepos);
328         if (bounce) *bounce *= 1.5f;
329                                 
330                                 
331         deflected= SoftBodyDetectCollision(s_actpos, s_futurepos, collisionpos,
332                                         facenormal, cf, force , 1,
333                                         G.scene->r.cfra, ob->lay, ob);
334         return(deflected);
335                                 
336 }
337
338 /* for future use (BM)
339 static int sb_deflect_edge_face(Object *ob,float *actpos, float *futurepos,float *collisionpos, float *facenormal,float *slip ,float *bounce)
340 {
341         int deflected;
342         float dummy[3],s_actpos[3], s_futurepos[3];
343         SoftBody *sb= ob->soft; // is supposed to be there
344         VECCOPY(s_actpos,actpos);
345         VECCOPY(s_futurepos,futurepos);
346         if (slip)   *slip   *= 0.98f;
347         if (bounce) *bounce *= 1.5f;
348                                 
349                                 
350         deflected= SoftBodyDetectCollision(s_actpos, s_futurepos, collisionpos,
351                                         facenormal, dummy, dummy , 2,
352                                         G.scene->r.cfra, ob->lay, ob);
353         return(deflected);
354                                 
355 }
356 */
357 // some functions removed here .. to help HOS on next merge (BM)
358
359 #define USES_FIELD              1
360 #define USES_DEFLECT    2
361 static int is_there_deflection(unsigned int layer)
362 {
363         Base *base;
364         int retval= 0;
365         
366         for(base = G.scene->base.first; base; base= base->next) {
367                 if( (base->lay & layer) && base->object->pd) {
368                         if(base->object->pd->forcefield) retval |= USES_FIELD;
369                         if(base->object->pd->deflect) retval |= USES_DEFLECT;
370                 }
371         }
372         return retval;
373 }
374
375 static void softbody_calc_forces(Object *ob, float forcetime)
376 {
377 /* rule we never alter free variables :bp->vec bp->pos in here ! 
378  * this will ruin adaptive stepsize AKA heun! (BM) 
379  */
380         SoftBody *sb= ob->soft; // is supposed to be there
381         BodyPoint  *bp;
382         BodyPoint *bproot;
383         BodySpring *bs; 
384         float iks, ks, kd, gravity, actspringlen, forcefactor, sd[3],
385         fieldfactor = 1000.0f, 
386         windfactor  = 250.0f;   
387         int a, b, do_effector;
388         
389         /* clear forces */
390         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
391                 bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
392         }
393         
394         gravity = sb->grav * sb_grav_force_scale(ob);   
395         /* check! */
396         do_effector= is_there_deflection(ob->lay);
397         iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
398         bproot= sb->bpoint; /* need this for proper spring addressing */
399         
400         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
401                 if(bp->goal < SOFTGOALSNAP){ // ommit this bp when i snaps
402                         float auxvect[3]; // aux unit vector  
403                         float velgoal[3];
404                         float absvel =0, projvel= 0;
405                         
406                         /* do goal stuff */
407                         if(ob->softflag & OB_SB_GOAL) {
408                                 /* true elastic goal */
409                                 VecSubf(auxvect,bp->origT,bp->pos);
410                                 ks  = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
411                                 bp->force[0]= ks*(auxvect[0]);
412                                 bp->force[1]= ks*(auxvect[1]);
413                                 bp->force[2]= ks*(auxvect[2]);
414                                 /* calulate damping forces generated by goals*/
415                                 VecSubf(velgoal,bp->origS, bp->origE);
416                                 kd =  sb->goalfrict * sb_fric_force_scale(ob) ;
417                                 
418                                 if (forcetime > 0.0 ) { // make sure friction does not become rocket motor on time reversal
419                                         bp->force[0]-= kd * (velgoal[0] + bp->vec[0]);
420                                         bp->force[1]-= kd * (velgoal[1] + bp->vec[1]);
421                                         bp->force[2]-= kd * (velgoal[2] + bp->vec[2]);
422                                 }
423                                 else {
424                                         bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
425                                         bp->force[1]-= kd * (velgoal[1] - bp->vec[1]);
426                                         bp->force[2]-= kd * (velgoal[2] - bp->vec[2]);
427                                 }
428                         }
429                         /* done goal stuff */
430                         
431                         
432                         /* gravitation */
433                         bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
434                         
435                         /* particle field & vortex */
436                         if(do_effector & USES_FIELD) {
437                                 float force[3]= {0.0f, 0.0f, 0.0f};
438                                 float speed[3]= {0.0f, 0.0f, 0.0f};
439                                 float eval_sb_fric_force_scale = sb_fric_force_scale(ob); // just for calling functio once
440                                 
441                                 pdDoEffector(bp->pos, force, speed, (float)G.scene->r.cfra, ob->lay,PE_WIND_AS_SPEED);
442                                 // note: now we have wind as motion of media, so we can do anisotropic stuff here, 
443                                 // if we had vertex normals here(BM)
444                                 /* apply forcefield*/
445                                 VecMulf(force,fieldfactor* eval_sb_fric_force_scale); 
446                                 VECADD(bp->force, bp->force, force);
447                                 
448                                 /* friction in moving media */  
449                                 kd= sb->mediafrict* eval_sb_fric_force_scale;  
450                                 bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale);
451                                 bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale);
452                                 bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale);
453                                 /* now we'll have nice centrifugal effect for vortex */
454                                 
455                         }
456                         else {
457                                 /* friction in media (not) moving*/
458                                 kd= sb->mediafrict* sb_fric_force_scale(ob);  
459                                 /* assume it to be proportional to actual velocity */
460                                 bp->force[0]-= bp->vec[0]*kd;
461                                 bp->force[1]-= bp->vec[1]*kd;
462                                 bp->force[2]-= bp->vec[2]*kd;
463                                 /* friction in media done */
464                         }
465                         
466                         /*other forces*/
467                         /* this is the place where other forces can be added
468                         yes, constraints and collision stuff should go here too (read baraff papers on that!)
469                         */
470                         /* try to match moving collision targets */
471                         /* master switch to turn collision off (BM)*/
472                         //if(0) {
473                         if(do_effector & USES_DEFLECT) {
474                                 /*sorry for decl. here i'll move 'em up when WIP is done (BM) */
475                                 float defforce[3];
476                                 float collisionpos[3],facenormal[3];
477                                 float cf = 1.0f;
478                                 float bounce = 0.5f;
479                                 kd = 1.0f;
480                                 defforce[0] = 0.0f;
481                                 defforce[1] = 0.0f;
482                                 defforce[2] = 0.0f;
483                                 
484                                 if (sb_deflect_face(ob,bp->pos, bp->pos, collisionpos, facenormal,defforce,&cf,&bounce)){
485                                         bp->force[0] += defforce[0]*kd;
486                                         bp->force[1] += defforce[1]*kd;
487                                         bp->force[2] += defforce[2]*kd;
488                                         bp->contactfrict = cf;
489                                 }
490                                 else{ 
491                                         bp->contactfrict = 0.0f;
492                                 }
493                                 
494                         }
495                         
496                         /*other forces done*/
497                         /* nice things could be done with anisotropic friction
498                         like wind/air resistance in normal direction
499                         --> having a piece of cloth sailing down 
500                         but this needs to have a *valid* vertex normal
501                         *valid* means to be calulated on time axis
502                         hrms .. may be a rough one could be used as well .. let's see 
503                         */
504                         
505                         if(ob->softflag & OB_SB_EDGES) {
506                                 if (sb->bspring){ // spring list exists at all ? 
507                                         for(b=bp->nofsprings;b>0;b--){
508                                                 bs = sb->bspring + bp->springs[b-1];
509                                                 if (( (sb->totpoint-a) == bs->v1) ){ 
510                                                         actspringlen= VecLenf( (bproot+bs->v2)->pos, bp->pos);
511                                                         VecSubf(sd,(bproot+bs->v2)->pos, bp->pos);
512                                                         Normalise(sd);
513                                                         
514                                                         // friction stuff V1
515                                                         VecSubf(velgoal,bp->vec,(bproot+bs->v2)->vec);
516                                                         kd = sb->infrict * sb_fric_force_scale(ob);
517                                                         absvel  = Normalise(velgoal);
518                                                         projvel = ABS(Inpf(sd,velgoal));
519                                                         kd *= absvel * projvel;
520                                                         Vec3PlusStVec(bp->force,-kd,velgoal);
521                                                         
522                                                         if(bs->len > 0.0) /* check for degenerated springs */
523                                                                 forcefactor = (bs->len - actspringlen)/bs->len * iks;
524                                                         else
525                                                                 forcefactor = actspringlen * iks;
526                                                         
527                                                         Vec3PlusStVec(bp->force,-forcefactor,sd);
528                                                         
529                                                 }
530                                                 
531                                                 if (( (sb->totpoint-a) == bs->v2) ){ 
532                                                         actspringlen= VecLenf( (bproot+bs->v1)->pos, bp->pos);
533                                                         VecSubf(sd,bp->pos,(bproot+bs->v1)->pos);
534                                                         Normalise(sd);
535                                                         
536                                                         // friction stuff V2
537                                                         VecSubf(velgoal,bp->vec,(bproot+bs->v1)->vec);
538                                                         kd = sb->infrict * sb_fric_force_scale(ob);
539                                                         absvel  = Normalise(velgoal);
540                                                         projvel = ABS(Inpf(sd,velgoal));
541                                                         kd *= absvel * projvel;
542                                                         Vec3PlusStVec(bp->force,-kd,velgoal);
543                                                         
544                                                         if(bs->len > 0.0)
545                                                                 forcefactor = (bs->len - actspringlen)/bs->len * iks;
546                                                         else
547                                                                 forcefactor = actspringlen * iks;
548                                                         Vec3PlusStVec(bp->force,+forcefactor,sd);                                                       
549                                                 }
550                                         }/* loop springs */
551                                 }/* existing spring list */ 
552                         }/*any edges*/
553                 }/*omit on snap */
554         }/*loop all bp's*/
555 }
556
557 static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err)
558 {
559         /* time evolution */
560         /* actually does an explicit euler step mode == 0 */
561         /* or heun ~ 2nd order runge-kutta steps, mode 1,2 */
562         SoftBody *sb= ob->soft; // is supposed to be there
563         BodyPoint *bp;
564         float dx[3],dv[3];
565         float timeovermass;
566         float maxerr = 0.0;
567         int a, do_effector;
568
569     forcetime *= sb_time_scale(ob);
570         /* check! */
571         do_effector= is_there_deflection(ob->lay);
572     
573         // claim a minimum mass for vertex 
574         if (sb->nodemass > 0.09999f) timeovermass = forcetime/sb->nodemass;
575         else timeovermass = forcetime/0.09999f;
576         
577         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
578                 if(bp->goal < SOFTGOALSNAP){
579                         
580                         /* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces  + more forces*/
581                         /* the ( ... )' operator denotes derivate respective time */
582                         /* the euler step for velocity then becomes */
583                         /* v(t + dt) = v(t) + a(t) * dt */ 
584                         bp->force[0]*= timeovermass; /* individual mass of node here */ 
585                         bp->force[1]*= timeovermass;
586                         bp->force[2]*= timeovermass;
587                         /* some nasty if's to have heun in here too */
588                         VECCOPY(dv,bp->force); 
589                         if (mode == 1){
590                                 VECCOPY(bp->prevvec, bp->vec);
591                                 VECCOPY(bp->prevdv, dv);
592                         }
593                         if (mode ==2){
594                                 /* be optimistic and execute step */
595                                 bp->vec[0] = bp->prevvec[0] + 0.5f * (dv[0] + bp->prevdv[0]);
596                                 bp->vec[1] = bp->prevvec[1] + 0.5f * (dv[1] + bp->prevdv[1]);
597                                 bp->vec[2] = bp->prevvec[2] + 0.5f * (dv[2] + bp->prevdv[2]);
598                                 /* compare euler to heun to estimate error for step sizing */
599                                 maxerr = MAX2(maxerr,ABS(dv[0] - bp->prevdv[0]));
600                                 maxerr = MAX2(maxerr,ABS(dv[1] - bp->prevdv[1]));
601                                 maxerr = MAX2(maxerr,ABS(dv[2] - bp->prevdv[2]));
602                         }
603                         else {VECADD(bp->vec, bp->vec, bp->force);}
604
605                         /* so here is (x)'= v(elocity) */
606                         /* the euler step for location then becomes */
607                         /* x(t + dt) = x(t) + v(t) * dt */ 
608                         
609                         VECCOPY(dx,bp->vec);
610                         dx[0]*=forcetime ; 
611                         dx[1]*=forcetime ; 
612                         dx[2]*=forcetime ; 
613                         
614                         /* again some nasty if's to have heun in here too */
615                         if (mode ==1){
616                                 VECCOPY(bp->prevpos,bp->pos);
617                                 VECCOPY(bp->prevdx ,dx);
618                         }
619                         
620                         if (mode ==2){
621                                 bp->pos[0] = bp->prevpos[0] + 0.5f * ( dx[0] + bp->prevdx[0]);
622                                 bp->pos[1] = bp->prevpos[1] + 0.5f * ( dx[1] + bp->prevdx[1]);
623                                 bp->pos[2] = bp->prevpos[2] + 0.5f* ( dx[2] + bp->prevdx[2]);
624                                 maxerr = MAX2(maxerr,ABS(dx[0] - bp->prevdx[0]));
625                                 maxerr = MAX2(maxerr,ABS(dx[1] - bp->prevdx[1]));
626                                 maxerr = MAX2(maxerr,ABS(dx[2] - bp->prevdx[2]));
627 /* kind of hack .. while inside collision target .. make movement more *viscous* */
628                                 if (bp->contactfrict > 0.0f){
629                                         bp->vec[0] *= (1.0 - bp->contactfrict);
630                                         bp->vec[1] *= (1.0 - bp->contactfrict);
631                                         bp->vec[2] *= (1.0 - bp->contactfrict);
632                                 }
633                         }
634                         else { VECADD(bp->pos, bp->pos, dx);}
635 // experimental particle collision suff was here .. just to help HOS on next merge (BM)
636                 }//snap
637         } //for
638         if (err){ /* so step size will be controlled by biggest difference in slope */
639                 *err = maxerr;
640         }
641 }
642
643 /* used by heun when it overshoots */
644 static void softbody_restore_prev_step(Object *ob)
645 {
646         SoftBody *sb= ob->soft; // is supposed to be there
647         BodyPoint *bp;
648         int a;
649         
650         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
651                 VECCOPY(bp->vec, bp->prevvec);
652                 VECCOPY(bp->pos, bp->prevpos);
653         }
654 }
655
656
657 static void softbody_apply_goalsnap(Object *ob)
658 {
659         SoftBody *sb= ob->soft; // is supposed to be there
660         BodyPoint *bp;
661         int a;
662         
663         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
664                 if (bp->goal >= SOFTGOALSNAP){
665                         VECCOPY(bp->prevpos,bp->pos);
666                         VECCOPY(bp->pos,bp->origT);
667                 }               
668         }
669 }
670
671 /* unused */
672 #if 0
673 static void softbody_force_goal(Object *ob)
674 {
675         SoftBody *sb= ob->soft; // is supposed to be there
676         BodyPoint *bp;
677         int a;
678         
679         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {           
680                 VECCOPY(bp->pos,bp->origT);
681                 bp->vec[0] = bp->origE[0] - bp->origS[0];
682                 bp->vec[1] = bp->origE[1] - bp->origS[1];
683                 bp->vec[2] = bp->origE[2] - bp->origS[2];               
684         }
685 }
686 #endif
687
688 /* expects full initialized softbody */
689 static void interpolate_exciter(Object *ob, int timescale, int time)
690 {
691         SoftBody *sb= ob->soft;
692         BodyPoint *bp;
693         float f;
694         int a;
695
696         // note: i removed Mesh usage here, softbody should remain generic! (ton)
697         
698         f = (float)time/(float)timescale;
699         
700         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {   
701                 bp->origT[0] = bp->origS[0] + f*(bp->origE[0] - bp->origS[0]); 
702                 bp->origT[1] = bp->origS[1] + f*(bp->origE[1] - bp->origS[1]); 
703                 bp->origT[2] = bp->origS[2] + f*(bp->origE[2] - bp->origS[2]); 
704                 if (bp->goal >= SOFTGOALSNAP){
705                         bp->vec[0] = bp->origE[0] - bp->origS[0];
706                         bp->vec[1] = bp->origE[1] - bp->origS[1];
707                         bp->vec[2] = bp->origE[2] - bp->origS[2];
708                 }
709         }
710         
711         if(ob->softflag & OB_SB_EDGES) {
712                 /* hrms .. do springs alter their lenght ?
713                 bs= ob->soft->bspring;
714                 bp= ob->soft->bpoint;
715                 for(a=0; (a<me->totedge && a < ob->soft->totspring ); a++, bs++) {
716                         bs->len= VecLenf( (bp+bs->v1)->origT, (bp+bs->v2)->origT);
717                 }
718                 */
719         }
720 }
721
722
723 /* ************ convertors ********** */
724
725 /*  for each object type we need;
726     - xxxx_to_softbody(Object *ob)      : a full (new) copy
727         - xxxx_update_softbody(Object *ob)  : update refreshes current positions
728     - softbody_to_xxxx(Object *ob)      : after simulation, copy vertex locations back
729 */
730
731 /* helper  call */
732 static int object_has_edges(Object *ob) 
733 {
734         if(ob->type==OB_MESH) {
735                 Mesh *me= ob->data;
736                 if(me->medge) return 1;
737         }
738         else if(ob->type==OB_LATTICE) {
739                 ;
740         }
741         
742         return 0;
743 }
744
745 /* helper  call */
746 static void set_body_point(Object *ob, BodyPoint *bp, float *vec)
747 {
748         
749         VECCOPY(bp->pos, vec);
750         Mat4MulVecfl(ob->obmat, bp->pos);  // yep, sofbody is global coords
751         VECCOPY(bp->origS, bp->pos);
752         VECCOPY(bp->origE, bp->pos);
753         VECCOPY(bp->origT, bp->pos);
754         
755         bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0;
756         bp->weight= 1.0;
757         bp->goal= ob->soft->defgoal;
758         
759         bp->nofsprings= 0;
760         bp->springs= NULL;
761         bp->contactfrict = 0.0f;
762 }
763
764
765 /* copy original (new) situation in softbody, as result of matrices or deform */
766 /* is assumed to enter function with ob->soft, but can be without points */
767 static void mesh_update_softbody(Object *ob)
768 {
769         Mesh *me= ob->data;
770         MVert *mvert= me->mvert;
771 /*      MEdge *medge= me->medge;  */ /*unused*/
772         BodyPoint *bp;
773         int a;
774         
775         /* possible after a file read... */
776         if(ob->soft->bpoint==NULL) sbObjectToSoftbody(ob);
777         
778         if(me->totvert) {
779         
780                 bp= ob->soft->bpoint;
781                 for(a=0; a<me->totvert; a++, mvert++, bp++) {
782                         VECCOPY(bp->origS, bp->origE);
783                         VECCOPY(bp->origE, mvert->co);
784                         Mat4MulVecfl(ob->obmat, bp->origE);
785                         VECCOPY(bp->origT, bp->origE);
786                 }
787                 
788                 if(ob->softflag & OB_SB_EDGES) {
789                         
790                         /* happens when in UI edges was set */
791                         if(ob->soft->bspring==NULL) 
792                                 if(object_has_edges(ob)) sbObjectToSoftbody(ob);
793                 
794                         /* hrms .. do springs alter their lenght ? (yes, mesh keys would (ton))
795                         if(medge) {
796                                 bs= ob->soft->bspring;
797                                 bp= ob->soft->bpoint;
798                                 for(a=0; (a<me->totedge && a < ob->soft->totspring ); a++, medge++, bs++) { 
799                                         bs->len= VecLenf( (bp+bs->v1)->origE, (bp+bs->v2)->origE);
800                                 }
801                         }
802                         */
803                 }
804         }
805 }
806
807
808 static void get_scalar_from_vertexgroup(Object *ob, int vertID, short groupindex, float *target)
809 /* result 0 on success, else indicates error number
810 -- kind of *inverse* result defintion,
811 -- but this way we can signal error condition to caller  
812 -- and yes this function must not be here but in a *vertex group module*
813 */
814 {
815         MDeformVert *dv;
816         int i;
817         
818         /* spot the vert in deform vert list at mesh */
819         if(ob->type==OB_MESH) {
820                 if (((Mesh *)ob->data)->dvert) {
821                         dv = ((Mesh*)ob->data)->dvert + vertID; 
822                         /* Lets see if this vert is in the weight group */
823                         for (i=0; i<dv->totweight; i++){
824                                 if (dv->dw[i].def_nr == groupindex){
825                                         *target= dv->dw[i].weight; /* got it ! */
826                                         break;
827                                 }
828                         }
829                 }
830         }
831
832
833 /* makes totally fresh start situation */
834 static void mesh_to_softbody(Object *ob)
835 {
836         SoftBody *sb;
837         Mesh *me= ob->data;
838         MVert *mvert= me->mvert;
839         MEdge *medge= me->medge;
840         BodyPoint *bp;
841         BodySpring *bs;
842         float goalfac;
843         int a, totedge;
844         
845         if (ob->softflag & OB_SB_EDGES) totedge= me->totedge;
846         else totedge= 0;
847         
848         /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
849         renew_softbody(ob, me->totvert, totedge);
850                 
851         /* we always make body points */
852         sb= ob->soft;   
853         bp= sb->bpoint;
854         goalfac= ABS(sb->maxgoal - sb->mingoal);
855         
856         for(a=me->totvert; a>0; a--, mvert++, bp++) {
857                 
858                 set_body_point(ob, bp, mvert->co);
859                 
860                 /* get scalar values needed  *per vertex* from vertex group functions,
861                 so we can *paint* them nicly .. 
862                 they are normalized [0.0..1.0] so may be we need amplitude for scale
863                 which can be done by caller but still .. i'd like it to go this way 
864                 */ 
865
866                 if(sb->vertgroup) {
867                         get_scalar_from_vertexgroup(ob, me->totvert - a, sb->vertgroup-1, &bp->goal);
868                         // do this always, regardless successfull read from vertex group
869                         bp->goal= sb->mingoal + bp->goal*goalfac;
870                 }
871                 /* a little ad hoc changing the goal control to be less *sharp* */
872                 bp->goal = (float)pow(bp->goal, 4.0f);
873                         
874                 /* to proove the concept
875                 this would enable per vertex *mass painting*
876                 strcpy(name,"SOFTMASS");
877                 error = get_scalar_from_named_vertexgroup(ob,name,me->totvert - a,&temp);
878                 if (!error) bp->mass = temp * ob->rangeofmass;
879                 */
880         }
881
882         /* but we only optionally add body edge springs */
883         if (ob->softflag & OB_SB_EDGES) {
884                 if(medge) {
885                         bs= sb->bspring;
886                         bp= sb->bpoint;
887                         for(a=me->totedge; a>0; a--, medge++, bs++) {
888                                 bs->v1= medge->v1;
889                                 bs->v2= medge->v2;
890                                 bs->strength= 1.0;
891                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
892                         }
893
894                 
895                         /* insert *diagonal* springs in quads if desired */
896                         if (ob->softflag & OB_SB_QUADS) {
897                                 add_mesh_quad_diag_springs(ob);
898                         }
899
900                         build_bps_springlist(ob); /* big mesh optimization */
901                 }
902         }
903         
904 }
905
906 /* copies current sofbody position in mesh, so do this within modifier stacks! */
907 static void softbody_to_mesh(Object *ob)
908 {
909         Mesh *me= ob->data;
910         MVert *mvert;
911         BodyPoint *bp;
912         int a;
913
914         bp= ob->soft->bpoint;
915         mvert= me->mvert;
916         for(a=me->totvert; a>0; a--, mvert++, bp++) {
917                 VECCOPY(mvert->co, bp->pos);
918                 Mat4MulVecfl(ob->imat, mvert->co);      // softbody is in global coords
919         }
920
921 }
922
923 /* makes totally fresh start situation */
924 static void lattice_to_softbody(Object *ob)
925 {
926         SoftBody *sb;
927         Lattice *lt= ob->data;
928         BodyPoint *bop;
929         BPoint *bp;
930         int a, totvert;
931         
932         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
933         
934         /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
935         renew_softbody(ob, totvert, 0);
936         
937         /* we always make body points */
938         sb= ob->soft;   
939         
940         for(a= totvert, bp= lt->def, bop= sb->bpoint; a>0; a--, bp++, bop++) {
941                 set_body_point(ob, bop, bp->vec);
942         }
943 }
944
945 /* copies current sofbody position */
946 static void softbody_to_lattice(Object *ob)
947 {
948         SoftBody *sb;
949         Lattice *lt= ob->data;
950         BodyPoint *bop;
951         BPoint *bp;
952         int a, totvert;
953         
954         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
955         sb= ob->soft;   
956         
957         for(a= totvert, bp= lt->def, bop= sb->bpoint; a>0; a--, bp++, bop++) {
958                 VECCOPY(bp->vec, bop->pos);
959                 Mat4MulVecfl(ob->imat, bp->vec);        // softbody is in global coords
960         }
961 }
962
963 /* copy original (new) situation in softbody, as result of matrices or deform */
964 /* is assumed to enter function with ob->soft, but can be without points */
965 static void lattice_update_softbody(Object *ob)
966 {
967         Lattice *lt= ob->data;
968         BodyPoint *bop;
969         BPoint *bp;
970         int a, totvert;
971         
972         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
973         
974         /* possible after a file read... */
975         if(ob->soft->bpoint==NULL) sbObjectToSoftbody(ob);
976         
977         for(a= totvert, bp= lt->def, bop= ob->soft->bpoint; a>0; a--, bp++, bop++) {
978                 VECCOPY(bop->origS, bop->origE);
979                 VECCOPY(bop->origE, bp->vec);
980                 Mat4MulVecfl(ob->obmat, bop->origE);
981                 VECCOPY(bop->origT, bop->origE);
982         }       
983 }
984
985
986 /* copies softbody result back in object */
987 /* only used in sbObjectStep() */
988 static void softbody_to_object(Object *ob)
989 {
990         
991         if(ob->soft==NULL) return;
992         
993         /* inverse matrix is not uptodate... */
994         Mat4Invert(ob->imat, ob->obmat);
995         
996         switch(ob->type) {
997         case OB_MESH:
998                 softbody_to_mesh(ob);
999                 break;
1000         case OB_LATTICE:
1001                 softbody_to_lattice(ob);
1002                 break;
1003         }
1004 }
1005
1006 /* copy original (new) situation in softbody, as result of matrices or deform */
1007 /* used in sbObjectStep() and sbObjectReset() */
1008 /* assumes to have ob->soft, but can be entered without points */
1009 static void object_update_softbody(Object *ob)
1010 {
1011         
1012         switch(ob->type) {
1013         case OB_MESH:
1014                 mesh_update_softbody(ob);
1015                 break;
1016         case OB_LATTICE:
1017                 lattice_update_softbody(ob);
1018                 break;
1019         }
1020         
1021 }
1022
1023 /* return 1 if succesfully baked and applied step */
1024 static int softbody_baked_step(Object *ob, float framenr)
1025 {
1026         SoftBody *sb= ob->soft;
1027         SBVertex *key0, *key1, *key2, *key3;
1028         BodyPoint *bp;
1029         float data[4], sfra, efra, cfra, dfra, fac;     // start, end, current, delta 
1030         int ofs1, a;
1031
1032         /* precondition check */
1033         if(sb==NULL || sb->keys==NULL || sb->totkey==0) return 0;
1034         /* so we got keys, but no bodypoints... even without simul we need it for the bake */
1035         if(sb->bpoint==NULL) sb->bpoint= MEM_callocN( sb->totpoint*sizeof(BodyPoint), "bodypoint");
1036         
1037         /* convert cfra time to system time */
1038         sfra= (float)sb->sfra;
1039         cfra= bsystem_time(ob, NULL, framenr, 0.0);
1040         efra= (float)sb->efra;
1041         dfra= (float)sb->interval;
1042
1043         /* offset in keys array */
1044         ofs1= floor( (cfra-sfra)/dfra );
1045
1046         if(ofs1 < 0) {
1047                 key0=key1=key2=key3= *sb->keys;
1048         }
1049         else if(ofs1 >= sb->totkey-1) {
1050                 key0=key1=key2=key3= *(sb->keys+sb->totkey-1);
1051         }
1052         else {
1053                 key1= *(sb->keys+ofs1);
1054                 key2= *(sb->keys+ofs1+1);
1055
1056                 if(ofs1>0) key0= *(sb->keys+ofs1-1);
1057                 else key0= key1;
1058                 if(ofs1<sb->totkey-1) key3= *(sb->keys+ofs1+1);
1059                 else key3= key2;
1060         }
1061         
1062         sb->ctime= cfra;        // needed?
1063         
1064         /* timing */
1065         fac= ((cfra-sfra)/dfra) - (float)ofs1;
1066         CLAMP(fac, 0.0, 1.0);
1067         set_four_ipo(fac, data, KEY_BSPLINE);
1068         
1069         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++, key0++, key1++, key2++, key3++) {
1070                 bp->pos[0]= data[0]*key0->vec[0] +  data[1]*key1->vec[0] + data[2]*key2->vec[0] + data[3]*key3->vec[0];
1071                 bp->pos[1]= data[0]*key0->vec[1] +  data[1]*key1->vec[1] + data[2]*key2->vec[1] + data[3]*key3->vec[1];
1072                 bp->pos[2]= data[0]*key0->vec[2] +  data[1]*key1->vec[2] + data[2]*key2->vec[2] + data[3]*key3->vec[2];
1073         }
1074         
1075         softbody_to_object(ob);
1076         
1077         return 1;
1078 }
1079
1080 /* only gets called after succesfully doing softbody_step */
1081 /* already checked for OB_SB_BAKE flag */
1082 static void softbody_baked_add(Object *ob, float framenr)
1083 {
1084         SoftBody *sb= ob->soft;
1085         SBVertex *key;
1086         BodyPoint *bp;
1087         float sfra, efra, cfra, dfra, fac1;     // start, end, current, delta 
1088         int ofs1, a;
1089         
1090         /* convert cfra time to system time */
1091         sfra= (float)sb->sfra;
1092         cfra= bsystem_time(ob, NULL, framenr, 0.0);
1093         efra= (float)sb->efra;
1094         dfra= (float)sb->interval;
1095         
1096         if(sb->totkey==0) {
1097                 if(sb->sfra >= sb->efra) return;                // safety, UI or py setting allows
1098                 if(sb->interval<1) sb->interval= 1;             // just be sure
1099                 
1100                 sb->totkey= 1 + (int)(ceil( (efra-sfra)/dfra ) );
1101                 sb->keys= MEM_callocN( sizeof(void *)*sb->totkey, "sb keys");
1102         }
1103         
1104         /* now find out if we have to store a key */
1105         
1106         /* offset in keys array */
1107         if(cfra==efra) {
1108                 ofs1= sb->totkey-1;
1109                 fac1= 0.0;
1110         }
1111         else {
1112                 ofs1= floor( (cfra-sfra)/dfra );
1113                 fac1= ((cfra-sfra)/dfra) - (float)ofs1;
1114         }       
1115         if( fac1 < 1.0/dfra ) {
1116                 
1117                 key= *(sb->keys+ofs1);
1118                 if(key == NULL) {
1119                         *(sb->keys+ofs1)= key= MEM_mallocN(sb->totpoint*sizeof(SBVertex), "softbody key");
1120                         
1121                         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++, key++) {
1122                                 VECCOPY(key->vec, bp->pos);
1123                         }
1124                 }
1125         }
1126 }
1127
1128 /* ************ Object level, exported functions *************** */
1129
1130 /* allocates and initializes general main data */
1131 SoftBody *sbNew(void)
1132 {
1133         SoftBody *sb;
1134         
1135         sb= MEM_callocN(sizeof(SoftBody), "softbody");
1136         
1137         sb->mediafrict= 0.5; 
1138         sb->nodemass= 1.0;
1139         sb->grav= 0.0; 
1140         sb->physics_speed= 1.0;
1141         sb->rklimit= 0.1;
1142
1143         sb->goalspring= 0.5; 
1144         sb->goalfrict= 0.0; 
1145         sb->mingoal= 0.0;  
1146         sb->maxgoal= 1.0;
1147         sb->defgoal= 0.7;
1148         
1149         sb->inspring= 0.5;
1150         sb->infrict= 0.5; 
1151         
1152         sb->interval= 10;
1153         sb->sfra= G.scene->r.sfra;
1154         sb->efra= G.scene->r.efra;
1155         
1156         return sb;
1157 }
1158
1159 /* frees all */
1160 void sbFree(SoftBody *sb)
1161 {
1162         free_softbody_intern(sb);
1163         MEM_freeN(sb);
1164 }
1165
1166
1167 /* makes totally fresh start situation */
1168 void sbObjectToSoftbody(Object *ob)
1169 {
1170
1171         switch(ob->type) {
1172         case OB_MESH:
1173                 mesh_to_softbody(ob);
1174                 break;
1175         case OB_LATTICE:
1176                 lattice_to_softbody(ob);
1177                 break;
1178         }
1179         
1180         if(ob->soft) ob->soft->ctime= bsystem_time(ob, NULL, (float)G.scene->r.cfra, 0.0);
1181         ob->softflag &= ~OB_SB_REDO;
1182 }
1183
1184 /* reset all motion */
1185 void sbObjectReset(Object *ob)
1186 {
1187         SoftBody *sb= ob->soft;
1188         BodyPoint *bp;
1189         int a;
1190         
1191         if(sb==NULL) return;
1192         sb->ctime= bsystem_time(ob, NULL, (float)G.scene->r.cfra, 0.0);
1193         
1194         object_update_softbody(ob);
1195         
1196         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
1197                 // origS is previous timestep
1198                 VECCOPY(bp->origS, bp->origE);
1199                 VECCOPY(bp->pos, bp->origE);
1200                 VECCOPY(bp->origT, bp->origE);
1201                 bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0f;
1202
1203                 // no idea about the Heun stuff! (ton)
1204                 VECCOPY(bp->prevpos, bp->pos);
1205                 VECCOPY(bp->prevvec, bp->vec);
1206                 VECCOPY(bp->prevdx, bp->vec);
1207                 VECCOPY(bp->prevdv, bp->vec);
1208         }
1209 }
1210
1211
1212 /* simulates one step. framenr is in frames */
1213 /* copies result back to object, displist */
1214 void sbObjectStep(Object *ob, float framenr)
1215 {
1216         SoftBody *sb;
1217         Base *base;
1218         float dtime;
1219         int timescale,t;
1220         float ctime, forcetime;
1221         float err;
1222
1223         /* this variable is set while transform(). with lattices also having a softbody, 
1224            it calls lattice_modifier() all the time... has no displist yet. Is temporal
1225            hack which should be resolved with proper depgraph usage + storage of deformed
1226            vertices in lattice (ton) */
1227         if(G.moving) return;
1228         
1229         /* baking works with global time */
1230         if(!(ob->softflag & OB_SB_BAKEDO) )
1231                 if(softbody_baked_step(ob, framenr) ) return;
1232         
1233         /* remake softbody if: */
1234         if( (ob->softflag & OB_SB_REDO) ||              // signal after weightpainting
1235                 (ob->soft==NULL) ||                                     // just to be nice we allow full init
1236                 (ob->soft->bpoint==NULL) )                      // after reading new file, or acceptable as signal to refresh
1237                         sbObjectToSoftbody(ob);
1238         
1239         sb= ob->soft;
1240
1241         /* still no points? go away */
1242         if(sb->totpoint==0) return;
1243         
1244         /* reset deflector cache, sumohandle is free, but its still sorta abuse... (ton) */
1245         for(base= G.scene->base.first; base; base= base->next) {
1246                 base->object->sumohandle= NULL;
1247         }
1248
1249         /* checking time: */
1250         ctime= bsystem_time(ob, NULL, framenr, 0.0);
1251         dtime= ctime - sb->ctime;
1252                 // bail out for negative or for large steps
1253         if(dtime<0.0 || dtime >= 9.9*G.scene->r.framelen) { // G.scene->r.framelen corrects for frame-mapping, so this is actually 10 frames for UI
1254                 sbObjectReset(ob);
1255                 return;
1256         }
1257         
1258         /* the simulator */
1259         
1260         if(dtime > 0.0) {       // note: what does this mean now? (ton)
1261                 //answer (BM) :
1262                 //dtime is still in [frames]
1263                 //we made sure dtime is >= 0.0
1264                 //but still need to handle dtime == 0.0 -> just return sb as is, just to be nice
1265                 object_update_softbody(ob);
1266                 
1267                 if (TRUE) {     // RSOL1 always true now (ton)
1268                         /* special case of 2nd order Runge-Kutta type AKA Heun */
1269                         float timedone =0.0; // how far did we get without violating error condition
1270                         /* loops = counter for emergency brake
1271                          * we don't want to lock up the system if physics fail
1272                          */
1273                         int loops =0 ; 
1274                         SoftHeunTol = sb->rklimit; // humm .. this should be calculated from sb parameters and sizes
1275
1276                         forcetime = dtime; /* hope for integrating in one step */
1277                         while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) )
1278                         {
1279                                 if (ABS(dtime) > 9.0 ){
1280                                         if(G.f & G_DEBUG) printf("SB_STEPSIZE \n");
1281                                         break; // sorry but i must assume goal movement can't be interpolated any more
1282                                 }
1283                                 //set goals in time
1284                                 interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime)));
1285                                 // do predictive euler step
1286                                 softbody_calc_forces(ob, forcetime);
1287                                 softbody_apply_forces(ob, forcetime, 1, NULL);
1288                                 // crop new slope values to do averaged slope step
1289                                 softbody_calc_forces(ob, forcetime);
1290                                 softbody_apply_forces(ob, forcetime, 2, &err);
1291                                 softbody_apply_goalsnap(ob);
1292
1293                                 if (err > SoftHeunTol){ // error needs to be scaled to some quantity
1294                                         softbody_restore_prev_step(ob);
1295                                         forcetime /= 2.0;
1296                                 }
1297                                 else {
1298                                         float newtime = forcetime * 1.1f; // hope for 1.1 times better conditions in next step
1299                                         
1300                                         if (err > SoftHeunTol/2.0){ // stay with this stepsize unless err really small
1301                                                 newtime = forcetime;
1302                                         }
1303                                         timedone += forcetime;
1304                                         if (forcetime > 0.0)
1305                                                 forcetime = MIN2(dtime - timedone,newtime);
1306                                         else 
1307                                                 forcetime = MAX2(dtime - timedone,newtime);
1308                                 }
1309                                 loops++;
1310                         }
1311                         // move snapped to final position
1312                         interpolate_exciter(ob, 2, 2);
1313                         softbody_apply_goalsnap(ob);
1314                         
1315                         if(G.f & G_DEBUG) {
1316                                 if (loops > HEUNWARNLIMIT) /* monitor high loop counts say 1000 after testing */
1317                                         printf("%d heun integration loops/frame \n",loops);
1318                         }
1319                 }
1320                 else{
1321                         /* do brute force explicit euler */
1322                         /* inner intagration loop */
1323                         /* */
1324                         // loop n times so that n*h = duration of one frame := 1
1325                         // x(t+h) = x(t) + h*v(t);
1326                         // v(t+h) = v(t) + h*f(x(t),t);
1327                         timescale = (int)(sb->rklimit * ABS(dtime)); 
1328                         for(t=1 ; t <= timescale; t++) {
1329                                 if (ABS(dtime) > 15 ) break;
1330                                 
1331                                 /* the *goal* mesh must use the n*h timing too !
1332                                 use *cheap* linear intepolation for that  */
1333                                 interpolate_exciter(ob,timescale,t);                    
1334                                 if (timescale > 0 ) {
1335                                         forcetime = dtime/timescale;
1336                                         
1337                                         /* does not fit the concept sloving ODEs :) */
1338                                         /*                      softbody_apply_goal(ob,forcetime );  */
1339                                         
1340                                         /* explicit Euler integration */
1341                                         /* we are not controling a nuclear power plant! 
1342                                         so rought *almost* physical behaviour is acceptable.
1343                                         in cases of *mild* stiffnes cranking up timscale -> decreasing stepsize *h*
1344                                         avoids instability      */
1345                                         softbody_calc_forces(ob,forcetime);
1346                                         softbody_apply_forces(ob,forcetime,0, NULL);
1347                                         softbody_apply_goalsnap(ob);
1348                                         
1349                                         //      if (0){
1350                                         /* ok here comes the ├╝berhammer
1351                                         use a semi implicit euler integration to tackle *all* stiff conditions 
1352                                         but i doubt the cost/benifit holds for most of the cases
1353                                         -- to be coded*/
1354                                         //      }
1355                                         
1356                                 }
1357                         }
1358                 }
1359                 
1360                 /* and apply to vertices */
1361                  softbody_to_object(ob);
1362                 
1363                 sb->ctime= ctime;
1364         } // if(ABS(dtime) > 0.0) 
1365         else {
1366                 // rule : you have asked for the current state of the softobject 
1367                 // since dtime= ctime - ob->soft->ctime== 0.0;
1368                 // and we were not notifified about any other time changes 
1369                 // so here it is !
1370                 softbody_to_object(ob);
1371         }
1372
1373         /* reset deflector cache */
1374         for(base= G.scene->base.first; base; base= base->next) {
1375                 if(base->object->sumohandle) {
1376                         MEM_freeN(base->object->sumohandle);
1377                         base->object->sumohandle= NULL;
1378                 }
1379         }
1380         
1381         if(ob->softflag & OB_SB_BAKEDO) softbody_baked_add(ob, framenr);
1382 }
1383