Index: trunk/psLib/src/dataManip/psStats.c
===================================================================
--- trunk/psLib/src/dataManip/psStats.c	(revision 643)
+++ trunk/psLib/src/dataManip/psStats.c	(revision 652)
@@ -3,6 +3,6 @@
 #include <string.h>
 #include <stdarg.h>
-#include "pslib.h"
-/******************************************************************************
+/******************************************************************************
+//#include "pslib.h"
 //#include "psMemory.h"
 //#include "psTrace.h"
@@ -64,10 +64,17 @@
     psHistogram *newHist = NULL;
     int numBins = 0;
+    psType dummyFloatType;
+    psType dummyIntType;
+
+    dummyFloatType.type = PS_TYPE_FLOAT;
+    dummyFloatType.dimen = PS_DIMEN_VECTOR;
+    dummyIntType.type = PS_TYPE_FLOAT;
+    dummyIntType.dimen = PS_DIMEN_VECTOR;
 
     numBins = 1 + ((upper - lower) / size);
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
-    newHist->lower = psVectorAlloc(numBins, PS_TYPE_FLOAT);
-    newHist->upper = psVectorAlloc(numBins, PS_TYPE_FLOAT);
-    newHist->nums =  psVectorAlloc(numBins, PS_TYPE_INT);
+    newHist->lower = psVectorAlloc(dummyFloatType, numBins);
+    newHist->upper = psVectorAlloc(dummyFloatType, numBins);
+    newHist->nums =  psVectorAlloc(dummyIntType, numBins);
     newHist->minVal = lower;
     newHist->maxVal = upper;
@@ -86,4 +93,11 @@
     int numBins = 0;
     int i;
+    psType dummyFloatType;
+    psType dummyIntType;
+
+    dummyFloatType.type = PS_TYPE_FLOAT;
+    dummyFloatType.dimen = PS_DIMEN_VECTOR;
+    dummyIntType.type = PS_TYPE_FLOAT;
+    dummyIntType.dimen = PS_DIMEN_VECTOR;
 
     if (lower->n != upper->n) {
@@ -92,11 +106,11 @@
     numBins = lower->n;
     newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
-    newHist->lower = psVectorAlloc(numBins, PS_TYPE_FLOAT);
-    newHist->upper = psVectorAlloc(numBins, PS_TYPE_FLOAT);
+    newHist->lower = psVectorAlloc(dummyFloatType, numBins);
+    newHist->upper = psVectorAlloc(dummyFloatType, numBins);
     for (i=0;i<numBins;i++) {
-        newHist->lower->vec.vecF[i] = lower->vec.vecF[i];
-        newHist->upper->vec.vecF[i] = upper->vec.vecF[i];
-    }
-    newHist->nums = psVectorAlloc(numBins, PS_TYPE_INT);
+        newHist->lower->vec.f[i] = lower->vec.f[i];
+        newHist->upper->vec.f[i] = upper->vec.f[i];
+    }
+    newHist->nums = psVectorAlloc(dummyIntType, numBins);
     newHist->minVal = minVal;
     newHist->maxVal = maxVal;
@@ -135,13 +149,13 @@
     binSize = (myHist->maxVal - myHist->minVal) / (float) myHist->nums->n;
     for (i=0;i<myVector->n;i++) {
-        if (myVector->vec.vecF[i] < myHist->minVal) {
+        if (myVector->vec.f[i] < myHist->minVal) {
             myHist->minNum++;
-        } else if (myVector->vec.vecF[i] > myHist->maxVal) {
+        } else if (myVector->vec.f[i] > myHist->maxVal) {
             myHist->maxNum++;
         } else {
-            binNum = (int) ((myVector->vec.vecF[i] - myHist->minVal) / binSize);
-            if ((myVector->vec.vecF[i] >= myHist->lower->vec.vecF[i]) &&
-                    (myVector->vec.vecF[i] <= myHist->upper->vec.vecF[i])) {
-                myHist->nums->vec.vecI[i]++;
+            binNum = (int) ((myVector->vec.f[i] - myHist->minVal) / binSize);
+            if ((myVector->vec.f[i] >= myHist->lower->vec.f[i]) &&
+                    (myVector->vec.f[i] <= myHist->upper->vec.f[i])) {
+                myHist->nums->vec.i[i]++;
             } else {
                 psError(__func__, "data value was not within the bounds of the bin it should habe been in.");
@@ -164,11 +178,11 @@
     if (maskVector != NULL) {
         for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.vecI[i])) {
-                mean+= maskVector->vec.vecI[i];
+            if (!(maskVal & maskVector->vec.i[i])) {
+                mean+= maskVector->vec.i[i];
             }
         }
     } else {
         for (i=0;i<myVector->n;i++) {
-            mean+= maskVector->vec.vecI[i];
+            mean+= maskVector->vec.i[i];
         }
     }
@@ -185,7 +199,7 @@
     if (maskVector != NULL) {
         for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.vecI[i])) {
-                if (myVector->vec.vecF[i] > max) {
-                    max = maskVector->vec.vecI[i];
+            if (!(maskVal & maskVector->vec.i[i])) {
+                if (myVector->vec.f[i] > max) {
+                    max = maskVector->vec.i[i];
                 }
             }
@@ -193,6 +207,6 @@
     } else {
         for (i=0;i<myVector->n;i++) {
-            if (myVector->vec.vecF[i] > max) {
-                max = maskVector->vec.vecI[i];
+            if (myVector->vec.f[i] > max) {
+                max = maskVector->vec.i[i];
             }
         }
@@ -210,7 +224,7 @@
     if (maskVector != NULL) {
         for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.vecI[i])) {
-                if (myVector->vec.vecF[i] < min) {
-                    min = maskVector->vec.vecI[i];
+            if (!(maskVal & maskVector->vec.i[i])) {
+                if (myVector->vec.f[i] < min) {
+                    min = maskVector->vec.i[i];
                 }
             }
@@ -218,6 +232,6 @@
     } else {
         for (i=0;i<myVector->n;i++) {
-            if (myVector->vec.vecF[i] < min) {
-                min = maskVector->vec.vecI[i];
+            if (myVector->vec.f[i] < min) {
+                min = maskVector->vec.i[i];
             }
         }
@@ -226,4 +240,7 @@
 }
 
+/******************************************************************************
+ 
+ *****************************************************************************/
 int p_psArrayNValues(const psVector *restrict myVector, // float
                      const psVector *restrict maskVector, // ints
@@ -235,5 +252,5 @@
     if (maskVector != NULL) {
         for (i=0;i<myVector->n;i++) {
-            if (!(maskVal & maskVector->vec.vecI[i])) {
+            if (!(maskVal & maskVector->vec.i[i])) {
                 numData++;
             }
@@ -243,4 +260,72 @@
     }
     return(numData);
+}
+
+
+
+#define MEDIAN_SIZE_THRESHOLD 10000
+float p_psArraySampleMedian(const psVector *restrict myVector,
+                            const psVector *restrict maskVector,
+                            unsigned int maskVal)
+{
+    psVector *unsortedVector = NULL;
+    psVector *sortedVector = NULL;
+    psType dummyFloatType;
+    int dataSize = 0;
+    int count = 0;
+    int i = 0;
+    float median = 0.0;
+
+    dummyFloatType.type = PS_TYPE_FLOAT;
+    dummyFloatType.dimen = PS_DIMEN_VECTOR;
+    dataSize = p_psArrayNValues(myVector, maskVector, maskVal);
+
+    if (dataSize < MEDIAN_SIZE_THRESHOLD) {
+        unsortedVector = psVectorAlloc(dummyFloatType, dataSize);
+        sortedVector   = psVectorAlloc(dummyFloatType, dataSize);
+
+        count = 0;
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->vec.i[i])) {
+                    unsortedVector->vec.f[count++] = maskVector->vec.f[i];
+                }
+            }
+            psSort(sortedVector, unsortedVector);
+        } else {
+            psSort(sortedVector, myVector);
+        }
+
+        if (0 == (dataSize % 2)) {
+            median = 0.5 * (sortedVector->vec.f[(dataSize/2)-1] +
+                            sortedVector->vec.f[dataSize/2]);
+        } else {
+            median = sortedVector->vec.f[dataSize/2];
+        }
+    } else {
+        // BROAD: Calculate the Robust Median
+        // Determine the LQ of the distribution.
+        // Determine the UQ of the distribution.
+        // Histogram the data with bin size (sigma_e = (UQ - LQ) / 1.34) / 10.0.
+        // Smooth the histogram with a Gaussian with sigma_s = sigma_e / 4
+        // Find the bin with the peak value between LQ and UQ (the MODE)
+        // dL = (UQ - LQ) / 8
+        // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
+        // The resulting fit parameters are the robust mean, mean_r, and sigma
+        psError(__func__, "GUS: p_psArraySampleMedian(), large data sample");
+    }
+
+    psVectorFree(unsortedVector);
+    psVectorFree(sortedVector);
+    return(median);
+}
+
+
+int p_psArrayXXX(const psVector *restrict myVector,
+                 const psVector *restrict maskVector,
+                 unsigned int maskVal)
+{
+    printf("ERROR: Don't call me: p_psArrayXXX()\n");
+    exit(1);
 }
 
