Home » 2012 » November

# Contents

## Free-Space Path Loss vs. Cable Loss

Click for more…! Consider a space link with 100 million kilometer distance and a transmit frequency of 2 GHz (S-Band). 1. Calculate the Free-Space Path Loss [dB].

```clear all;
clc;
close all;
F = 2;                                %transmit frequency (S-Band)[GHz]
d = 1e8;                              %[km]
L = 92.4 + 20*log10(F) + 20*log10(d); %[dB]
fprintf('Free-Space Path Loss L = %4.2f [dB] \n',L);
% Assume a typical coaxial cable with a loss of 0.3[dB/m] at S-Band.
% How long may this cable be for a total loss equal to the Free-Space
% Path Loss from above [m]?
c_l = 0.3;          %[dB/m]
l_cabel = L/c_l;    %[m]
fprintf('Equivalent Cabel Length = %4.2f [m] \n\n',l_cabel);```
```Free-Space Path Loss L = 258.42 [dB]
Equivalent Cabel Length = 861.40 [m]```

## Propagation Delays

1.Calculate the two-way propagation delays[min] between Earth and spacecrafts at diferent planets (from Mercury to Saturn; consider the following average distances between Sun-Planet: 0.3871AU for Mercury, 0.723AU for Venus, 1.524AU for Mars, 5.203AU for Jupiter, 9.582AU for Saturn). Assume conjunction Sun-Planet-Earth (Mercury, Venus), or Sun-Earth-Planet

```%(Mars,Jupiter, Saturn) for an easy calculation of the minimum distances
%(Note that this is not the worst case in terms of maximum distances to be
% considered for the actual link design).
v_light   = 300000;% Speed of light in vacuum [km/s]
R_mercury = 0.3871;    %[AU]
R_venus   = 0.7230;    %[AU]
R_earth   = 1.0000;    %[AU]
R_mars    = 1.5240;    %[AU]
R_jupiter = 5.2030;    %[AU]
R_saturn  = 9.5820;    %[AU]
AU        = 149597871; %[km]
mu_sun = 132712440018;          %(km^3s?2)```

Earth – Mercury

```p_d = 2*(R_earth - R_mercury)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Mercury = %4.2f [min]\n',p_d);
figure(1);
T_mercury  = 2*pi*((R_mercury*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:10*T_mercury;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_mercury = rem(time,T_mercury)/(T_mercury)*360;
phase = fi_earth - fi_mercury;
dist = (R_earth^2 + R_mercury^2 - 2*R_earth*R_mercury*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Mercury Comunication');
ylabel('Distance Earth - Mercury [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Mercury = 10.19 [min]`

Earth – Venus

```p_d = 2*(R_earth - R_venus)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Venus   = %4.2f  [min]\n',p_d);
figure(2);
T_venus  = 2*pi*((R_venus*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:7*T_venus;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_venus = rem(time,T_venus)/(T_venus)*360;
phase = fi_earth - fi_venus;
dist = (R_earth^2 + R_venus^2 - 2*R_earth*R_venus*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Venus Comunication');
ylabel('Distance Earth - Venus [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Venus   = 4.60  [min]`

Earth – Mars

```p_d = 2*(R_mars - R_earth)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Mars    = %4.2f  [min]\n',p_d);
figure(3);
T_mars  = 2*pi*((R_mars*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:6*T_mars;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_mars  = rem(time,T_mars)/(T_mars)*360;
phase = fi_earth - fi_mars;
dist = (R_earth^2 + R_mars^2 - 2*R_earth*R_mars*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Mars Comunication');
ylabel('Distance Earth - Mars [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Mars    = 8.71  [min]`

Earth – Jupiter

```p_d = 2*(R_jupiter- R_earth)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Jupiter    = %4.2f  [min]\n',p_d);
figure(4);
T_jupiter  = 2*pi*((R_mars*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:3*T_mars;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_jupiter = rem(time,T_jupiter)/(T_jupiter)*360;
phase = fi_earth - fi_jupiter;
dist = (R_earth^2 + R_jupiter^2 - 2*R_earth*R_jupiter*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title('Earth - Jupiter Comunication');
ylabel('Distance Earth - Jupiter [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Jupiter    = 69.86  [min]`

