/*****************************************************************************
** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
+** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: mathfuncs.cpp,v 1.7 2001/01/02 16:02:13 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
*
* 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)
+* 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)
+double
+integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
{
if (np < 2)
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
-
+ 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]);
}
-
- 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);
}
*
* SYNOPSIS
* t = normalizeAngle (theta)
-* double t Normalized angle
+* double t Normalized angle
* double theta Input angle
*/
-double
+double
normalizeAngle (double theta)
{
while (theta < 0.)
theta += TWOPI;
while (theta >= TWOPI)
theta -= TWOPI;
-
+
return (theta);
}
-void
+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];
mean += v;
}
mean /= nPoints;
-
+
static const int nbin = 1024;
int hist[ nbin ] = {0};
double spread = max - min;
stddev += diff * diff;
}
stddev = sqrt (stddev / nPoints);
-
+
int max_binindex = 0;
int max_bin = -1;
for (int ibin = 0; ibin < nbin; 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;
}
+
+