Index: /branches/ppTranslate/Makefile.am =================================================================== --- /branches/ppTranslate/Makefile.am (revision 28622) +++ /branches/ppTranslate/Makefile.am (revision 28622) @@ -0,0 +1,5 @@ +SUBDIRS = src + +CLEANFILES = *.pyc *~ core core.* + +ACLOCAL_AMFLAGS = -I m4 Index: /branches/ppTranslate/autogen.sh =================================================================== --- /branches/ppTranslate/autogen.sh (revision 28622) +++ /branches/ppTranslate/autogen.sh (revision 28622) @@ -0,0 +1,103 @@ +#!/bin/sh +# Run this to generate all the initial makefiles, etc. + +srcdir=`dirname $0` +test -z "$srcdir" && srcdir=. + +ORIGDIR=`pwd` +cd $srcdir + +PROJECT=ppTranslate +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal -I `pwd`/m4 $ACLOCAL_FLAGS" +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOLIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOLIZE failed" +$ACLOCAL || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /branches/ppTranslate/configure.ac =================================================================== --- /branches/ppTranslate/configure.ac (revision 28622) +++ /branches/ppTranslate/configure.ac (revision 28622) @@ -0,0 +1,44 @@ +dnl Process this file with autoconf to produce a configure script. +AC_PREREQ(2.61) + +AC_INIT([ppMops], [0.1.1], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) +AC_CONFIG_MACRO_DIR([m4]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +IPP_STDCFLAGS + +AC_LANG(C) +AC_GNU_SOURCE +AC_PROG_CC_C99 +AC_PROG_INSTALL +AC_PROG_LIBTOOL +AC_SYS_LARGEFILE + +PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) +PKG_CHECK_MODULES([PSMODULES], [psmodules >= 1.0.0]) +CHECK_ZLIB + +AC_PATH_PROG([ERRORCODES], [psParseErrorCodes], [missing]) +if test "$ERRORCODES" = "missing" ; then + AC_MSG_ERROR([psParseErrorCodes is required]) +fi + +dnl Set CFLAGS for build +IPP_STDOPTS +CFLAGS="${CFLAGS} -Wall -Werror" + +IPP_VERSION + +AC_SUBST([PPMOPS_CFLAGS]) +AC_SUBST([PPMOPS_LIBS]) + +AC_CONFIG_FILES([ + Makefile + src/Makefile +]) + +AC_OUTPUT Index: /branches/ppTranslate/doc/MOPS_ICD_lite.txt =================================================================== --- /branches/ppTranslate/doc/MOPS_ICD_lite.txt (revision 28622) +++ /branches/ppTranslate/doc/MOPS_ICD_lite.txt (revision 28622) @@ -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/ppTranslate/m4/check_zlib.m4 =================================================================== --- /branches/ppTranslate/m4/check_zlib.m4 (revision 28622) +++ /branches/ppTranslate/m4/check_zlib.m4 (revision 28622) @@ -0,0 +1,87 @@ +dnl @synopsis CHECK_ZLIB() +dnl +dnl This macro searches for an installed zlib library. If nothing was +dnl specified when calling configure, it searches first in /usr/local +dnl and then in /usr. If the --with-zlib=DIR is specified, it will try +dnl to find it in DIR/include/zlib.h and DIR/lib/libz.a. If +dnl --without-zlib is specified, the library is not searched at all. +dnl +dnl If either the header file (zlib.h) or the library (libz) is not +dnl found, the configuration exits on error, asking for a valid zlib +dnl installation directory or --without-zlib. +dnl +dnl The macro defines the symbol HAVE_LIBZ if the library is found. You +dnl should use autoheader to include a definition for this symbol in a +dnl config.h file. Sample usage in a C/C++ source is as follows: +dnl +dnl #ifdef HAVE_LIBZ +dnl #include +dnl #endif /* HAVE_LIBZ */ +dnl +dnl @category InstalledPackages +dnl @author Loic Dachary +dnl @version 2004-09-20 +dnl @license GPLWithACException + +AC_DEFUN([CHECK_ZLIB], +# +# Handle user hints +# +[AC_MSG_CHECKING(if zlib is wanted) +AC_ARG_WITH(zlib, +[ --with-zlib=DIR root directory path of zlib installation [defaults to + /usr/local or /usr if not found in /usr/local] + --without-zlib to disable zlib usage completely], +[if test "$withval" != no ; then + AC_MSG_RESULT(yes) + if test -d "$withval" + then + ZLIB_HOME="$withval" + else + AC_MSG_WARN([Sorry, $withval does not exist, checking usual places]) + fi +else + AC_MSG_RESULT(no) +fi]) + +ZLIB_HOME=/usr/local +if test ! -f "${ZLIB_HOME}/include/zlib.h" +then + ZLIB_HOME=/usr +fi + +# +# Locate zlib, if wanted +# +if test -n "${ZLIB_HOME}" +then + ZLIB_OLD_LDFLAGS=$LDFLAGS + ZLIB_OLD_CPPFLAGS=$LDFLAGS + LDFLAGS="$LDFLAGS -L${ZLIB_HOME}/lib" + CPPFLAGS="$CPPFLAGS -I${ZLIB_HOME}/include" + AC_LANG_SAVE + AC_LANG_C + AC_CHECK_LIB(z, inflateEnd, [zlib_cv_libz=yes], [zlib_cv_libz=no]) + AC_CHECK_HEADER(zlib.h, [zlib_cv_zlib_h=yes], [zlib_cv_zlib_h=no]) + AC_LANG_RESTORE + if test "$zlib_cv_libz" = "yes" -a "$zlib_cv_zlib_h" = "yes" + then + # + # If both library and header were found, use them + # + AC_CHECK_LIB(z, inflateEnd) + AC_MSG_CHECKING(zlib in ${ZLIB_HOME}) + AC_MSG_RESULT(ok) + else + # + # If either header or library was not found, revert and bomb + # + AC_MSG_CHECKING(zlib in ${ZLIB_HOME}) + LDFLAGS="$ZLIB_OLD_LDFLAGS" + CPPFLAGS="$ZLIB_OLD_CPPFLAGS" + AC_MSG_RESULT(failed) + AC_MSG_ERROR(either specify a valid zlib installation with --with-zlib=DIR or disable zlib usage with --without-zlib) + fi +fi + +]) Index: /branches/ppTranslate/src/Makefile.am =================================================================== --- /branches/ppTranslate/src/Makefile.am (revision 28622) +++ /branches/ppTranslate/src/Makefile.am (revision 28622) @@ -0,0 +1,56 @@ +bin_PROGRAMS = ppMops ppMonet + +if HAVE_SVNVERSION +PPTRANSLATE_VERSION=`$(SVNVERSION) ..` +else +PPTRANSLATE_VERSION="UNKNOWN" +endif + +if HAVE_SVN +PPTRANSLATE_BRANCH=`$(SVN) info .. | $(SED) -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'` +PPTRANSLATE_SOURCE=`$(SVN) info | $(SED) -n -e 's/Repository UUID: // p'` +else +PPTRANSLATE_BRANCH="UNKNOWN" +PPTRANSLATE_SOURCE="UNKNOWN" +endif + +# Force recompilation of ppMopsVersion.c, since it gets the version information +ppTranslateVersion.c: ppTranslateVersionDefinitions.h +ppTranslateVersionDefinitions.h: ppTranslateVersionDefinitions.h.in FORCE + -$(RM) ppTranslateVersionDefinitions.h + $(SED) -e "s|@PPTRANSLATE_VERSION@|\"$(PPTRANSLATE_VERSION)\"|" -e "s|@PPTRANSLATE_BRANCH@|\"$(PPTRANSLATE_BRANCH)\"|" -e "s|@PPTRANSLATE_SOURCE@|\"$(PPTRANSLATE_SOURCE)\"|" ppTranslateVersionDefinitions.h.in > ppTranslateVersionDefinitions.h +FORCE: ; + +ppMops_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULES_CFLAGS) $(PPMOPS_CFLAGS) +ppMops_LDFLAGS = $(PSLIB_LIBS) $(PPMOPS_LIBS) + +ppMops_SOURCES = \ + ppMops.c \ + ppTranslateVersion.c \ + ppMopsArguments.c \ + ppMopsDetections.c \ + ppMopsRead.c \ + ppMopsWrite.c \ + ppMopsMerge.c + +ppMonet_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULES_CFLAGS) $(PPMONET_CFLAGS) +ppMonet_LDFLAGS = $(PSLIB_LIBS) $(PPMONET_LIBS) + +ppMonet_SOURCES = \ + ppMonet.c \ + ppTranslateVersion.c \ + ppMonetArguments.c + + +noinst_HEADERS = \ + ppTranslateVersion.h \ + ppMops.h \ + ppMonet.h + + +clean-local: + -rm -f TAGS + +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /branches/ppTranslate/src/ppMonet.c =================================================================== --- /branches/ppTranslate/src/ppMonet.c (revision 28622) +++ /branches/ppTranslate/src/ppMonet.c (revision 28622) @@ -0,0 +1,165 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppTranslateVersion.h" +#include "ppMonet.h" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + + ppMonetArguments *args = ppMonetArgumentsParse(argc, argv); // Parsed arguments + if (!args) { + psErrorStackPrint(stderr, "Error parsing arguments"); + exit(PS_EXIT_CONFIG_ERROR); + } + + psTrace("ppMonet.read", 1, "Reading input detections\n"); + + psFits *fits = psFitsOpen(args->input, "r"); // FITS file + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open input %s", args->input); + exit(PS_EXIT_SYS_ERROR); + } + + // Set up output file + FILE *file = fopen(args->output, "w"); // File handle + if (!file) { + psErrorStackPrint(stderr, "Error opening file %s", args->output); + exit(PS_EXIT_SYS_ERROR); + } + gzFile gz = gzdopen(fileno(file), "w"); // Gzip file handle + gzprintf(gz, "# Pan-STARRS IPP detections for Dave Monet\n"); + gzprintf(gz, "#\n"); + gzprintf(gz, "# Translated from %s\n", args->input); + psString source = ppTranslateSource(), version = ppTranslateVersion(); + gzprintf(gz, "# S/W source = %s\n", source); + gzprintf(gz, "# S/W version = %s\n", version); + gzprintf(gz, "#\n"); + psFree(source); + psFree(version); + gzprintf(gz, "# exp_name = %s\n", args->exp_name); + gzprintf(gz, "# exp_id = %" PRId64 "\n", args->exp_id); + gzprintf(gz, "# chip_id = %" PRId64 "\n", args->chip_id); + gzprintf(gz, "# cam_id = %" PRId64 "\n", args->cam_id); + gzprintf(gz, "# ZP = %f\n", args->zp); + gzprintf(gz, "# ZP error = %f\n", args->zpErr); + gzprintf(gz, "# Astrometry rms = %f\n", args->rmsAstrom); + + psFree(args); + + psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header"); + exit(PS_EXIT_SYS_ERROR); + } + + double raBoresight = psMetadataLookupF64(NULL, header, "RA"); + double decBoresight = psMetadataLookupF64(NULL, header, "DEC"); + psString filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FILTERID")); + float airmass = psMetadataLookupF32(NULL, header, "AIRMASS"); + float exptime = psMetadataLookupF32(NULL, header, "EXPTIME"); + double posangle = RAD_TO_DEG(psMetadataLookupF64(NULL, header, "POSANGLE")); + double rotangle = psMetadataLookupF64(NULL, header, "ROTANGLE"); + double alt = psMetadataLookupF64(NULL, header, "ALT"); + double az = psMetadataLookupF64(NULL, header, "AZ"); + double mjd = psMetadataLookupF64(NULL, header, "MJD-OBS"); + + gzprintf(gz, "# RA = %.8lf\n", raBoresight); + gzprintf(gz, "# Dec = %.8lf\n", decBoresight); + gzprintf(gz, "# Filter = %s\n", filter); + gzprintf(gz, "# Airmass = %f\n", airmass); + gzprintf(gz, "# ExpTime = %f\n", exptime); + gzprintf(gz, "# PA = %lf\n", posangle); + gzprintf(gz, "# Rotator = %lf\n", rotangle); + gzprintf(gz, "# Alt = %.8lf\n", alt); + gzprintf(gz, "# Az = %.8lf\n", az); + gzprintf(gz, "# MJD = %.6lf\n", mjd); + gzprintf(gz, "#\n"); + gzprintf(gz, "# chip , x , y , RA , Dec , Mag , MagErr , Width , Flags\n"); + + psFree(header); + + int numHDU = psFitsGetSize(fits); // Number of HDUs + for (int i = 1; i < numHDU; i++) { + if (!psFitsMoveExtNum(fits, i, false)) { + psError(PS_ERR_IO, false, "Unable to move to HDU %d", i); + exit(PS_EXIT_SYS_ERROR); + } + + if (psFitsGetExtType(fits) != PS_FITS_TYPE_BINARY_TABLE) { + psTrace("ppMonet", 1, "Skipping extension %d: not a binary table", i); + continue; + } + + psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header %d\n", i); + exit(PS_EXIT_SYS_ERROR); + return false; + } + + const char *extname = psMetadataLookupStr(NULL, header, "EXTNAME"); + const char *exttype = extname + strlen(extname) - strlen(EXT_TYPE); + if (strcmp(exttype, EXT_TYPE) != 0) { + psTrace("ppMonet", 1, "Skipping extension %d (%s): not correct type", i, extname); + psFree(header); + continue; + } + + long size = psFitsTableSize(fits); // Size of table + if (size <= 0) { + psErrorStackPrint(stderr, "Unable to determine size of table %d", i); + psFree(header); + exit(PS_EXIT_SYS_ERROR); + } + + psString chipName = psStringNCopy(extname, strlen(extname) - strlen(EXT_TYPE)); // Name of chip + psFree(header); + + psTrace("ppMonet.read", 3, "Reading %ld rows from %s\n", size, extname); + psArray *table = psFitsReadTable(fits); // Table of interest + if (!table) { + psError(PS_ERR_IO, false, "Unable to read table %d", i); + psFree(chipName); + exit(PS_EXIT_SYS_ERROR); + } + + psTrace("ppMonet.read", 3, "Writing %ld rows for %s\n", size, chipName); + for (long j = 0; j < size; j++) { + psMetadata *row = table->data[j]; // Row of interest + + float x = psMetadataLookupF32(NULL, row, "X_PSF"); + float y = psMetadataLookupF32(NULL, row, "Y_PSF"); + double ra = psMetadataLookupF64(NULL, row, "RA_PSF"); + double dec = psMetadataLookupF64(NULL, row, "DEC_PSF"); + float mag = psMetadataLookupF32(NULL, row, "PSF_INST_MAG"); + float magErr = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG"); + float major = psMetadataLookupF32(NULL, row, "PSF_MAJOR"); + float minor = psMetadataLookupF32(NULL, row, "PSF_MINOR"); + psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS"); + + float width = 0.5 * (major + minor); + + gzprintf(gz, "%s , %.2f , %.2f , %.8lf , %.8lf , %.3f , %.3f , %.2f , %#08x\n", + chipName, x, y, ra, dec, mag, magErr, width, flags); + } + psFree(table); + psFree(chipName); + } + gzclose(gz); + psFitsClose(fits); + psFree(filter); + + psTrace("ppMonet.read", 1, "Done reading input detections\n"); + + psLibFinalize(); + + return PS_EXIT_SUCCESS; +} Index: /branches/ppTranslate/src/ppMonet.h =================================================================== --- /branches/ppTranslate/src/ppMonet.h (revision 28622) +++ /branches/ppTranslate/src/ppMonet.h (revision 28622) @@ -0,0 +1,23 @@ +#ifndef PP_MOPS_H +#define PP_MOPS_H + +#include + +#define EXT_TYPE ".psf" // Extension type to consider + +// Configuration data +typedef struct { + psString 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 + float zp, zpErr; // Magnitude zero point and error + float rmsAstrom; // Astrometric solution RMS + psString output; // Output filename +} ppMonetArguments; + +/// Parse arguments +ppMonetArguments *ppMonetArgumentsParse(int argc, char *argv[]); + +#endif Index: /branches/ppTranslate/src/ppMonetArguments.c =================================================================== --- /branches/ppTranslate/src/ppMonetArguments.c (revision 28622) +++ /branches/ppTranslate/src/ppMonetArguments.c (revision 28622) @@ -0,0 +1,88 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppMonet.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-Monet detection translator\n\n"); + fprintf(stderr, "Usage: %s INPUT_NAME OUTPUT_NAME\n", program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +static void monetArgumentsFree(ppMonetArguments *args) +{ + psFree(args->input); + psFree(args->exp_name); + psFree(args->output); + return; +} + +ppMonetArguments *ppMonetArgumentsAlloc(void) +{ + ppMonetArguments *args = psAlloc(sizeof(ppMonetArguments)); // Data to return + psMemSetDeallocator(args, (psFreeFunc)monetArgumentsFree); + + args->input = NULL; + args->exp_name = NULL; + args->exp_id = 0; + args->chip_id = 0; + args->cam_id = 0; + args->zp = NAN; + args->zpErr = NAN; + args->rmsAstrom = NAN; + args->output = NULL; + + return args; +} + + +ppMonetArguments *ppMonetArgumentsParse(int argc, char *argv[]) +{ + assert(argv); + + psTrace("ppMonet.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); + 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); + } + + ppMonetArguments *args = ppMonetArgumentsAlloc(); // Arguments, to return + + args->input = psStringCopy(argv[1]); + 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->zp = psMetadataLookupF32(NULL, arguments, "-zp"); + args->zpErr = psMetadataLookupF32(NULL, arguments, "-zp_error"); + args->rmsAstrom = psMetadataLookupF32(NULL, arguments, "-astrom_rms"); + + psTrace("ppMonet.args", 1, "Done parsing command-line arguments\n"); + + return args; +} Index: /branches/ppTranslate/src/ppMonetMerge.c =================================================================== --- /branches/ppTranslate/src/ppMonetMerge.c (revision 28622) +++ /branches/ppTranslate/src/ppMonetMerge.c (revision 28622) @@ -0,0 +1,171 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppMonet.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 ppMonetDetections *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); +} + + +ppMonetDetections *ppMonetMerge(const psArray *detections) +{ + PS_ASSERT_ARRAY_NON_NULL(detections, NULL); + + psTrace("ppMonet.merge", 1, "Merging detections from %ld inputs\n", detections->n); + + ppMonetDetections *merged = NULL; // Merged list + int num = 1; // Number of merged files + for (int i = 0; i < detections->n; i++) { + ppMonetDetections *det = detections->data[i]; // Detections of interest + if (!det) { + psTrace("ppMonet.merge", 3, "Ignoring NULL input %d\n", i); + continue; + } else if (det->num == 0) { + psTrace("ppMonet.merge", 3, "Ignoring empty input %d\n", i); + continue; + } + num++; + if (!merged) { + psTrace("ppMonet.merge", 3, "Accepting %ld detections from input %d\n", det->num, i); + merged = psMemIncrRefCounter(det); + continue; + } + psTrace("ppMonet.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("ppMonet.merge", 9, "No matches for source %d in input %d\n", j, i); + psFree(indices); + ppMonetDetectionsCopySingle(merged, det, j); + continue; + } + psTrace("ppMonet.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("ppMonet.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; + } + ppMonetDetectionsCopySingle(merged, det, j); + } else { + psTrace("ppMonet.merge", 6, "Old sources clobber new source\n"); + } + psFree(indices); + } + + psTrace("ppMonet.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num); + + psFree(tree); + ppMonetDetectionsPurge(merged); + } + + psTrace("ppMonet.merge", 2, "%ld sources in merged detections list\n", merged->num); + + merged->seeing /= num; + + return merged; +} + Index: /branches/ppTranslate/src/ppMops.c =================================================================== --- /branches/ppTranslate/src/ppMops.c (revision 28622) +++ /branches/ppTranslate/src/ppMops.c (revision 28622) @@ -0,0 +1,232 @@ +#include +#include + +#include "ppMops.h" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + + ppMopsArguments *args = ppMopsArgumentsParse(argc, argv); // Parsed arguments + if (!args) { + psErrorStackPrint(stderr, "Error parsing arguments"); + 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 (!ppMopsPurgeDuplicates(detections)) { + psErrorStackPrint(stderr, "Unable to merge detections"); + exit(PS_EXIT_SYS_ERROR); + } + + if (!ppMopsWrite(detections, args)) { + psErrorStackPrint(stderr, "Unable to write detections"); + exit(PS_EXIT_SYS_ERROR); + } + + psFree(detections); + psFree(args); + + psLibFinalize(); + + 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: /branches/ppTranslate/src/ppMops.h =================================================================== --- /branches/ppTranslate/src/ppMops.h (revision 28622) +++ /branches/ppTranslate/src/ppMops.h (revision 28622) @@ -0,0 +1,62 @@ +#ifndef PP_MOPS_H +#define PP_MOPS_H + +#include +#include + +#define IN_EXTNAME "SkyChip.psf" // Extension name for data in input +#define OBSERVATORY_CODE "F51" // IAU Observatory Code + +// Configuration data +typedef struct { + 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 +} ppMopsArguments; + +/// Parse arguments +ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]); + +typedef struct { + psString component; // Component name + psMetadata *header; // FITS header (extension names *.hdr) + 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 of exposure mid-point + float seeing; // Seeing of exposure + int naxis1, naxis2; // Size of image + long num; // Number of detections + psMetadata *psfHeader; // FITS header (extension names *.psf) + psMetadata *table; // Columns of data + psVector *x, *y; // Image coordinates + psVector *ra, *dec; // Sky coordinates + psMetadata *deteffHeader; // Detection efficiency header (extension names *.deteff) + psMetadata *deteffTable; // Detection efficiency table +} ppMopsDetections; + +// Allocator +ppMopsDetections *ppMopsDetectionsAlloc(void); + +/// Read detections +psArray *ppMopsRead(const ppMopsArguments *args); + +/// Purge duplicate detections +bool ppMopsPurgeDuplicates(const psArray *detections); + +/// Write detections +bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args); + +#endif Index: /branches/ppTranslate/src/ppMopsArguments.c =================================================================== --- /branches/ppTranslate/src/ppMopsArguments.c (revision 28622) +++ /branches/ppTranslate/src/ppMopsArguments.c (revision 28622) @@ -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/ppTranslate/src/ppMopsDetections.c =================================================================== --- /branches/ppTranslate/src/ppMopsDetections.c (revision 28622) +++ /branches/ppTranslate/src/ppMopsDetections.c (revision 28622) @@ -0,0 +1,55 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +static void mopsDetectionsFree(ppMopsDetections *det) +{ + psFree(det->component); + psFree(det->raBoresight); + psFree(det->decBoresight); + psFree(det->filter); + psFree(det->header); + psFree(det->table); + psFree(det->x); + psFree(det->y); + psFree(det->ra); + psFree(det->dec); + psFree(det->deteffHeader); + psFree(det->deteffTable); + + return; +} + +ppMopsDetections *ppMopsDetectionsAlloc(void) +{ + ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return + psMemSetDeallocator(det, (psFreeFunc)mopsDetectionsFree); + + det->component = NULL; + 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->header = NULL; + det->table = NULL; + det->x = NULL; + det->y = NULL; + det->ra = NULL; + det->dec = NULL; + det->deteffHeader = NULL; + det->deteffTable = NULL; + + return det; +} Index: /branches/ppTranslate/src/ppMopsMerge.c =================================================================== --- /branches/ppTranslate/src/ppMopsMerge.c (revision 28622) +++ /branches/ppTranslate/src/ppMopsMerge.c (revision 28622) @@ -0,0 +1,290 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppMops.h" + +#define LEAF_SIZE 2 // 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 + +#if 0 +#undef psTrace +#define psTrace(facil, level, ...) \ + if (level <= 5) { \ + fprintf(stderr, __VA_ARGS__); \ + } +#endif + +// 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 / 2.0; + float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0; + return PS_SQR(dx) + PS_SQR(dy); +} + + +bool ppMopsPurgeDuplicates(const psArray *detections) +{ + PS_ASSERT_ARRAY_NON_NULL(detections, NULL); + + int numInputs = detections->n; // Number of inputs + psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs); + + long total = 0; // Total number of sources + psArray *duplicates = psArrayAlloc(numInputs); // Vector of duplicate bits for each input + psVector *dupNum = psVectorAlloc(numInputs, PS_TYPE_U32); // Number of duplicates for each input + psVectorInit(dupNum, 0); + for (int i = 0; i < numInputs; i++) { + ppMopsDetections *det = detections->data[i]; // Detections from + if (!det || det->num == 0) { + continue; + } + psVector *dupes = duplicates->data[i] = psVectorAlloc(det->num, PS_TYPE_U8); + psVectorInit(dupes, 0); + total += det->num; + } + psTrace("ppMops.merge", 2, "Total detections: %ld\n", total); + + psVector *raMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged RAs + psVector *decMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged Decs + psVector *sourceMerged = psVectorAlloc(total, PS_TYPE_U16); // Source image of merged sources + psVector *indexMerged = psVectorAlloc(total, PS_TYPE_U32); // Index of merged sources + long num = 0; // Number of merged sources + + const char *raBoresight = NULL, *decBoresight = NULL; // Boresight coordinates + const char *filter = NULL; // Filter name + float airmass = NAN, exptime = NAN; // Exposure details + double posangle = NAN, alt = NAN, az = NAN; // Telescope pointing + double mjd = NAN; // Time of exposure + + for (int i = 0; i < numInputs; 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; + } + + // Check exposure characteristics + if (num == 0) { + raBoresight = det->raBoresight; + decBoresight = det->decBoresight; + filter = det->filter; + airmass = det->airmass; + exptime = det->exptime; + posangle = det->posangle; + alt = det->alt; + az = det->az; + } else { + if (strcmp(raBoresight, det->raBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s", + raBoresight, det->raBoresight); + return false; + } + if (strcmp(decBoresight, det->decBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s", + decBoresight, det->decBoresight); + return false; + } + if (strcmp(filter, det->filter) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s", + filter, det->filter); + return false; + } + if (fabsf(airmass - det->airmass) > AIRMASS_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f", + airmass, det->airmass); + return false; + } + if (fabsf(exptime - det->exptime) > EXPTIME_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f", + exptime, det->exptime); + return false; + } + if (fabs(posangle - det->posangle) > POSANGLE_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f", + posangle, det->posangle); + return false; + } + if (fabs(alt - det->alt) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf", + alt, det->alt); + return false; + } + if (fabs(az - det->az) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf", + az, det->az); + return false; + } + if (fabs(mjd - det->mjd) > MJD_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf", + mjd, det->mjd); + return false; + } + } + + psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i); + memcpy(&raMerged->data.F64[num], det->ra->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + memcpy(&decMerged->data.F64[num], det->dec->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + for (long j = 0, k = num; j < det->num; j++, k++) { + sourceMerged->data.U16[k] = i; + indexMerged->data.U32[k] = j; + } + num += det->num; + } + + psTrace("ppMops.merge", 3, "Generating kd-tree from %ld sources\n", num); + psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, raMerged, decMerged); // kd tree + if (!tree) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree"); + return false; + } + + for (int i = 0; i < detections->n; i++) { + ppMopsDetections *det = detections->data[i]; // Detections of interest + if (!det) { + continue; + } else if (det->num == 0) { + continue; + } + psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i); + + psVector *dupes = duplicates->data[i]; // Duplicates list + 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); + return false; + } + psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i); + psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i); + + if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) { + // It's myself + psFree(indices); + continue; + } + + // 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 mergeIndex = indices->data.S64[k]; // Index of point in merged list + int source = sourceMerged->data.U16[mergeIndex]; // Source image + if (source == i) { + continue; + } + long index = indexMerged->data.U32[mergeIndex]; // Index in source + psVector *dupes = duplicates->data[source]; // Duplicates list + if (dupes->data.U8[index]) { + continue; + } + + float distance = mergeDistance(detections->data[source], index); // Distance to centre of image + if (distance < bestDistance) { + bestDistance = distance; + bestIndex = index; + } + } + + if (bestIndex >= 0) { + float distance = mergeDistance(det, j); // Distance to centre of image + if (bestIndex >= 0 && distance < bestDistance) { + psTrace("ppMops.merge", 6, "New source clobbers old sources\n"); + for (int k = 0; k < indices->n; k++) { + long mergeIndex = indices->data.S64[k]; // Index of point + int source = sourceMerged->data.U16[mergeIndex]; // Source image + if (source == i) { + continue; + } + long index = indexMerged->data.U32[mergeIndex]; // Index in source + psVector *dupes = duplicates->data[source]; // Duplicates list + if (!dupes->data.U8[index]) { + dupes->data.U8[index] = 0xFF; + dupNum->data.U32[source]++; + } + } + } else { + psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); + dupes->data.U8[j] = 0xFF; + dupNum->data.U32[i]++; + } + } + psFree(indices); + } + psFree(coords); + } + psFree(tree); + psFree(raMerged); + psFree(decMerged); + psFree(sourceMerged); + psFree(indexMerged); + + // Remove duplicates + for (int i = 0; i < detections->n; i++) { + ppMopsDetections *det = detections->data[i]; // Detections of interest + if (!det) { + continue; + } else if (det->num == 0) { + continue; + } + psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i); + +#define VECTOR_PURGE_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + for (int i = 0, j = 0; i < vector->n; i++) { \ + if (!dupes->data.U8[i]) { \ + if (i == j) { \ + j++; \ + continue; \ + } \ + vector->data.TYPE[j] = vector->data.TYPE[i]; \ + } \ + } \ + break; + + psVector *dupes = duplicates->data[i]; // Duplicates + psArray *table = psListToArray(det->table->list); // Table of data + for (int j = 0; j < table->n; j++) { + psMetadataItem *item = table->data[j]; // Table item + psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type); + psVector *vector = item->data.V; // Vector to purge + switch (vector->type.type) { + VECTOR_PURGE_CASE(U8); + VECTOR_PURGE_CASE(U16); + VECTOR_PURGE_CASE(U32); + VECTOR_PURGE_CASE(U64); + VECTOR_PURGE_CASE(S8); + VECTOR_PURGE_CASE(S16); + VECTOR_PURGE_CASE(S32); + VECTOR_PURGE_CASE(S64); + VECTOR_PURGE_CASE(F32); + VECTOR_PURGE_CASE(F64); + default: + psAbort("Unrecognised vector type: %x", vector->type.type); + } + } + } + + return true; +} + Index: /branches/ppTranslate/src/ppMopsRead.c =================================================================== --- /branches/ppTranslate/src/ppMopsRead.c (revision 28622) +++ /branches/ppTranslate/src/ppMopsRead.c (revision 28622) @@ -0,0 +1,119 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#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++) { + const char *name = inNames->data[i]; // File name + 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; + } + + 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 = detections->data[i] = ppMopsDetectionsAlloc(); + det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension + det->header = header; + det->num = 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, "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")); + + det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); + det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); + + det->psfHeader = psFitsReadHeader(NULL, fits); + if (!det->psfHeader) { + psError(psErrorCodeLast(), false, "Unable to read header %d", i); + return false; + } + psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest + if (!table) { + psError(psErrorCodeLast(), false, "Unable to read table %d", i); + return false; + } + + psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); + psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); + + 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 false; + } + + psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); + + if (!psFitsMoveExtName(fits, "SkyChip.deteff")) { + psWarning("No detection efficiencies included in %s", det->component); + psErrorStackPrint(stderr, "No detection efficiencies included in %s", det->component); + psErrorClear(); + continue; + } + + det->deteffHeader = psFitsReadHeader(NULL, fits); + if (!det->deteffHeader) { + psWarning("Unable to read detection efficiency header in %s", det->component); + psErrorClear(); + continue; + } + det->deteffTable = psFitsReadTableAllColumns(fits); + if (!det->deteffTable) { + psFree(det->deteffHeader); + det->deteffHeader = NULL; + psWarning("Unable to read detection efficiency table in %s", det->component); + psErrorClear(); + continue; + } + psTrace("ppMops.read", 2, "Read detection efficiency from %s\n", det->component); + psFitsClose(fits); + } + + psTrace("ppMops.read", 1, "Done reading input detections\n"); + + return detections; +} Index: /branches/ppTranslate/src/ppMopsWrite.c =================================================================== --- /branches/ppTranslate/src/ppMopsWrite.c (revision 28622) +++ /branches/ppTranslate/src/ppMopsWrite.c (revision 28622) @@ -0,0 +1,86 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" +#include "ppTranslateVersion.h" + +bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args) +{ + psTrace("ppMops.write", 1, "Writing %ld extensions to %s", detections->n, args->output); + + psFits *fits = psFitsOpen(args->output, "w"); // FITS file + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open output file."); + return false; + } + + // Primary header + { + psMetadata *header = psMetadataAlloc(); // Header to write + ppTranslateVersionHeader(header); + 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); + 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); + psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); + + ppMopsDetections *det = detections->data[0]; // Representative set of detections + 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); + psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter); + psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass at boresight", det->airmass); + psMetadataAddF32(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime); + psMetadataAddF64(header, PS_LIST_TAIL, "POSANG", 0, "Position angle", det->posangle); + psMetadataAddF64(header, PS_LIST_TAIL, "ALT", 0, "Altitude of boresight", det->alt); + psMetadataAddF64(header, PS_LIST_TAIL, "AZ", 0, "Azimuth of boresight", det->az); + psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of observation", det->mjd); + + if (!psFitsWriteBlank(fits, header, NULL)) { + psError(psErrorCodeLast(), false, "Unable to write primary header"); + return false; + } + psFree(header); + } + + for (int i = 0; i < detections->n; i++) { + ppMopsDetections *det = detections->data[i]; // Detections for extension + psTrace("ppMops.write", 1, "Writing extension %d to %s", i, args->output); + psString hdrName = NULL, psfName = NULL, deteffName = NULL; + psStringAppend(&hdrName, "%s.hdr", det->component); + psStringAppend(&psfName, "%s.psf", det->component); + psStringAppend(&deteffName, "%s.deteff", det->component); + + if (!psFitsWriteBlank(fits, det->header, hdrName)) { + psError(psErrorCodeLast(), false, "Unable to write header %d", i); + return false; + } + if (!psFitsWriteTableAllColumns(fits, det->psfHeader, det->table, psfName)) { + psError(psErrorCodeLast(), false, "Unable to write table %d", i); + return false; + } + if (det->deteffHeader && det->deteffTable && + !psFitsWriteTableAllColumns(fits, det->deteffHeader, det->deteffTable, deteffName)) { + psError(psErrorCodeLast(), false, "Unable to write detection efficiency %d", i); + return false; + } + psFree(hdrName); + psFree(psfName); + psFree(deteffName); + } + psFitsClose(fits); + + psTrace("ppMops.write", 1, "Done writing to %s", args->output); + + return true; +} Index: /branches/ppTranslate/src/ppTranslateVersion.c =================================================================== --- /branches/ppTranslate/src/ppTranslateVersion.c (revision 28622) +++ /branches/ppTranslate/src/ppTranslateVersion.c (revision 28622) @@ -0,0 +1,101 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppTranslateVersionDefinitions.h" + +#ifndef PPTRANSLATE_VERSION +#error "PPTRANSLATE_VERSION is not set" +#endif +#ifndef PPTRANSLATE_BRANCH +#error "PPTRANSLATE_BRANCH is not set" +#endif +#ifndef PPTRANSLATE_SOURCE +#error "PPTRANSLATE_SOURCE is not set" +#endif + +psString ppTranslateVersion(void) +{ + char *value = NULL; + psStringAppend(&value, "%s@%s", PPTRANSLATE_BRANCH, PPTRANSLATE_VERSION); + return value; +} + +psString ppTranslateSource(void) +{ + return psStringCopy(PPTRANSLATE_SOURCE); +} + +psString ppTranslateVersionLong(void) +{ + psString version = ppTranslateVersion(); // Version, to return + psString source = ppTranslateSource(); // Source + + psStringPrepend(&version, "ppTranslate "); + psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__); + psFree(source); + +#ifdef __OPTIMIZE__ + psStringAppend(&version, " optimised"); +#else + psStringAppend(&version, " unoptimised"); +#endif + + return version; +}; + + +bool ppTranslateVersionHeader(psMetadata *header) +{ + PS_ASSERT_METADATA_NON_NULL(header, false); + + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psString history = NULL; // History string + psStringAppend(&history, "ppTranslate at %s", timeString); + psFree(timeString); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history); + psFree(history); + + psLibVersionHeader(header); + + psString version = ppTranslateVersion(); // Software version + psString source = ppTranslateSource(); // Software source + + psMetadataAddStr(header, PS_LIST_TAIL, "TRANSL_V", PS_META_REPLACE, NULL, PPTRANSLATE_VERSION); + + psStringPrepend(&version, "ppTranslate version: "); + psStringPrepend(&source, "ppTranslate source: "); + + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source); + + psFree(version); + psFree(source); + + return true; +} + +void ppTranslateVersionPrint(void) +{ + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psLogMsg("ppTranslate", PS_LOG_INFO, "ppTranslate at %s", timeString); + psFree(timeString); + + psString pslib = psLibVersionLong();// psLib version + psString ppTranslate = ppTranslateVersionLong(); // ppTranslate version + + psLogMsg("ppTranslate", PS_LOG_INFO, "%s", pslib); + psLogMsg("ppTranslate", PS_LOG_INFO, "%s", ppTranslate); + + psFree(pslib); + psFree(ppTranslate); + + return; +} Index: /branches/ppTranslate/src/ppTranslateVersion.h =================================================================== --- /branches/ppTranslate/src/ppTranslateVersion.h (revision 28622) +++ /branches/ppTranslate/src/ppTranslateVersion.h (revision 28622) @@ -0,0 +1,21 @@ +#ifndef PP_TRANSLATE_VERSION +#define PP_TRANSLATE_VERSION + +#include + +/// Return version +psString ppTranslateVersion(void); + +/// Return source +psString ppTranslateSource(void); + +/// Return detailed version information +psString ppTranslateVersionLong(void); + +/// Put version into header +bool ppTranslateVersionHeader(psMetadata *header); + +/// Print version information +void ppTranslateVersionPrint(void); + +#endif Index: /branches/ppTranslate/src/ppTranslateVersionDefinitions.h.in =================================================================== --- /branches/ppTranslate/src/ppTranslateVersionDefinitions.h.in (revision 28622) +++ /branches/ppTranslate/src/ppTranslateVersionDefinitions.h.in (revision 28622) @@ -0,0 +1,8 @@ +#ifndef PPTRANSLATE_VERSION_DEFINITIONS_H +#define PPTRANSLATE_VERSION_DEFINITIONS_H + +#define PPTRANSLATE_VERSION @PPTRANSLATE_VERSION@ // SVN version +#define PPTRANSLATE_BRANCH @PPTRANSLATE_BRANCH@ // SVN branch +#define PPTRANSLATE_SOURCE @PPTRANSLATE_SOURCE@ // SVN source + +#endif