Index: /branches/eam_branches/eam_branch_20080511/ppSim/.cvsignore =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/.cvsignore (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/.cvsignore (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/Makefile.am =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/Makefile.am (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/Makefile.am (revision 22153) @@ -0,0 +1,3 @@ +SUBDIRS = src + +CLEANFILES = *~ core core.* Index: /branches/eam_branches/eam_branch_20080511/ppSim/autogen.sh =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/autogen.sh (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/autogen.sh (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/configure.ac =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/configure.ac (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/configure.ac (revision 22153) @@ -0,0 +1,31 @@ +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([PSPHOT], [psphot >= 0.8.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_branches/eam_branch_20080511/ppSim/notes.txt =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/notes.txt (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/notes.txt (revision 22153) @@ -0,0 +1,128 @@ + + +2008.05.17 + + pmFPAfiles used in ppSim: + + PPSIM.INPUT : an input image loaded by ppSim; the simulated features + may be added to this image. In fact, the features are + generated in the output image, and the two images are + added (as images) together. + + PPSIM.OUTPUT : the output image. this file is used for the main + ppSimLoop function + + PPSIM.CHIP : chip-mosaiced version of the output image. only + generated if photometry is requested. Is this + actually used as an output image? + + PPSIM.REAL.SOURCES : real sources loaded by PPSIM. these sources + are subtracted from the image before testing + for the detectability of the fake sources and + before measuring the flux of fake or forced + soures. They should be the sources measured in + a previous psphot run, though any input list + could be used. Results are very ill-defined if + the sources do not correspond to actual image + sources! + + PSPHOT.PSF.LOAD : input PSF model used to measure the input sources + + PPSIM.SOURCES : output fake sources (injected fake sources with real values) + PPSIM.FAKE.SOURCES : output fake photometry (measured photometry for all fake sources) + PPSIM.FORCE.SOURCES : output force photometry + + ** used within psphot functions: + PSPHOT.BACKMDL : model background generated by psphotModelBackground + PSPHOT.BACKMDL.STDEV : model background error generated by psphotModelBackground + PSPHOT.BACKGND : full-scale model background generated by psphotSubtractBackground + PSPHOT.BACKSUB : background-subtracted image + +2008.05.15 + + For fake and force photometry, we need to load the known existing + sources, then perform a complete linear fit solution to the complete + set of both real (known) sources plus the fake and/or forced + photometry positions (as PSFs). The fake and forced photometry + positions need to be measured independently. This provides the real + photometry for these sources ** in the presence of the other real + sources **. In order to test the detectability of the fake sources + (or the forced sources, for that matter), we need to generate the + smoothed detection image, and determine the peak value at the + positions of the sources. + + ppSimPhotomReadout outline: + + * we have three source lists (loaded before the function is called): + * realSources : these are loaded from a CMF-style file (or equiv) + * fakeSources : these are defined internally, and need a link from + the input fake source to the measured fake source + * forceSources : these are loaded from a DVO database via getstar + + * we require a supplied PSF + + * need to subtract the background (before or after subtracting the + sources?) + + * subtract the real sources + * model the background + * build the detection image + * measure the detectability of the fake and force sources + * replace the real sources + + * merge real + fake sources + * linear fit to merged source list + * replace the sources (make 'no-subtract' option?) + + * merge real + forced sources + * linear fit to merged source list + * replace the sources (make 'no-subtract' option?) + + + +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_branches/eam_branch_20080511/ppSim/psf.txt =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/psf.txt (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/psf.txt (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/.cvsignore =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/.cvsignore (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/.cvsignore (revision 22153) @@ -0,0 +1,9 @@ +config.h +config.h.in +Makefile +Makefile.in +ppSim +stamp-h1 +.deps +.libs +ppSimSequence Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/Makefile.am =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 22153) @@ -0,0 +1,59 @@ +bin_PROGRAMS = ppSim ppSimSequence + +ppSim_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PSASTRO_CFLAGS) $(ppSim_CFLAGS) +ppSim_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSPHOT_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 \ + ppSimLoadSpots.c \ + ppSimPhotom.c \ + ppSimPhotomReadoutFake.c \ + ppSimPhotomReadoutForce.c \ + ppSimPhotomFiles.c \ + ppSimLoadForceSources.c \ + ppSimMergeReadouts.c \ + ppSimDetections.c \ + ppSimMergeSources.c \ + ppSimMosaicChip.c \ + ppSimRandomGaussian.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_branches/eam_branch_20080511/ppSim/src/ppSim.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSim.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSim.c (revision 22153) @@ -0,0 +1,42 @@ +# include "ppSim.h" + +int failure (pmConfig *config, psExit status, char *message) { + ppSimRandomGaussianFree (); + psErrorStackPrint(stderr, message); + psFree(config); + pmModelClassCleanup(); + psLibFinalize(); + exit (status); +} + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + if (!pmModelClassInit ()) abort(); + + pmConfig *config = pmConfigRead(&argc, argv, PPSIM_RECIPE); // Configuration + if (!config) { + failure (config, PS_EXIT_CONFIG_ERROR, "Unable to read configurations."); + } + + if (!ppSimArguments(argc, argv, config)) { + failure (config, PS_EXIT_CONFIG_ERROR, "Error parsing command-line arguments"); + } + + if (!ppSimCreate(config)) { + failure (config, PS_EXIT_CONFIG_ERROR, "Unable to create output file."); + } + + if (!ppSimLoop(config)) { + failure (config, PS_EXIT_SYS_ERROR, "Unable to generate data."); + } + + ppSimRandomGaussianFree (); + psFree(config); + pmModelClassCleanup(); + pmConfigDone(); + pmConceptsDone(); + psLibFinalize(); + + exit (PS_EXIT_SUCCESS); +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSim.h =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 22153) @@ -0,0 +1,201 @@ +#ifndef PP_SIM_H +#define PP_SIM_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include // for unlink +#include +#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 +bool 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 (bool *status, pmReadout *readout, pmConfig *config, const psRandom *rng) ; +bool ppSimMakeDark (pmReadout *readout, pmConfig *config); +bool ppSimMakeSky (pmReadout *readout, psImage *expCorr, ppSimType type, pmConfig *config); + +bool ppSimLoadSpots (pmFPA *fpa, 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); + +bool ppSimMosaicChip(pmConfig *config, const psMaskType blankMask, const pmFPAview *view, + const char *outFile, const char *inFile); + +bool ppSimPhotom (pmConfig *config, pmFPAview *view); + + +/// 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); + +float ppSimArgToRecipeF32(bool *status, + 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 + ); + +int ppSimArgToRecipeS32(bool *status, + 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 + ); + +char *ppSimArgToRecipeStr(bool *status, + 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 ppSimArgToRecipeBool(bool *status, + 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 + ); + +ppSimType ppSimTypeFromString (char *typeStr); +char *ppSimTypeToString (ppSimType type); + +float ppSimGetZeroPoint (psMetadata *recipe, char *filter); + +bool ppSimMergeReadouts (pmConfig *config, pmFPAview *view); + +double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma); +double ppSimRandomGaussianNorm (const psRandom *rnd); +void ppSimRandomGaussianFree(); + +bool ppSimPhotomFiles (pmConfig *config, pmFPAfile *fakeFile, pmFPAfile *forceFile); + +bool ppSimPhotomReadoutFake(pmConfig *config, const pmFPAview *view); +bool ppSimPhotomReadoutForce(pmConfig *config, const pmFPAview *view); + +psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view); +bool ppSimDetections (psImage *significance, psMetadata *recipe, psArray *sources); +psArray *ppSimMergeSources (psArray *in1, psArray *in2); + +psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename); +bool ppSimDefinePixels (psArray *sources, pmReadout *readout, psMetadata *recipe); + +#endif Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 22153) @@ -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]) * ppSimRandomGaussianNorm(rng); + signal->data.F32[y][x] /= gain; // Converting to ADU + } + } + + // XXX why return this?? + return true; +} + Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 22153) @@ -0,0 +1,42 @@ +#include "ppSim.h" + +bool ppSimAddOverscan (pmReadout *readout, pmConfig *config, psVector *biasCols, psVector *biasRows, psRandom *rng) { + + bool mdok; + + // skip this step if we did not generate a bias level + if (biasCols == NULL) return true; + if (biasRows == NULL) return true; + + 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_branches/eam_branch_20080511/ppSim/src/ppSimArguments.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 22153) @@ -0,0 +1,229 @@ +# 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); +} + +// this function supplements the RECIPE:OPTIONS folder with command-line options +bool ppSimArguments(int argc, char **argv, pmConfig *config) +{ + bool status; + + assert(config); + + // save the following command-line options in the arguments structure. these will later be + // parsed and moved to the config->recipes:PPSIM_RECIPE folder + + 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); + psMetadataAddBool(arguments, PS_LIST_TAIL, "+photom", 0, "Perform photometry on fake sources", false); + + 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"); + if (psErrorCodeLast () != PS_ERR_NONE) { psAbort ("problem with -psf or -psflist option"); } + + pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT.SOURCES", "-cmf", "-cmflist"); + if (psErrorCodeLast () != PS_ERR_NONE) { psAbort ("problem with -cmf or -cmflist option"); } + + // Only one of -camera and -file is needed or allowed. The -camera option would have been + // already parsed by pmConfigRead in ppSim.c and resulted in a value for config->camera + bool loadImage = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-input", "-inputlist"); + if (!config->camera && !loadImage) { + psError(PS_ERR_IO, true, "A camera name (-camera NAME) or an input image (-input NAME) must be specified"); + psFree(arguments); + return false; + } + if (config->camera && loadImage) { + psError(PS_ERR_IO, true, "Only one of (-camera NAME) and (-file NAME) may be specified"); + psFree(arguments); + return false; + } + + if (!psArgumentParse(arguments, &argc, argv)) { + psError(PS_ERR_IO, false, "error in command-line arguments"); + psFree(arguments); + return false; + } + + if (argc != 2) { + psError(PS_ERR_IO, true, "Missing output filename"); + psFree(arguments); + return false; + } + + // output filename + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + // save the 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) { + psError(PS_ERR_IO, false, "Unable to find recipe options for %s", PPSIM_RECIPE); + psFree(arguments); + return false; + } + + // these arguments may be used whether the input image is created or loaded + ppSimArgToRecipeF32(&status, options, "STARS.LUM", arguments, "-starslum"); + ppSimArgToRecipeF32(&status, options, "STARS.MAG", arguments, "-starsmag"); + ppSimArgToRecipeF32(&status, options, "STARS.DENSITY", arguments, "-starsdensity"); + ppSimArgToRecipeBool(&status, options, "PHOTOM", arguments, "+photom"); + + // if we are loading the input image (not creating it), then we can skip the remaining arguments + if (loadImage) { + // if we are supplying an input image, it is so we may supply fake stars; force it to + // be treated as an OBJECT image. + psMetadataAddStr(options, PS_LIST_TAIL, "IMAGE.TYPE", 0, "Exposure type", "OBJECT"); + + // check for these options as well + ppSimArgToRecipeStr(&status, options, "PSF.MODEL", arguments, "-psfclass"); // PSF model class + ppSimArgToRecipeStr(&status, options, "GALAXY.MODEL", arguments, "-galmodel"); // Galaxy model name + + // 'seeing' is not required: we can load a psf-model instead; but if not, it is allowed + ppSimArgToRecipeF32(&status, options, "SEEING", arguments, "-seeing"); // seeing (FWHM in arcsec) + + // 'scale' is not required: we can use the WCS instead; but if not, it is allowed + ppSimArgToRecipeF32(&status, options, "PIXEL.SCALE", arguments, "-scale"); // Plate scale + + // XXX we need to be more flexible: get this from the input image header or WCS + float ra0 = psMetadataLookupF32(NULL, arguments, "-ra"); // Right Ascension of boresight + float dec0 = psMetadataLookupF32(NULL, arguments, "-dec"); // Declination of boresight + 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); + + psFree (arguments); + return true; + } + + // 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) { + psError(PS_ERR_IO, false, "Unable to find FORMATS in camera configuration."); + psFree(arguments); + return false; + } + psMetadata *format = psMetadataLookupMetadata(NULL, formats, formatName); // Format of interest + if (!format) { + psError(PS_ERR_IO, false, "Unable to find format %s in camera FORMATS.", formatName); + psFree(arguments); + return false; + } + config->format = psMemIncrRefCounter(format); + } + + char *typeStr = ppSimArgToRecipeStr (&status, options, "IMAGE.TYPE", arguments, "-type"); // Requested exposure type + if (typeStr == NULL) { + psError(PS_ERR_IO, false, "An exposure type must be specified using -type"); + psFree(arguments); + return false; + } + ppSimType type = ppSimTypeFromString (typeStr); + + ppSimArgToRecipeStr(&status, options, "FILTER", arguments, "-filter"); // Filter name + + // set the exposure time + if (type == PPSIM_TYPE_BIAS) { + psMetadataAddF32(options, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", 0.0); + } else { + ppSimArgToRecipeF32(&status, options, "EXPTIME", arguments, "-exptime"); + if (!status) { + psError(PS_ERR_IO, false, "The exposure time must be specified using -exptime"); + psFree(arguments); + return false; + } + } + + // these values all get moved from arguments to RECIPE:OPTIONS + ppSimArgToRecipeF32(&status, options, "BIAS.LEVEL", arguments, "-biaslevel"); + ppSimArgToRecipeF32(&status, options, "BIAS.RANGE", arguments, "-biasrange"); + ppSimArgToRecipeF32(&status, options, "DARK.RATE", arguments, "-darkrate"); + ppSimArgToRecipeF32(&status, options, "FLAT.SIGMA", arguments, "-flatsigma"); + ppSimArgToRecipeF32(&status, options, "FLAT.RATE", arguments, "-flatrate"); + ppSimArgToRecipeF32(&status, options, "SHUTTER.TIME", arguments, "-shuttertime"); + ppSimArgToRecipeF32(&status, options, "SKY.RATE", arguments, "-skyrate"); + ppSimArgToRecipeS32(&status, options, "BINNING", arguments, "-bin"); + + 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)) { + psError(PS_ERR_IO, false, "-ra, -dec, -pa, -scale, -zp, -seeing must be specified for OBJECT type"); + psFree(arguments); + return false; + } + + ppSimArgToRecipeF32(&status, options, "PIXEL.SCALE", arguments, "-scale"); // Plate scale (arcsec / pixel) + ppSimArgToRecipeF32(&status, options, "SKY.MAGS", arguments, "-skymags"); // Plate scale + ppSimArgToRecipeStr(&status, options, "PSF.MODEL", arguments, "-psfclass"); // PSF model class + ppSimArgToRecipeStr(&status, options, "GALAXY.MODEL", arguments, "-galmodel"); // Galaxy model name + + // the user supplies FWHM in arcsec; here we convert to Sigma in pixels + psMetadataAddF32(options, PS_LIST_TAIL, "SEEING", 0, "Seeing FWHM (arcsec)", seeing); + 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); + + ppSimArgToRecipeF32(&status, options, "ZEROPOINT", arguments, "-zp"); // Zero point + } + + psFree(arguments); + return true; +} + +/* the following elements come from the config->arguments: + + PSPHOT.PSF + INPUT + OUTPUT + + all other values should come from the recipe +*/ Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimBadPixels.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimBadPixels.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimBadPixels.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimBounds.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimBounds.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimBounds.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimCreate.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 22153) @@ -0,0 +1,232 @@ +# include "ppSim.h" + +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, "PPSIM.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; + } + } else { + // If an image is supplied, we still generate a fake image and merge them together downstream + // (otherwise, we get the variance wrong). + simImage = false; + if (input->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.INPUT is not of type IMAGE"); + return NULL; + } + } + + // generate the fpa structure used by the output camera (determined from INPUT or specified) + assert (config->camera); + fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the observation + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + return NULL; + } + + // define the output image file -- this is the basis for the ppSimLoop + pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSIM.OUTPUT"); + if (!output) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from PPSIM.OUTPUT. Did you forget to specify the format?"); + return NULL; + } + if (output->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "PPSIM.OUTPUT type is not IMAGE"); + psFree(fpa); + return NULL; + } + // XXX we should not require the output image to be written + output->save = true; + + config->format = psMemIncrRefCounter (output->format); + config->formatName = psStringCopy (output->formatName); + + // the recipe is now fully realized for the desired camera + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate + + if (type == PPSIM_TYPE_OBJECT) { + // adjust the seeing by the scale + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); + if (isnan(seeing)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "seeing is not defined"); + psFree(fpa); + return NULL; + } + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); + if (isnan(scale)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "pixel scale is not defined"); + psFree(fpa); + return NULL; + } + psMetadataAddF32(recipe, PS_LIST_TAIL, "SEEING", PS_META_REPLACE, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); + + // if we have been supplied an input image, but no ra & dec, use the input image values + if (input) { + float ra = psMetadataLookupF32(&status, recipe, "RA"); + if (isnan(ra)) { + ra = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.RA"); + psMetadataAddF32(recipe, PS_LIST_TAIL, "RA", PS_META_REPLACE, "ra (radians)", ra); + } + float dec = psMetadataLookupF32(&status, recipe, "DEC"); + if (isnan(dec)) { + dec = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.DEC"); + psMetadataAddF32(recipe, PS_LIST_TAIL, "DEC", PS_META_REPLACE, "dec (radians)", dec); + } + } + } + + if ((type == PPSIM_TYPE_OBJECT) || (type == PPSIM_TYPE_FLAT)) { + // determine the zeropoint from the filter + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); + if (isnan(zp)) { + char *filter = psMetadataLookupStr(&status, recipe, "FILTER"); + float zp = ppSimGetZeroPoint (recipe, filter); + psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + } + } + + // For photometry, we operate on the chip-mosaicked image. we create a copy of the mosaicked + // image for psphot so we can write out a clean image + bool doPhotom = psMetadataLookupBool(&status, recipe, "PHOTOM"); // Density of fakes + if (doPhotom) { + + // XXX at the moment, we can perform photometry on the fake sources and the forced + // photometry positions, but not one or the other. Also, we only support photometry in + // the context of an image previously analysed by psphot. Add support for: + // * photometry of a pure fake image (requires peak detection and psf creation) + // * photometry of forced source positions without fake image (this might work, it just + // requires not generating any fake features). + + if (!input) { + psError(PS_ERR_UNKNOWN, false, "input image not found; currently required for photometry"); + return NULL; + } + + // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) + pmFPAfile *fakeImage = pmFPAfileDefineChipMosaic(config, output->fpa, "PPSIM.FAKE.CHIP"); + if (!fakeImage) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FAKE.CHIP")); + psFree(fpa); + return NULL; + } + if (fakeImage->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.FAKE.CHIP is not of type IMAGE"); + psFree(fpa); + return NULL; + } + + // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?) + pmFPAfile *forceImage = NULL; + if (input) { + forceImage = pmFPAfileDefineChipMosaic(config, input->fpa, "PPSIM.FORCE.CHIP"); + if (!forceImage) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FORCE.CHIP")); + psFree(fpa); + return NULL; + } + if (forceImage->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.FORCE.CHIP is not of type IMAGE"); + psFree(fpa); + return NULL; + } + } + + // define associated psphot input/output files + if (!ppSimPhotomFiles (config, fakeImage, forceImage)) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files for psphot"); + psFree(fpa); + return NULL; + } + } else { + // have we supplied a psf model? this happens in ppSimPhotomFiles if we request a photometry + // analysis. however, even if we do not, a psf model may be used to generate the fake + // sources. + if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { + // tie the psf file to the chipMosaic + pmFPAfileBindFromArgs(&status, output, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD"); + psFree(fpa); + return NULL; + } + } + } + + // PPSIM.SOURCES carries the constructed, fake sources with their true parameters + // XXX only invoke this code for OBJECT types of images? + pmFPAfile *simSources = pmFPAfileDefineOutput (config, output->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 derive certain values from the image, if possible + if (input) { + // we need to extract certain metadata from the image and populate the recipe. + // or else we need to set the fpa concepts based on the recipe options... + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + ppSimArgToRecipeF32(&status, recipe, "EXPTIME", input->fpa->concepts, "FPA.EXPOSURE"); + char *filter = ppSimArgToRecipeStr(&status, recipe, "FILTER", input->fpa->concepts, "FPA.FILTER"); + + float zp = ppSimGetZeroPoint (recipe, filter); + psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + } + + pmFPALevel phuLevel = pmFPAPHULevel(output->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, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } + } + + pmChip *chip; // Chip from FPA + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + if (phuLevel == PM_FPA_LEVEL_CHIP) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } + } + + pmCell *cell; // Cell from chip + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + if (phuLevel == PM_FPA_LEVEL_CELL) { + if (!pmFPAAddSourceFromView(fpa, "Simulation", view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } + } + } + } + + psFree(fpa); + psFree(view); + + return output; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimDetectionLimits.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimDetectionLimits.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimDetectionLimits.c (revision 22153) @@ -0,0 +1,47 @@ +# include "ppSim.h" + +// compare injected sources (PPSIM.SOURCES) and measured fake sources (PPSIM.FAKE.SOURCES) +// to determine detection limits and recovered magnitude errors +bool ppSimDetectionLimits (pmConfig *config, pmFPAview *view) { + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + // XXX do we need to ask if the analysis was performed to 1st or 2nd detection limit? + float NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + psArray *injectedSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.SOURCES"); + psArray *measuredSources = psMetadataLookupPtr (NULL, readout->analysis, "PPSIM.FAKE.SOURCES"); + psAssert (injectedSources->n == measuredSources->n, "mis-match between injected and measured sources"); + + psVector *mag = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *dmag = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *SN = psVectorAlloc (injectedSources->n, PS_TYPE_F32); + psVector *detect = psVectorAlloc (injectedSources->n, PS_TYPE_BOOL); + + // first measure, for each source, if it was detected (SN > limit) + // and the magnitude offset (dM = M_inject - M_measure) + for (int i = 0; i < injectedSources->n; i++) { + pmSource *injectSource = injectedSources->data[i]; + pmSource *measureSource = measuredSources->data[i]; + + SN->data.F32[i] = measureSource->peak->SN; + dmag->data.F32[i] = injectSource->psfMag - measureSource->psfMag; + mag->data.F32[i] = injectSource->psfMag; + detect->data.Bool[i] = (measureSource->peak->SN >= NSIGMA_PEAK); + } + + // generate a histogram consisting of the instrumental magnitude bin, the + // XXX what is resolution? (user parameter?) + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimDetections.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimDetections.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimDetections.c (revision 22153) @@ -0,0 +1,30 @@ +# include "ppSim.h" + +bool ppSimDetections (psImage *significance, psMetadata *recipe, psArray *sources) { + psAssert (sources, "programming error: ppSimDetections passed NULL sources"); + + bool status; + + // for each source, measure the significance of the peak at the given coordinate + // where does this get stored? + + // XXX need to get the effective Area from the PSF sigma + float SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "SIGMA_SMOOTH"); + psAssert (status, "SIGMA_SMOOTH missing: call psphotSignificanceImage first"); + float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH); + + int row0 = significance->row0; + int col0 = significance->col0; + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + pmPeak *peak = source->peak; + psAssert (peak, "peak is not defined for the source"); + + peak->value = significance->data.F32[peak->y-row0][peak->x-col0]; + peak->SN = sqrt(peak->value*effArea); + + } + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimInsertGalaxies.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimInsertGalaxies.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimInsertGalaxies.c (revision 22153) @@ -0,0 +1,189 @@ +# 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; } + + 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, "PIXEL.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, recipe, "BINNING"); // Binning in x and y + + // determine the galaxy model + char *modelName = psMetadataLookupStr(&mdok, recipe, "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_branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 22153) @@ -0,0 +1,140 @@ +# 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; } + + 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, "PIXEL.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, recipe, "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; + + // insert the source flux in the image + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY); + + // 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); + + // add the sources to the source array + psArrayAdd (sources, 100,source); + psFree(source); // Drop reference + } + + // NOTE: the pmFPAfile "PPSIM.OUTPUT" points at these sources + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "fake sources", sources); + psFree(sources); + + // XXX many leaks in here, i think + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c (revision 22153) @@ -0,0 +1,113 @@ +# include "ppSim.h" + +// Reset a pointer: free and set to NULL +#define RESET(PTR) \ + psFree(PTR); \ + PTR = NULL; + +psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view) { + + bool status; + + pmFPAfile *outFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); + psAssert(outFile, "missing PPSIM.OUTPUT"); + psAssert(outFile->fpa, "missing fpr for PPSIM.OUTPUT"); + + pmChip *outChip = pmFPAviewThisChip (view, outFile->fpa); + pmPSF *psf = psMetadataLookupPtr (&status, outChip->analysis, "PSPHOT.PSF"); + assert (psf); + + psArray *spots = psMetadataLookupPtr(&status, outFile->fpa->analysis, "FORCED.SPOTS"); + PS_ASSERT_PTR_NON_NULL (spots, NULL); + + // in this program, we are looping over the output image, rather than the input as in ppImage + pmFPAfile *srcFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.REAL.SOURCES"); + psAssert(srcFile, "missing PPSIM.REAL.SOURCES"); + + // select the fpa and chip from PPSIM.REAL.SOURCES, which carries the astrometry in the headers + pmFPA *fpa = srcFile->fpa; + pmChip *chip = pmFPAviewThisChip (view, fpa); + + // we use the readout for the PPSIM.FORCE.CHIP, which contains the chip-mosaicked image(s) + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.CHIP"); + + // XXX we should have only one cell and readout per chip, right? + + // read WCS from existing output from psphot & psastro + // XXX can we be more flexible? ie, use the header of PPSIM.FORCE.CHIP? + bool bilevelAstrometry = false; + status = psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); + psAssert (status, "failed to set astrometry"); + + // the pixel scale is only used to set the focal-plane coordinate system + // setting this value to 1 makes this system have the units of pixels + float pixelScale = 1.0; + status = psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry); + psAssert (status, "failed to set astrometry"); + + float minX = 0.0; + float maxX = readout->image->numCols; + float minY = 0.0; + float maxY = readout->image->numRows; + + // the sources is a subset within range of this chip + psArray *sources = psArrayAllocEmpty (100); + + // Select the spots within range of this readout. Project the spots to this chip + for (int i = 0; i < spots->n; i++) { + pmAstromObj *obj = spots->data[i]; + + psProject (obj->TP, obj->sky, fpa->toSky); + psPlaneTransformApply (obj->FP, fpa->fromTPA, obj->TP); + psPlaneTransformApply (obj->chip, chip->fromFPA, obj->FP); + + // limit the X,Y range of the objs to the selected chip + if (obj->chip->x < minX) continue; + if (obj->chip->x > maxX) continue; + if (obj->chip->y < minY) continue; + if (obj->chip->y > maxY) continue; + + // convert the pmAstromObj to pmSource + + // instantiate a model for the PSF at this location, Io = 1.0 + pmModel *model = pmModelFromPSFforXY (psf, obj->chip->x, obj->chip->y, 1.0); + + // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) + // XXX use a fixed radius?? + // float radius = model->modelRadius (model->params, roughNoise); + // radius = PS_MAX (radius, 1.0); + float radius = 5.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; + + // these are not really needed since we will be fitting for them + source->psfMag = NAN; + source->errMag = NAN; + + // insert the source flux in the image + // XXX not sure the offset is really 0,0 + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, 0.0, 0.0); + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, 0.0, 0.0); + + // 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); + + psArrayAdd (sources, 100, source); + psFree(source); // Drop local reference + } + psTrace ("psastro", 4, "Added %ld sources\n", sources->n); + + return sources; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadSpots.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadSpots.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadSpots.c (revision 22153) @@ -0,0 +1,71 @@ +# include "ppSim.h" + +// Load the relevant forced-photometry positions for this field. This function does not determine +// the pixel coordinates, merely the celestial coordinates for sources in the general vicinity. The +// sources are saved on the fpa->analysis metadata as FORCED.SPOTS. We always create an entry; it +// will be empty if no sources were selected or forced photometry is not desired. +bool ppSimLoadSpots (pmFPA *fpa, pmConfig *config) { + + bool status; + + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + bool forcedPhot = psMetadataLookupBool(&status, recipe, "FORCED.PHOT"); // Density of fakes + if (!forcedPhot) { + psArray *spots = psArrayAllocEmpty (1); + psMetadataAddArray (fpa->analysis, PS_LIST_TAIL, "FORCED.SPOTS", PS_META_REPLACE, "forced photometry positions", spots); + psFree (spots); // free the extra (local) reference; a copy remains on fpa->analysis + return true; + } + + // We read the catalogue stars using psastroLoadRefstars + psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE); + if (!astroRecipe) { + psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE); + return false; + } + + float ra0 = psMetadataLookupF32(NULL, recipe, "RA"); // Boresight RA (radians) + float dec0 = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians) + float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + + if (isnan(ra0) || isnan(dec0)) { + psError(PS_ERR_UNKNOWN, false, "image boresite coords not defined."); + return false; + } + + char *catdir = psMetadataLookupStr(NULL, recipe, "FORCED.CATDIR"); + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale; + psFree(bounds); + + // modify the PSASTRO recipe to use the values desired for FORCED.PHOT. we can use a view on + // the PSASTRO recipe because we are not calling psastro elsewhere in this program. XXX need to + // use the WCS to define the overlap region + psMetadataAddStr(astroRecipe, PS_LIST_TAIL, "DVO.CATDIR", PS_META_REPLACE, "", catdir); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", ra0 - radius); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "RA_MAX", PS_META_REPLACE, "", ra0 + radius); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", dec0 - radius); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", dec0 + radius); + + // tell psastroLoadRefstars to ignore photcode and maglim + psMetadataAddStr(astroRecipe, PS_LIST_TAIL, "DVO.GETSTAR.PHOTCODE", PS_META_REPLACE, "", "NONE"); + psMetadataAddF32(astroRecipe, PS_LIST_TAIL, "DVO.GETSTAR.MAG.MAX", PS_META_REPLACE, "", NAN); + + // load refstars from the catalog + psArray *spots = psastroLoadRefstars(config); + if (!spots || spots->n == 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars."); + psFree(spots); + return false; + } + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", spots->n); + + psMetadataAddArray (fpa->analysis, PS_LIST_TAIL, "FORCED.SPOTS", PS_META_REPLACE, "forced photometry positions", spots); + psFree (spots); // free the extra (local) reference; a copy remains on fpa->analysis + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadStars.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadStars.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoadStars.c (revision 22153) @@ -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; + } + + // relevant metadata + 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, "PIXEL.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) { + 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_branches/eam_branch_20080511/ppSim/src/ppSimLoop.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 22153) @@ -0,0 +1,220 @@ +#include "ppSim.h" + +# define ESCAPE(CODE,MSG) { \ + psError(CODE, false, MSG); \ + psFree (rng); \ + return false; } + +bool ppSimLoop(pmConfig *config) +{ + bool status; + + // XXX if we are supplying a PSF, then we should use that to specify the seeing. + // we will need to force the psf to be loaded here (deactivate everyone, activate psf, load + // it, then calculate seeing as appropriate). + + PS_ASSERT_PTR_NON_NULL(config, PS_EXIT_PROG_ERROR); + + // in this program, we are looping over the output image, rather than the input as in ppImage + pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); // Output file + assert(file); + + pmFPA *fpa = file->fpa; // FPA for file + assert(fpa); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate + int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y + + psArray *stars = psArrayAllocEmpty (1); + psArray *galaxies = psArrayAllocEmpty (1); + if (type == PPSIM_TYPE_OBJECT) { + // Load forced-photometry positions (these are placed on fpa->analysis for use in ppSimPhotomReadout) + if (!ppSimLoadSpots (fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load forced-photometry spots"); + + // Load catalogue stars + if (!ppSimLoadStars (stars, fpa, config)) ESCAPE (PS_ERR_UNKNOWN, "failed to load catalog stars"); + + // Add random stars + if (!ppSimMakeStars (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random stars"); + + // Add random galaxies + if (!ppSimMakeGalaxies (galaxies, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random galaxies"); + } + + 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); + if (fpa->hdu) { // XXX only do this if there is no INPUT image + if (!ppSimInitHeader(config, fpa, NULL, NULL)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); + } + + 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); + + pmReadout *readout; + while ((readout = pmFPAviewNextReadout (view, fpa, 1))) { + + // 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 (&status, readout, config, rng); + if (!status) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure"); + if (type == PPSIM_TYPE_BIAS) goto done; + + if (!ppSimMakeDark (readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating dark structure"); + if (type == PPSIM_TYPE_DARK) goto done; + + if (!ppSimMakeSky (readout, expCorr, type, config)) ESCAPE (PS_ERR_UNKNOWN, "problem generating sky background"); + if (type == PPSIM_TYPE_FLAT) goto done; + + if (type == PPSIM_TYPE_OBJECT) { + if (!ppSimInsertStars (readout, expCorr, stars, config)) ESCAPE (PS_ERR_UNKNOWN, "problem inserting stars"); + } + + if (type == PPSIM_TYPE_OBJECT) { + if (!ppSimInsertGalaxies (readout, expCorr, galaxies, config)) ESCAPE (PS_ERR_UNKNOWN, "problem inserting galaxies"); + } + + psFree(expCorr); + + done: + if (!ppSimAddNoise(readout->image, readout->weight, cell, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding noise"); + if (!ppSimSaturate(readout, config)) ESCAPE (PS_ERR_UNKNOWN, "problem setting saturation levels"); + + if (!ppSimBadPixels(readout, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding bad pixels"); + + if (!ppSimAddOverscan (readout, config, biasCols, biasRows, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding overscan region"); + psFree(biasRows); + + readout->data_exists = true; + readout->parent->data_exists = true; + readout->parent->parent->data_exists = true; + + // if there is an input image, merge it with the simulated image + if (!ppSimMergeReadouts (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image"); + } + psFree(biasCols); + + if (!ppSimUpdateConceptsCell (cell, config)) ESCAPE (PS_ERR_UNKNOWN, "problem updating cell concepts"); + + if (cell->hdu) { + // XXX only do this if there is no INPUT image? + if (!ppSimInitHeader(config, NULL, NULL, cell)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); + } + + 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; + } + } + + // XXX why no UpdateConceptsChip?? + + if (chip->hdu) { + // XXX only do this if there is no INPUT image + if (!ppSimInitHeader(config, NULL, chip, NULL)) ESCAPE (PS_ERR_UNKNOWN, "problem setting output header"); + } + + // we perform photometry on the readouts of this chip in the output + if (!ppSimPhotom (config, view)) ESCAPE (PS_ERR_UNKNOWN, "problem performing photometry"); + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_IO, false, "Unable to write file."); + psFree(rng); + psFree(view); + return false; + } + } + + psFree(stars); + psFree(galaxies); + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_IO, false, "Unable to write file."); + psFree(rng); + psFree(view); + return false; + } + + psFree(rng); + psFree(view); + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeBias.c (revision 22153) @@ -0,0 +1,63 @@ +# include "ppSim.h" + +psVector *ppSimMakeBias (bool *status, pmReadout *readout, pmConfig *config, const psRandom *rng) { + + bool mdok; + + if (status) *status = true; + + pmCell *cell = readout->parent; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + bool bias = psMetadataLookupBool(&mdok, recipe, "BIAS"); // Generate a Bias? + if (!bias) return NULL; + + 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."); + *status = false; + 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_branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c (revision 22153) @@ -0,0 +1,42 @@ +# include "ppSim.h" + +psVector *ppSimMakeBiassec (pmCell *cell, pmConfig *config) { + + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + bool bias = psMetadataLookupBool(&mdok, recipe, "BIAS"); // Generate a Bias? + if (!bias) return NULL; + + 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_branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeDark.c (revision 22153) @@ -0,0 +1,41 @@ +# 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 + + bool dark = psMetadataLookupBool(&mdok, recipe, "DARK"); // Generate a DARK? + if (!dark) return true; + + float darkRate = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate + if (isnan(darkRate)) { + psError(PS_ERR_UNKNOWN, false, "missing DARK.RATE\n"); + return false; + } + + float expTime = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time + if (isnan(expTime)) { + psError(PS_ERR_UNKNOWN, false, "missing EXPTIME\n"); + return false; + } + + psTrace("ppSim", 6, "darkRate: %f, expTime: %f\n", darkRate, expTime); + + 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_branches/eam_branch_20080511/ppSim/src/ppSimMakeGalaxies.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeGalaxies.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeGalaxies.c (revision 22153) @@ -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 true; + + 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, "PIXEL.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 true; + + // 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 true; + } + + // 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 true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 22153) @@ -0,0 +1,106 @@ +# include "ppSim.h" + +// this function sets the skyRate to a value for the night sky (SKY.RATE or SKY.MAGS) or for a +// flat-field image (FLAT.RATE). Include a shutter correction and a scattered light source + +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 + + bool sky = psMetadataLookupBool(&status, recipe, "SKY"); // Generate a SKY flux? + bool flat = psMetadataLookupBool(&status, recipe, "FLAT"); // Apply flat-field term? + + 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 scatterFrac = psMetadataLookupF32(&status, recipe, "SCATTER.FRAC"); // scattered light fraction (max) + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + if (!isnan(skyMags)) { + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + if (type == PPSIM_TYPE_FLAT) { + skyRate = flatRate; + } + + 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, recipe, "BINNING"); // Binning in x and y + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + int dXfpa = bounds->x1 - bounds->x0; + int dYfpa = bounds->y1 - bounds->y0; + + // Correct chip offsets so that boresight is in the middle of the FPA + x0Chip -= 0.5 * dXfpa; + y0Chip -= 0.5 * dYfpa; + + 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 with flatValue = 1.0 at the field center + float flatValue = 1.0; + if (flat) { + // we make the flat-field have a response of 1.0 at the field center (like a vignetting) + flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))); + } + + float scatterRate = 0.0; + + if (sky) { + // add a scattered light term to the flat-field images (no + if (type == PPSIM_TYPE_FLAT) { + float xF = 2.0*(xFPA / dXfpa) - 1.0; + scatterRate = scatterFrac * PS_SQR(xF); + } + + // Sky background + float skyFlux = (skyRate * (flatValue + scatterRate)) * realExpTime; // Flux from sky + signal->data.F32[y][x] += skyFlux; + variance->data.F32[y][x] += skyFlux; + } + + // used later to modify the star and galaxy photometry + if (expCorr) { + expCorr->data.F32[y][x] = realExpTime / expTime; + } + + // TO DO: Add fringes + + } + } + psFree(bounds); + return true; +} + Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeStars.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeStars.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMakeStars.c (revision 22153) @@ -0,0 +1,147 @@ +# include "ppSim.h" + +# define ESCAPE(MSG) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, MSG); \ + return false; } + +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, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + + if (isnan(darkRate)) darkRate = 0.0; + if (isnan(expTime)) ESCAPE("EXPTIME is not defined"); + if (isnan(zp)) ESCAPE("ZEROPOINT is not defined"); + if (isnan(seeing)) ESCAPE("SEEING is not defined"); + if (isnan(scale)) ESCAPE("PIXEL.SCALE is not defined"); + + 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 + + // if we tried and failed to load reference stars, set more artificial limits + if (!status) { + refMag = brightMag; + refSum = 1; + } + } 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_branches/eam_branch_20080511/ppSim/src/ppSimMergeReadouts.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMergeReadouts.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMergeReadouts.c (revision 22153) @@ -0,0 +1,48 @@ +# include "ppSim.h" + +// XXX add bounds to the inputs? +bool ppSimMergeReadouts (pmConfig *config, pmFPAview *view) { + + // bool mdok; + + // if we have an input image, we need to add the synthetic data on top of it below + // XXX we may potentially use this input file to skip missing elements + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSIM.INPUT"); + if (!input) return true; + + pmReadout *inReadout = pmFPAviewThisReadout (view, input->fpa); + if (!inReadout) return true; + + if (!inReadout->weight) { + if (!pmReadoutGenerateWeight(inReadout, true)) { + psError (PS_ERR_UNKNOWN, false, "trouble creating weight"); + return false; + } + } + + // output must exist or we made a programming error + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT"); // Output file + assert(output); + + // XXX require outReadout? + pmReadout *outReadout = pmFPAviewThisReadout (view, output->fpa); + if (!outReadout) return true; + + psImage *inSignal = inReadout->image; + psImage *inVariance = inReadout->weight; + + psImage *outSignal = outReadout->image; + psImage *outVariance = outReadout->weight; + + assert (inSignal->numRows == outSignal->numRows); + assert (inSignal->numCols == outSignal->numCols); + + for (int y = 0; y < inSignal->numRows; y++) { + for (int x = 0; x < inSignal->numCols; x++) { + outSignal->data.F32[y][x] += inSignal->data.F32[y][x]; + outVariance->data.F32[y][x] += inVariance->data.F32[y][x]; + } + } + return true; +} + Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMergeSources.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMergeSources.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMergeSources.c (revision 22153) @@ -0,0 +1,17 @@ +# include "ppSim.h" + +psArray *ppSimMergeSources (psArray *in1, psArray *in2) { + + psArray *out = psArrayAlloc (in1->n + in2->n); + + for (int i = 0; i < in1->n; i++) { + out->data[i] = psMemIncrRefCounter (in1->data[i]); + } + + int nOff = in1->n; + for (int i = nOff; i < out->n; i++) { + out->data[i] = psMemIncrRefCounter (in2->data[i-nOff]); + } + + return out; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMosaicChip.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMosaicChip.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimMosaicChip.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotom.c (revision 22153) @@ -0,0 +1,62 @@ +# 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; + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + bool doPhotom = psMetadataLookupBool(&status, recipe, "PHOTOM"); // Density of fakes + if (!doPhotom) return true; + + psphotModelClassInit (); + + int blankMask = 0; // XXX not sure what this should be set to... + ppSimMosaicChip(config, blankMask, view, "PPSIM.FORCE.CHIP", "PPSIM.INPUT"); + ppSimMosaicChip(config, blankMask, view, "PPSIM.FAKE.CHIP", "PPSIM.OUTPUT"); + + // use PPSIM.FAKE.CHIP as the base since it is guaranteed to exist (derived from PPSIM.OUTPUT) + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PPSIM.FAKE.CHIP"); + if (!status) { + psError(PSPHOT_ERR_CONFIG, false, "PPSIM.FAKE.CHIP I/O file is not defined"); + return false; + } + + // XXX If we want to be able to write out both the positive and residual images, we will + // need to define an additional pmFPAfile container, and then make a copy as is done for + // PPIMAGE.CHIP -> PSPHOT.INPUT in ppImage. At the moment, PPSIM.CHIP will write out the + // mosaicked image with the sources subtracted. + + // XXX if input -> PSPHOT.INPUT, which is tied to PPSIM.CHIP with pmFPAfileDefineFromFile, + // then the following code creates the chip copy: + // 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 (!ppSimPhotomReadoutFake (config, view)) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + return false; + } + if (!ppSimPhotomReadoutForce (config, view)) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + return false; + } + // ppSimDetectionLimits (config, view); + } + } + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomFiles.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomFiles.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomFiles.c (revision 22153) @@ -0,0 +1,83 @@ +# include "ppSim.h" + +// define the needed / desired I/O files +bool ppSimPhotomFiles (pmConfig *config, pmFPAfile *fakeFile, pmFPAfile *forceFile) { + + bool status = false; + + // select recipe options supplied on command line + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + + // optionally save the background model (small FITS image) + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) { + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, DX, DY, "PSPHOT.BACKMDL"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKMDL"); + return false; + } + output->save = true; + } + // optionally save the background model's standard deviation (small FITS image) + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL.STDEV")) { + int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN"); + int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN"); + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, DX, DY, "PSPHOT.BACKMDL.STDEV"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKMDL.STDEV"); + return false; + } + output->save = true; + } + // optionally save the full background image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKGND")) { + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, 1, 1, "PSPHOT.BACKGND"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKGND"); + return false; + } + output->save = true; + } + // optionally save the background-subtracted image + if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKSUB")) { + pmFPAfile *output = pmFPAfileDefineFromFile (config, fakeFile, 1, 1, "PSPHOT.BACKSUB"); + if (!output) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.BACKSUB"); + return false; + } + output->save = true; + } + + // the input sources are required for sensible measurements of the fake or force sources + pmFPAfile *cmffile = pmFPAfileDefineFromArgs (&status, config, "PPSIM.REAL.SOURCES", "INPUT.SOURCES"); + if (!cmffile) { + psError(PS_ERR_UNKNOWN, false, "input sources not found; required for photometry (or error with PPSIM.REAL.SOURCES filerule)"); + return false; + } + + // the fake sources are carried on the fakeFile->fpa structures (PPSIM.FAKE.CHIP) + pmFPAfile *outsources = pmFPAfileDefineOutput (config, fakeFile->fpa, "PPSIM.FAKE.SOURCES"); + if (!outsources) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PPSIM.FAKE.SOURCES"); + return false; + } + outsources->save = true; + + // the forced sources are carried on the forceFile->fpa structures (PPSIM.FORCE.CHIP) + pmFPAfile *outforced = pmFPAfileDefineOutput (config, forceFile->fpa, "PPSIM.FORCE.SOURCES"); + if (!outforced) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PPSIM.FORCE.SOURCES"); + return false; + } + outforced->save = true; + + // we bind the psffile to the chip-mosaicked image so we have a reference to the image + pmFPAfile *psffile = pmFPAfileBindFromArgs(&status, fakeFile, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + if (!psffile) { + psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.PSF.LOAD"); + return false; + } + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadout.c (revision 22153) @@ -0,0 +1,184 @@ +# include "ppSim.h" + +psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) { + + pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename); + PS_ASSERT_PTR_NON_NULL (readout, NULL); + + psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + return sources; +} + +bool ppSimDefinePixels (psArray *sources, pmReadout *readout, psMetadata *recipe) { + + bool status; + + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + if (!status) return NULL; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + } + return true; +} + +bool ppSimPhotomReadout(pmConfig *config, const pmFPAview *view) { + + psTimerStart ("psphotReadout"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + +# if 0 + // set the photcode for this image + if (!psphotAddPhotcode (recipe, config, view, "PPSIM.CHIP")) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } +# endif + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // the previously measured real sources are loaded into the PPSIM.REAL.SOURCES pmFPAfile + psArray *realSources = ppSimSelectSources (config, view, "PPSIM.REAL.SOURCES"); + PS_ASSERT_PTR_NON_NULL (realSources, false); + + // the newly injected fake sources are saved on the PPSIM.OUTPUT pmFPAfile + psArray *injectedSources = ppSimSelectSources (config, view, "PPSIM.OUTPUT"); + PS_ASSERT_PTR_NON_NULL (injectedSources, false); + + // XXX need to define the source pixels + ppSimDefinePixels (realSources, readout, recipe); + ppSimDefinePixels (injectedSources, readout, recipe); + + // make a copy of the fake-source parameters so these can be saved independently from the + // measured fake-source parameters + psArray *fakeSources = psArrayAlloc (injectedSources->n); + for (int i = 0; i < injectedSources->n; i++) { + fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]); + } + + // load the forced source lists + psArray *forceSources = ppSimLoadForceSources (config, view); + psAssert (forceSources, "failed to load force photometry sources"); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + // this function uses PSPHOT.PSF.LOAD as the pmFPAfile + pmPSF *psf = psphotLoadPSF (config, view, recipe); + assert (psf); + + // remove all sources + psphotRemoveAllSources (realSources, recipe); + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, "PPSIM.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PPSIM.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // add noise for real (subtracted) objects + psphotAddNoise (readout, realSources, recipe); + + // XXX fake sources should measure peak->x,y, force sources should not + psMaskType maskVal = 0xff; + psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); + ppSimDetections (significance, recipe, fakeSources); + ppSimDetections (significance, recipe, forceSources); + psFree (significance); + + // remove noise for subtracted objects (ie, return to normal noise level) + psphotSubNoise (readout, realSources, recipe); + + // replace all sources + psphotReplaceAllSources (realSources, recipe); + + // construct an initial model for each object + psphotGuessModels (readout, realSources, recipe, psf); + psphotGuessModels (readout, fakeSources, recipe, psf); + psphotGuessModels (readout, forceSources, recipe, psf); + + psArray *sources = NULL; + + // linear fit to real + fake sources + sources = ppSimMergeSources (realSources, fakeSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + + // linear fit to real + forced sources (if any were selected) + if (forceSources->n) { + sources = ppSimMergeSources (realSources, forceSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + } + + // XXX do we need to measure aperture photometry corrections? + // XXX do we store the pre and post correction magnitudes? + // XXX for the fake sources, these must be identically zero + // XXX for the force sources, need the apresid to put the mags on the correct system + if (0) { + if (!psphotApResid (readout, forceSources, recipe, psf)) { + psLogMsg ("psphot", 3, "failed on psphotApResid"); + return psphotReadoutCleanup (config, readout, recipe, NULL, psf, NULL); + } + } + + // calculate source magnitudes (for which set??) + pmReadout *background = psphotSelectBackground (config, view, false); + psphotMagnitudes(fakeSources, recipe, psf, background); + psphotMagnitudes(forceSources, recipe, psf, background); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (realSources); + psphotSourceFreePixels (fakeSources); + psphotSourceFreePixels (forceSources); + + // create the exported-metadata and free local data + // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? + // (or don't supply the sources, and do this with a different function) + psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); + + // add forceSources to the readout (real and fake already there) + // these are outputs: we need to generate the fpa structure + + pmCell *forceCell = pmFPAfileThisCell (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceCell, "no cell?"); + pmChip *forceChip = pmFPAfileThisChip (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceChip, "no chip?"); + pmReadout *forceReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.SOURCES"); + if (!forceReadout) { + forceReadout = pmReadoutAlloc (forceCell); + psFree (forceReadout); // there is a copy on 'cell' as well + } + psAssert (forceReadout, "no forceReadout?"); + pmChipSetDataStatus (forceChip, true); + psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "forced photometry ", forceSources); + psFree (forceSources); + + pmCell *fakeCell = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?"); + pmChip *fakeChip = pmFPAfileThisChip (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeChip, "no chip?"); + pmReadout *fakeReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.SOURCES"); + if (!fakeReadout) { + fakeReadout = pmReadoutAlloc (fakeCell); + psFree (fakeReadout); // there is a copy on 'cell' as well + } + psAssert (fakeReadout, "no fakeReadout?"); + pmChipSetDataStatus (fakeChip, true); + psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources); + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadoutFake.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadoutFake.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadoutFake.c (revision 22153) @@ -0,0 +1,124 @@ +# include "ppSim.h" + +bool ppSimPhotomReadoutFake(pmConfig *config, const pmFPAview *view) { + + psTimerStart ("psphotReadout"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + +# if 0 + // set the photcode for this image + if (!psphotAddPhotcode (recipe, config, view, "PPSIM.FAKE.CHIP")) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } +# endif + + // *** in this section, perform the photometry for real + fake sources on PPSIM.FAKE.CHIP *** + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // the previously measured real sources are loaded into the PPSIM.REAL.SOURCES pmFPAfile + psArray *realMeasuredSources = ppSimSelectSources (config, view, "PPSIM.REAL.SOURCES"); + PS_ASSERT_PTR_NON_NULL (realMeasuredSources, false); + + // the newly injected fake sources are saved on the PPSIM.OUTPUT pmFPAfile + psArray *injectedSources = ppSimSelectSources (config, view, "PPSIM.OUTPUT"); + PS_ASSERT_PTR_NON_NULL (injectedSources, false); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); + + // XXX need to define the source pixels + ppSimDefinePixels (realMeasuredSources, readout, recipe); + ppSimDefinePixels (injectedSources, readout, recipe); + + // make a copy of the injected fake-source parameters so these can be saved independently + // from the measured fake-source parameters + psArray *fakeSources = psArrayAlloc (injectedSources->n); + for (int i = 0; i < injectedSources->n; i++) { + fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]); + } + + // make a copy of the measured real-source parameters so these can be fitted here and in + // ppSimPhotomReadoutForce, if desired. + psArray *realSources = psArrayAlloc (realMeasuredSources->n); + for (int i = 0; i < realMeasuredSources->n; i++) { + realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]); + } + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + // this function uses PSPHOT.PSF.LOAD as the pmFPAfile + pmPSF *psf = psphotLoadPSF (config, view, recipe); + assert (psf); + + // remove all sources + psphotRemoveAllSources (realSources, recipe); + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, "PPSIM.FAKE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PPSIM.FAKE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // add noise for real (subtracted) objects + psphotAddNoise (readout, realSources, recipe); + + // XXX fake sources should measure peak->x,y, force sources should not + psMaskType maskVal = 0xff; + psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); + ppSimDetections (significance, recipe, fakeSources); + psFree (significance); + + // remove noise for subtracted objects (ie, return to normal noise level) + psphotSubNoise (readout, realSources, recipe); + + // replace all sources + psphotReplaceAllSources (realSources, recipe); + + // construct an initial model for each object + psphotGuessModels (readout, realSources, recipe, psf); + psphotGuessModels (readout, fakeSources, recipe, psf); + + // linear fit to real + fake sources + psArray *sources = ppSimMergeSources (realSources, fakeSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + + // calculate source magnitudes (for which set??) + pmReadout *background = psphotSelectBackground (config, view, false); + psphotMagnitudes(fakeSources, recipe, psf, background); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (realSources); + psphotSourceFreePixels (fakeSources); + + // create the exported-metadata and free local data + // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? + // (or don't supply the sources, and do this with a different function) + psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); + + pmCell *fakeCell = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?"); + pmChip *fakeChip = pmFPAfileThisChip (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeChip, "no chip?"); + pmReadout *fakeReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FAKE.SOURCES"); + if (!fakeReadout) { + fakeReadout = pmReadoutAlloc (fakeCell); + psFree (fakeReadout); // there is a copy on 'cell' as well + } + psAssert (fakeReadout, "no fakeReadout?"); + pmChipSetDataStatus (fakeChip, true); + psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources); + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadoutForce.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadoutForce.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimPhotomReadoutForce.c (revision 22153) @@ -0,0 +1,117 @@ +# include "ppSim.h" + +bool ppSimPhotomReadoutForce(pmConfig *config, const pmFPAview *view) { + + psTimerStart ("psphotReadout"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + +# if 0 + // set the photcode for this image + if (!psphotAddPhotcode (recipe, config, view, "PPSIM.FAKE.CHIP")) { + psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); + return false; + } +# endif + + // *** in this section, perform the photometry for real + force sources on PPSIM.FORCE.CHIP *** + + // find the currently selected readout. + // we always perform photometry on the mosaiced chip + pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.CHIP"); + PS_ASSERT_PTR_NON_NULL (readout, false); + + // the previously measured real sources are loaded into the PPSIM.REAL.SOURCES pmFPAfile + psArray *realMeasuredSources = ppSimSelectSources (config, view, "PPSIM.REAL.SOURCES"); + PS_ASSERT_PTR_NON_NULL (realMeasuredSources, false); + + // load the forced source lists + psArray *forceSources = ppSimLoadForceSources (config, view); + psAssert (forceSources, "failed to load force photometry sources"); + + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndWeight (config, readout, recipe); + + // XXX need to define the source pixels + ppSimDefinePixels (realMeasuredSources, readout, recipe); + ppSimDefinePixels (forceSources, readout, recipe); + + // make a copy of the measured real-source parameters so these can be fitted here and in + // ppSimPhotomReadoutForce, if desired. + psArray *realSources = psArrayAlloc (realMeasuredSources->n); + for (int i = 0; i < realMeasuredSources->n; i++) { + realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]); + } + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe + // this function uses PSPHOT.PSF.LOAD as the pmFPAfile + pmPSF *psf = psphotLoadPSF (config, view, recipe); + assert (psf); + + // remove all sources + psphotRemoveAllSources (realSources, recipe); + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, "PPSIM.FORCE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + if (!psphotSubtractBackground (config, view, "PPSIM.FORCE.CHIP")) { + return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + } + + // add noise for real (subtracted) objects + psphotAddNoise (readout, realSources, recipe); + + // XXX fake sources should measure peak->x,y, force sources should not + psMaskType maskVal = 0xff; + psImage *significance = psphotSignificanceImage (readout, recipe, 1, maskVal); + ppSimDetections (significance, recipe, forceSources); + psFree (significance); + + // remove noise for subtracted objects (ie, return to normal noise level) + psphotSubNoise (readout, realSources, recipe); + + // replace all sources + psphotReplaceAllSources (realSources, recipe); + + // construct an initial model for each object + psphotGuessModels (readout, realSources, recipe, psf); + psphotGuessModels (readout, forceSources, recipe, psf); + + // linear fit to real + force sources + psArray *sources = ppSimMergeSources (realSources, forceSources); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); // XXX option to NOT subtract + psphotReplaceAllSources (sources, recipe); + psFree (sources); // only frees the merged references + + // calculate source magnitudes (for which set??) + pmReadout *background = psphotSelectBackground (config, view, false); + psphotMagnitudes(forceSources, recipe, psf, background); + + // drop the references to the image pixels held by each source + psphotSourceFreePixels (realSources); + psphotSourceFreePixels (forceSources); + + // create the exported-metadata and free local data + // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? + // (or don't supply the sources, and do this with a different function) + psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); + + pmCell *forceCell = pmFPAfileThisCell (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceCell, "no cell?"); + pmChip *forceChip = pmFPAfileThisChip (config->files, view, "PPSIM.FORCE.SOURCES"); psAssert (forceChip, "no chip?"); + pmReadout *forceReadout = pmFPAfileThisReadout (config->files, view, "PPSIM.FORCE.SOURCES"); + if (!forceReadout) { + forceReadout = pmReadoutAlloc (forceCell); + psFree (forceReadout); // there is a copy on 'cell' as well + } + psAssert (forceReadout, "no forceReadout?"); + pmChipSetDataStatus (forceChip, true); + psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "force photometry ", forceSources); + + return true; +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimRandomGaussian.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimRandomGaussian.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimRandomGaussian.c (revision 22153) @@ -0,0 +1,102 @@ +# include "ppSim.h" + +static int Ngaussint = 0; +static double *gaussint = NULL; + +extern double drand48(); + +double p_ppSimGaussian (double x, double mean, double sigma) { + + double f; + + f = exp (-0.5 * PS_SQR(x - mean) / PS_SQR(sigma)) / sqrt(2 * M_PI * PS_SQR(sigma)); + + return (f); + +} + +void ppSimRandomGaussianFree() +{ + psFree (gaussint); + return; +} + +void ppSimRandomGaussianAlloc (int Nbin) { + + gaussint = (double *) psAlloc(Nbin*sizeof(double)); + return; +} + +/* integrate a gaussian from -5 sigma to +5 sigma */ +void p_ppSimRandomGaussianInit () { + + int i; + long A, B; + double val, x, dx, dx1, dx2, dx3, df; + double mean, sigma; + + /* no need to generate this if it already exists */ + if (gaussint) return; + + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + + Ngaussint = 0x1000; + ppSimRandomGaussianAlloc (Ngaussint + 1); + + val = 0; + dx = 1.0 / Ngaussint; + dx1 = dx / 3.0; + dx2 = 2.0*dx/3.0; + dx3 = dx; + mean = 0.0; + sigma = 1.0; + + for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx) { + df = (3.0*p_ppSimGaussian(x , mean, sigma) + + 9.0*p_ppSimGaussian(x+dx1, mean, sigma) + + 9.0*p_ppSimGaussian(x+dx2, mean, sigma) + + 3.0*p_ppSimGaussian(x+dx3, mean, sigma)) * (dx1/8.0); + val += df; + if (val > (i + 0.5) / (double) Ngaussint) { + gaussint[i] = x + dx / 2.0; + i++; + } + } +} + +// XXX we are using drand48() rather than the random var supplied by rnd +double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma) { + + int i; + double y; + + if (gaussint == NULL) { + p_ppSimRandomGaussianInit (); + } + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]*sigma + mean; + + return (y); + +} + +// XXX we are using drand48() rather than the random var supplied by rnd +double ppSimRandomGaussianNorm (const psRandom *rnd) { + + int i; + double y; + + if (gaussint == NULL) { + p_ppSimRandomGaussianInit (); + } + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]; + + return (y); +} Index: /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSaturate.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSaturate.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSaturate.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimSequence.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequence.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequence.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimSequence.h =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequence.h (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequence.h (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimSequenceBias.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceBias.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceBias.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimSequenceDark.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceDark.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceDark.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimSequenceFlat.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceFlat.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceFlat.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimSequenceObject.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceObject.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSequenceObject.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimSetPSF.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSetPSF.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimSetPSF.c (revision 22153) @@ -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, recipe, "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_branches/eam_branch_20080511/ppSim/src/ppSimStars.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimStars.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimStars.c (revision 22153) @@ -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_branches/eam_branch_20080511/ppSim/src/ppSimUtils.c =================================================================== --- /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimUtils.c (revision 22153) +++ /branches/eam_branches/eam_branch_20080511/ppSim/src/ppSimUtils.c (revision 22153) @@ -0,0 +1,301 @@ +# 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, "PIXEL.SCALE"); // plate scale in arcsec / pixel + scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel + + int binning = psMetadataLookupS32(NULL, recipe, "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; +} + +char *ppSimTypeToString (ppSimType type) { + + char *typeStr; + + switch (type) { + case PPSIM_TYPE_BIAS: typeStr = psStringCopy ("BIAS"); break; + case PPSIM_TYPE_DARK: typeStr = psStringCopy ("DARK"); break; + case PPSIM_TYPE_FLAT: typeStr = psStringCopy ("FLAT"); break; + case PPSIM_TYPE_OBJECT: typeStr = psStringCopy ("OBJECT"); break; + default: + psAbort("Should never get here."); + } + return (typeStr); +} + +ppSimType ppSimTypeFromString (char *typeStr) { + + if (!strcasecmp (typeStr, "BIAS")) return PPSIM_TYPE_BIAS; + if (!strcasecmp (typeStr, "DARK")) return PPSIM_TYPE_DARK; + if (!strcasecmp (typeStr, "FLAT")) return PPSIM_TYPE_FLAT; + if (!strcasecmp (typeStr, "OBJECT")) return PPSIM_TYPE_OBJECT; + psAbort("Should never get here."); +} + +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, recipe, "FILTER"); // Filter name + if (!filter) { + filter = "NONE"; + } + + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + + 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.EXPOSURE", 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, recipe, "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; +} + +bool ppSimRecipeValidation (pmConfig *config) { + + bool status; + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + int binning = psMetadataLookupS32(&status, recipe, "BINNING"); // 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); + } + return true; +} + +// Get a value from the command-line arguments and add it to recipe options +float ppSimArgToRecipeF32(bool *status, + 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 myStatus; + float value = psMetadataLookupF32(&myStatus, arguments, argName); // Value of interest + if (status) { *status = myStatus; } + if (isnan(value)) return value; + + psMetadataAddF32(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +// Get a value from the command-line arguments and add it to recipe options +int ppSimArgToRecipeS32(bool *status, + 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 myStatus; + int value = psMetadataLookupS32(&myStatus, arguments, argName); // Value of interest + if (status) { *status = myStatus; } + + psMetadataAddS32(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +// Get a value from the command-line arguments and add it to recipe options +bool ppSimArgToRecipeBool(bool *status, + 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 myStatus; + bool value = psMetadataLookupS32(&myStatus, arguments, argName); // Value of interest + if (status) { *status = myStatus; } + + psMetadataAddBool(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +// Get a value from the command-line arguments and add it to recipe options +char *ppSimArgToRecipeStr(bool *status, + 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 myStatus; + + char *value = psMetadataLookupStr(&myStatus, arguments, argName); // Value of interest + if (status) { + *status = myStatus; + } + psMetadataAddStr(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +float ppSimGetZeroPoint (psMetadata *recipe, char *filter) { + + bool mdok; + float 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 (&mdok, entry, "FILTER"); + assert (filterName); + + if (strcmp(filterName, filter)) continue; + + zp = psMetadataLookupF32 (&mdok, entry, "ZERO_PT"); + assert (mdok); + psFree (entries); + return zp; + } + psFree (entries); + return NAN; +} + +psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) { + + pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename); + PS_ASSERT_PTR_NON_NULL (readout, NULL); + + psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + return sources; +} + +bool ppSimDefinePixels (psArray *sources, pmReadout *readout, psMetadata *recipe) { + + bool status; + + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + if (!status) return NULL; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + } + return true; +} +