/*---------------------------------------------

Simpson integration technique for

evaluating  double integrals

-----------------------------------------------*/

 #include "math.h"

float fxy[50][50],fy[50];

float pi;

main() {

float f();

float llx, lly, ulx, uly,x,y;

float hx, hy ,ef,of,simpson;

int nosx, nosy,i,j;

/*----------------------------------------------

Simpson integration constants.

            nos -> number of strips

            ul  -> upper limit of integration

            ll  -> lower limit of integration

        h   -> incremental value per strip

-------------------------------------------------*/

            pi = 3.1415926;

            nosx = 30;  nosy = 40;

            llx = 0.1e-8;lly=0.1e-8;

            ulx = pi; uly=2.0*pi;

            hx = (ulx-llx)/nosx;

            hy = (uly-lly)/nosy;

            printf ("\n\nDouble integration parameters: \n");

            printf (" Step size hx and hy:  %f,%f\n", hx,hy);

            printf (" Number of steps (xy,y): %d,%d\n\n",nosx,nosy);

/*

            Calculate all the points within the integration domain

*/

              for (j=0; j <= nosy; j++) {

                       y = j*hy + lly;

                       for (i=0; i<= nosx; i++) {

                             x = i*hx + llx;

                             fxy[i][j] = f(x,y);

                       }

            }

/*

Now perform a Simpson integration along

the x axis and accumulate results using

the y axis variable as an index

*/

            for (j=0; j <= nosy; j++) {

                        of = fxy[1][j];

                        ef=0.0;

                        for (i=2; i <= nosx-2; i +=2) {

                                    ef += fxy[i][j];

                                    of += fxy[i+1][j];

                        }

                        fy[j]=fxy[0][j] + fxy[nosx][j]+2.0*ef+4*of;

            }

/*Lastly perform Simpson integration

along the y axis.

*/

            of=fy[1];

            ef=0.0;

            for (j=2; j <= nosy-2; j+=2) {

                        ef += fy[j];

                        of += fy[j+1];

            }

            simpson = (hx*hy/9.0)*(fy[0]+fy[nosy]+2.0*ef+4.0*of);

            printf ("Result = %f\n\n", simpson);

}         

/*Enter the function to be integrated here.

*/

float f(x,y) float x,y; {

double zd;

            zd = cos(pi*cos(x)/2.0)*cos(pi*(1.0-sin(x)*cos(y))/4.0);

            zd = pow(zd,2.0)/sin(x);

            return(fabs(zd));

}