Index: include/ode/collision.h
===================================================================
--- include/ode/collision.h	(revision 1641)
+++ include/ode/collision.h	(working copy)
@@ -936,6 +936,10 @@
 ODE_API void dGeomCapsuleGetParams (dGeomID ccylinder, dReal *radius, dReal *length);
 ODE_API dReal dGeomCapsulePointDepth (dGeomID ccylinder, dReal x, dReal y, dReal z);
 
+#define ODE_CAPSULE_CAN_SET_AXIS	//for testing purposes
+ODE_API void dGeomCapsuleSetAxis (dGeomID ccylinder, int axis);
+ODE_API int dGeomCapsuleGetAxis (dGeomID ccylinder);
+
 // For now we want to have a backwards compatible C-API, note: C++ API is not.
 #define dCreateCCylinder dCreateCapsule
 #define dGeomCCylinderSetParams dGeomCapsuleSetParams
Index: include/ode/odecpp_collision.h
===================================================================
--- include/ode/odecpp_collision.h	(revision 1641)
+++ include/ode/odecpp_collision.h	(working copy)
@@ -301,6 +301,11 @@
     { dGeomCapsuleSetParams (_id, radius, length); }
   void getParams (dReal *radius, dReal *length) const
     { dGeomCapsuleGetParams (_id,radius,length); }
+
+  void setAxis (int axis)
+    { dGeomCapsuleSetAxis (_id, axis); }
+  int getAxis () const
+    { return dGeomCapsuleGetAxis (_id); }  
 };
 
 
Index: ode/src/collision_std.h
===================================================================
--- ode/src/collision_std.h	(revision 1641)
+++ ode/src/collision_std.h	(working copy)
@@ -112,7 +112,8 @@
 
 
 struct dxCapsule : public dxGeom {
-  dReal radius,lz;	// radius, length along z axis
+  dReal radius,length;
+  int axis;
   dxCapsule (dSpaceID space, dReal _radius, dReal _length);
   void computeAABB();
 };
Index: ode/src/capsule.cpp
===================================================================
--- ode/src/capsule.cpp	(revision 1641)
+++ ode/src/capsule.cpp	(working copy)
@@ -46,12 +46,13 @@
 // capped cylinder public API
 
 dxCapsule::dxCapsule (dSpaceID space, dReal _radius, dReal _length) :
-  dxGeom (space,1)
+  dxGeom (space,1),
+  axis(2)
 {
   dAASSERT (_radius >= 0 && _length >= 0);
   type = dCapsuleClass;
   radius = _radius;
-  lz = _length;
+  length = _length;
   updateZeroSizedFlag(!_radius/* || !_length -- zero length capsule is not a zero sized capsule*/);
 }
 
@@ -61,9 +62,9 @@
   const dMatrix3& R = final_posr->R;
   const dVector3& pos = final_posr->pos;
   
-  dReal xrange = dFabs(R[2]  * lz) * REAL(0.5) + radius;
-  dReal yrange = dFabs(R[6]  * lz) * REAL(0.5) + radius;
-  dReal zrange = dFabs(R[10] * lz) * REAL(0.5) + radius;
+  dReal xrange = dFabs(R[axis]   * length) * REAL(0.5) + radius;
+  dReal yrange = dFabs(R[axis+4] * length) * REAL(0.5) + radius;
+  dReal zrange = dFabs(R[axis+8] * length) * REAL(0.5) + radius;
   aabb[0] = pos[0] - xrange;
   aabb[1] = pos[0] + xrange;
   aabb[2] = pos[1] - yrange;
@@ -85,21 +86,19 @@
   dAASSERT (radius >= 0 && length >= 0);
   dxCapsule *c = (dxCapsule*) g;
   c->radius = radius;
-  c->lz = length;
+  c->length = length;
   c->updateZeroSizedFlag(!radius/* || !length -- zero length capsule is not a zero sized capsule*/);
   dGeomMoved (g);
 }
 
-
 void dGeomCapsuleGetParams (dGeomID g, dReal *radius, dReal *length)
 {
   dUASSERT (g && g->type == dCapsuleClass,"argument not a ccylinder");
   dxCapsule *c = (dxCapsule*) g;
   *radius = c->radius;
-  *length = c->lz;
+  *length = c->length;
 }
 
