Space Communication
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
Source https://smallsat.wordpress.com/
% 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
Iterative Numerical Solution for Mach, P, and T distribution along the SSME Nozzle
Contents
Iterative Numerical Solution for Mach, P, and T distribution along the
SSME Nozzle
Source https://smallsat.wordpress.com
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
Small Satellites 