# include "psphotInternal.h"

# define PETROSIAN_RATIO 0.2
# define PETROSIAN_RADII 2.0

// generate the Petrosian radius and flux from the mean surface brightness (r_i)

float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);

bool psphotPetrosianStats (pmSource *source) {

    psAssert (source, "missing source");
    psAssert (source->extpars, "missing extpars");
    psAssert (source->extpars->petProfile, "missing petProfile");

    pmSourceRadialProfile *profile = source->extpars->petProfile;

    psVector *binSB      = profile->binSB;
    psVector *binSBstdev = profile->binSBstdev;
    psVector *binRad     = profile->radialBins;
    psVector *area       = profile->area;

    psVector *fluxSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    psVector *fluxSumErr2 = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    psVector *refRadius   = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    psVector *petRatio    = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    psVector *petRatioErr = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    psVector *meanSB      = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);

    float petRadius = NAN;
    float petFlux = NAN;

    bool anyPetro = false;
    bool manyPetro = false;
    bool above = true;
    float Asum = 0.0;
    float Fsum = 0.0;
    float dFsum2 = 0.0;

    float nSigma = 3.0;
    int lowestSignificantRadius = 0;
    float lowestSignificantRatio = 1.0;

    // find the Petrosian Radius and Petrosian Flux

    int nOut = 0;
    for (int i = 0; i < binSB->n; i++) {
	// skip nan bins (do not contribute to flux or area)
	if (!isfinite(binSB->data.F32[i])) continue;

	float Area = area->data.F32[i];
	Asum += Area;
	Fsum += binSB->data.F32[i] * Area;
	dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area);

	float areaInner = 0.5 * Area;
	float fluxInner = 0.5 * Area * binSB->data.F32[i];
	float fluxInnerErr2 = PS_SQR(binSBstdev->data.F32[i] * 0.5 * Area);
	if (nOut > 0) {
	    areaInner += areaSum->data.F32[nOut-1];
	    fluxInner += fluxSum->data.F32[nOut-1];
	    fluxInnerErr2 += fluxSumErr2->data.F32[nOut-1];
	}

	// ratio = binSB / meanSB
	// meanSB = flux / area
	// flux = sum(binSB(i) * area(i)
	// fluxErr^2 = sum(binSBerr(i)^2 area(i)^2)
	// meanSBerr^2 = fluxErr^2 / area^2
	// (ratioErr/ratio)^2 = (binSBerr/binSB)^2 + (meanSBerr/meanSB)^2

	psVectorAppend(meanSB, (fluxInner / areaInner));

	float ratio = binSB->data.F32[i] / meanSB->data.F32[nOut];
	psVectorAppend(petRatio, ratio);

	float meanSBerr = sqrt(fluxInnerErr2) / areaInner;
	float ratioErr = fabs(ratio) * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut]));

	psVectorAppend(petRatioErr, ratioErr);

	psVectorAppend(areaSum, Asum);
	psVectorAppend(fluxSum, Fsum);
	psVectorAppend(fluxSumErr2, dFsum2);
	psVectorAppend(refRadius, binRad->data.F32[i]);

	psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", 
		 i, refRadius->data.F32[nOut], 
		 binSB->data.F32[i], binSBstdev->data.F32[i], 
		 meanSB->data.F32[nOut], meanSBerr, 
		 petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], 
		 fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
    
	// anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
	// we will keep and report the last (largest radius) value
	if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) {
	    // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
	    if (i == 0) { 
		// assume Fmax @ R = 0.0
		petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
	    } else {
		petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
	    }
	    above = false;
	    if (anyPetro) manyPetro = true;
	    anyPetro = true;
	}
    
	// anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
	// we will keep and report the last (largest radius) value
	// find the last signficant measurement of the petrosian ratio
	if (above && (petRatio->data.F32[nOut] < lowestSignificantRatio) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) {
	    lowestSignificantRadius = nOut;
	    lowestSignificantRatio = petRatio->data.F32[nOut];
	}
    
	// reset on transitions up, but do not re-calculate rad_90, flux_90
	if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) {
	    above = true;
	}
	nOut ++;
    }

    if (!anyPetro) {
	// interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
	if (lowestSignificantRadius == 0) { 
	    // assume Fmax @ R = 0.0
	    petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
	} else {
	    petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
	}
    }

    // now measure the flux within PETROSIAN_RADII * petRadius 
    float apRadius = PETROSIAN_RADII * petRadius;
    for (int i = 0; i < refRadius->n; i++) {
	// XXX use bisection to do this faster:
	if (refRadius->data.F32[i] > apRadius) {
	    if (i == 0) {
		psWarning ("does this case make any sense? (refRadius[0] > apRadius)");
		continue;
	    } else {
		petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
		break;
	    }
	}
    }

    // now measure the radii R90 and R50 where flux = 0.9 (or 0.5) * petFlux;
    float flux90 = 0.9 * petFlux;
    float flux50 = 0.5 * petFlux;
    float R50 = NAN;
    float R90 = NAN;
    bool found50 = false;
    bool found90 = false;
    // XXX use bisection to do this faster:
    for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) {
	if (!found50 && (fluxSum->data.F32[i] > flux50)) {
	    if (i == 0) {
		psWarning ("does this case make any sense? (fluxSum[0] > flux50)");
		continue;
	    } else {
		R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
		found50 = true;
	    }
	}
	if (!found90 && (fluxSum->data.F32[i] > flux90)) {
	    if (i == 0) {
		psWarning ("does this case make any sense? (fluxSum[0] > flux90)");
		continue;
	    } else {
		R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
		found90 = true;
	    }
	}
    }


    // XXX save flags (anyPetro, manyPetro)
    source->extpars->petrosianRadius = petRadius;
    source->extpars->petrosianFlux   = petFlux;
    source->extpars->petrosianR50    = R50;
    source->extpars->petrosianR90    = R90;
    
    // XXX add the errors
    source->extpars->petrosianRadiusErr = NAN;
    source->extpars->petrosianFluxErr   = NAN;
    source->extpars->petrosianR50Err    = NAN;
    source->extpars->petrosianR90Err    = NAN;

    fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
    psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);

    psFree(fluxSum);
    psFree(fluxSumErr2);
    psFree(refRadius);
    psFree(petRatio);
    psFree(petRatioErr);
    psFree(meanSB);
    psFree(areaSum);

    return true;
}

float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
    float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0);
    return Y;
}

