/******************************************************************************
    params->data.F32[0] = So;
    params->data.F32[1] = Zo;
    params->data.F32[2] = Xo;
    params->data.F32[3] = Yo;
    params->data.F32[4] = sqrt(2.0) / SigmaX;
    params->data.F32[5] = sqrt(2.0) / SigmaY;
    params->data.F32[6] = Sxy;
*****************************************************************************/

psF32 pmModelFunc_TAUSS(psVector *deriv,
                        const psVector *params,
                        const psVector *x)
{
    psF32 X  = x->data.F32[0] - params->data.F32[2];
    psF32 Y  = x->data.F32[1] - params->data.F32[3];
    psF32 px = params->data.F32[4]*X;
    psF32 py = params->data.F32[5]*Y;
    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
    psF32 r  = exp(-z);
    psF32 q  = params->data.F32[1]*r;
    psF32 f  = q + params->data.F32[0];

    if (deriv != NULL) {
        deriv->data.F32[0] = +1.0;
        deriv->data.F32[1] = +r;
        deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
        deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
        deriv->data.F32[4] = -2.0*q*px*X;
        deriv->data.F32[5] = -2.0*q*py*Y;
        deriv->data.F32[6] = -q*X*Y;
    }
    return(f);
}

psF64 pmModelFlux_TAUSS(const psVector *params)
{
    psF64 A1   = PS_SQR(params->data.F32[4]);
    psF64 A2   = PS_SQR(params->data.F32[5]);
    psF64 A3   = PS_SQR(params->data.F32[6]);
    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
    // Area is equivalent to 2 pi sigma^2

    psF64 Flux = params->data.F32[1] * Area;

    return(Flux);
}

// return the radius which yields the requested flux
psF64 pmModelRadius_TAUSS  (const psVector *params, psF64 flux)
{

    if (flux <= 0)
        return (1.0);
    if (params->data.F32[1] <= 0)
        return (1.0);
    if (flux >= params->data.F32[1])
        return (1.0);

    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
    return (radius);
}

// define the parameter limits
bool pmModelLimits_TAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max)
{

    *beta_lim   = psVectorAlloc (7, PS_TYPE_F32);
    *params_min = psVectorAlloc (7, PS_TYPE_F32);
    *params_max = psVectorAlloc (7, PS_TYPE_F32);

    beta_lim[0][0].data.F32[0] = 1000;
    beta_lim[0][0].data.F32[1] = 10000;
    beta_lim[0][0].data.F32[2] = 5;
    beta_lim[0][0].data.F32[3] = 5;
    beta_lim[0][0].data.F32[4] = 0.5;
    beta_lim[0][0].data.F32[5] = 0.5;
    beta_lim[0][0].data.F32[6] = 0.5;

    params_min[0][0].data.F32[0] = -1000;
    params_min[0][0].data.F32[1] = 0;
    params_min[0][0].data.F32[2] = -100;
    params_min[0][0].data.F32[3] = -100;
    params_min[0][0].data.F32[4] = 0.01;
    params_min[0][0].data.F32[5] = 0.01;
    params_min[0][0].data.F32[6] = -5.0;

    params_max[0][0].data.F32[0] = 1e5;
    params_max[0][0].data.F32[1] = 1e6;
    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
    params_max[0][0].data.F32[4] = 2.0;
    params_max[0][0].data.F32[5] = 2.0;
    params_max[0][0].data.F32[6] = +5.0;

    return (TRUE);
}

// make an initial guess for parameters
bool pmModelGuess_TAUSS (pmModel *model, pmSource *source)
{

    pmMoments *moments = source->moments;
    psF32     *params  = model->params->data.F32;

    params[0] = moments->Sky;
    params[1] = moments->Peak - moments->Sky;
    params[2] = moments->x;
    params[3] = moments->y;
    params[4] = 1.2 / moments->Sx;
    params[5] = 1.2 / moments->Sy;
    params[6] = 0.0;
    return(true);
}

// construct the PSF model from the FLT model and the psf
bool pmModelFromPSF_TAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
{

    psF32 *out = modelPSF->params->data.F32;
    psF32 *in  = modelFLT->params->data.F32;

    out[0] = in[0];
    out[1] = in[1];
    out[2] = in[2];
    out[3] = in[3];

    for (int i = 4; i < 7; i++) {
        psPolynomial2D *poly = psf->params->data[i-4];
        // XXX: be wary of a bug here.  EAM had his own version of psPolynomial2DEval().
        // I think the reason was that hewas using nX=xOrder.  Not sure.  I changed this
        // without verifying.
        // out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
        out[i] = psPolynomial2DEval(poly, out[2], out[3]);
    }
    return(true);
}

// check the status of the fitted model
bool pmModelFitStatus_TAUSS (pmModel *model)
{

    psF32 dP;
    bool  status;

    psF32 *PAR  = model->params->data.F32;
    psF32 *dPAR = model->dparams->data.F32;

    dP = 0;
    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
    dP = sqrt (dP);

    status = true;
    status &= (dP < 0.5);
    status &= (PAR[PM_PAR_FLUX] > 0);
    status &= ((dPAR[PM_PAR_FLUX]/PAR[PM_PAR_FLUX]) < 0.5);

    if (status)
        return true;
    return false;
}
