Index: /trunk/ppTranslate/src/ppMops.c =================================================================== --- /trunk/ppTranslate/src/ppMops.c (revision 34278) +++ /trunk/ppTranslate/src/ppMops.c (revision 34279) @@ -67,5 +67,4 @@ } - psArray *detections = ppMopsRead(args); // Detections from each input if (!detections) { @@ -92,202 +91,9 @@ 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: /trunk/ppTranslate/src/ppMops.h =================================================================== --- /trunk/ppTranslate/src/ppMops.h (revision 34278) +++ /trunk/ppTranslate/src/ppMops.h (revision 34279) @@ -30,32 +30,4 @@ } 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[]); @@ -75,5 +47,7 @@ long numGood; // Number of "good" detections psS64 diffSkyfileId; // unique id for input skyfile - psMetadata *table; // Columns from the input file + 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) psVector *x, *y; // Image coordinates psVector *ra, *dec; // Sky coordinates @@ -103,4 +77,5 @@ /// @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: /trunk/ppTranslate/src/ppMopsArguments.c =================================================================== --- /trunk/ppTranslate/src/ppMopsArguments.c (revision 34278) +++ /trunk/ppTranslate/src/ppMopsArguments.c (revision 34279) @@ -34,5 +34,5 @@ fprintf(stderr, "\n"); psArgumentHelp(arguments); - fprintf(stderr, "\t\tCMF file version can be set to either 1 for PS1_DV1 or 2 for PS1_DV1\n"); + 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\tSee IPP-MOPS ICD for details\n"); psLibFinalize(); Index: /trunk/ppTranslate/src/ppMopsGetSkyChipPsfVersion.c =================================================================== --- /trunk/ppTranslate/src/ppMopsGetSkyChipPsfVersion.c (revision 34278) +++ /trunk/ppTranslate/src/ppMopsGetSkyChipPsfVersion.c (revision 34279) @@ -14,4 +14,7 @@ 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: /trunk/ppTranslate/src/ppMopsMerge.c =================================================================== --- /trunk/ppTranslate/src/ppMopsMerge.c (revision 34278) +++ /trunk/ppTranslate/src/ppMopsMerge.c (revision 34279) @@ -71,4 +71,5 @@ 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: /trunk/ppTranslate/src/ppMopsRead.c =================================================================== --- /trunk/ppTranslate/src/ppMopsRead.c (revision 34278) +++ /trunk/ppTranslate/src/ppMopsRead.c (revision 34279) @@ -8,201 +8,346 @@ #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; - } - - 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; + 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; + } + //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"); } - 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); - } - - psTrace("ppMops.read", 1, "Done reading input detections\n"); - - return detections; -} + } + } + + 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; +} + +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: /trunk/ppTranslate/src/ppMopsWrite.c =================================================================== --- /trunk/ppTranslate/src/ppMopsWrite.c (revision 34278) +++ /trunk/ppTranslate/src/ppMopsWrite.c (revision 34279) @@ -9,378 +9,419 @@ #include "ppTranslateVersion.h" -static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, char *outColName, char *inColName, bool convertTo32); +static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, int extension, 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; + 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); \ + } + 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; +} + +// 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) +{ + 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) { - continue; + 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); } - 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); + } + 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]; } - - 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; + } + } + } 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; } -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 = 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; }