BSplines.m

%% Get the B-spline from control points that are specified manually
%% Saikat Mukherjee Purdue University July 5 2022

%% Console
close all;
clear all;
format long;
clc;

R      = 461.5;
theta  = 300;
rho_l  = 996.513027468132;
rho_v  = 0.025589673682918975;
p_sat  = 0.003536806752274e6;
dp_l   = 22.3e6;
dp_v   = 1.38e6;
DeltaP = 1.0e2;
mu_sat = -5.36e3;

%% Control Point 1 : Vapor Saturation
rho1 = rho_v;
p1   = p_sat;

%% Control Point 6 : Liquid saturation
rho6 = rho_l;
p6   = p_sat;

%% Control Point 2 : Match the slope from the vapor side EoS
rho2 = rho1 + DeltaP/dp_v;
p2   = p1 + DeltaP;

%% Control Point 5 : Match the slope from liquid EoS
rho5 = rho6 - DeltaP/dp_l;
p5   = p6 - DeltaP;

%% Control Point 3
rho3 = 0.25*(rho_l+3*rho_v);
p3   = p2;

%% Control Point 4
rho4 = 0.25*(rho_v+3*rho_l);
p4   = p5;

% Control points
C = [rho1 rho2 rho3 rho4 rho5 rho6;  % x components
     p1 p2 p3 p4 p5 p6]; % y components

% B-Spline Degree
m = 2;

% B-Spline
s = f_Bspline(C, m, 0.01);

rho_global = s(1,:);
p_global   = s(2,:);

save test4.dat s -ascii

psi_rhov = rho_v*R*theta*-1.036970015019896;
psi_rhol = rho_l*R*theta*-0.038751812411856;
C_Integration = ((psi_rhol-psi_rhov) - rho_l*trapz( rho_global, s(2,:)./s(1,:).^2 ))/(rho_l-rho_v);
for i =2:400
    rho_subset = rho_global(1:i);
    integrand_subset = s(2,1:1:i)./s(1,1:1:i).^2;
    psi(i-1) = psi_rhov + rho_global(i)*trapz(rho_subset, integrand_subset) + C_Integration*(rho_global(i)-rho_v);
end
psi_global = [psi_rhov psi psi_rhol];

f   = sqrt(abs( 2*( psi_global-mu_sat*rho_global ) ));
Integral = trapz(rho_global,f);
lambda = (0.0728/Integral)^2;

%% Plot

figure(1);
plot(C(1,:), C(2,:),'o');
hold on;
plot(C(1,:), C(2,:),'--');
hold on;
plot(rho_global,p_global);
grid off;
xlabel('x');
ylabel('y');
legend({'Control points','Polygon','B-Spline-Approx.'}, 'Location', 'southeast');
ylim([min(C(2,:)) max(s(2,:))]);

figure(2);
plot(rho_global,psi_global,'r');
hold on;
plot(rho_global,mu_sat*rho_global,'k');

function s = f_Bspline(C, m, step)
% Calculates a B-Spline of degree m using the control points C.
% 
% C: 2-dimensional control points (x_0, ... , x_n; y_0, ... , y_n)
% m: B-Spline degree
% step: Step size of t.
%
% s: B-Spline. s(1,:) -> x component, s(2,:) -> y component

%% Parameters

% Control point's x and y components
x = C(1,:);
y = C(2,:);

% Number of control point - 1
n = size(x,2) - 1;

% Knot vector
T = f_BSpline_KnotVector(m,n);

% B-Spline intervall
t = 0:step:(n-m+1);

%% Calculate B-Spline

for z=1:1:size(t,2)
    ti = t(z);
    s(1,z) = 0;
    s(2,z) = 0;
    for i=0:1:n
        % Base B-Spline
        B = f_BSpline_BaseSpline(i, m, ti, T);
        % x component
        s(1,z) = s(1,z) + x(i+1) * B;
        % y component
        s(2,z) = s(2,z) + y(i+1) * B;
    end
end

end

function T = f_BSpline_KnotVector(m,n)
% Calculate knot sequence.
% m: Degree of B-Spline
% n: Number of control points - 1
%
% T = [t0, ... , t_n+m+1] : Knot sequence / vector

T = zeros(1, (n+m+2));
for j=0:1:(n+m+1)
    if j <= m
        Ti = 0;
    elseif j >= (m+1) && j <= n
        Ti = j - m;
    elseif j > n
        Ti = n - m + 1;
    end
    T(j+1) = Ti;
end

end

function B = f_BSpline_BaseSpline(i, k, t, T)
% Calculate Base B-Spline B_i,k
% k: B-Spline degree
% T: Knot sequence
% t: Current t parameter
%
% B: Base B-Spline at t.

% Index shift
j = i + 1;

if k == 0
    % Corrected end of recursion
    if (t >= T(j) && t < T(j+1))
        B = 1;
    % Exception: the last basis function is equal to 1 also at the last knot.
    % There might be a more elegant way of writing it; I am no Matlab expert.
    elseif ( T(j+1) == T(end) && t == T(end))
        B = 1;
    else
        B = 0;
    end
else
    % Check dividing by zero
    if T(j+k) ~= T(j)
        A = (t -T(j))/(T(j+k) - T(j));
    else
        A = 0;
    end
    if T(j+k+1) ~= T(j+1)
        B = (T(j+k+1) - t) / (T(j+k+1) - T(j+1));
    else
        B = 0;
    end
    % Calculate base B-Spline
    B1 = f_BSpline_BaseSpline(i,   k-1, t, T);
    B2 = f_BSpline_BaseSpline(i+1, k-1, t, T);
    B = A * B1 + B * B2;
end


end