X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsupport%2Fmathfuncs.cpp;h=4a96e8f617fd1aa55bbfa77de19698ffa8057d5b;hp=289c8b7f9f27589920c42435c6d50e79b8a4a77e;hb=f7ee98f7d964ed361068179f0e7ea4475ed1abdf;hpb=8a7697ce57b56cdc43698cd1241ad98d49f9b5ac diff --git a/libctsupport/mathfuncs.cpp b/libctsupport/mathfuncs.cpp index 289c8b7..4a96e8f 100644 --- a/libctsupport/mathfuncs.cpp +++ b/libctsupport/mathfuncs.cpp @@ -26,38 +26,38 @@ * * 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); } @@ -67,28 +67,28 @@ integrateSimpson (const double xmin, const double xmax, const double *y, const i * * 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 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]; @@ -102,7 +102,7 @@ vectorNumericStatistics (std::vector vec, const int nPoints, double& min mean += v; } mean /= nPoints; - + static const int nbin = 1024; int hist[ nbin ] = {0}; double spread = max - min; @@ -116,7 +116,7 @@ vectorNumericStatistics (std::vector vec, const int nPoints, double& min stddev += diff * diff; } stddev = sqrt (stddev / nPoints); - + int max_binindex = 0; int max_bin = -1; for (int ibin = 0; ibin < nbin; ibin++) { @@ -125,11 +125,11 @@ vectorNumericStatistics (std::vector vec, const int nPoints, double& min 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