#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <math.h>
#include <pslib.h>

#include "pmHDU.h"
#include "pmFPA.h"

#include "pmMoments.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmTrend2D.h"
#include "pmPSF.h"
#include "pmModel.h"
#include "pmModelClass.h"
#include "pmPeaks.h"
#include "pmSource.h"
#include "pmSourceUtils.h"
#include "pmModelUtils.h"

#include "pmReadoutFake.h"

#define MODEL_TYPE "PS_MODEL_RGAUSS"    // Type of model to use
#define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model


// Given an object model, circularise it by setting the axes to be identical
static bool circulariseModel(pmModel *model // Model to circularise
    )
{
    assert(model);

    psF32 *params = model->params->data.F32; // Model parameters
    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes
    // Curiously, the minor axis can be larger than the major axis, so need to check.
    if (axes.major >= axes.minor) {
        axes.minor = axes.major;
    } else {
        axes.major = axes.minor;
    }
    return pmPSF_AxesToModel(params, axes);
}

bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
                              const psVector *xOffset, const psVector *yOffset, const pmPSF *psf,
                              float minFlux, int radius, bool circularise)
{
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_INT_LARGER_THAN(numCols, 0, false);
    PS_ASSERT_INT_LARGER_THAN(numRows, 0, false);
    PS_ASSERT_ARRAY_NON_NULL(sources, false);

    if (xOffset || yOffset) {
        PS_ASSERT_VECTOR_NON_NULL(xOffset, false);
        PS_ASSERT_VECTOR_NON_NULL(yOffset, false);
        PS_ASSERT_VECTORS_SIZE_EQUAL(xOffset, yOffset, false);
        PS_ASSERT_VECTOR_TYPE(xOffset, PS_TYPE_S32, false);
        PS_ASSERT_VECTOR_TYPE_EQUAL(xOffset, yOffset, false);
        if (xOffset->n != sources->n) {
            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                    "Number of offset vectors (%ld) and sources (%ld) doesn't match",
                    xOffset->n, sources->n);
            return false;
        }
    }
    PS_ASSERT_PTR_NON_NULL(psf, false);
    if (radius > 0 && isfinite(minFlux) && minFlux > 0.0) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot define both minimum flux and fixed radius.");
        return false;
    }

    readout->image = psImageRecycle(readout->image, numCols, numRows, PS_TYPE_F32);
    psImageInit(readout->image, 0);

    int numSources = sources->n;          // Number of stars

    pmModel *fakeModel = pmModelFromPSFforXY(psf, (float)numCols / 2.0, (float)numRows / 2.0,
                                             1.0); // Fake model, with central intensity of 1.0

    float flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0

    if (circularise && !circulariseModel(fakeModel)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
        psFree(fakeModel);
        return false;
    }
    psFree(fakeModel);

    for (int i = 0; i < numSources; i++) {
        pmSource *source = sources->data[i]; // Source of interest
        if (!isfinite(source->psfMag)) {
            continue;
        }
        float x, y;                     // Coordinates of source
        if (source->modelPSF) {
            x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
            y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
        } else {
            x = source->peak->xf;
            y = source->peak->yf;
        }

        pmModel *fakeModel = pmModelFromPSFforXY(psf, x, y, powf(10.0, -0.4 * source->psfMag) / flux0);
        if (!fakeModel) {
            psErrorClear();
            continue;
        }
        if (circularise && !circulariseModel(fakeModel)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
            psFree(fakeModel);
            return false;
        }

        psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
                fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS],
                fakeModel->params->data.F32[PM_PAR_I0]);

        pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
        fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
        float fakeRadius = radius > 0 ? radius : fakeModel->modelRadius(fakeModel->params, minFlux); // Radius

        if (xOffset) {
            if (!pmSourceDefinePixels(fakeSource, readout, x + xOffset->data.S32[i],
                                      y + yOffset->data.S32[i], fakeRadius)) {
                psErrorClear();
                continue;
            }
            if (!pmModelAddWithOffset(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0,
                                      - xOffset->data.S32[i], - yOffset->data.S32[i])) {
                psErrorClear();
                continue;
            }
        } else {
            if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
                psErrorClear();
                continue;
            }
            if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) {
                psErrorClear();
                continue;
            }
        }
        psFree(fakeSource);
        psFree(fakeModel);
    }

    return true;
}
