Attribute VB_Name = "INTEGRATION_GAULEG_LIBR"

'-------------------------------------------------------------------------------------
'-------------------------------------------------------------------------------------

Option Explicit     'Requires that all variables to be declared explicitly.
Option Base 1       'The "Option Base" statement allows to specify 0 or 1 as the
                    'default first index of arrays.
'-------------------------------------------------------------------------------------
'-------------------------------------------------------------------------------------


'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : INTEGRATION_7_GAULEG_FUNC
'DESCRIPTION   : Find the area under the curvey = f (x), -1 = x = 1.
'weights for 7-point Gauss-Legendre integration
'(only 4 values out of 7 are given as they are symmetric)
'LIBRARY       : INTEGRATION
'GROUP         : GAULEG
'ID            : 001
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 21/01/2009
'************************************************************************************
'************************************************************************************

Function INTEGRATION_7_GAULEG_FUNC(ByVal FUNC_STR_NAME As String, _
ByVal LOWER_BOUND As Double, _
ByVal UPPER_BOUND As Double, _
ByVal tolerance As Double, _
Optional ByVal nLOOPS As Long = 30)
  
  ' weights for 7-point Gauss-Legendre integration
  ' (only 4 values out of 7 are given as they are symmetric)
  
  Dim i As Single
  Dim j As Single
  
  Dim k As Long

  Dim TEMP_SUM As Double
  Dim FUNC_VAL As Double

  Dim TEMP_HALF As Double
  Dim TEMP_CENTER As Double
  
  Dim ATEMP_VAL As Double 'ATEMP_VAL (abscissa) and f(ATEMP_VAL)
  Dim BTEMP_VAL As Double 'will be result of BTEMP_VAL integral
  Dim CTEMP_VAL As Double 'will be result of CTEMP_VAL integral

  Dim ATEMP_ARR() As Double
  Dim BTEMP_ARR() As Double
  Dim CTEMP_ARR() As Double
  
  On Error GoTo ERROR_LABEL
  
  k = 0
  ReDim ATEMP_ARR(0 To 3)
  
  ATEMP_ARR(0) = 0.417959183673469
  ATEMP_ARR(1) = 0.381830050505119
  ATEMP_ARR(2) = 0.279705391489277
  ATEMP_ARR(3) = 0.12948496616887
  
  ReDim BTEMP_ARR(0 To 7) ' weights for 15-point Gauss-Kronrod integration
  
  BTEMP_ARR(0) = 0.209482141084728
  BTEMP_ARR(1) = 0.204432940075298
  BTEMP_ARR(2) = 0.190350578064785
  BTEMP_ARR(3) = 0.169004726639267
  BTEMP_ARR(4) = 0.140653259715525
  BTEMP_ARR(5) = 0.10479001032225
  BTEMP_ARR(6) = 0.063092092629979
  BTEMP_ARR(7) = 0.022935322010529
  
  ReDim CTEMP_ARR(0 To 7) ' abscissae (evaluation points) for 15-point
  'Gauss-Kronrod integration
  
  CTEMP_ARR(0) = 0#
  CTEMP_ARR(1) = 0.207784955007898
  CTEMP_ARR(2) = 0.405845151377397
  CTEMP_ARR(3) = 0.586087235467691
  CTEMP_ARR(4) = 0.741531185599394
  CTEMP_ARR(5) = 0.864864423359769
  CTEMP_ARR(6) = 0.949107912342758
  CTEMP_ARR(7) = 0.991455371120813
  
  TEMP_HALF = (UPPER_BOUND - LOWER_BOUND) / 2
  TEMP_CENTER = (LOWER_BOUND + UPPER_BOUND) / 2
  FUNC_VAL = Excel.Application.Run(FUNC_STR_NAME, TEMP_CENTER)
  BTEMP_VAL = FUNC_VAL * ATEMP_ARR(0)
  CTEMP_VAL = FUNC_VAL * BTEMP_ARR(0)
  
  ' calculate BTEMP_VAL and half of CTEMP_VAL
  i = 2
  For j = 1 To 3
    ATEMP_VAL = TEMP_HALF * CTEMP_ARR(i)
    TEMP_SUM = Excel.Application.Run(FUNC_STR_NAME, (TEMP_CENTER - ATEMP_VAL)) + _
               Excel.Application.Run(FUNC_STR_NAME, (TEMP_CENTER + ATEMP_VAL))
    BTEMP_VAL = BTEMP_VAL + TEMP_SUM * ATEMP_ARR(j)
    CTEMP_VAL = CTEMP_VAL + TEMP_SUM * BTEMP_ARR(i)
    i = i + 2
  Next j
  
  ' calculate other half of CTEMP_VAL
  For i = 1 To 7 Step 2
  
      ATEMP_VAL = TEMP_HALF * CTEMP_ARR(i)
      TEMP_SUM = Excel.Application.Run(FUNC_STR_NAME, (TEMP_CENTER - ATEMP_VAL)) + _
        Excel.Application.Run(FUNC_STR_NAME, (TEMP_CENTER + ATEMP_VAL))
      CTEMP_VAL = CTEMP_VAL + TEMP_SUM * BTEMP_ARR(i)
  
  Next i
  
  ' multiply by (LOWER_BOUND - UPPER_BOUND) / 2
  BTEMP_VAL = TEMP_HALF * BTEMP_VAL
  CTEMP_VAL = TEMP_HALF * CTEMP_VAL
  
  ' 15 more function evaluations have been used
  k = k + 15
  
  ' error is <= CTEMP_VAL - BTEMP_VAL
  ' if error is larger than tolerance then split the interval
  ' in two and integrate recursively
  If (Abs(CTEMP_VAL - BTEMP_VAL) < tolerance) Then
     INTEGRATION_7_GAULEG_FUNC = CTEMP_VAL
     Exit Function
  Else
      If k + 30 > nLOOPS Then: GoTo ERROR_LABEL
      'maximum number of function evaluations exceeded"
     INTEGRATION_7_GAULEG_FUNC = _
     INTEGRATION_7_GAULEG_FUNC(FUNC_STR_NAME, LOWER_BOUND, TEMP_CENTER, tolerance / 2, nLOOPS) + _
     INTEGRATION_7_GAULEG_FUNC(FUNC_STR_NAME, TEMP_CENTER, UPPER_BOUND, tolerance / 2, nLOOPS)
  End If

