From the attached pictures of the equations, the paper simply states that amongst all the curves discussed, the Poly4_Pre curve gives the maximum value of H_e for a given L_e i.e. there is no need to maximize the value of H_e as it is simply the value of the polynomial Poly4_Pre evaluated at x=L_e. Further, this H_e value is used to fit the Bezier curve.
The confusion you had might be due to you assuming the last point At E: X=L_e, Y=H_e as a boundary condition where you assumed the value of H_e to be the value you got from the circular arc.
Omitting this boundary condition, we will get a curve which gives the value of H_e when evaluated at L_e.
This can be done as shown:
Constants
G=1.4;
Gp=G+1;
Gm=G-1;
M = 2;
nu = sqrt(Gp/Gm).*atand(sqrt(Gm*(M.^2-1)/Gp))-atand(sqrt(M.^2-1));
ThetaMax = nu/2;
y0 = 1;
L_e = 3.155*y0*sind(ThetaMax);
% Define x values for plotting
x_values = linspace(0, L_e, 100);
Circular Arc
% Define parameters
y0 = 1; % Initial height, example value
ThetaMax = 15; % Maximum angle in degrees
Le = 3.155 * y0 * sind(ThetaMax); % Length along symmetry axis
% Calculate the radius R based on the slope condition
R = Le / tand(ThetaMax);
% Define the circle center
center_y = y0 + R;
% Calculate corresponding y values for the circular arc
circle_y = center_y - sqrt(R^2 - x_values.^2);
Poly2
syms x a0 a1 a2
% Polynomial of degree 2
Y = a2*x^2 + a1*x + a0;
% Boundary Conditions
% At O: X = 0, Y = h_h
cond1 = subs(Y, x, 0) == y0;
% At O: dY/dX = 0
dY_dx = diff(Y, x);
cond2 = subs(dY_dx, x, 0) == 0;
% At E: dY/dX = tan(theta_max)
cond3 = subs(dY_dx, x, L_e) == tand(ThetaMax);
% Solve for coefficients
coeffs = solve([cond1, cond2, cond3], [a0, a1, a2]);
a0 = coeffs.a0;
a1 = coeffs.a1;
a2 = coeffs.a2;
% Define the polynomial function using the obtained coefficients
Y_poly2 = matlabFunction(a2*x.^2 + a1*x + a0);
poly2_y = Y_poly2(x_values);
Poly3
syms x a0 a1 a2 a3
% Polynomial of degree 4
Y = a3*x^3 + a2*x^2 + a1*x + a0;
% Boundary Conditions
% At O: X = 0, Y = h_h
cond1 = subs(Y, x, 0) == y0;
% At O: dY/dX = 0
dY_dx = diff(Y, x);
cond2 = subs(dY_dx, x, 0) == 0;
% At E: dY/dX = tan(theta_max)
cond3 = subs(dY_dx, x, L_e) == tand(ThetaMax);
% At E: d2Y/dX2 = 0
d2Y_dx2 = diff(dY_dx, x);
cond4 = subs(d2Y_dx2, x, L_e) == 0;
% Solve for coefficients
coeffs = solve([cond1, cond2, cond3, cond4], [a0, a1, a2, a3]);
a0 = coeffs.a0;
a1 = coeffs.a1;
a2 = coeffs.a2;
a3 = coeffs.a3;
% Define the polynomial function using the obtained coefficients
Y_poly3 = matlabFunction(a3*x.^3 + a2*x.^2 + a1*x + a0);
poly3_y = Y_poly3(x_values);
Poly4
syms x a0 a1 a2 a3 a4
% Polynomial of degree 4
Y = a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0;
% Boundary Conditions
% At O: X = 0, Y = h_h
cond1 = subs(Y, x, 0) == y0;
% At O: dY/dX = 0
dY_dx = diff(Y, x);
cond2 = subs(dY_dx, x, 0) == 0;
% At E: dY/dX = tan(theta_max)
cond3 = subs(dY_dx, x, L_e) == tand(ThetaMax);
% At O: d2Y/dX2 = 0
d2Y_dx2 = diff(dY_dx, x);
cond4 = subs(d2Y_dx2, x, L_e) == 0;
% At E: d2Y/dX2 = 0
cond5 = subs(d2Y_dx2, x, 0) == 0;
% Solve for coefficients
coeffs = solve([cond1, cond2, cond3, cond4, cond5], [a0, a1, a2, a3, a4]);
a0 = coeffs.a0;
a1 = coeffs.a1;
a2 = coeffs.a2;
a3 = coeffs.a3;
a4 = coeffs.a4;
Y_poly4 = matlabFunction(a4*x.^4 + a3*x.^3 + a2*x.^2 + a1*x + a0);
poly4_y = Y_poly4(x_values);
Poly4Pre
syms x a0 a1 a2 a3 a4
% Polynomial of degree 4
Y = a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0;
% Boundary Conditions
% At O: X = 0, Y = h_h
cond1 = subs(Y, x, 0) == y0;
% At O: dY/dX = 0
dY_dx = diff(Y, x);
cond2 = subs(dY_dx, x, 0) == 0;
% At E: dY/dX = tan(theta_max)
cond3 = subs(dY_dx, x, L_e) == tand(ThetaMax);
% At E: d2Y/dX2 = 0
d2Y_dx2 = diff(dY_dx, x);
cond4 = subs(d2Y_dx2, x, L_e) == 0;
% Solve for coefficients
coeffs = solve([cond1, cond2, cond3, cond4], [a0, a1, a2, a3, a4]);
a0 = coeffs.a0;
a1 = coeffs.a1;
a2 = coeffs.a2;
a3 = coeffs.a3;
a4 = coeffs.a4;
Y_poly4Pre = matlabFunction(a4*x.^4 + a3*x.^3 + a2*x.^2 + a1*x + a0);
poly4Pre_y = Y_poly4Pre(x_values);
H_e = poly4Pre_y(end)
Bézier
% Control points as symbolic variables
syms x1 y1 x2 y2 real
% Define control points
P0 = [0, y0];
P1 = [x1, y1];
P2 = [x2, y2];
P3 = [L_e, H_e];
% Bézier curve as a symbolic function of t
syms t
B = (1-t)^3 * P0 + 3*(1-t)^2 * t * P1 + 3*(1-t) * t^2 * P2 + t^3 * P3;
% Separate x and y components
Bx = B(1);
By = B(2);
% First derivatives
dBx_dt = diff(Bx, t);
dBy_dt = diff(By, t);
% Slope (dy/dx) at t
slope = dBy_dt / dBx_dt;
% double derivative
DD = diff(slope,t)/dBx_dt;
% Boundary conditions
eq1 = subs(By, t, 0) == y0; % At t = 0, y = y0
eq2 = subs(slope, t, 0) == 0; % At t = 0, dy/dx = 0
eq3 = subs(slope, t, 1) == tand(ThetaMax); % At t = 1, dy/dx = tan(thetaMax)
eq4 = subs(By, t, 1) == H_e; % At t = 1, y = He
eq5 = subs(DD,t,0) == 0;
eq6 = subs(DD,t,1) == 0;
% Solve for the unknowns
sol = solve([eq1, eq2, eq3, eq4, eq5, eq6], [x1, y1, x2, y2]);
% Display the solution
x1_sol = double(sol.x1);
y1_sol = double(sol.y1);
x2_sol = double(sol.x2);
y2_sol = double(sol.y2);
% Substitute the solutions back into the control points
P1_sol = [x1_sol, y1_sol];
P2_sol = [x2_sol, y2_sol];
% Plot the Bézier curve
t_vals = linspace(0, 1, 100);
B_vals = subs(B, {x1, y1, x2, y2}, {x1_sol, y1_sol, x2_sol, y2_sol});
B_vals = subs(B_vals,t,t_vals');
% Extract x and y coordinates
x_vals = B_vals(:, 1);
y_vals = B_vals(:, 2);
Plotting All
figure;
hold on;
plot(x_values, circle_y, 'DisplayName', 'Circular Arc');
plot(x_values, poly2_y, 'DisplayName', 'Poly2');
plot(x_values, poly3_y , 'DisplayName', 'Poly3');
plot(x_values, poly4_y , 'DisplayName', 'Poly4');
plot(x_values, poly4Pre_y , 'DisplayName', 'Poly4Pre');
plot(x_vals, y_vals, 'DisplayName', 'Bezier');
xlabel('X');
ylabel('Y');
title('Comparing Plots');
legend('show', 'Location', 'northwest');
grid on;
hold off;
This graph looks similar to the results shown in the paper which you were referring to.
