Index: branches/eam_branches/20090715/ppMops/ICDlite.txt =================================================================== --- branches/eam_branches/20090715/ppMops/ICDlite.txt (revision 25401) +++ branches/eam_branches/20090715/ppMops/ICDlite.txt (revision 25401) @@ -0,0 +1,192 @@ +== General notes == + +Though efforts have been made to avoid this from occurring, where +header keyword names are longer than the usual 8 character maximum, +the HIERARCH convention will be used. + +Celestial coordinates are in ICRS. + +Magnitudes are AB magnitudes (flux zero point is 3631 Jy). + +Unless otherwise noted, error values are statistical only, and include +no contribution by systematic errors. + +== Primary header == + +The primary header will not contain any information beyond that +written automatically by cfitsio, and is present only for compliance +with the FITS standard. + +== FITS table == + +The FITS table will include all detections within an exposure. + +Duplicates from overlapping skycells will have been filtered out by +IPP, with the measurements from the source closest to the centre of a +skycell included (and all others discarded). + +=== Header === + +The header will contain information relevant to the exposure as a whole. + + * Version information + * SWSOURCE (string): source of software (e.g., "60eb6cdc-a59c-4636-a4e0-dba66a9721fd") + * SWVERSN (string): version of software (e.g., "trunk/ppMops@24658") + + * Provenance information + * EXP_NAME (string): Exposure name (e.g., "o1234g5678") + * EXP_ID (S64): Exposure identifier + * CHIP_ID (S64): Chip stage identifier + * CAM_ID (S64): Camera stage identifier + * FAKE_ID (S64): Fake stage identifier + * WARP_ID (S64): Warp stage identifier + * DIFF_ID (S64): Diff stage identifier + * DIFF_POS (boolean): Sense of subtraction; T for forward, F for backward + + * Exposure details + * MJD-OBS (F64): TAI MJD of exposure mid-point + * RA (string): Reported Right Ascension of telescope boresight, sexagesimal hours (e.g., "12:34:56.789") + * DEC (string): Reported Declination of telescope boresight, sexagesimal degrees (e.g., "-12:34:56.78") + * TEL_ALT (F64): Reported telescope altitude, degrees + * TEL_AZ (F64): Reported telescope azimuth, degrees + * EXPTIME (F32): Exposure time, seconds + * ROTANGLE (F64): Rotator angle, degrees + * FILTER (string): Filter name (e.g., "r.00000") + * AIRMASS (F32): Airmass for exposure + * OBSCODE (string): IAU observatory code (i.e., "F51" for PS1) + * SEEING (F32): Measured seeing at diff stage, arcsec + * MAGZP (F32): Magnitude zero point + * MAGZPERR (F32): Error in magnitude zero point + * ASTRORMS (F32): RMS of astrometric fit, arcsec + + * Detection efficiency + * DE_MAGnn (F32): Magnitude (calibrated) for detection efficiency + * DE_EFFnn (F32): Detection efficiency (0..1) + +==== Current limitations ==== + +The IPP is not yet calculating detection efficiencies (it is still +being developed). Further, it is not yet clear how to merge the +detection efficiency measurements for different skycells. Until these +issues are resolved, the detection efficiency values will be fake. + +==== Differences from original ICD ==== + + * SWSOURCE, SWVERSN has replaced TABLEVER which was never defined + * All of the "Provenance information" has replaced FPA_ID, for better tracking of the source of each detection + * SEEING and MAGZP added as indicators of the quality of the exposure + * MAGZPERR added to make absolute (i.e., across exposures) magnitude errors more accurate + * ASTRORMS added to make absolute (i.e., across exposures) astrometry errors more accurate + * DE_MAGnn and DE_EFFnn replace DE1 through DE10, which were never well defined + * Removed LIMITMAG, which was never well defined, and is unnecessary given the DE_MAGnn and DE_EFFnn + +=== Table === + +The table will contain information relevant to the individual +detections within the exposure. + + * RA (F64): Right Ascension of detection centre, degrees + * RA_ERR (F64): Error in RA, degrees + * DEC (F64): Declination of detection centre, degrees + * DEC_ERR (F64): Error in DEC, degrees + * MAG (F32): Calibrated magnitude of detection + * MAG_ERR (F32): Error in MAG + * STARPSF (F32): A PSF/extended source separator + * ANGLE (F64): Angle of trail fit to source, degrees E of N + * ANGLE_ERR (F64): Error in ANGLE, degrees + * LENGTH (F32): Length of trail fit to source, degrees + * LENGTH_ERR (F32): Error in LENGTH, degrees + * FLAGS (S32): IPP detection flags, bit mask + * DIFF_SKYFILE_ID (S64): IPP diff_skyfile_id for source + +==== Current limitations ==== + +The IPP is not yet fitting trails to sources in the difference images. +Until this is being done, the ANGLE, ANGLE_ERR, LENGTH and LENGTH_ERR +values will be zero. + +The value being written as STARPSF is EXT_NSIGMA from the IPP CMF +files; it is not clear that this is what MOPS wants, but this is +probably irrelevant until the trail fitting has been implemented. + +==== Differences from original ICD ==== + + * RA_DEG, DEC_DEG renamed RA, DEC to match apparent naming policy + * *_SIG renamed *_ERR to be more clear + * ANG, ANG_SIG, LEN, LEN_SIG spelled out as ANGLE, ANGLE_ERR, LENGTH, LENGTH_ERR + * FLUX, FLUX_SIG renamed MAG, MAG_ERR since IPP writes magnitudes + * FLAGS added to allow additional weeding out of bad detections + * DIFF_SKYFILE_ID added to allow trace back to IPP diff skyfile, for postage stamps + + +=== Example === + +{{{ +XTENSION= 'BINTABLE' / binary table extension +BITPIX = 8 / 8-bit bytes +NAXIS = 2 / 2-dimensional binary table +NAXIS1 = 72 / width of table in bytes +NAXIS2 = 42032 / number of rows in table +PCOUNT = 0 / size of special data area +GCOUNT = 1 / one data group (required keyword) +TFIELDS = 13 / number of fields in each row +TTYPE1 = 'RA ' / label for field 1 +TFORM1 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE2 = 'RA_ERR ' / label for field 2 +TFORM2 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE3 = 'DEC ' / label for field 3 +TFORM3 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE4 = 'DEC_ERR ' / label for field 4 +TFORM4 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE5 = 'MAG ' / label for field 5 +TFORM5 = '1E ' / data format of field: 4-byte REAL +TTYPE6 = 'MAG_ERR ' / label for field 6 +TFORM6 = '1E ' / data format of field: 4-byte REAL +TTYPE7 = 'STARPSF ' / label for field 7 +TFORM7 = '1E ' / data format of field: 4-byte REAL +TTYPE8 = 'ANGLE ' / label for field 8 +TFORM8 = '1E ' / data format of field: 4-byte REAL +TTYPE9 = 'ANGLE_ERR' / label for field 9 +TFORM9 = '1E ' / data format of field: 4-byte REAL +TTYPE10 = 'LENGTH ' / label for field 10 +TFORM10 = '1E ' / data format of field: 4-byte REAL +TTYPE11 = 'LENGTH_ERR' / label for field 11 +TFORM11 = '1E ' / data format of field: 4-byte REAL +TTYPE12 = 'FLAGS ' / label for field 12 +TFORM12 = '1J ' / data format of field: 4-byte INTEGER +TZERO12 = 2147483648 / offset for unsigned integers +TSCAL12 = 1 / data are not scaled +TTYPE13 = 'DIFF_SKYFILE_ID' / label for field 13 +TFORM13 = '1K ' / data format of field: 8-byte INTEGER +SWSOURCE= '60eb6cdc-a59c-4636-a4e0-dba66a9721fd' / Software source +SWVERSN = 'branches/pap_mops/ppMops@25227' / Software version +HISTORY ppMops at 2009-09-02T03:48:46.695783 +HISTORY psLib version: branches/pap_mops/psLib@25227 +HISTORY psLib source: 60eb6cdc-a59c-4636-a4e0-dba66a9721fd +HISTORY ppMops version: branches/pap_mops/ppMops@25227 +HISTORY ppMops source: 60eb6cdc-a59c-4636-a4e0-dba66a9721fd +EXP_NAME= 'o4995g0129o' / Exposure name +EXP_ID = 77164 / Exposure identifier +CHIP_ID = 24019 / Chip stage identifier +CAM_ID = 17726 / Cam stage identifier +FAKE_ID = 10227 / Fake stage identifier +WARP_ID = 8842 / Warp stage identifier +DIFF_ID = 0 / Diff stage identifier +DIFF_POS= F / Positive subtraction? +MJD-OBS = 54995.4740598313 / MJD of exposure midpoint +RA = '18:25:01.988' / Right Ascension of boresight +DEC = '-17:20:40.069' / Declination of boresight +TEL_ALT = 51.951873 / Telescope altitude +TEL_AZ = 179.483883 / Telescope azimuth +EXPTIME = 38. / Exposure time (sec) +ROTANGLE= 333.1039 / Rotator position angle +FILTER = 'r.00000 ' / Filter name +AIRMASS = 1.269 / Airmass of exposure +OBSCODE = 'F51 ' / IAU Observatory code +SEEING = 1.678401 / Mean seeing +MAGZP = 28.65226 / Magnitude zero point +MAGZPERR= 0.353063 / Error in magnitude zero point +ASTRORMS= 0.3111496 / RMS of astrometric fit +EXTNAME = 'MOPS_TRANSIENT_DETECTIONS' +END +}}} Index: branches/eam_branches/20090715/ppMops/src/Makefile.am =================================================================== --- branches/eam_branches/20090715/ppMops/src/Makefile.am (revision 25022) +++ branches/eam_branches/20090715/ppMops/src/Makefile.am (revision 25401) @@ -28,5 +28,9 @@ ppMops.c \ ppMopsVersion.c \ - ppMopsData.c + ppMopsArguments.c \ + ppMopsDetections.c \ + ppMopsRead.c \ + ppMopsWrite.c \ + ppMopsMerge.c noinst_HEADERS = \ Index: branches/eam_branches/20090715/ppMops/src/ppMops.c =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMops.c (revision 25022) +++ branches/eam_branches/20090715/ppMops/src/ppMops.c (revision 25401) @@ -6,19 +6,42 @@ int main(int argc, char *argv[]) { - if (argc != 4) { - fprintf(stderr, "Insufficient arguments.\n"); - fprintf(stderr, "Usage: %s DETECTIONS ZP OUTPUT\n", argv[0]); + psLibInit(NULL); + + ppMopsArguments *args = ppMopsArgumentsParse(argc, argv); // Parsed arguments + if (!args) { + psErrorStackPrint(stderr, "Error parsing arguments"); exit(PS_EXIT_CONFIG_ERROR); } - ppMopsData *data = ppMopsDataAlloc(); // Configuration data - data->detections = psStringCopy(argv[1]); - data->zp = atof(argv[2]); - data->output = psStringCopy(argv[3]); - - if (!isfinite(data->zp)) { - psErrorStackPrint(stderr, "Zero point is unknown\n"); - exit(PS_EXIT_CONFIG_ERROR); - } + psArray *detections = ppMopsRead(args); // Detections from each input + if (!detections) { + psErrorStackPrint(stderr, "Unable to read detections"); + exit(PS_EXIT_SYS_ERROR); + } + + ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections + psFree(detections); + if (!merged) { + psErrorStackPrint(stderr, "Unable to merge detections"); + exit(PS_EXIT_SYS_ERROR); + } + + if (!ppMopsWrite(merged, args)) { + psErrorStackPrint(stderr, "Unable to write detections"); + exit(PS_EXIT_SYS_ERROR); + } + + psFree(merged); + psFree(args); + + psLibFinalize(); + + return PS_EXIT_SUCCESS; +} + + +#if 0 + ps + psArray *detections = NULL; // Detections @@ -120,13 +143,9 @@ double alt = psMetadataLookupF64(NULL, header, "FPA.ALT"); double az = psMetadataLookupF64(NULL, header, "FPA.AZ"); - int imageid = psMetadataLookupS32(NULL, header, "IMAGEID"); + 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")); - - // XXX This is wrong - int fpaid = psMetadataLookupS32(NULL, header, "IMAGEID"); - psMetadataAddStr(outHeader, PS_LIST_TAIL, "RA", 0, "Right ascension of boresight", ra); @@ -139,11 +158,13 @@ psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", alt); psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", az); - psMetadataAddS32(outHeader, PS_LIST_TAIL, "DIFFIMID", 0, "Difference image identifier", imageid); - psMetadataAddS32(outHeader, PS_LIST_TAIL, "FPA_ID", 0, "Exposure identifier", fpaid); + 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)", 25.0); + 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); @@ -210,4 +231,4 @@ psFree(data); - return PS_EXIT_SUCCESS; -} +#endif + Index: branches/eam_branches/20090715/ppMops/src/ppMops.h =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMops.h (revision 25022) +++ branches/eam_branches/20090715/ppMops/src/ppMops.h (revision 25401) @@ -11,14 +11,63 @@ PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SKY_FAILURE) // Flags to exclude - // Configuration data typedef struct { - psString detections; // Detections filename - float zp; // Magnitude zero point + psArray *input; // Input filenames + psString exp_name; // Exposure name + psS64 exp_id; // Exposure identifier + psS64 chip_id; // Chip stage identifier + psS64 cam_id; // Camera stage identifier + psS64 fake_id; // Fake stage identifier + psS64 warp_id; // Warp stage identifier + psS64 diff_id; // Diff stage identifier + bool positive; // Sense of subtraction, T=positive, F=negative + float zp, zpErr; // Magnitude zero point and error + float rmsAstrom; // Astrometric solution RMS psString output; // Output filename -} ppMopsData; +} ppMopsArguments; -// Allocator -ppMopsData *ppMopsDataAlloc(void); +/// Parse arguments +ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]); + +typedef struct { + psString raBoresight, decBoresight; // RA,Dec of telescope boresight + psString filter; // Filter for exposure + float airmass; // Airmass of exposure + float exptime; // Exposure time + double posangle; // Position angle + double alt, az; // Telescope altitude and azimuth + double mjd; // Modified Julian Date + float seeing; // Seeing of exposure + long num; // Number of detections + psVector *x, *y; // Image coordinates + psVector *ra, *dec; // Sky coordinates + psVector *raErr, *decErr; // Error in sky coordinates + psVector *mag, *magErr; // Magnitude and associated error + psVector *extended; // Measure of extendedness + psVector *angle, *angleErr; // Angle of trail and associated error + psVector *length, *lengthErr; // Length of trail and associated error + psVector *flags; // psphot flags + psVector *diffSkyfileId; // Identifier for source image + psVector *naxis1, *naxis2; // Size of image + psVector *mask; // Mask for detections +} ppMopsDetections; + +ppMopsDetections *ppMopsDetectionsAlloc(long num); + +/// Copy a detection +bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index); + +/// Purge the detections list of masked detections +bool ppMopsDetectionsPurge(ppMopsDetections *detections); + + +/// Read detections +psArray *ppMopsRead(const ppMopsArguments *args); + +/// Merge detections +ppMopsDetections *ppMopsMerge(const psArray *detections); + +/// Write detections +bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args); Index: branches/eam_branches/20090715/ppMops/src/ppMopsArguments.c =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMopsArguments.c (revision 25401) +++ branches/eam_branches/20090715/ppMops/src/ppMopsArguments.c (revision 25401) @@ -0,0 +1,107 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppMops.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments // Command-line arguments + ) +{ + fprintf(stderr, "\nPan-STARRS IPP-MOPS detection translator\n\n"); + fprintf(stderr, "Usage: %s INPUT_LIST OUTPUT_NAME\n", program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +static void mopsArgumentsFree(ppMopsArguments *args) +{ + psFree(args->input); + psFree(args->exp_name); + psFree(args->output); + return; +} + +ppMopsArguments *ppMopsArgumentsAlloc(void) +{ + ppMopsArguments *args = psAlloc(sizeof(ppMopsArguments)); // Data to return + psMemSetDeallocator(args, (psFreeFunc)mopsArgumentsFree); + + args->input = NULL; + args->exp_name = NULL; + args->exp_id = 0; + args->chip_id = 0; + args->cam_id = 0; + args->fake_id = 0; + args->warp_id = 0; + args->diff_id = 0; + args->zp = NAN; + args->positive = true; + args->zpErr = NAN; + args->rmsAstrom = NAN; + args->output = NULL; + + return args; +} + + +ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]) +{ + assert(argv); + + psTrace("ppMops.args", 1, "Parsing command-line arguments\n"); + + psArgumentVerbosity(&argc, argv); + + psMetadata *arguments = psMetadataAlloc(); // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-exp_name", 0, "Exposure name", NULL); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-exp_id", 0, "Exposure identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-chip_id", 0, "Chip stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-cam_id", 0, "Camera stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-fake_id", 0, "Fake stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-warp_id", 0, "Warp stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-diff_id", 0, "Diff stage identifier", 0); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Inverse subtraction?", false); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Magnitude zero point", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp_error", 0, "Error in magnitude zero point", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-astrom_rms", 0, "Astrometric solution RMS", NAN); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) { + usage(argv[0], arguments); + } + + ppMopsArguments *args = ppMopsArgumentsAlloc(); // Arguments, to return + + psString inList = psSlurpFilename(argv[1]); // List of filenames + args->input = psStringSplitArray(inList, "\n", false); + psFree(inList); + if (!args->input || args->input->n == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No inputs provided."); + return NULL; + } + args->output = psStringCopy(argv[2]); + + args->exp_name = psMetadataLookupStr(NULL, arguments, "-exp_name"); + args->exp_id = psMetadataLookupS64(NULL, arguments, "-exp_id"); + args->chip_id = psMetadataLookupS64(NULL, arguments, "-chip_id"); + args->cam_id = psMetadataLookupS64(NULL, arguments, "-cam_id"); + args->fake_id = psMetadataLookupS64(NULL, arguments, "-fake_id"); + args->warp_id = psMetadataLookupS64(NULL, arguments, "-warp_id"); + args->diff_id = psMetadataLookupS64(NULL, arguments, "-diff_id"); + args->positive = !psMetadataLookupBool(NULL, arguments, "-inverse"); // NOTE: negated + + args->zp = psMetadataLookupF32(NULL, arguments, "-zp"); + args->zpErr = psMetadataLookupF32(NULL, arguments, "-zp_error"); + args->rmsAstrom = psMetadataLookupF32(NULL, arguments, "-astrom_rms"); + + psTrace("ppMops.args", 1, "Done parsing command-line arguments\n"); + + return args; +} Index: branches/eam_branches/20090715/ppMops/src/ppMopsData.c =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMopsData.c (revision 25022) +++ (revision ) @@ -1,29 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include "ppMops.h" - -static void mopsDataFree(ppMopsData *data) -{ - psFree(data->detections); - psFree(data->output); - return; -} - -ppMopsData *ppMopsDataAlloc(void) -{ - ppMopsData *data = psAlloc(sizeof(ppMopsData)); // Data to return - psMemSetDeallocator(data, (psFreeFunc)mopsDataFree); - - data->detections = NULL; - data->zp = NAN; - data->output = NULL; - - return data; -} - - Index: branches/eam_branches/20090715/ppMops/src/ppMopsDetections.c =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMopsDetections.c (revision 25401) +++ branches/eam_branches/20090715/ppMops/src/ppMopsDetections.c (revision 25401) @@ -0,0 +1,204 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +static void mopsDetectionsFree(ppMopsDetections *det) +{ + psFree(det->raBoresight); + psFree(det->decBoresight); + psFree(det->filter); + psFree(det->x); + psFree(det->y); + psFree(det->ra); + psFree(det->dec); + psFree(det->raErr); + psFree(det->decErr); + psFree(det->mag); + psFree(det->magErr); + psFree(det->extended); + psFree(det->angle); + psFree(det->angleErr); + psFree(det->length); + psFree(det->lengthErr); + psFree(det->flags); + psFree(det->diffSkyfileId); + psFree(det->naxis1); + psFree(det->naxis2); + psFree(det->mask); + return; +} + +ppMopsDetections *ppMopsDetectionsAlloc(long num) +{ + ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return + psMemSetDeallocator(det, (psFreeFunc)mopsDetectionsFree); + + det->raBoresight = NULL; + det->decBoresight = NULL; + det->filter = NULL; + det->airmass = NAN; + det->exptime = NAN; + det->posangle = NAN; + det->alt = NAN; + det->az = NAN; + det->mjd = NAN; + det->seeing = NAN; + det->num = 0; + det->x = psVectorAllocEmpty(num, PS_TYPE_F32); + det->y = psVectorAllocEmpty(num, PS_TYPE_F32); + det->ra = psVectorAllocEmpty(num, PS_TYPE_F64); + det->dec = psVectorAllocEmpty(num, PS_TYPE_F64); + det->raErr = psVectorAllocEmpty(num, PS_TYPE_F64); + det->decErr = psVectorAllocEmpty(num, PS_TYPE_F64); + det->mag = psVectorAllocEmpty(num, PS_TYPE_F32); + det->magErr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->extended = psVectorAllocEmpty(num, PS_TYPE_F32); + det->angle = psVectorAllocEmpty(num, PS_TYPE_F32); + det->angleErr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->length = psVectorAllocEmpty(num, PS_TYPE_F32); + det->lengthErr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->flags = psVectorAllocEmpty(num, PS_TYPE_U32); + det->diffSkyfileId = psVectorAllocEmpty(num, PS_TYPE_S64); + det->naxis1 = psVectorAllocEmpty(num, PS_TYPE_S32); + det->naxis2 = psVectorAllocEmpty(num, PS_TYPE_S32); + det->mask = psVectorAllocEmpty(num, PS_TYPE_U8); + + return det; +} + + +ppMopsDetections *ppMopsDetectionsRealloc(ppMopsDetections *det, long num) +{ + det->x = psVectorRealloc(det->x, num); + det->y = psVectorRealloc(det->y, num); + det->ra = psVectorRealloc(det->ra, num); + det->dec = psVectorRealloc(det->dec, num); + det->raErr = psVectorRealloc(det->raErr, num); + det->decErr = psVectorRealloc(det->decErr, num); + det->mag = psVectorRealloc(det->mag, num); + det->magErr = psVectorRealloc(det->magErr, num); + det->extended = psVectorRealloc(det->extended, num); + det->angle = psVectorRealloc(det->angle, num); + det->angleErr = psVectorRealloc(det->angleErr, num); + det->length = psVectorRealloc(det->length, num); + det->lengthErr = psVectorRealloc(det->lengthErr, num); + det->flags = psVectorRealloc(det->flags, num); + det->diffSkyfileId = psVectorRealloc(det->diffSkyfileId, num); + det->naxis1 = psVectorRealloc(det->naxis1, num); + det->naxis2 = psVectorRealloc(det->naxis2, num); + det->mask = psVectorRealloc(det->mask, num); + + return det; +} + + +bool ppMopsDetectionsAdd(ppMopsDetections *det, float x, float y, double ra, double dec, + double raErr, double decErr, float mag, float magErr, float extended, + float angle, float angleErr, float length, float lengthErr, + psU32 flags, psS64 diffSkyfileId, int naxis1, int naxis2) +{ + psVectorAppend(det->x, x); + psVectorAppend(det->y, y); + psVectorAppend(det->ra, ra); + psVectorAppend(det->dec, dec); + psVectorAppend(det->raErr, raErr); + psVectorAppend(det->decErr, decErr); + psVectorAppend(det->mag, mag); + psVectorAppend(det->magErr, magErr); + psVectorAppend(det->extended, extended); + psVectorAppend(det->angle, angle); + psVectorAppend(det->angleErr, angleErr); + psVectorAppend(det->length, length); + psVectorAppend(det->lengthErr, lengthErr); + psVectorAppend(det->flags, flags); + psVectorAppend(det->diffSkyfileId, diffSkyfileId); + psVectorAppend(det->naxis1, naxis1); + psVectorAppend(det->naxis2, naxis2); + psVectorAppend(det->mask, 0); + return true; +} + + +bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index) +{ + psVectorAppend(target->x, source->x->data.F32[index]); + psVectorAppend(target->y, source->y->data.F32[index]); + psVectorAppend(target->ra, source->ra->data.F64[index]); + psVectorAppend(target->dec, source->dec->data.F64[index]); + psVectorAppend(target->raErr, source->raErr->data.F64[index]); + psVectorAppend(target->decErr, source->decErr->data.F64[index]); + psVectorAppend(target->mag, source->mag->data.F32[index]); + psVectorAppend(target->magErr, source->magErr->data.F32[index]); + psVectorAppend(target->extended, source->extended->data.F32[index]); + psVectorAppend(target->angle, source->angle->data.F32[index]); + psVectorAppend(target->angleErr, source->angleErr->data.F32[index]); + psVectorAppend(target->length, source->length->data.F32[index]); + psVectorAppend(target->lengthErr, source->lengthErr->data.F32[index]); + psVectorAppend(target->flags, source->flags->data.U32[index]); + psVectorAppend(target->diffSkyfileId, source->diffSkyfileId->data.S64[index]); + psVectorAppend(target->naxis1, source->naxis1->data.S32[index]); + psVectorAppend(target->naxis2, source->naxis2->data.S32[index]); + psVectorAppend(target->mask, 0); + target->num++; + return true; +} + + +bool ppMopsDetectionsPurge(ppMopsDetections *det) +{ + long num = 0; + for (long i = 0; i < det->num; i++) { + if (!det->mask->data.U8[i]) { + if (i == num) { + // No need to copy + num++; + continue; + } + det->x->data.F32[num] = det->x->data.F32[i]; + det->y->data.F32[num] = det->y->data.F32[i]; + det->ra->data.F64[num] = det->ra->data.F64[i]; + det->dec->data.F64[num] = det->dec->data.F64[i]; + det->raErr->data.F64[num] = det->raErr->data.F64[i]; + det->decErr->data.F64[num] = det->decErr->data.F64[i]; + det->mag->data.F32[num] = det->mag->data.F32[i]; + det->magErr->data.F32[num] = det->magErr->data.F32[i]; + det->extended->data.F32[num] = det->extended->data.F32[i]; + det->angle->data.F32[num] = det->angle->data.F32[i]; + det->angleErr->data.F32[num] = det->angleErr->data.F32[i]; + det->length->data.F32[num] = det->length->data.F32[i]; + det->lengthErr->data.F32[num] = det->lengthErr->data.F32[i]; + det->flags->data.U32[num] = det->flags->data.U32[i]; + det->diffSkyfileId->data.S64[num] = det->diffSkyfileId->data.S64[i]; + det->naxis1->data.S32[num] = det->naxis1->data.S32[i]; + det->naxis2->data.S32[num] = det->naxis2->data.S32[i]; + det->mask->data.U8[num] = 0; + num++; + } + } + det->x->n = num; + det->y->n = num; + det->ra->n = num; + det->dec->n = num; + det->raErr->n = num; + det->decErr->n = num; + det->mag->n = num; + det->magErr->n = num; + det->extended->n = num; + det->angle->n = num; + det->angleErr->n = num; + det->length->n = num; + det->lengthErr->n = num; + det->flags->n = num; + det->diffSkyfileId->n = num; + det->naxis1->n = num; + det->naxis2->n = num; + det->mask->n = num; + det->num = num; + return true; +} + Index: branches/eam_branches/20090715/ppMops/src/ppMopsMerge.c =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMopsMerge.c (revision 25401) +++ branches/eam_branches/20090715/ppMops/src/ppMopsMerge.c (revision 25401) @@ -0,0 +1,171 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppMops.h" + +#define LEAF_SIZE 4 // Size of leaf +#define MATCH_RADIUS SEC_TO_RAD(1.0) // Matching radius +#define MJD_TOL 1.0/3600.0/24.0 // Tolerance for MJD matching +#define BORESIGHT_TOL SEC_TO_RAD(1.0) // Tolerance for boresight matching +#define EXPTIME_TOL 1.0e-3 // Tolerance for exposure time matching +#define POSANGLE_TOL SEC_TO_RAD(1.0) // Tolerance for position angle matching +#define AIRMASS_TOL 1.0e-3 // Tolerance for airmass matching + +// Get distance from detection to centre of image +static float mergeDistance(const ppMopsDetections *detections, // Detections of interest + long index // Index for source of interest + ) +{ + float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0; + float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0; + return PS_SQR(dx) + PS_SQR(dy); +} + + +ppMopsDetections *ppMopsMerge(const psArray *detections) +{ + PS_ASSERT_ARRAY_NON_NULL(detections, NULL); + + psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n); + + ppMopsDetections *merged = NULL; // Merged list + int num = 1; // Number of merged files + for (int i = 0; i < detections->n; i++) { + ppMopsDetections *det = detections->data[i]; // Detections of interest + if (!det) { + psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i); + continue; + } else if (det->num == 0) { + psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i); + continue; + } + num++; + if (!merged) { + psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i); + merged = psMemIncrRefCounter(det); + continue; + } + psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i); + + // XXX compare exposure properties + if (strcmp(merged->raBoresight, det->raBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s", + merged->raBoresight, det->raBoresight); + return NULL; + } + if (strcmp(merged->decBoresight, det->decBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s", + merged->decBoresight, det->decBoresight); + return NULL; + } + if (strcmp(merged->filter, det->filter) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s", + merged->filter, det->filter); + return NULL; + } + + if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f", + merged->airmass, det->airmass); + return NULL; + } + if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f", + merged->exptime, det->exptime); + return NULL; + } + if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f", + merged->posangle, det->posangle); + return NULL; + } + if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf", + merged->alt, det->alt); + return NULL; + } + if (fabs(merged->az - det->az) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf", + merged->az, det->az); + return NULL; + } + if (fabs(merged->mjd - det->mjd) > MJD_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf", + merged->mjd, det->mjd); + return NULL; + } + + merged->seeing += det->seeing; // Taking average + + psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree + if (!tree) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree"); + psFree(merged); + return NULL; + } + + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest + for (int j = 0; j < det->num; j++) { + coords->data.F64[0] = det->ra->data.F64[j]; + coords->data.F64[1] = det->dec->data.F64[j]; + psVector *indices = psTreeAllWithin(tree, coords, MATCH_RADIUS); // Indices for matching sources + if (!indices) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to search for matches"); + psFree(coords); + psFree(tree); + psFree(merged); + return NULL; + } + if (indices->n == 0) { + psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i); + psFree(indices); + ppMopsDetectionsCopySingle(merged, det, j); + continue; + } + psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i); + + // Which one do we keep? + float bestDistance = INFINITY; // Best distance to centre + long bestIndex = -1; // Index with best distance + for (int k = 0; k < indices->n; k++) { + long index = indices->data.S64[k]; // Index of point + float distance = mergeDistance(merged, index); // Distance to centre of image + if (distance < bestDistance) { + bestDistance = distance; + bestIndex = index; + } + } + + float distance = mergeDistance(det, j); // Distance to centre of image + if (distance < bestDistance) { + psTrace("ppMops.merge", 6, "New source clobbers old sources\n"); + // Blow away existing sources + for (int k = 0; k < indices->n; k++) { + long index = indices->data.S64[k]; // Index of point + merged->mask->data.U8[index] = 0xFF; + } + ppMopsDetectionsCopySingle(merged, det, j); + } else { + psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); + } + psFree(indices); + } + + psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num); + + psFree(tree); + ppMopsDetectionsPurge(merged); + } + + psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num); + + merged->seeing /= num; + + return merged; +} + Index: branches/eam_branches/20090715/ppMops/src/ppMopsRead.c =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMopsRead.c (revision 25401) +++ branches/eam_branches/20090715/ppMops/src/ppMopsRead.c (revision 25401) @@ -0,0 +1,166 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +psArray *ppMopsRead(const 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++) { + psFits *fits = psFitsOpen(inNames->data[i], "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; + } + + 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); + + 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 = 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") + + psMetadataLookupF32(NULL, header, "FWHM_MIN")); + + int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows + + psFree(header); + + psArray *table = psFitsReadTable(fits); // Table of interest + if (!table) { + psError(PS_ERR_IO, false, "Unable to read table %d", i); + return false; + } + psFitsClose(fits); + + double plateScale = 0.0; // Plate scale + long numGood = 0; // Number of good rows + for (long j = 0; j < size; j++) { + psMetadata *row = table->data[j]; // Row of interest + + psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS"); + if (flags & SOURCE_MASK) { + continue; + } + + det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF"); + det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF"); + det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF")); + det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF")); + det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG"); + det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG"); + det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA"); + det->angle->data.F32[numGood] = 0.0; + det->angleErr->data.F32[numGood] = 0.0; + det->length->data.F32[numGood] = 0.0; + det->lengthErr->data.F32[numGood] = 0.0; + det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS"); + det->diffSkyfileId->data.F32[numGood] = diffSkyfileId; + det->naxis1->data.S32[numGood] = naxis1; + det->naxis2->data.S32[numGood] = naxis2; + + // Calculate error in RA, Dec + double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); + double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); + double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); + double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); + + if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) || + !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) || + !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) || + !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle) || + (det->flags->data.U32[numGood] & SOURCE_MASK)) { + 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[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + det->mask->data.U8[numGood] = 0; + plateScale += scale; + numGood++; + } + det->seeing *= plateScale / numGood; + + det->x->n = numGood; + det->y->n = numGood; + det->ra->n = numGood; + det->dec->n = numGood; + det->raErr->n = numGood; + det->decErr->n = numGood; + det->mag->n = numGood; + det->magErr->n = numGood; + det->extended->n = numGood; + det->angle->n = numGood; + det->angleErr->n = numGood; + det->length->n = numGood; + det->lengthErr->n = numGood; + det->flags->n = numGood; + det->diffSkyfileId->n = numGood; + det->naxis1->n = numGood; + det->naxis2->n = numGood; + det->mask->n = numGood; + + det->num = numGood; + + if (isfinite(args->zp) && numGood > 0) { + psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); + } + + psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]); + + psFree(table); + detections->data[i] = det; + } + + psTrace("ppMops.read", 1, "Done reading input detections\n"); + + return detections; +} Index: branches/eam_branches/20090715/ppMops/src/ppMopsWrite.c =================================================================== --- branches/eam_branches/20090715/ppMops/src/ppMopsWrite.c (revision 25401) +++ branches/eam_branches/20090715/ppMops/src/ppMopsWrite.c (revision 25401) @@ -0,0 +1,119 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args) +{ + psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output); + + 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 = ppMopsSource(), version = ppMopsVersion(); + psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source); + psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version); + ppMopsVersionHeader(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); + + 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); + psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); + psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing); + 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); + + if (det->num == 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, "STARPSF", 0, "EXT_NSIGMA", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)", 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); + if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write empty table."); + psFree(header); + psFree(row); + return false; + } + psFree(row); + } else { + psArray *table = psArrayAlloc(det->num); // Table to write + for (long i = 0; i < det->num; i++) { + psMetadata *row = psMetadataAlloc(); // Output row + psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", + RAD_TO_DEG(det->ra->data.F64[i])); + psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", + det->raErr->data.F64[i]); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", + RAD_TO_DEG(det->dec->data.F64[i])); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", + det->decErr->data.F64[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", det->extended->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)", + det->angle->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)", + det->angleErr->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)", + det->length->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)", + det->lengthErr->data.F32[i]); + psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.U32[i]); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", + det->diffSkyfileId->data.S64[i]); + table->data[i] = row; + } + if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write table."); + psFree(header); + psFree(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; +}