r117: *** empty log message ***

[ctsim.git] / libctsupport / mathfuncs.cpp

diff --git a/libctsupport/mathfuncs.cpp b/libctsupport/mathfuncs.cpp
new file mode 100644 (file)
index 0000000..bc3b383
--- /dev/null
+++ b/libctsupport/mathfuncs.cpp
@@ -0,0 +1,83 @@
+/*****************************************************************************
+**  This is part of the CTSim program
+**  Copyright (C) 1983-2000 Kevin Rosenberg
+**
+**  $Id: mathfuncs.cpp,v 1.1 2000/06/22 10:17:28 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
+**  published by the Free Software Foundation.
+**
+**  This program is distributed in the hope that it will be useful,
+**  but WITHOUT ANY WARRANTY; without even the implied warranty of
+**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+**  GNU General Public License for more details.
+**
+**  You should have received a copy of the GNU General Public License
+**  along with this program; if not, write to the Free Software
+**  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+******************************************************************************/
+
+#include "ctsupport.h"
+
+
+/* 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
+ */
+
+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
+
+  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 */
+    area += width * (y[np-2] + y[np-1]) / 2;
+  
+  return (area);
+}
+
+
+/* NAME
+ *    norm_angle       Normalize angle to 0 to 2 * PI range
+ *
+ * SYNOPSIS
+ *    t = norm_angle (theta)
+ *    double t        Normalized angle
+ *    double theta     Input angle
+ */
+
+double 
+normalizeAngle (double theta)
+{
+  while (theta < 0.)
+    theta += TWOPI;
+  while (theta >= TWOPI)
+    theta -= TWOPI;
+  
+  return (theta);
+}