1 /*
2  * Tux Racer
3  * Copyright (C) 1999-2001 Jasmin F. Patry
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
18  */
19 
20 
21 #include "tuxracer.h"
22 #include "course_load.h"
23 #include "course_render.h"
24 #include "hier.h"
25 #include "hier_util.h"
26 #include "tux.h"
27 #include "part_sys.h"
28 #include "phys_sim.h"
29 #include "nmrcl.h"
30 #include "audio.h"
31 #include "loop.h"
32 
33 /*
34  * Constants
35  */
36 
37 /* Tux's mass (kg) */
38 #define TUX_MASS 20
39 
40 /* Width of Tux (m) */
41 #define TUX_WIDTH 0.45
42 
43 /* Minimum speed (m/s) */
44 #define MIN_TUX_SPEED 1.4
45 
46 /* Minimum speed for friction to take effect (m/s) */
47 #define MIN_FRICTION_SPEED 2.8
48 
49 /* Initial Tux speed (m/s) */
50 #define INIT_TUX_SPEED 3.0
51 
52 /* Maximum frictional force */
53 #define MAX_FRICTIONAL_FORCE 800
54 
55 /* Maximum angle by which frictional force is rotated when turning*/
56 #define MAX_TURN_ANGLE 45
57 
58 /* Maximum turning force perpendicular to direction of movement */
59 #define MAX_TURN_PERPENDICULAR_FORCE 400
60 
61 /* Maximum increase in friction when turning */
62 #define MAX_TURN_PENALTY 0.15
63 
64 /* Force applied by Tux when braking (N) */
65 #define BRAKE_FORCE 200
66 
67 /* Maximum roll angle (degrees) */
68 #define MAX_ROLL_ANGLE 30
69 
70 /* Roll angle while braking (degrees) */
71 #define BRAKING_ROLL_ANGLE 55
72 
73 /* Baseline friction coeff for rolling */
74 #define IDEAL_ROLL_FRIC_COEFF 0.35
75 
76 /* Baseline speed for rolling (m/s) */
77 #define IDEAL_ROLL_SPEED 6.0
78 
79 /* Tux's orientation is updated using a first-order filter; these are the
80    time constants of the filter when tux is on ground and airborne (s) */
81 #define TUX_ORIENTATION_TIME_CONSTANT 0.14
82 #define TUX_ORIENTATION_AIRBORNE_TIME_CONSTANT 0.5
83 
84 /* Max particles generated by turning (particles/s) */
85 #define MAX_TURN_PARTICLES 1500
86 
87 /* Max particles generated by rolling (particles/s) */
88 #define MAX_ROLL_PARTICLES 3000
89 
90 /* Particles generated by braking (particles/s) */
91 #define BRAKE_PARTICLES 4000
92 
93 /* Number of particles is multiplied by ( speed / PARTICLES_SPEED_FACTOR ) */
94 #define PARTICLE_SPEED_FACTOR 40
95 
96 /* Maximum deflection angle for particles (out to side of player)  */
97 #define MAX_PARTICLE_ANGLE 80
98 
99 /* Speed at which maximum deflection occurs */
100 #define MAX_PARTICLE_ANGLE_SPEED 50
101 
102 /* Particle speed = min ( MAX_PARTICLE_SPEED, speed*PARTICLE_SPEED_MULTIPLIER )
103  */
104 #define PARTICLE_SPEED_MULTIPLIER 0.3
105 #define MAX_PARTICLE_SPEED 2
106 
107 /* The compressiblity of the first section of the spring modelling Tux's
108    tush (m) */
109 #define TUX_GLUTE_STAGE_1_COMPRESSIBILITY 0.05
110 
111 /* The spring coefficient of first section of aforementioned tush (N/m) */
112 #define TUX_GLUTE_STAGE_1_SPRING_COEFF 1500
113 
114 /* The damping coefficient of first section of aforementioned tush (N*s/m) */
115 #define TUX_GLUTE_STAGE_1_DAMPING_COEFF 100
116 
117 /* The compressiblity of the second section of the spring modelling Tux's
118    tush (m) */
119 #define TUX_GLUTE_STAGE_2_COMPRESSIBILITY 0.12
120 
121 /* The spring coefficient of second section of aforementioned tush (N/m) */
122 #define TUX_GLUTE_STAGE_2_SPRING_COEFF 3000
123 
124 /* The damping coefficient of second section of aforementioned tush (N*s/m) */
125 #define TUX_GLUTE_STAGE_2_DAMPING_COEFF 500
126 
127 /* The spring coefficient of third section of aforementioned tush (N/m) */
128 #define TUX_GLUTE_STAGE_3_SPRING_COEFF 10000
129 
130 /* The damping coefficient of third section of aforementioned tush (N*s/m) */
131 #define TUX_GLUTE_STAGE_3_DAMPING_COEFF 1000
132 
133 /* The maximum force exertedd by both sections of spring (N) */
134 #define TUX_GLUTE_MAX_SPRING_FORCE 3000
135 
136 /* Maximum distance that player can penetrate into terrain (m) */
137 #define MAX_SURFACE_PENETRATION 0.2
138 
139 /* Density of air @ -3 C (kg/m^3) */
140 #define AIR_DENSITY 1.308
141 
142 /* Diameter of sphere used to model Tux for air resistance calcs (m) */
143 #define TUX_DIAMETER TUX_WIDTH
144 
145 /* Viscosity of air @ -3 C (N*s/m^2) */
146 #define AIR_VISCOSITY 17.00e-6
147 
148 /* Drag coefficient vs. Reynolds number table (from Janna, _Introduction to
149    Fluid Mechanics_, 3rd Edition, PWS-Kent, p 308)
150    Values are for Re in the range of 10^-1 to 10^6, spaced logarithmically
151    (all values are log base 10) */
152 static const double air_log_re[] = { -1,
153 				     0,
154 				     1,
155 				     2,
156 				     3,
157 				     4,
158 				     5,
159 				     6 };
160 static const double air_log_drag_coeff[] = { 2.25,
161 					     1.35,
162 					     0.6,
163 					     0.0,
164 					     -0.35,
165 					     -0.45,
166 					     -0.33,
167 					     -0.9 };
168 
169 /* Minimum time step for ODE solver (s) */
170 #define MIN_TIME_STEP 0.01
171 
172 /* Maximum time step for ODE solver (s) */
173 #define MAX_TIME_STEP 0.10
174 
175 /* Maximum distance of step for ODE solver (m) */
176 #define MAX_STEP_DISTANCE 0.20
177 
178 /* Tolerance on error in Tux's position (m) */
179 #define MAX_POSITION_ERROR 0.005
180 
181 /* Tolerance on error in Tux's velocity (m/s) */
182 #define MAX_VELOCITY_ERROR 0.05
183 
184 /* To smooth out the terrain, the terrain normals are interpolated
185    near the edges of the triangles.  This parameter controls how much
186    smoothing is done, as well as the size of the triangle boundary
187    over which smoothing occurs.  A value of 0 will result in no
188    smoothing, while a value of 0.5 will smooth the normals over the
189    entire surface of the triangle (though the smoothing will increase
190    as the value approaches infinity).  (dimensionless) */
191 #define NORMAL_INTERPOLATION_LIMIT 0.05
192 
193 /* The distance to move Tux up by so that he doesn't sit too low on the
194    ground (m) */
195 #define TUX_Y_CORRECTION_ON_STOMACH 0.33
196 
197 /* The square of the distance that Tux must move before we recompute a
198    collision (m^2) */
199 #define COLLISION_TOLERANCE 0.04
200 
201 /* Duraction of paddling motion (s) */
202 #define PADDLING_DURATION 0.40
203 
204 /* Force applied against ground when paddling (N) */
205 #define MAX_PADDLING_FORCE 122.5
206 
207 /* Ideal paddling friction coefficient */
208 #define IDEAL_PADDLING_FRIC_COEFF 0.35
209 
210 /* G force applied for jump with charge of 0 */
211 #define BASE_JUMP_G_FORCE 1.5
212 
213 /* Maximum G force applied during jump */
214 #define MAX_JUMP_G_FORCE 3
215 
216 /* Magnitude of force before damage is incurred (N) */
217 #define DAMAGE_RESISTANCE ( 4.0 * TUX_MASS * EARTH_GRAV )
218 
219 /* Damage scaling factor (health/(N*s)) */
220 #define DAMAGE_SCALE 9e-8
221 
222 /* Amount to scale tree force by (so that it does more damage than regular
223    collisions) (dimensionless) */
224 #define COLLISION_BASE_DAMAGE 0.02
225 
226 /* Damage incurred by paddling exertion (health/s) */
227 #define PADDLING_DAMAGE 0.02
228 
229 
230 /* Wind velocity (hard-coded for now) */
231 static vector_t wind_vel = { 250.0/3.6, 0, 0 };
232 static scalar_t wind_scale = 0.5;
233 
234 /*
235  * Static variables
236  */
237 
238 /* Friction coefficients (dimensionless) */
239 static scalar_t  fricCoeff[3] = { 0.35, /* ice */
240                                   0.7,  /* rock */
241                                   0.45  /* snow */
242                                 };
243 
244 
245 /* Current time step in ODE solver */
246 static double ode_time_step = -1;
247 
248 
249 /*
250  * Function definitions
251  */
252 
increment_turn_fact(player_data_t * plyr,scalar_t amt)253 void increment_turn_fact( player_data_t *plyr, scalar_t amt )
254 {
255     plyr->control.turn_fact += amt;
256     plyr->control.turn_fact = min( 1.0, max( plyr->control.turn_fact, -1.0 ) );
257 }
258 
get_min_tux_speed()259 scalar_t get_min_tux_speed()
260 {
261     return MIN_TUX_SPEED;
262 }
263 
set_friction_coeff(scalar_t fric[3])264 void set_friction_coeff( scalar_t fric[3] )
265 {
266     fricCoeff[0] = fric[0];
267     fricCoeff[1] = fric[1];
268     fricCoeff[2] = fric[2];
269 }
270 
get_min_y_coord()271 scalar_t get_min_y_coord()
272 {
273     scalar_t courseWidth, courseLength;
274     scalar_t angle;
275     scalar_t minY;
276 
277     get_course_dimensions( &courseWidth, &courseLength );
278     angle = get_course_angle();
279 
280     minY = -courseLength * tan( ANGLES_TO_RADIANS( angle ) );
281     return minY;
282 }
283 
get_indices_for_point(scalar_t x,scalar_t z,int * x0,int * y0,int * x1,int * y1)284 void get_indices_for_point( scalar_t x, scalar_t z,
285 			    int *x0, int *y0, int *x1, int *y1 )
286 {
287     scalar_t courseWidth, courseLength;
288     scalar_t xidx, yidx;
289     int nx, ny;
290 
291     get_course_dimensions( &courseWidth, &courseLength );
292     get_course_divisions( &nx, &ny );
293 
294     xidx = x / courseWidth * ( (scalar_t) nx - 1. );
295     yidx = -z / courseLength * ( (scalar_t) ny - 1. );
296 
297     if (xidx < 0) {
298         xidx = 0;
299     } else if ( xidx > nx-1 ) {
300         xidx = nx-1;
301     }
302 
303     if (yidx < 0) {
304         yidx = 0;
305     } else if ( yidx > ny-1 ) {
306         yidx = ny-1;
307     }
308 
309     /* I found that ceil(3) was quite slow on at least some architectures,
310        so I've replace it with an approximation */
311     *x0 = (int) (xidx);              /* floor(xidx) */
312     *x1 = (int) ( xidx + 0.9999 );   /* ceil(xidx) */
313     *y0 = (int) (yidx);              /* floor(yidx) */
314     *y1 = (int) ( yidx + 0.9999 );   /* ceil(yidx) */
315 
316     if ( *x0 == *x1 ) {
317 	if ( *x1 < nx - 1 )
318 	    (*x1)++;
319 	else
320 	    (*x0)--;
321     }
322 
323     if ( *y0 == *y1 ) {
324 	if ( *y1 < ny - 1 )
325 	    (*y1)++;
326 	else
327 	    (*y0)--;
328     }
329 
330     check_assertion( *x0 >= 0, "invalid x0 index" );
331     check_assertion( *x0 < nx, "invalid x0 index" );
332     check_assertion( *x1 >= 0, "invalid x1 index" );
333     check_assertion( *x1 < nx, "invalid x1 index" );
334     check_assertion( *y0 >= 0, "invalid y0 index" );
335     check_assertion( *y0 < ny, "invalid y0 index" );
336     check_assertion( *y1 >= 0, "invalid y1 index" );
337     check_assertion( *y1 < ny, "invalid y1 index" );
338 }
339 
find_barycentric_coords(scalar_t x,scalar_t z,index2d_t * idx0,index2d_t * idx1,index2d_t * idx2,scalar_t * u,scalar_t * v)340 void find_barycentric_coords( scalar_t x, scalar_t z,
341 			      index2d_t *idx0,
342 			      index2d_t *idx1,
343 			      index2d_t *idx2,
344 			      scalar_t *u, scalar_t *v )
345 {
346     scalar_t xidx, yidx;
347     int nx, ny;
348     int x0, x1, y0, y1;
349     scalar_t dx, ex, dz, ez, qx, qz, invdet; /* to calc. barycentric coords */
350     scalar_t courseWidth, courseLength;
351     scalar_t *elevation;
352 
353     elevation = get_course_elev_data();
354     get_course_dimensions( &courseWidth, &courseLength );
355     get_course_divisions( &nx, &ny );
356 
357     get_indices_for_point( x, z, &x0, &y0, &x1, &y1 );
358 
359     xidx = x / courseWidth * ( (scalar_t) nx - 1. );
360     yidx = -z / courseLength * ( (scalar_t) ny - 1. );
361 
362 
363     /* The terrain is meshed as follows:
364 	      ...
365 	   +-+-+-+-+            x<---+
366 	   |\|/|\|/|                 |
367         ...+-+-+-+-+...              V
368 	   |/|\|/|\|                 y
369 	   +-+-+-+-+
370 	      ...
371 
372        So there are two types of squares: those like this (x0+y0 is even):
373 
374        +-+
375        |/|
376        +-+
377 
378        and those like this (x0+y0 is odd).
379 
380        +-+
381        |\|
382        +-+
383     */
384 
385 #define IDX(x,y) make_index2d(x,y)
386 
387     if ( (x0 + y0) % 2 == 0 ) {
388 	if ( yidx - y0 < xidx - x0 ) {
389 	    *idx0 = IDX(x0, y0);
390 	    *idx1 = IDX(x1, y0);
391 	    *idx2 = IDX(x1, y1);
392 	} else {
393 	    *idx0 = IDX(x1, y1);
394 	    *idx1 = IDX(x0, y1);
395 	    *idx2 = IDX(x0, y0);
396 	}
397     } else {
398 	/* x0 + y0 is odd */
399 	if ( yidx - y0 + xidx - x0 < 1 ) {
400 	    *idx0 = IDX(x0, y0);
401 	    *idx1 = IDX(x1, y0);
402 	    *idx2 = IDX(x0, y1);
403 	} else {
404 	    *idx0 = IDX(x1, y1);
405 	    *idx1 = IDX(x0, y1);
406 	    *idx2 = IDX(x1, y0);
407 	}
408     }
409 #undef IDX
410 
411     dx = idx0->i - idx2->i;
412     dz = idx0->j - idx2->j;
413     ex = idx1->i - idx2->i;
414     ez = idx1->j - idx2->j;
415     qx = xidx - idx2->i;
416     qz = yidx - idx2->j;
417 
418     invdet = 1./(dx * ez - dz * ex);
419     *u = (qx * ez - qz * ex) * invdet;
420     *v = (qz * dx - qx * dz) * invdet;
421 }
422 
find_course_normal(scalar_t x,scalar_t z)423 vector_t find_course_normal( scalar_t x, scalar_t z )
424 {
425     vector_t *course_nmls;
426     vector_t tri_nml;
427     scalar_t xidx, yidx;
428     int nx, ny;
429     int x0, x1, y0, y1;
430     point_t p0, p1, p2;
431     index2d_t idx0, idx1, idx2;
432     vector_t n0, n1, n2;
433     scalar_t course_width, course_length;
434     scalar_t u, v;
435     scalar_t min_bary, interp_factor;
436     vector_t smooth_nml;
437     vector_t interp_nml;
438     scalar_t *elevation;
439 
440     elevation = get_course_elev_data();
441     get_course_dimensions( &course_width, &course_length );
442     get_course_divisions( &nx, &ny );
443     course_nmls = get_course_normals();
444 
445     get_indices_for_point( x, z, &x0, &y0, &x1, &y1 );
446 
447     xidx = x / course_width * ( (scalar_t) nx - 1. );
448     yidx = -z / course_length * ( (scalar_t) ny - 1. );
449 
450     find_barycentric_coords( x, z, &idx0, &idx1, &idx2, &u, &v );
451 
452     n0 = course_nmls[ idx0.i + nx * idx0.j ];
453     n1 = course_nmls[ idx1.i + nx * idx1.j ];
454     n2 = course_nmls[ idx2.i + nx * idx2.j ];
455 
456     p0 = COURSE_VERTEX( idx0.i, idx0.j );
457     p1 = COURSE_VERTEX( idx1.i, idx1.j );
458     p2 = COURSE_VERTEX( idx2.i, idx2.j );
459 
460     smooth_nml = add_vectors(
461 	scale_vector( u, n0 ),
462 	add_vectors( scale_vector( v, n1 ), scale_vector( 1.-u-v, n2 ) ) );
463 
464     tri_nml = cross_product(
465 	subtract_points( p1, p0 ), subtract_points( p2, p0 ) );
466     normalize_vector( &tri_nml );
467 
468     min_bary = min( u, min( v, 1. - u - v ) );
469     interp_factor = min( min_bary / NORMAL_INTERPOLATION_LIMIT, 1.0 );
470 
471     interp_nml = add_vectors(
472 	scale_vector( interp_factor, tri_nml ),
473 	scale_vector( 1.-interp_factor, smooth_nml ) );
474     normalize_vector( &interp_nml );
475 
476     return interp_nml;
477 }
478 
find_y_coord(scalar_t x,scalar_t z)479 scalar_t find_y_coord( scalar_t x, scalar_t z )
480 {
481     scalar_t ycoord;
482     point_t p0, p1, p2;
483     index2d_t idx0, idx1, idx2;
484     scalar_t u, v;
485     int nx, ny;
486     scalar_t *elevation;
487     scalar_t course_width, course_length;
488 
489     /* cache last request */
490     static scalar_t last_x, last_z, last_y;
491     static bool_t cache_full = False;
492 
493     if ( cache_full && last_x == x && last_z == z ) {
494 	return last_y;
495     }
496 
497     get_course_divisions( &nx, &ny );
498     get_course_dimensions( &course_width, &course_length );
499     elevation = get_course_elev_data();
500 
501     find_barycentric_coords( x, z, &idx0, &idx1, &idx2, &u, &v );
502 
503     p0 = COURSE_VERTEX( idx0.i, idx0.j );
504     p1 = COURSE_VERTEX( idx1.i, idx1.j );
505     p2 = COURSE_VERTEX( idx2.i, idx2.j );
506 
507     ycoord = u * p0.y + v * p1.y + ( 1. - u - v ) * p2.y;
508 
509     last_x = x;
510     last_z = z;
511     last_y = ycoord;
512     cache_full = True;
513 
514     return ycoord;
515 }
516 
get_surface_type(scalar_t x,scalar_t z,scalar_t weights[])517 void get_surface_type( scalar_t x, scalar_t z, scalar_t weights[] )
518 {
519     terrain_t *terrain;
520     scalar_t courseWidth, courseLength;
521     int nx, ny;
522     index2d_t idx0, idx1, idx2;
523     scalar_t u, v;
524     int i;
525 
526     find_barycentric_coords( x, z, &idx0, &idx1, &idx2, &u, &v );
527 
528     terrain = get_course_terrain_data();
529     get_course_dimensions( &courseWidth, &courseLength );
530     get_course_divisions( &nx, &ny );
531 
532     for (i=0; i<NumTerrains; i++) {
533 	weights[i] = 0;
534 	if ( terrain[ idx0.i + nx*idx0.j ] == i ) {
535 	    weights[i] += u;
536 	}
537 	if ( terrain[ idx1.i + nx*idx1.j ] == i ) {
538 	    weights[i] += v;
539 	}
540 	if ( terrain[ idx2.i + nx*idx2.j ] == i ) {
541 	    weights[i] += 1.0 - u - v;
542 	}
543     }
544 }
545 
get_local_course_plane(point_t pt)546 plane_t get_local_course_plane( point_t pt )
547 {
548     plane_t plane;
549 
550     pt.y = find_y_coord( pt.x, pt.z );
551 
552     plane.nml = find_course_normal( pt.x, pt.z );
553     plane.d = -( plane.nml.x * pt.x +
554 		 plane.nml.y * pt.y +
555 		 plane.nml.z * pt.z );
556 
557     return plane;
558 }
559 
update_paddling(player_data_t * plyr)560 static void update_paddling( player_data_t *plyr )
561 {
562     if ( plyr->control.is_paddling ) {
563 	if ( g_game.time - plyr->control.paddle_time >= PADDLING_DURATION ) {
564 	    print_debug( DEBUG_CONTROL, "paddling off" );
565 	    plyr->control.is_paddling = False;
566 	}
567     }
568 }
569 
set_tux_pos(player_data_t * plyr,point_t new_pos)570 void set_tux_pos( player_data_t *plyr, point_t new_pos )
571 {
572     char *tuxRoot;
573     scalar_t playWidth, playLength;
574     scalar_t courseWidth, courseLength;
575     scalar_t boundaryWidth;
576     scalar_t disp_y;
577 
578     get_play_dimensions( &playWidth, &playLength );
579     get_course_dimensions( &courseWidth, &courseLength );
580     boundaryWidth = ( courseWidth - playWidth ) / 2;
581 
582 
583     if ( new_pos.x < boundaryWidth ) {
584         new_pos.x = boundaryWidth;
585     } else if ( new_pos.x > courseWidth - boundaryWidth ) {
586         new_pos.x = courseWidth - boundaryWidth;
587     }
588 
589     if ( new_pos.z > 0 ) {
590         new_pos.z = 0;
591     } else if ( -new_pos.z >= playLength ) {
592         new_pos.z = -playLength;
593         set_game_mode( GAME_OVER );
594     }
595 
596     plyr->pos = new_pos;
597 
598     disp_y = new_pos.y + TUX_Y_CORRECTION_ON_STOMACH;
599 
600     tuxRoot = get_tux_root_node();
601     reset_scene_node( tuxRoot );
602     translate_scene_node( tuxRoot,
603 			  make_vector( new_pos.x, disp_y, new_pos.z ) );
604 }
605 
check_tree_collisions(player_data_t * plyr,point_t pos,point_t * tree_loc,scalar_t * tree_diam)606 bool_t check_tree_collisions( player_data_t *plyr, point_t pos,
607 			      point_t *tree_loc, scalar_t *tree_diam )
608 {
609     tree_t *trees;
610     int num_trees, i;
611     scalar_t diam = 0.0;
612     scalar_t height;
613     polyhedron_t ph, ph2;
614     matrixgl_t mat;
615     point_t loc;
616     bool_t hit = False;
617     vector_t distvec;
618     char *tux_root;
619     scalar_t squared_dist;
620     int tree_type;
621 
622     /* These variables are used to cache collision detection results */
623     static bool_t last_collision = False;
624     static point_t last_collision_tree_loc = { -999, -999, -999 };
625     static scalar_t last_collision_tree_diam = 0;
626     static point_t last_collision_pos = { -999, -999, -999 };
627 
628     /* If we haven't moved very much since the last call, we re-use
629        the results of the last call (significant speed-up) */
630     vector_t dist_vec = subtract_points( pos, last_collision_pos );
631     if ( MAG_SQD( dist_vec ) < COLLISION_TOLERANCE ) {
632 	if ( last_collision ) {
633 	    if ( tree_loc != NULL ) {
634 		*tree_loc = last_collision_tree_loc;
635 	    }
636 	    if ( tree_diam != NULL ) {
637 		*tree_diam = last_collision_tree_diam;
638 	    }
639 	    return True;
640 	} else {
641 	    return False;
642 	}
643     }
644 
645     trees = get_tree_locs();
646     num_trees = get_num_trees();
647     tree_type = trees[0].tree_type;
648     ph = get_tree_polyhedron(tree_type);
649 
650     for (i=0; i<num_trees; i++) {
651         diam = trees[i].diam;
652         height = trees[i].height;
653         loc = trees[i].ray.pt;
654 
655         distvec = make_vector( loc.x - pos.x, 0.0, loc.z - pos.z );
656 
657 	/* check distance from tree; .6 is the radius of a bounding
658            sphere around tux (approximate) */
659 	squared_dist = ( diam/2. + 0.6 );
660 	squared_dist *= squared_dist;
661         if ( MAG_SQD( distvec ) > squared_dist ) {
662             continue;
663         }
664 
665 	/* have to look at polyhedron - switch to correct one if necessary */
666 	if ( tree_type != trees[i].tree_type )  {
667 	    tree_type = trees[i].tree_type;
668 	    ph = get_tree_polyhedron(tree_type);
669 	}
670 
671         ph2 = copy_polyhedron( ph );
672 
673         make_scaling_matrix( mat, diam, height, diam );
674         trans_polyhedron( mat, ph2 );
675         make_translation_matrix( mat, loc.x, loc.y, loc.z );
676         trans_polyhedron( mat, ph2 );
677 
678 	tux_root = get_tux_root_node();
679 	reset_scene_node( tux_root );
680 	translate_scene_node( tux_root,
681 			      make_vector( pos.x, pos.y, pos.z ) );
682         hit = collide( tux_root, ph2 );
683 
684         free_polyhedron( ph2 );
685 
686         if ( hit == True ) {
687 	    if ( tree_loc != NULL ) {
688 		*tree_loc = loc;
689 	    }
690 	    if ( tree_diam != NULL ) {
691 		*tree_diam = diam;
692 	    }
693 
694 
695             break;
696         }
697     }
698 
699     last_collision_tree_loc = loc;
700     last_collision_tree_diam = diam;
701     last_collision_pos = pos;
702 
703     if ( hit ) {
704 	last_collision = True;
705 
706 	/* Record collision in player data so that health can be adjusted */
707 	plyr->collision = True;
708     } else {
709 	last_collision = False;
710     }
711 
712     return hit;
713 }
714 
check_item_collection(player_data_t * plyr,point_t pos)715 void check_item_collection( player_data_t *plyr, point_t pos )
716 {
717     item_t *items;
718     int num_items, i;
719     scalar_t diam = 0.0;
720     scalar_t height;
721     point_t loc;
722     vector_t distvec;
723     scalar_t squared_dist;
724     int item_type;
725 
726     /* These variables are used to cache collision detection results */
727     static point_t last_collision_pos = { -999, -999, -999 };
728 
729     /* If we haven't moved very much since the last call, we re-use
730        the results of the last call (significant speed-up) */
731     vector_t dist_vec = subtract_points( pos, last_collision_pos );
732     if ( MAG_SQD( dist_vec ) < COLLISION_TOLERANCE ) {
733 	return ;
734     }
735 
736     items = get_item_locs();
737     num_items = get_num_items();
738     item_type = items[0].item_type;
739 
740     for (i=0; i<num_items; i++) {
741 
742 	if ( items[i].collectable != 1 ) {
743 	    continue;
744 	}
745 
746         diam = items[i].diam;
747         height = items[i].height;
748         loc = items[i].ray.pt;
749 
750         distvec = make_vector( loc.x - pos.x, 0.0, loc.z - pos.z );
751 
752 	/* check distance from tree; .6 is the radius of a bounding
753            sphere around tux (approximate) */
754 	squared_dist = ( diam/2. + 0.6 );
755 	squared_dist *= squared_dist;
756         if ( MAG_SQD( distvec ) > squared_dist ) {
757             continue;
758         }
759 
760 	if ( (pos.y - 0.6 >= loc.y && pos.y - 0.6 <= loc.y + height) ||
761 	     (pos.y + 0.6 >= loc.y && pos.y + 0.6 <= loc.y + height) ||
762 	     (pos.y - 0.6 <= loc.y && pos.y + 0.6 >= loc.y + height) ) {
763 	    /* close enough to hitting the flag */
764 	    /* play a noise? */
765 	    items[i].collectable = 0;
766 	    plyr->herring += 1;
767 	    play_sound( "item_collect", 0 );
768 	}
769 
770     }
771 
772 }
773 
get_compression_depth(terrain_t surf_type)774 scalar_t get_compression_depth( terrain_t surf_type )
775 {
776     switch ( surf_type ) {
777     case Ice:
778 	return 0.03;
779     case Rock:
780 	return 0.01;
781     case Snow:
782 	return 0.11;
783     default:
784 	code_not_reached();
785     }
786     return 0;
787 }
788 
789 /*
790  * Check for tree collisions and adjust position and velocity appropriately.
791  */
adjust_for_tree_collision(player_data_t * plyr,point_t pos,vector_t * vel)792 static void adjust_for_tree_collision( player_data_t *plyr,
793 				       point_t pos, vector_t *vel )
794 {
795     vector_t treeNml;
796     point_t treeLoc;
797     bool_t treeHit;
798     scalar_t speed;
799     scalar_t costheta;
800     scalar_t tree_diam;
801 
802     treeHit = check_tree_collisions( plyr, pos, &treeLoc, &tree_diam );
803     if (treeHit) {
804 	/*
805 	 * Calculate the normal vector to the tree; here we model the tree
806 	 * as a vertical cylinder.
807 	 */
808         treeNml.x = pos.x - treeLoc.x;
809         treeNml.y = 0;
810         treeNml.z = pos.z - treeLoc.z;
811         normalize_vector( &treeNml );
812 
813 	/* Reduce speed by a minimum of 30% */
814         speed = normalize_vector( vel );
815         speed *= 0.7;
816 
817 	/*
818 	 * If Tux is moving into the tree, reduce the speed further,
819 	 * and reflect the velocity vector off of the tree
820 	 */
821         costheta = dot_product( *vel, treeNml );
822         if (costheta < 0 ) {
823 	    /* Reduce speed */
824 
825             speed *= 1 + costheta;
826             speed *= 1 + costheta;
827 
828 	    /* Do the reflection */
829             *vel = add_vectors(
830 		scale_vector( -2. * dot_product( *vel, treeNml ) , treeNml ),
831 		*vel );
832 
833 	    normalize_vector( vel );
834 
835         }
836 
837 	speed = max( speed, get_min_tux_speed() );
838 
839         *vel = scale_vector( speed, *vel );
840     }
841 }
842 
843 /*
844  * Adjusts velocity so that his speed doesn't drop below the minimum
845  * speed
846  */
adjust_velocity(vector_t * vel,point_t pos,plane_t surf_plane,scalar_t dist_from_surface)847 scalar_t adjust_velocity( vector_t *vel, point_t pos,
848 			  plane_t surf_plane, scalar_t dist_from_surface )
849 {
850     vector_t surf_nml;
851     scalar_t speed;
852 
853     surf_nml = surf_plane.nml;
854 
855     speed = normalize_vector( vel );
856 
857     if ( speed < EPS )
858     {
859 	if ( fabs( surf_nml.x ) + fabs( surf_nml.z ) > EPS ) {
860 	    *vel = project_into_plane(
861 		surf_nml, make_vector( 0.0, -1.0, 0.0 ) );
862 	    normalize_vector( vel );
863 	} else {
864 	    *vel = make_vector( 0.0, 0.0, -1.0 );
865 	}
866     }
867 
868     speed = max( get_min_tux_speed(), speed );
869 
870     *vel = scale_vector( speed, *vel );
871     return speed;
872 }
873 
adjust_position(point_t * pos,plane_t surf_plane,scalar_t dist_from_surface)874 void adjust_position( point_t *pos, plane_t surf_plane,
875 		      scalar_t dist_from_surface )
876 {
877 
878     if ( dist_from_surface < -MAX_SURFACE_PENETRATION )
879     {
880 	*pos = move_point( *pos,
881 			   scale_vector( -MAX_SURFACE_PENETRATION -
882 					 dist_from_surface,
883 					 surf_plane.nml ) );
884     }
885     return;
886 }
887 
adjust_tux_zvec_for_roll(player_data_t * plyr,vector_t vel,vector_t zvec)888 static vector_t adjust_tux_zvec_for_roll( player_data_t *plyr,
889 					  vector_t vel, vector_t zvec )
890 {
891     matrixgl_t rot_mat;
892 
893     vel = project_into_plane( zvec, vel );
894 
895     normalize_vector( &vel );
896 
897     if ( plyr->control.is_braking ) {
898 	make_rotation_about_vector_matrix( rot_mat, vel,
899 					   plyr->control.turn_fact *
900 					   BRAKING_ROLL_ANGLE );
901     } else {
902 	make_rotation_about_vector_matrix( rot_mat, vel,
903 					   plyr->control.turn_fact *
904 					   MAX_ROLL_ANGLE );
905     }
906 
907     return transform_vector( rot_mat, zvec );
908 }
909 
adjust_surf_nml_for_roll(player_data_t * plyr,vector_t vel,scalar_t fric_coeff,vector_t nml)910 static vector_t adjust_surf_nml_for_roll( player_data_t *plyr,
911 					  vector_t vel,
912 					  scalar_t fric_coeff,
913 					  vector_t nml )
914 {
915     matrixgl_t rot_mat;
916     scalar_t angle;
917     scalar_t speed;
918     scalar_t roll_angle;
919 
920     speed = normalize_vector( &vel );
921 
922     vel = project_into_plane( nml, vel );
923 
924     normalize_vector( &vel );
925     if ( plyr->control.is_braking ) {
926 	roll_angle = BRAKING_ROLL_ANGLE;
927     } else {
928 	roll_angle = MAX_ROLL_ANGLE;
929     }
930 
931     angle = plyr->control.turn_fact *
932 	roll_angle *
933 	min( 1.0, max(0.0, fric_coeff)/IDEAL_ROLL_FRIC_COEFF )
934 	* min(1.0, max(0.0,speed-MIN_TUX_SPEED)/
935 	      (IDEAL_ROLL_SPEED-MIN_TUX_SPEED));
936 
937     make_rotation_about_vector_matrix( rot_mat, vel, angle );
938 
939     return transform_vector( rot_mat, nml );
940 }
941 
942 
adjust_orientation(player_data_t * plyr,scalar_t dtime,point_t pos,vector_t vel,scalar_t dist_from_surface,vector_t surf_nml)943 void adjust_orientation( player_data_t *plyr, scalar_t dtime,
944 			 point_t pos, vector_t vel,
945 			 scalar_t dist_from_surface, vector_t surf_nml )
946 {
947     vector_t new_x, new_y, new_z;
948     matrixgl_t cob_mat, inv_cob_mat;
949     matrixgl_t rot_mat;
950     quaternion_t new_orient;
951     char* tux_root;
952     scalar_t time_constant;
953     static vector_t minus_z_vec = { 0., 0., -1. };
954     static vector_t y_vec = { 0., 1., 0. };
955 
956     if ( dist_from_surface > 0 ) {
957 	new_y = scale_vector( 1., vel );
958 	normalize_vector( &new_y );
959 	new_z = project_into_plane( new_y, make_vector(0., -1., 0.) );
960 	normalize_vector( &new_z);
961 	new_z = adjust_tux_zvec_for_roll( plyr, vel, new_z );
962     } else {
963 	new_z = scale_vector( -1., surf_nml );
964 	new_z = adjust_tux_zvec_for_roll( plyr, vel, new_z );
965 	new_y = project_into_plane( surf_nml, scale_vector( 1., vel ) );
966 	normalize_vector(&new_y);
967     }
968 
969     new_x = cross_product( new_y, new_z );
970 
971     make_change_of_basis_matrix( cob_mat, inv_cob_mat, new_x, new_y, new_z );
972 
973     new_orient = make_quaternion_from_matrix( cob_mat );
974 
975     if ( !plyr->orientation_initialized ) {
976 	plyr->orientation_initialized = True;
977 	plyr->orientation = new_orient;
978     }
979 
980     time_constant = dist_from_surface > 0
981 	? TUX_ORIENTATION_AIRBORNE_TIME_CONSTANT
982 	: TUX_ORIENTATION_TIME_CONSTANT;
983 
984     plyr->orientation = interpolate_quaternions(
985 	plyr->orientation, new_orient,
986 	min( dtime / time_constant, 1.0 ) );
987 
988     plyr->plane_nml = rotate_vector( plyr->orientation, minus_z_vec );
989     plyr->direction = rotate_vector( plyr->orientation, y_vec );
990 
991     make_matrix_from_quaternion( cob_mat, plyr->orientation );
992 
993 
994     /* Trick rotations */
995     new_y = make_vector( cob_mat[1][0], cob_mat[1][1], cob_mat[1][2] );
996     make_rotation_about_vector_matrix( rot_mat, new_y,
997 				       ( plyr->control.barrel_roll_factor * 360 ) );
998     multiply_matrices( cob_mat, rot_mat, cob_mat );
999     new_x = make_vector( cob_mat[0][0], cob_mat[0][1], cob_mat[0][2] );
1000     make_rotation_about_vector_matrix( rot_mat, new_x,
1001 				       plyr->control.flip_factor * 360 );
1002     multiply_matrices( cob_mat, rot_mat, cob_mat );
1003 
1004 
1005 
1006     transpose_matrix( cob_mat, inv_cob_mat );
1007 
1008     tux_root = get_tux_root_node();
1009     transform_scene_node( tux_root, cob_mat, inv_cob_mat );
1010 }
1011 
adjust_particle_count(scalar_t particles)1012 scalar_t adjust_particle_count( scalar_t particles )
1013 {
1014     if ( particles < 1 ) {
1015 	if ( ( (scalar_t) rand() ) / RAND_MAX < particles ) {
1016 	    return 1.0;
1017 	} else {
1018 	    return 0.0;
1019 	}
1020     } else {
1021 	return particles;
1022     }
1023 }
1024 
generate_particles(player_data_t * plyr,scalar_t dtime,point_t pos,scalar_t speed)1025 void generate_particles( player_data_t *plyr, scalar_t dtime,
1026 			 point_t pos, scalar_t speed )
1027 {
1028     point_t left_part_pt, right_part_pt;
1029     scalar_t brake_particles;
1030     scalar_t turn_particles;
1031     scalar_t roll_particles;
1032     scalar_t surf_weights[NumTerrains];
1033     scalar_t surf_y;
1034     scalar_t left_particles, right_particles;
1035     vector_t left_part_vel, right_part_vel;
1036     matrixgl_t rot_mat;
1037     vector_t xvec;
1038 
1039     get_surface_type( pos.x, pos.z, surf_weights );
1040     surf_y = find_y_coord( pos.x, pos.z );
1041 
1042     if ( surf_weights[Snow] > 0.5 && pos.y < surf_y ) {
1043 	xvec = cross_product( plyr->direction, plyr->plane_nml );
1044 
1045         right_part_pt = left_part_pt = pos;
1046 
1047 	right_part_pt = move_point(
1048 	    right_part_pt,
1049 	    scale_vector( TUX_WIDTH/2.0, xvec ) );
1050 
1051 	left_part_pt = move_point(
1052 	    left_part_pt,
1053 	    scale_vector( -TUX_WIDTH/2.0, xvec ) );
1054 
1055         right_part_pt.y = left_part_pt.y  = surf_y;
1056 
1057 	brake_particles = dtime *
1058 	    BRAKE_PARTICLES * ( plyr->control.is_braking ? 1.0 : 0.0 )
1059 	    * min( speed / PARTICLE_SPEED_FACTOR, 1.0 );
1060 	turn_particles = dtime * MAX_TURN_PARTICLES
1061 	    * min( speed / PARTICLE_SPEED_FACTOR, 1.0 );
1062 	roll_particles = dtime * MAX_ROLL_PARTICLES
1063 	    * min( speed / PARTICLE_SPEED_FACTOR, 1.0 );
1064 
1065 	left_particles = turn_particles *
1066 	    fabs( min(plyr->control.turn_fact, 0.) ) +
1067 	    brake_particles +
1068 	    roll_particles * fabs( min(plyr->control.turn_animation, 0.) );
1069 
1070 	right_particles = turn_particles *
1071 	    fabs( max(plyr->control.turn_fact, 0.) ) +
1072 	    brake_particles +
1073 	    roll_particles * fabs( max(plyr->control.turn_animation, 0.) );
1074 
1075 	left_particles = adjust_particle_count( left_particles );
1076 	right_particles = adjust_particle_count( right_particles );
1077 
1078 	/* Create particle velocitites */
1079 	make_rotation_about_vector_matrix(
1080 	    rot_mat, plyr->direction,
1081 	    max( -MAX_PARTICLE_ANGLE,
1082 		 -MAX_PARTICLE_ANGLE * speed / MAX_PARTICLE_ANGLE_SPEED ) );
1083 	left_part_vel = transform_vector( rot_mat, plyr->plane_nml );
1084 	left_part_vel = scale_vector( min( MAX_PARTICLE_SPEED,
1085 					   speed * PARTICLE_SPEED_MULTIPLIER ),
1086 				      left_part_vel );
1087 
1088 	make_rotation_about_vector_matrix(
1089 	    rot_mat, plyr->direction,
1090 	    min( MAX_PARTICLE_ANGLE,
1091 		 MAX_PARTICLE_ANGLE * speed / MAX_PARTICLE_ANGLE_SPEED ) );
1092 	right_part_vel = transform_vector( rot_mat, plyr->plane_nml );
1093 	right_part_vel = scale_vector( min( MAX_PARTICLE_SPEED,
1094 					    speed * PARTICLE_SPEED_MULTIPLIER ),
1095 				       right_part_vel );
1096 
1097         create_new_particles( left_part_pt, left_part_vel,
1098 			      (int)left_particles );
1099         create_new_particles( right_part_pt, right_part_vel,
1100 			      (int)right_particles );
1101     }
1102 }
1103 
1104 
1105 /*
1106  * Calculate the magnitude of force due to air resistance (wind)
1107  */
calc_wind_force(vector_t player_vel)1108 vector_t calc_wind_force( vector_t player_vel )
1109 {
1110     vector_t total_vel;
1111     scalar_t wind_speed;
1112     scalar_t re;         /* Reynolds number */
1113     scalar_t df_mag;     /* magnitude of drag force */
1114     scalar_t drag_coeff; /* drag coefficient */
1115     int table_size;      /* size of drag coeff table */
1116 
1117     static scalar_t last_time_called = -1;
1118 
1119     total_vel = scale_vector( -1, player_vel );
1120 
1121     if ( g_game.race.windy ) {
1122 	/* adjust wind_scale with a random walk */
1123 	if ( last_time_called != g_game.time ) {
1124 	    wind_scale = wind_scale +
1125 		(rand()/(scalar_t)RAND_MAX-0.50) * 0.15;
1126 	    wind_scale = min( 1.0, max( 0.0, wind_scale ) );
1127 	}
1128 
1129 	total_vel = add_vectors( total_vel,
1130 				 scale_vector( wind_scale, wind_vel ) );
1131     }
1132 
1133     wind_speed = normalize_vector( &total_vel );
1134 
1135     re = AIR_DENSITY * wind_speed * TUX_DIAMETER / AIR_VISCOSITY;
1136 
1137     table_size = sizeof(air_log_drag_coeff) / sizeof(air_log_drag_coeff[0]);
1138 
1139     drag_coeff = pow( 10.0,
1140 		      lin_interp( air_log_re, air_log_drag_coeff,
1141 				  log10(re), table_size ) );
1142 
1143     df_mag = 0.5 * drag_coeff * AIR_DENSITY * ( wind_speed * wind_speed )
1144 	* ( M_PI * ( TUX_DIAMETER * TUX_DIAMETER ) * 0.25 );
1145 
1146     check_assertion( df_mag > 0, "Negative wind force" );
1147 
1148     last_time_called = g_game.time;
1149 
1150     return scale_vector( df_mag, total_vel );
1151 }
1152 
calc_spring_force(scalar_t compression,vector_t vel,vector_t surf_nml,vector_t * unclamped_f)1153 static vector_t calc_spring_force( scalar_t compression, vector_t vel,
1154 				   vector_t surf_nml, vector_t *unclamped_f )
1155 {
1156     scalar_t spring_vel; /* velocity perp. to surface (for damping) */
1157     scalar_t spring_f_mag; /* magnitude of force */
1158 
1159    check_assertion( compression >= 0,
1160 		    "spring can't have negative compression" );
1161 
1162     spring_vel = dot_product( vel, surf_nml );
1163 
1164     /* Stage 1 */
1165     spring_f_mag =
1166 	min( compression, TUX_GLUTE_STAGE_1_COMPRESSIBILITY )
1167 	* TUX_GLUTE_STAGE_1_SPRING_COEFF;
1168 
1169     /* Stage 2 */
1170     spring_f_mag +=
1171 	max( 0, min( compression - TUX_GLUTE_STAGE_1_COMPRESSIBILITY,
1172 		     TUX_GLUTE_STAGE_2_COMPRESSIBILITY ) )
1173 	* TUX_GLUTE_STAGE_2_SPRING_COEFF;
1174 
1175     /* Stage 3 */
1176     spring_f_mag +=
1177 	max( 0., compression - TUX_GLUTE_STAGE_2_COMPRESSIBILITY -
1178 	     TUX_GLUTE_STAGE_1_COMPRESSIBILITY )
1179 	* TUX_GLUTE_STAGE_3_SPRING_COEFF;
1180 
1181     /* Damping */
1182     spring_f_mag -= spring_vel * (
1183 	compression <= TUX_GLUTE_STAGE_1_COMPRESSIBILITY
1184 	? TUX_GLUTE_STAGE_1_SPRING_COEFF :
1185 	( compression <= TUX_GLUTE_STAGE_2_COMPRESSIBILITY
1186 	  ? TUX_GLUTE_STAGE_2_DAMPING_COEFF
1187 	  : TUX_GLUTE_STAGE_3_DAMPING_COEFF ) );
1188 
1189     /* Clamp to >= 0.0 */
1190     spring_f_mag = max( spring_f_mag, 0.0 );
1191 
1192     if ( unclamped_f != NULL ) {
1193 	*unclamped_f = scale_vector( spring_f_mag, surf_nml );
1194     }
1195 
1196     /* Clamp to <= TUX_GLUTE_MAX_SPRING_FORCE */
1197     spring_f_mag = min( spring_f_mag, TUX_GLUTE_MAX_SPRING_FORCE );
1198 
1199     return scale_vector( spring_f_mag, surf_nml );
1200 }
1201 
1202 
calc_net_force(player_data_t * plyr,point_t pos,vector_t vel)1203 static vector_t calc_net_force( player_data_t *plyr, point_t pos,
1204 				vector_t vel )
1205 {
1206     vector_t nml_f;      /* normal force */
1207     vector_t unclamped_nml_f; /* unclamped normal force (for damage calc) */
1208     vector_t fric_f;     /* frictional force */
1209     scalar_t fric_f_mag; /* frictional force magnitude */
1210     vector_t fric_dir;   /* direction of frictional force */
1211     vector_t grav_f;     /* gravitational force */
1212     vector_t air_f;      /* air resistance force */
1213     vector_t brake_f;    /* braking force */
1214     vector_t paddling_f; /* paddling force */
1215     vector_t net_force;  /* the net force (sum of all other forces) */
1216     scalar_t comp_depth; /* depth to which the terrain can be compressed */
1217     scalar_t speed;      /* speed (m/s) */
1218     vector_t orig_surf_nml; /* normal to terrain at current position */
1219     vector_t surf_nml;   /* normal to terrain w/ roll effect */
1220     scalar_t glute_compression; /* amt that Tux's tush has been compressed */
1221     scalar_t steer_angle; /* Angle to rotate fricitonal force for turning */
1222     matrixgl_t fric_rot_mat; /* Matrix to rotate frictional force */
1223     vector_t jump_f;
1224     plane_t surf_plane;
1225     scalar_t dist_from_surface;
1226     scalar_t surf_weights[NumTerrains];
1227     scalar_t surf_fric_coeff;
1228     int i;
1229 
1230     get_surface_type( pos.x, pos.z, surf_weights );
1231     surf_plane = get_local_course_plane( pos );
1232     orig_surf_nml = surf_plane.nml;
1233 
1234     surf_fric_coeff = 0;
1235     for (i=0; i<NumTerrains; i++) {
1236 	surf_fric_coeff += surf_weights[i] * fricCoeff[i];
1237     }
1238     surf_nml = adjust_surf_nml_for_roll( plyr, vel, surf_fric_coeff,
1239 					 orig_surf_nml );
1240 
1241     comp_depth = 0;
1242     for (i=0; i<NumTerrains; i++) {
1243 	comp_depth += surf_weights[i] * get_compression_depth( (terrain_t)i );
1244     }
1245 
1246 
1247     grav_f = make_vector( 0, -EARTH_GRAV * TUX_MASS, 0 );
1248 
1249     dist_from_surface = distance_to_plane( surf_plane, pos );
1250 
1251     if ( dist_from_surface <= 0 ) {
1252 	plyr->airborne = False;
1253     } else {
1254 	plyr->airborne = True;
1255     }
1256 
1257     /*
1258      * Calculate normal force
1259      */
1260     nml_f = make_vector( 0., 0., 0. );
1261     if ( dist_from_surface <= -comp_depth ) {
1262 	/*
1263 	 * Tux ended up below the surface.
1264 	 * Calculate the spring force exterted by his rear end. :-)
1265 	 */
1266 	glute_compression = -dist_from_surface - comp_depth;
1267 	check_assertion( glute_compression >= 0,
1268 			 "unexpected negative compression" );
1269 
1270 	nml_f = calc_spring_force( glute_compression, vel, surf_nml,
1271 				   &unclamped_nml_f );
1272     }
1273 
1274     /* Check if player is trying to jump */
1275     if ( plyr->control.begin_jump == True ) {
1276 	plyr->control.begin_jump = False;
1277 	if ( dist_from_surface <= 0 ) {
1278 	    plyr->control.jumping = True;
1279 	    plyr->control.jump_start_time = g_game.time;
1280 	} else {
1281 	    plyr->control.jumping = False;
1282 	}
1283     }
1284 
1285 
1286     /* Apply jump force in up direction for JUMP_FORCE_DURATION */
1287     if ( ( plyr->control.jumping ) &&
1288 	 ( g_game.time - plyr->control.jump_start_time <
1289 	   JUMP_FORCE_DURATION ) )
1290     {
1291 	jump_f = make_vector(
1292 	    0,
1293 	    BASE_JUMP_G_FORCE * TUX_MASS * EARTH_GRAV +
1294 	    plyr->control.jump_amt *
1295 	    (MAX_JUMP_G_FORCE-BASE_JUMP_G_FORCE) * TUX_MASS * EARTH_GRAV,
1296 	    0 );
1297 
1298     } else {
1299 	jump_f = make_vector( 0, 0, 0 );
1300 	plyr->control.jumping = False;
1301     }
1302 
1303     /* Use the unclamped normal force for damage calculation purposes */
1304     plyr->normal_force = unclamped_nml_f;
1305 
1306     /*
1307      * Calculate frictional force
1308      */
1309     fric_dir = vel;
1310     speed = normalize_vector( &fric_dir );
1311     fric_dir = scale_vector( -1.0, fric_dir );
1312 
1313     if ( dist_from_surface < 0 && speed > MIN_FRICTION_SPEED ) {
1314 	vector_t tmp_nml_f = nml_f;
1315 
1316 	fric_f_mag = normalize_vector( &tmp_nml_f ) * surf_fric_coeff;
1317 
1318 	fric_f_mag = min( MAX_FRICTIONAL_FORCE, fric_f_mag );
1319 
1320 	fric_f = scale_vector( fric_f_mag, fric_dir );
1321 
1322 
1323 	/*
1324 	 * Adjust friction for steering
1325 	 */
1326 	steer_angle = plyr->control.turn_fact * MAX_TURN_ANGLE;
1327 
1328 	if ( fabs( fric_f_mag * sin( ANGLES_TO_RADIANS( steer_angle ) ) ) >
1329 	     MAX_TURN_PERPENDICULAR_FORCE )
1330 	{
1331 	    check_assertion( fabs( plyr->control.turn_fact ) > 0,
1332 			     "steer angle is non-zero when player is not "
1333 			     "turning" );
1334 	    steer_angle = RADIANS_TO_ANGLES(
1335 		asin( MAX_TURN_PERPENDICULAR_FORCE / fric_f_mag ) ) *
1336 		plyr->control.turn_fact / fabs( plyr->control.turn_fact );
1337 	}
1338 	make_rotation_about_vector_matrix( fric_rot_mat, orig_surf_nml,
1339 					   steer_angle );
1340 	fric_f = transform_vector( fric_rot_mat, fric_f );
1341 	fric_f = scale_vector( 1.0 + MAX_TURN_PENALTY, fric_f );
1342 
1343 
1344 	/*
1345 	 * Calculate braking force
1346 	 */
1347 	if ( speed > get_min_tux_speed() && plyr->control.is_braking ) {
1348 	    brake_f = scale_vector(
1349 		surf_fric_coeff * BRAKE_FORCE, fric_dir );
1350 	} else {
1351 	    brake_f = make_vector( 0., 0., 0. );
1352 	}
1353 
1354     } else {
1355 	fric_f = brake_f = make_vector( 0., 0., 0. );
1356     }
1357 
1358     /*
1359      * Calculate air resistance
1360      */
1361     air_f = calc_wind_force( vel );
1362 
1363 
1364     /*
1365      * Calculate force from paddling
1366      */
1367     update_paddling( plyr );
1368     if ( plyr->control.is_paddling ) {
1369 	if ( plyr->airborne ) {
1370 	    paddling_f = make_vector( 0, 0, -TUX_MASS * EARTH_GRAV / 4.0 );
1371 	    paddling_f = rotate_vector( plyr->orientation, paddling_f );
1372 	} else {
1373 	    paddling_f = scale_vector(
1374 		-1 * min( MAX_PADDLING_FORCE,
1375 			  MAX_PADDLING_FORCE *
1376 			  ( MAX_PADDLING_SPEED - speed ) / MAX_PADDLING_SPEED *
1377 			  min(1.0,
1378 			      surf_fric_coeff/IDEAL_PADDLING_FRIC_COEFF)) ,
1379 		fric_dir );
1380 	}
1381     } else {
1382 	paddling_f = make_vector( 0., 0., 0. );
1383     }
1384 
1385 
1386     /*
1387      * Add all the forces
1388      */
1389     net_force = add_vectors( jump_f, add_vectors( grav_f,
1390 		  add_vectors( nml_f, add_vectors( fric_f,
1391 		   add_vectors( air_f, add_vectors( brake_f,
1392 		     paddling_f ))))));
1393 
1394     return net_force;
1395 }
1396 
adjust_time_step_size(scalar_t h,vector_t vel)1397 static scalar_t adjust_time_step_size( scalar_t h, vector_t vel )
1398 {
1399     scalar_t speed;
1400 
1401     speed = normalize_vector( &vel );
1402 
1403     h = max( h, MIN_TIME_STEP );
1404     h = min( h, MAX_STEP_DISTANCE / speed );
1405     h = min( h, MAX_TIME_STEP );
1406 
1407     return h;
1408 }
1409 
1410 /*
1411  * Solves the system of ordinary differential equations governing Tux's
1412  * movement.  Based on Matlab's ode23.m, (c) The MathWorks, Inc.
1413  */
solve_ode_system(player_data_t * plyr,scalar_t dtime)1414 void solve_ode_system( player_data_t *plyr, scalar_t dtime )
1415 {
1416     double t0, t, tfinal;
1417     ode_data_t *x, *y, *z, *vx, *vy, *vz; /* store estimates of derivs */
1418     ode_solver_t solver;
1419     double h;
1420     bool_t done = False;
1421     bool_t failed = False;
1422     point_t new_pos;
1423     vector_t new_vel, tmp_vel;
1424     scalar_t speed;
1425     vector_t new_f;
1426     point_t saved_pos;
1427     vector_t saved_vel, saved_f;
1428     double pos_err[3], vel_err[3], tot_pos_err, tot_vel_err;
1429     double err=0, tol=0;
1430     int i;
1431 
1432     /* Select a solver */
1433     switch ( getparam_ode_solver() ) {
1434     case EULER:
1435 	solver = new_euler_solver();
1436 	break;
1437     case ODE23:
1438 	solver = new_ode23_solver();
1439 	break;
1440     case ODE45:
1441 	solver = new_ode45_solver();
1442 	break;
1443     default:
1444 	setparam_ode_solver( ODE23 );
1445 	solver = new_ode23_solver();
1446     }
1447 
1448     /* Select an initial time step */
1449     h = ode_time_step;
1450 
1451     if ( h < 0 || solver.estimate_error == NULL ) {
1452 	h = adjust_time_step_size( dtime, plyr->vel );
1453     }
1454 
1455     t0 = 0;
1456     tfinal = dtime;
1457 
1458     t = t0;
1459 
1460     /* Create variables to store derivative estimates & other data */
1461     x  = solver.new_ode_data();
1462     y  = solver.new_ode_data();
1463     z  = solver.new_ode_data();
1464     vx = solver.new_ode_data();
1465     vy = solver.new_ode_data();
1466     vz = solver.new_ode_data();
1467 
1468     /* initialize state */
1469     new_pos = plyr->pos;
1470     new_vel = plyr->vel;
1471     new_f   = plyr->net_force;
1472 
1473     /* loop until we've integrated from t0 to tfinal */
1474     while (!done) {
1475 
1476 	if ( t >= tfinal ) {
1477 	    print_warning( CRITICAL_WARNING,
1478 			   "t >= tfinal in solve_ode_system()" );
1479 	    break;
1480 	}
1481 
1482 	/* extend h by up to 10% to reach tfinal */
1483 	if ( 1.1 * h > tfinal - t ) {
1484 	    h = tfinal-t;
1485 	    check_assertion( h >= 0., "integrated past tfinal" );
1486 	    done = True;
1487 	}
1488 
1489         print_debug( DEBUG_ODE, "h: %g", h );
1490 
1491 	saved_pos = new_pos;
1492 	saved_vel = new_vel;
1493 	saved_f = new_f;
1494 
1495 	/* Loop until error is acceptable */
1496 	failed = False;
1497 
1498 	for (;;) {
1499 
1500 	    /* Store initial conditions */
1501 	    solver.init_ode_data( x, new_pos.x, h );
1502 	    solver.init_ode_data( y, new_pos.y, h );
1503 	    solver.init_ode_data( z, new_pos.z, h );
1504 	    solver.init_ode_data( vx, new_vel.x, h );
1505 	    solver.init_ode_data( vy, new_vel.y, h );
1506 	    solver.init_ode_data( vz, new_vel.z, h );
1507 
1508 
1509 	    /* We assume that the first estimate in all ODE solvers is equal
1510 	       to the final value of the last iteration */
1511 	    solver.update_estimate( x, 0, new_vel.x );
1512 	    solver.update_estimate( y, 0, new_vel.y );
1513 	    solver.update_estimate( z, 0, new_vel.z );
1514 	    solver.update_estimate( vx, 0, new_f.x / TUX_MASS );
1515 	    solver.update_estimate( vy, 0, new_f.y / TUX_MASS );
1516 	    solver.update_estimate( vz, 0, new_f.z / TUX_MASS );
1517 
1518 	    /* Update remaining estimates */
1519 	    for ( i=1; i < solver.num_estimates(); i++ ) {
1520 		new_pos.x = solver.next_val( x, i );
1521 		new_pos.y = solver.next_val( y, i );
1522 		new_pos.z = solver.next_val( z, i );
1523 		new_vel.x = solver.next_val( vx, i );
1524 		new_vel.y = solver.next_val( vy, i );
1525 		new_vel.z = solver.next_val( vz, i );
1526 
1527 		solver.update_estimate( x, i, new_vel.x );
1528 		solver.update_estimate( y, i, new_vel.y );
1529 		solver.update_estimate( z, i, new_vel.z );
1530 
1531 		new_f = calc_net_force( plyr, new_pos, new_vel );
1532 
1533 		solver.update_estimate( vx, i, new_f.x / TUX_MASS );
1534 		solver.update_estimate( vy, i, new_f.y / TUX_MASS );
1535 		solver.update_estimate( vz, i, new_f.z / TUX_MASS );
1536 	    }
1537 
1538 	    /* Get final values */
1539 	    new_pos.x = solver.final_estimate( x );
1540 	    new_pos.y = solver.final_estimate( y );
1541 	    new_pos.z = solver.final_estimate( z );
1542 
1543 	    new_vel.x = solver.final_estimate( vx );
1544 	    new_vel.y = solver.final_estimate( vy );
1545 	    new_vel.z = solver.final_estimate( vz );
1546 
1547 	    /* If the current solver can provide error estimates, update h
1548 	       based on the error, and re-evaluate this step if the error is
1549 	       too large  */
1550 	    if ( solver.estimate_error != NULL ) {
1551 
1552 		/* Calculate the error */
1553 		pos_err[0] = solver.estimate_error( x );
1554 		pos_err[1] = solver.estimate_error( y );
1555 		pos_err[2] = solver.estimate_error( z );
1556 
1557 		vel_err[0] = solver.estimate_error( vx );
1558 		vel_err[1] = solver.estimate_error( vy );
1559 		vel_err[2] = solver.estimate_error( vz );
1560 
1561 		tot_pos_err = 0.;
1562 		tot_vel_err = 0.;
1563 		for ( i=0; i<3; i++ ) {
1564 		    pos_err[i] *= pos_err[i];
1565 		    tot_pos_err += pos_err[i];
1566 		    vel_err[i] *= vel_err[i];
1567 		    tot_vel_err += vel_err[i];
1568 		}
1569 		tot_pos_err = sqrt( tot_pos_err );
1570 		tot_vel_err = sqrt( tot_vel_err );
1571 
1572                 print_debug( DEBUG_ODE, "pos_err: %g, vel_err: %g",
1573 			     tot_pos_err, tot_vel_err );
1574 
1575 		if ( tot_pos_err / MAX_POSITION_ERROR >
1576 		     tot_vel_err / MAX_VELOCITY_ERROR )
1577 		{
1578 		    err = tot_pos_err;
1579 		    tol = MAX_POSITION_ERROR;
1580 		} else {
1581 		    err = tot_vel_err;
1582 		    tol = MAX_VELOCITY_ERROR;
1583 		}
1584 
1585 		/* Update h based on error */
1586 		if (  err > tol  && h > MIN_TIME_STEP + EPS  )
1587 		{
1588 		    done = False;
1589 		    if ( !failed ) {
1590 			failed = True;
1591 			h *=  max( 0.5, 0.8 *
1592 				   pow( tol/err,
1593 					solver.time_step_exponent() ) );
1594 		    } else {
1595 			h *= 0.5;
1596 		    }
1597 
1598 		    h = adjust_time_step_size( h, saved_vel );
1599 
1600 		    new_pos = saved_pos;
1601 		    new_vel = saved_vel;
1602 		    new_f = saved_f;
1603 
1604 		} else {
1605 		    /* Error is acceptable or h is at its minimum; stop */
1606 		    break;
1607 		}
1608 
1609 	    } else {
1610 		/* Current solver doesn't provide error estimates; stop */
1611 		break;
1612 	    }
1613 	}
1614 
1615 	/* Update time */
1616 	t = t + h;
1617 
1618 	tmp_vel = new_vel;
1619 	speed = normalize_vector( &tmp_vel );
1620 
1621 	/* only generate particles if we're drawing them */
1622 	if ( getparam_draw_particles() ) {
1623 	    generate_particles( plyr, h, new_pos, speed );
1624 	}
1625 
1626 	/* Calculate the final net force */
1627 	new_f = calc_net_force( plyr, new_pos, new_vel );
1628 
1629 	/* If no failures, compute a new h */
1630 	if ( !failed && solver.estimate_error != NULL ) {
1631 	    double temp = 1.25 * pow(err / tol, solver.time_step_exponent());
1632 	    if ( temp > 0.2 ) {
1633 		h = h / temp;
1634 	    } else {
1635 		h = 5.0 * h;
1636 	    }
1637 	}
1638 
1639 	/* Clamp h to constraints */
1640 	h = adjust_time_step_size( h, new_vel );
1641 
1642 	/* Important: to make trees "solid", we must manipulate the
1643 	   velocity here; if we don't and Tux is moving very quickly,
1644 	   he can pass through trees */
1645 	adjust_for_tree_collision( plyr, new_pos, &new_vel );
1646 
1647 	/* Try to collect items here */
1648 	check_item_collection( plyr, new_pos );
1649 
1650     }
1651 
1652     /* Save time step for next time */
1653     ode_time_step = h;
1654 
1655     plyr->vel = new_vel;
1656     plyr->pos = new_pos;
1657     plyr->net_force = new_f;
1658 
1659     free( x );
1660     free( y );
1661     free( z );
1662     free( vx );
1663     free( vy );
1664     free( vz );
1665 }
1666 
1667 
1668 /*
1669  * Updates Tux's position taking into account gravity, friction, tree
1670  * collisions, etc.  This is the main physics function.
1671  */
update_player_pos(player_data_t * plyr,scalar_t dtime)1672 void update_player_pos( player_data_t *plyr, scalar_t dtime )
1673 {
1674     vector_t surf_nml;   /* normal vector of terrain */
1675     vector_t tmp_vel;
1676     scalar_t speed;
1677     scalar_t paddling_factor;
1678     vector_t local_force;
1679     scalar_t flap_factor;
1680     plane_t surf_plane;
1681     scalar_t dist_from_surface;
1682 
1683     if ( dtime > 2. * EPS ) {
1684 	solve_ode_system( plyr, dtime );
1685     }
1686 
1687     tmp_vel = plyr->vel;
1688 
1689     /*
1690      * Set position, orientation, generate particles
1691      */
1692     surf_plane = get_local_course_plane( plyr->pos );
1693     surf_nml = surf_plane.nml;
1694     dist_from_surface = distance_to_plane( surf_plane, plyr->pos );
1695     adjust_velocity( &plyr->vel, plyr->pos, surf_plane,
1696 		     dist_from_surface );
1697     adjust_position( &plyr->pos, surf_plane, dist_from_surface );
1698 
1699     speed = normalize_vector( &tmp_vel );
1700 
1701     set_tux_pos( plyr, plyr->pos );
1702     adjust_orientation( plyr, dtime, plyr->pos, plyr->vel,
1703 			dist_from_surface, surf_nml );
1704 
1705     flap_factor = 0;
1706 
1707     if ( plyr->control.is_paddling ) {
1708 	scalar_t factor;
1709 	factor = (g_game.time - plyr->control.paddle_time) / PADDLING_DURATION;
1710 	if ( plyr->airborne ) {
1711 	    paddling_factor = 0;
1712 	    flap_factor = factor;
1713 	} else {
1714 	    paddling_factor = factor;
1715 	    flap_factor = 0;
1716 	}
1717     } else {
1718 	paddling_factor = 0.0;
1719     }
1720 
1721     /* calculate force in Tux's local coordinate system */
1722     local_force = rotate_vector( quaternion_conjugate( plyr->orientation ),
1723 				 plyr->net_force );
1724 
1725     if (plyr->control.jumping) {
1726 	flap_factor = (g_game.time - plyr->control.jump_start_time) /
1727 	    JUMP_FORCE_DURATION;
1728     }
1729 
1730     adjust_tux_joints( plyr->control.turn_animation, plyr->control.is_braking,
1731 		       paddling_factor, speed, local_force, flap_factor );
1732 }
1733 
1734 
init_physical_simulation()1735 void init_physical_simulation()
1736 {
1737     vector_t nml;
1738     matrixgl_t rotMat;
1739     scalar_t ycoord;
1740     player_data_t *plyr;
1741     vector_t init_vel;
1742     vector_t init_f;
1743     int i;
1744 
1745     for ( i=0; i<g_game.num_players; i++ ) {
1746 	plyr = get_player_data( i );
1747 
1748 	ycoord = find_y_coord( plyr->pos.x, plyr->pos.z );
1749 	nml = find_course_normal(plyr->pos.x, plyr->pos.z );
1750 	make_rotation_matrix( rotMat, -90., 'x' );
1751 	init_vel = transform_vector( rotMat, nml );
1752 	init_vel = scale_vector( INIT_TUX_SPEED, init_vel );
1753 	init_f = make_vector( 0., 0., 0. );
1754 
1755 	plyr->pos.y = ycoord;
1756 	plyr->vel = init_vel;
1757 	plyr->net_force = init_f;
1758 	plyr->control.turn_fact = 0.0;
1759 	plyr->control.turn_animation = 0.0;
1760 	plyr->control.barrel_roll_factor = 0.0;
1761 	plyr->control.flip_factor = 0.0;
1762 	plyr->control.is_braking = False;
1763 	plyr->orientation_initialized = False;
1764 	plyr->plane_nml = nml;
1765 	plyr->direction = init_vel;
1766 	plyr->normal_force = make_vector(0,0,0);
1767 	plyr->airborne = False;
1768 	plyr->collision = False;
1769 	plyr->control.jump_amt = 0;
1770 	plyr->control.is_paddling = False;
1771 	plyr->control.jumping = False;
1772 	plyr->control.jump_charging = False;
1773 	plyr->control.barrel_roll_left = False;
1774 	plyr->control.barrel_roll_right = False;
1775 	plyr->control.barrel_roll_factor = 0;
1776 	plyr->control.front_flip = False;
1777 	plyr->control.back_flip = False;
1778     }
1779 
1780     ode_time_step = -1;
1781 }
1782 
1783 
1784