clearvars;clc;close all;
clearvars;clc;
Ri = 1.6870;
Rstar = 8.0;
rstar = 1.0;
Rstar_rstar = 8.0;
theta = 25;
a = 2.5;
L = 10 ;
d = 3.0/2.0 * Rstar*tand(theta);
xplus = L-d;
rplus = (5.0*d^2) / (12.0*Rstar) +rstar;
b = 2.0*(xplus + (rplus-Ri)/tand(theta)-a );
e = (Ri - rstar)/tand(theta);
c = e- b/2.0 - (5.0/8.0) *d;
x_cylindrical = linspace(0,a,100);
r_cylindrical = Ri * ones(size(x_cylindrical));
plot(x_cylindrical,r_cylindrical,'.b')
hold on
x_upstream_q = linspace(a, a + b, 100);
r_upstream_q = Ri - (b*tand(theta))/2.0*((x_upstream_q-a)./b).^3.*(2-(x_upstream_q-a)./b);
plot(x_upstream_q,r_upstream_q,'.k')
hold on
x_conical = linspace(a+b,a+b+c,100);
r_conical = Ri + ((a+b)/2.0)*tand(theta) -x_conical.*tand(theta);
plot(x_conical,r_conical,'.r')
x_downstream_q = linspace(a+b+c,L,100);
r_downstream_q = ((L-x_downstream_q).^2)./(12.0*Rstar).*(6.0-((L-x_downstream_q)./d).^2) + rstar;
plot(x_downstream_q,r_downstream_q,'.g')
x_combined = [x_cylindrical, x_upstream_q, x_conical, x_downstream_q];
r_combined = [r_cylindrical, r_upstream_q, r_conical, r_downstream_q];
figure;
plot(x_combined, r_combined, 'LineWidth', 2);
xlabel('x');
ylabel('r(x)');
title('CQ-CQ Nozzle Curve');
grid on;
axis equal;
