
/******************************************************************************
    one component, two slopes:
    1 / (1 + z + z^Npow)
 
    params->data.F32[0] = So;
    params->data.F32[1] = Zo;
    params->data.F32[2] = Xo;
    params->data.F32[3] = Yo;
    params->data.F32[4] = 1 / SigmaX;
    params->data.F32[5] = 1 / SigmaY;
    params->data.F32[6] = Sxy;
    params->data.F32[7] = Npow
*****************************************************************************/

psF64 psModelFunc_RGAUSS(psVector *deriv,
                         const psVector *params,
                         const psVector *x)
{
    psF32 *PAR = params->data.F32;

    psF32 X  = x->data.F32[0] - PAR[2];
    psF32 Y  = x->data.F32[1] - PAR[3];
    psF32 px = PAR[4]*X;
    psF32 py = PAR[5]*Y;
    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y;

    // psF32 FACT = 1 + 5*exp(-5*PAR[7]);

    psF32 p  = pow(z, PAR[7] - 1.0);
    psF32 r  = 1.0 / (1 + z + z*p);
    psF32 f  = PAR[1]*r + PAR[0];

    if (deriv != NULL) {
        // note difference from a pure gaussian: q = params->data.F32[1]*r
        psF32 t = PAR[1]*r*r;
        psF32 q = t*(1 + PAR[7]*p);

        deriv->data.F32[0] = +1.0;
        deriv->data.F32[1] = +r;
        deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
        deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
        deriv->data.F32[4] = -q*px*X*2;
        deriv->data.F32[5] = -q*py*Y*2;
        deriv->data.F32[6] = -q*X*Y;
        deriv->data.F32[7] = -5.0*t*log(z)*p*z;

        // deriv->data.F32[4] = -1.8*q*px*X*2;
        // deriv->data.F32[5] = -1.8*q*py*Y*2;
        // deriv->data.F32[6] = -1.5*q*X*Y;
        // deriv->data.F32[7] = -5.0*t*log(z)*p*z;
    }
    return(f);
}

psF64 psModelFlux_RGAUSS(const psVector *params)
{
    float f, norm, z;

    psF32 *PAR = params->data.F32;

    psF64 A1   = PS_SQR(PAR[4]);
    psF64 A2   = PS_SQR(PAR[5]);
    psF64 A3   = PS_SQR(PAR[6]);
    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
    // Area is equivalent to 2 pi sigma^2

    // the area needs to be multiplied by the integral of f(z)
    norm = 0.0;
    for (z = 0.005; z < 50; z += 0.01) {
        f = 1.0 / (1 + z + pow(z, PAR[7]));
        norm += f;
    }
    norm *= 0.01;

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

    return(Flux);
}

// define this function so it never returns Inf or NaN
// return the radius which yields the requested flux
psF64 psModelRadius_RGAUSS  (const psVector *params, psF64 flux)
{
    psF64 z, f, p;
    psF32 *PAR = params->data.F32;

    if (flux <= 0)
        return (1.0);
    if (PAR[1] <= 0)
        return (1.0);
    if (flux >= PAR[1])
        return (1.0);

    // if Sx == Sy, sigma = Sx == Sy
    psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0);
    psF64 dz = 1.0 / (2.0 * sigma*sigma);
    psF64 limit = flux / PAR[1];

    // we can do this much better with intelligent choices here
    for (z = 0.0; z < 20.0; z += dz) {
        p = pow(z, PAR[7]);
        f = 1.0 / (1 + z + p);
        if (f < limit)
            break;
    }
    psF64 radius = sigma * sqrt (2.0 * z);
    if (isnan(radius)) {
        fprintf (stderr, "error in code\n");
    }
    return (radius);
}

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

    psVector *params = model->params;

    EllipseAxes axes;
    EllipseShape shape;
    EllipseMoments moments;

    moments.x2 = PS_SQR(source->moments->Sx);
    moments.y2 = PS_SQR(source->moments->Sy);
    moments.xy = source->moments->Sxy;

    axes = EllipseMomentsToAxes(moments);
    shape = EllipseAxesToShape(axes);

    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] = 1.0 / shape.sx;
    params->data.F32[5] = 1.0 / shape.sy;
    params->data.F32[6] = shape.sxy;
    params->data.F32[7] = 2.0;
    return(true);
}

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