Fix for ticket #71.

Problem was with the use of the haversine formula beyond it's valid
range of latitude/longitude pairs


git-svn-id: http://svn.osgeo.org/qgis/trunk@5296 c8812cc2-4d05-0410-92ff-de0c093fc19c

Author: g_j_m

src/core/qgsscalecalculator.cpp

@@ -79,23 +79,39 @@ double  QgsScaleCalculator::calculateGeographicDistance(QgsRect &mapExtent)
   // We'll use the middle latitude for the calculation
   // Note this is an approximation (although very close) but calculating scale
   // for geographic data over large extents is quasi-meaningless
-  double lat1 = (mapExtent.yMax() - mapExtent.yMin())/2 + mapExtent.yMin();
-  double lat2 = lat1;
-  double lon1 = mapExtent.xMin();
-  double lon2 = mapExtent.xMax();
-  double dlon = lon2 - lon1;
-  double dlat = lat2 - lat1;
-  double rads = (4 * atan(1.0))/180;
-  double a = pow((sin(dlat*rads/2)),2) + 
-    cos(lat1 *rads) * cos(lat2 *rads) *pow(sin(dlon*rads/2),2);
-  double c = 2 * atan2(sqrt(a), sqrt(1-a));
-  // calculate radius of earth
-  double ra = 6378;
-  // unused! double rb = 6357;
-  double e = .081082;
-  double R = ra* sqrt(1-pow(e,2))/(1 - pow(e,2)*pow(sin(lat1*rads),2));
-  double d = c *R; // kilometers;
-  double meters = d * 1000.0;
+
+  // The distance between two points on a sphere can be estimated
+  // using the Haversine formula. This gives the shortest distance
+  // between two points on the sphere. However, what we're after is
+  // the distance from the left of the given extent and the right of
+  // it. This is not necessarily the shortest distance between two
+  // points on a sphere.
+  //
+  // The code below uses the Haversine formula, but with some changes
+  // to cope with the above problem, and also to deal with the extent
+  // possibly extending beyond +/-180 degrees: 
+  //
+  // - Use the Halversine formula to calculate the distance from -90 to
+  //   +90 degrees at the mean latitude.
+  // - Scale this distance by the number of degrees between
+  //   mapExtent.xMin() and mapExtent.xMax();
+  // - For a slight improvemnt, allow for the ellipsoid shape of earth.
+
+
+  // For a longitude change of 180 degrees
+  double lat = (mapExtent.yMax() + mapExtent.yMin()) * 0.5;
+  const static double rads = (4.0 * atan(1.0))/180.0;
+  double a = pow(cos(lat * rads), 2);
+  double c = 2.0 * atan2(sqrt(a), sqrt(1.0-a));
+  const static double ra = 6378000; // [m]
+  const static double rb = 6357000; // [m]
+  // The eccentricity. This comes from sqrt(1.0 - rb*rb/(ra*ra));
+  const static double e = 0.0810820288;
+  double radius = ra * (1.0 - e*e) / 
+    pow(1.0 - e*e*sin(lat*rads)*sin(lat*rads), 1.5);
+  double meters = (mapExtent.xMax() - mapExtent.xMin()) / 180.0 * radius * c;
+
   QgsDebugMsg("Distance across map extent (m): " + QString::number(meters));
+
   return meters;
 }