arctan2 algorithm?

[Simon Brooke]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly low
level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

 

[Patricia Shanahan]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Are you sure you need atan2? Distances between arbitrary points can be
calculated, given x,y coordinates, using only add, subtract, multiply,
and square root.

Patricia

 

[Thomas Fritsch]

Are you sure you need atan2? Distances between arbitrary points can be
calculated, given x,y coordinates, using only add, subtract, multiply,
and square root.

(+ - * sqrt) are sufficient only for calculating distances on a /flat/
2-dimensional space.
But I think the OP talks about calculating large distances on the earth's
sphaerical surface from given longitudes/latitudes. For such calculations
one also needs: sin cos arccos. Because Java's Math has no arccos, one has
to calculate arccos somehow, probably by arctan2.

 

[Simon Brooke]

Are you sure you need atan2? Distances between arbitrary points can be
calculated, given x,y coordinates, using only add, subtract, multiply,
and square root.

Not on the surface of the Earth, they can't. The Earth is (roughly)
spherical.

I agree that a flat Earth would make life simpler.

 

[Simon Brooke]

(+ - * sqrt) are sufficient only for calculating distances on a /flat/
2-dimensional space.
But I think the OP talks about calculating large distances on the earth's
sphaerical surface from given longitudes/latitudes. For such calculations
one also needs: sin cos arccos. Because Java's Math has no arccos, one
has to calculate arccos somehow, probably by arctan2.

The algorithm I'm (haversine formula) using uses arctan2 rather than
arccos:

/**
* compute the distance in metres between this location and x.
*
* @param x the remote location
*
* @return the distance in metres
*/
public int distance( NMEALocation x )
{
int metres = 0;
double deltalat = this.latitude - x.latitude;
double deltalong = this.longitude - x.longitude;

double a =
( Math.sin( deltalat / 2 ) * Math.sin( deltalat / 2 ) ) +
( Math.cos( this.latitude ) * Math.cos( x.latitude ) *
Math.sin( deltalong / 2 ) * Math.sin( deltalong / 2 ) );
double c = 2 * Math.atan2( Math.sqrt( a ), Math.sqrt( 1 - a ) );
double kilometres = EARTHRADIUS * c;

return (int)(kilometres * 1000.0);
}

so the only thing I'm missing is arctan or arctan2.


--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/

;; MS Windows: A thirty-two bit extension ... to a sixteen bit
;; patch to an eight bit operating system originally coded for a
;; four bit microprocessor and sold by a two-bit company that
;; can't stand one bit of competition -- anonymous

 

[Fred Kleinschmidt]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly low
level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/
Iraq war: it's time for regime change...
... go now, Tony, while you can still go with dignity.
[update three years after this .sig was written: it's still
relevant]

Here's a quick approximation, accurate to about .05 radians (in C):
double arctan2( double y, double x ) {
const double ONEQTR_PI = M_PI / 4.0;
const double THRQTR_PI = 3.0 * M_PI / 4.0;
double r, angle;
double abs_y = fabs(y) + 1e-10f;
if ( x < 0.0f ) {
r = (x + abs_y) / (abs_y - x);
angle = THRQTR_PI;
}
else {
r = (x - abs_y) / (x + abs_y);
angle = ONEQTR_PI;
}
angle += (0.1963f * r * r - 0.9817f) * r;
if ( y < 0.0f )
return( -angle );
else
return( angle );
}

 

[tam]

Simon Brooke wrote:
....

The algorithm I'm (haversine formula) using uses arctan2 rather than
arccos:

/**
* compute the distance in metres between this location and x.
*
* @param x the remote location
*
* @return the distance in metres
*/
public int distance( NMEALocation x )
{
int metres = 0;
double deltalat = this.latitude - x.latitude;
double deltalong = this.longitude - x.longitude;

double a =
( Math.sin( deltalat / 2 ) * Math.sin( deltalat / 2 ) ) +
( Math.cos( this.latitude ) * Math.cos( x.latitude ) *
Math.sin( deltalong / 2 ) * Math.sin( deltalong / 2 ) );
double c = 2 * Math.atan2( Math.sqrt( a ), Math.sqrt( 1 - a ) );
double kilometres = EARTHRADIUS * c;

return (int)(kilometres * 1000.0);
}

You can use a formulation that needs only the arc sine.
E.g.,

c = 2*asin(sqrt(sin(deltalat/2)**2 +
cos(lat1)*cos(lat2)*sin(deltalon/2)**2))

where lat1 = this.latitude and lat2=x.latitude
and ** is used for exponentiation..

Regards,
Tom McGlynn

 

[Thomas Fritsch]

The algorithm I'm (haversine formula) using uses arctan2 rather than
arccos:

