Commit fc98486c authored by Matthias Steffen's avatar Matthias Steffen
Browse files

missing Python function added

parent 58b876bd
def mcubic(yinput, xinput, xout):
import numpy as np
#
nout=xout.size
nx = xinput.size
ny = yinput.size
if (nx == ny):
n=nx
else:
exit('xinput and yinput must have the same number of elements!')
#
# print('xinput: ',xinput )
# print('yinput: ',yinput )
# print('xout:',xout)
# --- we assume that xinput is a sequence at least 4 elements,
# composed of unique, monotonically increasing or decraesing values ---
if (n < 4):
exit('xinput, yinput must have at least 4 elements!')
if (xinput[0] < xinput[n-1]):
xin=xinput
yin=yinput
else:
xin=np.flip(xinput)
yin=np.flip(yinput)
#
# --- Compute x-step h, linear slope s and slope of parabola p ---
h= xin[1:n] - xin[0:n-1]
s=(yin[1:n] - yin[0:n-1])/h
p=(s[0:n-2]*h[1:n-1]+s[1:n-1]*h[0:n-2])/(h[0:n-2]+h[1:n-1])
q0=np.sign(s[0:n-2])+np.sign(s[1:n-1])
q1=np.array([ abs(s[0:n-2]), abs(s[1:n-1]), 0.5*abs(p) ])
q2=q0*q1.min(axis=0)
y1=np.empty(n,dtype='float64')
y1[1:n-1]=q2
# print('h:',h)
# print('s:',s)
# print('p:',p)
# print('q0:',q0)
# print('q1:',q1)
# print('q2:',q2,q2.dtype)
# print('y1:',y1,y1.dtype)
#
# --- Boundary handling d2y/dx2=0 [bcfac=1.5] ---
# --- y1(0) =1.5*s(0) -0.5*y1(1)
# --- y1(n-1)=1.5*s(n-2)-0.5*y1(n-2)
yleft=1.5*s[0] - 0.5*y1[1]
yrght=1.5*s[n-2] - 0.5*y1[n-2]
y1[0] =yleft
y1[n-1]=yrght
# --- Compute coefficients for cubic polynomial (inner points) ---
a=(y1[0:n-1]+y1[1:n]-2.0*s)/h**2
b=(3.0*s-2.0*y1[0:n-1]-y1[1:n])/h
c=y1[ 0:n-1]
d=yin[0:n-1]
# print('a:',a)
# print('b:',b)
# print('c:',c)
# print('d:',d)
#
# --- Linear extrapolation at upper boundary ---
a = np.append(a,0.0)
b = np.append(b,0.0)
c = np.append(c,y1[ n-1])
d = np.append(d,yin[n-1])
# print('a:',a)
# print('b:',b)
# print('c:',c)
# print('d:',d)
#
# --- Linear extrapolation at lower boundary ---
a = np.append(0.0,a)
b = np.append(0.0,b)
c = np.append(y1[ 0],c)
d = np.append(yin[0],d)
# print('a:',a)
# print('b:',b)
# print('c:',c)
# print('d:',d)
#
# --- Determine interval indices and relative x value ---
iout=np.empty(nout,dtype='int64')
xrel=np.empty(nout,dtype='float64')
for j in range(nout):
if (xout[j] < xin[0]):
iout[j]=0
xrel[j]=xout[j]-xin[0]
elif (xout[j] > xin[n-1]):
iout[j]=n
xrel[j]=xout[j]-xin[n-1]
else:
for i in range(1,n):
if (xin[i] >= xout[j]):
iout[j]=i
xrel[j]=xout[j]-xin[i-1]
break
#
# print('iout:',iout)
# print('xrel:',xrel)
# --- Interpolation ---
yout=((a[iout]*xrel+b[iout])*xrel+c[iout])*xrel+d[iout]
#
return yout
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment