%################################

%@@@@@@@@ Given information @@@@@

%Universal Gas Constant, J/(kmole. K)

%#################################
%@@@ The Propelllant properties @@
%#################################
Tc= 3370;
%The Flame Temperature (chamber temperature) (K)
Rho_p= 1772;
%The density of propellant ( kg/m^3)
M_w= 29.3;
%The molecular weight (g/mole)
R= 284;
%The exhaust Gas Constant (J/(kg.K))
Gamma= 1.17;
%Specific Heats Ratios of exhaust gases
n=0.4; k=4;
%Burning rate exponent and burnng rate constant
C_star= 1575;
%Characteristic Velocity (m/s)
%#############################################
%@@@ Properties of air at altitude of 25 m @@@
%#############################################
P0= 101026;
%Ambient Pressure (Pa)
T0 = 288.1;
%Ambient Temperature (K)
R_air = 286.9;
%Air Gas Constant (J/kg.K)
Rho_air= 1.222;
%the density of air (kg/m^3)
Gamma_air= 1.4;
%Specific Heats Ratios of air
%############################################
%@@@@ For the nozzle to be fully expanded @@@
%############################################
Pe = P0;
%An assumption for optimal solution
%###############################################
%@@@ Mission specifications and requirements @@@
%###############################################
M_No=0.85;
%Mach Number
minVelocity= 289.15;
% Minimum Velocity of flight (m/s)
minRange= 60000;
%Minimum target range (m)
%###############################################
%@@@@@@ Design Requirements and limitations @@@
%###############################################
maxMass= 800;
%Maximum Total mass of the missile (kg)
missilelength=6.0;
%Maximum length of the missile (m)
% given that:
warheadMass= 150;
guidaceMass= 40;

%The mass of the warhead (kg)

%###########################################
%@@@@@@@@@@ Methodoloy of designing @@@@@@@@
%###########################################
fprintf('\tNozzle Expansion Ratio Properties\n\t Me\t\t\t\t Epsilon\n');
%Epsilon is the expansion ratio (Epsolon= exitAera/ throatAera)
%For simplying the formulation: let:

a=Gamma-1;
b=Gamma+1;
for Me=1.2:0.05:3;
Epsilon = (((1+((a/2)*Me*Me))^(b/(2*a)))*((2/b)^(b/(2*a))))/Me;
fprintf('\t%0.2f \t\t\t \t%0.3f\n', Me, Epsilon);
end
%####################################################
%@@@@ chossing Initial value of Espilon @@@@@@@@@@@@@
%####################################################
Espilon= 7.440;
Me= 2.567;

%Corresponding Mach number to Espilon

fprintf('\nFor Expansion Ratio = %0.3f and Exit Mach Number = %0.3f: \n',
Epsilon, Me);
%##########################################################################
%@@@ Calculation of chamber pressure (Pc) for Me using isentropic formula
@@@@
%##########################################################################
Pc = P0*((1+(a/2)*Me*Me)^(Gamma/a));
fprintf('Chamber Pressure = %0.3f Pa\n', Pc);
%############################################
%@@@ Calculation of X-function
@@@@@@@@
%############################################
X_star= sqrt(Gamma)*(2/b)^(b/(2*a));
fprintf('X-function = %0.3f \n', X_star)
%#################################################
%@@@ Calculation of the coefficient of Thrust @@@@
%#################################################
Cf_optimum = sqrt((2*Gamma*Gamma/a)*((2/b)^(b/a))*(1-(Pe/Pc)^(a/Gamma)));
Cf_actual = 0.9*Cf_optimum;
fprintf('Thrust Coefficient = %0.3f \n', Cf_actual);
%####################################################
%@@@ Calculation of Specific Impulse(Isp) @@@@@@@@@@@
%####################################################
Isp= (C_star * Cf_actual)/g0;
fprintf('Specific Impulse = %0.5f s\n', Isp);
%##################################################
%@@@ Preliminary Design of Motor Dimensions
@@@@
%##################################################
D_motor= 0.83;
%Initial assumption of motor diameter (m)
L_charge= 2.63;
%Initial assumption of charge length (m)
t_casing_insulation= 0.025;
%Initial assumption of thickness of casing and
insulation

%charge diameter :
dc= D_motor-2.0*t_casing_insulation; %Determining diameter of charge
%Calculating burn aera (A_b):
A_b= 0.9*pi*dc*L_charge;
%Calculating throat aera and diameter :
A_star=((Rho_p*A_b*k*(1.0e-5)^n*sqrt(R*Tc))/(1000*X_star *(Pc)^(1-n)));
D_star= sqrt(4*A_star/pi);
fprintf('Throat Area = %0.5f m^2\nThroat Diameter = %0.5f m\n', A_star,
D_star);
%#############################################
%@@@@ Exhaust velocity determination @@@@@@@@@
%#############################################
Ve = sqrt((2*Gamma*R*Tc/a)*(1-(Pe/Pc)^(a/Gamma)));
fprintf('Exhaust Velocity = %0.5f m/s\n', Ve);
%###########################################
%@@@@ The Thrust (Fn) calculation @@@@@@@@@
%###########################################
Fn= Cf_actual* Pc * A_star;
fprintf('Thrust = %0.5f N or %0.5f kN\n', Fn, Fn/1000);
%@@@@ The mass flow rate (m_dot) calculation:
m_dot = Fn/Ve;
fprintf('Mass Flow Rate = %0.5f kg/s\n', m_dot);
%#######################################################
%@@@@ Calculation of the propellant mass (M_p) @@@@@@@@@
%#######################################################
%first calculating the graing volume (V0):
V0= pi*dc*dc*L_charge/4;
%Area of port (A_conduit) is to be at least four times the throat area (m^2:
A_conduit = 4*A_star;
D_conduit = sqrt(4*A_conduit/pi);
%volume of port (V_conduit) determination:
V_conduit = A_conduit*L_charge;
%The volume of grain (Vc):
Vc = V0 - V_conduit;
%The mass of propellant is:
M_p = Rho_p*Vc;
fprintf('\nPropellant volume = %0.5f m^3\n', Vc);
fprintf('\nPropellant Mass = %0.5f kg\n', M_p);
%########################################################################
%@@@@ The Burning rate (r_b)and burning time (t_b) calculation @@@@@@@@@
%########################################################################
r_b = k*((Pc*1.0e-5)^n)/1000;

t_b= M_p/m_dot;
fprintf('Burning Rate = %0.5f m/s\nburning time = %0.5f

%####################################################
%@@@ Calculation of Total Impulse(It) @@@@@@@@@@@
%####################################################
It=Fn*t_b;
fprintf('Total Impulse = %0.5f s\n', It);

s\n', r_b,t_b);