PBPK_public/pythonScripts/filionQuadrature.py

import numpy

def Taylor(x,coef):
    q=1
    s=0
    for c in coef:
        s+=c*q
        q*=x
    return s

def filionAlpha(th,th3,sth,cth):
    if abs(th)<0.1:
        return Taylor(th,[0,0,0,2/45,0,-2/315,0,2/4725])
    
    th2=th*th
    return (th2+sth*cth*th-2*sth*sth)/th3

def filionBeta(th,th3,sth,cth):
    if abs(th)<0.1:
        return Taylor(th,[2/3,0,2/15,0,-4/105,0,2/567,0,-4/22275])
    
    return 2*((1+cth*cth)*th-2*sth*cth)/th3

def filionGamma(th,th3,sth,cth):
    if abs(th)<0.1:
        return Taylor(th,[4/3,0,-2/15,0,1/210,0,-1/11340])
    
    return 4*(sth-th*cth)/th3

def cEven(qs,n):

   idx=numpy.arange(0,2*n+1,2)
   return numpy.sum(qs[idx])

def qFunc(qf,t,x0,h,n,func):

   qx=numpy.linspace(x0,x0+2*h*n,2*n+1)
   qs=[q*func(t*x) for x,q in zip(qx,qf)]
   qs[0]*=0.5
   qs[len(qs)-1]*=0.5
   return numpy.array(qs)
        
def cOdd(qs,n):
   idx=numpy.arange(1,2*n+1,2)
   return numpy.sum(qs[idx])

def filionQuadrature(qf,t,x0,h,n,func,debug=0):
   #calculate integral f(x) func(tx) in [x0,x0+2*n*h] 
   #where func is either cos or sine 
    cFunc=numpy.cos
    iFunc=numpy.sin
    sign=1
    if func=='sin':
        cFunc=numpy.sin
        iFunc=numpy.cos
        sign=-1
    th=t*h
    sth=numpy.sin(th)
    cth=numpy.cos(th)
    th3=th*th*th
    x2n=x0+2*n*h
    A=filionAlpha(th,th3,sth,cth)
    B=filionBeta(th,th3,sth,cth)
    C=filionGamma(th,th3,sth,cth)
    qs=qFunc(qf,t,x0,h,n,cFunc)
    CE=cEven(qs,n)
    CO=cOdd(qs,n)
    D=sign*(qf[2*n]*iFunc(t*x2n)-qf[0]*iFunc(t*x0))
    r=h*(A*D+B*CE+C*CO)
    if debug:
      print('fillon r={:.2f} A={:.2f} D={:.2f} B={:.2f} CE={:.2f} C={:.2f} CO={:.2f}'.format(r,A,D,B,CE,C,CO))
    return r
    
def tabulateFunction(f,x0,n,h):
   qx=[x0+i*h for i in range(2*n+1)]
   qy=[f(x) for x in qx]
   return qx,qy