%% Mission and System Design
% Main Source: Heister, Rocket Propulsion

% Orbital Mechanics
z=100; % km, target altitude
mu=4902.8; % km^3*s^-2, moon GM — https://ssd.jpl.nasa.gov/astro_par.html
r=1737; % km, moon radius
v_c=1000*sqrt(mu/(r+z)); % m/s, circular orbital velocity

% Design Parameter Calculations
dV_g_frac=0.20; % gravity loss fraction — Heister Problem 2.8, p.42 suggests 20% for lunar ascent
dV_g=dV_g_frac*2*v_c; % m/s, gravity loss
dV=2*v_c+dV_g; % m/s, delta V (required for lunar surface to orbit and back)
ISP_lit=336; % s, vacuum specific impulse for MMH/NTO from literature — http://www.astronautix.com/n/n2o4mmh.html
ISP_old=341.9974; % s, vacuum ISP found in second section after first design iteration
ISP_new=338.6933; % s, vacuum ISP found in seventh section after heat transfer analysis
of_lit=1.65; % http://www.astronautix.com/n/n2o4mmh.html#:~:text=Subtopics,First%20flight%201996.
of_new=2.2727; % O/F ratio found in second section after first iteration
g0=9.807; % m*s^-2, acceleration unit conversion constant
MR=exp(dV/(ISP_new*g0)); % mass ratio, derived from rocket equation

m_payload=20000; % kg, specified payload mass

m_inert=5000; % kg, inert mass estimation

m_prop=MR*(m_inert+m_payload)-m_inert-m_payload; % kg, propellant mass derived from MR definition

m_wet=m_payload+m_inert+m_prop;

rho_MMH=874; % kg*m^-3, density at room temp — https://pubchem.ncbi.nlm.nih.gov/compound/Methylhydrazine#section=Density
rho_NTO=1448; % kg*m^-3, density near boiling point
m_MMH=m_prop/(1+of_new); % kg, mass of MMH
m_NTO=of_new*(m_prop/(1+of_new)); % kg, mass of NTO
V_prop=m_MMH/rho_MMH+m_NTO/rho_NTO; % m^3, total volume of propellant
V_prop_tank=10; % m^3, tank volume estimation
n_tank=V_prop/V_prop_tank; % number of prop tanks

% Results
fprintf('\n--- Mission and System Design ---\n');
fprintf('Delta V                    = %.2f m/s\n',dV);
fprintf('Target Velocity            = %.2f m/s\n',v_c);
fprintf('Mass Ratio                 = %.2f\n',MR);
fprintf('Payload Mass               = %.2f kg\n',m_payload);
fprintf('Inert Mass                 = %.2f kg\n',m_inert);
fprintf('Propellant Mass            = %.2f kg\n',m_prop);
fprintf('Propellant Volume          = %.2f m^3\n',V_prop);
fprintf('Wet Mass                   = %.2f kg\n',m_wet);
fprintf('Number of Propellant Tanks = %.0f\n',n_tank);

%% Liquid Engine Design Overview

% Thrust
eta=1.3; % thrust-to-weight ratio
g_moon_surface=1000*mu/r^2; % m*s^-2, gravity at lunar surface
F_total=eta*m_wet*g_moon_surface;
n_engines=4;
F_engine=F_total/n_engines;

