%% FASE 2




%% 1. VULCAIN GEOMETRY


%% Berekenen Machgetal voor de Nozzle
%%

clear all
close all
clc

[nozzle(1).x, nozzle(1).A]=textread('Vulcain.txt','%f %f');     

ii = 1;
ii_max = 33;
ii2 = 34;
ii_max2 = 129;

nozzle(1).M = zeros(129, 1);
nozzle(1).M(1) = 0.01;                            % Zelfgekozen beginwaarde van de Machgetal
nozzle(1).M(34) = 1.01;

mm = 0.018015;                          % Molaire massa van water(damp) in Kg/mol
Ralg = 8.314472;                        % Algemene Gasconstante in J/(mol*K)
Rsp = Ralg/mm;                          % Speciefieke Gasconstante J/(Kg*K)
n = 6;                                  % Aantal vrijheidsgraden voor waterdamp
y = (n+2)/n;                            % Gamma
nozzle(1).Am = min(nozzle(1).A);                            % Kleinste Oppervlakte A* in m

while ii < ii_max;
                                                  
q = 1;

   while q < 5;
       F(ii) = nozzle(1).M(ii)*(1+((y-1)/2)*nozzle(1).M(ii)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((nozzle(1).Am)./nozzle(1).A(ii));
                                        % Functie van het Machgetal (M)
       f(ii) = (1+((y-1)/2)*nozzle(1).M(ii)^2)^(-((y+1)/(2*(y-1))))*(1-((nozzle(1).M(ii)^2*(y+1))/(2*(1+((y-1)/2)*nozzle(1).M(ii)^2))));
                                        % Afgeleide van de functie F(M)
   nozzle(1).M(ii) = nozzle(1).M(ii)-(F(ii)/f(ii));
 
   q = q+1;
   end
   
   nozzle(1).M(ii+1) = nozzle(1).M(ii); 
    
   ii = ii+1;
end


while ii2 < ii_max2;
                                                  
o = 1;

   while o < 5;
       G(ii2) = nozzle(1).M(ii2)*(1+((y-1)/2)*nozzle(1).M(ii2)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((nozzle(1).Am)./nozzle(1).A(ii2));
                                        % Functie van het Machgetal (M)
       g(ii2) = (1+((y-1)/2)*nozzle(1).M(ii2)^2)^(-((y+1)/(2*(y-1))))*(1-((nozzle(1).M(ii2)^2*(y+1))/(2*(1+((y-1)/2)*nozzle(1).M(ii2)^2))));
                                        % Afgeleide van de functie F(M)
   nozzle(1).M(ii2) = nozzle(1).M(ii2)-(G(ii2)/g(ii2));
 
   o = o+1;
   end
   
   nozzle(1).M(ii2+1) = nozzle(1).M(ii2); 
    
   ii2 = ii2+1;
end


%% Kort deel over de druk tov de hoogte
%%

p00 = 100000;                           % Druk op aardoppervlak in Pa

for z = 1:500;                          % Hoogte in Km
pp(z) = p00*exp(-0.119*z);

end


%% Berekening van de stuwkracht
%%

p0 = 11500000;                          % De begindruk in Pa
T0 = 3523;                              % De begintemperatuur in K
roh0 = p0/(Rsp*T0);                     % De begindichheid in Kg/(m^3)

nozzle(1).T = zeros(129, 1);
nozzle(1).roh = zeros(129, 1);
nozzle(1).p = zeros(129, 1);
nozzle(1).u = zeros(129, 1);

nozzle(1).T = T0./(1+((y-1)/2).*nozzle(1).M.^2);                               % Temperatuur op positie x van de nozzle in K
nozzle(1).roh = roh0./((1+((y-1)/2).*nozzle(1).M.^2).^(1/(y-1)));              % Dichtheid op positie x van de nozzle in Kg/(m^3)
nozzle(1).p = p0./((1+((y-1)/2).*nozzle(1).M.^2).^(y/(y-1)));                  % Druk op positie x van de nozzle in Pa
nozzle(1).u = sqrt(((y*Rsp*T0).*nozzle(1).M.^2)./(1+((y-1)/2).*nozzle(1).M.^2));         % Snelheid op positie x van de nozzle in m/s

Ae = nozzle(1).A(129);                            % De oppervlakte aan de uitgang van de nozzle in m
ue = nozzle(1).u(129);                            % Snelheid aan de uitgang van de nozzle in m/s
pe = nozzle(1).p(129);                            % Druk aan de uitgang van de nozzle in MPa
nozzle(1).ms = nozzle(1).roh.*nozzle(1).u.*nozzle(1).A;                         % Massastroom in Kg/s
mse = nozzle(1).ms(129);                          % Massastroom aan de uitgang van de nozzle in Kg/s
nozzle(1).Stuw = mse*ue+(pe-pp).*Ae;              % De Stuwkracht in N


%% Extra formules voor dimensieloos maken
%%

nozzle(1).rhorel = nozzle(1).roh/roh0;
nozzle(1).prel = nozzle(1).p/p0;




%% 2. VULCAIN BOUNDARY LAYER CORRECTED 1 GEOMETRY
% Alle variablen krijgen een V daarvoor

%% Berekenen Machgetal voor de Nozzle
%%

clear all
close all
clc

[nozzle(2).x, nozzle(2).A]=textread('Vulcaincorrected1.txt','%f %f');     

ii = 1;
ii_max = 33;
ii2 = 34;
ii_max2 = 129;

nozzle(2).M = zeros(129, 1);
nozzle(2).M(1) = 0.01;                           % Zelfgekozen beginwaarde van de Machgetal
nozzle(2).M(34) = 1.01;

mm = 0.018015;                          % Molaire massa van water(damp) in Kg/mol
Ralg = 8.314472;                        % Algemene Gasconstante in J/(mol*K)
Rsp = Ralg/mm;                          % Speciefieke Gasconstante J/(Kg*K)
n = 6;                                  % Aantal vrijheidsgraden voor waterdamp
y = (n+2)/n;                            % Gamma
nozzle(2).Am = min(nozzle(2).A);                          % Kleinste Oppervlakte A* in m

while ii < ii_max;
                                                  
q = 1;

   while q < 5;
       F(ii) = nozzle(2).M(ii)*(1+((y-1)/2)*nozzle(2).M(ii)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((nozzle(2).Am)./nozzle(2).A(ii));
                                        % Functie van het Machgetal (M)
       f(ii) = (1+((y-1)/2)*nozzle(2).M(ii)^2)^(-((y+1)/(2*(y-1))))*(1-((nozzle(2).M(ii)^2*(y+1))/(2*(1+((y-1)/2)*nozzle(2).M(ii)^2))));
                                        % Afgeleide van de functie F(M)
   nozzle(2).M(ii) = nozzle(2).M(ii)-(F(ii)/f(ii));
 
   q = q+1;
   end
   
   nozzle(2).M(ii+1) = nozzle(2).M(ii); 
    
   ii = ii+1;
end


while ii2 < ii_max2;
                                                  
o = 1;

   while o < 5;
       G(ii2) = nozzle(2).M(ii2)*(1+((y-1)/2)*nozzle(2).M(ii2)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((nozzle(2).Am)./nozzle(2).A(ii2));
                                        % Functie van het Machgetal (M)
       g(ii2) = (1+((y-1)/2)*nozzle(2).M(ii2)^2)^(-((y+1)/(2*(y-1))))*(1-((nozzle(2).M(ii2)^2*(y+1))/(2*(1+((y-1)/2)*nozzle(2).M(ii2)^2))));
                                        % Afgeleide van de functie F(M)
   nozzle(2).M(ii2) = nozzle(2).M(ii2)-(G(ii2)/g(ii2));
 
   o = o+1;
   end
   
   nozzle(2).M(ii2+1) = nozzle(2).M(ii2); 
    
   ii2 = ii2+1;
end


%% Kort deel over de druk tov de hoogte
%%

p00 = 100000;                           % Druk op aardoppervlak in Pa

for z = 1:500;                          % Hoogte in Km
pp(z) = p00*exp(-0.119*z);

end


%% Berekening van de stuwkracht
%%

p0 = 11500000;                          % De begindruk in Pa
T0 = 3523;                              % De begintemperatuur in K
roh0 = p0/(Rsp*T0);                     % De begindichheid in Kg/(m^3)

nozzle(2).T = zeros(129, 1);
nozzle(2).roh = zeros(129, 1);
nozzle(2).p = zeros(129, 1);
nozzle(2).u = zeros(129, 1);

nozzle(2).T = T0./(1+((y-1)/2).*nozzle(2).M.^2);                             % Temperatuur op positie x van de nozzle in K
nozzle(2).roh = roh0./((1+((y-1)/2).*nozzle(2).M.^2).^(1/(y-1)));            % Dichtheid op positie x van de nozzle in Kg/(m^3)
nozzle(2).p = p0./((1+((y-1)/2).*nozzle(2).M.^2).^(y/(y-1)));                % Druk op positie x van de nozzle in Pa
nozzle(2).u = sqrt(((y*Rsp*T0).*nozzle(2).M.^2)./(1+((y-1)/2).*nozzle(2).M.^2));      % Snelheid op positie x van de nozzle in m/s

Ae = nozzle(2).A(129);                          % De oppervlakte aan de uitgang van de nozzle in m
ue = nozzle(2).u(129);                          % Snelheid aan de uitgang van de nozzle in m/s
pe = nozzle(2).p(129);                          % Druk aan de uitgang van de nozzle in MPa
nozzle(2).ms = nozzle(2).roh.*nozzle(2).u.*nozzle(2).A;                     % Massastroom in Kg/s
mse = nozzle(2).ms(129);                        % Massastroom aan de uitgang van de nozzle in Kg/s
nozzle(2).Stuw = mse*ue+(pe-pp).*Ae;         % De Stuwkracht in N


%% Extra formules voor dimensieloos maken
%%

nozzle(2).rhorel = nozzle(2).roh/roh0;
nozzle(2).prel = nozzle(2).p/p0;




%% 3. VULCAIN BOUNDARY LAYER CORRECTED 2 GEOMETRY
% Alle variablen krijgen een V2 daarvoor

%% Berekenen Machgetal voor de Nozzle
%%

clear all
close all
clc

[nozzle(3).x,nozzle(3).A]=textread('Vulcaincorrected2.txt','%f %f');     

ii = 1;
ii_max = 33;
ii2 = 34;
ii_max2 = 129;

nozzle(3).M = zeros(129, 1);
nozzle(3).M(1) = 0.01;                          % Zelfgekozen beginwaarde van de Machgetal
nozzle(3).M(34) = 1.01;

mm = 0.018015;                          % Molaire massa van water(damp) in Kg/mol
Ralg = 8.314472;                        % Algemene Gasconstante in J/(mol*K)
Rsp = Ralg/mm;                          % Speciefieke Gasconstante J/(Kg*K)
n = 6;                                  % Aantal vrijheidsgraden voor waterdamp
y = (n+2)/n;                            % Gamma
nozzle(3).Am = min(nozzle(3).A);                        % Kleinste Oppervlakte A* in m

while ii < ii_max;
                                                  
q = 1;

   while q < 5;
       F(ii) = nozzle(3).M(ii)*(1+((y-1)/2)*nozzle(3).M(ii)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((nozzle(3).Am)./nozzle(3).A(ii));
                                        % Functie van het Machgetal (M)
       f(ii) = (1+((y-1)/2)*nozzle(3).M(ii)^2)^(-((y+1)/(2*(y-1))))*(1-((nozzle(3).M(ii)^2*(y+1))/(2*(1+((y-1)/2)*nozzle(3).M(ii)^2))));
                                        % Afgeleide van de functie F(M)
   nozzle(3).M(ii) = nozzle(3).M(ii)-(F(ii)/f(ii));
 
   q = q+1;
   end
   
   nozzle(3).M(ii+1) = nozzle(3).M(ii); 
    
   ii = ii+1;
end


while ii2 < ii_max2;
                                                  
o = 1;

   while o < 5;
       G(ii2) = nozzle(3).M(ii2)*(1+((y-1)/2)*nozzle(3).M(ii2)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((nozzle(3).Am)./nozzle(3).A(ii2));
                                        % Functie van het Machgetal (M)
       g(ii2) = (1+((y-1)/2)*nozzle(3).M(ii2)^2)^(-((y+1)/(2*(y-1))))*(1-((nozzle(3).M(ii2)^2*(y+1))/(2*(1+((y-1)/2)*nozzle(3).M(ii2)^2))));
                                        % Afgeleide van de functie F(M)
   nozzle(3).M(ii2) = nozzle(3).M(ii2)-(G(ii2)/g(ii2));
 
   o = o+1;
   end
   
   nozzle(3).M(ii2+1) = nozzle(3).M(ii2); 
    
   ii2 = ii2+1;
end


%% Kort deel over de druk tov de hoogte
%%

p00 = 100000;                           % Druk op aardoppervlak in Pa

for z = 1:500;                          % Hoogte in Km
pp(z) = p00*exp(-0.119*z);

end


%% Berekening van de stuwkracht
%%

p0 = 11500000;                          % De begindruk in Pa
T0 = 3523;                              % De begintemperatuur in K
roh0 = p0/(Rsp*T0);                     % De begindichheid in Kg/(m^3)

nozzle(3).T = zeros(129, 1);
nozzle(3).roh = zeros(129, 1);
nozzle(3).p = zeros(129, 1);
nozzle(3).u = zeros(129, 1);

nozzle(3).T = T0./(1+((y-1)/2).*nozzle(3).M.^2);                           % Temperatuur op positie x van de nozzle in K
nozzle(3).roh = roh0./((1+((y-1)/2).*nozzle(3).M.^2).^(1/(y-1)));          % Dichtheid op positie x van de nozzle in Kg/(m^3)
nozzle(3).p = p0./((1+((y-1)/2).*nozzle(3).M.^2).^(y/(y-1)));              % Druk op positie x van de nozzle in Pa
nozzle(3).u = sqrt(((y*Rsp*T0).*nozzle(3).M.^2)./(1+((y-1)/2).*nozzle(3).M.^2));   % Snelheid op positie x van de nozzle in m/s

Ae = nozzle(3).A(129);                        % De oppervlakte aan de uitgang van de nozzle in m
ue = nozzle(3).u(129);                        % Snelheid aan de uitgang van de nozzle in m/s
pe = nozzle(3).p(129);                        % Druk aan de uitgang van de nozzle in MPa
nozzle(3).ms = nozzle(3).roh.*nozzle(3).u.*nozzle(3).A;                 % Massastroom in Kg/s
mse = nozzle(3).ms(129);                      % Massastroom aan de uitgang van de nozzle in Kg/s
nozzle(3).Stuw = mse*ue+(pe-pp).*Ae;    % De Stuwkracht in N


%% Extra formules voor dimensieloos maken
%%

nozzle(3).rhorel = nozzle(3).roh/roh0;
nozzle(3).prel = nozzle(3).p/p0;




%% GRAFIEKEN PLOTTEN VAN DE DRIE NOZZLE`S
% Vulcain Geometry                            --> blauw
% Vulcain Boundary Layer Corrected 2 Geometry --> groen
% Vulcain Boundary Layer Corrected 2 Geometry --> rood

%%

axis equal;
ylim([0 2.5]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('A [m^2]');
hold on;
plot(nozzle(1).x, nozzle(1).A, 'b')
plot(nozzle(2).x, nozzle(2).A, 'g')
plot(nozzle(3).x, nozzle(3).A, 'r')
figure;

ylim([0 1]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('Relatieve Dichtheid');
hold on;
plot(x,rhorel, 'b')
plot(Vx, Vrhorel, 'g')
plot(V2x, V2rhorel, 'r')
title('Relatieve Dichtheid');
figure;

ylim([0 1]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('Relatieve Druk');
hold on;
plot(x,prel, 'b')
plot(Vx,Vprel, 'g')
plot(V2x,V2prel, 'r')
title('Relatieve Druk');
figure;

ylim([0 6]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('M [-]');
hold on;
plot(x,M, 'b')
plot(Vx,VM, 'g')
plot(V2x,V2M, 'r')
figure;

ylim([500 3600]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('T [K]');
hold on;
plot(x,T, 'b')
plot(Vx,VT, 'g')
plot(V2x,V2T, 'r')
figure;

ylim([0 7]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('roh [Kg/(m^3)]');
hold on;
plot(x,roh, 'b')
plot(Vx,Vroh, 'g')
plot(V2x,V2roh, 'r')
figure;

ylim([0 12000000]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('p [Pa]');
hold on;
plot(x,p, 'b')
plot(Vx,Vp, 'g')
plot(V2x,V2p, 'r')
figure;

ylim([0 3500]);
xlim([-0.5 2.5]);
xlabel('x [m]');
ylabel('u [m/s]');
hold on;
plot(x,u, 'b')
plot(Vx,Vu, 'g')
plot(V2x,V2u, 'r')
figure;

ylim([900000 1200000]);
xlim([-1000 100000]);
xlabel('Druk op hoogte: p [Pa]');
ylabel('Stuwkracht [N]');
hold on;
plot(pp,Stuw, 'b')
plot(pp,VStuw, 'g')
plot(pp,V2Stuw, 'r')
figure;

xlim([900000 1000000]);
ylim([0 1000]);
ylabel('Hoogte z [Km]');
xlabel('Stuwkracht [N]');
hold on;
plot(Stuw, z, 'b')
plot(VStuw, z, 'g')
plot(V2Stuw, z, 'r')