-
 dReal dGeomCapsulePointDepth (dGeomID g, dReal x, dReal y, dReal z)
 {
   dUASSERT (g && g->type == dCapsuleClass,"argument not a ccylinder");
@@ -113,18 +112,32 @@
   a[0] = x - pos[0];
   a[1] = y - pos[1];
   a[2] = z - pos[2];
-  dReal beta = dDOT14(a,R+2);
-  dReal lz2 = c->lz*REAL(0.5);
-  if (beta < -lz2) beta = -lz2;
-  else if (beta > lz2) beta = lz2;
-  a[0] = c->final_posr->pos[0] + beta*R[0*4+2];
-  a[1] = c->final_posr->pos[1] + beta*R[1*4+2];
-  a[2] = c->final_posr->pos[2] + beta*R[2*4+2];
+  dReal beta = dDOT14(a,R+c->axis);
+  dReal length2 = c->length*REAL(0.5);
+  if (beta < -length2) beta = -length2;
+  else if (beta > length2) beta = length2;
+  a[0] = c->final_posr->pos[0] + beta*R[c->axis];
+  a[1] = c->final_posr->pos[1] + beta*R[c->axis+4];
+  a[2] = c->final_posr->pos[2] + beta*R[c->axis+8];
   return c->radius -
     dSqrt ((x-a[0])*(x-a[0]) + (y-a[1])*(y-a[1]) + (z-a[2])*(z-a[2]));
 }
 
+void dGeomCapsuleSetAxis (dGeomID g, int axis)
+{
+  dUASSERT (g && g->type == dCapsuleClass,"argument not a capsule");
+  dAASSERT (axis >= 0 && axis <=2);
+  dxCapsule *c = (dxCapsule*) g;
+  c->axis = axis;
+  dGeomMoved (g);
+}
 
+int dGeomCapsuleGetAxis (dGeomID g)
+{
+  dUASSERT (g && g->type == dCapsuleClass,"argument not a capsule");
+  dxCapsule *c = (dxCapsule*) g;
+  return c->axis;
+}
 
 int dCollideCapsuleSphere (dxGeom *o1, dxGeom *o2, int flags,
 			     dContactGeom *contact, int skip)
@@ -144,18 +157,18 @@
 
   // find the point on the cylinder axis that is closest to the sphere
   dReal alpha = 
-    o1->final_posr->R[2]  * (o2->final_posr->pos[0] - o1->final_posr->pos[0]) +
-    o1->final_posr->R[6]  * (o2->final_posr->pos[1] - o1->final_posr->pos[1]) +
-    o1->final_posr->R[10] * (o2->final_posr->pos[2] - o1->final_posr->pos[2]);
-  dReal lz2 = ccyl->lz * REAL(0.5);
-  if (alpha > lz2) alpha = lz2;
-  if (alpha < -lz2) alpha = -lz2;
+    o1->final_posr->R[ccyl->axis]   * (o2->final_posr->pos[0] - o1->final_posr->pos[0]) +
+    o1->final_posr->R[ccyl->axis+4] * (o2->final_posr->pos[1] - o1->final_posr->pos[1]) +
+    o1->final_posr->R[ccyl->axis+8] * (o2->final_posr->pos[2] - o1->final_posr->pos[2]);
+  dReal length2 = ccyl->length * REAL(0.5);
+  if (alpha > length2) alpha = length2;
+  if (alpha < -length2) alpha = -length2;
 
   // collide the spheres
   dVector3 p;
-  p[0] = o1->final_posr->pos[0] + alpha * o1->final_posr->R[2];
-  p[1] = o1->final_posr->pos[1] + alpha * o1->final_posr->R[6];
-  p[2] = o1->final_posr->pos[2] + alpha * o1->final_posr->R[10];
+  p[0] = o1->final_posr->pos[0] + alpha * o1->final_posr->R[ccyl->axis];
+  p[1] = o1->final_posr->pos[1] + alpha * o1->final_posr->R[ccyl->axis+4];
+  p[2] = o1->final_posr->pos[2] + alpha * o1->final_posr->R[ccyl->axis+8];
   return dCollideSpheres (p,ccyl->radius,o2->final_posr->pos,sphere->radius,contact);
 }
 
