# include "psphot.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
  // XXX cache this?
  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");

  // first find flux at R0
  int firstBelow = -1;
  int lastAbove = -1;
  for (int i = 0; i < radius->n; i++) {
    if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i;
    if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i;
  }

  // average flux in this range
  float fluxR0 = 0.0;
  int fluxRn = 0;
  for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); 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?
  firstBelow = -1;
  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;
  }

  // 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;

  return true;

}
