/*
 *  ============================================================================
 *  Use Gauss-Legendre algorithm to compute locations and weighting coefficients
 *  for an arbitrary number of integration points (limit to 10 in test program).
 * 
 *  Written By : Mark Austin                                           July 1996 
 *  ============================================================================
 */

#include <stdio.h>
#include <math.h>

enum { False = 0, True = 1};

void GaussLegendre( int , double daX[], double daW[], double , double );

/*
 *  ============================================================================
 *  main() : Prompt user for lower and upper bounds on numerical integration,
 *           and the number of gauss quadrature points to be computed.
 *         : Call GaussLegendre() to compute quadrature points with method of
 *           Gauss-Legendre Integration.
 *  ============================================================================
 */

int main( void ) {
enum { MaxIntegrationPoints = 10};
double daX [ MaxIntegrationPoints ];  /* Location of integration points   */
double daW [ MaxIntegrationPoints ];  /* Weightings of integration points */
char caBuffer[8];   /* Buffer array to hold input from screen             */
int    iContinue;   /* Flag used continue/discontinue input               */
int    iNoPoints;   /* Number of Integration Points                       */
double  dMinValue, dMaxValue;   /* Lower and upper bound of integration [ dMinValue, dMaxValue ]  */
int           ii;

   /* [a] : Prompt user for desired number of integration points */

   iContinue = True;
   while( iContinue == (int) True ) {

     /* [a.1] : Ask user if input will be provided */

     printf("\n");
     printf("Would you like compute integration coefficients (yes/no/quit) ? ");
     scanf("%s%*c", caBuffer);

     if(strcmp(caBuffer,"yes") != 0)
        iContinue = False;
     else {

        /* [a.2] : Prompt User for number of integration points */

        printf("Enter lower bound of integration : ");
        scanf("%lf%*c", &dMinValue );
        printf("Enter upper bound of integration : ");
        scanf("%lf%*c", &dMaxValue );
        printf("Enter no integration points : ");
        scanf("%i%*c", &iNoPoints);

        /* [a.3] : Print summary of user input to screen */

        printf("Bounds of integration are : [ %f, %f ]\n", dMinValue, dMaxValue );
        printf("No of integration points is : %d\n", iNoPoints);

        /* [a.4] : Compute and print integration point locations and 
                   weighting coefficients */

        if (iNoPoints <= MaxIntegrationPoints ) {

            GaussLegendre( iNoPoints, daX, daW, dMinValue, dMaxValue );

            printf("\nNo Integration Points = %d\n", iNoPoints );
            for( ii = 1; ii <= iNoPoints ; ii = ii + 1 ) {
                 printf("X[ %2d ] = %10f : ", ii, daX[ii-1] );
                 printf("W[ %2d ] = %10f \n", ii, daW[ii-1] );
            }
        } else
            printf("\nRUN-TIME ERROR : Maximum allowable integration points is %d \n",
                   MaxIntegrationPoints );
     }
  }
}

/*
 *  =========================================================================
 *  GaussLegendre() : Given an interval [x1,x2] and the number of integration
 *  points n, this routine computes: 
 *
 *     daX [] = array containing locations of "n" integration points.
 *     daW [] = array containing weights for "n" integration points.
 *  
 *  Input :  iNoPoints : number of integration points (no limitation on n).
 *           daX[]     : array of local coordinates of integration points.
 *           daW[]     : array of weighting coefficients for integration.
 *           dMinValue : lower bound on integration.
 *           dMaxValue : upper bound on integration.
 *  Output : void :
 *
 *  Note that [ dMinValue1,dMaxValue ] need not be [-1,1]. The Gauss-Legendre
 *  algorithm is general purpose, and there is no limitation on the desired
 *  number of integration points n.
 *  =========================================================================
 */

#define eps 1.5e-16

void
GaussLegendre(int iNoPoints, double daX[], double daW[],
                  double dMinValue, double dMaxValue ) {
double dXmean, dXrange, dP1, dP2, dP3, dPp, dZ, dZ1;
int  ii, ij, im; 

   /* [a] : Roots are symmetric so only need to find half of them.   */

   im  = (iNoPoints + 1)/2;
   dXmean  = 0.5*( dMaxValue + dMinValue );
   dXrange = dMaxValue - dMinValue;

   /* [b] : Compute locations and weighting coefficients */

   for( ii = 1; ii <= im ; ii = ii + 1 ) {

      /* [b.1] : Compute "rough" approximation to root */

      dZ = cos( 3.141592654 * (ii - 0.25)/(iNoPoints + 0.5));

      /* [b.2] : Use Newton's method to refine approximation to root */

      do {
         dP1 = 1.0; dP2 = 0.0;
         for ( ij = 1; ij <= iNoPoints; ij = ij + 1 ) {
             dP3 = dP2; dP2 = dP1;
             dP1 = ((2.0*ij-1.0)*dZ*dP2-(ij-1.0)*dP3)/ij;
         }
         dPp = iNoPoints*(dZ*dP1-dP2)/(dZ*dZ-1.0);
         dZ1 = dZ;
         dZ  = dZ1-dP1/dPp;
      } while ( fabs(dZ-dZ1) > eps );

      /* [b.3] : Scale root to desired interval and put in symmetric counterpart */

      daX [ii-1]         = dXmean - 0.5*dXrange*dZ;
      daX [iNoPoints-ii] = dXmean + 0.5*dXrange*dZ;

      /* [b.4] : Compute weight and its symmetric counterpart */

      daW [ii-1]         = (dXrange/2.0) * 2.0/((1.0-dZ*dZ)*dPp*dPp);
      daW [iNoPoints-ii] = daW [ii-1];
   }
}