IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 6, 2012, 6:32:19 PM (14 years ago)
Author:
watersc1
Message:

merge from trunk. Preliminary versions of stacktool/warptool updates and regenerate_background.pl script. The warp code is finished and tested, and I still need to get the stack version resolved.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20120906/ppTranslate/src/ppMopsRead.c

    r34285 r34772  
     1#define DIRTY_BUG_FIX 0
     2
    13#ifdef HAVE_CONFIG_H
    24#include <config.h>
    … …  
    810#include "ppMops.h"
    911
    10 static void addDummyValues(psMetadata* md, long size);
     12static void addDummyValues(psMetadata* md, long size, int version);
    1113static void replaceDummyValuesF32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes);
     14static void replaceDummyValuesS32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes);
     15static void renameDummyValuesF32(const char* colNameSource, psMetadata* source, const char* colNameTarget, psMetadata* target, psVector* indexes);
    1216
    1317/*
    … …  
    2428  for (int i = 0; i < num; i++) {
    2529    const char *name = inNames->data[i];
     30    printf("%s\n", name);
    2631
    2732    psFits *fits = psFitsOpen(name,  "r"); // FITS file
    … …  
    6166    if (skyChipPsfVersion != args->version) {
    6267      if (skyChipPsfVersion > args->version) {
    63         psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
    64                   skyChipPsfVersion, args->version);
     68        psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
     69                  skyChipPsfVersion, args->version);
    6570      } else { // Necessarily: skyChipPsfVersion > args->version
    66         psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
    67                   skyChipPsfVersion, args->version);
     71        psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
     72                  skyChipPsfVersion, args->version);
    6873      }
    6974    }
    … …  
    95100    det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
    96101    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 
     102    //The global SEEING is the mean of all defined FWHM_MAJ
     103    det->seeing = psMetadataLookupF32(NULL, header, "FWHM_MAJ");
    101104    det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
    102105    det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    103 
    104106    psFree(header);
    105107
    … …  
    116118    if (skyChipPsfVersion >= 2) {
    117119      //First: append all the new columns that we want to the existing table
    118       addDummyValues(det->table, size);
     120      addDummyValues(det->table, size, args->version);
    119121      if (!psFitsMoveExtName(fits, "SkyChip.xfit")) {
    120         psTrace("ppMops.read", 3, "No fitted trails extension");
     122        psTrace("ppMops.read", 3, "No fitted trails extension");
    121123      } else {
    122         psTrace("ppMops.read", 3, "Fitted trails extension found\n");
    123         int fittedTrailsSize = psFitsTableSize(fits);
    124         if (fittedTrailsSize <= 0) {
    125           psErrorStackPrint(stderr, "Unable to determine size of fitted trails extension table %d", i);
    126           psTrace("ppMops.read", 3, "No entry in fitted trails extension!!!!\n");
    127           psErrorClear();
    128         } else {
    129           psMetadata* fittedTrails = psFitsReadTableAllColumns(fits); // Table of interest
    130           if (!fittedTrails) {
    131             psError(PS_ERR_IO, false, "Unable to read fittedTrails table in file %d", i);
    132             return NULL;
    133           }
    134           //Iterate on the different names and types expected in the fittedTrails parameters
    135           psVector* idet = psMetadataLookupVector(NULL, fittedTrails, "IPP_IDET");
    136           replaceDummyValuesF32("X_EXT", fittedTrails, det->table, idet);
    137           replaceDummyValuesF32("Y_EXT", fittedTrails, det->table, idet);
    138           replaceDummyValuesF32("X_EXT_SIG", fittedTrails, det->table, idet);
    139           replaceDummyValuesF32("Y_EXT_SIG", fittedTrails, det->table, idet);
    140           replaceDummyValuesF32("EXT_INST_MAG",  fittedTrails, det->table, idet);
    141           replaceDummyValuesF32("EXT_INST_MAG_SIG",  fittedTrails, det->table, idet);
    142           replaceDummyValuesF32("NPARAMS",  fittedTrails, det->table, idet);
    143           replaceDummyValuesF32("EXT_WIDTH_MAJ",  fittedTrails, det->table, idet);
    144           replaceDummyValuesF32("EXT_WIDTH_MIN",  fittedTrails, det->table, idet);
    145           replaceDummyValuesF32("EXT_THETA",  fittedTrails, det->table, idet);
    146           replaceDummyValuesF32("EXT_WIDTH_MAJ_ERR",  fittedTrails, det->table, idet);
    147           replaceDummyValuesF32("EXT_WIDTH_MIN_ERR",  fittedTrails, det->table, idet);
    148           replaceDummyValuesF32("EXT_THETA_ERR",  fittedTrails, det->table, idet);
    149         }
    150       }
    151     }
    152 
     124        psTrace("ppMops.read", 3, "Fitted trails extension found\n");
     125        psTrace("ppMops.read", 10, "Getting size?\n");
     126        int fittedTrailsSize = psFitsTableSize(fits);
     127        psTrace("ppMops.read", 10, "size = %d\n", fittedTrailsSize);
     128        if (fittedTrailsSize <= 0) {
     129          psErrorStackPrint(stderr, "Unable to determine size of fitted trails extension table %d", i);
     130          psTrace("ppMops.read", 3, "No entry in fitted trails extension!!!!\n");
     131          psErrorClear();
     132        } else {
     133          psTrace("ppMops.read", 10, "Reading table of interest?\n");
     134          psMetadata* fittedTrails = psFitsReadTableAllColumns(fits); // Table of interest
     135          psTrace("ppMops.read", 10, "OK for table of interest\n");
     136          if (!fittedTrails) {
     137            psError(PS_ERR_IO, false, "Unable to read fittedTrails table in file %d", i);
     138            return NULL;
     139          }
     140          //Iterate on the different names and types expected in the fittedTrails parameters
     141          psTrace("ppMops.read", 10, "Getting IPP_IDET\n");
     142          psVector* idet = psMetadataLookupVector(NULL, fittedTrails, "IPP_IDET");
     143          replaceDummyValuesF32("X_EXT", fittedTrails, det->table, idet);
     144          replaceDummyValuesF32("Y_EXT", fittedTrails, det->table, idet);
     145          replaceDummyValuesF32("X_EXT_SIG", fittedTrails, det->table, idet);
     146          replaceDummyValuesF32("Y_EXT_SIG", fittedTrails, det->table, idet);
     147          replaceDummyValuesF32("EXT_INST_MAG",  fittedTrails, det->table, idet);
     148          replaceDummyValuesF32("EXT_INST_MAG_SIG",  fittedTrails, det->table, idet);
     149          replaceDummyValuesS32("NPARAMS",  fittedTrails, det->table, idet);
     150          replaceDummyValuesF32("EXT_WIDTH_MAJ",  fittedTrails, det->table, idet);
     151          replaceDummyValuesF32("EXT_WIDTH_MIN",  fittedTrails, det->table, idet);
     152          replaceDummyValuesF32("EXT_THETA",  fittedTrails, det->table, idet);
     153          replaceDummyValuesF32("EXT_WIDTH_MAJ_ERR",  fittedTrails, det->table, idet);
     154          // EXT_WIDTH_MIN_ERR is actually undefined but set to 0. in CMF files
     155          // We explicitely let it set to NaN, hence the commented out
     156          // following line
     157          //replaceDummyValuesF32("EXT_WIDTH_MIN_ERR",  fittedTrails, det->table, idet);
     158          replaceDummyValuesF32("EXT_THETA_ERR",  fittedTrails, det->table, idet);
     159          psTrace("ppMops.read", 10, "Got all version 2 values\n");
     160          if (skyChipPsfVersion >= 3) {
     161            psTrace("ppMops.read", 10, "Getting PS1_DV3 data\n");
     162            replaceDummyValuesF32("RA_EXT",  fittedTrails, det->table, idet);
     163            // Note: RA_EXT_SIGMA is computed (see det->raExtErr)
     164            replaceDummyValuesF32("DEC_EXT",  fittedTrails, det->table, idet);
     165            // Note: DEC_EXT_SIGMA is computed (see det->decExtErr)
     166            //Note: POSANG_EXT is POSANGLE in xfit
     167            renameDummyValuesF32("POSANGLE",  fittedTrails, "POSANG_EXT", det->table, idet);
     168            //Not written anymore but still used in computation
     169            renameDummyValuesF32("PLTSCALE", fittedTrails, "PLTSCALE_EXT", det->table, idet);
     170            renameDummyValuesF32("EXT_INST_FLUX", fittedTrails, "EXT_FLUX", det->table, idet);
     171            replaceDummyValuesF32("EXT_CAL_MAG", fittedTrails, det->table, idet);
     172            renameDummyValuesF32("EXT_INST_MAG_SIG", fittedTrails, "EXT_MAG_SIG", det->table, idet);
     173            replaceDummyValuesF32("EXT_CHISQ", fittedTrails, det->table, idet);
     174            replaceDummyValuesS32("EXT_NDOF", fittedTrails, det->table, idet);
     175            psTrace("ppMops.read", 10, "Got all version 3 values\n");
     176          }
     177        }
     178      }
     179    }
    153180    psFitsClose(fits);
    154181
    155182    if (args->version == 0) {
    156183      if (skyChipPsfVersion < 2) {
    157         // XXX: TODO: Do we need to add dummy vectors for the missing columns?
     184        // XXX: TODO: Do we need to add dummy vectors for the missing columns?
    158185      }
    159186    }
    … …  
    161188    psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
    162189    psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
    163 
    164190    det->raErr = psVectorAlloc(size, PS_TYPE_F64);
    165191    det->decErr = psVectorAlloc(size, PS_TYPE_F64);
     192    if (skyChipPsfVersion >= 3) {
     193      det->raExtErr = psVectorAlloc(size, PS_TYPE_F64);
     194      det->decExtErr = psVectorAlloc(size, PS_TYPE_F64);
     195    } else {
     196      det->raExtErr = NULL;
     197      det->decExtErr = NULL;
     198    }
    166199    det->mask = psVectorAlloc(size, PS_TYPE_U8);
    167200
    … …  
    177210
    178211    // Add our new vectors to the table so that duplicates and masked items may be purged
     212    psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
    179213    psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
    180     psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
     214    if (skyChipPsfVersion >= 3) {
     215      psMetadataAddVector(table, PS_LIST_HEAD, "RA_EXT_SIGMA", 0, NULL, det->raExtErr);
     216      psMetadataAddVector(table, PS_LIST_HEAD, "DEC_EXT_SIGMA", 0, NULL, det->decExtErr);
     217    }
    181218
    182219    psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
    … …  
    189226    psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
    190227    psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
     228    psVector *xExtErrV = NULL;
     229    psVector *yExtErrV = NULL;
     230    psVector *scaleExtV = NULL;
     231    psVector *angleExtV = NULL;
     232    if (skyChipPsfVersion >= 3) {
     233      xExtErrV = psMetadataLookupVector(NULL, table, "X_EXT_SIG");
     234      yExtErrV = psMetadataLookupVector(NULL, table, "Y_EXT_SIG");
     235      scaleExtV = psMetadataLookupVector(NULL, table, "PLTSCALE_EXT");
     236      angleExtV = psMetadataLookupVector(NULL, table, "POSANG_EXT");
     237    }
    191238
    192239    double plateScale = 0.0;        // Plate scale
    193240    long numGood = 0;               // Number of good rows
    194241    for (long row = 0; row < size; row++) {
    195 
    196242      psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
    197243      if (flags & SOURCE_MASK) {
    198         psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
    199         det->mask->data.U8[row] = 0xFF;
    200         continue;
     244        psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
     245        det->mask->data.U8[row] = 0xFF;
     246        continue;
    201247      }
    202248
    … …  
    206252      double scale = scaleV->data.F32[row];
    207253      double angle = angleV->data.F32[row];
    208 
    209254      if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
    210           !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
    211           !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
    212           !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
    213         psTrace("ppMops.read", 10,
    214                 "Discarding row %ld from input %d because of non-finite values: "
    215                 "%f %f %lf %lf %f %f %f %f %f %f",
    216                 row, i,
    217                 det->x->data.F32[row], det->y->data.F32[row],
    218                 det->ra->data.F64[row], det->dec->data.F64[row],
    219                 mag->data.F32[row], magErr->data.F32[row],
    220                 xErr, yErr, scale, angle);
    221         det->mask->data.U8[row] = 0xFF;
    222         continue;
     255          !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
     256          !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
     257          !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
     258        psTrace("ppMops.read", 10,
     259                "Discarding row %ld from input %d because of non-finite values: "
     260                "%f %f %lf %lf %f %f %f %f %f %f",
     261                row, i,
     262                det->x->data.F32[row], det->y->data.F32[row],
     263                det->ra->data.F64[row], det->dec->data.F64[row],
     264                mag->data.F32[row], magErr->data.F32[row],
     265                xErr, yErr, scale, angle);
     266        det->mask->data.U8[row] = 0xFF;
     267        continue;
    223268      }
    224269
    … …  
    232277      det->mask->data.U8[row] = 0;
    233278      plateScale += scale;
     279      //Update the platescale value (should be a constant)
     280      if (isfinite(scale)) {
     281          det->platescale = scale;
     282      }
    234283      numGood++;
    235     }
    236     det->seeing *= ((float) plateScale) / ((float) numGood);
    237 
    238     // Are we using numGood for anything outside of this function?
     284
     285      // Same for EXT data if version permits
     286      if (skyChipPsfVersion >= 3) {
     287        xErr = xExtErrV->data.F32[row];
     288        yErr = yExtErrV->data.F32[row];
     289        scale = scaleExtV->data.F32[row]/3600.;
     290        angle = angleExtV->data.F32[row];
     291        cosAngle = cos(angle);
     292        sinAngle = sin(angle);
     293        cosAngle2 = PS_SQR(cosAngle);
     294        sinAngle2 = PS_SQR(sinAngle);
     295        xErr2 = PS_SQR(xErr);
     296        yErr2 = PS_SQR(yErr);
     297        errScale = scale / 3600.0;
     298        det->raExtErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     299        det->decExtErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     300      }
     301    }
    239302    det->numGood = numGood;
    240303
    … …  
    257320}
    258321
    259 static void addDummyValues(psMetadata* md, long size) {
     322static psVector* createDummyS32(long size) {
     323  psVector* dummy = psVectorAlloc(size, PS_TYPE_S32);
     324  psVectorInit(dummy, 0);
     325  return dummy;
     326}
     327
     328static void addDummyValues(psMetadata* md, long size, int version) {
    260329  psMetadataAdd(md, PS_LIST_TAIL, "X_EXT", PS_DATA_VECTOR, "EXT model x coordinate", createDummyF32(size));
    261330  psMetadataAdd(md, PS_LIST_TAIL, "Y_EXT", PS_DATA_VECTOR, "EXT model y coordinate", createDummyF32(size));
    … …  
    264333  psMetadataAdd(md, PS_LIST_TAIL, "EXT_INST_MAG",  PS_DATA_VECTOR, "EXT fit instrumental magnitude", createDummyF32(size));
    265334  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));
     335  psMetadataAdd(md, PS_LIST_TAIL, "NPARAMS",  PS_DATA_VECTOR, "Number of model parameters", createDummyS32(size));
    267336  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MAJ",  PS_DATA_VECTOR, "EXT width (major axis), length for trail", createDummyF32(size));
    268337  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MIN",  PS_DATA_VECTOR, "EXT width (minor axis), sigma for trail", createDummyF32(size));
    … …  
    271340  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",  PS_DATA_VECTOR, "EXT width error (minor axis)", createDummyF32(size));
    272341  psMetadataAdd(md, PS_LIST_TAIL, "EXT_THETA_ERR",  PS_DATA_VECTOR, "EXT orientation angle (error)", createDummyF32(size));
     342  if (version >= 3) {
     343    psTrace("ppMops.read", 10, "Adding columns for version PS1_DV%d\n", version);
     344    psMetadataAdd(md, PS_LIST_TAIL, "RA_EXT",  PS_DATA_VECTOR, "Fitted centroid RA", createDummyF32(size));
     345/*     psMetadataAdd(md, PS_LIST_TAIL, "RA_EXT_SIGMA",  PS_DATA_VECTOR, "Fitted RA sigma", createDummyF32(size)); *\/ */
     346    psMetadataAdd(md, PS_LIST_TAIL, "DEC_EXT",  PS_DATA_VECTOR, "Fitted centroid DEC", createDummyF32(size));
     347/*     psMetadataAdd(md, PS_LIST_TAIL, "DEC_EXT_SIGMA",  PS_DATA_VECTOR, "Fitted DEC sigma", createDummyF32(size)); *\/ */
     348    psMetadataAdd(md, PS_LIST_TAIL, "POSANG_EXT",  PS_DATA_VECTOR, "Fitted position angle", createDummyF32(size));
     349    psMetadataAdd(md, PS_LIST_TAIL, "PLTSCALE_EXT",  PS_DATA_VECTOR, "Plate scale at centroid", createDummyF32(size));
     350    psMetadataAdd(md, PS_LIST_TAIL, "EXT_FLUX",  PS_DATA_VECTOR, "Fitted flux", createDummyF32(size));
     351    psMetadataAdd(md, PS_LIST_TAIL, "EXT_CAL_MAG",  PS_DATA_VECTOR, "Calibrated mag", createDummyF32(size));
     352    psMetadataAdd(md, PS_LIST_TAIL, "EXT_MAG_SIG",  PS_DATA_VECTOR, "Mag sigma", createDummyF32(size));
     353    psMetadataAdd(md, PS_LIST_TAIL, "EXT_CHISQ",  PS_DATA_VECTOR, "Chi^2 of fit", createDummyF32(size));
     354    psMetadataAdd(md, PS_LIST_TAIL, "EXT_NDOF",  PS_DATA_VECTOR, "Fit degrees of freedom", createDummyS32(size));
     355  }
    273356}
    274357
    275358static void replaceDummyValuesF32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes) {
     359  psTrace("ppMops.read", 10, "Trying to get column %s\n", colName);
    276360  psVector* source_vector = psMetadataLookupVector(NULL, source, colName);
    277361  psVector* target_vector = psMetadataLookupVector(NULL, target, colName);
    … …  
    279363    target_vector->data.F32[indexes->data.S64[index]] = source_vector->data.F32[index];
    280364  }
    281 }
     365  psTrace("ppMops.read", 10, "OK for %s\n", colName);
     366}
     367
     368static void replaceDummyValuesS32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes) {
     369  psTrace("ppMops.read", 10, "Trying to get column %s\n", colName);
     370  psVector* source_vector = psMetadataLookupVector(NULL, source, colName);
     371  psVector* target_vector = psMetadataLookupVector(NULL, target, colName);
     372  for (long index = 0; index<indexes->n; ++index) {
     373    target_vector->data.S32[indexes->data.S64[index]] = source_vector->data.S32[index];
     374  }
     375  psTrace("ppMops.read", 10, "OK for %s\n", colName);
     376}
     377
     378static void renameDummyValuesF32(const char* colNameSource, psMetadata* source, const char* colNameTarget, psMetadata* target, psVector* indexes) {
     379  psTrace("ppMops.read", 10, "Trying to rename column %s to %s\n", colNameSource, colNameTarget);
     380  psVector* source_vector = psMetadataLookupVector(NULL, source, colNameSource);
     381  psVector* target_vector = psMetadataLookupVector(NULL, target, colNameTarget);
     382  for (long index = 0; index<indexes->n; ++index) {
     383    target_vector->data.F32[indexes->data.S64[index]] = source_vector->data.F32[index];
     384  }
     385  psTrace("ppMops.read", 10, "OK for %s (renamed by %s)\n", colNameSource, colNameTarget);
     386}
Note: See TracChangeset for help on using the changeset viewer.