Unable to get a continuous curvature for my nozzle contour

Greetings everyone, I am trying to fit a Bezier Curve as shown in the paper provide(fig1), but I am unable to get the desired plot(fig 2 and fig3). I have attached the codes as well as the paper and plot. Any help on this is highly appreciated!.
PS: The inflection point is the last point on xarc, i.e, xarc(end) and Nozzle_WithCircARc.m is the main file. It would be better to run the file cell by cell.
%function Nozzle_WithCircARc(G,Me,n,display)
%% Initialize datapoint matrices
clearvars;close all;clc;
G=1.4;
Me = 2.0;
n = 53; % speed index from pucketts paper
display = 0;
Km = zeros(n,n); % K- vlaues (Constant along right running characteristic lines)
Kp = zeros(n,n); % K+ vlaues (Constant along left running characteristic lines)
Theta = zeros(n,n); % Flow angles relative to the horizontal
Mu = zeros(n,n); % Mach angles
M = zeros(n,n); % Mach Numbers
x = zeros(n,n); % x-coordinates
y = zeros(n,n); % y-coordinates
%% Generate the convergent portion of a nozzle
% The inlet height/area and exit height/area(divergent) are same
% Therefore we have same A/A* = 1.6875
% Corresponding to Mach = 0.3722
%% Find NuMax (maximum expansion angle)
[~, NuMax, ~] = PMF(G,Me,0,0);
ThetaMax = NuMax/2;
%%
%ThetaMax0 = (1.0/1.687)^(2/9) * (NuMax/2);
%% Define some flow parameters of originating characteristic lines
dT = ThetaMax/n;
%no_ref = (ThetaMax – ThetaMax_ref)/dT;
ThetaArc(:,1) = (0:dT:ThetaMax);
NuArc = ThetaArc;
KmArc = ThetaArc + NuArc;
[~, ~, MuArc(:,1)] = PMF(G,0,NuArc(:,1),0);
%% Coordinates of wall along curve from throat
y0 = 1; % Define throat half-height
ThroatCurveRadius = 1.5*y0; % Radius of curvature just downstream of the throat
% for larger factors, ywall deviates from A/A* preferred value is 1.1
%L_e = 1.1 * y0 * sind(ThetaMax);
[xarc, yarc] = Arc(ThroatCurveRadius,ThetaArc); % Finds x- and y-coordinates for given theta-values
yarc(:,1) = yarc(:,1) + y0; % Defines offset due to arc being above horizontal
%% Fill in missing datapoint info along first C+ line
% First centerline datapoint done manually
Km(:,1) = KmArc(2:length(KmArc),1);
Theta(:,1) = ThetaArc(2:length(KmArc),1);
Nu(:,1) = Theta(:,1);
Kp(:,1) = Theta(:,1)-Nu(:,1);
M(1,1) = 1.0001;
Nu(1,1) = 0;
Mu(1,1) = 90;
y(1,1) = 0;
x(1,1) = xarc(2,1) + (y(1,1) – yarc(2,1))/tand((ThetaArc(2,1) – MuArc(2,1) – MuArc(2,1))/2);
% Finds the information at interior nodes along first C+ line
for i=2:n
[M(i,1), Nu(i,1), Mu(i,1)] = PMF(G,0,Nu(i,1),0);
s1 = tand((ThetaArc(i+1,1) – MuArc(i+1,1) + Theta(i,1) – Mu(i,1))/2);
s2 = tand((Theta(i-1,1) + Mu(i-1,1) + Theta(i,1) + Mu(i,1))/2);
x(i,1) = ((y(i-1,1)-x(i-1,1)*s2)-(yarc(i+1,1)-xarc(i+1,1)*s1))/(s1-s2);
y(i,1) = y(i-1,1) + (x(i,1)-x(i-1,1))*s2;

end
%% Find flow properties at remaining interior nodes
for j=2:n;
for i=1:n+1-j;

Km(i,j) = Km(i+1,j-1);

if i==1;

