r117: *** empty log message ***

[ctsim.git] / libctsupport / simpson.cpp

diff --git a/libctsupport/simpson.cpp b/libctsupport/simpson.cpp
deleted file mode 100644 (file)
index e7a5270..0000000
--- a/libctsupport/simpson.cpp
+++ /dev/null
@@ -1,62 +0,0 @@
-/*****************************************************************************
-**  This is part of the CTSim program
-**  Copyright (C) 1983-2000 Kevin Rosenberg
-**
-**  $Id: simpson.cpp,v 1.2 2000/06/19 19:04:05 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);
-}