Skip to main content
GameDev.net gamedev.net
🔒 Locked

Calculating Phi

Started by csuguy Sep 16, 2008 at 11:41 PM 10 replies 2.2k views
Original Post
csuguy
csuguy
Hey all, Ok, so I'm working on 3D rotations. Thanks to various online resources, I was able to understand these basic formulas: x = cos(phi)*sin(theta) y = sin(phi)*sin(theta) z = cos(theta) And then you times them by rho to scale of course. Anyways, I'm trying to take it a little further. I'm trying to do some rotations without phi and theta - and then calculate the phi and theta again. This is for rotating around the Y axis. Basically - I'm stuck when it comes to calculating phi. Using the above equations, I naturally got: phi = acos( x / sin(theta) ) Providing I'm not making some elementary mistake when reformatting the equation (wouldn't be the first time) that SHOULD calculate phi. HOWEVER, when I plug that into my c++ code - it, well... doesn't work. And I plug it in before I times it by Rho - so that's not the problem. Looking at it again, I thought that it obviously wouldn't work as it would only work well when sin(theta) was ~1 or ~-1. If its a small number, then its going to divide x multiple times - and thus the value would be greater than 1 (or less than -1) and thus acos wouldn't work. If sin(theta) returns 0 - I get a divide by zero error... So - how do I calculate this angle? Is there another equation - if so please describe in detail or provide good resource on the subject. Thank you, Ryan
Zipster
Zipster
Just think of the sin(theta) as a factor that makes the circumference circle smaller the closer you get to the poles of the sphere. It doesn't actually affect the shape of the circle, just makes it bigger or smaller. So in essence, if you look at it from the top down you still have a circle shape:

x = cos(phi) * arbitrary_scalar
y = sin(phi) * arbitrary_scalar

We know that sin(phi)/cos(phi) = tan(phi), so it follows that phi = atan( sin(phi) / cos(phi) ), where we make substitutions from the above equations to get phi = atan(y / x). You should use the atan2 library function to handle all the edge cases, particular when one or both of the parameters is 0.

To get theta, take a look at the circle from the side. You essentially have a right triangle, where the opposite side is the radius of the circumference circle at theta. The radius of the circle is just √(x2 + y2) -- also derived from the first two equations. The adjacent side of the right triangle is simply z, since that's the height of the circle from the equator. The properties of right triangles tells us that tan(theta) = √(x2 + y2) / z, or theta = atan( √(x2 + y2) / z). Once again, atan2 will be your friend here and handle all the edge cases!
csuguy
csuguy
Thank you very much! That makes perfect sense! Now to rotate, muahhahahhahahaha :D
Ryan
csuguy
csuguy
Ok - so here's my c++ code for setting phi:

ret.setPhi( atan2f(y,x) );

ret is an instance of my vector class. Anyways - when I use this, it works good for all values where y is positive. Once y hits 0 or negative - it freaks out and moves it to a seemingly random positive spot - where it works fine again. So how do I make this work for the top and bottom of my rotation?

Thanks again,
Ryan

[Edited by - csuguy on September 17, 2008 1:55:14 AM]
Zipster
Zipster
The atan2 function returns a value between -PI/2 and PI/2, and a negative Y value would make Phi negative. Is it possible your code doesn't handle that case?
csuguy
csuguy
Here's the code for the function - I'm doing something wrong :/

Vector Vector::convertToVector(float x,float y,float z)
{
Vector ret;

//Calculate Distance from origin(0,0,0) and then divide x,y, and z by it
float mag = Vector::getDistance(x,y,z,0,0,0);
ret.setMagnitude(mag);

x = x/mag;
y = y/mag;
z = z/mag;

ret.setTheta( toDeg(atanf(getDistance(x,y,0,0,0,0),z)));
ret.setPhi( toDeg(atan2f(y,x)) );

return ret;
}

It is being called after this function:

static void axisRotation(float* d1,float* d2, float angle)
{
float r = sqrt((*d1)*(*d1)+(*d2)*(*d2));

float cur_ang = toDeg(acos(*d1/r));
cur_ang += angle;

*d1 = cos(toRad(cur_ang))*r;
*d2 = sin(toRad(cur_ang))*r;
}

I'm passing in the currentt X and Y floats for d1 and d2 respectively. I'm increasing the angle by +/- 1.0f, depending upon whether the user presses up/down.

