' ******************************************************************
' Integration by trapezoidal rule
' ******************************************************************

FUNCTION TrapInt(X() AS DOUBLE, Y() AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Computes the integral of a function defined by point coordinates
' (X, Y) using the trapezoidal rule
' ------------------------------------------------------------------

  DIM AS DOUBLE  Sum
  DIM AS INTEGER I, J, Lb, Ub

  Lb = LBOUND(X)
  Ub = UBOUND(X)

  Sum = 0
  FOR I = Lb TO Ub - 1
    J = I + 1
    Sum = Sum + 0.5 * (X(J) - X(I)) * (Y(J) + Y(I))
  NEXT I
  TrapInt = Sum
END FUNCTION

SUB ConvTrap(Func1 AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
             Func2 AS FUNCTION(X AS DOUBLE) AS DOUBLE,_
             T()   AS DOUBLE,_
             Y()   AS DOUBLE)
' ------------------------------------------------------------------
' Computes the convolution product of 2 functions Func1 and Func2
' by the trapezoidal rule over an array T(0..N) of equally spaced
' abscissas, with T(0) = 0. The result is returned in Y(0..N)
' ------------------------------------------------------------------

  DIM AS INTEGER N = UBOUND(T)
  DIM AS DOUBLE  F2(0 TO N)
  DIM AS DOUBLE  F10, Sum
  DIM AS INTEGER I, J

  F10 = Func1(0)

  FOR I = 0 TO N
    F2(I) = Func2(T(I))
  NEXT I

  FOR I = 1 TO N
    Sum = 0.5 * (Func1(T(I)) * F2(0) + F10 * F2(I))
    FOR J = 1 TO I - 1
      Sum = Sum + Func1(T(I) - T(J)) * F2(J)
    NEXT J
    Y(I) = Sum * T(1)
  NEXT I
END SUB

