Changeset 1293

Timestamp:

Jul 23, 2004, 4:46:45 PM (22 years ago)

Author:

gusciora

Message:

...

Location:

trunk/psLib/src

Files:

• astro/psCoord.c (modified) (5 diffs)
• astro/psCoord.h (modified) (5 diffs)
• astronomy/psCoord.c (modified) (5 diffs)
• astronomy/psCoord.h (modified) (5 diffs)
• dataManip/psStats.c (modified) (17 diffs)
• math/psStats.c (modified) (17 diffs)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-23 03:13:39 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-24 02:42:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22 +22 @@
 #include "psCoord.h"
 #include "psMemory.h"
-
+#include "psAbort.h"
+#include <math.h>
+#include <float.h>
 
 psPlane *psPlaneTransformApply(psPlane *out,
@@ -127 +129 @@
 }
 
+// XXX: Is this the correct way to calculate this?
+float cot(float x)
+{
+    return(1.0 / atan(x));
+}
+
+float arg(float x, float y)
+{
+    if (x > 0) {
+        return(atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(-0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(M_PI + atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(-M_PI + atan(y/x));
+    }
+
+    psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
+    return(0.0);
+}
+
+
+psPlane *psProject(const psSphere *coord,
+                   const psProjection *projection)
+{
+    float R = 0.0;
+    float alpha = 0.0;
+    psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = cot(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = cos(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->x = coord->d;
+        tmp->y = coord->r;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->x = coord->d;
+        tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        alpha = 1.0 / ((180.0 / M_PI) *
+                       sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5)));
+
+        tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
+        tmp->y = alpha * sin(coord->d);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
+
+psSphere *psDeproject(const psPlane *coord,
+                      const psProjection *projection)
+{
+    float R = 0.0;
+    float chu = 0.0;
+    float chu1 = 0.0;
+    float chu2 = 0.0;
+    psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = atan(180.0 / (R * M_PI));
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = acos((R * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->d = coord->x;
+        tmp->r = coord->y;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->d = coord->x;
+        tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        chu1 = (coord->x * M_PI) / 720.0;
+        chu1*= chu1;
+        chu2 = (coord->y * M_PI) / 360.0;
+        chu2*= chu2;
+        chu = sqrt(1.0 - chu1 - chu2);
+        tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0,
+                           (coord->x * chu * M_PI)/360.0);
+        tmp->r = asin((coord->y * chu * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
 
 psSphere *psSphereGetOffset(const psSphere *restrict position1,
@@ -133 +248 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
-
-
+    psAbort(__func__, "This function has not be dedfined.\n");
+    return(NULL);
+}
+
+
+// XXX: I copied the algorithm from the ADD exactly.  Arguments mode and unit
+// are ignored.
 psSphere *psSphereSetOffset(const psSphere *restrict position,
                             const psSphere *restrict offset,
@@ -142 +260 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
+    psPlane lin;
+    psSphere *tmp;
+    psProjection proj;
+
+    if (mode == PS_LINEAR) {
+        proj.R = position->r;
+        proj.D = position->d;
+        proj.Xs = 0.0;
+        proj.Ys = 0.0;
+        proj.type = PS_PROJ_TAN;
+
+        lin.x = offset->r;
+        lin.y = offset->d;
+
+        tmp = psDeproject(&lin, &proj);
+        return(tmp);
+    } else if (mode == PS_SPHERICAL) {
+        psAbort(__func__, "Sperical offset modes are not defined.\n");
+    }
+    psAbort(__func__, "Unrecognized offset mode\n");
+    return(NULL);
+}

trunk/psLib/src/astro/psCoord.h

@@ -10 +10 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-23 03:13:39 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-24 02:42:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 typedef struct
 {
-    double x;    ///< x position
-    double y;    ///< y position
+    double x;      ///< x position
+    double y;      ///< y position
     double xErr;   ///< Error in x position
     double yErr;   ///< Error in y position
@@ -36 +36 @@
 typedef struct
 {
-    double r;    ///< RA
-    double d;    ///< Dec
+    double r;      ///< RA
+    double d;      ///< Dec
     double rErr;   ///< Error in RA
     double dErr;   ///< Error in Dec
@@ -63 +63 @@
     double sinNPlat;   ///< sin of North Pole latitude
     double cosNPlat;   ///< cos of North Pole latitude
-    double sinZP;   ///< sin of Forst PT os Ares lon
-    double cosZP;   ///< cos of Forst PT os Ares lon
+    double sinZP;      ///< sin of Forst PT os Ares lon
+    double cosZP;      ///< cos of Forst PT os Ares lon
 }
 psSphereTransform;
 
 typedef enum {
-    PS_PROJ_TAN,   ///< Tangent projection
-    PS_PROJ_SIN,   ///< Sine projection
-    PS_PROJ_AIT,   ///< Aitoff projection
-    PS_PROJ_PAR,   ///< Par projection
-    PS_PROJ_GLS,   ///< GLS projection
-    PS_PROJ_NTYPE   ///< Number of types; must be last.
+    PS_PROJ_TAN,       ///< Tangent projection
+    PS_PROJ_SIN,       ///< Sine projection
+    PS_PROJ_AIT,       ///< Aitoff projection
+    PS_PROJ_PAR,       ///< Par projection
+    PS_PROJ_GLS,       ///< GLS projection
+    PS_PROJ_CAR,       ///< CAR projection
+    PS_PROJ_MER,       ///< MER projection
+    PS_PROJ_NTYPE      ///< Number of types; must be last.
 } psProjectionType;
 
 typedef struct
 {
-    double R;    ///< Coordinates of projection center
-    double D;    ///< Coordinates of projection center
+    double R;     ///< Coordinates of projection center
+    double D;     ///< Coordinates of projection center
     double Xs;    ///< plate-scale in X direction
     double Ys;    ///< plate-scale in Y direction
@@ -128 +130 @@
                    const psProjection *projection);
 
-psSphere *psProject(const psPlane *coord,
-                    const psProjection *projection);
+psSphere *psDeproject(const psPlane *coord,
+                      const psProjection *projection);
 
 psSphere *psSphereGetOffset(const psSphere *restrict position1,

trunk/psLib/src/astronomy/psCoord.c

@@ -10 +10 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-23 03:13:39 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-24 02:42:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22 +22 @@
 #include "psCoord.h"
 #include "psMemory.h"
-
+#include "psAbort.h"
+#include <math.h>
+#include <float.h>
 
 psPlane *psPlaneTransformApply(psPlane *out,
@@ -127 +129 @@
 }
 
+// XXX: Is this the correct way to calculate this?
+float cot(float x)
+{
+    return(1.0 / atan(x));
+}
+
+float arg(float x, float y)
+{
+    if (x > 0) {
+        return(atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(-0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(M_PI + atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(-M_PI + atan(y/x));
+    }
+
+    psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
+    return(0.0);
+}
+
+
+psPlane *psProject(const psSphere *coord,
+                   const psProjection *projection)
+{
+    float R = 0.0;
+    float alpha = 0.0;
+    psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = cot(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = cos(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->x = coord->d;
+        tmp->y = coord->r;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->x = coord->d;
+        tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        alpha = 1.0 / ((180.0 / M_PI) *
+                       sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5)));
+
+        tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
+        tmp->y = alpha * sin(coord->d);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
+
+psSphere *psDeproject(const psPlane *coord,
+                      const psProjection *projection)
+{
+    float R = 0.0;
+    float chu = 0.0;
+    float chu1 = 0.0;
+    float chu2 = 0.0;
+    psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = atan(180.0 / (R * M_PI));
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = acos((R * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->d = coord->x;
+        tmp->r = coord->y;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->d = coord->x;
+        tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        chu1 = (coord->x * M_PI) / 720.0;
+        chu1*= chu1;
+        chu2 = (coord->y * M_PI) / 360.0;
+        chu2*= chu2;
+        chu = sqrt(1.0 - chu1 - chu2);
+        tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0,
+                           (coord->x * chu * M_PI)/360.0);
+        tmp->r = asin((coord->y * chu * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
 
 psSphere *psSphereGetOffset(const psSphere *restrict position1,
@@ -133 +248 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
-
-
+    psAbort(__func__, "This function has not be dedfined.\n");
+    return(NULL);
+}
+
+
+// XXX: I copied the algorithm from the ADD exactly.  Arguments mode and unit
+// are ignored.
 psSphere *psSphereSetOffset(const psSphere *restrict position,
                             const psSphere *restrict offset,
@@ -142 +260 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
+    psPlane lin;
+    psSphere *tmp;
+    psProjection proj;
+
+    if (mode == PS_LINEAR) {
+        proj.R = position->r;
+        proj.D = position->d;
+        proj.Xs = 0.0;
+        proj.Ys = 0.0;
+        proj.type = PS_PROJ_TAN;
+
+        lin.x = offset->r;
+        lin.y = offset->d;
+
+        tmp = psDeproject(&lin, &proj);
+        return(tmp);
+    } else if (mode == PS_SPHERICAL) {
+        psAbort(__func__, "Sperical offset modes are not defined.\n");
+    }
+    psAbort(__func__, "Unrecognized offset mode\n");
+    return(NULL);
+}

trunk/psLib/src/astronomy/psCoord.h

@@ -10 +10 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-23 03:13:39 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-24 02:42:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 typedef struct
 {
-    double x;    ///< x position
-    double y;    ///< y position
+    double x;      ///< x position
+    double y;      ///< y position
     double xErr;   ///< Error in x position
     double yErr;   ///< Error in y position
@@ -36 +36 @@
 typedef struct
 {
-    double r;    ///< RA
-    double d;    ///< Dec
+    double r;      ///< RA
+    double d;      ///< Dec
     double rErr;   ///< Error in RA
     double dErr;   ///< Error in Dec
@@ -63 +63 @@
     double sinNPlat;   ///< sin of North Pole latitude
     double cosNPlat;   ///< cos of North Pole latitude
-    double sinZP;   ///< sin of Forst PT os Ares lon
-    double cosZP;   ///< cos of Forst PT os Ares lon
+    double sinZP;      ///< sin of Forst PT os Ares lon
+    double cosZP;      ///< cos of Forst PT os Ares lon
 }
 psSphereTransform;
 
 typedef enum {
-    PS_PROJ_TAN,   ///< Tangent projection
-    PS_PROJ_SIN,   ///< Sine projection
-    PS_PROJ_AIT,   ///< Aitoff projection
-    PS_PROJ_PAR,   ///< Par projection
-    PS_PROJ_GLS,   ///< GLS projection
-    PS_PROJ_NTYPE   ///< Number of types; must be last.
+    PS_PROJ_TAN,       ///< Tangent projection
+    PS_PROJ_SIN,       ///< Sine projection
+    PS_PROJ_AIT,       ///< Aitoff projection
+    PS_PROJ_PAR,       ///< Par projection
+    PS_PROJ_GLS,       ///< GLS projection
+    PS_PROJ_CAR,       ///< CAR projection
+    PS_PROJ_MER,       ///< MER projection
+    PS_PROJ_NTYPE      ///< Number of types; must be last.
 } psProjectionType;
 
 typedef struct
 {
-    double R;    ///< Coordinates of projection center
-    double D;    ///< Coordinates of projection center
+    double R;     ///< Coordinates of projection center
+    double D;     ///< Coordinates of projection center
     double Xs;    ///< plate-scale in X direction
     double Ys;    ///< plate-scale in Y direction
@@ -128 +130 @@
                    const psProjection *projection);
 
-psSphere *psProject(const psPlane *coord,
-                    const psProjection *projection);
+psSphere *psDeproject(const psPlane *coord,
+                      const psProjection *projection);
 
 psSphere *psSphereGetOffset(const psSphere *restrict position1,

trunk/psLib/src/dataManip/psStats.c

@@ -131 +131 @@
 mean of the input vector.
 Inputs
+    myVector
+    maskVector
+    maskVal
+    stats
 Returns
     NULL
@@ -404 +408 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    psStats *stats2 = NULL;                     // Temporary stats structure
 
 
     // Determine if the number of data points exceed a threshold which will
     // cause to generate robust stats, as opposed to exact stats.
-
-    if (myVector->n > stats->sampleLimit) {
-        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
-
-        // Calculate the robust quartiles.
-        stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-
-        // Store the robust quartiles into the sample quartile members.
-        stats->sampleMedian = stats2->robustMedian;
-
-        // Free temporary data buffers.
-        psFree(stats2);
-
-        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
-        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
-
-        return;
-    }
+    // XXX: This code is no longer used.  We now calculate the median exactly
+    // regardless of the vector size.
+    /*
+        if (myVector->n > stats->sampleLimit) {
+            psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
+     
+            // Calculate the robust quartiles.
+            stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+            p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
+     
+            // Store the robust quartiles into the sample quartile members.
+            stats->sampleMedian = stats2->robustMedian;
+     
+            // Free temporary data buffers.
+            psFree(stats2);
+     
+            // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
+            stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
+     
+            return;
+        }
+    */
 
     // Determine how many data points fit inside this min/max range
@@ -483 +489 @@
 
     // Calculate the median exactly.
+    // XXX: Is this the correct action?
     if (0 == (nValues % 2)) {
         stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] +
@@ -498 +505 @@
     This routine smoothes the data in the input robustHistogram with a
     Gaussian of width sigma.
+ 
+    XXX: Must consult with IfA on the proper values for the Gaussian
+    smoothing.  Currently, the Gaussian coefficients are such that only
+    the center coefficient is non-zero.  This is because "e" is being
+    raised to a very large negative number for all points except the
+    center point.
+ 
+    XXX: use a static variable for gaussianCoefs[] and compute them once.
  *****************************************************************************/
 psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
@@ -541 +556 @@
     }
 
+    // Perform the actual smoothing.
     for(i=0;i<robustHistogram->nums->n;i++) {
         smooth->data.F32[i] = 0.0;
@@ -576 +592 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-
-    // Determine if the number of data points exceed a threshold which will
-    // cause to generate robust stats, as opposed to exact stats.
-    if (myVector->n > stats->sampleLimit) {
-        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
-        psStats *stats2 = NULL;
-        // Calculate the robust quartiles.
-        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
-        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-
-        // Store the robust quartiles into the sample quartile members.
-        stats->sampleUQ = stats2->robustUQ;
-        stats->sampleLQ = stats2->robustLQ;
-
-        // Free temporary data buffers.
-        psFree(stats2);
-
-        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
-        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
-
-        return;
-    }
 
     // Determine how many data points fit inside this min/max range
@@ -861 +855 @@
 
 /*****************************************************************************
+p_psNormalizeVector(myData): this is a private function which normalizes the
+elements of a vector to a range between 0.0 and 1.0.
  *****************************************************************************/
 void p_psNormalizeVector(psVector *myData)
@@ -885 +881 @@
 
 
+/*****************************************************************************
+p_psGaussian(myData, myParams): an internal function, used by robust stats,
+intended to compute the Gaussian with a specified sigma/mean, at the
+specified data point.
+ *****************************************************************************/
 float p_psGaussian(const psVector *restrict myData,
                    const psVector *restrict myParams)
@@ -894 +895 @@
     tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
 
-    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
+    //    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
     return(tmp);
 }
 
+/*****************************************************************************
+p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which
+calculates the specified partial derivative of the above Gaussian function.
+ *****************************************************************************/
 float p_psGaussianDeriv(const psVector *restrict myData,
                         const psVector *restrict myParams,
@@ -922 +927 @@
 
 
+/*****************************************************************************
+p_psQuadratic(myData, myParams): an internal function, used by robust stats,
+intended to compute a quadratic, with the specified parameters, at the
+specified data point.
+ *****************************************************************************/
 float p_psQuadratic(const psVector *restrict myParams,
                     const psVector *restrict myCoords)
@@ -935 +945 @@
 }
 
+/*****************************************************************************
+p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which
+calculates the specified partial derivative of the above quadratic function.
+ *****************************************************************************/
 float p_psQuadraticDeriv(const psVector *restrict myParams,
                          const psVector *restrict myCoords,
@@ -959 +973 @@
 [rangeLow, rangeHigh].  It determines the x-value of that polynomial such
 that f(x) == midpoint.  This functions uses a binary-search algorithm on the
-range and assumes that the polnomial is monotonically increasing within that
-range.
+range and assumes that the polynomial is monotonically increasing or
+decreasing within that range.
  *****************************************************************************/
 float p_ps1DPolyMedian(psPolynomial1D *myPoly,
@@ -998 +1012 @@
 
 /******************************************************************************
-p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
-as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
-polynomials here) and a range of x-values for which it is defined:
-[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
-that f(x) == midpoint.  This functions uses a binary-search algorithm on the
-range and assumes that the polnomial is monotonically increasing within that
-range.
- *****************************************************************************/
-float p_psFitQuadratic(psHistogram *robustHistogram,
+p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given
+the specified histogram, with the specified pre-calculated cumulativeSums,
+this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum",
+and then finds the point x such that f(x) == fitFloat.
+ 
+XXX: This function is currently not being used.
+ *****************************************************************************/
+float p_psFitQuadratic(psHistogram *histogram,
+                       psVector *cumulativeSums,
                        int binNum,
                        float fitFloat)
 {
-    /*
-        if ((binNum > 0) &&
-            (binNum < (robustHistogram->nums->n+1))) {
-            x->data.F64[0] = (double) 0.5 *
-                      (robustHistogram->bounds->data.F32[binNum-1] +
-                       robustHistogram->bounds->data.F32[binNum]);
-            x->data.F64[1] = (double) 0.5 *
-                      (robustHistogram->bounds->data.F32[binNum] +
-                       robustHistogram->bounds->data.F32[binNum+1]);
-            x->data.F64[2] = (double) 0.5 *
-                      (robustHistogram->bounds->data.F32[binNum+1] +
-                       robustHistogram->bounds->data.F32[binNum+2]);
-     
-            y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1];
-            y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum];
-            y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1];
-     
-            if (!((y->data.F64[0] <= fitFloat) &&
-                 (fitFloat <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
-            }
-     
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-     
-            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
-            return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat);
-         } else {
-            return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] +
-                          robustHistogram->bounds->data.F32[binNum]));
-        }
-    */
+    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
+    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
+    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
+    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
+
+    if ((binNum > 0) &&
+            (binNum < (histogram->nums->n+1))) {
+        x->data.F64[0] = (double) 0.5 *
+                         (histogram->bounds->data.F32[binNum-1] +
+                          histogram->bounds->data.F32[binNum]);
+        x->data.F64[1] = (double) 0.5 *
+                         (histogram->bounds->data.F32[binNum] +
+                          histogram->bounds->data.F32[binNum+1]);
+        x->data.F64[2] = (double) 0.5 *
+                         (histogram->bounds->data.F32[binNum+1] +
+                          histogram->bounds->data.F32[binNum+2]);
+
+        y->data.F64[0] = cumulativeSums->data.F32[binNum-1];
+        y->data.F64[1] = cumulativeSums->data.F32[binNum];
+        y->data.F64[2] = cumulativeSums->data.F32[binNum+1];
+
+        if (!((y->data.F64[0] <= fitFloat) &&
+                (fitFloat <= y->data.F64[2]))) {
+            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
+        }
+
+        yErr->data.F64[0] = 1.0;
+        yErr->data.F64[1] = 1.0;
+        yErr->data.F64[2] = 1.0;
+
+        myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
+        return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2],
+                                fitFloat));
+    } else {
+        return(0.5 * (histogram->bounds->data.F32[binNum+1] +
+                      histogram->bounds->data.F32[binNum]));
+    }
+
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+    psFree(myPoly);
     return(0.0);
 }
@@ -1049 +1072 @@
 p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
 calculates a variety of robust stat measures:
-PS_STAT_ROBUST_MEAN
-PS_STAT_ROBUST_MEDIAN
-PS_STAT_ROBUST_MODE
-PS_STAT_ROBUST_STDEV
-PS_STAT_ROBUST_QUARTILE
+    PS_STAT_ROBUST_MEAN
+    PS_STAT_ROBUST_MEDIAN
+    PS_STAT_ROBUST_MODE
+    PS_STAT_ROBUST_STDEV
+    PS_STAT_ROBUST_QUARTILE
 I have included all that computation in a single function, as opposed to
-breaking it across several functions for one primary reason: 
-they all
+breaking it across several functions for one primary reason:  they all
 require the same basic initial processing steps (calculate the histogram,
-etc.) 
-however there is no currently defined means, in the SDRS, to
-specify this computation as a separate step.  If the above robust stat
-measures were calcualted in separate functiosn, then much of the initial
-processing would be duplicated.
+etc.)  however there is no currently defined means, in the SDRS, to specify
+this computation as a separate step.  If the above robust stat measures were
+calcualted in separate functiosn, then much of the initial processing would
+be duplicated.
+ 
 Inputs
-myVector
-maskVector
-maskVal
-stats
+    myVector
+    maskVector
+    maskVal
+    stats
 Returns
-NULL
+    NULL
 *****************************************************************************/
 void p_psVectorRobustStats(const psVector *restrict myVector,
@@ -1224 +1246 @@
         sumN50+= (float) robustHistogram->nums->data.S32[i];
     }
+
     // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
     dL = (UQBinNum - LQBinNum) / 4;
@@ -1244 +1267 @@
                 robustHistogramVector->data.F32[i];
             cumulativeMedian+= robustHistogramVector->data.F32[i];
-            sumNfit+= (float) robustHistogramVector->data.S32[i];
+            sumNfit+= (float) robustHistogram->nums->data.S32[i];
         }
     }
@@ -1780 +1803 @@
     }
 
+    // XXX: Should we abort if (stats->min == stats->max)?
+    if (stats->min >= stats->max) {
+        psAbort(__func__, "psVectorStats() called with range: %f to %f\n",
+                stats->min, stats->max);
+    }
+
+
+
     PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
     if (mask != NULL) {

trunk/psLib/src/math/psStats.c

@@ -131 +131 @@
 mean of the input vector.
 Inputs
+    myVector
+    maskVector
+    maskVal
+    stats
 Returns
     NULL
@@ -404 +408 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-    psStats *stats2 = NULL;                     // Temporary stats structure
 
 
     // Determine if the number of data points exceed a threshold which will
     // cause to generate robust stats, as opposed to exact stats.
-
-    if (myVector->n > stats->sampleLimit) {
-        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
-
-        // Calculate the robust quartiles.
-        stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-
-        // Store the robust quartiles into the sample quartile members.
-        stats->sampleMedian = stats2->robustMedian;
-
-        // Free temporary data buffers.
-        psFree(stats2);
-
-        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
-        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
-
-        return;
-    }
+    // XXX: This code is no longer used.  We now calculate the median exactly
+    // regardless of the vector size.
+    /*
+        if (myVector->n > stats->sampleLimit) {
+            psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
+     
+            // Calculate the robust quartiles.
+            stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+            p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
+     
+            // Store the robust quartiles into the sample quartile members.
+            stats->sampleMedian = stats2->robustMedian;
+     
+            // Free temporary data buffers.
+            psFree(stats2);
+     
+            // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
+            stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
+     
+            return;
+        }
+    */
 
     // Determine how many data points fit inside this min/max range
@@ -483 +489 @@
 
     // Calculate the median exactly.
+    // XXX: Is this the correct action?
     if (0 == (nValues % 2)) {
         stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] +
@@ -498 +505 @@
     This routine smoothes the data in the input robustHistogram with a
     Gaussian of width sigma.
+ 
+    XXX: Must consult with IfA on the proper values for the Gaussian
+    smoothing.  Currently, the Gaussian coefficients are such that only
+    the center coefficient is non-zero.  This is because "e" is being
+    raised to a very large negative number for all points except the
+    center point.
+ 
+    XXX: use a static variable for gaussianCoefs[] and compute them once.
  *****************************************************************************/
 psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
@@ -541 +556 @@
     }
 
+    // Perform the actual smoothing.
     for(i=0;i<robustHistogram->nums->n;i++) {
         smooth->data.F32[i] = 0.0;
@@ -576 +592 @@
     float rangeMin = 0.0;                       // Exclude data below this
     float rangeMax = 0.0;                       // Exclude date above this
-
-    // Determine if the number of data points exceed a threshold which will
-    // cause to generate robust stats, as opposed to exact stats.
-    if (myVector->n > stats->sampleLimit) {
-        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
-        psStats *stats2 = NULL;
-        // Calculate the robust quartiles.
-        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
-        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-
-        // Store the robust quartiles into the sample quartile members.
-        stats->sampleUQ = stats2->robustUQ;
-        stats->sampleLQ = stats2->robustLQ;
-
-        // Free temporary data buffers.
-        psFree(stats2);
-
-        // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
-        stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
-
-        return;
-    }
 
     // Determine how many data points fit inside this min/max range
@@ -861 +855 @@
 
 /*****************************************************************************
+p_psNormalizeVector(myData): this is a private function which normalizes the
+elements of a vector to a range between 0.0 and 1.0.
  *****************************************************************************/
 void p_psNormalizeVector(psVector *myData)
@@ -885 +881 @@
 
 
+/*****************************************************************************
+p_psGaussian(myData, myParams): an internal function, used by robust stats,
+intended to compute the Gaussian with a specified sigma/mean, at the
+specified data point.
+ *****************************************************************************/
 float p_psGaussian(const psVector *restrict myData,
                    const psVector *restrict myParams)
@@ -894 +895 @@
     tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
 
-    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
+    //    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
     return(tmp);
 }
 
+/*****************************************************************************
+p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which
+calculates the specified partial derivative of the above Gaussian function.
+ *****************************************************************************/
 float p_psGaussianDeriv(const psVector *restrict myData,
                         const psVector *restrict myParams,
@@ -922 +927 @@
 
 
+/*****************************************************************************
+p_psQuadratic(myData, myParams): an internal function, used by robust stats,
+intended to compute a quadratic, with the specified parameters, at the
+specified data point.
+ *****************************************************************************/
 float p_psQuadratic(const psVector *restrict myParams,
                     const psVector *restrict myCoords)
@@ -935 +945 @@
 }
 
+/*****************************************************************************
+p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which
+calculates the specified partial derivative of the above quadratic function.
+ *****************************************************************************/
 float p_psQuadraticDeriv(const psVector *restrict myParams,
                          const psVector *restrict myCoords,
@@ -959 +973 @@
 [rangeLow, rangeHigh].  It determines the x-value of that polynomial such
 that f(x) == midpoint.  This functions uses a binary-search algorithm on the
-range and assumes that the polnomial is monotonically increasing within that
-range.
+range and assumes that the polynomial is monotonically increasing or
+decreasing within that range.
  *****************************************************************************/
 float p_ps1DPolyMedian(psPolynomial1D *myPoly,
@@ -998 +1012 @@
 
 /******************************************************************************
-p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
-as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
-polynomials here) and a range of x-values for which it is defined:
-[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
-that f(x) == midpoint.  This functions uses a binary-search algorithm on the
-range and assumes that the polnomial is monotonically increasing within that
-range.
- *****************************************************************************/
-float p_psFitQuadratic(psHistogram *robustHistogram,
+p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given
+the specified histogram, with the specified pre-calculated cumulativeSums,
+this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum",
+and then finds the point x such that f(x) == fitFloat.
+ 
+XXX: This function is currently not being used.
+ *****************************************************************************/
+float p_psFitQuadratic(psHistogram *histogram,
+                       psVector *cumulativeSums,
                        int binNum,
                        float fitFloat)
 {
-    /*
-        if ((binNum > 0) &&
-            (binNum < (robustHistogram->nums->n+1))) {
-            x->data.F64[0] = (double) 0.5 *
-                      (robustHistogram->bounds->data.F32[binNum-1] +
-                       robustHistogram->bounds->data.F32[binNum]);
-            x->data.F64[1] = (double) 0.5 *
-                      (robustHistogram->bounds->data.F32[binNum] +
-                       robustHistogram->bounds->data.F32[binNum+1]);
-            x->data.F64[2] = (double) 0.5 *
-                      (robustHistogram->bounds->data.F32[binNum+1] +
-                       robustHistogram->bounds->data.F32[binNum+2]);
-     
-            y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1];
-            y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum];
-            y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1];
-     
-            if (!((y->data.F64[0] <= fitFloat) &&
-                 (fitFloat <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
-            }
-     
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-     
-            myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
-            return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat);
-         } else {
-            return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] +
-                          robustHistogram->bounds->data.F32[binNum]));
-        }
-    */
+    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
+    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
+    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
+    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
+
+    if ((binNum > 0) &&
+            (binNum < (histogram->nums->n+1))) {
+        x->data.F64[0] = (double) 0.5 *
+                         (histogram->bounds->data.F32[binNum-1] +
+                          histogram->bounds->data.F32[binNum]);
+        x->data.F64[1] = (double) 0.5 *
+                         (histogram->bounds->data.F32[binNum] +
+                          histogram->bounds->data.F32[binNum+1]);
+        x->data.F64[2] = (double) 0.5 *
+                         (histogram->bounds->data.F32[binNum+1] +
+                          histogram->bounds->data.F32[binNum+2]);
+
+        y->data.F64[0] = cumulativeSums->data.F32[binNum-1];
+        y->data.F64[1] = cumulativeSums->data.F32[binNum];
+        y->data.F64[2] = cumulativeSums->data.F32[binNum+1];
+
+        if (!((y->data.F64[0] <= fitFloat) &&
+                (fitFloat <= y->data.F64[2]))) {
+            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
+        }
+
+        yErr->data.F64[0] = 1.0;
+        yErr->data.F64[1] = 1.0;
+        yErr->data.F64[2] = 1.0;
+
+        myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
+        return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2],
+                                fitFloat));
+    } else {
+        return(0.5 * (histogram->bounds->data.F32[binNum+1] +
+                      histogram->bounds->data.F32[binNum]));
+    }
+
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+    psFree(myPoly);
     return(0.0);
 }
@@ -1049 +1072 @@
 p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
 calculates a variety of robust stat measures:
-PS_STAT_ROBUST_MEAN
-PS_STAT_ROBUST_MEDIAN
-PS_STAT_ROBUST_MODE
-PS_STAT_ROBUST_STDEV
-PS_STAT_ROBUST_QUARTILE
+    PS_STAT_ROBUST_MEAN
+    PS_STAT_ROBUST_MEDIAN
+    PS_STAT_ROBUST_MODE
+    PS_STAT_ROBUST_STDEV
+    PS_STAT_ROBUST_QUARTILE
 I have included all that computation in a single function, as opposed to
-breaking it across several functions for one primary reason: 
-they all
+breaking it across several functions for one primary reason:  they all
 require the same basic initial processing steps (calculate the histogram,
-etc.) 
-however there is no currently defined means, in the SDRS, to
-specify this computation as a separate step.  If the above robust stat
-measures were calcualted in separate functiosn, then much of the initial
-processing would be duplicated.
+etc.)  however there is no currently defined means, in the SDRS, to specify
+this computation as a separate step.  If the above robust stat measures were
+calcualted in separate functiosn, then much of the initial processing would
+be duplicated.
+ 
 Inputs
-myVector
-maskVector
-maskVal
-stats
+    myVector
+    maskVector
+    maskVal
+    stats
 Returns
-NULL
+    NULL
 *****************************************************************************/
 void p_psVectorRobustStats(const psVector *restrict myVector,
@@ -1224 +1246 @@
         sumN50+= (float) robustHistogram->nums->data.S32[i];
     }
+
     // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
     dL = (UQBinNum - LQBinNum) / 4;
@@ -1244 +1267 @@
                 robustHistogramVector->data.F32[i];
             cumulativeMedian+= robustHistogramVector->data.F32[i];
-            sumNfit+= (float) robustHistogramVector->data.S32[i];
+            sumNfit+= (float) robustHistogram->nums->data.S32[i];
         }
     }
@@ -1780 +1803 @@
     }
 
+    // XXX: Should we abort if (stats->min == stats->max)?
+    if (stats->min >= stats->max) {
+        psAbort(__func__, "psVectorStats() called with range: %f to %f\n",
+                stats->min, stats->max);
+    }
+
+
+
     PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
     if (mask != NULL) {