# include "ppSim.h"

bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng) {

    bool mdok;
    assert (galaxies);

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

    bool galaxyFake = psMetadataLookupBool(&mdok, recipe, "GALAXY.FAKE"); // Density of fakes
    if (!galaxyFake) return NULL;

    float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope
    float galaxyMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude
    float galaxyDensity = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes

    bool galaxyGrid = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID"); // Density of fakes
    int galaxyGridDX = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes
    int galaxyGridDY = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes
    
    float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes
    float galaxyRmajorMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MIN"); // Density of fakes

    float galaxyARatioMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MAX"); // Density of fakes
    float galaxyARatioMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MIN"); // Density of fakes

    float galaxyThetaMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MAX"); // Density of fakes
    float galaxyThetaMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MIN"); // Density of fakes

    float galaxyIndexMin  = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MIN"); // Density of fakes
    float galaxyIndexMax  = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MAX"); // Density of fakes

    // XXX push these into the recipe...
    float darkRate = psMetadataLookupF32(&mdok, config->arguments, "DARK.RATE"); // Dark rate
    float expTime = psMetadataLookupF32(&mdok, config->arguments, "EXPTIME"); // Exposure time

    float zp = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); // Photometric zero point
    float seeing = psMetadataLookupF32(&mdok, config->arguments, "SEEING"); // Seeing sigma (pix)
    float scale = psMetadataLookupF32(&mdok, 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 (galaxyDensity <= 0) return NULL;

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

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

    // Peak fluxes: faintest and brightest levels for random galaxy
    float faint = 0.1 * 10.0 * sqrt(2.0*M_PI) * seeing * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
    float bright = powf(10.0, -0.4 * (galaxyMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
    if (bright < faint) {
	psLogMsg("ppSim", PS_LOG_INFO,
		 "Image noise is above brightest random galaxy --- no random galaxy added.");
	return NULL;
    }

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

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

    if (galaxyGrid) {
	num = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY);
    
	psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", num);

	float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin) / num;
	float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin) / num;
	float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin) / num;
	float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin) / num;

	int i = 0;

	for (long iy = 0.5*galaxyGridDY; iy < ySize; iy += galaxyGridDY) {
	    for (long ix = 0.5*galaxyGridDX; ix < xSize; ix += galaxyGridDX) {
		ppSimGalaxy *galaxy = ppSimGalaxyAlloc ();

		// make fpa center of distribution
		galaxy->x    = ix;
		galaxy->y    = iy;

		galaxy->peak = 20000;

		galaxy->index = (galaxyIndexMin  + i * galaxyIndexSlope);
		galaxy->Rmaj  = (galaxyRmajorMin + i * galaxyRmajorSlope);
		galaxy->Rmin  = (galaxyARatioMin + i * galaxyARatioSlope) * galaxy->Rmaj;
		galaxy->theta = (galaxyThetaMin  + i * galaxyThetaSlope);

		psArrayAdd (galaxies, 100, galaxy);
		i++;
	    }
	}
    } else {    

	float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin);
	float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin);
	float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin);
	float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin);

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

	for (long i = 0; i < num; i++) {
	    ppSimGalaxy *galaxy = ppSimGalaxyAlloc ();

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

	    galaxy->flux = pow ((double)((i + 1) / norm), (double) (1.0 / galaxyLum));

	    galaxy->index = (psRandomUniform(rng) * galaxyIndexSlope  + galaxyIndexMin);
	    galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope  + galaxyThetaMin);
	    galaxy->Rmaj  = (psRandomUniform(rng) * galaxyRmajorSlope + galaxyRmajorMin);
	    galaxy->Rmin  = (psRandomUniform(rng) * galaxyARatioSlope + galaxyARatioMin) * galaxy->Rmaj;

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

	    psArrayAdd (galaxies, 100, galaxy);
	}
    }
    return galaxies;
}
