What is the formula behind atan2?

[ETAkoG1T]

I stumbled upon this function but couldn't find some clear explanation of the formula used. I basically want to know how to make my own easy copy of the function. (I like doing things manually just for practice and to understand it better)

This is the source code for a working program were I used the function:

#include <iostream>
#include "stdafx.h"

struct polar
{
	double distance;
	double angle;
};
struct rect
{
	double x;
	double y;
};
// Function prototypes
polar rect_to_polar(rect xypos);
void show_polar(polar dapos);

int main()
{
	using namespace std;
	rect rplace;
	polar pplace;

	cout << "Enter the x and y values: ";
	while (cin >> rplace.x >> rplace.y)
	{
		pplace = rect_to_polar(rplace);
		show_polar(pplace);
		cout << "Next two numbers (q to quit) : ";
	}
	cout << "Done. \n";
	return 0;
}

//Rectangular to polar coordinates
polar rect_to_polar(rect xypos)
{
	using namespace std;
	polar answer;

	answer.distance = sqrt (xypos.x * xypos.x + xypos.y * xypos.y);
	answer.angle = atan2 (xypos.y, xypos.x);
	return answer;
}


//Show polar coordinates, converting angle to degrees
void show_polar (polar dapos)
{
	using namespace std;
	const double Rad_to_deg = 57.29577951;

	cout << "distance = " << dapos.distance;
	cout << ", angle = " << dapos.angle * Rad_to_deg;
	cout << " degrees\n";
}

It successfully compiled in visual studio 2012, but I would like to know exactly what the function does :)
Thanks in advance!

[Thumper]

atan stands for arctangent, or the inverse tangent. The function atan2() takes the ratio x/y and returns the value of the angle in radians.
http://www.cplusplus.com/reference/clibrary/cmath/atan2/

[Zhuge]

You could write your own function that approximates arctan() using the Taylor polynomial approximation, which IIRC is atan2(y,x) = atan(y/x) = (y/x) - (1/3)(y/x)^3 + (1/5)(y/x)^5 ...

It continues with odd powers of (y/x), odd denominators to the coefficent, and alternating + / - signs. Take the polynomial out until you have the precision you want.

There may also be a table of lookups to help speed the calculation for values where the series doesn't converge very quickly.

[ETAkoG1T]

thanks :)
But the thing im really interested in knowing is if I were to write a formula or just calculate it without using a function how would you convert point x and y to radian? for instance,
x = 40 y = 50
radians = 0.89055
how do you come to 0.89055 without using atan? I searched around but couldn't find anywere talking about the actual formula only an explanation of what it does. I don't want to dismiss it as simply being magic!

I am only 14 years old so keep that in mind, it might be obvious to every one except me :P

EDIT:
Thanks Zhuge im getting closer! I couldn't find that formula anywere!

[ETAkoG1T]

Ok, this is getting too anoying! Can you just show how you would do it out with my code(change atan2 and that stuff with the formula).
I best I could do was something like this:

answer.angle = (xypos.y/xypos.x) - ((1/3)*(xypos.y/xypos.x))*((1/3)*(xypos.y/xypos.x))*((1/3)*(xypos.y/xypos.x)) + ((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x));

Im not sure if you can use ^ in c++. The answer I got for x = 40 and y= 50 was 1.25 instead of 0.89055, and it looks like it has nothing to do with accuracy!

[D80DSL3A]

Try this:

#include <iostream>
#include<iomanip>// for setprecision, setw
#include<cmath>

using namespace std;

double PI = 3.141592654;

// unfortunately, 3 formulas apply to different ranges of x
// atan(x) = x - x^3/3 + x^5/5 - ...                                     for -1 < x < 1
// atan(x) =  PI/2 - x^-1 + (3*x^3)^-1 - (5*x^5)^-1 + ... for x > 1
// atan(x) = -PI/2 - x^-1 + ... as above                              for x < -1
double myAtan(double x, int n, double& r_err )
{
	double a = 0.0;// 1st term
	double sum = 0.0;

	// special cases
	if( x == 1.0 ) return PI/4.0;
	if( x == -1.0 ) return -PI/4.0;

	if(n > 0)
	{
	    if( (x < -1.0) || (x > 1.0) )
	    {
	        // constant term
	        if( x > 1.0 )
                sum = PI/2.0;
            else
                sum = -PI/2.0;
            // initial value of a
            a = -1.0/x;
            for(int j=1; j<=n; j++)
            {
                sum += a;
                a *= -1.0*(2.0*j-1)/((2.0*j+1)*x*x);// next term from last
            }
	    }
	    else// -1 < x < 1
	    {
	        // constant term
	        sum = 0.0;
	        // initial value of a
            a = x;
            for(int j=1; j<=n; j++)
            {
                sum += a;
                a *= -1.0*(2.0*j-1)*x*x/(2.0*j+1);// next term from last
            }
	    }
		r_err = a;// max. error = 1st term not taken for alternating series
	}

	return sum;
}// end of mySine()

