Index: branches/eam_branches/ipp-20151113/Ohana/src/gcompare/src/gc_compare.c
===================================================================
--- branches/eam_branches/ipp-20151113/Ohana/src/gcompare/src/gc_compare.c	(revision 39266)
+++ branches/eam_branches/ipp-20151113/Ohana/src/gcompare/src/gc_compare.c	(revision 39316)
@@ -2,5 +2,56 @@
 # define D_NMATCH 500;
 
-double gcdist (double r1, double d1, double r2, double d2) {
+typedef struct {
+  double sn_r;
+  double cs_r;
+  double sn_d;
+  double cs_d;
+} gc_data_type;
+
+gc_data_type *gcdist_init (data_type data) {
+
+  gc_data_type *gcdata;
+  ALLOCATE (gcdata, gc_data_type, data.Nvalues);
+
+  int i;
+  for (i = 0; i < data.Nvalues; i++) {
+    gcdata[i].sn_r = sin(RAD_DEG*data.values[i].X);
+    gcdata[i].cs_r = cos(RAD_DEG*data.values[i].X);
+    gcdata[i].sn_d = sin(RAD_DEG*data.values[i].Y);
+    gcdata[i].cs_d = cos(RAD_DEG*data.values[i].Y);
+  }
+  return gcdata;
+}
+
+double gcdist_new (gc_data_type v1, gc_data_type v2) {
+
+  double t1 = SQ(v2.cs_d*(v2.sn_r*v1.cs_r - v1.sn_r*v2.cs_r));
+
+  double t2 = v1.cs_d*v2.sn_d;
+
+  double cs_dr = (v1.cs_r*v2.cs_r + v1.sn_r*v2.sn_r); // cos(r2 - r1)
+  double t3 = v1.sn_d*v2.cs_d * cs_dr;
+
+  double t4 = SQ(t2 - t3);
+
+  double y = sqrt(t1 + t4);
+
+  double x = v1.sn_d * v2.sn_d + v1.cs_d * v2.cs_d * cs_dr;
+
+  /*
+  fprintf (stderr, "t1: %f\n", t1);
+  fprintf (stderr, "t2: %f\n", t2);
+  fprintf (stderr, "t3: %f\n", t3);
+  fprintf (stderr, "t4: %f\n", t4);
+  fprintf (stderr, "cs: %f\n", cs_dr);
+
+  fprintf (stderr, "x: %f\n", x);
+  fprintf (stderr, "y: %f\n", y);
+  */
+
+  return(DEG_RAD*atan2(y,x));
+}
+
+double gcdist_old (double r1, double d1, double r2, double d2) {
   double num,den;
   r1 *= M_PI / 180;
@@ -9,11 +60,36 @@
   d2 *= M_PI / 180;
 
-  num = sqrt(pow((cos(d2) * sin(r2 - r1)),2) +
+# if (0)  
+  double t1 = pow((cos(d2) * sin(r2 - r1)),2);
+
+  double t2 = cos(d1) * sin(d2);
+  
+  double cs_dr = cos(r2 - r1);
+  double t3 = sin(d1) * cos(d2) * cs_dr;
+
+  double t4 = pow((t2 - t3), 2);
+
+  double y = sqrt(t1 + t4);
+# endif
+
+  num = sqrt(pow((cos(d2) * sin(r2 - r1)),2) + 
 	     pow((cos(d1) * sin(d2) -
 		  sin(d1) * cos(d2) * cos(r2 - r1)),2));
   den = (sin(d1) * sin(d2) + cos(d1) * cos(d2) * cos(r2 - r1));
+
+# if (0)
+  fprintf (stderr, "t1: %f\n", t1);
+  fprintf (stderr, "t2: %f\n", t2);
+  fprintf (stderr, "t3: %f\n", t3);
+  fprintf (stderr, "t4: %f\n", t4);
+  fprintf (stderr, "cs: %f\n", cs_dr);
+
+  fprintf (stderr, "y: %f\n", y);
+
+  fprintf (stderr, "%f %f\n", num, den);
+# endif
+
   return(atan2(num,den) * (180 / M_PI));
 }
-  
 
 match_type *gc_compare (data_type data1, data_type data2, int *Nmatches, double radius, double DX, double DY, double noauto) {
@@ -28,4 +104,7 @@
   ALLOCATE (match, match_type, NMATCH);
 
+  gc_data_type *gcdata1 = gcdist_init (data1);
+  gc_data_type *gcdata2 = gcdist_init (data2);
+
   for (i = 0; i < data1.Nvalues ;i++) {
     if (!(i % 100))
@@ -33,6 +112,18 @@
 
     for (j = 0; j < data2.Nvalues ; j++) {
-      dR = gcdist(data1.values[i].X,data1.values[i].Y,
-		  data2.values[j].X,data2.values[j].Y);
+
+      if (1) {
+	dR = gcdist_new(gcdata1[i], gcdata2[j]);
+      } else {
+	dR = gcdist_old(data1.values[i].X, data1.values[i].Y, data2.values[j].X, data2.values[j].Y);
+      }
+
+      /* 
+      double dR1 = gcdist_new(gcdata1[i], gcdata2[j]);
+      double dR2 = gcdist_old(data1.values[i].X, data1.values[i].Y, data2.values[j].X, data2.values[j].Y);
+      fprintf (stderr, "%f %f\n", dR1, dR2);
+      dR = dR1;
+      */
+
 /*       fprintf(stderr,"%g %g %g %g => %g\n", */
 /* 	      data1.values[i].X,data1.values[i].Y,data2.values[j].X,data2.values[j].Y, */