# include "psphotInternal.h"

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

  bool status;

  assert (source->extpars);
  assert (source->extpars->profile);
  assert (source->extpars->profile->radius);
  assert (source->extpars->profile->flux);

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

  // flux at which to measure isophotal parameters
  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
  assert (status);

  // first find flux at R0
  int firstAbove = -1;
  int lastBelow = -1;
  for (int i = 0; i < radius->n; i++) {
    if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i;
    if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i;
  }
  // if we don't go out far enough, we have a problem...
  if (lastBelow == radius->n - 1) {
    psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius...");
    // XXX skip object? raise a flag ?
    return false;
  }
  if (firstAbove < 0) {
    psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius");
    // XXX raise a flag ?
    return false;
  }

  // average flux in this range
  float fluxR0 = 0.0;
  int fluxRn = 0;
  for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) {
    fluxR0 += flux->data.F32[i];
    fluxRn ++;
  }
  fluxR0 /= (float)(fluxRn);

  // target flux for petrosian radius
  float fluxRP = fluxR0 * PETROSIAN_RF;

  // find the first bin below the flux level and the last above the level
  // XXX can this be done faster with bisection?
  // XXX do I need to worry about crazy outliers?
  // XXX should i be smoothing or fitting the curve?
  int firstBelow = -1;
  int lastAbove = -1;
  for (int i = 0; i < flux->n; i++) {
    if (flux->data.F32[i] > fluxRP) lastAbove = i;
    if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
  }
  // if we don't go out far enough, we have a problem...
  if (lastAbove == radius->n - 1) {
    psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius...");
    // XXX skip object? raise a flag ?
    return false;
  }
  if (firstBelow < 0) {
    psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius");
    // XXX raise a flag ?
    return false;
  }

  // need to examine pixels in this vicinity
  float fluxFirst = 0;
  float fluxLast = 0;
  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
    if (i <= firstBelow) {
      fluxFirst += flux->data.F32[i];
    }
    if (i <= lastAbove) {
      fluxLast += flux->data.F32[i];
    }
  }
  float fluxRPSum    = 0.5*(fluxLast + fluxFirst);
  float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst);
  // XXX need to use the weight appropriately here...

  float rad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
  float radErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);

  if (!source->extpars->petrosian) {
    source->extpars->petrosian = pmSourcePetrosianValuesAlloc ();
  }

  // these are uncalibrated: instrumental mags and pixel units
  source->extpars->petrosian->mag    = -2.5*log10(fluxRPSum);
  source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum;

  source->extpars->petrosian->rad    = rad;
  source->extpars->petrosian->radErr = radErr;

  psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n",
           source->extpars->petrosian->mag, source->extpars->petrosian->magErr,
           source->extpars->petrosian->rad, source->extpars->petrosian->radErr,
           source->peak->xf, source->peak->yf);

  return true;

}