// lets compare some values with the "real" atan() function
int main()
{
    setprecision(8);

    int Nterms = 0;
    double errorAmt = 3.0;
    double x[] = { 0.1, 0.3, 0.5, 0.9, 1.1, 2.5, 5.0, 10.0 };// will use these values
    int Ntrials = sizeof(x)/sizeof(double);

    cout << "Nterms = "; cin >> Nterms;

    for(int i=0; i<Ntrials; ++i  )
    {
        cout << "atan(" << setprecision(3) << x[i] << ") = " << setprecision(8) << atan(x[i]);
        cout << "  myAtan(" << x[i] << ") = " << setprecision(8) << myAtan(x[i], Nterms, errorAmt );
        cout << "  err = " << errorAmt << endl;
    }

    cout << endl;
    return 0;
}

I get the following output for N = 50 terms in all cases:

Nterms = 50
atan(0.1) = 0.099668652  myAtan(0.1) = 0.099668652  err = 9.9009901e-104
atan(0.3) = 0.29145679  myAtan(0.3) = 0.29145679  err = 1.5308243e-055
atan(0.5) = 0.46364761  myAtan(0.5) = 0.46364761  err = 3.905252e-033
atan(0.9) = 0.7328151  myAtan(0.9) = 0.73281497  err = 2.3668573e-007
atan(1.1) = 0.83298127  myAtan(1.1) = 0.83298163  err = -6.5315676e-007
atan(2.5) = 1.1902899  myAtan(2.5) = 1.1902899  err = -6.3641111e-043
atan(5) = 1.3734008  myAtan(5) = 1.3734008  err = -2.5101992e-073
atan(10) = 1.4711277  myAtan(10) = 1.4711277  err = -9.9009901e-104


Process returned 0 (0x0)   execution time : 3.125 s
Press any key to continue.

I see the accuracy drops to 6 significant figures around x = 1.0, but it's pretty low elsewhere. Try varying the number of terms (Ntrems) to see how that affect accuracy.

EDIT: Thought I'd better check at negative values of x:

Nterms = 50
atan(-0.1) = -0.099668652  myAtan(-0.1) = -0.099668652  err = -9.9009901e-104
atan(-0.3) = -0.29145679  myAtan(-0.3) = -0.29145679  err = -1.5308243e-055
atan(-0.5) = -0.46364761  myAtan(-0.5) = -0.46364761  err = -3.905252e-033
atan(-0.9) = -0.7328151  myAtan(-0.9) = -0.73281497  err = -2.3668573e-007
atan(-1.1) = -0.83298127  myAtan(-1.1) = -0.83298163  err = 6.5315676e-007
atan(-2.5) = -1.1902899  myAtan(-2.5) = -1.1902899  err = 6.3641111e-043
atan(-5) = -1.3734008  myAtan(-5) = -1.3734008  err = 2.5101992e-073
atan(-10) = -1.4711277  myAtan(-10) = -1.4711277  err = 9.9009901e-104

EDIT2: That's right, you wanted atan2():

// finds atan(y/x) in the range -PI to +PI
double myAtan2(double y, double x, int n, double& r_err )
{
    double u = myAtan( y/x, n, r_err );
    if( x < 0.0 )// 2nd, 3rd quadrant
    {
        if( u > 0.0 )// will go to 3rd quadrant
            u -= PI;
        else
            u += PI;
    }
    return u;
}

[ETAkoG1T]

Thanks!

But I don't really understand how to find out
y - axis
5
4
3 x
2
1 2 3 4 5 - x axis

How do I find out the angle of x?
I tried this on my program and got 30.96 but this is with atan2 and I want to find out how to do it manually, what you have found looks like it is from 1 value and not 2. Would you care to explain how I can use the information I get from your program to find out e.g. (5, 3)? Thanks anyways looks like you put alot of effort into it. It would also be nice if you explained more about exactly what your program did! :)

[D80DSL3A]

Yes, I gave a function for atan() first.

I gave a function for atan2() at the end of my post - myAtan2. This function takes 2 arguments (y,x) same as the "real" atan2().

To use myAtan2() you give 4 numbers:
double myAtan2(double y, double x, int n, double& r_err )
y and x are the values used in the regular atan2().
n = # of terms to sum the power series. I used n = 50 in my examples.
Finally, r_err is there to get the value of the error in stopping the sum at n terms. You must pass a variable for this (not a value). See my program above.
Use:

double errorAmt = 0.0;
double y = 3.0, x = 2.0;