Exit Function
ERROR_LABEL:
INTEGRATION_7_GAULEG_FUNC = Err.number
End Function


'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : INTEGRATION_8_GAULEG_FUNC
'DESCRIPTION   : Gauss-Legendre performed over 8 subintervals
'LIBRARY       : INTEGRATION
'GROUP         : GAULEG
'ID            : 002
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 21/01/2009
'************************************************************************************
'************************************************************************************

Function INTEGRATION_8_GAULEG_FUNC(ByVal FORMULA_STRING As String, _
ByVal VARIABLE_STRING As String, _
ByVal MIN_VAL As Double, _
ByVal MAX_VAL As Double, _
Optional ByVal INT_VALUE As Long = 8)

'-------------------------------------------------------------------
'a = 0
'b = 1

'Debug.Print INTEGRATION_8_GAULEG_FUNC("sin(x)", "x", a, b, 8)
'Debug.Print " 0.4596976941318603 (exact value)"

'Debug.Print INTEGRATION_8_GAULEG_FUNC("exp(@x)", "@x", a, b, 8)
'Debug.Print " 1.7182818284590452 (exact value)"
'-------------------------------------------------------------------


' INT_VALUE--> SubInterval: change as you like
' to any positive integer

Dim i As Long   ' COUNTER
Dim j As Long   ' COUNTER
Dim k As Long   ' degrees

Dim ATEMP_VAL As Double      ' subinterval ATEMP ... BTEMP
Dim BTEMP_VAL As Double

Dim TEMP_BETA As Double      '   "
Dim TEMP_ALPHA As Double     ' bounds to give +- 1 by change of variables
Dim TEMP_VALUE As Double        ' integration variable
Dim TEMP_LENGTH As Double   ' of length (MAX_VAL-MIN_VAL)/pieces

Dim TEMP_SUM As Double  ' to sum up
Dim TEMP_TOTAL As Double  ' summing over subintervals

Dim ZERO_VECTOR() As Double
Dim WEIGHT_VECTOR() As Double

On Error GoTo ERROR_LABEL

'--------------------------------------------------------------------------
k = 16
'--------------------------------------------------------------------------
ReDim ZERO_VECTOR(1 To k) ' Gauss zeros
ReDim WEIGHT_VECTOR(1 To k) ' Gauss weights
'--------------------------------------------------------------------------
'----Pre-computed data for Gauss-Legendre integration, do not change this--
'--------------------------------------------------------------------------

