mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
120 lines
3.2 KiB
Python
120 lines
3.2 KiB
Python
import numpy as np
|
|
|
|
|
|
def linearize(x, f, tolerance=0.001):
|
|
"""Return a tabulated representation of a one-variable function
|
|
|
|
Parameters
|
|
----------
|
|
x : Iterable of float
|
|
Initial x values at which the function should be evaluated
|
|
f : Callable
|
|
Function of a single variable
|
|
tolerance : float
|
|
Tolerance on the interpolation error
|
|
|
|
Returns
|
|
-------
|
|
numpy.ndarray
|
|
Tabulated values of the independent variable
|
|
numpy.ndarray
|
|
Tabulated values of the dependent variable
|
|
|
|
"""
|
|
# Make sure x is a numpy array
|
|
x = np.asarray(x)
|
|
|
|
# Initialize output arrays
|
|
x_out = []
|
|
y_out = []
|
|
|
|
# Initialize stack
|
|
x_stack = [x[0]]
|
|
y_stack = [f(x[0])]
|
|
|
|
for i in range(x.shape[0] - 1):
|
|
x_stack.insert(0, x[i + 1])
|
|
y_stack.insert(0, f(x[i + 1]))
|
|
|
|
while True:
|
|
# Get the bounding points currently on the stack
|
|
x_high, x_low = x_stack[-2:]
|
|
y_high, y_low = y_stack[-2:]
|
|
|
|
# Evaluate the function at the midpoint
|
|
x_mid = 0.5*(x_low + x_high)
|
|
y_mid = f(x_mid)
|
|
|
|
# Linearly interpolate between the bounding points
|
|
y_interp = y_low + (y_high - y_low)/(x_high - x_low)*(x_mid - x_low)
|
|
|
|
# Check the error on the interpolated point and compare to tolerance
|
|
error = abs((y_interp - y_mid)/y_mid)
|
|
if error > tolerance:
|
|
x_stack.insert(-1, x_mid)
|
|
y_stack.insert(-1, y_mid)
|
|
else:
|
|
x_out.append(x_stack.pop())
|
|
y_out.append(y_stack.pop())
|
|
if len(x_stack) == 1:
|
|
break
|
|
|
|
x_out.append(x_stack.pop())
|
|
y_out.append(y_stack.pop())
|
|
|
|
return np.array(x_out), np.array(y_out)
|
|
|
|
def thin(x, y, tolerance=0.001):
|
|
"""Check for (x,y) points that can be removed.
|
|
|
|
Parameters
|
|
----------
|
|
x : numpy.ndarray
|
|
Independent variable
|
|
y : numpy.ndarray
|
|
Dependent variable
|
|
tolerance : float
|
|
Tolerance on interpolation error
|
|
|
|
Returns
|
|
-------
|
|
numpy.ndarray
|
|
Tabulated values of the independent variable
|
|
numpy.ndarray
|
|
Tabulated values of the dependent variable
|
|
|
|
"""
|
|
# Initialize output arrays
|
|
x_out = x.copy()
|
|
y_out = y.copy()
|
|
|
|
N = x.shape[0]
|
|
i_left = 0
|
|
i_right = 2
|
|
|
|
while i_left < N - 2 and i_right < N:
|
|
m = (y[i_right] - y[i_left])/(x[i_right] - x[i_left])
|
|
|
|
for i in range(i_left + 1, i_right):
|
|
# Determine error in interpolated point
|
|
y_interp = y[i_left] + m*(x[i] - x[i_left])
|
|
if abs(y[i]) > 0.:
|
|
error = abs((y_interp - y[i])/y[i])
|
|
else:
|
|
error = 2*tolerance
|
|
|
|
if error > tolerance:
|
|
for i_remove in range(i_left + 1, i_right - 1):
|
|
x_out[i_remove] = np.nan
|
|
y_out[i_remove] = np.nan
|
|
i_left = i_right - 1
|
|
i_right = i_left + 1
|
|
break
|
|
|
|
i_right += 1
|
|
|
|
for i_remove in range(i_left + 1, i_right - 1):
|
|
x_out[i_remove] = np.nan
|
|
y_out[i_remove] = np.nan
|
|
|
|
return x_out[np.isfinite(x_out)], y_out[np.isfinite(y_out)]
|