X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsupport%2Fmathfuncs.cpp;h=c9b18ac3da75ca61cf6c10e82bfe5689a4bca0a7;hp=2eaa55153c6f73c6325dd66e4335dad3c0fcd912;hb=9b2bb510160bdb56f04847f5b55ab61dd8a47976;hpb=f7d2b7144f32a7bd157b7689022e62944b82fcc1

diff --git a/libctsupport/mathfuncs.cpp b/libctsupport/mathfuncs.cpp
index 2eaa551..c9b18ac 100644
--- a/libctsupport/mathfuncs.cpp
+++ b/libctsupport/mathfuncs.cpp
@@ -2,7 +2,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: mathfuncs.cpp,v 1.4 2000/12/21 03:40:58 kevin Exp $
+**  $Id: mathfuncs.cpp,v 1.6 2001/01/02 05:34:57 kevin Exp $
 **
 **  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
@@ -22,18 +22,18 @@
 
 
 /* NAME
- *    integrateSimpson         Integrate array of data by Simpson's rule
- *
- * 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)
- *
- * RETURNS
- *    integrand of function
- */
+*    integrateSimpson         Integrate array of data by Simpson's rule
+*
+* 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)
+*
+* RETURNS
+*    integrand of function
+*/
 
 double 
 integrateSimpson (const double xmin, const double xmax, const double *y, const int np)	
@@ -42,16 +42,16 @@ integrateSimpson (const double xmin, const double xmax, const double *y, const i
     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
-
+  
   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]);
   }
   
@@ -63,13 +63,13 @@ integrateSimpson (const double xmin, const double xmax, const double *y, const i
 
 
 /* NAME
- *    normalizeAngle       Normalize angle to 0 to 2 * PI range
- *
- * SYNOPSIS
- *    t = normalizeAngle (theta)
- *    double t	       Normalized angle
- *    double theta     Input angle
- */
+*    normalizeAngle       Normalize angle to 0 to 2 * PI range
+*
+* SYNOPSIS
+*    t = normalizeAngle (theta)
+*    double t	       Normalized angle
+*    double theta     Input angle
+*/
 
 double 
 normalizeAngle (double theta)
@@ -86,52 +86,52 @@ normalizeAngle (double theta)
 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];
-	int i;
-    for (i = 0; i < nPoints; i++) {
-		double v = vec[i];
-		if (v > max)
-			max = v;
-		if (v < min)
-			min = v;
-		mean += v;
-    }
-    mean /= nPoints;
-	
-    static const int nbin = 1024;
-    int hist[ nbin ] = {0};
-    double spread = max - min;
-    mode = 0;
-    stddev = 0;
-    for (i = 0; i < nPoints; i++) {
-		double v = vec[i];
-		int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
-		hist[b]++;
-		double diff = (v - mean);
-		stddev += diff * diff;
-    }
-    stddev = sqrt (stddev / nPoints);
-	
-    int max_binindex = 0;
-    int max_bin = -1;
-    for (int ibin = 0; ibin < nbin; ibin++) {
-		if (hist[ibin] > max_bin) {
-			max_bin = hist[ibin];
-			max_binindex = ibin;
-		}
+  if (nPoints <= 0)
+    return;
+  
+  mean = 0;
+  min = vec[0];
+  max = vec[0];
+  int i;
+  for (i = 0; i < nPoints; i++) {
+    double v = vec[i];
+    if (v > max)
+      max = v;
+    if (v < min)
+      min = v;
+    mean += v;
+  }
+  mean /= nPoints;
+  
+  static const int nbin = 1024;
+  int hist[ nbin ] = {0};
+  double spread = max - min;
+  mode = 0;
+  stddev = 0;
+  for (i = 0; i < nPoints; i++) {
+    double v = vec[i];
+    int b = static_cast<int>((((v - min) / spread) * (nbin - 1)) + 0.5);
+    hist[b]++;
+    double diff = (v - mean);
+    stddev += diff * diff;
+  }
+  stddev = sqrt (stddev / nPoints);
+  
+  int max_binindex = 0;
+  int max_bin = -1;
+  for (int ibin = 0; ibin < nbin; ibin++) {
+    if (hist[ibin] > max_bin) {
+      max_bin = hist[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;
+  }
+  
+  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;
 }