Theta(i,j) = 0;
Kp(i,j) = -Km(i,j);
Nu(i,j) = Km(i,j);
[M(i,j), Nu(i,j), Mu(i,j)] = PMF(G,0,Nu(i,j),0);
s1 = tand((Theta(i+1,j-1)-Mu(i+1,j-1)+Theta(i,j)-Mu(i,j))/2);
x(i,j) = x(i+1,j-1) – y(i+1,j-1)/s1;
y(i,j) = 0;

else

Kp(i,j) = Kp(i-1,j);
Theta(i,j) = (Km(i,j)+Kp(i,j))/2;
Nu(i,j) = (Km(i,j)-Kp(i,j))/2;
[M(i,j), Nu(i,j), Mu(i,j)] = PMF(G,0,Nu(i,j),0);
s1 = tand((Theta(i+1,j-1)-Mu(i+1,j-1)+Theta(i,j)-Mu(i,j))/2);
s2 = tand((Theta(i-1,j)+Mu(i-1,j)+Theta(i,j)+Mu(i,j))/2);
x(i,j) = ((y(i-1,j)-x(i-1,j)*s2)-(y(i+1,j-1)-x(i+1,j-1)*s1))/(s1-s2);
y(i,j) = y(i-1,j) + (x(i,j)-x(i-1,j))*s2;

end

end
end
%% Find wall node information
xwall = zeros(2*n,1);
ywall = xwall;
ThetaWall = ywall;
xwall(1:n,1) = xarc(2:length(xarc),1);
ywall(1:n,1) = yarc(2:length(xarc),1);
ThetaWall(1:n,1) = ThetaArc(2:length(xarc),1);
for i=1:n-1
ThetaWall(n+i,1) = ThetaWall(n-i,1); % criteria for stopping the reflection from the wall
end
%% Location of wall points
for i=1:n

s1 = tand((ThetaWall(n+i-1,1) + ThetaWall(n+i,1))/2);
s2 = tand(Theta(n+1-i,i)+Mu(n+1-i,i));
xwall(n+i,1) = ((y(n+1-i,i)-x(n+1-i,i)*s2)-(ywall(n+i-1,1)-xwall(n+i-1,1)*s1))/(s1-s2);
ywall(n+i,1) = ywall(n+i-1,1) + (xwall(n+i,1)-xwall(n+i-1,1))*s1;

end
%% Provide wall geometry to user
assignin(‘caller’,’xwall’,xwall)
assignin(‘caller’,’ywall’,ywall)
assignin(‘caller’,’Coords’,[xwall ywall])
%% Generate the convergent portion of nozzle
H_in = ywall(end);
L_e = (xwall(end)*(1.0/3.0));
[xconv,yconv] = Convergent_new_3rd(y0,H_in,L_e,n);
%%
% Draw contour and characteristic web
if display == 1

plot(xwall,ywall,’-‘)
axis equal
axis([0 ceil(xwall(length(xwall),1)) 0 ceil(ywall(length(ywall),1))])
hold on
plot(xarc,yarc,’k-‘)

for i=1:n-1
plot(x(1:n+1-i,i),y(1:n+1-i,i))
end
for i=1:n
plot([xarc(i,1) x(i,1)],[yarc(i,1) y(i,1)])
plot([x(n+1-i,i) xwall(i+n,1)],[y(n+1-i,i) ywall(i+n,1)])
end
for c=1:n
for r=2:n+1-c
plot([x(c,r) x(c+1,r-1)],[y(c,r) y(c+1,r-1)])
end
end
%hold on
%contourf(x, y, M, 20, ‘LineColor’, ‘none’); % Draw Mach number contours
%colorbar;
xlabel(‘Length [x/y0]’)
ylabel(‘Height [y/y0]’)

