IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 24, 2012, 3:52:04 PM (14 years ago)
Author:
eugene
Message:

merge changes from trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120905/ppTranslate/src/ppMopsRead.c

    r34285 r34566  
     1/*
     2 *
     3 *
     4 *
     5 *
     6 *
     7 *
     8 *
     9 *
     10 *
     11 *
     12 *
     13 *
     14 *
     15 *
     16 *
     17 *
     18 *
     19 *
     20 *
     21 *
     22 *
     23 *
     24 *
     25 *
     26 *
     27 *
     28 *
     29 *
     30 * THIS CODE CONTAINS DIRTY BUG FIXES!
     31 *
     32 *
     33 *
     34 *
     35 *
     36 *
     37 *
     38 * Look for "DIRTY_BUG_FIX" to know where
     39 *
     40 * and set the following to 0 to change the code
     41 */
     42#define DIRTY_BUG_FIX 1
     43/*
     44 *
     45 *
     46 *
     47 *
     48 *
     49 *
     50 *
     51 *
     52 *
     53 *
     54 *
     55 *
     56 *
     57 *
     58 *
     59 *
     60 *
     61 *
     62 *
     63 *
     64 *
     65 *
     66 *
     67 *
     68 *
     69 *
     70 *
     71 *
     72 *
     73 *
     74 *
     75 *
     76 *
     77 */
     78
    179#ifdef HAVE_CONFIG_H
    280#include <config.h>
    … …  
    886#include "ppMops.h"
    987
    10 static void addDummyValues(psMetadata* md, long size);
     88static void addDummyValues(psMetadata* md, long size, int version);
    1189static void replaceDummyValuesF32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes);
     90static void replaceDummyValuesS32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes);
     91static void renameDummyValuesF32(const char* colNameSource, psMetadata* source, const char* colNameTarget, psMetadata* target, psVector* indexes);
     92#if DIRTY_BUG_FIX
     93static psVector* createDummyF32_default(long size, float default_value);
     94#endif
    1295
    1396/*
    … …  
    24107  for (int i = 0; i < num; i++) {
    25108    const char *name = inNames->data[i];
     109    printf("%s\n", name);
    26110
    27111    psFits *fits = psFitsOpen(name,  "r"); // FITS file
    … …  
    95179    det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
    96180    det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;
    97 
    98     det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
    99                                  psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    100 
     181    //The global SEEING is the mean of all defined FWHM_MAJ
     182    det->seeing = psMetadataLookupF32(NULL, header, "FWHM_MAJ");
    101183    det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
    102184    det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    103 
    104185    psFree(header);
    105186
    … …  
    116197    if (skyChipPsfVersion >= 2) {
    117198      //First: append all the new columns that we want to the existing table
    118       addDummyValues(det->table, size);
     199      addDummyValues(det->table, size, args->version);
    119200      if (!psFitsMoveExtName(fits, "SkyChip.xfit")) {
    120201        psTrace("ppMops.read", 3, "No fitted trails extension");
    121202      } else {
    122203        psTrace("ppMops.read", 3, "Fitted trails extension found\n");
     204        psTrace("ppMops.read", 10, "Getting size?\n");
    123205        int fittedTrailsSize = psFitsTableSize(fits);
     206        psTrace("ppMops.read", 10, "size = %d\n", fittedTrailsSize);
    124207        if (fittedTrailsSize <= 0) {
    125208          psErrorStackPrint(stderr, "Unable to determine size of fitted trails extension table %d", i);
    … …  
    127210          psErrorClear();
    128211        } else {
     212          psTrace("ppMops.read", 10, "Reading table of interest?\n");
    129213          psMetadata* fittedTrails = psFitsReadTableAllColumns(fits); // Table of interest
     214          psTrace("ppMops.read", 10, "OK for table of interest\n");
    130215          if (!fittedTrails) {
    131216            psError(PS_ERR_IO, false, "Unable to read fittedTrails table in file %d", i);
    … …  
    133218          }
    134219          //Iterate on the different names and types expected in the fittedTrails parameters
     220          psTrace("ppMops.read", 10, "Getting IPP_IDET\n");
    135221          psVector* idet = psMetadataLookupVector(NULL, fittedTrails, "IPP_IDET");
    136222          replaceDummyValuesF32("X_EXT", fittedTrails, det->table, idet);
    … …  
    140226          replaceDummyValuesF32("EXT_INST_MAG",  fittedTrails, det->table, idet);
    141227          replaceDummyValuesF32("EXT_INST_MAG_SIG",  fittedTrails, det->table, idet);
    142           replaceDummyValuesF32("NPARAMS",  fittedTrails, det->table, idet);
     228          replaceDummyValuesS32("NPARAMS",  fittedTrails, det->table, idet);
    143229          replaceDummyValuesF32("EXT_WIDTH_MAJ",  fittedTrails, det->table, idet);
    144230          replaceDummyValuesF32("EXT_WIDTH_MIN",  fittedTrails, det->table, idet);
    145231          replaceDummyValuesF32("EXT_THETA",  fittedTrails, det->table, idet);
    146232          replaceDummyValuesF32("EXT_WIDTH_MAJ_ERR",  fittedTrails, det->table, idet);
    147           replaceDummyValuesF32("EXT_WIDTH_MIN_ERR",  fittedTrails, det->table, idet);
     233          // EXT_WIDTH_MIN_ERR is actually undefined but set to 0. in CMF files
     234          // We explicitely let it set to NaN, hence the commented out
     235          // following line
     236          //replaceDummyValuesF32("EXT_WIDTH_MIN_ERR",  fittedTrails, det->table, idet);
    148237          replaceDummyValuesF32("EXT_THETA_ERR",  fittedTrails, det->table, idet);
     238          psTrace("ppMops.read", 10, "Got all version 2 values\n");
     239          if (skyChipPsfVersion >= 3) {
     240            psTrace("ppMops.read", 10, "Getting PS1_DV3 data\n");
     241            replaceDummyValuesF32("RA_EXT",  fittedTrails, det->table, idet);
     242            // Note: RA_EXT_SIGMA is computed (see det->raExtErr)
     243            replaceDummyValuesF32("DEC_EXT",  fittedTrails, det->table, idet);
     244            // Note: DEC_EXT_SIGMA is computed (see det->decExtErr)
     245            //Note: POSANG_EXT is POSANGLE in xfit
     246#if DIRTY_BUG_FIX
     247#else
     248            renameDummyValuesF32("POSANGLE",  fittedTrails, "POSANG_EXT", det->table, idet);
     249#endif   
     250            //Note: PLTSCALE_EXT in xfit is PLTSCALE
     251            renameDummyValuesF32("PLTSCALE", fittedTrails, "PLTSCALE_EXT", det->table, idet);
     252            renameDummyValuesF32("EXT_INST_FLUX", fittedTrails, "EXT_FLUX", det->table, idet);
     253            replaceDummyValuesF32("EXT_CAL_MAG", fittedTrails, det->table, idet);
     254            renameDummyValuesF32("EXT_INST_MAG_SIG", fittedTrails, "EXT_MAG_SIG", det->table, idet);
     255            replaceDummyValuesF32("EXT_CHISQ", fittedTrails, det->table, idet);
     256            replaceDummyValuesS32("EXT_NDOF", fittedTrails, det->table, idet);
     257            psTrace("ppMops.read", 10, "Got all version 3 values\n");
     258          }
    149259        }
    150260      }
    151261    }
    152 
    153262    psFitsClose(fits);
    154263
    … …  
    161270    psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
    162271    psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
    163 
    164272    det->raErr = psVectorAlloc(size, PS_TYPE_F64);
    165273    det->decErr = psVectorAlloc(size, PS_TYPE_F64);
     274    if (skyChipPsfVersion >= 3) {
     275      det->raExtErr = psVectorAlloc(size, PS_TYPE_F64);
     276      det->decExtErr = psVectorAlloc(size, PS_TYPE_F64);
     277    } else {
     278      det->raExtErr = NULL;
     279      det->decExtErr = NULL;
     280    }
    166281    det->mask = psVectorAlloc(size, PS_TYPE_U8);
    167282
    … …  
    177292
    178293    // Add our new vectors to the table so that duplicates and masked items may be purged
     294    psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
    179295    psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
    180     psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
     296    if (skyChipPsfVersion >= 3) {
     297      psMetadataAddVector(table, PS_LIST_HEAD, "RA_EXT_SIGMA", 0, NULL, det->raExtErr);
     298      psMetadataAddVector(table, PS_LIST_HEAD, "DEC_EXT_SIGMA", 0, NULL, det->decExtErr);
     299    }
    181300
    182301    psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
    … …  
    189308    psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
    190309    psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
     310    psVector *xExtErrV = NULL;
     311    psVector *yExtErrV = NULL;
     312    psVector *scaleExtV = NULL;
     313    psVector *angleExtV = NULL;
     314    if (skyChipPsfVersion >= 3) {
     315      xExtErrV = psMetadataLookupVector(NULL, table, "X_EXT_SIG");
     316      yExtErrV = psMetadataLookupVector(NULL, table, "Y_EXT_SIG");
     317      scaleExtV = psMetadataLookupVector(NULL, table, "PLTSCALE_EXT");
     318      angleExtV = psMetadataLookupVector(NULL, table, "POSANG_EXT");
     319    }
    191320
    192321    double plateScale = 0.0;        // Plate scale
    193322    long numGood = 0;               // Number of good rows
    194323    for (long row = 0; row < size; row++) {
    195 
    196324      psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
    197325      if (flags & SOURCE_MASK) {
    … …  
    205333      double yErr = yErrV->data.F32[row];
    206334      double scale = scaleV->data.F32[row];
     335#if DIRTY_BUG_FIX
     336      double angle = 0. * angleV->data.F32[row];
     337#else
    207338      double angle = angleV->data.F32[row];
     339#endif
    208340
    209341      if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
    … …  
    233365      plateScale += scale;
    234366      numGood++;
    235     }
    236     det->seeing *= ((float) plateScale) / ((float) numGood);
    237 
    238     // Are we using numGood for anything outside of this function?
     367
     368      // Same for EXT data if version permits
     369      if (skyChipPsfVersion >= 3) {
     370        xErr = xExtErrV->data.F32[row];
     371        yErr = yExtErrV->data.F32[row];
     372        scale = scaleExtV->data.F32[row]/3600.;
     373#if DIRTY_BUG_FIX
     374        angle = 0.;
     375#else
     376        angle = angleExtV->data.F32[row];
     377#endif
     378        cosAngle = cos(angle);
     379        sinAngle = sin(angle);
     380        cosAngle2 = PS_SQR(cosAngle);
     381        sinAngle2 = PS_SQR(sinAngle);
     382        xErr2 = PS_SQR(xErr);
     383        yErr2 = PS_SQR(yErr);
     384        errScale = scale / 3600.0;
     385        det->raExtErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     386        det->decExtErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     387      }
     388    }
    239389    det->numGood = numGood;
    240390
    … …  
    256406  return dummy;
    257407}
    258 
    259 static void addDummyValues(psMetadata* md, long size) {
     408#if DIRTY_BUG_FIX
     409static psVector* createDummyF32_default(long size, float default_value) {
     410  psVector* dummy = psVectorAlloc(size, PS_TYPE_F32);
     411  psVectorInit(dummy, default_value);
     412  return dummy;
     413}
     414#endif
     415static psVector* createDummyS32(long size) {
     416  psVector* dummy = psVectorAlloc(size, PS_TYPE_S32);
     417  psVectorInit(dummy, NAN);
     418  return dummy;
     419}
     420
     421static void addDummyValues(psMetadata* md, long size, int version) {
    260422  psMetadataAdd(md, PS_LIST_TAIL, "X_EXT", PS_DATA_VECTOR, "EXT model x coordinate", createDummyF32(size));
    261423  psMetadataAdd(md, PS_LIST_TAIL, "Y_EXT", PS_DATA_VECTOR, "EXT model y coordinate", createDummyF32(size));
    … …  
    264426  psMetadataAdd(md, PS_LIST_TAIL, "EXT_INST_MAG",  PS_DATA_VECTOR, "EXT fit instrumental magnitude", createDummyF32(size));
    265427  psMetadataAdd(md, PS_LIST_TAIL, "EXT_INST_MAG_SIG",  PS_DATA_VECTOR, "Sigma of PSF instrumental magnitude", createDummyF32(size));
    266   psMetadataAdd(md, PS_LIST_TAIL, "NPARAMS",  PS_DATA_VECTOR, "Number of model parameters", createDummyF32(size));
     428  psMetadataAdd(md, PS_LIST_TAIL, "NPARAMS",  PS_DATA_VECTOR, "Number of model parameters", createDummyS32(size));
    267429  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MAJ",  PS_DATA_VECTOR, "EXT width (major axis), length for trail", createDummyF32(size));
    268430  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MIN",  PS_DATA_VECTOR, "EXT width (minor axis), sigma for trail", createDummyF32(size));
    … …  
    271433  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",  PS_DATA_VECTOR, "EXT width error (minor axis)", createDummyF32(size));
    272434  psMetadataAdd(md, PS_LIST_TAIL, "EXT_THETA_ERR",  PS_DATA_VECTOR, "EXT orientation angle (error)", createDummyF32(size));
     435  if (version >= 3) {
     436    psTrace("ppMops.read", 10, "Adding columns for version PS1_DV%d\n", version);
     437    psMetadataAdd(md, PS_LIST_TAIL, "RA_EXT",  PS_DATA_VECTOR, "Fitted centroid RA", createDummyF32(size));
     438/*     psMetadataAdd(md, PS_LIST_TAIL, "RA_EXT_SIGMA",  PS_DATA_VECTOR, "Fitted RA sigma", createDummyF32(size)); *\/ */
     439    psMetadataAdd(md, PS_LIST_TAIL, "DEC_EXT",  PS_DATA_VECTOR, "Fitted centroid DEC", createDummyF32(size));
     440/*     psMetadataAdd(md, PS_LIST_TAIL, "DEC_EXT_SIGMA",  PS_DATA_VECTOR, "Fitted DEC sigma", createDummyF32(size)); *\/ */
     441#if DIRTY_BUG_FIX
     442    psMetadataAdd(md, PS_LIST_TAIL, "POSANG_EXT",  PS_DATA_VECTOR, "Fitted position angle", createDummyF32_default(size, 0.));
     443#else
     444    psMetadataAdd(md, PS_LIST_TAIL, "POSANG_EXT",  PS_DATA_VECTOR, "Fitted position angle", createDummyF32(size));
     445#endif
     446    psMetadataAdd(md, PS_LIST_TAIL, "PLTSCALE_EXT",  PS_DATA_VECTOR, "Plate scale at centroid", createDummyF32(size));
     447    psMetadataAdd(md, PS_LIST_TAIL, "EXT_FLUX",  PS_DATA_VECTOR, "Fitted flux", createDummyF32(size));
     448    psMetadataAdd(md, PS_LIST_TAIL, "EXT_CAL_MAG",  PS_DATA_VECTOR, "Calibrated mag", createDummyF32(size));
     449    psMetadataAdd(md, PS_LIST_TAIL, "EXT_MAG_SIG",  PS_DATA_VECTOR, "Mag sigma", createDummyF32(size));
     450    psMetadataAdd(md, PS_LIST_TAIL, "EXT_CHISQ",  PS_DATA_VECTOR, "Chi^2 of fit", createDummyF32(size));
     451    psMetadataAdd(md, PS_LIST_TAIL, "EXT_NDOF",  PS_DATA_VECTOR, "Fit degrees of freedom", createDummyS32(size));
     452  }
    273453}
    274454
    275455static void replaceDummyValuesF32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes) {
     456  psTrace("ppMops.read", 10, "Trying to get column %s\n", colName);
    276457  psVector* source_vector = psMetadataLookupVector(NULL, source, colName);
    277458  psVector* target_vector = psMetadataLookupVector(NULL, target, colName);
    … …  
    279460    target_vector->data.F32[indexes->data.S64[index]] = source_vector->data.F32[index];
    280461  }
    281 }
     462  psTrace("ppMops.read", 10, "OK for %s\n", colName);
     463}
     464
     465static void replaceDummyValuesS32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes) {
     466  psTrace("ppMops.read", 10, "Trying to get column %s\n", colName);
     467  psVector* source_vector = psMetadataLookupVector(NULL, source, colName);
     468  psVector* target_vector = psMetadataLookupVector(NULL, target, colName);
     469  for (long index = 0; index<indexes->n; ++index) {
     470    target_vector->data.S32[indexes->data.S64[index]] = source_vector->data.S32[index];
     471  }
     472  psTrace("ppMops.read", 10, "OK for %s\n", colName);
     473}
     474
     475static void renameDummyValuesF32(const char* colNameSource, psMetadata* source, const char* colNameTarget, psMetadata* target, psVector* indexes) {
     476  psTrace("ppMops.read", 10, "Trying to rename column %s to %s\n", colNameSource, colNameTarget);
     477  psVector* source_vector = psMetadataLookupVector(NULL, source, colNameSource);
     478  psVector* target_vector = psMetadataLookupVector(NULL, target, colNameTarget);
     479  for (long index = 0; index<indexes->n; ++index) {
     480    target_vector->data.F32[indexes->data.S64[index]] = source_vector->data.F32[index];
     481  }
     482  psTrace("ppMops.read", 10, "OK for %s (renamed by %s)\n", colNameSource, colNameTarget);
     483}
Note: See TracChangeset for help on using the changeset viewer.