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

%braeunig example 1.5
if(0)
syms q;     % Propellant mass flow rate = 100 kg/s         
k = 1.221;    % Specific heat ratio = 1.221            
M = 21.4;      % Exhaust gas molecular weight = 21.4              
Tc = 3470;   % Combustion chamber temperature =  K        
Pc = 50;      % Combustion chamber pressure = atm        
Pa = 1;   % Ambient pressure = 1 atm
R = 8314.5;  % Universal gas constant N*m/kmol*K
Pe = Pa;     % Ambient air pressure
end
% 1.5 validates

%braeunig 1.7
if(0)
q = 500;     % Propellant mass flow rate = 500 kg/s         
k = 1.221;    % Specific heat ratio = 1.221            
M = 21.4;      % Exhaust gas molecular weight = 21.4              
Tc = 3470;   % Combustion chamber temperature =  K        
Pc = 5.066;      % Combustion chamber pressure = Pa        
Pa = .1013;   % Ambient pressure = 1 atm
R = 8314.5;  % Universal gas constant N*m/kmol*K
Pe = Pa;     % Ambient air pressure
end
%1.7 validates

%braeunig 1.8
if(1)
q = 500;     % Propellant mass flow rate = 500 kg/s         
k = 1.221;    % Specific heat ratio = 1.221            
M = 21.4;      % Exhaust gas molecular weight = 21.4              
Tc = 3470;   % Combustion chamber temperature =  K        
Pc = 5.066;      % Combustion chamber pressure = Pa        
Pa = .0795;   % Ambient pressure = 1 atm
R = 8314.5;  % Universal gas constant N*m/kmol*K
Pe = Pa;     % Ambient air pressure
end
%1.8 validates

if(0)
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 = 3.447;      % 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 = 127;

%theta = 30;

%Vc = Lstar*At

%Lcapprox = exp(0.029*log(Dt*100)^2+.47*log(Dt*100)+1.94) 
    %in cm, suggested starting point for design from Braeunig
%Lcapproxm = Lcapprox*.01
    
%Dc = ((Dt^3 + 24/pi*tand(theta)*Vc)/(Dc+6*tand(theta)*Lc))^.5

%solve('Vc = (pi/24)*(6*Lcapproxm*Dc^2+(Dc^3-Dt^3)/tand(30))')
%solve('.0059 = (pi/24)*(6*.0616*x^2+(x^3-.0077^3)/tand(30))')
%'Dc = ((Dt^3+24/pi*tand(30)*Vc)/(Dc+6*tand(30)*Lc))^.5')