Index: /tags/ipp-20120802/ppTranslate/src/ppMops.c =================================================================== --- /tags/ipp-20120802/ppTranslate/src/ppMops.c (revision 34280) +++ /tags/ipp-20120802/ppTranslate/src/ppMops.c (revision 34281) @@ -67,4 +67,5 @@ } + psArray *detections = ppMopsRead(args); // Detections from each input if (!detections) { @@ -91,9 +92,202 @@ psLibFinalize(); -/* fprintf (stderr, "found %d leaks at %s\n", */ -/* psMemCheckLeaks2 (0, */ -/* NULL, stdout, false, 500), "ppMops"); */ + fprintf (stderr, "found %d leaks at %s\n", + psMemCheckLeaks2 (0, + NULL, stdout, false, 500), "ppMops"); return PS_EXIT_SUCCESS; } + +#if 0 + ps + + + psArray *detections = NULL; // Detections + psMetadata *header = NULL; // Header for detections + { + psFits *fits = psFitsOpen(data->detections, "r"); // FITS file + header = psFitsReadHeader(NULL, fits); + if (!header) { + psErrorStackPrint(stderr, "Unable to read header"); + psFitsClose(fits); + psFree(data); + exit(PS_EXIT_DATA_ERROR); + } + if (!psFitsMoveExtName(fits, IN_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to move to extension %s", IN_EXTNAME); + psFitsClose(fits); + psFree(header); + psFree(data); + exit(PS_EXIT_DATA_ERROR); + } + detections = psFitsReadTable(fits); + psFitsClose(fits); + if (!detections) { + psErrorStackPrint(stderr, "Unable to read detections"); + psFree(data); + psFree(header); + exit(PS_EXIT_DATA_ERROR); + } + } + + // Translate the columns + int numIn = detections->n; // Number of rows in input + int numOut = 0; // Number of rows in output + psArray *output = psArrayAllocEmpty(numIn); // Output table + double plateScale = 0.0; // Average plate scale + for (int i = 0; i < numIn; i++) { + psMetadata *inRow = detections->data[i]; // Input row + + double ra = psMetadataLookupF64(NULL, inRow, "RA_PSF"); + double dec = psMetadataLookupF64(NULL, inRow, "DEC_PSF"); + double mag = psMetadataLookupF64(NULL, inRow, "PSF_INST_MAG") + data->zp; + double magErr = psMetadataLookupF64(NULL, inRow, "PSF_INST_MAG_SIG"); + float ext = psMetadataLookupF32(NULL, inRow, "EXT_NSIGMA"); + double xErr = psMetadataLookupF64(NULL, inRow, "X_PSF_SIG"); + double yErr = psMetadataLookupF64(NULL, inRow, "Y_PSF_SIG"); + double scale = psMetadataLookupF64(NULL, inRow, "PLTSCALE"); + double angle = psMetadataLookupF64(NULL, inRow, "POSANGLE"); + psU32 flags = psMetadataLookupU32(NULL, inRow, "FLAGS"); + + if (!isfinite(mag) || !isfinite(magErr) || (flags & SOURCE_MASK) || + ((flags & PM_SOURCE_MODE_DEFECT) && !(flags & PM_SOURCE_MODE_MOMENTS_FAILURE)) + ) { + // DEFECT can be due to bad moments + continue; + } + + psMetadata *outRow = output->data[numOut] = psMetadataAlloc(); // Output row + + numOut++; + plateScale += scale; + + // XXX Not at all sure I've got the angles around the right way here... + double cosAngle = cos(angle), sinAngle = sin(angle); + double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle); + double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); + double errScale = scale / 3600.0; + double raErr = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + double decErr = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_DEG", 0, "Right ascension (degrees)", ra); + psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_DEG", 0, "Declination (degrees)", dec); + psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_SIG", 0, "Right ascension error (degrees)", raErr); + psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_SIG", 0, "Declination error (degrees)", decErr); + psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG", 0, "Magnitude", mag); + psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG_SIG", 0, "Magnitude error", magErr); + psMetadataAddU32(outRow, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", flags); + + // The below need fixing + psMetadataAddF64(outRow, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", ext); + psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG", 0, "Position angle of trail (degrees)", 0.0); + psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG_SIG", 0, "Position angle error (degrees)", 0.0); + psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN", 0, "Length of trail (arcsec)", 0.0); + psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN_SIG", 0, "Length error (arcsec)", 0.0); + } + output->n = numOut; + plateScale /= numOut; + psFree(detections); + + // Translate the header + psMetadata *outHeader = psMetadataAlloc(); // Output header + ppMopsVersionHeader(outHeader); + { + const char *ra = psMetadataLookupStr(NULL, header, "FPA.RA"); + const char *dec = psMetadataLookupStr(NULL, header, "FPA.DEC"); + const char *filter = psMetadataLookupStr(NULL, header, "FPA.FILTER"); + float airmass = psMetadataLookupF32(NULL, header, "FPA.AIRMASS"); + float exptime = psMetadataLookupF32(NULL, header, "EXPTIME"); + double angle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE"); + double alt = psMetadataLookupF64(NULL, header, "FPA.ALT"); + double az = psMetadataLookupF64(NULL, header, "FPA.AZ"); + psS64 imageid = psMetadataLookupS64(NULL, header, "IMAGEID"); + double mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + exptime / 2.0 / 3600 / 24; + + float psf = plateScale * 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") + + psMetadataLookupF32(NULL, header, "FWHM_MIN")); + + psMetadataAddStr(outHeader, PS_LIST_TAIL, "RA", 0, "Right ascension of boresight", ra); + psMetadataAddStr(outHeader, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", dec); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", airmass); + psMetadataAddF64(outHeader, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", mjd); + psMetadataAddStr(outHeader, PS_LIST_TAIL, "FILTER", 0, "Filter name", filter); + psMetadataAddF64(outHeader, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", exptime); + psMetadataAddF64(outHeader, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", angle); + psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", alt); + psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", az); + psMetadataAddS64(outHeader, PS_LIST_TAIL, "DIFFIMID", 0, "Difference image identifier", imageid); + psMetadataAddStr(outHeader, PS_LIST_TAIL, "FPA_ID", 0, "Exposure name", data->exp_name); + psMetadataAddS64(outHeader, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", data->exp_id); + psMetadataAddBool(outHeader, PS_LIST_TAIL, "POSITIVE", 0, "Positive subtraction?", data->direction); + psMetadataAddStr(outHeader, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "STARPSF", 0, "Stellar PSF (arcsec)", psf); + + // These are completely fake + psMetadataAddF32(outHeader, PS_LIST_TAIL, "LIMITMAG", 0, "Limiting magnitude (FAKE)", 99.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE1", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE2", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE3", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE4", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE5", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE6", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE7", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE8", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE9", 0, "Detection efficiency (FAKE)", 0.0); + psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE10", 0, "Detection efficiency (FAKE)", 0.0); + } + psFree(header); + + // Write the new table + { + psFits *fits = psFitsOpen(data->output, "w"); // FITS file + if (!fits) { + psErrorStackPrint(stderr, "Unable to open %s for writing", data->output); + psFree(outHeader); + psFree(output); + psFree(data); + exit(PS_EXIT_SYS_ERROR); + } + + if (numOut == 0) { + // Write dummy table + psMetadata *outRow = psMetadataAlloc(); // Dummy output row + psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_DEG", 0, "Right ascension (degrees)", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_DEG", 0, "Declination (degrees)", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_SIG", 0, "Right ascension error (degrees)", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_SIG", 0, "Declination error (degrees)", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG_SIG", 0, "Magnitude error", NAN); + psMetadataAddU32(outRow, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0); + + // The below need fixing + psMetadataAddF64(outRow, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG", 0, "Position angle of trail (degrees)", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG_SIG", 0, "Position angle error (degrees)", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN", 0, "Length of trail (arcsec)", NAN); + psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN_SIG", 0, "Length error (arcsec)", NAN); + if (!psFitsWriteTableEmpty(fits, outHeader, outRow, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write table."); + psFree(outHeader); + psFree(output); + psFree(outRow); + psFree(data); + exit(PS_EXIT_SYS_ERROR); + } + psFree(outRow); + } else if (!psFitsWriteTable(fits, outHeader, output, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write table."); + psFree(outHeader); + psFree(output); + psFree(data); + exit(PS_EXIT_SYS_ERROR); + } + psFitsClose(fits); + } + + psFree(outHeader); + psFree(output); + psFree(data); + +#endif + Index: /tags/ipp-20120802/ppTranslate/src/ppMops.h =================================================================== --- /tags/ipp-20120802/ppTranslate/src/ppMops.h (revision 34280) +++ /tags/ipp-20120802/ppTranslate/src/ppMops.h (revision 34281) @@ -30,4 +30,32 @@ } ppMopsArguments; +#if 0 +#warning "IS THERE ANYTHING TO BE MODIFIED HERE?" +TTYPE19 = 'PSF_CHISQ' / label for field 19 +TFORM19 = '1E ' / data format of field: 4-byte REAL +TTYPE20 = 'CR_NSIGMA' / label for field 20 +TFORM20 = '1E ' / data format of field: 4-byte REAL +TTYPE21 = 'EXT_NSIGMA' / label for field 21 +TFORM21 = '1E ' / data format of field: 4-byte REAL +TTYPE22 = 'PSF_MAJOR' / label for field 22 +TFORM22 = '1E ' / data format of field: 4-byte REAL +TTYPE23 = 'PSF_MINOR' / label for field 23 +TFORM23 = '1E ' / data format of field: 4-byte REAL +TTYPE24 = 'PSF_THETA' / label for field 24 +TFORM24 = '1E ' / data format of field: 4-byte REAL +TTYPE25 = 'PSF_QF ' / label for field 25 +TFORM25 = '1E ' / data format of field: 4-byte REAL +TTYPE26 = 'PSF_NDOF' / label for field 26 +TFORM26 = '1J ' / data format of field: 4-byte INTEGER +TTYPE27 = 'PSF_NPIX' / label for field 27 +TFORM27 = '1J ' / data format of field: 4-byte INTEGER +TTYPE28 = 'MOMENTS_XX' / label for field 28 +TFORM28 = '1E ' / data format of field: 4-byte REAL +TTYPE29 = 'MOMENTS_XY' / label for field 29 +TFORM29 = '1E ' / data format of field: 4-byte REAL +TTYPE30 = 'MOMENTS_YY' / label for field 30 +TFORM30 = '1E ' / data format of field: 4-byte REAL +#endif + /// Parse arguments ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]); @@ -47,7 +75,5 @@ long numGood; // Number of "good" detections psS64 diffSkyfileId; // unique id for input skyfile - psMetadata *table; // Columns from the input file (SkyChip.psf extension) - psMetadata *fittedTrails; // Columns from the fitted trails (SkyChip.xfit extension) - long fittedTrailsSize; // Size of the fitted trails (SkyChip.xfit extension) + psMetadata *table; // Columns from the input file psVector *x, *y; // Image coordinates psVector *ra, *dec; // Sky coordinates @@ -77,5 +103,4 @@ /// @returns 1 if EXTTYPE of "SkyChip.psf" is PS1_DV1 /// @returns 2 if EXTTYPE of "SkyChip.psf" is PS1_DV2 -/// @returns 3 if EXTTYPE of "SkyChip.psf" is PS1_DV3 /// @returns 0 otherwise int ppMopsGetSkyChipPsfVersion(const psFits* fits); Index: /tags/ipp-20120802/ppTranslate/src/ppMopsArguments.c =================================================================== --- /tags/ipp-20120802/ppTranslate/src/ppMopsArguments.c (revision 34280) +++ /tags/ipp-20120802/ppTranslate/src/ppMopsArguments.c (revision 34281) @@ -34,5 +34,5 @@ fprintf(stderr, "\n"); psArgumentHelp(arguments); - fprintf(stderr, "\t\tCMF file version can be set to either 1 for PS1_DV1, 2 for PS1_DV2, or 3 for PS1_DV3\n"); + fprintf(stderr, "\t\tCMF file version can be set to either 1 for PS1_DV1 or 2 for PS1_DV1\n"); fprintf(stderr, "\t\tSee IPP-MOPS ICD for details\n"); psLibFinalize(); Index: /tags/ipp-20120802/ppTranslate/src/ppMopsGetSkyChipPsfVersion.c =================================================================== --- /tags/ipp-20120802/ppTranslate/src/ppMopsGetSkyChipPsfVersion.c (revision 34280) +++ /tags/ipp-20120802/ppTranslate/src/ppMopsGetSkyChipPsfVersion.c (revision 34281) @@ -14,7 +14,4 @@ psFree(headerSkyChip); return 2; - } else if (strcmp(version, "PS1_DV3") == 0) { - psFree(headerSkyChip); - return 2; } psWarning("Unsupported EXTTYPE in SkyChip.psf table: [%s]", version); Index: /tags/ipp-20120802/ppTranslate/src/ppMopsMerge.c =================================================================== --- /tags/ipp-20120802/ppTranslate/src/ppMopsMerge.c (revision 34280) +++ /tags/ipp-20120802/ppTranslate/src/ppMopsMerge.c (revision 34281) @@ -71,5 +71,4 @@ for (int i = 0; i < numInputs; i++) { - psTrace("ppMops.merge", 1, "Processing batch %d\n", i); ppMopsDetections *det = detections->data[i]; // Detections of interest if (!det) { Index: /tags/ipp-20120802/ppTranslate/src/ppMopsRead.c =================================================================== --- /tags/ipp-20120802/ppTranslate/src/ppMopsRead.c (revision 34280) +++ /tags/ipp-20120802/ppTranslate/src/ppMopsRead.c (revision 34281) @@ -8,346 +8,201 @@ #include "ppMops.h" -static psMetadata* sort_by_increasing_idet(psMetadata *fittedTrail, long length, long fillLength); - /* ppMopsRead possibly modifies the args->version if the user did not set it explicitely. -*/ + */ psArray *ppMopsRead(ppMopsArguments *args) { - psTrace("ppMops.read", 1, "Reading input detections\n"); - - psArray *inNames = args->input; // Input names - long num = inNames->n; // Number of inputs - psArray *detections = psArrayAlloc(num); // Array of detections, to return - for (int i = 0; i < num; i++) { - const char *name = inNames->data[i]; - - psFits *fits = psFitsOpen(name, "r"); // FITS file - if (!fits) { - psError(PS_ERR_IO, false, "Unable to open input %d", i); - return false; + psTrace("ppMops.read", 1, "Reading input detections\n"); + + psArray *inNames = args->input; // Input names + long num = inNames->n; // Number of inputs + psArray *detections = psArrayAlloc(num); // Array of detections, to return + for (int i = 0; i < num; i++) { + const char *name = inNames->data[i]; + + psFits *fits = psFitsOpen(name, "r"); // FITS file + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open input %d", i); + return false; + } + + psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header %d", i); + return false; + } + + psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image + if (diffSkyfileId == 0) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i); + return false; + } + + if (!psFitsMoveExtName(fits, "SkyChip.psf")) { + psError(PS_ERR_IO, false, "Unable to move to HDU with detections"); + return false; + } + int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits); + if (args->version == 0) { + psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion); + args->version = (unsigned short) skyChipPsfVersion; + } + if (skyChipPsfVersion == 0) { + // Try to read with the user specified version? + skyChipPsfVersion = args->version; + } + /* Display a warning message if there are version + inconsistencies between the file and the flag (note that + those inconsistencies might be wanted) */ + if (skyChipPsfVersion != args->version) { + if (skyChipPsfVersion > args->version) { + psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n", + skyChipPsfVersion, args->version); + } else { // Necessarily: skyChipPsfVersion > args->version + psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n", + skyChipPsfVersion, args->version); + } + } + + long size = psFitsTableSize(fits); // Size of table + if (size <= 0) { + psErrorStackPrint(stderr, "Unable to determine size of table %d", i); + psErrorClear(); + psWarning("Ignoring input %d", i); + psFree(header); + psFitsClose(fits); + continue; + } + ppMopsDetections *det = ppMopsDetectionsAlloc(size); + detections->data[i] = det; + det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension + det->num = size; + det->diffSkyfileId = diffSkyfileId; + + psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); + + det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA")); + det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC")); + det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER")); + det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS"); + det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME"); + det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE"); + det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT"); + det->az = psMetadataLookupF64(NULL, header, "FPA.AZ"); + det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24; + + det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") + + psMetadataLookupF32(NULL, header, "FWHM_MIN")); + + det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns + det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows + + psFree(header); + + psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest + if (!table) { + psError(PS_ERR_IO, false, "Unable to read table %d", i); + return NULL; + } + det->table = table; + psFitsClose(fits); + if (args->version == 0) { + if (skyChipPsfVersion < 2) { + // XXX: TODO: Do we need to add dummy vectors for the missing columns? + } + } + + psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); + psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); + + det->raErr = psVectorAlloc(size, PS_TYPE_F64); + det->decErr = psVectorAlloc(size, PS_TYPE_F64); + det->mask = psVectorAlloc(size, PS_TYPE_U8); + + // convert ra and dec to radians for use in the purge duplicates function + det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); + det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); + det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF")); + det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF")); + if (!det->ra || !det->dec || !det->x || !det->y) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns"); + return NULL; + } + + // Add our new vectors to the table so that duplicates and masked items may be purged + psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr); + psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr); + + psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); + + psVector *mag = psMetadataLookupVector(NULL, table, "PSF_INST_MAG"); + psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG"); + psVector *xErrV = psMetadataLookupVector(NULL, table, "X_PSF_SIG"); + psVector *yErrV = psMetadataLookupVector(NULL, table, "Y_PSF_SIG"); + psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE"); + psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE"); + psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS"); + + double plateScale = 0.0; // Plate scale + long numGood = 0; // Number of good rows + for (long row = 0; row < size; row++) { + + psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS"); + if (flags & SOURCE_MASK) { + psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags); + det->mask->data.U8[row] = 0xFF; + continue; + } + + // Calculate error in RA, Dec + + double xErr = xErrV->data.F32[row]; + double yErr = yErrV->data.F32[row]; + double scale = scaleV->data.F32[row]; + double angle = angleV->data.F32[row]; + + if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) || + !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) || + !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) || + !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) { + psTrace("ppMops.read", 10, + "Discarding row %ld from input %d because of non-finite values: " + "%f %f %lf %lf %f %f %f %f %f %f", + row, i, + det->x->data.F32[row], det->y->data.F32[row], + det->ra->data.F64[row], det->dec->data.F64[row], + mag->data.F32[row], magErr->data.F32[row], + xErr, yErr, scale, angle); + det->mask->data.U8[row] = 0xFF; + continue; + } + + // XXX Not at all sure I've got the angles around the right way here... + double cosAngle = cos(angle), sinAngle = sin(angle); + double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle); + double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); + double errScale = scale / 3600.0; + det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + det->mask->data.U8[row] = 0; + plateScale += scale; + numGood++; + } + det->seeing *= ((float) plateScale) / ((float) numGood); + + // Are we using numGood for anything outside of this function? + det->numGood = numGood; + + if (isfinite(args->zp) && numGood > 0) { + psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); + } + + psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name); } - psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header - if (!header) { - psError(PS_ERR_IO, false, "Unable to read header %d", i); - return false; - } - - psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image - if (diffSkyfileId == 0) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i); - return false; - } - - if (!psFitsMoveExtName(fits, "SkyChip.psf")) { - psError(PS_ERR_IO, false, "Unable to move to HDU with detections"); - return false; - } - int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits); - if (args->version == 0) { - psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion); - args->version = (unsigned short) skyChipPsfVersion; - } - if (skyChipPsfVersion == 0) { - // Try to read with the user specified version? - skyChipPsfVersion = args->version; - } - /* Display a warning message if there are version - inconsistencies between the file and the flag (note that - those inconsistencies might be wanted) */ - if (skyChipPsfVersion != args->version) { - if (skyChipPsfVersion > args->version) { - psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n", - skyChipPsfVersion, args->version); - } else { // Necessarily: skyChipPsfVersion > args->version - psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n", - skyChipPsfVersion, args->version); - } - } - - long size = psFitsTableSize(fits); // Size of table - if (size <= 0) { - psErrorStackPrint(stderr, "Unable to determine size of table %d", i); - psErrorClear(); - psWarning("Ignoring input %d", i); - psFree(header); - psFitsClose(fits); - continue; - } - ppMopsDetections *det = ppMopsDetectionsAlloc(size); - detections->data[i] = det; - det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension - det->num = size; - det->diffSkyfileId = diffSkyfileId; - - psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); - - det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA")); - det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC")); - det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER")); - det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS"); - det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME"); - det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE"); - det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT"); - det->az = psMetadataLookupF64(NULL, header, "FPA.AZ"); - det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24; - - det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") + - psMetadataLookupF32(NULL, header, "FWHM_MIN")); - - det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns - det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows - - psFree(header); - - psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest - if (!table) { - psError(PS_ERR_IO, false, "Unable to read table %d", i); - return NULL; - } - //Sort table by idet if necessary - psMetadata* sortedTable = sort_by_increasing_idet(table, det->num, det->num); - if (sortedTable != table) { - psFree(table); - table = sortedTable; - } // else det->table is sorted - det->table = table; - psTrace("ppMops.read", 9, "PSF table:\n%s\n", psMetadataConfigFormat(det->table)); - psTrace("ppMops.read", 9, "End of PSF table\n"); - - //fittedTrail extension (Supported for releases 2 and above) - if (skyChipPsfVersion >= 2) { - if (!psFitsMoveExtName(fits, "SkyChip.xfit")) { - psTrace("ppMops.read", 3, "No fitted trails extension"); - } else { - psTrace("ppMops.read", 3, "Fitted trails extension found\n"); - det->fittedTrailsSize = psFitsTableSize(fits); - if (det->fittedTrailsSize <= 0) { - psErrorStackPrint(stderr, "Unable to determine size of fitted trails extension table %d", i); - psTrace("ppMops.read", 3, "No entry in fitted trails extension!!!!\n"); - psErrorClear(); - det->fittedTrails = NULL; - } else { - det->fittedTrails = psFitsReadTableAllColumns(fits); // Table of interest - if (!det->fittedTrails) { - psError(PS_ERR_IO, false, "Unable to read fittedTrails table %d", i); - return NULL; - } - psMetadata* sortedFittedTrails = sort_by_increasing_idet(det->fittedTrails, det->fittedTrailsSize, det->num); - if (sortedFittedTrails != det->fittedTrails) { - psFree(det->fittedTrails); - det->fittedTrails = sortedFittedTrails; - } // else det->fittedTrails is sorted - psTrace("ppMops.read", 9, "Fitted trail:\n%s\n", psMetadataConfigFormat(det->fittedTrails)); - psTrace("ppMops.read", 9, "End of Fitted trail\n"); - } - } - } - - psFitsClose(fits); - - if (args->version == 0) { - if (skyChipPsfVersion < 2) { - // XXX: TODO: Do we need to add dummy vectors for the missing columns? - } - } - - psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); - psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); - - det->raErr = psVectorAlloc(size, PS_TYPE_F64); - det->decErr = psVectorAlloc(size, PS_TYPE_F64); - det->mask = psVectorAlloc(size, PS_TYPE_U8); - - // convert ra and dec to radians for use in the purge duplicates function - det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); - det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); - det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF")); - det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF")); - if (!det->ra || !det->dec || !det->x || !det->y) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns"); - return NULL; - } - - // Add our new vectors to the table so that duplicates and masked items may be purged - psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr); - psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr); - - psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); - - psVector *mag = psMetadataLookupVector(NULL, table, "PSF_INST_MAG"); - psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG"); - psVector *xErrV = psMetadataLookupVector(NULL, table, "X_PSF_SIG"); - psVector *yErrV = psMetadataLookupVector(NULL, table, "Y_PSF_SIG"); - psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE"); - psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE"); - psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS"); - - double plateScale = 0.0; // Plate scale - long numGood = 0; // Number of good rows - for (long row = 0; row < size; row++) { - - psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS"); - if (flags & SOURCE_MASK) { - psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags); - det->mask->data.U8[row] = 0xFF; - continue; - } - - // Calculate error in RA, Dec - - double xErr = xErrV->data.F32[row]; - double yErr = yErrV->data.F32[row]; - double scale = scaleV->data.F32[row]; - double angle = angleV->data.F32[row]; - - if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) || - !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) || - !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) || - !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) { - psTrace("ppMops.read", 10, - "Discarding row %ld from input %d because of non-finite values: " - "%f %f %lf %lf %f %f %f %f %f %f", - row, i, - det->x->data.F32[row], det->y->data.F32[row], - det->ra->data.F64[row], det->dec->data.F64[row], - mag->data.F32[row], magErr->data.F32[row], - xErr, yErr, scale, angle); - det->mask->data.U8[row] = 0xFF; - continue; - } - - // XXX Not at all sure I've got the angles around the right way here... - double cosAngle = cos(angle), sinAngle = sin(angle); - double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle); - double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); - double errScale = scale / 3600.0; - det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); - det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); - det->mask->data.U8[row] = 0; - plateScale += scale; - numGood++; - } - det->seeing *= ((float) plateScale) / ((float) numGood); - - // Are we using numGood for anything outside of this function? - det->numGood = numGood; - - if (isfinite(args->zp) && numGood > 0) { - psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); - } - - psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name); - } - - psTrace("ppMops.read", 1, "Done reading input detections\n"); - - return detections; + psTrace("ppMops.read", 1, "Done reading input detections\n"); + + return detections; } - -static long* sorted_idets_indexes(psVector *idets); -static int compare (const void *a, const void *b); -static psVector* fill_F32_vector(psVector* in, int size, long* indexes, int indexes_length); - -// We want the entries in table/fittedTrails to be sorted by increasing idet -static psMetadata* sort_by_increasing_idet(psMetadata *aTable, long length, long fillLength) { - if (aTable == NULL) { - return aTable; - } - psVector *idets = psMetadataLookupVector(NULL, aTable, "IPP_IDET"); - if (!idets) { - psError(PS_ERR_PROGRAMMING, true, "Failed to find IPP_IDET column"); - return NULL; - } - bool is_sorted = true; - int count = 0; - long previous_value = -1; - while (is_sorted && (countdata.S64[count]); - previous_value = idets->data.S64[count]; - count += 1; - } - if (is_sorted && (length==fillLength)) { - psTrace("ppMops.read", 5, "Values are already sorted by idet and the table is complete\n"); - return aTable; - } - //otherwise sort (even if it should not happen) and fill - long* sortedIndexVector = sorted_idets_indexes(idets); - psMetadata* sortedTable = psMetadataAlloc(); - char* column_name = "X_EXT"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "Y_EXT"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "X_EXT_SIG"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "Y_EXT_SIG"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_INST_MAG"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_INST_MAG_SIG"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "NPARAMS"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_WIDTH_MAJ"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_WIDTH_MIN"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_THETA"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_WIDTH_MAJ_ERR"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_WIDTH_MIN_ERR"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - column_name = "EXT_THETA_ERR"; - psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, - fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) ); - return sortedTable; -} - -static psVector* fill_F32_vector(psVector* in, int length, long* indexes, int indexes_length) { - if (in == NULL) { - return NULL; - } - psVector* out = psVectorAlloc(length, in->type.type); - int current_index = 0; - for (int i = 0; i < indexes_length; i++) { - while ( (current_index < indexes[i]) && (current_index < length) ) { - out->data.F32[current_index] = NAN; - current_index += 1; - } - out->data.F32[current_index] = in->data.F32[i]; - current_index += 1; - } - while ( current_index < length ) { - out->data.F32[current_index] = NAN; - current_index += 1; - } - return out; -} - -static int compare (const void *a, const void *b) { - return ( *((long*) a)-*((long*) b)); -} - -static long* sorted_idets_indexes(psVector *idets) { - long* sortedIndexes = (long*) malloc(idets->n * sizeof(long)); - for (int i=0;in;i++) { - sortedIndexes[i] = idets->data.S64[i]; - } - qsort(sortedIndexes, idets->n, sizeof(long), (void *) compare); - return sortedIndexes; -} Index: /tags/ipp-20120802/ppTranslate/src/ppMopsWrite.c =================================================================== --- /tags/ipp-20120802/ppTranslate/src/ppMopsWrite.c (revision 34280) +++ /tags/ipp-20120802/ppTranslate/src/ppMopsWrite.c (revision 34281) @@ -9,419 +9,378 @@ #include "ppTranslateVersion.h" -static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, int extension, char *outColName, char *inColName, bool convertTo32); +static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, char *outColName, char *inColName, bool convertTo32); static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName); bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args) { - psFits *fits = psFitsOpen(args->output, "w"); // FITS file - if (!fits) { - psError(PS_ERR_IO, false, "Unable to open output file."); - return false; - } - - psMetadata *header = psMetadataAlloc(); // Header to write - psString source = ppTranslateSource(), version = ppTranslateVersion(); - psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source); - psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version); - ppTranslateVersionHeader(header); - psFree(source); - psFree(version); - - psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name); - psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id); - psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id); - psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id); - psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id); - psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id); - psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id); - psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive); - - // Get these header words from the first input non null input - ppMopsDetections *det = NULL; - for (int d = 0; d < detections->n && det == NULL; d++) { - det = detections->data[d]; - } - if (det != NULL) { - psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd); - psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight); - psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight); - psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt); - psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az); - psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime); - psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle); - psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter); - psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass); - psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing); - } else { - psWarning("no inputs with surviving detections. output header will be incomplete"); - } - psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); - psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp); - psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr); - psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom); - - //field in header that tells about the CMF version - char cmfVersion[8]; - sprintf(cmfVersion, PS1_DV_FORMAT, args->version); - psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion); - - // Find the total number of detections - - long total = 0; - for (long i=0; in; i++) { - ppMopsDetections *det = detections->data[i]; - if (!det) { - continue; - } - total += det->num; - } - - psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output); - - if (total == 0) { - // Write dummy table - psMetadata *row = psMetadataAlloc(); // Output row - psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN); - psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN); - psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN); - psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN); - psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0); - psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN); - psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0); - psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN); - psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0); - psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0); - - psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0); - psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN); - psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN); - - if (args->version == 2) { - // Write data of version 2 (see ICD) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", - NAN); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags (group 2)", - 0); - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", - 0); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", - 0); - if (args->version == 3) { - // TODO: Write data of version 3 (see ICD) - } - } - - if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) { - psErrorStackPrint(stderr, "Unable to write empty table."); - psFree(header); - psFree(row); - return false; - } - psFree(row); - } else { - -#define addColumn(_outName, _inName, _convertTo32) \ - if (!addOutputColumn(table, detections, total, 0, _outName, _inName, _convertTo32)) { \ - psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \ - return false; \ - } - - // Allocate the output table - psMetadata *table = psMetadataAlloc(); - addColumn("RA", "RA_PSF", 0); - addColumn("RA_ERR", NULL, 0); // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG - addColumn("DEC", "DEC_PSF", 0); - addColumn("DEC_ERR", NULL, 0); // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG - addColumn("MAG", "PSF_INST_MAG", 0); - addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0); - addColumn("PSF_CHI2", "PSF_CHISQ", 0); - addColumn("PSF_DOF", "PSF_NDOF", 1); - addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0); - addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0); - addColumn("PSF_MAJOR", NULL, 0); - addColumn("PSF_MINOR", NULL, 0); - addColumn("PSF_THETA", NULL, 0); - addColumn("PSF_QUALITY", "PSF_QF", 0); - addColumn("PSF_NPIX", NULL, 1); - addColumn("MOMENTS_XX", NULL, 0); - addColumn("MOMENTS_XY", NULL, 0); - addColumn("MOMENTS_YY", NULL, 0); - addColumn("N_POS", "DIFF_NPOS", 1); - addColumn("F_POS", "DIFF_FRATIO", 0); - addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0); - addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0); - addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0); - addColumn("FLAGS", "FLAGS", 1); - addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID"); - if (args->version >= 2) { - addColumn("IPP_IDET", NULL, 1); - addColumn("PSF_INST_FLUX", NULL, 0); - addColumn("PSF_INST_FLUX_SIG", NULL, 0); - addColumn("AP_MAG", NULL, 0); - addColumn("AP_MAG_RAW", NULL, 0); - addColumn("AP_MAG_RADIUS", NULL, 0); - addColumn("AP_FLUX", NULL, 0); - addColumn("AP_FLUX_SIG", NULL, 0); - addColumn("PEAK_FLUX_AS_MAG", NULL, 0); - addColumn("CAL_PSF_MAG", NULL, 0); - addColumn("CAL_PSF_MAG_SIG", NULL, 0); - addColumn("SKY", NULL, 0); - addColumn("SKY_SIGMA", NULL, 0); - addColumn("PSF_QF_PERFECT", NULL, 0); - addColumn("MOMENTS_R1", NULL, 0); - addColumn("MOMENTS_RH", NULL, 0); - addColumn("KRON_FLUX", NULL, 0); - addColumn("KRON_FLUX_ERR", NULL, 0); - addColumn("KRON_FLUX_INNER", NULL, 0); - addColumn("KRON_FLUX_OUTER", NULL, 0); - addColumn("DIFF_R_P", NULL, 0); - addColumn("DIFF_SN_P", NULL, 0); - addColumn("DIFF_R_M", NULL, 0); - addColumn("DIFF_SN_M", NULL, 0); - addColumn("FLAGS2", NULL, 1); - addColumn("IPP_IDET", NULL, 0); - addColumn("N_FRAMES", NULL, 0); - addColumn("PADDING", NULL, 0); - if (det->fittedTrails != NULL) { -#define addFittedTrailColumn(_outName, _inName, _convertTo32) \ - if (!addOutputColumn(table, detections, total, 1, _outName, _inName, _convertTo32)) { \ - psTrace("ppMops.write", 1, "Failed to add column %s -> Replaced with NAN", _outName); \ + + psFits *fits = psFitsOpen(args->output, "w"); // FITS file + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open output file."); + return false; + } + + psMetadata *header = psMetadataAlloc(); // Header to write + psString source = ppTranslateSource(), version = ppTranslateVersion(); + psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source); + psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version); + ppTranslateVersionHeader(header); + psFree(source); + psFree(version); + + psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name); + psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id); + psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id); + psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id); + psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id); + psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id); + psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id); + psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive); + + // Get these header words from the first input non null input + ppMopsDetections *det = NULL; + for (int d = 0; d < detections->n && det == NULL; d++) { + det = detections->data[d]; + } + if (det != NULL) { + psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd); + psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight); + psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight); + psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt); + psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az); + psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime); + psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle); + psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter); + psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass); + psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing); + } else { + psWarning("no inputs with surviving detections. output header will be incomplete"); + } + psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); + psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp); + psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr); + psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom); + + //field in header that tells about the CMF version + char cmfVersion[8]; + sprintf(cmfVersion, PS1_DV_FORMAT, args->version); + psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion); + + // Find the total number of detections + + long total = 0; + for (long i=0; in; i++) { + ppMopsDetections *det = detections->data[i]; + if (!det) { + continue; } - addFittedTrailColumn("X_EXT", NULL, 0); - addFittedTrailColumn("Y_EXT", NULL, 0); - addFittedTrailColumn("X_EXT_SIG", NULL, 0); - addFittedTrailColumn("Y_EXT_SIG", NULL, 0); - addFittedTrailColumn("EXT_INST_MAG", NULL, 0); - addFittedTrailColumn("EXT_INST_MAG_SIG", NULL, 0); - addFittedTrailColumn("NPARAMS", NULL, 0); - addFittedTrailColumn("EXT_WIDTH_MAJ", NULL, 0); - addFittedTrailColumn("EXT_WIDTH_MIN", NULL, 0); - addFittedTrailColumn("EXT_THETA", NULL, 0); - addFittedTrailColumn("EXT_WIDTH_MAJ_ERR", NULL, 0); - addFittedTrailColumn("EXT_WIDTH_MIN_ERR", NULL, 0); - addFittedTrailColumn("EXT_THETA_ERR", NULL, 0); - } - } - if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) { - psError(psErrorCodeLast(), false, "Unable to write table"); - return false; - } - psFree(table); - } - - psFree(header); - psFitsClose(fits); - - psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output); - - return true; + total += det->num; + } + + psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output); + + if (total == 0) { + // Write dummy table + psMetadata *row = psMetadataAlloc(); // Output row + psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN); + psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0); + psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN); + psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN); + psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0); + + psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0); + psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN); + + if (args->version == 2) { + // Write data of version 2 (see ICD) + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", + NAN); + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags (group 2)", + 0); + psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", + 0); + psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", + 0); + } + + if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write empty table."); + psFree(header); + psFree(row); + return false; + } + psFree(row); + } else { + +#define addColumn(_outName, _inName, _convertTo32) \ + if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \ + psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \ + return false; \ + } + + // Allocate the output table + psMetadata *table = psMetadataAlloc(); + addColumn("RA", "RA_PSF", 0); + addColumn("RA_ERR", NULL, 0); // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG + addColumn("DEC", "DEC_PSF", 0); + addColumn("DEC_ERR", NULL, 0); // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG + addColumn("MAG", "PSF_INST_MAG", 0); + addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0); + addColumn("PSF_CHI2", "PSF_CHISQ", 0); + addColumn("PSF_DOF", "PSF_NDOF", 1); + addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0); + addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0); + addColumn("PSF_MAJOR", NULL, 0); + addColumn("PSF_MINOR", NULL, 0); + addColumn("PSF_THETA", NULL, 0); + addColumn("PSF_QUALITY", "PSF_QF", 0); + addColumn("PSF_NPIX", NULL, 1); + addColumn("MOMENTS_XX", NULL, 0); + addColumn("MOMENTS_XY", NULL, 0); + addColumn("MOMENTS_YY", NULL, 0); + addColumn("N_POS", "DIFF_NPOS", 1); + addColumn("F_POS", "DIFF_FRATIO", 0); + addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0); + addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0); + addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0); + addColumn("FLAGS", "FLAGS", 1); + addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID"); + if (args->version == 2) { + addColumn("IPP_IDET", NULL, 1); + addColumn("PSF_INST_FLUX", NULL, 0); + addColumn("PSF_INST_FLUX_SIG", NULL, 0); + addColumn("AP_MAG", NULL, 0); + addColumn("AP_MAG_RAW", NULL, 0); + addColumn("AP_MAG_RADIUS", NULL, 0); + addColumn("AP_FLUX", NULL, 0); + addColumn("AP_FLUX_SIG", NULL, 0); + addColumn("PEAK_FLUX_AS_MAG", NULL, 0); + addColumn("CAL_PSF_MAG", NULL, 0); + addColumn("CAL_PSF_MAG_SIG", NULL, 0); + addColumn("SKY", NULL, 0); + addColumn("SKY_SIGMA", NULL, 0); + addColumn("PSF_QF_PERFECT", NULL, 0); + addColumn("MOMENTS_R1", NULL, 0); + addColumn("MOMENTS_RH", NULL, 0); + addColumn("KRON_FLUX", NULL, 0); + addColumn("KRON_FLUX_ERR", NULL, 0); + addColumn("KRON_FLUX_INNER", NULL, 0); + addColumn("KRON_FLUX_OUTER", NULL, 0); + addColumn("DIFF_R_P", NULL, 0); + addColumn("DIFF_SN_P", NULL, 0); + addColumn("DIFF_R_M", NULL, 0); + addColumn("DIFF_SN_M", NULL, 0); + addColumn("FLAGS2", NULL, 1); + addColumn("IPP_IDET", NULL, 0); + addColumn("N_FRAMES", NULL, 0); + addColumn("PADDING", NULL, 0); + } + if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) { + psError(psErrorCodeLast(), false, "Unable to write table"); + return false; + } + psFree(table); + } + + psFree(header); + psFitsClose(fits); + + psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output); + + return true; } -// extension parameter values: -// 0: SkyChip.psf -// 1: SkyChip.xfit -// Any other value is ignored -static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, int extension, char *outColumnName, char *inColumnName, bool convertTo32) +static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32) { - if (inColumnName == NULL) { - inColumnName = outColumnName; - } - psVector *out = NULL; - if (convertTo32) { - // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers - // We want to write 32 bits to the output. - int next = 0; - for (long i=0; in; i++) { - ppMopsDetections *det = detections->data[i]; - if (!det || det->num == 0) { - // no detections survived for this input - continue; - } - psVector *in = NULL; - if (extension == 0) { - in = psMetadataLookupVector(NULL, det->table, inColumnName); - } else if (extension == 1) { - in = psMetadataLookupVector(NULL, det->fittedTrails, inColumnName); - } else { - //Error? - } - if (!in) { - psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s (convertTo32 is true)", inColumnName); - return false; - } - if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) { - psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d", - inColumnName, in->type.type); - return false; - } - if (out == NULL) { - // First time through set up the output vector and the copy parameters - if (in->type.type == PS_TYPE_S64) { - out = psVectorAlloc(outputSize, PS_TYPE_S32); - } else { - out = psVectorAlloc(outputSize, PS_TYPE_U32); - } - } - for (long d=0; d < det->num; d++) { - if (in->type.type == PS_TYPE_S64) { - out->data.S32[next++] = in->data.S64[d]; - } else { - out->data.U32[next++] = in->data.U64[d]; - } - } - } - } else { - void *next = NULL; - int elementSize = 0; // size of elements in vector... We are making assumptions here about the organization of primitives in memory so we can use memcopy - for (long i=0; in; i++) { - ppMopsDetections *det = detections->data[i]; - if (!det || det->num == 0) { - // no detections survived for this input - continue; - } - psVector *in = NULL; - if (extension == 0) { - psTrace("ppMops.write", 1, "Using det->table for %s", inColumnName); - in = psMetadataLookupVector(NULL, det->table, inColumnName); - } else if (extension == 1) { - psTrace("ppMops.write", 1, "Using det->fittedTrails for %s", inColumnName); - in = psMetadataLookupVector(NULL, det->fittedTrails, inColumnName); - } else { - //Error? - } - if (!in) { - psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s (convertTo32 is false)", inColumnName); - out = psVectorAlloc(outputSize, PS_TYPE_F32); - psVectorInit(out, NAN); - psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out); - return false; - } - if (out == NULL) { - // First time through set up the output vector and the copy parameters - out = psVectorAlloc(outputSize, in->type.type); - next = (void *) out->data.U8; - switch (in->type.type) { - case PS_TYPE_S8: - elementSize = sizeof(psS8); - break; - case PS_TYPE_U8: - elementSize = sizeof(psU8); - break; - case PS_TYPE_S16: - elementSize = sizeof(psS16); - break; - case PS_TYPE_U16: - elementSize = sizeof(psU16); - break; - case PS_TYPE_S32: - elementSize = sizeof(psS32); - break; - case PS_TYPE_U32: - elementSize = sizeof(psU32); - break; - case PS_TYPE_S64: - elementSize = sizeof(psS64); - break; - case PS_TYPE_U64: - elementSize = sizeof(psU64); - break; - case PS_TYPE_F32: - elementSize = sizeof(psF32); - break; - case PS_TYPE_F64: - elementSize = sizeof(psF64); - break; - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type); - return false; - } - } - // We are doing nasty things here so we can use memcpy. - // It would be safer to do a proper loop over the elements. - long toCopy = det->num * elementSize; - memcpy(next, in->data.U8, toCopy); - next += toCopy; - } - } - // Finally add the new column to the output table - psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out); - psFree(out); // drop reference - return true; + if (inColumnName == NULL) { + inColumnName = outColumnName; + } + + psVector *out = NULL; + if (convertTo32) { + // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers + // We want to write 32 bits to the output. + int next = 0; + for (long i=0; in; i++) { + ppMopsDetections *det = detections->data[i]; + if (!det || det->num == 0) { + // no detections survived for this input + continue; + } + psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName); + if (!in) { + psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName); + return false; + } + if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) { + psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d", + inColumnName, in->type.type); + return false; + } + if (out == NULL) { + // First time through set up the output vector and the copy parameters + if (in->type.type == PS_TYPE_S64) { + out = psVectorAlloc(outputSize, PS_TYPE_S32); + } else { + out = psVectorAlloc(outputSize, PS_TYPE_U32); + } + } + for (long d=0; d < det->num; d++) { + if (in->type.type == PS_TYPE_S64) { + out->data.S32[next++] = in->data.S64[d]; + } else { + out->data.U32[next++] = in->data.U64[d]; + } + } + } + } else { + void *next = NULL; + int elementSize = 0; // size of elements in vector... We are making assumptions here about the organization of primitives in memory so we can use memcopy + for (long i=0; in; i++) { + ppMopsDetections *det = detections->data[i]; + if (!det || det->num == 0) { + // no detections survived for this input + continue; + } + psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName); + if (!in) { + psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName); + return false; + } + if (out == NULL) { + // First time through set up the output vector and the copy parameters + out = psVectorAlloc(outputSize, in->type.type); + next = (void *) out->data.U8; + switch (in->type.type) { + case PS_TYPE_S8: + elementSize = sizeof(psS8); + break; + case PS_TYPE_U8: + elementSize = sizeof(psU8); + break; + case PS_TYPE_S16: + elementSize = sizeof(psS16); + break; + case PS_TYPE_U16: + elementSize = sizeof(psU16); + break; + case PS_TYPE_S32: + elementSize = sizeof(psS32); + break; + case PS_TYPE_U32: + elementSize = sizeof(psU32); + break; + case PS_TYPE_S64: + elementSize = sizeof(psS64); + break; + case PS_TYPE_U64: + elementSize = sizeof(psU64); + break; + case PS_TYPE_F32: + elementSize = sizeof(psF32); + break; + case PS_TYPE_F64: + elementSize = sizeof(psF64); + break; + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type); + return false; + + } + } + // We are doing nasty things here so we can use memcpy. + // It would be safer to do a proper loop over the elements. + long toCopy = det->num * elementSize; + memcpy(next, in->data.U8, toCopy); + next += toCopy; + } + } + + // Finally add the new column to the output table + psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out); + psFree(out); // drop reference + + return true; } - static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName) { - psVector *out = psVectorAlloc(total, PS_TYPE_S64); - long next = 0; - for (long i = 0; in; i++) { - ppMopsDetections *det = detections->data[i]; - if (!det) { - continue; - } - psS64 diffSkyfileId = det->diffSkyfileId; - for (long j = 0; j < det->num; j++) { - out->data.S64[next++] = diffSkyfileId; - } - } - psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out); - psFree(out); - return true; + psVector *out = psVectorAlloc(total, PS_TYPE_S64); + long next = 0; + for (long i = 0; in; i++) { + ppMopsDetections *det = detections->data[i]; + if (!det) { + continue; + } + psS64 diffSkyfileId = det->diffSkyfileId; + for (long j = 0; j < det->num; j++) { + out->data.S64[next++] = diffSkyfileId; + } + } + psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out); + psFree(out); + return true; }