Index: trunk/ippToPsps/src/Makefile.am
===================================================================
--- trunk/ippToPsps/src/Makefile.am	(revision 27806)
+++ 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 <limits.h>
+#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; i<MAXDETECT;i++) assocDate[i] = (char*)calloc(20,sizeof(char));
+
+    // some stuff is the same for all detections so we can populate here
+    for (long s = 0; s<MAXDETECT; s++) {
+
+        filterIDs[s] = filterID;
+        strcpy(assocDate[s], timeStr);
+    }
+
+    long maxObjID = LONG_MIN; 
+    long minObjID = LONG_MAX;
+    short nOta = 0;
+    long i;
+    bool isDuplicate;
+    uint32_t numOfDuplicates;
+    uint32_t numInvalidFlux;
+    long numDetectionsOut;
+    long removeList[MAXDETECT];
+    long thisObjId;
+
+    // loop round all 60 chips
+    for (int x=3; x<4; x++) {
+        for (int y=3; y<4; y++) {
+
+            // dodge the corners
+            if (x==0 && y==0) continue;
+            if (x==0 && y==7) continue;
+            if (x==7 && y==0) continue;
+            if (x==7 && y==7) continue;
+
+            sprintf(ccdNumber, "%d%d", x, y);
+
+            // check we can move 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 this chip\n", extensionName);
+                continue;
+            }
+
+            // move to header extension
+            sprintf(extensionName, "XY%s.hdr", ccdNumber);
+            status=0;
+            if (fits_movnam_hdu(fitsIn, IMAGE_HDU, extensionName, 0, &status)) {
+                psError(PS_ERR_IO, false, "Can't move to extension: %s skipping this chip\n", extensionName);
+                continue;
+            }
+
+            // stuff to save from psf.hdr
+            status=0; fits_read_key(fitsIn, TFLOAT, "FWHM_MAJ", &fwhmMaj, NULL, &status);
+            status=0; fits_read_key(fitsIn, TFLOAT, "FWHM_MIN", &fwhmMin, NULL, &status);
+            status=0; fits_read_key(fitsIn, TFLOAT, "IQ_FW1", &momentMaj, NULL, &status);
+            status=0; fits_read_key(fitsIn, TFLOAT, "IQ_FW2", &momentMin, NULL, &status);
+            status=0; fits_read_key(fitsIn, TDOUBLE, "MJD-OBS", &obsTime, NULL, &status);
+
+            status=0; fits_read_key(fitsIn, TLONG, "IMAGEID", &imageId, NULL, &status);
+            status=0; fits_read_key(fitsIn, TLONG, "SOURCEID", &sourceId, NULL, &status);
+
+            // access DVO database
+            skylist = dvoSkyListByExternID(this->dvoConfig, 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; s<nDet; s++) {
+
+                isDuplicate = false;
+
+                // count jumps in smf file
+                if (s>0 && (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<totalDetections; i++) {
+                        if (thisObjId == objID[i]) {
+                            removeList[numOfDuplicates+numInvalidFlux] = s+1;
+                            numOfDuplicates++;
+                            isDuplicate=true;
+                            break;
+                        }
+                    }
+                    // populate some arrays to be shoved in table in bulk later
+                    objID[s] = thisObjId;
+                    detectID[s] = dvoDetections[d].meas.extID;
+                    ippObjID[s] = (uint64_t)dvoDetections[d].ave.catID*1000000000 + (uint64_t)dvoDetections[d].ave.objID;
+                    ippDetectID[s] = dvoDetections[d].meas.detID;
+                    imageID[s] = pspsImageId;
+                    obsTimes[s] = obsTime;
+
+                    ippToPsps_getFlux(zeroPoint, exposureTime, instMagErr[s], &instFluxErr[s]);
+
+                    // check for invalid flux values (if not already labelled as a duplicate)
+                    if ((!ippToPsps_getFlux(zeroPoint, exposureTime, peakMag[s], &peakFlux[s]) ||
+                                !ippToPsps_getFlux(zeroPoint, exposureTime, instMag[s], &instFlux[s])) && !isDuplicate) {
+                        removeList[numOfDuplicates+numInvalidFlux] = s+1;
+                        numInvalidFlux++;
+                    }
+
+                    // store max/min objID
+                    if (objID[s] > 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; i<MAXDETECT;i++) free(assocDate[i]);
+    free(assocDate);
+
+    // write number of images we have found into FrameMeta table
+    status=0; 
+    if (fits_movnam_hdu(this->fitsOut, 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;
+}
+
+