** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: mathfuncs.cpp,v 1.2 2000/08/31 08:38:58 kevin Exp $
+** $Id: mathfuncs.cpp,v 1.6 2001/01/02 05:34:57 kevin Exp $
**
** 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
- */
+* 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)
return (0.);
else if (np == 2)
return ((xmax - xmin) * (y[0] + y[1]) / 2);
-
+
double area = 0;
int nDiv = (np - 1) / 2; // number of divisions
double width = (xmax - xmin) / (double) (np - 1); // width of cells
-
+
for (int i = 1; i <= nDiv; i++) {
int xr = 2 * i;
int xl = xr - 2; // 2 * (i - 1) == 2 * i - 2 == xr - 2
int xm = xr - 1; // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
-
+
area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
}
/* NAME
- * normalizeAngle Normalize angle to 0 to 2 * PI range
- *
- * SYNOPSIS
- * t = normalizeAngle (theta)
- * double t Normalized angle
- * double theta Input angle
- */
+* 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)
return (theta);
}
+\r
+\r
+void \r
+vectorNumericStatistics (std::vector<double> vec, const int nPoints, double& min, double& max, double& mean, double& mode, double& median, double& stddev)\r
+{\r
+ if (nPoints <= 0)\r
+ return;\r
+ \r
+ mean = 0;\r
+ min = vec[0];\r
+ max = vec[0];\r
+ int i;\r
+ for (i = 0; i < nPoints; 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 /= nPoints;\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 < nPoints; 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 / nPoints);\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 (nPoints % 2) // Odd\r
+ median = vec[((nPoints - 1) / 2)];\r
+ else // Even\r
+ median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;\r
+}\r