Earth – Saturn

```p_d = 2*(R_saturn- R_earth)/v_light*AU/60;
fprintf('Min Two-way Propagation Delay, Earth - Saturn    = %4.2f  [min]\n',p_d);
figure(5);
T_saturn  = 2*pi*((R_saturn*AU)^3/mu_sun)^0.5/86400;   % solar days
T_earth = 2*pi*((R_earth*AU)^3/mu_sun)^0.5/86400; % solar days
time = 0:1:3*T_mars;
fi_earth = rem(time,T_earth)/(T_earth)*360;
fi_saturn = rem(time,T_saturn)/(T_saturn)*360;
phase = fi_earth - fi_saturn;
dist = (R_earth^2 + R_saturn^2 - 2*R_earth*R_saturn*cosd(phase)).^0.5;
subplot(2,1,1);
plot(time/T_earth,dist);
title ('Earth - Saturn Comunication');
ylabel('Distance Earth - Saturn [AU]');
p_d = 2*dist/v_light*AU/60;
subplot(2,1,2);
plot(time/T_earth,p_d,'k');
ylabel('Two-way Propagation Delay [min]');
xlabel('Time [Earth Period]');```
`Min Two-way Propagation Delay, Earth - Saturn    = 142.65  [min]`

```%2. For the various cases calculate the Free-Space Path Loss[dB], assuming
% an RF frequency of 6 GHz.
F = 6; %transmit frequency [GHz]
% Earth - Mercury
L = 92.4 + 20*log10(F) + 20*log10((R_earth - R_mercury)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Mercury L = %4.2f [dB] \n',L);
% Earth - Venus
L = 92.4 + 20*log10(F) + 20*log10((R_earth - R_venus)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Venus   L = %4.2f  [dB] \n',L);
% Earth - Mars
L = 92.4 + 20*log10(F) + 20*log10((R_mars - R_earth)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Mars    L = %4.2f [dB] \n',L);
% Earth - Jupiter
L = 92.4 + 20*log10(F) + 20*log10((R_jupiter - R_earth)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Jupiter L = %4.2f [dB] \n',L);
% Earth - Saturn
L = 92.4 + 20*log10(F) + 20*log10((R_saturn - R_earth)*AU); %[dB]
fprintf('Free-Space Path Loss, Earth - Saturn  L = %4.2f [dB] \n\n',L);```
```Free-Space Path Loss, Earth - Mercury L = 267.21 [dB]
Free-Space Path Loss, Earth - Venus   L = 260.31  [dB]
Free-Space Path Loss, Earth - Mars    L = 265.85 [dB]
Free-Space Path Loss, Earth - Jupiter L = 283.93 [dB]
Free-Space Path Loss, Earth - Saturn  L = 290.13 [dB]```

## Satellite Design

Your satellite designer wants to reduce the satellite transmitter output power from 50 W to 25 W to save weight. How much is this reduction expressed in dB scale?

```P1 = 50;                % [W]
P1_db = 10*log10(P1);   % 17 dBW
P2 = 25;                % [W]
P2_db = 10*log10(P2);   % 14 dBW
red = (P2_db - P1_db);  % - 3  dBW

fprintf('Reduction expressed in DB scale %4.0f\n',red);

% If you want to maintain the satellite-to-ground
% station data link at the same data rate, you could achieve this by
% modifying the antenna on ground: By which factor do you have to increase
% the diameter of a your dish antenna then?
% Pt1*A1r = Pt2*A2r  A2r/A1r = Pt1/Pt2 = 2
% Pt - Transmited Power ;
% Ar - Effective Area of receiving antena;

R_factor = sqrt(P1/P2);
fprintf('Diameter increase your dish antenna %4.2f \n',R_factor);```
```Reduction expressed in DB scale   -3
Diameter increase your dish antenna 1.41```

Published with MATLAB® 7.10

## Hohmann transfer orbit

Find more at https://smallsat.wordpress.com

