/** @file  pmModelUtils.c
 *
 *  Functions to manipulate object models
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2008-10-08 21:53:08 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

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

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"

#include "pmTrend2D.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmSpan.h"
#include "pmFootprintSpans.h"
#include "pmFootprint.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmModelFuncs.h"
#include "pmModel.h"
#include "pmModelUtils.h"
#include "pmModelClass.h"
#include "pmSourceMasks.h"
#include "pmSourceExtendedPars.h"
#include "pmSourceDiffStats.h"
#include "pmSource.h"
#include "pmSourceFitModel.h"
#include "pmPSF.h"

#include "pmErrorCodes.h"

/*****************************************************************************
pmModelFromPSF (*modelEXT, *psf):  use the model position parameters to
construct a realization of the PSF model at the object coordinates
 *****************************************************************************/
pmModel *pmModelFromPSF (pmModel *modelEXT, const pmPSF *psf)
{
    PS_ASSERT_PTR_NON_NULL(psf, NULL);
    PS_ASSERT_PTR_NON_NULL(modelEXT, NULL);

    // allocate a new pmModel to hold the PSF version
    pmModel *modelPSF = pmModelAlloc (psf->type);

    // set model parameters for this source based on PSF information
    if (!modelEXT->modelFromPSF (modelPSF, modelEXT, psf)) {
        psTrace ("psModules.objects", 3, "Failed to set model params from PSF");
        psFree(modelPSF);
        return NULL;
    }
    // note that model->residuals is just a reference
    modelPSF->residuals = psf->residuals;

    return (modelPSF);
}

// instantiate a model for the PSF at this location with peak flux
// NOTE: psf and (Xo,Yo) are defined wrt chip coordinates
pmModel *pmModelFromPSFforXY (const pmPSF *psf, float Xo, float Yo, float Io)
{
    PS_ASSERT_PTR_NON_NULL(psf, NULL);

    // allocate a new pmModel to hold the PSF version
    pmModel *modelPSF = pmModelAlloc (psf->type);

    // set model parameters for this source based on PSF information
    if (!modelPSF->modelParamsFromPSF (modelPSF, psf, Xo, Yo, Io)) {
        psFree(modelPSF);
        return NULL;
    }

    // note that model->residuals is just a reference
    modelPSF->residuals = psf->residuals;

    return (modelPSF);
}

// set this model to have the requested flux
bool pmModelSetFlux (pmModel *model, float flux) {
    PS_ASSERT_PTR_NON_NULL(model, NULL);
    PS_ASSERT_PTR_NON_NULL(model->params, NULL);

    // set Io to be 1.0
    model->params->data.F32[PM_PAR_I0] = 1.0;

    // determine the normalized flux
    float normFlux = model->modelFlux (model->params);
    assert (isfinite(normFlux));
    assert (normFlux > 0);

    // set the desired normalization
    model->params->data.F32[PM_PAR_I0] = flux / normFlux;

    return true;
}

bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments) {

    psEllipseMoments emoments;
    emoments.x2 = moments->Mxx;
    emoments.xy = moments->Mxy;
    emoments.y2 = moments->Myy;

    // force the axis ratio to be < 20.0
    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);

    if (!isfinite(axes.major)) return false;
    if (!isfinite(axes.minor)) return false;
    if (!isfinite(axes.theta)) return false;

    // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width 
    float scale = (isfinite(moments->Mrf) && (moments->Mrf > 0.0)) ? moments->Mrf / axes.major : 1.0;
    axes.major *= scale;
    axes.minor *= scale;

    psEllipseShape shape = psEllipseAxesToShape (axes);

    if (!isfinite(shape.sx))  return false;
    if (!isfinite(shape.sy))  return false;
    if (!isfinite(shape.sxy)) return false;

    // set the shape parameters
    *Sxx  = PS_MAX(0.5, M_SQRT2*shape.sx);
    *Syy  = PS_MAX(0.5, M_SQRT2*shape.sy);
    *Sxy  = shape.sxy;

    return true;
}

bool pmModelSetNorm (float *Io, pmSource *source) {

    *Io = source->peak->rawFlux;

#ifndef ALLOW_NONFINITE_PEAK
    // Gene says fail of peak !finite
    if (!isfinite(*Io)) return false;
#else 
    // This is the way it used to be. Somtimes an infinite value Io made it's way down the pipeline
    // causing assertion failures
    if (!isfinite(*Io) && !source->moments) return false;

    *Io = source->moments->Peak;
    if (!isfinite(*Io)) return false;
#endif

    return true;
}

bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source) {

    bool useMoments = pmSourcePositionUseMoments(source);
    
    if (useMoments) {
	*Xo = source->moments->Mx;
	*Yo = source->moments->My;
    } else {
	*Xo = source->peak->xf;
	*Yo = source->peak->yf;
    }

    return true;
}
