Index: /trunk/psLib/src/dataManip/Makefile
===================================================================
--- /trunk/psLib/src/dataManip/Makefile	(revision 562)
+++ /trunk/psLib/src/dataManip/Makefile	(revision 563)
@@ -3,5 +3,5 @@
 endif
 
-TARGET = libpsDataManip.a
+TARGET = psStats.o libpsDataManip.a
 
 all: $(TARGET)
@@ -13,5 +13,5 @@
 %.o: %.c
 	@echo "    Compiling $<. "
-	$(CC) $(CFLAGS) -I$(includedir) -c $< -o $@
+	$(CC) $(CFLAGS) -I../../include -c $< -o $@
 
 libpsDataManip.a: $(SRC_OBJS)
Index: /trunk/psLib/src/dataManip/psStats.c
===================================================================
--- /trunk/psLib/src/dataManip/psStats.c	(revision 563)
+++ /trunk/psLib/src/dataManip/psStats.c	(revision 563)
@@ -0,0 +1,333 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdarg.h>
+#include "pslib.h"
+//#include "psMemory.h"
+//#include "psTrace.h"
+//#include "psString.h"
+//#include "psError.h"
+#include "psArray.h"
+#include "psStats.h"
+
+/******************************************************************************
+    psStatsAlloc(): This routine must create a new psStats data structure.
+ *****************************************************************************/
+psStats *psStatsAlloc(psStatsOptions options)
+{
+    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->robustMeanNvalues = 0;
+    newStruct->robustMedian = 0.0;
+    newStruct->robustMedianNvalues = 0;
+    newStruct->robustMode = 0.0;
+    newStruct->robustModeNvalues = 0;
+    newStruct->robustStdev = 0.0;
+    newStruct->robustUQ = 0.0;
+    newStruct->robustLQ = 0.0;
+    newStruct->clippedMean = 0.0;
+    newStruct->clippedMeanNvalues = 0;
+    newStruct->clippedStdev = 0.0;
+    newStruct->clipSigma = 0.0;
+    newStruct->clipIter = 0;
+    newStruct->min = 0.0;
+    newStruct->max = 0.0;
+    newStruct->nValues = 0;
+    newStruct->options = options;
+
+    return(newStruct);
+}
+
+/******************************************************************************
+    psStatsFree(): This routine must free the psStats data structure.
+ *****************************************************************************/
+void psStatsFree(psStats *stats)
+{
+    psFree(stats);
+}
+
+/******************************************************************************
+    HISTOGRAM FUNCTIONS
+ *****************************************************************************/
+psHistogram *psHistogramAlloc(float lower,
+                              float upper,
+                              float size )
+{
+    psHistogram *newHist = NULL;
+    int numBins = 0;
+
+    numBins = 1 + ((upper - lower) / size);
+    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
+    newHist->lower = psFloatArrayAlloc(numBins);
+    newHist->upper = psFloatArrayAlloc(numBins);
+    newHist->nums = psIntArrayAlloc(numBins);
+    newHist->minVal = lower;
+    newHist->maxVal = upper;
+    newHist->minNum = 0;
+    newHist->maxNum = 0;
+
+    return(newHist);
+}
+
+psHistogram *psHistogramAllocGeneric(const psFloatArray *restrict lower,
+                                     const psFloatArray *restrict upper,
+                                     float minVal,
+                                     float maxVal)
+{
+    psHistogram *newHist = NULL;
+    int numBins = 0;
+    int i;
+
+    if (lower->n != upper->n) {
+        psError(__func__, "There is a different number of upper/lower values");
+    }
+    numBins = lower->n;
+    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
+    newHist->lower = psFloatArrayAlloc(numBins);
+    newHist->upper = psFloatArrayAlloc(numBins);
+    for (i=0;i<numBins;i++) {
+        newHist->lower->arr[i] = lower->arr[i];
+        newHist->upper->arr[i] = upper->arr[i];
+    }
+    newHist->nums = psIntArrayAlloc(numBins);
+    newHist->minVal = minVal;
+    newHist->maxVal = maxVal;
+    newHist->minNum = 0;
+    newHist->maxNum = 0;
+
+    return(newHist);
+}
+
+void psHistogramFree(psHistogram *restrict myHist)
+{
+    psFree(myHist);
+}
+
+
+/*****************************************************************************
+    NOTE: Can we assume fixed size bins in the myHist data structure?
+ 
+    NOTE: This procedure assumes that the following has been already set
+        myHist->lower
+        myHist->upper
+ myHist->minVal
+ myHist->maxVal
+    This procedure sets
+        myHist->nums
+        myHist->minNum
+        myHist->maxNum
+ *****************************************************************************/
+psHistogram *psGetArrayHistogram(       psHistogram *restrict myHist,
+                                        const psFloatArray *restrict myArray)
+{
+    int i = 0;
+    float binSize = 0.0;
+    int binNum = 0;
+
+    binSize = (myHist->maxVal - myHist->minVal) / (float) myHist->nums->n;
+    for (i=0;i<myArray->n;i++) {
+        if (myArray->arr[i] < myHist->minVal) {
+            myHist->minNum++;
+        } else if (myArray->arr[i] > myHist->maxVal) {
+            myHist->maxNum++;
+        } else {
+            binNum = (int) ((myArray->arr[i] - myHist->minVal) / binSize);
+            if ((myArray->arr[i] >= myHist->lower->arr[i]) &&
+                    (myArray->arr[i] <= myHist->upper->arr[i])) {
+                myHist->nums->arr[i]++;
+            } else {
+                psError(__func__, "data value was not within the bounds of the bin it should habe been in.");
+            }
+        }
+    }
+    return(myHist);
+}
+
+/******************************************************************************
+    MISC STATISTICAL FUNCTIONS
+ *****************************************************************************/
+float p_psArraySampleMean(const psFloatArray *restrict myArray,
+                          const psIntArray *restrict maskArray,
+                          unsigned int maskVal)
+{
+    int i = 0;
+    float mean = 0.0;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                mean+= maskArray->arr[i];
+            }
+        }
+    } else {
+        for (i=0;i<myArray->n;i++) {
+            mean+= maskArray->arr[i];
+        }
+    }
+    return(mean);
+}
+
+float p_psArrayMax(const psFloatArray *restrict myArray,
+                   const psIntArray *restrict maskArray,
+                   unsigned int maskVal)
+{
+    int i = 0;
+    float max = -1e99;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                if (myArray->arr[i] > max) {
+                    max = maskArray->arr[i];
+                }
+            }
+        }
+    } else {
+        for (i=0;i<myArray->n;i++) {
+            if (myArray->arr[i] > max) {
+                max = maskArray->arr[i];
+            }
+        }
+    }
+    return(max);
+}
+
+float p_psArrayMin(const psFloatArray *restrict myArray,
+                   const psIntArray *restrict maskArray,
+                   unsigned int maskVal)
+{
+    int i = 0;
+    float min = 1e99;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                if (myArray->arr[i] < min) {
+                    min = maskArray->arr[i];
+                }
+            }
+        }
+    } else {
+        for (i=0;i<myArray->n;i++) {
+            if (myArray->arr[i] < min) {
+                min = maskArray->arr[i];
+            }
+        }
+    }
+    return(min);
+}
+
+int p_psArrayNValues(const psFloatArray *restrict myArray,
+                     const psIntArray *restrict maskArray,
+                     unsigned int maskVal)
+{
+    int i = 0;
+    int numData = 0;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                numData++;
+            }
+        }
+    } else {
+        numData = myArray->n;
+    }
+    return(numData);
+}
+
+psStats *psArrayStats(const psFloatArray *restrict myArray,
+                      const psIntArray *restrict maskArray,
+                      unsigned int maskVal,
+                      psStats *stats)
+{
+    psStats *newStruct = NULL;
+
+    newStruct = psStatsAlloc(stats->options);
+    if (myArray == NULL) {
+        psError(__func__,
+                "Input data array (myArray) was NULL.");
+    }
+
+    if ((maskArray != NULL) &&
+            (myArray->n != maskArray->n)) {
+        psError(__func__, "Array data and array mask are of different sizes.");
+    }
+
+    switch (stats->options) {
+    case PS_STAT_SAMPLE_MEAN:
+        newStruct->max = p_psArraySampleMean(myArray, maskArray, maskVal);
+        break;
+    case PS_STAT_SAMPLE_MEDIAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_SAMPLE_STDEV:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_SAMPLE_UQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_SAMPLE_LQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEAN_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEDIAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEDIAN_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MODE:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MODE_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_STDEV:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_UQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_LQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_MEAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_MEAN_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_MEAN_NSIGMA:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_STDEV:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_MAX:
+        newStruct->max = p_psArrayMax(myArray, maskArray, maskVal);
+        break;
+    case PS_STAT_MIN:
+        newStruct->min = p_psArrayMin(myArray, maskArray, maskVal);
+        break;
+    case PS_STAT_NVALUES:
+        newStruct->nValues = p_psArrayNValues(myArray, maskArray, maskVal);
+        break;
+    default:
+        psError(__func__, "Unknown options 0x%x.\n", stats->options);
+        break;
+    }
+
+    return(newStruct);
+}
Index: /trunk/psLib/src/dataManip/psStats.h
===================================================================
--- /trunk/psLib/src/dataManip/psStats.h	(revision 563)
+++ /trunk/psLib/src/dataManip/psStats.h	(revision 563)
@@ -0,0 +1,125 @@
+#if !defined(PS_STATS_H)
+#define PS_STATS_H
+
+/** \file psStats.h
+ *  \brief basic statistical operations
+ *  \ingroup MathGroup
+ */
+
+/** statistics which may be calculated */
+typedef enum {
+    PS_STAT_SAMPLE_MEAN           = 0x000001,
+    PS_STAT_SAMPLE_MEDIAN         = 0x000002,
+    PS_STAT_SAMPLE_STDEV          = 0x000004,
+    PS_STAT_SAMPLE_UQ     = 0x000008,
+    PS_STAT_SAMPLE_LQ     = 0x000010,
+    PS_STAT_ROBUST_MEAN           = 0x000020,
+    PS_STAT_ROBUST_MEAN_NVALUES   = 0x000040,
+    PS_STAT_ROBUST_MEDIAN         = 0x000080,
+    PS_STAT_ROBUST_MEDIAN_NVALUES = 0x000100,
+    PS_STAT_ROBUST_MODE           = 0x000200,
+    PS_STAT_ROBUST_MODE_NVALUES   = 0x000400,
+    PS_STAT_ROBUST_STDEV          = 0x000800,
+    PS_STAT_ROBUST_UQ     = 0x001000,
+    PS_STAT_ROBUST_LQ     = 0x002000,
+    PS_STAT_CLIPPED_MEAN          = 0x004000,
+    PS_STAT_CLIPPED_MEAN_NVALUES  = 0x008000,
+    PS_STAT_CLIPPED_MEAN_NSIGMA   = 0x010000,
+    PS_STAT_CLIPPED_STDEV         = 0x020000,
+    PS_STAT_MAX         = 0x040000,
+    PS_STAT_MIN         = 0x080000,
+    PS_STAT_NVALUES     = 0x100000
+} psStatsOptions;
+
+
+/** This is the generic statistics structure */
+typedef struct
+{
+    double sampleMean;   ///< formal mean of sample
+    double sampleMedian;  ///< formal median of sample
+    double sampleStdev;   ///< standard deviation of sample
+    double sampleUQ;   ///< upper quartile of sample
+    double sampleLQ;   ///< lower quartile of sample
+    double robustMean;   ///< robust mean of array
+    int    robustMeanNvalues;  ///< number of measurements used for robust mean
+    double robustMedian;  ///< robust median of array
+    int    robustMedianNvalues;  ///< number of measurements used for robust median
+    double robustMode;   ///< Robust mode of array
+    int    robustModeNvalues;  ///< Number of measurements used for robust mode
+    double robustStdev;   ///< robust standard deviation of array
+    double robustUQ;   ///< robust upper quartile
+    double robustLQ;   ///< robust lower quartile
+    double clippedMean;   ///< Nsigma clipped mean
+    int    clippedMeanNvalues;  ///< number of data points used for clipped mean
+    double clippedStdev;  ///< standard deviation after clipping
+    double clipSigma;   ///< Nsigma used for clipping; user input
+    int    clipIter;   ///< Number of clipping iterations; user input
+    double min;    ///< minimum data value in array
+    double max;    ///< maximum data value in array
+    int    nValues;   ///< number of data values in array
+    psStatsOptions options;  ///< bitmask of calculated values
+}
+psStats;
+
+
+/** Do Statistics on an array.  Returns a status value. \ingroup MathGroup */
+psStats *
+psArrayStats(const psFloatArray *restrict myArray, ///< Array to be analysed
+             const psIntArray *restrict maskArray, ///< Ignore elements where (maskArray & maskVal) != 0
+             ///< May be NULL
+             unsigned int maskVal, ///< Only mask elements with one of these bits set in maskArray
+             psStats *stats  ///< stats structure defines stats to be calculated and how
+            );
+
+/** Constructor */
+psStats *
+psStatsAlloc(psStatsOptions options ///< Statistics to measure
+            );
+
+/** Destructor */
+void
+psStatsFree(psStats *restrict stats ///< Stats structure to destroy
+           );
+
+/***********************************************************************************************************/
+
+/** Histograms  */
+typedef struct
+{
+    const psFloatArray *restrict lower; ///< Lower bounds for the bins
+    const psFloatArray *restrict upper; ///< Upper bounds for the bins
+    psIntArray *nums;   ///< Number in each of the bins
+    float minVal, maxVal;  ///< Minimum and maximum values
+    int minNum, maxNum;   ///< Number below the minimum and above the maximum
+}
+psHistogram;
+
+/** Constructor \ingroup MathGroup */
+psHistogram *
+psHistogramAlloc(float lower,  ///< Lower limit for the bins
+                 float upper,  ///< Upper limit for the bins
+                 float size  ///< Size of the bins
+                );
+
+/** Generic constructor \ingroup MathGroup */
+psHistogram *
+psHistogramAllocGeneric(const psFloatArray *restrict lower, ///< Lower bounds for the bins
+                        const psFloatArray *restrict upper, ///< Upper bounds for the bins
+                        float minVal, ///< Minimum value
+                        float maxVal ///< Maximum value
+                       );
+
+/** Destructor \ingroup MathGroup **/
+void
+psHistogramFree(psHistogram *restrict myHist ///< Histogram to destroy
+               );
+
+
+/** Calculate a histogram \ingroup MathGroup **/
+psHistogram *
+psGetArrayHistogram(psHistogram *restrict myHist, ///< Histogram data
+                    const psFloatArray *restrict myArray ///< Array to analyse
+                   );
+
+#endif
+
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 563)
+++ /trunk/psLib/src/math/psStats.c	(revision 563)
@@ -0,0 +1,333 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdarg.h>
+#include "pslib.h"
+//#include "psMemory.h"
+//#include "psTrace.h"
+//#include "psString.h"
+//#include "psError.h"
+#include "psArray.h"
+#include "psStats.h"
+
+/******************************************************************************
+    psStatsAlloc(): This routine must create a new psStats data structure.
+ *****************************************************************************/
+psStats *psStatsAlloc(psStatsOptions options)
+{
+    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->robustMeanNvalues = 0;
+    newStruct->robustMedian = 0.0;
+    newStruct->robustMedianNvalues = 0;
+    newStruct->robustMode = 0.0;
+    newStruct->robustModeNvalues = 0;
+    newStruct->robustStdev = 0.0;
+    newStruct->robustUQ = 0.0;
+    newStruct->robustLQ = 0.0;
+    newStruct->clippedMean = 0.0;
+    newStruct->clippedMeanNvalues = 0;
+    newStruct->clippedStdev = 0.0;
+    newStruct->clipSigma = 0.0;
+    newStruct->clipIter = 0;
+    newStruct->min = 0.0;
+    newStruct->max = 0.0;
+    newStruct->nValues = 0;
+    newStruct->options = options;
+
+    return(newStruct);
+}
+
+/******************************************************************************
+    psStatsFree(): This routine must free the psStats data structure.
+ *****************************************************************************/
+void psStatsFree(psStats *stats)
+{
+    psFree(stats);
+}
+
+/******************************************************************************
+    HISTOGRAM FUNCTIONS
+ *****************************************************************************/
+psHistogram *psHistogramAlloc(float lower,
+                              float upper,
+                              float size )
+{
+    psHistogram *newHist = NULL;
+    int numBins = 0;
+
+    numBins = 1 + ((upper - lower) / size);
+    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
+    newHist->lower = psFloatArrayAlloc(numBins);
+    newHist->upper = psFloatArrayAlloc(numBins);
+    newHist->nums = psIntArrayAlloc(numBins);
+    newHist->minVal = lower;
+    newHist->maxVal = upper;
+    newHist->minNum = 0;
+    newHist->maxNum = 0;
+
+    return(newHist);
+}
+
+psHistogram *psHistogramAllocGeneric(const psFloatArray *restrict lower,
+                                     const psFloatArray *restrict upper,
+                                     float minVal,
+                                     float maxVal)
+{
+    psHistogram *newHist = NULL;
+    int numBins = 0;
+    int i;
+
+    if (lower->n != upper->n) {
+        psError(__func__, "There is a different number of upper/lower values");
+    }
+    numBins = lower->n;
+    newHist = (psHistogram *) psAlloc(sizeof(psHistogram));
+    newHist->lower = psFloatArrayAlloc(numBins);
+    newHist->upper = psFloatArrayAlloc(numBins);
+    for (i=0;i<numBins;i++) {
+        newHist->lower->arr[i] = lower->arr[i];
+        newHist->upper->arr[i] = upper->arr[i];
+    }
+    newHist->nums = psIntArrayAlloc(numBins);
+    newHist->minVal = minVal;
+    newHist->maxVal = maxVal;
+    newHist->minNum = 0;
+    newHist->maxNum = 0;
+
+    return(newHist);
+}
+
+void psHistogramFree(psHistogram *restrict myHist)
+{
+    psFree(myHist);
+}
+
+
+/*****************************************************************************
+    NOTE: Can we assume fixed size bins in the myHist data structure?
+ 
+    NOTE: This procedure assumes that the following has been already set
+        myHist->lower
+        myHist->upper
+ myHist->minVal
+ myHist->maxVal
+    This procedure sets
+        myHist->nums
+        myHist->minNum
+        myHist->maxNum
+ *****************************************************************************/
+psHistogram *psGetArrayHistogram(       psHistogram *restrict myHist,
+                                        const psFloatArray *restrict myArray)
+{
+    int i = 0;
+    float binSize = 0.0;
+    int binNum = 0;
+
+    binSize = (myHist->maxVal - myHist->minVal) / (float) myHist->nums->n;
+    for (i=0;i<myArray->n;i++) {
+        if (myArray->arr[i] < myHist->minVal) {
+            myHist->minNum++;
+        } else if (myArray->arr[i] > myHist->maxVal) {
+            myHist->maxNum++;
+        } else {
+            binNum = (int) ((myArray->arr[i] - myHist->minVal) / binSize);
+            if ((myArray->arr[i] >= myHist->lower->arr[i]) &&
+                    (myArray->arr[i] <= myHist->upper->arr[i])) {
+                myHist->nums->arr[i]++;
+            } else {
+                psError(__func__, "data value was not within the bounds of the bin it should habe been in.");
+            }
+        }
+    }
+    return(myHist);
+}
+
+/******************************************************************************
+    MISC STATISTICAL FUNCTIONS
+ *****************************************************************************/
+float p_psArraySampleMean(const psFloatArray *restrict myArray,
+                          const psIntArray *restrict maskArray,
+                          unsigned int maskVal)
+{
+    int i = 0;
+    float mean = 0.0;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                mean+= maskArray->arr[i];
+            }
+        }
+    } else {
+        for (i=0;i<myArray->n;i++) {
+            mean+= maskArray->arr[i];
+        }
+    }
+    return(mean);
+}
+
+float p_psArrayMax(const psFloatArray *restrict myArray,
+                   const psIntArray *restrict maskArray,
+                   unsigned int maskVal)
+{
+    int i = 0;
+    float max = -1e99;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                if (myArray->arr[i] > max) {
+                    max = maskArray->arr[i];
+                }
+            }
+        }
+    } else {
+        for (i=0;i<myArray->n;i++) {
+            if (myArray->arr[i] > max) {
+                max = maskArray->arr[i];
+            }
+        }
+    }
+    return(max);
+}
+
+float p_psArrayMin(const psFloatArray *restrict myArray,
+                   const psIntArray *restrict maskArray,
+                   unsigned int maskVal)
+{
+    int i = 0;
+    float min = 1e99;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                if (myArray->arr[i] < min) {
+                    min = maskArray->arr[i];
+                }
+            }
+        }
+    } else {
+        for (i=0;i<myArray->n;i++) {
+            if (myArray->arr[i] < min) {
+                min = maskArray->arr[i];
+            }
+        }
+    }
+    return(min);
+}
+
+int p_psArrayNValues(const psFloatArray *restrict myArray,
+                     const psIntArray *restrict maskArray,
+                     unsigned int maskVal)
+{
+    int i = 0;
+    int numData = 0;
+
+    if (maskArray != NULL) {
+        for (i=0;i<myArray->n;i++) {
+            if (!(maskVal & maskArray->arr[i])) {
+                numData++;
+            }
+        }
+    } else {
+        numData = myArray->n;
+    }
+    return(numData);
+}
+
+psStats *psArrayStats(const psFloatArray *restrict myArray,
+                      const psIntArray *restrict maskArray,
+                      unsigned int maskVal,
+                      psStats *stats)
+{
+    psStats *newStruct = NULL;
+
+    newStruct = psStatsAlloc(stats->options);
+    if (myArray == NULL) {
+        psError(__func__,
+                "Input data array (myArray) was NULL.");
+    }
+
+    if ((maskArray != NULL) &&
+            (myArray->n != maskArray->n)) {
+        psError(__func__, "Array data and array mask are of different sizes.");
+    }
+
+    switch (stats->options) {
+    case PS_STAT_SAMPLE_MEAN:
+        newStruct->max = p_psArraySampleMean(myArray, maskArray, maskVal);
+        break;
+    case PS_STAT_SAMPLE_MEDIAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_SAMPLE_STDEV:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_SAMPLE_UQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_SAMPLE_LQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEAN_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEDIAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MEDIAN_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MODE:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_MODE_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_STDEV:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_UQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_ROBUST_LQ:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_MEAN:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_MEAN_NVALUES:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_MEAN_NSIGMA:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_CLIPPED_STDEV:
+        printf("Code me!\n");
+        break;
+    case PS_STAT_MAX:
+        newStruct->max = p_psArrayMax(myArray, maskArray, maskVal);
+        break;
+    case PS_STAT_MIN:
+        newStruct->min = p_psArrayMin(myArray, maskArray, maskVal);
+        break;
+    case PS_STAT_NVALUES:
+        newStruct->nValues = p_psArrayNValues(myArray, maskArray, maskVal);
+        break;
+    default:
+        psError(__func__, "Unknown options 0x%x.\n", stats->options);
+        break;
+    }
+
+    return(newStruct);
+}
Index: /trunk/psLib/src/math/psStats.h
===================================================================
--- /trunk/psLib/src/math/psStats.h	(revision 563)
+++ /trunk/psLib/src/math/psStats.h	(revision 563)
@@ -0,0 +1,125 @@
+#if !defined(PS_STATS_H)
+#define PS_STATS_H
+
+/** \file psStats.h
+ *  \brief basic statistical operations
+ *  \ingroup MathGroup
+ */
+
+/** statistics which may be calculated */
+typedef enum {
+    PS_STAT_SAMPLE_MEAN           = 0x000001,
+    PS_STAT_SAMPLE_MEDIAN         = 0x000002,
+    PS_STAT_SAMPLE_STDEV          = 0x000004,
+    PS_STAT_SAMPLE_UQ     = 0x000008,
+    PS_STAT_SAMPLE_LQ     = 0x000010,
+    PS_STAT_ROBUST_MEAN           = 0x000020,
+    PS_STAT_ROBUST_MEAN_NVALUES   = 0x000040,
+    PS_STAT_ROBUST_MEDIAN         = 0x000080,
+    PS_STAT_ROBUST_MEDIAN_NVALUES = 0x000100,
+    PS_STAT_ROBUST_MODE           = 0x000200,
+    PS_STAT_ROBUST_MODE_NVALUES   = 0x000400,
+    PS_STAT_ROBUST_STDEV          = 0x000800,
+    PS_STAT_ROBUST_UQ     = 0x001000,
+    PS_STAT_ROBUST_LQ     = 0x002000,
+    PS_STAT_CLIPPED_MEAN          = 0x004000,
+    PS_STAT_CLIPPED_MEAN_NVALUES  = 0x008000,
+    PS_STAT_CLIPPED_MEAN_NSIGMA   = 0x010000,
+    PS_STAT_CLIPPED_STDEV         = 0x020000,
+    PS_STAT_MAX         = 0x040000,
+    PS_STAT_MIN         = 0x080000,
+    PS_STAT_NVALUES     = 0x100000
+} psStatsOptions;
+
+
+/** This is the generic statistics structure */
+typedef struct
+{
+    double sampleMean;   ///< formal mean of sample
+    double sampleMedian;  ///< formal median of sample
+    double sampleStdev;   ///< standard deviation of sample
+    double sampleUQ;   ///< upper quartile of sample
+    double sampleLQ;   ///< lower quartile of sample
+    double robustMean;   ///< robust mean of array
+    int    robustMeanNvalues;  ///< number of measurements used for robust mean
+    double robustMedian;  ///< robust median of array
+    int    robustMedianNvalues;  ///< number of measurements used for robust median
+    double robustMode;   ///< Robust mode of array
+    int    robustModeNvalues;  ///< Number of measurements used for robust mode
+    double robustStdev;   ///< robust standard deviation of array
+    double robustUQ;   ///< robust upper quartile
+    double robustLQ;   ///< robust lower quartile
+    double clippedMean;   ///< Nsigma clipped mean
+    int    clippedMeanNvalues;  ///< number of data points used for clipped mean
+    double clippedStdev;  ///< standard deviation after clipping
+    double clipSigma;   ///< Nsigma used for clipping; user input
+    int    clipIter;   ///< Number of clipping iterations; user input
+    double min;    ///< minimum data value in array
+    double max;    ///< maximum data value in array
+    int    nValues;   ///< number of data values in array
+    psStatsOptions options;  ///< bitmask of calculated values
+}
+psStats;
+
+
+/** Do Statistics on an array.  Returns a status value. \ingroup MathGroup */
+psStats *
+psArrayStats(const psFloatArray *restrict myArray, ///< Array to be analysed
+             const psIntArray *restrict maskArray, ///< Ignore elements where (maskArray & maskVal) != 0
+             ///< May be NULL
+             unsigned int maskVal, ///< Only mask elements with one of these bits set in maskArray
+             psStats *stats  ///< stats structure defines stats to be calculated and how
+            );
+
+/** Constructor */
+psStats *
+psStatsAlloc(psStatsOptions options ///< Statistics to measure
+            );
+
+/** Destructor */
+void
+psStatsFree(psStats *restrict stats ///< Stats structure to destroy
+           );
+
+/***********************************************************************************************************/
+
+/** Histograms  */
+typedef struct
+{
+    const psFloatArray *restrict lower; ///< Lower bounds for the bins
+    const psFloatArray *restrict upper; ///< Upper bounds for the bins
+    psIntArray *nums;   ///< Number in each of the bins
+    float minVal, maxVal;  ///< Minimum and maximum values
+    int minNum, maxNum;   ///< Number below the minimum and above the maximum
+}
+psHistogram;
+
+/** Constructor \ingroup MathGroup */
+psHistogram *
+psHistogramAlloc(float lower,  ///< Lower limit for the bins
+                 float upper,  ///< Upper limit for the bins
+                 float size  ///< Size of the bins
+                );
+
+/** Generic constructor \ingroup MathGroup */
+psHistogram *
+psHistogramAllocGeneric(const psFloatArray *restrict lower, ///< Lower bounds for the bins
+                        const psFloatArray *restrict upper, ///< Upper bounds for the bins
+                        float minVal, ///< Minimum value
+                        float maxVal ///< Maximum value
+                       );
+
+/** Destructor \ingroup MathGroup **/
+void
+psHistogramFree(psHistogram *restrict myHist ///< Histogram to destroy
+               );
+
+
+/** Calculate a histogram \ingroup MathGroup **/
+psHistogram *
+psGetArrayHistogram(psHistogram *restrict myHist, ///< Histogram data
+                    const psFloatArray *restrict myArray ///< Array to analyse
+                   );
+
+#endif
+