Average distance in latitude/longitude cell using Haversine formula

[caratacus]

Homework Statement

I've constructed a 5°x5° latitude/longitude cell, from 40-45° N and 85-90°W. This puts the center somewhere near the southern tip of Lake Michigan. I'm trying to find the average distance from the center of that cell (42.5°N, 87.5°W) to any other point in that cell.

Homework Equations

I know that the Haversine formula is used to calculate great circle distances between points. I formatted it nicely and placed it into Mathematica:

radians = Pi/180;
startLong = 87.5;
startLat = 42.5;
dlon = startLong - x;
dlat = startLat - y;
a = Sin[dlat/2*radians]^2 + 
   Cos[y*radians]*Cos[startLat*radians]*Sin[dlon/2*radians]^2;
c = 2*ArcTan[Sqrt[a], Sqrt[1 - a]];
d = 6370*c;

The Attempt at a Solution

My thought has been take the Haversine formula and take a double integral, with respect to x from -90° to -85°, then with respect to y from 40° to 45°, where x is longitude and y is latitude. However, this integral can crash Mathematica without breaking a sweat, and the NIntegrate function on Mathematica fails to approximate it correctly; I know because I've changed the starting point to a set of coordinates thousands of miles away, yet the integral approximation hardly changes its value.

In Mathematica,

NIntegrate[d, {y, 40*radians, 45*radians}, {x, 85*radians, 90*radians}]

yields 73.2622. This might be plausible, though when I change the startLat variable to 12.5, the same integral yields 72.7352. This means the average distance from somewhere near Central America to any point in the cell I've constructed is ~73km. Since this is obviously false, I have to conclude this double integral won't yield a proper answer, or there must be something I've overlooked. Any thoughts?

 

[voko]

The haversine formula does not involve arctan. There is a great-circle distance formula with arctan, but it is also different from the one you have.

http://en.wikipedia.org/wiki/Great-circle_distance

Note the integral won't give you the average; you need to divide the value of the integral by the area of your cell.

 

[caratacus]

Thanks for the reply! Here's a short modification to the equation, then, using arcsine:

radians = Pi/180;
startLong = 87.5;
startLat = 42.5;
dlon = startLong - x;
dlat = startLat - y;
a = Sin[dlat/2*radians]^2 + 
   Cos[y*radians]*Cos[startLat*radians]*Sin[dlon/2*radians]^2;
c = 2*ArcSin[Sqrt[a]];
d = 6370*c;

The same integral approximation yields 79.1375km, and substituting the same 12.5° in for the startLat and computing the integral has almost the exact same answer.

Is this integral the proper way of going about solving this problem? I'm not going to bother dividing the answer by ~200,000 sq. km. or whatever the area of the cell is (I'd have to check on that) if this method isn't yielding reasonable answers to begin with.

 

[voko]

The integral, in principle, is a correct approach. Can you check that the code, prior to integration, returns reasonable values?

 

[caratacus]

Using the second version with arcsine instead of arctangent,

x = 15; y = 10; d

yields 7825.55km (distance from 10°N, 15°W to 42°N, 87.5°W...why I have x for longitude and y for latitude is beyond me). This is only a few km off of an online great circle distance calculator (http://williams.best.vwh.net/gccalc.htm , for one). The equation looks fine...I'm starting to wonder if Mathematica doesn't approximate convoluted trigonometric integrals very well.

 

[voko]

Here is one potential problem: it seems you do integration with NIntegrate[d, {y, 40*radians, 45*radians}, {x, 85*radians, 90*radians}] - note that you specify the x/y limits in radians. Yet inside the haversine code you assume x and y are in degrees. Check that the angle measure is consistent throughout.

 

[caratacus]

Aha, thank you! That helps a lot; the integral has a reasonable amount of variation when I integrate with respect to different values of x and y. However, there still seems to be something wrong:

NIntegrate[d, {y, 40, 45}, {x, 85, 90}]

now returns 4650.41 km. This can't be the average distance to a point within the cell, as one of the longest distances (to 45°N, 90°W) is around 342km away. Any other ideas?

 

[voko]

As I remarked earlier, the integral gives you the "sum" of the all distances. To convert that to average, divide by the area of the integration rectangle; in this case 5x5 = 25, thus the average is 186 km.

 

[caratacus]

Awesome. Thank you so much for you help!