/**
* compute the distance in metres between this location and x.
*
* @param x the remote location
*
* @return the distance in metres
*/
public int distance( NMEALocation x )
{
int metres = 0;
double deltalat = this.latitude - x.latitude;
double deltalong = this.longitude - x.longitude;

double a =
( Math.sin( deltalat / 2 ) * Math.sin( deltalat / 2 ) ) +
( Math.cos( this.latitude ) * Math.cos( x.latitude ) *
Math.sin( deltalong / 2 ) * Math.sin( deltalong / 2 ) );
double c = 2 * Math.atan2( Math.sqrt( a ), Math.sqrt( 1 - a ) );
double kilometres = EARTHRADIUS * c;

return (int)(kilometres * 1000.0);
}

so the only thing I'm missing is arctan or arctan2.

arctan2(y,x) = + arctan(y/x)
or - arctan(y/x) (depending on the signs of x and y)

I did a quick google with: arctan +taylor

The Taylor series of arctan(x) (see
<http://en.wikipedia.org/wiki/Taylor_series> ) seems quite easy to
implement. But unfortunately it converges well only for small x
(i.e. for distances small compared to the earth radius).

For bigger x this funny page might help:
http://www.opensky.ca/~jdhildeb/arctan/arctan_taylor.html

 

[Patricia Shanahan]

Not on the surface of the Earth, they can't. The Earth is (roughly)
spherical.

I agree that a flat Earth would make life simpler.

But I can SEE the Earth is flat.

Well, if you insist on treating it as an approximate sphere, I agree you
do need at least one of the inverse trig functions. How accurate do you
need it to be?

Patricia

 

[tam]

Simon Brooke wrote:
...


You can use a formulation that needs only the arc sine.
E.g.,

c = 2*asin(sqrt(sin(deltalat/2)**2 +
cos(lat1)*cos(lat2)*sin(deltalon/2)**2))

where lat1 = this.latitude and lat2=x.latitude
and ** is used for exponentiation..

Regards,
Tom McGlynn

Hmmm... I had thought that I had seen that you had access to asin
but not atan, but I see that that's wrong. Still you may find that the
asin expression is easier to handle. Since asin has only a limited
domain, a series expansion might converge more rapidly than for
atan where you have to worry about the entire number line.
Regards,
Tom McGlynn

 

[Simon Brooke]

I'm doing stuff on JavaME for which needs to do distance between
arbitrary geographic locations, so I need Math.atan2( double, double),
which of course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly
low level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/
Iraq war: it's time for regime change...
... go now, Tony, while you can still go with dignity.
[update three years after this .sig was written: it's still
relevant]

Here's a quick approximation, accurate to about .05 radians (in C):
double arctan2( double y, double x ) {
const double ONEQTR_PI = M_PI / 4.0;
const double THRQTR_PI = 3.0 * M_PI / 4.0;
double r, angle;
double abs_y = fabs(y) + 1e-10f;
if ( x < 0.0f ) {
r = (x + abs_y) / (abs_y - x);
angle = THRQTR_PI;
}
else {
r = (x - abs_y) / (x + abs_y);
angle = ONEQTR_PI;
}
angle += (0.1963f * r * r - 0.9817f) * r;
if ( y < 0.0f )
return( -angle );
else
return( angle );
}

Thanks, that looks promising.

 

[William R. Frensley]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly low
level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

Look at M. Abramowitz and I. Segun, "Handbook of Mathematical Functions"
(Dover, 1965, or US Govt Printing Office back in the 1950s), page 81.
They give several series and rational approximations for arctan(z), for
-1<z<1, with different levels of complexity and accuracy.
To do atan2, let z = y/x or x/y, whichever gives a value in the
right range, and then use either the resulting angle, or its complement.

- Bill Frensley

 

[ducnbyu]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly low
level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/
Iraq war: it's time for regime change...
... go now, Tony, while you can still go with dignity.
[update three years after this .sig was written: it's still relevant]

This shows how it's done in hardware using formulas easy to reproduce
in java.

http://www.electronics.dit.ie/staff/aschwarzbacher/research/ecs99.pdf

The key is to store precomputed values of atan(1/(2^n)) (n starts at 0)
and binary search for the desired accuracy. I tested this in Excel and
it converges to 8 decimal places in about 30 iterations. There's some
low hanging fruit for optimization. I didn't find the explanation for
choosing "d" completely clear, but it's just saying if y > 0 then d is
1 if y < 0 then d is -1 and if y = 0 then d is zero because you hit on
an exact answer (based on the accuracy of your precomputed table) and
you don't want further computations to change the result. In other
words you'd just want to stop iterating. If y never gets to zero then
you stop iterating when you run out of precomputed values.

 

[ducnbyu]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly low
level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/
Iraq war: it's time for regime change...
... go now, Tony, while you can still go with dignity.
[update three years after this .sig was written: it's still relevant]

This shows how it's done in hardware using formulas easy to reproduce
in java.

http://www.electronics.dit.ie/staff/aschwarzbacher/research/ecs99.pdf

The key is to store precomputed values of atan(1/(2^n)) (n starts at 0)
and binary search for the desired accuracy. I tested this in Excel and
it converges to 8 decimal places in about 30 iterations. There's some
low hanging fruit for optimization. I didn't find the explanation for
choosing "d" completely clear, but it's just saying if y > 0 then d is
1 if y < 0 then d is -1 and if y = 0 then d is zero because you hit on
an exact answer (based on the accuracy of your precomputed table) and
you don't want further computations to change the result. In other
words you'd just want to stop iterating. If y never gets to zero then
you stop iterating when you run out of precomputed values.

BTW the above is for atan. You will need to query your values of x and
y to produce the correct result in -180 to 180 and -90 to 90 degree
spaces for geo coordinates for a true atan2 result.

 

[Jeff]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly low
level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/
Iraq war: it's time for regime change...
... go now, Tony, while you can still go with dignity.
[update three years after this .sig was written: it's still
relevant]

Here's a quick approximation, accurate to about .05 radians (in C):
double arctan2( double y, double x ) {
const double ONEQTR_PI = M_PI / 4.0;
const double THRQTR_PI = 3.0 * M_PI / 4.0;
double r, angle;
double abs_y = fabs(y) + 1e-10f;
if ( x < 0.0f ) {
r = (x + abs_y) / (abs_y - x);
angle = THRQTR_PI;
}
else {
r = (x - abs_y) / (x + abs_y);
angle = ONEQTR_PI;
}
angle += (0.1963f * r * r - 0.9817f) * r;
if ( y < 0.0f )
return( -angle );
else
return( angle );
}

Similar, written in Casl (a bit closer to basic, perhaps easier to
translate to Java). 3 related routines, arcsin, arcsinr (arcsin in
radians), and arctan, all you need then is a square root routine (I can
supply that too if needed):
function cslASINr(numeric val) as numeric;
# Function to calcuate arcsine of val, in radians
cslASINr =
1.570796327-cslSQRT(1-val)*(1.5707288-0.2121144*val+0.074261*val*val-0.0187293*val*val*val);
end;

function cslASIN(numeric val) as numeric;
# Function to calculate arcsine of val, in degrees
cslASIN = cslASINr(val)/0.01745329252;
end;

function cslATAN(numeric val) as numeric;
# Function to calculate arctangent of val, in degrees
cslATAN = cslASIN(val/cslSQRT(val*val+1));
end;

function cslATANr(numeric val) as numeric;
# Function to calculate arctangent of val, in radians
cslATANr = cslASINr(val/cslSQRT(val*val+1));
end;

 

[Simon Brooke]

Here's a quick approximation, accurate to about .05 radians (in C):

OK, I've gone with an adaptation of this. I have to admit I'm currently
getting wildly inaccurate results, but that's a matter of debugging. The
code I'm using, in case anyone wants to tell me how I've got it wrong, is:

/** One quarter of PI */
protected static final double ONEQTR_PI = Math.PI / 4.0;

/** Three quarters of PI */
protected static final double THRQTR_PI = ( 3.0 * Math.PI ) / 4.0;


/**
* Thanks to Fred Kleinschmidt (fred.l.kleinmschmidt[at]boeing.com>)
* in news post &lt;J9Ly5D.J5q[at]news.boeing.com&gt;
*
* @param y an angle, expressed as a double number of degrees
* @param x a different angle, expressed as a double number of degrees
*
* @return the arctangent between these two angles, as a double number
* of degrees
*/
private double arctan2( double y, double x )
{
y = Math.toRadians( y );
x = Math.toRadians( x );

double result = 0.0;
double r;
double theta;
double abs_y = Math.abs( y ) + 1e-10f;

if ( x < 0.0f )
{
r = ( x + abs_y ) / ( abs_y - x );
theta = THRQTR_PI;
}
else
{
r = ( x - abs_y ) / ( x + abs_y );
theta = ONEQTR_PI;
}

theta += ( ( ( 0.1963f * r * r ) - 0.9817f ) * r );

if ( y < 0.0f )
{
result = Math.toDegrees( 0 - theta);
}
else
{
result = Math.toDegrees( theta);
}

return result;
}


--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/
; gif ye hes forget our auld plane Scottis quhilk your mother lerit you,
; in tymes cuming I sall wryte to you my mind in Latin, for I am nocht
; acquyntit with your Southeron
;; Letter frae Ninian Winyet tae John Knox datit 27t October 1563

 

[Jeff]

I'm doing stuff on JavaME for which needs to do distance between arbitrary
geographic locations, so I need Math.atan2( double, double), which of
course isn't in JavaME.

Bother.

I've looked at Math.atan2 in src.zip, but it's doing a lot of horribly low
level stuff that I don't really understand. Can anyone point me to a
reasonably high-level algorithm for arctangent, in some reasonably high
level language or pseudo-language (doesn't have to be Java)?

--
(e-mail address removed) (Simon Brooke) http://www.jasmine.org.uk/~simon/
Iraq war: it's time for regime change...
... go now, Tony, while you can still go with dignity.
[update three years after this .sig was written: it's still relevant]

This shows how it's done in hardware using formulas easy to reproduce
in java.

http://www.electronics.dit.ie/staff/aschwarzbacher/research/ecs99.pdf

The key is to store precomputed values of atan(1/(2^n)) (n starts at 0)
and binary search for the desired accuracy. I tested this in Excel and
it converges to 8 decimal places in about 30 iterations. There's some
low hanging fruit for optimization. I didn't find the explanation for
choosing "d" completely clear, but it's just saying if y > 0 then d is
1 if y < 0 then d is -1 and if y = 0 then d is zero because you hit on
an exact answer (based on the accuracy of your precomputed table) and
you don't want further computations to change the result. In other
words you'd just want to stop iterating. If y never gets to zero then
you stop iterating when you run out of precomputed values.

BTW the above is for atan. You will need to query your values of x and
y to produce the correct result in -180 to 180 and -90 to 90 degree
spaces for geo coordinates for a true atan2 result.

I apologize for not knowing how to post code properly (sscce???), but I
put this together from a project from another language for PDAs. It is
pretty accurate, does not require anything more than +-*/. I tested
it and in a reasonable range (10-80 degrees tested) it was accurate to
3-4 places precision. It is quick-and-dirty, can be adjusted (iteration
count in the square root function, for example, can be adjusted up or
down).

/**
*
* @author Jeffrey Summers
*/
public class MiniMath {

/** Creates a new instance of MiniMath */
public MiniMath() {
}

public double mm_sqrt(double src) {
//method to calc sqrt using successive approximations
double hi, lo, midl;
midl = 0.0;
if (src > 1) {
lo = 1;
hi = src;
} else {
lo = src;
hi = 1;
}
for (int i = 1; i < 100; i++) {
midl = (hi + lo)/2.0;
if (midl*midl > src) {
hi = midl;
} else {
lo = midl;
}

}
return midl;
}
public double mm_asn(double src) {
// Calculate arcsin from degrees
return
(1.570796327-mm_sqrt(1-src)*(1.5707288-0.2121144*src+0.074261*src*src-0.0187293*src*src*src))/0.01745329252;
}
public double mm_atn(double src) {
// Calculate arctan from degrees
return (mm_asn(src/mm_sqrt(src*src+1)));
}
}

 

[Simon Brooke]

But I can SEE the Earth is flat.

Well, if you insist on treating it as an approximate sphere, I agree you
do need at least one of the inverse trig functions. How accurate do you
need it to be?

Ideally, to within +- 10 metres.

Currently I'm getting results which are NOTHING LIKE that good. Earlier
this morning it had me walking at 1207 Km/h, which I suspect was a trifle
inaccurate...

 

[Patricia Shanahan]

Ideally, to within +- 10 metres.

Currently I'm getting results which are NOTHING LIKE that good. Earlier
this morning it had me walking at 1207 Km/h, which I suspect was a trifle
inaccurate...

10 meters, for a circle the radius of the earth, is about 1.5e-6
radians. I think one of the suggestions was claimed to be accurate to
about 8 decimal places, which should be good enough.

I would test the atan2 implementation in a full Java environment,
comparing it to Math.atan2 which is far more accurate than you need.
Accept any answer that differs by no more than a couple of millionths of
a radian from the Math.atan2 answer.

Patricia

 

[Patricia Shanahan]

Ideally, to within +- 10 metres.

Currently I'm getting results which are NOTHING LIKE that good. Earlier
this morning it had me walking at 1207 Km/h, which I suspect was a trifle
inaccurate...

If your current problem is mainly local walking, you can do a temporary
flat-Earth approximation to use until you get your atan2 working.

Do spherical geometry, in a full Java environment, to get the surface
distance corresponding to a degree of latitude and longitude in your
area. If you approximate to a sphere, the latitude distance is a
constant, and only longitude depends on location.

Pick an arbitrary local origin, and calculate flat Earth cartesian
coordinates relative to it by multiplying the differences longitude and
latitude by the conversion factors.

Patricia