%a script for determining key liquid-fueled rocket motor parameters

%inits


clc
clear

R = 8314.5;         %universal gas constant N*m/kmol*K
g0 = 9.80665;       %gravitational attraction, ms^-2

q = .0917;            %mass flow, kg/s
k = 1.215;          %specific heat ratio, dimensionless
M = 22.0;           %exhaust gas molecular weight
MR = 2.35;              %mixture ratio
Tc = 3625;          %combustion chamber temperature, K
Pc = 12.67*10^6;     %combustion chamber pressure, Pa
Pa = 101325;         %ambient pressure, Pa

Isp = 360;          %specific impulse
Lstar = 120;        %combustion chamber design parameter
theta = 60;           %constriction half-angle
Lc_alt = 8;         %top-down chamber length, cm

s_y = 250*10^6;         %yield strength of structural steel
f_s = 2;                %safety factor for combustion chamber and other components

%maths

PcPSI = 1.45*10^-4 * Pc;
        %combustion chamber pressure in psi
Pt = Pc * (1 + (k - 1) / 2)^(-k/(k-1));
        %throat pressure
Tt = Tc * (1 / (1 + (k - 1) / 2));
        %throat temperature
At = (q / Pt) * ((R * Tt) / (M * k))^0.5;
        %throat area
Nm2 = (2 / (k - 1)) * ((Pc / Pa)^((k-1)/k) - 1);
        %mach number squared
Nm = Nm2^0.5;
        %mach number
Ae = (At / Nm) * ((1 + (k - 1) / 2 * Nm2)/((k + 1) / 2))^((k+1)/(2*(k-1)));
        %exhaust area
Ve = ( (2 * k / (k - 1)) * (R * Tc / M) * (1 - (Pa / Pc)^((k-1)/k)) )^0.5;
        %gas exit velocity
ER = Ae / At;
        %expansion ratio
F = Isp * q * g0;
        %developed force
FLBS = F / 4.44822162;
        %force in pounds
Dt = (4*At/pi())^0.5;
        %throat diameter
De = (4*Ae/pi())^0.5;
        %exhaust diameter
Decm = De * 100;
Dtcm = Dt * 100;
Vc = Lstar * At;
        %chamber volume
Vccm = Vc * 10^6;
        %chamber volume in cubic centimeters
Lccm = exp(0.029*log(Dt*100)^2+.47*log(Dt*100)+1.94);
        %approximate chamber length
Dccm = 10;
Dccmold = 0;
i = 0;
while abs(Dccm - Dccmold) > 1 * 10^(-6)
    Dccmold = Dccm;
    Dccm = ((Dtcm^3 + 24/pi*tand(theta)*Vccm)/(Dccmold+6*tand(theta)*Lccm))^.5;
    i = i + 1;
    if i > 50
        break
    end
end

j=0;

Dccm_override = 10;
Dccm_override_old = 0;

%manual chamber length override
while abs(Dccm_override - Dccm_override_old) > 1 * 10^(-6)
    Dccm_override_old = Dccm_override;
    Dccm_override = ((Dtcm^3 + 24/pi*tand(theta)*Vccm)/(Dccm_override_old+6*tand(theta)*Lc_alt))^.5;
    j = j + 1;
    if i > 50
        break
    end
end

Dc = Dccm * .01;

%disp('Dccm =')
%disp(Dccm)

%disp('Dccm_override = ')
%disp(Dccm_override)

%disp('Lcapprox = ')
%disp(Lccm)

%disp('Lc_alt = ')
%disp(Lc_alt)

%disp('Dtcm = ')
%disp(Dtcm)  

%disp('Decm = ')
%disp(Decm)


%a section to calculate wall thickness in the combustion chamber
s_all = s_y / f_s;
combustion_wall_thickness = (Pc - Pa) * Dc / s_all;

%disp('Combustion Chamber Wall Thickness, m')
%disp(combustion_wall_thickness)


%a section to calculate fuel and oxidizer orifice sizes and angles
%inits
Cd = 1;         %coefficient of discharge from hole. this is a lie and should be fixed ASAP.
rho_f = 819;     %kg/m^3
rho_o = 1146;      %kg/m^3
num_el = 4;     %number of injector elements

%calcs
q_f = q / (1 + MR);     %mass flow of fuel, kg/s
q_fi = q_f / num_el;
q_o = MR * q_f;         %mass flow of oxidizer, kg/s
q_oi = q_o / (num_el * 2); 
delta_P = Pc * 0.5;     %pressure drop across injector plate
A_f = (q_f / rho_f) / ((2 * delta_P/rho_f) ^ 0.5);  %total fuel injector hole area
A_o = (q_o / rho_o) / ((2 * delta_P/rho_o) ^ 0.5);  %total oxidizer injector hole area

A_fi = A_f/num_el;          %area of individual fuel holes
A_oi = A_o/(2 * num_el);    %area of individual oxidizer holes: twice that of fuel holes
D_fi = (4 * A_fi / pi) ^ .5;    %diameter of fuel holes
D_oi = (4 * A_oi / pi) ^ .5;    %diameter of oxidizer holes
plate_thickness = 80 * (D_fi + D_oi) / 2;

v_fi = q_f / (rho_f * A_fi);        %velocity of fuel leaving injector hole
v_oi = cosd(30) * q_o / (rho_o * A_oi);        %velocity of oxygen leaving injector hole

%gamma values are the angle between a stream and the chamber axis. we
%desire axial flow, ultimately.

gamma_f = 45;
gamma_o = asind((v_fi * sind(gamma_f))/(MR * v_oi));