end
%% Non-scaled/non-dimensionalized plot
figure (1)
plot(xwall,ywall,’.b’)
title("Non dimensionalized plot")
axis equal
hold on
plot(xarc,yarc,’-‘)
hold on
plot(xconv,yconv,’-‘)
xlabel(‘xwall’)
ylabel(‘ywall’)
axis equal
%%
%{
plot(xwall,ThetaWall)
xlabel(‘xwall’)
ylabel(‘ThetaWall’)
%%
x1 = linspace(min(xwall),max(xwall),1000);
%y1 = spline(xw,yw,x1);
y1 = interp1(xwall,ywall,x1,’spline’);
figure;
title(‘Subplot 2: Cubic spline interpolation’)
plot(x1,y1,’r’)
%}
%%
%{
for i = 1:length(M)
M_centerline(i) = M(1,i);
end
for i =1:length(x)
x_axial(i) = x(1,i);
end
plot(x_axial’, Mcenterline’,’r’,’LineWidth’,2)
%}
%% Find the scaling factor
Mexit = Me;
[~,~,~,~,area] = flowisentropic(G,Mexit);
Hexit = 11.2798; %mm
Full_Throat_height = Hexit/area;
half_yt = Full_Throat_height/2.0;
%% Scaling of factors
xarc = half_yt.*xarc;
yarc = half_yt.*yarc;
xwall = half_yt.*xwall;
ywall = half_yt.*ywall;
xconv = half_yt.*xconv;
yconv = half_yt.*yconv;
%% Scaled up coordinates in mm
%{
title("Dimensionalized Plot")
plot(xwall,ywall,’-‘)
hold on
plot(xarc,yarc,’-‘)
hold on
plot(xconv,yconv,’-‘)
xlabel(‘xwall’)
ylabel(‘ywall’)
axis equal
%}
%% Combine Coordinates
% Combine xconv and xwall
coords_new_x = [xconv; xwall];
coords_new_y = [yconv; ywall];
%%
figure (2)
hold on;
plot(xconv, yconv, ‘.b’, ‘LineWidth’, 1.5) % Convergent section
plot(xwall, ywall, ‘.r’, ‘LineWidth’, 1.5) % Wall section
plot(xarc,yarc,’.g’,’LineWidth’,1.5) % arc
plot(xconv, -1.*yconv, ‘.b’, ‘LineWidth’, 1.5) % Convergent section
plot(xwall, -1.*ywall, ‘.r’, ‘LineWidth’, 1.5) % Wall section
plot(xarc,-1.*yarc,’.g’,’LineWidth’,1.5) % arc
xlabel(‘x [mm]’);
ylabel(‘y [mm]’);
title(‘Scaled up CD nozzle’);
legend(‘Convergent Section’, ‘Straightening Section’, ‘Initial Expansion (Arc)’);
grid on;
axis equal
%% Export coordinates
x1 = linspace(min(coords_new_x ),max(coords_new_x),1000);
%y1 = spline(xw,yw,x1);
y1 = interp1(coords_new_x,coords_new_y,x1,’spline’);
x2 = x1′;
y2 = y1′;
%figure;
%plot(x1,y1,’r’)
%writematrix(x2,’Spline.xlsx’,’Sheet’,1,’Range’,’A1′);
%writematrix(y2,’Spline.xlsx’,’Sheet’,1,’Range’,’B1′);
%%
%{
plot(x2,y2)
%%
plot(xwall,ywall,’-‘)
axis equal
axis([0 ceil(xwall(length(xwall),1)) 0 ceil(ywall(length(ywall),1))])
hold on
plot(xarc,yarc,’k-‘)
%}
%% Coordinates of wall
%{
writematrix(coords_new_x,’coordinates.xlsx’,’Sheet’,1,’Range’,’A1′);
writematrix(coords_new_y,’coordinates.xlsx’,’Sheet’,1,’Range’,’B1′);
%}
%%
%{
plot(x2,y2,’-k’)
hold on
plot(x2,-1.*(y2),’-r’)
grid on
title("CD Nozzle using a Circular Arc")
axis equal
%}
%% Analyze the discontinuity
% Find second-order derivative of wall contour
dy2_dx2 = zeros(size(xwall));
% Loop over all points except the last two
for i = 1:length(xwall)-2
% Calculate the second-order derivative using the forward difference formula
dy2_dx2(i) = (ywall(i+2) – 2*ywall(i+1) + ywall(i)) / (xwall(i+1) – xwall(i))^2;
end
dy2_dx2(end-1:end) = NaN; % last two points
%% Bezier Curve
threshold = 0.05; % Define a threshold for discontinuity detection (can be adjusted)
discontinuities = find(abs(diff(dy2_dx2)) > threshold);
%%
if ~isempty(discontinuities)
for idx = 1:length(discontinuities)
point_idx = discontinuities(idx); % Index of discontinuity
% Skip a few points around the discontinuity
region_start = max(1, point_idx – 2); % Upstream region
region_end = min(length(xwall), point_idx + 2); % Downstream region
% Use Bézier curve fitting only in this region to smooth the discontinuity
P0 = [xwall(region_start), ywall(region_start)]; % Start point
P1 = [xwall(point_idx), ywall(point_idx)]; % Control point (inflection/discontinuity)
P2 = [xwall(region_end), ywall(region_end)]; % End point
% Parameter t varies between 0 and 1 to interpolate the Bézier curve
t = linspace(0, 1, region_end – region_start + 1)’;
xwall(region_start:region_end) = (1-t).^2 * P0(1) + 2*(1-t).*t * P1(1) + t.^2 * P2(1);
ywall(region_start:region_end) = (1-t).^2 * P0(2) + 2*(1-t).*t * P1(2) + t.^2 * P2(2);
end
end
%%
plot(xwall,ywall,’.b’)
axis equal
%% Checking for continuous second order derivative of wall contour
cont_dy2_dx2 = zeros(size(xwall));
% Loop over all points except the last two (since forward difference needs i+2)
for i = 1:length(xwall)-2
% Calculate the second-order derivative using the forward difference formula
cont_dy2_dx2(i) = (ywall(i+2) – 2*ywall(i+1) + ywall(i)) / (xwall(i+1) – xwall(i))^2;
end
cont_dy2_dx2(end-1:end) = NaN;
plot(xwall./xwall(end), dy2_dx2)Greetings everyone, I am trying to fit a Bezier Curve as shown in the paper provide(fig1), but I am unable to get the desired plot(fig 2 and fig3). I have attached the codes as well as the paper and plot. Any help on this is highly appreciated!.
PS: The inflection point is the last point on xarc, i.e, xarc(end) and Nozzle_WithCircARc.m is the main file. It would be better to run the file cell by cell.
%function Nozzle_WithCircARc(G,Me,n,display)
%% Initialize datapoint matrices
clearvars;close all;clc;
G=1.4;
Me = 2.0;
n = 53; % speed index from pucketts paper
display = 0;
Km = zeros(n,n); % K- vlaues (Constant along right running characteristic lines)
Kp = zeros(n,n); % K+ vlaues (Constant along left running characteristic lines)
Theta = zeros(n,n); % Flow angles relative to the horizontal
Mu = zeros(n,n); % Mach angles
M = zeros(n,n); % Mach Numbers
x = zeros(n,n); % x-coordinates
y = zeros(n,n); % y-coordinates
%% Generate the convergent portion of a nozzle
% The inlet height/area and exit height/area(divergent) are same
% Therefore we have same A/A* = 1.6875
% Corresponding to Mach = 0.3722
%% Find NuMax (maximum expansion angle)
[~, NuMax, ~] = PMF(G,Me,0,0);
ThetaMax = NuMax/2;
%%
%ThetaMax0 = (1.0/1.687)^(2/9) * (NuMax/2);
%% Define some flow parameters of originating characteristic lines
dT = ThetaMax/n;
%no_ref = (ThetaMax – ThetaMax_ref)/dT;
ThetaArc(:,1) = (0:dT:ThetaMax);
NuArc = ThetaArc;
KmArc = ThetaArc + NuArc;
[~, ~, MuArc(:,1)] = PMF(G,0,NuArc(:,1),0);
%% Coordinates of wall along curve from throat
y0 = 1; % Define throat half-height
ThroatCurveRadius = 1.5*y0; % Radius of curvature just downstream of the throat
% for larger factors, ywall deviates from A/A* preferred value is 1.1
%L_e = 1.1 * y0 * sind(ThetaMax);
[xarc, yarc] = Arc(ThroatCurveRadius,ThetaArc); % Finds x- and y-coordinates for given theta-values
yarc(:,1) = yarc(:,1) + y0; % Defines offset due to arc being above horizontal
%% Fill in missing datapoint info along first C+ line
% First centerline datapoint done manually
Km(:,1) = KmArc(2:length(KmArc),1);
Theta(:,1) = ThetaArc(2:length(KmArc),1);
Nu(:,1) = Theta(:,1);
Kp(:,1) = Theta(:,1)-Nu(:,1);
M(1,1) = 1.0001;
Nu(1,1) = 0;
Mu(1,1) = 90;
y(1,1) = 0;
x(1,1) = xarc(2,1) + (y(1,1) – yarc(2,1))/tand((ThetaArc(2,1) – MuArc(2,1) – MuArc(2,1))/2);
% Finds the information at interior nodes along first C+ line
for i=2:n
[M(i,1), Nu(i,1), Mu(i,1)] = PMF(G,0,Nu(i,1),0);
s1 = tand((ThetaArc(i+1,1) – MuArc(i+1,1) + Theta(i,1) – Mu(i,1))/2);
s2 = tand((Theta(i-1,1) + Mu(i-1,1) + Theta(i,1) + Mu(i,1))/2);
x(i,1) = ((y(i-1,1)-x(i-1,1)*s2)-(yarc(i+1,1)-xarc(i+1,1)*s1))/(s1-s2);
y(i,1) = y(i-1,1) + (x(i,1)-x(i-1,1))*s2;

