IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38612


Ignore:
Timestamp:
Jul 22, 2015, 8:13:53 AM (11 years ago)
Author:
eugene
Message:

relastro builds

Location:
branches/eam_branches/ipp-20150625/Ohana/src/relastro
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r38605 r38612  
    2626$(SRC)/FitChip.$(ARCH).o             \
    2727$(SRC)/FitMosaic.$(ARCH).o           \
     28$(SRC)/FitSimple.$(ARCH).o           \
     29$(SRC)/FitAstromOps.$(ARCH).o           \
    2830$(SRC)/FitPM.$(ARCH).o               \
    2931$(SRC)/FitPMandPar.$(ARCH).o         \
    3032$(SRC)/FitPosPMfixed.$(ARCH).o       \
    31 $(SRC)/FitSimple.$(ARCH).o           \
    3233$(SRC)/ImageOps.$(ARCH).o            \
    3334$(SRC)/MosaicOps.$(ARCH).o           \
    … …  
    99100RELASTRO_CLIENT = \
    100101$(SRC)/ConfigInit.$(ARCH).o          \
     102$(SRC)/FitSimple.$(ARCH).o           \
     103$(SRC)/FitAstromOps.$(ARCH).o           \
    101104$(SRC)/FitPM.$(ARCH).o               \
    102 $(SRC)/FitPar.$(ARCH).o              \
    103105$(SRC)/FitPMandPar.$(ARCH).o         \
    104 $(SRC)/FitSimple.$(ARCH).o           \
     106$(SRC)/FitPosPMfixed.$(ARCH).o       \
    105107$(SRC)/ImageOps.$(ARCH).o            \
    106108$(SRC)/MosaicOps.$(ARCH).o           \

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

    r38605 r38612  
    147147  double C_blue;
    148148  double C_red;
     149  int measure;
    149150} FitAstromPoint;
    150151
    … …  
    709710int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog);
    710711int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange);
    711 int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff);
     712int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int isStack);
    712713int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
    713714int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff);

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c

    r38605 r38612  
    44
    55FitStats *FitStatsInit (int Nmax, int Nboot) {
     6
     7  FitStats *fitStats = NULL;
     8  ALLOCATE (fitStats, FitStats, 1);
    69
    710  // counters to record successful fits or failures
    … …  
    2730  fitStats->NpointsAlloc = Nmax;
    2831  if (Nmax > 0) {
    29     ALLOCATE (fitStats->points, FitAstromPoints, Nmax);
    30     ALLOCATE (fitStats->sample, FitAstromPoints, Nmax);
     32    ALLOCATE (fitStats->points, FitAstromPoint, Nmax);
     33    ALLOCATE (fitStats->sample, FitAstromPoint, Nmax);
    3134  }
    3235
    … …  
    4851  // use J2000 as a reference time
    4952  fitStats->T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
    50   return;
     53  return fitStats;
    5154}
    5255
    … …  
    107110}
    108111
    109 void FitAstromPointInit (FitAstromObject *object) {
     112void FitAstromPointInit (FitAstromPoint *object) {
    110113  object->X      = 0.0;
    111114  object->Y      = 0.0;
    … …  
    122125  object->C_blue = 0.0;
    123126  object->C_red  = 0.0;
     127  object->measure= -1;
    124128  return;
    125129}

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c

    r38604 r38612  
    11# include "relastro.h"
    22
    3 int FitPMandPAR (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     3int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
    44
    55  int i;

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c

    r38606 r38612  
    6666  int k;
    6767
    68   int setRefColor = areImagesMatched();
    69 
    7068  /* calculate the average value of R,D for a single star */
    7169
    … …  
    9290
    9391  // select the measurements to be used in this analysis
    94   UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, cat, measOff);
     92  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, cat, measOff, FALSE);
    9593
    9694  // if there are no exposure detections, use the stack position
    … …  
    107105
    108106  double Tmean, Trange, parRange;
    109   FitAstromPoints_Project (fitStats->points, fitStats->Npoints, &Tmean, &Trange, &parRange);
     107  FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
    110108
    111109  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
    … …  
    120118      goto justPosition;
    121119    }
    122     if (Npoints <= PM_TOOFEW) {
     120    if (fitStats->Npoints <= PM_TOOFEW) {
    123121      mode = FIT_AVERAGE;
    124122      goto justPosition;
    … …  
    161159  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
    162160  if (mode == FIT_PM_AND_PAR) {
    163     if (N <= PAR_TOOFEW) {
     161    if (parRange < PAR_FACTOR_MIN) {
    164162      mode = FIT_PM_ONLY;
    165163      goto justPosition;
    166164    }
    167     if (parRange < PAR_FACTOR_MIN) {
     165    if (fitStats->Npoints <= PAR_TOOFEW) {
    168166      mode = FIT_PM_ONLY;
    169167      goto justPosition;
    … …  
    220218  for (k = 0; k < fitStats->Npoints; k++) {
    221219    int Nm = fitStats->points[k].measure;
     220    myAssert (Nm >= 0, "oops");
    222221    measure[Nm].dbFlags |= ID_MEAS_USED_OBJ;
    223222    if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; }
    … …  
    225224
    226225# if (0)
     226  int setRefColor = areImagesMatched();
    227227  if (setRefColor) {
    228228    // need to reassign here if isfinite()
    … …  
    345345int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) {
    346346
    347   off_t k;
    348 
    349347  // set the default values
    350348  average[0].Rstk  = NAN; // RA in degrees
    … …  
    354352
    355353  /* calculate the average value of R,D for a single star */
    356   PMFit fitAve;
    357   memset (&fitAve, 0, sizeof(fitAve));
    358   fitAve.chisq = NAN;
     354  FitAstromResult fitPos;
     355  FitAstromResultInit (&fitPos);
    359356
    360357  if (average[0].Nmeasure == 0) return TRUE;
    361 
    362   int N = 0;
    363358
    364359  int XVERB = FALSE;
    … …  
    366361  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
    367362
    368   // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
    369   for (k = 0; k < average[0].Nmeasure; k++) {
    370 
    371     if (XVERB) {
    372       char *date = ohana_sec_to_date (measure[k].t);
    373       int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
    374       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);
    375       free (date);
    376     }
    377 
    378     // SKIP everything except gpc1 stack data
    379     if (!isGPC1stack(measure[k].photcode)) continue;
    380 
    381     // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
    382     int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
    383     if (!keepMeasure) {
    384       continue;
    385     }
    386 
    387     R[N] = getMeanR (&measure[k], average, secfilt);
    388     D[N] = getMeanD (&measure[k], average, secfilt);
    389 
    390     // dX, dY : error in arcsec --
    391     dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
    392     dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
    393 
    394     // allow a given photcode or measurement to be
    395     // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
    396     if (isnan(dX[N])) continue;
    397     if (isnan(dY[N])) continue;
    398 
    399     // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
    400     // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
    401     // non-circular errors (different values for X and Y), then dR and dD will be
    402     // incorrect: they would need to be rotated to take out the position angle
    403     dR[N] = dX[N] / 3600.0;
    404     dD[N] = dY[N] / 3600.0;
    405 
    406     // XXX use a different flag for stack measurements?
    407     // measure[k].dbFlags |= ID_MEAS_USED_OBJ;
    408     // if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
    409 
    410     N++;
    411   } // loop over measurements : average[0].Nmeasure
    412 
    413   // if we have too few good detections for the desired fit, or too limited a
    414   // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
    415   // for even the average position, consider including the lower-quality detections?
     363  // select the measurements to be used in this analysis
     364  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, 0, 0, TRUE);
    416365
    417366  // too few measurements for average position (require 2 values)
    418   if (N < 1) return FALSE; // XXX ??
    419 
    420   // find the mean position
    421   StatType statsR, statsD;
    422   liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
    423   liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
    424 
    425   fitAve.Ro = statsR.mean;
    426   fitAve.dRo = 3600.0*statsR.sigma;
    427 
    428   fitAve.Do = statsD.mean;
    429   fitAve.dDo = 3600.0*statsD.sigma;
    430 
    431   fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq);
    432   fitAve.Nfit = N;
     367  if (fitStats->Npoints < 1) return FALSE; // XXX ??
     368 
     369  double Tmean, Trange, parRange;
     370  FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
     371
     372  FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
     373  FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
     374
     375  // project Ro, Do back to RA,DEC
     376  XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
    433377
    434378  // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE;
    … …  
    438382                      average[0].R,
    439383                      average[0].D,
    440                       fitAve.Ro, fitAve.Do,
    441                       3600*(average[0].R - fitAve.Ro),
    442                       3600*(average[0].D - fitAve.Do));
     384                      fitPos.Ro, fitPos.Do,
     385                      3600*(average[0].R - fitPos.Ro),
     386                      3600*(average[0].D - fitPos.Do));
    443387
    444388  // make sure that the fit succeeded
    445389  int status = TRUE;
    446   status &= finite(fitAve.Ro);
    447   status &= finite(fitAve.Do);
    448   status &= finite(fitAve.dRo);
    449   status &= finite(fitAve.dDo);
     390  status &= finite(fitPos.Ro);
     391  status &= finite(fitPos.Do);
     392  status &= finite(fitPos.dRo);
     393  status &= finite(fitPos.dDo);
    450394  if (!status) {
    451395    fitStats->Nskip ++;
    … …  
    458402
    459403  double dXoff, dYoff;
    460   RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &fitStats->coords);
     404  RD_to_XY (&dXoff, &dYoff, fitPos.Ro, fitPos.Do, &fitStats->coords);
    461405  float dPos = hypot (dXoff, dYoff);
    462406  if (dPos > MaxMeanOffset) {
    463407    if (fitStats->Noffset < 100) {
    464       fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff);
     408      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitPos.Ro, fitPos.Do, dXoff, dYoff);
    465409    }
    466410    fitStats->Noffset ++;
    … …  
    469413
    470414  // set the stack position values
    471   average[0].Rstk  = fitAve.Ro; // RA in degrees
    472   average[0].Dstk  = fitAve.Do; // DEC in degrees
    473   average[0].dRstk = fitAve.dRo; // RA scatter in arcsec
    474   average[0].dDstk = fitAve.dDo; // DEC scatter in arcsec
     415  average[0].Rstk  = fitPos.Ro; // RA in degrees
     416  average[0].Dstk  = fitPos.Do; // DEC in degrees
     417  average[0].dRstk = fitPos.dRo; // RA scatter in arcsec
     418  average[0].dDstk = fitPos.dDo; // DEC scatter in arcsec
    475419
    476420  return (TRUE);
    477421}
    478422
    479 int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff) {
     423int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int isStack) {
    480424
    481425  // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries
    … …  
    495439    }
    496440
    497     // SKIP gpc1 stack data
    498     if (isGPC1stack(measure[k].photcode)) continue;
    499 
    500441    // SKIP gpc1 forced-warp data
    501442    if (isGPC1warp(measure[k].photcode)) continue;
     443
     444    // SKIP gpc1 stack data
     445    if (isStack) {
     446      if (!isGPC1stack(measure[k].photcode)) continue;
     447    } else {
     448      if ( isGPC1stack(measure[k].photcode)) continue;
     449    }
    502450
    503451    // reset the bit to note that a detection was used (or not)
    … …  
    518466    // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
    519467    if (ExcludeBogus) {
    520       fitStats->coords.crval1 = average[0].R;
    521       fitStats->coords.crval2 = average[0].D;
     468      fit->coords.crval1 = average[0].R;
     469      fit->coords.crval2 = average[0].D;
    522470      double Xi, Yi;
    523       RD_to_XY (&Xi, &Yi, Ri, Di, &fitStats->coords);
     471      RD_to_XY (&Xi, &Yi, Ri, Di, &fit->coords);
    524472      double radius = hypot(Xi, Yi);
    525473      if (radius > ExcludeBogusRadius) {
    … …  
    565513    points[Npoints].dT = measure[k].dt;
    566514
     515    // XXX
     516    int setRefColor = 0;
    567517    if (setRefColor) {
    568518      points[Npoints].C_blue = getColorBlue (measOff+k, cat);
    … …  
    584534  int k;
    585535
     536  int Npoints = fitStats->Npoints;
    586537  FitAstromPoint *points = fitStats->points;
    587538
    … …  
    598549    Tmin = MIN(Tmin, points[k].T);
    599550    Tmax = MAX(Tmax, points[k].T);
    600     Tmean += points[k].T;
     551    *Tmean += points[k].T;
    601552
    602553    // at this point, T is in years since J2000
    … …  
    627578  for (k = 0; k < Npoints; k++) {
    628579    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords);
    629     points[k].T -= Tmean;
     580    points[k].T -= *Tmean;
    630581  }       
    631582  return TRUE;
    … …  
    636587  int i, j;
    637588
    638   Nmax = 0;
     589  int Nmax = 0;
    639590  for (i = 0; i < Ncatalog; i++) {
    640591    for (j = 0; j < catalog[i].Naverage; j++) {
    … …  
    642593    }
    643594  }
    644 
    645595  return Nmax;
    646596}
    … …  
    664614  ALLOCATE (values, double, Nfit);
    665615 
     616  int i;
     617
    666618  for (i = 0; i < Nfit; i++) {
    667619    switch (mode) {
    … …  
    745697int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) {
    746698
     699  int i;
     700
    747701  // add up the chi square for the fit
    748702  double chisq = 0.0;
    … …  
    755709  switch (mode) {
    756710    case FIT_AVERAGE:
    757       fit->chisq = chisq / (2.0*Npts - 2.0);
     711      fit->chisq = chisq / (2.0*Npoints - 2.0);
    758712      break;
    759713    case FIT_PM_ONLY:
    760       fit->chisq = chisq / (2.0*Npts - 4.0);
     714      fit->chisq = chisq / (2.0*Npoints - 4.0);
    761715      break;
    762716    case FIT_PM_AND_PAR:
    763       fit->chisq = chisq / (2.0*Npts - 5.0);
     717      fit->chisq = chisq / (2.0*Npoints - 5.0);
    764718      break;
    765719    default:
Note: See TracChangeset for help on using the changeset viewer.