Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 679)
+++ trunk/psLib/src/math/psStats.c	(revision 682)
@@ -3,5 +3,4 @@
 #include <string.h>
 #include <stdarg.h>
-/******************************************************************************
 //#include "pslib.h"
 #include "psMemory.h"
@@ -11,6 +10,9 @@
 //#include "psError.h"
 //#include "psArray.h"
- *****************************************************************************/
 #include "psStats.h"
+
+#include "float.h"
+#include <math.h>
+
 /******************************************************************************
     psStatsAlloc(): This routine must create a new psStats data structure.
@@ -19,26 +21,27 @@
 {
     psStats *newStruct = NULL;
+
     newStruct = (psStats *) psAlloc(sizeof(psStats));
-    newStruct->sampleMean = 0.0;
-    newStruct->sampleMedian = 0.0;
-    newStruct->sampleStdev = 0.0;
-    newStruct->sampleUQ = 0.0;
-    newStruct->sampleLQ = 0.0;
-    newStruct->robustMean = 0.0;
+    newStruct->sampleMean = NAN;
+    newStruct->sampleMedian = NAN;
+    newStruct->sampleStdev = NAN;
+    newStruct->sampleUQ = NAN;
+    newStruct->sampleLQ = NAN;
+    newStruct->robustMean = NAN;
     newStruct->robustMeanNvalues = 0;
-    newStruct->robustMedian = 0.0;
+    newStruct->robustMedian = NAN;
     newStruct->robustMedianNvalues = 0;
-    newStruct->robustMode = 0.0;
+    newStruct->robustMode = NAN;
     newStruct->robustModeNvalues = 0;
-    newStruct->robustStdev = 0.0;
-    newStruct->robustUQ = 0.0;
-    newStruct->robustLQ = 0.0;
-    newStruct->clippedMean = 0.0;
+    newStruct->robustStdev = NAN;
+    newStruct->robustUQ = NAN;
+    newStruct->robustLQ = NAN;
+    newStruct->clippedMean = NAN;
     newStruct->clippedMeanNvalues = 0;
-    newStruct->clippedStdev = 0.0;
-    newStruct->clipSigma = 0.0;
+    newStruct->clippedStdev = NAN;
+    newStruct->clipSigma = NAN;
     newStruct->clipIter = 0;
-    newStruct->min = 0.0;
-    newStruct->max = 0.0;
+    newStruct->min = NAN;
+    newStruct->max = NAN;
     newStruct->nValues = 0;
     newStruct->options = options;
@@ -156,7 +159,8 @@
     MISC STATISTICAL FUNCTIONS
  *****************************************************************************/
-float p_psArraySampleMean(const psVector *restrict myVector,
-                          const psVector *restrict maskVector,
-                          unsigned int maskVal)
+void p_psArraySampleMean(const psVector *restrict myVector,
+                         const psVector *restrict maskVector,
+                         unsigned int maskVal,
+                         psStats *newStruct)
 {
     int i = 0;
@@ -178,15 +182,15 @@
         mean/= (float) myVector->n;
     }
-    return(mean);
-}
-
-float p_psArrayMax(const psVector *restrict myVector,
-                   const psVector *restrict maskVector,
-                   unsigned int maskVal)
+    newStruct->sampleMean = mean;
+}
+
+void p_psArrayMax(const psVector *restrict myVector,
+                  const psVector *restrict maskVector,
+                  unsigned int maskVal,
+                  psStats *newStruct)
 {
     int i = 0;
     float max = -1e99;
 
-    printf("COOL: p_psArrayMax()\n");
     if (maskVector != NULL) {
         for (i=0;i<myVector->n;i++) {
@@ -204,10 +208,11 @@
         }
     }
-    return(max);
-}
-
-float p_psArrayMin(const psVector *restrict myVector, /* FLOATS */
-                   const psVector *restrict maskVector, /* INTS */
-                   unsigned int maskVal)
+    newStruct->max = max;
+}
+
+void p_psArrayMin(const psVector *restrict myVector, /* FLOATS */
+                  const psVector *restrict maskVector, /* INTS */
+                  unsigned int maskVal,
+                  psStats *newStruct)
 {
     int i = 0;
@@ -229,5 +234,5 @@
         }
     }
-    return(min);
+    newStruct->min = min;
 }
 
@@ -235,7 +240,8 @@
  
  *****************************************************************************/
