PVector rotatePointAboutAxis(final PVector p,final float theta,final PVector p1,final PVector p2){

 
 float d;
 
 PVector q1=PVector.sub(p,p1);
 PVector u=PVector.sub(p2,p1);
 u.normalize();
 d=sqrt(u.y*u.y+u.z*u.z);
 
 PVector q2=new PVector();
 
 if(d!=0){
   q2.x=q1.x;
   q2.y=q1.y*u.z/d-q1.z*u.y/d;
   q2.z=q1.y*u.y/d+q1.z*u.z/d;
 }
 else{
  q2.set(q1); 
 }
  
 q1.x=q2.x*d-q2.z*u.x;
 q1.y=q2.y;
 q1.z=q2.x*u.x+q2.z*d;
 
 float ct=cos(theta);
 float st=sin(theta);
 
 q2.x=q1.x*ct-q1.y*st;
 q2.y=q1.x*st+q1.y*ct;
 q2.z=q1.z;
 
 q1.x=q2.x*d+q2.z*u.x;
 q1.y=q2.y;
q1.z=-q2.x*u.x+q2.z*d;

if(d!=0){
   q2.x=q1.x;
   q2.y=q1.y*u.z/d+q1.z*u.y/d;
   q2.z=-q1.y*u.y/d+q1.z*u.z/d;
 }
 else{
  q2.set(q1); 
 }

q1.x=q2.x+p1.x;
q1.y=q2.y+p1.y;
q1.z=q2.z+p1.z;
  
 return q1;
}