Index: /branches/eam_branch_20080511/ppSim/.cvsignore =================================================================== --- /branches/eam_branch_20080511/ppSim/.cvsignore (revision 17629) +++ /branches/eam_branch_20080511/ppSim/.cvsignore (revision 17629) @@ -0,0 +1,15 @@ +aclocal.m4 +autom4te.cache +config.guess +config.log +config.status +config.sub +configure +depcomp +install-sh +libtool +ltmain.sh +Makefile +Makefile.in +missing +test Index: /branches/eam_branch_20080511/ppSim/Makefile.am =================================================================== --- /branches/eam_branch_20080511/ppSim/Makefile.am (revision 17629) +++ /branches/eam_branch_20080511/ppSim/Makefile.am (revision 17629) @@ -0,0 +1,3 @@ +SUBDIRS = src + +CLEANFILES = *~ core core.* Index: /branches/eam_branch_20080511/ppSim/autogen.sh =================================================================== --- /branches/eam_branch_20080511/ppSim/autogen.sh (revision 17629) +++ /branches/eam_branch_20080511/ppSim/autogen.sh (revision 17629) @@ -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=ppConfigDump +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal $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/eam_branch_20080511/ppSim/configure.ac =================================================================== --- /branches/eam_branch_20080511/ppSim/configure.ac (revision 17629) +++ /branches/eam_branch_20080511/ppSim/configure.ac (revision 17629) @@ -0,0 +1,30 @@ +AC_PREREQ(2.61) + +AC_INIT([ppConfigDump], [1.1.0], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) + +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([PSMODULE], [psmodules >= 1.0.0]) +PKG_CHECK_MODULES([PSASTRO], [psastro >= 0.9.0]) + +IPP_STDOPTS +CFLAGS="${CFLAGS=} -Wall -Werror" + +AC_CONFIG_FILES([ + Makefile + src/Makefile +]) +AC_OUTPUT Index: /branches/eam_branch_20080511/ppSim/notes.txt =================================================================== --- /branches/eam_branch_20080511/ppSim/notes.txt (revision 17629) +++ /branches/eam_branch_20080511/ppSim/notes.txt (revision 17629) @@ -0,0 +1,47 @@ + +ppSim development work still needed: + +* save sources in output files +* load input image + +* optionally apply detrend data (bias, dark, etc). when adding data + to a loaded image, we only need to inject fake stars. + +* examine the range of star fluxes. (too many are being generated). +* check the overlap between the real and random star fluxes + +* generate stars with / without known exposure correction (ie, stars + could be generated with known fluxes or know magnitudes) + +* generate galaxies (objects with random values for major, minor, theta) +* generate convolved galaxies + +* should not have to supply ZP, Seeing, PA, Scale on command line for all ppSim runs +* DVO.CATDIR path://SIMTEST is not being interpolated? +* -D PSASTRO:DVO.CATDIR is not being interpretted? + +* fake stars do not have random magnitudes? +* add generated stars to an output file (pmSource list?) +* make star source an option (random / dvo database) +* add the gain to the metadata + +ppSimLoop: + + * load catalogs stars (ra, dec, mag, x, y) + * generate random stars (ra, dec, mag, x, y) + * loop over chip,cell,readout + ** generate output images (signal, variance) + ** add the bias + ** add the dark + ** add the sky signal + ** add the stars (modified by the flat, shutter) + ** add the galaxies (modified by the flat, shutter) + ** add the poisson noise + ** add the overscan + readnoise + ** save the stars + ** save the galaxies + + ** update the concept info + ** update the header wcs info + ** write out the image data + Index: /branches/eam_branch_20080511/ppSim/psf.txt =================================================================== --- /branches/eam_branch_20080511/ppSim/psf.txt (revision 17629) +++ /branches/eam_branch_20080511/ppSim/psf.txt (revision 17629) @@ -0,0 +1,59 @@ + +Code to load a psf model and construct psf images. This assumes there +is an analysis loop which is performing pmFPAfileIOChecks at the +appropriate levels. + +** to supply a psf model from the command, define the command-line arguments (ppSimArguments.c) + + pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); + +** if a psf model is supplied, generate a pmFPAfile to carry it. bind + it to the pmFPAfile in your main loop (ppSimCreate.c): + + // have we supplied a psf model? + if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { + bool status = false; + + // tie the psf file to the chipMosaic + // file is the pmFPAfile used for the analysis loop + pmFPAfileBindFromArgs(&status, file, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); + // XXX free things to be freed... + return NULL; + } + } + +** generate source models from the psf and insert into a readout: + + // load the psf model from the corresponding chip. note: you may + // need to select the chip from the current readout, eg: (ppSimInsertSources.c) + // pmCell *cell = readout->parent; + // pmChip *chip = cell->parent; + pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF"); + + + // instantiate a model for the PSF at this location, set desired flux + + // note that pmModelFromPSFforXY takes the chip coordinate. you many + // need to convert the x,y coordinate in the readout to/from the x,y + // coordinate. (ppSimInsertSources.c) + pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0); + pmModelSetFlux (model, flux); + + // define the radius of valid pixels + float radius = model->modelRadius (model->params, noise); + radius = PS_MAX (radius, 1.0); + + // construct a source, with model flux pixels set, based on the model + pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR); + + // insert the source flux in the image. dX,dY is the coordinate of + // the readout 0,0 pixel in the chip frame. + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + +** given a pmSource, generate the its cached model flux: + + // this fills in the pixels in the image source->modelFlux + pmSourceCacheModel (source, maskVal); + Index: /branches/eam_branch_20080511/ppSim/src/.cvsignore =================================================================== --- /branches/eam_branch_20080511/ppSim/src/.cvsignore (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/.cvsignore (revision 17629) @@ -0,0 +1,9 @@ +config.h +config.h.in +Makefile +Makefile.in +ppSim +stamp-h1 +.deps +.libs +ppSimSequence Index: /branches/eam_branch_20080511/ppSim/src/Makefile.am =================================================================== --- /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17629) @@ -0,0 +1,48 @@ +bin_PROGRAMS = ppSim ppSimSequence + +ppSim_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS) $(ppSim_CFLAGS) +ppSim_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS) +ppSim_SOURCES = \ + ppSim.c \ + ppSimArguments.c \ + ppSimCreate.c \ + ppSimLoadStars.c \ + ppSimMakeStars.c \ + ppSimMakeGalaxies.c \ + ppSimMakeBiassec.c \ + ppSimMakeBias.c \ + ppSimMakeDark.c \ + ppSimMakeSky.c \ + ppSimInsertStars.c \ + ppSimInsertGalaxies.c \ + ppSimAddOverscan.c \ + ppSimAddNoise.c \ + ppSimSaturate.c \ + ppSimBounds.c \ + ppSimStars.c \ + ppSimSetPSF.c \ + ppSimUtils.c \ + ppSimLoop.c \ + ppSimBadPixels.c + +ppSimSequence_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS) $(ppSim_CFLAGS) +ppSimSequence_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS) +ppSimSequence_SOURCES = \ + ppSimSequence.c \ + ppSimSequenceBias.c \ + ppSimSequenceDark.c \ + ppSimSequenceFlat.c \ + ppSimSequenceObject.c + +noinst_HEADERS = \ + ppSim.h \ + ppSimSequence.h + + +CLEANFILES = *~ + +clean-local: + -rm -f TAGS +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /branches/eam_branch_20080511/ppSim/src/ppSim.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSim.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSim.c (revision 17629) @@ -0,0 +1,35 @@ +# include "ppSim.h" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + if (!pmModelClassInit ()) abort(); + + pmConfig *config = pmConfigRead(&argc, argv, PPSIM_RECIPE); // Configuration + if (!config) { + psErrorStackPrint(stderr, "Unable to read configurations."); + exit(PS_EXIT_CONFIG_ERROR); + } + + ppSimArguments(argc, argv, config); + + if (!ppSimCreate(config)) { + psErrorStackPrint(stderr, "Unable to create output file."); + psFree(config); + exit(PS_EXIT_CONFIG_ERROR); + } + + if (!ppSimLoop(config)) { + psErrorStackPrint(stderr, "Unable to generate data."); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + + psFree(config); + pmModelClassCleanup(); + pmConfigDone(); + pmConceptsDone(); + psLibFinalize(); + + return PS_EXIT_SUCCESS; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSim.h =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17629) @@ -0,0 +1,141 @@ +#ifndef PP_SIM_H +#define PP_SIM_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include // for unlink +#include +#include +#include + +#define PPSIM_RECIPE "PPSIM" +#define OUTPUT_FILE "PPSIM.OUTPUT" + +// Compare a value with minimum and maximum values, replacing where required. +#define COMPARE(VALUE,MIN,MAX) { \ + if (VALUE < MIN) { MIN = VALUE; } \ + if (VALUE > MAX) { MAX = VALUE; } \ + } + +// Return cell position, given an FPA position; calculations are all done in pixel units +#define PPSIM_FPA_TO_CELL(pos, cell0, cellParity, binning, chip0, chipParity) \ + (((pos) - (chip0))*(chipParity) - (cell0))*(cellParity) / (binning) + + +// Return FPA position, given a cell position; calculations are all done in pixel units +#define PPSIM_CELL_TO_FPA(pos, cell0, cellParity, binning, chip0, chipParity) \ + ((chip0) + (binning)*(chipParity)*((cell0) + (cellParity)*(pos))) + +// Type of image to simulate +typedef enum { + PPSIM_TYPE_NONE, // No type set + PPSIM_TYPE_BIAS, // Bias image + PPSIM_TYPE_DARK, // Dark image + PPSIM_TYPE_FLAT, // Flat-field image + PPSIM_TYPE_OBJECT // Object image +} ppSimType; + +typedef struct { + double ra; + double dec; + float mag; + float x; + float y; + float flux; + float peak; +} ppSimStar; + +typedef struct { + double ra; + double dec; + float mag; + float x; + float y; + float flux; + + float peak; + float Rmaj; + float Rmin; + float theta; + float index; +} ppSimGalaxy; + +ppSimStar *ppSimStarAlloc (); +ppSimGalaxy *ppSimGalaxyAlloc (); + +/// Parse command-line arguments +void ppSimArguments(int argc, char *argv[], ///< Command-line arguments + pmConfig *config ///< Configuration + ); + +/// Create output file +/// +/// Returns a borrowed pointer to the FPA file. +pmFPAfile *ppSimCreate(pmConfig *config ///< Configuration + ); + +// Return bounds of a chip, based on the concepts +psRegion *ppSimChipBounds(const pmChip *chip, // Chip for which to determine size + pmFPAview *view // View for chip + ); + +// Return bounds of an FPA, based on the concepts +psRegion *ppSimFPABounds(const pmFPA *fpa // FPA for which to determine size + ); + +/// Loop over the output file, generating simulated data +bool ppSimLoop(pmConfig *config ///< Configuration + ); + +psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config); +psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) ; +bool ppSimMakeDark (pmReadout *readout, pmConfig *config); +bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config); + +bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config); +bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng); +bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config); + +bool ppSimInitHeader(pmConfig *config, + pmFPA *fpa, + pmChip *chip, + pmCell *cell); + +bool ppSimSaturate(pmReadout *readout, // Image to apply saturation + const pmConfig *config); // Saturation level + +bool ppSimUpdateConceptsFPA (pmFPA *fpa, pmConfig *config); +bool ppSimUpdateConceptsCell (pmCell *cell, pmConfig *config); + +bool ppSimAddOverscan (pmReadout *readout, pmConfig *config, psVector *biasCols, psVector *biasRows, psRandom *rng); + +bool ppSimAddNoise(psImage *signal, + psImage *variance, + const pmCell *cell, + const pmConfig *config, + const psRandom *rng // Random number generator + ); + +bool ppSimSetPSF (pmChip *chip, pmConfig *config); + +bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng); +bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config); + +/// Add bad pixels to an image +bool ppSimBadPixels(pmReadout *readout, ///< Readout for which to generate bad pixels + const pmConfig *config, ///< Configuration + psRandom *rng ///< Random number generator + ); + +float ppSimStarSkyNoise (float skySigma, float seeingSigma); +float ppSimStarPeakToFlux (float peak, float seeingSigma); +float ppSimStarFluxToPeak (float flux, float seeingSigma); +float ppSimFluxToMag (float flux, float zp); +float ppSimMagToFlux (float mag, float zp); + +#endif Index: /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 17629) @@ -0,0 +1,42 @@ +# include "ppSim.h" + +// Add noise to an image +bool ppSimAddNoise(psImage *signal, // Signal image, modified and returned + psImage *variance, + const pmCell *cell, + const pmConfig *config, + const psRandom *rng // Random number generator + ) +{ + assert(signal->type.type == PS_TYPE_F32); + assert(signal->numCols == variance->numCols && signal->numRows == variance->numRows); + + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // CCD gain, e/ADU + if (isnan(gain)) { + psWarning("CELL.GAIN is not set; reverting to recipe value GAIN."); + gain = psMetadataLookupF32(&mdok, recipe, "GAIN"); + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find GAIN in recipe."); + return false; + } + } + + // Add the noise into the image + for (int y = 0; y < signal->numRows; y++) { + for (int x = 0; x < signal->numCols; x++) { + // XXX is psRandomGaussian doing this reasonally optimally? + // (generate a static array with the cumulative distribution, use the + // random number to select a bin from the histogram) + signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * psRandomGaussian(rng); + signal->data.F32[y][x] /= gain; // Converting to ADU + } + } + + // XXX why return this?? + return true; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 17629) @@ -0,0 +1,38 @@ +#include "ppSim.h" + +bool ppSimAddOverscan (pmReadout *readout, pmConfig *config, psVector *biasCols, psVector *biasRows, psRandom *rng) { + + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + float readnoise = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.READNOISE");// CCD read noise, e + if (isnan(readnoise)) { + psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE."); + readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + return NULL; + } + } + + // Add overscan + // XXX put this in a wrapper + for (int j = 0; j < biasCols->n; j++) { + psImage *signal = psImageAlloc(biasCols->data.S32[j], biasRows->n, PS_TYPE_F32); // Overscan + for (int y = 0; y < signal->numRows; y++) { + for (int x = 0; x < signal->numCols; x++) { + signal->data.F32[y][x] = biasRows->data.F32[y]; + } + } + psImage *variance = psImageAlloc(biasCols->data.S32[j], biasRows->n, PS_TYPE_F32); // Variance + psImageInit(variance, PS_SQR(readnoise)); + + ppSimAddNoise(signal, variance, readout->parent, config, rng); + psListAdd(readout->bias, PS_LIST_TAIL, signal); + + psFree(variance); + psFree(signal); // Drop reference + } + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 17629) @@ -0,0 +1,273 @@ +#include "ppSim.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config // Configuration + ) +{ + fprintf(stderr, "\nPan-STARRS data simulator\n\n"); + fprintf(stderr, "Usage: %s -camera CAMERA_NAME (output)\n", program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psFree(arguments); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +// Get a value from the command-line arguments and add it to recipe options +bool valueArgRecipe(psMetadata *options, // Target to which to add value + const char *recipeName, // Name for value in the recipe + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments + ) +{ + bool status; // Status of MD lookup + float value = psMetadataLookupF32(&status, arguments, argName); // Value of interest + if (isnan(value)) return true; + status = psMetadataAddF32(options, PS_LIST_TAIL, recipeName, 0, NULL, value); + return status; +} + +// this function supplements the RECIPE:OPTIONS folder with command-line options +void ppSimArguments(int argc, char *argv[], pmConfig *config) +{ + bool mdok; // Status of MD lookup + + assert(config); + + // save the following command-line options in the arguments structure + psMetadata *arguments = psMetadataAlloc(); // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-format", 0, "Camera format name", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Exposure type (BIAS|DARK|FLAT|OBJECT)", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-filter", 0, "Filter name", NULL); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-exptime", 0, "Exposure time (s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-biaslevel", 0, "Bias level (e)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-biasrange", 0, "Bias range (e)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-darkrate", 0, "Dark rate (e/s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-flatsigma", 0, "Flat sigma", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-flatrate", 0, "Flat rate (e/s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-shuttertime", 0, "Shutter time (s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-skyrate", 0, "Sky rate (e/s)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-skymags", 0, "Sky brightness in mags / square arcsec", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-ra", 0, "RA (degrees)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-dec", 0, "Dec (degrees)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-pa", 0, "Position angle (degrees)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-scale", 0, "Plate scale (arcsec/pixel)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Photometric zero point", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-seeing", 0, "Seeing FWHM (arcsec)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-starslum", 0, "Fake star luminosity function slope", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsmag", 0, "Brightest magnitude for fake stars", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsdensity", 0, "Density of fake stars at magnitude", NAN); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-psfclass", 0, "Type of PSF model", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-galmodel", 0, "Type of Galaxy model", NULL); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin", 0, "Binning in x and y", 1); + + if (psArgumentGet (argc, argv, "-h")) { usage(argv[0], arguments, config); } + if (psArgumentGet (argc, argv, "--h")) { usage(argv[0], arguments, config); } + if (psArgumentGet (argc, argv, "--help")) { usage(argv[0], arguments, config); } + + pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); + + // only one of -camera and -file is needed + bool status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); + if (!config->camera && !status) { + psErrorStackPrint(stderr, "A camera name (-camera NAME) or an image (-file NAME) must be specified"); + exit(PS_EXIT_CONFIG_ERROR); + } + + if (!psArgumentParse(arguments, &argc, argv)) { + exit(PS_EXIT_CONFIG_ERROR); + } + + if (argc != 2) { + psErrorStackPrint(stderr, "Missing output filename"); + exit(PS_EXIT_CONFIG_ERROR); + } + + // apply an alternate camera format + psString formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format + if (formatName) { + // XXX delay the config below until ppSimCreate? + config->formatName = psMemIncrRefCounter(formatName); + + psMetadata *formats = psMetadataLookupMetadata(NULL, config->camera, "FORMATS"); // The camera formats + if (!formats) { + psErrorStackPrint(stderr, "Unable to find FORMATS in camera configuration."); + psFree(arguments); + psFree(config); + exit(PS_EXIT_CONFIG_ERROR); + } + psMetadata *format = psMetadataLookupMetadata(NULL, formats, formatName); // Format of interest + if (!format) { + psErrorStackPrint(stderr, "Unable to find format %s in camera FORMATS.", formatName); + psFree(arguments); + psFree(config); + exit(PS_EXIT_CONFIG_ERROR); + } + config->format = psMemIncrRefCounter(format); + } + + // specify the type of simulated image to produce + // XXX this should not be required if we supplied INPUT + const char *typeStr = psMetadataLookupStr(NULL, arguments, "-type"); // Exposure type + if (!typeStr) { + psErrorStackPrint(stderr, "An exposure type must be specified using -type"); + exit(PS_EXIT_CONFIG_ERROR); + } + ppSimType type = PPSIM_TYPE_NONE; // Type to simulate + if (strcasecmp(typeStr, "BIAS") == 0) { + type = PPSIM_TYPE_BIAS; + } else if (strcasecmp(typeStr, "DARK") == 0) { + type = PPSIM_TYPE_DARK; + } else if (strcasecmp(typeStr, "FLAT") == 0) { + type = PPSIM_TYPE_FLAT; + } else if (strcasecmp(typeStr, "OBJECT") == 0) { + type = PPSIM_TYPE_OBJECT; + } else { + psErrorStackPrint(stderr, "Unrecognised exposure type: %s", typeStr); + exit(PS_EXIT_CONFIG_ERROR); + } + assert(type != PPSIM_TYPE_NONE); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "TYPE", 0, "Exposure type", type); + + const char *filter = psMetadataLookupStr(NULL, arguments, "-filter"); // Filter name + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "FILTER", 0, "Filter name", filter); + + // save the following additional recipe values based on command-line options + // these options override the PPSIM recipe values loaded from recipe files + psMetadata *options = pmConfigRecipeOptions (config, PPSIM_RECIPE); + if (!options) { + psErrorStackPrint(stderr, "Unable to find recipe options for %s", PPSIM_RECIPE); + psFree(arguments); + psFree(config); + exit(PS_EXIT_CONFIG_ERROR); + } + + float expTime; + if (type == PPSIM_TYPE_BIAS) { + expTime = 0.0; + } else { + expTime = psMetadataLookupF32(NULL, arguments, "-exptime"); // Exposure time + if (isnan(expTime)) { + psErrorStackPrint(stderr, "The exposure time must be specified using -exptime"); + exit(PS_EXIT_CONFIG_ERROR); + } + } + psMetadataAddF32(options, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", expTime); + + // these values all get moved from arguments to RECIPE:OPTIONS + valueArgRecipe(options, "BIAS.LEVEL", arguments, "-biaslevel"); + valueArgRecipe(options, "BIAS.RANGE", arguments, "-biasrange"); + valueArgRecipe(options, "DARK.RATE", arguments, "-darkrate"); + valueArgRecipe(options, "FLAT.SIGMA", arguments, "-flatsigma"); + valueArgRecipe(options, "FLAT.RATE", arguments, "-flatrate"); + valueArgRecipe(options, "SHUTTER.TIME", arguments, "-shuttertime"); + valueArgRecipe(options, "SKY.RATE", arguments, "-skyrate"); + valueArgRecipe(options, "STARS.LUM", arguments, "-starslum"); + valueArgRecipe(options, "STARS.MAG", arguments, "-starsmag"); + valueArgRecipe(options, "STARS.DENSITY", arguments, "-starsdensity"); + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + int biasOrder = psMetadataLookupS32(&mdok, recipe, "BIAS.ORDER"); // Overscan polynomial order + if (!mdok) { + psWarning("BIAS.ORDER(S32) is not set in the recipe %s --- assuming %d", PPSIM_RECIPE, biasOrder); + } + psMetadataAddS32(options, PS_LIST_TAIL, "BIAS.ORDER", 0, "Overscan polynomial order", biasOrder); + + int binning = psMetadataLookupS32(NULL, arguments, "-bin"); // Binning in x and y + if (binning <= 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Binning (%d) is non-positive.", binning); + exit(PS_EXIT_CONFIG_ERROR); + } + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BINNING", 0, "Binning in x and y", binning); + + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + if (type == PPSIM_TYPE_OBJECT) { + // Load values required for adding stars + float ra0 = psMetadataLookupF32(NULL, arguments, "-ra"); // Right Ascension of boresight + float dec0 = psMetadataLookupF32(NULL, arguments, "-dec"); // Declination of boresight + float pa = psMetadataLookupF32(NULL, arguments, "-pa"); // Position angle + float seeing = psMetadataLookupF32(NULL, arguments, "-seeing"); // seeing (FWHM in arcsec) + + // XXX scale and zp should be supplied by the config file (allow override, but this is camera-dependent) + if (isnan(ra0) || isnan(dec0) || isnan(pa) || isnan(seeing)) { + psErrorStackPrint(stderr, "-ra, -dec, -pa, -scale, -zp, -seeing must be specified for OBJECT type"); + exit(PS_EXIT_CONFIG_ERROR); + } + + float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point + if (isnan(zp)) { + // use the filter to get the zeropoint from the recipe + psMetadataItem *zpItem = psMetadataLookup (recipe, "ZEROPTS"); + // check that item is multi... + + psArray *entries = psListToArray (zpItem->data.list); + + // search for matching filter + for (int i = 0; i < entries->n; i++) { + psMetadataItem *item = entries->data[i]; + psMetadata *entry = item->data.V; + + char *filterName = psMetadataLookupStr (&status, entry, "FILTER"); + assert (filterName); + + if (strcmp(filterName, filter)) continue; + + zp = psMetadataLookupF32 (&status, entry, "ZERO_PT"); + assert (status); + break; + } + psFree (entries); + } + + float scale = psMetadataLookupF32(NULL, arguments, "-scale"); // Plate scale + if (isnan(scale)) { + scale = psMetadataLookupF32 (&status, recipe, "PIXEL.SCALE"); + } + + float skymags = psMetadataLookupF32(NULL, arguments, "-skymags"); // Position angle + if (isnan(skymags)) { + skymags = psMetadataLookupF32 (&status, recipe, "SKY.MAGS"); + } + + psMetadataAddF32(options, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0); + psMetadataAddF32(options, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)", dec0 * M_PI / 180.0); + psMetadataAddF32(options, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",pa * M_PI / 180.0); + + // the user supplies FWHM in arcsec; here we convert to Sigma in pixels + psMetadataAddF32(options, PS_LIST_TAIL, "SEEING", 0, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); + + psMetadataAddF32(options, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale); + psMetadataAddF32(options, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp); + psMetadataAddF32(options, PS_LIST_TAIL, "SKY.MAGS", 0, "sky surface brightness", skymags); + + const char *psfClass = psMetadataLookupStr(NULL, arguments, "-psfclass"); // PSF model class + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "PSF.MODEL", 0, "PSF model class", psfClass); + + const char *galModel = psMetadataLookupStr(NULL, arguments, "-galmodel"); // Galaxy model name + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "GALAXY.MODEL", 0, "Galaxy model", galModel); + } + + psFree(arguments); + return; +} + +/* the following elements come from the config->arguments: + + PSPHOT.PSF + INPUT + TYPE + FILTER + BIAS.ORDER + BINNING + OUTPUT + PSF.MODEL + GALAXY.MODEL + + all othr values should come from the recipe +*/ Index: /branches/eam_branch_20080511/ppSim/src/ppSimBadPixels.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimBadPixels.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimBadPixels.c (revision 17629) @@ -0,0 +1,64 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include + +#include "ppSim.h" + + +// The idea is that the pattern should be completely *deterministic* (despite the use of 'random' numbers --- +// with the seed specified, they should be deterministic) so that multiple calls of this function will result +// in the same pattern. These are set to (really) random values so that they do not dark-subtract or +// flat-field. + +bool ppSimBadPixels(pmReadout *readout, const pmConfig *config, psRandom *rng) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_IMAGE_NON_NULL(readout->image, false); + PS_ASSERT_PTR_NON_NULL(readout->parent, false); + PS_ASSERT_PTR_NON_NULL(config, false); + + bool mdok; // Status of MD lookup + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + if (!recipe) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSIM_RECIPE); + return false; + } + + psU64 seed = psMetadataLookupU64(&mdok, recipe, "BADPIX.SEED"); // Seed for RNG + if (seed == 0 || !mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "BADPIX.SEED not set, or zero."); + return false; + } + + float frac = psMetadataLookupF32(&mdok, recipe, "BADPIX.FRAC"); // Fraction of bad pixels + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "BADPIX.FRAC not set."); + return false; + } + if (frac == 0.0) { + // Nothing to do + return true; + } + + psRandom *pseudoRNG = psRandomAlloc(PS_RANDOM_TAUS, seed); // Pseudo-random number generator + + psImage *image = readout->image; // Image of interest + int numCols = image->numCols, numRows = image->numRows; // Size of image + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (psRandomUniform(pseudoRNG) < frac) { + image->data.F32[y][x] *= psRandomGaussian(rng); + } + } + } + + psFree(pseudoRNG); + + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimBounds.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimBounds.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimBounds.c (revision 17629) @@ -0,0 +1,98 @@ +#include "ppSim.h" + +// Return bounds of a chip, based on the concepts +psRegion *ppSimChipBounds(const pmChip *chip, // Chip for which to determine size + pmFPAview *view // View for chip + ) +{ + assert(chip); + assert(view); + + // Bounds of chip + int xMin = +INT_MAX; + int xMax = -INT_MAX; + int yMin = +INT_MAX; + int yMax = -INT_MAX; + + int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); + int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); + int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY"); + int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY"); + + if ((xParityChip != 1 && xParityChip != -1) || (yParityChip != 1 && yParityChip != -1)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Chip parities are not set."); + psFree(view); + return NULL; + } + + pmCell *cell; // Cell from chip + while ((cell = pmFPAviewNextCell(view, chip->parent, 1))) { + int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE"); + int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE"); + + if (xSize == 0 || ySize == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell sizes are not set."); + psFree(view); + return NULL; + } + + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + if ((xParityCell != 1 && xParityCell != -1) || (yParityCell != 1 && yParityCell != -1)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell parities are not set."); + psFree(view); + return NULL; + } + + // "Left", "Right", "Bottom" and "Top" don't take into account the parity + int cellLeft = PPSIM_CELL_TO_FPA(0, x0Cell, xParityCell, 1, x0Chip, xParityChip); + int cellRight = PPSIM_CELL_TO_FPA(xSize, x0Cell, xParityCell, 1, x0Chip, xParityChip); + int cellBottom = PPSIM_CELL_TO_FPA(0, y0Cell, yParityCell, 1, y0Chip, yParityChip); + int cellTop = PPSIM_CELL_TO_FPA(ySize, y0Cell, yParityCell, 1, y0Chip, yParityChip); + + COMPARE(cellLeft, xMin, xMax); + COMPARE(cellRight, xMin, xMax); + COMPARE(cellBottom, yMin, yMax); + COMPARE(cellTop, yMin, yMax); + } + + return psRegionAlloc(xMin, xMax, yMin, yMax); +} + +// Return bounds of an FPA, based on the concepts +psRegion *ppSimFPABounds(const pmFPA *fpa // FPA for which to determine size + ) +{ + assert(fpa); + + pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA + + // Bounds of focal plane + int xMin = +INT_MAX; + int xMax = -INT_MAX; + int yMin = +INT_MAX; + int yMax = -INT_MAX; + + pmChip *chip; // Chip from FPA + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + psRegion *bounds = ppSimChipBounds(chip, view); // Bounds for chip + COMPARE(bounds->x0, xMin, xMax); + COMPARE(bounds->x1, xMin, xMax); + COMPARE(bounds->y0, yMin, yMax); + COMPARE(bounds->y1, yMin, yMax); + psFree(bounds); + } + + psFree(view); + psRegion *bounds = psRegionAlloc(xMin, xMax, yMin, yMax); + + if (!bounds || (bounds->x0 == 0.0 && bounds->x1 == 0.0) || (bounds->y0 == 0.0 && bounds->y1 == 0.0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine bounds of FPA"); + return NULL; + } + + return bounds; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17629) @@ -0,0 +1,132 @@ +# include "ppSim.h" + +// XXX this function forces us to define the camera (on the command line) and format (via the +// PPSIM.OUTPUT entry). In this case, we need to set config->format,formatName based on these +// values. This will be a problem when we want to load an input image (in order to add fake +// stars). We will need to add some logic in ppSimArguments to distinguish the cases of 1) +// input image and 2) specified camera + +pmFPAfile *ppSimCreate(pmConfig *config) +{ + bool status; + bool simImage = false; + PS_ASSERT_PTR_NON_NULL(config, NULL); + pmFPA *fpa = NULL; + + // the input image defines the camera. if it is not supplied, the user must have + // supplied a camera and other metadata on the command line + pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPIMAGE.INPUT", "INPUT"); + if (!input) { + // if we have not specified the camera already, we need to interpolate the recipes associated with this camera, and read other command-line recipes + if (!pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) { + psError(PS_ERR_IO, false, "Error merging recipes from camera config for %s", config->cameraName); + return NULL; + } + + simImage = true; + fpa = pmFPAConstruct(config->camera); // FPA to contain the observation + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + return NULL; + } + + } else { + simImage = false; + if (input->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPIMAGE.INPUT is not of type IMAGE"); + return NULL; + } + fpa = input->fpa; + } + + // define the output image file + pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE); + if (!file) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s. " + "Did you forget to specify the format?", OUTPUT_FILE); + return NULL; + } + if (file->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "%s type is not IMAGE", OUTPUT_FILE); + psFree(fpa); + psFree(file); + return NULL; + } + file->save = true; + + config->format = psMemIncrRefCounter (file->format); + config->formatName = psStringCopy (file->formatName); + + // have we supplied a psf model? + if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { + bool status = false; + + // tie the psf file to the chipMosaic + pmFPAfileBindFromArgs(&status, file, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); + psFree(fpa); + psFree(file); + return NULL; + } + } + + // XXX only invoke this code for OBJECT types of images + // PPSIM.SOURCES carries the constructed, fake sources with their true parameters + pmFPAfile *simSources = pmFPAfileDefineOutput (config, file->fpa, "PPSIM.SOURCES"); + if (!simSources) { + psError(PS_ERR_UNKNOWN, false, "Cannot find a rule for PPSIM.SOURCES"); + return false; + } + simSources->save = true; + + // if we have loaded an input image, we do not need to populate the fpa + if (!simImage) { + return file; + } + + pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level at which PHU goes + + pmFPAview *view = pmFPAviewAlloc(0);// View for current level + + if (phuLevel == PM_FPA_LEVEL_FPA) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + psFree(file); + return NULL; + } + } + + pmChip *chip; // Chip from FPA + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + if (phuLevel == PM_FPA_LEVEL_CHIP) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + psFree(file); + return NULL; + } + } + + pmCell *cell; // Cell from chip + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + if (phuLevel == PM_FPA_LEVEL_CELL) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + psFree(file); + return NULL; + } + } + } + } + + psFree(fpa); + psFree(view); + + return file; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimInsertGalaxies.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimInsertGalaxies.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimInsertGalaxies.c (revision 17629) @@ -0,0 +1,192 @@ +# include "ppSim.h" +static char *defaultModel = "PS_MODEL_SERSIC"; + +bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config) { + + bool mdok; + + assert (readout); + assert (galaxies); + + if (!galaxies->n) { return true; } + + // XXX is this needed? + // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources + + pmCell *cell = readout->parent; + pmChip *chip = cell->parent; + + // XXX this is an estimate of the sky noise based on the inputs to the image simulation. + // XXX update this to allow the estimate based on the measured sky background + // XXX this is missing the gain. + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate + + float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e + if (isnan(readnoise)) { + psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE."); + readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + return false; + } + } + + float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate + if (isnan(skyRate)) { + float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok); + float scale = psMetadataLookupF32(&mdok, recipe, "SCALE"); assert (mdok); + float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + + // Rough noise estimate, appropriate for entire cell (use for source radius?) + float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + + int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); + int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); + int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY"); + int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY"); + + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + + // determine the galaxy model + char *modelName = psMetadataLookupStr(&mdok, config->arguments, "GALAXY.MODEL"); // galaxy model name + if (modelName == NULL) { + modelName = defaultModel; + } + pmModelType type = pmModelClassGetType (modelName); + if (type == -1) { + psError (PS_ERR_UNKNOWN, false, "invalid model name"); + return false; + } + + int nParam = pmModelClassParameterCount (type); + + pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF"); + assert (psf); + + int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning); + int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning); + + // psMetadataLookupPtr (readout->analysis, "PSPHOT.SOURCES", 0); + + psArray *sources = psArrayAllocEmpty (galaxies->n); + + + // add sources to the readout image & weight + for (long i = 0; i < galaxies->n; i++) { + ppSimGalaxy *galaxy = galaxies->data[i]; + + // galaxy->x,y are in fpa coordinates + + // Position on the cell and peak flux + float xChip = PM_FPA_TO_CHIP(galaxy->x, x0Chip, xParityChip); + float yChip = PM_FPA_TO_CHIP(galaxy->y, y0Chip, yParityChip); + + // Position on the cell and peak flux + float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning); + float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning); + + if (xCell < 0) continue; + if (yCell < 0) continue; + if (xCell > readout->image->numCols) continue; + if (yCell > readout->image->numRows) continue; + // XXX need to apply col0, row0 if readout is a subarray + + // XXX apply the expCorr to the galaxy->flux before setting the model flux + + // instantiate a model for the PSF at this location, set desired flux + pmModel *model = pmModelAlloc (type); + + psF32 *PAR = model->params->data.F32; + + PAR[PM_PAR_I0] = galaxy->peak; + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_XPOS] = xChip; + PAR[PM_PAR_YPOS] = yChip; + + psEllipseAxes axes; + axes.major = galaxy->Rmaj; + axes.minor = galaxy->Rmin; + axes.theta = galaxy->theta; + pmPSF_AxesToModel (PAR, axes); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + if (nParam == 8) { + PAR[PM_PAR_7] = galaxy->index; + } + + // set the normalization to get the desired flux + // XXX for now, use peak: pmModelSetFlux (model, galaxy->flux); + + // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) + float radius = model->modelRadius (model->params, roughNoise); + radius = PS_MAX (radius, 1.0); + // XXX the exp(-r^0.25) models can go way out if allowed... + radius = PS_MIN (radius, 150.0); + + // construct a source, with model flux pixels set, based on the model + pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_EXTENDED); + + // XXX set the mag & err values (should this be done in pmSourceFromModel?) + // XXX i should be applying the gain and the correct effective area + source->psfMag = -2.5*log10(galaxy->flux); + source->extMag = -2.5*log10(galaxy->flux); + source->errMag = sqrt(Area*PS_SQR(roughNoise) + galaxy->flux) / galaxy->flux; + + // XXX add the sources to a source array + +# if (0) + if (CONVOLVED_FIT) { + + // select the PSF (XXX : move out of loop) + pmPSF *psfModel = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF"); + assert (psf); + + // supply the psf to the source (normalized?) + source->modelPSF = pmModelFromPSFforXY (psfModel, xChip, yChip, 1.0); + + // instantiate the psf flux + pmSourceCachePSF (source, maskVal); + + // convert the cached cached psf model for this source to a psKernel + // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius) + // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius) + psKernel *psf = psphotKernelFromPSF (source, psfSize); + + // instantiate the source model flux + // XXX this may need to handle the offset? + pmSourceCacheModel (source, maskVal); + + // make a storage buffer for the output image + psImage *buffer = psImageCopy (source->modelFlux, PS_TYPE_F32); + + // convolve the psf image with the model image + psImageConvolveDirect (buffer, source->modelFlux, psf); + + // save the result back in the source + psFree (source->modelFlux); + source->modelFlux = buffer; + } +# endif + + // insert the source flux in the image + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + psArrayAdd (sources, 100,source); + } + + // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT" + // psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, "psphot sources", sources); + + // XXX many leaks in here, i think + return true; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 17629) @@ -0,0 +1,143 @@ +# include "ppSim.h" + +// Reset a pointer: free and set to NULL +#define RESET(PTR) \ + psFree(PTR); \ + PTR = NULL; + + +bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) { + + bool mdok; + + assert (readout); + assert (stars); + + if (!stars->n) { return true; } + + // XXX is this needed? + // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources + + pmCell *cell = readout->parent; + pmChip *chip = cell->parent; + + // XXX this is an estimate of the sky noise based on the inputs to the image simulation. + // XXX update this to allow the estimate based on the measured sky background + // XXX this is missing the gain. + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate + + float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e + if (isnan(readnoise)) { + psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE."); + readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + return false; + } + } + + float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate + if (isnan(skyRate)) { + float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok); + float scale = psMetadataLookupF32(&mdok, recipe, "SCALE"); assert (mdok); + float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + + // Rough noise estimate, appropriate for entire cell (use for source radius?) + float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + + int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); + int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); + int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY"); + int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY"); + + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + + pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF"); + assert (psf); + + int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning); + int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning); + + psArray *sources = psArrayAllocEmpty (stars->n); + + // add sources to the readout image & weight + psTrace("ppSim", 1, "Inserting %ld stars...\n", stars->n); + for (long i = 0; i < stars->n; i++) { + ppSimStar *star = stars->data[i]; + psTrace("ppSim", 10, "Inserting star at %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux); + + // star->x,y are in fpa coordinates + + // Position on the cell and peak flux + float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip); + float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip); + + // Position on the cell and peak flux + float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning); + float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning); + + // XXX Note, the below does not put the edges of stars on the readout if they fall slightly off + // This will be visible as cut-off stars at amplifier boundaries (e.g., Megacam) + if (xCell < 0) continue; + if (yCell < 0) continue; + if (xCell > readout->image->numCols) continue; + if (yCell > readout->image->numRows) continue; + // XXX need to apply col0, row0 if readout is a subarray + + // Apply the expCorr to the star->flux before setting the model flux + float flux = star->flux * expCorr->data.F32[(int)yCell][(int)xCell]; + + // instantiate a model for the PSF at this location, set desired flux + pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0); + pmModelSetFlux (model, flux); + + // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) + float radius = model->modelRadius (model->params, roughNoise); + radius = PS_MAX (radius, 1.0); + + // construct a source, with model flux pixels set, based on the model + pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR); + + // XXX set the mag & err values (should this be done in pmSourceFromModel?) + // XXX i should be applying the gain and the correct effective area + psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0); + psF64 Area = 2.0 * M_PI * axes.major * axes.minor; + + source->psfMag = -2.5*log10(flux); + source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux; + + // XXX add the sources to a source array + + // insert the source flux in the image + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + psArrayAdd (sources, 100,source); + psFree(source); // Drop reference + + // Blow away the image parts of the source, which makes the memory explode + RESET(source->pixels); + RESET(source->weight); + RESET(source->maskObj); + RESET(source->maskView); + RESET(source->modelFlux); + RESET(source->psfFlux); + RESET(source->blends); + } + + // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT" + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); + psFree(sources); + + // XXX many leaks in here, i think + return true; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoadSpots.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoadSpots.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoadSpots.c (revision 17629) @@ -0,0 +1,68 @@ +#include "ppSim.h" + +bool ppSimLoadSpots (pmFPA *fpa, pmConfig *config) { + + bool mdok; + assert (stars); + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + bool forcedPhot = psMetadataLookupBool(&mdok, recipe, "FORCED.PHOT"); // Density of fakes + if (!forcedPhot) return true; + + // Read catalogue stars using psastro + // XXX probably need to modify this... + psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE); + if (!astroRecipe) { + psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE); + return NULL; + } + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale; + psFree(bounds); + + float x0fpa = 0.5*(bounds->x0 + bounds->x1); + float y0fpa = 0.5*(bounds->y0 + bounds->y1); + + // XXX need to use the WCS to define the overlap region + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN", PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius); + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX", PS_DATA_F32 | PS_META_REPLACE, "", ra0 + radius); + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", dec0 - radius); + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius); + psArray *refStars = psastroLoadRefstars(config); + if (!refStars || refStars->n == 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars."); + psFree(refStars); + return NULL; + } + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n); + + // convert the pmAstromObj sources to the storage format used by the CMF files + + long oldSize = stars->n; + stars = psArrayRealloc (stars, refStars->n); + + // Conversion loop + for (long i = 0; i < refStars->n; i++) { + ppSimStar *star = ppSimStarAlloc (); + + pmAstromObj *ref = refStars->data[i]; // Reference star + star->ra = ref->sky->r; // RA of star + star->dec = ref->sky->d; // Dec of star + star->mag = ref->Mag; // Magnitude of star + + // Apply rotation, make FPA center of boresite + // convert the ra,dec to x,y using examples in psastro. + // star->x = cos(pa) * xi - sin(pa) * eta + x0fpa; + // star->y = sin(pa) * xi + cos(pa) * eta + y0fpa; + + stars->data[oldSize + i] = star; + + psTrace("ppSim", 10, "Adding catalogue star: %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux); + } + + // XXX these need to be saved on PSPHOT.INPUT.CMF + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoadStars.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoadStars.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoadStars.c (revision 17629) @@ -0,0 +1,98 @@ +#include "ppSim.h" + +bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config) { + + bool mdok; + assert (stars); + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + bool starsReal = psMetadataLookupBool(&mdok, recipe, "STARS.REAL"); // Density of fakes + if (!starsReal) return true; + + // Read catalogue stars using psastro + psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE); + if (!astroRecipe) { + psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE); + return NULL; + } + + // XXX push these into the recipe in ppSimArguments() + float zp = psMetadataLookupF32(NULL, recipe, "ZEROPOINT"); // Photometric zero point + float ra0 = psMetadataLookupF32(NULL, recipe, "RA"); // Boresight RA (radians) + float dec0 = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians) + float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) + float seeing = psMetadataLookupF32(NULL, recipe, "SEEING"); // Seeing SIGMA (pixels) + float scale = psMetadataLookupF32(NULL, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time (sec) + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale; + psFree(bounds); + + float x0fpa = 0.5*(bounds->x0 + bounds->x1); + float y0fpa = 0.5*(bounds->y0 + bounds->y1); + + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN", PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius); + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX", PS_DATA_F32 | PS_META_REPLACE, "", ra0 + radius); + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", dec0 - radius); + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius); + psArray *refStars = psastroLoadRefstars(config); + if (!refStars || refStars->n == 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars."); + psFree(refStars); + return NULL; + } + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n); + + long oldSize = stars->n; + stars = psArrayRealloc (stars, refStars->n); + + // Conversion loop + for (long i = 0; i < refStars->n; i++) { + ppSimStar *star = ppSimStarAlloc (); + + pmAstromObj *ref = refStars->data[i]; // Reference star + star->ra = ref->sky->r; // RA of star + star->dec = ref->sky->d; // Dec of star + star->mag = ref->Mag; // Magnitude of star + + // Convert to x,y position on tangent plane, in pixels + float div = (sin(star->ra) * sin(dec0) + cos(star->dec) * cos(dec0) * cos(star->ra - ra0)) * scale; + float xi = cos(star->dec) * sin(star->ra - ra0) / div; + float eta = (sin(star->dec) * cos(dec0) - cos(star->dec) * sin(dec0) * cos(star->ra - ra0)) / div; + + // Apply rotation, make FPA center of boresite + star->x = cos(pa) * xi - sin(pa) * eta + x0fpa; + star->y = sin(pa) * xi + cos(pa) * eta + y0fpa; + + // Convert magnitude to peak flux + star->flux = ppSimMagToFlux (star->mag, zp) * expTime; + star->peak = ppSimStarFluxToPeak (star->flux, seeing); + + stars->data[oldSize + i] = star; + + psTrace("ppSim", 10, "Adding catalogue star: %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux); + } + stars->n = oldSize + refStars->n; + + pmLumFunc *lumfunc = psastroLuminosityFunction (refStars); + psFree(refStars); + + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0); + + if (lumfunc) { + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MIN", PS_META_REPLACE, "min valid magnitude", lumfunc->mMin); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MAX", PS_META_REPLACE, "max valid magnitude", lumfunc->mMax); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.SLOPE", PS_META_REPLACE, "log-mag histogram slope", lumfunc->slope); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.OFFSET", PS_META_REPLACE, "log-mag histogram offset", lumfunc->offset); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.PEAK", PS_META_REPLACE, "magnitude of peak bin", lumfunc->mPeak); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.NUM.PEAK", PS_META_REPLACE, "number of stars in peak bin", lumfunc->nPeak); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.SUM.PEAK", PS_META_REPLACE, "sum of stars up to peak bin", lumfunc->sPeak); + psFree (lumfunc); + } + + return stars; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17629) @@ -0,0 +1,201 @@ +#include "ppSim.h" + +bool ppSimLoop(pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(config, PS_EXIT_PROG_ERROR); + + pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, OUTPUT_FILE); // Output file + assert(file); + pmFPA *fpa = file->fpa; // FPA for file + assert(fpa); + + ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + + int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + + // Load catalogue stars + psArray *stars = psArrayAllocEmpty (1); + if (type == PPSIM_TYPE_OBJECT) { + ppSimLoadStars (stars, fpa, config); + } + + // Add random stars + if (type == PPSIM_TYPE_OBJECT) { + ppSimMakeStars (stars, fpa, config, rng); + } + + // Add random galaxies + psArray *galaxies = psArrayAllocEmpty (1); + if (type == PPSIM_TYPE_OBJECT) { + ppSimMakeGalaxies (galaxies, fpa, config, rng); + } + + pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA + + // XXX if we include the psphot analysis, we will need to activate the correct pmFPAfiles + // at the correct times. ppSim operates on PPSIM.OUTPUT and PPSIM.SOURCES + + // load any needed files (eg, input image, PSF) + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for fpa in ppSim\n"); + psFree(view); + return false; + } + + ppSimUpdateConceptsFPA (fpa, config); + + pmChip *chip; // Chip from FPA + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + + // load any needed files (eg, input image, PSF) + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d in ppSim\n", view->chip); + psFree (view); + return false; + } + + if (type == PPSIM_TYPE_OBJECT) { + if (!ppSimSetPSF (chip, config)) { + psError(PS_ERR_UNKNOWN, false, "failed IO for chip %d in ppSim\n", view->chip); + psFree (view); + return false; + } + } + + pmCell *cell; // Cell from chip + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + + // check that we are able to work with this cell (readdir must be 1) + int readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); // Read direction + if (readdir != 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "CELL.READDIR = 1 is the only supported mode."); + psFree(rng); + return false; + } + + // if no readout exists, generate a single readout (upgrade to allow multiple + // readouts?) + if (cell->readouts->n == 0) { + // Size, position and orientation of cell + int numCols = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE") / binning; + int numRows = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE") / binning; + + // generate a new readout for this cell + pmReadout *readout = pmReadoutAlloc(cell); // Readout within cell + + // TO DO: Decide if cell is to be windowed, reduce numCols, numRows appropriately + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Signal in pixels + readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels + + psImageInit (readout->image, 0.0); + psImageInit (readout->weight, 0.0); + + psFree(readout); // Drop reference + } + + psVector *biasCols = ppSimMakeBiassec (cell, config); + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + assert (readout); + + // if we have not read in a weight or generated a fake image above, we need to + // build one here + if (!readout->weight) { + if (!pmReadoutGenerateWeight(readout, true)) { + psError (PS_ERR_UNKNOWN, false, "trouble creating weight"); + return false; + } + } + + psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects + if (type == PPSIM_TYPE_OBJECT) { + expCorr = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32); + } + + psVector *biasRows = ppSimMakeBias (readout, config, rng); + if (type == PPSIM_TYPE_BIAS) goto done; + + ppSimMakeDark (readout, config); + if (type == PPSIM_TYPE_DARK) goto done; + + ppSimMakeSky (readout, expCorr, type, config); + if (type == PPSIM_TYPE_FLAT) goto done; + + if (type == PPSIM_TYPE_OBJECT) { + ppSimInsertStars (readout, expCorr, stars, config); + } + + if (type == PPSIM_TYPE_OBJECT) { + ppSimInsertGalaxies (readout, expCorr, galaxies, config); + } + + psFree(expCorr); + + done: + ppSimAddNoise(readout->image, readout->weight, cell, config, rng); + ppSimSaturate(readout, config); + + ppSimBadPixels(readout, config, rng); + + ppSimAddOverscan (readout, config, biasCols, biasRows, rng); + psFree(biasRows); + + readout->data_exists = true; + readout->parent->data_exists = true; + readout->parent->parent->data_exists = true; + } + psFree(biasCols); + + ppSimUpdateConceptsCell (cell, config); + + if (cell->hdu) { + ppSimInitHeader(config, NULL, NULL, cell); + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_IO, false, "Unable to write file."); + psFree(rng); + psFree(view); + // return PS_EXIT_SYS_ERROR; + return false; + } + } + + if (chip->hdu) { + ppSimInitHeader(config, NULL, chip, NULL); + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_IO, false, "Unable to write file."); + psFree(rng); + psFree(view); + // return PS_EXIT_SYS_ERROR; + return false; + } + + } + + psFree(stars); + psFree(galaxies); + + if (fpa->hdu) { + ppSimInitHeader(config, fpa, NULL, NULL); + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_IO, false, "Unable to write file."); + psFree(rng); + psFree(view); + // return PS_EXIT_SYS_ERROR; + return false; + } + + psFree(rng); + psFree(view); + + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c (revision 17629) @@ -0,0 +1,57 @@ +# include "ppSim.h" + +psVector *ppSimMakeBias (pmReadout *readout, pmConfig *config, const psRandom *rng) { + + bool mdok; + + pmCell *cell = readout->parent; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + float biasLevel = psMetadataLookupF32(NULL, recipe, "BIAS.LEVEL"); // Bias level + float biasRange = psMetadataLookupF32(NULL, recipe, "BIAS.RANGE"); // Bias range + int biasOrder = psMetadataLookupS32(NULL, recipe, "BIAS.ORDER"); // Bias order + + float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e + if (isnan(readnoise)) { + psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE."); + readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + return NULL; + } + } + + psImage *signal = readout->image; + psImage *variance = readout->weight; + + int numRows = signal->numRows; + int numCols = signal->numCols; + + // Polynomial for bias + psPolynomial1D *biasPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, biasOrder); + for (int j = 0; j < biasOrder + 1; j++) { + biasPoly->coeff[j] = biasRange * psRandomGaussian(rng); + } + + psVector *biasRows = psVectorAlloc(numRows, PS_TYPE_F32); // Bias value, per row + int biasOffset = 0.5 * numRows * psRandomUniform(rng); // Offset to prevent common pattern + + for (int y = 0; y < numRows; y++) { + // Adjust bias level for this row + biasRows->data.F32[y] = psPolynomial1DEval(biasPoly, (float)(y + biasOffset) / + (float)numRows - 0.5) + biasLevel; + + for (int x = 0; x < numCols; x++) { + + // Bias level + signal->data.F32[y][x] += biasRows->data.F32[y]; + variance->data.F32[y][x] += PS_SQR(readnoise); + + } + } + psFree(biasPoly); + + return biasRows; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c (revision 17629) @@ -0,0 +1,39 @@ +# include "ppSim.h" + +psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config) { + + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + psList *biassec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); // Bias regions + int numBias = psListLength(biassec); // Number of bias sections + + psVector *biasCols; + if (numBias > 0) { + biasCols = psVectorAlloc(numBias, PS_TYPE_S32); + psListIterator *iter = psListIteratorAlloc(biassec, PS_LIST_HEAD, false); // Iterator + psRegion *bias; // Bias region, from iteration + int i = 0; // Counter + while ((bias = psListGetAndIncrement(iter))) { + biasCols->data.S32[i++] = bias->x1 = bias->x0; + } + } else { + biasCols = psVectorAlloc(1, PS_TYPE_S32); + biasCols->data.S32[0] = psMetadataLookupS32(&mdok, recipe, "OVERSCAN.SIZE"); + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find OVERSCAN.SIZE in recipe."); + psFree(biasCols); + return NULL; + } + psRegion *newBias = psRegionAlloc(NAN, NAN, NAN, NAN); // New bias region, to be set later + biassec = psListAlloc(newBias); + psFree(newBias); + psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE, + "Bias sections", biassec); + psFree(biassec); + numBias = 1; + } + + return biasCols; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c (revision 17629) @@ -0,0 +1,27 @@ +# include "ppSim.h" + +// XXX add bounds to the inputs? +bool ppSimMakeDark (pmReadout *readout, pmConfig *config) { + + bool mdok; + + psImage *signal = readout->image; + psImage *variance = readout->weight; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + + for (int y = 0; y < signal->numRows; y++) { + for (int x = 0; x < signal->numCols; x++) { + + // Dark current + float darkCurrent = darkRate * expTime; // Dark current accumulated + signal->data.F32[y][x] += darkCurrent; + variance->data.F32[y][x] += darkCurrent; + } + } + return true; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeGalaxies.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeGalaxies.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeGalaxies.c (revision 17629) @@ -0,0 +1,141 @@ +# include "ppSim.h" + +bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng) { + + bool mdok; + assert (galaxies); + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + bool galaxyFake = psMetadataLookupBool(&mdok, recipe, "GALAXY.FAKE"); // Density of fakes + if (!galaxyFake) return NULL; + + float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope + float galaxyMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude + float galaxyDensity = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes + + bool galaxyGrid = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID"); // Density of fakes + int galaxyGridDX = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes + int galaxyGridDY = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes + + float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes + float galaxyRmajorMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MIN"); // Density of fakes + + float galaxyARatioMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MAX"); // Density of fakes + float galaxyARatioMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MIN"); // Density of fakes + + float galaxyThetaMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MAX"); // Density of fakes + float galaxyThetaMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MIN"); // Density of fakes + + float galaxyIndexMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MIN"); // Density of fakes + float galaxyIndexMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MAX"); // Density of fakes + + float darkRate = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time + float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point + float seeing = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix) + float scale = psMetadataLookupF32(&mdok, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float skyRate = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate + if (isnan(skyRate)) { + float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + + if (galaxyDensity <= 0) return NULL; + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + // Size of focal plane (can't I get this from the config?) + int xSize = bounds->x1 - bounds->x0; + int ySize = bounds->y1 - bounds->y0; + psFree(bounds); + + // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it + // shouldn't be too bad. + float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + psFree(bounds); + return false; + } + + // Peak fluxes: faintest and brightest levels for random galaxy + float faint = 0.1 * 10.0 * sqrt(2.0*M_PI) * seeing * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + float bright = powf(10.0, -0.4 * (galaxyMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime; + if (bright < faint) { + psLogMsg("ppSim", PS_LOG_INFO, + "Image noise is above brightest random galaxy --- no random galaxy added."); + return NULL; + } + + // Normalisation, set by the specified stellar density at the specified bright magnitude + float norm = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) / + powf(bright, galaxyLum); + + // Total number of galaxy down to the faint flux end + long num = expf(logf(norm) + galaxyLum * logf(faint)) + 0.5; + + if (galaxyGrid) { + num = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY); + + psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", num); + + float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin) / num; + float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin) / num; + float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin) / num; + float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin) / num; + + int i = 0; + + for (long iy = 0.5*galaxyGridDY; iy < ySize; iy += galaxyGridDY) { + for (long ix = 0.5*galaxyGridDX; ix < xSize; ix += galaxyGridDX) { + ppSimGalaxy *galaxy = ppSimGalaxyAlloc (); + + // make fpa center of distribution + galaxy->x = ix; + galaxy->y = iy; + + galaxy->peak = 20000; + + galaxy->index = (galaxyIndexMin + i * galaxyIndexSlope); + galaxy->Rmaj = (galaxyRmajorMin + i * galaxyRmajorSlope); + galaxy->Rmin = (galaxyARatioMin + i * galaxyARatioSlope) * galaxy->Rmaj; + galaxy->theta = (galaxyThetaMin + i * galaxyThetaSlope); + + psArrayAdd (galaxies, 100, galaxy); + i++; + } + } + } else { + + float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin); + float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin); + float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin); + float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin); + + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld galaxies down to %f mag\n", num, + -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp); + + for (long i = 0; i < num; i++) { + ppSimGalaxy *galaxy = ppSimGalaxyAlloc (); + + // make fpa center of distribution + galaxy->x = psRandomUniform(rng) * xSize; // x position + galaxy->y = psRandomUniform(rng) * ySize; // y position + + galaxy->flux = pow ((double)((i + 1) / norm), (double) (1.0 / galaxyLum)); + + galaxy->index = (psRandomUniform(rng) * galaxyIndexSlope + galaxyIndexMin); + galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope + galaxyThetaMin); + galaxy->Rmaj = (psRandomUniform(rng) * galaxyRmajorSlope + galaxyRmajorMin); + galaxy->Rmin = (psRandomUniform(rng) * galaxyARatioSlope + galaxyARatioMin) * galaxy->Rmaj; + + galaxy->flux = expf((logf(i + 1) - logf(norm)) / galaxyLum); // Peak flux + galaxy->peak = galaxy->flux / (2.0*M_PI*PS_SQR(seeing)); + // galaxy->peak = 5000; + + psArrayAdd (galaxies, 100, galaxy); + } + } + return galaxies; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 17629) @@ -0,0 +1,85 @@ +# include "ppSim.h" + +// XXX add bounds to the inputs? +bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config) { + + bool status; + + psImage *signal = readout->image; + psImage *variance = readout->weight; + + pmCell *cell = readout->parent; + pmChip *chip = cell->parent; + pmFPA *fpa = chip->parent; + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time + float flatSigma = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient + float flatRate = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate + float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate + if (isnan(skyRate)) { + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); + float scale = psMetadataLookupF32(&status, recipe, "SCALE"); assert (status); + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + + int x0Chip = psMetadataLookupS32(&status, chip->concepts, "CHIP.X0"); + int y0Chip = psMetadataLookupS32(&status, chip->concepts, "CHIP.Y0"); + int xParityChip = psMetadataLookupS32(&status, chip->concepts, "CHIP.XPARITY"); + int yParityChip = psMetadataLookupS32(&status, chip->concepts, "CHIP.YPARITY"); + + int x0Cell = psMetadataLookupS32(&status, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(&status, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(&status, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(&status, cell->concepts, "CELL.YPARITY"); + + int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + + // Correct chip offsets so that boresight is in the middle of the FPA + x0Chip -= 0.5 * (bounds->x1 - bounds->x0); + y0Chip -= 0.5 * (bounds->y1 - bounds->y0); + + for (int y = 0; y < signal->numRows; y++) { + + float yFPA = PPSIM_CELL_TO_FPA(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 / + (bounds->y1 - bounds->y0); // Relative y position in FPA + + for (int x = 0; x < signal->numCols; x++) { + float xFPA = PPSIM_CELL_TO_FPA(x, x0Cell, xParityCell, binning, x0Chip, xParityChip) * 2.0 / + (bounds->x1 - bounds->x0); // Relative x position in FPA + + // Shutter: adjust exposure time + float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0; + + // Gaussian flat-field over the FPA + float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) / + flatSigma / sqrtf(2.0 * M_PI); + + if (type == PPSIM_TYPE_FLAT) { + float flatFlux = flatRate * flatValue * realExpTime; // Flux from flat-field + signal->data.F32[y][x] += flatFlux; + variance->data.F32[y][x] += flatFlux; + continue; + } + + expCorr->data.F32[y][x] = realExpTime / expTime; + + // Sky background + float skyFlux = skyRate * flatValue * realExpTime; // Flux from sky + signal->data.F32[y][x] += skyFlux; + variance->data.F32[y][x] += skyFlux; + + // TO DO: Add fringes + + } + } + psFree(bounds); + return true; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeStars.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeStars.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeStars.c (revision 17629) @@ -0,0 +1,132 @@ +# include "ppSim.h" + +bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) { + + bool status; + assert (stars); + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + bool starsFake = psMetadataLookupBool(&status, recipe, "STARS.FAKE"); // Density of fakes + if (!starsFake) return true; + + bool starsReal = psMetadataLookupBool(&status, recipe, "STARS.REAL"); // were real stars generated? + bool matchDensity = psMetadataLookupBool(&status, recipe, "MATCH.DENSITY"); // match real star density? + + float starsAlpha = psMetadataLookupF32(&status, recipe, "STARS.LUM"); // Star luminosity func slope + float starsDensity = psMetadataLookupF32(&status, recipe, "STARS.DENSITY"); // Density of fakes + float brightMag = psMetadataLookupF32(&status, recipe, "STARS.MAG"); // Star brightest magnitude + float nSigmaLim = psMetadataLookupF32(&status, recipe, "STARS.SIGMA.LIM"); // significance of faintest stars + + float darkRate = psMetadataLookupF32(&status, recipe, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) + float scale = psMetadataLookupF32(&status, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + + bool flatLum = psMetadataLookupBool(&status,recipe, "STARS.FLAT.LUM"); // were real stars generated? + float flatNum = psMetadataLookupS32(&status, recipe, "STARS.FLAT.NUM"); // were real stars generated? + + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate + if (isnan(skyRate)) { + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + // Size of focal plane (can't I get this from the config?) + int xSize = bounds->x1 - bounds->x0; + int ySize = bounds->y1 - bounds->y0; + psFree(bounds); + + // choose reference magnitude & density to set normalization + float refMag = 0; + float refSum = 0; + if (starsReal && matchDensity) { + refMag = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.MAG.PEAK"); // Star brightest magnitude + refSum = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.SUM.PEAK"); // Star brightest magnitude + assert (status); + } else { + refMag = brightMag; + refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + } + psTrace("ppSim", 6, "refMag: %f, refSum: %f\n", refMag, refSum); + + if (refSum <= 0) return true; + + // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it + // shouldn't be too bad. + float readnoise = psMetadataLookupF32(&status, recipe, "READNOISE"); // Default read noise + if (!status) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + psFree(bounds); + return false; + } + + // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars + + // faint limit is set by detection limit, in turn set by sky noise + float skySigma = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + float skyNoise = ppSimStarSkyNoise (skySigma, seeing); + + // total flux of faintest star: + float faintCounts = nSigmaLim * skyNoise; + float faintMag = ppSimFluxToMag ((faintCounts / expTime), zp); + + float brightCounts = ppSimMagToFlux (brightMag, zp) * expTime; // Bright limit is specified by user as mag + + psTrace("ppSim", 6, "Faint limit: %f counts = %f mags\n", faintCounts, faintMag); + psTrace("ppSim", 6, "Bright limit: %f counts = %f mags\n", brightCounts, brightMag); + if (brightCounts < faintCounts) { + psLogMsg("ppSim", PS_LOG_INFO, + "Image noise is above brightest random star --- no random stars added."); + return true; + } + + // given log_10 (dN / dmag) = alpha mag + beta + // or dN / dmag = No 10^(alpha mag), then: + // N(m < Mo) = alpha * No * log(10.0) * 10^(alpha*Mo) + + // XXX this needs to handle the case of a flat distribution for efficiency testing + + // Normalization, set by the specified stellar density at the specified bright magnitude + float norm = refSum / (starsAlpha * logf(10.0) * powf(10.0, (starsAlpha * refMag))); + float normScale = norm * starsAlpha * logf(10.0); + + // Total number of stars down to the faint flux end + long nTotal = 0; + if (flatLum) { + nTotal = flatNum; + } else { + nTotal = normScale * powf (10.0, (starsAlpha * faintMag)); + } + + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars between %f and %f mag\n", nTotal, brightMag, faintMag); + + long oldSize = stars->n; + psArrayRealloc (stars, stars->n + nTotal); + + for (long i = 0; i < nTotal; i++) { + ppSimStar *star = ppSimStarAlloc (); + + // make fpa center of distribution + star->x = psRandomUniform(rng) * xSize; // x position + star->y = psRandomUniform(rng) * ySize; // y position + + float starMag = 0; + if (flatLum) { + starMag = psRandomUniform(rng) * (faintMag - brightMag) + brightMag; + } else { + starMag = PS_MIN (faintMag, PS_MAX (brightMag, (log10((i + 1.0) / normScale) / starsAlpha))); + } + + star->flux = ppSimMagToFlux (starMag, zp) * expTime; + star->peak = ppSimStarFluxToPeak (star->flux, seeing); + + stars->data[oldSize + i] = star; + } + stars->n = oldSize + nTotal; + + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimMosaicChip.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMosaicChip.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMosaicChip.c (revision 17629) @@ -0,0 +1,36 @@ +#include "ppSim.h" + +// XXX this is essentially identical to ppImageMosaicChip +bool ppSimMosaicChip(pmConfig *config, const psMaskType blankMask, const pmFPAview *view, + const char *outFile, const char *inFile) +{ + bool status; // Status of MD lookup + + pmFPAfile *in = psMetadataLookupPtr(&status, config->files, inFile); // Input file + if (!status) { + psErrorStackPrint(stderr, "Can't find required I/O file!\n"); + exit(EXIT_FAILURE); + } + + pmFPAfile *out = psMetadataLookupPtr(&status, config->files, outFile); // Output file + if (!status) { + psErrorStackPrint(stderr, "Can't find required I/O file!\n"); + exit(EXIT_FAILURE); + } + + pmChip *outChip = pmFPAviewThisChip(view, out->fpa); + pmChip *inChip = pmFPAviewThisChip(view, in->fpa); + if (!outChip->hdu && !outChip->parent->hdu) { + const char *name = psMetadataLookupStr(&status, in->fpa->concepts, "FPA.NAME"); // Name of FPA + pmFPAAddSourceFromView(out->fpa, name, view, out->format); + } + + psTrace("pmChipMosaic", 5, "mosaic chip %s to %s (xbin,ybin: %d,%d to %d,%d)\n", + in->name, out->name, in->xBin, in->yBin, out->xBin, out->yBin); + + // XXX mosaic the chip, making a deep copy. this has the side effect of making the + // output image products pure trimmed images, but also increases the memory footprint. + status = pmChipMosaic(outChip, inChip, true, blankMask); + return status; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c (revision 17629) @@ -0,0 +1,43 @@ +# include "ppSim.h" + +// In this function, we perform the psphot analysis routine for the chip-mosaicked images +bool ppSimPhotom (pmConfig *config, pmFPAview *view) { + + bool status; + pmCell *cell; + pmReadout *readout; + + psphotModelClassInit (); + + // find or define a pmFPAfile PSPHOT.INPUT + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + if (!status) { + psError(PSPHOT_ERR_CONFIG, false, "PSPHOT.INPUT I/O file is not defined"); + return false; + } + + // we make a new copy of the output chip to keep psphot from modifying the output image + // XXX is this needed for ppSim? (yes, unless we do not do photometry and make a simulated image) + pmChip *oldChip = pmFPAviewThisChip (view, input->src); + pmChip *newChip = pmFPAviewThisChip (view, input->fpa); + pmChipCopy (newChip, oldChip); + + // iterate over the cells and readout for this chip + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("ppSimPhotom", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (! cell->process || ! cell->file_exists) { continue; } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // run the actual photometry analysis + if (!psphotReadout (config, view)) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + return false; + } + } + } + + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimSaturate.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSaturate.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSaturate.c (revision 17629) @@ -0,0 +1,34 @@ +# include "ppSim.h" + +// Apply saturation limit to image +bool ppSimSaturate(pmReadout *readout, // Image to apply saturation + const pmConfig *config // configuration data + ) +{ + bool mdok; + + psImage *image = readout->image; + pmCell *cell = readout->parent; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + float saturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION"); + if (isnan(saturation)) { + psWarning("CELL.SATURATION is not set; reverting to recipe value SATURATION."); + saturation = psMetadataLookupF32(&mdok, recipe, "SATURATION"); + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find SATURATION in recipe."); + return false; + } + } + + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + if (image->data.F32[y][x] > saturation) { + image->data.F32[y][x] = saturation; + } + } + } + return true; +} + Index: /branches/eam_branch_20080511/ppSim/src/ppSimSequence.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSequence.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSequence.c (revision 17629) @@ -0,0 +1,169 @@ +# include "ppSimSequence.h" +# include + +int main (int argc, char **argv) { + + bool status; + int argNum; + char line[1024]; + unsigned int nFail; + + psLibInit(NULL); + + char *dbname = NULL; + if ((argNum = psArgumentGet (argc, argv, "-dbname"))) { + psArgumentRemove(argNum, &argc, argv); + dbname = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } + + char *path = NULL; + if ((argNum = psArgumentGet (argc, argv, "-path"))) { + psArgumentRemove(argNum, &argc, argv); + path = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + sprintf (line, "mkdir -p %s", path); + system (line); + } + + char *workdir = NULL; + if ((argNum = psArgumentGet (argc, argv, "-workdir"))) { + psArgumentRemove(argNum, &argc, argv); + workdir = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } + + char *dvodb = NULL; + if ((argNum = psArgumentGet (argc, argv, "-dvodb"))) { + psArgumentRemove(argNum, &argc, argv); + dvodb = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } + + char *tess_id = NULL; + if ((argNum = psArgumentGet (argc, argv, "-tess_id"))) { + psArgumentRemove(argNum, &argc, argv); + tess_id = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } + + char *basename = NULL; + if ((argNum = psArgumentGet (argc, argv, "-basename"))) { + psArgumentRemove(argNum, &argc, argv); + basename = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } else { + basename = psStringCopy ("simtest"); + } + + char *label = NULL; + if ((argNum = psArgumentGet (argc, argv, "-label"))) { + psArgumentRemove(argNum, &argc, argv); + label = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } + + if (argc != 4) { + fprintf (stderr, "USAGE: ppSimSequence (sequence) (simulate) (inject) [options]\n"); + fprintf (stderr, "generates a set of simulated data defined by the sequence file\n"); + fprintf (stderr, " (sequence) : a mdc-file describing the desired image sequences\n"); + fprintf (stderr, " (simulate) : an output file with commands to generate the images\n"); + fprintf (stderr, " (inject) : an output file with commands to inject the images into the pipeline\n"); + fprintf (stderr, "options:\n"); + fprintf (stderr, " -dbname (dbname)\n"); + fprintf (stderr, " -path (path)\n"); + fprintf (stderr, " -workdir (workdir)\n"); + fprintf (stderr, " -basename (basename)\n"); + fprintf (stderr, " -label (label)\n"); + fprintf (stderr, " -dvodb (dvodb)\n"); + fprintf (stderr, " -tess_id (tess_id)\n"); + exit (2); + } + + // load the sequence description + psMetadata *config = psMetadataConfigRead (NULL, &nFail, argv[1], false); + if (!config) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "unable to read sequence description from %s", argv[1]); + exit (1); + } + + FILE *simfile = fopen (argv[2], "w"); + if (!simfile) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "unable to open %s for output", argv[2]); + exit (1); + } + + FILE *inject = fopen (argv[3], "w"); + if (!inject) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "unable to open %s for output", argv[3]); + exit (1); + } + + // build the base injectCommand string + psString injectCommand = psStringCopy ("ipp_inject_fileset.pl --telescope SIMTEST"); + if (dbname) psStringAppend (&injectCommand, " --dbname %s", dbname); + if (workdir) psStringAppend (&injectCommand, " --workdir %s", workdir); + if (label) psStringAppend (&injectCommand, " --label %s", label); + if (dvodb) psStringAppend (&injectCommand, " --dvodb %s", dvodb); + if (tess_id) psStringAppend (&injectCommand, " --tess_id %s", tess_id); + + unsigned long seed = psMetadataLookupS32 (&status, config, "RND_SEED"); + if (!status) seed = 0; + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed); + + psMetadataItem *item = psMetadataLookup (config, "SEQUENCE"); + if (item == NULL) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "missing SEQUENCE description"); + exit (1); + } + + psArray *sequences = NULL; + if (item->type == PS_DATA_METADATA) { + sequences = psArrayAlloc(1); + sequences->data[0] = psMemIncrRefCounter (item->data.V); + } else { + if (item->type != PS_DATA_METADATA_MULTI) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "SEQUENCE is not MULTI or METADATA"); + exit (1); + } + sequences = psListToArray (item->data.list); + } + + for (int i = 0; i < sequences->n; i++) { + + // XXX this is not obvious: the entry on the list is a metadata item + // containing the psMetadata :: why is this not just a metadata? + psMetadataItem *item = sequences->data[i]; + psMetadata *sequence = item->data.V; + // double check item type? + + // determine the sequence type + char *type = psMetadataLookupStr (&status, sequence, "OBSTYPE"); + if (!status) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "SEQUENCE %d is missing a type", i); + exit (1); + } + + if (!strcasecmp (type, "BIAS")) { + ppSimSequenceBias (simfile, inject, sequence, i, rng, path, basename, injectCommand); + continue; + } + if (!strcasecmp (type, "DARK")) { + ppSimSequenceDark (simfile, inject, sequence, i, rng, path, basename, injectCommand); + continue; + } + if (!strcasecmp (type, "FLAT")) { + ppSimSequenceFlat (simfile, inject, sequence, i, rng, path, basename, injectCommand); + continue; + } + if (!strcasecmp (type, "OBJECT")) { + ppSimSequenceObject (simfile, inject, sequence, i, rng, path, basename, injectCommand); + continue; + } + + psLogMsg ("ppSimSequence", PS_LOG_WARN, "SEQUENCE %d has an unknown type: %s", i, type); + exit (1); + } + + exit (0); +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimSequence.h =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSequence.h (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSequence.h (revision 17629) @@ -0,0 +1,19 @@ +#ifndef PP_SIM_SEQUENCE_H +#define PP_SIM_SEQUENCE_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +bool ppSimSequenceBias (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand); +bool ppSimSequenceDark (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand); +bool ppSimSequenceFlat (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand); +bool ppSimSequenceObject (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand); + +#endif Index: /branches/eam_branch_20080511/ppSim/src/ppSimSequenceBias.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSequenceBias.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSequenceBias.c (revision 17629) @@ -0,0 +1,50 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceBias (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand) { + + bool status, setLevel, setRange; + + char *camera = psMetadataLookupStr (&status, sequence, "CAMERA"); + + // optional details + float level = psMetadataLookupF32 (&setLevel, sequence, "BIAS.LEVEL"); + float range = psMetadataLookupF32 (&setRange, sequence, "BIAS.RANGE"); + + int nImages = psMetadataLookupS32 (&status, sequence, "NIMAGES"); + + // loop over the filters & exposure times + int nImage = 0; + for (int i = 0; i < nImages; i++) { + + psString command = NULL; + + psStringAppend (&command, "ppSim -type BIAS"); + psStringAppend (&command, " -camera %s", camera); + + if (setLevel) psStringAppend (&command, " -biaslevel %f", level); + if (setRange) psStringAppend (&command, " -biasrange %f", range); + + // XXX need to add output filename + psString filename = NULL; + if (path) { + psStringAppend (&filename, "%s/%s.%03d.%03d", path, basename, nSeq, nImage); + } else { + psStringAppend (&filename, "%s.%03d.%03d", basename, nSeq, nImage); + } + psStringAppend (&command, " %s", filename); + + fprintf (simfile, "%s\n", command); + psFree (command); + + // path should be dirname/filename + command = psStringCopy (injectCommand); + psStringAppend (&command, " --camera %s", camera); + psStringAppend (&command, " %s*.fits", filename); + fprintf (inject, "%s\n", command); + psFree (command); + psFree (filename); + + nImage ++; + } + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimSequenceDark.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSequenceDark.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSequenceDark.c (revision 17629) @@ -0,0 +1,67 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceDark (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand) { + + bool status, setRate; + float min, max = 0; + + char *camera = psMetadataLookupStr (&status, sequence, "CAMERA"); + + setRate = false; + min = psMetadataLookupF32 (&status, sequence, "DARK.MIN"); + if (status) { + max = psMetadataLookupF32 (&status, sequence, "DARK.MAX"); + setRate = true; + } + + psVector *exptimes = psMetadataLookupPtr (&status, sequence, "EXPTIMES"); + psVector *nImages = psMetadataLookupPtr (&status, sequence, "NIMAGES"); + + assert (exptimes->n == nImages->n); + + // loop over the filters & exposure times + int nImage = 0; + for (int i = 0; i < nImages->n; i++) { + + float exptime = exptimes->data.F32[i]; + float n = nImages->data.S32[i]; + + for (int j = 0; j < n; j++) { + psString command = NULL; + + psStringAppend (&command, "ppSim -type DARK"); + psStringAppend (&command, " -camera %s", camera); + + if (setRate) { + double frnd = psRandomUniform(rng); + float rate = min + (max - min)*frnd; + psStringAppend (&command, " -darkrate %f", rate); + } + + psStringAppend (&command, " -exptime %f", exptime); + + // XXX need to add output filename + psString filename = NULL; + if (path) { + psStringAppend (&filename, "%s/%s.%03d.%03d", path, basename, nSeq, nImage); + } else { + psStringAppend (&filename, "%s.%03d.%03d", basename, nSeq, nImage); + } + psStringAppend (&command, " %s", filename); + + fprintf (simfile, "%s\n", command); + psFree (command); + + // path should be dirname/filename + command = psStringCopy (injectCommand); + psStringAppend (&command, " --camera %s", camera); + psStringAppend (&command, " %s*.fits", filename); + fprintf (inject, "%s\n", command); + psFree (command); + psFree (filename); + + nImage ++; + } + } + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimSequenceFlat.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSequenceFlat.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSequenceFlat.c (revision 17629) @@ -0,0 +1,64 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceFlat (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand) { + + bool status; + + char *camera = psMetadataLookupStr (&status, sequence, "CAMERA"); + + // determine the filters & exposure times + char *filterList = psMetadataLookupStr (&status, sequence, "FILTERS"); + psArray *filters = psStringSplitArray (filterList, ",: ", false); + + psVector *exptimes = psMetadataLookupPtr (&status, sequence, "EXPTIMES"); + + if (filters->n != exptimes->n) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "mis-match in filter and exptime lists"); + exit (1); + } + + // number of images for each filter, exptime set + int nSetup = psMetadataLookupS32 (&status, sequence, "NSETUP"); + + int nImage = 0; + + // loop over the filters & exposure times + for (int i = 0; i < filters->n; i++) { + + // loop over the filters & exposure times + for (int j = 0; j < nSetup; j++) { + + psString command = NULL; + + psStringAppend (&command, "ppSim -type FLAT"); + psStringAppend (&command, " -camera %s", camera); + + psStringAppend (&command, " -filter %s", (char *) filters->data[i]); + psStringAppend (&command, " -exptime %f", exptimes->data.F32[i]); + + // XXX need to add output filename + psString filename = NULL; + if (path) { + psStringAppend (&filename, "%s/%s.%03d.%03d", path, basename, nSeq, nImage); + } else { + psStringAppend (&filename, "%s.%03d.%03d", basename, nSeq, nImage); + } + psStringAppend (&command, " %s", filename); + + + fprintf (simfile, "%s\n", command); + psFree (command); + + // path should be dirname/filename + command = psStringCopy (injectCommand); + psStringAppend (&command, " --camera %s", camera); + psStringAppend (&command, " %s*.fits", filename); + fprintf (inject, "%s\n", command); + psFree (command); + psFree (filename); + + nImage ++; + } + } + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimSequenceObject.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSequenceObject.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSequenceObject.c (revision 17629) @@ -0,0 +1,125 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceObject (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *injectCommand) { + + bool status; + + // generate ppSim lines that look like: + // ppSim -camera $camera -type OBJECT -filter $filter -exptime $exptime + // -skyrate $sky -ra $ra -dec $dec -pa $pa -scale $scale -zp $zp -seeing $seeing $filename", + + char *camera = psMetadataLookupStr (&status, sequence, "CAMERA"); + + // sequence reference coordinate + float Ro = psMetadataLookupF32 (&status, sequence, "CENTER.RA"); + float Do = psMetadataLookupF32 (&status, sequence, "CENTER.DEC"); + + // determine the filters & exposure times + char *filterList = psMetadataLookupStr (&status, sequence, "FILTERS"); + psArray *filters = psStringSplitArray (filterList, ",: ", false); + + psVector *exptimes = psMetadataLookupPtr (&status, sequence, "EXPTIMES"); + + psVector *skymags = psMetadataLookupPtr (&status, sequence, "SKYMAGS"); + + float IQmin = psMetadataLookupF32 (&status, sequence, "IQ_MIN"); + float IQmax = psMetadataLookupF32 (&status, sequence, "IQ_MAX"); + + if (filters->n != exptimes->n) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "mis-match in filter and exptime lists"); + exit (1); + } + + // track the number of files produced + int nImage = 0; + + // loop over the filters & exposure times + for (int i = 0; i < filters->n; i++) { + + // offset parameters + float dR = psMetadataLookupF32 (&status, sequence, "OFFSET.RA"); + float dD = psMetadataLookupF32 (&status, sequence, "OFFSET.DEC"); + + int nR = psMetadataLookupS32 (&status, sequence, "OFFSET.NR"); + int nD = psMetadataLookupS32 (&status, sequence, "OFFSET.ND"); + + // loop over the offset sequence + for (int iR = 0; iR < nR; iR++) { + for (int iD = 0; iD < nD; iD++) { + + // RA & DEC in degrees XXX (should be radians...) + // offsets are in arcseconds + float R = Ro + dR*(iR - 0.5*nR + 0.5) / cos (RAD_DEG*Do) / 3600.0; + float D = Do + dD*(iD - 0.5*nD + 0.5) / 3600.0; + + // dither parameters + float dr = psMetadataLookupF32 (&status, sequence, "DITHER.RA"); + float dd = psMetadataLookupF32 (&status, sequence, "DITHER.DEC"); + + int nr = psMetadataLookupS32 (&status, sequence, "DITHER.NR"); + int nd = psMetadataLookupS32 (&status, sequence, "DITHER.ND"); + + // loop over the dither sequence + for (int ir = 0; ir < nr; ir++) { + for (int id = 0; id < nd; id++) { + + // ra, dec in degrees; offsets in arcsec + float ra = R + dr*(ir - 0.5*nr + 0.5) / cos (RAD_DEG*D) / 3600.0; + float dec = D + dd*(id - 0.5*nd + 0.5) / 3600.0; + + // rotation sequence parameters + float pos_min = psMetadataLookupF32 (&status, sequence, "POS_MIN"); + float pos_max = psMetadataLookupF32 (&status, sequence, "POS_MAX"); + float pos_delta = psMetadataLookupF32 (&status, sequence, "POS_DELTA"); + assert (pos_delta > 0.0); + assert (pos_max >= pos_min); + + // loop over rotation sequence + for (float pos = pos_min; pos <= pos_max; pos += pos_delta) { + + psString command = NULL; + + psStringAppend (&command, "ppSim -type OBJECT"); + psStringAppend (&command, " -camera %s", camera); + psStringAppend (&command, " -filter %s", (char *) filters->data[i]); + psStringAppend (&command, " -exptime %f", exptimes->data.F32[i]); + psStringAppend (&command, " -skymags %f", skymags->data.F32[i]); + + psStringAppend (&command, " -ra %f", ra); + psStringAppend (&command, " -dec %f", dec); + psStringAppend (&command, " -pa %f", pos); + + double frnd = psRandomUniform(rng); + float seeing = IQmin + (IQmax - IQmin)*frnd; + + psStringAppend (&command, " -seeing %f", seeing); + + // XXX need to add output filename + psString filename = NULL; + if (path) { + psStringAppend (&filename, "%s/%s.%03d.%03d", path, basename, nSeq, nImage); + } else { + psStringAppend (&filename, "%s.%03d.%03d", basename, nSeq, nImage); + } + psStringAppend (&command, " %s", filename); + + fprintf (simfile, "%s\n", command); + psFree (command); + + // path should be dirname/filename + command = psStringCopy (injectCommand); + psStringAppend (&command, " --camera %s", camera); + psStringAppend (&command, " %s*.fits", filename); + fprintf (inject, "%s\n", command); + psFree (command); + psFree (filename); + + nImage ++; + } + } + } + } + } + } + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimSetPSF.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimSetPSF.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimSetPSF.c (revision 17629) @@ -0,0 +1,85 @@ +# include "ppSim.h" +static char *defaultModel = "PS_MODEL_QGAUSS"; + +bool ppSimSetPSF (pmChip *chip, pmConfig *config) { + + bool status, mdok; + pmPSF *psf = NULL; + pmTrend2D *param = NULL; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + // the pmPSF IO functions stores the PSF on the chip->analysis + psf = psMetadataLookupPtr (&status, chip->analysis, "PSPHOT.PSF"); + if (psf) { + return true; + } + + // no supplied PSF, build one using supplied value for seeing seeing is already corrected + // for the pixel scale, and is converted from FWHM to SIGMA (this is done in + // ppSimArguments) + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) + + char *psfModelName = psMetadataLookupStr(&status, config->arguments, "PSF.MODEL"); // Name of PSF model + if (psfModelName == NULL) { + psfModelName = defaultModel; + } + + // structure to store user options defining the psf + pmPSFOptions *options = pmPSFOptionsAlloc (); + options->type = pmModelClassGetType (psfModelName); + if (options->type == -1) { + psError (PS_ERR_UNKNOWN, false, "invalid model name"); + return false; + } + + // XXX this is messed up: CHIP.XSIZE and CHIP.YSIZE are not seta + int xSize = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XSIZE"); + int ySize = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YSIZE"); + xSize = 1000; + ySize = 1000; + + // no spatial variation + options->psfTrendMode = PM_TREND_POLY_ORD; + options->psfTrendNx = 0; + options->psfTrendNy = 0; + options->psfFieldNx = xSize; + options->psfFieldNy = ySize; + + // generate the psf + psf = pmPSFAlloc (options); + + psEllipseAxes axes; + psEllipsePol pol; + + // supply the semi-major axis (these are SIGMA values in PIXELS) + axes.major = seeing; + axes.minor = seeing; + axes.theta = 0.0; + + pol = psEllipseAxesToPol (axes); + + param = psf->params->data[PM_PAR_E0]; + param->poly->coeff[0][0] = pol.e0; + + param = psf->params->data[PM_PAR_E1]; + param->poly->coeff[0][0] = pol.e1; + + param = psf->params->data[PM_PAR_E2]; + param->poly->coeff[0][0] = pol.e2; + + if (!strcasecmp (psfModelName, "PS_MODEL_QGAUSS")) { + param = psf->params->data[PM_PAR_7]; + param->poly->coeff[0][0] = 1.0; + } + + if (!strcasecmp (psfModelName, "PS_MODEL_RGAUSS")) { + param = psf->params->data[PM_PAR_7]; + param->poly->coeff[0][0] = 1.0; + } + + psMetadataAdd (chip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); + psFree(psf); // Drop reference + + return true; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimStars.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimStars.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimStars.c (revision 17629) @@ -0,0 +1,59 @@ +# include "ppSim.h" + +void ppSimStarFree(ppSimStar *star) +{ + return; +} + +ppSimStar *ppSimStarAlloc () { + + ppSimStar *star = (ppSimStar *) psAlloc(sizeof(ppSimStar)); + psMemSetDeallocator(star, (psFreeFunc) ppSimStarFree); + + return star; +} + +void ppSimGalaxyFree(ppSimGalaxy *galaxy) +{ + return; +} + +ppSimGalaxy *ppSimGalaxyAlloc () { + + ppSimGalaxy *galaxy = (ppSimGalaxy *) psAlloc(sizeof(ppSimGalaxy)); + psMemSetDeallocator(galaxy, (psFreeFunc) ppSimGalaxyFree); + + return galaxy; +} + +float ppSimStarSkyNoise (float skySigma, float seeingSigma) { + + float skyNoise = skySigma * sqrt(4*M_PI*PS_SQR(seeingSigma)); + return skyNoise; +} + +float ppSimStarPeakToFlux (float peak, float seeingSigma) { + + float psfArea = 2.0*M_PI*PS_SQR(seeingSigma); + float flux = peak * psfArea; + return flux; +} + +float ppSimStarFluxToPeak (float flux, float seeingSigma) { + + float psfArea = 2.0*M_PI*PS_SQR(seeingSigma); + float peak = flux / psfArea; + return peak; +} + +float ppSimFluxToMag (float flux, float zp) { + + float mag = -2.5*log10(flux) + zp; + return mag; +} + +float ppSimMagToFlux (float mag, float zp) { + + float flux = powf (10.0, -0.4*(mag - zp)); + return flux; +} Index: /branches/eam_branch_20080511/ppSim/src/ppSimUtils.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimUtils.c (revision 17629) +++ /branches/eam_branch_20080511/ppSim/src/ppSimUtils.c (revision 17629) @@ -0,0 +1,154 @@ +# include "ppSim.h" + +// Generate a header containing WCS keywords +// this function is called with only one of fpa, chip, cell not NULL +bool ppSimInitHeader(pmConfig *config, + pmFPA *fpa, + pmChip *chip, + pmCell *cell) +{ + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + float ra0 = psMetadataLookupF32(NULL, recipe, "RA"); // Boresight RA (radians) + float dec0 = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians) + float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) + float scale = psMetadataLookupF32(NULL, recipe, "SCALE"); // plate scale in arcsec / pixel + scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel + + int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + + float x0 = 0.0, y0 = 0.0; + int xParity = 0, yParity = 0; + if (cell) { + int x0Chip = psMetadataLookupS32(NULL, cell->parent->concepts, "CHIP.X0"); + int y0Chip = psMetadataLookupS32(NULL, cell->parent->concepts, "CHIP.Y0"); + int xParityChip = psMetadataLookupS32(NULL, cell->parent->concepts, "CHIP.XPARITY"); + int yParityChip = psMetadataLookupS32(NULL, cell->parent->concepts, "CHIP.YPARITY"); + + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + + x0 = PPSIM_FPA_TO_CELL(0.0, x0Cell, xParityCell, binning, x0Chip, xParityChip); + y0 = PPSIM_FPA_TO_CELL(0.0, y0Cell, yParityCell, binning, y0Chip, yParityChip); + xParity = xParityCell * xParityChip; + yParity = yParityCell * yParityChip; + } + if (chip) { + int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0"); + int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0"); + int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY"); + int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY"); + + x0 = PPSIM_FPA_TO_CELL(0.0, 0, 1, binning, x0Chip, xParityChip); + y0 = PPSIM_FPA_TO_CELL(0.0, 0, 1, binning, y0Chip, yParityChip); + xParity = xParityChip; + yParity = yParityChip; + } + if (fpa) { + psRegion *bounds = ppSimFPABounds (fpa); + x0 = 0.5 * (bounds->x1 - bounds->x0); + y0 = 0.5 * (bounds->y1 - bounds->y0); + xParity = 1; + yParity = 1; + psFree (bounds); + } + assert(xParity != 0 && yParity != 0); + + psMetadata *header = psMetadataAlloc(); // Header, to return + pmAstromWCS *wcs = pmAstromWCSAlloc(1, 1); // WCS structure + wcs->toSky = psProjectionAlloc(ra0, dec0, scale * xParity, scale * yParity, PS_PROJ_TAN); + wcs->cdelt1 = scale * PM_DEG_RAD * xParity; + wcs->cdelt2 = scale * PM_DEG_RAD * yParity; + wcs->crpix1 = x0; + wcs->crpix2 = y0; + wcs->trans->x->coeff[1][0] = cos(pa) * wcs->cdelt1; + wcs->trans->x->coeff[0][1] = -sin(pa) * wcs->cdelt1; + wcs->trans->y->coeff[1][0] = sin(pa) * wcs->cdelt2; + wcs->trans->y->coeff[0][1] = cos(pa) * wcs->cdelt2; + + // These aren't used by pmAstromWCStoHeader, but set them anyway + wcs->trans->x->coeff[0][0] = ra0; + wcs->trans->y->coeff[0][0] = dec0; + wcs->trans->x->coeff[1][1] = 0.0; + wcs->trans->y->coeff[1][1] = 0.0; + + pmAstromWCStoHeader(header, wcs); + psFree(wcs); + + if (cell) { + cell->hdu->header = header; + cell->data_exists = true; + } + if (chip) { + chip->hdu->header = header; + chip->data_exists = true; + } + if (fpa) { + fpa->hdu->header = header; + } + + return true; +} + +bool ppSimUpdateConceptsFPA (pmFPA *fpa, pmConfig *config) { + + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + + const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name + if (!filter) { + filter = "NONE"; + } + + ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate + + // Update FPA concepts + const char *typeStr; // Exposure type String + switch (type) { + case PPSIM_TYPE_BIAS: typeStr = "BIAS"; break; + case PPSIM_TYPE_DARK: typeStr = "DARK"; break; + case PPSIM_TYPE_FLAT: typeStr = "FLAT"; break; + case PPSIM_TYPE_OBJECT: typeStr = "OBJECT"; break; + default: + psAbort("Should never get here."); + } + + psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.OBSTYPE", PS_META_REPLACE, + "Observation type", typeStr); + psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.OBJECT", PS_META_REPLACE, + "Observation name", typeStr); + psMetadataAddF32(fpa->concepts, PS_LIST_TAIL, "FPA.EXPTIME", PS_META_REPLACE, + "Exposure time (sec)", expTime); + psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.FILTERID", PS_META_REPLACE, "Filter name", filter); + psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.FILTER", PS_META_REPLACE, "Filter name", filter); + + return true; +} + +bool ppSimUpdateConceptsCell (pmCell *cell, pmConfig *config) { + + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + + psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, + "Exposure time (sec)", expTime); + psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.DARKTIME", PS_META_REPLACE, + "Dark time (sec)", expTime); + psMetadataAddS32(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE, + "Binning in x", binning); + psMetadataAddS32(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, + "Binning in y", binning); + + return true; +}