# include "psphot.h"

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

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

  bool status;

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

  // flux at which to measure isophotal parameters
  // XXX cache this?
  float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");

  // find the first bin below the flux level and the last above the level
  // XXX can this be done faster with disection?
  // 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] > ISOPHOT_FLUX) lastAbove = i;
    if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i;
  }

  // need to examine pixels in this vicinity
  float isophotalFluxFirst = 0;
  float isophotalFluxLast = 0;
  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
    if (i <= firstBelow) {
      isophotalFluxFirst += flux->data.F32[i];
    }
    if (i <= lastAbove) {
      isophotalFluxLast += flux->data.F32[i];
    }
  }
  float isophotalFlux    = 0.5*(isophotalFluxLast + isophotalFluxFirst);
  float isophotalFluxErr = 0.5*fabs(isophotalFluxLast - isophotalFluxFirst);

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

  if (!source->extpars->isophot) {
    source->extpars->isophot = pmSourceIsophotalValuesAlloc ();
  }
  
  // these are uncalibrated: instrumental mags and pixel units
  source->extpars->isophot->mag    = -2.5*log10(isophotalFlux);
  source->extpars->isophot->magErr = isophotalFluxErr / isophotalFlux;

  source->extpars->isophot->rad    = isophotalRad;
  source->extpars->isophot->radErr = isophotalRadErr;

  return true;

}