@@ -178,13 +191,13 @@
 
   // get p1,p2 = cylinder axis endpoints, get radius
   dVector3 p1,p2;
-  dReal clen = cyl->lz * REAL(0.5);
-  p1[0] = o1->final_posr->pos[0] + clen * o1->final_posr->R[2];
-  p1[1] = o1->final_posr->pos[1] + clen * o1->final_posr->R[6];
-  p1[2] = o1->final_posr->pos[2] + clen * o1->final_posr->R[10];
-  p2[0] = o1->final_posr->pos[0] - clen * o1->final_posr->R[2];
-  p2[1] = o1->final_posr->pos[1] - clen * o1->final_posr->R[6];
-  p2[2] = o1->final_posr->pos[2] - clen * o1->final_posr->R[10];
+  dReal clen = cyl->length * REAL(0.5);
+  p1[0] = o1->final_posr->pos[0] + clen * o1->final_posr->R[cyl->axis];
+  p1[1] = o1->final_posr->pos[1] + clen * o1->final_posr->R[cyl->axis+4];
+  p1[2] = o1->final_posr->pos[2] + clen * o1->final_posr->R[cyl->axis+8];
+  p2[0] = o1->final_posr->pos[0] - clen * o1->final_posr->R[cyl->axis];
+  p2[1] = o1->final_posr->pos[1] - clen * o1->final_posr->R[cyl->axis+4];
+  p2[2] = o1->final_posr->pos[2] - clen * o1->final_posr->R[cyl->axis+8];
   dReal radius = cyl->radius;
 
   // copy out box center, rotation matrix, and side array
@@ -221,17 +234,17 @@
   contact->side2 = -1;
 
   // copy out some variables, for convenience
-  dReal lz1 = cyl1->lz * REAL(0.5);
-  dReal lz2 = cyl2->lz * REAL(0.5);
+  dReal length1 = cyl1->length * REAL(0.5);
+  dReal length2 = cyl2->length * REAL(0.5);
   dReal *pos1 = o1->final_posr->pos;
   dReal *pos2 = o2->final_posr->pos;
   dReal axis1[3],axis2[3];
-  axis1[0] = o1->final_posr->R[2];
-  axis1[1] = o1->final_posr->R[6];
-  axis1[2] = o1->final_posr->R[10];
-  axis2[0] = o2->final_posr->R[2];
-  axis2[1] = o2->final_posr->R[6];
-  axis2[2] = o2->final_posr->R[10];
+  axis1[0] = o1->final_posr->R[cyl1->axis];
+  axis1[1] = o1->final_posr->R[cyl1->axis+4];
+  axis1[2] = o1->final_posr->R[cyl1->axis+8];
+  axis2[0] = o2->final_posr->R[cyl2->axis];
+  axis2[1] = o2->final_posr->R[cyl2->axis+4];
+  axis2[2] = o2->final_posr->R[cyl2->axis+8];
 
   // if the cylinder axes are close to parallel, we'll try to detect up to
   // two contact points along the body of the cylinder. if we can't find any
@@ -257,10 +270,10 @@
     dReal q[3];
     for (i=0; i<3; i++) q[i] = pos1[i]-pos2[i];
     dReal k = dDOT (axis1,q);
