IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 30, 2010, 9:31:50 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r28243 r28794  
    44
    55#include <stdio.h>
    6 #include <string.h>
    76#include <pslib.h>
    87
    … …  
    1716    psArray *detections = psArrayAlloc(num); // Array of detections, to return
    1817    for (int i = 0; i < num; i++) {
    19         const char *name = inNames->data[i]; // File name
    20         psFits *fits = psFitsOpen(name, "r"); // FITS file
     18        psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file
    2119        if (!fits) {
    2220            psError(PS_ERR_IO, false, "Unable to open input %d", i);
    … …  
    2624        if (!header) {
    2725            psError(PS_ERR_IO, false, "Unable to read header %d", i);
     26            return false;
     27        }
     28
     29        psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
     30        if (diffSkyfileId == 0) {
     31            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
    2832            return false;
    2933        }
    … …  
    4347            continue;
    4448        }
    45         ppMopsDetections *det = detections->data[i] = ppMopsDetectionsAlloc();
    46         det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
    47         det->header = header;
    48         det->num = size;
     49        ppMopsDetections *det = ppMopsDetectionsAlloc(size);
     50
     51        psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
    4952
    5053        det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
    … …  
    6164                             psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    6265
    63         det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1");
    64         det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2");
     66        int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
     67        int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    6568
    66         det->psfHeader = psFitsReadHeader(NULL, fits);
    67         if (!det->psfHeader) {
    68             psError(psErrorCodeLast(), false, "Unable to read header %d", i);
     69        psFree(header);
     70
     71        psArray *table = psFitsReadTable(fits); // Table of interest
     72        if (!table) {
     73            psError(PS_ERR_IO, false, "Unable to read table %d", i);
    6974            return false;
    7075        }
    71         psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest
    72         if (!table) {
    73             psError(psErrorCodeLast(), false, "Unable to read table %d", i);
    74             return false;
     76        psFitsClose(fits);
     77
     78        double plateScale = 0.0;        // Plate scale
     79        long numGood = 0;               // Number of good rows
     80        for (long j = 0; j < size; j++) {
     81            psMetadata *row = table->data[j]; // Row of interest
     82
     83            psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS");
     84            if (flags & SOURCE_MASK) {
     85                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags);
     86                continue;
     87            }
     88
     89            det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF");
     90            det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF");
     91            det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF"));
     92            det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF"));
     93            det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG");
     94            det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG");
     95            det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ");
     96            det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF");
     97            det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA");
     98            det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA");
     99            det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR");
     100            det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR");
     101            det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA");
     102            det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF");
     103            det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX");
     104            det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX");
     105            det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY");
     106            det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY");
     107            det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS");
     108            det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;
     109            det->naxis1->data.S32[numGood] = naxis1;
     110            det->naxis2->data.S32[numGood] = naxis2;
     111
     112            det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS");
     113            det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO");
     114            det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD");
     115            det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK");
     116            det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL");
     117
     118            // Calculate error in RA, Dec
     119            double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG");
     120            double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG");
     121            double scale = psMetadataLookupF64(NULL, row, "PLTSCALE");
     122            double angle = psMetadataLookupF64(NULL, row, "POSANGLE");
     123
     124            if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||
     125                !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||
     126                !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) ||
     127                !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
     128                psTrace("ppMops.read", 10,
     129                        "Discarding row %ld from input %d because of non-finite values: "
     130                        "%f %f %lf %lf %f %f %f %f %f %f",
     131                        j, i,
     132                        det->x->data.F32[numGood], det->y->data.F32[numGood],
     133                        det->ra->data.F64[numGood], det->dec->data.F64[numGood],
     134                        det->mag->data.F32[numGood], det->magErr->data.F32[numGood],
     135                        xErr, yErr, scale, angle);
     136                continue;
     137            }
     138
     139            // XXX Not at all sure I've got the angles around the right way here...
     140            double cosAngle = cos(angle), sinAngle = sin(angle);
     141            double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
     142            double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
     143            double errScale = scale / 3600.0;
     144            det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     145            det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     146
     147            det->mask->data.U8[numGood] = 0;
     148            plateScale += scale;
     149            numGood++;
     150        }
     151        det->seeing *= plateScale / numGood;
     152
     153        det->x->n = numGood;
     154        det->y->n = numGood;
     155        det->ra->n = numGood;
     156        det->dec->n = numGood;
     157        det->raErr->n = numGood;
     158        det->decErr->n = numGood;
     159        det->mag->n = numGood;
     160        det->magErr->n = numGood;
     161        det->chi2->n = numGood;
     162        det->dof->n = numGood;
     163        det->cr->n = numGood;
     164        det->extended->n = numGood;
     165        det->psfMajor->n = numGood;
     166        det->psfMinor->n = numGood;
     167        det->psfTheta->n = numGood;
     168        det->quality->n = numGood;
     169        det->numPix->n = numGood;
     170        det->xxMoment->n = numGood;
     171        det->xyMoment->n = numGood;
     172        det->yyMoment->n = numGood;
     173        det->flags->n = numGood;
     174        det->diffSkyfileId->n = numGood;
     175        det->naxis1->n = numGood;
     176        det->naxis2->n = numGood;
     177        det->mask->n = numGood;
     178        det->nPos->n = numGood;
     179        det->fPos->n = numGood;
     180        det->nRatioBad->n = numGood;
     181        det->nRatioMask->n = numGood;
     182        det->nRatioAll->n = numGood;
     183
     184        det->num = numGood;
     185
     186        if (isfinite(args->zp) && numGood > 0) {
     187            psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
    75188        }
    76189
    77         psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
    78         psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
     190        psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]);
    79191
    80         det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    81         det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    82         det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
    83         det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
    84         if (!det->ra || !det->dec || !det->x || !det->y) {
    85             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
    86             return false;
    87         }
    88 
    89         psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
    90 
    91         if (!psFitsMoveExtName(fits, "SkyChip.deteff")) {
    92             psWarning("No detection efficiencies included in %s", det->component);
    93             psErrorStackPrint(stderr, "No detection efficiencies included in %s", det->component);
    94             psErrorClear();
    95             continue;
    96         }
    97 
    98         det->deteffHeader = psFitsReadHeader(NULL, fits);
    99         if (!det->deteffHeader) {
    100             psWarning("Unable to read detection efficiency header in %s", det->component);
    101             psErrorClear();
    102             continue;
    103         }
    104         det->deteffTable = psFitsReadTableAllColumns(fits);
    105         if (!det->deteffTable) {
    106             psFree(det->deteffHeader);
    107             det->deteffHeader = NULL;
    108             psWarning("Unable to read detection efficiency table in %s", det->component);
    109             psErrorClear();
    110             continue;
    111         }
    112         psTrace("ppMops.read", 2, "Read detection efficiency from %s\n", det->component);
    113         psFitsClose(fits);
     192        psFree(table);
     193        detections->data[i] = det;
    114194    }
    115195
Note: See TracChangeset for help on using the changeset viewer.