```A spececraft is in a 300 km circular eart orbit. Calculate
(a)The orbital delta v required for a Hohmann to a 3000 km coplanar circular
earth orbit and
(b)The transfer orbit time```
```clear all;
clc;
close all;
Alt_a = 300;
Alt_b = 3000;
Re = 6378;
R_a = Alt_a + Re;
R_b = Alt_b + Re;
mu = 398600;
V_a = (mu/R_a)^0.5;                       %Velocity of Circular Orbit at A
fprintf('Alt_a = %8.2f  Alt_b = %8.2f \n',Alt_a,Alt_b);
fprintf('Velocity of Circular Orbit at A  = %8.4f km/s \n',V_a);
a = (R_a+R_b)/2;                          %Eliptical Orbit Semimajor Axis
fprintf('Eliptical Orbit Semimajor Axis = %8.4f km/s \n',a);
Vtp = (2*mu*(1/R_a - 1/(R_a + R_b)))^0.5; %Transfer orbit at Perigee
fprintf('Transfer orbit at Perigee = %8.4f km/s \n',Vtp);
dVa = Vtp - V_a;
fprintf('dV burn at Perigee = %8.4f km/s \n',dVa);
Vtb = (2*mu*(1/R_b - 1/(R_a + R_b)))^0.5;
fprintf('Transfer orbit at Apogee = %8.4f km/s \n',Vtp);
V_b = (mu/R_b)^0.5;                       %Velocity of Circular Orbit at B
fprintf('Velocity of Circular Orbit at B  = %8.4f km/s \n',V_b);
dVb = V_b - Vtb;
fprintf('dV burn at Apogee = %8.4f km/s \n\n',dVb);
dVt = dVa + dVb;
fprintf('Total delta Velocity required  = %8.4f km/s \n',dVt);
T = 2*pi*(a^3/mu)^0.5;                          %Orbital period
fprintf('Transfer orbit time  = %d m %2.1f s \n\n\n',floor(T/120),rem(T/2,60));```
```Alt_a =   300.00  Alt_b =  3000.00
Velocity of Circular Orbit at A  =   7.7258 km/s
Eliptical Orbit Semimajor Axis = 8028.0000 km/s
Transfer orbit at Perigee =   8.3502 km/s
dV burn at Perigee =   0.6244 km/s
Transfer orbit at Apogee =   8.3502 km/s
Velocity of Circular Orbit at B  =   6.5195 km/s
dV burn at Apogee =   0.5734 km/s

Total delta Velocity required  =   1.1977 km/s
Transfer orbit time  = 59 m 39.3 s```

Assuming the orbits of earth and Mars are circular and coplanar, calculate (a) the time required for a Hohmann transfer from earth orbit to Mars orbit (b) the initial position of Mars ( ? ) in its orbit relative to earth for interception to occur.Radius of earth orbit  1.496  10^8 km. Radius of Mars orbit  2.279*10^8 km.Sun 1.327  10 11 km 3 /s 2 .

```R_earth = 1.496E8;
R_mars = 2.279E8;
mu = 1.327E11;
V_earth = (mu/R_earth)^0.5;           %Velocity of Circular Orbit of erath
fprintf('Velocity of Circular Orbit of Earth  = %8.4f km/s \n',V_earth);
a = (R_earth + R_mars)/2;                  %Eliptical Orbit Semimajor Axis
fprintf('Eliptical Orbit Semimajor Axis = %8.4f km/s \n',a);
Vtp =(2*mu*(1/R_earth - 1/(R_earth + R_mars)))^0.5; %Trans.orbit at Perigee
fprintf('Transfer orbit at Periapsis = %8.4f km/s \n',Vtp);
dVp = Vtp - V_earth;
fprintf('dV burn at Periapsis = %8.4f km/s \n',dVp);

Vta = (2*mu*(1/R_mars - 1/(R_mars + R_earth)))^0.5;
V_mars = (mu/R_mars)^0.5;           %Velocity of Circular Mars Orbit
fprintf('Velocity of Circular Mars Orbit = %8.4f km/s \n',V_mars);
dVa = V_mars - Vta;
fprintf('dV burn at Apoapsis = %8.4f km/s \n',dVa);

dVt = dVa + dVp;
fprintf('Total delta Velocity required  = %8.4f km/s \n',dVt);
T = 2*pi*(a^3/mu)^0.5;                          %Orbital period

% The initial position of Mars ( ? ) in its orbit relative to earth for
% interception to occur
T_mars = 2*pi*(R_mars^3/mu)^0.5;
angle_mars = pi/T_mars*T;
angle_mars_t =( pi - angle_mars);
fprintf('Angle Mars Traveled    = %4.3f deg \n\n',angle_mars*180/pi);

fprintf('Transfer time   = %4.3f days \n',T/(120*60*24));
fprintf('The initial position of Mars = %4.3f deg \n',angle_mars_t*180/pi);```
```Velocity of Circular Orbit of Earth  =  29.7831 km/s
Eliptical Orbit Semimajor Axis = 188750000.0000 km/s
Transfer orbit at Periapsis =  32.7264 km/s
dV burn at Periapsis =   2.9433 km/s
Velocity of Circular Mars Orbit =  24.1303 km/s
dV burn at Apoapsis =   2.6478 km/s
Total delta Velocity required  =   5.5911 km/s
Angle Mars Traveled    = 135.671 deg

Transfer time   = 258.840 days
The initial position of Mars = 44.329 deg```

Two geocentric elliptical orbits have common apse lines and their perigees are on the same side of the earth. The fi rst orbit has a perigee radius of r p  7000 km and e  0.3, whereas for the second orbit r p  32,000 km and e  0.5. (a) Find the minimum total delta-v and the time of flight for a transfer from the perigee of the inner orbit to the apogee of the outer orbit. (b) Do part (a) for a transfer from the apogee of the inner orbit to the perigee of the outer orbit.

```fprintf('Part A\n');
mu = 398600;
R1_p = 7000;                %km
e1 = 0.3;
R2_p = 32000;               %km
e2 = 0.5;
R1_a = R1_p*(e1 + 1)/(1 - e1);
h1 = (2*mu*R1_a*R1_p/(R1_a + R1_p))^0.5;

V1_p = h1/R1_p;             %Velocity at Perigee 1st orbit
V1_a = h1/R1_a;
R2_a = R2_p*(e2 + 1)/(1 - e2);
h2 = (2*mu*R2_a*R2_p/(R2_a + R2_p))^0.5;

V2_a = h2/R2_a;             %Velocity at Apogee 2nd orbit
V2_p = h2/R2_p;
x_dist = R2_a + R1_p;
R3_p = R1_p;
R3_a = x_dist - R3_p;
e3 = (R3_a - R3_p)/(R3_a + R3_p);
h3 = (2*mu*R3_a*R3_p/(R3_a + R3_p))^0.5;

V3_p = h3/R3_p;
V3_a = h3/R3_a;

dV13_p = V3_p - V1_p;
fprintf('\ndV burn at Perigee Orbit 1 = %8.4f km/s \n',dV13_p);
dV23_a = V2_a - V3_a;
fprintf('dV burn at Apogee Orbit 2 = %8.4f km/s \n',dV23_a);
dV_total = dV13_p + dV23_a;
fprintf('\nTotal Velocity Increase  = %8.4f km/s \n',dV_total);
a3 = (R3_p + R3_a)/2;
T3 = 2*pi*(a3^3/mu)^0.5;
fprintf('Transfer time   = %4.3f hours \n\n',T3/(120*60));

% Transfer from the apogee of the inner orbit to the perigee of the
% outer orbit.
fprintf('Part B');
R4_p = R1_a;
R4_a = R2_p;
e4 = (R4_p + R4_a)/2;
h4 =(2*mu*R4_a*R4_p/(R4_a + R4_p))^0.5;
V4_p = h4/R4_p;
V4_a = h4/R4_a;

dV14_ap = abs(V4_p - V1_a);
fprintf('\n\ndV burn at Apogee Orbit 1 = %8.4f km/s \n',dV14_ap);
dV42_ap = abs(V4_a - V2_p);
fprintf('dV burn at Perigee Orbit 2 = %8.4f km/s \n',dV42_ap);
dV_total = dV14_ap + dV42_ap;
fprintf('\nTotal Velocity Increase  = %8.4f km/s \n',dV_total);
a4 = (R1_a + R2_p)/2;
T4 = 2*pi*(a4^3/mu)^0.5;
fprintf('Transfer time   = %4.3f hours \n',T4/(120*60));```
```Part A

dV burn at Perigee Orbit 1 =   1.6989 km/s
dV burn at Apogee Orbit 2 =   0.6896 km/s

Total Velocity Increase  =   2.3885 km/s
Transfer time   = 16.154 hours

Part B

dV burn at Apogee Orbit 1 =   1.9708 km/s
dV burn at Perigee Orbit 2 =   1.6398 km/s

Total Velocity Increase  =   3.6106 km/s
Transfer time   = 4.665 hours```