-int p_psArrayNValues(const psVector *restrict myVector, // float
-                     const psVector *restrict maskVector, // ints
-                     unsigned int maskVal)
+void p_psArrayNValues(const psVector *restrict myVector,
+                      const psVector *restrict maskVector,
+                      unsigned int maskVal,
+                      psStats *newStruct)
 {
     int i = 0;
@@ -251,5 +257,5 @@
         numData = myVector->n;
     }
-    return(numData);
+    newStruct->nValues = numData;
 }
 
@@ -257,7 +263,8 @@
 
 #define MEDIAN_SIZE_THRESHOLD 10000
-float p_psArraySampleMedian(const psVector *restrict myVector,
-                            const psVector *restrict maskVector,
-                            unsigned int maskVal)
+void p_psArraySampleMedian(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *newStruct)
 {
     psVector *unsortedVector = NULL;
@@ -268,5 +275,5 @@
     float median = 0.0;
 
-    dataSize = p_psArrayNValues(myVector, maskVector, maskVal);
+    p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
 
     if (dataSize < MEDIAN_SIZE_THRESHOLD) {
@@ -307,14 +314,56 @@
     psVectorFree(unsortedVector);
     psVectorFree(sortedVector);
-    return(median);
-}
-
-
-int p_psArrayXXX(const psVector *restrict myVector,
-                 const psVector *restrict maskVector,
-                 unsigned int maskVal)
+    newStruct->sampleMedian = median;
+}
+
+
+float p_psArrayXXX(const psVector *restrict myVector,
+                   const psVector *restrict maskVector,
+                   unsigned int maskVal)
 {
     printf("ERROR: Don't call me: p_psArrayXXX()\n");
     exit(1);
+}
+
+void p_psArraySampleStdev(const psVector *restrict myVector,
+                          const psVector *restrict maskVector,
+                          unsigned int maskVal,
+                          psStats *newStruct)
+{
+    int i            = 0;
+    int countInt     = 0;
+    float countFloat = 0.0;
+    float       mean = 0.0;
+    float       diff = 0.0;
+    float sumSquares = 0.0;
+    float sumDiffs = 0.0;
+
+    if (0 != isnan(newStruct->sampleMean)) {
+        p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
+    }
+    mean = newStruct->sampleMean;
+
+    if (maskVector != NULL) {
+        for (i=0;i<myVector->n;i++) {
+            if (!(maskVal & maskVector->vec.ui8[i])) {
+                diff = myVector->vec.f[i] - mean;
+                sumSquares+= (diff * diff);
+                sumDiffs+= diff;
+                countInt++;
+
+            }
+        }
+    } else {
+        for (i=0;i<myVector->n;i++) {
+            diff = myVector->vec.f[i] - mean;
+            sumSquares+= (diff * diff);
+            sumDiffs+= diff;
+            countInt++;
+        }
+        countInt = myVector->n;
+    }
+    countFloat = (float) countInt;
+    newStruct->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
+                                     sumDiffs/countFloat))/ (countFloat-1));
 }
 
@@ -333,5 +382,4 @@
     psStats *newStruct = NULL;
 
-    newStruct = psStatsAlloc(stats->options);
     if (myVector == NULL) {
         psAbort(__func__,
@@ -354,29 +402,26 @@
         }
     }
+    newStruct = psStatsAlloc(stats->options);
 
     if (stats->options & PS_STAT_SAMPLE_MEAN) {
-        newStruct->sampleMean = p_psArraySampleMean(myVector,
-                                maskVector,
-                                maskVal);
+        p_psArraySampleMean(myVector, maskVector, maskVal, newStruct);
     }
 
     if (stats->options & PS_STAT_MAX) {
-        newStruct->max = p_psArrayMax(myVector, maskVector, maskVal);
+        p_psArrayMax(myVector, maskVector, maskVal, newStruct);
     }
 
     if (stats->options & PS_STAT_MIN) {
-        newStruct->min = p_psArrayMin(myVector, maskVector, maskVal);
+        p_psArrayMin(myVector, maskVector, maskVal, newStruct);
     }
 
     if (stats->options & PS_STAT_NVALUES) {
-        newStruct->nValues = p_psArrayNValues(myVector, maskVector, maskVal);
-    }
-
-    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
-        newStruct->sampleMedian = p_psArrayXXX(myVector, maskVector, maskVal);
-    }
+        p_psArrayNValues(myVector, maskVector, maskVal, newStruct);
+    }
+
+    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {}
 
     if (stats->options & PS_STAT_SAMPLE_STDEV) {
-        newStruct->sampleStdev = p_psArrayXXX(myVector, maskVector, maskVal);
+        p_psArraySampleStdev(myVector, maskVector, maskVal, newStruct);
     }