Index: /trunk/psLib/src/types/psTree.c
===================================================================
--- /trunk/psLib/src/types/psTree.c	(revision 28231)
+++ /trunk/psLib/src/types/psTree.c	(revision 28232)
@@ -430,11 +430,10 @@
         switch (dim) {
           case 2: {
-              // Haversine formula
+              // Haversine formula, modulo a factor of 1/2
               double dphi = coords->data.F64[1] - tree->data->F64[index][1];
-              double sindphi = sin(dphi / 2.0);
+              double haverPhi = 1.0 - cos(dphi);
               double dlambda = coords->data.F64[0] - tree->data->F64[index][0];
-              double sindlambda = sin(dlambda / 2.0);
-              return PS_SQR(sindphi) +
-                  cos(coords->data.F64[1]) * cos(tree->data->F64[index][1]) * PS_SQR(sindlambda);
+              double haverLambda = 1.0 - cos(dlambda);
+              return haverPhi + cos(coords->data.F64[1]) * cos(tree->data->F64[index][1]) * haverLambda;
           }
           default:
@@ -446,10 +445,4 @@
 
     return NAN;
-}
-
-static inline double sin2diff(double coord1, double coord2)
-{
-    double sinDiff = sin(coord1 - coord2);
-    return PS_SQR(sinDiff);
 }
 
@@ -479,25 +472,19 @@
           double decMin = tree->min->F64[index][1], decMax = tree->max->F64[index][1]; // Dec bounds
 
-          double raDist = 0.0;
+          // Haversine formula, modulo a factor of 1/2
+          // This seems to deal with the wrap at RA=0=2pi, probably ascending the tree and descending on the
+          // other side of the wrap.
           if (ra < raMin || ra > raMax) {
-              double ra1 = sin2diff(ra, raMin);
-              double ra2 = sin2diff(ra + 2.0 * M_PI, raMin);
-              raDist = PS_MIN(ra1, ra2);
-              double ra3 = sin2diff(ra, raMax);
-              if (ra3 < raDist) {
-                  raDist = ra3;
-              }
-              double ra4 = sin2diff(ra + 2.0 * M_PI, raMax);
-              if (ra4 < raDist) {
-                  raDist = ra4;
-              }
+              double ra1 = cos(ra - raMin), ra2 = cos(ra - raMax); // Options for RA distance
+              double raDist = PS_MAX(ra1, ra2);                    // Will give the smallest distance
+              double dec0 = (fabs(dec - decMin) < fabs(dec - decMax)) ? decMin : decMax; // Closest Dec limit
+              distance += cos(dec0) * cos(dec) * (1.0 - raDist);
           }
 
-          double decDist = 0.0;
           if (dec < decMin || dec > decMax) {
-              double dec1 = sin2diff(dec, decMin);
-              double dec2 = sin2diff(dec, decMax);
-              decDist = PS_MIN(dec1, dec2);
+              double dec1 = cos(dec - decMin), dec2 = cos(dec - decMax); // Options for Dec distance
+              distance += 1.0 - PS_MAX(dec1, dec2);
           }
+
           break;
       }
@@ -619,9 +606,7 @@
         // Using the square of the distance as the distance measure
         return PS_SQR(radius);
-      case PS_TREE_SPHERICAL: {
-          // Using a rearrangement of the Haversine formula
-          double sindist = sin(radius / 2.0);
-          return PS_SQR(sindist);
-      }
+      case PS_TREE_SPHERICAL:
+        // Using Haversine formula, modulo a factor of 1/2
+        return 1.0 - cos(radius);
       default:
         psAbort("Unrecognised type: %x", tree->type);