\ Cauchy principal value integration using
\ 6-point Gauss-Legendre quadrature
\


\ ---------------------------------------------------
\     (c) Copyright 1999  Julian V. Noble.          \
\       Permission is granted by the author to      \
\       use this software for any application pro-  \
\       vided this copyright notice is preserved.   \
\ ---------------------------------------------------

\ This is an ANS Forth program requiring the
\   FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependencies:
\   Assumes independent floating point stack
\   Uses a FORmula TRANslator for clarity
\
\ Description of algorithm:
\   integrate f(x)/(x - x0) from x0-delta to x0+delta
\   using even-order Gauss-Legendre quadrature formula
\
\ Usage: use( fname delta x0 )pval  ( f: -- integral)

\ Example:
\   15 set-precision  ok
\   use( fexp 1e0 0e0 )pval fs. 2.11450175075134E0  ok
\
\   FVARIABLE t
\   : f1   ( f: x -- f1)   f" t="  f" (t-1)/(1-t^3)"  ;
\   use( f1  0.25e0 1e0 )pval fs. 1.67823854224122E-1  ok

MARKER  -pval

: undefined    BL WORD  FIND  NIP  0=  ;

undefined  f"  [IF]  include  ftran201.f  [THEN]

FVARIABLE wt1    FVARIABLE xi1      FVARIABLE y1    FVARIABLE y1m
FVARIABLE wt2    FVARIABLE xi2      FVARIABLE y2    FVARIABLE y2m
FVARIABLE wt3    FVARIABLE xi3      FVARIABLE y3    FVARIABLE y3m

FVARIABLE x0     FVARIABLE delta

f$" xi1 = 0.238619186083197"
f$" xi2 = 0.661209386466265"
f$" xi3 = 0.932469514203152"

f$" wt1 = 0.467913934572691"
f$" wt2 = 0.360761573048139"
f$" wt3 = 0.171324492379170"

v: fdummy

: )pval    ( xt -- )   ( f: delta x0 -- integral)
        defines fdummy          \ vector function name
        x0 F!   delta F!
        f" y1 = x0 + delta*xi1" \ scale interval
        f" y2 = x0 + delta*xi2"
        f" y3 = x0 + delta*xi3"
        f" y1m = x0 - delta*xi1"
        f" y2m = x0 - delta*xi2"
        f" y3m = x0 - delta*xi3"

        f" wt1 * ( fdummy(y1) - fdummy(y1m) ) / xi1 +
           wt2 * ( fdummy(y2) - fdummy(y2m) ) / xi2 +
           wt3 * ( fdummy(y3) - fdummy(y3m) ) / xi3 "
;