/****************************************************************/
/*                                                              */
/*          Integral package  -  version 27.01.89               */
/*                                                              */
/****************************************************************/

#include "comal/comal.h"

extern CmlFlt NumFunc();

int c2 = 2, c3 = 3, c6 = 6;

#define f(x)  (NumFunc(fnct,x))

void integral(wrkspc,wrktop)      /* integral(fnct(),a,b)       */
CmlFlt  *wrkspc,*wrktop;          /* pointers to workspace      */
{
  struct refparam *fnct;
  CmlFlt a,b,xi,dx,sum;           /* Floating point variables   */
  int i,n;
  ULONG *arg;                     /* Pointer to arguments       */

  b = *(--wrkspc); a = *(--wrkspc);
  fnct = (struct refparam *)wrkspc;
  (--fnct);
  if (fnct->Type!=fltuser && fnct->Type!=fltcode) RetStat(76);
  if (fnct->Type==fltuser) 
    if (fnct->Attributes & (IMPORTBIT | CLOSEDBIT | EXTERNALBIT)==0)
      RetStat(76);
  arg = (ULONG *)fnct->Data;
  if (*arg != 0x01002F) RetStat(76);

  WorkBot = (APTR)wrkspc;         /* Set work space bottom ..   */
  WorkTop = (APTR)wrktop;         /* .. and top                 */

  n = 32;
  dx = (b-a)/n;
  xi = a; sum = (f(a)-f(b))*dx/c6;
  for (i = 1; i<=n; i++)
    {
     xi = xi+dx;
     sum = sum+(f(xi)+c2*f(xi-dx/c2))*dx/c3;
    };
  RetValue(&sum);
}

                              /**********************************/
struct varnode intlnk = {     /* Identifier table for integral  */
  0L,                         /* Last link                      */
  "INTEGRAL",                 /* Name                           */
  fltcode,                    /* Type                           */
  3                      };   /* one argument                   */
UWORD intpar1 = fltref;       /* Float reference parameter      */
UWORD intpar2 = fltval;       /* Float value parameter          */
UWORD intpar3 = fltval;       /* Float value parameter          */
void (*introut)() = &integral;/* Routine address                */
                              /**********************************/

topptr topnode = &intlnk;     /* Top link points to cotlnk      */

void init(pck)
struct package *pck;
{
  pck->var = &topnode;        /* Initialise top-id link         */
}