X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsupport%2Fsimpson.cpp;fp=libctsupport%2Fsimpson.cpp;h=0000000000000000000000000000000000000000;hp=e7a52706152d855b96f9d36b5be2ce5d11bb86bf;hb=2d39e823ba389fc68e5317c422b55be006094252;hpb=a95e41ac40cd2f3a4401d921618604cf33f2a904 diff --git a/libctsupport/simpson.cpp b/libctsupport/simpson.cpp deleted file mode 100644 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); -}