Index: /branches/ccl_branches/ipponly-20191108/ppSim/Makefile.am =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/Makefile.am (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/Makefile.am (revision 41095) @@ -0,0 +1,3 @@ +SUBDIRS = src + +CLEANFILES = *~ core core.* Index: /branches/ccl_branches/ipponly-20191108/ppSim/autogen.sh =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/autogen.sh (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/autogen.sh (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/configure.ac =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/configure.ac (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/configure.ac (revision 41095) @@ -0,0 +1,33 @@ +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_STDLDFLAGS + +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_VERSION + +IPP_STDOPTS +IPP_STDCFLAGS + +AC_CONFIG_FILES([ + Makefile + src/Makefile +]) +AC_OUTPUT Index: /branches/ccl_branches/ipponly-20191108/ppSim/notes.txt =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/notes.txt (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/notes.txt (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/psf.txt =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/psf.txt (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/psf.txt (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/src/Makefile.am =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/Makefile.am (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/Makefile.am (revision 41095) @@ -0,0 +1,62 @@ +bin_PROGRAMS = ppSim ppSimSequence + +# Force recompilation of ppSimVersion.c, since it gets the version information +ppSimVersion.c: ppSimVersionDefinitions.h +ppSimVersionDefinitions.h: ppSimVersionDefinitions.h.in FORCE + pslib-setsvnversion.pl PPSIM ppSimVersionDefinitions.h.in ppSimVersionDefinitions.h +FORCE: ; + +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 \ + ppSimMakeStarGrid.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 \ + ppSimSmoothReadout.c \ + ppSimSetPSF.c \ + ppSimUtils.c \ + ppSimLoop.c \ + ppSimMergeReadouts.c \ + ppSimMergeSources.c \ + ppSimMosaicChip.c \ + ppSimRandomGaussian.c \ + ppSimBadPixels.c \ + ppSimBadCTE.c \ + ppSimVersion.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/ccl_branches/ipponly-20191108/ppSim/src/ppSim.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSim.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSim.c (revision 41095) @@ -0,0 +1,42 @@ +# include "ppSim.h" + +int failure (pmConfig *config, psExit status, char *message) { + ppSimRandomGaussianFree (); + psErrorStackPrint(stderr, "%s", 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/ccl_branches/ipponly-20191108/ppSim/src/ppSim.h =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSim.h (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSim.h (revision 41095) @@ -0,0 +1,217 @@ +#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(void); +ppSimGalaxy *ppSimGalaxyAlloc(void); + +/// 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 ppSimMakeStarGrid(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng); +bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config); + +bool ppSimSmoothReadout(pmReadout *input, psMetadata *recipe); + +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 psImageMaskType blankMask, const pmFPAview *view, const char *outFile, const char *inFile); + +bool ppSimPhotom (pmConfig *config, pmFPAview *view); + + +// add a bad CTE region +bool ppSimBadCTE(psImage *image, // Signal image, modified and returned + const pmConfig *config // configuration + ); + +/// 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, const char *filter); + +bool ppSimMergeReadouts (pmConfig *config, pmFPAview *view); + +double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma); +double ppSimRandomGaussianNorm (const psRandom *rnd); +void ppSimRandomGaussianFree(void); + +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); + +/// Return software version +psString ppSimVersion(void); + +/// Return software source +psString ppSimSource(void); + +/// Return long version information +psString ppSimVersionLong(void); + +#endif Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimAddNoise.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimAddNoise.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimAddNoise.c (revision 41095) @@ -0,0 +1,43 @@ +# 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 + + // the recipe should set GAIN to NAN, and only modify to override the concept value + float gain = psMetadataLookupF32(&mdok, recipe, "GAIN"); // CCD gain, e/ADU + if (isnan(gain)) { + gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); + if (!mdok) { + psWarning("CELL.GAIN is not set; assuming gain of 1.0."); + gain = 1.0; + } + } + + // 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 the variance probably should be scaled as well -- we only are OK since GAIN = 1 + } + } + + // XXX why return this?? + return true; +} + Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimAddOverscan.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimAddOverscan.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimAddOverscan.c (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/src/ppSimArguments.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimArguments.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimArguments.c (revision 41095) @@ -0,0 +1,253 @@ +# 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, "-obs_mode", 0, "observation mode", 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, "-dettemp", 0, "Detector Temperature (C)", 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, "-seeing-max", 0, "Seeing FWHM (arcsec) at end of ramp", NAN); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-seeing-ramp", 0, "Use a seeing ramp", false); + 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); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-nx", 0, "cell x-size in pixels", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-ny", 0, "cell y-size in pixels", 0); + 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"); + psString obs_mode = psMetadataLookupStr(&status, arguments, "-obs_mode"); + if (obs_mode) { + psMetadataAddStr(options, PS_LIST_TAIL, "OBS_MODE", 0, "observation mode", obs_mode); + } + + + // 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) + ppSimArgToRecipeF32(&status, options, "SEEING.MAX", arguments, "-seeing-max"); // seeing (FWHM in arcsec) + ppSimArgToRecipeBool(&status, options, "SEEING.RAMP", arguments, "-seeing-ramp"); // 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); + // XXX handle error + + 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, "DET.TEMP", arguments, "-dettemp"); + 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"); + + ppSimArgToRecipeS32(&status, config->arguments, "NX.CELL", arguments, "-nx"); + ppSimArgToRecipeS32(&status, config->arguments, "NY.CELL", arguments, "-ny"); + + 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; we convert to Sigma in pixels (in ppSimCreate.c:103) + 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); + + if (psMetadataLookupBool(NULL, arguments, "-seeing-ramp")) { + float seeingMax = psMetadataLookupF32(NULL, arguments, "-seeing-max"); // seeing at ramp end + psMetadataAddF32(options, PS_LIST_TAIL, "SEEING.MAX", 0, "Seeing FWHM (arcsec)", seeingMax); + psMetadataAddBool(options, PS_LIST_TAIL, "SEEING.RAMP", 0, "use seeing ramp", true); + } + + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimBadCTE.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBadCTE.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBadCTE.c (revision 41095) @@ -0,0 +1,34 @@ +# include "ppSim.h" + +// Add a "bad CTE" region to the image : it is not really modeled as bad CTE, rather it is +// simply smoothed by a Gaussian kernel. This has the same effect on the image variance as bad +// CTE, but not (quite) the same effect on stellar photometry + +bool ppSimBadCTE(psImage *image, // Signal image, modified and returned + const pmConfig *config // configuration + ) +{ + assert(image->type.type == PS_TYPE_F32); + + bool status; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + if (!psMetadataLookupBool (&status, recipe, "BADCTE")) { + return true; + } + + char *region = psMetadataLookupStr (&status, recipe, "BADCTE.REGION"); + psRegion badCTE = psRegionForImage (image, psRegionFromString (region)); + if (psRegionIsNaN (badCTE)) psAbort("analysis region mis-defined"); + + float sigma = psMetadataLookupF32 (&status, recipe, "BADCTE.SIGMA"); + int Nsigma = psMetadataLookupS32 (&status, recipe, "BADCTE.NSIGMA"); + + psImage *subset = psImageSubset (image, badCTE); + if (!subset) psAbort ("error in badCTE region?"); + + psImageSmooth (subset, sigma, Nsigma); + + return true; +} + Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBadPixels.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBadPixels.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBadPixels.c (revision 41095) @@ -0,0 +1,107 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include + +#include "ppSim.h" + +// Mode for bad pixels +typedef enum { + BADPIX_ERROR, // Error: mode not set + BADPIX_RANDOM, // Random values for bad pixels + BADPIX_NAN, // NAN values for bad pixels +} badpixMode; + + +// 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; + } + + int size = psMetadataLookupS32(&mdok, recipe, "BADPIX.SIZE"); // (Half-)size of bad pixels + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "BADPIX.SIZE not set."); + return false; + } + + psString modeString = psMetadataLookupStr(&mdok, recipe, "BADPIX.MODE"); // Mode for bad pixels + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "BADPIX.MODE not set."); + return false; + } + badpixMode mode = BADPIX_ERROR; + if (strcmp(modeString, "RANDOM") == 0) { + mode = BADPIX_RANDOM; + } else if (strcmp(modeString, "NAN") == 0) { + mode = BADPIX_NAN; + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "BADPIX.MODE not recognised."); + return false; + } + + psRandom *pseudoRNG = psRandomAllocSpecific(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++) { + int vMin = PS_MAX(y - size, 0), vMax = PS_MIN(y + size, numRows - 1); // Extent of CR + for (int x = 0; x < numCols; x++) { + if (psRandomUniform(pseudoRNG) < frac) { + int uMin = PS_MAX(x - size, 0), uMax = PS_MIN(x + size, numCols - 1); // Extent of CR + for (int v = vMin; v <= vMax; v++) { + for (int u = uMin; u <= uMax; u++) { + switch (mode) { + case BADPIX_RANDOM: + image->data.F32[v][u] *= psRandomGaussian(rng); + break; + case BADPIX_NAN: + image->data.F32[v][u] = NAN; + break; + default: + psAbort("Unrecognised mode: %x", mode); + } + } + } + } + } + } + + psFree(pseudoRNG); + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBounds.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBounds.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimBounds.c (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/src/ppSimCreate.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimCreate.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimCreate.c (revision 41095) @@ -0,0 +1,236 @@ +# include "ppSim.h" + +pmFPAfile *ppSimCreate(pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(config, NULL); + + bool status; + 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). + if (input->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSIM.INPUT is not of type IMAGE"); + return NULL; + } + } + + assert (config->camera); + + int nx = psMetadataLookupS32 (&status, config->arguments, "NX.CELL"); + int ny = psMetadataLookupS32 (&status, config->arguments, "NY.CELL"); + + if (nx || ny) { + // need to find format(s) + psMetadata *formats = psMetadataLookupPtr(&status, config->camera, "FORMATS"); + psAssert (formats, "FORMATS missing from camera config"); + + psMetadataIterator *iter = psMetadataIteratorAlloc(formats, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *format; // Item from iteration + while ((format = psMetadataGetAndIncrement(iter))) { + psAssert (format->type == PS_DATA_METADATA, "unexpected format block"); + + if (nx) { + psMetadata *defaults = psMetadataLookupPtr(&status, format->data.md, "DEFAULTS"); + psMetadataLookupS32 (&status, defaults, "CELL.XSIZE"); + psAssert (status, "CELL.XSIZE should be in DEFAULTS"); + psMetadataAddF32(defaults, PS_LIST_TAIL, "CELL.XSIZE", PS_META_REPLACE, "", nx); + } + if (ny) { + psMetadata *defaults = psMetadataLookupPtr(&status, format->data.md, "DEFAULTS"); + psMetadataLookupS32 (&status, defaults, "CELL.YSIZE"); + psAssert (status, "CELL.YSIZE should be in DEFAULTS"); + psMetadataAddF32(defaults, PS_LIST_TAIL, "CELL.YSIZE", PS_META_REPLACE, "", ny); + } + } + psFree(iter); + } + + // generate the fpa structure used by the output camera (determined from INPUT or specified) + 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); + + float seeingMax = psMetadataLookupF32(&status, recipe, "SEEING.MAX"); + psMetadataAddF32(recipe, PS_LIST_TAIL, "SEEING.MAX", PS_META_REPLACE, "Seeing SIGMA (pixels)", seeingMax / 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) { + psError(PS_ERR_UNKNOWN, false, "in-line photometry in ppSim had been deprecated"); + return NULL; + } + + // 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.FILTERID"); + + float zp = ppSimGetZeroPoint(recipe, filter); + if (!isfinite(zp)) { + psError(PS_ERR_UNKNOWN, false, "Unable to find zero point for filter %s", filter); + psFree(fpa); + return NULL; + } + 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, 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, 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, view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA."); + psFree(fpa); + psFree(view); + return NULL; + } + } + // XXX this is a hack, but I don't have a better way at the moment: assumes a single cell per chip + if (nx) { + psMetadataAddS32(chip->concepts, PS_LIST_TAIL, "CHIP.XSIZE", PS_META_REPLACE, "", nx); + } else { + nx = psMetadataLookupF32(&status, cell->concepts, "CELL.XSIZE"); + psMetadataAddS32(chip->concepts, PS_LIST_TAIL, "CHIP.XSIZE", PS_META_REPLACE, "", nx); + } + if (ny) { + psMetadataAddS32(chip->concepts, PS_LIST_TAIL, "CHIP.YSIZE", PS_META_REPLACE, "", ny); + } else { + ny = psMetadataLookupF32(&status, cell->concepts, "CELL.YSIZE"); + psMetadataAddS32(chip->concepts, PS_LIST_TAIL, "CHIP.YSIZE", PS_META_REPLACE, "", ny); + } + } + } + + psFree(fpa); + psFree(view); + + return output; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimDetectionLimits.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimDetectionLimits.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimDetectionLimits.c (revision 41095) @@ -0,0 +1,48 @@ +# include "ppSim.h" + +// compare injected sources (PPSIM.SOURCES) and measured fake sources (PPSIM.FAKE.SOURCES) +// to determine detection limits and recovered magnitude errors +// XXX this function is not used +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/ccl_branches/ipponly-20191108/ppSim/src/ppSimDetections.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimDetections.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimDetections.c (revision 41095) @@ -0,0 +1,29 @@ +# include "ppSim.h" + +// XXX this function is not used +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->detValue = significance->data.F32[peak->y-row0][peak->x-col0]; + } + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimInsertGalaxies.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimInsertGalaxies.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimInsertGalaxies.c (revision 41095) @@ -0,0 +1,192 @@ +# include "ppSim.h" +static char *defaultModel = "PS_MODEL_SERSIC"; + +// this is used to measure the inserted flux as a difference of 2 large numbers : must use +// 'double' or we get floating point errors and trends on a scale of ~0.005 mags at > -10 mags +// instrumental +double imageSum (psImage *image) { + double sum = 0.0; + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + sum += image->data.F32[iy][ix]; + } + } + return sum; +} + +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"); + psAssert (psf, "missing PSPHOT.PSF on chip.analysis"); + + int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning); + int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning); + + pmDetections *detections = psMetadataLookupPtr (&mdok, readout->analysis, "PSPHOT.DETECTIONS"); + if (!detections) { + detections = pmDetectionsAlloc(); + detections->allSources = psArrayAllocEmpty (galaxies->n); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY | PS_META_REPLACE, "psphot detections", detections); + } else { + psMemIncrRefCounter (detections); + } + psArray *sources = sources = detections->allSources; + + // output filename + char outname[1024]; + char *outroot = psMetadataLookupStr(&mdok, config->arguments, "OUTPUT"); + sprintf (outname, "%s.dat", outroot); + + FILE *outfile = fopen (outname, "a"); + + // 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; + + if (type == pmModelClassGetType ("PS_MODEL_TRAIL")) { + PAR[PM_PAR_LENGTH] = galaxy->Rmaj; + PAR[PM_PAR_SIGMA] = galaxy->Rmin; + PAR[PM_PAR_THETA] = galaxy->theta; + } else { + psEllipseAxes axes; + axes.major = galaxy->Rmaj; + axes.minor = galaxy->Rmin; + axes.theta = galaxy->theta; + pmPSF_AxesToModel (PAR, axes, model->class->useReff); + } + psF64 Area = 2.0 * M_PI * galaxy->Rmaj * galaxy->Rmin; + + if (nParam == 8) { + PAR[PM_PAR_7] = galaxy->index; + } + + // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0) + float radius = model->class->modelRadius (model->params, 0.1*roughNoise); + radius = PS_MAX (radius, 1.0); + + // XXX the exp(-r^0.25) models can go way out if allowed... + radius = PS_MIN (radius, 300.0); + + // construct a source, with model flux pixels set, based on the model + pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_EXTENDED); + + galaxy->flux = model->class->modelFlux (model->params); + + // 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->psfMagErr = sqrt(Area*PS_SQR(roughNoise) + galaxy->flux) / galaxy->flux; + + // insert the source flux in the image + double sum1 = imageSum(source->pixels); + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + double sum2 = imageSum(source->pixels); + float flux = sum2 - sum1; + + // fprintf (stderr, "flux: %f, sum1: %f, sum2: %f, flux2: %f, dflux: %f\n", galaxy->flux, sum1, sum2, sum2 - sum1, galaxy->flux - sum2 + sum1); + + float par8 = (model->params->n == 8) ? model->params->data.F32[7] : 0.0; + fprintf (outfile, "%8.3f %8.3f %10.2f %2d %7.3f %5.3f %5.3f %5.3f %5.3f %5.3f\n", galaxy->x, galaxy->y, flux, 1, source->psfMag, source->psfMagErr, galaxy->Rmaj, galaxy->Rmin, galaxy->theta, par8); + + psArrayAdd (sources, 100,source); + } + fclose (outfile); + + // XXX many leaks in here, i think + psFree (detections); + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimInsertStars.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimInsertStars.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimInsertStars.c (revision 41095) @@ -0,0 +1,185 @@ +# include "ppSim.h" + +// Reset a pointer: free and set to NULL +#define RESET(PTR) \ + psFree(PTR); \ + PTR = NULL; + +double imageSum (psImage *image); + +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 + + bool psfConvolve = psMetadataLookupBool(NULL, recipe, "PSF.CONVOLVE"); // smooth the image with the PSF? + + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate + float nSigmaLim = psMetadataLookupF32(NULL, recipe, "STARS.SIGMA.LIM"); // significance of faintest stars + + 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); + + // output filename + char outname[1024]; + char *outroot = psMetadataLookupStr(&mdok, config->arguments, "OUTPUT"); + sprintf (outname, "%s.dat", outroot); + + FILE *outfile = fopen (outname, "w"); + psAssert (outfile, "cannot write output"); + + float radius = -1.0; + + // add sources to the readout image & variance + 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); + + // we must set the radius on the first object, and leave it there for everyone else + // (inserting flux with a variable radius introduces a flux-dependent bias) + if (radius < 0) { + radius = model->class->modelRadius (model->params, nSigmaLim * roughNoise); + radius = PS_MAX (radius, 1.0); + radius = PS_MIN (radius, 50.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, model->class->useReff); + float Area = 2.0 * M_PI * axes.major * axes.minor; + + source->psfMag = -2.5*log10(star->flux); + source->psfMagErr = sqrt(Area*PS_SQR(roughNoise) + flux) / flux; + + // set the expected model errors + model->dparams->data.F32[PM_PAR_I0] = source->psfMagErr * model->params->data.F32[PM_PAR_I0]; + + float par8 = (model->params->n == 8) ? model->params->data.F32[7] : 0.0; + float starFlux = 0.0; + + // if psfConvolve is TRUE, we will (elsewhere) convolve the image we a PSF + // in this case, simply place delta functions in the image + if (psfConvolve) { + readout->image->data.F32[(int)(yCell)][(int)(xCell)] += flux; + } else { + // insert the source flux in the image and measure inserted flux + // this is a difference of 2 large numbers : must use 'double' or we get floating + // point errors and trends on a scale of ~0.005 mags at > -10 mags instrumental + double sum1 = imageSum(source->pixels); + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + double sum2 = imageSum(source->pixels); + starFlux = sum2 - sum1; + + // insert the source flux in the noise image + 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->variance); + RESET(source->maskObj); + RESET(source->maskView); + RESET(source->modelFlux); + RESET(source->psfImage); + RESET(source->blends); + } + fprintf (outfile, "%8.3f %8.3f %10.5f %2d %7.3f %5.3f %5.3f %5.3f %5.3f %5.3f : %f %f\n", star->x, star->y, starFlux, 0, source->psfMag, source->psfMagErr, axes.major, axes.minor, axes.theta, par8, model->params->data.F32[PM_PAR_I0], radius); + + // add the sources to the source array + psArrayAdd (sources, 100,source); + psFree(source); // Drop reference + } + fclose (outfile); + + pmDetections *detections = pmDetectionsAlloc(); + detections->allSources = sources; + + // save detections on the readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "fake sources", detections)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + psFree(detections); + + // XXX many leaks in here, i think + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimLoadStars.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimLoadStars.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimLoadStars.c (revision 41095) @@ -0,0 +1,106 @@ +# 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; + + float x0fpa = 0.5*(bounds->x0 + bounds->x1); + float y0fpa = 0.5*(bounds->y0 + bounds->y1); + + psFree(bounds); + + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN", PS_DATA_F32 | PS_META_REPLACE, "", + ra0 - radius / cos(dec0)); + psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX", PS_DATA_F32 | PS_META_REPLACE, "", + ra0 + radius / cos(dec0)); + 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, "PPSIM.OUTPUT"); + 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); + + psProjection *proj = psProjectionAlloc(ra0, dec0, scale, scale, PS_PROJ_TAN); // Projection + + // 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 + psPlane plane; // Plane (detector) coordinates + psSphere sphere; // Sphere (sky) coordinates + sphere.r = star->ra; + sphere.d = star->dec; + psProject(&plane, &sphere, proj); + + // Apply rotation, make FPA center of boresite + star->x = cos(pa) * plane.x - sin(pa) * plane.y + x0fpa; + star->y = sin(pa) * plane.x + cos(pa) * plane.y + 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; + psFree(proj); + + pmLumFunc *lumfunc = psastroLuminosityFunction (refStars, NULL); + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimLoop.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimLoop.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimLoop.c (revision 41095) @@ -0,0 +1,232 @@ +#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); // 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 + + bool psfConvolve = psMetadataLookupBool(NULL, recipe, "PSF.CONVOLVE"); // smooth the image with the PSF? + + 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"); + } + + 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 stars + if (!ppSimMakeStarGrid (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 + + // 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; + } + + 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->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels + + psImageInit (readout->image, 0.0); + psImageInit (readout->variance, 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 variance or generated a fake image above, we need to + // build one here + if (!readout->variance) { + if (!pmReadoutGenerateVariance(readout, NULL, true)) { + psError (PS_ERR_UNKNOWN, false, "trouble creating variance"); + 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); + + // we have two options for generating images which have a PSF: + // 1) lay down stars with the PSF model applied : in this case, galaxies do NOT have the PSF + // 2) lay down delta functions for stars and smooth the image with a PSF : in this case, the stars land at integer pixel locations + if (psfConvolve) { + if (!ppSimSmoothReadout(readout, recipe)) ESCAPE (PS_ERR_UNKNOWN, "problem smoothing image"); + } + + done: + if (!ppSimAddNoise(readout->image, readout->variance, cell, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "problem adding noise"); + + if (!ppSimBadCTE(readout->image, config)) ESCAPE (PS_ERR_UNKNOWN, "problem inducing bad cte"); + + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeBias.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeBias.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeBias.c (revision 41095) @@ -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->variance; + + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeBiassec.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeBiassec.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeBiassec.c (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeDark.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeDark.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeDark.c (revision 41095) @@ -0,0 +1,48 @@ +# include "ppSim.h" + +// XXX add bounds to the inputs? +bool ppSimMakeDark (pmReadout *readout, pmConfig *config) { + + bool mdok; + + psImage *signal = readout->image; + psImage *variance = readout->variance; + + 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 detTemp = psMetadataLookupF32(&mdok, recipe, "DET.TEMP"); // Detector Temperature + float darkTemp = psMetadataLookupF32(&mdok, recipe, "DARK.TEMP"); // Dark temperature trend + float darkTemp2 = psMetadataLookupF32(&mdok, recipe, "DARK.TEMP2"); // Dark temperature trend 2nd order term + if (!isfinite(detTemp)) detTemp = 0.0; + if (!isfinite(darkTemp)) darkTemp = 0.0; + if (!isfinite(darkTemp2)) darkTemp2 = 0.0; + + float expTime = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time + if (isnan(expTime)) { + psError(PS_ERR_UNKNOWN, false, "missing EXPTIME\n"); + return false; + } + + // Dark current + float darkCurrent = (darkRate + darkTemp*detTemp + darkTemp2*PS_SQR(detTemp)) * expTime; // Dark current accumulated + + psTrace("ppSim", 6, "dark current: %f, darkRate: %f, darkTempTrend: %f, darkTempSecondOrder: %f, detTemp: %f, expTime: %f\n", darkCurrent, darkRate, darkTemp, darkTemp2, detTemp, expTime); + + for (int y = 0; y < signal->numRows; y++) { + for (int x = 0; x < signal->numCols; x++) { + signal->data.F32[y][x] += darkCurrent; + variance->data.F32[y][x] += darkCurrent; + } + } + return true; +} + Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeGalaxies.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeGalaxies.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeGalaxies.c (revision 41095) @@ -0,0 +1,255 @@ +# include "ppSim.h" + +bool ppSimSetGalaxyPeak (ppSimGalaxy *galaxy, pmModelType type, float index, float seeing); + +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 galaxyAlpha = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope + float brightMag = 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 + bool galaxyGridRandom = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID.RANDOM"); // Density of fakes + // float galaxyGridPeak = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.PEAK"); // peak flux of fakes + float galaxyGridMag = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.MAG"); // peak flux of fakes + + float galaxyGridXoff = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.XOFF"); // centroid offset + float galaxyGridYoff = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.YOFF"); // centroid offset + + 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 + galaxyThetaMax *= PS_RAD_DEG; + galaxyThetaMin *= PS_RAD_DEG; + + 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); + } + + // determine the galaxy model + char *modelName = psMetadataLookupStr(&mdok, recipe, "GALAXY.MODEL"); // galaxy model name + pmModelType type = pmModelClassGetType (modelName); + if (type == -1) { + psError (PS_ERR_UNKNOWN, false, "invalid model name"); + return false; + } + + if (galaxyDensity <= 0) return true; + + if (galaxyGridRandom) { + long A, B; + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + } + + // 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; + // } + + // faintest and brightest magnitudes for random galaxies + // float brightFlux = ppSimMagToFlux (brightMag, zp) * expTime; + float faintMag = brightMag + 5.0; + if (brightMag > faintMag) { + 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 refSum = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + float normLum = refSum / (galaxyAlpha * logf(10.0) * powf(10.0, (galaxyAlpha * brightMag))); + float normScale = normLum * galaxyAlpha * logf(10.0); + + // Total number of galaxy down to the faint flux end + long nTotal = 0; + if (galaxyAlpha < 0.1) { + nTotal = 100; + } else { + nTotal = normScale * powf (10.0, (galaxyAlpha * faintMag)); + } + + if (galaxyGrid) { + fprintf (stderr, "galaxy grid inst mag: %f\n", galaxyGridMag - zp - 2.5*log10(expTime)); + nTotal = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY); + + psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", nTotal); + + float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin); + float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin); + float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin); + float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin); + + int i = 0; + // float refNorm = 1.0; + // float ourNorm = 1.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 + // center the galaxy peak on the center of a pixel + galaxy->x = ix + 0.5 + galaxyGridXoff; + galaxy->y = iy + 0.5 + galaxyGridYoff; + + float galaxyGridFlux = expTime * powf(10.0, -0.4 * (galaxyGridMag - zp)); + + // galaxyIndex from user should be for function of this form: exp(-r^(1/n)) + // galaxy->index = (1/2n) + + float rndValue; + + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; + float index = (galaxyIndexMin + rndValue * galaxyIndexSlope); + galaxy->index = 0.5/index; // factor of 0.5 because the Sersic model creates exp(-z^n), not exp(-r^n) + + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; + float scale = (galaxyRmajorMin + rndValue * galaxyRmajorSlope); + + galaxy->Rmaj = scale; + + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; + galaxy->Rmin = (galaxyARatioMin + rndValue * galaxyARatioSlope) * galaxy->Rmaj; + + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; + galaxy->theta = (galaxyThetaMin + rndValue * galaxyThetaSlope); + + galaxy->flux = galaxyGridFlux; + ppSimSetGalaxyPeak (galaxy, type, index, seeing); + + // if (0) { + // fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f, theta: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io, galaxy->theta); + // } + + psArrayAdd (galaxies, 100, galaxy); + i++; + } + } + } else { + + fprintf (stderr, "galaxy ref inst mag: %f\n", brightMag - zp - 2.5*log10(expTime)); + 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", nTotal, faintMag); + + for (long i = 0; i < nTotal; i++) { + ppSimGalaxy *galaxy = ppSimGalaxyAlloc (); + + // make fpa center of distribution + galaxy->x = psRandomUniform(rng) * xSize; // x position + galaxy->y = psRandomUniform(rng) * ySize; // y position + + float index = (sqrt(psRandomUniform(rng)) * galaxyIndexSlope + galaxyIndexMin); + galaxy->index = 0.5/index; // factor of 0.5 because the Sersic model creates exp(-z^n), not exp(-r^n) + + float scale = (psRandomUniform(rng) * galaxyRmajorSlope + galaxyRmajorMin); + galaxy->Rmaj = scale; + galaxy->Rmin = (PS_SQR(psRandomUniform(rng)) * galaxyARatioSlope + galaxyARatioMin) * galaxy->Rmaj; + galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope + galaxyThetaMin); + + // XXX probably should limit the allowed thinness of a galaxy + float galaxyMag = PS_MIN (faintMag, PS_MAX (brightMag, (log10((i + 1.0) / normScale) / galaxyAlpha))); + + galaxy->flux = ppSimMagToFlux (galaxyMag, zp) * expTime; + ppSimSetGalaxyPeak (galaxy, type, index, seeing); + + if (1) { + fprintf (stderr, "Rmaj: %f, Rmin: %f, theta: %f, scale: %f, index: %f, peak: %f, mag: %f\n", galaxy->Rmaj, galaxy->Rmin, galaxy->theta*180.0/M_PI, scale, index, galaxy->peak, -2.5*log10(galaxy->flux)); + } + + psArrayAdd (galaxies, 100, galaxy); + } + } + return true; +} + +bool ppSimSetGalaxyPeak (ppSimGalaxy *galaxy, pmModelType type, float index, float seeing) { + + bool isSersicType = false; + if (type == pmModelClassGetType ("PS_MODEL_DEV")) { + index = 4.0; + isSersicType = true; + } + if (type == pmModelClassGetType ("PS_MODEL_EXP")) { + index = 1.0; + isSersicType = true; + } + if (type == pmModelClassGetType ("PS_MODEL_SERSIC")) { + isSersicType = true; + } + + // exponential + // f = I0 exp (sqrt(z)) (index = 1.0) + + // f = I0 exp (-0.5*z) + + if (isSersicType) { + // for a sersic model, + // float bn = 1.9992*index - 0.3271; + // float Io = exp(bn); + + // the integral of a Sersic has an analytical form as follows: + // float logGamma = lgamma(2.0*index); + // float bnFactor = pow(bn, 2.0*index); + // float norm = 2.0 * M_PI * PS_SQR(galaxy->Rmaj) * index * Io * exp(logGamma) / bnFactor; + + // find the flux of a sersic with peak = 1.0: + float norm = pmSersicNorm (index); + float flux = 2.0 * M_PI * galaxy->Rmaj * galaxy->Rmin * norm; + + // XXX probably should limit the allowed thinness of a galaxy + galaxy->peak = galaxy->flux / flux; + return true; + } + + if (type == pmModelClassGetType ("PS_MODEL_TRAIL")) { + // galaxy->Rmaj -> PM_PAR_LENGTH + // seeing -> PM_PAR_SIGMA + galaxy->Rmin = seeing; // Rmin is used for sigma + galaxy->peak = galaxy->flux / (galaxy->Rmaj * galaxy->Rmin * 2.0 * sqrt(2.0 * M_PI)); + return true; + } + + galaxy->peak = galaxy->flux / (2 * M_PI * galaxy->Rmaj * galaxy->Rmin); + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeSky.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeSky.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeSky.c (revision 41095) @@ -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->variance; + + 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 + if (type == PPSIM_TYPE_FLAT) { + scatterRate = scatterFrac * PS_SQR(xFPA); + } + + // 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) { + // exposure correction is (effective exposure time) * (flatValue) + expCorr->data.F32[y][x] = flatValue * realExpTime / expTime; + } + + // TO DO: Add fringes + + } + } + psFree(bounds); + return true; +} + Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeStarGrid.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeStarGrid.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeStarGrid.c (revision 41095) @@ -0,0 +1,70 @@ +# include "ppSim.h" + +# define ESCAPE(MSG) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, MSG); \ + return false; } + +bool ppSimMakeStarGrid(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) { + + bool status; + assert (stars); + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + bool starsGrid = psMetadataLookupBool(&status, recipe, "STARS.GRID"); // Grid of stars + if (!starsGrid) return true; + + float starMag = psMetadataLookupF32(&status, recipe, "STARS.GRID.MAG"); // magnitude of grid stars + float dX = psMetadataLookupS32(&status, recipe, "STARS.GRID.DX"); // grid spacing in X + float dY = psMetadataLookupS32(&status, recipe, "STARS.GRID.DY"); // grid spacing in Y + + 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(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"); + + // 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); + + int nX = xSize / dX; + int oX = 0.5*dX; + + int nY = ySize / dY; + int oY = 0.5*dY; + + psLogMsg("ppSim", PS_LOG_INFO, "Adding %d stars with %f mag\n", nX * nY, starMag); + + long oldSize = stars->n; + psArrayRealloc (stars, stars->n + (nX * nY)); + + int i = 0; + for (int ix = 0; ix < nX; ix ++) { + for (int iy = 0; iy < nY; iy ++) { + + ppSimStar *star = ppSimStarAlloc (); + + // make fpa center of distribution + star->x = oX + dX*ix; // x position + star->y = oY + dY*iy; // y position + + star->flux = ppSimMagToFlux (starMag, zp) * expTime; + star->peak = ppSimStarFluxToPeak (star->flux, seeing); + + stars->data[oldSize + i] = star; + i++; + } + } + stars->n = oldSize + i; + + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeStars.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeStars.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMakeStars.c (revision 41095) @@ -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 = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + } + } 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimMergeReadouts.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMergeReadouts.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMergeReadouts.c (revision 41095) @@ -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->variance) { + if (!pmReadoutGenerateVariance(inReadout, NULL, true)) { + psError (PS_ERR_UNKNOWN, false, "trouble creating variance"); + 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->variance; + + psImage *outSignal = outReadout->image; + psImage *outVariance = outReadout->variance; + + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimMergeSources.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMergeSources.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMergeSources.c (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/src/ppSimMosaicChip.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMosaicChip.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimMosaicChip.c (revision 41095) @@ -0,0 +1,35 @@ +#include "ppSim.h" + +// XXX this is essentially identical to ppImageMosaicChip +bool ppSimMosaicChip(pmConfig *config, const psImageMaskType 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) { + pmFPAAddSourceFromView(out->fpa, 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimRandomGaussian.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimRandomGaussian.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimRandomGaussian.c (revision 41095) @@ -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 (void) { + + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimSaturate.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSaturate.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSaturate.c (revision 41095) @@ -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/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequence.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequence.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequence.c (revision 41095) @@ -0,0 +1,255 @@ +# include "ppSimSequence.h" +# include + +int main (int argc, char **argv) { + + bool status; + int argNum; + 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); + psString line = NULL; // Line to execute + psStringAppend(&line, "mkdir -p %s", path); + if (!system(line)) { + psWarning("Unable to create directory %s", path); + } + psFree(line); + } + + char *camera = NULL; + if ((argNum = psArgumentGet (argc, argv, "-camera"))) { + psArgumentRemove(argNum, &argc, argv); + camera = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } + + char *workdir = NULL; + if ((argNum = psArgumentGet (argc, argv, "-workdir"))) { + psArgumentRemove(argNum, &argc, argv); + workdir = psStringCopy (argv[argNum]); + psArgumentRemove(argNum, &argc, argv); + } + + char *ppsim_recipe = NULL; + if ((argNum = psArgumentGet (argc, argv, "-ppsim_recipe"))) { + psArgumentRemove(argNum, &argc, argv); + ppsim_recipe = 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, " -camera (camera) [otherwise must be set in sequences file]\n"); + fprintf (stderr, " -ppsim_recipe (recipe)\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); + + // build the base ppSimCommand string + psString ppSimCommand = psStringCopy ("ppSim"); + if (ppsim_recipe) psStringAppend (&ppSimCommand, " -recipe PPSIM %s", ppsim_recipe); + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + + // global camera option (if not set, we look in each sequence) + if (camera == NULL) { + camera = psMetadataLookupStr (&status, config, "CAMERA"); + } + + psArray *files = NULL; + psArray *sequences = NULL; + + { // find the FILERULE (if it exists) to set the file extension (.fits is default) + psMetadataItem *filerule = psMetadataLookup (config, "FILERULE"); + if (filerule == NULL) { + psLogMsg ("ppSimSequence", PS_LOG_INFO, "no FILERULE, assuming .fits ending"); + files = psArrayAlloc(1); + files->data[0] = psStringCopy("fits"); + goto sequence; + } + if (filerule->type == PS_DATA_METADATA_MULTI) { + psArray *rules = psListToArray (filerule->data.list); + psAssert (rules, "failed to get array from list?"); + psAssert (rules->n > 1, "supposed to be multiple entries in the list?"); + files = psArrayAllocEmpty(rules->n); + for (int i = 0; i < rules->n; i++) { + psMetadataItem *item = rules->data[i]; + if (item->type != PS_DATA_STRING) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "invalid FILERULE type"); + exit (PS_EXIT_CONFIG_ERROR); + } + psArrayAdd (files, 16, item->data.str); + } + goto sequence; + } + if (filerule->type == PS_DATA_STRING) { + files = psArrayAlloc(1); + files->data[0] = psStringCopy(filerule->data.str); + goto sequence; + } + psLogMsg ("ppSimSequence", PS_LOG_WARN, "invalid FILERULE type"); + exit (PS_EXIT_CONFIG_ERROR); + } + +sequence: + { // find the set of sequences which define the ppSim data to be produced + psMetadataItem *item = psMetadataLookup (config, "SEQUENCE"); + if (item == NULL) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "missing SEQUENCE description"); + exit (PS_EXIT_CONFIG_ERROR); + } + + 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); + } + + // determine the camera for the sequence and define the ppSim command + if (camera == NULL) { + camera = psMetadataLookupStr (&status, sequence, "CAMERA"); + } + if (!camera) { + psLogMsg ("ppSimSequence", PS_LOG_WARN, "CAMERA is not defined"); + exit (1); + } + + psString injectCommandReal = NULL; + psString ppSimCommandReal = NULL; + + psStringAppend (&injectCommandReal, "%s --camera %s", injectCommand, camera); + psStringAppend (&ppSimCommandReal, "%s -camera %s", ppSimCommand, camera); + + if (!strcasecmp (type, "BIAS")) { + ppSimSequenceBias (simfile, inject, sequence, i, rng, path, basename, ppSimCommandReal, injectCommandReal, files); + psFree (injectCommandReal); + psFree (ppSimCommandReal); + continue; + } + if (!strcasecmp (type, "DARK")) { + ppSimSequenceDark (simfile, inject, sequence, i, rng, path, basename, ppSimCommandReal, injectCommandReal, files); + psFree (injectCommandReal); + psFree (ppSimCommandReal); + continue; + } + if (!strcasecmp (type, "FLAT")) { + ppSimSequenceFlat (simfile, inject, sequence, i, rng, path, basename, ppSimCommandReal, injectCommandReal, files); + psFree (injectCommandReal); + psFree (ppSimCommandReal); + continue; + } + if (!strcasecmp (type, "OBJECT")) { + ppSimSequenceObject (simfile, inject, sequence, i, rng, path, basename, ppSimCommandReal, injectCommandReal, files); + psFree (injectCommandReal); + psFree (ppSimCommandReal); + continue; + } + + psLogMsg ("ppSimSequence", PS_LOG_WARN, "SEQUENCE %d has an unknown type: %s", i, type); + exit (1); + } + + exit (0); +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequence.h =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequence.h (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequence.h (revision 41095) @@ -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 *ppSimCommand, const char *injectCommand, psArray *files); +bool ppSimSequenceDark (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *ppSimCommand, const char *injectCommand, psArray *files); +bool ppSimSequenceFlat (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *ppSimCommand, const char *injectCommand, psArray *files); +bool ppSimSequenceObject (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *ppSimCommand, const char *injectCommand, psArray *files); + +#endif Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceBias.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceBias.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceBias.c (revision 41095) @@ -0,0 +1,53 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceBias (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *ppSimCommand, const char *injectCommand, psArray *files) { + + bool status, setLevel, setRange; + + // 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++) { + + // define the 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); + } + + // define the ppSim command + psString command = NULL; + + psStringAppend (&command, "%s -type BIAS", ppSimCommand); + + if (setLevel) psStringAppend (&command, " -biaslevel %f", level); + if (setRange) psStringAppend (&command, " -biasrange %f", range); + + psStringAppend (&command, " %s", filename); + + fprintf (simfile, "%s\n", command); + psFree (command); + + // define the inject command + // path should be dirname/filename + + // we use the filename above (really the file root) to construct the filenames + for (int i = 0; i < files->n; i++) { + command = psStringCopy (injectCommand); + psStringAppend (&command, " %s.%s", filename, (char *) files->data[i]); + fprintf (inject, "%s\n", command); + psFree (command); + } + + psFree (filename); + nImage ++; + } + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceDark.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceDark.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceDark.c (revision 41095) @@ -0,0 +1,71 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceDark (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *ppSimCommand, const char *injectCommand, psArray *files) { + + bool status, setRate; + float min, max = 0; + + 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++) { + + // 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); + } + + // define the ppSim command + psString command = NULL; + + psStringAppend (&command, "%s -type DARK", ppSimCommand); + + if (setRate) { + double frnd = psRandomUniform(rng); + float rate = min + (max - min)*frnd; + psStringAppend (&command, " -darkrate %f", rate); + } + + psStringAppend (&command, " -exptime %f", exptime); + + psStringAppend (&command, " %s", filename); + + fprintf (simfile, "%s\n", command); + psFree (command); + + // define the inject command + // path should be dirname/filename + + // we use the filename above (really the file root) to construct the filenames + for (int i = 0; i < files->n; i++) { + command = psStringCopy (injectCommand); + psStringAppend (&command, " %s.%s", filename, (char *) files->data[i]); + fprintf (inject, "%s\n", command); + psFree (command); + } + + psFree (filename); + nImage ++; + } + } + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceFlat.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceFlat.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceFlat.c (revision 41095) @@ -0,0 +1,66 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceFlat (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *ppSimCommand, const char *injectCommand, psArray *files) { + + bool status; + + // 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++) { + + // define the 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); + } + + // define the ppSim comand + psString command = NULL; + + psStringAppend (&command, "%s -type FLAT", ppSimCommand); + + psStringAppend (&command, " -filter %s", (char *) filters->data[i]); + psStringAppend (&command, " -exptime %f", exptimes->data.F32[i]); + + psStringAppend (&command, " %s", filename); + + fprintf (simfile, "%s\n", command); + psFree (command); + + // define the inject command + // path should be dirname/filename + + // we use the filename above (really the file root) to construct the filenames + for (int i = 0; i < files->n; i++) { + command = psStringCopy (injectCommand); + psStringAppend (&command, " %s.%s", filename, (char *) files->data[i]); + fprintf (inject, "%s\n", command); + psFree (command); + } + + psFree (filename); + nImage ++; + } + } + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceObject.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceObject.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSequenceObject.c (revision 41095) @@ -0,0 +1,129 @@ +# include "ppSimSequence.h" + +bool ppSimSequenceObject (FILE *simfile, FILE *inject, psMetadata *sequence, int nSeq, psRandom *rng, const char *path, const char *basename, const char *ppSimCommand, const char *injectCommand, psArray *files) { + + 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", + + // 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) { + + // define the 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); + } + + // define the ppSim command + psString command = NULL; + + psStringAppend (&command, "%s -type OBJECT", ppSimCommand); + 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); + psStringAppend (&command, " -obs_mode OBJECT.%s", (char *) filters->data[i]); + + double frnd = psRandomUniform(rng); + float seeing = IQmin + (IQmax - IQmin)*frnd; + + psStringAppend (&command, " -seeing %f", seeing); + + psStringAppend (&command, " %s", filename); + + fprintf (simfile, "%s\n", command); + psFree (command); + + // define the inject command + // path should be dirname/filename + + // we use the filename above (really the file root) to construct the filenames + for (int i = 0; i < files->n; i++) { + command = psStringCopy (injectCommand); + psStringAppend (&command, " %s.%s", filename, (char *) files->data[i]); + fprintf (inject, "%s\n", command); + psFree (command); + } + + psFree (filename); + nImage ++; + } + } + } + } + } + } + return true; +} Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSetPSF.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSetPSF.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSetPSF.c (revision 41095) @@ -0,0 +1,148 @@ +# 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) + float aRatio = psMetadataLookupF32(&status, recipe, "PSF.ARATIO"); // Seeing SIGMA (pixels) + float theta = psMetadataLookupF32(&status, recipe, "PSF.THETA"); // Seeing SIGMA (pixels) + + float seeingMax = psMetadataLookupF32(&status, recipe, "SEEING.MAX"); // Seeing SIGMA (pixels) + bool seeingRamp = psMetadataLookupBool(&status, recipe, "SEEING.RAMP"); // 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"); + + if (seeingRamp) { + psEllipseAxes axes; + psEllipsePol pol; + psEllipsePol polMax; + + // try spatial variation + options->psfTrendMode = PM_TREND_POLY_ORD; + options->psfTrendNx = 1; + options->psfTrendNy = 1; + options->psfFieldNx = xSize; + options->psfFieldNy = ySize; + + // generate the psf + psf = pmPSFAlloc (options); + + // supply the semi-major axis (these are SIGMA values in PIXELS) + axes.major = seeing; + axes.minor = aRatio * seeing; + axes.theta = theta * PS_RAD_DEG; + + pol = psEllipseAxesToPol (axes); + + axes.major = seeingMax; + axes.minor = aRatio * seeingMax; + polMax = psEllipseAxesToPol (axes); + + param = psf->params->data[PM_PAR_E0]; + param->poly->coeff[0][0] = pol.e0; + param->poly->coeff[1][0] = (polMax.e0 - pol.e0) / xSize; + param->poly->coeff[0][1] = (polMax.e0 - pol.e0) / ySize; + + param = psf->params->data[PM_PAR_E1]; + param->poly->coeff[0][0] = pol.e1; + param->poly->coeff[1][0] = 0.0; + param->poly->coeff[0][1] = 0.0; + + param = psf->params->data[PM_PAR_E2]; + param->poly->coeff[0][0] = pol.e2; + param->poly->coeff[1][0] = 0.0; + param->poly->coeff[0][1] = 0.0; + + if (!strcasecmp (psfModelName, "PS_MODEL_QGAUSS")) { + param = psf->params->data[PM_PAR_7]; + param->poly->coeff[0][0] = 1.0; + param->poly->coeff[1][0] = 0.0; + param->poly->coeff[0][1] = 0.0; + } + + if (!strcasecmp (psfModelName, "PS_MODEL_RGAUSS")) { + param = psf->params->data[PM_PAR_7]; + param->poly->coeff[0][0] = 1.0; + param->poly->coeff[1][0] = 0.0; + param->poly->coeff[0][1] = 0.0; + } + + } else { + psEllipseAxes axes; + psEllipsePol pol; + + // 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); + + // supply the semi-major axis (these are SIGMA values in PIXELS) + axes.major = seeing; + axes.minor = aRatio * seeing; + axes.theta = theta * PS_RAD_DEG; + + pol = psEllipseAxesToPol (axes); + psEllipseAxes testaxes = psEllipsePolToAxes (pol, 0.1); + fprintf (stderr, "psf in axes: %f x %f @ %f\n", axes.major, axes.minor, axes.theta*PS_DEG_RAD); + fprintf (stderr, "psf ot axes: %f x %f @ %f\n", testaxes.major, testaxes.minor, testaxes.theta*PS_DEG_RAD); + + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimSmoothReadout.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSmoothReadout.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimSmoothReadout.c (revision 41095) @@ -0,0 +1,49 @@ +#include "ppSim.h" + +bool ppSimSmoothReadout(pmReadout *input, psMetadata *recipe) +{ + bool status; + psTimerStart ("ppSmooth"); + + float nSigma = psMetadataLookupF32(&status, recipe, "CONVOLVE.NSIGMA"); // SIGMA convolutions (pixels) + if (!status) nSigma = 5.0; + + float sigma = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) + + char *modelName = psMetadataLookupStr(&status, recipe, "PSF.MODEL"); // Seeing SIGMA (pixels) + if (!strcmp (modelName, "PS_MODEL_GAUSS")) { + // bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve + // smooth the image in place, applying the mask as we go + psImageSmooth(input->image, sigma, nSigma); + psLogMsg("ppSmooth", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark("ppSmooth")); + return true; + } + + if (!strcmp (modelName, "PS_MODEL_PS1_V1")) { + // bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve + // smooth the image in place, applying the mask as we go + psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc(nSigma); + psImageSmooth2dCacheKernel_PS1_V1 (smdata, sigma, 0.2); + + psImageSmooth2dCache_F32 (input->image, smdata); + + psFree (smdata); + psLogMsg("ppSmooth", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark("ppSmooth")); + return true; + } + + // psImageConvolveSetThreads(oldThreads); + psLogMsg("ppSmooth", PS_LOG_MINUTIA, "failed to smooth image: %f sec\n", psTimerMark("ppSmooth")); + return false; +} + +// The variance is smoothed by the PSF^2, renomalized to maintain the input level of the +// variance. We achieve this by smoothing with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with +// unity normalization. Note that this process yields a smoothed image with correlated errors. +// The pixel-to-pixel variations in smooth_im will be decreased by a factor of +// 4*pi*SIGMA_SMTH^2, but for measurements based on apertures comparable to or larger than the +// smoothing kernel, the effective per-pixel variance is maintained. + + +// XXX this only allows for Gaussian PSFs. To extend this, we need to convolve with an image +// of the PSF model Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimStars.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimStars.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimStars.c (revision 41095) @@ -0,0 +1,59 @@ +# include "ppSim.h" + +void ppSimStarFree(ppSimStar *star) +{ + return; +} + +ppSimStar *ppSimStarAlloc(void) { + + ppSimStar *star = (ppSimStar *) psAlloc(sizeof(ppSimStar)); + psMemSetDeallocator(star, (psFreeFunc) ppSimStarFree); + + return star; +} + +void ppSimGalaxyFree(ppSimGalaxy *galaxy) +{ + return; +} + +ppSimGalaxy *ppSimGalaxyAlloc(void) { + + 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/ccl_branches/ipponly-20191108/ppSim/src/ppSimUtils.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimUtils.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimUtils.c (revision 41095) @@ -0,0 +1,303 @@ +# 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."); + // XXX raise error instead of aborting +} + +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 + float detTemp = psMetadataLookupF32(NULL, recipe, "DET.TEMP"); // Detector Temperature + + 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 + char *obs_mode = psMetadataLookupStr(NULL, recipe, "OBS_MODE"); + + 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); + psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.OBS.MODE", PS_META_REPLACE, "Observation mode", obs_mode); + psMetadataAddF32(fpa->concepts, PS_LIST_TAIL, "FPA.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", expTime); + psMetadataAddF32(fpa->concepts, PS_LIST_TAIL, "FPA.TEMP", PS_META_REPLACE, "Detector Temperature (C)", detTemp); + 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); + psMetadataAddF32(fpa->concepts, PS_LIST_TAIL, "FPA.AIRMASS", PS_META_REPLACE, "airmass (dummy value)", 1.0); + + // test time + psTime *obstime = psTimeFromMJD (55440.0); + psMetadataAddTime(fpa->concepts, PS_LIST_TAIL, "FPA.TIME", PS_META_REPLACE, "sample time", obstime); + + 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); + + // test time + psTime *obstime = psTimeFromMJD (55440.0); + psMetadataAddTime(cell->concepts, PS_LIST_TAIL, "CELL.TIME", PS_META_REPLACE, "sample time", obstime); + + 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 +// if it is not specified, do not override the existing recipe value +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; + } + if (!value) return NULL; + psMetadataAddStr(options, PS_LIST_TAIL, recipeName, PS_META_REPLACE, NULL, value); + return value; +} + +float ppSimGetZeroPoint(psMetadata *recipe, const char *filter) +{ + PS_ASSERT_METADATA_NON_NULL(recipe, NAN); + PS_ASSERT_STRING_NON_EMPTY(filter, NAN); + + // use the filter to get the zeropoint from the recipe + psMetadata *zeropoints = psMetadataLookupMetadata(NULL, recipe, "ZEROPTS"); + if (!zeropoints) { + psError(PS_ERR_UNKNOWN, false, "Unable to find ZEROPTS in recipe"); + return NAN; + } + + return psMetadataLookupF32(NULL, zeropoints, filter); +} + +psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) { + + bool status; + + pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename); + PS_ASSERT_PTR_NON_NULL (readout, NULL); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + 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; +} + Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimVersion.c =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimVersion.c (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimVersion.c (revision 41095) @@ -0,0 +1,42 @@ +#include "ppSim.h" +#include "ppSimVersionDefinitions.h" + +#ifndef PPSIM_VERSION +#error "PPSIM_VERSION is not set" +#endif +#ifndef PPSIM_BRANCH +#error "PPSIM_BRANCH is not set" +#endif +#ifndef PPSIM_SOURCE +#error "PPSIM_SOURCE is not set" +#endif + +psString ppSimVersion(void) +{ + char *value = NULL; + psStringAppend(&value, "%s@%s", PPSIM_BRANCH, PPSIM_VERSION); + return value; +} + +psString ppSimSource(void) +{ + return psStringCopy(PPSIM_SOURCE); +} + +psString ppSimVersionLong(void) +{ + psString version = ppSimVersion(); // Version, to return + psString source = ppSimSource(); // Source + + psStringPrepend(&version, "ppSim "); + psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__); + psFree(source); + +#ifdef __OPTIMIZE__ + psStringAppend(&version, " optimised"); +#else + psStringAppend(&version, " unoptimised"); +#endif + + return version; +}; Index: /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimVersionDefinitions.h.in =================================================================== --- /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimVersionDefinitions.h.in (revision 41095) +++ /branches/ccl_branches/ipponly-20191108/ppSim/src/ppSimVersionDefinitions.h.in (revision 41095) @@ -0,0 +1,8 @@ +#ifndef PPSIM_VERSION_DEFINITIONS_H +#define PPSIM_VERSION_DEFINITIONS_H + +#define PPSIM_VERSION @PPSIM_VERSION@ // SVN version +#define PPSIM_BRANCH @PPSIM_BRANCH@ // SVN branch +#define PPSIM_SOURCE @PPSIM_SOURCE@ // SVN source + +#endif