/******************************************************************************
 * this file defines the WAUSS source shape model (XXX need a better name!).  Note that these
 * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include'
 * statements of their own.  The models use a psVector to represent the set of parameters, with
 * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
 * may thus vary depending on the specifics of the model.  All models which are used a PSF
 * representations share a few parameters, for which # define names are listed in pmModel.h:

   power-law with fitted linear term
   1 / (1 + Az + Bz^2 + z^3/6)

   * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
   * PM_PAR_I0 1    - central intensity
   * PM_PAR_XPOS 2  - X center of object
   * PM_PAR_YPOS 3  - Y center of object
   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
   * PM_PAR_SXY 6   - X*Y term of elliptical contour
   * PM_PAR_7   7   - amplitude of the linear component (A)
   * PM_PAR_8   8   - amplitude of the quadratic component (B)
   *****************************************************************************/

/***
    XXXX the model in this file needs to be tested more carefully.
    fix up the code to follow conventions in the other model function files.
***/

XXX broken code

# define PM_MODEL_FUNC       pmModelFunc_WAUSS
# define PM_MODEL_FLUX       pmModelFlux_WAUSS
# define PM_MODEL_GUESS      pmModelGuess_WAUSS
# define PM_MODEL_LIMITS     pmModelLimits_WAUSS
# define PM_MODEL_RADIUS     pmModelRadius_WAUSS
# define PM_MODEL_FROM_PSF   pmModelFromPSF_WAUSS
# define PM_MODEL_FIT_STATUS pmModelFitStatus_WAUSS

psF64 PS_MODEL_FUNC (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[2];
    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 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
    psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
    psF32 f = params->data.F32[1]*r + params->data.F32[0];

    if (deriv != NULL) {
        // note difference from gaussian: q = params->data.F32[1]*r
        psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
        deriv->data.F32[0] = +1.0;
        deriv->data.F32[1] = +r;
        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
        deriv->data.F32[3] = q*(2.0*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;
        deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
        deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
        // The values of 100 dampen the swing of params->data.F32[7,8] */
    }
    return(f);
}

bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
{

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

    beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000;
    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
    beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;

    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
    params_min[0][0].data.F32[PM_PAR_I0] = 0;
    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
    params_min[0][0].data.F32[PM_PAR_7] = 0.1;

    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
    params_max[0][0].data.F32[PM_PAR_7] = 10.0;

    return (TRUE);
}

bool PS_MODEL_GUESS  (psModel *model, psSource *source)
{

    psVector *params = model->params;

    params->data.F32[0] = source->moments->Sky;
    params->data.F32[1] = source->peak->counts - source->moments->Sky;
    params->data.F32[2] = source->moments->x;
    params->data.F32[3] = source->moments->y;
    params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
    params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
    params->data.F32[6] = source->moments->Sxy;
    // XXX: What are these?
    // params->data.F32[7] = B2;
    // params->data.F32[8] = B3;
    return(true);
}

// this is probably wrong since it uses the gauss integral 2 pi sigma^2
psF64 PS_MODEL_FLUX (const psVector *params)
{
    psF64 A1   = 1 / PS_SQR(params->data.F32[4]);
    psF64 A2   = 1 / PS_SQR(params->data.F32[5]);
    psF64 A3   = params->data.F32[6];
    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(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 PS_MODEL_RADIUS   (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] <= 0)
        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);
}

bool PS_MODEL_FROM_PSF  (psModel *modelPSF, psModel *modelFLT, const 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 < 9; i++) {
        psPolynomial2D *poly = psf->params->data[i-4];
        out[i] = Polynomial2DEval (poly, out[2], out[3]);
    }
    return(true);
}

bool PM_MODEL_FIT_STATUS (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_I0] > 0);
    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);

    return status;
}

# undef PM_MODEL_FUNC
# undef PM_MODEL_FLUX
# undef PM_MODEL_GUESS
# undef PM_MODEL_LIMITS
# undef PM_MODEL_RADIUS
# undef PM_MODEL_FROM_PSF
# undef PM_MODEL_FIT_STATUS
