2ff4e0bcff85645ddb6c375b9dccb0bd57f78202
1 /*****************************************************************************
2 **  This is part of the CTSim program
3 **  Copyright (c) 1983-2001 Kevin Rosenberg
4 **
5 **  \$Id: mathfuncs.cpp,v 1.8 2001/01/28 19:10:18 kevin Exp \$
6 **
7 **  This program is free software; you can redistribute it and/or modify
8 **  it under the terms of the GNU General Public License (version 2) as
9 **  published by the Free Software Foundation.
10 **
11 **  This program is distributed in the hope that it will be useful,
12 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 **  GNU General Public License for more details.
15 **
16 **  You should have received a copy of the GNU General Public License
17 **  along with this program; if not, write to the Free Software
18 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
19 ******************************************************************************/
21 #include "ctsupport.h"
24 /* NAME
25 *    integrateSimpson         Integrate array of data by Simpson's rule
26 *
27 * SYNOPSIS
28 *    double integrateSimpson (xmin, xmax, y, np)
29 *    double xmin, xmax          Extent of integration
30 *    double y[]         Function values to be integrated
31 *    int np                     number of data points
32 *                               (must be an odd number and at least 3)
33 *
34 * RETURNS
35 *    integrand of function
36 */
38 double
39 integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
40 {
41   if (np < 2)
42     return (0.);
43   else if (np == 2)
44     return ((xmax - xmin) * (y + y) / 2);
46   double area = 0;
47   int nDiv = (np - 1) / 2;  // number of divisions
48   double width = (xmax - xmin) / (double) (np - 1);     // width of cells
50   for (int i = 1; i <= nDiv; i++) {
51     int xr = 2 * i;
52     int xl = xr - 2;       // 2 * (i - 1) == 2 * i - 2 == xr - 2
53     int xm = xr - 1;       // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
55     area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
56   }
58   if ((np & 1) == 0)            /* do last trapazoid */
59     area += width * (y[np-2] + y[np-1]) / 2;
61   return (area);
62 }
65 /* NAME
66 *    normalizeAngle       Normalize angle to 0 to 2 * PI range
67 *
68 * SYNOPSIS
69 *    t = normalizeAngle (theta)
70 *    double t          Normalized angle
71 *    double theta     Input angle
72 */
74 double
75 normalizeAngle (double theta)
76 {
77   while (theta < 0.)
78     theta += TWOPI;
79   while (theta >= TWOPI)
80     theta -= TWOPI;
82   return (theta);
83 }
86 void
87 vectorNumericStatistics (std::vector<double> vec, const int nPoints, double& min, double& max, double& mean, double& mode, double& median, double& stddev)
88 {
89   if (nPoints <= 0)
90     return;
92   mean = 0;
93   min = vec;
94   max = vec;
95   int i;
96   for (i = 0; i < nPoints; i++) {
97     double v = vec[i];
98     if (v > max)
99       max = v;
100     if (v < min)
101       min = v;
102     mean += v;
103   }
104   mean /= nPoints;
106   static const int nbin = 1024;
107   int hist[ nbin ] = {0};
108   double spread = max - min;
109   mode = 0;
110   stddev = 0;
111   for (i = 0; i < nPoints; i++) {
112     double v = vec[i];
113     int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
114     hist[b]++;
115     double diff = (v - mean);
116     stddev += diff * diff;
117   }
118   stddev = sqrt (stddev / nPoints);
120   int max_binindex = 0;
121   int max_bin = -1;
122   for (int ibin = 0; ibin < nbin; ibin++) {
123     if (hist[ibin] > max_bin) {
124       max_bin = hist[ibin];
125       max_binindex = ibin;
126     }
127   }
129   mode = (max_binindex * spread / (nbin - 1)) + min;
131   std::sort(vec.begin(), vec.end());
133   if (nPoints % 2)  // Odd
134     median = vec[((nPoints - 1) / 2)];
135   else        // Even
136     median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;
137 }