IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36528 for trunk


Ignore:
Timestamp:
Feb 18, 2014, 2:42:32 PM (12 years ago)
Author:
eugene
Message:

add support for PS1_V4 alt format; fix some dvomerge issues; split lists for psf, aper, kron mag averages

Location:
trunk/Ohana
Files:
21 edited

Legend:

Unmodified
Added
Removed
  • trunk/Ohana

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/dvomerge

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c

    r34277 r36528  
    167167  off_t i, oldID, newID;
    168168
     169  if (!IDmap->old) {
     170    fprintf (stderr, "input database has image IDs, but no Image table\n");
     171    return FALSE;
     172  }
     173
     174  off_t lastID = IDmap->old[IDmap->Nmap-1];
     175
    169176  for (i = 0; i < catalog[0].Nmeasure; i++) {
    170177    oldID = catalog[0].measure[i].imageID;
    171178    if (oldID == 0) continue;
    172179
    173     if (!IDmap->old) {
    174       fprintf (stderr, "input database has image IDs, but no Image table\n");
    175     }
    176 
    177180    newID = dvo_map_image_ID (IDmap, oldID);
    178181    if (newID == 0) {
     182      if (oldID > lastID) {
     183        fprintf (stderr, "problem with image IDs : input out of range\n");
     184        fprintf (stderr, "old ID: "OFF_T_FMT", last ID: "OFF_T_FMT"\n", oldID, lastID);
     185        exit (2);
     186      }
    179187      if (!IDmap->notFound[oldID]) {
    180188        fprintf (stderr, "cannot find image ID "OFF_T_FMT"\n",  oldID);
  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c

    r35765 r36528  
    111111  SetPhotcodeTable(NULL);
    112112
    113   dvomergeUpdate_catalogs (input, output, outsky, inlist, NsecfiltInput, NsecfiltOutput, &IDmap, secfiltMap);
     113  int status = dvomergeUpdate_catalogs (input, output, outsky, inlist, NsecfiltInput, NsecfiltOutput, &IDmap, secfiltMap);
    114114
    115115  // save the output sky table copy
     
    123123  gettimeofday (&stop, NULL);
    124124  dtime = DTIME (stop, start);
     125
     126  if (!status) {
     127    fprintf (stderr, "ERROR: elapsed time %9.4f sec\n", dtime);
     128    exit (3);
     129  }
     130
    125131  fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime);
    126 
    127132  exit (0);
    128133}
  • trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c

    r36486 r36528  
    140140    }
    141141
     142    dvo_update_image_IDs (IDmap, &incatalog);
     143
    142144    // merge input into the appropriate output tables
    143145    for (j = 0; j < outlist[0].Nregions; j++) {
     
    156158      LoadCatalog (&outcatalog, outlist[0].regions[j], outcatalog.filename, "w", NsecfiltOutput);
    157159
    158       dvo_update_image_IDs (IDmap, &incatalog);
    159160      merge_catalogs_old (outlist[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap);
    160161
    161162      outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    162163
    163       if (outstat[i].missed) {
     164      if (outstat[j].missed) {
    164165        dmhObjectAdd (outstat[j].history, &outcatalog.header, inStats);
    165166      }
     
    280281  }
    281282  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
    282     HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     283    int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     284    if (!status) {
     285      fprintf (stderr, "error running one of the remote clients\n");
     286      return status;
     287    }
    283288  }
    284289
  • trunk/Ohana/src/libdvo/include/dvo.h

    r36490 r36528  
    412412} MeasureTiny;
    413413
     414/** STRUCT DEFINITION **/
     415typedef struct {
     416  double           R;                    // RA (decimal degrees )
     417  double           D;                    // DEC (decimal degrees )
     418  float            dR;                   // RA error (arcsec)
     419  float            dD;                   // DEC error (arcsec)
     420  float            uR;                   // RA*cos(D) proper-motion (arcsec/year)
     421  float            uD;                   // DEC proper-motion (arcsec/year)
     422  float            duR;                  // RA*cos(D) p-m error (arcsec/year)
     423  float            duD;                  // DEC p-m error (arcsec/year)
     424  float            P;                    // parallax (arcsec)
     425  float            dP;                   // parallax error (arcsec)
     426  float            ChiSqAve;             // astrometry analysis chisq
     427  float            ChiSqPM;              // astrometry analysis chisq
     428  float            ChiSqPar;             // astrometry analysis chisq
     429  int              Tmean;                // mean epoch (PM,PAR ref) (unix time seconds)
     430  int              Trange;               // mean epoch (PM,PAR ref) (unix time seconds)
     431  float            Xp;                   // unused
     432  unsigned short   Npos;                 // number of detections used for astrometry
     433  unsigned short   Nmeasure;             // number of psf measurements
     434  unsigned short   Nmissing;             // number of missings
     435  unsigned short   Nextend;              // number of extended measurements
     436  uint32_t         measureOffset;        // offset to first psf measurement
     437  uint32_t         missingOffset;        // offset to first missing obs
     438  uint32_t         extendOffset;         // offset to first extended measurement
     439  uint32_t         flags;                // average object flags (star; ghost; etc)
     440  uint32_t         photFlagsUpper;       // upper bit of 2 bit summary of per-measure photflags
     441  uint32_t         photFlagsLower;       // lower bit of 2 bit summary of per-measure photflags
     442  unsigned int     objID;                // unique ID for object in table
     443  unsigned int     catID;                // unique ID for table in which object was first realized
     444  uint64_t         extID;                // external ID for object (eg PSPS objID)
     445} Average_PS1_V4alt;
     446
     447Average_PS1_V4alt *gfits_table_get_Average_PS1_V4alt (FTable *table, off_t *Ndata, char *swapped);
     448Average *Average_PS1_V4alt_ToInternal (Average_PS1_V4alt *in, off_t Nvalues);
    414449
    415450/** STRUCT DEFINITION **/
     
    592627
    593628float PhotAperInst (Measure *measure);
     629float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt);
    594630float PhotKronInst (Measure *measure);
    595631float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt);
  • trunk/Ohana/src/libdvo/src/dbExtractAverages.c

    r36084 r36528  
    243243          break;
    244244
     245        case MAG_APER:
     246          value.Flt = PhotAperAve (field->photcode, average, secfilt);
     247          break;
     248
    245249        case MAG_KRON:
    246250          value.Flt = PhotKronAve (field->photcode, average, secfilt);
  • trunk/Ohana/src/libdvo/src/dvo_convert.c

    r36490 r36528  
    6363    fprintf (stderr, "EXTNAME missing for average table\n");
    6464    return (FALSE);
     65  }
     66
     67  // block to convert broken tables (PS1_V4 made before the Xfix addition)
     68  if (!strcmp (extname, "DVO_AVERAGE_PS1_V4") && (ftable[0].header[0].Naxis[0] == 120)) {
     69    Average_PS1_V4alt *tmpAverage;
     70    tmpAverage = gfits_table_get_Average_PS1_V4alt (ftable, Naverage, NULL);
     71    if (!tmpAverage) {
     72      fprintf (stderr, "ERROR: failed to read averages\n");
     73      exit (2);
     74    }
     75    average = Average_PS1_V4alt_ToInternal (tmpAverage, *Naverage);
     76    free (tmpAverage);
     77    *format = DVO_FORMAT_PS1_V4;
     78    return (average);
    6579  }
    6680
  • trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c

    r36490 r36528  
    647647}
    648648
     649// 'primary' is needed to conform with the API for Loneos and Elixir, but is not used
     650Average *Average_PS1_V4alt_ToInternal (Average_PS1_V4alt *in, off_t Nvalues) {
     651
     652  off_t i;
     653  Average *out;
     654
     655  ALLOCATE_ZERO (out, Average, Nvalues);
     656
     657  for (i = 0; i < Nvalues; i++) {
     658    dvo_average_init (&out[i]);
     659
     660    out[i].R             = in[i].R;     
     661    out[i].D             = in[i].D;     
     662    out[i].dR            = in[i].dR;
     663    out[i].dD            = in[i].dD;
     664    out[i].uR            = in[i].uR;
     665    out[i].uD            = in[i].uD;
     666    out[i].duR           = in[i].duR;
     667    out[i].duD           = in[i].duD;
     668    out[i].P             = in[i].P;
     669    out[i].dP            = in[i].dP;
     670    out[i].ChiSqAve      = in[i].ChiSqAve;     
     671    out[i].ChiSqPM       = in[i].ChiSqPM;     
     672    out[i].ChiSqPar      = in[i].ChiSqPar;     
     673    out[i].Tmean         = in[i].Tmean;     
     674    out[i].Trange        = in[i].Trange;     
     675    out[i].Npos          = in[i].Npos;     
     676    out[i].Nmeasure      = in[i].Nmeasure;     
     677    out[i].Nmissing      = in[i].Nmissing;     
     678    out[i].Nextend       = in[i].Nextend;     
     679    out[i].measureOffset = in[i].measureOffset;
     680    out[i].missingOffset = in[i].missingOffset;
     681    out[i].extendOffset  = in[i].extendOffset;
     682    out[i].flags         = in[i].flags;   
     683    out[i].photFlagsUpper = in[i].photFlagsUpper;   
     684    out[i].photFlagsLower = in[i].photFlagsLower;   
     685    out[i].objID         = in[i].objID;
     686    out[i].catID         = in[i].catID;
     687    out[i].extID         = in[i].extID;
     688  }
     689  return (out);
     690}
     691
     692/******/
     693
     694int gfits_convert_Average_PS1_V4alt (Average_PS1_V4alt *data, off_t size, off_t nitems) {
     695
     696  off_t i;
     697  unsigned char *byte, tmp;
     698
     699  if (size != 120) {
     700    fprintf (stderr, "WARNING: mismatch in data types Average_PS1_V4alt: "OFF_T_FMT" vs %d\n",  size,  120);
     701    return (FALSE);
     702  }
     703
     704  /* provide initial values to avoid compiler warnings for non-BYTE_SWAP arch */
     705  i = tmp = 0;
     706  byte = NULL;
     707
     708# ifdef BYTE_SWAP
     709  byte = (unsigned char *) data;
     710  for (i = 0; i < nitems; i++, byte += 129) {
     711    /** BYTE SWAP **/
     712    SWAP_DBLE (0); // RA
     713    SWAP_DBLE (8); // DEC
     714    SWAP_WORD (16); // RA_ERR
     715    SWAP_WORD (20); // DEC_ERR
     716    SWAP_WORD (24); // U_RA
     717    SWAP_WORD (28); // U_DEC
     718    SWAP_WORD (32); // V_RA_ERR
     719    SWAP_WORD (36); // V_DEC_ERR
     720    SWAP_WORD (40); // PAR
     721    SWAP_WORD (44); // PAR_ERR
     722    SWAP_WORD (48); // CHISQ_POS
     723    SWAP_WORD (52); // CHISQ_PM
     724    SWAP_WORD (56); // CHISQ_PAP
     725    SWAP_WORD (60); // MEAN_EPOCH
     726    SWAP_WORD (64); // TIME_RANGE
     727    SWAP_WORD (68); // SIGMA
     728    SWAP_BYTE (72); // NUMBER_POS
     729    SWAP_BYTE (74); // NMEASURE
     730    SWAP_BYTE (76); // NMISSING
     731    SWAP_BYTE (78); // NEXTEND
     732    SWAP_WORD (80); // OFF_MEASURE
     733    SWAP_WORD (84); // OFF_MISSING
     734    SWAP_WORD (88); // OFF_EXTEND
     735    SWAP_WORD (92); // FLAGS
     736    SWAP_WORD (96); // PHOTFLAGS_U
     737    SWAP_WORD (100); // PHOTFLAGS_L
     738    SWAP_WORD (104); // OBJ_ID
     739    SWAP_WORD (108); // CAT_ID
     740    SWAP_DBLE (112); // EXT_ID
     741  }
     742# endif 
     743
     744  return (TRUE);
     745}
     746
     747/*** add test of EXTNAME and header-defined columns? ***/
     748/* return internal structure representation */
     749Average_PS1_V4alt *gfits_table_get_Average_PS1_V4alt (FTable *ftable, off_t *Ndata, char *swapped) {
     750
     751  int Ncols;
     752  Average_PS1_V4alt *data;
     753
     754  Ncols = ftable[0].header[0].Naxis[0];
     755  if (Ncols != 120) {
     756    fprintf (stderr, "ERROR: mis-match in table size: width is %d but should be %d bytes\n", Ncols, 176);
     757    return NULL;
     758  }
     759
     760  *Ndata = ftable[0].header[0].Naxis[1];
     761  data = (Average_PS1_V4alt *) ftable[0].buffer;
     762  if ((swapped == NULL) || (*swapped == FALSE)) {
     763    if (!gfits_convert_Average_PS1_V4alt (data, sizeof (Average_PS1_V4alt), *Ndata)) {
     764      return NULL;
     765    }
     766    gfits_table_scale_data (ftable);
     767    if (swapped != NULL) *swapped = TRUE;
     768  }
     769  return (data);
     770}
     771
  • trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c

    r34844 r36528  
    560560}
    561561
     562float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt) {
     563
     564  int Ns;
     565  float Maper;
     566
     567  if (code == NULL) return NAN;
     568
     569  Ns = photcodes[0].hashNsec[code[0].code];
     570  Maper = (Ns == -1) ? NAN : secfilt[Ns].Map;
     571  return (Maper);
     572}
     573
    562574float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt) {
    563575
  • trunk/Ohana/src/libohana/src

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/opihi

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/opihi/cmd.astro

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/opihi/cmd.data

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/opihi/dvo/dvo_client.c

    r35263 r36528  
    9090  }
    9191
     92  // these are set in 'startup.c' for readline-based programs
     93  set_variable ("PID", getpid());
     94  set_str_variable ("KAPA", "kapa");
     95  set_int_variable ("UNSIGN", 0);
     96  gfits_set_unsign_mode (FALSE);
     97
    9298  return;
    9399}
  • trunk/Ohana/src/opihi/lib.shell/stack_math.c

    r36489 r36528  
    11251125  if (!strcmp (op, "floor"))  S_FUNC(floor (M1), ST_SCALAR_FLT);
    11261126  if (!strcmp (op, "ceil"))   S_FUNC(ceil (M1), ST_SCALAR_FLT);
    1127   if (!strcmp (op, "rint"))   S_FUNC(nearbyint (M1), ST_SCALAR_FLT);
     1127  // if (!strcmp (op, "rint"))   S_FUNC(nearbyint (M1), ST_SCALAR_FLT);
    11281128  if (!strcmp (op, "exp"))    S_FUNC(exp (M1), ST_SCALAR_FLT);
    11291129  if (!strcmp (op, "ten"))    S_FUNC(pow (10.0,M1), ST_SCALAR_FLT);
     
    12071207  if (!strcmp (op, "floor"))  V_FUNC(floor (*M1), ST_SCALAR_FLT);
    12081208  if (!strcmp (op, "ceil"))   V_FUNC(ceil (*M1), ST_SCALAR_FLT);
    1209   if (!strcmp (op, "rint"))   V_FUNC(nearbyint (*M1), ST_SCALAR_FLT);
     1209  // if (!strcmp (op, "rint"))   V_FUNC(nearbyint (*M1), ST_SCALAR_FLT);
    12101210  if (!strcmp (op, "exp"))    V_FUNC(exp(*M1), ST_SCALAR_FLT);
    12111211  if (!strcmp (op, "ten"))    V_FUNC(pow(10.0,*M1), ST_SCALAR_FLT);
     
    12821282  if (!strcmp (op, "floor")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = floor (*M1); }}
    12831283  if (!strcmp (op, "ceil"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = ceil (*M1); }}
    1284   if (!strcmp (op, "rint"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = nearbyint (*M1); }}
     1284  // if (!strcmp (op, "rint"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = nearbyint (*M1); }}
    12851285
    12861286  if (!strcmp (op, "exp"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = exp(*M1);          }}
  • trunk/Ohana/src/opihi/pcontrol/CheckBusyJob.c

    r28242 r36528  
    9999
    100100  // XXX runaway job if output too large?
    101   if (job[0].stdout_buf.size > 0x1000000) abort();
    102   if (job[0].stderr_buf.size > 0x1000000) abort();
     101  if (job[0].stdout_buf.size > 0x10000000) abort();
     102  if (job[0].stderr_buf.size > 0x10000000) abort();
    103103
    104104  // job has exited : move to DONE stack
  • trunk/Ohana/src/relastro/src

    • Property svn:mergeinfo deleted
  • trunk/Ohana/src/relphot/include/relphot.h

    r36491 r36528  
    8686  // others are only used on the final output steps.
    8787
    88   double *list;               // list of measure.mag values for a given star
    89   double *dlist;              // mag errors for a star
    90   double *wlist;              // weights to use for mean mags
    91   double *aplist;             // ap mags for a star
    92   double *kronlist;           // kron mags for a star
     88  double *Mpsflist;           // list of measure.mag values for a given star
     89  double *dpsflist;           // mag errors for a star
     90  double *wpsflist;           // weights to use for mean mags
     91
     92  double *Maplist;            // ap mags for a star
     93  double *daplist;            // ap mags for a star
     94  double *waplist;            // ap mags for a star
     95
     96  double *Mkronlist;          // kron mags for a star
    9397  double *dkronlist;          // kron mag errors
     98  double *wkronlist;          // kron mag errors
    9499
    95100  double *psfqf_list;         // psfqf for all filters
     
    97102  double *stargal_list;       // stargal for all filters
    98103
    99   StatType stats;
     104  StatType psfstats;
    100105  StatType apstats;
    101106  StatType kronstats;
  • trunk/Ohana/src/relphot/src/StarOps.c

    r36491 r36528  
    6363  results->Ngrid = 0;
    6464  if (allocLists) {
    65     ALLOCATE (results->list, double, Nmax);
    66     ALLOCATE (results->dlist, double, Nmax);
    67     ALLOCATE (results->wlist, double, Nmax);
     65    ALLOCATE (results->Mpsflist, double, Nmax);
     66    ALLOCATE (results->dpsflist, double, Nmax);
     67    ALLOCATE (results->wpsflist, double, Nmax);
    6868  }
    6969}
    7070
    7171void SetMrelInfoFree (SetMrelInfo *results) {
    72   free (results->list);
    73   free (results->dlist);
    74   free (results->wlist);
     72  free (results->Mpsflist);
     73  free (results->dpsflist);
     74  free (results->wpsflist);
    7575}
    7676
    7777void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
    7878  summary->Nfew  += results->Nfew ;
    79   summary->Ncode  += results->Ncode ;
     79  summary->Ncode += results->Ncode ;
    8080  summary->Nsys  += results->Nsys ;
    8181  summary->Nbad  += results->Nbad ;
     
    158158  SetMrelInfoInit (&summary, FALSE);
    159159  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
    160   ALLOCATE (results.aplist, double, Nmax);
    161   ALLOCATE (results.kronlist, double, Nmax);
     160
     161  ALLOCATE (results.Maplist,  double, Nmax);
     162  ALLOCATE (results.daplist,  double, Nmax);
     163  ALLOCATE (results.waplist,  double, Nmax);
     164
     165  ALLOCATE (results.Mkronlist, double, Nmax);
    162166  ALLOCATE (results.dkronlist, double, Nmax);
     167  ALLOCATE (results.wkronlist, double, Nmax);
    163168
    164169  ALLOCATE (results.psfqf_list, double, Nmax);
     
    183188
    184189  SetMrelInfoFree (&results);
    185   free (results.aplist);
    186   free (results.kronlist);
     190  free (results.Maplist);
     191  free (results.daplist);
     192  free (results.waplist);
     193
     194  free (results.Mkronlist);
    187195  free (results.dkronlist);
     196  free (results.wkronlist);
    188197
    189198  free (results.psfqf_list);
     
    331340  exit (3);
    332341
     342# if (0)
     343
    333344  off_t j, k, m, ID;
    334345  int N;
     
    336347
    337348  StatType stats, apstats, kronstats;
    338   liststats_setmode (&stats, STATMODE);
     349  liststats_setmode (&psfstats, STATMODE);
    339350  liststats_setmode (&apstats, STATMODE);
    340351  liststats_setmode (&kronstats, STATMODE);
     
    837848  if (primaryCell) free (primaryCell);
    838849  return (TRUE);
     850# endif
    839851}
    840852
  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r36491 r36528  
    3434  off_t j;
    3535
    36   liststats_setmode (&results->stats, STATMODE);
     36  liststats_setmode (&results->psfstats, STATMODE);
    3737  liststats_setmode (&results->apstats, STATMODE);
    3838  liststats_setmode (&results->kronstats, STATMODE);
     
    6767
    6868  off_t k;
    69   int N;
    7069  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
    7170
    72   double *list      = results->list;
    73   double *dlist     = results->dlist;
    74   double *wlist     = results->wlist;
    75   double *aplist    = results->aplist;
    76   double *kronlist  = results->kronlist;
     71  // we are measuring means for 3 types of mags: psf, ap, kron.  I am using the psf mag
     72  // error for the ap mags, but krons have their own errors.  it is an open question if I
     73  // should be doing weighted or unweighted fits (this is a user option)
     74  double *Mpsflist  = results->Mpsflist;
     75  double *dpsflist  = results->dpsflist;
     76  double *wpsflist  = results->wpsflist;
     77
     78  double *Maplist   = results->Maplist;
     79  double *daplist   = results->daplist;
     80  double *waplist   = results->waplist;
     81
     82  double *Mkronlist = results->Mkronlist;
    7783  double *dkronlist = results->dkronlist;
    78 
    79   StatType *stats = &results->stats;
    80   StatType *apstats = &results->apstats;
     84  double *wkronlist = results->wkronlist;
     85
     86  StatType *psfstats  = &results->psfstats;
     87  StatType *apstats   = &results->apstats;
    8188  StatType *kronstats = &results->kronstats;
    8289
     
    125132    off_t meas = measureOffset;
    126133
    127     N = 0;
     134    int Nap = 0;
     135    int Npsf = 0;
     136    int Nkron = 0;
    128137    for (k = 0; k < averageT[0].Nmeasure; k++, meas++) {
    129138
     
    164173      int myUbercalDist = getUbercalDist(meas, cat);
    165174      minUbercalDist = MIN(minUbercalDist, myUbercalDist);
     175
     176      int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
     177      float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
     178
     179      if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure
    166180
    167181      if (isSetMrelFinal) {
     
    231245          if (measureT[k].M < MaxMagForceSynth) {
    232246            forceSynth = TRUE;
    233             forceSynthEntry = N;
     247            forceSynthEntry = Npsf;
    234248          } else {
    235249            if (pass < 4) {
     
    240254        }
    241255
    242         // Map (aplist) and Mkron (kronlist,dkronlist) are used to calculate mean mags per filter
     256        // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter
    243257        float Map = PhotAper (&measure[k]);
    244         aplist[N] = Map - Mcal - Mmos - Mgrid;
     258        if (!isnan(Map)) {
     259            Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
     260            daplist[Nap] = dMpsf;
     261            waplist[Nap] = isUbercal ? 10.0 : 1.0;
     262            Nap ++;
     263        }
    245264
    246265        float Mkron = PhotKron (&measure[k]);
    247         kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
    248         dkronlist[N] = measure[k].dMkron;
    249       }
     266        if (!isnan(Mkron)) {
     267          Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
     268          dkronlist[Nkron] = measure[k].dMkron;
     269          wkronlist[Nkron] = isUbercal ? 10.0 : 1.0;
     270          Nkron ++;
     271        }
     272      } // if (isSetMrelFinal)
    250273
    251274      // dlist gives the error per measurement, wlist gives the weight
     
    255278      // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
    256279      // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    257       dlist[N] = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
    258       wlist[N] = 1.0;
    259       list[N]  = Msys - Mcal - Mmos - Mgrid;
    260 
    261       // up-weight the ubercal values (or convergence can take a long time...)
    262       if (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    263         haveUbercal = TRUE;
    264         wlist[N] = 10.0;
    265       }
     280      Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
     281      dpsflist[Npsf] = dMpsf;
     282      wpsflist[Npsf] = isUbercal ? 10.0 : 1.0;
     283
    266284
    267285      // tie down reference photometry if the -refcode (code) option is selected
    268286      // eg, -refcode g_SDSS
    269287      // this probably makes no sense in the context of multifilter analysis
     288      // XXX probably need to use the photocde table to assign reference mag weights.
    270289      if (refPhotcode) {
    271290        if (code->code == refPhotcode->code) {
    272           wlist[N] = 100.0;
    273         }
    274       }
    275       N++;
     291          wpsflist[Npsf] = 100.0;
     292        }
     293      }
     294      Npsf ++;
    276295    }
    277296
     
    288307
    289308    // when performing the grid analysis, STAR_TOOFEW should be set to 1;
    290     if (N < Nminmeas) { /* too few measurements */
    291       // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, N, STAR_TOOFEW);
     309    if (Npsf < Nminmeas) { /* too few measurements */
     310      // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW);
    292311      secfilt[Nsec].flags |= ID_STAR_FEW;
    293312      if (Ncode == 0) {
     
    303322    if (forceSynth) {
    304323      // use the single SYNTH value instead of the other mags here
    305       myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
    306       list[0]  = list[forceSynthEntry];
    307       dlist[0] = dlist[forceSynthEntry];
    308       wlist[0] = wlist[forceSynthEntry];
    309       N = 1;
    310     }
    311     liststats (list, dlist, wlist, N, stats);
    312 
    313     secfilt[Nsec].M  = stats->mean;
    314     secfilt[Nsec].dM = stats->error;
    315     secfilt[Nsec].Xm = (stats->Nmeas > 1) ? 100.0*log10(stats->chisq + 1e-4) : NAN_S_SHORT;
     324      myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error");
     325      Mpsflist[0] = Mpsflist[forceSynthEntry];
     326      dpsflist[0] = dpsflist[forceSynthEntry];
     327      wpsflist[0] = wpsflist[forceSynthEntry];
     328      Npsf = 1;
     329    }
     330    liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
     331
     332    secfilt[Nsec].M  = psfstats->mean;
     333    secfilt[Nsec].dM = psfstats->error;
     334    secfilt[Nsec].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT;
    316335
    317336    // when running -averages, we have no information about the images, so we cannot set this
     
    323342      found[Nsec] = TRUE;
    324343
    325       secfilt[Nsec].Mstdev = 1000.0*stats->sigma; // Mstdev is in millimags (not enough space for more precision)
     344      secfilt[Nsec].Mstdev = 1000.0*psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    326345      // secfilt[Nsec].Ncode = Ncode;
    327       secfilt[Nsec].Nused = stats->Nmeas;
    328 
    329       secfilt[Nsec].M_80 = 1000 * stats->Upper80;
    330       secfilt[Nsec].M_20 = 1000 * stats->Lower20;
     346      secfilt[Nsec].Nused = psfstats->Nmeas;
     347
     348      secfilt[Nsec].M_80 = 1000 * psfstats->Upper80;
     349      secfilt[Nsec].M_20 = 1000 * psfstats->Lower20;
    331350
    332351      // NOTE : use the modified weight for apmags as well as psf mags
    333       liststats (aplist, dlist, wlist, N, apstats);
    334       secfilt[Nsec].Map  = apstats->mean;
    335 
    336       liststats (kronlist, dkronlist, wlist, N, kronstats);
    337       secfilt[Nsec].Mkron  = kronstats->mean;
    338       secfilt[Nsec].dMkron = kronstats->error;
     352      liststats (Maplist, daplist, waplist, Nap, apstats);
     353      secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN;
     354
     355      liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
     356      secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN;
     357      secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN;
    339358
    340359      // NOTE: for 2MASS measurements, Next should be 1, as should N
    341       if ((Next > 0) && (Next > 0.5*N)) {
     360      if ((Next > 0) && (Next > 0.5*Npsf)) {
    342361        secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
    343362      }
Note: See TracChangeset for help on using the changeset viewer.