-    dReal a1lo = -lz1;
-    dReal a1hi = lz1;
-    dReal a2lo = -lz2 - k;
-    dReal a2hi = lz2 - k;
+    dReal a1lo = -length1;
+    dReal a1hi = length1;
+    dReal a2lo = -length2 - k;
+    dReal a2hi = length2 - k;
     dReal lo = (a1lo > a2lo) ? a1lo : a2lo;
     dReal hi = (a1hi < a2hi) ? a1hi : a2hi;
     if (lo <= hi) {
@@ -301,18 +314,18 @@
 	  
   // use the closest point algorithm
   dVector3 a1,a2,b1,b2;
-  a1[0] = o1->final_posr->pos[0] + axis1[0]*lz1;
-  a1[1] = o1->final_posr->pos[1] + axis1[1]*lz1;
-  a1[2] = o1->final_posr->pos[2] + axis1[2]*lz1;
-  a2[0] = o1->final_posr->pos[0] - axis1[0]*lz1;
-  a2[1] = o1->final_posr->pos[1] - axis1[1]*lz1;
-  a2[2] = o1->final_posr->pos[2] - axis1[2]*lz1;
-  b1[0] = o2->final_posr->pos[0] + axis2[0]*lz2;
-  b1[1] = o2->final_posr->pos[1] + axis2[1]*lz2;
-  b1[2] = o2->final_posr->pos[2] + axis2[2]*lz2;
-  b2[0] = o2->final_posr->pos[0] - axis2[0]*lz2;
-  b2[1] = o2->final_posr->pos[1] - axis2[1]*lz2;
-  b2[2] = o2->final_posr->pos[2] - axis2[2]*lz2;
+  a1[0] = o1->final_posr->pos[0] + axis1[0]*length1;
+  a1[1] = o1->final_posr->pos[1] + axis1[1]*length1;
+  a1[2] = o1->final_posr->pos[2] + axis1[2]*length1;
+  a2[0] = o1->final_posr->pos[0] - axis1[0]*length1;
+  a2[1] = o1->final_posr->pos[1] - axis1[1]*length1;
+  a2[2] = o1->final_posr->pos[2] - axis1[2]*length1;
+  b1[0] = o2->final_posr->pos[0] + axis2[0]*length2;
+  b1[1] = o2->final_posr->pos[1] + axis2[1]*length2;
+  b1[2] = o2->final_posr->pos[2] + axis2[2]*length2;
+  b2[0] = o2->final_posr->pos[0] - axis2[0]*length2;
+  b2[1] = o2->final_posr->pos[1] - axis2[1]*length2;
+  b2[2] = o2->final_posr->pos[2] - axis2[2]*length2;
 
   dClosestLineSegmentPoints (a1,a2,b1,b2,sphere1,sphere2);
   return dCollideSpheres (sphere1,cyl1->radius,sphere2,cyl2->radius,contact);
@@ -331,11 +344,11 @@
   dxPlane *plane = (dxPlane*) o2;
 
   // collide the deepest capping sphere with the plane
-  dReal sign = (dDOT14 (plane->p,o1->final_posr->R+2) > 0) ? REAL(-1.0) : REAL(1.0);
+  dReal sign = (dDOT14 (plane->p,o1->final_posr->R+ccyl->axis) > 0) ? REAL(-1.0) : REAL(1.0);
   dVector3 p;
-  p[0] = o1->final_posr->pos[0] + o1->final_posr->R[2]  * ccyl->lz * REAL(0.5) * sign;
-  p[1] = o1->final_posr->pos[1] + o1->final_posr->R[6]  * ccyl->lz * REAL(0.5) * sign;
-  p[2] = o1->final_posr->pos[2] + o1->final_posr->R[10] * ccyl->lz * REAL(0.5) * sign;
+  p[0] = o1->final_posr->pos[0] + o1->final_posr->R[ccyl->axis]   * ccyl->length * REAL(0.5) * sign;
+  p[1] = o1->final_posr->pos[1] + o1->final_posr->R[ccyl->axis+4] * ccyl->length * REAL(0.5) * sign;
+  p[2] = o1->final_posr->pos[2] + o1->final_posr->R[ccyl->axis+8] * ccyl->length * REAL(0.5) * sign;
 
   dReal k = dDOT (p,plane->p);
   dReal depth = plane->p[3] - k + ccyl->radius;
@@ -351,9 +364,9 @@
   int ncontacts = 1;
   if ((flags & NUMC_MASK) >= 2) {
     // collide the other capping sphere with the plane
-    p[0] = o1->final_posr->pos[0] - o1->final_posr->R[2]  * ccyl->lz * REAL(0.5) * sign;
-    p[1] = o1->final_posr->pos[1] - o1->final_posr->R[6]  * ccyl->lz * REAL(0.5) * sign;
-    p[2] = o1->final_posr->pos[2] - o1->final_posr->R[10] * ccyl->lz * REAL(0.5) * sign;
+    p[0] = o1->final_posr->pos[0] - o1->final_posr->R[ccyl->axis]   * ccyl->length * REAL(0.5) * sign;
+    p[1] = o1->final_posr->pos[1] - o1->final_posr->R[ccyl->axis+4] * ccyl->length * REAL(0.5) * sign;
+    p[2] = o1->final_posr->pos[2] - o1->final_posr->R[ccyl->axis+8] * ccyl->length * REAL(0.5) * sign;
 
     k = dDOT (p,plane->p);
     depth = plane->p[3] - k + ccyl->radius;
Index: ode/src/collision_trimesh_ccylinder.cpp
===================================================================
--- ode/src/collision_trimesh_ccylinder.cpp	(revision 1641)
+++ ode/src/collision_trimesh_ccylinder.cpp	(working copy)
@@ -188,9 +188,11 @@
 	dMatrix3  m_mCapsuleRotation;
 	dVector3   m_vCapsulePosition;
 	dVector3   m_vCapsuleAxis;
+	
 	// static data
 	dReal      m_vCapsuleRadius;
 	dReal      m_fCapsuleSize;
+	int        m_iCapsuleAxisId;
 
 	// mesh data
 	// dMatrix4  mHullDstPl;
@@ -215,10 +217,7 @@
 	int m_iStride;
 };
 
-// Capsule lie on axis number 3 = (Z axis)
-static const int nCAPSULE_AXIS = 2;
 
-
 #if OPTIMIZE_CONTACTS
 
 // Use to classify contacts to be "near" in position
@@ -945,11 +944,12 @@
 
 	const dVector3* pDst = (const dVector3*)dGeomGetPosition(Capsule);
 	memcpy(m_vCapsulePosition, pDst, sizeof(dVector3));
+	
+	m_iCapsuleAxisId=dGeomCapsuleGetAxis(Capsule);
+	m_vCapsuleAxis[0] = m_mCapsuleRotation[m_iCapsuleAxisId];
+	m_vCapsuleAxis[1] = m_mCapsuleRotation[m_iCapsuleAxisId+4];
+	m_vCapsuleAxis[2] = m_mCapsuleRotation[m_iCapsuleAxisId+8];
 
-	m_vCapsuleAxis[0] = m_mCapsuleRotation[0*4 + nCAPSULE_AXIS];
-	m_vCapsuleAxis[1] = m_mCapsuleRotation[1*4 + nCAPSULE_AXIS];
-	m_vCapsuleAxis[2] = m_mCapsuleRotation[2*4 + nCAPSULE_AXIS];
-
 	// Get size of Capsule
 	dGeomCapsuleGetParams(Capsule, &m_vCapsuleRadius, &m_fCapsuleSize);
 	m_fCapsuleSize += 2*m_vCapsuleRadius;
@@ -1008,7 +1008,8 @@
 	const dVector3 &vCapsulePosition = cData.m_vCapsulePosition;
 
 	Point cCenter(/*(float)*/ vCapsulePosition[0], /*(float)*/ vCapsulePosition[1], /*(float)*/ vCapsulePosition[2]);
-	Point cExtents(/*(float)*/ cData.m_vCapsuleRadius, /*(float)*/ cData.m_vCapsuleRadius,/*(float)*/ cData.m_fCapsuleSize/2);
+	Point cExtents(/*(float)*/ cData.m_vCapsuleRadius, /*(float)*/ cData.m_vCapsuleRadius,/*(float)*/ cData.m_vCapsuleRadius);
+	cExtents[cData.m_iCapsuleAxisId]=cData.m_fCapsuleSize/2;
 
 	Matrix3x3 obbRot;
 
@@ -1139,7 +1140,6 @@
 
 // GIMPACT version
 #if dTRIMESH_GIMPACT
-#define nCAPSULE_AXIS 2
 // capsule - trimesh  By francisco leon
 int dCollideCCTL(dxGeom *o1, dxGeom *o2, int flags, dContactGeom *contact, int skip)
 {
@@ -1162,9 +1162,10 @@
 	dVector3* pDst = (dVector3*)dGeomGetPosition(gCylinder);
 	memcpy(vCapsulePosition,pDst,sizeof(dVector3));
 	//Axis
-	vCapsuleAxis[0] = mCapsuleRotation[0*4 + nCAPSULE_AXIS];
-	vCapsuleAxis[1] = mCapsuleRotation[1*4 + nCAPSULE_AXIS];
-	vCapsuleAxis[2] = mCapsuleRotation[2*4 + nCAPSULE_AXIS];
+    int capsule_axis=dGeomCapsuleGetAxis(gCylinder);
+	vCapsuleAxis[0] = mCapsuleRotation[capsule_axis];
+	vCapsuleAxis[1] = mCapsuleRotation[capsule_axis+4];
+	vCapsuleAxis[2] = mCapsuleRotation[capsule_axis+8];
 	// Get size of CCylinder
 	dGeomCCylinderGetParams(gCylinder,&vCapsuleRadius,&fCapsuleSize);
 	fCapsuleSize*=0.5f;
Index: ode/src/ray.cpp
===================================================================
--- ode/src/ray.cpp	(revision 1641)
+++ ode/src/ray.cpp	(working copy)
@@ -364,7 +364,8 @@
   contact->side1 = -1;
   contact->side2 = -1;
   
-  dReal lz2 = ccyl->lz * REAL(0.5);
+  dReal lz2 = ccyl->length * REAL(0.5);
+  int axis = ccyl->axis;
 
   // compute some useful info
   dVector3 cs,q,r;
@@ -372,10 +373,10 @@
   cs[0] = ray->final_posr->pos[0] - ccyl->final_posr->pos[0];
   cs[1] = ray->final_posr->pos[1] - ccyl->final_posr->pos[1];
   cs[2] = ray->final_posr->pos[2] - ccyl->final_posr->pos[2];
-  k = dDOT41(ccyl->final_posr->R+2,cs);	// position of ray start along ccyl axis
-  q[0] = k*ccyl->final_posr->R[0*4+2] - cs[0];
-  q[1] = k*ccyl->final_posr->R[1*4+2] - cs[1];
-  q[2] = k*ccyl->final_posr->R[2*4+2] - cs[2];
+  k = dDOT41(ccyl->final_posr->R+axis,cs);	// position of ray start along ccyl axis
+  q[0] = k*ccyl->final_posr->R[axis] - cs[0];
+  q[1] = k*ccyl->final_posr->R[axis+4] - cs[1];
+  q[2] = k*ccyl->final_posr->R[axis+8] - cs[2];
   C = dDOT(q,q) - ccyl->radius*ccyl->radius;
   // if C < 0 then ray start position within infinite extension of cylinder
 
@@ -384,9 +385,9 @@
   if (C < 0) {
     if (k < -lz2) k = -lz2;
     else if (k > lz2) k = lz2;
-    r[0] = ccyl->final_posr->pos[0] + k*ccyl->final_posr->R[0*4+2];
-    r[1] = ccyl->final_posr->pos[1] + k*ccyl->final_posr->R[1*4+2];
-    r[2] = ccyl->final_posr->pos[2] + k*ccyl->final_posr->R[2*4+2];
+    r[0] = ccyl->final_posr->pos[0] + k*ccyl->final_posr->R[axis];
+    r[1] = ccyl->final_posr->pos[1] + k*ccyl->final_posr->R[axis+4];
+    r[2] = ccyl->final_posr->pos[2] + k*ccyl->final_posr->R[axis+8];
     if ((ray->final_posr->pos[0]-r[0])*(ray->final_posr->pos[0]-r[0]) +
 	(ray->final_posr->pos[1]-r[1])*(ray->final_posr->pos[1]-r[1]) +
 	(ray->final_posr->pos[2]-r[2])*(ray->final_posr->pos[2]-r[2]) < ccyl->radius*ccyl->radius) {
@@ -402,10 +403,10 @@
     if (k < 0) k = -lz2; else k = lz2;
   }
   else {
-    dReal uv = dDOT44(ccyl->final_posr->R+2,ray->final_posr->R+2);
-    r[0] = uv*ccyl->final_posr->R[0*4+2] - ray->final_posr->R[0*4+2];
-    r[1] = uv*ccyl->final_posr->R[1*4+2] - ray->final_posr->R[1*4+2];
-    r[2] = uv*ccyl->final_posr->R[2*4+2] - ray->final_posr->R[2*4+2];
+    dReal uv = dDOT44(ccyl->final_posr->R+axis,ray->final_posr->R+axis);
+    r[0] = uv*ccyl->final_posr->R[axis]   - ray->final_posr->R[axis];
+    r[1] = uv*ccyl->final_posr->R[axis+4] - ray->final_posr->R[axis+4];
+    r[2] = uv*ccyl->final_posr->R[axis+8] - ray->final_posr->R[axis+8];
     dReal A = dDOT(r,r);
     dReal B = 2*dDOT(q,r);
     k = B*B-4*A*C;
@@ -428,21 +429,21 @@
 
       // the ray intersects the infinite cylinder. check to see if the
       // intersection point is between the caps
-      contact->pos[0] = ray->final_posr->pos[0] + alpha*ray->final_posr->R[0*4+2];
-      contact->pos[1] = ray->final_posr->pos[1] + alpha*ray->final_posr->R[1*4+2];
-      contact->pos[2] = ray->final_posr->pos[2] + alpha*ray->final_posr->R[2*4+2];
+      contact->pos[0] = ray->final_posr->pos[0] + alpha*ray->final_posr->R[axis];
+      contact->pos[1] = ray->final_posr->pos[1] + alpha*ray->final_posr->R[axis+4];
+      contact->pos[2] = ray->final_posr->pos[2] + alpha*ray->final_posr->R[axis+8];
       q[0] = contact->pos[0] - ccyl->final_posr->pos[0];
       q[1] = contact->pos[1] - ccyl->final_posr->pos[1];
       q[2] = contact->pos[2] - ccyl->final_posr->pos[2];
-      k = dDOT14(q,ccyl->final_posr->R+2);
+      k = dDOT14(q,ccyl->final_posr->R+axis);
       dReal nsign = inside_ccyl ? REAL(-1.0) : REAL(1.0);
       if (k >= -lz2 && k <= lz2) {
 	contact->normal[0] = nsign * (contact->pos[0] -
-				      (ccyl->final_posr->pos[0] + k*ccyl->final_posr->R[0*4+2]));
+				      (ccyl->final_posr->pos[0] + k*ccyl->final_posr->R[axis]));
 	contact->normal[1] = nsign * (contact->pos[1] -
-				      (ccyl->final_posr->pos[1] + k*ccyl->final_posr->R[1*4+2]));
+				      (ccyl->final_posr->pos[1] + k*ccyl->final_posr->R[axis+4]));
 	contact->normal[2] = nsign * (contact->pos[2] -
-				      (ccyl->final_posr->pos[2] + k*ccyl->final_posr->R[2*4+2]));
+				      (ccyl->final_posr->pos[2] + k*ccyl->final_posr->R[axis+8]));
 	dNormalize3 (contact->normal);
 	contact->depth = alpha;
 	return 1;
@@ -456,9 +457,9 @@
 
   // check for ray intersection with the caps. k must indicate the cap
   // position to check
-  q[0] = ccyl->final_posr->pos[0] + k*ccyl->final_posr->R[0*4+2];
-  q[1] = ccyl->final_posr->pos[1] + k*ccyl->final_posr->R[1*4+2];
-  q[2] = ccyl->final_posr->pos[2] + k*ccyl->final_posr->R[2*4+2];
+  q[0] = ccyl->final_posr->pos[0] + k*ccyl->final_posr->R[axis];
+  q[1] = ccyl->final_posr->pos[1] + k*ccyl->final_posr->R[axis+4];
+  q[2] = ccyl->final_posr->pos[2] + k*ccyl->final_posr->R[axis+8];
   return ray_sphere_helper (ray,q,ccyl->radius,contact, inside_ccyl);
 }
 
