1 /*****************************************************************************
2 ** This is part of the CTSim program
3 ** Copyright (C) 1983-2000 Kevin Rosenberg
5 ** $Id: mathfuncs.cpp,v 1.2 2000/08/31 08:38:58 kevin Exp $
7 ** This program is free software; you can redistribute it and/or modify
8 ** it under the terms of the GNU General Public License (version 2) as
9 ** published by the Free Software Foundation.
11 ** This program is distributed in the hope that it will be useful,
12 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 ** GNU General Public License for more details.
16 ** You should have received a copy of the GNU General Public License
17 ** along with this program; if not, write to the Free Software
18 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 ******************************************************************************/
21 #include "ctsupport.h"
25 * integrateSimpson Integrate array of data by Simpson's rule
28 * double integrateSimpson (xmin, xmax, y, np)
29 * double xmin, xmax Extent of integration
30 * double y[] Function values to be integrated
31 * int np number of data points
32 * (must be an odd number and at least 3)
35 * integrand of function
39 integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
44 return ((xmax - xmin) * (y[0] + y[1]) / 2);
47 int nDiv = (np - 1) / 2; // number of divisions
48 double width = (xmax - xmin) / (double) (np - 1); // width of cells
50 for (int i = 1; i <= nDiv; i++) {
52 int xl = xr - 2; // 2 * (i - 1) == 2 * i - 2 == xr - 2
53 int xm = xr - 1; // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
55 area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
58 if ((np & 1) == 0) /* do last trapazoid */
59 area += width * (y[np-2] + y[np-1]) / 2;
66 * normalizeAngle Normalize angle to 0 to 2 * PI range
69 * t = normalizeAngle (theta)
70 * double t Normalized angle
71 * double theta Input angle
75 normalizeAngle (double theta)
79 while (theta >= TWOPI)