end
%% Find flow properties at remaining interior nodes
for j=2:n;
for i=1:n+1-j;

Km(i,j) = Km(i+1,j-1);

if i==1;

else

end

plot(xwall,ywall,’-‘)
axis equal
axis([0 ceil(xwall(length(xwall),1)) 0 ceil(ywall(length(ywall),1))])
hold on
plot(xarc,yarc,’k-‘)

end
%% Non-scaled/non-dimensionalized plot
figure (1)
plot(xwall,ywall,’.b’)
title("Non dimensionalized plot")
axis equal
hold on
plot(xarc,yarc,’-‘)
hold on
plot(xconv,yconv,’-‘)
xlabel(‘xwall’)
ylabel(‘ywall’)
axis equal
%%
%{
plot(xwall,ThetaWall)
xlabel(‘xwall’)
ylabel(‘ThetaWall’)
%%
x1 = linspace(min(xwall),max(xwall),1000);
%y1 = spline(xw,yw,x1);
y1 = interp1(xwall,ywall,x1,’spline’);
figure;
title(‘Subplot 2: Cubic spline interpolation’)
plot(x1,y1,’r’)
%}
%%
%{
for i = 1:length(M)
M_centerline(i) = M(1,i);
end
for i =1:length(x)
x_axial(i) = x(1,i);
end
plot(x_axial’, Mcenterline’,’r’,’LineWidth’,2)
%}
%% Find the scaling factor
Mexit = Me;
[~,~,~,~,area] = flowisentropic(G,Mexit);
Hexit = 11.2798; %mm
Full_Throat_height = Hexit/area;
half_yt = Full_Throat_height/2.0;
%% Scaling of factors
xarc = half_yt.*xarc;
yarc = half_yt.*yarc;
xwall = half_yt.*xwall;
ywall = half_yt.*ywall;
xconv = half_yt.*xconv;
yconv = half_yt.*yconv;
%% Scaled up coordinates in mm
%{
title("Dimensionalized Plot")
plot(xwall,ywall,’-‘)
hold on
plot(xarc,yarc,’-‘)
hold on
plot(xconv,yconv,’-‘)
xlabel(‘xwall’)
ylabel(‘ywall’)
axis equal
%}
%% Combine Coordinates
% Combine xconv and xwall
coords_new_x = [xconv; xwall];
coords_new_y = [yconv; ywall];
%%
figure (2)
hold on;
plot(xconv, yconv, ‘.b’, ‘LineWidth’, 1.5) % Convergent section
plot(xwall, ywall, ‘.r’, ‘LineWidth’, 1.5) % Wall section
plot(xarc,yarc,’.g’,’LineWidth’,1.5) % arc
plot(xconv, -1.*yconv, ‘.b’, ‘LineWidth’, 1.5) % Convergent section
plot(xwall, -1.*ywall, ‘.r’, ‘LineWidth’, 1.5) % Wall section
plot(xarc,-1.*yarc,’.g’,’LineWidth’,1.5) % arc
xlabel(‘x [mm]’);
ylabel(‘y [mm]’);
title(‘Scaled up CD nozzle’);
legend(‘Convergent Section’, ‘Straightening Section’, ‘Initial Expansion (Arc)’);
grid on;
axis equal
%% Export coordinates
x1 = linspace(min(coords_new_x ),max(coords_new_x),1000);
%y1 = spline(xw,yw,x1);
y1 = interp1(coords_new_x,coords_new_y,x1,’spline’);
x2 = x1′;
y2 = y1′;
%figure;
%plot(x1,y1,’r’)
%writematrix(x2,’Spline.xlsx’,’Sheet’,1,’Range’,’A1′);
%writematrix(y2,’Spline.xlsx’,’Sheet’,1,’Range’,’B1′);
%%
%{
plot(x2,y2)
%%
plot(xwall,ywall,’-‘)
axis equal
axis([0 ceil(xwall(length(xwall),1)) 0 ceil(ywall(length(ywall),1))])
hold on
plot(xarc,yarc,’k-‘)
%}
%% Coordinates of wall
%{
writematrix(coords_new_x,’coordinates.xlsx’,’Sheet’,1,’Range’,’A1′);
writematrix(coords_new_y,’coordinates.xlsx’,’Sheet’,1,’Range’,’B1′);
%}
%%
%{
plot(x2,y2,’-k’)
hold on
plot(x2,-1.*(y2),’-r’)
grid on
title("CD Nozzle using a Circular Arc")
axis equal
%}
%% Analyze the discontinuity
% Find second-order derivative of wall contour
dy2_dx2 = zeros(size(xwall));
% Loop over all points except the last two
for i = 1:length(xwall)-2
% Calculate the second-order derivative using the forward difference formula
dy2_dx2(i) = (ywall(i+2) – 2*ywall(i+1) + ywall(i)) / (xwall(i+1) – xwall(i))^2;
end
dy2_dx2(end-1:end) = NaN; % last two points
%% Bezier Curve
threshold = 0.05; % Define a threshold for discontinuity detection (can be adjusted)
discontinuities = find(abs(diff(dy2_dx2)) > threshold);
%%
if ~isempty(discontinuities)
for idx = 1:length(discontinuities)
point_idx = discontinuities(idx); % Index of discontinuity
% Skip a few points around the discontinuity
region_start = max(1, point_idx – 2); % Upstream region
region_end = min(length(xwall), point_idx + 2); % Downstream region
% Use Bézier curve fitting only in this region to smooth the discontinuity
P0 = [xwall(region_start), ywall(region_start)]; % Start point
P1 = [xwall(point_idx), ywall(point_idx)]; % Control point (inflection/discontinuity)
P2 = [xwall(region_end), ywall(region_end)]; % End point
% Parameter t varies between 0 and 1 to interpolate the Bézier curve
t = linspace(0, 1, region_end – region_start + 1)’;
xwall(region_start:region_end) = (1-t).^2 * P0(1) + 2*(1-t).*t * P1(1) + t.^2 * P2(1);
ywall(region_start:region_end) = (1-t).^2 * P0(2) + 2*(1-t).*t * P1(2) + t.^2 * P2(2);
end
end
%%
plot(xwall,ywall,’.b’)
axis equal
%% Checking for continuous second order derivative of wall contour
cont_dy2_dx2 = zeros(size(xwall));
% Loop over all points except the last two (since forward difference needs i+2)
for i = 1:length(xwall)-2
% Calculate the second-order derivative using the forward difference formula
cont_dy2_dx2(i) = (ywall(i+2) – 2*ywall(i+1) + ywall(i)) / (xwall(i+1) – xwall(i))^2;
end
cont_dy2_dx2(end-1:end) = NaN;
plot(xwall./xwall(end), dy2_dx2) Greetings everyone, I am trying to fit a Bezier Curve as shown in the paper provide(fig1), but I am unable to get the desired plot(fig 2 and fig3). I have attached the codes as well as the paper and plot. Any help on this is highly appreciated!.
PS: The inflection point is the last point on xarc, i.e, xarc(end) and Nozzle_WithCircARc.m is the main file. It would be better to run the file cell by cell.
%function Nozzle_WithCircARc(G,Me,n,display)
%% Initialize datapoint matrices
clearvars;close all;clc;
G=1.4;
Me = 2.0;
n = 53; % speed index from pucketts paper
display = 0;
Km = zeros(n,n); % K- vlaues (Constant along right running characteristic lines)
Kp = zeros(n,n); % K+ vlaues (Constant along left running characteristic lines)
Theta = zeros(n,n); % Flow angles relative to the horizontal
Mu = zeros(n,n); % Mach angles
M = zeros(n,n); % Mach Numbers
x = zeros(n,n); % x-coordinates
y = zeros(n,n); % y-coordinates
%% Generate the convergent portion of a nozzle
% The inlet height/area and exit height/area(divergent) are same
% Therefore we have same A/A* = 1.6875
% Corresponding to Mach = 0.3722
%% Find NuMax (maximum expansion angle)
[~, NuMax, ~] = PMF(G,Me,0,0);
ThetaMax = NuMax/2;
%%
%ThetaMax0 = (1.0/1.687)^(2/9) * (NuMax/2);
%% Define some flow parameters of originating characteristic lines
dT = ThetaMax/n;
%no_ref = (ThetaMax – ThetaMax_ref)/dT;
ThetaArc(:,1) = (0:dT:ThetaMax);
NuArc = ThetaArc;
KmArc = ThetaArc + NuArc;
[~, ~, MuArc(:,1)] = PMF(G,0,NuArc(:,1),0);
%% Coordinates of wall along curve from throat
y0 = 1; % Define throat half-height
ThroatCurveRadius = 1.5*y0; % Radius of curvature just downstream of the throat
% for larger factors, ywall deviates from A/A* preferred value is 1.1
%L_e = 1.1 * y0 * sind(ThetaMax);
[xarc, yarc] = Arc(ThroatCurveRadius,ThetaArc); % Finds x- and y-coordinates for given theta-values
yarc(:,1) = yarc(:,1) + y0; % Defines offset due to arc being above horizontal
%% Fill in missing datapoint info along first C+ line
% First centerline datapoint done manually
Km(:,1) = KmArc(2:length(KmArc),1);
Theta(:,1) = ThetaArc(2:length(KmArc),1);
Nu(:,1) = Theta(:,1);
Kp(:,1) = Theta(:,1)-Nu(:,1);
M(1,1) = 1.0001;
Nu(1,1) = 0;
Mu(1,1) = 90;
y(1,1) = 0;
x(1,1) = xarc(2,1) + (y(1,1) – yarc(2,1))/tand((ThetaArc(2,1) – MuArc(2,1) – MuArc(2,1))/2);
% Finds the information at interior nodes along first C+ line
for i=2:n
[M(i,1), Nu(i,1), Mu(i,1)] = PMF(G,0,Nu(i,1),0);
s1 = tand((ThetaArc(i+1,1) – MuArc(i+1,1) + Theta(i,1) – Mu(i,1))/2);
s2 = tand((Theta(i-1,1) + Mu(i-1,1) + Theta(i,1) + Mu(i,1))/2);
x(i,1) = ((y(i-1,1)-x(i-1,1)*s2)-(yarc(i+1,1)-xarc(i+1,1)*s1))/(s1-s2);
y(i,1) = y(i-1,1) + (x(i,1)-x(i-1,1))*s2;

end
%% Find flow properties at remaining interior nodes
for j=2:n;
for i=1:n+1-j;

Km(i,j) = Km(i+1,j-1);

if i==1;

else

end

plot(xwall,ywall,’-‘)
axis equal
axis([0 ceil(xwall(length(xwall),1)) 0 ceil(ywall(length(ywall),1))])
hold on
plot(xarc,yarc,’k-‘)

Cart

Cart