# include "psphotInternal.h"

# define COMPARE_RADIUS(A,B) (radius->data.F32[A] < radius->data.F32[B])
# define SWAP_RADIUS(TYPE,A,B) { \
  float tmp; \
  if (A != B) { \
    tmp = radius->data.F32[A]; \
    radius->data.F32[A] = radius->data.F32[B]; \
    radius->data.F32[B] = tmp; \
    tmp = flux->data.F32[A]; \
    flux->data.F32[A] = flux->data.F32[B]; \
    flux->data.F32[B] = tmp; \
    tmp = weight->data.F32[A]; \
    weight->data.F32[A] = weight->data.F32[B]; \
    weight->data.F32[B] = tmp; \
  } \
}

bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psMaskType maskVal) {

    // allocate pmSourceExtendedParameters, if not already defined
    if (!source->extpars) {
        source->extpars = pmSourceExtendedParsAlloc ();
    }

    if (!source->extpars->profile) {
        source->extpars->profile = pmSourceRadialProfileAlloc ();
    }

    int nPts = source->pixels->numRows * source->pixels->numCols;
    source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    source->extpars->profile->flux   = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    source->extpars->profile->weight = psVectorAllocEmpty (nPts, PS_TYPE_F32);

    psVector *radius = source->extpars->profile->radius;
    psVector *flux   = source->extpars->profile->flux;
    psVector *weight = source->extpars->profile->weight;

    // XXX use the extended source model here for Xo, Yo?
    // XXX define a radius scaled to the elliptical contour?

    int n = 0;
    
    float Xo = 0.0;
    float Yo = 0.0;

    if (source->modelEXT) {
      Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
      Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
    } else {
      Xo = source->peak->xf - source->pixels->col0;
      Yo = source->peak->yf - source->pixels->row0;
    }
    for (int iy = 0; iy < source->pixels->numRows; iy++) {
        for (int ix = 0; ix < source->pixels->numCols; ix++) {
            if (source->maskObj->data.U8[iy][ix]) continue;
            radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ;
            flux->data.F32[n]   = source->pixels->data.F32[iy][ix];
            weight->data.F32[n] = source->weight->data.F32[iy][ix];
            n++;
        }
    }
    radius->n = n;
    weight->n = n;
    flux->n = n;

    // sort the vector set by the radius
    PSSORT (radius->n, COMPARE_RADIUS, SWAP_RADIUS, NONE);

    return true;
}