My idea to do rotations around an axis is to simply isolate the two other dimensions and do 2D rotations with those 2 dimensions while leaving the third one be. The two dimensions would keep the same scale and thus it would remain a unit sphere. This maybe a bad way of doing it - but it SEEMS like it should work. At the very least I should be able to get the basic XY rotation down using this method. I did something wrong somewhere... If you can find it please tell me and explain.

Thank you!
Ryan
Zipster
Zipster
A few things I noticed (in the comments). I believe your problem is in axisRotation on the acos call:
Vector Vector::convertToVector(float x,float y,float z){/* Is this vector initialized?  If not, then setting the magnitude below doesn't make sense */Vector ret;//Calculate Distance from origin(0,0,0) and then divide x,y, and z by it/* The distance from the origin is just the length of the vector [x,y,z], so that should eliminate a subtraction (hopefully the compiler optimizes it anyway but it would be more clear if you just did it here) */float mag = Vector::getDistance(x,y,z,0,0,0);ret.setMagnitude(mag);x = x/mag;y = y/mag;z = z/mag;/* Again, probably more clear just to do sqrt(x*x + y*y) */ret.setTheta( toDeg(atanf(getDistance(x,y,0,0,0,0),z)));ret.setPhi( toDeg(atan2f(y,x)) );return ret;}static void axisRotation(float* d1,float* d2, float angle){float r = sqrt((*d1)*(*d1)+(*d2)*(*d2));/* Here is likely your problem.  acos is not a good function to use because it only returns a value in the range [0,PI], which is positive Y.  atan2(*d2,*d1) should work a lot better since it disambiguates all quadrants and handles zeros properly */float cur_ang = toDeg(acos(*d1/r));cur_ang += angle;*d1 = cos(toRad(cur_ang))*r;*d2 = sin(toRad(cur_ang))*r;}
Cornstalks
Cornstalks
Your convertToVector function looks wrong. I don't quite get the goal of the axisRotation function so I'm not going to really look at it, plus I trust Zipster [smile]. Using the formula

α = atan2(|u×v|, u⋅v)

where |u×v| is the cross product of u and v and u⋅v is the dot product of u and v, and α is the angle between the vectors u and v.

To find φ, the angle between the z-axis and the vector (x, y, z), we simply do the following:

u = (0, 0, 1)
v = (x, y, z)

φ = atan2(|u×v|, u⋅v)

φ = atan2(√(x2 + y2), z)

Of course, this is exactly what Zipster said in his original post. But we can also use this method to find θ, the angle between the x-axis and the vector (x, y, 0). Now we just have to do

u = (1, 0, 0)
v = (x, y, 0)

θ = atan2(|u×v|, u⋅v)

θ = atan2(y, x)

So finally the function convertToVecotr becomes this:

Vector Vector::convertToVector(float x,float y,float z){    Vector vector;        vector.setMagnitude(std::sqrt(x * x + y * y + z * z));    vector.setPhi(std::atan2(std::sqrt(x * x + y * y), z));    vector.setTheta(std::atan2(y, x));        return vector;}


Hopefully I didn't make any mistakes when writing this post... I'm not going to guarantee it's all perfect.
csuguy
csuguy
Quote:
Original post by Zipster
A few things I noticed (in the comments). I believe your problem is in axisRotation on the acos call:
*** Source Snippet Removed ***


Thank you Zipster - I made some of your recommended changes, particularly with the acos and atan. Can't believe I didn't catch that one myself. Unfortunately that didn't solve the problem :/ - I get some pretty random movements. I'm going to take what you and Mike have said and see if I can't solve it before posting my new code.

Again - thanks for all your help!
Ryan
csuguy
csuguy
Ah! Lots of math lol. Anyways - I tried plugging in your formula (with a few corrections) but it didn't work :/. If I pass it a point like: (50,0,0) it turns it into (0,0,50)?! lol. I noticed that your equation for setPhi and setTheta were reversed so I switched them - but I just got the same results...

In order to work this out I'm going to take what you and Zipster have said and I'm going to make a sample program to test everything out. No custom classes - just the math.

If I can't get that to work then I'll post my sample code to see if guys can't solve it.

Thanks so much for your help!
Ryan

Cornstalks
Cornstalks
Quote:
Original post by csuguy
Ah! Lots of math lol. Anyways - I tried plugging in your formula (with a few corrections) but it didn't work :/. If I pass it a point like: (50,0,0) it turns it into (0,0,50)?! lol. I noticed that your equation for setPhi and setTheta were reversed so I switched them - but I just got the same results...

In order to work this out I'm going to take what you and Zipster have said and I'm going to make a sample program to test everything out. No custom classes - just the math.

If I can't get that to work then I'll post my sample code to see if guys can't solve it.

Thanks so much for your help!
Ryan

Yeah, I switched phi and theta like that because it's a standard math convention for phi to be the angle between the z-axis and the point (plus I hadn't realized you were using the opposite). Clicky. If you switch phi and theta, you get exactly what you had before (except for normalizing the vector, because that step is unnecessary). Anyway, I'll leave it up to you to choose what theta and phi represent. Just realize that if you are using the standard convention, then

