IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33422


Ignore:
Timestamp:
Mar 6, 2012, 7:07:53 AM (14 years ago)
Author:
eugene
Message:

mods from Measure to MeasureTiny

Location:
branches/eam_branches/ipp-20111122/Ohana/src/relastro
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/parallel.txt

    r33419 r33422  
    1515 * relastro_objects, high_speed_catalogs, and standard relastro all call UpdateObjects, one way or another.
    1616
    17  The big question is if we can use
     17 The big question is if / how we can use MeasureTiny (or if we need a separate astrometry version)
    1818
    19  UpdateObjects:
    20  *       if (catalog[i].average[j].Nmeasure == 0) {
     19  * average.d fields we are NOT using in UpdateObjects (28 / 120 bytes)
     20  FIELD Nmissing,       NMISSING,    unsigned short,  number of missings
     21  FIELD Nextend,        NEXTEND,     unsigned short,  number of extended measurements
     22  FIELD missingOffset,  OFF_MISSING, int,             offset to first missing obs
     23  FIELD extendOffset,   OFF_EXTEND,  int,             offset to first extended measurement
     24  FIELD photFlagsUpper, PHOTFLAGS_U, uint32_t,        upper bit of 2 bit summary of per-measure photflags
     25  FIELD photFlagsLower, PHOTFLAGS_L, uint32_t,        lower bit of 2 bit summary of per-measure photflags
     26  FIELD extID,          EXT_ID,      uint64_t,        external ID for object (eg PSPS objID)
     27
     28  * measure.d fields in use in UpdateObjects
     29* measure[m].t
     30* measure[m].dbFlags
     31* measure[m].dt
     32* measure[0].dR
     33* measure[0].dD
     34* measure[0].M
     35* measure[0].dM
     36* measure[0].photcode
     37* measure[0].photFlags
     38* measure[m].Xccd;
     39* measure[m].Yccd;
     40* measure[m].averef;
     41* measure[0].imageID;
     42* measure[0].airmass
     43* measure[0].Mcal (only used in PhotRelTiny)
     44+ measure[0].dXccd
     45+ measure[0].dYccd
     46+ measure[m].dRsys
     47
     48  float          dR;
     49  float          dD;
     50  float          M;
     51  float          Mcal;
     52  float          dM;
     53  float          airmass;
     54  float          Xccd;
     55  float          Yccd;
     56  float          dt;
     57  int            t;
     58  unsigned int   averef;
     59  unsigned int   imageID;
     60  unsigned int   dbFlags;
     61  unsigned int   photFlags;
     62- int            catID; // unsigned int?
     63  unsigned short photcode;
     64
     65maybe
     66FIELD Map,            M_APER,       float,          aperture mag,                   mag
     67FIELD psfQual,        PSF_QF,       float,          psf coverage/quality factor
     68
     69not used
     70FIELD dMcal,          MAG_CAL_ERR,  float,          systematic calibration error,   mag
     71FIELD az,             AZ,           float,          telescope azimuth
     72FIELD Sky,            SKY_FLUX,     float,          local estimate of sky flux,     counts/sec
     73FIELD dSky,           SKY_FLUX_ERR, float,          local estimate of sky flux,     counts/sec
     74FIELD t_msec,         TIME_MSEC,    unsigned short, time fraction of second,        milliseconds
     75FIELD detID,          DET_ID,       unsigned int,   detection ID
     76FIELD objID,          OBJ_ID,       unsigned int,   unique ID for object in table
     77FIELD catID,          CAT_ID,       unsigned int,   unique ID for table in which object was first realized
     78FIELD psfChisq,       PSF_CHISQ,    float,          psf fit chisq
     79FIELD psfNdof,        PSF_NDOF,     int,            psf degrees of freedom
     80FIELD psfNpix,        PSF_NPIX,     int,            psf number of pixels
     81FIELD crNsigma,       CR_NSIGMA,    float,          Nsigma deviation towards CR
     82FIELD extNsigma,      EXT_NSIGMA,   float,          Nsigma deviation towards EXT
     83FIELD FWx,            FWHM_MAJOR,   short,          object fwhm major axis,         1/100 of pixels
     84FIELD FWy,            FWHM_MINOR,   short,          object fwhm minor axis,         1/100 of pixels
     85FIELD theta,          PSF_THETA,    short,          angle wrt ccd X dir,            (0xffff/360) deg
     86FIELD Mxx,            MXX,          short,          second moments in pixel coords, 1/100 of pixels^2
     87FIELD Mxy,            MXY,          short,          second moments in pixel coords, 1/100 of pixels^2
     88FIELD Myy,            MYY,          short,          second moments in pixel coords, 1/100 of pixels^2
     89FIELD posangle,       POSANGLE,     short,          position angle sky to chip,     (0xffff/360) deg
     90FIELD pltscale,       PLTSCALE,     float,          plate scale,                    arcsec/pixel
     91

  • branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FixProblemImages.c

    r32346 r33422  
    55
    66int FixProblemImages (SkyList *skylist) {
     7
     8    fprintf (stderr, "this function is currently not used : check on change from measure -> measureT in load_catalogs\n");
     9    abort();
    710
    811  int Nbad;

  • branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c

    r32695 r33422  
    181181
    182182  off_t idx, ID;
    183   Measure *measure;
     183  MeasureTiny *measure;
    184184  int i, found;
    185185
    186   measure = &catalog[cat].measure[meas];
     186  measure = &catalog[cat].measureT[meas];
    187187
    188188  ID = measure[0].imageID;
    … …  
    241241
    242242  off_t i;
    243   Measure *measure;
    244 
    245   measure = &catalog[cat].measure[meas];
     243  MeasureTiny *measure;
     244
     245  measure = &catalog[cat].measureT[meas];
    246246
    247247  /* find the image that supplied this measurement */
    … …  
    333333}
    334334
    335 void dump_measures(Average *average, Measure *measure) {
     335void dump_measures(Average *average, MeasureTiny *measure) {
    336336
    337337  off_t j, off;
    … …  
    406406    c = clist[im][i];
    407407
    408     X = catalog[c].measure[m].Xccd;
    409     Y = catalog[c].measure[m].Yccd;
    410     n = catalog[c].measure[m].averef;
     408    MeasureTiny *measure = &catalog[c].measureT[m];
     409
     410    X = measure[0].Xccd;
     411    Y = measure[0].Yccd;
     412    n = measure[0].averef;
    411413
    412414    if (moscoords == NULL) {
    … …  
    429431    if (fabs(dR) > 3.0*ADDSTAR_RADIUS) {
    430432      fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
    431       dump_measures (&catalog[c].average[n], catalog[c].measure);
     433      dump_measures (&catalog[c].average[n], catalog[c].measureT);
    432434      // abort ();
    433435    }
    434436    if (fabs(dD) > 3.0*ADDSTAR_RADIUS) {
    435437      fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
    436       dump_measures (&catalog[c].average[n], catalog[c].measure);
     438      dump_measures (&catalog[c].average[n], catalog[c].measureT);
    437439      // abort ();
    438440    }
    439441
    440442    // complain if the new location is far from the old location
    441     if (fabs(catalog[c].measure[m].dR - dR) > DPOS_MAX_ASEC) {
     443    if (fabs(catalog[c].measureT[m].dR - dR) > DPOS_MAX_ASEC) {
    442444      fprintf (stderr, "measurement is far from original location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
    443       dump_measures (&catalog[c].average[n], catalog[c].measure);
     445      dump_measures (&catalog[c].average[n], catalog[c].measureT);
    444446      // abort();
    445447    }
    446     if (fabs(catalog[c].measure[m].dD - dD) > DPOS_MAX_ASEC) {
     448    if (fabs(measure[0].dD - dD) > DPOS_MAX_ASEC) {
    447449      fprintf (stderr, "measurement is far from original location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
    448450      dump_measures (&catalog[c].average[n], catalog[c].measure);
    … …  
    450452    }
    451453
    452     dPos += SQ(catalog[c].measure[m].dR - dR) + SQ(catalog[c].measure[m].dD - dD);
     454    dPos += SQ(measure[0].dR - dR) + SQ(measure[0].dD - dD);
    453455    nPos ++;
    454456
    455     catalog[c].measure[m].dR = dR;
    456     catalog[c].measure[m].dD = dD;
    457    
    458     if (catalog[c].measure[m].dR > +180.0*3600.0) {
     457    measure[0].dR = dR;
     458    measure[0].dD = dD;
     459   
     460    if (measure[0].dR > +180.0*3600.0) {
    459461      // average on high end of boundary, move star up
    460462      R += 360.0;
    461       catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    462     }
    463     if (catalog[c].measure[m].dR < -180.0*3600.0) {
     463      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
     464    }
     465    if (measure[0].dR < -180.0*3600.0) {
    464466      // average on low end of boundary, move star down
    465467      R -= 360.0;
    466       catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     468      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
    467469    }
    468470
    469471    // set the systematic error for this image:
    470     catalog[c].measure[m].dRsys = ToShortPixels(dPosSys);
     472    measure[0].dRsys = ToShortPixels(dPosSys);
    471473  }
    472474
    … …  
    487489  Mosaic *mosaic;
    488490  Coords *moscoords, *imcoords, *oldcoords;
     491
     492  fprintf (stderr, "fix resetImageRaw wrt MeasureTiny\n");
     493  abort();
    489494
    490495  // check if this image is bad and should be skipped
    … …  
    511516    c = clist[im][i];
    512517
    513     X = catalog[c].measure[m].Xccd;
    514     Y = catalog[c].measure[m].Yccd;
    515     n = catalog[c].measure[m].averef;
     518    // XXX unclear if this should use Measure or MeasureTiny; it is only called by FixProblemImages, which has not been updated
     519    Measure *measure = &catalog[c].measure[m];
     520
     521    X = measure[0].Xccd;
     522    Y = measure[0].Yccd;
     523    n = measure[0].averef;
    516524
    517525    dR = dD = 0.0;
    … …  
    527535    }
    528536
    529     catalog[c].measure[m].dR = dR;
    530     catalog[c].measure[m].dD = dD;
    531 
    532     if (catalog[c].measure[m].dR > +180.0*3600.0) {
     537    measure[0].dR = dR;
     538    measure[0].dD = dD;
     539
     540    if (measure[0].dR > +180.0*3600.0) {
    533541      // average on high end of boundary, move star up
    534542      R += 360.0;
    535       catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    536     }
    537     if (catalog[c].measure[m].dR < -180.0*3600.0) {
     543      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
     544    }
     545    if (measure[0].dR < -180.0*3600.0) {
    538546      // average on low end of boundary, move star down
    539547      R -= 360.0;
    540       catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     548      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
    541549    }
    542550  }
    … …  
    573581    c = clist[im][i];
    574582
     583    MeasureTiny *measure = &catalog[c].measure[m];
     584
    575585    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    576     raw[i].X = catalog[c].measure[m].Xccd;
    577     raw[i].Y = catalog[c].measure[m].Yccd;
    578 
    579     raw[i].Mag  = catalog[c].measure[m].M;
    580     raw[i].dMag = catalog[c].measure[m].dM;
    581     raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
    582 
    583     n = catalog[c].measure[m].averef;
     586    raw[i].X = measure[0].Xccd;
     587    raw[i].Y = measure[0].Yccd;
     588
     589    raw[i].Mag  = measure[0].M;
     590    raw[i].dMag = measure[0].dM;
     591    raw[i].dPos = GetAstromError (&measure[0], ERROR_MODE_POS);
     592
     593    n = measure[0].averef;
    584594
    585595    // an object with only one detection provides no information about the image calibration
    … …  
    589599      raw[i].mask |= 0x0001;
    590600    }
    591     if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
     601    if (!finite(measure[0].dR) || !finite(measure[0].dD)) {
    592602      raw[i].mask |= 0x0002;
    593603    }
    … …  
    646656    n = catalog[c].measure[m].averef;
    647657
     658    MeasureTiny *measure = &catalog[c].measure[m];
     659
    648660    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    649661    ref[i].R = catalog[c].average[n].R;
    650662    ref[i].D = catalog[c].average[n].D;
    651663
    652     ref[i].Mag  = catalog[c].measure[m].M;
    653     ref[i].dMag = catalog[c].measure[m].dM;
    654     ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
     664    ref[i].Mag  = measure[0].M;
     665    ref[i].dMag = measure[0].dM;
     666    ref[i].dPos = GetAstromError (&measure[0], ERROR_MODE_POS);
    655667
    656668    ref[i].mask = FALSE;
    … …  
    683695
    684696void FlagOutliers (Catalog *catalog) {
    685   // XXX FlagOutliers is just using FlagOutliers2D
     697
     698  // XXX FlagOutliers is now just using FlagOutliers2D
    686699  FlagOutliers2D(catalog);
    687700  return;
    … …  
    716729  for (j = 0; j < catalog[0].Naverage; j++) {
    717730   
     731    // pointer to this set of measurements
     732    m = catalog[0].average[j].measureOffset;
     733    MeasureTiny *measure = &catalog[0].measureT[m];
     734
    718735    /* accumulate list of valid measurements */
    719     m = catalog[0].average[j].measureOffset;
    720736    N = 0;
    721     for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
    722      
     737    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    723738      // skip measurements based on user selected criteria
    724       if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
    725       R[N] = catalog[0].measure[m].dR;
    726       D[N] = catalog[0].measure[m].dD;
    727       dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
    728       dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
     739      if (!MeasFilterTest(&measure[k], FALSE)) continue;
     740      R[N] = measure[k].dR;
     741      D[N] = measure[k].dD;
     742      dR[N] = GetAstromError (&measure[k], ERROR_MODE_RA);
     743      dD[N] = GetAstromError (&measure[k], ERROR_MODE_DEC);
    729744      if (isnan(R[N]) || isnan(D[N])) continue;
    730745      N++;
    … …  
    741756   
    742757    /* compare per-object distance to this standard deviation, and flag outliers*/
    743     m = catalog[0].average[j].measureOffset;
    744758    N = 0;
    745     for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
     759    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    746760      //reset flag on each invocation
    747       catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     761      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    748762
    749763      // skip measurements based on user selected criteria
    750       if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
     764      if (!MeasFilterTest(&measure[k], FALSE)) continue;
    751765     
    752       x = catalog[0].measure[m].dR - statsR.median;
    753       y = catalog[0].measure[m].dD - statsD.median;
     766      x = measure[k].dR - statsR.median;
     767      y = measure[k].dD - statsD.median;
    754768      theta = atan2(y,x);
    755769      if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) +
    756770                         SQR(statsD.sigma * Ns * sin(theta)))) {   
    757         catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
     771        measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
    758772        Ndel++;
    759773      }
    … …  
    809823  for (j = 0; j < catalog[0].Naverage; j++) {
    810824   
     825    // pointer to this set of measurements
     826    m = catalog[0].average[j].measureOffset;
     827    MeasureTiny *measure = &catalog[0].measureT[m];
     828
    811829    /* accumulate list of valid measurements */
    812     m = catalog[0].average[j].measureOffset;
    813830    N = 0;
    814    
    815     for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
    816      
    817       //reset flag on each invocation
    818       catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     831    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
     832
     833      // reset flag on each invocation
     834      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    819835     
    820836      // skip measurements based on user selected criteria
    821       if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
    822       R[N] = catalog[0].measure[m].dR;
    823       D[N] = catalog[0].measure[m].dD;
    824       dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
    825       dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
     837      if (!MeasFilterTest(&measure[k], FALSE)) continue;
     838      R[N] = measure[k].dR;
     839      D[N] = measure[k].dD;
     840      dR[N] = GetAstromError(&measure[k], ERROR_MODE_RA);
     841      dD[N] = GetAstromError(&measure[k], ERROR_MODE_DEC);
    826842      if (isnan(R[N]) || isnan(D[N])) continue;
    827843      N++;
    … …  
    837853   
    838854    /* calculate deviations of all points*/
    839     m = catalog[0].average[j].measureOffset;
    840855    N = 0;
    841     for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
    842       //skip bad measurements
    843       if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; 
    844       x = catalog[0].measure[m].dR - statsR.median;
    845       y = catalog[0].measure[m].dD - statsD.median;
     856    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
     857      // skip bad measurements
     858      if (!MeasFilterTest(&measure[k], FALSE)) continue; 
     859      x = measure[k].dR - statsR.median;
     860      y = measure[k].dD - statsD.median;
    846861      theta = atan2(y,x);
    847862      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) +
    848863                             SQR(statsD.sigma * Ns * sin(theta)));     
    849       index[N] = m;
     864      index[N] = k;
    850865      N++;
    851866    }
    … …  
    856871
    857872    // recalculate image center, sigma based on closest 50% of points
    858     for(k = 0;  k < N; k++) {
     873    for (k = 0;  k < N; k++) {
    859874      off_t ind = (off_t) index[k];
    860       R[k] = catalog[0].measure[ind].dR;
    861       D[k] = catalog[0].measure[ind].dD;
    862       dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA);
    863       dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC);
     875      R[k] = measure[ind].dR;
     876      D[k] = measure[ind].dD;
     877      dR[k] = GetAstromError(&measure[ind], ERROR_MODE_RA);
     878      dD[k] = GetAstromError(&measure[ind], ERROR_MODE_DEC);
    864879    }
    865880    liststats (R, dR, N, &statsR);
    … …  
    869884   
    870885    // use these new statistics to flag outliers
    871     m = catalog[0].average[j].measureOffset;
    872886    N = 0;
    873     for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
     887    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
    874888      //skip bad measurements
    875       if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; 
    876       x = catalog[0].measure[m].dR - statsR.median;
    877       y = catalog[0].measure[m].dD - statsD.median;
     889      if (!MeasFilterTest(&measure[k], FALSE)) continue; 
     890      x = measure[k].dR - statsR.median;
     891      y = measure[k].dD - statsD.median;
    878892      theta = atan2(y,x);
    879893      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) +
    880894                             SQR(statsD.sigma * Ns * sin(theta)));     
    881895      if ((d2[N]) > 1) {
    882         catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
     896        measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
    883897        Ndel ++;
    884898      }
    885899      N++;
    886900      Nave++;
    887     }  //done rejecting outliers
     901    }  // done rejecting outliers
    888902
    889903    // examine results
    890904    // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
    891905   
    892   } //done looping over objects
     906  } // done looping over objects
    893907 
    894908  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
    … …  
    905919/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */
    906920// we only optionally apply the sigma limit: for object averages, this should not be used (should it?)
    907 int MeasFilterTest(Measure *measure, int applySigmaLim) {
     921int MeasFilterTest(MeasureTiny *measure, int applySigmaLim) {
    908922  int found, k;
    909923  long mask;
    … …  
    958972  /* select measurements by mag limit */
    959973  if (ImagSelect) {
    960     mag = PhotInst (measure);
     974    mag = PhotInstTiny (measure);
    961975    if (mag < ImagMin || mag > ImagMax) return FALSE;
    962976  }

  • branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/StarMaps.c

    r32346 r33422  
    6868  for (i = 0; i < catalog[0].Nmeasure; i++) {
    6969   
    70     N = getImageByID(catalog[0].measure[i].imageID);
     70    MeasureTiny *measure = &catalog[0].measureT[i]
     71
     72    N = getImageByID(measure[0].imageID);
    7173    if (N < 0) continue;
    7274
    73     xbin = catalog[0].measure[i].Xccd / starmap[N].Nx;
    74     ybin = catalog[0].measure[i].Yccd / starmap[N].Ny;
     75    xbin = measure[0].Xccd / starmap[N].Nx;
     76    ybin = measure[0].Yccd / starmap[N].Ny;
    7577
    7678    xbin = MAX(0, MIN(NX_MAP-1, xbin));

  • branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjectOffsets.c

    r30616 r33422  
    4040    UpdateMeasures (&catalog, 1);
    4141
     42  XXX: note that this version of UpdateObjects is / should be expecting catalog.measure, not catalog.measureT
    4243    UpdateObjects (&catalog, 1);
    4344

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

    r32854 r33422  
    8686      XVERB = FALSE;
    8787
    88       // skip objects which are known to be problematic
    89       // XXX include this code or not?
    90 # if (0)
    91       if (catalog[i].average[j].code & STAR_BAD) {
    92         Nskip ++;
    93         continue; 
    94       }
    95 # endif
    96 
    9788      if (catalog[i].average[j].Nmeasure == 0) {
    9889          continue;
    … …  
    10192      N = 0;
    10293      m = catalog[i].average[j].measureOffset;
    103       Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25);
     94      MeasureTiny *measure = &catalog[i].measureT[m];
     95
     96      Tmin = Tmax = (measure[0].t - T2000) / (86400*365.25);
    10497      mode = FIT_MODE;
    10598
    10699      // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
    107100      Tmean = 0;
    108       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    109 
    110         //does the measurement pass the supplied filtering constraints?
    111         if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) {
    112           catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
     101      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
     102
     103        // does the measurement pass the supplied filtering constraints?
     104        if (!MeasFilterTest(&measure[k], FALSE)) {
     105          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    113106          continue;
    114107        }
    115108
    116109        //outlier rejection
    117         if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
    118           catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
     110        if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
     111          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    119112          continue;
    120113        }
    121114
    122         // exclude measurements by previous outlier detection
    123         // XXX include this code or not?
    124 # if (0)
    125         if (catalog[i].measure[m].dbFlags & MEAS_BAD) {
    126           catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
    127           continue;
    128         }
    129 # endif
    130 
    131         catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
    132 
    133         R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    134         D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     115        measure[k].dbFlags |= ID_MEAS_USED_OBJ;
     116
     117        R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     118        D[N] = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    135119
    136120        // XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects
    137121        // to get Time in years relative to UNIX Tzero (1970/01/01)
    138         T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
     122        T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
    139123
    140124        Tmin = MIN(Tmin, T[N]);
    … …  
    143127
    144128        // dX, dY : error in arcsec --
    145         dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
    146         dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
     129        dX[N] = GetAstromError (&measure[k], ERROR_MODE_RA);
     130        dY[N] = GetAstromError (&measure[k], ERROR_MODE_DEC);
    147131
    148132        // add systematic error in quadrature, if desired
    149133        // only do this after the fit has converged (or you will never improve the poor images)
    150134        // if (INCLUDE_SYS_ERR) {
    151         // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys);
     135        // float dRsys = FromShortPixels(measure[k].dRsys);
    152136        // dX[N] = hypot(dX[N], dRsys);
    153137        // dY[N] = hypot(dY[N], dRsys);
    … …  
    157141        // dY[N] = 0.1;
    158142
    159         dT[N] = catalog[i].measure[m].dt;
     143        dT[N] = measure[k].dt;
    160144
    161145        // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
    … …  
    204188      XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
    205189      XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
    206       // XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
    207190
    208191      // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
    … …  
    311294      // the measure fields must be updated before the average fields
    312295      m = catalog[i].average[j].measureOffset;
     296      MeasureTiny *measure = &catalog[i].measureT[m];
    313297      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    314         setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    315         setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     298        setMeanR (fit.Ro, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
     299        setMeanD (fit.Do, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    316300      }     
    317301

  • branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c

    r31664 r33422  
    1212
    1313  /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
     14  memset(subcatalog, 0, sizeof(Catalog));
    1415  NAVERAGE = 50;
    1516  NMEASURE = 1000;
    16   ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
    17   ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
    18   ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
     17  ALLOCATE (subcatalog[0].average,  Average,    NAVERAGE);
     18  ALLOCATE (subcatalog[0].secfilt,  SecFilt,    NAVERAGE*Nsecfilt);
     19  ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
    1920  Nmeasure = Naverage = 0;
    2021
    … …  
    7071      // allowed.
    7172
    72       // CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
    73 
    74       subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset];
    75       subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
    76       subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
    77       subcatalog[0].measure[Nmeasure].averef   = Naverage;
     73      CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
     74      // subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset];
     75      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
     76      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
     77      subcatalog[0].measureT[Nmeasure].averef   = Naverage;
    7878      if (RESET) {
    79         subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    80         subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
     79        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     80        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
    8181      }
    8282
    … …  
    8585      if (Nmeasure == NMEASURE) {
    8686        NMEASURE += 1000;
    87         REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
     87        REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
    8888      }
    8989    }
    … …  
    9696    }
    9797  }
    98   REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
    99   REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
    100   REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt);
     98  REALLOCATE (subcatalog[0].average,  Average,    MAX (Naverage, 1));
     99  REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1));
     100  REALLOCATE (subcatalog[0].secfilt,  SecFilt,    MAX (Naverage, 1)*Nsecfilt);
    101101  subcatalog[0].Naverage = Naverage;
    102102  subcatalog[0].Nmeasure = Nmeasure;
    … …  
    168168
    169169  // allocate the output data
    170   ALLOCATE (tmpcatalog.average, Average, Nmax);
    171   ALLOCATE (tmpcatalog.measure, Measure, NMEASURE);
    172   ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt);
     170  ALLOCATE (tmpcatalog.average,  Average,    Nmax);
     171  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
     172  ALLOCATE (tmpcatalog.secfilt,  SecFilt,    Nmax * Nsecfilt);
    173173
    174174  off_t Nmeasure = 0;
    … …  
    181181    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
    182182      off_t offset = subcatalog[0].average[ave].measureOffset + j;
    183       tmpcatalog.measure[Nmeasure] = subcatalog[0].measure[offset];
    184       tmpcatalog.measure[Nmeasure].averef = i;
     183      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
     184      tmpcatalog.measureT[Nmeasure].averef = i;
    185185      Nmeasure ++;
     186    }
     187    for (j = 0; j < Nsecfilt; j++) {
     188      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
    186189    }
    187190  }
    … …  
    193196
    194197  free (subcatalog[0].average);
    195   free (subcatalog[0].measure);
     198  free (subcatalog[0].measureT);
    196199  free (subcatalog[0].secfilt);
    197200
    198201  subcatalog[0].average = tmpcatalog.average;
    199   subcatalog[0].measure = tmpcatalog.measure;
     202  subcatalog[0].measureT = tmpcatalog.measureT;
    200203  subcatalog[0].secfilt = tmpcatalog.secfilt;
    201204  subcatalog[0].Naverage = Nmax;
Note: See TracChangeset for help on using the changeset viewer.