% CEA Simulation Properties & Respective Core Calculations
load('thermo_lib.mat')
T_prop=298.15; % K, propellant tanks at room temp
% sweep over chamber pressure values
Pc_vals(1,:)=linspace(100,200); % psia
Pc_Pa_vals=Pc_vals*6894.76; % Pa
C_star_vals=zeros(2,length(Pc_vals));
ISP_vals=zeros(2,length(Pc_vals));
m_dot_vals=zeros(2,length(Pc_vals));
m_dot_engine_vals=zeros(2,length(Pc_vals));
At_vals=zeros(2,length(Pc_vals));
Cf_vals=zeros(2,length(Pc_vals));
t_burn_vals=zeros(2,length(Pc_vals));
for i=1:length(Pc_vals)
    CEA_data=CEA('problem','rocket','equilibrium','o/f',1.65,'p,psia',Pc_vals(i),'supson',50,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','CSTAR','Ivac','end');
    C_star_vals(1,i)=CEA_data.output.eql.cstar(2); % m/s, Cstar
    ISP_vals(1,i)=CEA_data.output.eql.isp_vac(3); % s, vacuum ISP
    m_dot_vals(1,i)=F_total/(ISP_vals(1,i)*g0); % kg/s, total mass flow rate
    m_dot_engine_vals(1,i)=m_dot_vals(1,i)/n_engines; % kg/s, mass flow rate per engine
    At_vals(1,i)=m_dot_engine_vals(1,i)*C_star_vals(1,i)/Pc_Pa_vals(1,i); % m^2, throat area
    Cf_vals(1,i)=F_engine/(Pc_Pa_vals(1,i)*At_vals(1,i)); % coefficient of thrust
    t_burn_vals(1,i)=(m_prop/m_dot_vals(1,i))/60; % min, burn time
end
figure
subplot(2,3,1); plot(Pc_vals,C_star_vals(1,:),'LineWidth',2); xlabel('Chamber Pressure (psia)'); ylabel('Cstar (m/s)');
subplot(2,3,2); plot(Pc_vals,ISP_vals(1,:),'LineWidth',2); xlabel('Chamber Pressure (psia)'); ylabel('Vacuum ISP (s)');
subplot(2,3,3); plot(Pc_vals,m_dot_engine_vals(1,:),'LineWidth',2); xlabel('Chamber Pressure (psia)'); ylabel('M_{dot} per Engine (kg/s)');
subplot(2,3,4); plot(Pc_vals,At_vals(1,:),'LineWidth',2); xlabel('Chamber Pressure (psia)'); ylabel('Throat Area (m^2)');
subplot(2,3,5); plot(Pc_vals,Cf_vals(1,:),'LineWidth',2);  xlabel('Chamber Pressure (psia)'); ylabel('Coefficient of Thrust');
subplot(2,3,6); plot(Pc_vals,t_burn_vals(1,:),'LineWidth',2); xlabel('Chamber Pressure (psia)'); ylabel('Burn Time (min)');
% sweep over chamber pressures and O/F ratios via a contour plot
Pc_grid_vals=linspace(100, 200, 20);
of_grid_vals=linspace(1, 2.5, 20);
ISP_Matrix = zeros(length(Pc_grid_vals), length(of_grid_vals));
for i = 1:length(Pc_grid_vals)
    for j = 1:length(of_grid_vals)
        CEA_data=CEA('problem','rocket','equilibrium','o/f',of_grid_vals(j),'p,psia',Pc_grid_vals(i),'supson',50,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','CSTAR','Ivac','end');
        ISP_Matrix(i,j)=CEA_data.output.eql.isp_vac(3);
    end
end
figure
contourf(of_grid_vals,Pc_grid_vals,ISP_Matrix,50,'LineColor','none');
colorbar;
xlabel('O/F Ratio'); ylabel('Chamber Pressure (psia)'); title('Vacuum ISP Contour Map');
% set pressure chamber
Pc=100; % psia
Pc_Pa=Pc*6894.76;
% sweep over O/F ratios
of_vals=linspace(1,2.5,100);
for i=1:length(Pc_vals)
    CEA_data=CEA('problem','rocket','equilibrium','o/f',of_vals(i),'p,psia',Pc,'supson',50,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','CSTAR','Ivac','end');
    C_star_vals(2,i)=CEA_data.output.eql.cstar(2); % m/s, Cstar
    ISP_vals(2,i)=CEA_data.output.eql.isp_vac(3); % s, vacuum ISP
    m_dot_vals(2,i)=F_total/(ISP_vals(2,i)*g0); % kg/s, total mass flow rate
    m_dot_engine_vals(2,i)=m_dot_vals(2,i)/n_engines; % kg/s, mass flow rate per engine
    At_vals(2,i)=m_dot_engine_vals(2,i)*C_star_vals(2,i)/Pc_Pa; % m^2, throat area
    Cf_vals(2,i)=F_engine/(Pc_Pa*At_vals(2,i)); % coefficient of thrust
    t_burn_vals(2,i)=(m_prop/m_dot_vals(2,i))/60; % min, burn time
end
figure
subplot(2,3,1); plot(of_vals,C_star_vals(2,:),'LineWidth',2); xlabel('O/F Ratio'); ylabel('Cstar (m/s)');
subplot(2,3,2); plot(of_vals,ISP_vals(2,:),'LineWidth',2); xlabel('O/F Ratio'); ylabel('Vacuum ISP (s)');
subplot(2,3,3); plot(of_vals,m_dot_engine_vals(2,:),'LineWidth',2); xlabel('O/F Ratio'); ylabel('M_{dot} per Engine (kg/s)');
subplot(2,3,4); plot(of_vals,At_vals(2,:),'LineWidth',2); xlabel('O/F Ratio'); ylabel('Throat Area (m^2)');
subplot(2,3,5); plot(of_vals,Cf_vals(2,:),'LineWidth',2);  xlabel('O/F Ratio'); ylabel('Coefficient of Thrust');
subplot(2,3,6); plot(of_vals,t_burn_vals(2,:),'LineWidth',2); xlabel('O/F Ratio'); ylabel('Burn Time (min)');
% set O/F ratio & ISP for iterations with section 1
[~,of_index]=max(ISP_vals(2,:)); % s, vacuum ISP and its maximizing index
of=of_vals(of_index); % O/F ratio at ISP maximizing index
% final calculations
ISP=ISP_new; % from section 7 heat transfer analysis
C_star=C_star_vals(2,of_index); % m/s, Cstar
m_dot=F_total/(ISP*g0); % kg/s, total mass flow rate
m_dot_engine=m_dot/n_engines; % kg/s, mass flow rate per engine
At=m_dot_engine*C_star/Pc_Pa; % m^2, throat area
Cf=F_engine/(Pc_Pa*At); % coefficient of thrust
t_burn=(m_prop/m_dot)/60; % min, burn time

% Throttling
eta_surface=1; % set thrust to weight ratio equal to one for safe, hover maneuver landing
m_wet_surface=(m_inert+m_payload)*exp(dV/(2*ISP*g0)); % kg, wet mass at surface -- rocket equation using half ∆V
F_min=m_wet_surface*g_moon_surface*eta_surface; % N, minimum total thrust
F_min_engine=F_min/n_engines;
F_ratio=F_min/F_total;
% shut two engines off to get to 50%, throttle each engine to 86% of max

epsilon=50; % expansion ratio, standard for vacuum conditions

% Results
fprintf('\n--- Liquid Engine Design Overview ---\n');
fprintf('Maximum Total Thrust           = %.2f N\n',F_total);
fprintf('Minimum Total Thrust           = %.2f N\n',F_min);
fprintf('Min:Max Thrust Ratio           = %.2f\n',F_ratio);
fprintf('Max thrust Per Engine          = %.2f N\n',F_engine);
fprintf('Min thrust Per Engine          = %.2f N\n',F_min_engine);
fprintf('O/F Ratio                      = %.2f\n',of);
fprintf('Number of Engines              = %.0f\n',n_engines);
fprintf('Specific (Vacuum) Impulse, ISP = %.0f s\n',ISP);
fprintf('Total Mass Flow Rate           = %.2f kg/s\n',m_dot);
fprintf('Mass Flow Rate per Engine      = %.2f kg/s\n',m_dot_engine);
fprintf('Total Burn Time                = %.2f min\n',t_burn);
fprintf('Chamber Pressure               = %.0f Pa = %.0f psi\n',Pc_Pa,Pc);
fprintf('Characteristic Velocity, C*    = %.0f m/s\n',C_star)
fprintf('Throat Area                    = %.4f m^2\n',At);
fprintf('Coefficient of Thrust          = %.2f\n',Cf);
fprintf('Expansion Ratio                = %.0f\n',epsilon);

%% Pressurization / Propellant

% Helium Sizing
% assume ideal gas, isentropic
dP_inj=0.20*Pc_Pa; % Pa, assume 20% pressure drop across injector
dP_lines=0.05*Pc_Pa; % Pa, assume 5% pressure drop across feed lines
dP_misc=0.05*Pc_Pa;
P_prop_tank=Pc_Pa+dP_inj+dP_misc+dP_lines;
R_He=2077.1; % J*kg^-1*K^-1, R constant for Helium
gamma_He=5/3; % Cp/Cv of Helium
m_He=P_prop_tank*V_prop/(R_He*T_prop); % kg, mass Helium required to fill prop tank
rho_He=0.1785; % kg*m^-3, density of helium at 0˚C
P_He_Pa=4500*6894.76; % Pa, pressure of helium tank
V_He_tank=m_He*R_He*T_prop/P_He_Pa; % m^3, helium tank volume
r_He_tank=(V_He_tank*3/4)^(1/3); % m, radius of helium tank
sig_y_tank=880; % MPa, yield strength of Ti-6Al-4V, material of helium tank - https://www.amardeepsteel.com/blog/titanium-alloys-6-4.html
t_He_tank=P_He_Pa*10^-6*r_He_tank*1000/(2*sig_y_tank); % mm, thickness of helium tank

% Propellant Conditions
P_injector_inlet=Pc_Pa+dP_inj; % Pa, pressure right before injector

% Results
fprintf('\n--- Pressurization and Pressurant Conditions ---\n');
fprintf('Helium Tank Pressure     = %.0f psi\n',P_He_Pa/6894.76);
fprintf('Helium Tank Volume       = %.2f m^3\n',V_He_tank);
fprintf('Helium Tank Diameter     = %.2f m\n',r_He_tank*2);
fprintf('Helium Tank Thickness    = %.2f mm\n',t_He_tank);
fprintf('Helium Mass              = %.2f kg\n',m_He)
fprintf('Propellant Tank Pressure = %.2f psi\n',P_prop_tank/6894.76);

fprintf('\n--- Propellant Conditions Entering Injector ---\n');
fprintf('Phase                    = Liquid\n')
fprintf('Temperature              = %.2f K\n',T_prop)
fprintf('Inlet Pressure           = %.2f Psi\n',P_injector_inlet/6894.76)

%% Injector Design

% Inconel 625
% NTO/MMH --> unlike doublet impinging element injector because:
    % simliar volume flows
    % large interfacial contact, beneficial for hypergolics

m_dot_MMH_engine=m_dot_engine/(1+of); % kg/s, mass flow of fuel per engine
m_dot_NTO_engine=m_dot_engine*of/(1+of); % kg/s, mass flow of oxidizer per engine

% Orifice Area Calculations
% main orifices - this is the main fuel stream
percent_fuel_film=0.05;
percent_fuel_core=1-percent_fuel_film;
m_dot_MMH_core=percent_fuel_core*m_dot_MMH_engine;
m_dot_NTO_core=percent_fuel_core*m_dot_NTO_engine;
Cd=0.7; % assume somewhat sharp corners on orifices to stay conservative
N_main_doublets=50; % assume 50 main doublets
m_dot_NTO_orifice_main=m_dot_NTO_core/N_main_doublets; % mass flow NTO per main orifice
m_dot_MMH_orifice_main=m_dot_MMH_core/N_main_doublets; % mass flow MMH per main orifice
A_NTO_main=m_dot_NTO_orifice_main/(Cd*sqrt(2*rho_NTO*dP_inj)); % m^2, NTO main orifice area
A_MMH_main=m_dot_MMH_orifice_main/(Cd*sqrt(2*rho_MMH*dP_inj)); % m^2, MMH main orifice area
D_NTO_main=sqrt(A_NTO_main*4/pi); % m, NTO orifice diameter
D_MMH_main=sqrt(A_MMH_main*4/pi); % m, MMH orifice diameter
% cooling orifices - for film cooling
m_dot_film=percent_fuel_film*(m_dot_MMH_engine+m_dot_NTO_engine); % kg/s, total mass flow of propellant for film cooling
of_film=0.5; % calculated in heat xfr section
m_dot_MMH_film=m_dot_film/(1+of_film); % kg/s, mass flow of fuel per engine
m_dot_NTO_film=m_dot_film*of_film/(1+of_film); % kg/s, mass flow of oxidizer per engine
N_film_doublets=10; % assume 10 film doublets 
m_dot_NTO_orifice_film=m_dot_NTO_film/N_film_doublets; % mass flow NTO per film orifice
m_dot_MMH_orifice_film=m_dot_MMH_film/N_film_doublets; % mass flow MMH per film orifice
A_NTO_film=m_dot_NTO_orifice_film/(Cd*sqrt(2*rho_NTO*dP_inj)); % m^2, NTO main orifice area
A_MMH_film=m_dot_MMH_orifice_film/(Cd*sqrt(2*rho_MMH*dP_inj)); % m^2, MMH main orifice area
D_NTO_film=sqrt(A_NTO_film*4/pi); % m, NTO orifice diameter
D_MMH_film=sqrt(A_MMH_film*4/pi); % m, MMH orifice diameter

% Orifice Angle Calculations
% main angles
v_NTO_main=m_dot_NTO_orifice_main/(rho_NTO*A_NTO_main); % m/s, NTO injector velocity
v_MMH_main=m_dot_MMH_orifice_main/(rho_MMH*A_MMH_main); % m/s, MMH injector velocity
sin_ratio_main=m_dot_NTO_orifice_main*v_NTO_main/(m_dot_MMH_orifice_main*v_MMH_main); % ratio of the sine of the angle of the MMH orifice to that of the NTO orifice
% sum_a_vals=a_MMH+a_NTO --> a_MMH=sum_a_vals-a_NTO (1)
% sind(a_MMH)=sin_ratio*sind(a_NTO) (2)
% (1) into (2), solve for a_NTO --> trig identities -->
% a_NTO=atand(sind(sum_a_vals)/(sin_ratio+cos(sum_a_vals)))
% choose sum_a = 45˚ to keep the collisions furthest away from injector
% plate and mitigate blowback
sum_a=45;
a_NTO_main=atand(sind(sum_a)/(sin_ratio_main+cosd(sum_a))); % degrees, chosen angle of the NTO orifice
a_MMH_main=sum_a-a_NTO_main; % degrees, chosen angle of the MMH orifice
% cooling angles
v_NTO_film=m_dot_NTO_film/(rho_NTO*A_NTO_film); % m/s, NTO injector velocity
v_MMH_film=m_dot_MMH_film/(rho_MMH*A_MMH_film); % m/s, MMH injector velocity
sin_ratio_film=m_dot_NTO_film*v_NTO_film/(m_dot_MMH_film*v_MMH_film); % ratio of the sine of the angle of the MMH orifice to that of the NTO orifice
a_NTO_film=atand(sind(sum_a)/(sin_ratio_film+cosd(sum_a))); % degrees, chosen angle of the NTO orifice
a_MMH_film=sum_a-a_NTO_film; % degrees, chosen angle of the MMH orifice

% Details from CAD
% main
dist_main=9; % mm, distance from center of MMH orifice to that of NTO orifice
dist_impinge_main=10.71; % mm, distance to impingement
% cooling
dist_film=9; % mm, distance from center of MMH orifice to that of NTO orifice
dist_impinge_film=10.46; % mm, distance to impingement

% Results
fprintf('\n--- Injector Design ---\n');
fprintf('O/F Ratio                      = %.2f\n',of);
fprintf('MMH Mass Flow Per Engine       = %.2f kg/s\n',m_dot_MMH_engine);
fprintf('NTO Mass Flow Per Engine       = %.2f kg/s\n', m_dot_NTO_engine);
fprintf('Percent of Propellant for Film = %.2f\n',percent_fuel_film);
fprintf('MMH Main Orifice Diameter      = %.2f mm\n',D_MMH_main*10^3);
fprintf('NTO Main Orifice Diameter      = %.2f mm\n',D_NTO_main*10^3);
fprintf('MMH Main Orifice Angle         = %.2f degrees\n',a_MMH_main);
fprintf('NTO Main Orifice Angle         = %.2f degrees\n',a_NTO_main);
fprintf('MMH Film Orifice Diameter      = %.2f mm\n',D_MMH_film*10^3);
fprintf('NTO Film Orifice Diameter      = %.2f mm\n',D_NTO_film*10^3);
fprintf('MMH Film Orifice Angle         = %.2f degrees\n',a_MMH_film);
fprintf('NTO Film Orifice Angle         = %.2f degrees\n',a_NTO_film);
fprintf('Distance between Orifices      = %.0f mm\n',dist_film);
fprintf('Main Impingement Distance      = %.2f mm\n', dist_impinge_main);
fprintf('Film Impingement Distance      = %.2f mm\n',dist_impinge_film);
fprintf('Number of Main Doublets        = %.0f\n',N_main_doublets);
fprintf('Number of Film Doublets        = %.0f\n',N_film_doublets);
%% Combustion Modeling
% CEA analysis integrated into each part
CEA_model=CEA('problem','rocket','eql','o/f',of,'p,psia',Pc,'supson',epsilon,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','end');
gamma_model=CEA_model.output.eql.gamma(1); % specific heat ratio at chamber
MW_model=CEA_model.output.eql.mw(1); % kg/kmol, combined molecular weight of propellant at chamber
T_af_model=CEA_model.output.eql.temperature(1); % K, adiabatic flame temp

fprintf('\n--- Combustion Properties ---\n')
fprintf('Specific Heat Ratio  = %.2f\n',gamma_model);
fprintf('Molecular Weight     = %.2f\n',MW_model);
fprintf('Adiabatic Flame Temp = %.2f\n',T_af_model);
%% Chamber + Nozzle Design

L_star=35*0.0254; % m, assumption for characteristic length - Huzel Table 4-1 p. 87
Vc=L_star*At; % m^2, chamber volume
CR=2.5; % chamber contraction ratio used throughout Huzel
Ac=CR*At; % m^2, chamber cross-sectional area
Dc=sqrt(4*Ac/pi); % m, chamber diameter
Dt=sqrt(4*At/pi); % m, throat diameter
Lc=Vc/Ac; % m, chamber length

% 80% Bell Nozzle Design (Rao parabolic approximation) - Huzel & http://www.aspirespace.org.uk/downloads/Thrust%20optimised%20parabolic%20nozzle.pdf
Lf=0.8; % length fraction, 80% of standard conical nozzle
De=sqrt(epsilon)*Dt; % m, exit diameter
Ln=0.8*(sqrt(epsilon)-1)*(Dt/2)/tand(15); % m, nozzle length
Rt=Dt/2; % m, throat radius
Re=De/2; % m, exit radius
% Figure 4-14 Huzel, for Lf=0.8 yields:
theta_n=32.5; % degrees, initial parabola angle
theta_e=7.5; % degrees, final parabola angle
% throat entrant
thetas_ent=linspace(-135,-90,100); % parameter range
x_ent=1.5*Rt*cosd(thetas_ent); % eqn 1.1
y_ent=1.5*Rt*sind(thetas_ent)+1.5*Rt+Rt; % eqn 1.2
% chamber
x_c=linspace(1.5*Rt*cosd(-135)-Lc,1.5*Rt*cosd(-135),100);
y_c=Dc/2*ones(1,100);
% throat exit
thetas_ex=linspace(-90,theta_n-90,100); % parameter range
x_ex=0.382*Rt*cosd(thetas_ex); % eqn 2.1
y_ex=0.382*Rt*sind(thetas_ex)+0.382*Rt+Rt; % eqn 2.2
% Bell (Bezier curve)
t_bell=linspace(0,1,100); % parameter range
Nx=0.382*Rt*cosd(theta_n-90); % point N x-coordinate, set eqn 2.1 parameter to theta_n-90
Ny=0.382*Rt*sind(theta_n-90)+0.382*Rt+Rt; % point N y-coordinate, set eqn 2.2 parameter to theta_n-90
Ex=Ln; % exit x-coordinate
Ey=Re; % exit y-coordinate
m1=tand(theta_n); % gradient 1
m2=tand(theta_e); % gradient 2
C1=Ny-m1*Nx; % intercept 1
C2=Ey-m2*Ex; % intercept 2
Qx=(C2-C1)/(m1-m2); % point Q x-coordinate
Qy=(m1*C2-m2*C1)/(m1-m2); % point Q y-coordinate
x_bell=(1-t_bell).^2.*Nx+2.*(1-t_bell).*t_bell.*Qx+t_bell.^2.*Ex;
y_bell=(1-t_bell).^2.*Ny+2.*(1-t_bell).*t_bell.*Qy+t_bell.^2.*Ey;

% plot nozzle
x_nozzle=[x_c,x_ent,x_ex,x_bell];
y_nozzle_upper=[y_c,y_ent,y_ex,y_bell];
y_nozzle_lower=-y_nozzle_upper;
figure; hold on; axis equal; axis tight; grid on; grid minor;
plot(x_nozzle,y_nozzle_upper,x_nozzle,y_nozzle_lower,'Color','blue','LineWidth',2);
title('Nozzle Geometry'); xlabel('Length (m)'); ylabel('Radius (m)')

% local flow properties
R_local=y_nozzle_upper; % m, local radius
A_local=pi.*R_local.^2; % m^2, local area
eps_local=A_local./At; % local area ratio
eps_index=find(x_nozzle==0,1); % find index of throat
eps_subson=eps_local(1:eps_index); % subsonic area ratios
eps_supson=eps_local(eps_index+1:end); % supersonic area ratios
M_local=zeros(1,length(eps_local));
P_local=zeros(1,length(eps_local));
T_local=zeros(1,length(eps_local));
a_local=zeros(1,length(eps_local));
for i=1:length(eps_subson)
    CEA_nozzle=CEA('problem','rocket','eql','o/f',of,'p,psia',Pc,'subson',eps_subson(i),'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','end'); % CEA simulation
    M_local(i)=CEA_nozzle.output.eql.mach(3); % local Mach number
    P_local(i)=CEA_nozzle.output.eql.pressure(3)*14.5038; % psi, local pressure
    T_local(i)=CEA_nozzle.output.eql.temperature(3); % K, local temp
    a_local(i)=CEA_nozzle.output.eql.sonvel(3); % m/s, local sonic velocity
end
for i=1:length(eps_supson)
    CEA_nozzle=CEA('problem','rocket','eql','o/f',of,'p,psia',Pc,'supson',eps_supson(i),'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','end'); % CEA simulation
    M_local(i+eps_index)=CEA_nozzle.output.eql.mach(3); % local Mach number
    P_local(i+eps_index)=CEA_nozzle.output.eql.pressure(3)*14.5038; % psi, local pressure
    T_local(i+eps_index)=CEA_nozzle.output.eql.temperature(3); % K, local temp
    a_local(i+eps_index)=CEA_nozzle.output.eql.sonvel(3); % m/s, local sonic velocity
end
v_local=M_local.*a_local; % m/s, local velocity

% plot local flow properties
figure
subplot(4,1,1);
plot(x_nozzle,M_local,'Color','Blue','LineWidth',2);
xlabel('Length (m)'); ylabel('Mach Number')
subplot(4,1,2)
plot(x_nozzle,P_local,'Color','Blue','LineWidth',2);
xlabel('Length (m)'); ylabel('Pressure (psia)')
subplot(4,1,3)
plot(x_nozzle,T_local,'Color','Blue','LineWidth',2);
xlabel('Length (m)'); ylabel('Temperature (K)');
subplot(4,1,4)
plot(x_nozzle,v_local,'Color','Blue','LineWidth',2);
xlabel('Length (m)'); ylabel('Velocity (m/s)');

% hoop stress, thin-wall assumption
sig_hoop_local=P_local.*6894.76.*R_local/0.001; % 20 mm thickness

figure; hold on
plot(x_nozzle,sig_hoop_local*10^-6,'LineWidth',2,'DisplayName','Hoop Stress');
skirt_idx=find(eps_local>=20,1,'first');
xline(x_nozzle(skirt_idx),'LineWidth',2,'Color','r','DisplayName','Nozzle Start');
xlabel('Length (m)'); ylabel('Hoop Stress (MPa)'); legend('show')

fprintf('\n--- Nozzle Design ---\n');
fprintf('80 Percent Bell Nozzle\n');
fprintf('\n--- Chamber Design ---\n');
fprintf('Characteristic Length = %.2f cm\n',L_star*100);
fprintf('Contraction Ratio     = %.2f\n',CR);
fprintf('Chamber Length        = %.2f cm\n',Lc*100);
fprintf('Chamber Diameter      = %.2f cm\n',Dc*100);

%% Heat Transfer Analysis

% Tw calculation function, incorporating ISP tradeoffs - adapted from HW 4
function [Tw,ISP_true]=Tw_calc(of_film,Pc,eps_local,T_prop,D_local,emissivity,percent_fuel_film,m_dot_MMH_engine,m_dot_NTO_engine)
    % CEA on film cooling stream
    CEA_Tw=CEA('problem','rocket','eql','o/f',of_film,'p,psia',Pc,'supson',eps_local,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','transport','end'); % CEA simulation
    if eps_local==1 % workaround to solve for Tw at throat
        idx=2;
    else
        idx=3;
    end
    Pr=CEA_Tw.output.eql.prandtl.eql(idx); % Prandtl number
    Cp=CEA_Tw.output.eql.cp(idx)*1000; % J*kg^-1*K^-1, specific heat constant pressure
    mu=CEA_Tw.output.eql.viscosity(idx)*10^-4; % Pa*s, viscosity
    rho=CEA_Tw.output.eql.density(idx);
    M=CEA_Tw.output.eql.mach(idx); % Mach number
    a=CEA_Tw.output.eql.sonvel(idx); % m/s, sonic velocity
    v=M*a; % m/s, velocity
    mu_0=CEA_Tw.output.eql.viscosity(1)*10^-4; % Pa*s, chamber viscosity
    T=CEA_Tw.output.eql.temperature(idx); % K, temperature
    gamma=CEA_Tw.output.eql.gamma(idx); % specific heat ratio
    r=Pr^(1/3); % recovery constant
    Tr=T*(1+(gamma-1)/2*r*M^2); % K, recovery temp

    % Newton-Raphson solver for wall temp
    Tw=T-10; % initial guess, wall temp 10 K cooler than film cooling stream temp
    for n=1:30
        rho_am=0.5*rho*(1+T/Tw); % arithmetic mean density
        mu_w=mu*(Tw/T)^0.68125; % wall viscosity
        mu_am=(mu+mu_w)/2; % arithmetic mean viscosity
        hg=0.026/D_local^0.2 * Cp*mu^0.2/Pr^0.6 * (rho*v)^0.8 * rho_am/rho * (mu_am/mu_0)^0.2; % Bartz equation for heat transfer coefficient
        f=hg*(Tr-Tw)-emissivity*5.67*10^-8*Tw^4; % set convection minus radiation equal to f; we want f to eventually be zero 
        dfdT=-hg-4*emissivity*5.67*10^-8*Tw^3;
        Tw_new=Tw+0.5*(-f/dfdT);
        Tw_new=max(Tw_new,300);
        Tw_new=min(Tw_new,Tr);
        if abs(Tw_new-Tw)<0.1
            Tw=Tw_new;
            break
        end
        Tw=Tw_new;
    end

    % ISP trade-off
    percent_fuel_core=1-percent_fuel_film;
    
    % mass fraction of core
    m_dot_NTO_core=percent_fuel_core*m_dot_NTO_engine;
    m_dot_MMH_core=percent_fuel_core*m_dot_MMH_engine;
    m_dot_core=m_dot_NTO_core+m_dot_MMH_core;
    m_dot=m_dot_NTO_engine+m_dot_MMH_engine;
    mass_frac_core=m_dot_core/m_dot;

    % mass fraction of film
    m_dot_NTO_film=percent_fuel_film*m_dot_NTO_engine;
    m_dot_MMH_film=percent_fuel_film*m_dot_MMH_engine;
    m_dot_film=m_dot_NTO_film+m_dot_MMH_film;
    mass_frac_film=m_dot_film/m_dot;

    CEA_core=CEA('problem','rocket','eql','o/f',2.27,'p,psia',Pc,'supson',eps_local,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','transport','end'); % CEA core simulation
    ISP_core=CEA_core.output.eql.isp_vac(3);
    ISP_film=CEA_Tw.output.eql.isp_vac(3);
    ISP_true=mass_frac_core*ISP_core+mass_frac_film*ISP_film;
end

% chamber Tw function
function Tw_actual=Tw_chamber_calc(of_film,Pc,eps_local,T_prop,D_local,k_phenol,t_phenol,k_carbon,t_carbon,CR)
    % CEA on film cooling stream
    CEA_Tw=CEA('problem','rocket','eql','o/f',of_film,'p,psia',Pc,'subson',eps_local,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','transport','end'); % CEA simulation
    if eps_local==CR % workaround to solve for Tw at chamber
        idx=1;
    elseif eps_local==1
        idx=2;
    else
        idx=3;
    end
    Pr=CEA_Tw.output.eql.prandtl.eql(idx); % Prandtl number
    Cp=CEA_Tw.output.eql.cp(idx)*1000; % J*kg^-1*K^-1, specific heat constant pressure
    mu=CEA_Tw.output.eql.viscosity(idx)*10^-4; % Pa*s, viscosity
    rho=CEA_Tw.output.eql.density(idx);
    M=CEA_Tw.output.eql.mach(idx); % Mach number
    a=CEA_Tw.output.eql.sonvel(idx); % m/s, sonic velocity
    v=M*a; % m/s, velocity
    mu_0=CEA_Tw.output.eql.viscosity(1)*10^-4; % Pa*s, chamber viscosity
    T=CEA_Tw.output.eql.temperature(idx); % K, temperature
    gamma=CEA_Tw.output.eql.gamma(idx); % specific heat ratio
    r=Pr^(1/3); % recovery constant
    Tr=T*(1+(gamma-1)/2*r*M^2); % K, recovery temp (in chamber, becomes stagnation temp)

    % Newton-Raphson solver for wall temp
    Tw=T-10; % initial guess, wall temp 10 K cooler than film cooling stream temp
    T_ext=400; % K, keep outside surface of chamber at 400 K;
    for n=1:30
        rho_am=0.5*rho*(1+T/Tw); % arithmetic mean density
        mu_w=mu*(Tw/T)^0.68125; % wall viscosity
        mu_am=(mu+mu_w)/2; % arithmetic mean viscosity
        hg=0.026/D_local^0.2 * Cp*mu^0.2/Pr^0.6 * (rho*v)^0.8 * rho_am/rho * (mu_am/mu_0)^0.2; % Bartz equation for heat transfer coefficient
        f=hg*(Tr-Tw)-(Tw-T_ext)/(t_phenol/k_phenol+t_carbon/k_carbon); % set convection minus conduction into ablative material (phenolic resin) equal to f; we want f to eventually be zero 
        dfdT=-hg-1/(t_phenol/k_phenol+t_carbon/k_carbon);
        Tw_new=Tw+0.5*(-f/dfdT);
        Tw_new=max(Tw_new,300);
        Tw_new=min(Tw_new,Tr);
        if abs(Tw_new-Tw)<0.1
            Tw=Tw_new;
            break
        end
        Tw=Tw_new;
    end
    % but wait, Tw is the temp of the ablative material, which burns off
    % calculate the temperature at the inside face of the carbon fiber
    q_cond=(Tw-T_ext)/(t_phenol/k_phenol+t_carbon/k_carbon); % conduction
    Tw_actual=Tw-(q_cond*(t_phenol/k_phenol)); % derived from fourier's law
end

% nozzle skirt material: R512E slurry silicide coated C-103 columbium
% run Tw function over various O/F ratios and supersonic area ratios
of_film=linspace(0.5,of,4); % iterate over O/F ratios up to main stream O/F ratio
eps_supson_Tw=linspace(1,50,100); % supersonic area ratios
D_local_Tw=sqrt(eps_supson_Tw).*Dt; % m, diameters
emiss_alloy=0.7; % lower bound value for alloy - https://niobium.tech/-/media/niobiumtech/attachments-biblioteca-tecnica/nt_thermomechanical-processing-and-properties-of-niobium-alloys.pdf
percent_fuel_film=0.05; % assumption of percent fuel used for film cooling
Tw=zeros(length(eps_supson_Tw),length(of_film));
ISP_true=zeros(length(eps_supson_Tw),length(of_film));
eta_film=zeros(length(eps_supson_Tw),length(of_film));
for i=1:length(of_film)
    for j=1:length(eps_supson_Tw)
        [Tw(i,j),ISP_true(i,j)]=Tw_calc(of_film(i),Pc,eps_supson_Tw(j),T_prop,D_local_Tw(j),emiss_alloy,percent_fuel_film,m_dot_MMH_engine,m_dot_NTO_engine);
    end
end
% chamber and nozzle material: carbon phenolic (ablative)
% run Tw_chamber function over various O/F ratios and subsonic area ratios
eps_subson_Tw=linspace(CR,1,100); % subsonic area ratios
D_local_Tw_chamber=sqrt(eps_subson_Tw).*Dt; % m, diameters
k_phenol=0.2; % W*m^-1*K^-1, thermal conductivity of phenolic - when phenolic resin chars, its k value should go up, however we will leave it as is for conservative design https://pmc.ncbi.nlm.nih.gov/articles/PMC9404460/#:~:text=The%20thermal%20conductivity%20and%20density%20of%20phenolic,about%200.2%20W/mK%20and%201.1%20g/cm3%2C%20respectively.
t_phenol=10*0.33*10^-3; % m, thickness of ten ply phenolic laminate on T-300 carbon fiber - https://www.azistausa.com/composites/prepregs/phenolic-prepregs#:~:text=%25%20(170%C2%B0C%2C%201,Strength%20(ILSS):%2044%20MPa.
k_carbon=9; % W*m^-1*K^-1, thermal conductivity of carbon fiber - https://dragonplate.com/carbon-fiber-thermal-conductivity-and-the-coefficient-of-thermal-expansion
t_carbon=0.001+t_phenol; % m, thickness of carbon fiber, assuming total thickness of 1 mm + 10 ply
% note: to find actual thickness, must look at ablation regression rate (cite https://pmc.ncbi.nlm.nih.gov/articles/PMC10003895/#:~:text=This%20indicates%20that%20increased%20surface%20ablation%20and,heat%20transfer%20to%20the%2030%C2%B0%20carbon%20phenolic),
% which depends on temperature and heat flux, among other variables.
% however, assume above thickness.
Tw_chamber=zeros(length(eps_subson_Tw),length(of_film));
for i=1:length(of_film)
    for j=1:length(eps_subson_Tw)
        Tw_chamber(i,j)=Tw_chamber_calc(of_film(i),Pc,eps_subson_Tw(j),T_prop,D_local_Tw_chamber(j),k_phenol,t_phenol,k_carbon,t_carbon,CR);
    end
end
% plot wall temp and ISP values vs. supersonic area ratios
colors=lines(length(of_film));
% Tw
figure; yyaxis left; hold on; grid on; grid minor
for i=1:length(of_film)
    legend_name=sprintf('Wall Temp (O/F = %.2f)',of_film(i));
    plot(eps_supson_Tw,Tw(i,:),'-','LineWidth',2,'Color',colors(i,:),'DisplayName',legend_name)
end
ylabel('Wall Temp (K)'); legend('show')
ax=gca; ax.YColor='k';
% add core temp values too
% since we are not iterating over the length of the nozzle, must calculate
% core temps
T_core=zeros(1,length(eps_supson_Tw));
for i=1:length(eps_supson_Tw)
    T_core(i)=CEA('problem','rocket','eql','o/f',of,'p,psia',Pc,'supson',eps_supson_Tw(i),'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','end').output.eql.temperature(3);
end
plot(eps_supson_Tw,T_core,'-','LineWidth',4,'Color','k','DisplayName','Core Temp');
% add failure temp to determine where film cooling ends
T_fail=(1500+273.15)-100; % K, 1500˚C survivability ceiling adjusted by 100 degrees for conservative calculations - https://www.sciencedirect.com/science/article/pii/S0257897224001403#:~:text=Introduction,components%20%5B4%2C5%5D.
yline(T_fail,'LineWidth',4,'Color','r','DisplayName','Failure Temp of Skirt');
% add vertical line to show where skirt begins
xline(20,'-','LineWidth',4,'Color','b','DisplayName','Skirt Begins');
% ISP_true
yyaxis right; hold on;
for i=1:length(of_film)
    legend_name=sprintf('ISP (O/F = %.2f)',of_film(i));
    plot(eps_supson_Tw,ISP_true(i,:),'--','LineWidth',2,'Color',colors(i,:),'DisplayName',legend_name)
end
ylabel('True Vacuum ISP (s)');
ax.YColor='k';
xlabel('Supersonic Area Ratios');
legend('show')
title('Film Cooling in Nozzle at Various Mixture Ratios');

% plot wall temp in chamber
figure; hold on; grid on; grid minor;
% Flip the x-axis direction
set(gca, 'XDir', 'reverse');
for i=1:length(of_film)
    legend_name=sprintf('Wall Temp (O/F = %.2f)',of_film(i));
    plot(eps_subson_Tw,Tw_chamber(i,:),'-','LineWidth',2,'Color',colors(i,:),'DisplayName',legend_name)
end
xlabel('Subsonic Area Ratios'); ylabel('Wall Temp (K)'); legend('show');
title('Film Cooling in Chamber at Various Mixture Ratios');

%based on graph analysis, will use film cooling/ablation at O/F=0.5 until area ratio
%of 20:1, where the nozzle skirt will be attached
ISP_current=ISP_true(1,100); % new true ISP value of engine!

T_gas_eff=CEA('problem','rocket','eql','o/f',0.5,'p,psia',Pc,'supson',20,'reactants','fuel','CH6N2(L)','wt%',100,'t,k',T_prop,'oxid','N2O4(L)','wt%',100,'t,k',T_prop,'output','gamma','transport','end').output.eql.temperature(3);

film_efficiency=(Tw(skirt_idx)-T_gas_eff)/(T_prop-T_gas_eff); % film cooling efficiency at nozzle skirt start

fprintf('\nNew Vacuum ISP = %.2f\n',ISP_current);
fprintf('Film Cooling Efficiency = %.2f percent\n',film_efficiency);
%% CAD Data Export
% CAD time! Export x and y data as MILLIMETERS
z_nozzle=zeros(length(x_nozzle),1); % column vector of zeros for SolidWorks z column
all_coords=[x_nozzle',y_nozzle_upper',z_nozzle]*1000; % to mm
unique_coords=unique(all_coords,'rows','stable');

nozzle_idx=find(eps_local>=20,1,'first'); % find coordinates up to area ratio 20 for nozzle
new_nozzle_idx=find(unique_coords(:,1)>=x_nozzle(nozzle_idx)*1000,1,'first'); % update index for unique coords

nozzle_coords=unique_coords(1:new_nozzle_idx,:); % up to area ratio 20
skirt_coords=unique_coords(new_nozzle_idx:end,:); % after area ratio 20
writematrix(nozzle_coords,'nozzle_coords.txt','Delimiter','tab');
writematrix(skirt_coords,'skirt_coords.txt','Delimiter','tab');
writematrix(unique_coords,'all_coords.txt','Delimiter','tab');

% P_local as function of x (from throat)
P_local_Pa=P_local*6894.76; % PASCALS
P_eqn=polyfit(x_nozzle(101:end),P_local_Pa(101:end),6); % higher order, the better
P_vals=polyval(P_eqn,x_nozzle);
plot(x_nozzle*1000,P_local_Pa,x_nozzle*1000,P_vals)

% P_local as function of x simple b/c SolidWorks keeps crashing.
P_start=600000;
P_end=1000;
x_end=x_nozzle(nozzle_idx)*1000+(-x_nozzle(1)*1000); % in mm, offset since the coordinates on CAD is set at the chamber beginning
coeff=(P_end-P_start)/x_end
fprintf('Function: P = %.3fx + %.0f\n',coeff,P_start)
%% 
k_eq=(t_carbon+t_phenol)/(t_carbon/k_carbon+t_phenol/k_phenol);