# 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)

    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");
    xSize = 1000;
    ySize = 1000;

    // 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);

    psEllipseAxes axes;
    psEllipsePol pol;

    // 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);
    
    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;
}
