Data is stored in a MATLAB readable text file in seven rows for p, q, r, u, v, w, a or b.

To plot a surface in MATLAB, we use surf(x,y,z) where \(x, y, z\) are matrices of second order and dimensions \(m \times n\).

Sphere as mesh of various densities

A “sheet” of points, used by MATLAB to create a mesh of a surface for plotting, consists of columns and rows that form intersecting curves which create a net that MATLAB fills.

If the data is stored in MATLAB-readable format in a text file, such as curves.txt, we can use the load() function to extract all the data into a variable and then distribute it among various variables.

impdata = load("curves.txt", "~ascii");
p      = impdata(1,:);
q      = impdata(2,:);
r      = impdata(3,:);
u      = impdata(4,:);
v      = impdata(5,:);
w      = impdata(6,:);
minor  = impdata(7,:);
clear impdata;
                

To simplify the mathematics of a channel, we will make an analogy to one of its special cases, a right cylindrical surface.

A right cylindrical surface is a channel formed by a series of spheres, with their directrix being a line. Alternatively, a solid cylinder can be thought of as the result of integrating circular discs oriented perpendicular to the directrix. A cylindrical surface, on the other hand, is the result of integrating circles or the perimeters of the discs.

Cylinder as Integration of Circular Sheets

We can visualize a channel as a cylinder with varying cross-sectional radii. This allows us to calculate the position of points on elliptical cross-sections to create a mesh of the channel.

Global polynomial fitting involves approximating a dataset using a single polynomial equation of degree k that best represents the data. Once fitted, the polynomial can be differentiated to compute derivatives (e.g., slopes) at any given data point.

Calculating the slope of a curve can be quite tedious to code. First, we need to determine the coefficients of the fitted polynomial and then use these coefficients to calculate the derivative of the polynomial at each point. MATLAB’s built-in functions polyfit(x,y,k), polyder(func), and polyval(func,val) simplify this process.

% Degree of the polynomial fit
polynomial_order = 3;  % Choose based on the complexity of the data

% Number of points in the dataset
count = size(x, 2);

% Fit polynomials globally for y and z
p_y = polyfit(x, y, polynomial_order);  % Polynomial coefficients for y
p_z = polyfit(x, z, polynomial_order);  % Polynomial coefficients for z

% Compute the derivative coefficients
dp_y = polyder(p_y);  % Derivative coefficients for y
dp_z = polyder(p_z);  % Derivative coefficients for z

% Initialize arrays to store tangents (derivatives) at each point
tangent_y = zeros(1, count);
tangent_z = zeros(1, count);

% Evaluate the derivative polynomials at each point in x
for n = 1:count
tangent_y(n) = polyval(dp_y, x(n));  % Derivative of y at x(n)
tangent_z(n) = polyval(dp_z, x(n));  % Derivative of z at x(n)
end
                

Operational Complexity:

• Polynomial Fitting (polyfit): \( O(n \cdot k + k^3) \) where \( n \) is the number of points and \( k \) is the degree.
• Derivative Computation (polyder): \( O(k) \) for each polynomial.
• Polynomial Evaluation (polyval): \( O(n \cdot k) \).

Total Complexity: \(O(n \cdot k + k^3)\)

Key Observations:

• For large \( n \), the term \( O(n \cdot k) \) dominates.
• For large \( k \), the cubic cost \( O(k^3) \) dominates.
• The code is efficient for small \( k \) but becomes computationally expensive for high-degree polynomials.

Limitations: Overestimation or underestimation may occur if the proper degree is not used.

Spline interpolation is a piecewise polynomial method used to approximate a dataset by constructing smooth polynomials between consecutive points. It avoids overfitting and instabilities seen in global polynomial fitting.

window_size = 5;    % Number of points in the fitting window
max_order   = 3;    % Maximum degree of the polynomial fit
count       = size(x, 2); % Total number of points

tangent_y = zeros(1, count);
tangent_z = zeros(1, count);

for n = 1:count
% Determine the indices for the fitting window
half_window = floor(window_size / 2);
start_idx = max(1, n - half_window);
end_idx = min(count, n + half_window);

x_window = x(start_idx:end_idx);
y_window = y(start_idx:end_idx);
z_window = z(start_idx:end_idx);

% Determine the polynomial order based on the window size
current_order = min(max_order, length(x_window) - 1);

p_y = polyfit(x_window, y_window, current_order);
p_z = polyfit(x_window, z_window, current_order);

dp_y = polyder(p_y);
dp_z = polyder(p_z);

tangent_y(n) = polyval(dp_y, x(n));
tangent_z(n) = polyval(dp_z, x(n));
end
                

Operational Complexity:

• Polynomial Fitting (polyfit): \( O(m \cdot k^2 + k^3) \) where \( m \) is the window size and \( k \) is the degree.
• Derivative Computation (polyder): \( O(k) \).
• Polynomial Evaluation (polyval): \( O(k) \).

Total Complexity per iteration: \(O(m \cdot k^2 + k^3)\)

Limitations: MATLAB’s polyfit may error if consecutive data points are too close (e.g., differences less than \(0.01\)).

Data that we use as test data is very dense and smooth. Therefore, we will use the Central Finite Difference.

Channel Directrix is the curve that passes through the centers of all ellipses.

\(D = \frac{F_{1} + F_{2}}{2}\)

Each ellipse must satisfy:

• Major Axis must pass through the Focal Points.
• Ellipses must be perpendicular to the Directrix; that is, the normal vector of the ellipse is tangential to the directrix.
• Major and Minor Axes must be mutually perpendicular.

In practice, it is challenging to satisfy all conditions. Therefore, we choose to enforce that the Major Axis always passes through the Focal Points.

Let \(\hat{n}\) be a unit tangent vector to the Directrix. Then:

\(\hat{a} = \frac{1}{2}\frac{\begin{bmatrix} F_{x}^{2} - F_{x}^{1} & F_{y}^{2} - F_{y}^{1} & F_{z}^{2} - F_{z}^{1} \end{bmatrix}}{\sqrt{(F_{x}^{2} - F_{x}^{1})^{2} + (F_{y}^{2} - F_{y}^{1})^{2} + (F_{z}^{2} - F_{z}^{1})^{2}}}\)

\(\hat{b} = \hat{n} \times \hat{a}\)
\(\hat{a} = \hat{b} \times \hat{n}\)

Thus, the minor axis is always perpendicular to the tangent of the directrix.

To generate each elliptical cross-section, MATLAB computes the major and minor axes for each ellipse, then shifts the points by the ellipse’s center coordinates.

t      = 0:0.01:2*pi+0.1;
cyl_x  = zeros(size(t'*x));  % Container matrices for coordinates
cyl_y  = cyl_x;
cyl_z  = cyl_x;

for n = 1:size(x,2)
norm_vec  = [1 tangent_y(n) tangent_z(n)];
norm_vec  = norm_vec / norm(norm_vec);  % Unit vector
focal_vec = [(p(n)-u(n))/2 (q(n)-v(n))/2 (r(n)-w(n))/2];
minor_vec = cross(norm_vec, focal_vec);
minor_vec = minor_vec / norm(minor_vec);
major_vec = focal_vec / norm(focal_vec);
point_vec = (major(n) * (cos(t)' * major_vec) + minor(n) * (sin(t)' * minor_vec));

% Each column corresponds to one elliptical cross-section
cyl_x(:,n) = point_vec(:,1);
cyl_y(:,n) = point_vec(:,2);
cyl_z(:,n) = point_vec(:,3);
end

k       = ones(size(t));
shift_x = k' * x;
shift_y = k' * y;
shift_z = k' * z;
cyl_x   = cyl_x + shift_x;
cyl_y   = cyl_y + shift_y;
cyl_z   = cyl_z + shift_z;