x = rho * sin(phi) * cos(theta)
y = rho * sin(phi) * sin(theta)
z = rho * cos(theta)

Using my original code and not changing the code in your original post would cause some weirdness when converting to cartesian coordinates.

The only other thing I can think of is to make sure you are passing in the proper coordinates for the rotation axis. i.e. if you are rotating about the z-axis you do axisRotation(&x, &y, zangle); if you were rotating about the y-axis then you do axisRotation(&z, &x, yangle); if you were rotating about the x-axis you would do axisRotation(&y, &z, zangle).
csuguy
csuguy
SUCCESS!!!! I have successfully written a program (code below) which will rotate around a sphere. It not only uses the arrow keys to change the phi and theta values, but it allows for axis rotation using the '-' and '+/=' keys! You can set the rate at which it rotates using the number keys 1-9! Oh - and you can switch the axis your rotating around by hitting x,y, or z respectively.

There are a few problems I'm not sure how to solve, though they are pretty small errors. 1) Sometimes it will randomly go back to a previous position. 2) Sometimes there is a slight jump when you switch from axis rotation to normal phi theta rotation. There are pretty small errors - but if you have any know-how as to how to fix these do tell.

Well, compile it and tell me what you think:


/* Basic Rotation Program
Written By: Ryan Williams
Copyright 2008

Special Thanks:
Zipster &
MikeTacular
*/

#include <GL/glut.h>
#include <stdlib.h>
#include <math.h>
#include <fstream>

#define PI 3.14159265358979323846

float toDeg(float rad);
float toRad(float deg);

/*Variables*/
float rho = 25.0f,phi = toRad(90.0f), theta = 0.0f;
float x,y,z;
float cx = 0.0f,cy = 0.0f,cz = 100.0f;
char axis = 'z'; //Which axis to rotate around
float deg_chg = 1.0f;

std::ofstream out("diagnostics.txt",std::ios::out);

void axisRotation(float* d1,float* d2, float deg)
{
out << "axisRotation()\n";

float m = sqrt( (*d1)*(*d1) + (*d2)*(*d2) );
out << "Calculated Magnitude: " << m << "\n";

float a = toDeg(atan2f( (*d2)/m,(*d1)/m )) + deg;
out << "Calculated Angle: " << a << "\n";

(*d1) = cos(toRad(a))*m;
(*d2) = sin(toRad(a))*m;

out << "Calculated D1: " << (*d1) << "\n";
out << "Calculated D2: " << (*d2) << "\n";
out << "\n"<< "\n"<< "\n";
}

float toDeg(float rad)
{
return rad*180.0f/PI;
}

float toRad(float deg)
{
return deg*PI/180.0f;
}

void toVector()
{
out << "toVector() \n";
out << "X: " << x << "\n";
out << "Y: " << y << "\n";
out << "Z: " << z << "\n";
out << "OLD RHO: "<<rho<< "\n";
out << "OLD PHI: "<<toDeg(phi)<<"\n";
out << "OLD THETA: "<<toDeg(theta)<< "\n\n";


rho = sqrt( x*x + y*y + z*z );
out << "NEW RHO: " << rho << "\n";
float z2 = (z>0)? 1-(z/rho) : 1+(z/rho);
phi = atan2(z2,z/rho);
out << "NEW PHI: "<< toDeg(phi) << "\n";
theta = atan2(y/rho/sin(phi), x/rho/sin(phi));
out << "NEW THETA: " << toDeg(theta) << "\n";
out << "\n"<< "\n"<< "\n";
}

void toCoords()
{
out << "toCoords() \n";
out << "OLD X: " << x << "\n";
out << "OLD Y: " << y << "\n";
out << "OLD Z: " << z << "\n";
out << "RHO: "<<rho<< "\n";
out << "PHI: "<<toDeg(phi)<<"\n";
out << "THETA: "<<toDeg(theta)<< "\n\n";
x = rho*cos(theta)*sin(phi);
y = rho*sin(theta)*sin(phi);
z = rho*cos(phi);
out << "NEW X: " << x << "\n";
out << "NEW Y: " << y << "\n";
out << "NEW Z: " << z << "\n";
out << "\n"<< "\n"<< "\n";
}


