# include "ppSim.h"

#define FAINT_FUDGE_FACTOR   0.4        // Fraction of the noise in a (boxcar) PSF for the faint limit


bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) {

    bool mdok;
    assert (stars);

    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe

    bool starsFake = psMetadataLookupBool(&mdok, recipe, "STARS.FAKE"); // Density of fakes
    if (!starsFake) return true;

    float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
    float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude
    float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes

    float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
    float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time

    float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)

    float skyRate = psMetadataLookupF32(&mdok, config->arguments, "SKY.RATE"); // Sky rate
    if (isnan(skyRate)) {
	float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS");  assert (mdok);
	skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp));
    }

    if (starsDensity <= 0) return true;

    // Size of FPA
    psRegion *bounds = ppSimFPABounds (fpa);
    // Size of focal plane (can't I get this from the config?)
    int xSize = bounds->x1 - bounds->x0;
    int ySize = bounds->y1 - bounds->y0;
    psFree(bounds);

    // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it
    // shouldn't be too bad.
    float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise
    if (!mdok) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
        psFree(bounds);
        return false;
    }

    // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars
    float faint = FAINT_FUDGE_FACTOR * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime) *
        (2.0 * sqrt(2.0 * log(2.0)) * seeing); // Faint limit is related to the noise in a (boxcar) PSF
    float bright = powf(10.0, -0.4 * (starsMag - zp)) * expTime; // Bright limit is specified by user as mag
    psTrace("ppSim", 6, "Faint limit: %f\n", faint);
    psTrace("ppSim", 6, "Bright limit: %f\n", bright);
    if (bright < faint) {
        psLogMsg("ppSim", PS_LOG_INFO,
                 "Image noise is above brightest random star --- no random stars added.");
        return true;
    }

    // Normalisation, set by the specified stellar density at the specified bright magnitude
    float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
        powf(bright, starsLum);

    // Total number of stars down to the faint flux end
    long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;

    psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n",
             num, -2.5 * log10(faint / expTime) + zp);

    long oldSize = stars->n;
    psArrayRealloc (stars, stars->n + num);

    for (long i = 0; i < num; i++) {
        ppSimStar *star = ppSimStarAlloc ();

        // make fpa center of distribution
        star->x    = psRandomUniform(rng) * xSize; // x position
        star->y    = psRandomUniform(rng) * ySize; // y position

        star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
        star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));

        stars->data[oldSize + i] = star;
    }
    stars->n = oldSize + num;

    return true;
}
