IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38297


Ignore:
Timestamp:
May 21, 2015, 9:23:17 AM (11 years ago)
Author:
eugene
Message:

relphot rework compiles

Location:
branches/eam_branches/ipp-20150419/Ohana/src/relphot
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h

    r38293 r38297  
    136136  float *psfQfPerfMax;        // max psfQfperf value for this secfilt
    137137
     138  int *Nmeas;                 // count of PS1 exposure (chip) measurements for this secfilt
    138139  int *Next;                  // count of PS1 exposure (chip) measurements for this secfilt
    139140  int *NexpPS1;               // count of PS1 exposure (chip) measurements for this secfilt
    … …  
    141142
    142143  float *minUbercalDist;
    143 
    144   int NextPS1;                // count of PS1 exposure (chip) measurements for this secfilt
    145   int NpsfPS1;                // count of PS1 exposure (chip) measurements for this secfilt
    146144
    147145  StatType psfstats;
    … …  
    396394void          liststats_setmode   PROTO((StatType *stats, char *strmode));
    397395int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
     396int           liststats_init      PROTO((StatType *stats));
    398397Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile));
    399398Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile));
    … …  
    433432int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    434433void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage));
    435 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr));
     434int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    436435int           setMave             PROTO((Catalog *catalog, int Ncatalog));
    437436void          set_ZP              PROTO((double ZERO));
    … …  
    495494
    496495int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure);
    497 int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
    498 int setMrelAverageExposure (off_t measureOffset, int cat, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT);
    499 int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure);
    500 int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure);
     496int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
     497int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results);
     498int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results);
     499int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results);
    501500
    502501int setGlobalObjStats (Average *average, Measure *measure);
    … …  
    562561int whichGPC1filter (int photcode);
    563562int is2MASS (int photcode);
     563int isTYCHO (int photcode);
     564
     565int magStatsByRanking (StatDataSet *dataset, StatType *stats);
    564566
    565567SynthZeroPoints *SynthZeroPointsLoad (char *filename);
    … …  
    570572int relphot_synthphot_catalog (Catalog *catalog, SynthZeroPoints *zpts);
    571573int relphot_synthphot_average (Average *average, SecFilt *secfilt, Measure *measure, SynthZeroPoints *zpts);
     574
     575void setMeasureRank (Catalog *catalog);
     576int getImageFlags (off_t meas, int cat);
     577int getMosaicFlags (off_t meas, int cat);
     578

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c

    r38293 r38297  
    111111    results->psfQfPerfMax[i]  = NAN;
    112112
     113    results->Nmeas[i] = 0;
    113114    results->Next[i] = 0;
    114115    results->NexpPS1[i] = 0;
    … …  
    121122    results->kronData[i].Nlist = 0;
    122123  }
    123   results->NextPS1 = 0;
    124   results->NpsfPS1 = 0;
    125124}
    126125
    … …  
    149148    ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt);
    150149
     150    ALLOCATE (results->Nmeas,          int, results->Nsecfilt);
    151151    ALLOCATE (results->Next,           int, results->Nsecfilt);
    152152    ALLOCATE (results->NexpPS1,        int, results->Nsecfilt);
    … …  
    156156    ALLOCATE (results->minUbercalDist, float, results->Nsecfilt);
    157157
    158     SetMrelInfoResetObject (results, results->Nsecfilt);
     158    SetMrelInfoResetObject (results);
    159159
    160160    for (i = 0; i < NsynthNames; i++) {
    … …  
    259259        break;
    260260      case 1:
    261         setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
     261        setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
    262262        break;
    263263      default:
    … …  
    274274}
    275275
    276 int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {
     276int setMrelOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    277277
    278278  int i;
    … …  
    293293      // break;
    294294      case 1:
    295         setMrel_catalog_alt (catalog, i, pass, flatcorr, &results, Nsecfilt);
     295        setMrel_catalog_alt (catalog, i, TRUE, flatcorr, &results, Nsecfilt);
    296296        break;
    297297      default:
    … …  
    400400        break;
    401401      case 1:
    402         setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
     402        setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
    403403        break;
    404404      default:

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/extra.c

    r38293 r38297  
    7777}
    7878
     79int isTYCHO (int photcode) {
     80
     81  if ((photcode == 2020) || (photcode == 2021)) return TRUE;
     82  return FALSE;
     83}
     84

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c

    r38293 r38297  
    11# include "relphot.h"
    22# define UBERCAL_WEIGHT 100.0
     3void sort_StatDataSet (StatDataSet *dataset);
    34
    45# define SKIP_THIS_MEAS(REASON) {                       \
    … …  
    1617# define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > 0.0) && (MAG < 30.0))
    1718
     19# if (0)
    1820static float MagToFlux (float Mag) {
    1921  float Flux = pow(10.0, -0.4*(Mag));
    2022  return (Flux);
    2123}
     24# endif
    2225
    2326int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) {
    … …  
    3841}
    3942
    40 int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
     43int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
    4144
    4245  off_t j;
    … …  
    4851  SetMrelInfoReset (results); // reset the countesrs
    4952
    50   int isSetMrelFinal = (pass >= 0);
    51 
    5253  for (j = 0; j < catalog[Nc].Naverage; j++) {
    53     // we are guaranteed to have averageT, but not average
    54     off_t m = catalog[Nc].averageT[j].measureOffset;
    55     Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL;
    56     Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL;
    57     SecFilt *secfilt = catalog[Nc].secfilt ? &catalog[Nc].secfilt[j*Nsecfilt] : NULL;
    58     setMrelAverageExposure (m, Nc, isSetMrelFinal, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);
     54
     55    setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, flatcorr, results);
    5956
    6057    // only apply Stack operation on setMrelFinal in first pass
    61     if (isSetMrelFinal && (pass == 0) && !IS_DIFF_DB) {
    62       setMrelAverageStack (m, Nc, flatcorr, results, average, secfilt, &catalog[Nc].measure[m], catalog[Nc].Nmeasure);
    63       setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
     58    if (isSetMrelFinal && !IS_DIFF_DB) {
     59      setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results);
    6460    }
    6561
    6662    // only measure force-warp mean values if issetMrelFinal (make it optional?)
    6763    if (isSetMrelFinal) {
    68       setMrelAverageForcedWarp (m, Nc, pass, results, average, secfilt, measure, &catalog[Nc].foundWarp_t[Nsecfilt*j], catalog[Nc].Nmeasure);
    69       setGlobalObjStats (average, measure);
     64      setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results);
    7065    }
    7166  }
    … …  
    9186// 4) some reference photcode of some kind can be specified as fixed and have a high weight
    9287
     88static SynthZeroPoints *synthzpts = NULL;
     89
    9390// set mean of chip measurements (selected by photcode range for now):
    94 int setMrelAverageExposure (off_t measureOffset, int cat, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT) {
    95 
     91int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
     92
     93  int Nsec;
    9694  off_t k;
     95
     96  // we are guaranteed to have averageT, but not average
     97  AverageTiny *averageT = &catalog[cat].averageT[ave];
     98
     99  off_t measureOffset = averageT->measureOffset;
     100  MeasureTiny *measureT = &catalog[cat].measureT[measureOffset];
     101
     102  // we are NOT guaranteed to have average, measure, secfilt
     103  Average *average  = catalog[cat].average     ? &catalog[cat].average[ave]               : NULL;
     104  Measure *measure  = catalog[cat].measure     ? &catalog[cat].measure[measureOffset]     : NULL;
     105  SecFilt *secfilt  = catalog[cat].secfilt     ? &catalog[cat].secfilt[ave*Nsecfilt]      : NULL;
     106  char *measureRank = catalog[cat].measureRank ? &catalog[cat].measureRank[measureOffset] : NULL;
    97107
    98108  // we are measuring means for 3 types of mags: psf, ap, kron.
    … …  
    110120
    111121  SetMrelInfoResetObject (results); // reset the per-object arrays
     122
     123  int NextPS1 = 0;
     124  int NpsfPS1 = 0;
     125  int GoodPS1 = FALSE;
     126  int Good2MASS = FALSE;
     127  int Galaxy2MASS = FALSE;
     128  int haveTYCHO = FALSE;
    112129
    113130  // assign measurements to the photcode lists
    … …  
    124141    if (isGPC1warp(measureT[k].photcode)) continue;
    125142   
    126     if (isTycho(measureT[k].photcode)) { haveTYCHO = TRUE; }
     143    if (isTYCHO(measureT[k].photcode)) { haveTYCHO = TRUE; }
    127144
    128145    if (is2MASS(measureT[k].photcode)) {
    … …  
    138155    int Nsec = GetPhotcodeNsec (code->equiv);
    139156    if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode
     157
     158    results->Nmeas[Nsec] ++;
    140159
    141160    // force the use of SYN even if we have PS1 mags?
    … …  
    148167
    149168    // assign the max psfQF values for each Nsec value
    150     if (measureT[k].psfQF     > results->psfQfMax[Nsec])     results->psfQfMax[Nsec]     = measureT[k].psfQF;
    151     if (measureT[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measureT[k].psfQFperf;
    152    
     169    if (isSetMrelFinal) {
     170      if (measure[k].psfQF     > results->psfQfMax[Nsec])     results->psfQfMax[Nsec]     = measure[k].psfQF;
     171      if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf;
     172    }
     173
    153174    // are we a PS1 exposure photcode? (hard-wired photcodes)
    154175    if (isGPC1chip(measureT[k].photcode)) { results->NexpPS1[Nsec] ++; results->havePS1[Nsec] = TRUE; }
    … …  
    156177    // ** Choose the calibration (depends on the mode : do I have an image reference or not?)
    157178    float Mcal = 0, Mmos = 0, Mgrid = 0;
    158     if (getImageEntry (measureOffset + k, cat) < 0) {
     179    off_t meas = measureOffset + k;
     180    if (getImageEntry (meas, cat) < 0) {
    159181      // measurements without an image are either external reference photometry or
    160182      // data for which the associated image has not been loaded (probably because of
    … …  
    173195    // XXX ubercal dist should be by Nsec
    174196    int myUbercalDist = getUbercalDist(meas, cat);
    175     result->minUbercalDist[Nsec] = MIN(result->minUbercalDist[Nsec], myUbercalDist);
     197    results->minUbercalDist[Nsec] = MIN(results->minUbercalDist[Nsec], myUbercalDist);
    176198
    177199    int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
    … …  
    184206    if (isSetMrelFinal) {
    185207      Map = PhotCat (&measure[k], MAG_CLASS_APER);
    186       dMap = MAX (hypot(measureT[k].dMap,   code->photomErrSys), MIN_ERROR);
     208      dMap = MAX (hypot(measure[k].dMap,   code->photomErrSys), MIN_ERROR);
    187209      if (CHECK_VALID_MAG(Map)) {
    188         // XXX need to carry a counter for each list set
    189         int Nap = results->apData[Nsec].Nlist;
    190         results->apData[Nsec].flxlist[Nap] = Map - Mcal - Mmos - Mgrid;
    191         results->apData[Nsec].errlist[Nap] = dMpsf; // no other option
    192         results->apData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
    193         results->apData[Nsec].ranking[Nap] = measureRank[k];
    194         results->apData[Nsec].measSeq[Nap] = k;
    195         results->apData[Nsec].Nlist ++;
     210        int Nap = results->aperData[Nsec].Nlist;
     211        results->aperData[Nsec].flxlist[Nap] = Map - Mcal - Mmos - Mgrid;
     212        results->aperData[Nsec].errlist[Nap] = dMap;
     213        results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
     214        results->aperData[Nsec].ranking[Nap] = measureRank[k];
     215        results->aperData[Nsec].measSeq[Nap] = k;
     216        results->aperData[Nsec].Nlist ++;
    196217      }
    197218    }
    … …  
    200221    if (isSetMrelFinal) {
    201222      Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
    202       dMkron = MAX (hypot(measureT[k].dMkron, code->photomErrSys), MIN_ERROR);
     223      dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR);
    203224      if (CHECK_VALID_MAG(Mkron)) {
    204225        int Nkron = results->kronData[Nsec].Nlist;
    205226        results->kronData[Nsec].flxlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
    206         results->kronData[Nsec].errlist[Nkron] = measure[k].dMkron;
     227        results->kronData[Nsec].errlist[Nkron] = dMkron;
    207228        results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    208229        results->kronData[Nsec].ranking[Nkron] = measureRank[k];
    … …  
    237258      if (dMagAp > SigmaAp) {
    238259        results->Next[Nsec] ++;
    239         results->NextPS1 ++;
     260        NextPS1 ++;
    240261      } else {
    241         results->NpsfPS1 ++;
     262        NpsfPS1 ++;
    242263      }
    243264    }
    … …  
    246267  // now calculate the mean stats for the Nsec bands.
    247268  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    248     dvo_init_secfilt (&secfilt[Nsec]);
     269    dvo_secfilt_init (&secfilt[Nsec]); // XXX does this reset astrometry bits? no, right?
    249270
    250271    // XXX hardwired
    251     if ((thisCode < 6) || (thisCode == 9)) {
     272    if ((Nsec < 4) || (Nsec == 8)) {
    252273      secfilt[Nsec].Ncode = results->NexpPS1[Nsec];
    253274    } else {
    254       secfilt[Nsec].Ncode = results->Ncode[Nsec]; // 2MASS data if it exists
     275      secfilt[Nsec].Ncode = results->Nmeas[Nsec]; // 2MASS data if it exists
     276    }
     277
     278    if (haveTYCHO) {
     279      // secfilt[Nsec].flags |= ID_SECF_HAVE_TYCHO;
    255280    }
    256281
    … …  
    263288      float ZP = 0.0;
    264289
    265       if (!haveTycho) {
     290      // if we have loaded the synth zpt correction, and do NOT have a tycho measurement, apply the correction
     291      if (!haveTYCHO && synthzpts) {
    266292        // need to look up the (X,Y) coords from this (R,D) location
    267293        double X, Y;
    268294        double R = ohana_normalize_angle_to_midpoint (measure->R, 0.0); // XXX 0.0 or 180.0 for center?
    269         RD_to_XY (&X, &Y, R, measure->D, &zpts->coords);
     295        RD_to_XY (&X, &Y, R, measure->D, &synthzpts->coords);
    270296        if (X < 0) return FALSE;
    271297        if (Y < 0) return FALSE;
    272         if (X >= zpts->Nx) return FALSE;
    273         if (Y >= zpts->Ny) return FALSE;
     298        if (X >= synthzpts->Nx) return FALSE;
     299        if (Y >= synthzpts->Ny) return FALSE;
    274300       
    275301        int Xpix = X;
    276302        int Ypix = Y;
    277         int Npix = Xpix + Ypix*zpts->Nx;
     303        int Npix = Xpix + Ypix*synthzpts->Nx;
    278304       
    279         float *value = (float *) zpts->matrix[Nsec].buffer;
     305        float *value = (float *) synthzpts->matrix[Nsec].buffer;
    280306        ZP = !isnan(value[Npix]) ? value[Npix] : 0.0;
    281         secfilt[Nsec].flags |= ID_SECF_HAVE_TYCHO;
    282307      }
    283308      secfilt[Nsec].M      = Mpsf + ZP;
    284309      secfilt[Nsec].dM     = 0.6;
    285       secfilt[Nsec].chipsq = 0.0;
     310      secfilt[Nsec].Mchisq = 0.0;
    286311      secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
    287312      continue;
    … …  
    289314
    290315    // if too few valid measurements meet the minimum criteria, go to the next entry
    291     int Nranking = magStatsByRanking (&results->psfData[Nsec], &results->psfstats);
     316    StatType *psfstats = &results->psfstats;
     317    int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats);
    292318    if (Nranking < Nminmeas) {
    293319      secfilt[Nsec].flags |= ID_STAR_FEW;
    … …  
    300326
    301327    // when running -averages, we have no information about the images, so we cannot set this
    302     if (minUbercalDist > -1) {
    303       secfilt[Nsec].ubercalDist = minUbercalDist;
     328    if (results->minUbercalDist[Nsec] > -1) {
     329      secfilt[Nsec].ubercalDist = results->minUbercalDist[Nsec];
    304330    }
    305331
    306332    if (isSetMrelFinal) {
    307       secfilt[Nsec].Mstdev = psfstats->sigma;
    308       secfilt[Nsec].Nused  = psfstats->Nmeas;
    309       secfilt[Nsec].Mmax = psfstats->max;
    310       secfilt[Nsec].Mmin = psfstats->min;
    311 
    312       secfilt[Nsec].psfQfMax     = psfQfMax;
    313       secfilt[Nsec].psfQfPerfMax = psfQfPerfMax;
    314 
    315       Nranking = magStatsByRanking (&results->aperData[Nsec], results->apstats);
    316       secfilt[Nsec].Map     = apstats->mean;
    317       secfilt[Nsec].dMap    = apstats->error;
    318       secfilt[Nsec].sMap    = apstats->sigma;
    319       secfilt[Nsec].NusedAp = Nranking;
    320 
     333      if (Nranking) {
     334        secfilt[Nsec].Mstdev = psfstats->sigma;
     335        secfilt[Nsec].Nused  = psfstats->Nmeas;
     336        secfilt[Nsec].Mmax   = psfstats->max;
     337        secfilt[Nsec].Mmin   = psfstats->min;
     338      }
     339
     340      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
     341      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
     342
     343      StatType *apstats = &results->apstats;
     344      Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
     345      if (Nranking) {
     346        secfilt[Nsec].Map     = apstats->mean;
     347        secfilt[Nsec].dMap    = apstats->error;
     348        secfilt[Nsec].sMap    = apstats->sigma;
     349        secfilt[Nsec].NusedAp = Nranking;
     350      }
     351
     352      StatType *kronstats = &results->kronstats;
    321353      Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
    322       secfilt[Nsec].Mkron     = kronstats->mean;
    323       secfilt[Nsec].dMkron    = kronstats->error;
    324       secfilt[Nsec].sMkron    = kronstats->sigma;
    325       secfilt[Nsec].NusedKron = Nranking;
     354      if (Nranking) {
     355        secfilt[Nsec].Mkron     = kronstats->mean;
     356        secfilt[Nsec].dMkron    = kronstats->error;
     357        secfilt[Nsec].sMkron    = kronstats->sigma;
     358        secfilt[Nsec].NusedKron = Nranking;
     359      }
    326360
    327361      // does this object appear extended in > 50% of measurements?
    328       if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*results->Npsf[Nsec])) {
     362      if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*psfstats->Nmeas)) {
    329363        secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
    330364      }
    … …  
    362396
    363397  if (isSetMrelFinal) {
    364     if (results->NextPS1 && (results->NextPS1 > results->NpsfPS1)) {
     398    if (NextPS1 && (NextPS1 > NpsfPS1)) {
    365399      average[0].flags |= ID_OBJ_EXT;
    366400    }
    … …  
    378412}
    379413
     414# if (0)
    380415// set mean of chip measurements (selected by photcode range for now):
    381416int setMrelAverageExposure_old (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
    … …  
    750785  return (TRUE);
    751786}
     787# endif
    752788
    753789// only apply Stack operation on setMrelFinal in first pass
    … …  
    756792// 2) select the BEST detections per filter (regardless of PRIMARY)
    757793// 3) apply the zero point and AB->Jy transformations
    758 int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure) {
     794int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
     795
     796  // we are guaranteed to have average, measure, secfilt
     797  Average *average = &catalog[cat].average[ave];
     798  off_t measureOffset = average->measureOffset;
     799  int Nmeasure = average->Nmeasure;
     800
     801  Measure *measure = &catalog[cat].measure[measureOffset];
     802  SecFilt *secfilt = &catalog[cat].secfilt[ave*Nsecfilt];
    759803
    760804  off_t k;
    761805
    762   float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0;
     806  float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0;
    763807
    764808  // set the primary projection cell and skycell for this coordinate
    … …  
    801845
    802846    off_t meas = measureOffset;
    803     for (k = 0; k < average[0].Nmeasure; k++, meas++) {
     847    for (k = 0; k < Nmeasure; k++, meas++) {
    804848
    805849      // only examine gpc1 stack data
    … …  
    9621006    average[0].flags |= ID_OBJ_BAD_STACK;
    9631007  }
     1008
     1009  setGlobalObjStats (average, measure);
     1010
    9641011  return (TRUE);
    9651012}
    … …  
    9741021// somewhat simplified relative to chip-photometry:
    9751022// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
    976 int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure) {
     1023int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
     1024
     1025# if (0)
     1026
     1027  off_t measureOffset = average->measureOffset;
     1028
     1029  // we are guaranteed to have average, measure, secfilt
     1030  Average *average = &catalog[cat].average[ave];
     1031  Measure *measure = &catalog[cat].measure[measureOffset];
     1032  SecFilt *secfilt = &catalog[cat].secfilt[ave*Nsecfilt];
    9771033
    9781034  off_t k;
    979   float Mcal= 0;
     1035  float Mcal = 0;
    9801036
    9811037  // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
    … …  
    11801236    secfilt[Nsec].NwarpGood = NwarpGood;
    11811237  }
     1238
     1239  setGlobalObjStats (average, measure);
     1240
     1241# endif
     1242
    11821243  return (TRUE);
    11831244}
    … …  
    12201281  // we have a list of measurements for this Nsec value, along with errors, weights, and their rank
    12211282  // first, sort by the rank (increasing)
    1222   sort_measures_by_ranking (dataset->ranking, dataset->measSeq, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->Nlist);
     1283  sort_StatDataSet (dataset);
    12231284
    12241285  // find the values with the same minimum rank:
    … …  
    12291290  for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);
    12301291
    1231   // if no valid measurements meet the minimum criteria, go to the next entry
    1232   if (Nranking < Nminmeas) {
    1233     return 0;
    1234   }
    1235 
    12361292  liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);
    12371293  return (Nranking);
    12381294}
     1295
     1296/* sort a coordinate pair (X,Y) and the associated index (S) */
     1297void sort_StatDataSet (StatDataSet *dataset) {
     1298 
     1299# define SWAPFUNC(A,B){ double dtmp; int itmp; \
     1300    dtmp = dataset->flxlist[A]; dataset->flxlist[A] = dataset->flxlist[B]; dataset->flxlist[B] = dtmp;     \
     1301    dtmp = dataset->errlist[A]; dataset->errlist[A] = dataset->errlist[B]; dataset->errlist[B] = dtmp;     \
     1302    dtmp = dataset->wgtlist[A]; dataset->wgtlist[A] = dataset->wgtlist[B]; dataset->wgtlist[B] = dtmp;     \
     1303    itmp = dataset->ranking[A]; dataset->ranking[A] = dataset->ranking[B]; dataset->ranking[B] = itmp;     \
     1304    itmp = dataset->measSeq[A]; dataset->measSeq[A] = dataset->measSeq[B]; dataset->measSeq[B] = itmp;     \
     1305  }
     1306# define COMPARE(A,B)(dataset->ranking[A] < dataset->ranking[B])
     1307 
     1308  OHANA_SORT (dataset->Nlist, COMPARE, SWAPFUNC);
     1309 
     1310# undef SWAPFUNC
     1311# undef COMPARE
     1312
     1313}
     1314

  • branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelFinal.c

    r38293 r38297  
    102102
    103103  ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure);
    104   setMeasurementRank (catalog);
    105   setMrelOutput (catalog);
     104  setMeasureRank (catalog);
     105  setMrelOutput (catalog, 1, flatcorr);
    106106  if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr);
    107107
    … …  
    311311      // measurements ranked by inst mag limit (not REF mags)
    312312      if (ImagSelect) {
    313         mag = PhotInst (&measure[i], MAG_CLASS_PSF);
     313        float mag = PhotInst (&measure[i], MAG_CLASS_PSF);
    314314        if (mag < ImagMin)                           { measureRank[i] = 8; continue; }
    315315        if (mag > ImagMax)                           { measureRank[i] = 8; continue; }
Note: See TracChangeset for help on using the changeset viewer.