@@ -265,56 +350,56 @@
     switch (stats->options) {
     case PS_STAT_SAMPLE_MEAN:
-        newStruct->max = p_psArraySampleMean(myVector, maskVector, maskVal);
+        newStruct->sampleMean = p_psArraySampleMean(myVector, maskVector, maskVal);
         break;
     case PS_STAT_SAMPLE_MEDIAN:
-        printf("Code me!\n");
+        newStruct->sampleMedian = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_SAMPLE_STDEV:
-        printf("Code me!\n");
+        newStruct->sampleStdev = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_SAMPLE_UQ:
-        printf("Code me!\n");
+        newStruct->sampleUQ = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_SAMPLE_LQ:
-        printf("Code me!\n");
+        newStruct->sampleLQ = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_MEAN:
-        printf("Code me!\n");
+        newStruct->robustMean = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_MEAN_NVALUES:
-        printf("Code me!\n");
+        newStruct->robustMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_MEDIAN:
-        printf("Code me!\n");
+        newStruct->robustMedian = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_MEDIAN_NVALUES:
-        printf("Code me!\n");
+        newStruct->robustMedianNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_MODE:
-        printf("Code me!\n");
+        newStruct->robustMode = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_MODE_NVALUES:
-        printf("Code me!\n");
+        newStruct->robustModeNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_STDEV:
-        printf("Code me!\n");
+        newStruct->robustStdev = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_UQ:
-        printf("Code me!\n");
+        newStruct->robustUQ = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_ROBUST_LQ:
-        printf("Code me!\n");
+        newStruct->robustLQ = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_CLIPPED_MEAN:
-        printf("Code me!\n");
+        newStruct->clippedMean = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_CLIPPED_MEAN_NVALUES:
-        printf("Code me!\n");
+        newStruct->clippedMeanNvalues = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_CLIPPED_MEAN_NSIGMA:
-        printf("Code me!\n");
+        newStruct->clipSigma = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_CLIPPED_STDEV:
-        printf("Code me!\n");
+        newStruct->clippedStdev = p_psArrayXXX(myVector, maskVector, maskVal);
         break;
     case PS_STAT_MAX: