Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 31602) +++ trunk/Ohana/src/relastro/Makefile (revision 31635) @@ -61,4 +61,5 @@ $(SRC)/high_speed_catalogs.$(ARCH).o \ $(SRC)/high_speed_objects.$(ARCH).o \ +$(SRC)/high_speed_utils.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 31602) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 31635) @@ -131,4 +131,6 @@ int NphotcodesGroupA, NphotcodesGroupB; PhotCode **photcodesGroupA, **photcodesGroupB; +char WHERE_A[10000], WHERE_B[10000]; +SkyRegionSelection SELECTION; int AreaSelect; @@ -361,2 +363,8 @@ int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit); + +int initializeConstraints(); +int applyConstraintsA(); +int applyConstraintsB(); +void setupAreaSelection(SkyRegion *region); + Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 31602) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 31635) @@ -48,4 +48,8 @@ ScanConfig(config, "PHOTCODE_FILE", "%s", 0, MasterPhotcodeFile); + // Defaults for WHERE_A are established in db_utils.c: intializeConstraints + ScanConfig(config, "WHERE_A", "%s", 0, WHERE_A); + ScanConfig(config, "WHERE_B", "%s", 0, WHERE_B); + sprintf (ImageCat, "%s/Images.dat", CATDIR); Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 31602) +++ trunk/Ohana/src/relastro/src/args.c (revision 31635) @@ -307,5 +307,5 @@ fprintf (stderr, " -update-chips\n"); fprintf (stderr, " -update-mosaics\n"); - fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius)\n\n"); + fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius) (whereA) (whereB) (output catdir)\n\n"); fprintf (stderr, " additional options: \n"); fprintf (stderr, " -time (start)(stop)\n"); Index: trunk/Ohana/src/relastro/src/high_speed_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/high_speed_catalogs.c (revision 31602) +++ trunk/Ohana/src/relastro/src/high_speed_catalogs.c (revision 31635) @@ -9,4 +9,7 @@ Catalog catalog; + // tell libdvo the CATDIR + dvo_set_catdir(CATDIR); + // load the current sky table (layout of all SkyRegions) sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE); @@ -16,7 +19,11 @@ if (UserCatalog) { skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC); + setupAreaSelection(skylist->regions[0]); } else { skylist = SkyListByPatch (sky, -1, &UserPatch); + setupAreaSelection(&UserPatch); } + + initializeConstraints(); // load data from each region file, only use bright stars Index: trunk/Ohana/src/relastro/src/high_speed_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 31602) +++ trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 31635) @@ -4,9 +4,4 @@ # define NEXT_J { j++; continue; } -// XXX these were used to spot-check specific objects -# 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) { @@ -14,5 +9,6 @@ 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; + int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj; + int foundA, foundB; double *X1, *Y1; double dX, dY, dR, RADIUS2; @@ -20,9 +16,9 @@ Catalog catalog1; - int zcode, zNsec, ycode, yNsec, jcode, jNsec, hcode, hNsec, kcode, kNsec, USNO_R, USNO_N, Nsecfilt; + int Nsecfilt; char filename[1024]; Noff = strlen(CATDIR); - sprintf (filename, "%s/%s", HIGH_SPEED_DIR, &catalog[0].filename[Noff]); - printf("%s\n",filename); + sprintf (filename, "%s%s", HIGH_SPEED_DIR, &catalog[0].filename[Noff]); + fprintf(stderr, "%s\n",filename); dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/ catalog1.filename = strcreate(filename); @@ -33,7 +29,4 @@ 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; @@ -56,22 +49,4 @@ REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt); - ycode = GetPhotcodeCodebyName("y"); - yNsec = GetPhotcodeNsec(ycode); - - zcode = GetPhotcodeCodebyName("z"); - zNsec = GetPhotcodeNsec(zcode); - - jcode = GetPhotcodeCodebyName("J"); - jNsec = GetPhotcodeNsec(jcode); - - hcode = GetPhotcodeCodebyName("H"); - hNsec = GetPhotcodeNsec(hcode); - - kcode = GetPhotcodeCodebyName("K"); - kNsec = GetPhotcodeNsec(kcode); - - USNO_R = GetPhotcodeCodebyName("USNO_RED"); - USNO_N = GetPhotcodeCodebyName("USNO_N"); - // high-speed between different surveys (easier case): // we have two sets of photcodes (A) and (B) and are looking for objects @@ -94,5 +69,5 @@ memset (groupB, 0, catalog[0].Naverage*sizeof(int)); - fprintf (stderr, "checking "OFF_T_FMT" objects\n", catalog[0].Naverage); + if (VERBOSE) fprintf (stderr, "checking "OFF_T_FMT" objects\n", catalog[0].Naverage); Nslow = 0; Ninvalid = 0; @@ -100,21 +75,6 @@ NgroupBbad = 0; - 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 - // mark (exclude) objects with both sets of target photcodes for (i = 0; i < catalog[0].Naverage; i++) { - - // 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? m = catalog[0].average[i].measureOffset; @@ -132,5 +92,4 @@ for (j = 0; !foundB && (j < catalog[0].average[i].Nmeasure); j++, m++) { - if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupB, NphotcodesGroupB)) { foundB = TRUE; @@ -138,4 +97,5 @@ } + // object found in both - mark as slow moving if (foundA && foundB) { slowMoving[i] = TRUE; @@ -144,102 +104,30 @@ } - // additional constraints: - // * group A : require - - // 2MASS detections: + // Apply additional constraints: if (foundA && !foundB) { - // average-based tests: - - - valid = TRUE; - 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; - 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) || isnan(catalog[0].measure[m].M)); - } - if (catalog[0].measure[m].photcode == USNO_N) { - valid &= ((catalog[0].measure[m].M > 18.0) || isnan(catalog[0].measure[m].M)); - } - }*/ - - // ((objflags & 524288) == 524288) && - // ((objflags & 262144) == 0) && - // (z<1.0) && - // (y<1.0) && - // ((yp[$i]-J)>1.25) && ((yp[$i]-J)<5.0) && - // ((USNO_RED<1.0)||(USNO_RED>20.0)) && - // ((USNO_N>18.0)||(USNO_N<1.0)) - - if (valid) { - groupA[i] = TRUE; - NgroupA ++; - continue; - } else { - NgroupAbad ++; - // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n", i); - } - } - - // PS1 detections: - if (foundB && !foundA) { - - // average-based tests: - valid = TRUE; - 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)||(catalog[0].secfilt[i*Nsecfilt + zNsec].Nused > 1)); - valid &= ((catalog[0].secfilt[i*Nsecfilt + yNsec].dM < 0.2)||(catalog[0].secfilt[i*Nsecfilt + zNsec].dM < 0.2)); - - /*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) || isnan(catalog[0].measure[m].M)); - } - if (catalog[0].measure[m].photcode == USNO_N) { - valid &= ((catalog[0].measure[m].M > 18.5) || isnan(catalog[0].measure[m].M)); - } - }*/ - - // (abs(glat)>15.0) && - // ((objflags & 65536) == 0) && - // ((objflags & 131072) == 131072) && - // (J<1.0) && - // (y:nphot>1) && - // (y:err<0.2) && - // ((USNO_RED>20.0)||(USNO_RED<1.0)) && - // ((USNO_N>18.5)||(USNO_N<1.0)) && - // ((z<1.0)||(z:nphot == 1)||((z:nphot>1)&&((z-y)>0.2))) - - if (valid) { - groupB[i] = TRUE; - NgroupB ++; - continue; - } else { - NgroupBbad ++; - // fprintf (stderr, "skip "OFF_T_FMT" (group B)\n", i); - } - } - - // this object does not have a detection from either groupA or groupB -- skip as if slow + if (applyConstraintsA(catalog, i)) { + groupA[i] = TRUE; + NgroupA++; + continue; + } else { + NgroupAbad ++; + // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n", i); + } + + } else if (foundB && !foundA) { + + if (applyConstraintsB(catalog, i)) { + groupB[i] = TRUE; + NgroupB++; + continue; + } else { + NgroupBbad++; + } + } + // this object does not have a detection matching the contraints from either gorupA or groupB -- skip as if slow slowMoving[i] = TRUE; Ninvalid ++; } + fprintf (stderr, OFF_T_FMT" slow, "OFF_T_FMT" invalid, ("OFF_T_FMT" group A, "OFF_T_FMT" group B), "OFF_T_FMT" total objects; "OFF_T_FMT" in group A, "OFF_T_FMT" in group B\n", Nslow, Ninvalid, NgroupAbad, NgroupBbad, catalog[0].Naverage, NgroupA, NgroupB); @@ -284,10 +172,4 @@ ni = N1[i]; nj = N1[j]; - // if (ni == 131030) { - // fprintf (stderr, "got 2mass\n"); - // } - // if (nj == 53581) { - // fprintf (stderr, "got ps1\n"); - // } if (slowMoving[ni]) NEXT_I; @@ -319,10 +201,4 @@ 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) ***/ @@ -330,5 +206,4 @@ 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; @@ -363,21 +238,4 @@ REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt); } - // for the moment, just write the matches to stdout: - //float pmx = -dX; // proper-motion displacement in ra direction (B - A) [arcsec] - //float pmy = -dY; // proper-motion displacement in dec direction (B - A) [arcsec] - /* fprintf (stdout, "%11.7f %11.7f %11.7f %11.7f : %7.2f %7.2f %7.2f: %7.3f %6.3f %7.3f %6.3f %7.3f %6.3f %7.3f %6.3f %7.3f %6.3f\n", - catalog[0].average[ni].R, catalog[0].average[ni].D, - catalog[0].average[nj].R, catalog[0].average[nj].D, - dX, dY, sqrt(dR), - catalog[0].secfilt[nj*Nsecfilt + zNsec].M, - catalog[0].secfilt[nj*Nsecfilt + zNsec].dM, - catalog[0].secfilt[nj*Nsecfilt + yNsec].M, - catalog[0].secfilt[nj*Nsecfilt + yNsec].dM, - catalog[0].secfilt[ni*Nsecfilt + jNsec].M, - catalog[0].secfilt[ni*Nsecfilt + jNsec].dM, - catalog[0].secfilt[ni*Nsecfilt + hNsec].M, - catalog[0].secfilt[ni*Nsecfilt + hNsec].dM, - catalog[0].secfilt[ni*Nsecfilt + kNsec].M, - catalog[0].secfilt[ni*Nsecfilt + kNsec].dM);*/ } i++; Index: trunk/Ohana/src/relastro/src/high_speed_utils.c =================================================================== --- trunk/Ohana/src/relastro/src/high_speed_utils.c (revision 31635) +++ trunk/Ohana/src/relastro/src/high_speed_utils.c (revision 31635) @@ -0,0 +1,139 @@ +# include "relastro.h" +# include "kapa.h" + + // !ID_OBJ_EXT_ALT ID_OBJ_GOOD_ALT +#define DEFAULT_WHERE_A "((flags & 0x02000000 == 0) && (flags & 0x08000000))" + + // !ID_OBJ_EXT ID_OBJ_GOOD +#define DEFAULT_WHERE_B "((flags & 0x01000000) == 0) && ((flags & 0x04000000) && ((y:nphot > 1) || (z:nphot > 1)) &&((y:err < 0.2) || (z:err < 0.2)))" + + +extern char **isolate_elements(int, char **, int *); + +int print_error(); + +static dbField *fieldsA = NULL; +static int NfieldsA = 0; +static dbStack *stackA = NULL; +static int NstackA = 0; +static dbField *fieldsB = NULL; +static int NfieldsB = 0; +static dbStack *stackB = NULL; +static int NstackB = 0; +static dbValue *valuesA; +// static int NvaluesA; +static dbValue *valuesB; +// static int NvaluesB; + +void setupAreaSelection(SkyRegion *region) { + SELECTION.name = NULL; + SELECTION.list = NULL; + SELECTION.useDisplay = FALSE; + SELECTION.useSkyregion = TRUE; + set_skyregion(region->Rmin, region->Rmax, region->Dmin, region->Dmax); +} + + +// split string into an array of words using space as delimiter by a space. +static char ** splitToWords(char *string, int *n) { + char **words; + int NwordsAllocated = 20; + ALLOCATE(words, char *, NwordsAllocated); + + char *copy_of_str = strdup(string); + + char *saveptr = NULL; + char *p = strtok_r(copy_of_str, " ", &saveptr); + *n = 0; + while (p) { + if (*n == NwordsAllocated) { + NwordsAllocated += 20; + REALLOCATE (words, char *, NwordsAllocated); + } + words[(*n)++] = p; + p = strtok_r(NULL, " ", &saveptr); + } + + return words; +} + +// Create the data strucutres required to apply avextract like cuts +int setupConstraints(char *whereString, dbField **pFields, int *pNfields, dbStack **pStack, int *pNstack) { + + // split the string into words + int nWords; + char **words = splitToWords(whereString, &nWords); + if (!nWords) { + fprintf(stderr, "unable to parse constraint\n"); + return FALSE; + } + + // parse it into elements of the where condition + char **cstack = NULL; + int Ncstack; + cstack = isolate_elements(nWords, words, &Ncstack); + free(words); + + *pStack = dbRPN(Ncstack, cstack, pNstack); + if (Ncstack && !*pNstack) { + print_error(); + return FALSE; + } + dbAstroRegionLimits(pStack, pNstack, &SELECTION, DVO_TABLE_AVERAGE); + + *pFields = NULL; + *pNfields = 0; + + return dbCheckStack(*pStack, *pNstack, DVO_TABLE_AVERAGE, pFields, pNfields); +} + +int initializeConstraints() { + // if constaints haven't been set by the user use 2mass + if (!WHERE_A[0]) { + strcpy(WHERE_A, DEFAULT_WHERE_A); + } + if (!WHERE_B[0]) { + strcpy(WHERE_B, DEFAULT_WHERE_B); + } + + if (! setupConstraints(WHERE_A, &fieldsA, &NfieldsA, &stackA, &NstackA)) { + print_error(); + fprintf(stderr, "failed to set up constratints for group A: %s\n", WHERE_A); + exit(1); + } + ALLOCATE(valuesA, dbValue, NfieldsA); + if (! setupConstraints(WHERE_B, &fieldsB, &NfieldsB, &stackB, &NstackB)) { + print_error(); + fprintf(stderr, "failed to set up constratints for group A: %s\n", WHERE_A); + exit(1); + } + ALLOCATE(valuesB, dbValue, NfieldsB); + return TRUE; +} + +// XXX: TODO: We should be able to combine fieldsA and fieldsB and do a single extract + +int applyConstraintsA(Catalog *catalog, off_t i) { + off_t m = catalog[0].average[i].measureOffset; + off_t n; + int Nsecfilt = GetPhotcodeNsecfilt(); + + dbExtractAveragesInit(); + for (n = 0; n < NfieldsA; n++) { + valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsA[n]); + } + return dbBooleanCond(stackA, NstackA, valuesA); +} + +int applyConstraintsB(Catalog *catalog, off_t i) { + off_t m = catalog[0].average[i].measureOffset; + int Nsecfilt = GetPhotcodeNsecfilt(); + + dbExtractAveragesInit(); + + off_t n; + for (n = 0; n < NfieldsB; n++) { + valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsB[n]); + } + return dbBooleanCond(stackB, NstackB, valuesB); +} Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 31602) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 31635) @@ -63,8 +63,8 @@ /* parse the comma-separated list of photcodesGroupA */ + fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); list = PHOTCODE_A_LIST; while ((codename = strtok_r (list, ",", &ptr)) != NULL) { list = NULL; // pass NULL on successive strtok_r calls - fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); fprintf (stderr, "codename: %s\n", codename); if ((photcodesGroupA[NphotcodesGroupA] = GetPhotcodebyName (codename)) == NULL) { @@ -83,9 +83,9 @@ ALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES); + fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); /* parse the comma-separated list of photcodesGroupB */ list = PHOTCODE_B_LIST; while ((codename = strtok_r (list, ",", &ptr)) != NULL) { list = NULL; // pass NULL on successive strtok_r calls - fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); fprintf (stderr, "codename: %s\n", codename); if ((photcodesGroupB[NphotcodesGroupB] = GetPhotcodebyName (codename)) == NULL) { Index: trunk/Ohana/src/relastro/src/relastroVisual.c =================================================================== --- trunk/Ohana/src/relastro/src/relastroVisual.c (revision 31602) +++ trunk/Ohana/src/relastro/src/relastroVisual.c (revision 31635) @@ -3,5 +3,4 @@ */ -#include "kapa.h" #include "relastro.h"