Newton Excel Bach, not (just) an Excel Blog

• Lets look at some details of the coding.

The original code (available here) includes routines “intdeini” to generate an array of points and weights of the quadrature formula, and “intde” to carry out the integration of any supplied function “f”.

Modifications to intdeini (renamed getaw) were straightforward:

• an “implicit none” statement was added
• !f2py lines were added for the input and output
• The eps and tinyln constants were hard coded

      subroutine getaw(lenaw, aw)
        implicit none
        integer lenaw
!f2py intent(in) lenaw
        real*8 aw(0:lenaw)
!f2py intent(out)  aw
        real*8 eps
        real*8 efs, hoff
        integer noff, nk, k, j
        real*8 pi2, tinyln, epsln, h0, ehp, ehm, h, t, ep, em, xw, wg
    ! ---- adjustable parameter ----
          efs = 0.1d0
          hoff = 8.45d0
    ! ------------------------------
          eps =   2.22044604925031E-16
          pi2 = 2 * atan(1.0d0)
          tinyln =  708.396418532264
...

The intde routine (renamed intde3) needed a little more work. The code below shows the f2py comment lines required to allow a python function to be called from the Fortran code:

      subroutine intde3(a, b, aw, x, lenaw)
      implicit none
!f2py intent(callback) f
      external f
      real*8 f
!     The following lines define the signature of func for F2PY:
!f2py real*8 y
!f2py y = fun(y)

      real*8 x(3)
!f2py intent(out) x 

      integer lenaw
!f2py intent(in) lenaw
      real*8 aw(0:lenaw)
!f2py intent(in)  aw
      real*8 a, b
!f2py intent(in) a, b

...

The other major change was:

•  Where the calls to the target function included a calculation, this was replaced with a variable.

So lines in the original code such as:

...
        xa = ba * aw(j)
          fa = f(a + xa)
          fb = f(b - xa)
...

were replaced with:

...
          xa = ba * aw(j)
          c = a + xa
          fa = f(c)
          c = b - xa
          fb = f(c)
...

Finally 11 short functions were added to the Fortran code for testing purposes:

•  ...       real*8 function f1(x)
  real*8 x
  !f2py intent(in)  x
  f1 = 1 / sqrt(x)
  end function f1
  ... 

  The Fortran code in tsc2.f90 was then compiled in a form that could be accessed from Python using F2PY:
  Open a Command Prompt window (right-click the Windows Icon and select Command Prompt)

• If necessary, go to the directory containing your Fortran code
• Enter f2py -c -m tsc2 tsc2.f90

F2PY generates a large volume of text during the compilation process, but will terminate either with an error message, or if compilation was successful, with the message shown below:

The resulting compiled file (tsc2.pyd) can be called from Python, in the same way as any other Python module:

To call from Excel(or from another Python function) it is convenient to create an interface function, to first call the getaw function, then pass the function to be integrated in alternative ways:

def Quad_TS(Func, limits, ftype, lenaw, var):
    res = np.zeros(6)
    stime = time.clock()
    x = np.zeros(3)
    a = limits[0]
    b = limits[1]
  
    try:
        aw[0] != 0.
    except:
        aw = tsc2.getaw(lenaw)
        
    stime2 = time.clock()
    if ftype == 0:
        c = tsc2.intde3(a, b, aw, tsc2.f1)
    elif ftype > 0:
        Func = getattr(tsc2, Func)
        c = tsc2.intde3(a, b, aw, Func)
    else:
        if ftype == -2: Func = 'lambda ' + var + ': ' + Func
        c = tsc2.intde3(a, b, aw, eval(Func))
    etime = time.clock()
    res[0] = c[0]
    res[1] = c[1]
    res[2] = c[2]
    res[3] = etime - stime2
    res[4] = etime - stime    
    return res

The ftype argument currently has four options for the function to be integrated:

• ftype = 0 uses the hard coded function tsc2.f1
• ftype > 0 converts the string Func argument into the named function from the tsc2 module.  In future versions this will be extended to allow the module name to be passed, as well as the function name.
• ftype = -2 converts the passed text string into a Python lambda function, with the variable symbol passed in the var argument.
• ftype = any other negative number converts the string Func to a Python Function object, using Eval.

Finally the Python interface function (Quad-TS) is called from Excel, using Excel-Python:

Function Quad_TS(Func As String, Limits As Variant, Optional FType As Long = 1, Optional LenAW As Long = 11150, Optional VarSym As String = "x")
Dim Methods As Variant, Result As Variant, STime As Double, Rtn As Variant, IFunc As String, RunAW As Boolean
    On Error GoTo rtnerr
    STime = MicroTimer
    
    GetArray Limits, 1
        
    Set Methods = Py.Module(modname1)
    Set Result = Py.Call(Methods, "Quad_TS", Py.Tuple(Func, Limits, FType, LenAW, VarSym))
    
    Rtn = Py.Var(Result)
    Rtn(4) = MicroTimer - STime
    Quad_TS = Rtn
    Exit Function

rtnerr:
    Quad_TS = Err.Description

End Function