IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38536


Ignore:
Timestamp:
Jun 24, 2015, 2:11:27 PM (11 years ago)
Author:
eugene
Message:

basic working version of addstar multi-file

Location:
branches/eam_branches/ipp-20150616/Ohana/src/addstar
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/LoadData.c

    r38533 r38536  
    1313  dvo_catalog_init (catalog, TRUE);
    1414  ALLOCATE (catalog->average, Average, 1);
    15   ALLOCATE (catalog->secfilt, SecFilt, 1);
    1615  ALLOCATE (catalog->measure, Measure, 1);
    1716  ALLOCATE (catalog->lensing, Lensing, 1);
    … …  
    117116   
    118117    AddstarClientOptions matchOptions = *options;
    119     matchOptions.radius = 0.1; // tight radius at this stage
     118    matchOptions.radius = 0.4; // tight radius at this stage
    120119    matchOptions.calibrate = FALSE;
    121120    matchOptions.only_match = FALSE;
    122121    matchOptions.nosort = FALSE;
     122    matchOptions.photcode = 0; // use an invalid photcode to avoid touching secfilt
    123123
    124124    find_matches_closest (&region, newcat, catalog, matchOptions);

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/LoadStars.c

    r38525 r38536  
    1212
    1313  Catalog *catalog = NULL;
    14   // ALLOCATE (catalog, Catalog, 1);
    15   // dvo_catalog_init (catalog, TRUE);
     14  ALLOCATE (catalog, Catalog, 1);
     15  dvo_catalog_init (catalog, TRUE);
     16  ALLOCATE (catalog->average, Average, 1);
     17  ALLOCATE (catalog->measure, Measure, 1);
     18  ALLOCATE (catalog->lensing, Lensing, 1);
    1619
    1720  for (i = 0; i < Nfile; i++) {
    … …  
    6972    }
    7073
     74    off_t NimagesStart = *Nimages;
    7175    Catalog *newcat = LoadData (f, &file[i], images, Nimages, headers, extsize, headerSets, NheaderSets, options);
    72     catalog = newcat;
    7376
    74     // XXX I need to use find_matches_closest here, but for testing:
    75     // find_matches_closest (region, catalog, newcat, options);
    76     // dvo_catalog_free (newcat);
    77     // free (newcat);
     77    // if we added data from an image, then we can merge it in
     78    if (*Nimages > NimagesStart) {
     79      double Rref, Dref;
     80      XY_to_RD (&Rref, &Dref, 0.0, 0.0, &images[0][*Nimages-1].coords);
     81      SkyRegion region;
     82      region.Rmin = Rref - 180.0;
     83      region.Rmax = Rref + 180.0;
     84      region.Dmin = Dref -  90.0;
     85      region.Dmax = Dref +  90.0;
     86     
     87      AddstarClientOptions matchOptions = *options;
     88      matchOptions.radius = 0.4; // tight radius at this stage
     89      matchOptions.calibrate = FALSE;
     90      matchOptions.only_match = FALSE;
     91      matchOptions.nosort = FALSE;
     92      matchOptions.photcode = 0; // use an invalid photcode to avoid touching secfilt
     93     
     94      find_matches_closest (&region, newcat, catalog, matchOptions);
     95    }
     96    dvo_catalog_free (newcat);
     97    free (newcat);
    7898
    7999  next_file:
    … …  
    95115    exit (0);
    96116  }
    97   ohana_memcheck (TRUE);
    98117
    99118  AddstarFileFree (file, Nfile);

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/addstar.c

    r38525 r38536  
    196196  SkyListFree (skylist);
    197197
    198   if (options.mode != ADDSTAR_MODE_REFCAT) {
    199     dvo_catalog_free (newcat);
    200     free (newcat);
    201   }
    202 
    203198  // We only measure a single value for the entire mosaic (add all images to this function)
    204199  if (options.calibrate) { FindCalibration (&images[0]); }
    … …  
    236231  float dtime = DTIME (stopAddstar, startAddstar);
    237232  fprintf (stderr, "SUCCESS: elapsed time %9.4f sec for "OFF_T_FMT" stars (%5d matches), "OFF_T_FMT" average, "OFF_T_FMT" measure, "OFF_T_FMT" lensing\n", dtime, newcat->Naverage, Nmatch,  Naverage,  Nmeasure, Nlensing);
     233
     234  if (options.mode != ADDSTAR_MODE_REFCAT) {
     235    dvo_catalog_free (newcat);
     236    free (newcat);
     237  }
    238238
    239239  // XXX test

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/find_matches_closest.c

    r38533 r38536  
    11# include "addstar.h"
    22
    3 int find_matches_closest (SkyRegion *region, Catalog *newcat, Catalog *catalog, AddstarClientOptions options) {
    4  
    5   off_t i, j, n, N, J, Jmin, status, Nstars;
     3int find_matches_closest (SkyRegion *region, Catalog *srccat, Catalog *tgtcat, AddstarClientOptions options) {
     4
     5  off_t i, j, n, N, J, Jmin, status, Nstars, mSrc;
    66  double RADIUS, RADIUS2, Rmin;
    77  double *X1, *Y1, *X2, *Y2;
    88  double dX, dY, dR;
    9   off_t *N1, *N2, *next_meas, *next_lens;
     9  off_t *N1, *N2;
    1010  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS;
    1111  int Nsecfilt, Nsec;
    … …  
    1313  Coords tcoords;
    1414
    15   // XXX this probably is no longer needed (absorbed by catalog->average[j].Nmeasure
    16   if (NSTAR_GROUP <= 0) {
    17       fprintf (stderr, "ERROR: NSTAR_GROUP NOT SET!\n");
    18       exit (1);
    19   }
    20 
    2115  /* photcode data - must by of type DEP; options.photcode is equiv photcode for all input
    2216     images this function requires incoming stars to have the same photcode.equiv value.  if
    … …  
    2721  Nsec     = GetPhotcodeNsec (options.photcode);
    2822
    29   /** allocate local arrays (newcat) **/
    30   ALLOCATE (X1, double, newcat->Naverage);
    31   ALLOCATE (Y1, double, newcat->Naverage);
    32   ALLOCATE (N1, off_t,  newcat->Naverage);
    33 
    34   if (!newcat->found_t) {
    35     ALLOCATE (newcat->found_t, off_t, newcat->Naverage);
    36     for (i = 0; i < newcat->Naverage; i++) {
    37       newcat->found_t[i] = -1;
    38     }
    39   }
    40 
    41   /** allocate local arrays (catalog) **/
    42   NAVE = Nave = catalog[0].Naverage;
     23  /** allocate local arrays (srccat) **/
     24  ALLOCATE (X1, double, srccat->Naverage);
     25  ALLOCATE (Y1, double, srccat->Naverage);
     26  ALLOCATE (N1, off_t,  srccat->Naverage);
     27
     28  if (!srccat->found_t) {
     29    ALLOCATE (srccat->found_t, off_t, srccat->Naverage);
     30    for (i = 0; i < srccat->Naverage; i++) {
     31      srccat->found_t[i] = -1;
     32    }
     33  }
     34
     35  /** allocate local arrays (tgtcat) **/
     36  NAVE = Nave = tgtcat[0].Naverage;
    4337  ALLOCATE (X2, double, NAVE);
    4438  ALLOCATE (Y2, double, NAVE);
    4539  ALLOCATE (N2, off_t,  NAVE);
    46   ALLOCATE (catalog[0].found_t, off_t, NAVE);
     40
     41  if (!tgtcat->found_t) {
     42    ALLOCATE (tgtcat[0].found_t, off_t, NAVE);
     43    for (i = 0; i < NAVE; i++) {
     44      tgtcat->found_t[i] = -1;
     45    }
     46  }
    4747  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
    4848
    4949  /* internal counters */
    5050  Nmatch = 0;
    51   NMEAS = Nmeas = catalog[0].Nmeasure;
    52   NLENS = Nlens = catalog[0].Nlensing;
    53 
    54   // current max obj ID for this catalog
    55   objID = catalog[0].objID;
    56   catID = catalog[0].catID;
     51  NMEAS = Nmeas = tgtcat[0].Nmeasure;
     52  NLENS = Nlens = tgtcat[0].Nlensing;
     53
     54  off_t *next_meas = NULL;
     55  off_t *next_lens = NULL;
     56
     57  // current max obj ID for this tgtcat
     58  objID = tgtcat[0].objID;
     59  catID = tgtcat[0].catID;
    5760
    5861  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has
    5962   * the advantage that every point in R,D has a mapping to a unique X,Y.  However, note
    6063   * that not all possible X,Y points map back to R,D and the local plate scale changes
    61    * far from the projection pole. We use the center of the region (catalog) for crval1,2.
     64   * far from the projection pole. We use the center of the region (tgtcat) for crval1,2.
    6265   */
    6366  InitCoords (&tcoords, "DEC--ARC");
    … …  
    7275  /* build spatial index (RA sort) referencing input array sequence */
    7376  Nstars = 0;
    74   for (i = 0; i < newcat->Naverage; i++) {
    75     status = RD_to_XY (&X1[Nstars], &Y1[Nstars], newcat->average[i].R, newcat->average[i].D, &tcoords);
     77  for (i = 0; i < srccat->Naverage; i++) {
     78    status = RD_to_XY (&X1[Nstars], &Y1[Nstars], srccat->average[i].R, srccat->average[i].D, &tcoords);
    7679    if (!status) continue;
    7780    N1[Nstars] = i;
    … …  
    7982  }
    8083  if (Nstars < 1) {
    81     if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
     84    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", tgtcat[0].filename);
    8285    free (X1);
    8386    free (Y1);
    … …  
    8790    free (N2);
    8891    return (0);
     92    // XXX check on this
    8993  }
    9094  if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars);
    … …  
    9296  /* build spatial index (RA sort) */
    9397  for (i = 0; i < Nave; i++) {
    94     RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
     98    RD_to_XY (&X2[i], &Y2[i], tgtcat[0].average[i].R, tgtcat[0].average[i].D, &tcoords);
    9599    N2[i] = i;
    96     catalog[0].found_t[N2[i]] = -1;
     100    tgtcat[0].found_t[N2[i]] = -1;
    97101  }
    98102  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    99103
    100104  /* set up pointers for linked list of measure */
    101   if (catalog[0].sorted && (catalog[0].Nmeasure == catalog[0].Nmeasure_disk)) {
    102     // this version is only valid if we have done a full catalog load, and if the catalog
     105  if (tgtcat[0].sorted && (tgtcat[0].Nmeasure == tgtcat[0].Nmeasure_disk)) {
     106    // this version is only valid if we have done a full tgtcat load, and if the tgtcat
    103107    // is sorted while processed
    104     next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
    105     next_lens = init_lensing_links (catalog[0].average, Nave, catalog[0].lensing, Nlens);
     108    next_meas = init_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
     109    next_lens = init_lensing_links (tgtcat[0].average, Nave, tgtcat[0].lensing, Nlens);
    106110  } else {
    107     next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
    108     next_lens = build_lensing_links (catalog[0].average, Nave, catalog[0].lensing, Nlens);
     111    next_meas = build_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
     112    next_lens = build_lensing_links (tgtcat[0].average, Nave, tgtcat[0].lensing, Nlens);
    109113  }   
    110114
    … …  
    113117  RADIUS2 = RADIUS*RADIUS;
    114118
    115   /** find matched stars **/
    116   // XXX could use NSTAR_GROUP to do this match more quicky if NSTAR_GROUP > 1
     119  /****************** find matched stars ********************/
     120
    117121  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    118122    if (!finite(X1[i]) || !finite(Y1[i])) {
    … …  
    138142
    139143    // XXX check that this is allocated
    140     if (newcat->found_t[N1[i]] != -1) {
    141         /* this star has already been assigned to an object in this or another catalog */
     144    if (srccat->found_t[N1[i]] != -1) {
     145        /* this star has already been assigned to an object in this or another tgtcat */
    142146        i++;
    143147        continue;
    … …  
    178182
    179183    /* make sure there is space for next entry */
    180     if (Nmeas >= NMEAS) {
    181       NMEAS = Nmeas + 1000;
     184    if (Nmeas >= NMEAS - srccat->average[N].Nmeasure) {
     185      NMEAS = Nmeas + 1000 + srccat->average[N].Nmeasure;
    182186      REALLOCATE (next_meas, off_t, NMEAS);
    183       REALLOCATE (catalog[0].measure, Measure, NMEAS);
    184     }
    185     if (Nlens >= NLENS) {
    186       NLENS = Nlens + 1000;
     187      REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
     188    }
     189    if (Nlens >= NLENS - srccat->average[N].Nlensing) {
     190      NLENS = Nlens + 1000 + srccat->average[N].Nlensing;
    187191      REALLOCATE (next_lens, off_t, NLENS);
    188       REALLOCATE (catalog[0].lensing, Lensing, NLENS);
    189     }
    190 
    191     /* add to end of measurement list */
    192     add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
    193        
     192      REALLOCATE (tgtcat[0].lensing, Lensing, NLENS);
     193    }
     194
    194195    /** add measurements for this star **/
    195 
    196     // set the new measurements
    197     catalog[0].measure[Nmeas]          = newcat->measure[N];
    198 
    199     // measure now carries R,D (not dR,dD)
    200     catalog[0].measure[Nmeas].dbFlags  = 0;
    201     catalog[0].measure[Nmeas].averef   = n;
    202     catalog[0].measure[Nmeas].objID    = catalog[0].average[n].objID;
    203     catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    204 
    205     float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]);
    206 
    207     // rationalize dR
    208     if (dRoff > +180.0*3600.0) {
    209       // average on high end of boundary, move star up
    210       catalog[0].measure[Nmeas].R += 360.0;
    211       dRoff -= 360.0*3600.0;
    212     }
    213     if (dRoff < -180.0*3600.0) {
    214       // average on low end of boundary, move star down
    215       catalog[0].measure[Nmeas].R -= 360.0;
    216       dRoff += 360.0*3600.0;
    217     }
    218     if (fabs(dRoff) > 10*RADIUS) {
     196    off_t mOff = srccat->average[N].measureOffset;
     197    for (mSrc = 0; mSrc < srccat->average[N].Nmeasure; mSrc++) {
     198      /* add to end of measurement list */
     199      add_meas_link (&tgtcat[0].average[n], next_meas, Nmeas, NMEAS);
     200
     201      // set the new measurements
     202      tgtcat[0].measure[Nmeas]          = srccat->measure[mSrc + mOff];
     203
     204      // measure now carries R,D (not dR,dD)
     205      tgtcat[0].measure[Nmeas].dbFlags  = 0;
     206      tgtcat[0].measure[Nmeas].averef   = n;
     207      tgtcat[0].measure[Nmeas].objID    = tgtcat[0].average[n].objID;
     208      tgtcat[0].measure[Nmeas].catID    = tgtcat[0].catID;
     209
     210      float dRoff = dvoOffsetR(&tgtcat[0].measure[Nmeas], &tgtcat[0].average[n]);
     211
     212      // rationalize dR
     213      if (dRoff > +180.0*3600.0) {
     214        // average on high end of boundary, move star up
     215        tgtcat[0].measure[Nmeas].R += 360.0;
     216        dRoff -= 360.0*3600.0;
     217      }
     218      if (dRoff < -180.0*3600.0) {
     219        // average on low end of boundary, move star down
     220        tgtcat[0].measure[Nmeas].R -= 360.0;
     221        dRoff += 360.0*3600.0;
     222      }
     223      if (fabs(dRoff) > 10*RADIUS) {
    219224        // take declination into account and check again.
    220         double cosD = cos(RAD_DEG*catalog[0].average[n].D);
     225        double cosD = cos(RAD_DEG*tgtcat[0].average[n].D);
    221226        if (fabs(dRoff*cosD) > 10*RADIUS) {
    222             fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n",
    223                      catalog[0].average[n].R, catalog[0].average[n].D,
    224                      newcat->average[N].R, newcat->average[N].D,
    225                      X1[i], X2[Jmin],
    226                      Y1[i], Y2[Jmin]);
    227             // XXX abort on this? -- this is a bad failure...
     227          fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n",
     228                   tgtcat[0].average[n].R, tgtcat[0].average[n].D,
     229                   srccat->average[N].R, srccat->average[N].D,
     230                   X1[i], X2[Jmin],
     231                   Y1[i], Y2[Jmin]);
     232          // XXX abort on this? -- this is a bad failure...
    228233        }
    229     }
    230 
    231     // add the lensing values if they exist
    232     if (newcat->lensing) {
    233       add_lens_link (&catalog[0].average[n], next_lens, Nlens, NLENS); // ?
    234       catalog[0].lensing[Nlens] = newcat->lensing[N];
    235      
    236       catalog[0].lensing[Nlens].averef = n;
    237       catalog[0].lensing[Nlens].objID = catalog[0].average[n].objID;
    238       catalog[0].lensing[Nlens].catID = catalog[0].catID;
    239       catalog[0].average[n].Nlensing ++;
    240       Nlens ++;
    241     }
    242 
    243     /* adds the measurement to the calibration if appropriate color terms are found */
    244     /* we call this before (optionally) setting the average magnitude to avoid auto-correlations */
    245     if (options.calibrate) {
    246       AddToCalibration (&catalog[0].average[n], &catalog[0].secfilt[n*Nsecfilt], catalog[0].measure, &catalog[0].measure[Nmeas], next_meas, N);
    247     }
    248 
    249     /* set the average magnitude if not already set and the photcode.equiv is not 0 */
    250     /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    251     if (Nsec > -1) {
    252       if (isnan(catalog[0].secfilt[n*Nsecfilt+Nsec].M)) {
    253         catalog[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF);
    254       }
    255     }
    256 
    257     /* if we choose to flag close encounters, see find_matches.c */
    258     /* if we choose to calculate RA,DEC averages, see update_coords.c */
    259 
    260     /* Nm is updated, but not written out in -update mode (for existing entries)
    261        Nm is recalculated in build_meas_links if loaded table is not sorted */
    262     newcat->found_t[N] = Nmeas;
    263     catalog[0].found_t[n] = Nmeas;
    264     catalog[0].average[n].Nmeasure ++;
    265     Nmeas ++;
    266     i++;
    267   }
    268 
    269   /** incorporate unmatched image stars, if this star is in field of this catalog **/
    270   /* these new entries are all written out in UPDATE mode */
    271   for (i = 0; (i < Nstars) && !options.only_match; i += NSTAR_GROUP) {
    272     fprintf (stderr, "loop %d\n", (int) i);
    273     /* make sure there is space for next entry */
    274     if (Nmeas >= NMEAS - NSTAR_GROUP) {
    275       NMEAS = Nmeas + 1000;
    276       REALLOCATE (next_meas, off_t, NMEAS);
    277       REALLOCATE (catalog[0].measure, Measure, NMEAS);
    278     }
    279     if (Nlens >= NLENS - NSTAR_GROUP) {
    280       NLENS = Nlens + 1000;
    281       REALLOCATE (next_lens, off_t, NLENS);
    282       REALLOCATE (catalog[0].lensing, Lensing, NLENS);
    283     }
    284     if (Nave >= NAVE) {
    285       NAVE = Nave + 1000;
    286       REALLOCATE (catalog[0].average, Average, NAVE);
    287       REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
    288     }
    289 
    290     if (newcat->found_t[i] != -1) continue;
    291     if (!IN_REGION (newcat->average[i].R, newcat->average[i].D)) continue;
    292 
    293     dvo_average_init (&catalog[0].average[Nave]);
    294     catalog[0].average[Nave].R             = newcat->average[i].R;
    295     catalog[0].average[Nave].D             = newcat->average[i].D;
    296 
    297     catalog[0].average[Nave].Nmeasure      = NSTAR_GROUP;
    298     catalog[0].average[Nave].measureOffset = Nmeas;
    299     catalog[0].average[Nave].objID         = objID;
    300     catalog[0].average[Nave].catID         = catID;
    301 
    302     ohana_memcheck(1);
    303 
    304     if (PSPS_ID) {
    305         catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D);
    306     }
    307 
    308     objID ++;
    309 
    310     for (j = 0; j < Nsecfilt; j++) {
    311       // XXXX inconsistent with my alloc above
    312       dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
    313     }
    314 
    315     ohana_memcheck(1);
    316 
    317     for (j = 0; j < NSTAR_GROUP; j++) {
    318       // supply the measurments from this detection
    319       dvo_measure_init (&catalog[0].measure[Nmeas]);
    320       catalog[0].measure[Nmeas] = newcat->measure[i + j];
    321 
    322       // the following measure elements cannot be set until here:
    323       catalog[0].measure[Nmeas].dbFlags  = 0;
    324       catalog[0].measure[Nmeas].averef   = Nave;
    325       catalog[0].measure[Nmeas].objID    = catalog[0].average[Nave].objID;
    326       catalog[0].measure[Nmeas].catID    = catalog[0].catID;
     234      }
     235
     236      /* adds the measurement to the calibration if appropriate color terms are found */
     237      /* we call this before (optionally) setting the average magnitude to avoid auto-correlations */
     238      if (options.calibrate) {
     239        AddToCalibration (&tgtcat[0].average[n], &tgtcat[0].secfilt[n*Nsecfilt], tgtcat[0].measure, &tgtcat[0].measure[Nmeas], next_meas, N);
     240      }
    327241
    328242      /* set the average magnitude if not already set and the photcode.equiv is not 0 */
    329243      /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    330244      if (Nsec > -1) {
    331         catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF);
    332       }
    333 
    334       /* next[Nmeas] should always be -1 in this context (it is always the only
    335          measurement for the star) */
    336       newcat->found_t[i+j] = Nmeas;
    337       next_meas[Nmeas] = -1;  // inital value here update below
     245        if (isnan(tgtcat[0].secfilt[n*Nsecfilt+Nsec].M)) {
     246          tgtcat[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     247        }
     248      }
     249
     250      /* Nm is updated, but not written out in -update mode (for existing entries)
     251         Nm is recalculated in build_meas_links if loaded table is not sorted */
     252      tgtcat[0].average[n].Nmeasure ++;
    338253      Nmeas ++;
    339254    }
    340     for (j = 0; j < NSTAR_GROUP - 1; j++) {
    341       next_meas[Nmeas - NSTAR_GROUP + j] = Nmeas - NSTAR_GROUP + j + 1;
    342     }
    343 
    344     ohana_memcheck(1);
     255
     256    // add the lensing values if they exist
     257    if (srccat->lensing) {
     258      /** add measurements for this star **/
     259      off_t mOff = srccat->average[N].lensingOffset;
     260      for (mSrc = 0; mSrc < srccat->average[N].Nlensing; mSrc++) {
     261        /* add to end of lensing list */
     262
     263        add_lens_link (&tgtcat[0].average[n], next_lens, Nlens, NLENS); // ?
     264        tgtcat[0].lensing[Nlens] = srccat->lensing[mSrc + mOff];
     265     
     266        tgtcat[0].lensing[Nlens].averef = n;
     267        tgtcat[0].lensing[Nlens].objID = tgtcat[0].average[n].objID;
     268        tgtcat[0].lensing[Nlens].catID = tgtcat[0].catID;
     269        tgtcat[0].average[n].Nlensing ++;
     270        Nlens ++;
     271      }
     272    }
     273
     274    /* if we choose to flag close encounters, see find_matches.c */
     275    /* if we choose to calculate RA,DEC averages, see update_coords.c */
     276
     277    srccat->found_t[N] = n;
     278    tgtcat[0].found_t[n] = -1;
     279    i++;
     280  }
     281
     282  /*************** add unmatched stars *************************/
     283
     284  /** incorporate unmatched image stars, if this star is in field of this tgtcat **/
     285  /* these new entries are all written out in UPDATE mode */
     286  for (i = 0; (i < Nstars) && !options.only_match; i++) {
     287    /* make sure there is space for next entry */
     288    if (Nmeas >= NMEAS - srccat->average[i].Nmeasure) {
     289      NMEAS = Nmeas + 1000 + srccat->average[i].Nmeasure;
     290      REALLOCATE (next_meas, off_t, NMEAS);
     291      REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
     292    }
     293    if (Nlens >= NLENS - srccat->average[i].Nlensing) {
     294      NLENS = Nlens + 1000 + srccat->average[i].Nlensing;
     295      REALLOCATE (next_lens, off_t, NLENS);
     296      REALLOCATE (tgtcat[0].lensing, Lensing, NLENS);
     297    }
     298    if (Nave >= NAVE) {
     299      NAVE = Nave + 1000;
     300      REALLOCATE (tgtcat[0].average, Average, NAVE);
     301      REALLOCATE (tgtcat[0].found_t, off_t, NAVE);
     302      if (tgtcat[0].secfilt) {
     303        // we only update the secfilt table if it has been allocated for output
     304        REALLOCATE (tgtcat[0].secfilt, SecFilt, NAVE*tgtcat[0].Nsecfilt);
     305      }
     306    }
     307
     308    if (srccat->found_t[i] != -1) continue;
     309    if (!IN_REGION (srccat->average[i].R, srccat->average[i].D)) continue;
     310
     311    dvo_average_init (&tgtcat[0].average[Nave]);
     312    tgtcat[0].average[Nave].R              = srccat->average[i].R;
     313    tgtcat[0].average[Nave].D              = srccat->average[i].D;
     314
     315    tgtcat[0].average[Nave].Nmeasure       = srccat->average[i].Nmeasure;
     316    tgtcat[0].average[Nave].measureOffset  = Nmeas;
     317    tgtcat[0].average[Nave].objID          = objID;
     318    tgtcat[0].average[Nave].catID          = catID;
     319
     320    // this star has not yet been merged into another tgtcat
     321    tgtcat[0].found_t[Nave]                = -1;
     322
     323    if (PSPS_ID) {
     324        tgtcat[0].average[Nave].extID = CreatePSPSObjectID(tgtcat[0].average[Nave].R, tgtcat[0].average[Nave].D);
     325    }
     326
     327    objID ++;
     328
     329    // we only update the secfilt table if it has been allocated for output
     330    for (j = 0; tgtcat[0].secfilt && (j < Nsecfilt); j++) {
     331      dvo_secfilt_init (&tgtcat[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
     332    }
     333
     334    off_t mOff = srccat->average[i].measureOffset;
     335    for (mSrc = 0; mSrc < srccat->average[i].Nmeasure; mSrc++) {
     336      // supply the measurements from this detection
     337      dvo_measure_init (&tgtcat[0].measure[Nmeas]);
     338      tgtcat[0].measure[Nmeas] = srccat->measure[mSrc + mOff];
     339
     340      // the following measure elements cannot be set until here:
     341      tgtcat[0].measure[Nmeas].dbFlags  = 0;
     342      tgtcat[0].measure[Nmeas].averef   = Nave;
     343      tgtcat[0].measure[Nmeas].objID    = tgtcat[0].average[Nave].objID;
     344      tgtcat[0].measure[Nmeas].catID    = tgtcat[0].catID;
     345
     346      /* set the average magnitude if not already set and the photcode.equiv is not 0 */
     347      /* in UPDATE mode, this value is not saved; use relphot to recalculate */
     348      if (Nsec > -1) {
     349        tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
     350      }
     351
     352      Nmeas ++;
     353
     354      // update the next_meas pointer for this entry (last one for this star is -1)
     355      next_meas[Nmeas-1] = (mSrc < srccat->average[i].Nmeasure - 1) ? Nmeas : -1;
     356    }
     357    srccat->found_t[i] = Nave;
    345358
    346359    // if we have lensing data, insert that as well
    347     if (newcat->lensing) {
    348       catalog[0].average[Nave].Nlensing = NSTAR_GROUP;
    349       catalog[0].average[Nave].lensingOffset = Nlens;
    350       for (j = 0; j < NSTAR_GROUP; j++) {
    351         // add the lensing values if they exist
    352         if (newcat->lensing) {
    353           dvo_lensing_init (&catalog[0].lensing[Nlens]);
    354           catalog[0].lensing[Nlens] = newcat->lensing[i + j];
     360    if (srccat->lensing) {
     361      tgtcat[0].average[Nave].Nlensing = srccat->average[i].Nlensing;
     362      tgtcat[0].average[Nave].lensingOffset = Nlens;
     363
     364      off_t mOff = srccat->average[i].lensingOffset;
     365      for (mSrc = 0; mSrc < srccat->average[i].Nlensing; mSrc++) {
     366        dvo_lensing_init (&tgtcat[0].lensing[Nlens]);
     367        tgtcat[0].lensing[Nlens] = srccat->lensing[mSrc + mOff];
    355368         
    356           catalog[0].lensing[Nlens].averef = Nave;
    357           catalog[0].lensing[Nlens].objID = catalog[0].average[Nave].objID;
    358           catalog[0].lensing[Nlens].catID = catalog[0].catID;
    359           next_lens[Nlens] = -1;
    360           Nlens ++;
    361         }
    362       }
    363       for (j = 0; j < NSTAR_GROUP - 1; j++) {
    364         next_lens[Nlens - NSTAR_GROUP + j] = Nlens - NSTAR_GROUP + j + 1;
    365       }
    366     }
    367     ohana_memcheck(1);
    368 
     369        tgtcat[0].lensing[Nlens].averef = Nave;
     370        tgtcat[0].lensing[Nlens].objID = tgtcat[0].average[Nave].objID;
     371        tgtcat[0].lensing[Nlens].catID = tgtcat[0].catID;
     372        Nlens ++;
     373        next_lens[Nlens-1] = (mSrc < srccat->average[i].Nlensing - 1) ? Nlens : -1;
     374      }
     375    }
    369376    Nave ++;
    370377  }
    371      
    372   REALLOCATE (catalog[0].average, Average, Nave);
    373   REALLOCATE (catalog[0].measure, Measure, Nmeas);
    374   REALLOCATE (catalog[0].lensing, Lensing, Nlens);
     378
     379  REALLOCATE (tgtcat[0].average, Average, Nave);
     380  REALLOCATE (tgtcat[0].measure, Measure, Nmeas);
     381  REALLOCATE (tgtcat[0].lensing, Lensing, Nlens);
    375382 
    376383  if (options.nosort) {
    377     catalog[0].sorted = FALSE;
     384    tgtcat[0].sorted = FALSE;
    378385  } else {
    379     catalog[0].sorted = TRUE;
    380     catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
    381     catalog[0].lensing = sort_lensing (catalog[0].average, Nave, catalog[0].lensing, Nlens, next_lens);
    382   }
    383 
    384   /* note stars which have been found in this catalog */
    385   for (i = 0; i < newcat->Naverage; i++) {
    386     if (newcat->found_t[i] > -1) {
    387       newcat->found_t[i] = -2;
    388     }
    389   }
    390 
    391   /* check if the catalog has changed?  if no change, no need to write */
    392   catalog[0].objID    = objID; // new max value, save on catalog close
    393   catalog[0].Naverage = Nave;
    394   catalog[0].Nmeasure = Nmeas;
    395   catalog[0].Nlensing = Nlens;
    396   catalog[0].Nsecfilt_mem = Nave*Nsecfilt;
     386    tgtcat[0].sorted = TRUE;
     387    tgtcat[0].measure = sort_measure (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas, next_meas);
     388    tgtcat[0].lensing = sort_lensing (tgtcat[0].average, Nave, tgtcat[0].lensing, Nlens, next_lens);
     389  }
     390
     391  /* note stars which have been found in this tgtcat */
     392  for (i = 0; i < srccat->Naverage; i++) {
     393    if (srccat->found_t[i] > -1) {
     394      srccat->found_t[i] = -2;
     395    }
     396  }
     397
     398  /* check if the tgtcat has changed?  if no change, no need to write */
     399  tgtcat[0].objID    = objID; // new max value, save on tgtcat close
     400  tgtcat[0].Naverage = Nave;
     401  tgtcat[0].Nmeasure = Nmeas;
     402  tgtcat[0].Nlensing = Nlens;
     403  tgtcat[0].Nsecfilt_mem = tgtcat[0].secfilt ? Nave*Nsecfilt : 0;
    397404  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nlens, Nmatch);
    398405
    … …  
    412419   notes:
    413420   
    414    for finding if a catalog star is in an image or an image star is in the catalog:
     421   for finding if a tgtcat star is in an image or an image star is in the tgtcat:
    415422   
    416    catalogs have boundaries defined by RA and DEC, but they may curve in projection
     423   tgtcats have boundaries defined by RA and DEC, but they may curve in projection
    417424   images have boundaries which are lines in pixels coords, but curve in RA and DEC
    418425   
    419    catalog[0].found[Ncat] but stars[Nstar].found
     426   tgtcat[0].found[Ncat] but stars[Nstar].found
    420427   
    421428*/

  • branches/eam_branches/ipp-20150616/Ohana/src/addstar/test/simple.dvo

    r38524 r38536  
    88
    99macro test.one
    10   test.fields PS1_V5    PS1_V4
     10  test.fields PS1_V5 PS1_V5
    1111end
    1212
Note: See TracChangeset for help on using the changeset viewer.