IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 17557


Ignore:
Timestamp:
May 7, 2008, 10:55:21 AM (18 years ago)
Author:
eugene
Message:

fix errors with star density normalizations and flux consistency

Location:
trunk/ppSim/src
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSim/src/ppSim.h

    r16493 r17557  
    133133    );
    134134
     135float ppSimStarSkyNoise (float skySigma, float seeingSigma);
     136float ppSimStarPeakToFlux (float peak, float seeingSigma);
     137float ppSimStarFluxToPeak (float flux, float seeingSigma);
     138float ppSimFluxToMag (float flux, float zp);
     139float ppSimMagToFlux (float mag, float zp);
    135140
    136141#endif

  • trunk/ppSim/src/ppSimArguments.c

    r16494 r17557  
    1818}
    1919
    20 // Get a value from the command-line arguments or the recipe, and add it to the target
    21 bool valueArgRecipe(pmConfig *config,   // Configuration
    22                     psMetadata *arguments, // Command-line arguments
    23                     const char *argName, // Argument name in the command-line arguments
    24                     const psMetadata *recipe, // Recipe
     20// Get a value from the command-line arguments and add it to recipe options
     21bool valueArgRecipe(psMetadata *options,    // Target to which to add value
    2522                    const char *recipeName, // Name for value in the recipe
    26                     psMetadata *target // Target to which to add value
     23                    psMetadata *arguments,  // Command-line arguments
     24                    const char *argName     // Argument name in the command-line arguments
    2725    )
    2826{
    29     float value = psMetadataLookupF32(NULL, arguments, argName); // Value of interest
    30     if (isnan(value)) {
    31         bool mdok;                      // Status of MD lookup
    32         value = psMetadataLookupF32(&mdok, recipe, recipeName);
    33         if (!mdok) {
    34             psErrorStackPrint(stderr, "Unable to find %s in recipe %s", recipeName, PPSIM_RECIPE);
    35             psFree((psPtr)arguments);
    36             psFree(config);
    37             exit(PS_EXIT_CONFIG_ERROR);
    38         }
    39     }
    40     return psMetadataAddF32(target, PS_LIST_TAIL, recipeName, 0, NULL, value);
     27    bool status;                                                    // Status of MD lookup
     28    float value = psMetadataLookupF32(&status, arguments, argName); // Value of interest
     29    if (isnan(value)) return true;
     30    status = psMetadataAddF32(options, PS_LIST_TAIL, recipeName, 0, NULL, value);
     31    return status;
    4132}
    4233
     34// this function supplements the RECIPE:OPTIONS folder with command-line options
    4335void ppSimArguments(int argc, char *argv[], pmConfig *config)
    4436{
     37    bool mdok;                          // Status of MD lookup
     38
    4539    assert(config);
    4640
     41    // save the following command-line options in the arguments structure
    4742    psMetadata *arguments = psMetadataAlloc(); // Command-line arguments
    4843    psMetadataAddStr(arguments, PS_LIST_TAIL, "-format", 0, "Camera format name", NULL);
    … …  
    6863    psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsdensity", 0, "Density of fake stars at magnitude", NAN);
    6964    psMetadataAddStr(arguments, PS_LIST_TAIL, "-psfclass", 0, "Type of PSF model", NULL);
     65    psMetadataAddStr(arguments, PS_LIST_TAIL, "-galmodel", 0, "Type of Galaxy model", NULL);
    7066    psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin", 0, "Binning in x and y", 1);
    7167
    … …  
    9288    }
    9389
     90    // apply an alternate camera format
    9491    psString formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format
    9592    if (formatName) {
    … …  
    114111    }
    115112
     113    // specify the type of simulated image to produce
     114    // XXX this should not be required if we supplied INPUT
    116115    const char *typeStr = psMetadataLookupStr(NULL, arguments, "-type"); // Exposure type
    117116    if (!typeStr) {
    … …  
    138137    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "FILTER", 0, "Filter name", filter);
    139138
     139    // save the following additional recipe values based on command-line options
     140    // these options override the PPSIM recipe values loaded from recipe files
     141    psMetadata *options = pmConfigRecipeOptions (config, PPSIM_RECIPE);
     142    if (!options) {
     143        psErrorStackPrint(stderr, "Unable to find recipe options for %s", PPSIM_RECIPE);
     144        psFree(arguments);
     145        psFree(config);
     146        exit(PS_EXIT_CONFIG_ERROR);
     147    }
     148
    140149    float expTime;
    141150    if (type == PPSIM_TYPE_BIAS) {
    … …  
    148157        }
    149158    }
    150     psMetadataAddF32(config->arguments, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", expTime);
    151 
    152     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe for ppSim
    153     if (!recipe) {
    154         psErrorStackPrint(stderr, "Unable to find recipe %s", PPSIM_RECIPE);
    155         psFree(arguments);
    156         psFree(config);
    157         exit(PS_EXIT_CONFIG_ERROR);
    158     }
    159 
    160     valueArgRecipe(config, arguments, "-biaslevel",    recipe, "BIAS.LEVEL",    config->arguments);
    161     valueArgRecipe(config, arguments, "-biasrange",    recipe, "BIAS.RANGE",    config->arguments);
    162     valueArgRecipe(config, arguments, "-darkrate",     recipe, "DARK.RATE",     config->arguments);
    163     valueArgRecipe(config, arguments, "-flatsigma",    recipe, "FLAT.SIGMA",    config->arguments);
    164     valueArgRecipe(config, arguments, "-flatrate",     recipe, "FLAT.RATE",     config->arguments);
    165     valueArgRecipe(config, arguments, "-shuttertime",  recipe, "SHUTTER.TIME",  config->arguments);
    166     valueArgRecipe(config, arguments, "-skyrate",      recipe, "SKY.RATE",      config->arguments);
    167     valueArgRecipe(config, arguments, "-starslum",     recipe, "STARS.LUM",     config->arguments);
    168     valueArgRecipe(config, arguments, "-starsmag",     recipe, "STARS.MAG",     config->arguments);
    169     valueArgRecipe(config, arguments, "-starsdensity", recipe, "STARS.DENSITY", config->arguments);
    170 
    171     bool mdok;                          // Status of MD lookup
     159    psMetadataAddF32(options, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", expTime);
     160
     161    // these values all get moved from arguments to RECIPE:OPTIONS
     162    valueArgRecipe(options, "BIAS.LEVEL",    arguments, "-biaslevel");
     163    valueArgRecipe(options, "BIAS.RANGE",    arguments, "-biasrange");
     164    valueArgRecipe(options, "DARK.RATE",     arguments, "-darkrate");
     165    valueArgRecipe(options, "FLAT.SIGMA",    arguments, "-flatsigma");
     166    valueArgRecipe(options, "FLAT.RATE",     arguments, "-flatrate");
     167    valueArgRecipe(options, "SHUTTER.TIME",  arguments, "-shuttertime");
     168    valueArgRecipe(options, "SKY.RATE",      arguments, "-skyrate");
     169    valueArgRecipe(options, "STARS.LUM",     arguments, "-starslum");
     170    valueArgRecipe(options, "STARS.MAG",     arguments, "-starsmag");
     171    valueArgRecipe(options, "STARS.DENSITY", arguments, "-starsdensity");
     172
     173    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
     174
    172175    int biasOrder = psMetadataLookupS32(&mdok, recipe, "BIAS.ORDER"); // Overscan polynomial order
    173176    if (!mdok) {
    174177        psWarning("BIAS.ORDER(S32) is not set in the recipe %s --- assuming %d", PPSIM_RECIPE, biasOrder);
    175178    }
    176     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BIAS.ORDER", 0,
    177                      "Overscan polynomial order", biasOrder);
     179    psMetadataAddS32(options, PS_LIST_TAIL, "BIAS.ORDER", 0, "Overscan polynomial order", biasOrder);
    178180
    179181    int binning = psMetadataLookupS32(NULL, arguments, "-bin"); // Binning in x and y
    … …  
    190192        float ra0     = psMetadataLookupF32(NULL, arguments, "-ra"); // Right Ascension of boresight
    191193        float dec0    = psMetadataLookupF32(NULL, arguments, "-dec"); // Declination of boresight
    192         float pa      = psMetadataLookupF32(NULL, arguments, "-pa"); // Position angle
    193         float seeing  = psMetadataLookupF32(NULL, arguments, "-seeing"); // Zero point
     194        float pa      = psMetadataLookupF32(NULL, arguments, "-pa");  // Position angle
     195        float seeing  = psMetadataLookupF32(NULL, arguments, "-seeing"); // seeing (FWHM in arcsec)
    194196
    195197        // XXX scale and zp should be supplied by the config file (allow override, but this is camera-dependent)
    … …  
    199201        }
    200202
    201         float zp      = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point
     203        float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point
    202204        if (isnan(zp)) {
    203205            // use the filter to get the zeropoint from the recipe
    … …  
    234236        }
    235237
    236         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0);
    237         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)", dec0 * M_PI / 180.0);
    238         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",pa * M_PI / 180.0);
    239         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SEEING", 0, "Seeing sigma (pix)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale);
    240 
    241         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale);
    242         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp);
    243         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SKY.MAGS", 0, "sky surface brightness", skymags);
    244 
    245         const char *psfClass = psMetadataLookupStr(NULL, arguments, "-psfclass"); // Filter name
     238        psMetadataAddF32(options, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0);
     239        psMetadataAddF32(options, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)", dec0 * M_PI / 180.0);
     240        psMetadataAddF32(options, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",pa * M_PI / 180.0);
     241
     242        // the user supplies FWHM in arcsec; here we convert to Sigma in pixels
     243        psMetadataAddF32(options, PS_LIST_TAIL, "SEEING", 0, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale);
     244
     245        psMetadataAddF32(options, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale);
     246        psMetadataAddF32(options, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp);
     247        psMetadataAddF32(options, PS_LIST_TAIL, "SKY.MAGS", 0, "sky surface brightness", skymags);
     248
     249        const char *psfClass = psMetadataLookupStr(NULL, arguments, "-psfclass"); // PSF model class
    246250        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "PSF.MODEL", 0, "PSF model class", psfClass);
     251
     252        const char *galModel = psMetadataLookupStr(NULL, arguments, "-galmodel"); // Galaxy model name
     253        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "GALAXY.MODEL", 0, "Galaxy model", galModel);
    247254    }
    248255
    … …  
    250257    return;
    251258}
     259
     260/* the following elements come from the config->arguments:
     261   
     262   PSPHOT.PSF
     263   INPUT
     264   TYPE
     265   FILTER
     266   BIAS.ORDER
     267   BINNING
     268   OUTPUT
     269   PSF.MODEL
     270   GALAXY.MODEL
     271
     272   all othr values should come from the recipe
     273*/

  • trunk/ppSim/src/ppSimCreate.c

    r14667 r17557  
    1818    pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPIMAGE.INPUT", "INPUT");
    1919    if (!input) {
     20        // if we have not specified the camera already, we need to interpolate the recipes associated with this camera, and read other command-line recipes
     21        if (!pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) {
     22            psError(PS_ERR_IO, false, "Error merging recipes from camera config for %s", config->cameraName);
     23            return NULL;
     24        }
     25
    2026        simImage = true;
    2127        fpa = pmFPAConstruct(config->camera); // FPA to contain the observation
    … …  
    2430            return NULL;
    2531        }
     32       
    2633    } else {
    2734        simImage = false;
    … …  
    3340    }
    3441
     42    // define the output image file
    3543    pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE);
    3644    if (!file) {

  • trunk/ppSim/src/ppSimInsertGalaxies.c

    r16496 r17557  
    2222    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
    2323
    24     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
    25     float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
     24    float expTime   = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time
     25    float darkRate  = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate
     26
    2627    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
    2728    if (isnan(readnoise)) {
    … …  
    3435    }
    3536
    36     float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
     37    float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate
    3738    if (isnan(skyRate)) {
    38         float zp      = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); assert (mdok);
    39         float scale   = psMetadataLookupF32(&mdok, config->arguments, "SCALE");     assert (mdok);
    40         float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS");  assert (mdok);
    41         skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp));
     39        float zp      = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok);
     40        float scale   = psMetadataLookupF32(&mdok, recipe, "SCALE");     assert (mdok);
     41        float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
     42        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    4243    }
    4344   
    … …  
    5859
    5960    // determine the galaxy model
    60     char *modelName = psMetadataLookupStr(&mdok, config->arguments, "GALAXY.MODEL"); // Seeing sigma (pix)
     61    char *modelName = psMetadataLookupStr(&mdok, config->arguments, "GALAXY.MODEL"); // galaxy model name
    6162    if (modelName == NULL) {
    6263        modelName = defaultModel;

  • trunk/ppSim/src/ppSimInsertStars.c

    r16496 r17557  
    2727    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
    2828
    29     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
    30     float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
     29    float expTime   = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time
     30    float darkRate  = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate
     31
    3132    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
    3233    if (isnan(readnoise)) {
    … …  
    3940    }
    4041
    41     float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
     42    float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate
    4243    if (isnan(skyRate)) {
    43         float zp      = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); assert (mdok);
    44         float scale   = psMetadataLookupF32(&mdok, config->arguments, "SCALE");     assert (mdok);
    45         float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS");  assert (mdok);
    46         skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp));
     44        float zp      = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok);
     45        float scale   = psMetadataLookupF32(&mdok, recipe, "SCALE");     assert (mdok);
     46        float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
     47        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    4748    }
    4849

  • trunk/ppSim/src/ppSimLoadStars.c

    r14813 r17557  
    1818    }
    1919
    20     float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
    21     float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians)
    22     float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
    23     float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)
    24     float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
    25     float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
    26     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     20    // XXX push these into the recipe in ppSimArguments()
     21    float zp      = psMetadataLookupF32(NULL, recipe, "ZEROPOINT"); // Photometric zero point
     22    float ra0     = psMetadataLookupF32(NULL, recipe, "RA");        // Boresight RA (radians)
     23    float dec0    = psMetadataLookupF32(NULL, recipe, "DEC");       // Boresight Dec (radians)
     24    float pa      = psMetadataLookupF32(NULL, recipe, "PA");        // Position angle (radians)
     25    float seeing  = psMetadataLookupF32(NULL, recipe, "SEEING");    // Seeing SIGMA (pixels)
     26    float scale   = psMetadataLookupF32(NULL, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
     27    float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME");   // Exposure time (sec)
    2728
    2829    // Size of FPA
    … …  
    6869
    6970        // Convert magnitude to peak flux
    70         star->flux = powf(10.0, -0.4 * (star->mag - zp)) * expTime;
    71         star->peak = star->flux / (2.0*M_PI * PS_SQR(seeing));
     71        star->flux = ppSimMagToFlux (star->mag, zp) * expTime;
     72        star->peak = ppSimStarFluxToPeak (star->flux, seeing);
     73
    7274        stars->data[oldSize + i] = star;
    7375
    … …  
    7577    }
    7678    stars->n = oldSize + refStars->n;
     79
     80    pmLumFunc *lumfunc = psastroLuminosityFunction (refStars);
    7781    psFree(refStars);
    7882
    7983    psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0);
    8084    psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0);
     85   
     86    if (lumfunc) {
     87        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MIN",   PS_META_REPLACE, "min valid magnitude",             lumfunc->mMin);
     88        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MAX",   PS_META_REPLACE, "max valid magnitude",             lumfunc->mMax);
     89        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.SLOPE",  PS_META_REPLACE, "log-mag histogram slope",         lumfunc->slope);
     90        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.OFFSET", PS_META_REPLACE, "log-mag histogram offset",        lumfunc->offset);
     91        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.PEAK",  PS_META_REPLACE, "magnitude of peak bin",           lumfunc->mPeak);
     92        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.NUM.PEAK",  PS_META_REPLACE, "number of stars in peak bin", lumfunc->nPeak);
     93        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.SUM.PEAK",  PS_META_REPLACE, "sum of stars up to peak bin", lumfunc->sPeak);
     94        psFree (lumfunc);
     95    }
    8196
    8297    return stars;

  • trunk/ppSim/src/ppSimMakeBias.c

    r14667 r17557  
    99    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
    1010
    11     float biasLevel = psMetadataLookupF32(NULL, config->arguments, "BIAS.LEVEL"); // Bias level
    12     float biasRange = psMetadataLookupF32(NULL, config->arguments, "BIAS.RANGE"); // Bias range
    13     int biasOrder = psMetadataLookupS32(NULL, config->arguments, "BIAS.ORDER"); // Bias order
     11    float biasLevel = psMetadataLookupF32(NULL, recipe, "BIAS.LEVEL"); // Bias level
     12    float biasRange = psMetadataLookupF32(NULL, recipe, "BIAS.RANGE"); // Bias range
     13    int biasOrder   = psMetadataLookupS32(NULL, recipe, "BIAS.ORDER"); // Bias order
    1414
    1515    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e

  • trunk/ppSim/src/ppSimMakeDark.c

    r14657 r17557  
    44bool ppSimMakeDark (pmReadout *readout, pmConfig *config) {
    55
     6    bool mdok;
     7
    68    psImage *signal = readout->image;
    79    psImage *variance = readout->weight;
    810
    9     float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
    10     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     11    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
     12
     13    float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate
     14    float expTime  = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time
    1115
    1216    for (int y = 0; y < signal->numRows; y++) {

  • trunk/ppSim/src/ppSimMakeGalaxies.c

    r16498 r17557  
    1111    if (!galaxyFake) return NULL;
    1212
    13     float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope
    14     float galaxyMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude
    15     float galaxyDensity = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes
     13    float galaxyLum       = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope
     14    float galaxyMag       = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude
     15    float galaxyDensity   = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes
    1616
    17     bool galaxyGrid = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID"); // Density of fakes
    18     int galaxyGridDX = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes
    19     int galaxyGridDY = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes
     17    bool galaxyGrid       = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID"); // Density of fakes
     18    int galaxyGridDX      = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes
     19    int galaxyGridDY      = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes
    2020   
    2121    float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes
    … …  
    2525    float galaxyARatioMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MIN"); // Density of fakes
    2626
    27     float galaxyThetaMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MAX"); // Density of fakes
    28     float galaxyThetaMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MIN"); // Density of fakes
     27    float galaxyThetaMax  = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MAX"); // Density of fakes
     28    float galaxyThetaMin  = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MIN"); // Density of fakes
    2929
    3030    float galaxyIndexMin  = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MIN"); // Density of fakes
    3131    float galaxyIndexMax  = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MAX"); // Density of fakes
    3232
    33     // XXX push these into the recipe...
    34     float darkRate = psMetadataLookupF32(&mdok, config->arguments, "DARK.RATE"); // Dark rate
    35     float expTime = psMetadataLookupF32(&mdok, config->arguments, "EXPTIME"); // Exposure time
    36 
    37     float zp = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); // Photometric zero point
    38     float seeing = psMetadataLookupF32(&mdok, config->arguments, "SEEING"); // Seeing sigma (pix)
    39     float scale = psMetadataLookupF32(&mdok, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
    40 
    41     float skyRate = psMetadataLookupF32(&mdok, config->arguments, "SKY.RATE"); // Sky rate
     33    float darkRate        = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate
     34    float expTime         = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time
     35    float zp              = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point
     36    float seeing          = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix)
     37    float scale           = psMetadataLookupF32(&mdok, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
     38    float skyRate         = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate
    4239    if (isnan(skyRate)) {
    43         float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS");  assert (mdok);
    44         skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp));
     40        float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
     41        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    4542    }
    4643

  • trunk/ppSim/src/ppSimMakeSky.c

    r16497 r17557  
    1313    pmFPA  *fpa  = chip->parent;
    1414
    15     float expTime     = psMetadataLookupF32(&status, config->arguments, "EXPTIME"); // Exposure time
    16     float flatSigma   = psMetadataLookupF32(&status, config->arguments, "FLAT.SIGMA"); // Flat-field coefficient
    17     float flatRate    = psMetadataLookupF32(&status, config->arguments, "FLAT.RATE"); // Flat-field rate
    18     float shutterTime = psMetadataLookupF32(&status, config->arguments, "SHUTTER.TIME"); // Shutter time
     15    psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe
    1916
    20     float skyRate     = psMetadataLookupF32(&status, config->arguments, "SKY.RATE"); // Sky rate
     17    float expTime     = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time
     18    float flatSigma   = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient
     19    float flatRate    = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate
     20    float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time
     21    float skyRate     = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate
    2122    if (isnan(skyRate)) {
    22         float zp      = psMetadataLookupF32(&status, config->arguments, "ZEROPOINT"); assert (status);
    23         float scale   = psMetadataLookupF32(&status, config->arguments, "SCALE");     assert (status);
    24         float skyMags = psMetadataLookupF32(&status, config->arguments, "SKY.MAGS");  assert (status);
    25         skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp));
     23        float zp      = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status);
     24        float scale   = psMetadataLookupF32(&status, recipe, "SCALE");     assert (status);
     25        float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS");  assert (status);
     26        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    2627    }
    2728

  • trunk/ppSim/src/ppSimMakeStars.c

    r16498 r17557  
    11# include "ppSim.h"
    2 
    3 #define FAINT_FUDGE_FACTOR   0.4        // Fraction of the noise in a (boxcar) PSF for the faint limit
    4 
    52
    63bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
    74
    8     bool mdok;
     5    bool status;
    96    assert (stars);
    107
    11     psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
     8    psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe
    129
    13     bool starsFake = psMetadataLookupBool(&mdok, recipe, "STARS.FAKE"); // Density of fakes
     10    bool starsFake = psMetadataLookupBool(&status, recipe, "STARS.FAKE"); // Density of fakes
    1411    if (!starsFake) return true;
    1512
    16     float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
    17     float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude
    18     float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes
     13    bool starsReal     = psMetadataLookupBool(&status, recipe, "STARS.REAL"); // were real stars generated?
     14    bool matchDensity  = psMetadataLookupBool(&status, recipe, "MATCH.DENSITY"); // match real star density?
    1915
    20     float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
    21     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     16    float starsAlpha   = psMetadataLookupF32(&status, recipe, "STARS.LUM"); // Star luminosity func slope
     17    float starsDensity = psMetadataLookupF32(&status, recipe, "STARS.DENSITY"); // Density of fakes
     18    float brightMag    = psMetadataLookupF32(&status, recipe, "STARS.MAG"); // Star brightest magnitude
     19    float nSigmaLim    = psMetadataLookupF32(&status, recipe, "STARS.SIGMA.LIM"); // significance of faintest stars
    2220
    23     float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
    24     float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
    25     float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
     21    float darkRate     = psMetadataLookupF32(&status, recipe, "DARK.RATE"); // Dark rate
     22    float expTime      = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time
     23    float zp           = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point
     24    float seeing       = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels)
     25    float scale        = psMetadataLookupF32(&status, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
    2626
    27     float skyRate = psMetadataLookupF32(&mdok, config->arguments, "SKY.RATE"); // Sky rate
     27    float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate
    2828    if (isnan(skyRate)) {
    29         float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS");  assert (mdok);
    30         skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp));
     29        float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS");  assert (status);
     30        skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
    3131    }
    32 
    33     if (starsDensity <= 0) return true;
    3432
    3533    // Size of FPA
    … …  
    4038    psFree(bounds);
    4139
     40    // choose reference magnitude & density to set normalization
     41    float refMag = 0;
     42    float refSum = 0;
     43    if (starsReal && matchDensity) {
     44        refMag = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.MAG.PEAK"); // Star brightest magnitude
     45        refSum = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.SUM.PEAK"); // Star brightest magnitude
     46        assert (status);
     47    } else {
     48        refMag = brightMag;
     49        refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI);
     50    }
     51    psTrace("ppSim", 6, "refMag: %f, refSum: %f\n", refMag, refSum);
     52
     53    if (refSum <= 0) return true;
     54
    4255    // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it
    4356    // shouldn't be too bad.
    44     float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise
    45     if (!mdok) {
     57    float readnoise = psMetadataLookupF32(&status, recipe, "READNOISE"); // Default read noise
     58    if (!status) {
    4659        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
    4760        psFree(bounds);
    … …  
    5063
    5164    // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars
    52     float faint = FAINT_FUDGE_FACTOR * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime) *
    53         (2.0 * sqrt(2.0 * log(2.0)) * seeing); // Faint limit is related to the noise in a (boxcar) PSF
    54     float bright = powf(10.0, -0.4 * (starsMag - zp)) * expTime; // Bright limit is specified by user as mag
    55     psTrace("ppSim", 6, "Faint limit: %f\n", faint);
    56     psTrace("ppSim", 6, "Bright limit: %f\n", bright);
    57     if (bright < faint) {
     65
     66    // faint limit is set by detection limit, in turn set by sky noise
     67    float skySigma = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
     68    float skyNoise = ppSimStarSkyNoise (skySigma, seeing);
     69
     70    // total flux of faintest star:
     71    float faintCounts = nSigmaLim * skyNoise;
     72    float faintMag = ppSimFluxToMag ((faintCounts / expTime), zp);
     73
     74    float brightCounts = ppSimMagToFlux (brightMag, zp) * expTime; // Bright limit is specified by user as mag
     75
     76    psTrace("ppSim", 6, "Faint limit: %f counts = %f mags\n", faintCounts, faintMag);
     77    psTrace("ppSim", 6, "Bright limit: %f counts = %f mags\n", brightCounts, brightMag);
     78    if (brightCounts < faintCounts) {
    5879        psLogMsg("ppSim", PS_LOG_INFO,
    5980                 "Image noise is above brightest random star --- no random stars added.");
    … …  
    6182    }
    6283
    63     // Normalisation, set by the specified stellar density at the specified bright magnitude
    64     float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
    65         powf(bright, starsLum);
     84    // given log_10 (dN / dmag) = alpha mag + beta
     85    // or dN / dmag = No 10^(alpha mag), then:
     86    // N(m < Mo) = alpha * No * log(10.0) * 10^(alpha*Mo)
     87
     88    // XXX this needs to handle the case of a flat distribution for efficiency testing
     89
     90    // Normalization, set by the specified stellar density at the specified bright magnitude
     91    float norm = refSum / (starsAlpha * logf(10.0) * powf(10.0, (starsAlpha * refMag)));
     92    float normScale = norm * starsAlpha * logf(10.0);
    6693
    6794    // Total number of stars down to the faint flux end
    68     long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
     95    long nTotal = normScale * powf (10.0, (starsAlpha * faintMag));
    6996
    70     psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n",
    71              num, -2.5 * log10(faint / expTime) + zp);
     97    psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", nTotal, faintMag);
    7298
    7399    long oldSize = stars->n;
    74     psArrayRealloc (stars, stars->n + num);
     100    psArrayRealloc (stars, stars->n + nTotal);
    75101
    76     for (long i = 0; i < num; i++) {
     102    for (long i = 0; i < nTotal; i++) {
    77103        ppSimStar *star = ppSimStarAlloc ();
    78104
    … …  
    81107        star->y    = psRandomUniform(rng) * ySize; // y position
    82108
    83         star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
    84         star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));
     109        // XXX this needs to respect the bright limit
     110        float starMag = log10((i + 1.0) / normScale) / starsAlpha;
     111       
     112        star->flux = ppSimMagToFlux (starMag, zp) * expTime;
     113        star->peak = ppSimStarFluxToPeak (star->flux, seeing);
    85114
    86115        stars->data[oldSize + i] = star;
    87116    }
    88     stars->n = oldSize + num;
     117    stars->n = oldSize + nTotal;
    89118
    90119    return true;

  • trunk/ppSim/src/ppSimSequence.c

    r16499 r17557  
    11# include "ppSimSequence.h"
     2# include <sys/stat.h>
    23
    34int main (int argc, char **argv) {
    45
     6    bool status;
    57    int argNum;
    6     bool status;
     8    char line[1024];
    79    unsigned int nFail;
    810
    … …  
    1719
    1820    char *path = NULL;
    19     if (psArgumentGet (argc, argv, "-path")) {
     21    if ((argNum = psArgumentGet (argc, argv, "-path"))) {
    2022        psArgumentRemove(argNum, &argc, argv);
    2123        path = psStringCopy (argv[argNum]);
    2224        psArgumentRemove(argNum, &argc, argv);
     25        sprintf (line, "mkdir -p %s", path);
     26        system (line);
    2327    }
    2428
    2529    char *workdir = NULL;
    26     if (psArgumentGet (argc, argv, "-workdir")) {
     30    if ((argNum = psArgumentGet (argc, argv, "-workdir"))) {
    2731        psArgumentRemove(argNum, &argc, argv);
    2832        workdir = psStringCopy (argv[argNum]);
    … …  
    3135
    3236    char *basename = NULL;
    33     if (psArgumentGet (argc, argv, "-basename")) {
     37    if ((argNum = psArgumentGet (argc, argv, "-basename"))) {
    3438        psArgumentRemove(argNum, &argc, argv);
    3539        basename = psStringCopy (argv[argNum]);
    … …  
    130134        exit (1);
    131135    }
     136
    132137    exit (0);
    133138}

  • trunk/ppSim/src/ppSimSetPSF.c

    r15005 r17557  
    11# include "ppSim.h"
    2 static char *defaultModel = "PS_MODEL_GAUSS";
     2static char *defaultModel = "PS_MODEL_QGAUSS";
    33
    44bool ppSimSetPSF (pmChip *chip, pmConfig *config) {
    55
    6     bool status;
     6    bool status, mdok;
    77    pmPSF *psf = NULL;
    88    pmTrend2D *param = NULL;
     9
     10    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
    911
    1012    // the pmPSF IO functions stores the PSF on the chip->analysis
    … …  
    1416    }
    1517
    16     // no supplied PSF, build one using supplied value for seeing
    17     // XXX need to correct for the pixel scale
    18     float seeing   = psMetadataLookupF32(&status, config->arguments, "SEEING"); // Seeing sigma (pix)
     18    // no supplied PSF, build one using supplied value for seeing seeing is already corrected
     19    // for the pixel scale, and is converted from FWHM to SIGMA (this is done in
     20    // ppSimArguments)
     21    float seeing   = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels)
    1922
    2023    char *psfModelName = psMetadataLookupStr(&status, config->arguments, "PSF.MODEL"); // Name of PSF model
    … …  
    5053    psEllipsePol pol;
    5154
    52     // supply the semi-major axis
     55    // supply the semi-major axis (these are SIGMA values in PIXELS)
    5356    axes.major = seeing;
    5457    axes.minor = seeing;

  • trunk/ppSim/src/ppSimStars.c

    r14667 r17557  
    2626    return galaxy;
    2727}
     28
     29float ppSimStarSkyNoise (float skySigma, float seeingSigma) {
     30
     31    float skyNoise = skySigma * sqrt(4*M_PI*PS_SQR(seeingSigma));
     32    return skyNoise;
     33}
     34
     35float ppSimStarPeakToFlux (float peak, float seeingSigma) {
     36
     37    float psfArea = 2.0*M_PI*PS_SQR(seeingSigma);
     38    float flux = peak * psfArea;
     39    return flux;
     40}
     41
     42float ppSimStarFluxToPeak (float flux, float seeingSigma) {
     43
     44    float psfArea = 2.0*M_PI*PS_SQR(seeingSigma);
     45    float peak = flux / psfArea;
     46    return peak;
     47}
     48
     49float ppSimFluxToMag (float flux, float zp) {
     50
     51    float mag = -2.5*log10(flux) + zp;
     52    return mag;
     53}
     54
     55float ppSimMagToFlux (float mag, float zp) {
     56
     57    float flux = powf (10.0, -0.4*(mag - zp));
     58    return flux;
     59}

  • trunk/ppSim/src/ppSimUtils.c

    r14798 r17557  
    88                     pmCell *cell)
    99{
     10    bool mdok;
    1011
    11     float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians)
    12     float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
    13     float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)
    14     float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // plate scale in arcsec / pixel
     12    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
     13
     14    float ra0   = psMetadataLookupF32(NULL, recipe, "RA");  // Boresight RA (radians)
     15    float dec0  = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians)
     16    float pa    = psMetadataLookupF32(NULL, recipe, "PA");  // Position angle (radians)
     17    float scale = psMetadataLookupF32(NULL, recipe, "SCALE"); // plate scale in arcsec / pixel
    1518    scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel
    1619
    … …  
    9497bool ppSimUpdateConceptsFPA (pmFPA *fpa, pmConfig *config) {
    9598
    96     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     99    bool mdok;
     100
     101    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
     102
     103    float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time
    97104
    98105    const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name
    … …  
    128135bool ppSimUpdateConceptsCell (pmCell *cell, pmConfig *config) {
    129136
     137    bool mdok;
     138
     139    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
     140
    130141    int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y
    131     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     142    float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time
    132143
    133144    psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,
Note: See TracChangeset for help on using the changeset viewer.