# include <stdio.h>
# include <math.h>


/*
This is my function to rotate a helix. I think it also translates
it centering the helix on the origin. I think to remember also
that I still had to solve the problem of specular images or
something like that, and I suspect there is something in the
function which does it but not quite. Thank you for the new
alignment. I have now a whole set of matrices of neighbors
with different measures and some ideas about a flexible
desriptor of the environment at each position and how to
choose it, which I would like to discuss with you sometimes,
maybe between tennis games.
*/

/*
   vector = CAs of helix coming in 
   newvector = CAs of helix coming out 
   n = number of CAs
   newc = rotation matrix
*/
double transform(double vector[][3],int n,double newvector[][3],double newc[][3])
{
  int i,j;
  double	
    error=0.0,segno,
    A,B,C,K,Xx,Xy,Xz,Yx,Yy,Yz,Zx,Zy,Zz,distm,dist,
    mx,my,mz,x0,y0,z0,d1,d2,
    Sx=0.,Sy=0.,Sz=0.,Sxq=0.,Syq=0.,Szq=0.,Sxz=0.,Syz=0.;
  
  for(i=0;i<n;++i) {
    Sx+=vector[i][0];
    Sy+=vector[i][1];
    Sz+=vector[i][2];
  }
  mx = Sx/(double)n;
  my = Sy/(double)n;
  mz = Sz/(double)n;
  
  Sx = Sy = Sz = 0.0;
  
  for(i=0;i<n;++i) {
    vector[i][0]-=mx;
    vector[i][1]-=my;
    vector[i][2]-=mz;
    Sx+=vector[i][0];
    Sy+=vector[i][1];
    Sz+=vector[i][2];
    Sxq+=vector[i][0]*vector[i][0];
    Syq+=vector[i][1]*vector[i][1];
    Szq+=vector[i][2]*vector[i][2];
    Sxz+=vector[i][0]*vector[i][2];
    Syz+=vector[i][1]*vector[i][2];
  }
  A = (Szq-Sxq+sqrt((Sxq-Szq)*(Sxq-Szq)+4.*Sxz*Sxz))/(2.*Sxz);
  C = (Szq-Sxq-sqrt((Sxq-Szq)*(Sxq-Szq)+4.*Sxz*Sxz))/(2.*Sxz);
  d1 = (Szq+A*A*Sxq-2.*A*Sxz)/((1.+A*A)*(1.+A*A));
  d2 = (Szq+C*C*Sxq-2.*C*Sxz)/((1.+C*C)*(1.+C*C));
  if(d1<d2) A = 1./A;
  else A = 1./C;
  
  B = (Szq-Syq+sqrt((Syq-Szq)*(Syq-Szq)+4.*Syz*Syz))/(2.*Syz);
  C = (Szq-Syq-sqrt((Syq-Szq)*(Syq-Szq)+4.*Syz*Syz))/(2.*Syz);
  d1 = (Szq+B*B*Syq-2.*B*Syz)/((1.+B*B)*(1.+B*B));
  d2 = (Szq+C*C*Sxq-2.*C*Sxz)/((1.+C*C)*(1.+C*C));
  if(d1<d2) B = 1./B;
  else B = 1./C;
  
  K=sqrt(1.+A*A+B*B);

  Xz = newc[0][2] = A/K;
  Yz = newc[1][2] = B/K;
  Zz = newc[2][2] = 1./K;
  
  for(i=0;i<n;++i)
    {
      segno = newvector[i][2] = Xz*vector[i][0]+Yz*vector[i][1]+Zz*vector[i][2];
      if(segno>0. && !i)
	{
	  newvector[i][2] = -newvector[i][2];
	  Xz = -Xz;
	  Yz = -Yz;
	  Zz = -Zz;
	}
    }
  
  A = -(Xz+Yz)/Zz;
  K=sqrt(1.+1.+A*A);
  
  Xx = newc[0][0] = 1./K;
  Yx = newc[1][0] = 1./K;
  Zx = newc[2][0] = A/K;
  
  
  for(i=0;i<n;++i)
    newvector[i][0] = Xx*vector[i][0]+Yx*vector[i][1]+Zx*vector[i][2];
  
  
  A = (Xz*Zx-Xx*Zz)/(Yx*Zz-Yz*Zx);
  B = -(Xz+A*Yz)/Zz;
  K=sqrt(1.+A*A+B*B);
  
  Xy = newc[0][1] = 1./K;
  Yy = newc[1][1] = A/K;
  Zy = newc[2][1] = B/K;
  
  for(i=0;i<n;++i)
    {
      newvector[i][1] = Xy*vector[i][0]+Yy*vector[i][1]+Zy*vector[i][2];
      /*   fprintf(stdout,"%6.2f %6.2f %6.2f\n",newvector[i][0],newvector[i][1],newvector[i][2]); */
    }
  
  
  distm=0.0;
  for(i=0;i<n;++i)
    {
      dist = 0.0;
      for(j=0;j<2;++j) dist += newvector[i][j]*newvector[i][j];
      distm += sqrt(dist);
    }
  distm = distm/(double)n;
  
  return(distm);
}