/*****************************************************************************
** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
+** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: mathfuncs.cpp,v 1.3 2000/12/20 20:08:48 kevin Exp $
+** $Id$
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
/* NAME
- * integrateSimpson Integrate array of data by Simpson's rule
- *
- * SYNOPSIS
- * double integrateSimpson (xmin, xmax, y, np)
- * double xmin, xmax Extent of integration
- * double y[] Function values to be integrated
- * int np number of data points
- * (must be an odd number and at least 3)
- *
- * RETURNS
- * integrand of function
- */
-
-double
-integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
+* integrateSimpson Integrate array of data by Simpson's rule
+*
+* SYNOPSIS
+* double integrateSimpson (xmin, xmax, y, np)
+* double xmin, xmax Extent of integration
+* double y[] Function values to be integrated
+* int np number of data points
+* (must be an odd number and at least 3)
+*
+* RETURNS
+* integrand of function
+*/
+
+double
+integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
{
if (np < 2)
return (0.);
double area = 0;
int nDiv = (np - 1) / 2; // number of divisions
- double width = (xmax - xmin) / (double) (np - 1); // width of cells
+ double width = (xmax - xmin) / (double) (np - 1); // width of cells
for (int i = 1; i <= nDiv; i++) {
int xr = 2 * i;
area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
}
-
- if ((np & 1) == 0) /* do last trapazoid */
+
+ if ((np & 1) == 0) /* do last trapazoid */
area += width * (y[np-2] + y[np-1]) / 2;
-
+
return (area);
}
/* NAME
- * normalizeAngle Normalize angle to 0 to 2 * PI range
- *
- * SYNOPSIS
- * t = normalizeAngle (theta)
- * double t Normalized angle
- * double theta Input angle
- */
-
-double
+* normalizeAngle Normalize angle to 0 to 2 * PI range
+*
+* SYNOPSIS
+* t = normalizeAngle (theta)
+* double t Normalized angle
+* double theta Input angle
+*/
+
+double
normalizeAngle (double theta)
{
while (theta < 0.)
theta += TWOPI;
while (theta >= TWOPI)
theta -= TWOPI;
-
+
return (theta);
}
-\r
-\r
-void \r
-vectorNumericStatistics (std::vector<double> vec, double& min, double& max, double& mean, double& mode, double& median, double& stddev)\r
-{\r
- int n = vec.size();\r
- if (n <= 0)\r
- return;\r
-\r
- mean = 0;\r
- min = vec[0];\r
- max = vec[0];\r
- int i;\r
- for (i = 0; i < n; i++) {\r
- double v = vec[i];\r
- if (v > max)\r
- max = v;\r
- if (v < min)\r
- min = v;\r
- mean += v;\r
- }\r
- mean /= n;\r
- \r
- static const int nbin = 1024;\r
- int hist[ nbin ] = {0};\r
- double spread = max - min;\r
- mode = 0;\r
- stddev = 0;\r
- for (i = 0; i < n; i++) {\r
- double v = vec[i];\r
- int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);\r
- hist[b]++;\r
- double diff = (v - mean);\r
- stddev += diff * diff;\r
- }\r
- stddev = sqrt (stddev / n);\r
- \r
- int max_binindex = 0;\r
- int max_bin = -1;\r
- for (int ibin = 0; ibin < nbin; ibin++) {\r
- if (hist[ibin] > max_bin) {\r
- max_bin = hist[ibin];\r
- max_binindex = ibin;\r
- }\r
- }\r
- \r
- mode = (max_binindex * spread / (nbin - 1)) + min;\r
- \r
- std::sort(vec.begin(), vec.end());\r
- \r
- if (n % 2) // Odd\r
- median = vec[((n - 1) / 2)];\r
- else // Even\r
- median = (vec[ (n / 2) - 1 ] + vec[ n / 2 ]) / 2;\r
-}\r
+
+
+void
+vectorNumericStatistics (std::vector<double> vec, const int nPoints, double& min, double& max, double& mean, double& mode, double& median, double& stddev)
+{
+ if (nPoints <= 0)
+ return;
+
+ mean = 0;
+ min = vec[0];
+ max = vec[0];
+ int i;
+ for (i = 0; i < nPoints; i++) {
+ double v = vec[i];
+ if (v > max)
+ max = v;
+ if (v < min)
+ min = v;
+ mean += v;
+ }
+ mean /= nPoints;
+
+ static const int nbin = 1024;
+ int hist[ nbin ] = {0};
+ double spread = max - min;
+ mode = 0;
+ stddev = 0;
+ for (i = 0; i < nPoints; i++) {
+ double v = vec[i];
+ int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
+ hist[b]++;
+ double diff = (v - mean);
+ stddev += diff * diff;
+ }
+ stddev = sqrt (stddev / nPoints);
+
+ int max_binindex = 0;
+ int max_bin = -1;
+ for (int ibin = 0; ibin < nbin; ibin++) {
+ if (hist[ibin] > max_bin) {
+ max_bin = hist[ibin];
+ max_binindex = ibin;
+ }
+ }
+
+ mode = (max_binindex * spread / (nbin - 1)) + min;
+
+ std::sort(vec.begin(), vec.end());
+
+ if (nPoints % 2) // Odd
+ median = vec[((nPoints - 1) / 2)];
+ else // Even
+ median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;
+}
+
+