How to rotate an object around a local axis in OpenGL?

I am working on an ongoing project where I want to align the chain links so that it matches the contours of the Bezier curve. I am currently following these steps.

• Drawing a curve.
• Use the display list to create a single link in the chain.
• Use the FOR loop to repeatedly call a function that calculates the angle between two points of the curve, returns the angle and the axis around which the link should be rotated.
• Turn the corner "a" and move to a new position, place the link in a new position.

Edit: I must also say that the centers of two one-and-a-half tori must lie on the Bezier curve. I also know that the method that I use to draw the torus is tiring, later I will use TRIANGLE_FAN or QUAD_STRIP to draw the torus more efficiently.

While at first glance this logic looks like it displays the chain correctly, the final result is not what I imagined. Here's a picture of how the chain looks.

I read that before translating you need to translate the object to the origin? Would I just call glTranslate (0,0,0) and then follow step 4 from above?

I have included the appropriate code from what I have done so far, I would appreciate any suggestions that will help me work correctly with the code.

/* this function calculates the angle between two vectors oldPoint and new point contain the x,y,z coordinates of the two points,axisOfRot is used to return the x,y,z coordinates of the rotation axis*/
double getAxisAngle(pointType oldPoint[],
pointType newPoint[],pointType axisOfRot[]){

float tmpPoint[3];   

float normA = 0.0,normB = 0.0,AB = 0.0,angle=0.0;
int i;

axisOfRot->x= oldPoint->y * newPoint->z - oldPoint->z * newPoint->y;
axisOfRot->y= oldPoint->z * newPoint->x - oldPoint->x * newPoint->z;
axisOfRot->z= oldPoint->x * newPoint->y - oldPoint->y * newPoint->x; 


normA=sqrt(oldPoint->x * oldPoint->x + oldPoint->y * oldPoint->y + oldPoint->z * 
oldPoint->z);
normB=sqrt(newPoint->x * newPoint->x + newPoint->y * newPoint->y + newPoint->z *    
newPoint->z);

tmpPoint[0] =  oldPoint->x * newPoint->x;
tmpPoint[1] =  oldPoint->y * newPoint->y;
tmpPoint[2] =  oldPoint->z * newPoint->z;   

for(i=0;i<=2;i++)
 AB+=tmpPoint[i];  

AB /= (normA * normB);  


    return angle = (180/PI)*acos(AB);  
 }

 /* this function calculates and returns the next point on the curve give the 4 initial points for the curve, t is the tension of the curve */
    void bezierInterpolation(float t,pointType cPoints[],
    pointType newPoint[]){

newPoint->x = pow(1 - t, 3) * cPoints[0].x +3 * pow(1 - t , 2) * t * cPoints[1].x + 3 
* pow(1 - t, 1) * pow(t, 2) * cPoints[2].x + pow(t, 3) * cPoints[3].x;

newPoint->y = pow(1 - t, 3) * cPoints[0].y +3 * pow(1 - t , 2) * t * cPoints[1].y + 3 
* pow(1 - t, 1) * pow(t, 2) * cPoints[2].y + pow(t, 3) * cPoints[3].y;

newPoint->z = pow(1 - t, 3) * cPoints[0].z +3 * pow(1 - t , 2) * t * cPoints[1].z + 3 
* pow(1 - t, 1) * pow(t, 2) * cPoints[2].z + pow(t, 3) * cPoints[3].z;

}

/* the two lists below are used to create a single link in a chain, I realize that creating a half torus using cylinders is a bad idea, I will use GL_STRIP or TRIANGLE_FAN once I get the alignment right 
*/
torusList=glGenLists(1);   

glNewList(torusList,GL_COMPILE);
for (i=0; i<=180; i++)
{
  degInRad = i*DEG2RAD;
  glPushMatrix();
  glTranslatef(cos(degInRad)*radius,sin(degInRad)*radius,0);
  glRotated(90,1,0,0);
  gluCylinder(quadric,Diameter/2,Diameter/2,Height/5,10,10);
  glPopMatrix();      
}

 glEndList();

    /*! create a list for the link , 2 half torus and 2 columns */

 linkList = glGenLists(1);

 glNewList(linkList, GL_COMPILE); 
 glPushMatrix();
 glCallList(torusList); 
 glRotatef(90,1,0,0);
 glTranslatef(radius,0,0);
 gluCylinder(quadric, Diameter/2, Diameter/2, Height,10,10); 
 glTranslatef(-(radius*2),0,0);
 gluCylinder(quadric, Diameter/2, Diameter/2, Height,10,10); 
 glTranslatef(radius,0, Height);
 glRotatef(90,1,0,0);
 glCallList(torusList); 
 glPopMatrix();
 glEndList();

Finally, here is the code to create three links in the chain

t=0.031; 
bezierInterpolation(t,cPoints,newPoint);  
a=getAxisAngle(oldPoint,newPoint,axisOfRot);
 glTranslatef(newPoint->x,newPoint->y,newPoint->z); 
   glRotatef(a,axisOfRot->x,axisOfRot->y,axisOfRot->z);     
    glCallList(DLid);
   glRotatef(-a,axisOfRot->x,axisOfRot->y,axisOfRot->z); 
 glTranslatef(-newPoint->x,-newPoint->y,-newPoint->z);  

oldPoint[0]=newPoint[0];
bezierInterpolation(t+=GAP,cPoints,newPoint);  
a=getAxisAngle(oldPoint,newPoint,axisOfRot);
 glTranslatef(newPoint->x,newPoint->y,newPoint->z); 
  glRotatef(90,0,1,0);
   glRotatef(a,axisOfRot->x,axisOfRot->y,axisOfRot->z);     
    glCallList(DLid);
   glRotatef(-a,axisOfRot->x,axisOfRot->y,axisOfRot->z); 
  glRotatef(90,0,1,0);  
 glTranslatef(-newPoint->x,-newPoint->y,-newPoint->z);  

 oldPoint[0]=newPoint[0];
 bezierInterpolation(t+=GAP,cPoints,newPoint);    
  a=getAxisAngle(oldPoint,newPoint,axisOfRot); 
   glTranslatef(newPoint->x,newPoint->y,newPoint->z);       
    glRotatef(-a,axisOfRot->x,axisOfRot->y,axisOfRot->z); 
     glCallList(DLid);
    glRotatef(a,axisOfRot->x,axisOfRot->y,axisOfRot->z); 
   glTranslatef(-newPoint->x,-newPoint->y,newPoint->z);