IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37116


Ignore:
Timestamp:
Jul 24, 2014, 5:10:56 PM (12 years ago)
Author:
eugene
Message:

merge relphot changes from eam_branches/ipp-20140717

Location:
trunk/Ohana
Files:
1 deleted
48 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/addstar/Makefile

    r35263 r37116  
    7878$(SRC)/parse_time.$(ARCH).o \
    7979$(SRC)/replace_match.$(ARCH).o \
    80 $(SRC)/resort_catalog.$(ARCH).o \
     80$(SRC)/resort_catalogs.$(ARCH).o \
    8181$(SRC)/resort_threaded.$(ARCH).o \
    8282$(SRC)/resort_unthreaded.$(ARCH).o \
    83 $(SRC)/resort_unthreaded_catalogs.$(ARCH).o \
     83$(SRC)/resort_catalog.$(ARCH).o \
    8484$(SRC)/StarOps.$(ARCH).o \
    8585$(SRC)/ReadStarsFITS.$(ARCH).o \
    … …  
    105105$(SRC)/addstar_client.$(ARCH).o \
    106106$(SRC)/args_parallel_client.$(ARCH).o \
    107 $(SRC)/resort_catalog.$(ARCH).o \
     107$(SRC)/resort_catalogs.$(ARCH).o \
    108108$(SRC)/resort_threaded.$(ARCH).o \
    109109$(SRC)/resort_unthreaded.$(ARCH).o \
    110 $(SRC)/resort_unthreaded_catalogs.$(ARCH).o
     110$(SRC)/resort_catalog.$(ARCH).o
    111111
    112112ADDSTARD = \

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

    r37036 r37116  
    237237Stars     *rd_gsc                 PROTO((char *filename, unsigned int *nstars));
    238238int        replace_match          PROTO((Average *average, Measure *measure, Stars *star));
    239 int        resort_threaded        PROTO((AddstarClientOptions *options, SkyTable *sky));
    240 int        resort_unthreaded      PROTO((AddstarClientOptions *options, SkyTable *sky, int hostID, char *hostpath));
    241 int        resort_unthreaded_catalogs PROTO((AddstarClientOptions *options, SkyList *skylist, int hostID, char *hostpath));
     239int        resort_catalogs        PROTO((AddstarClientOptions *options, SkyTable *sky));
     240int        resort_catalogs_parallel PROTO((AddstarClientOptions *options, SkyList *sky));
     241int        resort_threaded        PROTO((SkyList *skylist, int ForceSort));
     242int        resort_unthreaded      PROTO((SkyList *skylist, int ForceSort));
    242243void       resort_catalog         PROTO((Catalog *catalog));
    243244void       resort_catalog_old     PROTO((Catalog *catalog));
     245
    244246Stars     *ReadStarsFITS          PROTO((FILE *f, Header *header, Header *in_theader, unsigned int *nstars));
    245247Stars     *ReadStarsTEXT          PROTO((FILE *f, unsigned int *nstars));
    … …  
    253255off_t     *init_lensing_links     PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
    254256off_t     *init_missing_links     PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing));
    255 off_t      add_meas_link          PROTO((Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE));
    256 off_t      add_miss_link          PROTO((Average *average, off_t *next_miss, off_t Nmissing));
    257 off_t      add_lens_link          PROTO((Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING));
     257int        add_meas_link          PROTO((Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE));
     258int        add_miss_link          PROTO((Average *average, off_t *next_miss, off_t Nmissing));
     259int        add_lens_link          PROTO((Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING));
    258260off_t     *build_measure_links    PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
    259261off_t     *build_lensing_links    PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
    … …  
    316318void *memrchr(const void *s, int c, size_t n);
    317319int addstar_create_ID ();
     320int strextend (char *input, char *format,...);
    318321
    319322/**

  • trunk/Ohana/src/addstar/src/BoundaryTreeIO.c

    r34260 r37116  
    11# include "addstar.h"
     2
     3/** this file is deprecated (functions moved to libdvo) ***/
    24
    35# define GET_COLUMN_NEW(OUT,NAME,TYPE)                                  \

  • trunk/Ohana/src/addstar/src/addstar.c

    r37036 r37116  
    4040 
    4141  if (options.mode == ADDSTAR_MODE_RESORT) {
    42     if (NTHREADS == 0) {
    43       resort_unthreaded (&options, sky, 0, NULL);
    44     } else {
    45       resort_threaded (&options, sky);
    46     }
     42    resort_catalogs (&options, sky);
    4743    exit (0);
    4844  }
    … …  
    7167      skylist = SkyListForStars (sky, -1, stars, Nstars);
    7268      break;
    73     case ADDSTAR_MODE_RESORT:
    7469    case ADDSTAR_MODE_REFCAT:
    7570      skylist = SkyListByPatch (sky, -1, &UserPatch);
    … …  
    161156        if (Nsubset) free (subset);
    162157        break;
    163       case ADDSTAR_MODE_RESORT:
    164         if (options.nosort == 3) catalog.sorted = FALSE;
    165 
    166         // no need to resort empty catalogs
    167         if (catalog.Naves_disk == 0) {
    168           dvo_catalog_unlock (&catalog);
    169           dvo_catalog_free (&catalog);
    170           continue;
    171         }
    172 
    173         if (OLD_RESORT) {
    174           resort_catalog_old (&catalog);
    175         } else {
    176           resort_catalog (&catalog);
    177         }
    178         Nsubset = 1;
    179         break;
    180158    }
    181159    /* report total updated values */

  • trunk/Ohana/src/addstar/src/addstar_client.c

    r33963 r37116  
    2525
    2626  if (options.mode == ADDSTAR_MODE_RESORT) {
    27     resort_unthreaded (&options, sky, HOST_ID, HOSTDIR);
     27    resort_catalogs (&options, sky);
    2828    exit (0);
    2929  }

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

    r35760 r37116  
    3434    remove_argument (N, &argc, argv);
    3535  }
     36  OLD_RESORT = FALSE;
     37  if ((N = get_argument (argc, argv, "-old-resort"))) {
     38    remove_argument (N, &argc, argv);
     39    OLD_RESORT = TRUE;
     40  }
     41
    3642  if ((N = get_argument (argc, argv, "-create-id"))) {
    3743    options.mode = ADDSTAR_MODE_CREATE_ID;
    … …  
    350356  // XXX for the moment, make this selection manual.  it needs to be automatic
    351357  // based on the state of the SkyTable
     358  HOST_ID = 0;
     359  HOSTDIR = NULL;
    352360  PARALLEL = FALSE;
    353361  if ((N = get_argument (argc, argv, "-parallel"))) {
    … …  
    466474}
    467475
     476int strextend (char *input, char *format,...) {
     477
     478  char tmpextra[1024], tmpline[1024];
     479  va_list argp;
     480
     481  va_start (argp, format);
     482  vsnprintf (tmpextra, 1024, format, argp);
     483  snprintf (tmpline, 1024, "%s %s", input, tmpextra);
     484  strcpy (input, tmpline);
     485
     486  return TRUE;
     487}
     488
    468489/** addstar modes:
    469490 

  • trunk/Ohana/src/addstar/src/args_parallel_client.c

    r33963 r37116  
    8787  }
    8888
     89  OLD_RESORT = FALSE;
     90  if ((N = get_argument (argc, argv, "-old-resort"))) {
     91    remove_argument (N, &argc, argv);
     92    OLD_RESORT = TRUE;
     93  }
     94
    8995  /* extra error messages */
    9096  VERBOSE = FALSE;
    … …  
    118124  exit (2);
    119125}
     126
     127int strextend (char *input, char *format,...) {
     128
     129  char tmpextra[1024], tmpline[1024];
     130  va_list argp;
     131
     132  va_start (argp, format);
     133  vsnprintf (tmpextra, 1024, format, argp);
     134  snprintf (tmpline, 1024, "%s %s", input, tmpextra);
     135  strcpy (input, tmpline);
     136
     137  return TRUE;
     138}
     139

  • trunk/Ohana/src/addstar/src/build_links.c

    r37036 r37116  
    102102
    103103/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */
    104 off_t add_meas_link (Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE) {
     104int add_meas_link (Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE) {
    105105
    106106  off_t k, m;
    … …  
    133133  return (TRUE);
    134134}
    135 
    136 /* Missing does not carry enough information to reconstruct the links
    137    we must always save the missing table, if it exists */
    138135
    139136Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next_meas) {
    … …  
    160157}
    161158
     159/*******************************************************************************************/
     160
    162161/* build the initial links assuming the table is sorted */
    163162off_t *init_missing_links (Average *average, off_t Naverage, Missing *missing, off_t Nmissing) {
    … …  
    186185}
    187186
    188 off_t add_miss_link (Average *average, off_t *next_miss, off_t Nmissing) {
     187int add_miss_link (Average *average, off_t *next_miss, off_t Nmissing) {
    189188
    190189  off_t k, m;
    … …  
    205204}
    206205
     206/* Missing does not carry enough information to reconstruct the links
     207   we must always save the missing table, if it exists */
     208
    207209Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss) {
    208210
    … …  
    307309
    308310/* average[].lensingOffset, average[].Nlensing are valid within an addstar run */
    309 off_t add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) {
     311int add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) {
    310312
    311313  off_t k, m;

  • trunk/Ohana/src/addstar/src/findskycell.c

    r37036 r37116  
    2727int apply_tree (char *treefile, char *datafile);
    2828
    29 float SCALE = 1.0;
    30 int NX_SUB = 1;
    31 int NY_SUB = 1;
     29float SCALE  = 1.0;
     30float NX_SUB = 1.0;
     31float NY_SUB = 1.0;
    3232
    3333double R_MIN = 0.0;
    … …  
    217217  ALLOCATE (tree.Xo,   double *, tree.Nzone);
    218218  ALLOCATE (tree.Yo,   double *, tree.Nzone);
    219   ALLOCATE (tree.dX,      int *, tree.Nzone);
    220   ALLOCATE (tree.dY,      int *, tree.Nzone);
     219  ALLOCATE (tree.dX,    float *, tree.Nzone);
     220  ALLOCATE (tree.dY,    float *, tree.Nzone);
    221221  ALLOCATE (tree.cell,    int *, tree.Nzone);
    222222  ALLOCATE (tree.name,  char **, tree.Nzone);
    … …  
    251251    ALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
    252252    ALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
    253     ALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
    254     ALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
     253    ALLOCATE (tree.dX[zone],    float, tree.NBAND[zone]);
     254    ALLOCATE (tree.dY[zone],    float, tree.NBAND[zone]);
    255255    ALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
    256256    ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
    … …  
    260260      tree.Xo[zone][band] = NAN;
    261261      tree.Yo[zone][band] = NAN;
    262       tree.dX[zone][band] = -1;
    263       tree.dY[zone][band] = -1;
     262      tree.dX[zone][band] = NAN;
     263      tree.dY[zone][band] = NAN;
    264264      tree.cell[zone][band] = -1;
    265265      ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH);
    … …  
    286286      REALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
    287287      REALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
    288       REALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
    289       REALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
     288      REALLOCATE (tree.dX[zone],    float, tree.NBAND[zone]);
     289      REALLOCATE (tree.dY[zone],    float, tree.NBAND[zone]);
    290290      REALLOCATE (tree.cell[zone],    int, tree.NBAND[zone]);
    291291      REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
    … …  
    295295        tree.Xo[zone][band] = NAN;
    296296        tree.Yo[zone][band] = NAN;
    297         tree.dX[zone][band] = -1;
    298         tree.dY[zone][band] = -1;
     297        tree.dX[zone][band] = NAN;
     298        tree.dY[zone][band] = NAN;
    299299        tree.cell[zone][j] = -1;
    300300        ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);

  • trunk/Ohana/src/addstar/src/resort_catalog.c

    r37036 r37116  
    77void resort_catalog_old (Catalog *catalog) {
    88
    9   off_t *next_meas;
    10   off_t Naves, Nmeas;
     9  off_t *next_meas, *next_lens;
     10  off_t Naves, Nmeas, Nlens;
    1111
    1212  if (catalog[0].sorted == TRUE) return;
    … …  
    1616  /* internal counters */
    1717  Nmeas = catalog[0].Nmeasure;
     18  Nlens = catalog[0].Nlensing;
    1819  Naves = catalog[0].Naverage;
    1920 
    2021  /* set up pointers for linked list of measure, missing */
    2122  next_meas = build_measure_links (catalog[0].average, Naves, catalog[0].measure, Nmeas);
     23  next_lens = build_lensing_links (catalog[0].average, Naves, catalog[0].lensing, Nlens);
    2224
    2325  catalog[0].sorted = TRUE;
    2426  catalog[0].measure = sort_measure (catalog[0].average, Naves, catalog[0].measure, Nmeas, next_meas);
     27  catalog[0].lensing = sort_lensing (catalog[0].average, Naves, catalog[0].lensing, Nlens, next_lens);
    2528
    2629  MARKTIME ("  match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeas, (long long) Naves);

  • trunk/Ohana/src/addstar/src/resort_threaded.c

    r37036 r37116  
    1313  int iThread;
    1414  int state;
    15   int nosort;
     15  int forcesort;
    1616  char *filename;
    1717  SkyRegion *region;
    … …  
    3737    }
    3838
     39    // chose the catalog file (local or remote?)
     40    char hostfile[1024];
     41    snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, threadData->region->name);
     42    catalog.filename  = HOST_ID ? hostfile : threadData->filename;
     43
    3944    // set the parameters which guide catalog open/load/create
    40     catalog.filename  = threadData->filename;
    4145    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
    4246    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
    … …  
    6468
    6569    // this is an overloaded value to mean 'force sort'
    66     if (threadData->nosort == 3) catalog.sorted = FALSE;
     70    if (threadData->forcesort) catalog.sorted = FALSE;
    6771
    6872    if (OLD_RESORT) {
    … …  
    8892}
    8993
    90 int resort_threaded (AddstarClientOptions *options, SkyTable *sky) {
     94int resort_threaded (SkyList *skylist, int ForceSort) {
    9195
    9296  int j, launched, done, Nwait, Nrun;
    … …  
    97101  ThreadData *threadData;
    98102  pthread_t *threads;
    99 
    100   double dtime;
    101   struct timeval start, stop;
    102   gettimeofday (&start, NULL);
    103 
    104   SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch);
    105  
    106   // in these cases, limit the sky catalogs to an existing subset
    107   if (options->existing_regions) {
    108     SkyList *tmp;
    109     tmp = SkyListExistingSubset (skylist, CATDIR);
    110     SkyListFree (skylist);
    111     skylist = tmp;
    112   }
    113   if (VERBOSE) fprintf (stderr, "writing to "OFF_T_FMT" regions\n",  skylist[0].Nregions);
    114 
    115   if (options->only_images) {
    116     fprintf (stderr, "-image (only images) makes no sense with -resort\n");
    117     exit (2);
    118   }
    119   if (options->only_match) {
    120     fprintf (stderr, "-only-match makes no sense with -resort\n");
    121     exit (2);
    122   }
    123   if (options->nosort == 1) {
    124     fprintf (stderr, "-nosort makes no sense with -resort\n");
    125     exit (2);
    126   }
    127   if (options->update) {
    128     fprintf (stderr, "-update makes no sense with -resort\n");
    129     exit (2);
    130   }
    131   if (options->calibrate) {
    132     fprintf (stderr, "-cal (calibrate) makes no sense with -resort\n");
    133     exit (2);
    134   }
    135103
    136104  ALLOCATE(threads, pthread_t, NTHREADS);
    … …  
    148116  Naverage = Nmeasure = 0;
    149117  for (i = 0; i < skylist[0].Nregions; i++) {
     118
     119    // does this host ID match the desired location for the table?
     120    if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue;
    150121
    151122    launched = FALSE;
    … …  
    166137        if (threadData[j].state == TS_FAIL) goto failure;
    167138        if (threadData[j].state != TS_WAIT) continue;
    168         threadData[j].nosort = options->nosort;
     139        threadData[j].forcesort = ForceSort;
    169140        threadData[j].filename = skylist[0].filename[i];
    170141        threadData[j].region = skylist[0].regions[i];
    … …  
    203174  }
    204175
    205   gettimeofday (&stop, NULL);
    206   dtime = DTIME (stop, start);
    207   fprintf (stderr, "SUCCESS: elapsed time %9.4f sec for, "OFF_T_FMT" average, "OFF_T_FMT" measure\n", dtime, Naverage, Nmeasure);
    208 
    209176  // for (i = 0; i < NTHREADS; i++) {
    210177  //   pthread_cancel(threads[i]);

  • trunk/Ohana/src/addstar/src/resort_unthreaded.c

    r35142 r37116  
    11# include "addstar.h"
    22
    3 int resort_unthreaded (AddstarClientOptions *options, SkyTable *sky, int hostID, char *hostpath) {
     3// pass in options and skylist?
     4int resort_unthreaded (SkyList *skylist, int ForceSort) {
    45
    5   double dtime;
    6   struct timeval start, stop;
     6  off_t i;
     7  off_t Naverage, Nmeasure;
     8  Catalog catalog;
    79
    8   gettimeofday (&start, NULL);
     10  /* match stars to existing catalog data (or otherwise manipulate catalog data) */
     11  Naverage = Nmeasure = 0;
     12  for (i = 0; i < skylist[0].Nregions; i++) {
    913
    10   SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch);
    11  
    12   // in these cases, limit the sky catalogs to an existing subset
    13   if (options->existing_regions) {
    14     SkyList *tmp;
    15     tmp = SkyListExistingSubset (skylist, CATDIR);
    16     SkyListFree (skylist);
    17     skylist = tmp;
     14    // does this host ID match the desired location for the table?
     15    if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue;
     16
     17    // chose the catalog file (local or remote?)
     18    char hostfile[1024];
     19    snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name);
     20    catalog.filename  = HOST_ID ? hostfile : skylist[0].filename[i];
     21
     22    // set the parameters which guide catalog open/load/create
     23    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
     24    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
     25    catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_LENSING;
     26    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
     27
     28    // an error exit status here is a significant error (disk I/O or file access)
     29    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
     30      fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename);
     31      exit (2);
     32    }
     33
     34    // Naves_disk == 0 implies an empty catalog file, skip empty catalogs
     35    if (catalog.Naves_disk == 0) {
     36      dvo_catalog_unlock (&catalog);
     37      dvo_catalog_free (&catalog);
     38      continue;
     39    }
     40
     41    // this is an overloaded value to mean 'force sort'
     42    if (ForceSort) catalog.sorted = FALSE;
     43
     44    if (OLD_RESORT) {
     45      resort_catalog_old (&catalog);
     46    } else {
     47      resort_catalog (&catalog);
     48    }
     49
     50    // report total updated values
     51    Naverage += catalog.Naverage;
     52    Nmeasure += catalog.Nmeasure;
     53
     54    // write out catalog, if appropriate
     55    SetProtect (TRUE);
     56    dvo_catalog_save (&catalog, VERBOSE);
     57    SetProtect (FALSE);
     58    dvo_catalog_unlock (&catalog);
     59    dvo_catalog_free (&catalog);
    1860  }
    19   if (VERBOSE) fprintf (stderr, "writing to "OFF_T_FMT" regions\n",  skylist[0].Nregions);
    20 
    21   if (options->only_images) {
    22     fprintf (stderr, "-image (only images) makes no sense with -resort\n");
    23     exit (2);
    24   }
    25   if (options->only_match) {
    26     fprintf (stderr, "-only-match makes no sense with -resort\n");
    27     exit (2);
    28   }
    29   if (options->nosort == 1) {
    30     fprintf (stderr, "-nosort makes no sense with -resort\n");
    31     exit (2);
    32   }
    33   if (options->update) {
    34     fprintf (stderr, "-update makes no sense with -resort\n");
    35     exit (2);
    36   }
    37   if (options->calibrate) {
    38     fprintf (stderr, "-cal (calibrate) makes no sense with -resort\n");
    39     exit (2);
    40   }
    41 
    42   resort_unthreaded_catalogs (options, skylist, hostID, hostpath);
    43 
    44   gettimeofday (&stop, NULL);
    45   dtime = DTIME (stop, start);
    46   fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime);
    47 
    4861  return TRUE;
    4962}
    5063
     64

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

    r35765 r37116  
    129129int        merge_catalogs_old     PROTO((SkyRegion *region, Catalog *output, Catalog *input, double RADIUS, int *secflitMap));
    130130
     131off_t     *build_measure_links    PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
    131132off_t     *init_measure_links     PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
     133int        add_meas_link          PROTO((Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE));
     134Measure   *sort_measure           PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next));
     135
     136off_t     *build_lensing_links    PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
     137off_t     *init_lensing_links     PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
     138int        add_lens_link          PROTO((Average *average, off_t *next, off_t Nlensing, off_t NLENSING));
     139Lensing   *sort_lensing           PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next));
     140
    132141off_t     *init_missing_links     PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing));
    133 off_t     *build_measure_links    PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
    134 int        add_meas_link          PROTO((Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE));
    135142int        add_miss_link          PROTO((Average *average, off_t *next, off_t Nmissing));
    136 Measure   *sort_measure           PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next));
    137143Missing   *sort_missing           PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss));
    138144

  • trunk/Ohana/src/dvomerge/src/LoadCatalog.c

    r29938 r37116  
    88
    99    // always load all of the data (if any exists)
    10     catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     10    catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING | LOAD_LENSOBJ;
    1111
    1212    catalog[0].catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data

  • trunk/Ohana/src/dvomerge/src/build_links.c

    r27435 r37116  
    11# include "dvomerge.h"
     2
     3/*
     4
     5There are two modes for the measure table: sorted and unsorted.
     6
     7In sorted mode, all measures associated with a given average are in a single block.  The block
     8is pointed to by average->measureOffset and the range is average->Nmeasure
     9
     10In unsorted mode, it is not possible to go directly from average to measure without scanning.
     11In this case, it is necessary to use the value measure->averef to find the corresponding
     12average entry. 
     13
     14Note that average->measureOffset and measure->averef are only valid for a given load of the
     15data: they refer to the sequence number in the data blocks.
     16
     17next_meas is a list of the equivalent sequence of the measure block as if it were sorted.
     18
     19to find the sequence of measurements for a given average:
     20n_0 = average->measureOffset
     21n_1 = next_meas[n_0]
     22n_i = next_meas[n_i-1]
     23
     24*/
    225
    326/* build the initial links assuming the table is sorted,
    … …  
    629
    730  off_t i, j, N;
    8   off_t *next;
     31  off_t *next_meas;
    932
    1033  N = 0;
    1134
    12   ALLOCATE (next, off_t, Nmeasure);
     35  ALLOCATE (next_meas, off_t, Nmeasure);
    1336  for (i = 0; i < Naverage; i++, N++) {
    1437    for (j = 0; j < average[i].Nmeasure - 1; j++, N++) {
    15       next[N] = N + 1;
     38      next_meas[N] = N + 1;
    1639      if (N >= Nmeasure) {
    1740        fprintf (stderr, "WARNING: N out of bounds (1)\n");
    1841      }
    1942    }
    20     next[N] = -1;
     43    next_meas[N] = -1;
    2144    if (N >= Nmeasure) {
    2245      fprintf (stderr, "WARNING: N out of bounds (2)\n");
    … …  
    2851    }
    2952  }
    30   return (next);
    31 }
    32 
    33 /* average[].measureOffset, average[].Nmeasure are valid within an addstar run */
    34 int add_meas_link (Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE) {
    35 
    36   off_t k, m;
    37 
    38   /* if we have trouble, check validity of next[m] : m < Nmeasure */
    39   m = average[0].measureOffset; 
    40 
    41   for (k = 0; k < average[0].Nmeasure - 1; k++)  {
    42     m = next[m];
    43     if (m >= NMEASURE) {
    44       fprintf (stderr, "WARNING: m out of bounds (3)\n");
    45     }
    46   }
    47 
    48   /* set up references */
    49   next[Nmeasure] = -1;
    50   if (Nmeasure >= NMEASURE) {
    51     fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n");
    52   }
    53 
    54   if (m == -1) {
    55     average[0].measureOffset = Nmeasure;
    56   } else {
    57     next[m] = Nmeasure;
    58     if (m >= NMEASURE) {
    59       fprintf (stderr, "WARNING: m out of bounds (4)\n");
    60     }
    61   }
    62 
    63   return (TRUE);
    64 }
    65 
    66 Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next) {
    67 
    68   off_t i, k, n, N;
    69   Measure *tmpmeasure;
    70 
    71   /* fix order of Measure (memory intensive, but fast) */
    72   N = 0;
    73   ALLOCATE (tmpmeasure, Measure, Nmeasure);
    74   for (i = 0; i < Naverage; i++) {
    75     n = average[i].measureOffset;
    76     average[i].measureOffset = N;
    77     for (k = 0; k < average[i].Nmeasure; k++, N++) {
    78       tmpmeasure[N] = measure[n];
    79       tmpmeasure[N].averef = i;
    80       n = next[n];
    81     }
    82   }
    83   free (measure);
    84   return (tmpmeasure);
     53  return (next_meas);
    8554}
    8655
    … …  
    9564
    9665  off_t i, m, k, Nm, averef;
    97   off_t *next;
    98 
    99   ALLOCATE (next, off_t, Nmeasure);
     66  off_t *next_meas;
     67
     68  ALLOCATE (next_meas, off_t, Nmeasure);
    10069
    10170  /* reset the Nm, offset values for average */
    … …  
    10877    averef = measure[Nm].averef;
    10978    m = average[averef].measureOffset; 
    110     next[Nm] = -1;
     79    next_meas[Nm] = -1;
    11180
    11281    if (m == -1) { /* no links yet for source */
    … …  
    11685    }
    11786
    118     for (k = 0; next[m] != -1; k++) {
    119       m = next[m];
     87    for (k = 0; next_meas[m] != -1; k++) {
     88      m = next_meas[m];
    12089      if (m >= Nmeasure) {
    12190        fprintf (stderr, "WARNING: m out of bounds (1)\n");
    … …  
    12493
    12594    average[averef].Nmeasure = k + 2;
    126     next[m] = Nm;
     95    next_meas[m] = Nm;
    12796    if (m >= Nmeasure) {
    12897      fprintf (stderr, "WARNING: m out of bounds (2)\n");
    12998    }
    13099  }
    131   return (next);
    132 }
     100  return (next_meas);
     101}
     102
     103/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */
     104int add_meas_link (Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE) {
     105
     106  off_t k, m;
     107
     108  /* if we have trouble, check validity of next_meas[m] : m < Nmeasure */
     109  m = average[0].measureOffset; 
     110
     111  for (k = 0; k < average[0].Nmeasure - 1; k++)  {
     112    m = next_meas[m];
     113    if (m >= NMEASURE) {
     114      fprintf (stderr, "WARNING: m out of bounds (3)\n");
     115    }
     116  }
     117
     118  /* set up references */
     119  next_meas[Nmeasure] = -1;
     120  if (Nmeasure >= NMEASURE) {
     121    fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n");
     122  }
     123
     124  if (m == -1) {
     125    average[0].measureOffset = Nmeasure;
     126  } else {
     127    next_meas[m] = Nmeasure;
     128    if (m >= NMEASURE) {
     129      fprintf (stderr, "WARNING: m out of bounds (4)\n");
     130    }
     131  }
     132
     133  return (TRUE);
     134}
     135
     136Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next_meas) {
     137
     138  off_t i, k, n, N;
     139  Measure *tmpmeasure;
     140
     141  /* fix order of Measure (memory intensive, but fast) */
     142  N = 0;
     143  ALLOCATE (tmpmeasure, Measure, Nmeasure);
     144  for (i = 0; i < Naverage; i++) {
     145    n = average[i].measureOffset;
     146    average[i].measureOffset = N;
     147    for (k = 0; k < average[i].Nmeasure; k++, N++) {
     148      if (n == -1) abort();
     149      tmpmeasure[N] = measure[n];
     150      if (measure[n].averef != i) abort();
     151      tmpmeasure[N].averef = i;
     152      n = next_meas[n];
     153    }
     154  }
     155  free (measure);
     156  return (tmpmeasure);
     157}
     158
     159/*******************************************************************************************/
    133160
    134161/* build the initial links assuming the table is sorted */
    … …  
    136163
    137164  off_t i, j, N;
    138   off_t *next;
     165  off_t *next_miss;
    139166
    140167  N = 0;
    141168
    142   ALLOCATE (next, off_t, Nmissing);
     169  ALLOCATE (next_miss, off_t, Nmissing);
    143170  for (i = 0; i < Naverage; i++) {
    144171    for (j = 0; j < average[i].Nmissing - 1; j++, N++) {
    145       next[N] = N + 1;
     172      next_miss[N] = N + 1;
    146173    }
    147174    if (average[i].Nmissing > 0) {
    148       next[N] = -1;
     175      next_miss[N] = -1;
    149176      if (N >= Nmissing) {
    150177        fprintf (stderr, "overflow in init_missing_links");
    … …  
    155182
    156183  }
    157   return (next);
    158 }
    159 
    160 int add_miss_link (Average *average, off_t *next, off_t Nmissing) {
     184  return (next_miss);
     185}
     186
     187int add_miss_link (Average *average, off_t *next_miss, off_t Nmissing) {
    161188
    162189  off_t k, m;
    … …  
    165192  if (average[0].Nmissing < 1) {
    166193    average[0].missingOffset = Nmissing;
    167     next[Nmissing] = -1;
     194    next_miss[Nmissing] = -1;
    168195    return (TRUE);
    169196  }
    170197
    171198  m = average[0].missingOffset; 
    172   for (k = 0; k < average[0].Nmissing - 1; k++) m = next[m];
     199  for (k = 0; k < average[0].Nmissing - 1; k++) m = next_miss[m];
    173200  /* set up references */
    174   next[Nmissing] = -1;
    175   next[m] = Nmissing;
     201  next_miss[Nmissing] = -1;
     202  next_miss[m] = Nmissing;
    176203  return (TRUE);
    177204}
    … …  
    180207   we must always save the missing table, if it exists */
    181208
    182 Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next) {
     209Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss) {
    183210
    184211  off_t i, k, n, N;
    … …  
    193220    for (k = 0; k < average[i].Nmissing; k++, N++) {
    194221      tmpmissing[N] = missing[n];
    195       n = next[n];
     222      n = next_miss[n];
    196223    }
    197224  }
    … …  
    199226  return (tmpmissing);
    200227}
     228
     229/*******************************************************************************************/
     230
     231/* build the initial links assuming the table is sorted,
     232   not partial, and has a correct set of average[].lensingOffset,Nlensing values */
     233off_t *init_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) {
     234
     235  off_t i, j, N;
     236  off_t *next_lens;
     237
     238  N = 0;
     239
     240  ALLOCATE (next_lens, off_t, Nlensing);
     241  for (i = 0; i < Naverage; i++) {
     242    if (!average[i].Nlensing) continue;
     243    for (j = 0; j < average[i].Nlensing - 1; j++, N++) {
     244      next_lens[N] = N + 1;
     245      if (N >= Nlensing) {
     246        fprintf (stderr, "WARNING: N out of bounds (1)\n");
     247      }
     248    }
     249    next_lens[N] = -1;
     250    if (N >= Nlensing) {
     251      fprintf (stderr, "WARNING: N out of bounds (2)\n");
     252    }
     253
     254    if (N >= Nlensing) {
     255      fprintf (stderr, "overflow in init_lensing_links\n");
     256      abort ();
     257    }
     258    N++;
     259  }
     260  return (next_lens);
     261}
     262
     263/* construct lensing links which are valid FOR THIS LOAD
     264 * - if we have a full load, we will get links which can
     265 *   be used by other programs (eg, relphot, etc)
     266 * - if we have a partial load, the links are only valid
     267 *   for that partial load
     268 */
     269
     270off_t *build_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) {
     271
     272  off_t i, m, k, Nm, averef;
     273  off_t *next_lens;
     274
     275  ALLOCATE (next_lens, off_t, Nlensing);
     276
     277  /* reset the Nm, offset values for average */
     278  for (i = 0; i < Naverage; i++) {
     279    average[i].lensingOffset = -1;
     280    average[i].Nlensing     =  0;
     281  }
     282
     283  for (Nm = 0; Nm < Nlensing; Nm++) {
     284    averef = lensing[Nm].averef;
     285    m = average[averef].lensingOffset; 
     286    next_lens[Nm] = -1;
     287
     288    if (m == -1) { /* no links yet for source */
     289      average[averef].lensingOffset = Nm;
     290      average[averef].Nlensing     = 1;
     291      continue;
     292    }
     293
     294    for (k = 0; next_lens[m] != -1; k++) {
     295      m = next_lens[m];
     296      if (m >= Nlensing) {
     297        fprintf (stderr, "WARNING: m out of bounds (1)\n");
     298      }
     299    }
     300
     301    average[averef].Nlensing = k + 2;
     302    next_lens[m] = Nm;
     303    if (m >= Nlensing) {
     304      fprintf (stderr, "WARNING: m out of bounds (2)\n");
     305    }
     306  }
     307  return (next_lens);
     308}
     309
     310/* average[].lensingOffset, average[].Nlensing are valid within an addstar run */
     311int add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) {
     312
     313  off_t k, m;
     314
     315  /* if we have trouble, check validity of next_lens[m] : m < Nlensing */
     316  m = average[0].lensingOffset; 
     317
     318  for (k = 0; k < average[0].Nlensing - 1; k++)  {
     319    m = next_lens[m];
     320    if (m >= NLENSING) {
     321      fprintf (stderr, "WARNING: m out of bounds (3)\n");
     322    }
     323  }
     324
     325  /* set up references */
     326  next_lens[Nlensing] = -1;
     327  if (Nlensing >= NLENSING) {
     328    fprintf (stderr, "WARNING: Nlensing out of bounds (1)\n");
     329  }
     330
     331  if (m == -1) {
     332    average[0].lensingOffset = Nlensing;
     333  } else {
     334    next_lens[m] = Nlensing;
     335    if (m >= NLENSING) {
     336      fprintf (stderr, "WARNING: m out of bounds (4)\n");
     337    }
     338  }
     339
     340  return (TRUE);
     341}
     342
     343Lensing *sort_lensing (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next_lens) {
     344
     345  off_t i, k, n, np, N;
     346  Lensing *tmplensing;
     347
     348  /*
     349  for (i = 0; i < Naverage; i++) {
     350    if (average[i].Nlensing != 4) {
     351      fprintf (stderr, "check %d %d %d\n", (int) i, (int) average[i].Nlensing, (int) average[i].lensingOffset);
     352    }
     353  }
     354  */
     355
     356  /* fix order of Lensing (memory intensive, but fast) */
     357  N = 0;
     358  ALLOCATE (tmplensing, Lensing, Nlensing);
     359  for (i = 0; i < Naverage; i++) {
     360    if (!average[i].Nlensing) continue;
     361    n = average[i].lensingOffset;
     362    average[i].lensingOffset = N;
     363    for (k = 0; k < average[i].Nlensing; k++, N++) {
     364      if (n == -1) {
     365        fprintf (stderr, "entry after %d has a problem\n", (int) np);
     366        abort();
     367      }
     368      tmplensing[N] = lensing[n];
     369      if (lensing[n].averef != i) abort();
     370      tmplensing[N].averef = i;
     371      np = n;
     372      n = next_lens[n];
     373    }
     374  }
     375  free (lensing);
     376  return (tmplensing);
     377}
     378
     379

  • trunk/Ohana/src/dvomerge/src/dvomergeContinue.c

    r34260 r37116  
    131131      merge_catalogs_old (&outsky[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap);
    132132
    133       outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    134 
    135133      // if we receive a signal which would cause us to exit, wait until the full catalog is written
    136134      SetProtect (TRUE);

  • trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c

    r34260 r37116  
    7474      merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap);
    7575
    76       outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    77 
    7876      // if we receive a signal which would cause us to exit, wait until the full catalog is written
    7977      SetProtect (TRUE);

  • trunk/Ohana/src/dvomerge/src/dvomergeCreate.c

    r33963 r37116  
    178178    SkyListFree (inlist);
    179179
    180     outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    181180    dvo_catalog_save (&outcatalog, VERBOSE);
    182181    dvo_catalog_unlock (&outcatalog);

  • trunk/Ohana/src/dvomerge/src/dvomergeFromList.c

    r35590 r37116  
    143143    merge_catalogs_old (&skyregion, &outcatalog, &incatalog, RADIUS, secfiltMap);
    144144   
    145     outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    146 
    147145    dmhObjectAdd (outstat[0].history, &outcatalog.header, inStats);
    148146

  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c

    r36528 r37116  
    1313  int NsecfiltInput, NsecfiltOutput;
    1414
    15   double dtime;
    16   struct timeval start, stop;
    17   gettimeofday (&start, NULL);
     15  INITTIME;
    1816
    1917  CONTINUE = FALSE;
    … …  
    121119  }
    122120
    123   gettimeofday (&stop, NULL);
    124   dtime = DTIME (stop, start);
    125 
    126121  if (!status) {
    127     fprintf (stderr, "ERROR: elapsed time %9.4f sec\n", dtime);
     122    MARKTIME ("ERROR: elapsed time %9.4f sec\n", dtime);
    128123    exit (3);
    129124  }
    130125
    131   fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime);
     126  MARKTIME ("SUCCESS: elapsed time %9.4f sec\n", dtime);
    132127  exit (0);
    133128}

  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c

    r37046 r37116  
    163163      LoadCatalog (&outcatalog, outlist[0].regions[j], outcatalog.filename, "w", NsecfiltOutput);
    164164
     165      // if no catalog already exists, use the input catalog to define the format
     166      if (outcatalog.Naves_disk == 0) {
     167        outcatalog.catformat = incatalog.catformat;
     168      }
     169
    165170      merge_catalogs_old (outlist[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap);
    166 
    167       outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    168171
    169172      if (outstat[j].missed) {

  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c

    r35765 r37116  
    7272      merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap);
    7373
    74       outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    75 
    7674      // if we receive a signal which would cause us to exit, wait until the full catalog is written
    7775      SetProtect (TRUE);
    … …  
    115113  int NsecfiltInput, NsecfiltOutput;
    116114
    117   double dtime;
    118   struct timeval start, stop;
    119   gettimeofday (&start, NULL);
     115  INITTIME;
    120116
    121117  CONTINUE = FALSE;
    … …  
    134130
    135131  fprintf (stderr, "WARNING / ERROR : multi-threaded dvomerge does not handle merge tracking yet\n");
     132  fprintf (stderr, "WARNING / ERROR : multi-threaded dvomerge does work in a parallel context\n");
    136133  exit (2);
    137134
    … …  
    293290  }
    294291
    295   gettimeofday (&stop, NULL);
    296   dtime = DTIME (stop, start);
    297   fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime);
    298 
     292  MARKTIME ("SUCCESS: elapsed time %9.4f sec\n", dtime);
    299293  exit (0);
    300294
    … …  
    318312    }
    319313  }
     314  MARKTIME ("ERROR: elapsed time %9.4f sec\n", dtime);
    320315  exit (1);
    321316}

  • trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c

    r28855 r37116  
    88// input entries always define new objects
    99
    10 #define notyet 1
    11 
    1210int merge_catalogs_new (SkyRegion *region, Catalog *output, Catalog *input, int *secfiltMap) {
    1311 
    1412  off_t i, j, offset;
    15   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, NsecfiltIn, NsecfiltOut, Nm;
     13  off_t NAVERAGE, NMEASURE, NLENSING, Naverage, Nmeasure, Nlensing, NsecfiltIn, NsecfiltOut, Nm;
    1614
    1715  Naverage = output[0].Naverage;
    … …  
    6361      }
    6462      output[0].average[Naverage].Nmeasure = Nm;
     63
     64      Nm = 0;
     65      for (j = 0; j < input[0].average[i].Nlensing; j++) {
     66          offset = input[0].average[i].lensingOffset + j;
     67
     68          output[0].lensing[Nlensing] = input[0].lensing[offset];
     69          output[0].lensing[Nlensing].averef = Naverage;
     70
     71          Nlensing ++;
     72          Nm ++;
     73          if (Nlensing == NLENSING) {
     74              NLENSING += 1000;
     75              REALLOCATE (output[0].lensing, Lensing, NLENSING);
     76          }
     77      }
     78      output[0].average[Naverage].Nlensing = Nm;
     79
    6580      Naverage ++;
    6681      if (Naverage == NAVERAGE) {
    … …  
    7287  REALLOCATE (output[0].average, Average, MAX (Naverage, 1));
    7388  REALLOCATE (output[0].measure, Measure, MAX (Nmeasure, 1));
     89  REALLOCATE (output[0].lensing, Lensing, MAX (Nlensing, 1));
    7490  REALLOCATE (output[0].secfilt, SecFilt, NsecfiltOut*MAX (Naverage, 1));
    7591  output[0].Naverage = Naverage;
    76   output[0].Nmeasure = Nmeasure;
     92  output[0].Nlensing = Nlensing;
    7793  output[0].Nsecf_mem = Naverage * NsecfiltOut;
    7894
    … …  
    84100 
    85101  if (VERBOSE) {
    86       fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures) for catalog\n",
     102      fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures, "OFF_T_FMT" lensing) for catalog\n",
    87103                i,
    88104                output[0].Naverage,
    89                 output[0].Nmeasure);
     105                output[0].Nmeasure,
     106                output[0].Nlensing);
    90107  }
    91108  return (TRUE);

  • trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c

    r37046 r37116  
    1414  double *X1, *Y1, *X2, *Y2;
    1515  double dX, dY, dR;
    16   off_t *N1, *N2,  *next_meas;
    17   off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
     16  off_t *N1, *N2, *next_meas, *next_lens;
     17  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS;
    1818  int NsecfiltIn;
    1919  int NsecfiltOut;
    … …  
    5252  Nmatch = 0;
    5353  NMEAS = Nmeas = output[0].Nmeasure;
     54  NLENS = Nlens = output[0].Nlensing;
    5455
    5556  // current max obj ID for this catalog
    … …  
    113114  }
    114115
    115 
    116116  /* set up pointers for linked list of measure */
    117117  if (output[0].sorted && (output[0].Nmeasure >= output[0].Nmeas_disk)) {
    … …  
    119119    // is sorted while processed
    120120    next_meas = init_measure_links (output[0].average, Nave, output[0].measure, Nmeas);
     121    next_lens = init_lensing_links (output[0].average, Nave, output[0].lensing, Nlens);
    121122  } else {
    122123    next_meas = build_measure_links (output[0].average, Nave, output[0].measure, Nmeas);
     124    next_lens = build_lensing_links (output[0].average, Nave, output[0].lensing, Nlens);
    123125  }   
    124126
    … …  
    182184      REALLOCATE (output[0].measure, Measure, NMEAS);
    183185    }
     186    if (Nlens + input[0].average[N].Nmeasure >= NLENS) {
     187      NLENS = Nlens + input[0].average[N].Nmeasure + 1000;
     188      REALLOCATE (next_lens, off_t, NLENS);
     189      REALLOCATE (output[0].lensing, Lensing, NLENS);
     190    }
    184191
    185192    // 4) average properties from the input and the output db need to be properly merged.
    … …  
    191198      if (REPLACE_BY_PHOTCODE) {
    192199        // index to first measure for this object
     200        // XXX this does not support lensing measurements
    193201        int Mout = output[0].average[n].measureOffset; 
    194202        if (replace_match(&output[0].average[n], &output[0].measure[Mout], &input[0].average[N], &input[0].measure[offset])) {
    … …  
    248256    }
    249257
     258    // if lensing measurements exist, add them too
     259    if (output[0].lensing) {
     260      for (Nin = 0; Nin < input[0].average[N].Nlensing; Nin++) {
     261        /* add to end of lensing list */
     262        add_lens_link (&output[0].average[n], next_lens, Nlens, NLENS);
     263       
     264        // set the new lensing
     265        off_t lensoff = input[0].average[N].lensingOffset + Nin;
     266        output[0].lensing[Nlens] = input[0].lensing[lensoff];
     267
     268        output[0].lensing[Nlens].averef   = n;
     269        output[0].lensing[Nlens].objID    = output[0].average[n].objID;
     270        output[0].lensing[Nlens].catID    = output[0].catID;
     271        output[0].average[n].Nlensing ++;
     272        Nlens ++;
     273      }
     274    }
     275
    250276    // update the average properties to reflect the incoming entries:
    251277    // if the original value is NAN but the input value is not, accept the input:
    … …  
    288314      REALLOCATE (output[0].measure, Measure, NMEAS);
    289315    }
     316    if (Nlens + input[0].average[N].Nlensing >= NLENS) {
     317      NLENS = Nlens + input[0].average[N].Nlensing + 1000;
     318      REALLOCATE (next_lens, off_t, NLENS);
     319      REALLOCATE (output[0].lensing, Lensing, NLENS);
     320    }
    290321    if (Nave >= NAVE) {
    291322      NAVE = Nave + 1000;
    … …  
    355386      next_meas[Nmeas - Ngroup + j] = Nmeas - Ngroup + j + 1;
    356387    }
     388
     389    /** add lensing for this input average object **/
     390    if (output[0].lensing) {
     391      output[0].average[Nave].lensingOffset  = Nlens;
     392      for (Nin = 0; Nin < input[0].average[N].Nlensing; Nin ++) {
     393        // supply the lensing values from this detection
     394        off_t lensoff = input[0].average[N].lensingOffset + Nin;
     395        output[0].lensing[Nlens]           = input[0].lensing[lensoff];
     396
     397        // the following lensing elements cannot be set until here:
     398        output[0].lensing[Nlens].averef   = Nave;
     399        output[0].lensing[Nlens].objID    = output[0].average[Nave].objID;
     400        output[0].lensing[Nlens].catID    = output[0].catID;
     401
     402        // as we add lensing, update Nlensing to match
     403        output[0].average[Nave].Nlensing ++;
     404
     405        /* we set next[Nlens] to -1 here, and update correctly below */
     406        next_lens[Nlens] = -1;
     407        Nlens ++;
     408      }
     409      int Ngroup = input[0].average[N].Nlensing;
     410      for (j = 0; j < Ngroup - 1; j++) {
     411        next_lens[Nlens - Ngroup + j] = Nlens - Ngroup + j + 1;
     412      }
     413    }
     414
    357415    Nave ++;
    358416  }
    … …  
    362420  REALLOCATE (output[0].average, Average, Nave);
    363421  REALLOCATE (output[0].measure, Measure, Nmeas);
     422  REALLOCATE (output[0].lensing, Lensing, Nlens);
    364423 
    365424# define NOSORT 0
    … …  
    369428    output[0].sorted = TRUE;
    370429    output[0].measure = sort_measure (output[0].average, Nave, output[0].measure, Nmeas, next_meas);
     430    output[0].lensing = sort_lensing (output[0].average, Nave, output[0].lensing, Nlens, next_lens);
    371431  }
    372432
    … …  
    384444  output[0].Naverage = Nave;
    385445  output[0].Nmeasure = Nmeas;
     446  output[0].Nlensing = Nlens;
    386447  output[0].Nsecf_mem = Nave*NsecfiltOut;
    387   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);
     448  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nlens, Nmatch);
    388449
    389450  free (next_meas);
     451  free (next_lens);
    390452
    391453  free (X2);

  • trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c

    r37048 r37116  
    352352    PrintIOBuffer (sec_buffer, "%hd, ", secfilt->ubercalDist); 
    353353    PrintIOBuffer (sec_buffer, "%.6f, ", secfilt->Mstdev);       
    354     PrintIOBuffer (sec_buffer, "%lu  ", secfilt->stackDetectID);
     354    PrintIOBuffer (sec_buffer, "%lu  ", 0L); // secfilt->stackDetectID
    355355    PrintIOBuffer (sec_buffer, "),\n");
    356356    secfilt ++;

  • trunk/Ohana/src/getstar

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • trunk/Ohana/src/libautocode/def/secfilt-ps1-v4.d

    r37034 r37116  
    2323FIELD  ubercalDist,   UBERCAL_DIST,    short,      number of images from an ubercal-image
    2424FIELD  Mstdev,        MAG_STDEV,       short,      standard deviation of measurements,    millimags
    25 FIELD  stackDetectID, STACK_DETECT_ID, uint64_t,   detection ID on stack used for fluxes (if any)
     25FIELD  stackPrmryOff, STACK_PRIMARY_OFF, int,      measure entry which is primary stack detection
     26FIELD  stackBestOff,  STACK_BEST_OFF,    int,      measure entry which is best stack detection

  • trunk/Ohana/src/libautocode/def/secfilt-ps1-v5.d

    r37034 r37116  
    4242FIELD  dFapStk,       FLUX_AP_STK_ERR,   float,    mean flux ap error
    4343
    44 ## need an 8-byte boundary here ##
    45 FIELD  stackDetectID, STACK_DETECT_ID,   uint64_t, detection ID on stack used for fluxes (if any)
     44FIELD  stackPrmryOff, STACK_PRIMARY_OFF,  int,     measure entry which is primary stack detection
     45FIELD  stackBestOff,  STACK_BEST_OFF,     int,     measure entry which is best stack detection
    4646
    4747## *** this section is for forced-warp mean values ***

  • trunk/Ohana/src/libautocode/def/secfilt.d

    r37034 r37116  
    4242FIELD  dFapStk,       FLUX_AP_STK_ERR,   float,    mean flux ap error
    4343
    44 ## need an 8-byte boundary here ##
    45 FIELD  stackDetectID, STACK_DETECT_ID, uint64_t,   detection ID on stack used for fluxes (if any)
     44FIELD  stackPrmryOff, STACK_PRIMARY_OFF,  int,     measure entry which is primary stack detection
     45FIELD  stackBestOff,  STACK_BEST_OFF,     int,     measure entry which is best stack detection
    4646
    4747## *** this section is for forced-warp mean values ***

  • trunk/Ohana/src/libdvo/include/dvo.h

    r37035 r37116  
    128128  ID_MEAS_BLEND_MEAS     = 0x00000100,  // detection is within radius of multiple objects
    129129  ID_MEAS_BLEND_OBJ      = 0x00000200,  // multiple detections within radius of object
    130   ID_MEAS_UNDEF_3        = 0x00000400,  // unused
     130  ID_MEAS_WARP_USED      = 0x00000400,  // measurement used to find mean warp photometry
    131131  ID_MEAS_UNDEF_4        = 0x00000800,  // unused
    132132  ID_MEAS_BLEND_MEAS_X   = 0x00001000,  // detection is within radius of multiple objects across catalogs                   
    … …  
    383383  char  ***name; // projection cell name
    384384 
    385   int NX_SUB;
    386   int NY_SUB;
     385  float NX_SUB;
     386  float NY_SUB;
    387387  double dPix;
    388388
    389389  double **Xo;
    390390  double **Yo;
    391   int **dX;
    392   int **dY;
     391  float **dX;
     392  float **dY;
    393393} BoundaryTree;
    394394
    … …  
    411411  double Do;
    412412  double dPix;
    413   int dX;
    414   int dY;
    415 
    416   int NX_SUB;
    417   int NY_SUB;
     413  float dX;
     414  float dY;
     415
     416  float NX_SUB;
     417  float NY_SUB;
    418418
    419419  char *basename;
    … …  
    710710float PhotFluxInstErr (Measure *measure, dvoMagClassType class);
    711711float PhotFluxCatErr (Measure *measure, dvoMagClassType class);
     712float PhotFluxAveErr (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source);
    712713
    713714float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt);
    714715float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt);
    715716
    716 float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt);
     717int   PhotNphot (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source);
     718float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source);
    717719float PhotMmin (PhotCode *code, Average *average, SecFilt *secfilt);
    718720float PhotMmax (PhotCode *code, Average *average, SecFilt *secfilt);

  • trunk/Ohana/src/libdvo/include/dvodb.h

    r37035 r37116  
    5454  MAG_OPTION_NPHOT, // Nused
    5555  MAG_OPTION_UC_DIST,
    56   MAG_OPTION_STACK_DET_ID,
    5756  MAG_OPTION_FLAGS,
    5857} dvoMagOptionType;
     58
     59//  MAG_OPTION_STACK_PRIMARY_OFF,
     60//  MAG_OPTION_STACK_BEST_OFF,
    5961
    6062typedef enum {

  • trunk/Ohana/src/libdvo/src/BoundaryTree.c

    r36490 r37116  
    8484  gfits_scan (headerPHU, "DEC_ORI",  "%lf", 1, &tree->DEC_origin);
    8585  gfits_scan (headerPHU, "DEC_OFF",  "%lf", 1, &tree->DEC_offset);
    86   gfits_scan (headerPHU, "NX_SUB",   "%d",  1, &tree->NX_SUB);
    87   gfits_scan (headerPHU, "NY_SUB",   "%d",  1, &tree->NY_SUB);
     86  gfits_scan (headerPHU, "NX_SUB",   "%f",  1, &tree->NX_SUB);
     87  gfits_scan (headerPHU, "NY_SUB",   "%f",  1, &tree->NY_SUB);
    8888  gfits_scan (headerPHU, "PIXSCALE", "%lf", 1, &tree->dPix);
    8989
    … …  
    111111  ALLOCATE (tree->Xo,   double *, tree->Nzone);
    112112  ALLOCATE (tree->Yo,   double *, tree->Nzone);
    113   ALLOCATE (tree->dX,      int *, tree->Nzone);
    114   ALLOCATE (tree->dY,      int *, tree->Nzone);
     113  ALLOCATE (tree->dX,    float *, tree->Nzone);
     114  ALLOCATE (tree->dY,    float *, tree->Nzone);
    115115  ALLOCATE (tree->cell,    int *, tree->Nzone);
    116116  ALLOCATE (tree->projID,  int *, tree->Nzone);
    … …  
    121121    ALLOCATE (tree->Xo[i],   double, tree->Nband[i]);
    122122    ALLOCATE (tree->Yo[i],   double, tree->Nband[i]);
    123     ALLOCATE (tree->dX[i],      int, tree->Nband[i]);
    124     ALLOCATE (tree->dY[i],      int, tree->Nband[i]);
     123    ALLOCATE (tree->dX[i],    float, tree->Nband[i]);
     124    ALLOCATE (tree->dY[i],    float, tree->Nband[i]);
    125125    ALLOCATE (tree->cell[i],    int, tree->Nband[i]);
    126126    ALLOCATE (tree->projID[i],  int, tree->Nband[i]);
    … …  
    147147  GET_COLUMN_NEW(&headerCell, &ftableCell, Xo,    "X_CENT",      double);
    148148  GET_COLUMN_NEW(&headerCell, &ftableCell, Yo,    "Y_CENT",      double);
    149   GET_COLUMN_NEW(&headerCell, &ftableCell, dX,    "X_GRID",      int);
    150   GET_COLUMN_NEW(&headerCell, &ftableCell, dY,    "Y_GRID",      int);
     149  GET_COLUMN_NEW(&headerCell, &ftableCell, dX,    "X_GRID",      float);
     150  GET_COLUMN_NEW(&headerCell, &ftableCell, dY,    "Y_GRID",      float);
    151151  GET_COLUMN_NEW(&headerCell, &ftableCell, name,  "NAME",        char); // XXX how is this done?
    152152  gfits_free_header (&headerCell);
    … …  
    214214  gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset);
    215215
    216   gfits_modify (&header, "NX_SUB", "%d", 1, tree->NX_SUB);
    217   gfits_modify (&header, "NY_SUB", "%d", 1, tree->NY_SUB);
     216  gfits_modify (&header, "NX_SUB", "%f", 1, tree->NX_SUB);
     217  gfits_modify (&header, "NY_SUB", "%f", 1, tree->NY_SUB);
    218218  gfits_modify (&header, "PIXSCALE", "%lf", 1, tree->dPix);
    219219
    … …  
    289289    gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0);
    290290    gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0);
    291     gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0);
    292     gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0);
     291    gfits_define_bintable_column (&theader, "E", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0);
     292    gfits_define_bintable_column (&theader, "E", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0);
    293293    gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0);
    294294
    … …  
    310310    double *Xo            ; ALLOCATE (Xo,    double, Ncell);
    311311    double *Yo            ; ALLOCATE (Yo,    double, Ncell);
    312     int    *dX            ; ALLOCATE (dX,    int,    Ncell);
    313     int    *dY            ; ALLOCATE (dY,    int,    Ncell);
     312    float  *dX            ; ALLOCATE (dX,    float,  Ncell);
     313    float  *dY            ; ALLOCATE (dY,    float,  Ncell);
    314314    char   *name          ; ALLOCATE (name,  char,   Ncell*BOUNDARY_TREE_NAME_LENGTH);
    315315

  • trunk/Ohana/src/libdvo/src/TessellationTable.c

    r36490 r37116  
    2828    tess[i].dPix = NAN;
    2929
    30     tess[i].dX = 0;
    31     tess[i].dY = 0;
    32 
    33     tess[i].NX_SUB = 0;
    34     tess[i].NY_SUB = 0;
     30    tess[i].dX = NAN;
     31    tess[i].dY = NAN;
     32
     33    tess[i].NX_SUB = NAN;
     34    tess[i].NY_SUB = NAN;
    3535
    3636    tess[i].basename = NULL;
    … …  
    123123  GET_COLUMN_NEW(Xo,     "X_CENT", double);
    124124  GET_COLUMN_NEW(Yo,     "Y_CENT", double);
    125   GET_COLUMN_NEW(dX,     "X_GRID", int);
    126   GET_COLUMN_NEW(dY,     "Y_GRID", int);
     125  GET_COLUMN_NEW(dX,     "X_GRID", float);
     126  GET_COLUMN_NEW(dY,     "Y_GRID", float);
    127127  GET_COLUMN_NEW(Rmin,   "R_MIN",  double);
    128128  GET_COLUMN_NEW(Rmax,   "R_MAX",  double);
    … …  
    130130  GET_COLUMN_NEW(Dmax,   "D_MAX",  double);
    131131  GET_COLUMN_NEW(dPix,   "SCALE",  double);
    132   GET_COLUMN_NEW(NX_SUB, "NX_SUB", int);
    133   GET_COLUMN_NEW(NY_SUB, "NY_SUB", int);
     132  GET_COLUMN_NEW(NX_SUB, "NX_SUB", float);
     133  GET_COLUMN_NEW(NY_SUB, "NY_SUB", float);
    134134  GET_COLUMN_NEW(TYPE,   "TYPE",   int);
    135135
    … …  
    280280    gfits_modify (&header, "DEC_OFF", "%lf", 1, rings->DEC_offset);
    281281
    282     gfits_modify (&header, "NX_SUB", "%d", 1, rings->NX_SUB);
    283     gfits_modify (&header, "NY_SUB", "%d", 1, rings->NY_SUB);
     282    gfits_modify (&header, "NX_SUB", "%f", 1, rings->NX_SUB);
     283    gfits_modify (&header, "NY_SUB", "%f", 1, rings->NY_SUB);
    284284    gfits_modify (&header, "PIXSCALE", "%lf", 1, rings->dPix);
    285285  }
    … …  
    300300  gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0);
    301301  gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0);
    302   gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0);
    303   gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0);
     302  gfits_define_bintable_column (&theader, "E", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0);
     303  gfits_define_bintable_column (&theader, "E", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0);
    304304  gfits_define_bintable_column (&theader, "D", "R_MIN",  "RA limit (lower)", "none", 1.0, 0.0);
    305305  gfits_define_bintable_column (&theader, "D", "R_MAX",  "RA limit (upper)", "none", 1.0, 0.0);
    … …  
    307307  gfits_define_bintable_column (&theader, "D", "D_MAX",  "DEC limit (upper)", "none", 1.0, 0.0);
    308308  gfits_define_bintable_column (&theader, "D", "SCALE",  "pixel scale for projection cell", "none", 1.0, 0.0);
    309   gfits_define_bintable_column (&theader, "J", "NX_SUB", "skycell subdivision in x", "none", 1.0, 0.0);
    310   gfits_define_bintable_column (&theader, "J", "NY_SUB", "skycell subdivision in y", "none", 1.0, 0.0);
     309  gfits_define_bintable_column (&theader, "E", "NX_SUB", "skycell subdivision in x", "none", 1.0, 0.0);
     310  gfits_define_bintable_column (&theader, "E", "NY_SUB", "skycell subdivision in y", "none", 1.0, 0.0);
    311311  gfits_define_bintable_column (&theader, "J", "TYPE",   "type of tessellation", "none", 1.0, 0.0);
    312312
    … …  
    330330  double *Xo            ; ALLOCATE (Xo    ,    double, Nout);
    331331  double *Yo            ; ALLOCATE (Yo    ,    double, Nout);
    332   int    *dX            ; ALLOCATE (dX    ,    int,    Nout);
    333   int    *dY            ; ALLOCATE (dY    ,    int,    Nout);
     332  float  *dX            ; ALLOCATE (dX    ,    float,  Nout);
     333  float  *dY            ; ALLOCATE (dY    ,    float,  Nout);
    334334  double *Rmin          ; ALLOCATE (Rmin  ,    double, Nout);
    335335  double *Rmax          ; ALLOCATE (Rmax  ,    double, Nout);
    … …  
    337337  double *Dmax          ; ALLOCATE (Dmax  ,    double, Nout);
    338338  double *dPix          ; ALLOCATE (dPix  ,    double, Nout);
    339   int    *NX_SUB        ; ALLOCATE (NX_SUB,    int,    Nout);
    340   int    *NY_SUB        ; ALLOCATE (NY_SUB,    int,    Nout);
     339  float  *NX_SUB        ; ALLOCATE (NX_SUB,    float,  Nout);
     340  float  *NY_SUB        ; ALLOCATE (NY_SUB,    float,  Nout);
    341341  int    *TYPE          ; ALLOCATE (TYPE,      int,    Nout);
    342342
    … …  
    559559    int xi = x / tess[myTess].tree->dX[zone][band];
    560560    int yi = y / tess[myTess].tree->dY[zone][band];
     561
     562    xi = MAX(MIN(xi, tess[myTess].tree->NX_SUB - 1.0), 0);
     563    yi = MAX(MIN(yi, tess[myTess].tree->NY_SUB - 1.0), 0);
     564
    561565    int N = xi + tess[myTess].tree->NX_SUB * yi;
    562566 

  • trunk/Ohana/src/libdvo/src/dbExtractAverages.c

    r37035 r37116  
    200200      }
    201201
     202      /* in some cases below we need a measurement which matches the requested photcode.
     203         here is some sample code to find this.
     204
     205          for (i = 0; i < average[0].Nmeasure; i++) {
     206            if (field->photcode->code != measure[i].photcode) continue;
     207            value.Flt = measure[i].M;
     208            break;
     209          }
     210      */
     211
    202212      switch (field->magOption) {
    203213        case MAG_OPTION_MAG:
    … …  
    207217              break;
    208218            case MAG_LEVEL_REF:
    209               // why is measure needed here?
     219              // which measure is needed here?
    210220              // PhotRef  (field->photcode, average, secfilt, measure, field->magClass, field->magSource);
    211221              break;
    … …  
    242252              break;
    243253            case MAG_LEVEL_REF:
    244               // value.Flt = PhotFluxRef  (equiv, average, secfilt, measure, field->magClass, field->magSource);
     254              // XXX which measure is needed here?
     255              // value.Flt = PhotFluxRef  (field->photcode, average, secfilt, measure, field->magClass, field->magSource);
    245256              break;
    246257            case MAG_LEVEL_INST:
    … …  
    258269            case MAG_LEVEL_AVE:
    259270            case MAG_LEVEL_REF:
    260               // value.Flt = PhotFluxAveErr (equiv, average, secfilt, field->magClass, field->magSource); 
     271              value.Flt = PhotFluxAveErr (field->photcode, average, secfilt, field->magClass, field->magSource); 
    261272              break;
    262273            case MAG_LEVEL_INST:
    … …  
    271282
    272283        case MAG_OPTION_STDEV:
     284          value.Flt = PhotMstdev (field->photcode, average, secfilt, field->magClass, field->magSource);
     285          break;
    273286        case MAG_OPTION_CHISQ:
     287          value.Flt = PhotXm (field->photcode, average, secfilt);
     288          break;
    274289        case MAG_OPTION_MIN:
     290          value.Flt = PhotMmin (field->photcode, average, secfilt);
     291          break;
    275292        case MAG_OPTION_MAX:
     293          value.Flt = PhotMmax (field->photcode, average, secfilt);
     294          break;
    276295        case MAG_OPTION_NCODE: {
    277296          int Nsec = GetPhotcodeNsec (field->photcode->code);
    … …  
    281300        }
    282301        case MAG_OPTION_NPHOT: {
    283           int Nsec = GetPhotcodeNsec (field->photcode->code);
    284           if (Nsec == -1) break;
    285           value.Int = secfilt[Nsec].Nused;
     302          value.Int = PhotNphot (field->photcode, average, secfilt, field->magClass, field->magSource);
    286303          break;
    287304        }
    288305        case MAG_OPTION_UC_DIST:
    289         case MAG_OPTION_STACK_DET_ID:
    290           value.Int = PhotStackID (field->photcode, average, secfilt);
    291           break;
     306          value.Flt = PhotUCdist (field->photcode, average, secfilt);
     307          break;
     308        // XX case MAG_OPTION_STACK_DET_ID:
     309        // XX   value.Int = PhotStackID (field->photcode, average, secfilt);
     310        // XX   break;
    292311        case MAG_OPTION_FLAGS: {
    293312          int Nsec = GetPhotcodeNsec (field->photcode->code);
    … …  
    301320      break;
    302321
    303 # if (0)
    304       switch (field->magMode) {
    305         case MAG_AVE:
    306           value.Flt = PhotAve  (field->photcode, average, secfilt);
    307           break;
    308         case MAG_REF:
    309           value.Flt = PhotRef  (field->photcode, average, secfilt, measure);
    310           break;
    311         case MAG_INST:
    312         case MAG_CAT:
    313         case MAG_SYS:
    314         case MAG_REL:
    315         case MAG_CAL:
    316           // XXX need to code this correctly: this returns just the first matching value
    317           value.Flt = NAN;
    318           for (i = 0; i < average[0].Nmeasure; i++) {
    319             if (field->photcode->code != measure[i].photcode) continue;
    320             value.Flt = measure[i].M;
    321             break;
    322           }
    323           break;
    324         case MAG_CHISQ:
    325           value.Flt = PhotXm (field->photcode, average, secfilt);
    326           break;
    327         case MAG_ERR:
    328           if ((field->photcode->type == PHOT_REF) || (field->photcode->type == PHOT_DEP)) {
    329             for (i = 0; i < average[0].Nmeasure; i++) {
    330               if (field->photcode->code != measure[i].photcode) continue;
    331               value.Flt = measure[i].dM;
    332               break;
    333             }
    334           } else {
    335             value.Flt = PhotAveErr (field->photcode, average, secfilt);
    336           }
    337           break;
    338         case MAG_STDEV:
    339           value.Flt = PhotMstdev (field->photcode, average, secfilt);
    340           break;
    341         case MAG_PHOT_FLAGS:
    342           if ((field->photcode->type == PHOT_REF) || (field->photcode->type == PHOT_DEP)) {
    343             for (i = 0; i < average[0].Nmeasure; i++) {
    344               if (field->photcode->code != measure[i].photcode) continue;
    345               value.Int = measure[i].photFlags;
    346               break;
    347             }
    348           } else {
    349             int Nsec = GetPhotcodeNsec (field->photcode->code);
    350             if (Nsec == -1) break;
    351             value.Int = secfilt[Nsec].flags;
    352           }
    353           break;
    354         case MAG_NCODE:
    355           // XXX push these into dvo_photcode_ops APIs
    356           // XXX do I need to allow for conversion to equiv?
    357           int Nsec = GetPhotcodeNsec (field->photcode->code);
    358           if (Nsec == -1) break;
    359           value.Int = secfilt[Nsec].Ncode;
    360           break;
    361         case MAG_NPHOT:
    362           int Nsec = GetPhotcodeNsec (field->photcode->code);
    363           if (Nsec == -1) break;
    364           value.Int = secfilt[Nsec].Nused;
    365           break;
    366 
    367         case MAG_APER:
    368         case MAG_APER_AVE:
    369           value.Flt = PhotAperAve (field->photcode, average, secfilt);
    370           break;
    371 
    372         case MAG_APER_REF:
    373           value.Flt = PhotAperRef (field->photcode, average, secfilt, measure);
    374           break;
    375 
    376         case MAG_KRON:
    377         case MAG_KRON_AVE:
    378           value.Flt = PhotKronAve (field->photcode, average, secfilt);
    379           break;
    380 
    381         case MAG_KRON_REF:
    382           value.Flt = PhotKronRef (field->photcode, average, secfilt, measure);
    383           break;
    384 
    385         case MAG_KRON_ERR:
    386           value.Flt = PhotKronAveErr (field->photcode, average, secfilt);
    387           break;
    388 
    389         case MAG_MIN:
    390           value.Flt = PhotMmin (field->photcode, average, secfilt);
    391           break;
    392         case MAG_MAX:
    393           value.Flt = PhotMmax (field->photcode, average, secfilt);
    394           break;
    395         case MAG_UC_DIST:
    396           value.Flt = PhotUCdist (field->photcode, average, secfilt);
    397           break;
    398         case MAG_STACK_DET_ID:
    399           value.Int = PhotStackID (field->photcode, average, secfilt);
    400           break;
    401 
    402         case MAG_FLUX_PSF:
    403           value.Flt = PhotAveFluxPSF (field->photcode, average, secfilt);
    404           break;
    405         case MAG_FLUX_PSF_ERR:
    406           value.Flt = PhotAvedFluxPSF (field->photcode, average, secfilt);
    407           break;
    408         case MAG_FLUX_KRON:
    409           value.Flt = PhotAveFluxKron (field->photcode, average, secfilt);
    410           break;
    411         case MAG_FLUX_KRON_ERR:
    412           value.Flt = PhotAvedFluxKron (field->photcode, average, secfilt);
    413           break;
    414       }
    415       break;
    416     case AVE_dMAG:
    417       value.Flt = PhotAveErr (field->photcode, average, secfilt);
    418       break;
    419     case AVE_Xm:
    420       value.Flt = PhotXm (field->photcode, average, secfilt);
    421       break;
    422 # endif
    423322    case AVE_TYPE:
    424323      break;

  • trunk/Ohana/src/libdvo/src/dbExtractMeasures.c

    r37035 r37116  
    248248        case MAG_OPTION_NPHOT:
    249249        case MAG_OPTION_UC_DIST:
    250         case MAG_OPTION_STACK_DET_ID:
    251250        case MAG_OPTION_FLAGS:
    252251        case MAG_OPTION_NONE:

  • trunk/Ohana/src/libdvo/src/dbFields.c

    r37035 r37116  
    7878  if (!strcasecmp (string, "nphot"))        return (MAG_OPTION_NPHOT);
    7979  if (!strcasecmp (string, "uc_dist"))      return (MAG_OPTION_UC_DIST);
    80   if (!strcasecmp (string, "stack_det_id")) return (MAG_OPTION_STACK_DET_ID);
    8180  if (!strcasecmp (string, "flags"))        return (MAG_OPTION_FLAGS);
    8281
    … …  
    145144  if (!strcasecmp (string, "max"))         return (MAG_MAX);
    146145  if (!strcasecmp (string, "ucdist"))      return (MAG_UC_DIST);
    147   if (!strcasecmp (string, "stackDetectID")) return (MAG_STACK_DET_ID);
    148146  if (!strcasecmp (string, "fluxpsf"))     return (MAG_FLUX_PSF);
    149147  if (!strcasecmp (string, "fluxpsferr"))  return (MAG_FLUX_PSF_ERR);
    … …  
    296294
    297295  switch (field->magOption) {
    298     case MAG_OPTION_STACK_DET_ID:
    299296    case MAG_OPTION_NCODE:
    300297    case MAG_OPTION_NPHOT:

  • trunk/Ohana/src/libdvo/src/dvo_catalog.c

    r37035 r37116  
    197197  secfilt->dFapStk     = NAN;
    198198
    199   secfilt->stackDetectID = 0;
     199  secfilt->stackPrmryOff = -1;
     200  secfilt->stackBestOff  = -1;
    200201
    201202  secfilt->MpsfWrp     = NAN;
    … …  
    860861    catalog[0].lensobj = NULL;
    861862  }
    862   // if (catalog[0].X_t)      free (catalog[0].X_t);
    863   // if (catalog[0].Y_t)      free (catalog[0].Y_t);
    864   if (catalog[0].nOwn_t)      free (catalog[0].nOwn_t);
    865   if (catalog[0].found_t)     free (catalog[0].found_t);
    866   // if (catalog[0].image_t)  free (catalog[0].image_t);
    867   // if (catalog[0].mosaic_t) free (catalog[0].mosaic_t);
    868   if (catalog[0].foundWarp_t) free (catalog[0].foundWarp_t);
     863  if (catalog[0].nOwn_t)      { free (catalog[0].nOwn_t); catalog[0].nOwn_t = NULL; }
     864  if (catalog[0].found_t)     { free (catalog[0].found_t); catalog[0].found_t = NULL; }
     865  if (catalog[0].foundWarp_t) { free (catalog[0].foundWarp_t); catalog[0].foundWarp_t = NULL; }
    869866}
    870867

  • trunk/Ohana/src/libdvo/src/dvo_convert.c

    r37035 r37116  
    392392    return (lensing); }
    393393
     394# define SKIPPING_FORMAT(NAME, FORMAT, TYPE)            \
     395  if (!strcmp (extname, NAME)) {                                        \
     396    fprintf (stderr, "ERROR: format %s not defined for lensing, skipping\n", NAME); \
     397    *Nlensing = 0; \
     398    return NULL; \
     399  }
     400
    394401  if (!strcmp (extname, "DVO_LENSING")) {
    395402    lensing = gfits_table_get_Lensing (ftable, Nlensing, NULL);
    … …  
    402409  }
    403410
    404 // CONVERT_FORMAT ("DVO_LENSING_PS1_REF",         PS1_REF,         PS1_REF);
    405 // CONVERT_FORMAT ("DVO_LENSING_ELIXIR",          ELIXIR,          Elixir);
    406 // CONVERT_FORMAT ("DVO_LENSING_LONEOS",          LONEOS,          Loneos);
    407 // CONVERT_FORMAT ("DVO_LENSING_PANSTARRS_DEV_0", PANSTARRS_DEV_0, Panstarrs_DEV_0);
    408 // CONVERT_FORMAT ("DVO_LENSING_PANSTARRS_DEV_1", PANSTARRS_DEV_1, Panstarrs_DEV_1);
    409 // CONVERT_FORMAT ("DVO_LENSING_PS1_DEV_1",       PS1_DEV_1,       PS1_DEV_1);
    410 // CONVERT_FORMAT ("DVO_LENSING_PS1_DEV_2",       PS1_DEV_2,       PS1_DEV_2);
    411 // CONVERT_FORMAT ("DVO_LENSING_PS1_V1",          PS1_V1,          PS1_V1);
    412 // CONVERT_FORMAT ("DVO_LENSING_PS1_V2",          PS1_V2,          PS1_V2);
    413 // CONVERT_FORMAT ("DVO_LENSING_PS1_V3",          PS1_V3,          PS1_V3);
    414 // CONVERT_FORMAT ("DVO_LENSING_PS1_V4",          PS1_V4,          PS1_V4);
    415   CONVERT_FORMAT ("DVO_LENSING_PS1_V5",          PS1_V5,          PS1_V5);
     411  SKIPPING_FORMAT ("DVO_LENSING_PS1_REF",         PS1_REF,         PS1_REF);
     412  SKIPPING_FORMAT ("DVO_LENSING_ELIXIR",          ELIXIR,          Elixir);
     413  SKIPPING_FORMAT ("DVO_LENSING_LONEOS",          LONEOS,          Loneos);
     414  SKIPPING_FORMAT ("DVO_LENSING_PANSTARRS_DEV_0", PANSTARRS_DEV_0, Panstarrs_DEV_0);
     415  SKIPPING_FORMAT ("DVO_LENSING_PANSTARRS_DEV_1", PANSTARRS_DEV_1, Panstarrs_DEV_1);
     416  SKIPPING_FORMAT ("DVO_LENSING_PS1_DEV_1",       PS1_DEV_1,       PS1_DEV_1);
     417  SKIPPING_FORMAT ("DVO_LENSING_PS1_DEV_2",       PS1_DEV_2,       PS1_DEV_2);
     418  SKIPPING_FORMAT ("DVO_LENSING_PS1_V1",          PS1_V1,          PS1_V1);
     419  SKIPPING_FORMAT ("DVO_LENSING_PS1_V2",          PS1_V2,          PS1_V2);
     420  SKIPPING_FORMAT ("DVO_LENSING_PS1_V3",          PS1_V3,          PS1_V3);
     421  SKIPPING_FORMAT ("DVO_LENSING_PS1_V4",          PS1_V4,          PS1_V4);
     422  CONVERT_FORMAT  ("DVO_LENSING_PS1_V5",          PS1_V5,          PS1_V5);
    416423# undef CONVERT_FORMAT
    417424

  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c

    r37035 r37116  
    251251    out[i].Mstdev        = in[i].Mstdev;     
    252252    out[i].ubercalDist   = in[i].ubercalDist;     
    253     out[i].stackDetectID = in[i].stackDetectID;     
     253    out[i].stackPrmryOff = in[i].stackPrmryOff;     
     254    out[i].stackBestOff  = in[i].stackBestOff;     
    254255 }
    255256  return (out);
    … …  
    281282    out[i].Mstdev        = in[i].Mstdev;     
    282283    out[i].ubercalDist   = in[i].ubercalDist;     
    283     out[i].stackDetectID = in[i].stackDetectID;
     284    out[i].stackPrmryOff = in[i].stackPrmryOff;     
     285    out[i].stackBestOff  = in[i].stackBestOff;     
    284286  }
    285287  return (out);

  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c

    r37035 r37116  
    319319    out[i].dFapStk       = in[i].dFapStk;
    320320
    321     out[i].stackDetectID = in[i].stackDetectID;     
     321    out[i].stackPrmryOff = in[i].stackPrmryOff;     
     322    out[i].stackBestOff  = in[i].stackBestOff;     
    322323
    323324    out[i].MpsfWrp       = in[i].MpsfWrp;
    … …  
    387388    out[i].dFapStk       = in[i].dFapStk;
    388389
    389     out[i].stackDetectID = in[i].stackDetectID;     
     390    out[i].stackPrmryOff = in[i].stackPrmryOff;     
     391    out[i].stackBestOff  = in[i].stackBestOff;     
    390392
    391393    out[i].MpsfWrp       = in[i].MpsfWrp;

  • trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c

    r37035 r37116  
    687687}
    688688
    689 float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt) {
    690 
    691   int Ns;
    692   float Mstdev;
    693 
    694   if (code == NULL) return NAN;
    695 
    696   Ns = photcodes[0].hashNsec[code[0].code];
    697   Mstdev  = (Ns == -1) ? NAN : 0.001*secfilt[Ns].Mstdev;
     689float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source) {
     690
     691  if (code == NULL) return NAN;
     692
     693  int Ns = photcodes[0].hashNsec[code[0].code];
     694  if (Ns == -1) return (NAN);
     695
     696  float Mstdev = NAN;
     697  switch (source) {
     698    case MAG_SRC_CHP:
     699      switch (class) {
     700        case MAG_CLASS_PSF:
     701          Mstdev = secfilt[Ns].Mstdev;
     702          break;
     703        case MAG_CLASS_KRON:
     704          Mstdev = secfilt[Ns].sMkron;
     705          break;
     706        case MAG_CLASS_APER:
     707          Mstdev = secfilt[Ns].sMap;
     708          break;
     709        default:
     710          break;
     711      }
     712      break;
     713    case MAG_SRC_WRP:
     714      switch (class) {
     715        case MAG_CLASS_PSF:
     716          Mstdev = secfilt[Ns].sFpsfWrp;
     717          break;
     718        case MAG_CLASS_KRON:
     719          Mstdev = secfilt[Ns].sFkronWrp;
     720          break;
     721        case MAG_CLASS_APER:
     722          Mstdev = secfilt[Ns].sFapWrp;
     723          break;
     724        default:
     725          break;
     726      }
     727      break;
     728    case MAG_SRC_STK:
     729    default:
     730      break;
     731  }
    698732  return (Mstdev);
     733}
     734
     735int PhotNphot (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source) {
     736
     737  if (code == NULL) return 0;
     738
     739  int Ns = photcodes[0].hashNsec[code[0].code];
     740  if (Ns == -1) return 0;
     741
     742  int Nphot = 0;
     743  switch (source) {
     744    case MAG_SRC_CHP:
     745      switch (class) {
     746        case MAG_CLASS_PSF:
     747          Nphot = secfilt[Ns].Nused;
     748          break;
     749        case MAG_CLASS_KRON:
     750          Nphot = secfilt[Ns].NusedKron;
     751          break;
     752        case MAG_CLASS_APER:
     753          Nphot = secfilt[Ns].NusedAp;
     754          break;
     755        default:
     756          break;
     757      }
     758      break;
     759    case MAG_SRC_WRP:
     760      switch (class) {
     761        case MAG_CLASS_PSF:
     762          Nphot = secfilt[Ns].NusedWrp;
     763          break;
     764        case MAG_CLASS_KRON:
     765          Nphot = secfilt[Ns].NusedKronWrp;
     766          break;
     767        case MAG_CLASS_APER:
     768          Nphot = secfilt[Ns].NusedApWrp;
     769          break;
     770        default:
     771          break;
     772      }
     773      break;
     774    case MAG_SRC_STK:
     775    default:
     776      break;
     777  }
     778  return (Nphot);
    699779}
    700780
    … …  
    783863}
    784864
    785 unsigned int PhotStackID (PhotCode *code, Average *average, SecFilt *secfilt) {
    786 
    787   int Ns;
    788   unsigned int ID;
    789 
    790   if (code == NULL) return 0;
    791 
    792   Ns = photcodes[0].hashNsec[code[0].code];
    793   ID = (Ns == -1) ? 0 : secfilt[Ns].stackDetectID;
    794   return (ID);
    795 }
     865// XX unsigned int PhotStackID (PhotCode *code, Average *average, SecFilt *secfilt) {
     866// XX
     867// XX   int Ns;
     868// XX   unsigned int ID;
     869// XX
     870// XX   if (code == NULL) return 0;
     871// XX
     872// XX   Ns = photcodes[0].hashNsec[code[0].code];
     873// XX   ID = (Ns == -1) ? 0 : secfilt[Ns].stackDetectID;
     874// XX   return (ID);
     875// XX }
    796876
    797877// Xm is now (2014.07.03) stored as the chisq except in dvo formats which use as short
    … …  
    11121192  }
    11131193  return (Fave);
     1194}
     1195
     1196float PhotFluxAveErr (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source) {
     1197
     1198  if (code == NULL) return NAN;
     1199
     1200  int Ns = photcodes[0].hashNsec[code[0].code];
     1201  if (Ns == -1) return (NAN);
     1202
     1203  float dFave = NAN;
     1204  switch (source) {
     1205    case MAG_SRC_CHP:
     1206      switch (class) {
     1207        case MAG_CLASS_PSF:
     1208          dFave = secfilt[Ns].dM * MagToFlux(secfilt[Ns].M);
     1209          break;
     1210        case MAG_CLASS_KRON:
     1211          dFave = secfilt[Ns].dMkron * MagToFlux(secfilt[Ns].Mkron);
     1212          break;
     1213        case MAG_CLASS_APER:
     1214          dFave = secfilt[Ns].dMap * MagToFlux(secfilt[Ns].Map);
     1215          break;
     1216        default:
     1217          break;
     1218      }
     1219      break;
     1220    case MAG_SRC_WRP:
     1221      switch (class) {
     1222        case MAG_CLASS_PSF:
     1223          dFave = secfilt[Ns].dFpsfWrp;
     1224          break;
     1225        case MAG_CLASS_KRON:
     1226          dFave = secfilt[Ns].dFkronWrp;
     1227          break;
     1228        case MAG_CLASS_APER:
     1229          dFave = secfilt[Ns].dFapWrp;
     1230          break;
     1231        default:
     1232          break;
     1233      }
     1234      break;
     1235    case MAG_SRC_STK:
     1236      switch (class) {
     1237        case MAG_CLASS_PSF:
     1238          dFave = secfilt[Ns].dFpsfStk;
     1239          break;
     1240        case MAG_CLASS_KRON:
     1241          dFave = secfilt[Ns].dFkronStk;
     1242          break;
     1243        case MAG_CLASS_APER:
     1244          dFave = secfilt[Ns].dFapStk;
     1245          break;
     1246        default:
     1247          break;
     1248      }
     1249      break;
     1250    default:
     1251      break;
     1252  }
     1253  return (dFave);
    11141254}
    11151255

  • trunk/Ohana/src/opihi/dvo/gstar.c

    r37049 r37116  
    2323  GSTAR_STACK_FLUX_KRON,
    2424  GSTAR_STACK_FLUX_KRON_ERR,
     25  GSTAR_STACK_FLUX_APER,
     26  GSTAR_STACK_FLUX_APER_ERR,
     27  GSTAR_WARP_FLUX_PSF,
     28  GSTAR_WARP_FLUX_PSF_ERR,
     29  GSTAR_WARP_FLUX_KRON,
     30  GSTAR_WARP_FLUX_KRON_ERR,
     31  GSTAR_WARP_FLUX_APER,
     32  GSTAR_WARP_FLUX_APER_ERR,
     33  GSTAR_STACK_MAG_PSF,
     34  GSTAR_STACK_MAG_PSF_ERR,
     35  GSTAR_STACK_MAG_KRON,
     36  GSTAR_STACK_MAG_KRON_ERR,
     37  GSTAR_STACK_MAG_APER,
     38  GSTAR_STACK_MAG_APER_ERR,
     39  GSTAR_WARP_MAG_PSF,
     40  GSTAR_WARP_MAG_PSF_ERR,
     41  GSTAR_WARP_MAG_KRON,
     42  GSTAR_WARP_MAG_KRON_ERR,
     43  GSTAR_WARP_MAG_APER,
     44  GSTAR_WARP_MAG_APER_ERR,
    2545} GSTAR_SECF_CODES;
    2646
    … …  
    3151void sort_index_float (float *X, int *IDX, int N) {
    3252
    33 # define SWAPFUNC(A,B){ float tmp; int itmp;    \
    34   tmp = X[A]; X[A] = X[B]; X[B] = tmp; \
    35   itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \
    36 }
     53# define SWAPFUNC(A,B){ float tmp; int itmp;            \
     54    tmp = X[A]; X[A] = X[B]; X[B] = tmp;                \
     55    itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp;      \
     56  }
    3757# define COMPARE(A,B)(X[A] < X[B])
    3858
    … …  
    5272  off_t i, Nstars, *N1;
    5373  off_t j, k, m, N, Nlo, Nhi;
    54   int Nsecfilt, NPTS, QUIET, FULL_OUTPUT, STACK_OUTPUT, INST, SHOW_MASKS;
     74  int Nsecfilt, NPTS, QUIET, INST, SHOW_MASKS;
    5575  int found, GetMeasures;
    5676  int SaveVectors;
    … …  
    109129  }
    110130
    111   FULL_OUTPUT = FALSE;
     131  int FULL_OUTPUT = FALSE;
    112132  if ((N = get_argument (argc, argv, "-full"))) {
    113133    FULL_OUTPUT = TRUE;
    … …  
    115135  }
    116136
    117   STACK_OUTPUT = FALSE;
    118   if ((N = get_argument (argc, argv, "-stack"))) {
     137  int STACK_OUTPUT = FALSE;
     138  if ((N = get_argument (argc, argv, "+stack"))) {
    119139    STACK_OUTPUT = TRUE;
     140    remove_argument (N, &argc, argv);
     141  }
     142
     143  int EXTRA_OUTPUT = FALSE;
     144  if ((N = get_argument (argc, argv, "+extras"))) {
     145    EXTRA_OUTPUT = TRUE;
     146    remove_argument (N, &argc, argv);
     147  }
     148
     149  int WARP_OUTPUT = FALSE;
     150  if ((N = get_argument (argc, argv, "+warp"))) {
     151    WARP_OUTPUT = TRUE;
     152    remove_argument (N, &argc, argv);
     153  }
     154
     155  int FLUX_OUTPUT = FALSE;
     156  if ((N = get_argument (argc, argv, "+flux"))) {
     157    FLUX_OUTPUT = TRUE;
     158    remove_argument (N, &argc, argv);
     159  }
     160
     161  int MAG_OUTPUT = TRUE;
     162  if ((N = get_argument (argc, argv, "-mag"))) {
     163    MAG_OUTPUT = FALSE;
     164    remove_argument (N, &argc, argv);
     165  }
     166
     167  int KRON_OUTPUT = FALSE;
     168  if ((N = get_argument (argc, argv, "+kron"))) {
     169    KRON_OUTPUT = TRUE;
     170    remove_argument (N, &argc, argv);
     171  }
     172
     173  int APER_OUTPUT = FALSE;
     174  if ((N = get_argument (argc, argv, "+aper"))) {
     175    APER_OUTPUT = TRUE;
    120176    remove_argument (N, &argc, argv);
    121177  }
    … …  
    185241  // an error exit status here is a significant error
    186242  if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) {
    187       fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename);
    188       exit (2);
     243    fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename);
     244    exit (2);
    189245  }
    190246  dvo_catalog_unlock (&catalog);
    … …  
    248304      k = N1[i];
    249305      if (!QUIET) {
    250         gprint (GP_LOG, "star: "OFF_T_FMT"\n",  i);
     306        gprint (GP_LOG, "star: "OFF_T_FMT"\n",  k);
    251307        gprint (GP_LOG, "%11.7f ", catalog.average[k].R);
    252308        gprint (GP_LOG, "%11.7f ", catalog.average[k].D);
    … …  
    259315       
    260316        if (FULL_OUTPUT) {
    261             gprint (GP_LOG, "%f ",     catalog.average[k].dR);
    262             gprint (GP_LOG, "%f ",     catalog.average[k].dD);
    263             gprint (GP_LOG, "%f ",     catalog.average[k].uR);
    264             gprint (GP_LOG, "%f ",     catalog.average[k].uD);
    265             gprint (GP_LOG, "%f ",     catalog.average[k].duR);
    266             gprint (GP_LOG, "%f ",     catalog.average[k].duD);
    267             gprint (GP_LOG, "%f ",     catalog.average[k].P);
    268             gprint (GP_LOG, "%f   ",     catalog.average[k].dP);
    269 
    270             gprint (GP_LOG, "%f   ",     catalog.average[k].ChiSqPM);
    271             gprint (GP_LOG, "%f   ",     catalog.average[k].ChiSqPar);
    272 
    273             date = ohana_sec_to_date (catalog.average[k].Tmean);
    274             gprint (GP_LOG, "%20s ",     date);
    275             gprint (GP_LOG, "%f   ",     catalog.average[k].Trange / 86400.0);
     317          gprint (GP_LOG, "%f ",     catalog.average[k].dR);
     318          gprint (GP_LOG, "%f ",     catalog.average[k].dD);
     319          gprint (GP_LOG, "%f ",     catalog.average[k].uR);
     320          gprint (GP_LOG, "%f ",     catalog.average[k].uD);
     321          gprint (GP_LOG, "%f ",     catalog.average[k].duR);
     322          gprint (GP_LOG, "%f ",     catalog.average[k].duD);
     323          gprint (GP_LOG, "%f ",     catalog.average[k].P);
     324          gprint (GP_LOG, "%f   ",     catalog.average[k].dP);
     325
     326          gprint (GP_LOG, "%f   ",     catalog.average[k].ChiSqPM);
     327          gprint (GP_LOG, "%f   ",     catalog.average[k].ChiSqPar);
     328
     329          date = ohana_sec_to_date (catalog.average[k].Tmean);
     330          gprint (GP_LOG, "%20s ",     date);
     331          gprint (GP_LOG, "%f   ",     catalog.average[k].Trange / 86400.0);
    276332        }
    277333
    278         gprint (GP_LOG, "\n");
     334        gprint (GP_LOG, "\n\n");
    279335     
    280336        /* filter names */
     337        gprint (GP_LOG, "filter     : ");
    281338        for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_FILTER_NAMES);
    282339        gprint (GP_LOG, "\n");
    283340
    284341        /* average mags */
     342        gprint (GP_LOG, "chp_psf_ave: ");
    285343        for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG);
    286344        gprint (GP_LOG, "\n");
    287345
    288346        /* average mag errors */
     347        gprint (GP_LOG, "chp_psf_err: ");
    289348        for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_ERR);
    290349        gprint (GP_LOG, "\n");
    291350
    292351        /* average mag chisq */
     352        gprint (GP_LOG, "chp_psf_chi: ");
    293353        for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_CHISQ);
    294354        gprint (GP_LOG, "\n");
    295355
    296         if (FULL_OUTPUT) {
     356       
     357        if (APER_OUTPUT) {
    297358          /* Map */
     359          gprint (GP_LOG, "chp_ap__ave: ");
    298360          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_AP_MAG);
    299361          gprint (GP_LOG, "\n");
    300 
     362        }
     363
     364        if (KRON_OUTPUT) {
    301365          /* Mkron */
     366          gprint (GP_LOG, "chp_krn_ave: ");
    302367          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_KRON_MAG);
    303368          gprint (GP_LOG, "\n");
    304369
    305370          /* dMkron */
     371          gprint (GP_LOG, "chp_krn_err: ");
    306372          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_KRON_MAG_ERR);
    307373          gprint (GP_LOG, "\n");
    308 
     374        }
     375
     376        if (EXTRA_OUTPUT) {
    309377          /* Mmin */
     378          gprint (GP_LOG, "chp_psf_min: ");
    310379          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_MIN);
    311380          gprint (GP_LOG, "\n");
    312381
    313382          /* Mmax */
     383          gprint (GP_LOG, "chp_psf_max: ");
    314384          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_MAX);;
    315385          gprint (GP_LOG, "\n");
    316386
    317387          /* UCDIST */
     388          gprint (GP_LOG, "ubercal_dst: ");
    318389          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_UCDIST);;
    319390          gprint (GP_LOG, "\n");
    320391
    321392          /* secfilt flags */
     393          gprint (GP_LOG, "filtflags:   ");
    322394          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_SECF_FLAGS);
    323395          gprint (GP_LOG, "\n");
    324 
    325396        }
    326         if (STACK_OUTPUT) {
     397        if (STACK_OUTPUT && MAG_OUTPUT) {
    327398          /* FluxPSF */
     399          gprint (GP_LOG, "stk_psf_ave: ");
     400          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_PSF);
     401          gprint (GP_LOG, "\n");
     402 
     403          /* dFluxPSF */
     404          gprint (GP_LOG, "stk_psf_err: ");
     405          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_PSF_ERR);
     406          gprint (GP_LOG, "\n");
     407 
     408          if (KRON_OUTPUT) {
     409            /* MagKron */
     410            gprint (GP_LOG, "stk_krn_ave: ");
     411            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_KRON);
     412            gprint (GP_LOG, "\n");
     413 
     414            /* dMagKron */
     415            gprint (GP_LOG, "stk_krn_err: ");
     416            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_KRON_ERR);
     417            gprint (GP_LOG, "\n");
     418          }
     419
     420          if (APER_OUTPUT) {
     421            /* MagAper */
     422            gprint (GP_LOG, "stk_ap__ave: ");
     423            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_APER);
     424            gprint (GP_LOG, "\n");
     425 
     426            /* dMagAper */
     427            gprint (GP_LOG, "stk_ap__err: ");
     428            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_APER_ERR);
     429            gprint (GP_LOG, "\n");
     430          }
     431        }
     432        if (WARP_OUTPUT && MAG_OUTPUT) {
     433          /* FluxPSF */
     434          gprint (GP_LOG, "wrp_psf_ave: ");
     435          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_PSF);
     436          gprint (GP_LOG, "\n");
     437 
     438          /* dFluxPSF */
     439          gprint (GP_LOG, "wrp_psf_err: ");
     440          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_PSF_ERR);
     441          gprint (GP_LOG, "\n");
     442 
     443          if (KRON_OUTPUT) {
     444            /* MagKron */
     445            gprint (GP_LOG, "wrp_krn_ave: ");
     446            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_KRON);
     447            gprint (GP_LOG, "\n");
     448 
     449            /* dMagKron */
     450            gprint (GP_LOG, "wrp_krn_err: ");
     451            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_KRON_ERR);
     452            gprint (GP_LOG, "\n");
     453
     454          }
     455
     456          if (APER_OUTPUT) {
     457            /* MagAper */
     458            gprint (GP_LOG, "wrp_ap__ave: ");
     459            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_APER);
     460            gprint (GP_LOG, "\n");
     461 
     462            /* dMagAper */
     463            gprint (GP_LOG, "wrp_ap__err: ");
     464            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_APER_ERR);
     465            gprint (GP_LOG, "\n");
     466          }
     467        }
     468        if (STACK_OUTPUT && FLUX_OUTPUT) {
     469          /* FluxPSF */
     470          gprint (GP_LOG, "stk_psf_ave: ");
    328471          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_PSF);
    329472          gprint (GP_LOG, "\n");
    330 
     473 
    331474          /* dFluxPSF */
     475          gprint (GP_LOG, "stk_psf_err: ");
    332476          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_PSF_ERR);
    333477          gprint (GP_LOG, "\n");
    334 
    335           /* FluxKron */
    336           for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON);
    337           gprint (GP_LOG, "\n");
    338 
    339           /* dFluxKron */
    340           for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON_ERR);
    341           gprint (GP_LOG, "\n");
    342 
     478 
     479          if (KRON_OUTPUT) {
     480            /* FluxKron */
     481            gprint (GP_LOG, "stk_krn_ave: ");
     482            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON);
     483            gprint (GP_LOG, "\n");
     484 
     485            /* dFluxKron */
     486            gprint (GP_LOG, "stk_krn_err: ");
     487            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON_ERR);
     488            gprint (GP_LOG, "\n");
     489          }
     490
     491          if (APER_OUTPUT) {
     492            /* FluxAper */
     493            gprint (GP_LOG, "stk_ap_ave: ");
     494            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_APER);
     495            gprint (GP_LOG, "\n");
     496 
     497            /* dFluxAper */
     498            gprint (GP_LOG, "stk_ap_err: ");
     499            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_APER_ERR);
     500            gprint (GP_LOG, "\n");
     501          }
    343502        }
    344       }
     503        if (WARP_OUTPUT && FLUX_OUTPUT) {
     504          /* FluxPSF */
     505          gprint (GP_LOG, "wrp_psf_ave: ");
     506          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_PSF);
     507          gprint (GP_LOG, "\n");
     508 
     509          /* dFluxPSF */
     510          gprint (GP_LOG, "wrp_psf_ave: ");
     511          for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_PSF_ERR);
     512          gprint (GP_LOG, "\n");
     513 
     514          if (KRON_OUTPUT) {
     515            /* FluxKron */
     516            gprint (GP_LOG, "wrp_krn_ave: ");
     517            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_KRON);
     518            gprint (GP_LOG, "\n");
     519 
     520            /* dFluxKron */
     521            gprint (GP_LOG, "wrp_krn_err: ");
     522            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_KRON_ERR);
     523            gprint (GP_LOG, "\n");
     524          }
     525
     526          if (APER_OUTPUT) {
     527            /* FluxAper */
     528            gprint (GP_LOG, "wrp_ap__ave: ");
     529            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_APER);
     530            gprint (GP_LOG, "\n");
     531 
     532            /* dFluxAper */
     533            gprint (GP_LOG, "wrp_ap__err: ");
     534            for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_APER_ERR);
     535            gprint (GP_LOG, "\n");
     536          }
     537        }
     538      }
     539      gprint (GP_LOG, "\n");
     540     
    345541      if (GetMeasures || SaveVectors) {
    346542
    … …  
    413609
    414610            if (FULL_OUTPUT) {
    415                 gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mcal);
    416                 gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Map);
    417                 gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mkron);
    418                 gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].dMkron);
    419                 gprint (GP_LOG, "%5.1f ", pow(10.0, 0.4*catalog.measure[Nv].dt));
    420                 gprint (GP_LOG, "%5.3f ", catalog.measure[Nv].airmass);
    421                 gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].az);
    422                 gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Xccd);
    423                 gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Yccd);
    424                 gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dXccd));
    425                 gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dYccd));
    426                 gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Sky);
    427                 gprint (GP_LOG, "%5.1f ", catalog.measure[Nv].dSky);
    428                 gprint (GP_LOG, "%8d ", catalog.measure[Nv].averef);
    429                 gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].detID);
    430                 gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].imageID);
    431                 gprint (GP_LOG, "%.3f ", catalog.measure[Nv].psfQF);
    432                 gprint (GP_LOG, "%.1f ", catalog.measure[Nv].psfChisq);
    433                 gprint (GP_LOG, "%.1f ", catalog.measure[Nv].crNsigma);
    434                 gprint (GP_LOG, "%.1f ", catalog.measure[Nv].extNsigma);
    435                 gprint (GP_LOG, "%.1f ", FromShortDegrees(catalog.measure[Nv].theta));
     611              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mcal);
     612              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Map);
     613              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mkron);
     614              gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].dMkron);
     615              gprint (GP_LOG, "%5.1f ", pow(10.0, 0.4*catalog.measure[Nv].dt));
     616              gprint (GP_LOG, "%5.3f ", catalog.measure[Nv].airmass);
     617              gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].az);
     618              gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Xccd);
     619              gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Yccd);
     620              gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dXccd));
     621              gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dYccd));
     622              gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Sky);
     623              gprint (GP_LOG, "%5.1f ", catalog.measure[Nv].dSky);
     624              gprint (GP_LOG, "%8d ", catalog.measure[Nv].averef);
     625              gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].detID);
     626              gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].imageID);
     627              gprint (GP_LOG, "%5.3f ", catalog.measure[Nv].psfQF);
     628              gprint (GP_LOG, "%7.1f ", catalog.measure[Nv].psfChisq);
     629              gprint (GP_LOG, "%3.1f ", catalog.measure[Nv].crNsigma);
     630              gprint (GP_LOG, "%4.1f ", catalog.measure[Nv].extNsigma);
     631              gprint (GP_LOG, "%5.1f ", FromShortDegrees(catalog.measure[Nv].theta));
    436632            }
    437             if (STACK_OUTPUT) {
    438                 gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxPSF);
    439                 gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxPSF);
    440                 gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxKron);
    441                 gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxKron);
     633            if (FLUX_OUTPUT) {
     634              gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxPSF);
     635              gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxPSF);
     636              gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxKron);
     637              gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxKron);
    442638            }
    443639            gprint (GP_LOG, "\n");
    … …  
    547743
    548744  switch (type) {
    549     case GSTAR_AVE_MAG: /* average mags */
    550       if (seq == -1) {
    551         print_double (NAN);
    552       } else {
    553         print_double (secfilt[seq].M);
    554       }
    555       break;
    556 
    557     case GSTAR_AVE_MAG_ERR: /* average mags errors */
    558       if (seq == -1) {
    559         print_double (NAN);
    560       } else {
    561         print_double (secfilt[seq].dM);
     745    case GSTAR_AVE_MAG_MIN: /* average ap mags */
     746      if (seq == -1) {
     747        print_double (NAN);
     748      } else {
     749        print_double (secfilt[seq].Mmin);
     750      }
     751      break;
     752
     753    case GSTAR_AVE_MAG_MAX: /* average ap mags */
     754      if (seq == -1) {
     755        print_double (NAN);
     756      } else {
     757        print_double (secfilt[seq].Mmax);
    562758      }
    563759      break;
    … …  
    592788      break;
    593789
     790      /*** CHIP : MAG ***/
     791
     792    case GSTAR_AVE_MAG: /* average mags */
     793      if (seq == -1) {
     794        print_double (NAN);
     795      } else {
     796        print_double (secfilt[seq].M);
     797      }
     798      break;
     799
     800    case GSTAR_AVE_MAG_ERR: /* average mags errors */
     801      if (seq == -1) {
     802        print_double (NAN);
     803      } else {
     804        print_double (secfilt[seq].dM);
     805      }
     806      break;
     807
    594808    case GSTAR_AVE_AP_MAG: /* average ap mags */
    595809      if (seq == -1) {
    … …  
    600814      break;
    601815
    602     case GSTAR_AVE_MAG_MIN: /* average ap mags */
    603       if (seq == -1) {
    604         print_double (NAN);
    605       } else {
    606         print_double (secfilt[seq].Mmin);
    607       }
    608       break;
    609 
    610     case GSTAR_AVE_MAG_MAX: /* average ap mags */
    611       if (seq == -1) {
    612         print_double (NAN);
    613       } else {
    614         print_double (secfilt[seq].Mmax);
    615       }
    616       break;
    617 
    618816    case GSTAR_AVE_KRON_MAG: /* average ap mags */
    619817      if (seq == -1) {
    … …  
    632830      break;
    633831
     832      /******************* STACK : FLUX **********************/
     833
    634834    case GSTAR_STACK_FLUX_PSF: /* average ap mags */
    635835      if (seq == -1) {
    … …  
    664864      break;
    665865
     866    case GSTAR_STACK_FLUX_APER: /* average ap mags */
     867      if (seq == -1) {
     868        print_double (NAN);
     869      } else {
     870        print_double_exp (secfilt[seq].FapStk);
     871      }
     872      break;
     873
     874    case GSTAR_STACK_FLUX_APER_ERR: /* average ap mags */
     875      if (seq == -1) {
     876        print_double (NAN);
     877      } else {
     878        print_double_exp (secfilt[seq].dFapStk);
     879      }
     880      break;
     881
     882      /******************* STACK : MAG **********************/
     883
     884    case GSTAR_STACK_MAG_PSF: /* average ap mags */
     885      if (seq == -1) {
     886        print_double (NAN);
     887      } else {
     888        print_double (secfilt[seq].MpsfStk);
     889      }
     890      break;
     891
     892    case GSTAR_STACK_MAG_PSF_ERR: /* average ap mags */
     893      if (seq == -1) {
     894        print_double (NAN);
     895      } else {
     896        print_double (secfilt[seq].dFpsfStk / secfilt[seq].FpsfStk);
     897      }
     898      break;
     899
     900    case GSTAR_STACK_MAG_KRON: /* average ap mags */
     901      if (seq == -1) {
     902        print_double (NAN);
     903      } else {
     904        print_double (secfilt[seq].MkronStk);
     905      }
     906      break;
     907
     908    case GSTAR_STACK_MAG_KRON_ERR: /* average ap mags */
     909      if (seq == -1) {
     910        print_double (NAN);
     911      } else {
     912        print_double (secfilt[seq].dFkronStk / secfilt[seq].FkronStk);
     913      }
     914      break;
     915
     916    case GSTAR_STACK_MAG_APER: /* average ap mags */
     917      if (seq == -1) {
     918        print_double (NAN);
     919      } else {
     920        print_double (secfilt[seq].MapStk);
     921      }
     922      break;
     923
     924    case GSTAR_STACK_MAG_APER_ERR: /* average ap mags */
     925      if (seq == -1) {
     926        print_double (NAN);
     927      } else {
     928        print_double (secfilt[seq].dFapStk / secfilt[seq].FapStk);
     929      }
     930      break;
     931
     932      /******************* WARP : FLUX **********************/
     933
     934    case GSTAR_WARP_FLUX_PSF: /* average ap mags */
     935      if (seq == -1) {
     936        print_double (NAN);
     937      } else {
     938        print_double_exp (secfilt[seq].FpsfWrp);
     939      }
     940      break;
     941
     942    case GSTAR_WARP_FLUX_PSF_ERR: /* average ap mags */
     943      if (seq == -1) {
     944        print_double (NAN);
     945      } else {
     946        print_double_exp (secfilt[seq].dFpsfWrp);
     947      }
     948      break;
     949
     950    case GSTAR_WARP_FLUX_KRON: /* average ap mags */
     951      if (seq == -1) {
     952        print_double (NAN);
     953      } else {
     954        print_double_exp (secfilt[seq].FkronWrp);
     955      }
     956      break;
     957
     958    case GSTAR_WARP_FLUX_KRON_ERR: /* average ap mags */
     959      if (seq == -1) {
     960        print_double (NAN);
     961      } else {
     962        print_double_exp (secfilt[seq].dFkronWrp);
     963      }
     964      break;
     965
     966    case GSTAR_WARP_FLUX_APER: /* average ap mags */
     967      if (seq == -1) {
     968        print_double (NAN);
     969      } else {
     970        print_double_exp (secfilt[seq].FapWrp);
     971      }
     972      break;
     973
     974    case GSTAR_WARP_FLUX_APER_ERR: /* average ap mags */
     975      if (seq == -1) {
     976        print_double (NAN);
     977      } else {
     978        print_double_exp (secfilt[seq].dFapWrp);
     979      }
     980      break;
     981
     982      /******************* WARP : MAG **********************/
     983
     984    case GSTAR_WARP_MAG_PSF: /* average ap mags */
     985      if (seq == -1) {
     986        print_double (NAN);
     987      } else {
     988        print_double (secfilt[seq].MpsfWrp);
     989      }
     990      break;
     991
     992    case GSTAR_WARP_MAG_PSF_ERR: /* average ap mags */
     993      if (seq == -1) {
     994        print_double (NAN);
     995      } else {
     996        print_double (secfilt[seq].dFpsfWrp / secfilt[seq].FpsfWrp);
     997      }
     998      break;
     999
     1000    case GSTAR_WARP_MAG_KRON: /* average ap mags */
     1001      if (seq == -1) {
     1002        print_double (NAN);
     1003      } else {
     1004        print_double (secfilt[seq].MkronWrp);
     1005      }
     1006      break;
     1007
     1008    case GSTAR_WARP_MAG_KRON_ERR: /* average ap mags */
     1009      if (seq == -1) {
     1010        print_double (NAN);
     1011      } else {
     1012        print_double (secfilt[seq].dFkronWrp / secfilt[seq].FkronWrp);
     1013      }
     1014      break;
     1015
     1016    case GSTAR_WARP_MAG_APER: /* average ap mags */
     1017      if (seq == -1) {
     1018        print_double (NAN);
     1019      } else {
     1020        print_double (secfilt[seq].MapWrp);
     1021      }
     1022      break;
     1023
     1024    case GSTAR_WARP_MAG_APER_ERR: /* average ap mags */
     1025      if (seq == -1) {
     1026        print_double (NAN);
     1027      } else {
     1028        print_double (secfilt[seq].dFapWrp / secfilt[seq].FapWrp);
     1029      }
     1030      break;
     1031
    6661032  }
    6671033}

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

    r37037 r37116  
    440440int MatchImageName (off_t meas, int cat, char *name);
    441441int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID);
     442int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID);
    442443
    443444int load_tess (char *treefile);
    … …  
    451452int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
    452453int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found);
    453 int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell);
     454int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure);
    454455int setGlobalObjStats (Average *average, Measure *measure);
    455456

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r37037 r37116  
    491491  if (projectID[i] != myProjectionID) return FALSE;
    492492  if (skycellID[i] != mySkycellID) return FALSE;
     493 
     494  return TRUE;
     495}
     496
     497// returns image.Mcal - ff(x,y)
     498int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID) {
     499
     500  off_t i;
     501
     502  if (!MeasureToImage) return FALSE;
     503
     504  if (meas < 0) return FALSE;
     505
     506  i = MeasureToImage[cat][meas];
     507  if (i == -1) return FALSE;
     508
     509  if (tessID[i]    == -1) return FALSE;
     510  if (projectID[i] == -1) return FALSE;
     511  if (skycellID[i] == -1) return FALSE;
     512
     513  *myTessID = tessID[i];
     514  *myProjectionID = projectID[i];
     515  *mySkycellID =  skycellID[i];
    493516 
    494517  return TRUE;

  • trunk/Ohana/src/relphot/src/StarOps.c

    r37037 r37116  
    10861086          if (ecode != thisCode) { continue; }
    10871087
     1088          // SKIP gpc1 stack data
     1089          if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
     1090         
     1091          // SKIP gpc1 forced-warp data
     1092          if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
     1093         
    10881094          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    10891095          Mcal  = getMcal  (m, i, flatcorr, catalog);
    … …  
    11311137          if (ecode != thisCode) { continue; }
    11321138
     1139          // SKIP gpc1 stack data
     1140          if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
     1141         
     1142          // SKIP gpc1 forced-warp data
     1143          if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
     1144         
    11331145          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
    11341146          Mcal  = getMcal  (m, i, flatcorr, catalog);

  • trunk/Ohana/src/relphot/src/bcatalog.c

    r37037 r37116  
    7878        int Nsec = GetPhotcodeNsec(thisCode);
    7979
    80         fprintf (stderr, "really use dvo_secfilt_init?");
     80        // fprintf (stderr, "really use dvo_secfilt_init?");
    8181        dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]);
    8282

  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r37037 r37116  
    5050  int isSetMrelFinal = (pass >= 0);
    5151
    52   char *primaryCell = NULL;
    53   if (isSetMrelFinal) {
    54     ALLOCATE (primaryCell, char, DVO_MAX_PATH);
    55   }
     52  // XX char *primaryCell = NULL;
     53  // XX if (isSetMrelFinal) {
     54  // XX   ALLOCATE (primaryCell, char, DVO_MAX_PATH);
     55  // XX }
    5656
    5757  for (j = 0; j < catalog[Nc].Naverage; j++) {
    … …  
    6464    // only apply Stack operation on setMrelFinal in first pass
    6565    if (isSetMrelFinal && (pass == 0)) {
    66       setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], primaryCell);
     66      setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m]);
    6767      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
    6868    }
    … …  
    7474    }
    7575  }
    76   if (primaryCell) free (primaryCell);
     76  // if (primaryCell) free (primaryCell);
    7777  return (TRUE);
    7878}
    … …  
    441441
    442442// only apply Stack operation on setMrelFinal in first pass
    443 int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell) {
    444 
    445   off_t k, ID;
    446 
    447   float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
     443// this function has 3 goals, not to be confused:
     444// 1) find and mark the PRIMARY detections (regardless of the quality of the detection)
     445// 2) select the BEST detections per filter (regardless of PRIMARY)
     446// 3) apply the zero point and AB->Jy transformations
     447int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure) {
     448
     449  off_t k;
     450
     451  float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0;
    448452
    449453  // set the primary projection cell and skycell for this coordinate
    … …  
    467471    int haveStack = FALSE;
    468472
    469     // need to find the measurement closest to the center of its skycell, as well as the
    470     // closest for the subset of primary projection cells
    471 
    472     float stackCenterOffsetMin = 1e9;
    473     off_t stackCenterMeasureMin = -1;
    474 
    475     float stackPrimaryOffsetMin = 1e9;
    476     off_t stackPrimaryMeasureMin = -1;
     473    float psfQFbest = 0.0;
     474
     475    off_t stackBestMeasure = -1;
     476    off_t stackPrimaryMeasure = -1;
    477477
    478478    int isBad = FALSE;
    … …  
    481481    off_t meas = measureOffset;
    482482    for (k = 0; k < average[0].Nmeasure; k++, meas++) {
     483
     484      // only examine gpc1 stack data
     485      if (!isGPC1stack(measure[k].photcode)) continue;
    483486
    484487      // skip measurements that do not match the current photcode
    … …  
    489492      // clear this bit for all measurements
    490493      measure[k].dbFlags &= ~ID_MEAS_STACK_PRIMARY;
    491 
    492       // only examine gpc1 stack data
    493       if (!isGPC1stack(measure[k].photcode)) continue;
     494      measure[k].dbFlags &= ~ID_MEAS_STACK_PHOT_SRC;
    494495
    495496      haveStack = TRUE;
    496497      haveStackObject = TRUE;
    497498     
     499      // match measurement to its image (this is just a check, right?)
     500      if (getImageEntry (meas, cat) < 0) {
     501        fprintf (stderr, "*** warning : this should be impossible\n");
     502        abort();
     503        continue;
     504      }
     505
     506      // find the PRIMARY measurement
     507
     508      // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the
     509      // primary skycell for this position
     510      if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
     511        stackPrimaryMeasure = k;
     512        measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
     513        secfilt[Nsec].stackPrmryOff = meas;
     514      }
     515
     516      // now choose the BEST measurements (may also be PRIMARY)
     517
    498518      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad);
    499519
    500       if (getImageEntry (meas, cat) < 0) continue;
    501520      // measurements without an image are either external reference photometry or
    502521      // data for which the associated image has not been loaded (probably because of
    … …  
    510529      if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid);
    511530
    512       Msys = PhotSys (&measure[k], &average[0], &secfilt[0], MAG_CLASS_PSF);
    513       if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys);
    514 
    515       unsigned int stackImageID;
    516 
    517       // which stack image should we use for the mean value?
    518       // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell
    519       if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
    520         float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
    521         if (stackPrimaryOffset < stackPrimaryOffsetMin) {
    522           stackPrimaryOffsetMin = stackPrimaryOffset;
    523           // stackPrimaryIDmin = stackImageID;
    524           stackPrimaryMeasureMin = k;
    525         }
    526       }
    527 
    528       // get the center distance for the generic case:
    529       float stackCenterOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
    530       if (stackCenterOffset < stackCenterOffsetMin) {
    531         stackCenterOffsetMin = stackCenterOffset;
    532         // stackCenterIDmin = stackImageID;
    533         stackCenterMeasureMin = k;
     531      // NOTE: negative and insignificant vlues values are allowed, but not NAN flux values
     532      Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF);
     533      if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst);
     534
     535      if (measure[k].psfQFperf > psfQFbest) {
     536        psfQFbest = measure[k].psfQFperf;
     537        stackBestMeasure = k;
    534538      }
    535539
    … …  
    554558
    555559    // measurements which are bad will not have a valid stack entry and are skipped
    556     k = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
     560    k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure;
    557561    if (k < 0) continue;
    558 
    559562
    560563    // we are now populating stackDetectID not stack Image ID in secfilt
    561564    // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
    562     ID = measure[k].extID; // for the stack, this is the stackDetectID
     565    // ID = measure[k].extID; // for the stack, this is the stackDetectID
    563566
    564567    // get the zero point for the selected image
    … …  
    605608    secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
    606609
     610    // Jy to AB mags
    607611    secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN;
    608612    secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
    609613    secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN;
    610614
    611     secfilt[Nsec].stackDetectID = ID;
    612 
     615    secfilt[Nsec].stackBestOff = k + measureOffset;
     616
     617    // this is the measurement used by secfilt[]
    613618    measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
    614     if (stackPrimaryMeasureMin >= 0) {
    615       measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
     619    if (k == stackPrimaryMeasure) {
    616620      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
    617621    }
    … …  
    642646
    643647  off_t k;
    644   float Fsys = 0, Mcal= 0;
     648  float Mcal= 0;
    645649
    646650  // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
    … …  
    664668
    665669  // option for a test print
    666   if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
     670  if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) {
    667671    fprintf (stderr, "test obj\n");
    668672    print_measure_set_alt (average, secfilt, measure);
    … …  
    680684    if (found[Nsec]) continue; 
    681685     
     686    // I need to restrict the measurements coming only from the same skycell as the secfilt values
     687    // if we cannot make that association, skip this set of warps
     688    int stkTessID, stkProjID, stkSkycellID;
     689    if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue;
     690
    682691    off_t meas = measureOffset;
    683692
    … …  
    694703      if (!code) continue;
    695704      if (code->equiv != thisCode) { continue; }
     705
     706      measure[k].dbFlags &= ~ID_MEAS_WARP_USED;
     707
     708      if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue;
    696709
    697710      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad);
    … …  
    720733      // are small (and not measured)
    721734      float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
    722       if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys);
     735      if (isnan(Fpsf)) SKIP_THIS_MEAS(Nsys);
    723736      // if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
    724737      // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
    725738
    726739      float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
     740      if (isnan(dFpsf)) SKIP_THIS_MEAS (Nsys);
     741
     742      int isBad = (measure[k].psfQF < 0.85);
     743      isBad |= isnan(measure[k].psfQF);
     744      if (isBad) SKIP_THIS_MEAS (Nsys);
     745
     746      int isPoor = (measure[k].psfQFperf < 0.85);
     747      if ((pass == 0) && isPoor) SKIP_THIS_MEAS (Nsys);
     748
    727749      dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error
    728750      Fpsflist[Npsf] = Fpsf * Fcal;
    … …  
    730752      wpsflist[Npsf] = 1.0;
    731753      Npsf ++;
     754
     755      measure[k].dbFlags |= ID_MEAS_WARP_USED;
    732756
    733757      float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
Note: See TracChangeset for help on using the changeset viewer.