Back to index

extremetuxracer  0.5beta
phys_sim.cpp
Go to the documentation of this file.
00001 /* 
00002  * PPRacer 
00003  * Copyright (C) 2004-2005 Volker Stroebel <volker@planetpenguin.de>
00004  *
00005  * Copyright (C) 1999-2001 Jasmin F. Patry
00006  * 
00007  * This program is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU General Public License
00009  * as published by the Free Software Foundation; either version 2
00010  * of the License, or (at your option) any later version.
00011  * 
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  * 
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00020  */
00021 
00022 #include "course_load.h"
00023 #include "course_render.h"
00024 #include "hier.h"
00025 #include "hier_util.h"
00026 #include "model_hndl.h"
00027 #include "part_sys.h"
00028 #include "snow.h"
00029 #include "phys_sim.h"
00030 #include "nmrcl.h"
00031 #include "game_config.h"
00032 #include "course_mgr.h"
00033 #include "stuff.h"
00034 
00035 
00036 #include "ppgltk/audio/audio.h"
00037 #include "ppgltk/alg/defs.h"
00038 
00039 #include "loop.h"
00040 
00041 #include "game_mgr.h"
00042 
00043 
00044 
00045 /*
00046  * Constants
00047  */
00048  
00049 extern terrain_tex_t terrain_texture[NUM_TERRAIN_TYPES];
00050 extern unsigned int num_terrains;
00051 
00052 /* Tux's mass (kg) */
00053 #define TUX_MASS 20
00054 
00055 /* Width of Tux (m) */
00056 #define TUX_WIDTH 0.45
00057 
00058 /* Minimum speed (m/s) */
00059 #define MIN_TUX_SPEED 1.4
00060 
00061 /* Minimum speed for friction to take effect (m/s) */
00062 #define MIN_FRICTION_SPEED 2.8
00063 
00064 /* Initial Tux speed (m/s) */
00065 #define INIT_TUX_SPEED 3.0
00066 
00067 /* Maximum frictional force */
00068 #define MAX_FRICTIONAL_FORCE 800
00069 
00070 /* Maximum angle by which frictional force is rotated when turning*/
00071 #define MAX_TURN_ANGLE 45
00072 
00073 /* Maximum turning force perpendicular to direction of movement */
00074 #define MAX_TURN_PERPENDICULAR_FORCE 400
00075 
00076 /* Maximum increase in friction when turning */
00077 #define MAX_TURN_PENALTY 0.15
00078 
00079 /* Force applied by Tux when braking (N) */
00080 #define BRAKE_FORCE 200
00081 
00082 /* Maximum roll angle (degrees) */
00083 #define MAX_ROLL_ANGLE 30
00084 
00085 /* Roll angle while braking (degrees) */
00086 #define BRAKING_ROLL_ANGLE 55
00087 
00088 /* Baseline friction coeff for rolling */
00089 #define IDEAL_ROLL_FRIC_COEFF 0.35
00090 
00091 /* Baseline speed for rolling (m/s) */
00092 #define IDEAL_ROLL_SPEED 6.0
00093 
00094 /* Tux's orientation is updated using a first-order filter; these are the 
00095    time constants of the filter when tux is on ground and airborne (s) */
00096 #define TUX_ORIENTATION_TIME_CONSTANT 0.14
00097 #define TUX_ORIENTATION_AIRBORNE_TIME_CONSTANT 0.5
00098 
00099 /* Max particles generated by turning (particles/s) */
00100 #define MAX_TURN_PARTICLES 1500
00101 
00102 /* Max particles generated by rolling (particles/s) */
00103 #define MAX_ROLL_PARTICLES 3000
00104 
00105 /* Particles generated by braking (particles/s) */
00106 #define BRAKE_PARTICLES 4000
00107 
00108 /* Number of particles is multiplied by ( speed / PARTICLES_SPEED_FACTOR ) */
00109 #define PARTICLE_SPEED_FACTOR 40
00110 
00111 /* Maximum deflection angle for particles (out to side of player)  */
00112 #define MAX_PARTICLE_ANGLE 80
00113 
00114 /* Speed at which maximum deflection occurs */
00115 #define MAX_PARTICLE_ANGLE_SPEED 50
00116 
00117 /* Particle speed = min ( MAX_PARTICLE_SPEED, speed*PARTICLE_SPEED_MULTIPLIER )
00118  */
00119 #define PARTICLE_SPEED_MULTIPLIER 0.3
00120 #define MAX_PARTICLE_SPEED 2
00121 
00122 /* The compressiblity of the first section of the spring modelling Tux's
00123    tush (m) */
00124 #define TUX_GLUTE_STAGE_1_COMPRESSIBILITY 0.05
00125 
00126 /* The spring coefficient of first section of aforementioned tush (N/m) */
00127 #define TUX_GLUTE_STAGE_1_SPRING_COEFF 1500
00128 
00129 /* The damping coefficient of first section of aforementioned tush (N*s/m) */
00130 #define TUX_GLUTE_STAGE_1_DAMPING_COEFF 100
00131 
00132 /* The compressiblity of the second section of the spring modelling Tux's
00133    tush (m) */
00134 #define TUX_GLUTE_STAGE_2_COMPRESSIBILITY 0.12
00135 
00136 /* The spring coefficient of second section of aforementioned tush (N/m) */
00137 #define TUX_GLUTE_STAGE_2_SPRING_COEFF 3000
00138 
00139 /* The damping coefficient of second section of aforementioned tush (N*s/m) */
00140 #define TUX_GLUTE_STAGE_2_DAMPING_COEFF 500
00141 
00142 /* The spring coefficient of third section of aforementioned tush (N/m) */
00143 #define TUX_GLUTE_STAGE_3_SPRING_COEFF 10000
00144 
00145 /* The damping coefficient of third section of aforementioned tush (N*s/m) */
00146 #define TUX_GLUTE_STAGE_3_DAMPING_COEFF 1000
00147 
00148 /* The maximum force exertedd by both sections of spring (N) */
00149 #define TUX_GLUTE_MAX_SPRING_FORCE 3000
00150 
00151 /* Maximum distance that player can penetrate into terrain (m) */
00152 #define MAX_SURFACE_PENETRATION 0.2
00153 
00154 /* Density of air @ -3 C (kg/m^3) */
00155 #define AIR_DENSITY 1.308
00156 
00157 /* Diameter of sphere used to model Tux for air resistance calcs (m) */
00158 #define TUX_DIAMETER TUX_WIDTH
00159 
00160 /* Viscosity of air @ -3 C (N*s/m^2) */
00161 #define AIR_VISCOSITY 17.00e-6
00162 
00163 /* Drag coefficient vs. Reynolds number table (from Janna, _Introduction to
00164    Fluid Mechanics_, 3rd Edition, PWS-Kent, p 308) 
00165    Values are for Re in the range of 10^-1 to 10^6, spaced logarithmically
00166    (all values are log base 10) */
00167 static const double air_log_re[] = { -1,
00168                                  0,
00169                                  1,
00170                                  2,
00171                                  3,
00172                                  4,
00173                                  5,
00174                                  6 };
00175 static const double air_log_drag_coeff[] = { 2.25,
00176                                         1.35,
00177                                         0.6,
00178                                         0.0,
00179                                         -0.35,
00180                                         -0.45,
00181                                         -0.33,
00182                                         -0.9 };
00183 
00184 /* Minimum time step for ODE solver (s) */
00185 #define MIN_TIME_STEP 0.01
00186 
00187 /* Maximum time step for ODE solver (s) */
00188 #define MAX_TIME_STEP 0.10
00189 
00190 /* Maximum distance of step for ODE solver (m) */
00191 #define MAX_STEP_DISTANCE 0.20
00192 
00193 /* Tolerance on error in Tux's position (m) */
00194 #define MAX_POSITION_ERROR 0.005
00195 
00196 /* Tolerance on error in Tux's velocity (m/s) */
00197 #define MAX_VELOCITY_ERROR 0.05
00198 
00199 /* To smooth out the terrain, the terrain normals are interpolated
00200    near the edges of the triangles.  This parameter controls how much
00201    smoothing is done, as well as the size of the triangle boundary
00202    over which smoothing occurs.  A value of 0 will result in no
00203    smoothing, while a value of 0.5 will smooth the normals over the
00204    entire surface of the triangle (though the smoothing will increase
00205    as the value approaches infinity).  (dimensionless) */
00206 #define NORMAL_INTERPOLATION_LIMIT 0.05
00207 
00208 /* The distance to move Tux up by so that he doesn't sit too low on the 
00209    ground (m) */
00210 #define TUX_Y_CORRECTION_ON_STOMACH 0.33
00211 
00212 /* The square of the distance that Tux must move before we recompute a 
00213    collision (m^2) */
00214 #define COLLISION_TOLERANCE 0.04
00215 
00216 /* Duraction of paddling motion (s) */
00217 #define PADDLING_DURATION 0.40
00218 
00219 /* Force applied against ground when paddling (N) */
00220 #define MAX_PADDLING_FORCE 122.5
00221 
00222 /* Ideal paddling friction coefficient */
00223 #define IDEAL_PADDLING_FRIC_COEFF 0.35
00224 
00225 /* G force applied for jump with charge of 0 */
00226 #define BASE_JUMP_G_FORCE 1.5
00227 
00228 /* Maximum G force applied during jump */
00229 #define MAX_JUMP_G_FORCE 3
00230 
00231 /* Magnitude of force before damage is incurred (N) */
00232 #define DAMAGE_RESISTANCE ( 4.0 * TUX_MASS * EARTH_GRAV )
00233 
00234 /* Damage scaling factor (health/(N*s)) */
00235 #define DAMAGE_SCALE 9e-8
00236 
00237 /* Amount to scale tree force by (so that it does more damage than regular
00238    collisions) (dimensionless) */
00239 #define COLLISION_BASE_DAMAGE 0.02
00240 
00241 /* Damage incurred by paddling exertion (health/s) */
00242 #define PADDLING_DAMAGE 0.02
00243 
00244 
00245 
00246 //WIND
00247 //wind_vel (the wind vector) and wind_scale (its random variations in intensity) are declared (and set) in snow.h
00248 
00249 /*
00250  * Static variables
00251  */
00252 
00253 /* Friction coefficients (dimensionless) */
00254 static double  fricCoeff[4] = { 0.35, /* ice */
00255                                   0.7,  /* rock */
00256                                   0.45,  /* snow */
00257                                                           -0.50  /* ramp */
00258                                 };
00259 
00260 
00261 /* Current time step in ODE solver */
00262 static double ode_time_step = -1;
00263 
00264 class Index2d{
00265 public:
00266        Index2d(){};
00267        Index2d(const int i, const int j)
00268        {
00269               this->i=i;
00270               this->j=j;
00271        }      
00272        int i;   
00273     int j;
00274 };     
00275 
00276 void
00277 set_wind_velocity(pp::Vec3d velocity, float scale)
00278 {
00279        wind_vel = velocity;
00280        wind_scale = scale;                                            
00281 }
00282 
00283 void increment_turn_fact( Player& plyr, double amt )
00284 {
00285     plyr.control.turn_fact += amt;
00286     plyr.control.turn_fact = MIN( 1.0, MAX( plyr.control.turn_fact, -1.0 ) );
00287 }
00288 
00289 double get_min_tux_speed()
00290 {
00291     return MIN_TUX_SPEED;   
00292 }
00293 
00294 void set_friction_coeff( const float fric[3] ) 
00295 {
00296     fricCoeff[0] = fric[0];
00297     fricCoeff[1] = fric[1];
00298     fricCoeff[2] = fric[2];
00299        fricCoeff[3] = fric[3];
00300 } 
00301 
00302 float get_min_y_coord() 
00303 {
00304     float courseWidth, courseLength;
00305     float angle;
00306     float minY;
00307 
00308     get_course_dimensions( &courseWidth, &courseLength );
00309     angle = get_course_angle();
00310 
00311     minY = -courseLength * tan( ANGLES_TO_RADIANS( angle ) );
00312     return minY;
00313 } 
00314 
00315 void get_indices_for_point( double x, double z, 
00316                          int *x0, int *y0, int *x1, int *y1 )
00317 {
00318     float courseWidth, courseLength;
00319     float xidx, yidx;
00320     int nx, ny;
00321 
00322     get_course_dimensions( &courseWidth, &courseLength );
00323     get_course_divisions( &nx, &ny );
00324     
00325     xidx = x / courseWidth * ( (float) nx - 1. );
00326     yidx = -z / courseLength * ( (float) ny - 1. );
00327 
00328     if (xidx < 0) {
00329         xidx = 0;
00330     } else if ( xidx > nx-1 ) {
00331         xidx = nx-1;
00332     } 
00333 
00334     if (yidx < 0) {
00335         yidx = 0;
00336     } else if ( yidx > ny-1 ) {
00337         yidx = ny-1;
00338     } 
00339 
00340     /* I found that ceil(3) was quite slow on at least some architectures, 
00341        so I've replace it with an approximation */
00342     *x0 = (int) (xidx);              /* floor(xidx) */
00343     *x1 = (int) ( xidx + 0.9999 );   /* ceil(xidx) */
00344     *y0 = (int) (yidx);              /* floor(yidx) */
00345     *y1 = (int) ( yidx + 0.9999 );   /* ceil(yidx) */
00346 
00347     if ( *x0 == *x1 ) {
00348        if ( *x1 < nx - 1 ) 
00349            (*x1)++;
00350        else 
00351            (*x0)--;
00352     } 
00353 
00354     if ( *y0 == *y1 ) {
00355        if ( *y1 < ny - 1 )
00356            (*y1)++;
00357        else 
00358            (*y0)--;
00359     } 
00360 
00361     check_assertion( *x0 >= 0, "invalid x0 index" );
00362     check_assertion( *x0 < nx, "invalid x0 index" );
00363     check_assertion( *x1 >= 0, "invalid x1 index" );
00364     check_assertion( *x1 < nx, "invalid x1 index" );
00365     check_assertion( *y0 >= 0, "invalid y0 index" );
00366     check_assertion( *y0 < ny, "invalid y0 index" );
00367     check_assertion( *y1 >= 0, "invalid y1 index" );
00368     check_assertion( *y1 < ny, "invalid y1 index" );
00369 }
00370 
00371 
00372 void find_barycentric_coords( float x, float z, 
00373                            Index2d *idx0, 
00374                            Index2d *idx1, 
00375                            Index2d *idx2,
00376                            float *u, float *v )
00377 {
00378     float xidx, yidx;
00379     int nx, ny;
00380     int x0, x1, y0, y1;
00381     float dx, ex, dz, ez, qx, qz, invdet; // to calc. barycentric coords 
00382     float courseWidth, courseLength;
00383     float *elevation;
00384 
00385     elevation = get_course_elev_data();
00386     get_course_dimensions( &courseWidth, &courseLength );
00387     get_course_divisions( &nx, &ny );
00388 
00389     get_indices_for_point( x, z, &x0, &y0, &x1, &y1 );
00390     
00391     xidx = x / courseWidth * ( (float) nx - 1. );
00392     yidx = -z / courseLength * ( (float) ny - 1. );
00393 
00394 
00395     /* The terrain is meshed as follows:
00396              ...
00397           +-+-+-+-+            x<---+
00398           |\|/|\|/|                 |
00399         ...+-+-+-+-+...              V
00400           |/|\|/|\|                 y
00401           +-+-+-+-+
00402              ...
00403        
00404        So there are two types of squares: those like this (x0+y0 is even):
00405 
00406        +-+
00407        |/|
00408        +-+
00409 
00410        and those like this (x0+y0 is odd).
00411        
00412        +-+
00413        |\|
00414        +-+
00415     */
00416 
00417 
00418     if ( (x0 + y0) % 2 == 0 ) {
00419        if ( yidx - y0 < xidx - x0 ) {
00420            *idx0 = Index2d(x0, y0); 
00421            *idx1 = Index2d(x1, y0); 
00422            *idx2 = Index2d(x1, y1); 
00423        } else {
00424            *idx0 = Index2d(x1, y1); 
00425            *idx1 = Index2d(x0, y1); 
00426            *idx2 = Index2d(x0, y0); 
00427        } 
00428     } else {
00429        // x0 + y0 is odd 
00430        if ( yidx - y0 + xidx - x0 < 1 ) {
00431            *idx0 = Index2d(x0, y0); 
00432            *idx1 = Index2d(x1, y0); 
00433            *idx2 = Index2d(x0, y1); 
00434        } else {
00435            *idx0 = Index2d(x1, y1); 
00436            *idx1 = Index2d(x0, y1); 
00437            *idx2 = Index2d(x1, y0); 
00438        } 
00439     }
00440 
00441     dx = idx0->i - idx2->i;
00442     dz = idx0->j - idx2->j;
00443     ex = idx1->i - idx2->i;
00444     ez = idx1->j - idx2->j;
00445     qx = xidx - idx2->i;
00446     qz = yidx - idx2->j;
00447 
00448     invdet = 1./(dx * ez - dz * ex);
00449     *u = (qx * ez - qz * ex) * invdet;
00450     *v = (qz * dx - qx * dz) * invdet;
00451 }
00452 
00453 
00454 pp::Vec3d find_course_normal( const float x, const float z )
00455 {
00456     pp::Vec3d *course_nmls;
00457     pp::Vec3d tri_nml;
00458     float xidx, yidx;
00459     int nx, ny;
00460     int x0, x1, y0, y1;
00461     pp::Vec3d p0, p1, p2;
00462     Index2d idx0, idx1, idx2;
00463     pp::Vec3d n0, n1, n2;
00464     float course_width, course_length;
00465     float u, v;
00466     float min_bary, interp_factor;
00467     pp::Vec3d smooth_nml;
00468     pp::Vec3d interp_nml;
00469     float *elevation;
00470 
00471     elevation = get_course_elev_data();
00472     get_course_dimensions( &course_width, &course_length );
00473     get_course_divisions( &nx, &ny );
00474     course_nmls = get_course_normals();
00475 
00476     get_indices_for_point( x, z, &x0, &y0, &x1, &y1 );
00477     
00478     xidx = x / course_width * ( (float) nx - 1. );
00479     yidx = -z / course_length * ( (float) ny - 1. );
00480 
00481     find_barycentric_coords( x, z, &idx0, &idx1, &idx2, &u, &v );
00482 
00483     n0 = course_nmls[ idx0.i + nx * idx0.j ];
00484     n1 = course_nmls[ idx1.i + nx * idx1.j ];
00485     n2 = course_nmls[ idx2.i + nx * idx2.j ];
00486 
00487     p0 = COURSE_VERTEX( idx0.i, idx0.j );
00488     p1 = COURSE_VERTEX( idx1.i, idx1.j );
00489     p2 = COURSE_VERTEX( idx2.i, idx2.j );
00490     
00491     smooth_nml = (u*n0)+(v*n1+( 1.-u-v)*n2);
00492 
00493     tri_nml = (p1-p0)^(p2-p0);
00494     tri_nml.normalize();
00495 
00496     min_bary = MIN( u, MIN( v, 1. - u - v ) );
00497     interp_factor = MIN( min_bary / NORMAL_INTERPOLATION_LIMIT, 1.0 );
00498 
00499     interp_nml = (interp_factor*tri_nml)+
00500        ((1.-interp_factor)*smooth_nml);
00501     interp_nml.normalize();
00502 
00503     return interp_nml;
00504 }
00505 
00506 float find_y_coord( float x, float z )
00507 {
00508     float ycoord;
00509     pp::Vec3d p0, p1, p2;
00510     Index2d idx0, idx1, idx2;
00511     float u, v;
00512     int nx, ny;
00513     float *elevation;
00514     float course_width, course_length;
00515 
00516     /* cache last request */
00517     static float last_x, last_z, last_y;
00518     static bool cache_full = false;
00519 
00520     if ( cache_full && last_x == x && last_z == z ) {
00521        return last_y;
00522     }
00523 
00524     get_course_divisions( &nx, &ny );
00525     get_course_dimensions( &course_width, &course_length );
00526     elevation = get_course_elev_data();
00527 
00528     find_barycentric_coords( x, z, &idx0, &idx1, &idx2, &u, &v );
00529 
00530     p0 = COURSE_VERTEX( idx0.i, idx0.j );
00531     p1 = COURSE_VERTEX( idx1.i, idx1.j );
00532     p2 = COURSE_VERTEX( idx2.i, idx2.j );
00533 
00534     ycoord = u * p0.y + v * p1.y + ( 1. - u - v ) * p2.y;
00535 
00536     last_x = x;
00537     last_z = z;
00538     last_y = ycoord;
00539     cache_full = true;
00540 
00541     return ycoord;
00542 } 
00543 
00544 void get_surface_type( float x, float z, float weights[] )
00545 {
00546     int *terrain;
00547     float courseWidth, courseLength;
00548     int nx, ny;
00549     Index2d idx0, idx1, idx2;
00550     float u, v;
00551     unsigned int i;
00552 
00553     find_barycentric_coords( x, z, &idx0, &idx1, &idx2, &u, &v );
00554 
00555     terrain = get_course_terrain_data();
00556     get_course_dimensions( &courseWidth, &courseLength );
00557     get_course_divisions( &nx, &ny );
00558 
00559     for (i=0; i<num_terrains; i++) {
00560        weights[i] = 0;
00561        if ( terrain[ idx0.i + nx*idx0.j ] == int(i) ) {
00562            weights[i] += u;
00563        }
00564        if ( terrain[ idx1.i + nx*idx1.j ] == int(i) ) {
00565            weights[i] += v;
00566        }
00567        if ( terrain[ idx2.i + nx*idx2.j ] == int(i) ) {
00568            weights[i] += 1.0 - u - v;
00569        }
00570     }
00571 } 
00572 
00573 pp::Plane get_local_course_plane( pp::Vec3d pt )
00574 {
00575     pp::Plane plane;
00576 
00577     pt.y = find_y_coord( pt.x, pt.z );
00578 
00579     plane.nml = find_course_normal( pt.x, pt.z );
00580     plane.d = -( plane.nml.x * pt.x + 
00581                plane.nml.y * pt.y +
00582                plane.nml.z * pt.z );
00583 
00584     return plane;
00585 }
00586 
00587 static void update_paddling( Player& plyr )
00588 {
00589     if ( plyr.control.is_paddling ) {
00590               if ( gameMgr->time - plyr.control.paddle_time >= PADDLING_DURATION ) {
00591                   print_debug( DEBUG_CONTROL, "paddling off" );
00592                   plyr.control.is_paddling = false;
00593               }
00594     }
00595 }
00596 
00597 void set_tux_pos( Player& plyr, pp::Vec3d new_pos )
00598 {
00599     char *tuxRoot;
00600     float playWidth, playLength;
00601     float courseWidth, courseLength;
00602     float boundaryWidth;
00603     float disp_y;
00604 
00605     get_play_dimensions( &playWidth, &playLength );
00606     get_course_dimensions( &courseWidth, &courseLength );
00607     boundaryWidth = ( courseWidth - playWidth ) / 2; 
00608 
00609 
00610     if ( new_pos.x < boundaryWidth ) {
00611         new_pos.x = boundaryWidth;
00612     } else if ( new_pos.x > courseWidth - boundaryWidth ) {
00613         new_pos.x = courseWidth - boundaryWidth;
00614     } 
00615 
00616     if ( new_pos.z > 0 ) {
00617         new_pos.z = 0;
00618     } else if ( -new_pos.z >= playLength ) {
00619         new_pos.z = -playLength;
00620         set_game_mode( GAME_OVER );
00621     } 
00622 
00623     plyr.pos = new_pos;
00624 
00625     disp_y = new_pos.y + TUX_Y_CORRECTION_ON_STOMACH; 
00626 
00627     tuxRoot = ModelHndl->get_tux_root_node();
00628     reset_scene_node( tuxRoot );
00629     translate_scene_node( tuxRoot, 
00630                        pp::Vec3d( new_pos.x, disp_y, new_pos.z ) );
00631 } 
00632 
00633 bool check_tree_collisions( Player& plyr, pp::Vec3d pos, 
00634                            pp::Vec3d *tree_loc, float *tree_diam )
00635 {
00636     if(getparam_disable_collision_detection()){
00637               return false;
00638        }      
00639        
00640        Tree *trees;
00641     int num_trees, i;
00642     float diam = 0.0; 
00643     float height;
00644        pp::Polyhedron *ph, ph2;
00645     pp::Matrix mat;
00646     pp::Vec3d loc;
00647     bool hit = false;
00648     pp::Vec3d distvec;
00649     char *tux_root;
00650     float squared_dist;
00651     int tree_type;
00652     
00653     /* These variables are used to cache collision detection results */
00654     static bool last_collision = false;
00655     static pp::Vec3d last_collision_tree_loc( -999, -999, -999 );
00656     static double last_collision_tree_diam = 0;
00657     static pp::Vec3d last_collision_pos( -999, -999, -999 );
00658 
00659     /* If we haven't moved very much since the last call, we re-use
00660        the results of the last call (significant speed-up) */
00661           
00662          
00663     pp::Vec3d dist_vec = pos - last_collision_pos;
00664     if ( dist_vec.length2() < COLLISION_TOLERANCE ) {
00665        if ( last_collision ) {
00666            if ( tree_loc != NULL ) {
00667               *tree_loc = last_collision_tree_loc;
00668            }
00669            if ( tree_diam != NULL ) {
00670               *tree_diam = last_collision_tree_diam;
00671            }
00672            return true;
00673        } else {
00674            return false;
00675        }
00676     }
00677 
00678     trees = get_tree_locs();
00679     num_trees = get_num_trees();
00680     tree_type = trees[0].type;
00681     ph = get_tree_polyhedron(tree_type);
00682 
00683     for (i=0; i<num_trees; i++) {
00684         diam = trees[i].diam;
00685         height = trees[i].height;
00686         loc = trees[i].ray.pt;
00687 
00688         distvec = pp::Vec3d( loc.x - pos.x, 0.0, loc.z - pos.z );
00689 
00690        /* check distance from tree; .6 is the radius of a bounding
00691            sphere around tux (approximate) */
00692        squared_dist = ( diam/2. + 0.6 );
00693        squared_dist *= squared_dist;
00694         if ( distvec.length2() > squared_dist ) {
00695             continue;
00696         } 
00697 
00698        /* have to look at polyhedron - switch to correct one if necessary */
00699        if ( tree_type != trees[i].type )  {
00700            tree_type = trees[i].type;
00701            ph = get_tree_polyhedron(tree_type);
00702               
00703               
00704               
00705        }
00706         ph2 = copy_polyhedron( *ph );
00707        
00708         //mat.makeScaling( diam, height, diam );
00709               mat.makeScaling( diam, diam, diam );
00710        
00711         trans_polyhedron( mat, ph2 );
00712         mat.makeTranslation( loc.x, loc.y, loc.z );
00713         trans_polyhedron( mat, ph2 );
00714 
00715        tux_root = ModelHndl->get_tux_root_node();
00716        reset_scene_node( tux_root );
00717        translate_scene_node( tux_root, 
00718                            pp::Vec3d( pos.x, pos.y, pos.z ) );
00719         hit = collide( tux_root, ph2 );
00720               
00721         free_polyhedron( ph2 );
00722 
00723         if ( hit == true ) {
00724            if ( tree_loc != NULL ) {
00725               *tree_loc = loc;
00726            }
00727            if ( tree_diam != NULL ) {
00728               *tree_diam = diam;
00729            }
00730 
00731 
00732             break;
00733         } 
00734     } 
00735 
00736     last_collision_tree_loc = loc;
00737     last_collision_tree_diam = diam;
00738     last_collision_pos = pos;
00739 
00740     if ( hit ) {
00741        last_collision = true;
00742 
00743        /* Record collision in player data so that health can be adjusted */
00744        plyr.collision = true;
00745     } else {
00746        last_collision = false;
00747     }
00748 
00749     return hit;
00750 } 
00751 
00752 
00753 
00754 
00755 
00756 void check_item_collection( Player& plyr, pp::Vec3d pos )
00757 {
00758     Item *items;
00759     int num_items, i;
00760     float diam = 0.0; 
00761     float height;
00762     pp::Vec3d loc;
00763     pp::Vec3d distvec;
00764     float squared_dist;
00765     int item_type;
00766 
00767     /* These variables are used to cache collision detection results */
00768     static pp::Vec3d last_collision_pos( -999, -999, -999 );
00769 
00770     /* If we haven't moved very much since the last call, we re-use
00771        the results of the last call (significant speed-up) */
00772     pp::Vec3d dist_vec = pos - last_collision_pos;
00773     if ( dist_vec.length2() < COLLISION_TOLERANCE ) {
00774               return ;
00775     }
00776 
00777     items = get_item_locs();
00778     num_items = get_num_items();
00779     item_type = items[0].type;
00780 
00781     for (i=0; i<num_items; i++) {
00782 
00783               /*
00784               if ( items[i].isCollectable()==false) {
00785                   continue;
00786               }
00787               */
00788               if ( items[i].getType() == Item::UNCOLLECTABLE ||
00789                       items[i].isCollected() )
00790               {
00791                   continue;
00792               }
00793 
00794         diam = items[i].diam;
00795         height = items[i].height;
00796         loc = items[i].ray.pt;
00797 
00798         distvec = pp::Vec3d( loc.x - pos.x, 0.0, loc.z - pos.z );
00799 
00800        /* check distance from tree; .6 is the radius of a bounding
00801            sphere around tux (approximate) */
00802        squared_dist = ( diam/2.0 + 0.6 );
00803        squared_dist *= squared_dist;
00804         if ( distvec.length2() > squared_dist ) {
00805             continue;
00806         } 
00807 
00808        if ( (pos.y - 0.6 >= loc.y && pos.y - 0.6 <= loc.y + height) ||
00809             (pos.y + 0.6 >= loc.y && pos.y + 0.6 <= loc.y + height) ||
00810             (pos.y - 0.6 <= loc.y && pos.y + 0.6 >= loc.y + height) ) {
00811            /* close enough to hitting the flag */
00812            /* play a noise? */
00813            items[i].setCollected();
00814               items[i].setDrawable(false);        
00815                       
00816               if(items[i].getType()==Item::HERRING){
00817                      plyr.herring += items[i].getScore();
00818                      if (plyr.herring <0) plyr.herring=0;
00819               }else{
00820                      plyr.incLives();                   
00821               }
00822 
00823            play_sound( "item_collect", 0 );
00824        }
00825 
00826     } 
00827 
00828 } 
00829 
00830 float get_compression_depth( const int terrain ) 
00831 {
00832        return terrain_texture[terrain].compression;
00833 }
00834 
00835 /*
00836  * Check for tree collisions and adjust position and velocity appropriately.
00837  */
00838 static void adjust_for_tree_collision( Player& plyr, 
00839                                    pp::Vec3d pos, pp::Vec3d *vel )
00840 {
00841     pp::Vec3d treeNml;
00842     pp::Vec3d treeLoc;
00843     bool treeHit;
00844     float speed;
00845     float costheta;
00846     float tree_diam;
00847 
00848     treeHit = check_tree_collisions( plyr, pos, &treeLoc, &tree_diam );
00849     if (treeHit) {
00850        /*
00851         * Calculate the normal vector to the tree; here we model the tree
00852         * as a vertical cylinder.
00853         */
00854         treeNml.x = treeLoc.x - pos.x;
00855         treeNml.y = 0;
00856         treeNml.z = treeLoc.z - pos.z;
00857         //Calculate distance from the player location to the center of the tree
00858         float lenNml = treeNml.normalize();
00859         
00860         float speed = vel->length();
00861 
00862         //velocity projected on the tree normal : t
00863         pp::Vec3d vel_t(treeNml);
00864         vel_t = vel_t * (*vel * treeNml);
00865 
00866         //perpendicular component of the velocity vector : n = vel - t
00867         pp::Vec3d vel_n(*vel);
00868         vel_n = *vel + (-1.0)*vel_t ;
00869 
00870         double n_factor = 1.0;
00871         double t_factor = 0.3;
00872         double random_t_factor = 1;//TODO
00873         //random_t_factor = static_cast<double>(std::rand())/static_cast<double>(RAND_MAX);;//random value between 0 and 1
00874         //random_t_factor *=0.3;
00875         //random_t_factor +=0.7;
00876         //We now have a random factor between 0.7 and 1.0
00877         double f_z = 1.0;//TODO
00878 
00879         vel_n=vel_n*n_factor;
00880         vel_t=vel_t*t_factor*random_t_factor*f_z;
00881 
00882         *vel = vel_t + vel_n;
00883         vel->normalize();
00884         *vel = *vel * speed;
00885     } 
00886 }
00887 
00888 /*
00889  * Adjusts velocity so that his speed doesn't drop below the minimum 
00890  * speed
00891  */
00892 double adjust_velocity( pp::Vec3d *vel, pp::Vec3d pos, 
00893                        pp::Plane surf_plane, float dist_from_surface )
00894 {
00895     pp::Vec3d surf_nml;
00896     float speed;
00897 
00898     surf_nml = surf_plane.nml;
00899 
00900     speed = vel->normalize();
00901 
00902     if ( speed < EPS )
00903     {
00904        if ( fabs( surf_nml.x ) + fabs( surf_nml.z ) > EPS ) {
00905            *vel = projectIntoPlane( 
00906               surf_nml, pp::Vec3d( 0.0, -1.0, 0.0 ) );
00907            vel->normalize();
00908        } else {
00909            *vel = pp::Vec3d( 0.0, 0.0, -1.0 );
00910        }
00911     }
00912 
00913     speed = MAX( get_min_tux_speed(), speed );
00914 
00915     *vel = speed*(*vel);
00916     return speed;
00917 }
00918 
00919 void adjust_position( pp::Vec3d *pos, pp::Plane surf_plane, 
00920                     float dist_from_surface )
00921 {
00922 
00923     if ( dist_from_surface < -MAX_SURFACE_PENETRATION )
00924     {
00925        *pos = ( *pos+ (-MAX_SURFACE_PENETRATION - dist_from_surface)* surf_plane.nml);
00926     }
00927     return;
00928 }
00929 
00930 static pp::Vec3d adjust_tux_zvec_for_roll( Player& plyr, 
00931                                      pp::Vec3d vel, pp::Vec3d zvec )
00932 {
00933     pp::Matrix rot_mat; 
00934 
00935     vel = projectIntoPlane( zvec, vel );
00936 
00937     vel.normalize();
00938 
00939     if ( plyr.control.is_braking ) {
00940        rot_mat.makeRotationAboutVector( vel, 
00941                                       plyr.control.turn_fact *
00942                                       BRAKING_ROLL_ANGLE );
00943     } else {
00944        rot_mat.makeRotationAboutVector( vel, 
00945                                       plyr.control.turn_fact *
00946                                       MAX_ROLL_ANGLE );
00947     }
00948 
00949     return rot_mat.transformVector(zvec);
00950 }
00951 
00952 static pp::Vec3d adjust_surf_nml_for_roll( Player& plyr, 
00953                                      pp::Vec3d vel, 
00954                                      float fric_coeff,
00955                                      pp::Vec3d nml )
00956 {
00957     pp::Matrix rot_mat; 
00958     float angle;
00959     float speed;
00960     float roll_angle;
00961 
00962     speed = vel.normalize();
00963 
00964     vel = projectIntoPlane( nml, vel );
00965 
00966     vel.normalize();
00967     if ( plyr.control.is_braking ) {
00968        roll_angle = BRAKING_ROLL_ANGLE;
00969     } else {
00970        roll_angle = MAX_ROLL_ANGLE;
00971     }
00972 
00973        angle = plyr.control.turn_fact *
00974        roll_angle *
00975        MIN( 1.0, MAX(0.0, fric_coeff)/IDEAL_ROLL_FRIC_COEFF )
00976        * MIN(1.0, MAX(0.0,speed-MIN_TUX_SPEED)/
00977              (IDEAL_ROLL_SPEED-MIN_TUX_SPEED));
00978     
00979     rot_mat.makeRotationAboutVector( vel, angle );
00980 
00981     return rot_mat.transformVector(nml);
00982 }
00983 
00984 
00985 void adjust_orientation( Player& plyr, float dtime, 
00986                       pp::Vec3d pos, pp::Vec3d vel,
00987                       float dist_from_surface, pp::Vec3d surf_nml )
00988 {
00989     pp::Vec3d new_x, new_y, new_z; 
00990     pp::Matrix inv_cob_mat;
00991     pp::Matrix rot_mat;
00992     pp::Quat new_orient;
00993     char* tux_root;
00994     float time_constant;
00995     static pp::Vec3d minus_z_vec(0., 0., -1.);
00996     static pp::Vec3d y_vec(0., 1., 0.);
00997 
00998     if ( dist_from_surface > 0 ) {
00999        new_y = 1.*vel;
01000        new_y.normalize();
01001        new_z = projectIntoPlane( new_y, pp::Vec3d(0., -1., 0.) );
01002        new_z.normalize();
01003        new_z = adjust_tux_zvec_for_roll( plyr, vel, new_z );
01004     } else { 
01005        new_z = -1.*surf_nml;
01006        new_z = adjust_tux_zvec_for_roll( plyr, vel, new_z );
01007        new_y = projectIntoPlane( surf_nml,1.*vel);
01008        new_y.normalize();
01009     }
01010 
01011     new_x = new_y^new_z;
01012 
01013        {
01014               pp::Matrix cob_mat;    
01015               pp::Matrix::makeChangeOfBasisMatrix( cob_mat, inv_cob_mat, new_x, new_y, new_z );
01016            new_orient = pp::Quat(cob_mat);
01017        }
01018     if ( !plyr.orientation_initialized ) {
01019        plyr.orientation_initialized = true;
01020        plyr.orientation = new_orient;
01021     }
01022 
01023     time_constant = dist_from_surface > 0
01024        ? TUX_ORIENTATION_AIRBORNE_TIME_CONSTANT
01025        : TUX_ORIENTATION_TIME_CONSTANT;
01026 
01027     plyr.orientation = pp::Quat::interpolate( 
01028        plyr.orientation, new_orient, 
01029        MIN( dtime / time_constant, 1.0 ) );
01030 
01031     plyr.plane_nml = plyr.orientation.rotate( minus_z_vec );
01032     plyr.direction = plyr.orientation.rotate( y_vec );
01033 
01034     pp::Matrix cob_mat( plyr.orientation );
01035 
01036 
01037     /* Trick rotations */
01038     new_y = pp::Vec3d( cob_mat.data[1][0], cob_mat.data[1][1], cob_mat.data[1][2] ); 
01039     rot_mat.makeRotationAboutVector( new_y, 
01040                                    ( plyr.control.barrel_roll_factor * 360 ) );
01041     cob_mat=rot_mat*cob_mat;
01042     new_x = pp::Vec3d( cob_mat.data[0][0], cob_mat.data[0][1], cob_mat.data[0][2] ); 
01043     rot_mat.makeRotationAboutVector( new_x, 
01044                                    plyr.control.flip_factor * 360 );
01045     cob_mat=rot_mat*cob_mat;
01046 
01047 
01048 
01049     inv_cob_mat.transpose(cob_mat);
01050 
01051     tux_root = ModelHndl->get_tux_root_node();
01052     transform_scene_node( tux_root, cob_mat, inv_cob_mat ); 
01053 }
01054 
01055 float adjust_particle_count( float particles ) 
01056 {
01057     if ( particles < 1 ) {
01058        if ( ( (float) rand() ) / RAND_MAX < particles ) {
01059            return 1.0;
01060        } else {
01061            return 0.0;
01062        }
01063     } else {
01064        return particles;
01065     }
01066 }
01067 
01068 void generate_particles( Player& plyr, float dtime, 
01069                       pp::Vec3d pos, float speed ) 
01070 {
01071     pp::Vec3d left_part_pt, right_part_pt;
01072     float brake_particles;
01073     float turn_particles;
01074     float roll_particles;
01075     float surf_weights[NUM_TERRAIN_TYPES]; 
01076     float surf_y;
01077     float left_particles, right_particles;
01078     pp::Vec3d left_part_vel, right_part_vel;
01079     pp::Matrix rot_mat;
01080     pp::Vec3d xvec;
01081        bool particles_type;
01082        GLuint particle_binding = 0;
01083        
01084        unsigned int i;
01085        
01086     get_surface_type( pos.x, pos.z, surf_weights );
01087     surf_y = find_y_coord( pos.x, pos.z );
01088 
01089        
01090        particles_type=false;
01091        for (i=0;i<num_terrains;i++){
01092               if (terrain_texture[i].partbind != 0){    
01093                      if (surf_weights[i] > 0.5) {
01094                             particles_type=true;
01095                             particle_binding = terrain_texture[i].partbind;
01096               }
01097               }
01098        }
01099        
01100        
01101        
01102        
01103     if ( particles_type== true && pos.y < surf_y ) {
01104        xvec =plyr.direction^plyr.plane_nml;
01105 
01106         right_part_pt = left_part_pt = pos;
01107 
01108        right_part_pt = right_part_pt + ((TUX_WIDTH/2.0)*xvec );
01109 
01110        left_part_pt = left_part_pt + ((-TUX_WIDTH/2.0)* xvec);
01111 
01112         right_part_pt.y = left_part_pt.y  = surf_y;
01113 
01114        brake_particles = dtime *
01115            BRAKE_PARTICLES * ( plyr.control.is_braking ? 1.0 : 0.0 )
01116            * MIN( speed / PARTICLE_SPEED_FACTOR, 1.0 );
01117        turn_particles = dtime * MAX_TURN_PARTICLES 
01118            * MIN( speed / PARTICLE_SPEED_FACTOR, 1.0 );
01119        roll_particles = dtime * MAX_ROLL_PARTICLES 
01120            * MIN( speed / PARTICLE_SPEED_FACTOR, 1.0 );
01121 
01122        left_particles = turn_particles * 
01123            fabs( MIN(plyr.control.turn_fact, 0.) ) + 
01124            brake_particles +
01125            roll_particles * fabs( MIN(plyr.control.turn_animation, 0.) );
01126 
01127        right_particles = turn_particles * 
01128            fabs( MAX(plyr.control.turn_fact, 0.) ) + 
01129            brake_particles +
01130            roll_particles * fabs( MAX(plyr.control.turn_animation, 0.) );
01131 
01132        left_particles = adjust_particle_count( left_particles );
01133        right_particles = adjust_particle_count( right_particles );
01134 
01135        /* Create particle velocitites */
01136        rot_mat.makeRotationAboutVector( plyr.direction,
01137            MAX( -MAX_PARTICLE_ANGLE, 
01138                -MAX_PARTICLE_ANGLE * speed / MAX_PARTICLE_ANGLE_SPEED ) );
01139        left_part_vel = rot_mat.transformVector(plyr.plane_nml);
01140        left_part_vel = MIN( MAX_PARTICLE_SPEED, 
01141                                       speed * PARTICLE_SPEED_MULTIPLIER )*
01142                                   left_part_vel;
01143        rot_mat.makeRotationAboutVector( plyr.direction,
01144            MIN( MAX_PARTICLE_ANGLE, 
01145                MAX_PARTICLE_ANGLE * speed / MAX_PARTICLE_ANGLE_SPEED ) );
01146        right_part_vel = rot_mat.transformVector( plyr.plane_nml );
01147        right_part_vel = MIN( MAX_PARTICLE_SPEED,
01148                                        speed * PARTICLE_SPEED_MULTIPLIER )*
01149                                    right_part_vel;
01150 
01151         create_new_particles( left_part_pt, left_part_vel, 
01152                            (int)left_particles, particle_binding );
01153         create_new_particles( right_part_pt, right_part_vel, 
01154                            (int)right_particles, particle_binding );
01155     } 
01156 }
01157 
01158 
01159 /*
01160  * Calculate the magnitude of force due to air resistance (wind)
01161  */
01162 pp::Vec3d calc_wind_force( pp::Vec3d player_vel )
01163 {
01164     pp::Vec3d total_vel;
01165     float wind_speed;
01166     float re;         /* Reynolds number */
01167     float df_mag;     /* magnitude of drag force */
01168     float drag_coeff; /* drag coefficient */
01169     int table_size;      /* size of drag coeff table */
01170 
01171     total_vel = -1*player_vel;
01172     
01173     if ( gameMgr->getCurrentRace().windy ) {
01174         total_vel = total_vel+(wind_scale*wind_vel);
01175     }
01176 
01177     wind_speed = total_vel.normalize();
01178 
01179     re = AIR_DENSITY * wind_speed * TUX_DIAMETER / AIR_VISCOSITY;
01180 
01181     table_size = sizeof(air_log_drag_coeff) / sizeof(air_log_drag_coeff[0]);
01182 
01183     drag_coeff = pow( 10.0, 
01184                     lin_interp( air_log_re, air_log_drag_coeff, 
01185                               log10(re), table_size ) );
01186 
01187     df_mag = 0.5 * drag_coeff * AIR_DENSITY * ( wind_speed * wind_speed )
01188        * ( M_PI * ( TUX_DIAMETER * TUX_DIAMETER ) * 0.25 );
01189 
01190     check_assertion( df_mag > 0, "Negative wind force" );
01191 
01192     return df_mag*total_vel;
01193 }
01194 
01195 static pp::Vec3d calc_spring_force( float compression, pp::Vec3d vel, 
01196                                pp::Vec3d surf_nml, pp::Vec3d *unclamped_f )
01197 {
01198     float spring_vel; /* velocity perp. to surface (for damping) */
01199     float spring_f_mag; /* magnitude of force */
01200 
01201    check_assertion( compression >= 0, 
01202                   "spring can't have negative compression" );
01203     
01204     spring_vel =vel*surf_nml;
01205 
01206     /* Stage 1 */
01207     spring_f_mag = 
01208        MIN( compression, TUX_GLUTE_STAGE_1_COMPRESSIBILITY ) 
01209        * TUX_GLUTE_STAGE_1_SPRING_COEFF;
01210 
01211     /* Stage 2 */
01212     spring_f_mag += 
01213        MAX( 0, MIN( compression - TUX_GLUTE_STAGE_1_COMPRESSIBILITY, 
01214                    TUX_GLUTE_STAGE_2_COMPRESSIBILITY ) )
01215        * TUX_GLUTE_STAGE_2_SPRING_COEFF;
01216 
01217     /* Stage 3 */
01218     spring_f_mag += 
01219        MAX( 0., compression - TUX_GLUTE_STAGE_2_COMPRESSIBILITY - 
01220             TUX_GLUTE_STAGE_1_COMPRESSIBILITY ) 
01221        * TUX_GLUTE_STAGE_3_SPRING_COEFF;
01222 
01223     /* Damping */
01224     spring_f_mag -= spring_vel * ( 
01225        compression <= TUX_GLUTE_STAGE_1_COMPRESSIBILITY 
01226        ? TUX_GLUTE_STAGE_1_SPRING_COEFF : 
01227        ( compression <= TUX_GLUTE_STAGE_2_COMPRESSIBILITY 
01228          ? TUX_GLUTE_STAGE_2_DAMPING_COEFF
01229          : TUX_GLUTE_STAGE_3_DAMPING_COEFF ) );
01230 
01231     /* Clamp to >= 0.0 */
01232     spring_f_mag = MAX( spring_f_mag, 0.0 );
01233 
01234     if ( unclamped_f != NULL ) {
01235        *unclamped_f = spring_f_mag*surf_nml;
01236     }
01237 
01238     /* Clamp to <= TUX_GLUTE_MAX_SPRING_FORCE */
01239     spring_f_mag = MIN( spring_f_mag, TUX_GLUTE_MAX_SPRING_FORCE );
01240 
01241     return spring_f_mag*surf_nml;
01242 }
01243 
01244 
01245 static pp::Vec3d calc_net_force( Player& plyr, pp::Vec3d pos, 
01246                             pp::Vec3d vel )
01247 {
01248     pp::Vec3d nml_f;      /* normal force */
01249     pp::Vec3d unclamped_nml_f; /* unclamped normal force (for damage calc) */
01250     pp::Vec3d fric_f;     /* frictional force */
01251     float fric_f_mag; /* frictional force magnitude */
01252     pp::Vec3d fric_dir;   /* direction of frictional force */
01253     pp::Vec3d grav_f;     /* gravitational force */
01254     pp::Vec3d air_f;      /* air resistance force */
01255     pp::Vec3d brake_f;    /* braking force */
01256     pp::Vec3d paddling_f; /* paddling force */
01257     pp::Vec3d net_force;  /* the net force (sum of all other forces) */
01258     float comp_depth; /* depth to which the terrain can be compressed */
01259     float speed;      /* speed (m/s) */
01260     pp::Vec3d orig_surf_nml; /* normal to terrain at current position */
01261     pp::Vec3d surf_nml;   /* normal to terrain w/ roll effect */
01262     float glute_compression; /* amt that Tux's tush has been compressed */
01263     float steer_angle; /* Angle to rotate fricitonal force for turning */
01264     pp::Matrix fric_rot_mat; /* Matrix to rotate frictional force */
01265     pp::Vec3d jump_f;
01266     pp::Plane surf_plane;
01267     float dist_from_surface;
01268     float surf_weights[NUM_TERRAIN_TYPES];
01269     float surf_fric_coeff;
01270     unsigned int i;
01271 
01272     get_surface_type( pos.x, pos.z, surf_weights );
01273     surf_plane = get_local_course_plane( pos );
01274     orig_surf_nml = surf_plane.nml;
01275 
01276     surf_fric_coeff = 0;
01277     for (i=0; i<num_terrains; i++) {
01278        /*
01279        surf_fric_coeff += surf_weights[i] * fricCoeff[terrain_texture[i].type];
01280        */
01281        surf_fric_coeff += surf_weights[i] * terrain_texture[i].friction;
01282     }
01283     surf_nml = adjust_surf_nml_for_roll( plyr, vel, surf_fric_coeff,
01284                                     orig_surf_nml );
01285     
01286     comp_depth = 0;
01287     for (i=0; i<num_terrains; i++) {
01288        comp_depth += surf_weights[i] * get_compression_depth(i);
01289     }
01290 
01291 
01292     grav_f = pp::Vec3d( 0, -EARTH_GRAV * TUX_MASS, 0 );
01293 
01294     dist_from_surface = surf_plane.distance( pos );
01295 
01296     if ( dist_from_surface <= 0 ) {
01297        plyr.airborne = false;
01298     } else {
01299        plyr.airborne = true;
01300     }
01301 
01302     /*
01303      * Calculate normal force
01304      */
01305     nml_f = pp::Vec3d( 0., 0., 0. );
01306     if ( dist_from_surface <= -comp_depth ) {
01307        /*
01308         * Tux ended up below the surface.
01309         * Calculate the spring force exterted by his rear end. :-)
01310         */
01311        glute_compression = -dist_from_surface - comp_depth;
01312        check_assertion( glute_compression >= 0, 
01313                       "unexpected negative compression" );
01314 
01315        nml_f = calc_spring_force( glute_compression, vel, surf_nml,
01316                                &unclamped_nml_f );
01317     }
01318 
01319     /* Check if player is trying to jump */
01320     if ( plyr.control.begin_jump == true ) {
01321        plyr.control.begin_jump = false;
01322        if ( dist_from_surface <= 0 ) {
01323            plyr.control.jumping = true;
01324            plyr.control.jump_start_time = gameMgr->time;
01325        } else {
01326            plyr.control.jumping = false;
01327        }
01328     }
01329 
01330 
01331     /* Apply jump force in up direction for JUMP_FORCE_DURATION */
01332     if ( ( plyr.control.jumping ) &&
01333         ( gameMgr->time - plyr.control.jump_start_time < 
01334           JUMP_FORCE_DURATION ) ) 
01335     {
01336        jump_f = pp::Vec3d( 
01337            0, 
01338            BASE_JUMP_G_FORCE * TUX_MASS * EARTH_GRAV + 
01339            plyr.control.jump_amt * 
01340            (MAX_JUMP_G_FORCE-BASE_JUMP_G_FORCE) * TUX_MASS * EARTH_GRAV, 
01341            0 );
01342 
01343     } else {
01344        jump_f = pp::Vec3d( 0, 0, 0 );
01345        plyr.control.jumping = false;
01346     }
01347 
01348     /* Use the unclamped normal force for damage calculation purposes */
01349     plyr.normal_force = unclamped_nml_f;
01350 
01351     /* 
01352      * Calculate frictional force
01353      */
01354     fric_dir = vel;
01355     speed = fric_dir.normalize();
01356     fric_dir = -1.0*fric_dir;
01357     
01358     if ( dist_from_surface < 0 && speed > MIN_FRICTION_SPEED ) {
01359        pp::Vec3d tmp_nml_f = nml_f;
01360 
01361        fric_f_mag = tmp_nml_f.normalize() * surf_fric_coeff;
01362 
01363        fric_f_mag = MIN( MAX_FRICTIONAL_FORCE, fric_f_mag );
01364 
01365        fric_f = fric_f_mag*fric_dir;
01366 
01367 
01368        /*
01369         * Adjust friction for steering
01370         */
01371        steer_angle = plyr.control.turn_fact * MAX_TURN_ANGLE;
01372 
01373        if ( fabs( fric_f_mag * sin( ANGLES_TO_RADIANS( steer_angle ) ) ) >
01374             MAX_TURN_PERPENDICULAR_FORCE ) 
01375        {
01376            //check_assertion( fabs( plyr->control.turn_fact ) > 0,
01377               //          "steer angle is non-zero when player is not "
01378               //          "turning" );
01379            steer_angle = RADIANS_TO_ANGLES( 
01380               asin( MAX_TURN_PERPENDICULAR_FORCE / fric_f_mag ) ) * 
01381               plyr.control.turn_fact / fabs( plyr.control.turn_fact );
01382        }
01383        fric_rot_mat.makeRotationAboutVector( orig_surf_nml, steer_angle );
01384        fric_f = fric_rot_mat.transformVector( fric_f );
01385        fric_f = (1.0 + MAX_TURN_PENALTY)*fric_f;
01386 
01387 
01388        /*
01389         * Calculate braking force
01390         */
01391        if ( speed > get_min_tux_speed() && plyr.control.is_braking ) {
01392            brake_f = (surf_fric_coeff * BRAKE_FORCE)*fric_dir; 
01393        } else {
01394            brake_f = pp::Vec3d( 0., 0., 0. );
01395        }
01396        
01397     } else {
01398        fric_f = brake_f = pp::Vec3d( 0., 0., 0. );
01399     }
01400 
01401     /*
01402      * Calculate air resistance
01403      */
01404     air_f = calc_wind_force( vel );
01405 
01406 
01407     /*
01408      * Calculate force from paddling
01409      */
01410     update_paddling( plyr );
01411     if ( plyr.control.is_paddling ) {
01412        if ( plyr.airborne ) {
01413            paddling_f = pp::Vec3d( 0, 0, -TUX_MASS * EARTH_GRAV / 4.0 );
01414            paddling_f = plyr.orientation.rotate( paddling_f );
01415        } else {
01416            paddling_f = ( 
01417               -1 * MIN( MAX_PADDLING_FORCE,  
01418                        MAX_PADDLING_FORCE * 
01419                        ( MAX_PADDLING_SPEED - speed ) / MAX_PADDLING_SPEED * 
01420                        MIN(1.0, 
01421                            surf_fric_coeff/IDEAL_PADDLING_FRIC_COEFF)))*
01422               fric_dir;
01423        }
01424     } else {
01425        paddling_f = pp::Vec3d( 0., 0., 0. );
01426     }
01427 
01428     
01429     /*
01430      * Add all the forces 
01431      */
01432     net_force = jump_f+grav_f+nml_f+fric_f+air_f+brake_f+paddling_f;
01433 
01434     return net_force;
01435 }
01436 
01437 static float adjust_time_step_size( float h, pp::Vec3d vel )
01438 {
01439     float speed;
01440 
01441     speed = vel.normalize();
01442 
01443     h = MAX( h, MIN_TIME_STEP );
01444     h = MIN( h, MAX_STEP_DISTANCE / speed );
01445     h = MIN( h, MAX_TIME_STEP );
01446 
01447     return h;
01448 }
01449 
01450 /*
01451  * Solves the system of ordinary differential equations governing Tux's 
01452  * movement.  Based on Matlab's ode23.m, (c) The MathWorks, Inc.
01453  */
01454 void solve_ode_system( Player& plyr, float dtime ) 
01455 {
01456     float t0, t, tfinal;
01457     ode_data_t *x, *y, *z, *vx, *vy, *vz; /* store estimates of derivs */
01458     ode_solver_t solver;
01459     float h;
01460     bool done = false;
01461     bool failed = false;
01462     pp::Vec3d new_pos;
01463     pp::Vec3d new_vel, tmp_vel;
01464     float speed;
01465     pp::Vec3d new_f;
01466     pp::Vec3d saved_pos;
01467     pp::Vec3d saved_vel, saved_f;
01468     float pos_err[3], vel_err[3], tot_pos_err, tot_vel_err;
01469     float err=0, tol=0;
01470     int i;
01471 
01472     /* Select a solver */
01473     switch ( getparam_ode_solver() ) {
01474     case EULER:
01475        solver = new_euler_solver();
01476        break;
01477     case ODE23:
01478        solver = new_ode23_solver();
01479        break;
01480     case ODE45:
01481        solver = new_ode45_solver();
01482        break;
01483     default:
01484        setparam_ode_solver( ODE23 );
01485        solver = new_ode23_solver();
01486     }
01487 
01488     /* Select an initial time step */
01489     h = ode_time_step;
01490 
01491     if ( h < 0 || solver.estimate_error == NULL ) {
01492        h = adjust_time_step_size( dtime, plyr.vel );
01493     }
01494 
01495     t0 = 0;
01496     tfinal = dtime;
01497 
01498     t = t0;
01499 
01500     /* Create variables to store derivative estimates & other data */
01501     x  = solver.new_ode_data();
01502     y  = solver.new_ode_data();
01503     z  = solver.new_ode_data();
01504     vx = solver.new_ode_data();
01505     vy = solver.new_ode_data();
01506     vz = solver.new_ode_data();
01507 
01508     /* initialize state */
01509     new_pos = plyr.pos;
01510     new_vel = plyr.vel;
01511     new_f   = plyr.net_force;
01512 
01513     /* loop until we've integrated from t0 to tfinal */
01514     while (!done) {
01515 
01516        if ( t >= tfinal ) {
01517            print_warning( CRITICAL_WARNING, 
01518                         "t >= tfinal in solve_ode_system()" );
01519            break;
01520        }
01521 
01522        /* extend h by up to 10% to reach tfinal */
01523        if ( 1.1 * h > tfinal - t ) {
01524            h = tfinal-t;
01525            check_assertion( h >= 0., "integrated past tfinal" );
01526            done = true;
01527        }
01528 
01529         print_debug( DEBUG_ODE, "h: %g", h );
01530 
01531        saved_pos = new_pos;
01532        saved_vel = new_vel;
01533        saved_f = new_f;
01534 
01535        /* Loop until error is acceptable */
01536        failed = false;
01537 
01538        for (;;) {
01539 
01540            /* Store initial conditions */
01541            solver.init_ode_data( x, new_pos.x, h );
01542            solver.init_ode_data( y, new_pos.y, h );
01543            solver.init_ode_data( z, new_pos.z, h );
01544            solver.init_ode_data( vx, new_vel.x, h );
01545            solver.init_ode_data( vy, new_vel.y, h );
01546            solver.init_ode_data( vz, new_vel.z, h );
01547 
01548 
01549            /* We assume that the first estimate in all ODE solvers is equal 
01550               to the final value of the last iteration */
01551            solver.update_estimate( x, 0, new_vel.x );
01552            solver.update_estimate( y, 0, new_vel.y );
01553            solver.update_estimate( z, 0, new_vel.z );
01554            solver.update_estimate( vx, 0, new_f.x / TUX_MASS );
01555            solver.update_estimate( vy, 0, new_f.y / TUX_MASS );
01556            solver.update_estimate( vz, 0, new_f.z / TUX_MASS );
01557 
01558            /* Update remaining estimates */
01559            for ( i=1; i < solver.num_estimates(); i++ ) {
01560               new_pos.x = solver.next_val( x, i );
01561               new_pos.y = solver.next_val( y, i );
01562               new_pos.z = solver.next_val( z, i );
01563               new_vel.x = solver.next_val( vx, i );
01564               new_vel.y = solver.next_val( vy, i );
01565               new_vel.z = solver.next_val( vz, i );
01566 
01567               solver.update_estimate( x, i, new_vel.x );
01568               solver.update_estimate( y, i, new_vel.y );
01569               solver.update_estimate( z, i, new_vel.z );
01570 
01571               new_f = calc_net_force( plyr, new_pos, new_vel );
01572 
01573               solver.update_estimate( vx, i, new_f.x / TUX_MASS );
01574               solver.update_estimate( vy, i, new_f.y / TUX_MASS );
01575               solver.update_estimate( vz, i, new_f.z / TUX_MASS );
01576            }
01577 
01578            /* Get final values */
01579            new_pos.x = solver.final_estimate( x );
01580            new_pos.y = solver.final_estimate( y );
01581            new_pos.z = solver.final_estimate( z );
01582 
01583            new_vel.x = solver.final_estimate( vx );
01584            new_vel.y = solver.final_estimate( vy );
01585            new_vel.z = solver.final_estimate( vz );
01586 
01587            /* If the current solver can provide error estimates, update h
01588               based on the error, and re-evaluate this step if the error is 
01589               too large  */
01590            if ( solver.estimate_error != NULL ) {
01591 
01592               /* Calculate the error */
01593               pos_err[0] = solver.estimate_error( x );
01594               pos_err[1] = solver.estimate_error( y );
01595               pos_err[2] = solver.estimate_error( z );
01596 
01597               vel_err[0] = solver.estimate_error( vx );
01598               vel_err[1] = solver.estimate_error( vy );
01599               vel_err[2] = solver.estimate_error( vz );
01600 
01601               tot_pos_err = 0.;
01602               tot_vel_err = 0.;
01603               for ( i=0; i<3; i++ ) {
01604                   pos_err[i] *= pos_err[i];
01605                   tot_pos_err += pos_err[i];
01606                   vel_err[i] *= vel_err[i];
01607                   tot_vel_err += vel_err[i];
01608               }
01609               tot_pos_err = sqrt( tot_pos_err );
01610               tot_vel_err = sqrt( tot_vel_err );
01611 
01612                 print_debug( DEBUG_ODE, "pos_err: %g, vel_err: %g", 
01613                           tot_pos_err, tot_vel_err );
01614 
01615               if ( tot_pos_err / MAX_POSITION_ERROR >
01616                    tot_vel_err / MAX_VELOCITY_ERROR )
01617               {
01618                   err = tot_pos_err;
01619                   tol = MAX_POSITION_ERROR;
01620               } else {
01621                   err = tot_vel_err;
01622                   tol = MAX_VELOCITY_ERROR;
01623               }
01624 
01625               /* Update h based on error */
01626               if (  err > tol  && h > MIN_TIME_STEP + EPS  )
01627               {
01628                   done = false;
01629                   if ( !failed ) {
01630                      failed = true;
01631                      h *=  MAX( 0.5, 0.8 *
01632                                pow( tol/err, 
01633                                    float(solver.time_step_exponent()) ) );
01634                   } else {
01635                      h *= 0.5;
01636                   }
01637 
01638                   h = adjust_time_step_size( h, saved_vel );
01639 
01640                   new_pos = saved_pos;
01641                   new_vel = saved_vel;
01642                   new_f = saved_f;
01643               
01644               } else {
01645                   /* Error is acceptable or h is at its minimum; stop */
01646                   break;
01647               }
01648 
01649            } else {
01650               /* Current solver doesn't provide error estimates; stop */
01651               break;
01652            }
01653        }
01654 
01655        /* Update time */
01656        t = t + h;
01657 
01658        tmp_vel = new_vel;
01659        speed = tmp_vel.normalize();
01660 
01661        /* only generate particles if we're drawing them */
01662        if ( getparam_draw_particles() ) {
01663            generate_particles( plyr, h, new_pos, speed );
01664        }
01665 
01666        /* Calculate the final net force */
01667        new_f = calc_net_force( plyr, new_pos, new_vel );
01668 
01669        /* If no failures, compute a new h */
01670        if ( !failed && solver.estimate_error != NULL ) {
01671               double temp = 1.25 * pow((double)err / tol,(double)solver.time_step_exponent());
01672               if ( temp > 0.2 ) {
01673               h = h / temp;
01674            } else {
01675               h = 5.0 * h;
01676            }
01677        }
01678 
01679        /* Clamp h to constraints */
01680        h = adjust_time_step_size( h, new_vel );
01681 
01682        /* Important: to make trees "solid", we must manipulate the 
01683           velocity here; if we don't and Tux is moving very quickly,
01684           he can pass through trees */
01685        adjust_for_tree_collision( plyr, new_pos, &new_vel );
01686 
01687        /* Try to collect items here */
01688        check_item_collection( plyr, new_pos );
01689 
01690     }
01691 
01692     /* Save time step for next time */
01693     ode_time_step = h;
01694 
01695     plyr.vel = new_vel;
01696     plyr.pos = new_pos;
01697     plyr.net_force = new_f;
01698 
01699     free( x );
01700     free( y );
01701     free( z );
01702     free( vx );
01703     free( vy );
01704     free( vz );
01705 }
01706 
01707 
01708 /* 
01709  * Updates Tux's position taking into account gravity, friction, tree 
01710  * collisions, etc.  This is the main physics function.
01711  */
01712 void update_player_pos( Player& plyr, float dtime )
01713 {
01714     pp::Vec3d surf_nml;   /* normal vector of terrain */
01715     pp::Vec3d tmp_vel;
01716     float speed;
01717     float paddling_factor; 
01718     pp::Vec3d local_force;
01719     float flap_factor;
01720     pp::Plane surf_plane;
01721     float dist_from_surface;
01722 
01723     if ( dtime > 2. * EPS ) {
01724        solve_ode_system( plyr, dtime );
01725     }
01726 
01727     tmp_vel = plyr.vel;
01728 
01729     /*
01730      * Set position, orientation, generate particles
01731      */
01732     surf_plane = get_local_course_plane( plyr.pos );
01733     surf_nml = surf_plane.nml;
01734     dist_from_surface = surf_plane.distance( plyr.pos );
01735     adjust_velocity( &plyr.vel, plyr.pos, surf_plane, 
01736                    dist_from_surface );
01737     adjust_position( &plyr.pos, surf_plane, dist_from_surface );
01738 
01739     speed = tmp_vel.normalize();
01740 
01741     set_tux_pos( plyr, plyr.pos );
01742     adjust_orientation( plyr, dtime, plyr.pos, plyr.vel, 
01743                      dist_from_surface, surf_nml );
01744 
01745     flap_factor = 0;
01746 
01747     if ( plyr.control.is_paddling ) {
01748        double factor;
01749        factor = (gameMgr->time - plyr.control.paddle_time) / PADDLING_DURATION;
01750        if ( plyr.airborne ) {
01751            paddling_factor = 0;
01752            flap_factor = factor;
01753        } else {
01754            paddling_factor = factor;
01755            flap_factor = 0;
01756        }
01757     } else {
01758        paddling_factor = 0.0;
01759     }
01760 
01761     /* calculate force in Tux's local coordinate system */
01762     local_force = plyr.orientation.conjugate().rotate(plyr.net_force);
01763 
01764     if (plyr.control.jumping) {
01765        flap_factor = (gameMgr->time - plyr.control.jump_start_time) / 
01766            JUMP_FORCE_DURATION;
01767     } 
01768 
01769     ModelHndl->adjust_tux_joints( plyr.control.turn_animation, plyr.control.is_braking,
01770                      paddling_factor, speed, local_force, flap_factor );
01771 }
01772 
01773 
01774 void init_physical_simulation()
01775 {
01776     pp::Vec3d nml;
01777     pp::Matrix rotMat;
01778     float ycoord;
01779     Player* plyr;
01780     pp::Vec3d init_vel;
01781     pp::Vec3d init_f;
01782     int i;
01783 
01784     for ( i=0; i<gameMgr->numPlayers; i++ ) {
01785               plyr = &players[i];
01786 
01787        ycoord = find_y_coord( plyr->pos.x, plyr->pos.z );
01788        nml = find_course_normal(plyr->pos.x, plyr->pos.z );
01789        rotMat.makeRotation( -90., 'x' );
01790        init_vel = rotMat.transformVector( nml );
01791        init_vel = INIT_TUX_SPEED*init_vel;
01792        init_f = pp::Vec3d( 0., 0., 0. );
01793 
01794        plyr->pos.y = ycoord;
01795        plyr->vel = init_vel;
01796        plyr->net_force = init_f;
01797        plyr->control.turn_fact = 0.0;
01798        plyr->control.turn_animation = 0.0;
01799        plyr->control.barrel_roll_factor = 0.0;
01800        plyr->control.flip_factor = 0.0;
01801        plyr->control.is_braking = false;
01802        plyr->orientation_initialized = false;
01803        plyr->plane_nml = nml;
01804        plyr->direction = init_vel;
01805        plyr->normal_force = pp::Vec3d(0,0,0);
01806        plyr->airborne = false;
01807        plyr->collision = false;
01808        plyr->control.jump_amt = 0;
01809        plyr->control.is_paddling = false;
01810        plyr->control.jumping = false;
01811        plyr->control.jump_charging = false;
01812        plyr->control.barrel_roll_left = false;
01813        plyr->control.barrel_roll_right = false;
01814        plyr->control.barrel_roll_factor = 0;
01815        plyr->control.front_flip = false;
01816        plyr->control.back_flip = false;
01817     }
01818 
01819     ode_time_step = -1;
01820 }