Changeset 28204

Timestamp:

Jun 3, 2010, 2:21:32 PM (16 years ago)

Author:

eugene

Message:

hard-wired rules for tdwarfs

File:

• trunk/Ohana/src/relastro/src/high_speed_objects.c (modified) (10 diffs)

trunk/Ohana/src/relastro/src/high_speed_objects.c

@@ -4 +4 @@
 # define NEXT_J { j++; continue; }
 
+# define REF_RA1 192.455781 
+# define REF_DEC1 2.563102 
+# define REF_RA2 192.444302113 
+# define REF_DEC2 2.5786559112 
 int high_speed_objects (SkyRegion *region, Catalog *catalog) {
 
   // XXX seems silly to require region here just to find the catalog bounds / center
 
-  off_t i, j, m, J, Ni, ni, nj, Nj, *N1, Nslow, NgroupA, NgroupB;
+  off_t i, j, m, J, ni, nj, *N1, Nslow, NgroupA, NgroupB;
   int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid;
   double *X1, *Y1;
-  double dX, dY, dR, RADIUS2;
+  double dX, dY, dR, RADIUS2, radius, yJ;
   Coords tcoords;
 
@@ -21 +25 @@
   yNsec = GetPhotcodeNsec(ycode);
   
-  zcode = GetPhotcodeCodebyName("y");
+  zcode = GetPhotcodeCodebyName("z");
   zNsec = GetPhotcodeNsec(zcode);
   
@@ -57 +61 @@
 
     if (i % 100 == 0) fprintf (stderr, ".");
+
+    radius = hypot((REF_DEC1 - catalog[0].average[i].D), (REF_RA1 - catalog[0].average[i].R));
+    if (radius < 0.0001) {
+      fprintf (stderr, "found it\n");
+    }
+    radius = hypot((REF_DEC2 - catalog[0].average[i].D), (REF_RA2 - catalog[0].average[i].R));
+    if (radius < 0.0001) {
+      fprintf (stderr, "found it\n");
+    }
 
     // do any of the measures for this object match group A?
@@ -89 +102 @@
       // average-based tests:
       valid = TRUE;
-      valid &= !(catalog[0].average[i].flags & 0x40000);
-      valid &= (catalog[0].average[i].flags & 0x80000);
+      valid &= ((catalog[0].average[i].flags & 0x40000) == 0);
+      valid &= ((catalog[0].average[i].flags & 0x80000)  > 0);
       valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0);
       valid &= (catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0);
@@ -99 +112 @@
       for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
         if (catalog[0].measure[m].photcode == USNO_R) {
-          valid &= ((catalog[0].measure[m].M > 20.0) || (catalog[0].measure[m].M < 1.0));
+          valid &= ((catalog[0].measure[m].M > 20.0) || isnan(catalog[0].measure[m].M));
         }
         if (catalog[0].measure[m].photcode == USNO_N) {
-          valid &= ((catalog[0].measure[m].M > 18.5) || (catalog[0].measure[m].M < 1.0));
+          valid &= ((catalog[0].measure[m].M > 18.0) || isnan(catalog[0].measure[m].M));
         }
       }
@@ -118 +131 @@
         NgroupA ++;
         continue;
+      } else {
+        // fprintf (stderr, "skip %lld (group A)\n", (long long) i);
       }
     }
@@ -126 +141 @@
       // average-based tests:
       valid = TRUE;
-      valid &= !(catalog[0].average[i].flags & 0x10000);
-      valid &= (catalog[0].average[i].flags & 0x20000);
+      valid &= ((catalog[0].average[i].flags & 0x10000) == 0);
+      valid &= ((catalog[0].average[i].flags & 0x20000)  > 0);
       valid &= (catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0);
       valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].Nused > 1);
       valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].dM < 0.2);
-      // XXX ignoring the z-band restriction for now
-      
+
+      if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) {
+        valid &= TRUE;
+      } else {
+        valid &= ((catalog[0].secfilt[i*Nsecfilt + zNsec].M - catalog[0].secfilt[i*Nsecfilt + yNsec].M) > 0.2);
+      }
+
       // measure-based tests:
       m = catalog[0].average[i].measureOffset;
       for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
         if (catalog[0].measure[m].photcode == USNO_R) {
-          valid &= ((catalog[0].measure[m].M > 20.0) || (catalog[0].measure[m].M < 1.0));
+          valid &= ((catalog[0].measure[m].M > 20.0) || isnan(catalog[0].measure[m].M));
         }
         if (catalog[0].measure[m].photcode == USNO_N) {
-          valid &= ((catalog[0].measure[m].M > 18.5) || (catalog[0].measure[m].M < 1.0));
+          valid &= ((catalog[0].measure[m].M > 18.5) || isnan(catalog[0].measure[m].M));
         }
       }
@@ -158 +178 @@
         NgroupB ++;
         continue;
+      } else {
+        // fprintf (stderr, "skip %lld (group B)\n", (long long) i);
       }
     }
@@ -207 +229 @@
     nj = N1[j];
 
+    if (ni == 131030) {
+      fprintf (stderr, "got 2mass\n");
+    }
+    if (nj == 53581) {
+      fprintf (stderr, "got ps1\n");
+    }
+
     if (slowMoving[ni]) NEXT_I;
     if (slowMoving[nj]) NEXT_J;
@@ -226 +255 @@
     for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {
       if (J == i) continue;  // avoid auto-matches
+
       dX = X1[i] - X1[J];
+
+      nj = N1[J];
+      if (!groupB[nj]) continue;
+
       dY = Y1[i] - Y1[J];
       dR = dX*dX + dY*dY;
       if (dR > RADIUS2) continue;
 
+      // y_groupA - J_groupB 
+      yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M;
+      if (yJ < 1.25) continue;
+      if (yJ > 5.0) continue;
+
       /*** a match is found (just print it for the moment) ***/
       Nmatch ++;
-      Ni = N1[i];
-      Nj = N1[J];
 
       if (1) fprintf (stderr, "match %d = %lld %lld : %f %f : %f %f : %f\n", 
-               Nmatch, (long long) i, (long long) J,
-               catalog[0].average[Ni].R, catalog[0].average[Ni].D,
-               catalog[0].average[Nj].R, catalog[0].average[Nj].D,
-               sqrt(dR));
+                      Nmatch, (long long) i, (long long) J,
+                      catalog[0].average[ni].R, catalog[0].average[ni].D,
+                      catalog[0].average[nj].R, catalog[0].average[nj].D,
+                      sqrt(dR));
     }
     i++;