rotation - Inertia Tensor (OpenGL) -

i have been working on rigid body mechanics system week , not sure doing wrong. have rotation matrix computes rotations in yxz order. rotations computed based off of previous rotation , updated offset. inertia tensor computed every time detect collision using rotation matrix i.e. = rot*i0*rotnot. i0 constant inertial frame when rotations lined principle axis. trying compute resultant angular , linear velocity of object after collision. code have below messy commented , based on off of wikipedia page http://en.wikipedia.org/wiki/collision_response scalar jr computation , angular , linear velocity computations.

//calculate inertia tensor based on rotation matrix     matrix rotmat = rectangle1->getrotmat();     matrix trotmat = rotmat.transpose();     i1 = rotmat*i1*trotmat;      //inverse of i1     //determinant of i1     float det1 = i1.getdeterminant();      //transpose of i1     matrix i1t = i1.transpose();      //cofactors     float *a = i1t.getmatrix();     matrix i1not(3,3);     i1not.setmatrix((a[4]*a[8] - a[5]*a[7]), -(a[3]*a[8] - a[5]*a[6]), (a[3]*a[7] - a[4]*a[6]),         -(a[1]*a[8] - a[2]*a[7]), (a[0]*a[8] - a[2]*a[6]), -(a[0]*a[7] - a[1]*a[6]),         (a[1]*a[5] - a[2]*a[4]), -(a[0]*a[5] - a[2]*a[3]), (a[0]*a[4] - a[1]*a[3]));      //scale determinant find inverse     i1not = i1not/det1;      //compute impluse     float e = 1.0f;     vector r1 = hitpoint - rectangle1->getlinpos();      vector r2 = hitpoint - rectangle2->getlinpos();      vector n = tempnorm;     n = n*-1.0f;     float jrnumer = ((rectangle2velocity-rectangle1velocity)*-(1.0f + e)).dotproduct(n);       //compute (inot*(r x n) x r)     vector tr1 = r1;     // tr1 = (r x n)     tr1.crossproduct(n);     float *i1a = i1not.getmatrix();     //(i1ir1n = inot * tr1)     vector i1ir1n(i1a[0]*tr1.getx()+i1a[1]*tr1.gety()+i1a[2]*tr1.getz(),         i1a[3]*tr1.getx()+i1a[4]*tr1.gety()+i1a[5]*tr1.getz(),         i1a[6]*tr1.getx()+i1a[7]*tr1.gety()+i1a[8]*tr1.getz());     //(i1ir1n x r)     i1ir1n.crossproduct(r1);     //compute (inot*(r x n) x r)dot(n) + 1/m1     float jrdenom = i1ir1n.dotproduct(n) + (1.0f/m1);      //jr magnitude of impulse     float jr = jrnumer/jrdenom;     vector v1 = rectangle1->getlinvel();     vector v1prime = v1 - n*(jr/m1);      //calculate angular velocity     vector r1 = r1;     //r1 = (r x n)     r1.crossproduct(n);     //ir1 = (inot * r1)     vector ir1(i1a[0]*r1.getx()+i1a[1]*r1.gety()+i1a[2]*r1.getz(),         i1a[3]*r1.getx()+i1a[4]*r1.gety()+i1a[5]*r1.getz(),         i1a[6]*r1.getx()+i1a[7]*r1.gety()+i1a[8]*r1.getz());      vector w1 = rectangle1->getangvel();     //w1' = w1 - jr*inot     vector w1prime = w1 - ir1*jr;