ZERO_VECTOR(1) = -0.98940093499165
ZERO_VECTOR(2) = -0.944575023073233
ZERO_VECTOR(3) = -0.865631202387832
ZERO_VECTOR(4) = -0.755404408355003
ZERO_VECTOR(5) = -0.617876244402644
ZERO_VECTOR(6) = -0.458016777657227
ZERO_VECTOR(7) = -0.281603550779259
ZERO_VECTOR(8) = -9.50125098376374E-02
ZERO_VECTOR(9) = 9.50125098376374E-02
ZERO_VECTOR(10) = 0.281603550779259
ZERO_VECTOR(11) = 0.458016777657227
ZERO_VECTOR(12) = 0.617876244402644
ZERO_VECTOR(13) = 0.755404408355003
ZERO_VECTOR(14) = 0.865631202387832
ZERO_VECTOR(15) = 0.944575023073233
ZERO_VECTOR(16) = 0.98940093499165
'-------------------------------------------------------------------------
WEIGHT_VECTOR(1) = 2.71524594117541E-02
WEIGHT_VECTOR(2) = 6.22535239386479E-02
WEIGHT_VECTOR(3) = 9.51585116824928E-02
WEIGHT_VECTOR(4) = 0.124628971255534
WEIGHT_VECTOR(5) = 0.149595988816577
WEIGHT_VECTOR(6) = 0.169156519395003
WEIGHT_VECTOR(7) = 0.182603415044924
WEIGHT_VECTOR(8) = 0.189450610455069
WEIGHT_VECTOR(9) = 0.189450610455069
WEIGHT_VECTOR(10) = 0.182603415044924
WEIGHT_VECTOR(11) = 0.169156519395003
WEIGHT_VECTOR(12) = 0.149595988816577
WEIGHT_VECTOR(13) = 0.124628971255534
WEIGHT_VECTOR(14) = 9.51585116824928E-02
WEIGHT_VECTOR(15) = 6.22535239386479E-02
WEIGHT_VECTOR(16) = 2.71524594117541E-02
'-------------------------------------------------------------------------

TEMP_LENGTH = (MAX_VAL - MIN_VAL) / INT_VALUE

ATEMP_VAL = MIN_VAL - TEMP_LENGTH
BTEMP_VAL = ATEMP_VAL + TEMP_LENGTH

TEMP_TOTAL = 0
For j = 1 To INT_VALUE
  ATEMP_VAL = ATEMP_VAL + TEMP_LENGTH
  BTEMP_VAL = ATEMP_VAL + TEMP_LENGTH

  TEMP_ALPHA = (BTEMP_VAL - ATEMP_VAL) / 2
  TEMP_BETA = (BTEMP_VAL + ATEMP_VAL) / 2
  TEMP_SUM = 0
  For i = 1 To k
    TEMP_VALUE = TEMP_ALPHA * ZERO_VECTOR(i) + _
        TEMP_BETA  ' change of variables for integration bounds
    TEMP_SUM = TEMP_SUM + WEIGHT_VECTOR(i) * _
    CALL_8_GAULEG_OBJ_FUNC(TEMP_VALUE, FORMULA_STRING, VARIABLE_STRING)
  Next i
  TEMP_TOTAL = TEMP_TOTAL + TEMP_SUM * TEMP_ALPHA _
  'change of variables for integration bounds
Next j

INTEGRATION_8_GAULEG_FUNC = TEMP_TOTAL

Exit Function
ERROR_LABEL:
INTEGRATION_8_GAULEG_FUNC = Err.number
End Function


'************************************************************************************
'************************************************************************************
'© Copyright NicoSystem 2009. All rights reserved by Rafael Nicolas Fermin Cota,
'San Francisco, CA. USA  www.rnfc.org
'nfermincota.hba2005@ivey.ca
'************************************************************************************
'************************************************************************************
'FUNCTION      : CALL_8_GAULEG_OBJ_FUNC
'DESCRIPTION   : Parse and evaluate Gauss-Legendre integration formula
'LIBRARY       : INTEGRATION
'GROUP         : GAULEG
'ID            : 003
'AUTHOR        : RAFAEL NICOLAS FERMIN COTA
'LAST UPDATE   : 21/01/2009
'************************************************************************************
'************************************************************************************

Private Function CALL_8_GAULEG_OBJ_FUNC(ByVal TEMP_VALUE As Double, _
ByVal FORMULA_STRING As String, _
Optional ByVal VARIABLE_STRING As String = "@x")

Dim TEMP_STR As String

On Error GoTo ERROR_LABEL

'Debug.Print CALL_8_GAULEG_OBJ_FUNC(1, "sin(x)", "x")
'Debug.Print CALL_8_GAULEG_OBJ_FUNC(0.5, "anyFct(@X)", "@X")

TEMP_STR = Replace(CStr(TEMP_VALUE), ",", ".") ' xlDecimalSeparator
TEMP_STR = Replace(FORMULA_STRING, VARIABLE_STRING, TEMP_STR)

CALL_8_GAULEG_OBJ_FUNC = Evaluate(Evaluate(TEMP_STR))

Exit Function
ERROR_LABEL:
CALL_8_GAULEG_OBJ_FUNC = Err.number
End Function