/* -------------------------------------------------------------
  this sample program produces a command which performs a 
  Simpson's rule numerical integration on either sin, cos
  or exp for a user specified range.  It demonstrates 
  parsing command line arguments, arrays of pointers to 
  functions, symbolic constants, array initialization, etc.
  It should be compiled with the "-lm" flag to include the
  functions for sin, cos and exp . 
  ----------------------------------------------------------------*/

#include <ctype.h>
#include <math.h>

#define SIN 0
#define COS 1
#define EXP 2

 /* the second dimension in the following declaration is
    7 to account for the null terminator */
static char labels[][7] = {"sine","cosine","exp"};

#define TOOBIG  1.e307

/*  argc is the count of arguments passed through the command line;
    argv is a pointer to an array of strings which contain the 
    arguments.  argv[0] is the command name, argv[1] is the first 
    argument, etc. You can also use an additional argument to main,
    envv, which is a pointer to an array of strings containing 
    information about environmental variables.    */

main(int argc,char **argv)

{ 
 int funcode,i,nseg;
 double start,stop;
 /* simp must have a prototype, even though it's defined in this file */
 double simp(double (*func)(),double start,double stop,long n);

 /* this declares funcs as an array of pointers to 
    functions returning double.  The function 
    declarations are in math.h             */
 static double (*funcs[])() = {sin,cos,exp};

 /* set up default values */

  funcode = -1;
  nseg = 20;
  start = stop = TOOBIG;

 /* --------  parse the command line -------*/
 /* the following flags may be used on the command line
   one of the following must be chosen:
      -c  for cosine, -s for sine or -e for exponential
   the start and stop values for integration must be given as
      -b begin_value    -f finish_value 
   by default, the integration is done using 10 segments; 
   to change the number of segments (it must be an even number) use
      -n  number_of_segments                              */

  i = 1;
  while(i < argc)
    {if(argv[i][0] == '-')goto doopt;
     i++;
     continue;

  /* the function atoi takes a pointer to a string, and 
     returns the integer equivalent; atof does the same
     for a double.  type "man 3 atof" for more information */
  
  doopt:
    if(argv[i][1] == 's')funcode = SIN;
    else if(argv[i][1] == 'c')funcode = COS;
    else if(argv[i][1] == 'e')funcode = EXP;
    else if(argv[i][1] == 'b') /* starting value */
      start = atof(argv[++i]);
    else if(argv[i][1] == 'f') /* finish value */
      stop = atof(argv[++i]);
    else if(argv[i][1] == 'n') /* number of iterations */
      nseg = atoi(argv[++i]);

  /* argv[0] points to the name of the command which is being
     executed.  It's often a good idea to print out the name
     of the command in an error message, in case it's not 
     clear where the program was called from */
 
    else 
      printf("(%s) %s unknown option.\n",argv[0],argv[i]);
    i++;  /*  advance to the next command line argument */
   }

  /* check for validity of parameters */
   
   if(funcode == -1)
     {printf("(%s) %s",argv[0],
       "One of -e (exp) -c (cosine) or -s (sine) must be specified.\n");
      exit(0); }

   if(start == TOOBIG || stop == TOOBIG)
     {printf("(%s) %s",argv[0],
       "You must specify both a begin (-b) and a finish (-f) value.\n");
      exit(0); }
    
   if(nseg % 2)   /* nseg must be even */
     {printf("(%s) %s",argv[0],
       "Number of segments must be even.\n");
      exit(0); }

    /* print the answer */

   printf("Integral of %s from %g to %g is %g.\n",
          labels[funcode],start,stop,
          simp(funcs[funcode],start,stop,nseg));
}
 

double simp(double (*func)(),double start,double stop,long n)

{
  double mult,x,t,t1,inc;
  long i;

  /*------------------------------------------------------------- 
    simpson's rule for evaluating an integral is:

   integral from start to stop of f(x) =
     (inc/3) * (f(start) + 4*f(start+inc) + 2*f(start+2*inc)
                      + ... + 4*f(start+(n-1)*inc) + f(stop))

   where inc = (stop - start) / n ,
   and n is the number of segments used to evaluate the integral
   (n must be even)
    -------------------------------------------------------------*/ 

  inc = (stop - start) / (double)n;

  t = func(x = start);
  mult = 4.;
  for(i=1;i<n;i++)
     {x += inc;
      t += mult * func(x);
      mult = mult == 4. ? 2. : 4.; }
  t += func(stop);

  return(t * inc / 3.);
}