1 /*****************************************************************************
2 ** This is part of the CTSim program
3 ** Copyright (c) 1983-2009 Kevin Rosenberg
5 ** This program is free software; you can redistribute it and/or modify
6 ** it under the terms of the GNU General Public License (version 2) as
7 ** published by the Free Software Foundation.
9 ** This program is distributed in the hope that it will be useful,
10 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ** GNU General Public License for more details.
14 ** You should have received a copy of the GNU General Public License
15 ** along with this program; if not, write to the Free Software
16 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 ******************************************************************************/
19 #include "ctsupport.h"
23 * integrateSimpson Integrate array of data by Simpson's rule
26 * double integrateSimpson (xmin, xmax, y, np)
27 * double xmin, xmax Extent of integration
28 * double y[] Function values to be integrated
29 * int np number of data points
30 * (must be an odd number and at least 3)
33 * integrand of function
37 integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
42 return ((xmax - xmin) * (y[0] + y[1]) / 2);
45 int nDiv = (np - 1) / 2; // number of divisions
46 double width = (xmax - xmin) / (double) (np - 1); // width of cells
48 for (int i = 1; i <= nDiv; i++) {
50 int xl = xr - 2; // 2 * (i - 1) == 2 * i - 2 == xr - 2
51 int xm = xr - 1; // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
53 area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
56 if ((np & 1) == 0) /* do last trapazoid */
57 area += width * (y[np-2] + y[np-1]) / 2;
64 * normalizeAngle Normalize angle to 0 to 2 * PI range
67 * t = normalizeAngle (theta)
68 * double t Normalized angle
69 * double theta Input angle
73 normalizeAngle (double theta)
77 while (theta >= TWOPI)
85 vectorNumericStatistics (std::vector<double> vec, const int nPoints, double& min, double& max, double& mean, double& mode, double& median, double& stddev)
94 for (i = 0; i < nPoints; i++) {
104 static const int nbin = 1024;
105 int hist[ nbin ] = {0};
106 double spread = max - min;
109 for (i = 0; i < nPoints; i++) {
111 int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
113 double diff = (v - mean);
114 stddev += diff * diff;
116 stddev = sqrt (stddev / nPoints);
118 int max_binindex = 0;
120 for (int ibin = 0; ibin < nbin; ibin++) {
121 if (hist[ibin] > max_bin) {
122 max_bin = hist[ibin];
127 mode = (max_binindex * spread / (nbin - 1)) + min;
129 std::sort(vec.begin(), vec.end());
131 if (nPoints % 2) // Odd
132 median = vec[((nPoints - 1) / 2)];
134 median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;