clc;
clear;

if(0)
q = 100;     % Propellant mass flow rate = 100 kg/s         
k = 1.20;    % Specific heat ratio = 1.20            
M = 24;      % Exhaust gas molecular weight = 24              
Tc = 3600;   % Combustion chamber temperature = 3600 K        
Pc = 5;      % Combustion chamber pressure = 5 MPa          
Pa = 0.05;   % Ambient pressure = 0.05 MPa
R = 8314.5;  % Universal gas constant N*m/kmol*K
Pe = Pa;     % Ambient air pressure
end

if(1)
q = 0.0917;  % Propellant mass flow rate kg/s         
k = 1.227;   % Specific heat ratio            
M = 20.9;    % Exhaust gas molecular weight              
Tc = 3280;   % Combustion chamber temperature K        
Pc = 2;      % Combustion chamber pressure MPa          
Pa = 0.1014; % Ambient pressure MPa
R = 8314.5;  % Universal gas constant N*m/kmol*K
Pe = Pa;     % Ambient air pressure
end

Pt = Pc * (1 + (k - 1) / 2)^(-k/(k-1)) * 10^6

Tt = Tc * (1 / (1 + (k - 1) / 2))

At = (q / Pt) * ((R * Tt) / (M * k))^0.5

Nm2 = (2 / (k - 1)) * ((Pc / Pa)^((k-1)/k) - 1)

Nm = Nm2^0.5

Ae = (At / Nm) * ((1 + (k - 1) / 2 * Nm2)/((k + 1) / 2))^((k+1)/(2*(k-1)))

Ve = ( (2 * k / (k - 1)) * (R * Tc / M) * (1 - (Pe / Pc)^((k-1)/k)) )^0.5

ExpansionRatio = Ae / At

F = q * Ve + (Pe - Pa) * Ae

ForceLBS = F / 4.44822162

Dt = (4*At/pi())^0.5

De = (4*Ae/pi())^0.5

%Isp = F / (q * 9.80665)

Lstar = 102;

theta = 30;

Vc = Lstar*At;

Lcapprox = exp(0.029*log(Dt)^2+.47*log(Dt)+1.94)

%Dc = ((Dt^3 + 24/pi*tand(theta)*Vc)/(Dc+6*tand(theta)*Lc))^.5

solve('Vc = (pi/24)*(6*Lc*Dc^2+(Dc^3-Dt^3)/tand(30))','Dc = ((Dt^3+24/pi*tand(30)*Vc)/(Dc+6*tand(30)*Lc))^.5')