%.................................................
%
%   Terminal Descent Engine Analysis
%
%   Type: Monopropellant Hydrazine Thruster
%         with bell nozzle
%
%.................................................

global game
global pe
global pa
global p0
global Arat
global Tav
global nang
global isp
global Cstar

global m_trans
global t_trans
global i_trans

%..Notes.........................................
%
%   Monopropellant Hyrazine Thruster
%   30% ammonia decomposition
%   Mol ratio leaving catalyst: 2 NH3, 1N2, 1H2
%   16 g/mol avg molar mass
%   255 isp
%   C* = 1350
%   1450K adiabatic reaction temperature
%
%   Ref: Rocket Propulsion Elements
%        George Sutton and Oscar Biblarz
%        7th edition
%        John Wiley and Sons, Inc.
%        New York
%        Chart: p.260
%
%   Note: Phoenix uses Hydrazine thruster for TD
%................................................

%..Calculate thrust coefficient

Ctau = sqrt(((2*game^2)/(game-1))*((2/(game+1))^((game+1)/(game-1)))*(1-(pe/p0)^((game-1)/game))+(pe/p0-pa/p0)*Arat);

%..Calculate nozzle throat and exit areas, m^2

At   = Tav/(Ctau*p0);

Ae   = At*Arat;

%..Calculate nozzle throat and exit diameters, m^2

Dt = 2*sqrt(At/pi);

De = 2*sqrt(Ae/pi);

%..Calculate nozzle length

L  = 0.7*(De/2 - Dt/2)/tan(nang);

%..Calculte const thrust massflow

mdot_c = Tav/(Ctau*Cstar);

%..Output data

format short g

disp(sprintf('-------------------------------------------------'))
disp(sprintf('Monopropellant Hydrazine Thruster Analysis'))
disp(sprintf('-------------------------------------------------'))
disp(sprintf('                                                 '))
disp(sprintf('Tank pressure       %g         atm',p0/101325))
disp(sprintf('Area ratio          %g           ',Arat))
disp(sprintf('Max Thrust          %g        N'  ,Tav))
disp(sprintf('Isp                 %g        sec',isp))
disp(sprintf('Char Velocity (C*)  %g       m/s',Cstar))
disp(sprintf('Thrust Coeff (Ct)   %g           ',Ctau))
disp(sprintf('Throat Diameter     %g  m^2',Dt))
disp(sprintf('Exit Diameter       %g  m^2',De))
disp(sprintf('Bell Nozzle Length  %g   m',L))
disp(sprintf('Const T mdot        %g   kg/s',mdot_c))
disp(sprintf('Var T mdot (AVG)    %g   kg/s', (m_trans(1)-m_trans(i_trans))/t_trans(i_trans) ))
disp(sprintf('-------------------------------------------------'))
Error using ==> mrdivide
Matrix dimensions must agree.