c07445253d3355a5f92df91b67842f7b44d6d904
[ctsim.git] / libctsupport / mathfuncs.cpp
1 /*****************************************************************************
2 **  This is part of the CTSim program
3 **  Copyright (C) 1983-2000 Kevin Rosenberg
4 **
5 **  $Id: mathfuncs.cpp,v 1.3 2000/12/20 20:08:48 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 ******************************************************************************/
20
21 #include "ctsupport.h"
22
23
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  */
37
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[0] + y[1]) / 2);
45
46   double area = 0;
47   int nDiv = (np - 1) / 2;  // number of divisions
48   double width = (xmax - xmin) / (double) (np - 1);     // width of cells
49
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
54
55     area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
56   }
57   
58   if ((np & 1) == 0)            /* do last trapazoid */
59     area += width * (y[np-2] + y[np-1]) / 2;
60   
61   return (area);
62 }
63
64
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  */
73
74 double 
75 normalizeAngle (double theta)
76 {
77   while (theta < 0.)
78     theta += TWOPI;
79   while (theta >= TWOPI)
80     theta -= TWOPI;
81   
82   return (theta);
83 }
84 \r
85 \r
86 void \r
87 vectorNumericStatistics (std::vector<double> vec, double& min, double& max, double& mean, double& mode, double& median, double& stddev)\r
88 {\r
89         int n = vec.size();\r
90         if (n <= 0)\r
91                 return;\r
92 \r
93     mean = 0;\r
94     min = vec[0];\r
95     max = vec[0];\r
96         int i;\r
97     for (i = 0; i < n; i++) {\r
98                 double v = vec[i];\r
99                 if (v > max)\r
100                         max = v;\r
101                 if (v < min)\r
102                         min = v;\r
103                 mean += v;\r
104     }\r
105     mean /= n;\r
106         \r
107     static const int nbin = 1024;\r
108     int hist[ nbin ] = {0};\r
109     double spread = max - min;\r
110     mode = 0;\r
111     stddev = 0;\r
112     for (i = 0; i < n; i++) {\r
113                 double v = vec[i];\r
114                 int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);\r
115                 hist[b]++;\r
116                 double diff = (v - mean);\r
117                 stddev += diff * diff;\r
118     }\r
119     stddev = sqrt (stddev / n);\r
120         \r
121     int max_binindex = 0;\r
122     int max_bin = -1;\r
123     for (int ibin = 0; ibin < nbin; ibin++) {\r
124                 if (hist[ibin] > max_bin) {\r
125                         max_bin = hist[ibin];\r
126                         max_binindex = ibin;\r
127                 }\r
128     }\r
129         \r
130     mode = (max_binindex * spread / (nbin - 1)) + min;\r
131         \r
132         std::sort(vec.begin(), vec.end());\r
133         \r
134         if (n % 2)  // Odd\r
135                 median = vec[((n - 1) / 2)];\r
136         else        // Even\r
137                 median = (vec[ (n / 2) - 1 ] + vec[ n / 2 ]) / 2;\r
138 }\r