%% FASE 2




%% 1. VULCAIN GEOMETRY


%% Berekenen Machgetal voor de Nozzle
%%

clear all
close all
clc

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

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

M = zeros(129, 1);
M(1) = 0.01;                            % Zelfgekozen beginwaarde van de Machgetal
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
Am = min(A);                            % Kleinste Oppervlakte A* in m

while ii < ii_max;
                                                  
q = 1;

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


while ii2 < ii_max2;
                                                  
o = 1;

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

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

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

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


%% Extra formules voor dimensieloos maken
%%

rhorel = roh/roh0;
prel = 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

[Vx,VA]=textread('Vulcaincorrected1.txt','%f %f');     

Vii = 1;
Vii_max = 33;
Vii2 = 34;
Vii_max2 = 129;

VM = zeros(129, 1);
VM(1) = 0.01;                           % Zelfgekozen beginwaarde van de Machgetal
VM(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
VAm = min(VA);                          % Kleinste Oppervlakte A* in m

while Vii < Vii_max;
                                                  
Vq = 1;

   while Vq < 5;
       VF(Vii) = VM(Vii)*(1+((y-1)/2)*VM(Vii)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((VAm)./VA(Vii));
                                        % Functie van het Machgetal (M)
       Vf(Vii) = (1+((y-1)/2)*VM(Vii)^2)^(-((y+1)/(2*(y-1))))*(1-((VM(Vii)^2*(y+1))/(2*(1+((y-1)/2)*VM(Vii)^2))));
                                        % Afgeleide van de functie F(M)
   VM(Vii) = VM(Vii)-(VF(Vii)/Vf(Vii));
 
   Vq = Vq+1;
   end
   
   VM(Vii+1) = VM(Vii); 
    
   Vii = Vii+1;
end


while Vii2 < Vii_max2;
                                                  
Vo = 1;

   while Vo < 5;
       VG(Vii2) = VM(Vii2)*(1+((y-1)/2)*VM(Vii2)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((VAm)./VA(Vii2));
                                        % Functie van het Machgetal (M)
       Vg(Vii2) = (1+((y-1)/2)*VM(Vii2)^2)^(-((y+1)/(2*(y-1))))*(1-((VM(Vii2)^2*(y+1))/(2*(1+((y-1)/2)*VM(Vii2)^2))));
                                        % Afgeleide van de functie F(M)
   VM(Vii2) = VM(Vii2)-(VG(Vii2)/Vg(Vii2));
 
   Vo = Vo+1;
   end
   
   VM(Vii2+1) = VM(Vii2); 
    
   Vii2 = Vii2+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)

VT = zeros(129, 1);
Vroh = zeros(129, 1);
Vp = zeros(129, 1);
Vu = zeros(129, 1);

VT = T0./(1+((y-1)/2).*VM.^2);                             % Temperatuur op positie x van de nozzle in K
Vroh = roh0./((1+((y-1)/2).*VM.^2).^(1/(y-1)));            % Dichtheid op positie x van de nozzle in Kg/(m^3)
Vp = p0./((1+((y-1)/2).*VM.^2).^(y/(y-1)));                % Druk op positie x van de nozzle in Pa
Vu = sqrt(((y*Rsp*T0).*VM.^2)./(1+((y-1)/2).*VM.^2));      % Snelheid op positie x van de nozzle in m/s

VAe = VA(129);                          % De oppervlakte aan de uitgang van de nozzle in m
Vue = Vu(129);                          % Snelheid aan de uitgang van de nozzle in m/s
Vpe = Vp(129);                          % Druk aan de uitgang van de nozzle in MPa
Vms = Vroh.*Vu.*VA;                     % Massastroom in Kg/s
Vmse = Vms(129);                        % Massastroom aan de uitgang van de nozzle in Kg/s
VStuw = Vmse*Vue+(Vpe-pp).*VAe;         % De Stuwkracht in N


%% Extra formules voor dimensieloos maken
%%

Vrhorel = Vroh/roh0;
Vprel = Vp/p0;




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

%% Berekenen Machgetal voor de Nozzle
%%

clear all
close all
clc

[V2x,V2A]=textread('Vulcaincorrected2.txt','%f %f');     

V2ii = 1;
V2ii_max = 33;
V2ii2 = 34;
V2ii_max2 = 129;

V2M = zeros(129, 1);
V2M(1) = 0.01;                          % Zelfgekozen beginwaarde van de Machgetal
V2M(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
V2Am = min(V2A);                        % Kleinste Oppervlakte A* in m

while V2ii < V2ii_max;
                                                  
V2q = 1;

   while V2q < 5;
       V2F(V2ii) = V2M(V2ii)*(1+((y-1)/2)*V2M(V2ii)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((V2Am)./V2A(V2ii));
                                        % Functie van het Machgetal (M)
       V2f(V2ii) = (1+((y-1)/2)*V2M(V2ii)^2)^(-((y+1)/(2*(y-1))))*(1-((V2M(V2ii)^2*(y+1))/(2*(1+((y-1)/2)*V2M(V2ii)^2))));
                                        % Afgeleide van de functie F(M)
   V2M(V2ii) = V2M(V2ii)-(V2F(V2ii)/V2f(V2ii));
 
   V2q = V2q+1;
   end
   
   V2M(V2ii+1) = V2M(V2ii); 
    
   V2ii = V2ii+1;
end


while V2ii2 < V2ii_max2;
                                                  
V2o = 1;

   while V2o < 5;
       V2G(V2ii2) = V2M(V2ii2)*(1+((y-1)/2)*V2M(V2ii2)^2)^(-(0.5*((y+1)/(y-1))))-(((y+1)/2)^(-(0.5*((y+1)/(y-1)))))*((V2Am)./V2A(V2ii2));
                                        % Functie van het Machgetal (M)
       V2g(V2ii2) = (1+((y-1)/2)*V2M(V2ii2)^2)^(-((y+1)/(2*(y-1))))*(1-((V2M(V2ii2)^2*(y+1))/(2*(1+((y-1)/2)*V2M(V2ii2)^2))));
                                        % Afgeleide van de functie F(M)
   V2M(V2ii2) = V2M(V2ii2)-(V2G(V2ii2)/V2g(V2ii2));
 
   V2o = V2o+1;
   end
   
   V2M(V2ii2+1) = V2M(V2ii2); 
    
   V2ii2 = V2ii2+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)

V2T = zeros(129, 1);
V2roh = zeros(129, 1);
V2p = zeros(129, 1);
V2u = zeros(129, 1);

V2T = T0./(1+((y-1)/2).*V2M.^2);                           % Temperatuur op positie x van de nozzle in K
V2roh = roh0./((1+((y-1)/2).*V2M.^2).^(1/(y-1)));          % Dichtheid op positie x van de nozzle in Kg/(m^3)
V2p = p0./((1+((y-1)/2).*V2M.^2).^(y/(y-1)));              % Druk op positie x van de nozzle in Pa
V2u = sqrt(((y*Rsp*T0).*V2M.^2)./(1+((y-1)/2).*V2M.^2));   % Snelheid op positie x van de nozzle in m/s

V2Ae = V2A(129);                        % De oppervlakte aan de uitgang van de nozzle in m
V2ue = V2u(129);                        % Snelheid aan de uitgang van de nozzle in m/s
V2pe = V2p(129);                        % Druk aan de uitgang van de nozzle in MPa
V2ms = V2roh.*V2u.*V2A;                 % Massastroom in Kg/s
V2mse = V2ms(129);                      % Massastroom aan de uitgang van de nozzle in Kg/s
V2Stuw = V2mse*V2ue+(V2pe-pp).*V2Ae;    % De Stuwkracht in N


%% Extra formules voor dimensieloos maken
%%

V2rhorel = V2roh/roh0;
V2prel = V2p/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(x,A, 'b')
plot(Vx, VA, 'g')
plot(V2x, V2A, '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')