Index: /trunk/ippToPsps/config/test/map.xml
===================================================================
--- /trunk/ippToPsps/config/test/map.xml (revision 27809)
+++ /trunk/ippToPsps/config/test/map.xml (revision 27809)
@@ -0,0 +1,95 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: /trunk/ippToPsps/config/test/tables.xml
===================================================================
--- /trunk/ippToPsps/config/test/tables.xml (revision 27809)
+++ /trunk/ippToPsps/config/test/tables.xml (revision 27809)
@@ -0,0 +1,109 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: /trunk/ippToPsps/src/Makefile.am
===================================================================
--- /trunk/ippToPsps/src/Makefile.am (revision 27808)
+++ /trunk/ippToPsps/src/Makefile.am (revision 27809)
@@ -30,4 +30,5 @@
ippToPspsVersion.c \
ippToPspsConfig.c \
+ ippToPspsBatchTest.c \
ippToPspsBatchInit.c \
ippToPspsBatchDetection.c \
Index: /trunk/ippToPsps/src/ippToPspsBatchTest.c
===================================================================
--- /trunk/ippToPsps/src/ippToPspsBatchTest.c (revision 27809)
+++ /trunk/ippToPsps/src/ippToPspsBatchTest.c (revision 27809)
@@ -0,0 +1,347 @@
+/** @file ippToPspsBatchDetection.c
+ *
+ * @ingroup ippToPsps
+ *
+ * @author IfA
+ * Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include
+#include "ippToPsps.h"
+#include "ippToPspsConfig.h"
+#include "ippToPspsTestEnums.h"
+#include "fitsio.h"
+
+// Gets flux from magnitude
+static __inline bool ippToPsps_getFlux(const float zeroPoint, const float exposureTime, const float magnitude, float* flux) {
+
+ // use uncalibrated instrumental flux
+ *flux = powf(10.0, -0.4*magnitude) / exposureTime;
+
+ // use calibrated flux in Janskys, where 3631 Jy is the zero point flux for the filter: constant over all filters because we're using AB magnitudes
+ //float flux = 3631 * powf(10.0, -0.4*(magnitude + zeroPoint)) / exposureTime;
+
+ // printf("mag=%f\texpTime=%f\tzeroPoint=%f\tflux=%f\n", magnitude, exposureTime, zeroPoint, flux);
+ return (!isfinite(*flux) || *flux < 0.000001) ? false : true;
+}
+
+/**
+ \brief test data for checking PSPS values correspond with IPP
+
+ */
+int ippToPsps_batchTest(IppToPsps *this) {
+
+ if (this->numOfInputFiles < 1) return PS_EXIT_DATA_ERROR;
+
+ int status = 0;
+ fitsfile *fitsIn;
+ int nKeys;
+
+ if (fits_open_file(&fitsIn, this->inputFiles[0], READONLY, &status)) {
+
+ fits_report_error(stderr, status);
+ return PS_EXIT_SYS_ERROR;
+ }
+
+ // get primary header and pull stuff out for later
+ fits_get_hdrspace(fitsIn, &nKeys, NULL, &status);
+
+ float zptObs, zeroPoint, exposureTime;
+ char filterType[20];
+ status=0; fits_read_key(fitsIn, TFLOAT, "ZPT_OBS", &zptObs, NULL, &status);
+ status=0; fits_read_key(fitsIn, TFLOAT, "ZPT_REF", &zeroPoint, NULL, &status);
+ status=0; fits_read_key(fitsIn, TFLOAT, "EXPREQ", &exposureTime, NULL, &status);
+ status=0; fits_read_key(fitsIn, TSTRING, "FILTERID", filterType, NULL, &status);
+
+ ippToPspsConfig_writeTable(this->config, fitsIn, this->fitsOut, 1, "FrameMeta", true);
+
+ // FrameMeta values
+ fits_write_col(this->fitsOut, TLONG, FRAMEMETA_FRAMEID, 1, 1, 1, &this->expId, &status);
+
+ int8_t filterID = -1;
+ ippToPspsConfig_getFilterId(this->config, filterType, &filterID);
+ fits_write_col(this->fitsOut, TBYTE, FRAMEMETA_FILTERID, 1, 1, 1, &filterID, &status);
+
+ // stuff to keep from psf.hdr header
+ double obsTime = 0.0;
+ int sourceId = -1;
+ int imageId = -1;
+ float fwhmMaj;
+ float fwhmMin;
+ float momentMaj;
+ float momentMin;
+ float psfFwhm;
+ float momentFwhm;
+ long pspsImageId = -1;
+
+ // DVO variables
+ SkyList *skylist = NULL;
+ dvoDetection *dvoDetections = NULL;
+ int maxDvoDetId = -1;
+ int numDvoDetections = -1;
+ Image *pImage = NULL;
+
+ // stuff for detections table
+ time_t now;
+ struct tm *ts;
+ char timeStr[12];
+ now = time(NULL);
+ ts = gmtime(&now); // TODO should be UTC
+ strftime(timeStr, sizeof(timeStr), "%Y-%m-%d", ts);
+
+ uint32_t totalDetections = 0;
+ uint32_t s,d;
+ uint32_t invalidDvoRows;
+ uint32_t smfJumps;
+ uint32_t unmatched;
+
+ long longnull = -999;
+ float floatnull = -999.0;
+ int anynull = 0;
+
+ char ccdNumber[3];
+ char extensionName[15];
+
+ long ippIDet[MAXDETECT];
+ float instMag[MAXDETECT];
+ float instMagErr[MAXDETECT];
+ float peakMag[MAXDETECT];
+ long objID[MAXDETECT];
+ long detectID[MAXDETECT];
+ long ippObjID[MAXDETECT];
+ long ippDetectID[MAXDETECT];
+ long imageID[MAXDETECT];
+ float obsTimes[MAXDETECT];
+ float instFlux[MAXDETECT];
+ float instFluxErr[MAXDETECT];
+ float peakFlux[MAXDETECT];
+ int8_t filterIDs[MAXDETECT];
+
+ char** assocDate = (char**)calloc(MAXDETECT, sizeof(char**));
+ for (uint32_t i=0; idvoConfig, sourceId, imageId, &pImage);
+ if (skylist == NULL) {
+ psError(PS_ERR_IO, false, "DVO: can't find SkyList for sourceId='%d' imageId='%d' (CCD = XY%s): skipping\n", sourceId, imageId, ccdNumber);
+ continue;
+ }
+
+ // create unique int from 'frameID' (aka exposure ID) and ccd number
+ pspsImageId = (this->expId*100) + pImage->ccdnum;
+
+ // now get DVO detections
+ dvoDetections = NULL;
+ numDvoDetections = dvoGetDetections(skylist, pImage->imageID, &dvoDetections, &maxDvoDetId);
+
+ // TODO check nDet < MAXDETECT
+
+ // create ImageMeta
+ ippToPspsConfig_writeTable(this->config, fitsIn, this->fitsOut, 1, "ImageMeta", true);
+ psfFwhm = (fwhmMaj+fwhmMin)/2;
+ momentFwhm = (momentMaj+momentMin)/2;
+ fits_write_col(this->fitsOut, TLONG, IMAGEMETA_PHOTOCALID, 1, 1, 1, &pImage->photcode, &status);
+ fits_write_col(this->fitsOut, TBYTE, IMAGEMETA_FILTERID, 1, 1, 1, &filterID, &status);
+ fits_write_col(this->fitsOut, TFLOAT, IMAGEMETA_PHOTOSCAT, 1, 1, 1, &zptObs, &status);
+ fits_write_col(this->fitsOut, TFLOAT, IMAGEMETA_PSFFWHM, 1, 1, 1, &psfFwhm, &status);
+ fits_write_col(this->fitsOut, TFLOAT, IMAGEMETA_PHOTOZERO, 1, 1, 1, &zptObs, &status);
+
+ // now move BACK to detections table in smf
+ sprintf(extensionName, "XY%s.psf", ccdNumber);
+ status=0;
+ if (fits_movnam_hdu(fitsIn, BINARY_TBL, extensionName, 0, &status)) {
+ psError(PS_ERR_IO, false, "Can't move to extension: %s skipping\n", extensionName);
+ continue;
+ }
+
+ // keep a running count of 'images' we find in order to write total to FrameMeta at the end
+ nOta++;
+
+ long nDet = 0;
+ if (fits_get_num_rows(fitsIn, &nDet, &status)) {
+ fits_report_error(stderr, status);
+ }
+
+ // loop round detections to populate some extra stuff
+ s = d = totalDetections = invalidDvoRows = smfJumps = unmatched = 0;
+
+ anynull = 0;
+ fits_read_col(fitsIn, TLONG, 1, 1, 1, nDet, &longnull, ippIDet, &anynull, &status);
+ fits_read_col(fitsIn, TFLOAT, 8, 1, 1, nDet, &floatnull, instMag, &anynull, &status);
+ fits_read_col(fitsIn, TFLOAT, 9, 1, 1, nDet, &floatnull, instMagErr, &anynull, &status);
+ fits_read_col(fitsIn, TFLOAT, 12, 1, 1, nDet, &floatnull, peakMag, &anynull, &status);
+
+ // DVO detections are ordered by IPP_IDET, which increments from 0 in SMF table
+ // both DVO detections and smf detections should increment. however, there may be gaps in smf and 'invalid' rows in dvo
+ numOfDuplicates = 0;
+ numInvalidFlux = 0;
+ for (long s = 0; s0 && (ippIDet[s] != (ippIDet[s-1] + 1))) smfJumps++;
+
+ while (ippIDet[s] > dvoDetections[d].meas.detID && dvoDetections[d].meas.detID <= maxDvoDetId) {
+
+ d++;
+ if (!dvoDetections[d].valid) invalidDvoRows++;
+ }
+
+ if (ippIDet[s] == dvoDetections[d].meas.detID ) {
+
+ // TODO HACK duplicates bug - don't include any duplicates
+ thisObjId = dvoDetections[d].ave.extID;
+ for (i=0; i maxObjID) maxObjID = objID[s];
+ if (objID[s] < minObjID) minObjID = objID[s];
+
+ totalDetections++;
+ }
+ // if we get here then we cannot find this smf detection in DVO
+ else {
+
+ unmatched++;
+ continue;
+ }
+
+ if (dvoDetections[d].meas.detID == maxDvoDetId) break;
+ }
+
+ // write number of rows (detections) to ImageMeta
+ numDetectionsOut = totalDetections - numOfDuplicates - numInvalidFlux;
+ fits_write_col(this->fitsOut, TLONG, IMAGEMETA_NDETECT, 1, 1, 1, &numDetectionsOut, &status);
+
+ // detections
+ ippToPspsConfig_writeTable(this->config, fitsIn, this->fitsOut, nDet, "Detection", false);
+ fits_write_col(this->fitsOut, TLONGLONG, DETECTION_OBJID, 1, 1, nDet, objID, &status);
+ fits_write_col(this->fitsOut, TLONGLONG, DETECTION_DETECTID, 1, 1, nDet, detectID, &status);
+ fits_write_col(this->fitsOut, TLONGLONG, DETECTION_IPPOBJID, 1, 1, nDet, ippObjID, &status);
+ fits_write_col(this->fitsOut, TLONGLONG, DETECTION_IPPDETECTID, 1, 1, nDet, ippDetectID, &status);
+ fits_write_col(this->fitsOut, TBYTE, DETECTION_FILTERID, 1, 1, nDet, filterIDs, &status);
+ fits_write_col(this->fitsOut, TFLOAT, DETECTION_INSTFLUX, 1, 1, nDet, instFlux, &status);
+ fits_write_col(this->fitsOut, TFLOAT, DETECTION_INSTFLUXERR, 1, 1, nDet, instFluxErr, &status);
+ fits_write_col(this->fitsOut, TFLOAT, DETECTION_PEAKADU, 1, 1, nDet, peakFlux, &status);
+ if (numOfDuplicates||numInvalidFlux) fits_delete_rowlist(this->fitsOut, removeList, numOfDuplicates+numInvalidFlux, &status);
+
+
+ psLogMsg("ippToPsps", PS_LOG_INFO,
+ "+-----------+---------+----------+------------------+---------------+--------------+\n"
+ "| Extension | Rows in | Rows out | Missing from DVO | Duplicate IDs | Invalid Flux |\n"
+ "| %s | %5ld | %5ld | %5d | %5d | %5d |\n",
+ extensionName, nDet, numDetectionsOut, unmatched, numOfDuplicates, numInvalidFlux);
+
+ if (dvoDetections!= NULL) dvoFree(dvoDetections);
+ SkyListFree(skylist);
+ }
+ }
+
+ for (uint32_t i=0; ifitsOut, BINARY_TBL, "FrameMeta", 0, &status))
+ psError(PS_ERR_IO, false, "Can't move back to FrameMeta extension to write number of OTAs\n");
+ else
+ fits_write_col(this->fitsOut, TSHORT, FRAMEMETA_NOTA, 1, 1, 1, &nOta, &status);
+
+ status=0;
+ if (fits_close_file(fitsIn, &status)) fits_report_error(stderr, status);
+
+ // write results
+ if (this->resultsFile) {
+
+ if (fprintf(this->resultsFile, "%ld\n", minObjID) < 1 )
+ psError(PS_ERR_IO, false, "Unable to write min Object ID to'%s'\n", this->resultsPath);
+ if (fprintf(this->resultsFile, "%ld\n", maxObjID) < 1 )
+ psError(PS_ERR_IO, false, "Unable to write max Object ID to'%s'\n", this->resultsPath);
+ }
+
+ psLogMsg("ippToPsps", PS_LOG_INFO, "Data written for a total of %d chips/OTAs", nOta);
+
+ return PS_EXIT_SUCCESS;
+}
+
+