Changeset 28232

Timestamp:

Jun 4, 2010, 4:13:32 PM (16 years ago)

Author:

Paul Price

Message:

Dealing with the wrap at RA=0=2pi. Using 2*hav(theta) = 1-cos(theta) instead of hav(theta)=sin^{2(theta/2) which involves more multiplications.}

File:

• trunk/psLib/src/types/psTree.c (modified) (4 diffs)

trunk/psLib/src/types/psTree.c

@@ -430 +430 @@
         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 +445 @@
 
     return NAN;
-}
-
-static inline double sin2diff(double coord1, double coord2)
-{
-    double sinDiff = sin(coord1 - coord2);
-    return PS_SQR(sinDiff);
 }
 
@@ -479 +472 @@
           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 +606 @@
         // 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);