IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29818


Ignore:
Timestamp:
Nov 24, 2010, 12:04:14 PM (15 years ago)
Author:
eugene
Message:

write results to separate catalog

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101103/Ohana/src/relastro/src/high_speed_objects.c

    r28811 r29818  
    1313  // XXX seems silly to require region here just to find the catalog bounds / center
    1414
    15   off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad;
    16   int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid;
     15  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;
    1717  double *X1, *Y1;
    18   double dX, dY, dR, RADIUS2, yJ;
     18  double dX, dY, dR, RADIUS2;
    1919  Coords tcoords;
     20  Catalog catalog1;
    2021
    2122  int zcode, zNsec, ycode, yNsec, jcode, jNsec, hcode, hNsec, kcode, kNsec, USNO_R, USNO_N, Nsecfilt;
     23  char filename[1024];
     24  char outdir[]="/data/ipp022.0/ndeacon/hispeedzy";
     25  Noff = strlen(CATDIR);
     26  sprintf (filename, "%s/%s", outdir, &catalog[0].filename[Noff]);
     27  dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/
     28  catalog1.filename = strcreate(filename);
     29
     30  SIGMA_LIM = 0.0;
    2231
    2332  Nsecfilt = GetPhotcodeNsecfilt();
     33  Naverage=catalog[0].Naverage;
     34  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);
     38  catalog1.filename = filename; // based on the input name, need to keep everything below the catdir portion
     39  catalog1.Nsecfilt  = Nsecfilt;
     40  // always load all of the data (if any exists)
     41  catalog1.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     42 
     43  catalog1.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
     44  catalog1.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
     45
     46  if (!dvo_catalog_open (&catalog1, region, VERBOSE,"w")) {
     47    fprintf (stderr, "ERROR: failure to open catalog file %s\n",
     48             filename);
     49    exit (2);
     50  }
     51
     52  NAVERAGE = 1000;
     53  NMEASURE = 10000;
     54  REALLOCATE (catalog1.average, Average, NAVERAGE);
     55  REALLOCATE (catalog1.measure, Measure, NMEASURE);
     56  REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
    2457
    2558  ycode = GetPhotcodeCodebyName("y");
    … …  
    67100  NgroupBbad = 0;
    68101
    69   fprintf (stdout, "#      RA_A  DEC_A              RA_B  DEC_B      :     pmx     pmy      dR  :    z      dz       y      dy       J      dJ       H      dH       K      dK\n");
     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");
    70103  //................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
    71104
    … …  
    88121    foundA = FALSE;
    89122    for (j = 0; !foundA && (j < catalog[0].average[i].Nmeasure); j++, m++) {
     123      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))
     124        {
     125          printf("Hello");
     126        }
     127
    90128      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupA, NphotcodesGroupA)) {
    91129        foundA = TRUE;
    … …  
    97135    foundB = FALSE;
    98136    for (j = 0; !foundB && (j < catalog[0].average[i].Nmeasure); j++, m++) {
     137                                                   
     138      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))
     139        {
     140          printf("Hello");
     141        }
     142
    99143      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupB, NphotcodesGroupB)) {
    100144        foundB = TRUE;
    … …  
    114158    if (foundA && !foundB) {
    115159      // average-based tests:
     160                                                   
     161      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))
     162        {
     163          printf("Hello");
     164        }
     165
    116166      valid = TRUE;
    117       valid &= ((catalog[0].average[i].flags & 0x40000) == 0);
    118       valid &= ((catalog[0].average[i].flags & 0x80000)  > 0);
    119       valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0);
    120       valid &= (catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0);
     167      valid &= ((catalog[0].average[i].flags & 0x02000000) == 0);
     168      valid &= ((catalog[0].average[i].flags & 0x08000000)  > 0);
     169      valid &= ((catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + yNsec].M)));
     170      valid &= (((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0))||(isnan(catalog[0].secfilt[i*Nsecfilt + zNsec].M)));
    121171      // XXX color restrictions are applied in the pair-wise matching below
    122172     
    123173      // measure-based tests:
    124       m = catalog[0].average[i].measureOffset;
     174      /*      m = catalog[0].average[i].measureOffset;
    125175      for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
    126176        if (catalog[0].measure[m].photcode == USNO_R) {
    … …  
    130180          valid &= ((catalog[0].measure[m].M > 18.0) || isnan(catalog[0].measure[m].M));
    131181        }
    132       }
     182        }*/
    133183
    134184      // ((objflags & 524288) == 524288) &&
    … …  
    154204
    155205      // average-based tests:
     206      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))
     207        {
     208          printf("Hello");
     209        }
    156210      valid = TRUE;
    157       valid &= ((catalog[0].average[i].flags & 0x10000) == 0);
    158       valid &= ((catalog[0].average[i].flags & 0x20000)  > 0);
    159       valid &= (catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0);
     211      valid &= ((catalog[0].average[i].flags & 0x01000000) == 0);
     212
     213      valid &= ((catalog[0].average[i].flags & 0x04000000)  > 0);
     214
     215      valid &= ((catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + jNsec].M)));
    160216      valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].Nused > 1);
    161217      valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].dM < 0.2);
    162218
    163       if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) {
     219      /*if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) {
    164220        valid &= TRUE;
    165221      } else {
    166222        valid &= ((catalog[0].secfilt[i*Nsecfilt + zNsec].M - catalog[0].secfilt[i*Nsecfilt + yNsec].M) > 0.2);
    167       }
     223        }*/
    168224
    169225      // measure-based tests:
    170       m = catalog[0].average[i].measureOffset;
     226      /*m = catalog[0].average[i].measureOffset;
    171227      for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
    172228        if (catalog[0].measure[m].photcode == USNO_R) {
    … …  
    176232          valid &= ((catalog[0].measure[m].M > 18.5) || isnan(catalog[0].measure[m].M));
    177233        }
    178       }
     234        }*/
    179235
    180236      // (abs(glat)>15.0) &&
    … …  
    236292
    237293  Nmatch = 0;
     294  Nmatchmeas = 0;
    238295  RADIUS2 = SQ(RADIUS);
    239296
    … …  
    243300    ni = N1[i];
    244301    nj = N1[j];
    245 
    246302    // if (ni == 131030) {
    247303    //   fprintf (stderr, "got 2mass\n");
    … …  
    268324
    269325    // within match range; look for valid matches
    270     for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {
     326    for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {     
    271327      if (J == i) continue;  // avoid auto-matches
    272328
    … …  
    281337
    282338      // y_groupA - J_groupB
    283       yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M;
     339
     340      /*      yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M;
    284341      if (yJ < 1.25) continue;
    285       if (yJ > 5.0) continue;
     342      if (yJ > 5.0) continue;*/
    286343
    287344      /*** a match is found (just print it for the moment) ***/
     345      /*Define a vector of matched indices and set averages in new catalogue*/
     346      Nepoch=2;
     347      FIT_MODE = FIT_PM_ONLY;
     348/*nv=(int *)malloc(Nepoch*sizeof(int));*/
     349      nv[0]=ni; /*THESE SHOULD BE CHANGED FOR MULTIPLE EPOCHS AS SHOULD nv*/
     350      nv[1]=nj;
     351
     352      catalog1.average[Nmatch]=catalog[0].average[nv[0]];
     353      /*Loop over index values and set measurements in new catalogue*/
     354      Nmatchmeasobj=0;
     355      catalog1.average[Nmatch].measureOffset=Nmatchmeas;
     356      for(l=0;l<Nepoch;l++) {
     357          m = catalog[0].average[nv[l]].measureOffset;
     358          for(i1=0;i1<catalog[0].average[nv[l]].Nmeasure;i1++) {
     359              catalog1.measure[Nmatchmeas]=catalog[0].measure[m+i1];
     360              /*Set offset RA and Dec wrt correct average value*/
     361              catalog1.measure[Nmatchmeas].dR=catalog[0].measure[m+i1].dR+3600.0*(catalog[0].average[nv[0]].R-catalog[0].average[nv[l]].R);
     362              catalog1.measure[Nmatchmeas].dD=catalog[0].measure[m+i1].dD+3600.0*(catalog[0].average[nv[0]].D-catalog[0].average[nv[l]].D);
     363              catalog1.measure[Nmatchmeas].averef = Nmatch;
     364              Nmatchmeasobj++;
     365              Nmatchmeas++;
     366
     367              if (Nmatchmeas == NMEASURE - 1) {
     368                NMEASURE += 10000;
     369                REALLOCATE (catalog1.measure, Measure, NMEASURE);
     370              }
     371          }
     372      }
     373      catalog1.average[Nmatch].Nmeasure=Nmatchmeasobj;
    288374      Nmatch ++;
    289375
     376      if (Nmatch == NAVERAGE - 1) {
     377        NAVERAGE += 1000;
     378        REALLOCATE (catalog1.average, Average, NAVERAGE);
     379        REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
     380      }
    290381      // for the moment, just write the matches to stdout:
    291       float pmx = -dX; // proper-motion displacement in ra  direction (B - A) [arcsec]
    292       float pmy = -dY; // proper-motion displacement in dec direction (B - A) [arcsec]
    293       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",
     382      //float pmx = -dX; // proper-motion displacement in ra  direction (B - A) [arcsec]
     383      //float pmy = -dY; // proper-motion displacement in dec direction (B - A) [arcsec]
     384      /*           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",
    294385               catalog[0].average[ni].R, catalog[0].average[ni].D,
    295386               catalog[0].average[nj].R, catalog[0].average[nj].D,
    296                pmx, pmy, sqrt(dR),
     387               dX, dY, sqrt(dR),
    297388               catalog[0].secfilt[nj*Nsecfilt + zNsec].M,
    298389               catalog[0].secfilt[nj*Nsecfilt + zNsec].dM,
    … …  
    304395               catalog[0].secfilt[ni*Nsecfilt + hNsec].dM,
    305396               catalog[0].secfilt[ni*Nsecfilt + kNsec].M,
    306                catalog[0].secfilt[ni*Nsecfilt + kNsec].dM);
     397               catalog[0].secfilt[ni*Nsecfilt + kNsec].dM);*/
    307398    }
    308399    i++;
    309400  }
     401  catalog1.Naverage=Nmatch;
     402  catalog1.Nmeasure=Nmatchmeas;
     403  catalog1.Nsecfilt=Nsecfilt;
     404  catalog1.Nsecf_mem=Nmatch*Nsecfilt;
     405
     406  UpdateObjects(&catalog1, 1);
     407
    310408  fprintf (stderr, "found %d matches\n", Nmatch);
    311  
     409  dvo_catalog_save (&catalog1, VERBOSE);
     410  dvo_catalog_unlock (&catalog1);
     411  dvo_catalog_free (&catalog1);
    312412  free (slowMoving);
    313413  free (groupA);
Note: See TracChangeset for help on using the changeset viewer.