comments done !

src/rotation.c

@@ -34,17 +34,43 @@ struct pdb_values* rotate_sophie(struct pdb_values* pdb, struct pdb_values* pdbR
     float x2i = 0, y2i = 0, z2i = 0;
     float x3i = 0, y3i = 0, z3i = 0;
 
-    // cosinuses and sinuses of d_angle
+    /* cosinuses and sinuses of the rotation angles dalpha, dbeta and dgamma
+     * In this function, alpha = dalpha, beta = dbeta, gamma = dgamma
+    */
     float c_a=cos(alpha),s_a=sin(alpha);
     float c_b=cos(beta),s_b=sin(beta);
     float c_g=cos(gamma),s_g=sin(gamma);
 
+    /* The angles alpha0, beta0, gamma0 corresponds to the angles for the axis 
+     * (origin,CA1 rotated by alpha), (origin, CA1 rotated by beta) and
+     * (origin, CA5 rotated by alpha then beta) respectively 
+     *
+     * We name these axis :
+     *  N : (origin,CA1 rotated by alpha) 
+     *  Z : (origin, CA5 rotated by alpha then beta) 
+     *
+     * The axis N therefore has the coordinates :
+     *  xN : c_a0
+     *  yN : s_a0
+     *  zN : 0
+     * This axis is fixed when performing the rotation beta
+     *
+     */
     float c_a0=cos(dock_res->alpha[nConf]),s_a0=sin(dock_res->alpha[nConf]);
     float c_b0=cos(dock_res->beta[nConf]), s_b0=sin(dock_res->beta[nConf]);
 
-    /* These values correspond to the final rotation axis obtained after rotating
-     * the Y axis by an angle of alpha0 then beta0*/
-    float rx = 0, ry = 0, rz = 0;
+    float xN = c_a0;
+    float yN = s_a0;
+    float zN = 0;
+    
+    /* The coordinates of the unit vector for the axis Z can be computed as such :
+     *  xZ : -s_a0 * c_b0
+     *  yZ :  c_a0 * c_b0
+     *  zZ :  s_b0
+     * */
+    float xZ = -s_a0*c_b0;
+    float yZ = c_a0*c_b0;
+    float zZ = s_b0;
 
     if(newPDB == NULL){
         newPDB = allocate_pdb(pdb->nbRes);
@@ -71,21 +97,14 @@ struct pdb_values* rotate_sophie(struct pdb_values* pdb, struct pdb_values* pdbR
             z1i=zi;
 
             /* beta rotation of the residue */
-            x2i=(c_a0*c_a0+(1.0-c_a0*c_a0)*c_b)*x1i+c_a0*s_a0*(1.0-c_b)*y1i+s_a0*s_b*z1i;
-            y2i=c_a0*s_a0*(1.0-c_b)*x1i+(s_a0*s_a0+(1.0-s_a0*s_a0)*c_b)*y1i-c_a0*s_b*z1i;
-            z2i=-s_a0*s_b*x1i+c_a0*s_b*y1i+c_b*z1i;
-            
-            /* We compute the value for the last rotation axis
-             * from the angles alpha0, beta0 and gamma0 
-             */
-            rx = -s_a0*c_b0;
-            ry = c_a0*c_b0;
-            rz = s_b0;
+            x2i = (xN*xN+(1.0-xN*xN)*c_b)*x1i + xN*yN*(1.0-c_b)*y1i         + yN*s_b*z1i;
+            y2i = xN*yN*(1.0-c_b)*x1i         + (yN*yN+(1.0-yN*yN)*c_b)*y1i - xN*s_b*z1i;
+            z2i = -yN*s_b*x1i                 + xN*s_b*y1i                  + c_b*z1i;
 
             // gamma rotation of the residue
-            x3i=(rx*rx+(1.0-rx*rx)*c_g)*x2i+(rx*ry*(1.0-c_g)-rz*s_g)*y2i+(rx*rz*(1.0-c_g)+ry*s_g)*z2i;
-            y3i=(rx*ry*(1.0-c_g)+rz*s_g)*x2i+(ry*ry+(1.0-ry*ry)*c_g)*y2i+(ry*rz*(1.0-c_g)-rx*s_g)*z2i;
-            z3i=(rx*rz*(1.0-c_g)-ry*s_g)*x2i+(ry*rz*(1.0-c_g)+rx*s_g)*y2i+(rz*rz+(1.0-rz*rz)*c_g)*z2i;
+            x3i = (xZ*xZ+(1.0-xZ*xZ)*c_g) *x2i + (xZ*yZ*(1.0-c_g)-zZ*s_g)*y2i + (xZ*zZ*(1.0-c_g)+yZ*s_g)*z2i;
+            y3i = (xZ*yZ*(1.0-c_g)+zZ*s_g)*x2i + (yZ*yZ+(1.0-yZ*yZ)*c_g) *y2i + (yZ*zZ*(1.0-c_g)-xZ*s_g)*z2i;
+            z3i = (xZ*zZ*(1.0-c_g)-yZ*s_g)*x2i + (yZ*zZ*(1.0-c_g)+xZ*s_g)*y2i + (zZ*zZ+(1.0-zZ*zZ)*c_g) *z2i;
 
 
             /* Now the residue has been rotated, and we need to translate it 
@@ -118,9 +137,11 @@ struct pdb_values* rotate_sophie(struct pdb_values* pdb, struct pdb_values* pdbR
 
 
 void getAnglesFromCA1andCA5_sophie(struct pdb_values* pdb, float* alpha, float* beta, float* gamma){
-    /* The goal of this function is to find the angles alpha, beta and gamma that correspond 
-     * to the angles alpha and beta of the first CA of the protein
-     * We then have to rotate the fifth CA by the angle alpha, then beta in order to compute the gamma angle
+    /* In this function, we first need to find the initial angles alpha and beta that correspond
+     * to the angles alpha and beta for the first CA of the protein. 
+     * Then, we need to apply a rotation alpha around the axis Z, then a rotation beta 
+     * around the axis X on the fifth atom CA of the protein.
+     * The coordinates of the then rotated atom CA 5 will give us the initial angle gamma 
      */
 
     *alpha = 0;