IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 27, 2015, 4:49:06 PM (11 years ago)
Author:
eugene
Message:

extensive work on relphot, relastro, uniphot, dvomerge aiming to the construction and calibration of PV3

Location:
trunk/Ohana
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r37807 r38986  
    22# define PAR_TOOFEW 5
    33
    4 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset);
    5 int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
    6 
    7 static off_t Nmax;
    8 static double *X, *dX;
    9 static double *Y, *dY;
    10 static double *R, *dR;
    11 static double *D, *dD;
    12 static double *pX;
    13 static double *pY;
    14 static double *T;
    15 static double *dT;
    16 static double *C_blue;
    17 static double *C_red;
    18 
    19 static Coords coords;
    20 
    21 static time_t T2000;
    22 
    23 void initFitStats (FitStats *fitStats) {
    24   fitStats->Nave = 0; 
    25   fitStats->Npm = 0;   
    26   fitStats->Npar = 0; 
    27   fitStats->Nskip = 0;
    28   fitStats->Noffset = 0;
    29   return;
    30 }
    31 
    32 void sumFitStats (FitStats *srcFitStats, FitStats *tgtFitStats) {
    33   tgtFitStats->Nave    += srcFitStats->Nave    ; 
    34   tgtFitStats->Npm     += srcFitStats->Npm     ;   
    35   tgtFitStats->Npar    += srcFitStats->Npar    ; 
    36   tgtFitStats->Nskip   += srcFitStats->Nskip   ;
    37   tgtFitStats->Noffset += srcFitStats->Noffset ;
    38   return;
    39 }
    40 
    41 void initObjectData (Catalog *catalog, int Ncatalog) {
    42 
    43   off_t i, j;
    44  
    45   Nmax = 0;
    46   for (i = 0; i < Ncatalog; i++) {
    47     for (j = 0; j < catalog[i].Naverage; j++) {
    48       Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
    49     }
    50   }
    51 
    52   ALLOCATE (R, double, MAX (1, Nmax));
    53   ALLOCATE (D, double, MAX (1, Nmax));
    54   ALLOCATE (T, double, MAX (1, Nmax));
    55   ALLOCATE (X, double, MAX (1, Nmax));
    56   ALLOCATE (Y, double, MAX (1, Nmax));
    57 
    58   ALLOCATE (dR, double, MAX (1, Nmax));
    59   ALLOCATE (dD, double, MAX (1, Nmax));
    60   ALLOCATE (dT, double, MAX (1, Nmax));
    61   ALLOCATE (dX, double, MAX (1, Nmax));
    62   ALLOCATE (dY, double, MAX (1, Nmax));
    63 
    64   ALLOCATE (pX, double, MAX (1, Nmax));
    65   ALLOCATE (pY, double, MAX (1, Nmax));
    66 
    67   ALLOCATE (C_blue, double, MAX (1, Nmax));
    68   ALLOCATE (C_red,  double, MAX (1, Nmax));
    69 
    70   /* project coordinates to a plane centered on the object with units of arcsec */
    71   InitCoords (&coords, "DEC--SIN");
    72   coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    73 
    74   // use J2000 as a reference time
    75   T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
    76 } 
    77 
    78 void freeObjectData () {
    79 
    80   free (R);
    81   free (D);
    82   free (T);
    83   free (X);
    84   free (Y);
    85 
    86   free (dR);
    87   free (dD);
    88   free (dT);
    89   free (dX);
    90   free (dY);
    91 
    92   free (pX);
    93   free (pY);
    94 
    95   free (C_blue);
    96   free (C_red);
    97 } 
     4int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog);
    985
    996// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    1007// be called with just MeasureTiny set and Measure == NULL
    1018int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop) {
    102 
    103   initObjectData (catalog, Ncatalog);
    1049
    10510  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
    … …  
    10914  }
    11015
    111   FitStats sumStatsChips; initFitStats (&sumStatsChips);
    112   FitStats sumStatsStack; initFitStats (&sumStatsStack);
     16  int NmeasureMax = CatalogMaxNmeasure (catalog, Ncatalog);
     17
     18  // allocate summary stats with Nmax = 0, Nboot = 0
     19  FitStats *sumStatsChips = FitStatsInit (0, 0);
     20  FitStats *sumStatsStack = FitStatsInit (0, 0);
     21
     22  FitStats *fitStatsChips = FitStatsInit (NmeasureMax, N_BOOTSTRAP_SAMPLES);
     23  FitStats *fitStatsStack = FitStatsInit (NmeasureMax, N_BOOTSTRAP_SAMPLES);
     24
     25  AstromErrorSetLoop (Nloop, FALSE);
    11326
    11427  int i;
    … …  
    11730    if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
    11831
    119     FitStats fitStatsChips; initFitStats (&fitStatsChips);
    120     FitStats fitStatsStack; initFitStats (&fitStatsStack);
     32    FitStatsReset (fitStatsChips);
     33    FitStatsReset (fitStatsStack);
    12134
    12235    off_t j;
    … …  
    12942      SecFilt *secfilt = &catalog[i].secfilt[j*Nsecfilt];
    13043
    131       UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack);
    132       UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m, Nloop);
    133     }
    134     if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
    135     if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
    136     sumFitStats (&fitStatsChips, &sumStatsChips);
    137     sumFitStats (&fitStatsStack, &sumStatsStack);
    138   }
    139   freeObjectData ();
    140 
    141   if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
    142   if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
     44      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, fitStatsStack);
     45      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, fitStatsChips, i, m);
     46    }
     47    if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips->Nave,  fitStatsChips->Npm,  fitStatsChips->Npar,  fitStatsChips->Nskip, fitStatsChips->Noffset);
     48    if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack->Nave,  fitStatsStack->Npm,  fitStatsStack->Npar,  fitStatsStack->Nskip, fitStatsStack->Noffset);
     49    FitStatsSum (fitStatsChips, sumStatsChips);
     50    FitStatsSum (fitStatsStack, sumStatsStack);
     51  }
     52
     53  // DumpObjectsWith2MASS (catalog, Ncatalog);
     54
     55  FitStatsFree (fitStatsChips);
     56  FitStatsFree (fitStatsStack);
     57
     58  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips->Nave + sumStatsChips->Npm + sumStatsChips->Npar),  sumStatsChips->Nave,  sumStatsChips->Npm,  sumStatsChips->Npar,  sumStatsChips->Nskip, sumStatsChips->Noffset);
     59  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack->Nave + sumStatsStack->Npm + sumStatsStack->Npar),  sumStatsStack->Nave,  sumStatsStack->Npm,  sumStatsStack->Npar,  sumStatsStack->Nskip, sumStatsStack->Noffset);
     60
     61  FitStatsFree (sumStatsChips);
     62  FitStatsFree (sumStatsStack);
     63
    14364  return (TRUE);
    14465}
    … …  
    14667// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    14768// be called with just MeasureTiny set and Measure == NULL
    148 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset) {
    149 
    150   int setRefColor = areImagesMatched();
     69int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog) {
     70
     71  int i;
     72  for (i = 0; i < Ncatalog; i++) {
     73    off_t j;
     74    for (j = 0; j < catalog[i].Naverage; j++) {
     75      /* calculate the average value of R,D for a single star */
     76      off_t m = catalog[i].average[j].measureOffset;
     77
     78      off_t k;
     79      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
     80        MeasureTiny *measure = &catalog[i].measureT[m+k];
     81        if (measure->dbFlags & ID_MEAS_OBJECT_HAS_2MASS) {
     82          fprintf (stderr, "0x%08x 0x%08x : %12.8f %12.8f %5d\n", catalog[i].average[j].objID, catalog[i].average[j].catID, measure->R, measure->D, measure->photcode);
     83        }
     84      }
     85    }
     86  }
     87  return (TRUE);
     88}
     89
     90// This function operates on both Measure and MeasureTiny.  In the big stages, this should
     91// be called with just MeasureTiny set and Measure == NULL
     92int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff) {
     93
     94  int k;
    15195
    15296  /* calculate the average value of R,D for a single star */
    15397
    154   PMFit fit;    memset (&fit,    0, sizeof(fit));
    155   PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN;
    156   PMFit fitPM;  memset (&fitPM,  0, sizeof(fitPM));  fitPM.chisq = NAN;
    157   PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN;
     98  FitAstromResult fitPos, fitPM, fitPar;
     99  FitAstromResultInit (&fitPos);
     100  FitAstromResultInit (&fitPM);
     101  FitAstromResultInit (&fitPar);
    158102
    159103  // if we fail to fit the astrometry for some reason, we need to set/reset these
    … …  
    167111  if (average[0].Nmeasure == 0) return TRUE;
    168112
    169   int NcBlue = 0;
    170   int NcRed = 0;
    171   int N = 0;
    172 
    173113  int mode = FIT_MODE; // start with the globally-defined fit mode
    174114
    … …  
    177117  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
    178118
    179   // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
    180   off_t k;
    181   for (k = 0; k < average[0].Nmeasure; k++) {
    182 
    183     if (XVERB) {
    184       char *date = ohana_sec_to_date (measure[k].t);
    185       int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
    186       fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
    187       free (date);
    188     }
    189 
    190     // SKIP gpc1 stack data
    191     if (isGPC1stack(measure[k].photcode)) continue;
    192 
    193     // SKIP gpc1 forced-warp data
    194     if (isGPC1warp(measure[k].photcode)) continue;
    195 
    196     // reset the bit to note that a detection was used (or not)
    197     measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    198     if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
    199 
    200     // does the measurement pass the supplied filtering constraints?
    201     // MeasFilterTestTiny does not test psfQF
    202     // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
    203     int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
    204     if (!keepMeasure) {
    205       continue;
    206     }
    207 
    208     double Ri = getMeanR (&measure[k], average, secfilt);
    209     double Di = getMeanD (&measure[k], average, secfilt);
    210 
    211     // if we are correcting for the Galaxy Motion Model, only should apply it here
    212     // (a) when we are working to correct the images (mean R,D assumed to be at J2000) and
    213     // (b) if we think the measure R,D is already at the image epoch position
    214     if (USE_GALAXY_MODEL && applyGalaxyOffset) {
    215       Ri -= measure[k].RoffGAL / 3600.0;
    216       Di -= measure[k].DoffGAL / 3600.0;
    217     }
    218 
    219     // XXX add in dR,dD GAL here
    220 
    221     // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
    222     if (ExcludeBogus) {
    223       coords.crval1 = average[0].R;
    224       coords.crval2 = average[0].D;
    225       double Xi, Yi;
    226       RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
    227       double radius = hypot(Xi, Yi);
    228       if (radius > ExcludeBogusRadius) {
    229         measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
    230         if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
    231         continue;
    232       }
    233       measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    234       if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
    235     }
    236 
    237     // outlier rejection
    238     if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
    239       continue;
    240     }
    241 
    242     R[N] = Ri;
    243     D[N] = Di;
    244 
    245     // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
    246     // T[] is time in years since J2000 (jd = 2451545)
    247     T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
    248 
    249     // dX, dY : error in arcsec --
    250     dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
    251     dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
    252 
    253     // allow a given photcode or measurement to be
    254     // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
    255     if (isnan(dX[N])) continue;
    256     if (isnan(dY[N])) continue;
    257 
    258     // add systematic error in quadrature, if desired
    259     // only do this after the fit has converged (or you will never improve the poor images)
    260     // if (INCLUDE_SYS_ERR) {
    261     // float dRsys = FromShortPixels(measure[k].dRsys);
    262     // dX[N] = hypot(dX[N], dRsys);
    263     // dY[N] = hypot(dY[N], dRsys);
    264     // }
    265 
    266     // dX[N] = 0.1;
    267     // dY[N] = 0.1;
    268 
    269     dT[N] = measure[k].dt;
    270 
    271     // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
    272     // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
    273     // non-circular errors (different values for X and Y), then dR and dD will be
    274     // incorrect: they would need to be rotated to take out the position angle
    275     dR[N] = dX[N] / 3600.0;
    276     dD[N] = dY[N] / 3600.0;
    277 
    278     if (setRefColor) {
    279       float colorBlue = getColorBlue (m+k, i);
    280       if (!isnan(colorBlue)) {
    281         C_blue[NcBlue] = colorBlue;
    282         NcBlue++;
    283       }
    284       float colorRed = getColorRed (m+k, i);
    285       if (!isnan(colorRed)) {
    286         C_red[NcRed] = colorRed;
    287         NcRed++;
    288       }
    289     }
    290 
    291     measure[k].dbFlags |= ID_MEAS_USED_OBJ;
    292     if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
    293 
    294     N++;
    295   } // loop over measurements : average[0].Nmeasure
    296 
    297   if (N < 1) {
     119  // select the measurements to be used in this analysis
     120  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, FALSE);
     121
     122  // if there are no exposure detections, use the stack position
     123  if (fitStats->Npoints < 1) {
    298124    if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) {
    299125      average[0].R  = average[0].Rstk;
    … …  
    301127      average[0].dR = average[0].dRstk;
    302128      average[0].dD = average[0].dDstk;
     129      average[0].flags |= ID_STACK_ASTROM;
    303130    }
    304131    return FALSE;
    305132  }
    306133
    307   // if we have too few good detections for the desired fit, or too limited a
    308   // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
    309   // for even the average position, consider including the lower-quality detections?
    310 
    311   // find Tmin & Tmax from the list of accepted measurements
    312   double Tmean = 0.0;
    313   double Tmin = T[0];
    314   double Tmax = T[0];
    315   for (k = 0; k < N; k++) {
    316     Tmin = MIN(Tmin, T[k]);
    317     Tmax = MAX(Tmax, T[k]);
    318     Tmean += T[k];
    319   }
    320   double Trange = Tmax - Tmin;
    321 
    322   if (RELASTRO_OP == OP_HIGH_SPEED) {
    323     Tmean = 0.5*(Tmax - Tmin);
    324   } else {
    325     Tmean /= (float) N;
    326   }
    327 
    328   /* we need to do the fit in a locally linear space; choose a ref coordinate */
    329   coords.crval1 = R[0];
    330   coords.crval2 = D[0];
     134  double Tmean, Trange, parRange;
     135  FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
    331136
    332137  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
    333138
    334   // project all of the R,D coordinates to a plane centered on this coordinate. set
    335   // the times to be relative to Tmean (this is required for parallax as well)
    336   for (k = 0; k < N; k++) {
    337     RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
    338     T[k] -= Tmean;
    339     if (XVERB) {
    340       fprintf (stderr, OFF_T_FMT" %f %f %f  %f %f +/- %f %f\n",  k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
    341     }
    342   }       
     139  // if we have too few good detections for the desired fit, or too limited a baseline,
     140  // use a fit with fewer parameters.
    343141
    344142  // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
    … …  
    346144    if (Trange < PM_DT_MIN) {
    347145      mode = FIT_AVERAGE;
    348       goto skipPM;
    349     }
    350     if (N <= PM_TOOFEW) {
     146      goto justPosition;
     147    }
     148    if (fitStats->Npoints <= PM_TOOFEW) {
    351149      mode = FIT_AVERAGE;
    352       goto skipPM;
    353     }
    354 
    355     FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
    356 
    357     if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
     150      goto justPosition;
     151    }
     152
     153    if (fitStats->NfitAlloc == 1) {
     154      // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
     155      FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
     156    } else {
     157      fitStats->Nfit = 0;
     158      for (k = 0; k < fitStats->NfitAlloc; k++) {
     159        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     160        if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
     161        fitStats->Nfit ++;
     162      }
     163      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     164      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     165      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     166      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     167    }
     168    FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY);
    358169
    359170    // project Ro, Do back to RA,DEC
    360     XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
    361     if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
     171    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
    362172    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    363     // XXX : does this make sense at 0,360 boundary?
    364173
    365174    fitPM.p  = fitPM.dp  = 0.0;
    … …  
    374183  }
    375184 
    376 skipPM:
    377185  // fit the parallax + proper-motion model
    378186  // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
    379187  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
    380188  if (mode == FIT_PM_AND_PAR) {
    381     if (Trange < PM_DT_MIN) {
    382       mode = FIT_PM_ONLY;
    383       goto skipPAR;
    384     }
    385     if (N <= PAR_TOOFEW) {
    386       mode = FIT_PM_ONLY;
    387       goto skipPAR;
    388     }
    389     float pXmin = +2.0;
    390     float pXmax = -2.0;
    391     float pYmin = +2.0;
    392     float pYmax = -2.0;
    393     for (k = 0; k < N; k++) {
    394       ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
    395       pXmin = MIN (pXmin, pX[k]);
    396       pXmax = MAX (pXmax, pX[k]);
    397       pYmin = MIN (pYmin, pY[k]);
    398       pYmax = MAX (pYmax, pY[k]);
    399     }
    400     float dXRange = pXmax - pXmin;
    401     float dYRange = pYmax - pYmin;
    402     float parRange = hypot (dXRange, dYRange);
    403        
    404189    if (parRange < PAR_FACTOR_MIN) {
    405190      mode = FIT_PM_ONLY;
    406       goto skipPAR;
    407     }
    408 
    409     FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
    410     if (XVERB) fprintf (stderr, "fitted PM+PAR:  %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq);
    411 
    412     XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
     191      goto justPosition;
     192    }
     193    if (fitStats->Npoints <= PAR_TOOFEW) {
     194      mode = FIT_PM_ONLY;
     195      goto justPosition;
     196    }
     197
     198    if (fitStats->NfitAlloc == 1) {
     199      // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
     200      FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
     201    } else {
     202      fitStats->Nfit = 0;
     203      for (k = 0; k < fitStats->NfitAlloc; k++) {
     204        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     205        FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
     206        fitStats->Nfit ++;
     207      }
     208      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     209      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     210      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     211      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     212      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
     213    }
     214    FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR);
     215
     216    // project Ro, Do back to RA,DEC
     217    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
     218    if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     219
    413220    average[0].flags |= ID_STAR_FIT_PAR;
    414221    fitStats->Npar ++;
    415222
    416     if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    417 
    418223    // XXX a hard-wired hack...
    419     if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
     224    if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) {
    420225      mode = FIT_PM_ONLY;
    421226    }
    422227  }       
    423228
    424 skipPAR:
     229justPosition:
    425230  {
    426     // ALWAYS fit the average model
    427     StatType statsR, statsD;
    428     liststats_pos (X, dX, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
    429     liststats_pos (Y, dY, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
     231    // use bootstrap resampling to check the error distribution
     232    // if we only have one point, this is silly...
     233   
     234    if (fitStats->NfitAlloc == 1) {
     235      FitAstromResultSetPM (&fitPos, 1, average);
     236      FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
     237    } else {
     238      fitStats->Nfit = 0;
     239      FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
     240      for (k = 0; k < fitStats->NfitAlloc; k++) {
     241        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     242        FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints);
     243        fitStats->Nfit ++;
     244      }
     245      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     246      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     247    }
     248    FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
    430249
    431250    // project Ro, Do back to RA,DEC
    432     XY_to_RD (&fitAve.Ro, &fitAve.Do, statsR.mean, statsD.mean, &coords);
    433     if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);
    434 
    435     fitAve.dRo = statsR.sigma;
    436     fitAve.dDo = statsD.sigma;
    437 
    438     fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
    439     fitAve.Nfit = N;
    440 
    441     fitAve.uR = fitAve.duR = 0.0;
    442     fitAve.uD = fitAve.duD = 0.0;
    443     fitAve.p  = fitAve.dp  = 0.0;
     251    XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
    444252    average[0].flags |= ID_STAR_FIT_AVE;
    445253    fitStats->Nave ++;
    446254  }
    447255
     256  // update the bit flags of which points were used
     257  for (k = 0; k < fitStats->Npoints; k++) {
     258    int Nm = fitStats->points[k].measure;
     259    myAssert (Nm >= 0, "oops");
     260    measure[Nm].dbFlags |= ID_MEAS_USED_OBJ;
     261    if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; }
     262  }
     263
     264  // we can set the star reference-image color only if we have loaded the image data
     265  int setRefColor = areImagesMatched();
    448266  if (setRefColor) {
     267    float *C_blue = NULL;
     268    float *C_red = NULL;
     269    ALLOCATE (C_blue, float, fitStats->Npoints);
     270    ALLOCATE (C_red, float, fitStats->Npoints);
     271
     272    int NcBlue = 0;
     273    int NcRed = 0;
     274
     275    for (k = 0; k < fitStats->Npoints; k++) {
     276      int Nm = fitStats->points[k].measure;
     277      float colorBlue = getColorBlue (measOff + Nm, cat);
     278      if (!isnan(colorBlue)) {
     279        C_blue[NcBlue] = colorBlue;
     280        NcBlue++;
     281      }
     282      float colorRed = getColorRed (measOff + Nm, cat);
     283      if (!isnan(colorRed)) {
     284        C_red[NcRed] = colorRed;
     285        NcRed++;
     286      }
     287    }
     288
     289    // need to reassign here if isfinite()
    449290    float colorMedian;
    450     dsort (C_blue, NcBlue);
     291    fsort (C_blue, NcBlue);
    451292    colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN;
    452293    average[0].refColorBlue = colorMedian;
    453     dsort (C_red, NcRed);
     294    fsort (C_red, NcRed);
    454295    colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN;
    455296    average[0].refColorRed = colorMedian;
     297
     298    free (C_blue);
     299    free (C_red);
    456300  }
    457301
    … …  
    459303  // XXXX for now, just use the mode as the result:
    460304  int result = mode;
     305  FitAstromResult fit;
     306  FitAstromResultInit (&fit);
    461307
    462308  switch (result) {
    463309    case FIT_AVERAGE:
    464310      average[0].flags |= ID_STAR_USE_AVE;
    465       fit = fitAve;
     311      fit = fitPos;
    466312      break;
    467313    case FIT_PM_ONLY:
    … …  
    471317    case FIT_PM_AND_PAR:
    472318      average[0].flags |= ID_STAR_USE_PAR;
    473       fit = fitPAR;
     319      fit = fitPar;
    474320      break;
    475321  }
    … …  
    503349
    504350  // what is the offset relative to the mean fit position?
    505   coords.crval1 = average[0].R;
    506   coords.crval2 = average[0].D;
    507   if (isnan(coords.crval1)) {
     351  fitStats->coords.crval1 = average[0].R;
     352  fitStats->coords.crval2 = average[0].D;
     353  if (isnan(fitStats->coords.crval1)) {
    508354    return (FALSE);
    509355  }
    510   if (isnan(coords.crval2)) {
     356  if (isnan(fitStats->coords.crval2)) {
    511357    return (FALSE);
    512358  }
    513359
    514360  double dXoff, dYoff;
    515   RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords);
     361  RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &fitStats->coords);
    516362  float dPos = hypot (dXoff, dYoff);
    517363  if (dPos > MaxMeanOffset) {
    518364    if (fitStats->Noffset < 100) {
    519       fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
     365      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
    520366    }
    521367    fitStats->Noffset ++;
    … …  
    531377                      average[0].uR,
    532378                      average[0].uD,
    533                       fitAve.chisq, fitPM.chisq, fitPAR.chisq);
     379                      fitPos.chisq, fitPM.chisq, fitPar.chisq);
    534380
    535381  average[0].R          = fit.Ro; // RA in degrees
    … …  
    546392  average[0].dP         = fit.dp; // parallax error in arcsec
    547393
    548   average[0].ChiSqAve   = fitAve.chisq;
     394  average[0].ChiSqAve   = fitPos.chisq;
    549395  average[0].ChiSqPM    = fitPM.chisq;
    550   average[0].ChiSqPar   = fitPAR.chisq;
    551 
    552   average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
     396  average[0].ChiSqPar   = fitPar.chisq;
     397
     398  average[0].Tmean      = (Tmean * 86400 * 365.25) + fitStats->T2000;
    553399  average[0].Trange     = (Trange * 86400 * 365.25);
    554400  average[0].Npos       = fit.Nfit;
    … …  
    563409// be called with just MeasureTiny set and Measure == NULL
    564410int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) {
    565 
    566   off_t k;
    567411
    568412  // set the default values
    … …  
    573417
    574418  /* calculate the average value of R,D for a single star */
    575   PMFit fitAve;
    576   memset (&fitAve, 0, sizeof(fitAve));
    577   fitAve.chisq = NAN;
     419  FitAstromResult fitPos;
     420  FitAstromResultInit (&fitPos);
    578421
    579422  if (average[0].Nmeasure == 0) return TRUE;
    580 
    581   int N = 0;
    582423
    583424  int XVERB = FALSE;
    … …  
    585426  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
    586427
     428  // select the measurements to be used in this analysis
     429  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE);
     430
     431  // too few measurements for average position (require 2 values)
     432  if (fitStats->Npoints < 1) return FALSE; // XXX ??
     433 
     434  double Tmean, Trange, parRange;
     435  FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
     436
     437  FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
     438  FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
     439
     440  // project Ro, Do back to RA,DEC
     441  XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
     442
     443  // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE;
     444  fitStats->Nave ++;
     445
     446  if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
     447                      average[0].R,
     448                      average[0].D,
     449                      fitPos.Ro, fitPos.Do,
     450                      3600*(average[0].R - fitPos.Ro),
     451                      3600*(average[0].D - fitPos.Do));
     452
     453  // make sure that the fit succeeded
     454  int status = TRUE;
     455  status &= finite(fitPos.Ro);
     456  status &= finite(fitPos.Do);
     457  status &= finite(fitPos.dRo);
     458  status &= finite(fitPos.dDo);
     459  if (!status) {
     460    fitStats->Nskip ++;
     461    return FALSE;
     462  }
     463
     464  // what is the offset relative to the mean fit position?
     465  fitStats->coords.crval1 = average[0].R;
     466  fitStats->coords.crval2 = average[0].D;
     467
     468  double dXoff, dYoff;
     469  RD_to_XY (&dXoff, &dYoff, fitPos.Ro, fitPos.Do, &fitStats->coords);
     470  float dPos = hypot (dXoff, dYoff);
     471  if (dPos > MaxMeanOffset) {
     472    if (fitStats->Noffset < 100) {
     473      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitPos.Ro, fitPos.Do, dXoff, dYoff);
     474    }
     475    fitStats->Noffset ++;
     476    return FALSE;
     477  }
     478
     479  // set the stack position values
     480  average[0].Rstk  = fitPos.Ro; // RA in degrees
     481  average[0].Dstk  = fitPos.Do; // DEC in degrees
     482  average[0].dRstk = fitPos.dRo; // RA scatter in arcsec
     483  average[0].dDstk = fitPos.dDo; // DEC scatter in arcsec
     484
     485  return (TRUE);
     486}
     487
     488int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int isStack) {
     489
     490  // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries
     491
     492  int has2MASS = FALSE;
     493
     494  int Npoints = fit->Npoints = 0;
     495  FitAstromPoint *points = fit->points;
     496
     497  int TESTPT2 = FALSE;
     498  TESTPT2 |= CAT_ID_SRC && OBJ_ID_SRC && (average[0].catID == CAT_ID_SRC) && (average[0].objID == OBJ_ID_SRC);
     499  TESTPT2 |= CAT_ID_DST && OBJ_ID_DST && (average[0].catID == CAT_ID_DST) && (average[0].objID == OBJ_ID_DST);
     500  if (TESTPT2) {
     501    fprintf (stderr, "got test det\n");
     502  }
     503
    587504  // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
     505  off_t k;
    588506  for (k = 0; k < average[0].Nmeasure; k++) {
    589507
    590     if (XVERB) {
     508    if (0) {
    591509      char *date = ohana_sec_to_date (measure[k].t);
    592510      int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
    593       fprintf (stderr, "stack: "OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
     511      fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
    594512      free (date);
    595513    }
    596514
    597     // SKIP everything except gpc1 stack data
    598     if (!isGPC1stack(measure[k].photcode)) continue;
    599 
     515    // SKIP gpc1 forced-warp data
     516    if (isGPC1warp(measure[k].photcode)) continue;
     517
     518    // SKIP gpc1 stack data
     519    if (isStack) {
     520      if (!isGPC1stack(measure[k].photcode)) continue;
     521    } else {
     522      if ( isGPC1stack(measure[k].photcode)) continue;
     523    }
     524
     525    // reset the bit to note that a detection was used (or not)
     526    measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
     527    if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
     528
     529    // does the measurement pass the supplied filtering constraints?
     530    // MeasFilterTestTiny does not test psfQF
    600531    // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
    601532    int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
    … …  
    604535    }
    605536
    606     R[N] = getMeanR (&measure[k], average, secfilt);
    607     D[N] = getMeanD (&measure[k], average, secfilt);
    608 
    609     // dX, dY : error in arcsec --
    610     dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
    611     dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
     537    double Ri = measure[k].R;
     538    double Di = measure[k].D;
     539
     540    // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
     541    if (ExcludeBogus) {
     542      fit->coords.crval1 = average[0].R;
     543      fit->coords.crval2 = average[0].D;
     544      double Xi, Yi;
     545      RD_to_XY (&Xi, &Yi, Ri, Di, &fit->coords);
     546      double radius = hypot(Xi, Yi);
     547      if (radius > ExcludeBogusRadius) {
     548        measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
     549        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
     550        continue;
     551      }
     552      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     553      if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
     554    }
     555
     556    // outlier rejection
     557    if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
     558      continue;
     559    }
     560
     561    FitAstromPointInit (&points[Npoints]);
     562
     563    points[Npoints].R = Ri;
     564    points[Npoints].D = Di;
     565
     566    // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
     567    // T[] is time in years since J2000 (jd = 2451545)
     568    points[Npoints].T = (measure[k].t - fit->T2000) / (86400*365.25) ; // time relative to J2000 in years
     569
     570    // add measured systematic error in quadrature?  only do this after the fit has
     571    // converged (or you will never improve the poor images)
     572
     573    // dX,dY are the X and Y direction errors in arcseconds.  dR, dD are the errors in
     574    // those directions in degrees.  IF we have non-circular errors (different values for
     575    // X and Y), then dR and dD will be incorrect: they would need to be rotated to take
     576    // out the position angle
     577
     578    // dX, dY : error in arcsec:
     579    points[Npoints].dX = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
     580    points[Npoints].dY = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
    612581
    613582    // allow a given photcode or measurement to be
    614583    // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
    615     if (isnan(dX[N])) continue;
    616     if (isnan(dY[N])) continue;
    617 
    618     // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
    619     // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
    620     // non-circular errors (different values for X and Y), then dR and dD will be
    621     // incorrect: they would need to be rotated to take out the position angle
    622     dR[N] = dX[N] / 3600.0;
    623     dD[N] = dY[N] / 3600.0;
    624 
    625     // XXX use a different flag for stack measurements?
    626     // measure[k].dbFlags |= ID_MEAS_USED_OBJ;
    627     // if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
    628 
    629     N++;
     584    if (isnan(points[Npoints].dX)) continue;
     585    if (isnan(points[Npoints].dY)) continue;
     586
     587    points[Npoints].dT = measure[k].dt;
     588
     589    points[Npoints].measure = k;
     590    Npoints++;
     591
     592    if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) {
     593      has2MASS = TRUE;
     594    }
     595
     596    myAssert (Npoints <= fit->NpointsAlloc, "oops");
    630597  } // loop over measurements : average[0].Nmeasure
    631598
    632   // if we have too few good detections for the desired fit, or too limited a
    633   // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
    634   // for even the average position, consider including the lower-quality detections?
    635 
    636   // too few measurements for average position (require 2 values)
    637   if (N < 1) return FALSE; // XXX ??
    638 
    639   // find the mean position
    640   StatType statsR, statsD;
    641   liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
    642   liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
    643 
    644   fitAve.Ro = statsR.mean;
    645   fitAve.dRo = 3600.0*statsR.sigma;
    646 
    647   fitAve.Do = statsD.mean;
    648   fitAve.dDo = 3600.0*statsD.sigma;
    649 
    650   fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq);
    651   fitAve.Nfit = N;
    652 
    653   // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE;
    654   fitStats->Nave ++;
    655 
    656   if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
    657                       average[0].R,
    658                       average[0].D,
    659                       fitAve.Ro, fitAve.Do,
    660                       3600*(average[0].R - fitAve.Ro),
    661                       3600*(average[0].D - fitAve.Do));
    662 
    663   // make sure that the fit succeeded
    664   int status = TRUE;
    665   status &= finite(fitAve.Ro);
    666   status &= finite(fitAve.Do);
    667   status &= finite(fitAve.dRo);
    668   status &= finite(fitAve.dDo);
    669   if (!status) {
    670     fitStats->Nskip ++;
    671     return FALSE;
    672   }
    673 
    674   // what is the offset relative to the mean fit position?
    675   coords.crval1 = average[0].R;
    676   coords.crval2 = average[0].D;
    677 
    678   double dXoff, dYoff;
    679   RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &coords);
    680   float dPos = hypot (dXoff, dYoff);
    681   if (dPos > MaxMeanOffset) {
    682     if (fitStats->Noffset < 100) {
    683       fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff);
    684     }
    685     fitStats->Noffset ++;
    686     return FALSE;
    687   }
    688 
    689   // set the stack position values
    690   average[0].Rstk  = fitAve.Ro; // RA in degrees
    691   average[0].Dstk  = fitAve.Do; // DEC in degrees
    692   average[0].dRstk = fitAve.dRo; // RA scatter in arcsec
    693   average[0].dDstk = fitAve.dDo; // DEC scatter in arcsec
    694 
     599  int TESTPT = FALSE;
     600  TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (average[0].catID == CAT_ID_SRC) && (average[0].objID == OBJ_ID_SRC);
     601  TESTPT |= CAT_ID_DST && OBJ_ID_DST && (average[0].catID == CAT_ID_DST) && (average[0].objID == OBJ_ID_DST);
     602  if (TESTPT) {
     603    fprintf (stderr, "got test det\n");
     604  }
     605 
     606  // XXX flag measurements from stars with 2MASS
     607  for (k = 0; k < average[0].Nmeasure; k++) {
     608    // reset the bit to note that a detection was used (or not)
     609    if (has2MASS) {
     610      measure[k].dbFlags |=  ID_MEAS_OBJECT_HAS_2MASS;
     611    } else {
     612      measure[k].dbFlags &= ~ID_MEAS_OBJECT_HAS_2MASS;
     613    }
     614  }
     615
     616  fit->Npoints = Npoints;
     617  return TRUE;
     618}
     619
     620int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) {
     621
     622  int k;
     623
     624  int Npoints = fitStats->Npoints;
     625  FitAstromPoint *points = fitStats->points;
     626
     627  // find Tmin & Tmax from the list of accepted measurements
     628  double Tmin  = points[0].T;
     629  double Tmax  = points[0].T;
     630  double pRmin = +2.0;
     631  double pRmax = -2.0;
     632  double pDmin = +2.0;
     633  double pDmax = -2.0;
     634
     635  *Tmean = 0.0;
     636
     637  double Tsum = 0.0;
     638  double Wsum = 0.0;
     639  for (k = 0; k < Npoints; k++) {
     640    Tmin = MIN(Tmin, points[k].T);
     641    Tmax = MAX(Tmax, points[k].T);
     642
     643    float wx = 1.0 / SQ(points[k].dX);
     644
     645    Tsum += points[k].T * wx;
     646    Wsum += wx;
     647
     648    // at this point, T is in years since J2000
     649    ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T);
     650    pRmin = MIN (pRmin, points[k].pR);
     651    pRmax = MAX (pRmax, points[k].pR);
     652    pDmin = MIN (pDmin, points[k].pD);
     653    pDmax = MAX (pDmax, points[k].pD);
     654  }
     655  *Trange = Tmax - Tmin;
     656
     657  // mean epoch
     658  *Tmean = Tsum / Wsum;
     659
     660  // for HIGH_SPEED, just use the center of the range
     661  if (RELASTRO_OP == OP_HIGH_SPEED) {
     662    *Tmean = 0.5*(Tmax - Tmin);
     663  }
     664
     665  *parRange = hypot (pRmax - pRmin, pDmax - pDmin);
     666
     667  /* we need to do the fit in a locally linear space; choose a ref coordinate */
     668  fitStats->coords.crval1 = points[0].R;
     669  fitStats->coords.crval2 = points[0].D;
     670
     671  // project all of the R,D coordinates to a plane centered on this coordinate. set
     672  // the times to be relative to Tmean
     673  for (k = 0; k < Npoints; k++) {
     674    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords);
     675    points[k].T -= *Tmean;
     676  }       
     677  return TRUE;
     678}
     679
     680int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) {
     681
     682  int i, j;
     683
     684  int Nmax = 0;
     685  for (i = 0; i < Ncatalog; i++) {
     686    for (j = 0; j < catalog[i].Naverage; j++) {
     687      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
     688    }
     689  }
     690  return Nmax;
     691}
     692
     693int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints) {
     694  int i;
     695
     696  // I need to draw Npoints random entries from 'points' with replacement:
     697  for (i = 0; i < Npoints; i++) {
     698    int N = Npoints * drand48();
     699    sample[i] = points[N];
     700  }
     701  return TRUE;
     702}
     703
     704// calculate mean and sigma points for the 5 fit parameter
     705int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) {
     706
     707  // generate a histogram for the selected element
     708  double *values = NULL;
     709  ALLOCATE (values, double, Nfit);
     710 
     711  int i;
     712
     713  for (i = 0; i < Nfit; i++) {
     714    switch (mode) {
     715      case FIT_RESULT_RA:
     716        values[i] = fit[i].Ro;
     717        break;
     718      case FIT_RESULT_DEC:
     719        values[i] = fit[i].Do;
     720        break;
     721      case FIT_RESULT_uR:
     722        values[i] = fit[i].uR;
     723        break;
     724      case FIT_RESULT_uD:
     725        values[i] = fit[i].uD;
     726        break;
     727      case FIT_RESULT_PLX:
     728        values[i] = fit[i].p;
     729        break;
     730      default:
     731        myAbort ("invalid option");
     732    }
     733  }
     734
     735  dsort (values, Nfit);
     736
     737  double median;
     738  if (Nfit % 2) {
     739    int Ncenter = Nfit / 2;
     740    median = values[Ncenter];
     741  } else {
     742    int Ncenter = Nfit / 2 - 1;
     743    median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
     744  }
     745
     746  double Slo = VectorFractionInterpolate (values, 0.158655, Nfit);
     747  double Shi = VectorFractionInterpolate (values, 0.841345, Nfit);
     748  double sigma = (Shi - Slo) / 2.0;
     749
     750  switch (mode) {
     751    case FIT_RESULT_RA:
     752      result->Ro = median;
     753      result->dRo = sigma;
     754      break;
     755    case FIT_RESULT_DEC:
     756      result->Do = median;
     757      result->dDo = sigma;
     758      break;
     759    case FIT_RESULT_uR:
     760      result->uR = median;
     761      result->duR = sigma;
     762      break;
     763    case FIT_RESULT_uD:
     764      result->uD = median;
     765      result->duD = sigma;
     766      break;
     767    case FIT_RESULT_PLX:
     768      result->p = median;
     769      result->dp = sigma;
     770      break;
     771    default:
     772      myAbort ("invalid option");
     773  }
     774
     775  return TRUE;
     776}
     777
     778double VectorFractionInterpolate (double *values, float fraction, int Npts) {
     779
     780  float F = fraction * Npts;
     781  int   N = fraction * Npts;
     782
     783  if (N < 0        ) return NAN;
     784  if (N >= Npts - 2) return NAN;
     785
     786  // interpolate between N,N+1
     787   
     788  double S = (F - N) * (values[N+1] - values[N]) + values[N];
     789  return S;
     790}
     791
     792int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) {
     793
     794  int i;
     795
     796  // add up the chi square for the fit
     797  double chisq = 0.0;
     798  for (i = 0; i < Npoints; i++) {
     799    double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR;
     800    double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD;
     801    chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
     802    chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
     803  }
     804  switch (mode) {
     805    case FIT_AVERAGE:
     806      fit->chisq = chisq / (2.0*Npoints - 2.0);
     807      break;
     808    case FIT_PM_ONLY:
     809      fit->chisq = chisq / (2.0*Npoints - 4.0);
     810      break;
     811    case FIT_PM_AND_PAR:
     812      fit->chisq = chisq / (2.0*Npoints - 5.0);
     813      break;
     814    default:
     815      myAbort ("invalid mode");
     816  }
     817  fit->Nfit = Npoints;
    695818  return (TRUE);
    696819}
    697820
    698 
    699 
    700 /* fitting proper-motion and parallax:
    701 
    702    given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
    703 
    704    let x,y be the coordinate in the linearized frame with y parallel to DEC lines
    705 
    706    L,B are the ecliptic longitude and latitude of the object,
    707    dL and dB are the offsets in the L and B directions
    708 
    709    dL = sin(t - topp)
    710    dB = cos(t - topp)*sin(B)
    711 
    712    these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
    713    for the source in the x,y frame is:
    714 
    715    x = Ro + uR * (t - to) + p * pR
    716    y = Do + uD * (t - to) + p * pD
    717 
    718    the unknowns in these equations are Ro, uR, Do, uD, and p
    719 
    720    XXX think through the concepts for the pole a bit better.  all objects near the pole
    721    move the same way with the same phase.  choose a projection center and define dL,dB relative
    722    to that center point coordinate system?
    723 
    724 */
     821int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average) {
     822
     823  int i;
     824
     825  if (USE_GALAXY_MODEL) {
     826    for (i = 0; i < Nfit; i++) {
     827      fit->uR = average->uRgal;
     828      fit->uD = average->uDgal;
     829    }
     830  } else {
     831    for (i = 0; i < Nfit; i++) {
     832      fit->uR = 0.0;
     833      fit->uD = 0.0;
     834    }
     835  }
     836
     837  return TRUE;
     838}
Note: See TracChangeset for help on using the changeset viewer.