Need help writing a formula in c++

[muhmuuh]

Hi!

First I know this is autoit forum but I've seen many users that are very skilled c++ programmers and also this is the best forum community I've ever seen so I think this the best place to find some help.

To the problem. I'm wirting an application in c++ that calculates distance between cities. The coordinates are given in degrees with latiudes and longitudes. I'm using this formula to calculate the distance.

d=2 R sin^-1 (sin² (a1-a2)/2 + cos a1 cos a2 sin² (b1-b2)/2)^1/2

Here, d is the distance between the towns, a1 and a2 are their latitudes,

b1 and b2 are their longitudes, and R is the radius of the earth.

The latitudes and longitudes are given in degrees but the functions sin and cos work with radians.

I've been working on this for two days and I'm stuck. I can't write this formula correctly in c++.

Here is what I've done

double d;
d=sin((Pi*((a1-a2)/2.0))/180.0);
d*=d;
double dd;
dd=sin((Pi*((b1-b2)/2.0))/180.0);
dd*=dd;
dd*=cos((Pi*a1)/180.0);
dd*=cos((Pi*a2)/180.0);
d+=dd;
d=pow(d, 0.5);
d=2*R/d;
cout << d << endl;

and also

d=2*R/sin(Pi*pow( sin((Pi*(a1-a2)/2.0)/180.0)*sin((Pi*(a1-a2)/2.0)/180.0)+cos(Pi*a1/180.0)*cos(Pi*a2/180.0)*sin((Pi*(b1-b2)/2.0)/180.0)*sin((Pi*(b1-b2)/2.0)/180.0), 0.5)/180.0);

but those lines don't give accurate results.

Can you tell me what I'm doing wrong?

Thank you for your time

p.s. Sorry for my English. I'm still learning it

Edited August 23, 2008 by muhmuuh

[monoceres]

Hi!

You're converting to radians wrong.

Here's how you should do it:

#define PI 3.14159

double DegToRad(double value)
{
       return 180.0/(PI*value);
}

And remember, x^(1/2)=Sqrt(x), I think it is more readable that way

[muhmuuh]

Hi!

You're converting to radians wrong.

Here's how you should do it:

#define PI 3.14159

double DegToRad(double value)
{
       return 180.0/(PI*value);
}

And remember, x^(1/2)=Sqrt(x), I think it is more readable that way

one degree is equivalent to π/180 radians

so x degrees should be equal to x*Pi/180 radians

right?

[monoceres]

so x degrees should be equal to x*Pi/180 radians

right?

Nope.

From the wikipedia page about Radians:

The radian is a unit of plane angle, equal to 180/π degrees, or about 57.2958 degrees.

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

Edited August 23, 2008 by monoceres

[muhmuuh]

I still think I'm right

The quote from wikipedia you posted shows how to convert radians to degrees.

I think I now how I'll prove I'm right.

Call you function with 1 degree. It will return 57.2958 radians so you say that 1degree=572958 radians which is not corret.

But if you call my function with 57.2958 it will return 1 radian so 57.2958 degrees = 1 radian which is correct

p.s.

just saw the wiki page about radians

http://en.wikipedia.org/wiki/Radian#Conversions

check after "Conversely, to convert from degrees to radians, multiply by π/180. For example,"

Edited August 23, 2008 by muhmuuh

[monoceres]

I still think I'm right

The quote from wikipedia you posted shows how to convert radians to degrees.

I think I now how I'll prove I'm right.

Call you function with 1 degree. It will return 57.2958 radians so you say that 1degree=572958 radians which is not corret.

But if you call my function with 57.2958 it will return 1 radian so 57.2958 degrees = 1 radian which is correct

p.s.

just saw the wiki page about radians

http://en.wikipedia.org/wiki/Radian#Conversions

check after "Conversely, to convert from degrees to radians, multiply by π/180. For example,"

Ok you're right.

Here's my attempt btw:

#include <iostream>
#include <cmath>

#define PI 3.14159

using namespace std;

double DegToRad(double);
double CalculateDistance(double,double,double,double,double);
int main ()
{
    double d;
    d=CalculateDistance(10,10,20,20,10000);
    cout << d;
    cin.get();  
    return 0;
}

double CalculateDistance(double a1,double a2,double b1,double b2,double radius)
{
       double dist;
       a1=DegToRad(a1);
       a2=DegToRad(a2);
       b1=DegToRad(b1);
       b2=DegToRad(b2);
       dist=2*radius*sqrt(pow(sin(pow(sin(a1-a2),2)/2+cos(a1)*cos(a2)*pow(sin(b1-b2),2)),-1));
       return dist;    
}

double DegToRad(double value)
{
       return value*(PI/180);
}

The weird thing is that it gives #INF when giving same coords.

When changing the -1 at the end to a 1 it gives zero so that could be right..

Edited August 23, 2008 by monoceres

[-Ultima-]

sin^-1(x) is arcsine... It's not the multiplicative inverse of sin(x). That's asin(x) in C++, not pow(sin(x), -1) or 1/sin(x). Assuming your formula is correct, the rough translation would be...

2*R*sqrt(asin( pow(sin(a1-a2), 2)/2 + cos(a1)*cos(a2)*pow(sin(b1-b2), 2)/2 ))

Thats after you change from degrees into radians for the input. Which reminds me: you might want to take a bit of effort to check to make sure inputs are within the domains of the respective functions before calling...

asin(x), -1 <= x <= 1

sqrt(x), 0 <= x

@monoceres: Well, if x^1 = 0, then x = 0, in which case x^(-1) = 1/0 = undefined, so it wouldn't really be all that weird, given that for the same point, sin(x-x) = sin(y-y) = sin(0) = 0, which would zero out the first parameter of your your pow(x, -1) call.

Edited August 23, 2008 by -Ultima-

[muhmuuh]

Hi!

I got it working

Here is the line:

d=2*R*asin(sqrt(sin((a1-a2)/2.0)*sin((a1-a2)/2.0)+cos(a1)*cos(a2)*sin((b1-b2)/2.0)*sin((b1-b2)/2.0)));

@Ultima

It is absolutely the same as yours after the edit

Thank you very much for the help

[JayBee]

Hi Muhmuuh,

You did a webmail cleaner app for me once morevthan few years ago (I attach screenshot of what it was), need you for quote on new little thing, and.....

......possibly bigger C programming  project later on, which we can also discuss. 

I need minor sorting and the difference computed of a 2 column tsv file and stop at the end, it would be similar small size project as that webmail cleaner you did for me.

Hope you can help by getting in contact with me soon, it's important thing, thx, Jerry

batchj2011 at gmail dot com

Only way I found you is by your ID on lower right corner of face of that software and googling this ID of yours. 

Thx, Truly, Jerry