static void resize(int w, int h)
{
if(h == 0)
h = 1;


float ratio = 1.0* w / h;

// Reset the coordinate system before modifying
glMatrixMode(GL_PROJECTION);
glLoadIdentity();

// Set the viewport to be the entire window
glViewport(0, 0, w, h);

// Set the correct perspective.
gluPerspective(45,ratio,1,1000);

glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
//cm.view();
}

static void display(void)
{
glMatrixMode(GL_MODELVIEW);
glClearColor (0.0,0.0,0.0,1.0);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glLoadIdentity();
gluLookAt(cx,cy,cz,0,0,0,0,1,0);

glPushMatrix();
glColor4f(1.0, 1.0, 1.0,.2);
glBegin(GL_LINES);
glVertex3f(-100.0,0.0,0.0);
glVertex3f(100.0,0.0,0.0);
glVertex3f(0.0,-100.0,0.0);
glVertex3f(0.0,100.0,0.0);
glVertex3f(0.0,0.0,-100.0);
glVertex3f(0.0,0.0,0.0);
glVertex3f(0.0,0.0,100.0);
glVertex3f(0.0,0.0,0.0);
glEnd();
glPopMatrix();

glPushMatrix();
glColor3d(1,0,0);
glTranslatef(x,y,z);
glutSolidCube(5);
glPopMatrix();

glutSwapBuffers();
}


static void key(unsigned char key, int kx, int ky)
{
switch (key)
{
case 27 :
case 'q':
out.close();
exit(0);
break;

case '+':
case '=':
if(axis == 'x')
axisRotation(&z,&y,deg_chg);
else if(axis == 'y')
axisRotation(&x,&z,deg_chg);
else
axisRotation(&x,&y,deg_chg);
toVector();
break;

case '-':
if(axis == 'x')
axisRotation(&z,&y,-deg_chg);
else if(axis == 'y')
axisRotation(&x,&z,-deg_chg);
else
axisRotation(&x,&y,-deg_chg);
toVector();
break;

case 'x':
axis = 'x';
break;

case 'y':
axis = 'y';
break;

case 'z':
axis = 'z';
break;

case '1':
deg_chg = 1.0f;
break;

case '2':
deg_chg = 2.0f;
break;

case '3':
deg_chg = 3.0f;
break;

case '4':
deg_chg = 4.0f;
break;

case '5':
deg_chg = 5.0f;
break;

case '6':
deg_chg = 6.0f;
break;

case '7':
deg_chg = 7.0f;
break;

case '8':
deg_chg = 8.0f;
break;

case '9':
deg_chg = 9.0f;
break;
}

display();
}

void key_special(int key, int x, int y)
{
switch(key)
{
case GLUT_KEY_UP:
//Increase Phi
out<< "UP!\n";
phi = toRad(toDeg(phi)+deg_chg);
toCoords();
break;

case GLUT_KEY_DOWN:
//Decrease Phi
out<< "DOWN!\n";
phi = toRad(toDeg(phi)-deg_chg);
toCoords();
break;

case GLUT_KEY_LEFT:
//Decrease Theta
out<< "LEFT!\n";
theta = toRad(toDeg(theta)-deg_chg);
toCoords();
break;

case GLUT_KEY_RIGHT:
//Increase Theta
out<< "RIGHT!\n";
theta = toRad(toDeg(theta)+deg_chg);
toCoords();
break;
}

if(phi < 0) phi += 360;
if(theta > 360) theta += 360;
if(phi > 360) phi -= 360;
if(theta > 360) theta -= 360;

display();
}

void init()
{
out << "INIT()\n";
toCoords();
}

int main(int argc, char *argv[])
{
init();
glutInit(&argc,argv);

glutInitDisplayMode(GLUT_DEPTH|GLUT_DOUBLE|GLUT_RGBA); /*Glut_Double provides a double buffer for smooth animation*/

glutInitWindowPosition(100,100); //Let it decide where to put it
glutInitWindowSize(800,600);
glutCreateWindow("Pocket Billiards");

glutDisplayFunc(display);
glutIdleFunc(display);
glutReshapeFunc(resize);

glutKeyboardFunc(key);
glutSpecialFunc(key_special);

glEnable(GL_DEPTH_TEST);
glutMainLoop();

return EXIT_SUCCESS;
}

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.

Sign in to reply to this topic.