#include <iostream.h>
#include <stdlib.h>
#include <math.h>

using namespace std;
/*     ***************************************************************** */
/*     INTEGRATION BY NEWTON-COTES RULE (N=6) */

/*     A, B = ARE THE INTEGRATION LIMITS */
/*     SUM  = THE RESULT OF THE INTEGRATION */
/*     THE FUNCTION TO BE INTEGRATED MUST BE DEFINED IN */
/*     A SUBPROGRAM FUNCTION FUN(ARGUMENTS) */

/*     SPECIFICALLY, THIS PROGRAM EVALUATES BESSEL FUNCTION */
/*     JN(X), OF ORDER N AND ARGUMENT X */
/*     ***************************************************************** */

const double pi = 3.1415927;
// method declarations
void integration(double *aa, double *bb, int *no, double *x, double *sum);
double fun(double *theta, int *n, double *x);

int main()
{
    /* Local variables */
    double a, b, x, sum;
    int i, n;

    a = 0.0;
    b = pi;
    n = 4;
    for (i = 1; i <= 20; i++) {
        x = static_cast<double>(i) * 0.1;
        integration(&a, &b, &n, &x, &sum);
        cout << "x = " << x << '\t' << "j0 = " << sum << endl;
    }

    system("PAUSE");
    return 0;
}

/*     ***************************************************************** */
/*     SUBROUTINE FOR INTEGRATION USING NEWTON-COTES RULE (N=6) */
/*     A IS THE LOWER LIMIT OF INTEGRATION */
/*     B IS THE UPPER LIMIT OF INTEGRATION */
/*     H IS THE INTERVAL SIZE */
/*     ***************************************************************** */
void integration(double *aa, double *bb, int *no, double *x, double *sum)
{
    /* Initialized data */
    double c[7] = { 41.0,216.0,27.0,272.0,27.0,216.0,41.0 };
    int nc = 840;
    /* Local variables */
    double a, h, t;
    int i, j, m, n, nn;

/*     START COMPUTATION */
    n = 6;
    m = 240;
/* M MUST BE A MULTIPLE OF N */
    h = (*bb - *aa) / static_cast<double>(m);
    *sum = 0.0;
    nn = m / n;
    t = *aa;
    for (i = 1; i <= nn; i++) {
	    for (j = 0; j <= n; j++) {
	        a = c[j] * fun(&t, no, x);
	        *sum += a;
            t += h;
        }
	    t -= h;
    }
    *sum *= h * static_cast<double>(n) / static_cast<double>(nc);
} /* integration */

/*     ***************************************************************** */
/*     THE FUNCTION TO BE INTEGRATED */
/*     ***************************************************************** */
double fun(double *theta, int *n, double *x)
{
    double ret_val;

    ret_val = cos(*x * sin(*theta) - static_cast<double>(*n) * *theta) / pi;
    return ret_val;
} /* fun */