IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 16, 2011, 11:18:00 AM (15 years ago)
Author:
bills
Message:

move the functions that implement the dvo "database lookups" (avextract, mextract, imextract) etc from opihi to libdvo so that they can be used in relastro -high-speed. Restructure relastro -high-speed to use configuration variables for the cuts on the two sets of photcodes. Default to ps1 versus 2mass as was previoulsy hard coded

Location:
trunk/Ohana/src/relastro
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r29001 r31635  
    6161$(SRC)/high_speed_catalogs.$(ARCH).o  \
    6262$(SRC)/high_speed_objects.$(ARCH).o  \
     63$(SRC)/high_speed_utils.$(ARCH).o  \
    6364$(SRC)/relastroVisual.$(ARCH).o
    6465

  • trunk/Ohana/src/relastro/include/relastro.h

    r31160 r31635  
    131131int           NphotcodesGroupA,  NphotcodesGroupB;
    132132PhotCode     **photcodesGroupA, **photcodesGroupB;
     133char          WHERE_A[10000],          WHERE_B[10000];
     134SkyRegionSelection SELECTION;
    133135
    134136int AreaSelect;
    … …  
    361363
    362364int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit);
     365
     366int initializeConstraints();
     367int applyConstraintsA();
     368int applyConstraintsB();
     369void setupAreaSelection(SkyRegion *region);
     370

  • trunk/Ohana/src/relastro/src/ConfigInit.c

    r30616 r31635  
    4848  ScanConfig(config, "PHOTCODE_FILE",          "%s",  0, MasterPhotcodeFile);
    4949
     50  // Defaults for WHERE_A are established in db_utils.c: intializeConstraints
     51  ScanConfig(config, "WHERE_A",                "%s",  0, WHERE_A);
     52  ScanConfig(config, "WHERE_B",                "%s",  0, WHERE_B);
     53
    5054  sprintf (ImageCat, "%s/Images.dat", CATDIR);
    5155

  • trunk/Ohana/src/relastro/src/args.c

    r31160 r31635  
    307307  fprintf (stderr, "  -update-chips\n");
    308308  fprintf (stderr, "  -update-mosaics\n");
    309   fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius)\n\n");
     309  fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius) (whereA) (whereB) (output catdir)\n\n");
    310310  fprintf (stderr, "  additional options: \n");
    311311  fprintf (stderr, "  -time (start)(stop)\n");

  • trunk/Ohana/src/relastro/src/high_speed_catalogs.c

    r31602 r31635  
    99  Catalog catalog;
    1010
     11  // tell libdvo the CATDIR
     12  dvo_set_catdir(CATDIR);
     13
    1114  // load the current sky table (layout of all SkyRegions)
    1215  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE);
    … …  
    1619  if (UserCatalog) {
    1720    skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC);
     21    setupAreaSelection(skylist->regions[0]);
    1822  } else {
    1923    skylist = SkyListByPatch (sky, -1, &UserPatch);
     24    setupAreaSelection(&UserPatch);
    2025  }
     26
     27  initializeConstraints();
    2128
    2229  // load data from each region file, only use bright stars

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

    r31160 r31635  
    44# define NEXT_J { j++; continue; }
    55
    6 // XXX these were used to spot-check specific objects
    7 # define REF_RA1 192.455781
    8 # define REF_DEC1 2.563102
    9 # define REF_RA2 192.444302113
    10 # define REF_DEC2 2.5786559112
    116int high_speed_objects (SkyRegion *region, Catalog *catalog) {
    127
    … …  
    149
    1510  off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad,l,i1,Noff,Nmeasure,Naverage, NAVERAGE, NMEASURE;
    16   int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid,Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj;
     11  int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj;
     12  int foundA, foundB;
    1713  double *X1, *Y1;
    1814  double dX, dY, dR, RADIUS2;
    … …  
    2016  Catalog catalog1;
    2117
    22   int zcode, zNsec, ycode, yNsec, jcode, jNsec, hcode, hNsec, kcode, kNsec, USNO_R, USNO_N, Nsecfilt;
     18  int Nsecfilt;
    2319  char filename[1024];
    2420  Noff = strlen(CATDIR);
    25   sprintf (filename, "%s/%s", HIGH_SPEED_DIR, &catalog[0].filename[Noff]);
    26   printf("%s\n",filename);
     21  sprintf (filename, "%s%s", HIGH_SPEED_DIR, &catalog[0].filename[Noff]);
     22  fprintf(stderr, "%s\n",filename);
    2723  dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/
    2824  catalog1.filename = strcreate(filename);
    … …  
    3329  Naverage=catalog[0].Naverage;
    3430  Nmeasure=catalog[0].Nmeasure;
    35   // ALLOCATE (catalog1.average, Average,Naverage);
    36   // ALLOCATE (catalog1.measure, Measure,Nmeasure);
    37   // ALLOCATE(catalog1.secfilt, SecFilt,catalog[0].Naverage*Nsecfilt);
    3831  catalog1.filename = filename; // based on the input name, need to keep everything below the catdir portion
    3932  catalog1.Nsecfilt  = Nsecfilt;
    … …  
    5649  REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
    5750
    58   ycode = GetPhotcodeCodebyName("y");
    59   yNsec = GetPhotcodeNsec(ycode);
    60  
    61   zcode = GetPhotcodeCodebyName("z");
    62   zNsec = GetPhotcodeNsec(zcode);
    63  
    64   jcode = GetPhotcodeCodebyName("J");
    65   jNsec = GetPhotcodeNsec(jcode);
    66  
    67   hcode = GetPhotcodeCodebyName("H");
    68   hNsec = GetPhotcodeNsec(hcode);
    69  
    70   kcode = GetPhotcodeCodebyName("K");
    71   kNsec = GetPhotcodeNsec(kcode);
    72  
    73   USNO_R = GetPhotcodeCodebyName("USNO_RED");
    74   USNO_N = GetPhotcodeCodebyName("USNO_N");
    75 
    7651  // high-speed between different surveys (easier case):
    7752  // we have two sets of photcodes (A) and (B) and are looking for objects
    … …  
    9469  memset (groupB, 0, catalog[0].Naverage*sizeof(int));
    9570
    96   fprintf (stderr, "checking "OFF_T_FMT" objects\n",  catalog[0].Naverage);
     71  if (VERBOSE) fprintf (stderr, "checking "OFF_T_FMT" objects\n",  catalog[0].Naverage);
    9772  Nslow = 0;
    9873  Ninvalid = 0;
    … …  
    10075  NgroupBbad = 0;
    10176
    102   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");
    103   //................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
    104 
    10577  // mark (exclude) objects with both sets of target photcodes
    10678  for (i = 0; i < catalog[0].Naverage; i++) {
    107 
    108     // if (i % 100 == 0) fprintf (stderr, ".");
    109 
    110     // radius = hypot((REF_DEC1 - catalog[0].average[i].D), (REF_RA1 - catalog[0].average[i].R));
    111     // if (radius < 0.0001) {
    112     //   fprintf (stderr, "found it\n");
    113     // }
    114     // radius = hypot((REF_DEC2 - catalog[0].average[i].D), (REF_RA2 - catalog[0].average[i].R));
    115     // if (radius < 0.0001) {
    116     //   fprintf (stderr, "found it\n");
    117     // }
    118 
    11979    // do any of the measures for this object match group A?
    12080    m = catalog[0].average[i].measureOffset;
    … …  
    13292    for (j = 0; !foundB && (j < catalog[0].average[i].Nmeasure); j++, m++) {
    13393         
    134 
    13594      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupB, NphotcodesGroupB)) {
    13695        foundB = TRUE;
    … …  
    13897    }
    13998
     99    // object found in both - mark as slow moving
    140100    if (foundA && foundB) {
    141101      slowMoving[i] = TRUE;
    … …  
    144104    }
    145105
    146     // additional constraints:
    147     // * group A : require
    148 
    149     // 2MASS detections:
     106    // Apply additional constraints:
    150107    if (foundA && !foundB) {
    151       // average-based tests:
    152          
    153 
    154       valid = TRUE;
    155       valid &= ((catalog[0].average[i].flags & 0x02000000) == 0);
    156       valid &= ((catalog[0].average[i].flags & 0x08000000)  > 0);
    157       valid &= ((catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + yNsec].M)));
    158       valid &= (((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0))||(isnan(catalog[0].secfilt[i*Nsecfilt + zNsec].M)));
    159       // XXX color restrictions are applied in the pair-wise matching below
    160      
    161       // measure-based tests:
    162       /*      m = catalog[0].average[i].measureOffset;
    163       for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
    164         if (catalog[0].measure[m].photcode == USNO_R) {
    165           valid &= ((catalog[0].measure[m].M > 20.0) || isnan(catalog[0].measure[m].M));
    166         }
    167         if (catalog[0].measure[m].photcode == USNO_N) {
    168           valid &= ((catalog[0].measure[m].M > 18.0) || isnan(catalog[0].measure[m].M));
    169         }
    170         }*/
    171 
    172       // ((objflags & 524288) == 524288) &&
    173       // ((objflags & 262144) == 0) &&
    174       // (z<1.0) &&
    175       // (y<1.0) &&
    176       // ((yp[$i]-J)>1.25) && ((yp[$i]-J)<5.0) &&
    177       // ((USNO_RED<1.0)||(USNO_RED>20.0)) &&
    178       // ((USNO_N>18.0)||(USNO_N<1.0))
    179 
    180       if (valid) {
    181         groupA[i] = TRUE;
    182         NgroupA ++;
    183         continue;
    184       } else {
    185         NgroupAbad ++;
    186         // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n",  i);
    187       }
    188     }
    189 
    190     // PS1 detections:
    191     if (foundB && !foundA) {
    192 
    193       // average-based tests:
    194       valid = TRUE;
    195       valid &= ((catalog[0].average[i].flags & 0x01000000) == 0);
    196 
    197       valid &= ((catalog[0].average[i].flags & 0x04000000)  > 0);
    198 
    199       valid &= ((catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + jNsec].M)));
    200       valid &= ((catalog[0].secfilt[i*Nsecfilt + yNsec].Nused > 1)||(catalog[0].secfilt[i*Nsecfilt + zNsec].Nused > 1));
    201       valid &= ((catalog[0].secfilt[i*Nsecfilt + yNsec].dM < 0.2)||(catalog[0].secfilt[i*Nsecfilt + zNsec].dM < 0.2));
    202 
    203       /*if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) {
    204         valid &= TRUE;
    205       } else {
    206         valid &= ((catalog[0].secfilt[i*Nsecfilt + zNsec].M - catalog[0].secfilt[i*Nsecfilt + yNsec].M) > 0.2);
    207         }*/
    208 
    209       // measure-based tests:
    210       /*m = catalog[0].average[i].measureOffset;
    211       for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
    212         if (catalog[0].measure[m].photcode == USNO_R) {
    213           valid &= ((catalog[0].measure[m].M > 20.0) || isnan(catalog[0].measure[m].M));
    214         }
    215         if (catalog[0].measure[m].photcode == USNO_N) {
    216           valid &= ((catalog[0].measure[m].M > 18.5) || isnan(catalog[0].measure[m].M));
    217         }
    218         }*/
    219 
    220       // (abs(glat)>15.0) &&
    221       // ((objflags & 65536) == 0) &&
    222       // ((objflags & 131072) == 131072) &&
    223       // (J<1.0) &&
    224       // (y:nphot>1) &&
    225       // (y:err<0.2) &&
    226       // ((USNO_RED>20.0)||(USNO_RED<1.0)) &&
    227       // ((USNO_N>18.5)||(USNO_N<1.0)) &&
    228       // ((z<1.0)||(z:nphot == 1)||((z:nphot>1)&&((z-y)>0.2)))
    229 
    230       if (valid) {
    231         groupB[i] = TRUE;
    232         NgroupB ++;
    233         continue;
    234       } else {
    235         NgroupBbad ++;
    236         // fprintf (stderr, "skip "OFF_T_FMT" (group B)\n",  i);
    237       }
    238     }
    239 
    240     // this object does not have a detection from either groupA or groupB -- skip as if slow
     108        if (applyConstraintsA(catalog, i)) {
     109            groupA[i] = TRUE;
     110            NgroupA++;
     111            continue;
     112        } else {
     113            NgroupAbad ++;
     114            // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n",  i);
     115        }
     116
     117    } else if (foundB && !foundA) {
     118
     119        if (applyConstraintsB(catalog, i)) {
     120            groupB[i] = TRUE;
     121            NgroupB++;
     122            continue;
     123        } else {
     124            NgroupBbad++;
     125        }
     126    }
     127    // this object does not have a detection matching the contraints from either gorupA or groupB -- skip as if slow
    241128    slowMoving[i] = TRUE;
    242129    Ninvalid ++;
    243130  }
     131
    244132  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);
    245133   
    … …  
    284172    ni = N1[i];
    285173    nj = N1[j];
    286     // if (ni == 131030) {
    287     //   fprintf (stderr, "got 2mass\n");
    288     // }
    289     // if (nj == 53581) {
    290     //   fprintf (stderr, "got ps1\n");
    291     // }
    292174
    293175    if (slowMoving[ni]) NEXT_I;
    … …  
    319201      dR = dX*dX + dY*dY;
    320202      if (dR > RADIUS2) continue;
    321 
    322       // y_groupA - J_groupB
    323 
    324       /*      yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M;
    325       if (yJ < 1.25) continue;
    326       if (yJ > 5.0) continue;*/
    327203
    328204      /*** a match is found (just print it for the moment) ***/
    … …  
    330206      Nepoch=2;
    331207      FIT_MODE = FIT_PM_ONLY;
    332 /*nv=(int *)malloc(Nepoch*sizeof(int));*/
    333208      nv[0]=ni; /*THESE SHOULD BE CHANGED FOR MULTIPLE EPOCHS AS SHOULD nv*/
    334209      nv[1]=nj;
    … …  
    363238        REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
    364239      }
    365       // for the moment, just write the matches to stdout:
    366       //float pmx = -dX; // proper-motion displacement in ra  direction (B - A) [arcsec]
    367       //float pmy = -dY; // proper-motion displacement in dec direction (B - A) [arcsec]
    368       /*           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",
    369                catalog[0].average[ni].R, catalog[0].average[ni].D,
    370                catalog[0].average[nj].R, catalog[0].average[nj].D,
    371                dX, dY, sqrt(dR),
    372                catalog[0].secfilt[nj*Nsecfilt + zNsec].M,
    373                catalog[0].secfilt[nj*Nsecfilt + zNsec].dM,
    374                catalog[0].secfilt[nj*Nsecfilt + yNsec].M,
    375                catalog[0].secfilt[nj*Nsecfilt + yNsec].dM,
    376                catalog[0].secfilt[ni*Nsecfilt + jNsec].M,
    377                catalog[0].secfilt[ni*Nsecfilt + jNsec].dM,
    378                catalog[0].secfilt[ni*Nsecfilt + hNsec].M,
    379                catalog[0].secfilt[ni*Nsecfilt + hNsec].dM,
    380                catalog[0].secfilt[ni*Nsecfilt + kNsec].M,
    381                catalog[0].secfilt[ni*Nsecfilt + kNsec].dM);*/
    382240    }
    383241    i++;

  • trunk/Ohana/src/relastro/src/initialize.c

    r29001 r31635  
    6363
    6464    /* parse the comma-separated list of photcodesGroupA */
     65    fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
    6566    list = PHOTCODE_A_LIST;
    6667    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
    6768      list = NULL; // pass NULL on successive strtok_r calls
    68       fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
    6969      fprintf (stderr, "codename: %s\n", codename);
    7070      if ((photcodesGroupA[NphotcodesGroupA] = GetPhotcodebyName (codename)) == NULL) {
    … …  
    8383    ALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES);
    8484
     85    fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
    8586    /* parse the comma-separated list of photcodesGroupB */
    8687    list = PHOTCODE_B_LIST;
    8788    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
    8889      list = NULL; // pass NULL on successive strtok_r calls
    89       fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
    9090      fprintf (stderr, "codename: %s\n", codename);
    9191      if ((photcodesGroupB[NphotcodesGroupB] = GetPhotcodebyName (codename)) == NULL) {

  • trunk/Ohana/src/relastro/src/relastroVisual.c

    r30616 r31635  
    33 */
    44
    5 #include "kapa.h"
    65#include "relastro.h"
    76
Note: See TracChangeset for help on using the changeset viewer.