Getting a point on a sphere to sit at the north pole with all other points rotating the same amount

[adrem7]

Say I have 10 points all on a sphere (randomly distributed) and I want to rotate the entire system to make sure one point sits at the north pole. How would I do this?

I went about it by looking at 3D rotation matrices:

http://en.wikipedia.org/wiki/Rotation_matrix

I rotate my point around the x-axis until the y component is zero and then rotate around the y-axis until the x component is zero. This should leave the point in question at either the north or south pole right?

My code looks like so:

#include <stdio.h> 
#include <string.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <time.h>
#include <stdlib.h>
#include <sstream>
using namespace std;
#define PI 3.14159265358979323846

int main()
{

  int a,b,c,f,i,j,k,m,n,s;
  double p,Time,Averagetime,Energy,energy,Distance,Length,DotProdForce,Forcemagnitude,
         ForceMagnitude[101],Force[101][4],E[1000001],En[501],x[101][4],y[101][4],
         Dist[101][101],Epsilon,z[101],theta,phi;

    /*  set the number of points */
    n=10;

    /* check that there are no more than 100 points */
    if(n>100){
      cout << n << " is too many points for me :-( \n";
      exit(0);
    }
  
    /* reset the random number generator */
    srand((unsigned)time(0));  
  
    for (i=1;i<=n;i++){
      x[i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
    
      Length=sqrt(pow(x[i][1],2)+pow(x[i][2],2)+pow(x[i][3],2));
    
      for (k=1;k<=3;k++){
        x[i][k]=x[i][k]/Length;
      }
    }

    /* calculate the energy */
    Energy=0.0;
  
    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
		            +pow(x[i][3]-x[j][3],2));
      
        Energy=Energy+1.0/Distance;
      }
    }
  
   cout << fixed << setprecision(10) << "energy=" << Energy << "\n";  
  
  /* Save Values so far */
  
  for(i=1;i<=n;i++){
    for(j=1;j<=3;j++){
      y[i][j]=x[i][j];
    }
  }
 
  /* Choose each point in turn and make it the north pole note this is what the while loop os for, but have set it to a<2 to just look at first point */
  
  a=1;
  b=0;
  c=0;
  
  while(a<2){
  
  /* Find theta and phi to rotate points by */
  
  cout << fixed << setprecision(5) << "x[" << a << "][1]=" << x[a][1] << 
  " x[" << a << "][2]=" << x[a][2] << " x[" << a << "][3]=" << x[a][3] << "\n";
  
  theta=x[a][2]/x[a][3];
  theta=b*PI+atan(theta);
 
  /* Rotate Points by theta around x axis and then by phi around y axis */
  
  for(i=1;i<=n;i++){
    x[i][1]=x[i][1];
    x[i][2]=x[i][2]*cos(theta)-x[i][3]*sin(theta);
    x[i][3]=x[i][2]*sin(theta)+x[i][3]*cos(theta);
  }
 
  phi=x[a][1]/x[a][3];
  phi=c*PI+atan(phi);
  
  for(i=1;i<=n;i++){
    x[i][1]=x[i][1]*cos(phi)-x[i][3]*sin(phi);
    x[i][2]=x[i][2];
    x[i][3]=x[i][1]*sin(phi)+x[i][3]*cos(phi);
  }
  
  cout << fixed << setprecision(5) << "x[" << a << "][1]=" << x[a][1] << 
  " x[" << a << "][2]=" << x[a][2] << " x[" << a << "][3]=" << x[a][3] << "\n";
    
   if(x[a][3]==1.0 && x[a][2]==x[a][3]==0)
    a=a+1;
  else if(b==0 && c==0)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        c=1;
      }
    }
  else if(b==0 && c==1)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        b=1;
        c=0;
      }
    }
  else if(b==1 && c==0)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        c=1;
      }
    }
  else if(b==1 && c==1)
    break;
 
  }
 
  energy=0.0;
  
  for(i=1;i<=n;i++){
    for(j=i+1;j<=n;j++){
      
      Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                    +pow(x[i][3]-x[j][3],2));
      energy=energy+1.0/Distance;
    }
  }  
    
  cout << fixed << setprecision(10) << "ENERGY=" << energy << "\n";  
  cout << fixed << setprecision(5) << "x[" << a << "][1]=" << x[a][1] << 
  " x[" << a << "][2]=" << x[a][2] << " x[" << a << "][3]=" << x[a][3] << "\n";
  
  /* Output to help with gnuin.txt */
  ofstream File4 ("mypoints");
  for(i=1;i<=n;i++){
    File4 << x[i][1] << " " <<   x[i][2] << " " << x[i][3] << "\n";
  }
  File4.close(); 
  
  return 0;

}

Ok, so there are loads of issues here, like the if statement (line 103) shouldn't really have a condition for equal to a double as it will never work, but I can sort that out later using indirect comparison (some epsilon stuff). My real query is why does the rotation even though it is acting on all the points take the points off the sphere?

(As you can see the values have been normalised to make them all on the unit sphere in line 38).

I am not sure sure if this a maths issue or a c++ issue, but I find that most c++ users understand maths, but not necessarily the other way around which is why I have asked here. (Note: the energy of the system should not change by just rotating points around the sphere as a whole).

Thanks, A.

[adrem7]

So I screwed up on line 85 by using the same variable. Should look like this:

 for(i=1;i<=n;i++){
    new_x[i][1]=x[i][1];
    new_x[i][2]=x[i][2]*cos(theta)-x[i][3]*sin(theta);
    new_x[i][3]=x[i][2]*sin(theta)+x[i][3]*cos(theta);
  }
 
  phi=new_x[a][1]/new_x[a][3];
  phi=c*PI+atan(phi);
  
  for(i=1;i<=n;i++){
    x[i][1]=new_x[i][1]*cos(phi)-new_x[i][3]*sin(phi);
    x[i][2]=new_x[i][2];
    x[i][3]=new_x[i][1]*sin(phi)+new_x[i][3]*cos(phi);
  }
  

My bad.