cout << "atan2(" << y <<  ", "  <<  x << ") = " << setprecision(8) << atan2(y,x);
cout << "  myAtan2(" << y <<  ", "  <<  x << ") = " << setprecision(8) << myAtan2(y, x, 50, errorAmt );
cout << "  err = " << errorAmt << endl

The function myAtan() simply finds each term in the power series for arctan(x) and adds them up.

[Imadatobanisa]

@fun2code
I don't know what the function 'atan' and 'atan2' are used for? Can you explain them and give me some examples?
Thanks you.

[D80DSL3A]

atan and atan2 are the function names for the inverse trigonometric function arctangent, or tan^-1. See:
http://en.wikipedia.org/wiki/Inverse_trigonometric_functions

See our reference section for info on the cmath header functions.
http://www.cplusplus.com/reference/clibrary/cmath/

[ETAkoG1T]

Yay, Thank you it actually worked! I don't really understand what n our & r_err does, but it worked now :) Ty, the case is solved, but I still have no idea how you did it! :)

[D80DSL3A]

You're welcome.
I was looking over one of your past posts and I saw that you nearly got it right yourself.
Where you tried:

answer.angle = (xypos.y/xypos.x) - ((1/3)*(xypos.y/xypos.x))*((1/3)*(xypos.y/xypos.x))*((1/3)*(xypos.y/xypos.x)) + ((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x))*((1/5)*(xypos.y/xypos.x));

You divided by 3 too much (and 5 in the 2nd term).
More simply:

double u = y/x;// for convenient substitution
double atanX = u - u*u*u/3.0 + u*u*u*u*u/5.0;// good for u*u < 1 

That gets it to 3 terms.
When we pass a value of 50 for n when calling myAtan2() we are getting our value to 50 terms!

[ETAkoG1T]

That looks like it solved some of the problem but the thing is that it is only accurate if the x axis number is higher than the y axis number! Or else it doesn't work properly.

This is all the code I needed to make it pretty accurate, only problem is that it only works if x > y, I got no idea why!

double x = (pxy->y) / (pxy->x);
	pda->angle = x - x*x*x/3.0 + x*x*x*x*x/5.0 - x*x*x*x*x*x*x/7.0 + x*x*x*x*x*x*x*x*x/9.0;   

I don't need to handle negative numbers. Do you know why and perhaps a way to fix this? ty!

[rcor]

Another option is employing a rational approximation as described in

http://ieeexplore.ieee.org/xpl/articleDetails.jsp?tp=&arnumber=6375931&contentType=Journals+%26+Magazines&queryText%3Dfull+quadrant+approximation


The following code uses the rational approximation to get the arctangent normalized to the [0 1) interval (you can multiply it by Pi/2 to get the real arctangent)

normalized_atan(x) ~ (B x + x^2) / (1 + 2 B x + x^2)

where B = 0.596227

The maximum error is 0.1620 º

#include <stdint.h>
#include <math.h>

// Approximates atan(x) normalized to the [-1,1] range
// with a maximum error of 0.1620 degrees.

float norm_atan( float x )
{
	static const uint32_t sign_mask = 0x80000000;
	static const float b = 0.596227f;
	
	// Extract the sign bit
	uint32_t ux_s  = sign_mask & (uint32_t &)x;

	// Calculate the arctangent in the first quadrant
	float bx_a = ::fabs( b * x );
	float num = bx_a + x * x;
	float atan_1q = num / ( 1.f + bx_a + num );
	
	// Restore the sign bit
	uint32_t atan_2q = ux_s | (uint32_t &)atan_1q;
	return (float &)atan_2q;
}

// Approximates atan2(y, x) normalized to the [0,4) range
// with a maximum error of 0.1620 degrees

float norm_atan2( float y, float x )
{
	static const uint32_t sign_mask = 0x80000000;
	static const float b = 0.596227f;

	// Extract the sign bits
	uint32_t ux_s  = sign_mask & (uint32_t &)x;
	uint32_t uy_s  = sign_mask & (uint32_t &)y;
	
	// Determine the quadrant offset
	float q = (float)( ( ~ux_s & uy_s ) >> 29 | ux_s >> 30 ); 

	// Calculate the arctangent in the first quadrant
	float bxy_a = ::fabs( b * x * y );
	float num = bxy_a + y * y;
	float atan_1q =  num / ( x * x + bxy_a + num );
	
	// Translate it to the proper quadrant
	uint32_t uatan_2q = (ux_s ^ uy_s) | (uint32_t &)atan_1q;
	return q + (float &)uatan_2q;
}

In case you need more precision, there is a 3rd order rational function:

normalized_atan(x) ~ ( C x + x^2 + x^3) / ( 1 + (C + 1) x + (C + 1) x^2 + x^3)

where C = (1 + sqrt(17)) / 8

which has a maximum approximation error of 0.00811º