Published with MATLAB® 7.10

## Best Fit Solution,Straight Line Fit

%https://smallsat.wordpress.com/

```clear all
clc
T = [0,20,100,100,200,400,400,800,1000,1200,1400,1600];
R = [10.96,10.72,14.1,14.85,17.9,25.4,26,40.3,47,52.7,58,63];
plot(T,R,'r*');
title('Straight Line Fit');
xlabel('Temperature (deg C)');
ylabel('Resistivity (Ohm cm)');
sy = 0;
sx = 0;
sxy = 0;
sx2 = 0;
n = 12;
for i=1:n
sy = sy +R(i);
sx = sx + T(i);
sxy = sxy + T(i)*R(i);
sx2 = sx2 + T(i)*T(i);
end
a = (sy*sx2-sx*sxy)/(n*sx2-sx^2);
b = (n*sxy-sx*sy)/(n*sx2-sx^2);
fprintf('\n\n Best Fit Coefficients')
fprintf('\n x = %6.4f T + %6.4f \n',b,a)
cm_fit =a+b*T;
hold on;
plot(T,cm_fit,'r');
legend('Data Points','Best Fit Solution','location','northwest');```
``` Best Fit Coefficients
x = 0.0337 T + 11.4695```

Published with MATLAB® 7.10

## Type K Thermocouple, Voltage vs Temperature

```% Polynomial Coefficients 0-500 °C[4]
% Type K thermocouple
c= [0,25.08355,7.860106e-2,-2.503131e-1,8.315270e-2,-1.228034e-2,...
9.804036e-4,-4.413030e-5,1.057734e-6,-1.052755e-8];
E= 0;
% Measured Voltage difference in mv
V = [1,E,E^2,E^3,E^4,E^5,E^6,E^7,E^8,E^9];
T_delta = c*V';
fprintf('Cold Bath Temp.  = %0.3f C\n' ,T_delta );
fprintf('Cold Bath Temp.  = %0.3f F\n\n' ,T_delta*1.8+32 );
E= 0.74;
% Measured Voltage difference in mv
V = [1,E,E^2,E^3,E^4,E^5,E^6,E^7,E^8,E^9];
T_delta = c*V';
fprintf('Room Bath Temp.  = %0.3f C\n' ,T_delta );
fprintf('Room Bath Temp.  = %0.3f F\n\n' ,T_delta*1.8+32 );

E= 3;
% Measured Voltage difference in mv
V = [1,E,E^2,E^3,E^4,E^5,E^6,E^7,E^8,E^9];
T_delta = c*V';
fprintf('Hot Bath Temp.  = %0.3f C\n' ,T_delta );
fprintf('Hot Bath Temp.  = %0.3f F\n\n' ,T_delta*1.8+32 );```
```Cold Bath Temp.  = 0.000 C
Cold Bath Temp.  = 32.000 F

Room Bath Temp.  = 18.526 C
Room Bath Temp.  = 65.346 F

Hot Bath Temp.  = 73.576 C
Hot Bath Temp.  = 164.436 F```

Published with MATLAB® 7.10

Contents

## SSME Nozzle

