# 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;
        }
    }

    // generate the fpa structure used by the output camera (determined from INPUT or specified)
    assert (config->camera);
    fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the observation
    if (!fpa) {
        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
        return NULL;
    }

    // define the output image file -- this is the basis for the ppSimLoop
    pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSIM.OUTPUT");
    if (!output) {
        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from PPSIM.OUTPUT. Did you forget to specify the format?");
        return NULL;
    }
    if (output->type != PM_FPA_FILE_IMAGE) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "PPSIM.OUTPUT type is not IMAGE");
        psFree(fpa);
        return NULL;
    }
    // XXX we should not require the output image to be written
    output->save = true;

    config->format = psMemIncrRefCounter (output->format);
    config->formatName = psStringCopy (output->formatName);

    // the recipe is now fully realized for the desired camera
    psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe

    char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate
    ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate

    if (type == PPSIM_TYPE_OBJECT) {
        // adjust the seeing by the scale
        float seeing = psMetadataLookupF32(&status, recipe, "SEEING");
        if (isnan(seeing)) {
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "seeing is not defined");
            psFree(fpa);
            return NULL;
        }
        float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE");
        if (isnan(scale)) {
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "pixel scale is not defined");
            psFree(fpa);
            return NULL;
        }
        psMetadataAddF32(recipe, PS_LIST_TAIL, "SEEING", PS_META_REPLACE, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale);

        // if we have been supplied an input image, but no ra & dec, use the input image values
        if (input) {
            float ra = psMetadataLookupF32(&status, recipe, "RA");
            if (isnan(ra)) {
                ra = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.RA");
                psMetadataAddF32(recipe, PS_LIST_TAIL, "RA", PS_META_REPLACE, "ra (radians)", ra);
            }
            float dec = psMetadataLookupF32(&status, recipe, "DEC");
            if (isnan(dec)) {
                dec = psMetadataLookupF64(&status, input->fpa->concepts, "FPA.DEC");
                psMetadataAddF32(recipe, PS_LIST_TAIL, "DEC", PS_META_REPLACE, "dec (radians)", dec);
            }
        }
    }

    if ((type == PPSIM_TYPE_OBJECT) || (type == PPSIM_TYPE_FLAT)) {
        // determine the zeropoint from the filter
        float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT");
        if (isnan(zp)) {
            char *filter = psMetadataLookupStr(&status, recipe, "FILTER");
            float zp = ppSimGetZeroPoint (recipe, filter);
            psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp);
        }
    }

    // For photometry, we operate on the chip-mosaicked image.  we create a copy of the mosaicked
    // image for psphot so we can write out a clean image
    bool doPhotom = psMetadataLookupBool(&status, recipe, "PHOTOM"); // Density of fakes
    if (doPhotom) {

        // XXX at the moment, we can perform photometry on the fake sources and the forced
        // photometry positions, but not one or the other.  Also, we only support photometry in
        // the context of an image previously analysed by psphot.  Add support for:
        // * photometry of a pure fake image (requires peak detection and psf creation)
        // * photometry of forced source positions without fake image (this might work, it just
        // requires not generating any fake features).

        if (!input) {
            psError(PS_ERR_UNKNOWN, false, "input image not found; currently required for photometry");
            return NULL;
        }

        // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?)
        pmFPAfile *fakeImage = pmFPAfileDefineChipMosaic(config, output->fpa, "PPSIM.FAKE.CHIP");
        if (!fakeImage) {
            psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FAKE.CHIP"));
            psFree(fpa);
            return NULL;
        }
        if (fakeImage->type != PM_FPA_FILE_IMAGE) {
            psError(PS_ERR_IO, true, "PPSIM.FAKE.CHIP is not of type IMAGE");
            psFree(fpa);
            return NULL;
        }

        // we need a chip image if we perform photometry (is it necessary to build it if we don't use it?)
        pmFPAfile *forceImage = NULL;
        if (input) {
            forceImage = pmFPAfileDefineChipMosaic(config, input->fpa, "PPSIM.FORCE.CHIP");
            if (!forceImage) {
                psError(PS_ERR_IO, false, _("Unable to generate new file from PPSIM.FORCE.CHIP"));
                psFree(fpa);
                return NULL;
            }
            if (forceImage->type != PM_FPA_FILE_IMAGE) {
                psError(PS_ERR_IO, true, "PPSIM.FORCE.CHIP is not of type IMAGE");
                psFree(fpa);
                return NULL;
            }
        }

        // define associated psphot input/output files
        if (!ppSimPhotomFiles (config, fakeImage, forceImage)) {
            psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files for psphot");
            psFree(fpa);
            return NULL;
        }
    } else {
        // have we supplied a psf model?  this happens in ppSimPhotomFiles if we request a photometry
        // analysis.  however, even if we do not, a psf model may be used to generate the fake
        // sources.
        if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) {
            // tie the psf file to the chipMosaic
            pmFPAfileBindFromArgs(&status, output, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF");
            if (!status) {
                psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD");
                psFree(fpa);
                return NULL;
            }
        }
    }

    // PPSIM.SOURCES carries the constructed, fake sources with their true parameters
    // XXX only invoke this code for OBJECT types of images?
    pmFPAfile *simSources = pmFPAfileDefineOutput (config, output->fpa, "PPSIM.SOURCES");
    if (!simSources) {
        psError(PS_ERR_UNKNOWN, false, "Cannot find a rule for PPSIM.SOURCES");
        return false;
    }
    simSources->save = true;

    // if we have loaded an input image, we derive certain values from the image, if possible
    if (input) {
        // we need to extract certain metadata from the image and populate the recipe.
        // or else we need to set the fpa concepts based on the recipe options...

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

        ppSimArgToRecipeF32(&status, recipe, "EXPTIME", input->fpa->concepts, "FPA.EXPOSURE");
        char *filter = ppSimArgToRecipeStr(&status, recipe, "FILTER", input->fpa->concepts, "FPA.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;
                }
            }
        }
    }

    psFree(fpa);
    psFree(view);

    return output;
}
