IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39375 for trunk


Ignore:
Timestamp:
Feb 23, 2016, 4:48:13 PM (10 years ago)
Author:
eugene
Message:

clarify rules for astrometry failures and small number cases; set stack based on median of stack points (not statistically independed); add stack tsets

Location:
trunk/Ohana/src/relastro
Files:
1 added
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r39370 r39375  
    1 default: relastro relastro_client
     1default: relastro relastro_client fitobj fitobj2 fitpm
    22help:
    33        @echo "make options: relastro (default)"
    … …  
    3333$(SRC)/FitPMandPar.$(ARCH).o         \
    3434$(SRC)/FitPosPMfixed.$(ARCH).o       \
     35$(SRC)/FitPosStack.$(ARCH).o       \
    3536$(SRC)/BootstrapOps.$(ARCH).o        \
    3637$(SRC)/ImageOps.$(ARCH).o            \
    … …  
    108109$(SRC)/FitPMandPar.$(ARCH).o         \
    109110$(SRC)/FitPosPMfixed.$(ARCH).o       \
     111$(SRC)/FitPosStack.$(ARCH).o       \
    110112$(SRC)/BootstrapOps.$(ARCH).o        \
    111113$(SRC)/ImageOps.$(ARCH).o            \
    … …  
    171173$(SRC)/FitPM.$(ARCH).o               \
    172174$(SRC)/FitPMandPar.$(ARCH).o               \
     175$(SRC)/FitPosStack.$(ARCH).o       \
    173176$(SRC)/BootstrapOps.$(ARCH).o        \
    174177$(SRC)/mkpolyterm.$(ARCH).o            \
    … …  
    187190$(SRC)/FitPM.$(ARCH).o               \
    188191$(SRC)/FitPMandPar.$(ARCH).o               \
     192$(SRC)/FitPosStack.$(ARCH).o       \
    189193$(SRC)/UpdateObjects.$(ARCH).o               \
    190194$(SRC)/BootstrapOps.$(ARCH).o        \
    … …  
    209213$(SRC)/FitPM.$(ARCH).o               \
    210214$(SRC)/FitPMandPar.$(ARCH).o               \
     215$(SRC)/FitPosStack.$(ARCH).o       \
    211216$(SRC)/UpdateObjects.$(ARCH).o               \
    212217$(SRC)/BootstrapOps.$(ARCH).o        \
    … …  
    221226$(FITOBJ2): $(INC)/relastro.h
    222227$(BIN)/fitobj2.$(ARCH): $(FITOBJ2)
     228
     229fitstk: $(BIN)/fitstk.$(ARCH)
     230
     231FITSTK = \
     232$(SRC)/fitstk.$(ARCH).o               \
     233$(SRC)/ParFactor.$(ARCH).o           \
     234$(SRC)/FitAstromOps.$(ARCH).o        \
     235$(SRC)/FitPosPMfixed.$(ARCH).o               \
     236$(SRC)/FitPM.$(ARCH).o               \
     237$(SRC)/FitPMandPar.$(ARCH).o               \
     238$(SRC)/FitPosStack.$(ARCH).o       \
     239$(SRC)/UpdateObjects.$(ARCH).o               \
     240$(SRC)/BootstrapOps.$(ARCH).o        \
     241$(SRC)/mkpolyterm.$(ARCH).o            \
     242$(SRC)/GetAstromError.$(ARCH).o            \
     243$(SRC)/ImageOps.$(ARCH).o            \
     244$(SRC)/MosaicOps.$(ARCH).o            \
     245$(SRC)/extra.$(ARCH).o            \
     246$(SRC)/plotstuff.$(ARCH).o            \
     247$(SRC)/fitpoly.$(ARCH).o
     248
     249$(FITSTK): $(INC)/relastro.h
     250$(BIN)/fitstk.$(ARCH): $(FITSTK)
    223251
    224252# $(SRC)/FitPM_IRLS.$(ARCH).o          \

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

    r39371 r39375  
    173173  int Nfit;
    174174  int NfitAlloc;
     175
     176  double * Xstack;
     177  double *dXstack;
     178  double * Ystack;
     179  double *dYstack;
    175180
    176181  FitAstromPoint *points;
    … …  
    755760int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average);
    756761void AstromErrorSetLoop (int Nloop, int isImageMode);
     762void FitPointsClearMasks (FitAstromPoint *points, int Npoints);
    757763
    758764int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    … …  
    760766int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    761767int FitPosPMfixed_Single (FitAstromResult *fit, FitAstromPoint *points, int Npoints);
     768int FitPosStack (FitAstromResult *fit, FitStats *fitstats);
    762769
    763770double MedianAbsDeviation(FitAstromPoint *points, int Npoints);

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

    r39371 r39375  
    8484  SetZeroPoint (25.0);
    8585
     86  if (USE_GALAXY_MODEL) {
     87    if (!InitGalaxyModel ("FEAST-HIPPARCOS")) {
     88      fprintf (stderr, "failed to init galaxy model\n");
     89      exit (2);
     90    }
     91  }
     92
    8693  free (config);
    8794  free (file);

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

    r39364 r39375  
    8585
    8686  // containers to hold the measurements for a given star
     87  fitStats-> Xstack = NULL; // pre-allocated array for median
     88  fitStats->dXstack = NULL; // pre-allocated array for median
     89  fitStats-> Ystack = NULL; // pre-allocated array for median
     90  fitStats->dYstack = NULL; // pre-allocated array for median
    8791  fitStats->points = NULL;
    8892  fitStats->sample = NULL;
    … …  
    9498    ALLOCATE (fitStats->sample, FitAstromPoint, Nmax);
    9599    ALLOCATE (fitStats->nomask, FitAstromPoint, Nmax);
     100    ALLOCATE (fitStats-> Xstack, double, Nmax);
     101    ALLOCATE (fitStats->dXstack, double, Nmax);
     102    ALLOCATE (fitStats-> Ystack, double, Nmax);
     103    ALLOCATE (fitStats->dYstack, double, Nmax);
    96104  }
    97105
    … …  
    143151  FREE (fitStats->nomask);
    144152
     153  FREE (fitStats-> Xstack);
     154  FREE (fitStats->dXstack);
     155  FREE (fitStats-> Ystack);
     156  FREE (fitStats->dYstack);
     157
    145158  FitAstromDataFree (fitStats->fitdataPos);
    146159  FitAstromDataFree (fitStats->fitdataPM);
    … …  
    168181
    169182  return fit;
     183}
     184
     185void FitPointsClearMasks (FitAstromPoint *points, int Npoints) {
     186
     187  int i;
     188  for (i = 0; i < Npoints; i++) {
     189    points[i].mask = FALSE;
     190  }
     191  return;
    170192}
    171193

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

    r39371 r39375  
    2626    fit->duD = 0.0;
    2727    fit->dp  = 0.0;
     28    fit->Nfit = 1;
    2829
    2930    found = i;

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

    r37807 r39375  
    6363
    6464    for (i = 0; i < catalog[c].Naverage; i++) {
    65       if (!(catalog[c].average[i].flags & ID_ICRF_QSO)) continue;
     65      if (!(catalog[c].average[i].flags & ID_OBJ_ICRF_QSO)) continue;
    6666
    6767      // only save a single value

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

    r39371 r39375  
    116116
    117117  // if we fail to fit the astrometry for some reason, we need to set/reset these
    118   average[0].flags |= ID_STAR_NO_ASTROM;
     118  average[0].flags |= ID_OBJ_NO_MEAN_ASTROM;
    119119  average[0].ChiSqAve  = NAN;
    120120  average[0].ChiSqPM   = NAN;
    … …  
    142142      average[0].dR = average[0].dRstk;
    143143      average[0].dD = average[0].dDstk;
    144       average[0].flags |= ID_STACK_ASTROM;
     144      average[0].flags |= ID_OBJ_STACK_FOR_MEAN;
    145145    }
    146146    return FALSE;
    … …  
    207207    if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    208208
    209     average[0].flags |= ID_STAR_FIT_PAR;
     209    average[0].flags |= ID_OBJ_FIT_PAR;
    210210    fitStats->Npar ++;
    211211
    … …  
    225225      mode = FIT_PM_ONLY;
    226226    } else {
    227       average[0].flags |= ID_STAR_USE_PAR;
     227      average[0].flags |= ID_OBJ_USE_PAR;
    228228    }
    229229  }       
    … …  
    242242    }
    243243
    244     if (average[0].flags & ID_STAR_USE_PAR) {
     244    if (average[0].flags & ID_OBJ_USE_PAR) {
    245245      // get the best pm fit and chisq given the set of points (mask is respected)
    246246      if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
    247         average[0].flags |= ID_STAR_BAD_PM;
     247        average[0].flags |= ID_OBJ_BAD_PM;
    248248        goto justPosition;
    249249      }
    … …  
    277277    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    278278
    279     average[0].flags |= ID_STAR_FIT_PM;
     279    average[0].flags |= ID_OBJ_FIT_PM;
    280280    fitStats->Npm ++;
    281281
    … …  
    290290    if (!valid) {
    291291      mode = FIT_AVERAGE;
    292       average[0].flags |= ID_STAR_BAD_PM;
     292      average[0].flags |= ID_OBJ_BAD_PM;
    293293    } else {
    294       if (!(average[0].flags & ID_STAR_USE_PAR)) {
    295         average[0].flags |= ID_STAR_USE_PM;
     294      if (!(average[0].flags & ID_OBJ_USE_PAR)) {
     295        average[0].flags |= ID_OBJ_USE_PM;
    296296      }
    297297    }
    … …  
    305305    FitAstromResultSetPM (&fitPos, 1, average);
    306306    // fprintf (stderr, "fit 1: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD);
    307     if (average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) {
     307    if (average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) {
    308308      if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     309        // if this fails, stick with the PM and/or PAR fit from above
    309310        goto doneWithFit;
    310311      }
    … …  
    313314        // I will not try to outlier-reject, just calculate the weighted average
    314315        if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     316          // if we have tried this, we have not masked any points; we will find a single unmasked point below
    315317          goto doneWithFit;
    316318        }
    317         // XXX Maybe add a flag here?
    318         average[0].flags |= ID_STAR_USE_AVE;
     319        average[0].flags |= ID_OBJ_USE_AVE;
     320        average[0].flags |= ID_OBJ_RAW_AVE;
    319321        goto useBasic;
    320322      }
    321323      if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    322         goto doneWithFit;
     324        // if the above fails, we need to clear the masks and try again below
     325        FitPointsClearMasks (fitStats->points, fitStats->Npoints);
     326        // just calculate the weighted average
     327        if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     328          // if this fails, find a single unmasked point below
     329          goto doneWithFit;
     330        }
     331        average[0].flags |= ID_OBJ_USE_AVE;
     332        average[0].flags |= ID_OBJ_RAW_AVE;
     333        goto useBasic;
    323334      }
    324335      // fprintf (stderr, "fit 2: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD);
    … …  
    327338        int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    328339        if (Nnomask < POS_MIN_NPTS_BOOT) {
     340          // if the above fails, we need to clear the masks and try again below
     341          FitPointsClearMasks (fitStats->points, fitStats->Npoints);
    329342          if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     343            // if this fails, find a single unmasked point below
    330344            goto doneWithFit;
    331345          }
    332           // XXX Maybe add a flag here?
    333           average[0].flags |= ID_STAR_USE_AVE;
     346          average[0].flags |= ID_OBJ_USE_AVE;
     347          average[0].flags |= ID_OBJ_RAW_AVE;
    334348          goto useBasic;
    335349        }
    … …  
    343357        BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    344358      }
    345       average[0].flags |= ID_STAR_USE_AVE;
     359      average[0].flags |= ID_OBJ_USE_AVE;
    346360    }
    347361
    … …  
    351365    if (fabs(fitPos.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    352366
    353     average[0].flags |= ID_STAR_FIT_AVE;
     367    average[0].flags |= ID_OBJ_FIT_AVE;
    354368    fitStats->Nave ++;
    355369  }
    … …  
    357371 doneWithFit:
    358372  // if no valid fit has been found, try to use a single unmasked point:
    359   if (!(average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM | ID_STAR_USE_AVE))) {
     373  if (!(average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM | ID_OBJ_USE_AVE))) {
    360374    if (!FitPosPMfixed_Single (&fitPos, fitStats->points, fitStats->Npoints)) {
    361375      fitStats->Nskip ++;
    … …  
    366380    if (fabs(fitPos.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    367381
    368     average[0].flags |= ID_STAR_USE_AVE;
    369     average[0].flags |= ID_STAR_FIT_AVE;
     382    average[0].flags |= ID_OBJ_USE_AVE;
     383    average[0].flags |= ID_OBJ_FIT_AVE;
    370384    fitStats->Nave ++;
    371385  }
    … …  
    422436  FitAstromResultInit (&fit);
    423437
    424   if (average[0].flags & ID_STAR_USE_PAR) {
    425     myAssert ((average[0].flags & (ID_STAR_USE_PM | ID_STAR_USE_AVE)) == 0, "programming error");
     438  if (average[0].flags & ID_OBJ_USE_PAR) {
     439    myAssert ((average[0].flags & (ID_OBJ_USE_PM | ID_OBJ_USE_AVE)) == 0, "programming error");
    426440    fit = fitPar;
    427441  }
    428   if (average[0].flags & ID_STAR_USE_PM) {
    429     myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_AVE)) == 0, "programming error");
     442  if (average[0].flags & ID_OBJ_USE_PM) {
     443    myAssert ((average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_AVE)) == 0, "programming error");
    430444    fit = fitPM;
    431445  }
    432   if (average[0].flags & ID_STAR_USE_AVE) {
    433     myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) == 0, "programming error");
     446  if (average[0].flags & ID_OBJ_USE_AVE) {
     447    myAssert ((average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) == 0, "programming error");
    434448    fit = fitPos;
    435449  }
    … …  
    516530
    517531  // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
    518   average[0].flags &= ~ID_STAR_NO_ASTROM;
     532  average[0].flags &= ~ID_OBJ_NO_MEAN_ASTROM;
    519533
    520534  return (TRUE);
    … …  
    554568    average[0].dRstk = average[0].dR; // RA scatter in arcsec
    555569    average[0].dDstk = average[0].dD; // DEC scatter in arcsec
     570    average[0].flags |= ID_OBJ_MEAN_FOR_STACK;
    556571  }
    557572
    … …  
    563578
    564579  // the stack positions are not statistically independent...
    565   FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
    566   FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
     580  FitPosStack (&fitPos, fitStats);
    567581
    568582  // project Ro, Do back to RA,DEC
    569583  XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
    570584
    571   // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE;
     585  // XXX choose stack flag? average[0].flags |= ID_OBJ_FIT_AVE;
    572586  fitStats->Nave ++;
    573587
    … …  
    610624  average[0].dRstk = fitPos.dRo; // RA scatter in arcsec
    611625  average[0].dDstk = fitPos.dDo; // DEC scatter in arcsec
     626  average[0].Tmean = (Tmean * 86400 * 365.25) + fitStats->T2000;
     627  average[0].flags &= ~ID_OBJ_MEAN_FOR_STACK;
    612628
    613629  return (TRUE);

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

    r39287 r39375  
    2525  int Nsecfilt;
    2626  Coords coords;
     27
     28  CoordTransform *transform = NULL;
     29  if (USE_GALAXY_MODEL) {
     30    transform = InitTransform (COORD_CELESTIAL, COORD_GALACTIC);
     31  }
    2732
    2833  /* for outlier rejection, project coordinates to a plane centered on the object with units of arcsec */
    … …  
    4348  Nmeasure = Naverage = 0;
    4449
    45   DVOAverageFlags averageBits =
    46     ID_STAR_FIT_AVE   |
    47     ID_STAR_FIT_PM    |
    48     ID_STAR_FIT_PAR   |
    49     ID_STAR_USE_AVE   |
    50     ID_STAR_USE_PM    |
    51     ID_STAR_USE_PAR   |
    52     ID_STAR_NO_ASTROM ;
     50  DVOAverageFlags astromBits =
     51    ID_OBJ_LARGE_PM        | // star with large proper motion
     52    ID_OBJ_RAW_AVE         | // simple weighted average position was used (no IRLS fitting)
     53    ID_OBJ_FIT_AVE         | // average position was fitted
     54    ID_OBJ_FIT_PM          | // proper motion model was fitted
     55    ID_OBJ_FIT_PAR         | // parallax model was fitted
     56    ID_OBJ_USE_AVE         | // average position used (not PM or PAR)
     57    ID_OBJ_USE_PM          | // proper motion used (not AVE or PAR)
     58    ID_OBJ_USE_PAR         | // parallax used (not AVE or PM)
     59    ID_OBJ_NO_MEAN_ASTROM  | // mean astrometry could not be measured
     60    ID_OBJ_STACK_FOR_MEAN  | // stack position used for mean astrometry
     61    ID_OBJ_MEAN_FOR_STACK  | // mean astrometry could not be measured
     62    ID_OBJ_BAD_PM;           // failure to measure proper-motion model
    5363
    5464  if (VERBOSE2 && ExcludeBogus && (fbogus == NULL)) {
    … …  
    7888
    7989    if (RESET) {
    80       subcatalog[0].average[Naverage].flags &= ~averageBits;
     90      subcatalog[0].average[Naverage].flags &= ~astromBits;
    8191    }
    8292
    8393    // if we want to do static correction to ICRF, keep the ICRF QSOs (regardless of other data quality)
    84     if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_ICRF_QSO)) {
     94    if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_OBJ_ICRF_QSO)) {
    8595      Nm = 0;
    8696      for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
    … …  
    115125    }
    116126
    117     // if we want to correct to the galay model, we need to use objects with distances determined by Green & Schlafly
     127    // if we want to correct to the galaxy model, we need to use objects with distances determined by Green & Schlafly
    118128    if (USE_GALAXY_MODEL) {
    119129      // if we use the galaxy model, we need to require the model pm to exist
    120       if (!catalog[0].average[i].Nstarpar) continue;
    121       if (!isfinite(catalog[0].average[i].uRgal)) continue;
    122       if (!isfinite(catalog[0].average[i].uDgal)) continue;
     130
     131      int gotStarpar = TRUE;
     132      gotStarpar = gotStarpar && (catalog[0].average[i].Nstarpar > 0);
     133      gotStarpar = gotStarpar && isfinite(catalog[0].average[i].uRgal);
     134      gotStarpar = gotStarpar && isfinite(catalog[0].average[i].uDgal);
     135      if (!gotStarpar) {
     136        // if we do not have starpar for some objects, generate an estimated starpar using d = 500pc
     137
     138        // NOTE: distance is in kiloparsec
     139        double distance = 0.5; // pow(10.0, 0.2*(starpar->DistMag + 5.0)) / 1000.0;
     140       
     141        double galLon, galLat;
     142        ApplyTransform (&galLon, &galLat, catalog[0].average[i].R, catalog[0].average[i].D, transform);
     143
     144        double Lrad = galLon * RAD_DEG;
     145        double Brad = galLat * RAD_DEG;
     146       
     147        double uL_gal, uB_gal;
     148        GalaxyMotionModel_radians(&uL_gal, &uB_gal, Lrad, Brad);
     149       
     150        double uL_sol, uB_sol;
     151        SolarMotionModel_radians(&uL_sol, &uB_sol, Lrad, Brad, distance);
     152       
     153        double uL = uL_gal + uL_sol;
     154        double uB = uB_gal + uB_sol;
     155       
     156        double uR, uD;
     157        TransformProperMotion (&uR, &uD, uL, uB, catalog[0].average[i].R, catalog[0].average[i].D, transform);
     158
     159        catalog[0].average[i].uRgal = uR;
     160        catalog[0].average[i].uDgal = uD;
     161      } else {
    123162     
    124       // filter stars which are not well-handled by the galactic motion model
    125       int mStarpar = catalog[0].average[i].starparOffset;
    126       myAssert (mStarpar < catalog[0].Nstarpar, "oops");
     163        // filter stars which are not well-handled by the galactic motion model
     164        int mStarpar = catalog[0].average[i].starparOffset;
     165        myAssert (mStarpar < catalog[0].Nstarpar, "oops");
    127166     
    128       StarPar *starpar = &catalog[0].starpar[mStarpar];
    129       if (starpar->DistMag  < MIN_DISTANCE_MOD) continue;
    130       if (starpar->DistMag  > MAX_DISTANCE_MOD) continue;
    131       if (starpar->dDistMag > MAX_DISTANCE_MOD_ERR) continue;
     167        StarPar *starpar = &catalog[0].starpar[mStarpar];
     168        if (starpar->DistMag  < MIN_DISTANCE_MOD) continue;
     169        if (starpar->DistMag  > MAX_DISTANCE_MOD) continue;
     170        if (starpar->dDistMag > MAX_DISTANCE_MOD_ERR) continue;
     171      }
    132172    }
    133173
    … …  
    264304  }
    265305
     306  if (transform) free (transform);
    266307  return (TRUE);
    267308}

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

    r39371 r39375  
    136136
    137137    double plx = 0.5*(drand48_cnt() + 0.0);
    138     if (FIT_MODE == FIT_PM_AND_PAR) plx = 0.0;
     138    if (FIT_MODE != FIT_PM_AND_PAR) plx = 0.0;
    139139
    140140    int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points
    141141    int Noutliers_max = Npoints * F_OUTLIERS;
    142142    int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points
     143
     144    if (Npoints + Noutliers == 0) continue;
    143145
    144146    // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
    … …  
    189191  catalog->average->catID = 1;
    190192
    191   catalog->average->R = Ro;
    192   catalog->average->D = Do;
    193   catalog->average->uR = uR;
    194   catalog->average->uD = uD;
    195   catalog->average->uRgal = uR;
    196   catalog->average->uDgal = uD;
     193  // catalog->average->R = Ro;
     194  // catalog->average->D = Do;
     195  // catalog->average->uR = uR;
     196  // catalog->average->uD = uD;
     197  // catalog->average->uRgal = uR;
     198  // catalog->average->uDgal = uD;
    197199
    198200  // int T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
    … …  
    250252    catalog->measureT[i].dt = 2.5*log10(30.0);
    251253  }
     254
     255  if (Ntotal > 0) {
     256    catalog->average->R = catalog->measureT[0].R;
     257    catalog->average->D = catalog->measureT[0].D;
     258    catalog->average->uRgal = uR;
     259    catalog->average->uDgal = uD;
     260  }
     261
    252262  catalog->average->Nmeasure = Ntotal;
    253263

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

    r38986 r39375  
    264264    good |= (testcat.average[0].ChiSqPM < MAX_CHISQ_PM);
    265265    good |= (testcat.average[0].Npos > MIN_NPOS_FRAC * testcat.average[0].Nmeasure);
    266     good &= ((testcat.average[0].flags & ID_STAR_FIT_PM) > 0);
     266    good &= ((testcat.average[0].flags & ID_OBJ_FIT_PM) > 0);
    267267
    268268    if (good) {

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

    r38986 r39375  
    1515    ID_MEAS_USED_CHIP;
    1616
    17   DVOAverageFlags averageBits =
    18     ID_STAR_FIT_AVE   |
    19     ID_STAR_FIT_PM    |
    20     ID_STAR_FIT_PAR   |
    21     ID_STAR_USE_AVE   |
    22     ID_STAR_USE_PM    |
    23     ID_STAR_USE_PAR   |
    24     ID_STAR_NO_ASTROM ;
     17  DVOAverageFlags astromBits =
     18    ID_OBJ_LARGE_PM        | // star with large proper motion
     19    ID_OBJ_RAW_AVE         | // simple weighted average position was used (no IRLS fitting)
     20    ID_OBJ_FIT_AVE         | // average position was fitted
     21    ID_OBJ_FIT_PM          | // proper motion model was fitted
     22    ID_OBJ_FIT_PAR         | // parallax model was fitted
     23    ID_OBJ_USE_AVE         | // average position used (not PM or PAR)
     24    ID_OBJ_USE_PM          | // proper motion used (not AVE or PAR)
     25    ID_OBJ_USE_PAR         | // parallax used (not AVE or PM)
     26    ID_OBJ_NO_MEAN_ASTROM  | // mean astrometry could not be measured
     27    ID_OBJ_STACK_FOR_MEAN  | // stack position used for mean astrometry
     28    ID_OBJ_MEAN_FOR_STACK  | // mean astrometry could not be measured
     29    ID_OBJ_BAD_PM;           // failure to measure proper-motion model
    2530
    2631  // XXX need to decide how to determine PARALLEL mode...
    … …  
    6974    if (RESET) {
    7075      for (j = 0; j < catalog.Naverage; j++) {
    71         catalog.average[j].flags &= ~averageBits;
     76        catalog.average[j].flags &= ~astromBits;
    7277        m = catalog.average[j].measureOffset;
    7378        for (k = 0; k < catalog.average[j].Nmeasure; k++) {

  • trunk/Ohana/src/relastro/test/dvo.photcodes

    r39363 r39375  
    3232
    3333  10233 GPC1.i.XY33          dep  24.611 -0.044 0.000     -     - 0.0000     0     3   0.000 1.000 0.000  0.000   0x0000 0x0000  0xe0440130 0x1003bc88
     34  11200 GPC1.i.SkyChip       dep  25.000  0.000 0.000     -     - 0.0000     0     3   0.000 1.000 0.000  0.000   0x0000 0x0000  0x0000 0x0000

  • trunk/Ohana/src/relastro/test/mana.sh

    r39371 r39375  
    2828  testnboot pm-boot
    2929  testnboot plx-boot
     30end
     31
     32macro teststk
     33  if ($0 != 3)
     34    echo "USAGE: testobj2 (Nstars) (foutliers)"
     35    break
     36  end
     37
     38  exec ../bin/fitstk.lin64 plx -Nstars $1 -foutliers $2 > test.stk.dat
     39  load.stack test.stk.dat
     40end
     41
     42macro testobj2
     43  if ($0 != 3)
     44    echo "USAGE: testobj2 (Nstars) (foutliers)"
     45    break
     46  end
     47
     48  exec ../bin/fitobj2.lin64 plx -Nstars $1 -foutliers $2 > test.plx.dat
     49  checkstats test.plx.dat
    3050end
    3151
    … …  
    191211end
    192212
     213macro load.stack
     214  if ($0 != 2)
     215    echo "USAGE: load.stack (file)"
     216    break
     217  end
     218
     219  data $1
     220  read Npts 1 Nout 2 Nstk 3 Ro 4 Do 5 uRo 6 uDo 7 Po 8 Rx 10 Dx 11 uRx 12 uDx 3 Px 14 Tx 16 dRo 18 dDo 19 duRo 20 duDo 21 dPo 22 Nfit 23 Rstk 25 Dstk 26 dRstk 27 dDstk 28 chisq 32
     221  # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
     222  # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22
     223
     224  foreach f uD uR P
     225    set d$f = $f\o - $f\x
     226  end
     227  set sdR = 3600*(Ro - Rstk)*dcos(Do)
     228  set sdD = 3600*(Do - Dstk)
     229
     230  set Rf = Rx - uRx*(Tx - $Tref)/3600.0/dcos(Dx)
     231  set Df = Dx - uDx*(Tx - $Tref)/3600.0
     232
     233  set dRf = 3600*(Ro - Rf)*dcos(Dx)
     234  set dDf = 3600*(Do - Df)
     235
     236  # echo "fit scatter"
     237  vstat -q dRf;  $value_dRf =  $SIGMA*1000
     238  vstat -q dDf;  $value_dDf =  $SIGMA*1000
     239  vstat -q duR;  $value_duR =  $SIGMA*1000
     240  vstat -q duD;  $value_duD =  $SIGMA*1000
     241  vstat -q dP ;  $value_dP  =  $SIGMA*1000
     242
     243  vstat -q sdR;  $value_sdR =  $SIGMA*1000
     244  vstat -q sdD;  $value_sdD =  $SIGMA*1000
     245
     246  # echo "fit errors"
     247  vstat -q dRo;  $value_dRo  = $MEAN*1000; $sigma_dRo  = $SIGMA*1000
     248  vstat -q dDo;  $value_dDo  = $MEAN*1000; $sigma_dDo  = $SIGMA*1000
     249  vstat -q duRo; $value_duRo = $MEAN*1000; $sigma_duRo = $SIGMA*1000
     250  vstat -q duDo; $value_duDo = $MEAN*1000; $sigma_duDo = $SIGMA*1000
     251  vstat -q dPo ; $value_dPo  = $MEAN*1000; $sigma_dPo  = $SIGMA*1000
     252
     253  # echo "Nfit"
     254  vstat -q Nfit; $mean_Nfit = $MEAN
     255
     256  # echo chisq
     257  vstat -q chisq; $mean_chisq = $MEAN
     258
     259  fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
     260 
     261  delete dRbin oRbin
     262  for i 0 100
     263    subset tmp = dRf if (Npts == $i)
     264    vstat tmp -q
     265    concat $SIGMA dRbin
     266    concat $MEAN oRbin
     267  end
     268
     269  create n 0 dRbin[]
     270  set dRmod = 0.010 / sqrt (n)
     271  lim n -0.001 0.015; clear; box; plot n dRbin
     272  plot -c red -pt 7 n dRmod
     273end
     274
    193275# determine the parallax factor (time) (ra) (dec) (parfR) (parfD) are vectors,
    194276macro getParFactor_vectors
Note: See TracChangeset for help on using the changeset viewer.