```clear all;
clc;
close all;
Pe = 16.8727;       %kPa
P0 = 20.4*1000;     %kPa
T0 = 3500;          %K

%Gas Properties
MW = 22;
gamma = 1.22;
%X Axis
X = [0.00 4.00 8.00 12.00 16.00 20.00 25.00 30.00 50.00 70.00 90.00 100.00...
120.00 140.00 160.00 180.00 200.00 220.00 240.00 260.00 280.00 300.00 305.00];
%Z Axis
Z = [17.50 15.50 13.00 12.25 13.00 15.50 18.50 22.00 33.00 41.50 50.50 54.50...
61.00 68.00 74.50 80.50 86.00 91.00 97.00 101.00 105.00 107.50 107.85];
%Diameter
Dia = 2*Z;

% Compute Mach Number at Exit
Me = (2/(gamma - 1)*(((P0/Pe)^((gamma - 1)/gamma)) - 1))^0.5;

fprintf('Mach Number at Exit = %8.6f \n',Me);
%Compute ratio of geometric throat area to Sonic throat area
AeAt = 77.5;                %SSME Nozzle geometry;
AtAs = (1/AeAt)*(1/Me*(2/(gamma + 1)*(1 + (gamma - 1)/2*Me^2))...
^((gamma + 1)/(2*(gamma - 1))));
fprintf('At/A* = %8.6f \n',AtAs);
% At/A* = 1, Nozzle is Choked

%Using  Iterative Solver at each Nozzle station
%Computing A/A* for each section
D_exit = 215.7;         %Exit Diameter
A_As = 77.5*(Dia/D_exit).^2;

oldMach = 4;
eps = 0.001;

for i =1:23
Mx(i) = Solver( oldMach,A_As(i),eps,gamma);
end
plot(X,Mx,'b');
hold on;
plot(X,Mx,'*r');
title('Mach Number Distribution in Nozzle');
xlabel('X axis [cm]');
ylabel('Mach Number');
grid on;

%Solving for Pressure Distribution in Nozzle

Px = P0./(( 1 + ((gamma - 1)/2)*Mx.*Mx).^(gamma/(gamma - 1)));
figure(2);
plot(X,Px,'b');
hold on;
plot(X,Px,'*r');
title('Pressure Distribution in Nozzle');
xlabel('X axis [cm]');
ylabel('Pressure[kPa]');
grid on;

%Pressure Distribution in Nozzle(normalized)
figure(3);
Pnx = 1./( 1 + ((gamma - 1)/2)*Mx.^2).^(gamma/(gamma - 1));
plot(X,Pnx,'b');
hold on;
plot(X,Pnx,'*r');
title('Pressure Distribution in Nozzle (normalized)');
xlabel('X axis [cm]');
ylabel('P/P_0');
grid on;

%Temperature Distribution in Nozzle
figure(4);
Tx = T0./( 1 + ((gamma - 1)/2)*Mx.^2);
plot(X,Tx,'b');
hold on;
plot(X,Tx,'*r');
title('Temperature Distribution in Nozzle');
xlabel('X axis [cm]');
ylabel('T  [ K ]');
grid on;

%Computing Exit Temperature, velocity based on Exit mach
T_exit = T0/( 1 + ((gamma - 1)/2)*Me^2);
fprintf('T_exit = %8.2f K\n',T_exit);

V_exit = (T_exit*gamma*8314.4126/MW)^0.5*Me;
fprintf('V_exit = %8.2f m/s \n',V_exit);```
```Mach Number at Exit = 4.858225
At/A* = 1.000152
T_exit =   973.23 K
V_exit =  3254.40 m/s```

function Mach = Solver( oldMach,A_As,eps,gamma)
maxit = 100000;
for i= 1: maxit
newMach = oldMach -(FOFM(oldMach,gamma)- A_As)/DFDM(oldMach,gamma);
oldMach = newMach;
fofm = FOFM(oldMach,gamma);

if( abs(fofm – A_As)/A_As <= eps)
%display(‘break’);
break;
end
end
Mach = oldMach;
end

function FOFM = FOFM(Mach,gamma)

FOFM = 1./Mach.*(2/(gamma + 1)*(1 + (gamma – 1)/2.*Mach.^2)).^…
((gamma + 1)/(2*(gamma – 1)));
end

function DFDM = DFDM(Mach,gamma)

DFDM = (2^((1 – 3*gamma)/(2 – 2*gamma)))*(Mach.^2 – 1)./…
(Mach.^2.*( 2 + Mach.^2*(gamma – 1))).*…
( 1 + (gamma – 1)/2*Mach.^2/(gamma + 1)).^((gamma + 1)/(2*(gamma -1)));
end