IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38593


Ignore:
Timestamp:
Jul 15, 2015, 8:03:07 PM (11 years ago)
Author:
eugene
Message:

updates to replace Stars with Catalog

Location:
branches/eam_branches/ipp-20150625/Ohana/src/addstar
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/Makefile

    r38553 r38593  
    4040# programs in 'SERVER' use the client-server concept and are out of date
    4141
    42 INSTALL = addstar sedstar loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell
    43 OTHERS  = load2mass loadwise loadtycho loadsupercos
     42INSTALL = addstar sedstar loadstarpar loadstarpar_client loadICRF loadICRF_client load2mass skycells mkcmf dumpskycells findskycell
     43OTHERS  = loadwise loadtycho loadsupercos
    4444SERVER  = addstar_client addstarc addstard addstart
    4545
    … …  
    5151# all: addstar sedstar load2mass skycells
    5252
    53 # $(SRC)/find_matches.$(ARCH).o \
    54 # $(SRC)/find_matches_closest.$(ARCH).o \
    55 # $(SRC)/find_matches_refstars.$(ARCH).o \
    56 # $(SRC)/find_matches_closest_refstars.$(ARCH).o \
    5753# $(SRC)/find_subset.$(ARCH).o \
    58 # $(SRC)/getgsc.$(ARCH).o \
    59 # $(SRC)/getusno.$(ARCH).o \
    60 # $(SRC)/getusnob.$(ARCH).o \
    61 # $(SRC)/greference.$(ARCH).o \
    62 # $(SRC)/grefstars.$(ARCH).o \
    6354# $(SRC)/dump.$(ARCH).o \
    6455# $(SRC)/ReadStarsTEXT.$(ARCH).o \
    6556# $(SRC)/StarOps.$(ARCH).o \
    66 # $(SRC)/replace_match.$(ARCH).o \
    67 #
    6857
    6958ADDSTAR = \
    … …  
    10493$(SRC)/ReadStarsUKIRT.$(ARCH).o \
    10594$(SRC)/ReadXradFITS.$(ARCH).o \
     95$(SRC)/find_matches.$(ARCH).o \
     96$(SRC)/find_matches_closest.$(ARCH).o \
     97$(SRC)/find_matches_refstars.$(ARCH).o \
     98$(SRC)/find_matches_closest_refstars.$(ARCH).o \
     99$(SRC)/getgsc.$(ARCH).o \
     100$(SRC)/getusno.$(ARCH).o \
     101$(SRC)/getusnob.$(ARCH).o \
     102$(SRC)/greference.$(ARCH).o \
     103$(SRC)/grefstars.$(ARCH).o \
     104$(SRC)/replace_match.$(ARCH).o \
    106105$(SRC)/FilterStars.$(ARCH).o \
    107106$(SRC)/ImageOptions.$(ARCH).o \
    … …  
    243242$(SRC)/load2mass.$(ARCH).o \
    244243$(SRC)/load2mass_as_rawdata.$(ARCH).o \
    245 $(SRC)/load2mass_catalog.$(ARCH).o \
    246244$(SRC)/get2mass_ops.$(ARCH).o \
    247245$(SRC)/get2mass_full.$(ARCH).o \
    … …  
    252250$(SRC)/update_coords.$(ARCH).o \
    253251$(SRC)/SkyRegionUtils.$(ARCH).o \
    254 $(SRC)/StarOps.$(ARCH).o \
    255 $(SRC)/ConfigInit.$(ARCH).o \
    256 $(SRC)/Shutdown.$(ARCH).o \
    257 $(SRC)/SetSignals.$(ARCH).o \
    258 $(SRC)/psps_ids.$(ARCH).o
     252$(SRC)/ConfigInit.$(ARCH).o \
     253$(SRC)/Shutdown.$(ARCH).o \
     254$(SRC)/SetSignals.$(ARCH).o \
     255$(SRC)/psps_ids.$(ARCH).o
     256
     257# $(SRC)/load2mass_catalog.$(ARCH).o
    259258
    260259LOAD-GALPHOT = \

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/2mass.h

    r24977 r38593  
    2222
    2323SkyTable *get2mass_acc (SkyRegion *patch, char *path, char *accel);
    24 Stars    *get2mass_2DR_data (SkyRegion *region, char *filename, SkyRegion *patch, int photcode, int *nstars);
    25 Stars    *get2mass_AS_data (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars);
    26 Stars    *get2mass_AS_rawdata (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars);
     24
     25// Stars    *get2mass_2DR_data (SkyRegion *region, char *filename, SkyRegion *patch, int photcode, int *nstars);
     26// Stars    *get2mass_AS_data (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars);
     27// Stars    *get2mass_AS_rawdata (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars);
    2728
    2829SkyTable *scan2mass_acc (char *path, char *accel);
    … …  
    3536int       load2mass_as_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options);
    3637SkyTable *load2mass_acc (char *path, char *accel);
    37 int       get2mass_3star (Stars *star, char *line, int Nmax);
    38 int       load2mass_catalog (Catalog *catalog, Stars *stars, int Nstars);
     38
     39// int       get2mass_3star (Stars *star, char *line, int Nmax);
     40// int       load2mass_catalog (Catalog *catalog, Stars *stars, int Nstars);
    3941
    4042int       get2mass_setup (int photcode);
    4143int       get2mass_coords (char *line, double *R, double *D, int Nmax);
    42 int       get2mass_star (Stars *star, char *line, int Nmax);
    43 int       get2mass_3star (Stars *star, char *line, int Nmax);
    4444
    45 int get2mass_3star_full (Stars **star, char *line, int Nmax);
     45// int       get2mass_star (Stars *star, char *line, int Nmax);
     46// int       get2mass_3star (Stars *star, char *line, int Nmax);
     47
     48int get2mass_3star_full (Measure *measure, char *line, int Nmax);
    4649char *next2MASSfield (char *line);
    47 int set2MASS_ph_qual (Stars *star, char qual);
    48 int set2MASS_rd_flag (Stars *star, char qual);
    49 int set2MASS_cc_flag (Stars *star, char qual);
    50 int set2MASS_bl_flag (Stars *star, char qual);
    51 int set2MASS_gal_flag (Stars *star, char qual);
    52 int set2MASS_mp_flag (Stars *star, char qual);
    53 int set2MASS_dup_flag (Stars *star, char qual);
    54 int set2MASS_use_flag (Stars *star, char qual);
     50int set2MASS_ph_qual (Measure *measure, char qual);
     51int set2MASS_rd_flag (Measure *measure, char qual);
     52int set2MASS_cc_flag (Measure *measure, char qual);
     53int set2MASS_bl_flag (Measure *measure, char qual);
     54int set2MASS_gal_flag (Measure *measure, char qual);
     55int set2MASS_mp_flag (Measure *measure, char qual);
     56int set2MASS_dup_flag (Measure *measure, char qual);
     57int set2MASS_use_flag (Measure *measure, char qual);
    5558int get2mass_sortStars (TMStars *tstars, int Ntstars);

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/addstar.h

    r38553 r38593  
    4949  int extnum_xrad;
    5050} HeaderSet;
    51 
    52 # if (0)
    53 typedef struct {
    54   Average  average;
    55   Measure  measure;
    56   Lensing *lensing; // optionally carry out the lensing measurements
    57   int found;
    58 } Stars;
    59 
    60 typedef struct {
    61   unsigned int Nstars;
    62   Stars *stars;
    63   off_t Nimages;
    64   Image *images;
    65   Coords *mosaic;
    66   SkyRegion *patch;
    67   AddstarClientOptions *options;
    68   char *refcat;
    69 } DVO_DATA;
    70 # endif
    7151
    7252# define IDTYPE int
    … …  
    222202int        find_matches           PROTO((SkyRegion *region, Catalog *newcat, Catalog *catalog, AddstarClientOptions options));
    223203
    224 # if (0)
    225 int        find_matches_refstars         PROTO((SkyRegion *region, Stars **stars, unsigned int Nstars, Catalog *catalog, AddstarClientOptions options));
    226 int        find_matches_closest_refstars PROTO((SkyRegion *region, Stars **stars, unsigned int Nstars, Catalog *catalog, AddstarClientOptions options));
    227 
    228 int        dump_rawstars          PROTO((Stars *stars, unsigned int Nstars));
    229 
    230 Stars    **find_subset            PROTO((SkyRegion *region, Stars *stars, unsigned int Nstars, unsigned int *NSTARS));
     204int        find_matches_refstars         PROTO((SkyRegion *region, Catalog *srccat, Catalog *tgtcat, AddstarClientOptions options));
     205int        find_matches_closest_refstars PROTO((SkyRegion *region, Catalog *srccat, Catalog *tgtcat, AddstarClientOptions options));
     206
     207Catalog   *greference             PROTO((char *Refcat, SkyRegion *catstats, int photcode));
     208
    231209int        gcatalog               PROTO((Catalog *catalog));
    232 Stars     *get2mass               PROTO((SkyRegion *patch, int photcode, int mode, unsigned int *NSTARS));
    233 
    234 Stars     *greference             PROTO((char *Refcat, SkyRegion *catstats, int photcode, unsigned int *nstars));
    235 Stars     *getgsc                 PROTO((SkyRegion *patch, unsigned int *NSTARS));
    236 Stars     *getusno                PROTO((SkyRegion *catstats, int photcode, unsigned int *Nstars));
    237 Stars     *getusnob               PROTO((SkyRegion *catstats, int photcode, unsigned int *Nstars));
     210Catalog   *getgsc                 PROTO((SkyRegion *patch));
     211Catalog   *getusno                PROTO((SkyRegion *catstats));
     212Catalog   *getusnob               PROTO((SkyRegion *catstats));
    238213
    239214// load text-based stars (REF only in the sense of REF photcodes)
    240 Stars     *grefstars              PROTO((char *file, int photcode, unsigned int *Nstars));
    241 
    242 Stars     *rd_gsc                 PROTO((char *filename, unsigned int *nstars));
    243 
    244 int        replace_match          PROTO((Average *average, Measure *measure, Stars *star));
    245 # endif
     215Catalog   *grefstars              PROTO((char *file, int photcode));
     216
     217Catalog   *rd_gsc                 PROTO((char *filename));
     218
     219int        replace_match          PROTO((Average *average, Measure *measure, Measure *newmeas, off_t *found));
    246220
    247221Catalog   *LoadStars              PROTO((char *file, Image **images, off_t *Nimages, AddstarClientOptions *options));
    … …  
    297271
    298272SkyList   *SkyListExistingSubset  PROTO((SkyList *input, char *path));
     273SkyList   *SkyListForStars        PROTO((SkyTable *table, int depth, Catalog *catalog));
    299274
    300275// these are all for the addstar client/server which has not been maintained
    … …  
    311286int        UpdateDatabase_Reflist PROTO((AddstarClientOptions *options, Stars *stars, unsigned int Nstars));
    312287int        UpdateDatabase_Refcat  PROTO((AddstarClientOptions *options, SkyRegion *UserPatch, char *refcat));
    313 SkyList   *SkyListForStars        PROTO((SkyTable *table, int depth, Stars *stars, unsigned int Nstars));
    314288int        SkyListSetPath         PROTO((SkyList *list, char *path));
    315289int        InitDataset            PROTO((void));

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/SkyListForStars.c

    r38553 r38593  
    11# include "addstar.h"
    22
    3 # if (0)
    4 void sort_stars_ra (Stars *stars, int N) {
     3void sort_stars_ra (Catalog *catalog) {
    54
    6 # define SWAPFUNC(A,B){ Stars tmp; tmp = stars[A]; stars[A] = stars[B]; stars[B] = tmp; }
    7 # define COMPARE(A,B)(stars[A].average.R < stars[B].average.R)
     5# define SWAPFUNC(A,B){ Average tmp; tmp = catalog->average[A]; catalog->average[A] = catalog->average[B]; catalog->average[B] = tmp; }
     6# define COMPARE(A,B)(catalog->average[A].R < catalog->average[B].R)
    87
    9   OHANA_SORT (N, COMPARE, SWAPFUNC);
     8  OHANA_SORT (catalog->Naverage, COMPARE, SWAPFUNC);
    109
    1110# undef SWAPFUNC
    … …  
    1413}
    1514
    16 SkyList *SkyListForStars (SkyTable *table, int depth, Stars *stars, unsigned int Nstars) {
     15SkyList *SkyListForStars (SkyTable *table, int depth, Catalog *catalog) {
    1716 
    1817  int i, j, Nr, NR;
    … …  
    2827  list[0].ownElements = FALSE; // free these elements when freeing the list
    2928
    30   sort_stars_ra (stars, Nstars); /* sort by RA */
    31   for (i = 0; i < Nstars; i++) {
    32     if (stars[i].found == -2) continue;
    33     here = SkyRegionByPoint (table, depth, stars[i].average.R, stars[i].average.D);
    34     stars[i].found = -2;
     29  sort_stars_ra (catalog); /* sort by RA */
     30
     31  ALLOCATE (catalog->found_t, off_t, catalog->Naverage);
     32  for (i = 0; i < catalog->Naverage; i++) {
     33    catalog->found_t[i] = 0;
     34  }
     35
     36  for (i = 0; i < catalog->Naverage; i++) {
     37    if (catalog->found_t[i] == -2) continue;
     38    here = SkyRegionByPoint (table, depth, catalog->average[i].R, catalog->average[i].D);
     39    catalog->found_t[i] = -2;
    3540    /* search forward for all contained stars */
    36     for (j = i; j < Nstars; j++) {
    37       if (stars[j].average.R >= here[0].regions[0][0].Rmax) break;
    38       if (stars[j].average.R <  here[0].regions[0][0].Rmin) break;
    39       if (stars[j].average.D <  here[0].regions[0][0].Dmin) continue;
    40       if (stars[j].average.D >= here[0].regions[0][0].Dmax) continue;
    41       stars[j].found = -2;
     41    for (j = i; j < catalog->Naverage; j++) {
     42      if (catalog->average[j].R >= here[0].regions[0][0].Rmax) break;
     43      if (catalog->average[j].R <  here[0].regions[0][0].Rmin) break;
     44      if (catalog->average[j].D <  here[0].regions[0][0].Dmin) continue;
     45      if (catalog->average[j].D >= here[0].regions[0][0].Dmax) continue;
     46      catalog->found_t[j] = -2;
    4247    }
    4348    list[0].regions[Nr] = here[0].regions[0];
    … …  
    5459
    5560  /* reset to -1 for all stars: required start for find_match_refstars */
    56   for (i = 0; i < Nstars; i++) stars[i].found = -1;
     61  free (catalog->found_t);
     62  catalog->found_t = NULL;
    5763  return (list);
    5864}
    … …  
    7076
    7177*/
    72 # endif

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

    r38553 r38593  
    5757      UpdateImageIDs (newcat, images, Nimages);
    5858
    59       // if ((DUMP != NULL) && !strcmp (DUMP, "rawstars")) dump_rawstars (stars, Nstars);
    6059      for (i = 0; i < Nimages; i++) {
    6160        newlist = SkyListByImage (sky, -1, &images[i]);
    … …  
    6665      ImageOptions (&options, images, Nimages);
    6766      break;
    68 # if (0)
    6967    case ADDSTAR_MODE_REFLIST:
    70       stars = grefstars (argv[1], options.photcode, &Nstars);
    71       skylist = SkyListForStars (sky, -1, stars, Nstars);
     68      newcat = grefstars (argv[1], options.photcode);
     69      skylist = SkyListForStars (sky, -1, newcat);
    7270      break;
    7371    case ADDSTAR_MODE_REFCAT:
    … …  
    8078      skylist[0].ownElements = FALSE;
    8179      break;
    82 # endif
    8380
    8481    default:
    … …  
    146143    switch (options.mode) {
    147144      case ADDSTAR_MODE_IMAGE:
    148         // Nsubset = Nstars;
    149145        if (options.closest) {
    150146          Nmatch += find_matches_closest (skylist[0].regions[i], newcat, &catalog, options);
    … …  
    153149        }
    154150        break;
    155 # if (0)
     151
    156152      case ADDSTAR_MODE_REFCAT:
    157         stars = greference (argv[1], skylist[0].regions[i], options.photcode, &Nstars);
     153        newcat = greference (argv[1], skylist[0].regions[i], options.photcode);
    158154      case ADDSTAR_MODE_REFLIST: {
    159         Stars **subset = find_subset (skylist[0].regions[i], stars, Nstars, &Nsubset);
    160155        if (options.closest) {
    161           Nmatch += find_matches_closest_refstars (skylist[0].regions[i], subset, Nsubset, &catalog, options);
     156          Nmatch += find_matches_closest_refstars (skylist[0].regions[i], newcat, &catalog, options);
    162157        } else {
    163           Nmatch += find_matches_refstars (skylist[0].regions[i], subset, Nsubset, &catalog, options);
     158          Nmatch += find_matches_refstars (skylist[0].regions[i], newcat, &catalog, options);
    164159        }
    165         if (Nsubset) free (subset);
    166160        break;
    167161      }
    168 # endif
    169162      default:
    170163        abort();

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches.c

    r38553 r38593  
    33int find_matches (SkyRegion *region, Catalog *srccat, Catalog *tgtcat, AddstarClientOptions options) {
    44
    5   off_t i, j, n, N, J, Nstars, mSrc;
     5  off_t i, j, J, Nstars, mSrc;
    66  double RADIUS, RADIUS2;
    77  double *X1, *Y1, *X2, *Y2;
    … …  
    99  off_t *N1, *N2;
    1010  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
    11   int Nsecfilt, Nsec, status;
    1211  unsigned int objID, catID;
    1312  Coords tcoords;
     13  int Nsecfilt, Nsec;
    1414
    1515  /* photcode data - must by of type DEP; options.photcode is equiv photcode for all input
    … …  
    5858  catID = tgtcat[0].catID;
    5959
    60   /* project onto rectilinear grid with 1 arcsec pixels. the choice of ZEA projection has the
    61    * advantage that every point in R,D has a mapping to a unique X,Y.  However, note that not all
    62    * possible X,Y points map back to R,D and the local plate scale changes substantially far from
    63    * the projection pole.  a better mapping might be ARC, not yet implemented (see
    64    * coordops.update.c).  We use the center of the region (tgtcat) for crval1,2.
     60  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has
     61   * the advantage that every point in R,D has a mapping to a unique X,Y.  However, note
     62   * that not all possible X,Y points map back to R,D and the local plate scale changes
     63   * substantially far from the projection pole.  We use the center of the region (tgtcat)
     64   * for crval1,2.
    6565   */
    6666  InitCoords (&tcoords, "DEC--ARC");
    6767  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
    68   if (region[0].Dmax < 90) {
     68  if ((region[0].Dmax < 90) && (region[0].Dmin > -90)) {
    6969    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
    7070  } else {
    71     tcoords.crval2 = 90.0;
     71    tcoords.crval2 = (region[0].Dmax >= 90) ? 90.0 : -90.0;
    7272  }
    7373  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    … …  
    7676  Nstars = 0;
    7777  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);
     78    int status = RD_to_XY (&X1[Nstars], &Y1[Nstars], srccat->average[i].R, srccat->average[i].D, &tcoords);
    7979    if (!status) continue;
    8080    N1[Nstars] = i;
    … …  
    101101  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    102102
    103   /* set up pointers for linked list of measure */
     103  /* set up pointers for linked list of measurements */
    104104  if (tgtcat[0].sorted && (tgtcat[0].Nmeasure == tgtcat[0].Nmeasure_disk)) {
     105    // this version is only valid if we have done a full catalog load, and if the catalog
     106    // is sorted while processed
    105107    next_meas = init_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
    106108  } else {
    … …  
    109111
    110112  /* choose a radius for matches (defined in args.c or ImageOptions.c) */
    111   RADIUS = options.radius;
     113  if (options.radius == 0) {
     114    RADIUS = 2.0; /* hardwired default for refstars */
     115  } else {
     116    RADIUS = options.radius; /* provided by config */
     117  }
    112118  RADIUS2 = RADIUS*RADIUS;
    113119
    114   /** find matched stars **/
     120  /****************** find matched stars ********************/
     121
    115122  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    116123    if (!finite(X1[i]) || !finite(Y1[i])) {
    … …  
    144151      /*** a match is found, add to average, measure ***/
    145152      Nmatch ++;
    146       n = N2[J];
    147       N = N1[i];
     153      off_t n = N2[J];
     154      off_t N = N1[i];
    148155
    149156      /* make sure there is space for next entry */
    … …  
    161168       
    162169        /** add measurements for this star **/
    163        
     170
    164171        // set the new measurements
    165172        tgtcat[0].measure[Nmeas]          = srccat->measure[mSrc + mOff];
    … …  
    231238          tgtcat[0].found_t[n] = Nmeas;
    232239        }
    233         /* Nm is updated, but not written out in -update mode (for existing entries)
    234            Nm is recalculated in build_meas_links if loaded table is not sorted */
    235240        tgtcat[0].average[n].Nmeasure ++;
    236241        Nmeas ++;
    … …  
    245250  }
    246251
    247 /* incorporate unmatched image stars, if this star is in field of this tgtcat */
    248 /* these new entries are all written out in UPDATE mode */
    249   for (i = 0; (i < Nstars) && !options.only_match; i ++) {
     252  /* incorporate unmatched image stars, if this star is in field of this tgtcat */
     253  /* these new entries are all written out in UPDATE mode */
     254  for (i = 0; (i < Nstars) && !options.only_match; i++) {
    250255    /* make sure there is space for next entry */
    251256    if (Nmeas >= NMEAS - srccat->average[i].Nmeasure) {
    … …  
    324329  }
    325330
    326 /* note stars which have been found in this tgtcat */
     331  /* note stars which have been found in this tgtcat */
    327332  for (i = 0; i < srccat->Nmeasure; i++) {
    328333    if (srccat->found_t[i] > -1) {
    … …  
    331336  }
    332337
    333 /* check if the tgtcat has changed?  if no change, no need to write */
     338  /* check if the tgtcat has changed?  if no change, no need to write */
    334339  tgtcat[0].objID    = objID; // new max value, save on tgtcat close
    335340  tgtcat[0].Naverage = Nave;
    … …  
    341346  free (Y1);
    342347  free (N1);
    343   free (N2);
    344348  free (X2);
    345349  free (Y2);
     350  free (N2);
    346351  free (next_meas);
    347352

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

    r38553 r38593  
    33int find_matches_closest (SkyRegion *region, Catalog *srccat, Catalog *tgtcat, AddstarClientOptions options) {
    44
    5   off_t i, j, n, N, J, Jmin, status, Nstars, mSrc;
     5  off_t i, j, J, Jmin, status, Nstars, mSrc;
    66  double RADIUS, RADIUS2, Rmin;
    77  double *X1, *Y1, *X2, *Y2;
    … …  
    99  off_t *N1, *N2;
    1010  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS;
    11   int Nsecfilt, Nsec;
    1211  unsigned int objID, catID;
    1312  Coords tcoords;
     13  int Nsecfilt, Nsec;
    1414
    1515  /* photcode data - must by of type DEP; options.photcode is equiv photcode for all input
    … …  
    6262   * the advantage that every point in R,D has a mapping to a unique X,Y.  However, note
    6363   * that not all possible X,Y points map back to R,D and the local plate scale changes
    64    * far from the projection pole. We use the center of the region (tgtcat) for crval1,2.
     64   * substantially far from the projection pole. We use the center of the region (tgtcat)
     65   * for crval1,2.
    6566   */
    6667  InitCoords (&tcoords, "DEC--ARC");
    6768  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
    68   if (region[0].Dmax < 90) {
     69  if ((region[0].Dmax < 90) && (region[0].Dmin > -90)) {
    6970    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
    7071  } else {
    71     tcoords.crval2 = 90.0;
     72    tcoords.crval2 = (region[0].Dmax >= 90) ? 90.0 : -90.0;
    7273  }
    7374  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    … …  
    114115
    115116  /* choose a radius for matches */
    116   RADIUS = options.radius; /* provided by config */
     117  if (options.radius == 0) {
     118    RADIUS = 2.0; /* hardwired default for refstars */
     119  } else {
     120    RADIUS = options.radius; /* provided by config */
     121  }
    117122  RADIUS2 = RADIUS*RADIUS;
    118123
    … …  
    178183    /*** a match is found, add to average, measure ***/
    179184    Nmatch ++;
    180     n = N2[Jmin];
    181     N = N1[i];
     185    off_t n = N2[Jmin];
     186    off_t N = N1[i];
    182187
    183188    /* make sure there is space for next entry */
    … …  
    355360      next_meas[Nmeas-1] = (mSrc < srccat->average[i].Nmeasure - 1) ? Nmeas : -1;
    356361    }
    357     srccat->found_t[i] = Nave;
    358 
    359362    // if we have lensing data, insert that as well
    360363    if (srccat->lensing) {
    … …  
    374377      }
    375378    }
     379    srccat->found_t[i] = Nave;
    376380    Nave ++;
    377381  }
    … …  
    407411  free (Y1);
    408412  free (N1);
    409   free (N2);
    410413  free (X2);
    411414  free (Y2);
     415  free (N2);
    412416  free (next_meas);
    413417  free (next_lens);

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_closest_refstars.c

    r38467 r38593  
    11# include "addstar.h"
    22
    3 /* differences with find_matches_closest():
    4  * accepts Stars **stars (vs Stars *stars)
    5  * NREFSTAR_GROUP vs NSTAR_GROUP
    6  * does not update lensing fields
    7  * allows option to replace existing entries
    8  * allows multiple photcodes
    9 
    10  */
    11 
    12 int find_matches_closest_refstars (SkyRegion *region, Stars **stars, unsigned int NstarsIn, Catalog *catalog, AddstarClientOptions options) {
    13 
    14   off_t i, j, n, N, J, Jmin, status;
     3int find_matches_closest_refstars (SkyRegion *region, Catalog *srccat, Catalog *tgtcat, AddstarClientOptions options) {
     4
     5  off_t i, j, J, Jmin, status, mSrc;
    156  double RADIUS, RADIUS2, Rmin;
    167  double *X1, *Y1, *X2, *Y2;
    178  double dX, dY, dR;
    18   off_t *N1, *N2, *next_meas;
     9  off_t *N1, *N2;
    1910  off_t Nave, NAVE, Nmeas, NMEAS, Nstars, Nmatch;
    2011  unsigned int objID, catID;
    … …  
    2213  int Nsecfilt;
    2314
    24   switch (NREFSTAR_GROUP) {
    25     case 1:
    26     case 3:
    27     case 4:
    28     case 6:
    29       break;
    30     default:
    31       fprintf (stderr, "ERROR: NREFSTAR_GROUP NOT SET!\n");
    32       exit (1);
    33   }
    34 
    3515  /* photcode data -- should not have to modify secfilt / average */
    3616  Nsecfilt = GetPhotcodeNsecfilt ();
    3717
    38   /** allocate local arrays (stars) **/
    39   ALLOCATE (X1, double, NstarsIn);
    40   ALLOCATE (Y1, double, NstarsIn);
    41   ALLOCATE (N1, off_t,  NstarsIn);
    42 
    43   /** allocate local arrays (catalog) **/
    44   NAVE = Nave = catalog[0].Naverage;
     18  /** allocate local arrays (srccat) **/
     19  ALLOCATE (X1, double, srccat->Naverage);
     20  ALLOCATE (Y1, double, srccat->Naverage);
     21  ALLOCATE (N1, off_t,  srccat->Naverage);
     22
     23  if (!srccat->found_t) {
     24    ALLOCATE (srccat->found_t, off_t, srccat->Naverage);
     25    for (i = 0; i < srccat->Naverage; i++) {
     26      srccat->found_t[i] = -1;
     27    }
     28  }
     29
     30  /** allocate local arrays (tgtcat) **/
     31  NAVE = Nave = tgtcat[0].Naverage;
    4532  ALLOCATE (X2, double, NAVE);
    4633  ALLOCATE (Y2, double, NAVE);
    4734  ALLOCATE (N2, off_t,  NAVE);
    48   ALLOCATE (catalog[0].found_t, off_t, NAVE);
     35
     36  if (!tgtcat->found_t) {
     37    ALLOCATE (tgtcat[0].found_t, off_t, NAVE);
     38    for (i = 0; i < NAVE; i++) {
     39      tgtcat->found_t[i] = -1;
     40    }
     41  }
    4942  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
    5043
    5144  /* internal counters */
    5245  Nmatch = 0;
    53   NMEAS = Nmeas = catalog[0].Nmeasure;
     46  NMEAS = Nmeas = tgtcat[0].Nmeasure;
    5447 
    55   // current max obj ID for this catalog
    56   objID = catalog[0].objID;
    57   catID = catalog[0].catID;
     48  off_t *next_meas = NULL;
     49
     50  // current max obj ID for this tgtcat
     51  objID = tgtcat[0].objID;
     52  catID = tgtcat[0].catID;
    5853
    5954  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has
    6055   * the advantage that every point in R,D has a mapping to a unique X,Y.  However, note
    6156   * that not all possible X,Y points map back to R,D and the local plate scale changes
    62    * far from the projection pole. We use the center of the region (catalog) for crval1,2.
     57   * substantially far from the projection pole. We use the center of the region (tgtcat)
     58   * for crval1,2.
    6359   */
    6460  InitCoords (&tcoords, "DEC--ARC");
    6561  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
    66   if (region[0].Dmax < 90) {
     62  if ((region[0].Dmax < 90) && (region[0].Dmin > -90)) {
    6763    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
    6864  } else {
    69     tcoords.crval2 = 90.0;
     65    tcoords.crval2 = (region[0].Dmax >= 90) ? 90.0 : -90.0;
    7066  }
    7167  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    … …  
    7369  /* build spatial index (RA sort) referencing input array sequence */
    7470  Nstars = 0;
    75   for (i = 0; i < NstarsIn; i++) {
    76     status = RD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i][0].average.R, stars[i][0].average.D, &tcoords);
     71  for (i = 0; i < srccat->Naverage; i++) {
     72    status = RD_to_XY (&X1[Nstars], &Y1[Nstars], srccat->average[i].R, srccat->average[i].D, &tcoords);
    7773    if (!status) continue;
    7874    N1[Nstars] = i;
    … …  
    8076  }
    8177  if (Nstars < 1) {
    82     if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
     78    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", tgtcat[0].filename);
    8379    free (X1);
    8480    free (Y1);
    … …  
    9389  /* build spatial index (RA sort) */
    9490  for (i = 0; i < Nave; i++) {
    95     RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
     91    RD_to_XY (&X2[i], &Y2[i], tgtcat[0].average[i].R, tgtcat[0].average[i].D, &tcoords);
    9692    N2[i] = i;
    97     catalog[0].found_t[N2[i]] = -1;
     93    tgtcat[0].found_t[N2[i]] = -1;
    9894  }
    9995  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    10096
    10197  /* set up pointers for linked list of measure */
    102   if (catalog[0].sorted && (catalog[0].Nmeasure == catalog[0].Nmeasure_disk)) {
    103     // this version is only valid if we have done a full catalog load, and if the catalog
     98  if (tgtcat[0].sorted && (tgtcat[0].Nmeasure == tgtcat[0].Nmeasure_disk)) {
     99    // this version is only valid if we have done a full tgtcat load, and if the tgtcat
    104100    // is sorted while processed
    105     next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
     101    next_meas = init_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
    106102  } else {
    107     next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
     103    next_meas = build_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
    108104  }   
    109105
    … …  
    116112  RADIUS2 = RADIUS*RADIUS;
    117113
    118   /** find matched stars **/
    119   // XXX could use NREFSTAR_GROUP to do this match more quicky
     114  /****************** find matched stars ********************/
     115
    120116  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    121117    if (!finite(X1[i]) || !finite(Y1[i])) {
    … …  
    140136    }
    141137
    142     if (stars[N1[i]][0].found != -1) {
    143         /* this star has already been assigned to an object in this or another catalog */
     138    // XXX check that this is allocated
     139    if (srccat->found_t[N1[i]] != -1) {
     140        /* this star has already been assigned to an object in this or another tgtcat */
    144141        i++;
    145142        continue;
    … …  
    169166    /*** a match is found, add to average, measure ***/
    170167    Nmatch ++;
    171     n = N2[Jmin];
    172     N = N1[i];
    173 
    174     /** in replace mode, search for entry and replace values M, dM, R, D */
    175     // XXX this fails for unsorted catalogs, right?
    176     if (options.replace && replace_match (&catalog[0].average[n], catalog[0].measure, stars[N])) {
    177         i++;
    178         continue;
    179     }
     168    off_t n = N2[Jmin];
     169    off_t N = N1[i];
    180170
    181171    /* make sure there is space for next entry */
    182     if (Nmeas >= NMEAS) {
    183       NMEAS = Nmeas + 1000;
     172    if (Nmeas >= NMEAS - srccat->average[N].Nmeasure) {
     173      NMEAS = Nmeas + 1000 + srccat->average[N].Nmeasure;
    184174      REALLOCATE (next_meas, off_t, NMEAS);
    185       REALLOCATE (catalog[0].measure, Measure, NMEAS);
    186     }
    187 
    188     /* add to end of measurement list */
    189     add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
    190        
     175      REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
     176    }
     177
    191178    /** add measurements for this star **/
    192 
    193     // set the new measurements
    194     catalog[0].measure[Nmeas]          = stars[N][0].measure;
    195 
    196     catalog[0].measure[Nmeas].dbFlags  = 0;
    197     catalog[0].measure[Nmeas].averef   = n;
    198     catalog[0].measure[Nmeas].objID    = catalog[0].average[n].objID;
    199     catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    200 
    201     float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]);
    202 
    203     // rationalize dR:
    204     if (dRoff > +180.0*3600.0) {
    205       // average on high end of boundary, move star up
    206       catalog[0].measure[Nmeas].R += 360.0;
    207       dRoff -= 360.0*3600.0;
    208     }
    209     if (dRoff < -180.0*3600.0) {
    210       // average on low end of boundary, move star down
    211       catalog[0].measure[Nmeas].R -= 360.0;
    212       dRoff += 360.0*3600.0;
    213     }
    214 
    215     // the reference data may not have per-object time specified
    216     catalog[0].measure[Nmeas].t        = (TIMEREF == 0) ? stars[N][0].measure.t      : TIMEREF; /** careful : time_t vs e_time **/
    217     catalog[0].measure[Nmeas].t_msec   = (TIMEREF == 0) ? stars[N][0].measure.t_msec : TIMEREF; /** careful : time_t vs e_time **/
    218 
    219     // we can choose to accept the proper-motion and parallax from the reference catalog
    220     if (ACCEPT_MOTION) {
    221       catalog[0].average[n].dR         = stars[N][0].average.dR;
    222       catalog[0].average[n].dD         = stars[N][0].average.dD;
    223       catalog[0].average[n].uR         = stars[N][0].average.uR;
    224       catalog[0].average[n].uD         = stars[N][0].average.uD;
    225       catalog[0].average[n].duR        = stars[N][0].average.duR;
    226       catalog[0].average[n].duD        = stars[N][0].average.duD;
    227       catalog[0].average[n].P          = stars[N][0].average.P;
    228       catalog[0].average[n].dP         = stars[N][0].average.dP;
    229       catalog[0].average[n].Tmean      = stars[N][0].average.Tmean;
    230     }
     179    off_t mOff = srccat->average[N].measureOffset;
     180    for (mSrc = 0; mSrc < srccat->average[N].Nmeasure; mSrc++) {
     181      /** in replace mode, search for entry and replace values M, dM, R, D */
     182      // XXX this fails for unsorted tgtcats, right?
     183      if (options.replace && replace_match (&tgtcat[0].average[n], tgtcat[0].measure, &srccat->measure[mSrc + mOff], &srccat->found_t[N])) continue;
     184
     185      /* add to end of measurement list */
     186      add_meas_link (&tgtcat[0].average[n], next_meas, Nmeas, NMEAS);
     187     
     188      // set the new measurements
     189      tgtcat[0].measure[Nmeas]          = srccat->measure[mSrc + mOff];
     190
     191      // measure now carries R,D (not dR,dD)
     192      tgtcat[0].measure[Nmeas].dbFlags  = 0;
     193      tgtcat[0].measure[Nmeas].averef   = n;
     194      tgtcat[0].measure[Nmeas].objID    = tgtcat[0].average[n].objID;
     195      tgtcat[0].measure[Nmeas].catID    = tgtcat[0].catID;
     196
     197      float dRoff = dvoOffsetR(&tgtcat[0].measure[Nmeas], &tgtcat[0].average[n]);
     198
     199      // rationalize dR:
     200      if (dRoff > +180.0*3600.0) {
     201        // average on high end of boundary, move star up
     202        tgtcat[0].measure[Nmeas].R += 360.0;
     203        dRoff -= 360.0*3600.0;
     204      }
     205      if (dRoff < -180.0*3600.0) {
     206        // average on low end of boundary, move star down
     207        tgtcat[0].measure[Nmeas].R -= 360.0;
     208        dRoff += 360.0*3600.0;
     209      }
     210      if (fabs(dRoff) > 10*RADIUS) {
     211        // take declination into account and check again.
     212        double cosD = cos(RAD_DEG*tgtcat[0].average[n].D);
     213        if (fabs(dRoff*cosD) > 10*RADIUS) {
     214          fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n",
     215                   tgtcat[0].average[n].R, tgtcat[0].average[n].D,
     216                   srccat->average[N].R, srccat->average[N].D,
     217                   X1[i], X2[Jmin],
     218                   Y1[i], Y2[Jmin]);
     219          // XXX abort on this? -- this is a bad failure...
     220        }
     221      }
     222
     223      // the reference data may not have per-object time specified
     224      if (TIMEREF) {
     225        tgtcat[0].measure[Nmeas].t        = TIMEREF;
     226        tgtcat[0].measure[Nmeas].t_msec   = 0;
     227      }
     228
     229      // we can choose to accept the proper-motion and parallax from the reference tgtcat
     230      if (ACCEPT_MOTION) {
     231        tgtcat[0].average[n].dR         = srccat->average[N].dR;
     232        tgtcat[0].average[n].dD         = srccat->average[N].dD;
     233        tgtcat[0].average[n].uR         = srccat->average[N].uR;
     234        tgtcat[0].average[n].uD         = srccat->average[N].uD;
     235        tgtcat[0].average[n].duR        = srccat->average[N].duR;
     236        tgtcat[0].average[n].duD        = srccat->average[N].duD;
     237        tgtcat[0].average[n].P          = srccat->average[N].P;
     238        tgtcat[0].average[n].dP         = srccat->average[N].dP;
     239        tgtcat[0].average[n].Tmean      = srccat->average[N].Tmean;
     240      }
    231241
    232242    /** don't update average / secfilt values for REF photcodes **/
    233243
    234     /*** handle multiple stars */
    235     /* this catalog star matches more than one image star */
    236     // XXX should this be an average flag?
    237     if (catalog[0].found_t[n] > -1) {
    238       catalog[0].measure[catalog[0].found_t[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
    239       catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
    240     } else {
    241       catalog[0].found_t[n] = Nmeas;
    242     }
    243 
    244     stars[N][0].found = Nmeas;
    245     catalog[0].average[n].Nmeasure ++;
    246     Nmeas ++;
    247 
    248     // make this optional:
    249     // update_coords (&catalog[0].average[n], &catalog[0].measure[0], next_meas);
     244      tgtcat[0].average[n].Nmeasure ++;
     245      Nmeas ++;
     246    }
     247
     248    srccat->found_t[N] = n;
     249    tgtcat[0].found_t[n] = -1;
    250250    i++;
    251251  }
     252
     253  /*************** add unmatched stars *************************/
    252254
    253255  /* Incorporate unmatched refcat stars.  Skip this step if we want to
    … …  
    255257     reference up-to-date with known stars only */
    256258
    257   for (i = 0; (i < Nstars) && !options.only_match; i+=NREFSTAR_GROUP) {
     259  for (i = 0; (i < Nstars) && !options.only_match; i++) {
    258260    /* make sure there is space for next entry */
    259     if (Nmeas >= NMEAS - NREFSTAR_GROUP) {
    260       NMEAS = Nmeas + 1000;
     261    if (Nmeas >= NMEAS - srccat->average[i].Nmeasure) {
     262      NMEAS = Nmeas + 1000 + srccat->average[i].Nmeasure;
    261263      REALLOCATE (next_meas, off_t, NMEAS);
    262       REALLOCATE (catalog[0].measure, Measure, NMEAS);
     264      REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
    263265    }
    264266    if (Nave >= NAVE) {
    265267      NAVE = Nave + 1000;
    266       REALLOCATE (catalog[0].average, Average, NAVE);
    267       REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
    268     }
    269 
    270     N = N1[i];
    271     if (stars[N][0].found >= 0) continue;
    272     if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue;
    273 
    274     dvo_average_init (&catalog[0].average[Nave]);
    275     catalog[0].average[Nave].R             = stars[N][0].average.R;
    276     catalog[0].average[Nave].D             = stars[N][0].average.D;
    277 
    278     catalog[0].average[Nave].Nmeasure      = NREFSTAR_GROUP;
    279     catalog[0].average[Nave].measureOffset = Nmeas;
    280     catalog[0].average[Nave].objID         = objID;
    281     catalog[0].average[Nave].catID         = catID;
     268      REALLOCATE (tgtcat[0].average, Average, NAVE);
     269      REALLOCATE (tgtcat[0].found_t, off_t, NAVE);
     270      if (tgtcat[0].secfilt) {
     271        // we only update the secfilt table if it has been allocated for output
     272        REALLOCATE (tgtcat[0].secfilt, SecFilt, NAVE*tgtcat[0].Nsecfilt);
     273      }
     274    }
     275
     276    if (srccat->found_t[i] != -1) continue;
     277    if (!IN_REGION (srccat->average[i].R, srccat->average[i].D)) continue;
     278
     279    dvo_average_init (&tgtcat[0].average[Nave]);
     280    tgtcat[0].average[Nave].R              = srccat->average[i].R;
     281    tgtcat[0].average[Nave].D              = srccat->average[i].D;
     282
     283    tgtcat[0].average[Nave].Nmeasure       = srccat->average[i].Nmeasure;
     284    tgtcat[0].average[Nave].measureOffset  = Nmeas;
     285    tgtcat[0].average[Nave].objID          = objID;
     286    tgtcat[0].average[Nave].catID          = catID;
    282287
    283288    if (PSPS_ID) {
    284         catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D);
     289        tgtcat[0].average[Nave].extID = CreatePSPSObjectID(tgtcat[0].average[Nave].R, tgtcat[0].average[Nave].D);
    285290    }
    286291
    287292    if (ACCEPT_MOTION) {
    288       catalog[0].average[Nave].dR          = stars[N][0].average.dR;
    289       catalog[0].average[Nave].dD          = stars[N][0].average.dD;
    290       catalog[0].average[Nave].uR          = stars[N][0].average.uR;
    291       catalog[0].average[Nave].uD          = stars[N][0].average.uD;
    292       catalog[0].average[Nave].duR         = stars[N][0].average.duR;
    293       catalog[0].average[Nave].duD         = stars[N][0].average.duD;
    294       catalog[0].average[Nave].P           = stars[N][0].average.P;
    295       catalog[0].average[Nave].dP          = stars[N][0].average.dP;
    296       catalog[0].average[Nave].Tmean       = stars[N][0].average.Tmean;
     293      tgtcat[0].average[Nave].dR           = srccat->average[i].dR;
     294      tgtcat[0].average[Nave].dD           = srccat->average[i].dD;
     295      tgtcat[0].average[Nave].uR           = srccat->average[i].uR;
     296      tgtcat[0].average[Nave].uD           = srccat->average[i].uD;
     297      tgtcat[0].average[Nave].duR          = srccat->average[i].duR;
     298      tgtcat[0].average[Nave].duD          = srccat->average[i].duD;
     299      tgtcat[0].average[Nave].P            = srccat->average[i].P;
     300      tgtcat[0].average[Nave].dP           = srccat->average[i].dP;
     301      tgtcat[0].average[Nave].Tmean        = srccat->average[i].Tmean;
    297302    }
    298303
    … …  
    300305
    301306    for (j = 0; j < Nsecfilt; j++) {
    302       dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
    303     }
    304 
    305     for (j = 0; j < NREFSTAR_GROUP; j++) {
    306       N = N1[i + j];
    307 
    308       catalog[0].measure[Nmeas]          = stars[N][0].measure;
    309       catalog[0].measure[Nmeas].dbFlags  = 0;
    310       catalog[0].measure[Nmeas].averef   = Nave;
    311       catalog[0].measure[Nmeas].objID    = catalog[0].average[Nave].objID;
    312       catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    313 
    314       catalog[0].measure[Nmeas].t        = (stars[N][0].measure.t == 0) ? TIMEREF : stars[N][0].measure.t;      /** careful : time_t vs e_time **/
    315       catalog[0].measure[Nmeas].t_msec   = (stars[N][0].measure.t == 0) ?       0 : stars[N][0].measure.t_msec; /** careful : time_t vs e_time **/
    316 
    317       stars[N][0].found = Nmeas;
    318       next_meas[Nmeas] = -1;
     307      dvo_secfilt_init (&tgtcat[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
     308    }
     309
     310    off_t mOff = srccat->average[i].measureOffset;
     311    for (mSrc = 0; mSrc < srccat->average[i].Nmeasure; mSrc++) {
     312      // supply the measurements from this detection
     313      dvo_measure_init (&tgtcat[0].measure[Nmeas]);
     314      tgtcat[0].measure[Nmeas] = srccat->measure[mSrc + mOff];
     315
     316      // the following measure elements cannot be set until here:
     317      tgtcat[0].measure[Nmeas].dbFlags  = 0;
     318      tgtcat[0].measure[Nmeas].averef   = Nave;
     319      tgtcat[0].measure[Nmeas].objID    = tgtcat[0].average[Nave].objID;
     320      tgtcat[0].measure[Nmeas].catID    = tgtcat[0].catID;
     321
     322      if (TIMEREF) {
     323        tgtcat[0].measure[Nmeas].t      = TIMEREF;
     324        tgtcat[0].measure[Nmeas].t_msec = 0;
     325      }
     326
    319327      Nmeas ++;
    320     }
    321     for (j = 0; j < NREFSTAR_GROUP - 1; j++) {
    322         next_meas[Nmeas - NREFSTAR_GROUP + j] = Nmeas - NREFSTAR_GROUP + j + 1;
    323     }
     328
     329      // update the next_meas pointer for this entry (last one for this star is -1)
     330      next_meas[Nmeas-1] = (mSrc < srccat->average[i].Nmeasure - 1) ? Nmeas : -1;
     331    }
     332    srccat->found_t[i] = Nave;
    324333    Nave ++;
    325334  }
    326335     
    327   REALLOCATE (catalog[0].average, Average, Nave);
    328   REALLOCATE (catalog[0].measure, Measure, Nmeas);
     336  REALLOCATE (tgtcat[0].average, Average, Nave);
     337  REALLOCATE (tgtcat[0].measure, Measure, Nmeas);
    329338
    330339  // XXX allow for unsorted output?
    331   catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
    332   catalog[0].sorted = TRUE;
    333 
    334   /* note stars which have been found in this catalog */
    335   for (i = 0; i < Nstars; i++) {
    336     if (stars[i][0].found > -1) {
    337       stars[i][0].found = -2;
    338     }
    339   }
    340 
    341   catalog[0].objID    = objID; // new max value, save on catalog close
    342   catalog[0].Naverage = Nave;
    343   catalog[0].Nmeasure = Nmeas;
    344   catalog[0].Nsecfilt_mem = Nave*Nsecfilt;
     340  tgtcat[0].measure = sort_measure (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas, next_meas);
     341  tgtcat[0].sorted = TRUE;
     342
     343  /* note stars which have been found in this tgtcat */
     344  for (i = 0; i < srccat->Naverage; i++) {
     345    if (srccat->found_t[i] > -1) {
     346      srccat->found_t[i] = -2;
     347    }
     348  }
     349
     350  /* check if the tgtcat has changed?  if no change, no need to write */
     351  tgtcat[0].objID    = objID; // new max value, save on tgtcat close
     352  tgtcat[0].Naverage = Nave;
     353  tgtcat[0].Nmeasure = Nmeas;
     354  tgtcat[0].Nsecfilt_mem = Nave*Nsecfilt;
    345355  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
    346356
    … …  
    348358  free (Y1);
    349359  free (N1);
    350   free (N2);
    351360  free (X2);
    352361  free (Y2);
     362  free (N2);
    353363  free (next_meas);
    354364
    355365  return (Nmatch);
    356366}
     367
     368/* differences with find_matches_closest():
     369 * accepts Stars **stars (vs Stars *stars)
     370 * NREFSTAR_GROUP vs NSTAR_GROUP
     371 * does not update lensing fields
     372 * allows option to replace existing entries
     373 * allows multiple photcodes
     374
     375 */
     376

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_refstars.c

    r38467 r38593  
    11# include "addstar.h"
    22
    3 int find_matches_refstars (SkyRegion *region, Stars **stars, unsigned int Nstars, Catalog *catalog, AddstarClientOptions options) {
    4 
    5   off_t i, j, n, N, J;
     3// see bottom of this file for differences between find_matches_refstars.c and find_matches.c
     4int find_matches_refstars (SkyRegion *region, Catalog *srccat, Catalog *tgtcat, AddstarClientOptions options) {
     5
     6  off_t i, j, J, Nstars, mSrc;
    67  double RADIUS, RADIUS2;
    78  double *X1, *Y1, *X2, *Y2;
    89  double dX, dY, dR;
    9   off_t *N1, *N2,  *next_meas;
     10  off_t *N1, *N2;
    1011  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
    1112  unsigned int objID, catID;
    … …  
    1314  int Nsecfilt;
    1415
    15   switch (NREFSTAR_GROUP) {
    16     case 1:
    17     case 3:
    18     case 4:
    19       break;
    20     default:
    21       fprintf (stderr, "ERROR: NREFSTAR_GROUP NOT SET!\n");
    22       exit (1);
    23   }
    24 
    2516  /* photcode data -- should not have to modify secfilt / average */
    2617  Nsecfilt = GetPhotcodeNsecfilt ();
    2718
    2819  /** allocate local arrays (stars) **/
    29   ALLOCATE (X1, double, Nstars);
    30   ALLOCATE (Y1, double, Nstars);
    31   ALLOCATE (N1, off_t,  Nstars);
    32 
    33   /** allocate local arrays (catalog) **/
    34   Nave = catalog[0].Naverage;
    35   NAVE = Nave + 1000;
     20  ALLOCATE (X1, double, srccat->Naverage);
     21  ALLOCATE (Y1, double, srccat->Naverage);
     22  ALLOCATE (N1, off_t,  srccat->Naverage);
     23
     24  myAbort ("figure out if found_t needs to be tracking Naverage or Nmeasure");
     25  if (!srccat->found_t) {
     26    ALLOCATE (srccat->found_t, off_t, srccat->Naverage);
     27    for (i = 0; i < srccat->Naverage; i++) {
     28      srccat->found_t[i] = -1;
     29    }
     30  }
     31
     32  /** allocate local arrays (tgtcat) **/
     33  NAVE = Nave = tgtcat[0].Naverage;
    3634  ALLOCATE (X2, double, NAVE);
    3735  ALLOCATE (Y2, double, NAVE);
    3836  ALLOCATE (N2, off_t,  NAVE);
    39   ALLOCATE (catalog[0].found_t, off_t, NAVE);
    40   REALLOCATE (catalog[0].average, Average, NAVE);
    41   REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
     37
     38  if (!tgtcat->found_t) {
     39    ALLOCATE (tgtcat[0].found_t, off_t, NAVE);
     40    for (i = 0; i < NAVE; i++) {
     41      tgtcat->found_t[i] = -1;
     42    }
     43  }
    4244  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
    4345
     46  /* internal counters */
    4447  Nmatch = 0;
    45   NMEAS = Nmeas = catalog[0].Nmeasure;
     48  NMEAS = Nmeas = tgtcat[0].Nmeasure;
    4649 
    47   // current max obj ID for this catalog
    48   objID = catalog[0].objID;
    49   catID = catalog[0].catID;
    50 
    51   /* project onto rectilinear grid with 1 arcsec pixels, sort by X */
    52   /* reference for coords is catalog center */
    53   InitCoords (&tcoords, "DEC--ZEA");
     50  off_t *next_meas = NULL;
     51
     52  // current max obj ID for this tgtcat
     53  objID = tgtcat[0].objID;
     54  catID = tgtcat[0].catID;
     55
     56  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has
     57   * the advantage that every point in R,D has a mapping to a unique X,Y.  However, note
     58   * that not all possible X,Y points map back to R,D and the local plate scale changes
     59   * substantially far from the projection pole.  We use the center of the region (tgtcat)
     60   * for crval1,2.
     61   */
     62  InitCoords (&tcoords, "DEC--ARC");
    5463  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
    55   tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
     64  if ((region[0].Dmax < 90) && (region[0].Dmin > -90)) {
     65    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
     66  } else {
     67    tcoords.crval2 = (region[0].Dmax >= 90) ? 90.0 : -90.0;
     68  }
    5669  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    57  
    58   for (i = 0; i < Nstars; i++) {
    59     RD_to_XY (&X1[i], &Y1[i], stars[i][0].average.R, stars[i][0].average.D, &tcoords);
    60     N1[i] = i;
     70
     71  /* build spatial index (RA sort) referencing input array sequence */
     72  Nstars = 0;
     73  for (i = 0; i < srccat->Naverage; i++) {
     74    int status = RD_to_XY (&X1[Nstars], &Y1[Nstars], srccat->average[i].R, srccat->average[i].D, &tcoords);
     75    if (!status) continue;
     76    N1[Nstars] = i;
     77    Nstars ++;
     78  }
     79  if (Nstars < 1) {
     80    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", tgtcat[0].filename);
     81    free (X1);
     82    free (Y1);
     83    free (N1);
     84    free (X2);
     85    free (Y2);
     86    free (N2);
     87    return (0);
    6188  }
    6289  if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars);
    6390 
     91  /* build spatial index (RA sort) */
    6492  for (i = 0; i < Nave; i++) {
    65     RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
     93    RD_to_XY (&X2[i], &Y2[i], tgtcat[0].average[i].R, tgtcat[0].average[i].D, &tcoords);
    6694    N2[i] = i;
    67     catalog[0].found_t[N2[i]] = -1;
     95    tgtcat[0].found_t[N2[i]] = -1;
    6896  }
    6997  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
    7098
    7199  /* set up pointers for linked list of measurements */
    72   if (catalog[0].sorted && (catalog[0].Nmeasure == catalog[0].Nmeasure_disk)) {
     100  if (tgtcat[0].sorted && (tgtcat[0].Nmeasure == tgtcat[0].Nmeasure_disk)) {
    73101    // this version is only valid if we have done a full catalog load, and if the catalog
    74102    // is sorted while processed
    75     next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
     103    next_meas = init_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
    76104  } else {
    77     next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
     105    next_meas = build_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
    78106  }   
    79107
    80   /* choose a radius for matches */
     108  /* choose a radius for matches (defined in args.c or ImageOptions.c) */
    81109  if (options.radius == 0) {
    82110    RADIUS = 2.0; /* hardwired default for refstars */
    … …  
    86114  RADIUS2 = RADIUS*RADIUS;
    87115
    88   /** find matched stars **/
    89   // XXX could use NREFSTAR_GROUP to do this match more quicky
     116  /****************** find matched stars ********************/
     117
    90118  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    91119    if (!finite(X1[i]) || !finite(Y1[i])) {
    … …  
    97125      continue;
    98126    }
    99    
     127
     128    /* negative dX: j is too large */
    100129    dX = X1[i] - X2[j];
    101130    if (dX <= -1.02*RADIUS) {
    … …  
    103132      continue;
    104133    }
     134    /* positive dX, i is too large */
    105135    if (dX >= 1.02*RADIUS) {
    106136      j++;
    … …  
    108138    }
    109139
    110     /* negative dX: j is too large, positive dX, i is too large */
     140    /* within match range; look for matches */
    111141    for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) {
    112142      dX = X1[i] - X2[J];
    … …  
    117147      /*** a match is found, add to average, measure ***/
    118148      Nmatch ++;
    119       n = N2[J];
    120       N = N1[i];
    121 
    122       /** in replace mode, search for entry and replace values M, dM, R, D */
    123       if (options.replace && replace_match (&catalog[0].average[n], catalog[0].measure, stars[N])) continue;
    124 
    125       if (Nmeas == NMEAS) {
    126         NMEAS = Nmeas + 1000;
     149      off_t n = N2[J];
     150      off_t N = N1[i];
     151
     152      /* make sure there is space for next entry */
     153      if (Nmeas >= NMEAS - srccat->average[N].Nmeasure) {
     154        NMEAS = Nmeas + 1000 + srccat->average[N].Nmeasure;
    127155        REALLOCATE (next_meas, off_t, NMEAS);
    128         REALLOCATE (catalog[0].measure, Measure, NMEAS);
     156        REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
    129157      }
    130158
    131       /* add to end of measurement list */
    132       add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
     159      /** add measurements for this star **/
     160      off_t mOff = srccat->average[N].measureOffset;
     161      for (mSrc = 0; mSrc < srccat->average[N].Nmeasure; mSrc++) {
     162        /** in replace mode, search for entry and replace values M, dM, R, D */
     163        if (options.replace && replace_match (&tgtcat[0].average[n], tgtcat[0].measure, &srccat->measure[mSrc + mOff], &srccat->found_t[N])) continue;
     164
     165        /* add to end of measurement list */
     166        add_meas_link (&tgtcat[0].average[n], next_meas, Nmeas, NMEAS);
    133167       
    134       /** add measurements for this star **/
    135 
    136       // set the new measurements
    137       catalog[0].measure[Nmeas]          = stars[N][0].measure;
    138 
    139       /** measure now stores R,D **/
    140       catalog[0].measure[Nmeas].dbFlags  = 0;
    141       catalog[0].measure[Nmeas].averef   = n;
    142       catalog[0].measure[Nmeas].objID    = catalog[0].average[n].objID;
    143       catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    144 
    145       float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]);
    146 
    147       // rationalize dR:
    148       if (dRoff > +180.0*3600.0) {
     168        /** add measurements for this star **/
     169
     170        // set the new measurements
     171        tgtcat[0].measure[Nmeas]          = srccat->measure[mSrc + mOff];
     172       
     173        // the following measure elements cannot be set until here:
     174        tgtcat[0].measure[Nmeas].dbFlags  = 0;
     175        tgtcat[0].measure[Nmeas].averef   = n; // this must be an absolute sequence number, if partial average is loaded
     176        tgtcat[0].measure[Nmeas].objID    = tgtcat[0].average[n].objID;
     177        tgtcat[0].measure[Nmeas].catID    = tgtcat[0].catID;
     178
     179        float dRoff = dvoOffsetR(&tgtcat[0].measure[Nmeas], &tgtcat[0].average[n]);
     180
     181        // rationalize R:
     182        if (dRoff > +180.0*3600.0) {
    149183          // average on high end of boundary, move star up
    150           catalog[0].measure[Nmeas].R += 360.0;
     184          tgtcat[0].measure[Nmeas].R += 360.0;
    151185          dRoff -= 360.0*3600.0;
     186        }
     187        if (dRoff < -180.0*3600.0) {
     188          // average on low end of boundary, move star down
     189          tgtcat[0].measure[Nmeas].R -= 360.0;
     190          dRoff += 360.0*3600.0;
     191        }
     192        if (fabs(dRoff) > 10*RADIUS) {
     193          // take declination into account and check again.
     194          double cosD = cos(RAD_DEG*tgtcat[0].average[n].D);
     195          if (fabs(dRoff*cosD) > 10*RADIUS) {
     196            fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n",
     197                     tgtcat[0].average[n].R, tgtcat[0].average[n].D,
     198                     srccat->average[N].R, srccat->average[N].D,
     199                     X1[i], X2[J],
     200                     Y1[i], Y2[J]);
     201          }
     202        }
     203
     204        if (TIMEREF) {
     205          tgtcat[0].measure[Nmeas].t        = TIMEREF; /** careful : time_t vs e_time **/
     206          tgtcat[0].measure[Nmeas].t_msec   = 0;
     207        }
     208
     209        if (ACCEPT_MOTION) {
     210          tgtcat[0].average[n].uR         = srccat->average[N].uR;
     211          tgtcat[0].average[n].uD         = srccat->average[N].uD;
     212          tgtcat[0].average[n].duR        = srccat->average[N].duR;
     213          tgtcat[0].average[n].duD        = srccat->average[N].duD;
     214          tgtcat[0].average[n].P          = srccat->average[N].P;
     215          tgtcat[0].average[n].dP         = srccat->average[N].dP;
     216        }
     217
     218        /** don't update average / secfilt values for REF photcodes **/
     219
     220        myAbort ("fix logic on multiple detections");
     221
     222        /*** flag multiple stars */
     223        /* this image star matches more than one tgtcat star */
     224        if (srccat->found_t[N] > -1) {
     225          tgtcat[0].measure[srccat->found_t[N]].dbFlags |= ID_MEAS_BLEND_MEAS;
     226          tgtcat[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_MEAS;
     227        }
     228        if (srccat->found_t[N] == -2) { /* this image star matches a tgtcat star on a neighboring tgtcat */
     229          tgtcat[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_MEAS_X;
     230        }
     231        if (srccat->found_t[N] == -1) { /* this image star matches only this tgtcat star */
     232          srccat->found_t[N] = Nmeas;  /* save first match, in case coincidences are found */
     233        }
     234        /* this tgtcat star matches more than one image star */
     235        if (tgtcat[0].found_t[n] > -1) {
     236          tgtcat[0].measure[tgtcat[0].found_t[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
     237          tgtcat[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
     238        } else {
     239          tgtcat[0].found_t[n] = Nmeas;
     240        }
     241        tgtcat[0].average[n].Nmeasure ++;
     242        Nmeas ++;
    152243      }
    153       if (dRoff < -180.0*3600.0) {
    154           // average on low end of boundary, move star down
    155           catalog[0].measure[Nmeas].R -= 360.0;
    156           dRoff += 360.0*3600.0;
    157       }
    158 
    159       catalog[0].measure[Nmeas].t        = (TIMEREF == 0) ? stars[N][0].measure.t      : TIMEREF; /** careful : time_t vs e_time **/
    160       catalog[0].measure[Nmeas].t_msec   = (TIMEREF == 0) ? stars[N][0].measure.t_msec : TIMEREF; /** careful : time_t vs e_time **/
    161 
    162       if (ACCEPT_MOTION) {
    163         catalog[0].average[n].uR         = stars[N][0].average.uR;
    164         catalog[0].average[n].uD         = stars[N][0].average.uD;
    165         catalog[0].average[n].duR        = stars[N][0].average.duR;
    166         catalog[0].average[n].duD        = stars[N][0].average.duD;
    167         catalog[0].average[n].P          = stars[N][0].average.P;
    168         catalog[0].average[n].dP         = stars[N][0].average.dP;
    169       }
    170 
    171       /** don't update average / secfilt values for REF photcodes **/
    172 
    173       /*** handle multiple stars */
    174       /* this image star matches more than one catalog star */
    175       if (stars[N][0].found > -1) {
    176         catalog[0].measure[stars[N][0].found].dbFlags |= ID_MEAS_BLEND_MEAS;
    177         catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_MEAS;
    178       }
    179       if (stars[N][0].found == -2) { /* this image star matches a catalog star on a neighboring catalog */
    180         catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_MEAS_X;
    181       }
    182       if (stars[N][0].found == -1) { /* this image star matches only this catalog star */
    183         stars[N][0].found = Nmeas;  /* save first match, in case coincidences are found */
    184       }
    185       /* this catalog star matches more than one image star */
    186       if (catalog[0].found_t[n] > -1) {
    187         catalog[0].measure[catalog[0].found_t[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
    188         catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
    189       } else {
    190         catalog[0].found_t[n] = Nmeas;
    191       }
    192 
    193       catalog[0].average[n].Nmeasure ++;
    194       Nmeas ++;
    195 
    196       // update_coords (&catalog[0].average[n], &catalog[0].measure[0], next_meas);
     244     
     245      // update_coords (&tgtcat[0].average[n], &tgtcat[0].measure[0], next_meas);
    197246    }
    198247    i++;
    199248  }
    200 
    201   /* we don't add missed entries for refcat
    202      (already in database, not refcat) */
    203249
    204250  /* incorporate unmatched refcat stars */
    … …  
    207253     the reference up-to-date with known stars only */
    208254
    209   for (i = 0; (i < Nstars) && !options.only_match; i+=NREFSTAR_GROUP) {
    210     if (Nmeas >= NMEAS - NREFSTAR_GROUP) {
    211       NMEAS = Nmeas + 1000;
     255  for (i = 0; (i < Nstars) && !options.only_match; i++) {
     256    /* make sure there is space for next entry */
     257    if (Nmeas >= NMEAS - srccat->average[i].Nmeasure) {
     258      NMEAS = Nmeas + 1000 + srccat->average[i].Nmeasure;
    212259      REALLOCATE (next_meas, off_t, NMEAS);
    213       REALLOCATE (catalog[0].measure, Measure, NMEAS);
    214     }
    215     if (Nave == NAVE) {
     260      REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
     261    }
     262    if (Nave >= NAVE) {
    216263      NAVE = Nave + 1000;
    217       REALLOCATE (catalog[0].average, Average, NAVE);
    218       REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
    219     }
    220 
    221     N = i;
    222     if (stars[N][0].found > -1) continue;
    223     if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue;
    224 
    225     dvo_average_init (&catalog[0].average[Nave]);
    226     catalog[0].average[Nave].R             = stars[N][0].average.R;
    227     catalog[0].average[Nave].D             = stars[N][0].average.D;
    228 
    229     catalog[0].average[Nave].Nmeasure      = NREFSTAR_GROUP;
    230     catalog[0].average[Nave].measureOffset = Nmeas;
    231     catalog[0].average[Nave].objID         = objID;
    232     catalog[0].average[Nave].catID         = catID;
     264      REALLOCATE (tgtcat[0].average, Average, NAVE);
     265      if (tgtcat[0].secfilt) {
     266        // we only update the secfilt table if it has been allocated for output
     267        REALLOCATE (tgtcat[0].secfilt, SecFilt, NAVE*tgtcat[0].Nsecfilt);
     268      }
     269    }
     270
     271    if (srccat->found_t[i] >= 0) continue;
     272    if (!IN_REGION (srccat->average[i].R, srccat->average[i].D)) continue;
     273
     274    dvo_average_init (&tgtcat[0].average[Nave]);
     275    tgtcat[0].average[Nave].R              = srccat->average[i].R;
     276    tgtcat[0].average[Nave].D              = srccat->average[i].D;
     277
     278    tgtcat[0].average[Nave].Nmeasure       = srccat->average[i].Nmeasure;
     279    tgtcat[0].average[Nave].measureOffset  = Nmeas;
     280    tgtcat[0].average[Nave].objID          = objID;
     281    tgtcat[0].average[Nave].catID          = catID;
    233282
    234283    if (PSPS_ID) {
    235         catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D);
     284      tgtcat[0].average[Nave].extID = CreatePSPSObjectID(tgtcat[0].average[Nave].R, tgtcat[0].average[Nave].D);
    236285    }
    237286
    238287    if (ACCEPT_MOTION) {
    239       catalog[0].average[Nave].dR          = stars[N][0].average.dR;
    240       catalog[0].average[Nave].dD          = stars[N][0].average.dD;
    241       catalog[0].average[Nave].uR          = stars[N][0].average.uR;
    242       catalog[0].average[Nave].uD          = stars[N][0].average.uD;
    243       catalog[0].average[Nave].duR         = stars[N][0].average.duR;
    244       catalog[0].average[Nave].duD         = stars[N][0].average.duD;
    245       catalog[0].average[Nave].P           = stars[N][0].average.P;
    246       catalog[0].average[Nave].dP          = stars[N][0].average.dP;
     288      tgtcat[0].average[Nave].dR           = srccat->average[i].dR;
     289      tgtcat[0].average[Nave].dD           = srccat->average[i].dD;
     290      tgtcat[0].average[Nave].uR           = srccat->average[i].uR;
     291      tgtcat[0].average[Nave].uD           = srccat->average[i].uD;
     292      tgtcat[0].average[Nave].duR          = srccat->average[i].duR;
     293      tgtcat[0].average[Nave].duD          = srccat->average[i].duD;
     294      tgtcat[0].average[Nave].P            = srccat->average[i].P;
     295      tgtcat[0].average[Nave].dP           = srccat->average[i].dP;
    247296    }
    248297
    249298    objID ++;
    250299
    251     for (j = 0; j < Nsecfilt; j++) {
    252       dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
    253     }
    254 
    255     for (j = 0; j < NREFSTAR_GROUP; j++) {
    256       N = i + j;
    257 
    258       catalog[0].measure[Nmeas]          = stars[N][0].measure;
    259 
    260       catalog[0].measure[Nmeas].t        = (stars[N][0].measure.t == 0) ? TIMEREF : stars[N][0].measure.t;      /** careful : time_t vs e_time **/
    261       catalog[0].measure[Nmeas].t_msec   = (stars[N][0].measure.t == 0) ?       0 : stars[N][0].measure.t_msec; /** careful : time_t vs e_time **/
    262       catalog[0].measure[Nmeas].dbFlags  = 0;
    263       catalog[0].measure[Nmeas].averef   = Nave;
    264       catalog[0].measure[Nmeas].objID    = catalog[0].average[Nave].objID;
    265       catalog[0].measure[Nmeas].catID    = catalog[0].catID;
    266 
    267       stars[N][0].found = Nmeas;
    268       next_meas[Nmeas] = -1;  // initial value here, update below
    269       Nmeas ++;
    270     }
    271 
    272     for (j = 0; j < NREFSTAR_GROUP - 1; j++) {
    273         next_meas[Nmeas - NREFSTAR_GROUP + j] = Nmeas - NREFSTAR_GROUP + j + 1;
     300    for (j = 0; tgtcat[0].secfilt && (j < Nsecfilt); j++) {
     301      dvo_secfilt_init (&tgtcat[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
     302    }
     303
     304    off_t mOff = srccat->average[i].measureOffset;
     305    for (mSrc = 0; mSrc < srccat->average[i].Nmeasure; mSrc++) {
     306      // supply the measurments from this detection
     307      dvo_measure_init (&tgtcat[0].measure[Nmeas]);
     308      tgtcat[0].measure[Nmeas] = srccat->measure[mSrc + mOff];
     309
     310      // the following measure elements cannot be set until here:
     311      tgtcat[0].measure[Nmeas].dbFlags   = 0;
     312      tgtcat[0].measure[Nmeas].averef    = Nave; // XXX EAM : must be absolute Nave if partial read
     313      tgtcat[0].measure[Nmeas].objID     = tgtcat[0].average[Nave].objID;
     314      tgtcat[0].measure[Nmeas].catID     = tgtcat[0].catID;
     315
     316      if (TIMEREF) {
     317        tgtcat[0].measure[Nmeas].t        = TIMEREF;
     318        tgtcat[0].measure[Nmeas].t_msec   =       0;
     319      }
     320
     321      /* next[Nmeas] should always be -1 in this context (it is always the only
     322         measurement for the star) */
     323      srccat->found_t[i] = Nmeas;
     324      Nmeas ++;
     325
     326      // update the next_meas pointer for this entry (last one for this star is -1)
     327      next_meas[Nmeas-1] = (mSrc < srccat->average[i].Nmeasure - 1) ? Nmeas : -1;
    274328    }
    275329    Nave ++;
    276330  }
    277331     
    278   REALLOCATE (catalog[0].average, Average, Nave);
     332  REALLOCATE (tgtcat[0].average, Average, Nave);
     333  REALLOCATE (tgtcat[0].measure, Measure, Nmeas);
    279334
    280335  // XXX allow for unsorted output?
    281   catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
    282 
    283   /* note stars which have been found in this catalog */
    284   for (i = 0; i < Nstars; i++) {
    285     if (stars[i][0].found > -1) {
    286       stars[i][0].found = -2;
    287     } else {
    288       stars[i][0].found = -3;
    289     }
    290   }
    291 
    292   catalog[0].objID    = objID; // new max value, save on catalog close
    293   catalog[0].Naverage = Nave;
    294   catalog[0].Nmeasure = Nmeas;
    295   catalog[0].Nsecfilt_mem = Nave*Nsecfilt;
    296   if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %d "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
     336  tgtcat[0].measure = sort_measure (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas, next_meas);
     337
     338  /* note stars which have been found in this tgtcat */
     339  for (i = 0; i < srccat->Nmeasure; i++) {
     340    if (srccat->found_t[i] > -1) {
     341      srccat->found_t[i] = -2;
     342    }
     343  }
     344
     345  tgtcat[0].objID    = objID; // new max value, save on tgtcat close
     346  tgtcat[0].Naverage = Nave;
     347  tgtcat[0].Nmeasure = Nmeas;
     348  tgtcat[0].Nsecfilt_mem = Nave*Nsecfilt;
     349  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
    297350
    298351  free (X1);
    … …  
    306359  return (Nmatch);
    307360}
     361
     362/** differences between find_matches_refstars.c and find_matches.c:
     363 * refstars has replace_match
     364 * refstars has TIMEREF and ACCEPT_MOTION
     365 * refstars does NOT have AddToCalibration
     366 * refstars does NOT have Nsec
     367 * refstars does NOT do update_coords
     368 **/

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/get2mass_full.c

    r34088 r38593  
    66// fill in the data for a JHK triplet star.  takes a pointer to the start of the line the
    77// RA and DEC have already been set
    8 int get2mass_3star_full (Stars **star, char *line, int Nmax) {
     8int get2mass_3star_full (Measure *measure, char *line, int Nmax) {
    99
    1010  char *ptr;
    … …  
    1717  ptr = next2MASSfield (ptr); // dec (assumed to be already set)
    1818  ptr = next2MASSfield (ptr); // err_maj
    19   star[0][0].measure.FWx = ToShortPixels(strtod (ptr, NULL));
     19  measure[0].FWx = ToShortPixels(strtod (ptr, NULL));
    2020  ptr = next2MASSfield (ptr); // err_min
    21   star[0][0].measure.FWy = ToShortPixels(strtod (ptr, NULL));
     21  measure[0].FWy = ToShortPixels(strtod (ptr, NULL));
    2222  ptr = next2MASSfield (ptr); // err_ang
    23   star[0][0].measure.theta = ToShortPixels(strtod (ptr, NULL));
    24 
    25   star[2][0].measure.FWx   = star[1][0].measure.FWx   = star[0][0].measure.FWx;
    26   star[2][0].measure.FWy   = star[1][0].measure.FWy   = star[0][0].measure.FWy;
    27   star[2][0].measure.theta = star[1][0].measure.theta = star[0][0].measure.theta;
     23  measure[0].theta = ToShortPixels(strtod (ptr, NULL));
     24
     25  measure[2].FWx   = measure[1].FWx   = measure[0].FWx;
     26  measure[2].FWy   = measure[1].FWy   = measure[0].FWy;
     27  measure[2].theta = measure[1].theta = measure[0].theta;
    2828
    2929  ptr = next2MASSfield (ptr); // designation (skip)
    3030
    3131  ptr = next2MASSfield (ptr); // j_m
    32   star[0][0].measure.M  = strtod (ptr, NULL);
     32  measure[0].M  = strtod (ptr, NULL);
    3333  ptr = next2MASSfield (ptr); // j_cmsig
    34   star[0][0].measure.dM = strtod (ptr, NULL);
     34  measure[0].dM = strtod (ptr, NULL);
    3535  ptr = next2MASSfield (ptr); // j_msigcom
    3636  dMfull = strtod (ptr, NULL);
    37   star[0][0].measure.dMcal = sqrt (SQ(dMfull) - SQ(star[0][0].measure.dM));
     37  measure[0].dMcal = sqrt (SQ(dMfull) - SQ(measure[0].dM));
    3838  ptr = next2MASSfield (ptr); // j_snr (skip)
    3939
    4040  ptr = next2MASSfield (ptr); // h_m
    41   star[1][0].measure.M  = strtod (ptr, NULL);
     41  measure[1].M  = strtod (ptr, NULL);
    4242  ptr = next2MASSfield (ptr); // h_cmsig
    43   star[1][0].measure.dM = strtod (ptr, NULL);
     43  measure[1].dM = strtod (ptr, NULL);
    4444  ptr = next2MASSfield (ptr); // h_msigcom
    4545  dMfull = strtod (ptr, NULL);
    46   star[1][0].measure.dMcal = sqrt (SQ(dMfull) - SQ(star[1][0].measure.dM));
     46  measure[1].dMcal = sqrt (SQ(dMfull) - SQ(measure[1].dM));
    4747  ptr = next2MASSfield (ptr); // h_snr (skip)
    4848
    4949  ptr = next2MASSfield (ptr); // k_m
    50   star[2][0].measure.M  = strtod (ptr, NULL);
     50  measure[2].M  = strtod (ptr, NULL);
    5151  ptr = next2MASSfield (ptr); // k_cmsig
    52   star[2][0].measure.dM = strtod (ptr, NULL);
     52  measure[2].dM = strtod (ptr, NULL);
    5353  ptr = next2MASSfield (ptr); // k_msigcom
    5454  dMfull = strtod (ptr, NULL);
    55   star[2][0].measure.dMcal = sqrt (SQ(dMfull) - SQ(star[2][0].measure.dM));
     55  measure[2].dMcal = sqrt (SQ(dMfull) - SQ(measure[2].dM));
    5656  ptr = next2MASSfield (ptr); // k_snr (skip)
    5757
    58   star[2][0].measure.photFlags = star[1][0].measure.photFlags = star[0][0].measure.photFlags = 0;
     58  measure[2].photFlags = measure[1].photFlags = measure[0].photFlags = 0;
    5959
    6060  ptr = next2MASSfield (ptr); // ph_qual
    61   set2MASS_ph_qual (star[0], ptr[0]);
    62   set2MASS_ph_qual (star[1], ptr[1]);
    63   set2MASS_ph_qual (star[2], ptr[2]);
     61  set2MASS_ph_qual (&measure[0], ptr[0]);
     62  set2MASS_ph_qual (&measure[1], ptr[1]);
     63  set2MASS_ph_qual (&measure[2], ptr[2]);
    6464
    6565  ptr = next2MASSfield (ptr); // rd_flg
    66   set2MASS_rd_flag (star[0], ptr[0]);
    67   set2MASS_rd_flag (star[1], ptr[1]);
    68   set2MASS_rd_flag (star[2], ptr[2]);
     66  set2MASS_rd_flag (&measure[0], ptr[0]);
     67  set2MASS_rd_flag (&measure[1], ptr[1]);
     68  set2MASS_rd_flag (&measure[2], ptr[2]);
    6969
    7070  ptr = next2MASSfield (ptr); // bl_flg
    71   set2MASS_bl_flag (star[0], ptr[0]);
    72   set2MASS_bl_flag (star[1], ptr[1]);
    73   set2MASS_bl_flag (star[2], ptr[2]);
     71  set2MASS_bl_flag (&measure[0], ptr[0]);
     72  set2MASS_bl_flag (&measure[1], ptr[1]);
     73  set2MASS_bl_flag (&measure[2], ptr[2]);
    7474
    7575  ptr = next2MASSfield (ptr); // cc_flg
    76   set2MASS_cc_flag (star[0], ptr[0]);
    77   set2MASS_cc_flag (star[1], ptr[1]);
    78   set2MASS_cc_flag (star[2], ptr[2]);
     76  set2MASS_cc_flag (&measure[0], ptr[0]);
     77  set2MASS_cc_flag (&measure[1], ptr[1]);
     78  set2MASS_cc_flag (&measure[2], ptr[2]);
    7979
    8080  ptr = next2MASSfield (ptr); // ndet (skip for now, XXX use somehow?)
    … …  
    8484
    8585  ptr = next2MASSfield (ptr); // gal_contam (one flag for all filters)
    86   set2MASS_gal_flag (star[0], ptr[0]);
    87   set2MASS_gal_flag (star[1], ptr[0]);
    88   set2MASS_gal_flag (star[2], ptr[0]);
     86  set2MASS_gal_flag (&measure[0], ptr[0]);
     87  set2MASS_gal_flag (&measure[1], ptr[0]);
     88  set2MASS_gal_flag (&measure[2], ptr[0]);
    8989
    9090  ptr = next2MASSfield (ptr); // mp_flg (one flag for all filters)
    91   set2MASS_mp_flag (star[0], ptr[0]);
    92   set2MASS_mp_flag (star[1], ptr[0]);
    93   set2MASS_mp_flag (star[2], ptr[0]);
     91  set2MASS_mp_flag (&measure[0], ptr[0]);
     92  set2MASS_mp_flag (&measure[1], ptr[0]);
     93  set2MASS_mp_flag (&measure[2], ptr[0]);
    9494
    9595  ptr = next2MASSfield (ptr); // pts_key (skip for now, XXX use somehow?)
    … …  
    101101
    102102  ptr = next2MASSfield (ptr); // x_scan
    103   star[0][0].measure.Xccd = strtod (ptr, NULL);
    104   star[2][0].measure.Xccd = star[1][0].measure.Xccd = star[0][0].measure.Xccd;
     103  measure[0].Xccd = strtod (ptr, NULL);
     104  measure[2].Xccd = measure[1].Xccd = measure[0].Xccd;
    105105
    106106  ptr = next2MASSfield (ptr); // jdate (julian date)
    107107  jd = strtod (ptr, NULL);
    108   star[0][0].measure.t = ohana_jd_to_sec (jd);
    109   star[2][0].measure.t = star[1][0].measure.t = star[0][0].measure.t;
     108  measure[0].t = ohana_jd_to_sec (jd);
     109  measure[2].t = measure[1].t = measure[0].t;
    110110
    111111  ptr = next2MASSfield (ptr); // j_psfchi
    112   star[0][0].measure.psfChisq = strtod (ptr, NULL);
     112  measure[0].psfChisq = strtod (ptr, NULL);
    113113  ptr = next2MASSfield (ptr); // h_psfchi
    114   star[1][0].measure.psfChisq = strtod (ptr, NULL);
     114  measure[1].psfChisq = strtod (ptr, NULL);
    115115  ptr = next2MASSfield (ptr); // k_psfchi
    116   star[2][0].measure.psfChisq = strtod (ptr, NULL);
     116  measure[2].psfChisq = strtod (ptr, NULL);
    117117
    118118  ptr = next2MASSfield (ptr); // j_m_stdap
    119   star[0][0].measure.Map = strtod (ptr, NULL);
     119  measure[0].Map = strtod (ptr, NULL);
    120120  ptr = next2MASSfield (ptr); // j_msig_stdap (skip?)
    121121
    122122  ptr = next2MASSfield (ptr); // h_m_stdap
    123   star[1][0].measure.Map = strtod (ptr, NULL);
     123  measure[1].Map = strtod (ptr, NULL);
    124124  ptr = next2MASSfield (ptr); // h_msig_stdap (skip?)
    125125
    126126  ptr = next2MASSfield (ptr); // k_m_stdap
    127   star[2][0].measure.Map = strtod (ptr, NULL);
     127  measure[2].Map = strtod (ptr, NULL);
    128128  ptr = next2MASSfield (ptr); // k_msig_stdap (skip?)
    129129
    … …  
    133133
    134134  ptr = next2MASSfield (ptr); // dup_src (one flag for all filters)
    135   set2MASS_dup_flag (star[0], ptr[0]);
    136   set2MASS_dup_flag (star[1], ptr[0]);
    137   set2MASS_dup_flag (star[2], ptr[0]);
     135  set2MASS_dup_flag (&measure[0], ptr[0]);
     136  set2MASS_dup_flag (&measure[1], ptr[0]);
     137  set2MASS_dup_flag (&measure[2], ptr[0]);
    138138
    139139  ptr = next2MASSfield (ptr); // use_src (one flag for all filters)
    140   set2MASS_use_flag (star[0], ptr[0]);
    141   set2MASS_use_flag (star[1], ptr[0]);
    142   set2MASS_use_flag (star[2], ptr[0]);
    143 
    144   star[0][0].measure.photcode  = TM_J;
    145   star[0][0].measure.detID   = 0;
    146   star[0][0].measure.imageID = 0;
    147 
    148   star[1][0].measure.photcode  = TM_H;
    149   star[1][0].measure.detID   = 0;
    150   star[1][0].measure.imageID = 0;
    151 
    152   star[2][0].measure.photcode  = TM_K;
    153   star[2][0].measure.detID   = 0;
    154   star[2][0].measure.imageID = 0;
     140  set2MASS_use_flag (&measure[0], ptr[0]);
     141  set2MASS_use_flag (&measure[1], ptr[0]);
     142  set2MASS_use_flag (&measure[2], ptr[0]);
     143
     144  measure[0].photcode  = TM_J;
     145  measure[0].detID   = 0;
     146  measure[0].imageID = 0;
     147
     148  measure[1].photcode  = TM_H;
     149  measure[1].detID   = 0;
     150  measure[1].imageID = 0;
     151
     152  measure[2].photcode  = TM_K;
     153  measure[2].detID   = 0;
     154  measure[2].imageID = 0;
    155155
    156156  return TRUE;
    … …  
    170170}
    171171
    172 int set2MASS_ph_qual (Stars *star, char qual) {
    173 
    174   switch (qual) {
    175     case 'A': star[0].measure.photFlags |= 0x00000001; break; // was: 0x0004
    176     case 'B': star[0].measure.photFlags |= 0x00000002; break; // was: 0x0005
    177     case 'C': star[0].measure.photFlags |= 0x00000004; break; // was: 0x0006
    178     case 'D': star[0].measure.photFlags |= 0x00000008; break; // was: 0x0007
    179     case 'E': star[0].measure.photFlags |= 0x00000010; break; // was: 0x0003
    180     case 'F': star[0].measure.photFlags |= 0x00000020; break; // was: 0x0002
    181     case 'U': star[0].measure.photFlags |= 0x00000040; break; // was: 0x0001
    182     case 'X': star[0].measure.photFlags |= 0x00000080; break; // was: 0x0000
     172int set2MASS_ph_qual (Measure *measure, char qual) {
     173
     174  switch (qual) {
     175    case 'A': measure[0].photFlags |= 0x00000001; break; // was: 0x0004
     176    case 'B': measure[0].photFlags |= 0x00000002; break; // was: 0x0005
     177    case 'C': measure[0].photFlags |= 0x00000004; break; // was: 0x0006
     178    case 'D': measure[0].photFlags |= 0x00000008; break; // was: 0x0007
     179    case 'E': measure[0].photFlags |= 0x00000010; break; // was: 0x0003
     180    case 'F': measure[0].photFlags |= 0x00000020; break; // was: 0x0002
     181    case 'U': measure[0].photFlags |= 0x00000040; break; // was: 0x0001
     182    case 'X': measure[0].photFlags |= 0x00000080; break; // was: 0x0000
    183183    default:
    184184      fprintf (stderr, "error!\n");
    … …  
    188188}
    189189
    190 int set2MASS_rd_flag (Stars *star, char qual) {
    191 
    192   switch (qual) {
    193     case '0': star[0].measure.photFlags |= 0x00000100; break; // was: 0x0000
    194     case '1': star[0].measure.photFlags |= 0x00000200; break; // was: 0x0010
    195     case '2': star[0].measure.photFlags |= 0x00000400; break; // was: 0x0020
    196     case '3': star[0].measure.photFlags |= 0x00000800; break; // was: 0x0030
    197     case '4': star[0].measure.photFlags |= 0x00001000; break; // was: 0x0040
    198     case '6': star[0].measure.photFlags |= 0x00002000; break; // was: 0x0050
    199     case '9': star[0].measure.photFlags |= 0x00004000; break; // was: 0x0060
     190int set2MASS_rd_flag (Measure *measure, char qual) {
     191
     192  switch (qual) {
     193    case '0': measure[0].photFlags |= 0x00000100; break; // was: 0x0000
     194    case '1': measure[0].photFlags |= 0x00000200; break; // was: 0x0010
     195    case '2': measure[0].photFlags |= 0x00000400; break; // was: 0x0020
     196    case '3': measure[0].photFlags |= 0x00000800; break; // was: 0x0030
     197    case '4': measure[0].photFlags |= 0x00001000; break; // was: 0x0040
     198    case '6': measure[0].photFlags |= 0x00002000; break; // was: 0x0050
     199    case '9': measure[0].photFlags |= 0x00004000; break; // was: 0x0060
    200200    default:
    201201      fprintf (stderr, "error!\n");
    … …  
    205205}
    206206
    207 int set2MASS_cc_flag (Stars *star, char qual) {
    208 
    209   switch (qual) {
    210     case 'p': star[0].measure.photFlags |= 0x00010000; break; // was: 0x0000
    211     case 'c': star[0].measure.photFlags |= 0x00020000; break; // was: 0x0100
    212     case 'd': star[0].measure.photFlags |= 0x00040000; break; // was: 0x0200
    213     case 's': star[0].measure.photFlags |= 0x00080000; break; // was: 0x0300
    214     case 'b': star[0].measure.photFlags |= 0x00010000; break; // was: 0x0400
    215     case '0': star[0].measure.photFlags |= 0x00020000; break; // was: 0x0500
     207int set2MASS_cc_flag (Measure *measure, char qual) {
     208
     209  switch (qual) {
     210    case 'p': measure[0].photFlags |= 0x00010000; break; // was: 0x0000
     211    case 'c': measure[0].photFlags |= 0x00020000; break; // was: 0x0100
     212    case 'd': measure[0].photFlags |= 0x00040000; break; // was: 0x0200
     213    case 's': measure[0].photFlags |= 0x00080000; break; // was: 0x0300
     214    case 'b': measure[0].photFlags |= 0x00010000; break; // was: 0x0400
     215    case '0': measure[0].photFlags |= 0x00020000; break; // was: 0x0500
    216216    default:
    217217      fprintf (stderr, "error!\n");
    … …  
    221221}
    222222
    223 int set2MASS_bl_flag (Stars *star, char qual) {
    224 
    225   switch (qual) {
    226     case '0': star[0].measure.photFlags &= ~0x00300000; break; // was: ~0x0008
    227     case '1': star[0].measure.photFlags |=  0x00100000; break; // was: ~0x0008
    228     default:  star[0].measure.photFlags |=  0x00200000; break; // was:  0x0008
    229   }     
    230   return (TRUE);
    231 }
    232 
    233 int set2MASS_gal_flag (Stars *star, char qual) {
    234 
    235   switch (qual) {
    236     case '0': star[0].measure.photFlags &= ~0x00c00000; break; // was: ~0x0080
    237     case '1': star[0].measure.photFlags |=  0x00400000; break; // was: ~0x0080
    238     default:  star[0].measure.photFlags |=  0x00800000;        // was:  0x0080
    239       star[0].measure.extNsigma = 100.0;
     223int set2MASS_bl_flag (Measure *measure, char qual) {
     224
     225  switch (qual) {
     226    case '0': measure[0].photFlags &= ~0x00300000; break; // was: ~0x0008
     227    case '1': measure[0].photFlags |=  0x00100000; break; // was: ~0x0008
     228    default:  measure[0].photFlags |=  0x00200000; break; // was:  0x0008
     229  }     
     230  return (TRUE);
     231}
     232
     233int set2MASS_gal_flag (Measure *measure, char qual) {
     234
     235  switch (qual) {
     236    case '0': measure[0].photFlags &= ~0x00c00000; break; // was: ~0x0080
     237    case '1': measure[0].photFlags |=  0x00400000; break; // was: ~0x0080
     238    default:  measure[0].photFlags |=  0x00800000;             // was:  0x0080
     239      measure[0].extNsigma = 100.0;
    240240      break;
    241241  }     
    … …  
    243243}
    244244
    245 int set2MASS_mp_flag (Stars *star, char qual) {
    246 
    247   switch (qual) {
    248     case '0': star[0].measure.photFlags &= ~0x03000000; break; // was: ~0x0800
    249     case '1': star[0].measure.photFlags |=  0x01000000; break; // was: ~0x0800
    250     default:  star[0].measure.photFlags |=  0x02000000; break; // was:  0x0800
    251   }     
    252   return (TRUE);
    253 }
    254 
    255 int set2MASS_dup_flag (Stars *star, char qual) {
    256 
    257   switch (qual) {
    258     case '0': star[0].measure.photFlags &= ~0x0c000000; break; // was: ~0x1000
    259     case '1': star[0].measure.photFlags |=  0x04000000; break; // was: ~0x1000
    260     default:  star[0].measure.photFlags |=  0x08000000; break; // was:  0x1000
    261   }     
    262   return (TRUE);
    263 }
    264 
    265 int set2MASS_use_flag (Stars *star, char qual) {
    266 
    267   switch (qual) {
    268     case '0': star[0].measure.photFlags &= ~0x10000000; break; // was: ~0x2000
    269     case '1': star[0].measure.photFlags |=  0x10000000; break; // was:  0x2000
     245int set2MASS_mp_flag (Measure *measure, char qual) {
     246
     247  switch (qual) {
     248    case '0': measure[0].photFlags &= ~0x03000000; break; // was: ~0x0800
     249    case '1': measure[0].photFlags |=  0x01000000; break; // was: ~0x0800
     250    default:  measure[0].photFlags |=  0x02000000; break; // was:  0x0800
     251  }     
     252  return (TRUE);
     253}
     254
     255int set2MASS_dup_flag (Measure *measure, char qual) {
     256
     257  switch (qual) {
     258    case '0': measure[0].photFlags &= ~0x0c000000; break; // was: ~0x1000
     259    case '1': measure[0].photFlags |=  0x04000000; break; // was: ~0x1000
     260    default:  measure[0].photFlags |=  0x08000000; break; // was:  0x1000
     261  }     
     262  return (TRUE);
     263}
     264
     265int set2MASS_use_flag (Measure *measure, char qual) {
     266
     267  switch (qual) {
     268    case '0': measure[0].photFlags &= ~0x10000000; break; // was: ~0x2000
     269    case '1': measure[0].photFlags |=  0x10000000; break; // was:  0x2000
    270270    default:  abort();
    271271  }     

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/get2mass_ops.c

    r34088 r38593  
    3636}
    3737
     38# if (0)
    3839// fill in the data for a single star.  takes a pointer to the start of the line
    3940int get2mass_star (Stars *star, char *line, int Nmax) {
    … …  
    136137  return TRUE;
    137138}
     139# endif
    138140
    139141// fill in the coords for a single star.  takes a pointer to the start of the line

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/getgsc.c

    r38467 r38593  
    44# define DNSTARS 1000
    55
    6 static short GSC_M;
    76static SkyTable *sky = NULL;
    87
    9 Stars *getgsc (SkyRegion *patch, unsigned int *NSTARS) {
     8static short int GSC_M;
     9
     10Catalog *getgsc (SkyRegion *patch) {
    1011 
    11   unsigned int i, Ngsc, Nstars;
    12   Stars *gsc;
    13   Stars *stars;
     12  unsigned int i;
    1413  SkyList *skylist;
    1514
    … …  
    2322  skylist = SkyListByPatch (sky, -1, patch);
    2423 
    25   Nstars = 0;
    26   ALLOCATE (stars, Stars, 1);
    27  
     24  Catalog *catalog = NULL;
     25  ALLOCATE (catalog, Catalog, 1);
     26  dvo_catalog_init (catalog, TRUE);
     27  ALLOCATE (catalog->average, Average, 1);
     28  ALLOCATE (catalog->measure, Measure, 1);
     29  ALLOCATE (catalog->lensing, Lensing, 1);
     30
     31  SkyRegion region;
     32  region.Rmin =   0.0;
     33  region.Rmax = 360.0;
     34  region.Dmin = -90.0;
     35  region.Dmax = +90.0;
     36
    2837  for (i = 0; i < skylist[0].Nregions; i++) {
    29     gsc = rd_gsc (skylist[0].filename[i], &Ngsc);
     38    Catalog *newcat = rd_gsc (skylist[0].filename[i]);
    3039
    31     REALLOCATE (stars, Stars, MAX (1, Nstars + Ngsc));
    32     memcpy (&stars[Nstars], gsc, Ngsc*sizeof(Stars));
    33     Nstars += Ngsc;
    34 
    35     free (gsc);
     40    AddstarClientOptions matchOptions;
     41    matchOptions.radius = 0.4; // tight radius at this stage
     42    matchOptions.calibrate = FALSE;
     43    matchOptions.only_match = FALSE;
     44    matchOptions.nosort = FALSE;
     45    matchOptions.photcode = 0; // use an invalid photcode to avoid touching secfilt
     46   
     47    find_matches_closest (&region, newcat, catalog, matchOptions);
     48   
     49    if (newcat) {
     50      dvo_catalog_free (newcat);
     51      free (newcat);
     52    }
    3653  }
    3754 
    38   // XXX don't free because we only load it once and save it
    3955  SkyListFree (skylist);
    40   // SkyTableFree (sky);
    4156
    42   if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", Nstars);
    43   *NSTARS = Nstars;
    44   return (stars);
     57  if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", (int) catalog->Naverage);
     58  return (catalog);
    4559} 
    4660
    47 Stars *rd_gsc (char *filename, unsigned int *nstars) {
     61Catalog *rd_gsc (char *filename) {
    4862 
    49   int i, NSTARS, Nstars, Nbytes, nbytes, Nline, Nbyte;
     63  int i, Nbytes, nbytes, Nline, Nbyte;
    5064  char *buffer;
    5165  FILE *f;
    52   Stars *stars;
    5366
    54   Nstars = 0;
    55   NSTARS = 1000;
    56   ALLOCATE (stars, Stars, NSTARS);
     67  int Nave = 0;
     68  int Nmeas = 0;
     69  int NAVE = 1000;
     70  int NMEAS = 1000;
     71
     72  Catalog *catalog = NULL;
     73  ALLOCATE (catalog, Catalog, 1);
     74  dvo_catalog_init (catalog, TRUE);
     75  ALLOCATE (catalog->average, Average, NAVE);
     76  ALLOCATE (catalog->measure, Measure, NMEAS);
    5777
    5878  f = fopen (filename, "r");
    … …  
    6888    for (i = 0; i < Nline; i++) {
    6989
    70       InitStar (&stars[Nstars]);
     90      dvo_average_init (&catalog->average[Nave]);
     91      dvo_measure_init (&catalog->measure[Nmeas]);
    7192
    7293      Nbyte = i*BYTES_STAR;
    73       dparse (&stars[Nstars].average.R, 1, &buffer[Nbyte]);
    74       dparse (&stars[Nstars].average.D, 2, &buffer[Nbyte]);
     94      dparse (&catalog->average[Nave].R, 1, &buffer[Nbyte]);
     95      dparse (&catalog->average[Nave].D, 2, &buffer[Nbyte]);
    7596
    76       if (stars[Nstars].average.R < UserPatch.Rmin) continue;
    77       if (stars[Nstars].average.R > UserPatch.Rmax) continue;
    78       if (stars[Nstars].average.D < UserPatch.Dmin) continue;
    79       if (stars[Nstars].average.D > UserPatch.Dmax) continue;
     97      if (catalog->average[Nave].R < UserPatch.Rmin) continue;
     98      if (catalog->average[Nave].R > UserPatch.Rmax) continue;
     99      if (catalog->average[Nave].D < UserPatch.Dmin) continue;
     100      if (catalog->average[Nave].D > UserPatch.Dmax) continue;
    80101
    81       stars[Nstars].measure.R = stars[Nstars].average.R;
    82       stars[Nstars].measure.D = stars[Nstars].average.D;
     102      catalog->measure[Nmeas].R = catalog->average[Nave].R;
     103      catalog->measure[Nmeas].D = catalog->average[Nave].D;
    83104
    84       fparse (&stars[Nstars].measure.M, 3, &buffer[Nbyte]);
    85       stars[Nstars].measure.dM        = NAN;
    86       stars[Nstars].measure.photcode = GSC_M;
    87       stars[Nstars].measure.t        = 0;
     105      fparse (&catalog->measure[Nmeas].M, 3, &buffer[Nbyte]);
     106      catalog->measure[Nmeas].dM = NAN;
     107      catalog->measure[Nmeas].photcode = GSC_M;
     108      catalog->measure[Nmeas].t = 0;
    88109
    89       Nstars ++;
    90       CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars, 1000);
     110      catalog->average[Nave].Nmeasure = 1;
     111      catalog->average[Nave].measureOffset = Nmeas;
     112
     113      Nave ++;
     114      Nmeas ++;
     115
     116      CHECK_REALLOCATE (catalog->average, Average, NAVE,  Nave,  1000);
     117      CHECK_REALLOCATE (catalog->measure, Measure, NMEAS, Nmeas, 1000);
    91118    }
    92119  }
    93120  fclose (f);
    94121  free (buffer);
    95   REALLOCATE (stars, Stars, Nstars);
    96   *nstars = Nstars;
    97   return (stars);
     122  return (catalog);
    98123}

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/getusno.c

    r38467 r38593  
    44# define NELEM  3
    55
    6 Stars *getusno (SkyRegion *catstats, int photcode, unsigned int *nstars) {
     6Catalog *getusno (SkyRegion *catstats) {
    77
    88  off_t offset;
    … …  
    1515  double dec;
    1616  int spd, spd_start, spd_end;
    17   int NUSNO, Nusno;
    18   Stars *stars;
    1917  short int USNO_RED, USNO_BLUE;
    2018
    … …  
    2220  NAMED_PHOTCODE (USNO_RED, "USNO_RED");
    2321  NAMED_PHOTCODE (USNO_BLUE, "USNO_BLUE");
    24   if (photcode == USNO_RED) goto good_code;
    25   if (photcode == USNO_BLUE) goto good_code;
    26   Shutdown ("USNO photcode not specified");
    27 good_code:
    2822
    2923  /* identify ra & dec range of interest */
    … …  
    4337  }
    4438
    45   Nusno = 0;
    46   NUSNO = 5000;
    47   ALLOCATE (stars, Stars, NUSNO);
     39  int Nave = 0;
     40  int Nmeas = 0;
     41  int NAVE = 1000;
     42  int NMEAS = 1000;
     43
     44  Catalog *catalog = NULL;
     45  ALLOCATE (catalog, Catalog, 1);
     46  dvo_catalog_init (catalog, TRUE);
     47  ALLOCATE (catalog->average, Average, NAVE);
     48  ALLOCATE (catalog->measure, Measure, NMEAS);
    4849
    4950  for (spd = spd_start; spd < spd_end; spd += 75) {
    … …  
    110111      if (buf[1] > iDEC1) continue;
    111112
    112       InitStar (&stars[Nusno]);
     113      dvo_average_init (&catalog->average[Nave]);
     114      dvo_measure_init (&catalog->measure[Nmeas+0]);
     115      dvo_measure_init (&catalog->measure[Nmeas+1]);
    113116
    114       stars[Nusno].average.R     = buf[0]/360000.0;
    115       stars[Nusno].average.D     = buf[1]/360000.0 - 90.0;
    116       stars[Nusno].measure.R     = stars[Nusno].average.R;
    117       stars[Nusno].measure.D     = stars[Nusno].average.D;
    118       stars[Nusno].measure.dM    = NAN;
     117      catalog->average[Nave].R     = buf[0]/360000.0;
     118      catalog->average[Nave].D     = buf[1]/360000.0 - 90.0;
    119119
    120       /* one pass of addstar does either r or b */
    121       if (photcode == USNO_RED) {
    122         stars[Nusno].measure.photcode = USNO_RED;
    123         stars[Nusno].measure.M        = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000))));
    124       }
    125       if (photcode == USNO_BLUE) {     
    126         stars[Nusno].measure.photcode = USNO_BLUE;
    127         stars[Nusno].measure.M        = fabs (0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000));
    128       }
    129       Nusno ++;
    130       CHECK_REALLOCATE (stars, Stars, NUSNO, Nusno, 5000);
     120      catalog->measure[Nmeas+0].R  = catalog->average[Nave].R;
     121      catalog->measure[Nmeas+0].D  = catalog->average[Nave].D;
     122      catalog->measure[Nmeas+0].dM = NAN;
     123      catalog->measure[Nmeas+1].R  = catalog->average[Nave].R;
     124      catalog->measure[Nmeas+1].D  = catalog->average[Nave].D;
     125      catalog->measure[Nmeas+1].dM = NAN;
     126
     127      catalog->measure[Nmeas+0].photcode = USNO_RED;
     128      catalog->measure[Nmeas+0].M        = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000))));
     129      catalog->measure[Nmeas+1].photcode = USNO_BLUE;
     130      catalog->measure[Nmeas+1].M        = fabs (0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000));
     131
     132      catalog->average[Nave].Nmeasure = 2;
     133      catalog->average[Nave].measureOffset = Nmeas;
     134
     135      Nave ++;
     136      Nmeas += 2;
     137
     138      CHECK_REALLOCATE (catalog->average, Average, NAVE,  Nave,  1000);
     139      CHECK_REALLOCATE (catalog->measure, Measure, NMEAS, Nmeas, 1000);
    131140    }
    132141    free (buffer);
    133142    fclose (f);
    134143  }
    135   REALLOCATE (stars, Stars, Nusno);
     144  catalog->Naverage = Nave;
     145  catalog->Nmeasure = Nmeas;
    136146
    137   *nstars = Nusno;
    138   if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nusno);
    139   return (stars);
     147  if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nave);
     148  return (catalog);
    140149}
    141150

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/getusnob.c

    r38467 r38593  
    33# define NELEM  20
    44
    5 Stars *getusnob (SkyRegion *catstats, int photcode, unsigned int *nstars) {
     5Catalog *getusnob (SkyRegion *catstats) {
    66
    77  off_t offset;
    8   int i, bin, first, last, nitems, Nitems, Nbins, Nitemsum;
     8  int i, j, bin, first, last, nitems, Nitems, Nbins, Nitemsum;
    99  float hours[100];
    1010  int start[100], number[100], *buffer, *buf;
    … …  
    1515  int iDEC0, iDEC1, iRA0, iRA1;
    1616  int spd, spd_start, spd_end;
    17   int NUSNO, Nusno, Nstars;
    1817  short int USNO_RED, USNO_BLUE, USNO_N;
    1918  e_time USNOepoch;
    20   Stars *stars;
    2119
    2220  m1 = m2 = 0.0;
    … …  
    2624  NAMED_PHOTCODE (USNO_BLUE, "USNO_BLUE");
    2725  NAMED_PHOTCODE (USNO_N,    "USNO_N");
    28   if (photcode == USNO_RED)  goto good_code;
    29   if (photcode == USNO_BLUE) goto good_code;
    30   if (photcode == USNO_N)    goto good_code;
    31   Shutdown ("USNO photcode not specified");
    32 good_code:
    3326
    3427  fprintf (stderr, "loading USNO-B 1.0\n");
    … …  
    5245
    5346  Nitemsum = 0;
    54   Nusno = 0;
    55   NUSNO = 5000;
    56   ALLOCATE (stars, Stars, NUSNO);
     47
     48  int Nave = 0;
     49  int Nmeas = 0;
     50  int NAVE = 1000;
     51  int NMEAS = 1000;
     52
     53  Catalog *catalog = NULL;
     54  ALLOCATE (catalog, Catalog, 1);
     55  dvo_catalog_init (catalog, TRUE);
     56  ALLOCATE (catalog->average, Average, NAVE);
     57  ALLOCATE (catalog->measure, Measure, NMEAS);
    5758
    5859  for (spd = spd_start; spd < spd_end; spd ++) {
    … …  
    9697
    9798    /* sum the number of stars in data segment of interest */
    98     Nstars = 0;
     99    int Nstars = 0;
    99100    for (bin = first; bin < last; bin++) {
    100101      Nstars += number[bin];
    … …  
    125126      /* USNO coords are reported for J2000 / epoch 2000.0 */
    126127      /* extract the basic stellar data */
    127       InitStar (&stars[Nusno]);
     128      dvo_average_init (&catalog->average[Nave]);
     129      dvo_measure_init (&catalog->measure[Nmeas+0]);
     130      dvo_measure_init (&catalog->measure[Nmeas+1]);
    128131
    129       stars[Nusno].average.R   = buf[0]/360000.0;
    130       stars[Nusno].average.D   = buf[1]/360000.0 - 90.0;
     132      catalog->average[Nave].R   = buf[0]/360000.0;
     133      catalog->average[Nave].D   = buf[1]/360000.0 - 90.0;
    131134
    132135      /* XXX uR cos(D) or just uR ??? */
    133       stars[Nusno].average.uR  = 2.0 * ((buf[2]       % 10000) - 5000);
    134       stars[Nusno].average.uD  = 2.0 * ((buf[2]/10000 % 10000) - 5000);
     136      catalog->average[Nave].uR  = 2.0 * ((buf[2]       % 10000) - 5000);
     137      catalog->average[Nave].uD  = 2.0 * ((buf[2]/10000 % 10000) - 5000);
    135138
    136       stars[Nusno].average.duR = (buf[3]      % 1000);
    137       stars[Nusno].average.duD = (buf[3]/1000 % 1000);
    138 
    139       stars[Nusno].average.dR  = 0.001 * (buf[4]      % 1000);
    140       stars[Nusno].average.dD  = 0.001 * (buf[4]/1000 % 1000);
    141 
    142       stars[Nusno].average.P   = 0;
    143       stars[Nusno].average.dP  = 0;
    144 
    145       stars[Nusno].measure.R   = stars[Nusno].average.R;
    146       stars[Nusno].measure.D   = stars[Nusno].average.D;
    147 
    148       /* USNO magnitude errors are reported as a fixed 0.3 mag */
    149       stars[Nusno].measure.dM  = 0.3;
     139      catalog->average[Nave].duR = (buf[3]      % 1000);
     140      catalog->average[Nave].duD = (buf[3]/1000 % 1000);
     141      catalog->average[Nave].dR  = 0.001 * (buf[4]      % 1000);
     142      catalog->average[Nave].dD  = 0.001 * (buf[4]/1000 % 1000);
     143      catalog->average[Nave].P   = 0;
     144      catalog->average[Nave].dP  = 0;
    150145
    151146      /* USNO-B uses J2000 equinox and 2000.0 epoch for coordinates */
    152147      /* the magnitudes have no temporal information */
    153       stars[Nusno].measure.t   = USNOepoch;
    154148
    155       /* one pass of addstar does r, b, or n */
    156       if (photcode == USNO_BLUE) {
    157           m1 = fabs(0.01 * (buf[5] % 10000)); /* 1st blue mag */
    158           m2 = fabs(0.01 * (buf[7] % 10000)); /* 1st blue mag */
     149      for (j = 0; j < 3; j++) {
     150        catalog->measure[Nmeas+j].R  = catalog->average[Nave].R;
     151        catalog->measure[Nmeas+j].D  = catalog->average[Nave].D;
     152        catalog->measure[Nmeas+j].dM = 0.3; /* USNO magnitude errors are reported as a fixed 0.3 mag */
     153        catalog->measure[Nmeas+j].t  = USNOepoch;
    159154      }
    160       if (photcode == USNO_RED) {
    161           m1 = fabs(0.01 * (buf[6] % 10000)); /* 1st red mag */
    162           m2 = fabs(0.01 * (buf[8] % 10000)); /* 2nd red mag */
    163       }
    164       if (photcode == USNO_N) {
    165           m1 = fabs(0.01 * (buf[9] % 10000)); /* N (IR) survey */
    166           m2 = 0;
    167       }
    168       stars[Nusno].measure.photcode = photcode;
    169155
    170       /* if two mags are available, get an average */
     156      catalog->measure[Nmeas+0].photcode = USNO_BLUE;
     157      m1 = fabs(0.01 * (buf[5] % 10000)); /* 1st blue mag */
     158      m2 = fabs(0.01 * (buf[7] % 10000)); /* 2nd blue mag */
    171159      if (m1 && m2) {
    172160        mag = 0.5*(m1 + m2);
    … …  
    174162        mag = (m1) ? m1 : m2;
    175163      }
    176       stars[Nusno].measure.M = (mag == 0.0) ? 32.0 : mag;
     164      catalog->measure[Nmeas+0].M = (mag == 0.0) ? 32.0 : mag;
    177165
    178       Nusno ++;
    179       CHECK_REALLOCATE (stars, Stars, NUSNO, Nusno, 5000);
     166      catalog->measure[Nmeas+1].photcode = USNO_RED;
     167      m1 = fabs(0.01 * (buf[6] % 10000)); /* 1st red mag */
     168      m2 = fabs(0.01 * (buf[8] % 10000)); /* 2nd red mag */
     169      if (m1 && m2) {
     170        mag = 0.5*(m1 + m2);
     171      } else {
     172        mag = (m1) ? m1 : m2;
     173      }
     174      catalog->measure[Nmeas+1].M = (mag == 0.0) ? 32.0 : mag;
     175     
     176      catalog->measure[Nmeas+2].photcode = USNO_N;
     177      catalog->measure[Nmeas+2].M = fabs(0.01 * (buf[9] % 10000)); /* N (IR) survey */
     178
     179      catalog->average[Nave].Nmeasure = 3;
     180      catalog->average[Nave].measureOffset = Nmeas;
     181
     182      Nave ++;
     183      Nmeas += 3;
     184
     185      CHECK_REALLOCATE (catalog->average, Average, NAVE,  Nave,  1000);
     186      CHECK_REALLOCATE (catalog->measure, Measure, NMEAS, Nmeas, 1000);
    180187    }
    181188    free (buffer);
    182189    Nitemsum += Nitems;
    183     if (VERBOSE) fprintf (stderr, "Nusno: %d, Nitems: %d, Nitemsum : %d, Nbytes stars: %d\n", Nusno, Nitems, Nitemsum, (unsigned int) (Nusno*sizeof(Stars)));
     190    if (VERBOSE) fprintf (stderr, "Naverage: %d, Nitems: %d, Nitemsum : %d\n", Nave, Nitems, Nitemsum);
    184191    fclose (f);
    185192  }
     193  catalog->Naverage = Nave;
     194  catalog->Nmeasure = Nmeas;
    186195
    187   REALLOCATE (stars, Stars, MAX (1, Nusno));
    188 
    189   *nstars = Nusno;
    190   if (VERBOSE) fprintf (stderr, "%d stars from USNO-B 1.0\n", Nusno);
    191   return (stars);
     196  if (VERBOSE) fprintf (stderr, "%d stars from USNO-B 1.0\n", Nave);
     197  return (catalog);
    192198}

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/greference.c

    r38467 r38593  
    33# define LOAD_DR2    1
    44
    5 Stars *greference (char *Refcat, SkyRegion *region, int photcode, unsigned int *nstars) {
     5Catalog *greference (char *Refcat, SkyRegion *region, int photcode) {
    66
    7   unsigned int Nstars;
    8   Stars *stars;
     7  Catalog *catalog;
    98
    109  if (VERBOSE) fprintf (stderr, "loading reference catalog data from %s\n", Refcat);
    1110  if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n", region[0].Rmin, region[0].Rmax, region[0].Dmin, region[0].Dmax);
    1211
    13   Nstars = 0;
    14   stars = NULL;
    15 
     12  /* get stars from HST GSC for the given region */
     13  if (!strcasecmp (Refcat, "GSC")) {
     14    catalog = getgsc (region);
     15  }
     16 
    1617  /* get stars from USNO for the given region */
    1718  if (!strcasecmp (Refcat, "USNO")) {
    18     stars = getusno (region, photcode, &Nstars);
     19    catalog = getusno (region);
    1920  }
    2021
    2122  /* get stars from the USNO B catalog for the given region */
    2223  if (!strcasecmp (Refcat, "USNOB")) {
    23     stars = getusnob (region, photcode, &Nstars);
     24    catalog = getusnob (region);
    2425  }
    2526
    2627  /* get stars from the Tycho catalog for the given region */
    2728  if (!strcasecmp (Refcat, "TYCHO")) {
    28     // stars = gettycho2 (region, photcode, &Nstars);
    2929    fprintf (stderr, "Tycho load via addstar is deprecated: use loadtycho\n");
    3030    exit (2);
    … …  
    3333  /* get stars from the Tycho catalog for the given region (old ingest) */
    3434  if (!strcasecmp (Refcat, "TYCHO_OLD")) {
    35     // stars = gettycho (region, photcode, &Nstars);
    3635    fprintf (stderr, "Tycho load via addstar is deprecated: use loadtycho\n");
    3736    exit (2);
    3837  }
    3938
    40   /* get stars from HST GSC for the given region */
    41   if (!strcasecmp (Refcat, "GSC")) {
    42     stars = getgsc (region, &Nstars);
    43   }
    44  
    4539  /* get stars from 2MASS for the given region */
    4640  if (!strcasecmp (Refcat, "2MASS")) {
    47     // stars = get2mass (region, photcode, LOAD_ALLSKY, &Nstars);
    4841    fprintf (stderr, "2MASS load via addstar is deprecated: use load2mass\n");
    4942    exit (2);
    … …  
    5245  /* get stars from 2MASS for the given region */
    5346  if (!strcasecmp (Refcat, "2MASS-ALLSKY")) {
    54     // stars = get2mass (region, photcode, LOAD_ALLSKY, &Nstars);
    5547    fprintf (stderr, "2MASS load via addstar is deprecated: use load2mass\n");
    5648    exit (2);
    … …  
    5951  /* get stars from 2MASS for the given region */
    6052  if (!strcasecmp (Refcat, "2MASS-DR2")) {
    61     // stars = get2mass (region, photcode, LOAD_DR2, &Nstars);
    6253    fprintf (stderr, "2MASS load via addstar is deprecated: use load2mass\n");
    6354    exit (2);
    6455  }
    6556 
    66   if (VERBOSE && (Nstars == 0)) fprintf (stderr, "no ref objects in region %s\n", region[0].name);
     57  if (VERBOSE && (catalog->Naverage == 0)) fprintf (stderr, "no ref objects in region %s\n", region[0].name);
    6758
    68   *nstars = Nstars;
    69   return (stars);
    70 
     59  return (catalog);
    7160}

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/grefstars.c

    r38467 r38593  
    22
    33/* read ASCII file with ref star data */
    4 Stars *grefstars (char *file, int photcode, unsigned int *Nstars) {
     4Catalog *grefstars (char *file, int photcode) {
    55
    66  FILE *f;
    7   int N, NSTARS;
    8   Stars *stars;
     7  int NSTARS;
    98  char line[256];
    109
    … …  
    1312  if (f == NULL) Shutdown ("can't read data from %s", file);
    1413
    15   NSTARS = 100;
    16   ALLOCATE (stars, Stars, NSTARS);
     14  NSTARS = 1000;
     15
     16  Catalog *catalog = NULL;
     17  ALLOCATE (catalog, Catalog, 1);
     18  dvo_catalog_init (catalog, TRUE);
     19  ALLOCATE (catalog->average, Average, NSTARS);
     20  ALLOCATE (catalog->measure, Measure, NSTARS);
    1721
    1822  /* read in stars line-by-line */
    19   for (N = 0; scan_line (f, line) != EOF; N++) {
     23  int N = 0;
     24  while (scan_line (f, line) != EOF) {
    2025    stripwhite (line);
    2126    if (line[0] == 0) continue;
    2227    if (line[0] == '#') continue;
    2328
    24     InitStar (&stars[N]);
     29    dvo_average_init (&catalog->average[N]);
     30    dvo_measure_init (&catalog->measure[N]);
    2531
    2632    int status = sscanf (line, "%lf %lf %f %f %x",
    27                      &stars[N].average.R,  &stars[N].average.D, 
    28                      &stars[N].measure.M, &stars[N].measure.dM,
    29                      &stars[N].measure.photFlags);
     33                         &catalog->average[N].R, &catalog->average[N].D, 
     34                         &catalog->measure[N].M, &catalog->measure[N].dM,
     35                         &catalog->measure[N].photFlags);
     36
    3037    if ((status != 4) && (status != 5)) {
    3138      fprintf (stderr, "error reading line: %s\n", line);
    32       N --;
    3339      continue;
    3440    }
    35     stars[N].average.R = ohana_normalize_angle (stars[N].average.R);
     41    catalog->average[N].R = ohana_normalize_angle (catalog->average[N].R);
    3642
    37     stars[N].measure.R = stars[N].average.R;
    38     stars[N].measure.D = stars[N].average.D;
     43    catalog->measure[N].R = catalog->average[N].R;
     44    catalog->measure[N].D = catalog->average[N].D;
    3945
    40     stars[N].measure.photcode = photcode;
     46    catalog->measure[N].photcode = photcode;
    4147
    42     CHECK_REALLOCATE (stars, Stars, NSTARS, N+1, 100);
     48    catalog->average[N].Nmeasure = 1;
     49    catalog->average[N].measureOffset = N;
     50
     51    N++;
     52    CHECK_REALLOCATE (catalog->average, Average, NSTARS, N, 1000);
     53    CHECK_REALLOCATE (catalog->measure, Measure, NSTARS, N, 1000);
    4354  }
    44   *Nstars = N;
    45   return (stars);
     55  catalog->Naverage = N;
     56  catalog->Nmeasure = N;
     57  return (catalog);
    4658}

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/load2mass_as_rawdata.c

    r38471 r38593  
    1313 
    1414  int i, j, verbose;
    15   int Nstars, NSTARS, Ntstars, NTSTARS;
     15  int Ntstars, NTSTARS;
    1616  int Nbyte, Nextra, offset;
    1717
    … …  
    2121  char *buffer, *p, *q;
    2222
    23   Stars **stars; // this is an array of pointers to be consistent with input to find_match_refstars
    2423  TMStars *tstars;
    2524  SkyList *skylist;
    … …  
    117116      if (DEBUG) fprintf (stderr, "writing to %s\n", skylist[0].filename[0]);
    118117
    119       // collect array of (Stars *) stars in a new output catalog
    120       Nstars = 0;
    121       NSTARS = 3000;
    122       ALLOCATE (stars, Stars *, NSTARS);
     118      int Nave = 0;
     119      int Nmeas = 0;
     120      int NAVE = 1000;
     121      int NMEAS = 3000;
     122
     123      Catalog *newcat = NULL;
     124      ALLOCATE (newcat, Catalog, 1);
     125      dvo_catalog_init (newcat, TRUE);
     126      ALLOCATE (newcat->average, Average, NAVE);
     127      ALLOCATE (newcat->measure, Measure, NMEAS);
    123128
    124129      // loop over stars in this 2mass region that are also in this output region
    … …  
    140145        offset = tstars[j].offset;
    141146
    142         ALLOCATE (stars[Nstars+0], Stars, 1);
    143         ALLOCATE (stars[Nstars+1], Stars, 1);
    144         ALLOCATE (stars[Nstars+2], Stars, 1);
    145 
    146         InitStar (stars[Nstars+0]);
    147         InitStar (stars[Nstars+1]);
    148         InitStar (stars[Nstars+2]);
    149 
    150         stars[Nstars+0][0].average.R = tstars[j].R;
    151         stars[Nstars+0][0].average.D = tstars[j].D;
    152         stars[Nstars+1][0].average.R = tstars[j].R;
    153         stars[Nstars+1][0].average.D = tstars[j].D;
    154         stars[Nstars+2][0].average.R = tstars[j].R;
    155         stars[Nstars+2][0].average.D = tstars[j].D;
    156         get2mass_3star_full (&stars[Nstars], &buffer[offset], Nbyte - offset);
     147        dvo_average_init (&newcat->average[Nave]);
     148        dvo_measure_init (&newcat->measure[Nmeas+0]);
     149        dvo_measure_init (&newcat->measure[Nmeas+1]);
     150        dvo_measure_init (&newcat->measure[Nmeas+2]);
     151
     152        newcat->average[Nave].R = tstars[j].R;
     153        newcat->average[Nave].D = tstars[j].D;
     154        get2mass_3star_full (&newcat->measure[Nmeas], &buffer[offset], Nbyte - offset);
    157155
    158156        tstars[j].flag = TRUE;
    159157
    160         Nstars += 3;
    161         if (Nstars >= NSTARS - 3) {
    162           NSTARS += 3000;
    163           REALLOCATE (stars, Stars *, NSTARS);
    164         }
    165       }
    166 
    167       if (!Nstars) {
    168         free (stars);
     158        newcat->average[Nave].Nmeasure = 3;
     159        newcat->average[Nave].measureOffset = Nmeas;
     160
     161        Nave ++;
     162        Nmeas += 3;
     163
     164        CHECK_REALLOCATE (newcat->average, Average, NAVE,  Nave,  1000);
     165        CHECK_REALLOCATE (newcat->measure, Measure, NMEAS, Nmeas, 3000);
     166      }
     167      newcat->Naverage = Nave;
     168      newcat->Nmeasure = Nmeas;
     169
     170      if (!newcat->Naverage) {
     171        dvo_catalog_free (newcat);
     172        free (newcat);
    169173        continue;
    170174      }
    171175
    172       if (DEBUG) fprintf (stderr, "selected %d stars (%10.6f - %10.6f, %10.6f - %10.6f)\n", Nstars,
     176      if (DEBUG) fprintf (stderr, "selected %d stars (%10.6f - %10.6f, %10.6f - %10.6f)\n", (int) newcat->Naverage,
    173177                          region[0].Rmin, region[0].Rmax, region[0].Dmin, region[0].Dmax);
    174178
    … …  
    193197        // assume no input star matches an existing star
    194198        // simply add to the existing table
    195         // load2mass_catalog (&catalog, stars, Nstars);
    196         find_matches_refstars (skylist[0].regions[0], stars, Nstars, &catalog, options);
    197         // load2mass_catalog (&catalog, stars, Nstars);
     199        find_matches_refstars (skylist[0].regions[0], newcat, &catalog, options);
    198200
    199201        dvo_catalog_save (&catalog, VERBOSE);
    200202        dvo_catalog_unlock (&catalog);
    201203        dvo_catalog_free (&catalog);
    202         // free (catalog.filename);
    203         // XXX don't free this! it points to an element of the skytable
    204204      }
    205205
    206206      SkyListFree (skylist);
    207       for (j = 0; j < Nstars; j++) free (stars[j]);
    208       free (stars);
     207      dvo_catalog_free (newcat);
     208      free (newcat);
    209209      VERBOSE = verbose;
    210210    }

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadwise_star_full.c

    r38062 r38593  
    22# include "WISE.h"
    33
    4 int setWISE_ph_qual (Stars *star, char qual);
    5 int setWISE_var_flag_allsky (Stars *star, char qual);
    6 int setWISE_ext_flag_allsky (Stars **star, char ptr);
    7 int setWISE_var_flag_prelim (Stars *star, char qual);
    8 int setWISE_ext_flag_prelim (Stars **star, char ptr);
    9 int setWISE_cc_flag (Stars *star, char qual);
    10 int setWISE_sat_flag (Stars *star, char *ptr);
    11 int setWISE_blend_flag (Stars **star, char *ptr);
     4int setWISE_ph_qual (Measure *measure, char qual);
     5int setWISE_var_flag_allsky (Measure *measure, char qual);
     6int setWISE_ext_flag_allsky (Measure *measure, char ptr);
     7int setWISE_var_flag_prelim (Measure *measure, char qual);
     8int setWISE_ext_flag_prelim (Measure *measure, char ptr);
     9int setWISE_cc_flag (Measure *measure, char qual);
     10int setWISE_sat_flag (Measure *measure, char *ptr);
     11int setWISE_blend_flag (Measure *measure, char *ptr);
    1212
    1313// fill in the data for a WISE quad star.  takes a pointer to the start of the line the
    1414// RA and DEC have already been set
    15 int loadwise_star_allwise (Stars **star, char *line, int Nmax) {
     15int loadwise_star_allwise (Measure *measure, char *line, int Nmax) {
    1616
    1717  int i;
    … …  
    2323
    2424  // I can assign dRA and dDEC to dX and dY if I can consistently set posangle and pltscale
    25   star[0][0].measure.posangle = 0.0;
    26   star[0][0].measure.pltscale = 1.0;
     25  measure[0].posangle = 0.0;
     26  measure[0].pltscale = 1.0;
    2727
    2828  ptr = skipNbounds (ptr, '|', 3, Nmax); // skip: desig, ra, dec,
    29   star[0][0].measure.dXccd = ToShortPixels((strtod (ptr, NULL) + 0.005)); // sig_ra -- add 5/1000 to avoid truncation
     29  measure[0].dXccd = ToShortPixels((strtod (ptr, NULL) + 0.005)); // sig_ra -- add 5/1000 to avoid truncation
    3030  ptr = nextWISEfield (ptr);
    31   star[0][0].measure.dYccd = ToShortPixels((strtod (ptr, NULL) + 0.005)); // sig_dec
     31  measure[0].dYccd = ToShortPixels((strtod (ptr, NULL) + 0.005)); // sig_dec
    3232  ptr = nextWISEfield (ptr);
    3333
    3434  // we only know a single set of values for all 4 bands
    3535  for (i = 1; i < 4; i++) {
    36       star[i][0].measure.dXccd    = star[0][0].measure.dXccd;
    37       star[i][0].measure.dYccd    = star[0][0].measure.dYccd;
    38       star[i][0].measure.posangle = star[0][0].measure.posangle;
    39       star[i][0].measure.pltscale = star[0][0].measure.pltscale;
     36      measure[i].dXccd    = measure[0].dXccd;
     37      measure[i].dYccd    = measure[0].dYccd;
     38      measure[i].posangle = measure[0].posangle;
     39      measure[i].pltscale = measure[0].pltscale;
    4040  }
    4141  ptr = skipNbounds (ptr, '|', 5, Nmax); // skip: sig_radec, glon, glat, elon, elat
    4242
    43   star[0][0].measure.Xccd = strtod (ptr, NULL); // wx
     43  measure[0].Xccd = strtod (ptr, NULL); // wx
    4444  ptr = nextWISEfield (ptr); // skip wx
    45   star[0][0].measure.Yccd = strtod (ptr, NULL); // wy
     45  measure[0].Yccd = strtod (ptr, NULL); // wy
    4646  ptr = nextWISEfield (ptr); // skip wy
    4747
    4848  for (i = 1; i < 4; i++) {
    49       star[i][0].measure.Xccd    = star[0][0].measure.Xccd;
    50       star[i][0].measure.Yccd    = star[0][0].measure.Yccd;
     49      measure[i].Xccd    = measure[0].Xccd;
     50      measure[i].Yccd    = measure[0].Yccd;
    5151  }
    5252  ptr = skipNbounds (ptr, '|', 4, Nmax); // skip: cntr, source_id, coadd_id, src
    … …  
    5555  for (i = 0; i < 4; i++) {
    5656    char *endpoint;
    57       star[i][0].measure.M  = strtod (ptr, &endpoint); // w?mpro
     57      measure[i].M  = strtod (ptr, &endpoint); // w?mpro
    5858      if (endpoint == ptr) {
    59         star[i][0].measure.M  = NAN;
     59        measure[i].M  = NAN;
    6060      }
    6161      ptr = nextWISEfield (ptr);
    62       star[i][0].measure.dM = strtod (ptr, &endpoint); // w?sigmpro
     62      measure[i].dM = strtod (ptr, &endpoint); // w?sigmpro
    6363      if (endpoint == ptr) {
    64         star[i][0].measure.dM  = NAN;
     64        measure[i].dM  = NAN;
    6565      }
    6666      ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w?sigmpro, w?snr
    67       star[i][0].measure.psfChisq = strtod (ptr, NULL); // w?rchi2
    68       star[i][0].measure.psfQF = strtod (ptr, NULL); // w?rchi2
     67      measure[i].psfChisq = strtod (ptr, NULL); // w?rchi2
     68      measure[i].psfQF = strtod (ptr, NULL); // w?rchi2
    6969      ptr = nextWISEfield (ptr); // skip : w1rchi2
    7070
    7171      // init the photFlags field
    72       star[i][0].measure.photFlags = 0;
     72      measure[i].photFlags = 0;
    7373  }
    7474
    … …  
    7676
    7777  // set blend flags for all 4 measures
    78   setWISE_blend_flag (star, ptr); // nb & na both used here
     78  setWISE_blend_flag (measure, ptr); // nb & na both used here
    7979  ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: nb, na
    8080
    8181  for (i = 0; i < 4; i++) {
    82       setWISE_sat_flag (star[i], ptr); // w1sat
     82      setWISE_sat_flag (&measure[i], ptr); // w1sat
    8383      ptr = nextWISEfield (ptr);
    8484  }
    … …  
    8888
    8989  for (i = 0; i < 4; i++) {
    90       setWISE_cc_flag (star[i], ptr[i]); // cc_flags
     90      setWISE_cc_flag (&measure[i], ptr[i]); // cc_flags
    9191  }
    9292  ptr = nextWISEfield (ptr); // skip cc_flags
    … …  
    9494
    9595  // set ext flags for all 4 measures
    96   setWISE_ext_flag_allsky (star, ptr[0]); // ext_flg
     96  setWISE_ext_flag_allsky (measure, ptr[0]); // ext_flg
    9797  ptr = nextWISEfield (ptr); // skip ext_flags
    9898
    9999  for (i = 0; i < 4; i++) {
    100       setWISE_var_flag_allsky (star[i], ptr[i]); // var_flg
     100      setWISE_var_flag_allsky (&measure[i], ptr[i]); // var_flg
    101101  }
    102102  ptr = nextWISEfield (ptr); // skip var_flags
    103103
    104104  for (i = 0; i < 4; i++) {
    105       setWISE_ph_qual (star[i], ptr[i]); // ph_qual
     105      setWISE_ph_qual (&measure[i], ptr[i]); // ph_qual
    106106  }
    107107
    … …  
    123123    // the release is based on data taken in the period 14 January 2010 to 29 April 2010
    124124    if (mjdmean == 0.0) {
    125       star[i][0].measure.t = 0;
     125      measure[i].t = 0;
    126126    } else {
    127       star[i][0].measure.t = ohana_mjd_to_sec (mjdmean);
     127      measure[i].t = ohana_mjd_to_sec (mjdmean);
    128128    }
    129129  }
    130130
    131   star[0][0].measure.photcode  = WISE_W1;
    132   star[0][0].measure.detID   = 0;
    133   star[0][0].measure.imageID = 0;
    134 
    135   star[1][0].measure.photcode  = WISE_W2;
    136   star[1][0].measure.detID   = 0;
    137   star[1][0].measure.imageID = 0;
    138 
    139   star[2][0].measure.photcode  = WISE_W3;
    140   star[2][0].measure.detID   = 0;
    141   star[2][0].measure.imageID = 0;
    142 
    143   star[3][0].measure.photcode  = WISE_W4;
    144   star[3][0].measure.detID   = 0;
    145   star[3][0].measure.imageID = 0;
     131  measure[0].photcode  = WISE_W1;
     132  measure[0].detID   = 0;
     133  measure[0].imageID = 0;
     134
     135  measure[1].photcode  = WISE_W2;
     136  measure[1].detID   = 0;
     137  measure[1].imageID = 0;
     138
     139  measure[2].photcode  = WISE_W3;
     140  measure[2].detID   = 0;
     141  measure[2].imageID = 0;
     142
     143  measure[3].photcode  = WISE_W4;
     144  measure[3].detID   = 0;
     145  measure[3].imageID = 0;
    146146
    147147  return TRUE;
    … …  
    156156// fill in the data for a WISE quad star.  takes a pointer to the start of the line the
    157157// RA and DEC have already been set
    158 int loadwise_star_allsky (Stars **star, char *line, int Nmax) {
     158int loadwise_star_allsky (Measure *measure, char *line, int Nmax) {
    159159
    160160  int i;
    … …  
    166166
    167167  // I can assign dRA and dDEC to dX and dY if I can consistently set posangle and pltscale
    168   star[0][0].measure.posangle = 0.0;
    169   star[0][0].measure.pltscale = 1.0;
     168  measure[0].posangle = 0.0;
     169  measure[0].pltscale = 1.0;
    170170
    171171  ptr = skipNbounds (ptr, '|', 3, Nmax); // skip: desig, ra, dec,
    172   star[0][0].measure.dXccd = ToShortPixels(strtod (ptr, NULL)); // sig_ra
     172  measure[0].dXccd = ToShortPixels(strtod (ptr, NULL)); // sig_ra
    173173  ptr = nextWISEfield (ptr);
    174   star[0][0].measure.dYccd = ToShortPixels(strtod (ptr, NULL)); // sig_dec
     174  measure[0].dYccd = ToShortPixels(strtod (ptr, NULL)); // sig_dec
    175175  ptr = nextWISEfield (ptr);
    176176
    177177  // we only know a single set of values for all 4 bands
    178178  for (i = 1; i < 4; i++) {
    179       star[i][0].measure.dXccd    = star[0][0].measure.dXccd;
    180       star[i][0].measure.dYccd    = star[0][0].measure.dYccd;
    181       star[i][0].measure.posangle = star[0][0].measure.posangle;
    182       star[i][0].measure.pltscale = star[0][0].measure.pltscale;
     179      measure[i].dXccd    = measure[0].dXccd;
     180      measure[i].dYccd    = measure[0].dYccd;
     181      measure[i].posangle = measure[0].posangle;
     182      measure[i].pltscale = measure[0].pltscale;
    183183  }
    184184  ptr = skipNbounds (ptr, '|', 5, Nmax); // skip: sig_radec, glon, glat, elon, elat
    185185
    186   star[0][0].measure.Xccd = strtod (ptr, NULL); // wx
     186  measure[0].Xccd = strtod (ptr, NULL); // wx
    187187  ptr = nextWISEfield (ptr); // skip wx
    188   star[0][0].measure.Yccd = strtod (ptr, NULL); // wy
     188  measure[0].Yccd = strtod (ptr, NULL); // wy
    189189  ptr = nextWISEfield (ptr); // skip wy
    190190
    191191  for (i = 1; i < 4; i++) {
    192       star[i][0].measure.Xccd    = star[0][0].measure.Xccd;
    193       star[i][0].measure.Yccd    = star[0][0].measure.Yccd;
     192      measure[i].Xccd    = measure[0].Xccd;
     193      measure[i].Yccd    = measure[0].Yccd;
    194194  }
    195195  ptr = skipNbounds (ptr, '|', 4, Nmax); // skip: cntr, source_id, coadd_id, src
    … …  
    198198  for (i = 0; i < 4; i++) {
    199199    char *endpoint;
    200       star[i][0].measure.M  = strtod (ptr, &endpoint); // w?mpro
     200      measure[i].M  = strtod (ptr, &endpoint); // w?mpro
    201201      if (endpoint == ptr) {
    202         star[i][0].measure.M  = NAN;
     202        measure[i].M  = NAN;
    203203      }
    204204      ptr = nextWISEfield (ptr);
    205       star[i][0].measure.dM = strtod (ptr, &endpoint); // w?sigmpro
     205      measure[i].dM = strtod (ptr, &endpoint); // w?sigmpro
    206206      if (endpoint == ptr) {
    207         star[i][0].measure.dM  = NAN;
     207        measure[i].dM  = NAN;
    208208      }
    209209      ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w?sigmpro, w?snr
    210       star[i][0].measure.psfChisq = strtod (ptr, NULL); // w?rchi2
     210      measure[i].psfChisq = strtod (ptr, NULL); // w?rchi2
    211211      ptr = nextWISEfield (ptr); // skip : w1rchi2
    212212
    213213      // init the photFlags field
    214       star[i][0].measure.photFlags = 0;
     214      measure[i].photFlags = 0;
    215215  }
    216216
    … …  
    218218
    219219  // set blend flags for all 4 measures
    220   setWISE_blend_flag (star, ptr); // nb & na both used here
     220  setWISE_blend_flag (measure, ptr); // nb & na both used here
    221221  ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: nb, na
    222222
    223223  for (i = 0; i < 4; i++) {
    224       setWISE_sat_flag (star[i], ptr); // w1sat
     224      setWISE_sat_flag (&measure[i], ptr); // w1sat
    225225      ptr = nextWISEfield (ptr);
    226226  }
    … …  
    228228
    229229  for (i = 0; i < 4; i++) {
    230       setWISE_cc_flag (star[i], ptr[i]); // cc_flg
     230      setWISE_cc_flag (&measure[i], ptr[i]); // cc_flg
    231231  }
    232232  ptr = nextWISEfield (ptr); // skip cc_flags
    233233
    234234  // set ext flags for all 4 measures
    235   setWISE_ext_flag_allsky (star, ptr[0]); // ext_flg
     235  setWISE_ext_flag_allsky (measure, ptr[0]); // ext_flg
    236236  ptr = nextWISEfield (ptr); // skip ext_flags
    237237
    238238  for (i = 0; i < 4; i++) {
    239       setWISE_var_flag_allsky (star[i], ptr[i]); // var_flg
     239      setWISE_var_flag_allsky (&measure[i], ptr[i]); // var_flg
    240240  }
    241241  ptr = nextWISEfield (ptr); // skip var_flags
    … …  
    261261    // the release is based on data taken in the period 14 January 2010 to 29 April 2010
    262262    if (mjdmean == 0.0) {
    263       star[i][0].measure.t = 0;
     263      measure[i].t = 0;
    264264    } else {
    265       star[i][0].measure.t = ohana_mjd_to_sec (mjdmean);
     265      measure[i].t = ohana_mjd_to_sec (mjdmean);
    266266    }
    267267  }
    268268
    269   star[0][0].measure.photcode  = WISE_W1;
    270   star[0][0].measure.detID   = 0;
    271   star[0][0].measure.imageID = 0;
    272 
    273   star[1][0].measure.photcode  = WISE_W2;
    274   star[1][0].measure.detID   = 0;
    275   star[1][0].measure.imageID = 0;
    276 
    277   star[2][0].measure.photcode  = WISE_W3;
    278   star[2][0].measure.detID   = 0;
    279   star[2][0].measure.imageID = 0;
    280 
    281   star[3][0].measure.photcode  = WISE_W4;
    282   star[3][0].measure.detID   = 0;
    283   star[3][0].measure.imageID = 0;
     269  measure[0].photcode  = WISE_W1;
     270  measure[0].detID   = 0;
     271  measure[0].imageID = 0;
     272
     273  measure[1].photcode  = WISE_W2;
     274  measure[1].detID   = 0;
     275  measure[1].imageID = 0;
     276
     277  measure[2].photcode  = WISE_W3;
     278  measure[2].detID   = 0;
     279  measure[2].imageID = 0;
     280
     281  measure[3].photcode  = WISE_W4;
     282  measure[3].detID   = 0;
     283  measure[3].imageID = 0;
    284284
    285285  return TRUE;
    … …  
    288288// fill in the data for a WISE quad star.  takes a pointer to the start of the line the
    289289// RA and DEC have already been set
    290 int loadwise_star_prelim (Stars **star, char *line, int Nmax) {
     290int loadwise_star_prelim (Measure *measure, char *line, int Nmax) {
    291291
    292292  int i;
    … …  
    298298
    299299  // I can assign dRA and dDEC to dX and dY if I can consistently set posangle and pltscale
    300   star[0][0].measure.posangle = 0.0;
    301   star[0][0].measure.pltscale = 1.0;
     300  measure[0].posangle = 0.0;
     301  measure[0].pltscale = 1.0;
    302302
    303303  ptr = skipNbounds (ptr, '|', 3, Nmax); // skip: desig, ra, dec,
    304   star[0][0].measure.dXccd = ToShortPixels(strtod (ptr, NULL)); // sig_ra
     304  measure[0].dXccd = ToShortPixels(strtod (ptr, NULL)); // sig_ra
    305305  ptr = nextWISEfield (ptr);
    306   star[0][0].measure.dYccd = ToShortPixels(strtod (ptr, NULL)); // sig_dec
     306  measure[0].dYccd = ToShortPixels(strtod (ptr, NULL)); // sig_dec
    307307  ptr = nextWISEfield (ptr);
    308308
    309309  // we only know a single set of values for all 4 bands
    310310  for (i = 1; i < 4; i++) {
    311       star[i][0].measure.dXccd    = star[0][0].measure.dXccd;
    312       star[i][0].measure.dYccd    = star[0][0].measure.dYccd;
    313       star[i][0].measure.posangle = star[0][0].measure.posangle;
    314       star[i][0].measure.pltscale = star[0][0].measure.pltscale;
     311      measure[i].dXccd    = measure[0].dXccd;
     312      measure[i].dYccd    = measure[0].dYccd;
     313      measure[i].posangle = measure[0].posangle;
     314      measure[i].pltscale = measure[0].pltscale;
    315315  }
    316316  ptr = skipNbounds (ptr, '|', 5, Nmax); // skip: sig_radec, glon, glat, elon, elat
    317317
    318   star[0][0].measure.Xccd = strtod (ptr, NULL); // wx
     318  measure[0].Xccd = strtod (ptr, NULL); // wx
    319319  ptr = nextWISEfield (ptr); // skip wx
    320   star[0][0].measure.Yccd = strtod (ptr, NULL); // wy
     320  measure[0].Yccd = strtod (ptr, NULL); // wy
    321321  ptr = nextWISEfield (ptr); // skip wy
    322322
    323323  for (i = 1; i < 4; i++) {
    324       star[i][0].measure.Xccd    = star[0][0].measure.Xccd;
    325       star[i][0].measure.Yccd    = star[0][0].measure.Yccd;
     324      measure[i].Xccd    = measure[0].Xccd;
     325      measure[i].Yccd    = measure[0].Yccd;
    326326  }
    327327  ptr = skipNbounds (ptr, '|', 4, Nmax); // skip: cntr, source_id, coadd_id, src
    … …  
    329329  // W1
    330330  for (i = 0; i < 4; i++) {
    331       star[i][0].measure.M  = strtod (ptr, NULL); // w?mpro
     331      measure[i].M  = strtod (ptr, NULL); // w?mpro
    332332      ptr = nextWISEfield (ptr);
    333       star[i][0].measure.dM = strtod (ptr, NULL); // w?sigmpro
     333      measure[i].dM = strtod (ptr, NULL); // w?sigmpro
    334334      ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w?sigmpro, w?snr
    335       star[i][0].measure.psfChisq = strtod (ptr, NULL); // w?rchi2
     335      measure[i].psfChisq = strtod (ptr, NULL); // w?rchi2
    336336      ptr = nextWISEfield (ptr); // skip : w1rchi2
    337337
    338338      // init the photFlags field
    339       star[i][0].measure.photFlags = 0;
     339      measure[i].photFlags = 0;
    340340  }
    341341
    … …  
    372372  double jd = 2455263.0; // NOTE : WISE prelim release does not contain per-detection time info.
    373373  // the release is based on data taken in the period 14 January 2010 to 29 April 2010
    374   star[0][0].measure.t = ohana_jd_to_sec (jd);
    375   star[1][0].measure.t = star[0][0].measure.t;
    376   star[2][0].measure.t = star[0][0].measure.t;
    377   star[3][0].measure.t = star[0][0].measure.t;
    378 
    379   star[0][0].measure.photcode  = WISE_W1;
    380   star[0][0].measure.detID   = 0;
    381   star[0][0].measure.imageID = 0;
    382 
    383   star[1][0].measure.photcode  = WISE_W2;
    384   star[1][0].measure.detID   = 0;
    385   star[1][0].measure.imageID = 0;
    386 
    387   star[2][0].measure.photcode  = WISE_W3;
    388   star[2][0].measure.detID   = 0;
    389   star[2][0].measure.imageID = 0;
    390 
    391   star[3][0].measure.photcode  = WISE_W4;
    392   star[3][0].measure.detID   = 0;
    393   star[3][0].measure.imageID = 0;
     374  measure[0].t = ohana_jd_to_sec (jd);
     375  measure[1].t = measure[0].t;
     376  measure[2].t = measure[0].t;
     377  measure[3].t = measure[0].t;
     378
     379  measure[0].photcode  = WISE_W1;
     380  measure[0].detID   = 0;
     381  measure[0].imageID = 0;
     382
     383  measure[1].photcode  = WISE_W2;
     384  measure[1].detID   = 0;
     385  measure[1].imageID = 0;
     386
     387  measure[2].photcode  = WISE_W3;
     388  measure[2].detID   = 0;
     389  measure[2].imageID = 0;
     390
     391  measure[3].photcode  = WISE_W4;
     392  measure[3].detID   = 0;
     393  measure[3].imageID = 0;
    394394
    395395  return TRUE;
    … …  
    436436# define FLAG_VARIABLE_LEVEL3 0x40000000 // var_flg == 8 or 9
    437437
    438 int setWISE_blend_flag (Stars **star, char *ptr) {
     438int setWISE_blend_flag (Measure *measure, char *ptr) {
    439439
    440440    int nb = atoi (ptr);
    … …  
    445445
    446446    if (na == 0) {
    447         star[0][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
    448         star[1][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
    449         star[2][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
    450         star[3][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
     447        measure[0].photFlags |= FLAG_BLEND_ACTIVE;
     448        measure[1].photFlags |= FLAG_BLEND_ACTIVE;
     449        measure[2].photFlags |= FLAG_BLEND_ACTIVE;
     450        measure[3].photFlags |= FLAG_BLEND_ACTIVE;
    451451    } else {
    452         star[0][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
    453         star[1][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
    454         star[2][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
    455         star[3][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
     452        measure[0].photFlags |= FLAG_BLEND_PASSIVE;
     453        measure[1].photFlags |= FLAG_BLEND_PASSIVE;
     454        measure[2].photFlags |= FLAG_BLEND_PASSIVE;
     455        measure[3].photFlags |= FLAG_BLEND_PASSIVE;
    456456    }
    457457    return TRUE;
    458458}
    459459
    460 int setWISE_sat_flag (Stars *star, char *ptr) {
     460int setWISE_sat_flag (Measure *measure, char *ptr) {
    461461
    462462    float sat = strtod (ptr, NULL);
    463463    if (sat > 0.0) {
    464         star[0].measure.photFlags |= FLAG_SATURATED_PIX;
     464        measure[0].photFlags |= FLAG_SATURATED_PIX;
    465465    }
    466466    return TRUE;
    467467}
    468468
    469 int setWISE_cc_flag (Stars *star, char qual) {
     469int setWISE_cc_flag (Measure *measure, char qual) {
    470470
    471471  switch (qual) {
    472     case 'p': star[0].measure.photFlags |= FLAG_CC_PERSIST; break;
    473     case 'h': star[0].measure.photFlags |= FLAG_CC_HALO;    break;
    474     case 'd': star[0].measure.photFlags |= FLAG_CC_SPIKE;   break;
    475     case 'o': star[0].measure.photFlags |= FLAG_CC_GHOST;   break;
    476 
    477     case 'P': star[0].measure.photFlags |= FLAG_CC_PERSIST | FLAG_CC_SPURIOUS; break;
    478     case 'H': star[0].measure.photFlags |= FLAG_CC_HALO    | FLAG_CC_SPURIOUS; break;
    479     case 'D': star[0].measure.photFlags |= FLAG_CC_SPIKE   | FLAG_CC_SPURIOUS; break;
    480     case 'O': star[0].measure.photFlags |= FLAG_CC_GHOST   | FLAG_CC_SPURIOUS; break;
     472    case 'p': measure[0].photFlags |= FLAG_CC_PERSIST; break;
     473    case 'h': measure[0].photFlags |= FLAG_CC_HALO;    break;
     474    case 'd': measure[0].photFlags |= FLAG_CC_SPIKE;   break;
     475    case 'o': measure[0].photFlags |= FLAG_CC_GHOST;   break;
     476
     477    case 'P': measure[0].photFlags |= FLAG_CC_PERSIST | FLAG_CC_SPURIOUS; break;
     478    case 'H': measure[0].photFlags |= FLAG_CC_HALO    | FLAG_CC_SPURIOUS; break;
     479    case 'D': measure[0].photFlags |= FLAG_CC_SPIKE   | FLAG_CC_SPURIOUS; break;
     480    case 'O': measure[0].photFlags |= FLAG_CC_GHOST   | FLAG_CC_SPURIOUS; break;
    481481
    482482    case '0': break;
    … …  
    487487}
    488488
    489 int setWISE_ext_flag_allsky (Stars **star, char value) {
     489int setWISE_ext_flag_allsky (Measure *measure, char value) {
    490490
    491491  switch (value) {
    … …  
    493493      return TRUE;
    494494    case '1':
    495       star[0][0].measure.photFlags |= FLAG_EXTENDED;
    496       star[1][0].measure.photFlags |= FLAG_EXTENDED;
    497       star[2][0].measure.photFlags |= FLAG_EXTENDED;
    498       star[3][0].measure.photFlags |= FLAG_EXTENDED;
     495      measure[0].photFlags |= FLAG_EXTENDED;
     496      measure[1].photFlags |= FLAG_EXTENDED;
     497      measure[2].photFlags |= FLAG_EXTENDED;
     498      measure[3].photFlags |= FLAG_EXTENDED;
    499499      return TRUE;
    500500    case '2':
    501       star[0][0].measure.photFlags |= FLAG_EXT_IN_XSC;
    502       star[1][0].measure.photFlags |= FLAG_EXT_IN_XSC;
    503       star[2][0].measure.photFlags |= FLAG_EXT_IN_XSC;
    504       star[3][0].measure.photFlags |= FLAG_EXT_IN_XSC;
     501      measure[0].photFlags |= FLAG_EXT_IN_XSC;
     502      measure[1].photFlags |= FLAG_EXT_IN_XSC;
     503      measure[2].photFlags |= FLAG_EXT_IN_XSC;
     504      measure[3].photFlags |= FLAG_EXT_IN_XSC;
    505505      return TRUE;
    506506    case '3':
    507       star[0][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
    508       star[1][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
    509       star[2][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
    510       star[3][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
     507      measure[0].photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
     508      measure[1].photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
     509      measure[2].photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
     510      measure[3].photFlags |= FLAG_EXTENDED | FLAG_EXT_IN_XSC;
    511511      return TRUE;
    512512    case '4':
    513       star[0][0].measure.photFlags |= FLAG_EXT_BY_XSC;
    514       star[1][0].measure.photFlags |= FLAG_EXT_BY_XSC;
    515       star[2][0].measure.photFlags |= FLAG_EXT_BY_XSC;
    516       star[3][0].measure.photFlags |= FLAG_EXT_BY_XSC;
     513      measure[0].photFlags |= FLAG_EXT_BY_XSC;
     514      measure[1].photFlags |= FLAG_EXT_BY_XSC;
     515      measure[2].photFlags |= FLAG_EXT_BY_XSC;
     516      measure[3].photFlags |= FLAG_EXT_BY_XSC;
    517517      return TRUE;
    518518    case '5':
    519       star[0][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
    520       star[1][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
    521       star[2][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
    522       star[3][0].measure.photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
     519      measure[0].photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
     520      measure[1].photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
     521      measure[2].photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
     522      measure[3].photFlags |= FLAG_EXTENDED | FLAG_EXT_BY_XSC;
    523523      return TRUE;
    524524    default:
    … …  
    529529}
    530530
    531 int setWISE_ext_flag_prelim (Stars **star, char value) {
     531int setWISE_ext_flag_prelim (Measure *measure, char value) {
    532532
    533533  switch (value) {
    … …  
    535535      return TRUE;
    536536    case '1':
    537       star[0][0].measure.photFlags |= FLAG_EXTENDED;
    538       star[1][0].measure.photFlags |= FLAG_EXTENDED;
    539       star[2][0].measure.photFlags |= FLAG_EXTENDED;
    540       star[3][0].measure.photFlags |= FLAG_EXTENDED;
     537      measure[0].photFlags |= FLAG_EXTENDED;
     538      measure[1].photFlags |= FLAG_EXTENDED;
     539      measure[2].photFlags |= FLAG_EXTENDED;
     540      measure[3].photFlags |= FLAG_EXTENDED;
    541541      return TRUE;
    542542    default:
    … …  
    548548
    549549// NOTE: var flag definition changed slightly between prelim & allsky
    550 int setWISE_var_flag_allsky (Stars *star, char qual) {
     550int setWISE_var_flag_allsky (Measure *measure, char qual) {
    551551
    552552  switch (qual) {
    … …  
    559559    case '4':
    560560    case '5':
    561         star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL1;
     561        measure[0].photFlags |= FLAG_VARIABLE_LEVEL1;
    562562        return TRUE;
    563563    case '6':
    564564    case '7':
    565         star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL2;
     565        measure[0].photFlags |= FLAG_VARIABLE_LEVEL2;
    566566        return TRUE;
    567567    case '8':
    568568    case '9':
    569       star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL3;
     569      measure[0].photFlags |= FLAG_VARIABLE_LEVEL3;
    570570      return TRUE;
    571571    default:
    … …  
    575575}
    576576
    577 int setWISE_var_flag_prelim (Stars *star, char qual) {
     577int setWISE_var_flag_prelim (Measure *measure, char qual) {
    578578
    579579  switch (qual) {
    … …  
    584584    case '3':
    585585    case '4':
    586         star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL1;
     586        measure[0].photFlags |= FLAG_VARIABLE_LEVEL1;
    587587        return TRUE;
    588588    case '5':
    589589    case '6':
    590590    case '7':
    591         star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL2;
     591        measure[0].photFlags |= FLAG_VARIABLE_LEVEL2;
    592592        return TRUE;
    593593    case '8':
    594594    case '9':
    595       star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL3;
     595      measure[0].photFlags |= FLAG_VARIABLE_LEVEL3;
    596596      return TRUE;
    597597    default:
    … …  
    601601}
    602602
    603 int setWISE_ph_qual (Stars *star, char qual) {
     603int setWISE_ph_qual (Measure *measure, char qual) {
    604604
    605605  switch (qual) {
    606     case 'A': star[0].measure.photFlags |= FLAG_PH_A; break;
    607     case 'B': star[0].measure.photFlags |= FLAG_PH_B; break;
    608     case 'C': star[0].measure.photFlags |= FLAG_PH_C; break;
    609     case 'U': star[0].measure.photFlags |= FLAG_PH_U; break;
    610     case 'X': star[0].measure.photFlags |= FLAG_PH_X; break;
    611     case 'Z': star[0].measure.photFlags |= FLAG_PH_Z; break;
     606    case 'A': measure[0].photFlags |= FLAG_PH_A; break;
     607    case 'B': measure[0].photFlags |= FLAG_PH_B; break;
     608    case 'C': measure[0].photFlags |= FLAG_PH_C; break;
     609    case 'U': measure[0].photFlags |= FLAG_PH_U; break;
     610    case 'X': measure[0].photFlags |= FLAG_PH_X; break;
     611    case 'Z': measure[0].photFlags |= FLAG_PH_Z; break;
    612612    default:
    613613      fprintf (stderr, "error in ph_flag: %c\n", qual);

  • branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/replace_match.c

    r37036 r38593  
    11# include "addstar.h"
    22
    3 int replace_match (Average *average, Measure *measure, Stars *star) {
     3int replace_match (Average *average, Measure *measure, Measure *newmeas, off_t *found) {
    44
    55  int i, j, m;
    … …  
    1111  for (i = 0; i < average[0].Nmeasure; i++) {
    1212    j = i + m;
    13     if (measure[j].photcode != star[0].measure.photcode) continue;
    14     measure[j].R = star[0].average.R;
    15     measure[j].D = star[0].average.D;
    16     measure[j].M  = star[0].measure.M;
    17     measure[j].dM = star[0].measure.dM;
    18     star[0].found = average[0].measureOffset + i;
     13    if (measure[j].photcode != newmeas->photcode) continue;
     14    measure[j].R  = newmeas->R;
     15    measure[j].D  = newmeas->D;
     16    measure[j].M  = newmeas->M;
     17    measure[j].dM = newmeas->dM;
     18    *found = average[0].measureOffset + i;
    1919    return (TRUE);
    2020  }
Note: See TracChangeset for help on using the changeset viewer.