Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c	(revision 25203)
+++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c	(revision 25204)
@@ -18,4 +18,5 @@
 
     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);
@@ -27,6 +28,7 @@
     bool manyPetro = false;
     bool above = true;
+    float Asum = 0.0;
     float Fsum = 0.0;
-    float Asum = 0.0;
+    float dFsum2 = 0.0;
 
     int nOut = 0;
@@ -36,19 +38,37 @@
 
 	float Area = area->data.F32[i];
+	Asum += Area;
 	Fsum += binSB->data.F32[i] * Area;
-	Asum += Area;
+	dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area);
 
+	float areaInner = 0.5 * Area;
 	float fluxInner = 0.5 * Area * binSB->data.F32[i];
-	float areaInner = 0.5 * Area;
+	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];
-	    areaInner += areaSum->data.F32[nOut-1];
+	    fluxInnerErr2 = PS_SQR(fluxSumErr2->data.F32[nOut-1] * areaSum->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));
-	psVectorAppend(petRatio, binSB->data.F32[i] / meanSB->data.F32[nOut]);
-	psVectorAppend(petRatioErr, binSBstdev->data.F32[i] / meanSB->data.F32[nOut]);
+
+	float ratio = binSB->data.F32[i] / meanSB->data.F32[nOut];
+	psVectorAppend(petRatio, ratio);
+
+	float meanSBerr = sqrt(fluxInnerErr2) / areaInner;
+	float ratioErr = 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(areaSum, Asum);
+	psVectorAppend(fluxSumErr2, dFsum2);
 	psVectorAppend(refRadius, binRad->data.F32[i]);