#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 "pmSpan.h"
#include "pmFootprint.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
#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models


// 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 pmReadoutFakeFromVectors(pmReadout *readout, int numCols, int numRows,
                              const psVector *x, const psVector *y, const psVector *mag,
                              const psVector *xOffset, const psVector *yOffset,
                              const pmPSF *psf, float minFlux, int radius,
                              bool circularise, bool normalisePeak)
{
    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_VECTOR_NON_NULL(x, false);
    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
    PS_ASSERT_VECTOR_NON_NULL(y, false);
    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
    PS_ASSERT_VECTOR_NON_NULL(mag, false);
    PS_ASSERT_VECTOR_TYPE(mag, PS_TYPE_F32, false);
    PS_ASSERT_VECTORS_SIZE_EQUAL(mag, x, false);
    long numSources = x->n;              // Number of sources
    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(yOffset, PS_TYPE_S32, false);
        if (xOffset->n != numSources) {
            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                    "Number of offset vectors (%ld) and sources (%ld) doesn't match",
                    xOffset->n, numSources);
            return false;
        }
    }
    PS_ASSERT_PTR_NON_NULL(psf, false);

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

    for (long i = 0; i < numSources; i++) {
        float flux = powf(10.0, -0.4 * mag->data.F32[i]); // Flux of source
        float xSrc = x->data.F32[i], ySrc = y->data.F32[i]; // Coordinates of source

        if (normalisePeak) {
            // Normalise flux
            pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation
            if (!normModel || (normModel->flags & MODEL_MASK)) {
                psFree(normModel);
                continue;
            }
	    // check that all params are valid:
	    bool validParams = true;
	    for (int n = 0; validParams && (n < normModel->params->n); n++) {
		if (n == PM_PAR_SKY) continue;
		if (n == PM_PAR_I0) continue;
		if (n == PM_PAR_XPOS) continue;
		if (n == PM_PAR_YPOS) continue;
		if (!isfinite(normModel->params->data.F32[n])) validParams = false;
	    }
	    if (!validParams) {
                psFree(normModel);
		continue;
	    }		
            if (circularise && !circulariseModel(normModel)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
                psFree(normModel);
                return false;
            }

            flux /= normModel->modelFlux(normModel->params);
            psFree(normModel);
        }

        pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux);
        if (!fakeModel || (fakeModel->flags & MODEL_MASK)) {
            psFree(fakeModel);
            continue;
        }
	// check that all params are valid:
	bool validParams = true;
	for (int n = 0; validParams && (n < fakeModel->params->n); n++) {
	    if (n == PM_PAR_SKY) continue;
	    if (n == PM_PAR_I0) continue;
	    if (n == PM_PAR_XPOS) continue;
	    if (n == PM_PAR_YPOS) continue;
	    if (!isfinite(fakeModel->params->data.F32[n])) validParams = false;
	}
	if (!validParams) {
	    psFree(fakeModel);
	    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(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
        float fakeRadius = 1.0;         // Radius of fake source
        if (isfinite(minFlux)) {
            fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
        }
        if (radius > 0) {
            fakeRadius = PS_MAX(fakeRadius, radius);
        }

        if (xOffset) {
            if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[i],
                                      ySrc + 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, xSrc, ySrc, fakeRadius)) {
                psErrorClear();
                continue;
            }
            if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) {
                psErrorClear();
                continue;
            }
        }
        psFree(fakeSource);
        psFree(fakeModel);
    }

    return true;
}


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

    int numSources = sources->n;          // Number of stars
    psVector *x = psVectorAllocEmpty(numSources, PS_TYPE_F32);
    psVector *y = psVectorAllocEmpty(numSources, PS_TYPE_F32);
    psVector *mag = psVectorAllocEmpty(numSources, PS_TYPE_F32);

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

        x->data.F32[numGood] = xSrc;
        y->data.F32[numGood] = ySrc;
        mag->data.F32[numGood] = source->psfMag;
        numGood++;
    }
    x->n = numGood;
    y->n = numGood;
    mag->n = numGood;

    bool status = pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf,
                                           minFlux, radius, circularise, normalisePeak);
    psFree(x);
    psFree(y);
    psFree(mag);

    return status;
}
