Index: /branches/eam_branches/ipp-20101103/Ohana/src/relastro/src/high_speed_objects.c =================================================================== --- /branches/eam_branches/ipp-20101103/Ohana/src/relastro/src/high_speed_objects.c (revision 29817) +++ /branches/eam_branches/ipp-20101103/Ohana/src/relastro/src/high_speed_objects.c (revision 29818) @@ -13,13 +13,46 @@ // XXX seems silly to require region here just to find the catalog bounds / center - off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad; - int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid; + off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad,l,i1,Noff,Nmeasure,Naverage, NAVERAGE, NMEASURE; + int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid,Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj; double *X1, *Y1; - double dX, dY, dR, RADIUS2, yJ; + double dX, dY, dR, RADIUS2; Coords tcoords; + Catalog catalog1; int zcode, zNsec, ycode, yNsec, jcode, jNsec, hcode, hNsec, kcode, kNsec, USNO_R, USNO_N, Nsecfilt; + char filename[1024]; + char outdir[]="/data/ipp022.0/ndeacon/hispeedzy"; + Noff = strlen(CATDIR); + sprintf (filename, "%s/%s", outdir, &catalog[0].filename[Noff]); + dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/ + catalog1.filename = strcreate(filename); + + SIGMA_LIM = 0.0; Nsecfilt = GetPhotcodeNsecfilt(); + Naverage=catalog[0].Naverage; + Nmeasure=catalog[0].Nmeasure; + // ALLOCATE (catalog1.average, Average,Naverage); + // ALLOCATE (catalog1.measure, Measure,Nmeasure); + // ALLOCATE(catalog1.secfilt, SecFilt,catalog[0].Naverage*Nsecfilt); + catalog1.filename = filename; // based on the input name, need to keep everything below the catdir portion + catalog1.Nsecfilt = Nsecfilt; + // always load all of the data (if any exists) + catalog1.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + catalog1.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog1.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + + if (!dvo_catalog_open (&catalog1, region, VERBOSE,"w")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", + filename); + exit (2); + } + + NAVERAGE = 1000; + NMEASURE = 10000; + REALLOCATE (catalog1.average, Average, NAVERAGE); + REALLOCATE (catalog1.measure, Measure, NMEASURE); + REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt); ycode = GetPhotcodeCodebyName("y"); @@ -67,5 +100,5 @@ NgroupBbad = 0; - fprintf (stdout, "# RA_A DEC_A RA_B DEC_B : pmx pmy dR : z dz y dy J dJ H dH K dK\n"); + fprintf (stdout, "# RA_A DEC_A RA_B DEC_B : pmx pmy dR dt: z dz y dy J dJ H dH K dK\n"); //................270.2346670 20.2471540 270.2170434 20.2717396 : -59.11 88.78 106.66 : 0.000 0.000 19.582 0.112 16.041 0.110 15.388 0.098 14.858 0.001 @@ -88,4 +121,9 @@ foundA = FALSE; for (j = 0; !foundA && (j < catalog[0].average[i].Nmeasure); j++, m++) { + if((catalog[0].average[i].R>204.1923)&&(catalog[0].average[i].R<204.1924)&&(catalog[0].average[i].D>11.376)&&(catalog[0].average[i].D<11.377)) + { + printf("Hello"); + } + if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupA, NphotcodesGroupA)) { foundA = TRUE; @@ -97,4 +135,10 @@ foundB = FALSE; for (j = 0; !foundB && (j < catalog[0].average[i].Nmeasure); j++, m++) { + + if((catalog[0].average[i].R>204.192)&&(catalog[0].average[i].R<204.1925)&&(catalog[0].average[i].D>11.376)&&(catalog[0].average[i].D<11.377)) + { + printf("Hello"); + } + if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupB, NphotcodesGroupB)) { foundB = TRUE; @@ -114,13 +158,19 @@ if (foundA && !foundB) { // average-based tests: + + if((catalog[0].average[i].R>204.1923)&&(catalog[0].average[i].R<204.1924)&&(catalog[0].average[i].D>11.376)&&(catalog[0].average[i].D<11.377)) + { + printf("Hello"); + } + valid = TRUE; - 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); + valid &= ((catalog[0].average[i].flags & 0x02000000) == 0); + valid &= ((catalog[0].average[i].flags & 0x08000000) > 0); + valid &= ((catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + yNsec].M))); + valid &= (((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0))||(isnan(catalog[0].secfilt[i*Nsecfilt + zNsec].M))); // XXX color restrictions are applied in the pair-wise matching below // measure-based tests: - m = catalog[0].average[i].measureOffset; + /* 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) { @@ -130,5 +180,5 @@ valid &= ((catalog[0].measure[m].M > 18.0) || isnan(catalog[0].measure[m].M)); } - } + }*/ // ((objflags & 524288) == 524288) && @@ -154,19 +204,25 @@ // average-based tests: + if((catalog[0].average[i].R>204.192)&&(catalog[0].average[i].R<204.193)&&(catalog[0].average[i].D>11.372)&&(catalog[0].average[i].D<11.373)) + { + printf("Hello"); + } valid = TRUE; - 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].average[i].flags & 0x01000000) == 0); + + valid &= ((catalog[0].average[i].flags & 0x04000000) > 0); + + valid &= ((catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + jNsec].M))); valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].Nused > 1); valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].dM < 0.2); - if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) { + /*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; + /*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) { @@ -176,5 +232,5 @@ valid &= ((catalog[0].measure[m].M > 18.5) || isnan(catalog[0].measure[m].M)); } - } + }*/ // (abs(glat)>15.0) && @@ -236,4 +292,5 @@ Nmatch = 0; + Nmatchmeas = 0; RADIUS2 = SQ(RADIUS); @@ -243,5 +300,4 @@ ni = N1[i]; nj = N1[j]; - // if (ni == 131030) { // fprintf (stderr, "got 2mass\n"); @@ -268,5 +324,5 @@ // within match range; look for valid matches - for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) { + for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) { if (J == i) continue; // avoid auto-matches @@ -281,18 +337,53 @@ // y_groupA - J_groupB - yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M; + + /* 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; + if (yJ > 5.0) continue;*/ /*** a match is found (just print it for the moment) ***/ + /*Define a vector of matched indices and set averages in new catalogue*/ + Nepoch=2; + FIT_MODE = FIT_PM_ONLY; +/*nv=(int *)malloc(Nepoch*sizeof(int));*/ + nv[0]=ni; /*THESE SHOULD BE CHANGED FOR MULTIPLE EPOCHS AS SHOULD nv*/ + nv[1]=nj; + + catalog1.average[Nmatch]=catalog[0].average[nv[0]]; + /*Loop over index values and set measurements in new catalogue*/ + Nmatchmeasobj=0; + catalog1.average[Nmatch].measureOffset=Nmatchmeas; + for(l=0;l