Index: /branches/pap_mops/ppMops/src/Makefile.am =================================================================== --- /branches/pap_mops/ppMops/src/Makefile.am (revision 25161) +++ /branches/pap_mops/ppMops/src/Makefile.am (revision 25162) @@ -28,5 +28,9 @@ ppMops.c \ ppMopsVersion.c \ - ppMopsData.c + ppMopsArguments.c \ + ppMopsDetections.c \ + ppMopsRead.c \ + ppMopsWrite.c \ + ppMopsMerge.c noinst_HEADERS = \ Index: /branches/pap_mops/ppMops/src/ppMops.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMops.c (revision 25161) +++ /branches/pap_mops/ppMops/src/ppMops.c (revision 25162) @@ -6,22 +6,19 @@ int main(int argc, char *argv[]) { - if (argc != 7) { - fprintf(stderr, "Insufficient arguments.\n"); - fprintf(stderr, "Usage: %s DETECTIONS ZP EXP_ID EXP_NAME DIRECTION OUTPUT\n", argv[0]); + 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->exp_id = atoll(argv[3]); - data->exp_name = psStringCopy(argv[4]); - data->direction = atoi(argv[5]); - data->output = psStringCopy(argv[6]); - - 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); + } + +#if 0 + ps + psArray *detections = NULL; // Detections @@ -211,4 +208,6 @@ psFree(data); +#endif + return PS_EXIT_SUCCESS; } Index: /branches/pap_mops/ppMops/src/ppMops.h =================================================================== --- /branches/pap_mops/ppMops/src/ppMops.h (revision 25161) +++ /branches/pap_mops/ppMops/src/ppMops.h (revision 25162) @@ -11,17 +11,54 @@ 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 - psString exp_name; // Exposure name - bool direction; // Direction of subtraction, 1=positive, 0=negative + 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 +} ppMopsDetections; + +ppMopsDetections *ppMopsDetectionsAlloc(long num); + +/// Read detections +psArray *ppMopsRead(const ppMopsArguments *args); + +/// Merge detections +ppMopsDetections *ppMopsDetectionsMerge(const psArray *detections); + +/// Write detections +bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args); Index: /branches/pap_mops/ppMops/src/ppMopsArguments.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsArguments.c (revision 25162) +++ /branches/pap_mops/ppMops/src/ppMopsArguments.c (revision 25162) @@ -0,0 +1,103 @@ +#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); + + 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", 1, "Done reading command-line arguments\n"); + + return args; +} Index: anches/pap_mops/ppMops/src/ppMopsData.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsData.c (revision 25161) +++ (revision ) @@ -1,32 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include "ppMops.h" - -static void mopsDataFree(ppMopsData *data) -{ - psFree(data->detections); - psFree(data->exp_name); - 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->exp_id = 0; - data->exp_name = NULL; - data->output = NULL; - - return data; -} - - Index: /branches/pap_mops/ppMops/src/ppMopsDetections.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsDetections.c (revision 25162) +++ /branches/pap_mops/ppMops/src/ppMopsDetections.c (revision 25162) @@ -0,0 +1,68 @@ +#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); + 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 = num; + det->x = psVectorAlloc(num, PS_TYPE_F32); + det->y = psVectorAlloc(num, PS_TYPE_F32); + det->ra = psVectorAlloc(num, PS_TYPE_F64); + det->dec = psVectorAlloc(num, PS_TYPE_F64); + det->raErr = psVectorAlloc(num, PS_TYPE_F64); + det->decErr = psVectorAlloc(num, PS_TYPE_F64); + det->mag = psVectorAlloc(num, PS_TYPE_F32); + det->magErr = psVectorAlloc(num, PS_TYPE_F32); + det->extended = psVectorAlloc(num, PS_TYPE_F32); + det->angle = psVectorAlloc(num, PS_TYPE_F64); + det->angleErr = psVectorAlloc(num, PS_TYPE_F64); + det->length = psVectorAlloc(num, PS_TYPE_F32); + det->lengthErr = psVectorAlloc(num, PS_TYPE_F32); + det->flags = psVectorAlloc(num, PS_TYPE_U32); + det->diffSkyfileId = psVectorAlloc(num, PS_TYPE_S64); + + return det; +} + + Index: /branches/pap_mops/ppMops/src/ppMopsMerge.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsMerge.c (revision 25162) +++ /branches/pap_mops/ppMops/src/ppMopsMerge.c (revision 25162) @@ -0,0 +1,14 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +ppMopsDetections *ppMopsMerge(const psArray *detections) +{ + + return NULL; +} Index: /branches/pap_mops/ppMops/src/ppMopsRead.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsRead.c (revision 25162) +++ /branches/pap_mops/ppMops/src/ppMopsRead.c (revision 25162) @@ -0,0 +1,111 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +psArray *ppMopsRead(const ppMopsArguments *args) +{ + 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); + + 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, "FPA.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")); + + 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 + for (long j = 0; j < size; j++) { + psMetadata *row = table->data[j]; // Row of interest + det->x->data.F32[j] = psMetadataLookupF32(NULL, row, "X_PSF"); + det->y->data.F32[j] = psMetadataLookupF32(NULL, row, "Y_PSF"); + det->ra->data.F32[j] = psMetadataLookupF32(NULL, row, "RA_PSF"); + det->dec->data.F32[j] = psMetadataLookupF32(NULL, row, "DEC_PSF"); + det->mag->data.F32[j] = psMetadataLookupF32(NULL, row, "Y_PSF"); + det->magErr->data.F32[j] = psMetadataLookupF32(NULL, row, "Y_PSF"); + det->extended->data.F32[j] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA"); + det->angle->data.F32[j] = 0.0; + det->angleErr->data.F32[j] = 0.0; + det->length->data.F32[j] = 0.0; + det->lengthErr->data.F32[j] = 0.0; + det->flags->data.F32[j] = psMetadataLookupU32(NULL, row, "FLAGS"); + det->diffSkyfileId->data.F32[j] = diffSkyfileId; + + // 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"); + + // 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.F32[j] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + det->decErr->data.F32[j] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + plateScale += scale; + } + det->seeing *= plateScale / size; + + psFree(table); + detections->data[i] = det; + } + + return detections; +} Index: /branches/pap_mops/ppMops/src/ppMopsWrite.c =================================================================== --- /branches/pap_mops/ppMops/src/ppMopsWrite.c (revision 25162) +++ /branches/pap_mops/ppMops/src/ppMopsWrite.c (revision 25162) @@ -0,0 +1,14 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args) +{ + + return true; +}